Plane-Wave Implementation and Performance of Ã -la-Carte Coulomb

Subscriber access provided by UNIV OF NEW ENGLAND ARMIDALE

Spectroscopy and Excited States

Plane-wave implementation and performance of à-la-carte Coulomb-attenuated exchange-correlation functionals for predicting optical excitation energies in some notorious cases Martin Peter Bircher, and Ursula Rothlisberger J. Chem. Theory Comput., Just Accepted Manuscript • Publication Date (Web): 16 Apr 2018 Downloaded from http://pubs.acs.org on April 16, 2018

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 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 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.

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 41 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

Plane-wave implementation and performance of à-la-carte Coulomb-attenuated exchange-correlation functionals for predicting optical excitation energies in some notorious cases Martin P. Bircher and Ursula Rothlisberger∗ Laboratoire de Chimie et Biochimie Computationnelles, Ecole polytechnique fédérale de Lausanne, Lausanne, Switzerland E-mail: [email protected]

Abstract Linear-response time-dependent density functional theory (LR-TDDFT) has become a valuable tool in the calculation of excited states of molecules of various sizes. However, standard generalised gradient approximation (GGA) and hybrid exchange-correlation (xc) functionals often fail to correctly predict charge-transfer (CT) excitations with low orbital overlap, thus limiting the scope of the method. The Coulomb-attenuation method (CAM) in the form of the CAM-B3LYP functional has been shown to reliably remedy this problem in many CT systems, making accurate predictions possible. However, in spite of a rather consistent performance across different orbital overlap regimes, some pitfalls remain. Here, we present a fully flexible and adaptable implementation of the CAM for Γ-point calculations within the plane-wave pseudopotential molecular

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

dynamics package CPMD and explore how customised xc functionals can improve the optical spectra of some notorious cases. We find that results obtained using plane waves agree well with those from all-electron calculations employing atom-centred bases, and that it is possible to construct a new Coulomb-attenuated xc functional based on simple considerations. We show that such a functional is able to outperform CAM-B3LYP in some notorious cases, while retaining similar accuracy in systems where CAM-B3LYP performs well.

1

Introduction

Be it for the vital conversion of sunlight to chemical energy in a leaf, for photochemical reactions causing harmful DNA damage to skin, for the blue fluorescence of scorpions, or for energy conversion in man-made solar cells: electronically excited states are of crucial importance to fundamental processes in Nature, and in scientific fields ranging from biology over chemistry to solid state physics. The theoretical description of the excitations which are at the base of these phenomena makes it possible to ultimately gain an improved understanding of these key events. A fully correlated description at the wave-function level is, unfortunately, prohibitively expensive for many, if not most systems of relevant size. It is thanks to the considerable progress in the field of Kohn-Sham time-dependent density functional theory (KS-TD-DFT) 1–3 that such excited state processes can these days be described at a comparably moderate computational cost. The application of the linear-response formulation of TD-DFT 4 is routinely used by many a computational chemist, and the ever increasing availability of computational resources has made it possible to describe larger and larger systems. However, like ground-state DFT and even more so, the choice of an appropriate exchange-correlation (xc) functional, or respectively xc kernel, is crucial, and often decides between results in good agreement with high-level reference data, and spectra which are considerably red-shifted and exhibit an incorrect ordering of states. 2


Page 2 of 41

Page 3 of 41 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 accuracy of the calculated transitions is mainly hampered by an inaccurate longrange decay of the xc potential, which leads to larger errors in the prediction of the excitation energies. 5 Hence, in analogy to their tendency to underestimate HOMO-LUMO gaps, generalised-gradient approximation (GGA) 6 xc functionals typically shift excitation energies by a considerable factor and may fail to recover the correct ordering of states. 7 Hybrid functionals, 8 which include a fixed proportion of exact exchange, may partly alleviate this problem for states in which there is an overlap between the orbitals involved in the transitions. However, even hybrid functionals will inevitably fail to describe low-overlap charge-transfer (CT) and Rydberg states, where the 1/r decay of the Coulomb operator is an important constituent in the correct description of the interaction between spatially distant orbitals. 9 But it is the inclusion of exact exchange that is most vital for these transitions: While an asymptotic correction of the (GGA) xc potential alone recovers the proper 1/r dependency and improves the description of Rydberg states, it cannot succesfully capture the effects of pronounced charge separation. 7 A promising remedy to this problem has been found in an appropriate splitting of the Coulomb operator, making it possible to ensure a correct decay of the xc potential at long-range, while keeping the accuracy and simplicity of a standard local formulation for the short-range components. The long-range correction (LC) scheme 5,10 and its generalisation, the Coulomb-attenuation method (CAM), 11 separate the Coulomb operator using an error function. While the short-range components are described using the GGA, the long-range components are taken into account using the exact exchange operator. This splitting captures the essentials of charge-transfer transitions and considerably increases the accuracy for such states. LC- and CAM-based xc functionals have therefore become a standard tool for the calculation of molecular excitation energies with pronounced charge-transfer character within TD-DFT. 12 The correspondingly modified Coulomb operator is easily implemented in a Gaussian basis. CAM-B3LYP is not only the first functional that was constructed using the Coulombattenuation method, it has also become the probably most prominent and abundantly used

3



functional of this family. Peach et al. 9 have assessed its performance on a diverse test set of main-group molecules shortly after the original paper presenting CAM-B3LYP, and many succesful applications of the functional followed: Among the substantial improvements documented with respect to B3LYP 8,13 or PBE, 14 Peach et al. have found that CAM-B3LYP accurately describes excitations in (poly)acenes including naphthalene, yielding correct state ordering. They could also show that the spectra of model dipeptides improve substantially with respect to non-Coulomb attenuated functionals. 7 Subsequently, other studies have found that the use of CAM-B3LYP can predict spectra of AT nucleobase dimers 15 and indoles 16 of varying substitution patterns with increased accuracy. All these systems exhibit chargetransfer or Rydberg transitions between spatially separated states, which results in the inevitable failure of GGA and conventional hybrid functionals. The orbital overlap in a transition can be quantitatively characterised by the parameter Λ introduced by Peach et al.; 9 a small value of Λ indicating a small overlap between the involved orbitals, a large value a substantial one. In their study covering excitations of different character, they have found that both GGA and conventional hybrids suffer from an inconsistent performance over the complete range of Λ values: For overlap values Λ < 0.4 in the case of PBE and Λ < 0.3 in the case of B3LYP, the errors in the excitation energies become substantially bigger. For CAM-B3LYP, no such correlation was found over the whole range 0 < Λ < 1. CAM-B3LYP was shown to fare particularly well for charge transfer excitations, especially in the ‘low overlap’ regime, Λ < 0.3; although some cases with charge transfer character in the ‘intermediate Λ regime’ are accurately described, too (e.g. the retinal protonated Schiff base 17,18 ). However, in systems where there is significant overlap between the orbitals involved in the transition, conventional hybrid functionals such as PBE0 19,20 often fare better, and CAM-B3LYP tends to red-shift the excitation energies. This is notably the case for the doubly fluorescent dye DMABN, for which CAM-B3LYP tends to overestimate the excitation energies of the S2 CT state with Λ = 0.72. 9 Fully longrange corrected functionals 21 such as LC-BLYP 5 or LC-PBE0 22 were not evaluated in the

4


Page 4 of 41

Page 5 of 41 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


study by Peach et al., but they have since been shown to perform well in certain systems with very low overlap that cannot be accurately described with CAM-B3LYP. 15,16,23 This is attributed to the absence of any GGA exchange at longest range, which benefits the description of Rydberg states, 9,21 but comes at the cost of an inferior average performance especially for local excitations. 11,16,21 Despite the absence of correlation between the predicted excitation energies and their Λ-values, there exist some systems even in the low-overlap range (Λ < 0.3) for which CAMB3LYP fails to deliver an accurate description. In the case of p-nitroaniline, the excitation energies are reasonably predicted, but solvatochromic shifts cannot be reproduced since the ground- and excited state dipoles are grossly over- and underestimated, respectively. 24 Whereas CAM-B3LYP predicts a correct state ordering and reasonable energetics for the excitations in the AT base pair, 15 the HOMO is predicted to lie on the wrong base when compared to higher-level wavefunction methods and basic considerations based on the ionization potential of the isolated bases. 25 Seemingly reasonable results may therefore be obtained based on the wrong physical reason. Similarly, the ordering of close-lying, low energy excitations may be inversed in some systems; this is the case for the β-dipeptide model system introduced by Serrano-Andrés et al. 26 and subsequently popularised in the aforementioned benchmark set by Peach et al. 9 Other studies have found the same problem to occur in the case of 7-azaindole, 16 even though the state ordering for other substituted indoles could be correctly predicted. A wrong ordering of states may be especially detrimental for excited state molecular dynamics, 27 where the forces excerted on the nuclei may substantially differ between the two swapped states, leading to a quantitatively as well as qualitatively wrong propagation of the system. For systems such as the β-dipeptide and 7-azaindole, 16 the use of an LC functional may yield a qualitatively correct ordering of the low-lying excitations, but the energies often remain too low. Alternatively, the range separation parameter µ in CAM-B3LYP may be tuned in order to ameliorate the performance of the functional. 28 This process known as γ-tuning adjusts the range separation parameter to a value that accurately

5



Page 6 of 41

reproduces ionization potentials. This approach often permits for an accurate description of the excitations of interest, 16 but it constitutes a rather system-specific remedy, lacking portability and thus predictive power in comparing various systems. All of these notorious systems are included in figure 2 and have been chosen here to test the performance of an ‘à-la-carte’ combination of range-separated functionals.

dipeptide

DMABN

naphthalene

β-dipeptide

p-nitroaniline

7-azaindole

retinal protonated Schiff base

AT base pair

Figure 1: Molecules explicitly discussed in this study, all of which contain excitations that are difficult to describe in TD-DFT when using a GGA or hybrid functional as the xc kernel. In other applications, a splitting of the Coulomb operator opposite to the LC and CAM scheme may be beneficial. This has been proposed in screened hybrid functionals for solid state applications, 29 where the exact exchange is limited to short- and the GGA exchange to long range. Functionals such as HSE03 29 yield results superior to those obtained with the GGA for many systems. Screened hybrids are especially beneficial in combination with a plane wave/pseudopotential approach, since they conveniently eliminate the divergence of the Coulomb operator at the G = 0 component of the plane wave basis in sufficiently large 6


Page 7 of 41 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


simulation cells. Since the splitting adopted in the LC and CAM schemes does not offer any particular advantage in solid state applications and does not eliminate this problematic divergence term, these methods have received much less attention in plane wave codes. In the following, we present the Coulomb-attenuation method applied to Γ-point calculations in a plane-wave pseudopotential framework. The implementation allows for the CAM to be combined with any exchange functional of choice, offering maximum flexibility. This makes the construction of customised ‘à-la-carte’ Coulomb-attenuated xc functionals possible, which can be tailored to any system of choice, thus maximising the performance of the method. In order to gain maximum flexibility, all CAM parameters can be chosen freely. Our implementation of the CAM in the molecular dynamics package CPMD 30 targets applications where Γ-point sampling is routinely used, and makes the sampling of large systems possible via the fully Hamiltonian QM/MM-scheme implemented in the CPMD code. Simulations of charge-transfer systems in the gas phase as well as in condensed matter therefore become feasible using a plane wave/pseudopotential approach. To facilitate the calculation of the necessary terms, a new driver for the calculation of the exchange-correlation energy has been implemented in the CPMD code. The paper is organised as follows: First, we give a short summary of the Coulombattenuation method, followed by a description of the implementation. We then give a more detailled account of the test systems used to benchmark both the implementation and a new Coulomb-attenuated xc functional. We discuss the performance of the CAM in plane waves with respect to the choice of pseudopotential, and in comparison to all-electron calculations with atom-centred basis sets. Finally, we will show how a flexible choice of the underlying GGA exchange-functional can improve accuracy in systems where standard functionals yield unsatisfactory results by comparing a customised ‘à-la-carte’ CAM-xc functional constructed based on simple considerations to the well-established CAM-B3LYP.

7



2

Page 8 of 41

Theory

2.1

The Coulomb-Attenuation Method

ˆ is split into two domains In the Coulomb-attenuation method, the Coulomb operator W dominated by long-range (lr) and short-range (sr) components each:

ˆ =W ˆ sr + W ˆ lr W 0 ˆ sr |ri = 1 − [α + β erf(µ |r − r |)] hr0 | W |r − r0 | 0 ˆ lr |ri = α + β erf (µ |r − r |) hr0 | W |r − r0 |

(1) (2) (3)

where α, β and µ are adjustable parameters 11 and α = 0 and β = 1 in the original LC method. 5,10 The first term is treated using a GGA expression for the exchange functional and becomes smaller for larger Coulomb distances, whereas the second term grows with increasing |r − r0 | and is treated using Fock’s expression for the exchange energy. The effective Coulomb operator in the exchange integrals then becomes:

ExHFX

1 XXX = 2 σ i j

ZZ

α ∗ ∗ ψiσ (r)ψjσ (r0 )

+ β erf(µ |r − r0 |) ψjσ (r)ψiσ (r0 )drdr0 0 |r − r |

(4)

Correspondingly, the GGA enhancement factor has to be adapted to the screened Coulomb operator. The adaptation is based on the LDA for a short-range Coulomb operator and apropriately generalised: 5

ExGGA

√ 1 X 4/3 8 1 = ρ K σ 1 − α − β aσ πerf + 2aσ (bσ − cσ ) 2 σ σ 3 2aσ

(5)

where Kσ is the spin-dependent formulation of the exchange enhancement factor and the

8


Page 9 of 41 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


terms due to the correction read: 1/2

µKσ aσ = √ 1/3 6 πρσ 1 bσ = exp − 2 − 1 4aσ 1 cσ = 2a2σ bσ + 2

2.2

(6) (7) (8)

ˆ in reciprocal space Singularity of the Coulomb operator W

At the Γ-point, the exact exchange energy can be efficiently computed with a mixed real-space reciprocal-space algorithm 31 after introducing a resolution of identity in G and rearranging the terms due to the complex conjugate of the Kohn-Sham (KS) orbitals (where hr|iσ i = hiσ |ri at the Γ-point):

ExHFX

=

XXXZ σ

=

i

i

0

0

0

0

0

Z

dr hjσ |r i hr |iσ i

ˆ |ri hiσ |ri hr|jσ i dr hr0 | W

(9)

j≥i

XXXZ σ

0


Z

0

0

0

dG hr |G i

Z

ˆ |Gi dG hG | W 0

Z dr hG|ri hiσ |ri hr|jσ i

j≥i

(10) Here, i, j index Kohn-Sham orbitals. The matrix elements of the Coulomb operator in ˆ |Gi = reciprocal space, hG0 | W

1 4π δ(G − G0 ), Ω G2

exhibit an integrable divergence at G = 0. 32

In practice, the plane wave/pseudopotential formalism relies on a discrete representation of points in direct and reciprocal space, and the integrals become sums associated to discrete Fourier transforms. The divergence term becomes problematic in this discrete form, and the

9



Page 10 of 41

ˆ Coulomb operator has to be replaced by a suitable generalisation Φ:

ExHFX

=

XXXZ σ

i

0

0

0


Z

0

0

Z

0

dG hr |G i

ˆ |Gi dG hG | Φ 0

Z dr hG|ri hiσ |ri hr|jσ i

j≥i

(11) ˆ |Gi, the offending divergence is screened by a suitIn the generalised matrix element hG0 | Φ able function χ:     1 4π δ(G − G0 ) for G 6= 0 2 0 ˆ hG | Φ |Gi = Ω G   χ(0)δ(G0 ) for G = 0

(12)

ˆ like W, ˆ is diagonal in G. where Ω denotes the supercell volume and Φ, In the initial treatment proposed by Gygi and Baldereschi, 32 the screening term χ is obtained numerically. An auxiliary function which exhibits the same singularity as the problematic term is added to and subtracted from the Coulomb operator, and the screening is given by the difference of the discrete representation of the auxiliary function as a sum over G and its analytical integral over a continuous range Q. Due to their particular choice of χ, the approach could not be applied to Γ-point sampling due to its poor convergence with respect to the number of k-points and simulation supercell size. In the following, we base our treatment on the scheme subsequently developed by Broqvist et al . 33 which - in contrast to the initial approach by Gygi and Baldereschi - selects an auxiliary function f (Q) that converges rapidly towards 1/Q2 : 2

e−γQ f (Q) = Q2

(13)

The G = 0 term is then given by the residual difference between the integral and discrete

10


Page 11 of 41 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


sum over the auxiliary function: 1 χ(0) = 2 2π

Z

4π X f (G) Ω G6=0 Q # " 2 4π X e−γG 1 , = lim √ − γ→0 πγ Ω G G2 f (Q)dQ −

(14)

(15)

where the second term (the discrete sum) is obtained numerically. This results in a more rapid convergence with respect to the size of the supercell, and the scheme can therefore be applied in Γ-point only calculations. Although a non-divergent analytical description for the G = 0 component is found for screened hybrids (it is simply π/Ω2 µ2 ), its accuracy depends on the size of the periodic supercell. If µ is small, the potential does not decay sufficiently rapidly with respect to the number of G-vectors and a treatment analogous to the non-screened Coulomb operator has to be used. 33 Broqvist et al have therefore suggested to treat the short range exact exchange analogously to the full Coulomb operator:  h i    1 4π 1 − e−G2 /4µ2 δ(G − G0 ) for G 6= 0 ˆ sr |Gi = Ω G2 hG0 | Φ   0 χ(µ)δ(G ˜ ) for G = 0

(16)

Using the same auxiliary function as for the full Coulomb potential, the resulting G = 0 term for the short-range screened exchange then reads: χ(µ) ˜ = χ(0) − χ

1 4µ2

(17)

where µ is the range separation parameter. Based on this treatment, the singularity correction for the Coulomb-attenuation method is easily found by using the identity erf(x) + erfc(x) = 1 and introducing the parameters α

11



Page 12 of 41

and β. We write for the long-range components:  h i    1 4π α + βe−G2 /4µ2 δ(G − G0 ) for G 6= 0 ˆ lr |Gi = Ω G2 hG0 | Φ   χ(µ, ¯ α, β)δ(G0 ) for G = 0

(18)

where the G = 0 term is simply a sum of the terms due to the full and the screened Coulomb potential, weighted by the attenuation parameters α and β. χ(µ, ¯ α, β) = αχ(0) + βχ

3

1 4µ2

(19)

Implementation

Some exchange-correlation functionals make use of a linear combination of different GGA exchange contributions, such as XLYP 34 (using 72.2% B88 35 and 34.7% PW91 36 exchange) or the well-known B3LYP (using 80% LDA, 72% of the B88 gradient correction term and 20% exact exchange). Accordingly, the usage of the CAM does not have to be intrinsically limited to a single type of exchange functional (as it is the case for the most prominent CAMB3LYP and LC-PBE0, where B88 and PBEx are attenuated, respectively). With respect to all possible combinations, our implementation achieves maximum flexibility in the choice of exchange functional by writing the exchange enhancement factor as a sum over individual contributions:

Kσx =

N X

cf Kσf

(20)

f

where Kσf denotes any GGA exchange-functional and cf are the corresponding weights of a P x total of N exchange functionals; N f cf = 1. The attenuation is then applied to Kσ after all N contributions have been added up. This makes it possible to consistently apply the CAM to any arbitrary combination of LDA and GGA exchange functionals, creating a custom

12


Page 13 of 41 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


‘à-la-carte’ xc functional. Three terms are needed for the propagation of a Kohn-Sham wavefunction with a GGA description: Z 1X E[ρ] = − η[Kσx ]Kσx [ρσ , ∇ρσ ]ρ4/3 σ (r)dr 2 σ Z 1X [ρσ (r)]dr =− 2 σ δ[ρ] δρσ δ[ρ] V2x = δ |∇ρσ | V1x =

(21) (22) (23) (24)

where we have introduced η as the attenuation function. For functionals that are derived for the closed-shell case or use a different definition of the enhancement factor, e.g. E[ρ] = R ρ(r)F [ρ]dr, the spin-dependent exchange enhancement factor Kσx is easily obtained from x using the spin-scaling relations. F [ρ] or Kαβ

All derivatives can be efficiently calculated by making extensive use of the chain rule and by transiently storing frequently used terms (notably δKx /δρ, δKx /δ∇ρ) on the stack. In order to make further performance gains, certain powers of the density and the gradient (ρ, ρ4/3 , ρ1/3 , |∇ρ|, ∇ρ2 ) are precomputed on a per-grid-point basis and reused in the calculation of Kσf and the attenuation function, thus avoiding repetitive, unnecessary operations. The implementation makes use of procedure pointers in order to facilitate the choice of functional. At the beginning of every run, the procedure pointer denoting Kσf is set to the selected exchange functional, no (explicit) if-construct is therefore necessary when looping over grid points. The new algorithm reaches an asymptotic speedup of 20% when calculating the LDA and GGA contributions with respect to the previous standard implementation, which obtained LDA and GGA contributions for each functional from seperate drivers and used extensive if-constructs to select among the functionals at every individual gridpoint.

13



Page 14 of 41

SUBROUTINE calculate_xc_energy(rho,grad,E_x) real, dimension(:), & intent(in) :: rho, grad real, intent(out) :: E_xc ! ! Stores rho, rho**1/3, rho**4/3, |\nabla rho|, |\nabla rho|**2 ! type(storage_t) :: storage ! ! Stores K, dK/drho, dK/d_nabla_rho, epsilon_x ! type(functional_t) :: functional forall p in gridpts CALL storage%precalculate_reused_terms_from(rho(p),grad(p)) forall f in K_x CALL functional%K(f)%calculate_energy_and_derivatives(storage,functional) if (functional_is_attenuated) CALL functional%attenuation(storage,functional) epsilon_x(p) = functional%epsilon_x E_x = sum(epsilon_x(:)) END SUBROUTINE calculate_xc_energy

14


Page 15 of 41 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

Test Set and Computational Details

4.1

Description of the Test Set

Basic convergence tests with respect to basis set and supercell size were carried out on a single, isolated water molecule. A test set for excitation energies was assembled comprising the molecules mentioned before and depicted in figure 2. The set includes the original test systems by Peach et al., 9 augmented with other systems where range-separation is known to be of importance: p-nitroaniline, 7-azaindole, the AT nucleobase dimer and the retinal protonated Schiff base. For the sake of comparison with localised basis sets, 5 of those systems were selected. The model dipeptide represents a typical system where CAM-B3LYP can successfully be applied: Peach et al found that a GGA (PBE) yields an inaccurate state ordering of the charge-transfer excitations, whereas B3LYP and CAM-B3LYP reproduce the ordering of the reference values. However, B3LYP is known to underestimate the energies of the n1 → π2∗ and π1 → π2∗ charge transfer excitations by up to 1.7eV, which are both reproduced by CAM-B3LYP with a reasonable accuracy of about 0.2eV. The spectrum of p-nitroaniline has also been reported to be reasonably predicted using CAM-B3LYP, 24 and although the charge separation in the first CT state was reported to be overestimated, we have used pnitroaniline as another probe known to benefit from range-separation. Naphthalene served as the most simple example of an acene and another notable example of the influence of range separation: Whereas the 1 B3u and 1 B3u states are inverted when using both PBE and B3LYP, only CAM-B3LYP recovers the correct state ordering. However, the reported excitation energies for the two lowest optically allowed singlet transitions deviate from the reference by about +0.21eV for 1 B3u and −0.16eV for 1 B2u , which results in a considerable underestimation of the state separation. The separation is reported to improve for larger acenes (anthracene, tetracene etc.); we have therefore chosen naphthalene as the most critical and sensitive compound to assess our implementation. An example where CAM-B3LYP

15



(narrowly) fails to deliver a quantitatively correct description of states is given by 7-azaindole. Whereas a comprehensive study of various substituted indoles found CAM-B3LYP to be in good agreement with reference values, in the case of 7-azaindole, 16 the 1 La and 1 Lb states have been reported to be swapped. Analogously to naphthalene, the states also lie much too close in energy, but are now also incorrectly ordered. This molecule can therefore serve as a representative of excitations with CT character and low orbital overlap where CAM-B3LYP surprisingly fails. The final molecule used in the selected test set is DMABN, with a Λ value of 0.72 for the S2 CT state it is a typical usage case for a non-range separated hybrid or even a simple GGA. 9 Indeed, whereas the errors for the 1 A and 1 B states are smaller than 0.2eV in conjecture with B3LYP, it is about doubled when using CAM-B3LYP. It serves as an example of an excitation with considerable overlap, thus completing the range of excitations covered here. The chosen test suite therefore includes both systems that are well described using CAM-B3LYP, as well as some notorious cases. In order to assess a possible gain of accuracy by using a customised functional rather than the established CAM-B3LYP, three additional systems were studied. The retinal protonated Schiff base constitutes a system where CAM-B3LYP has been succesfully used to predict both simple optical 17,37 and two-photon absorption spectra, 38,39 improving over conventional hybrids. In an investigation on the GC and AT nucleobase dimers, 15 it was found that CAMB3LYP can predict accurate excitation energies for the AT base pair, but the functional localises the HOMO and the LUMO on the wrong moieties, respectively. 25 According to the ionization potential of the isolated base, the HOMO should be localized on adenine, but it is predicted to lie on thymine. The AT base pair can therefore serve as a probe for the correct orbital localization obtained with a given functional. Similarly as in 7-azaindole, in the β-dipeptide studied by Peach et al, all PBE, B3LYP and CAM-B3LYP fail to describe the ordering of the π1 → π2∗ and n1 → π2∗ transitions, with errors being larger than 0.75eV for CAM-B3LYP and reaching a maximum value of about 4.5eV when using PBE. The βdipeptide therefore serves as yet another example in which CAM-B3LYP does not work even

16


Page 16 of 41

Page 17 of 41 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 a quantitative level, probably constituting the hardest test case for new range-separated functionals. Finally, in order to verify a possible overall gain of accuracy, the MAE of both CAM-B3LYP and CAM-O3LYP is compared for the complete test set by Peach et al. 9

4.2

Computational Setup

Calculations using Gaussian basis sets were either carried out using DALTON 2016 40 (pnitroaniline, the retinal protonated Schiff base and naphthalene) or the Gaussian09 suite of programs 41 (the remaining molecules) using Dunning’s correlation-consistent basis sets. 42 All calculations employing a Slater basis set 43 were carried out using ADF. 44–46 The structures for the dipeptides, DMABN and naphthalene were taken from the database published by Peach et al.. The structure from the retinal protonated Schiff base is an arbitrarily chosen snapshot extracted from molecular dynamics simulations at 300K, whereas the structure of p-nitroaniline was optimised using the aug-cc-pVTZ basis and the B3LYP xc functional. Structural optimisations for the remaining molecules were performed using the Gaussian suite of programs, following the published protocols of existing benchmarks for the AT base pair 15 and and 7-azaindole. 16 Excitation energies, where not quoted from the literature, were calculated using the Tamm-Dancoff approximation to TD-DFT 47 and the cc-pVDZ, aug-cc-pVDZ and d-aug-cc-pVDZ basis sets, respectively. The new exchange-correlation driver was implemented in a development version of CPMD 30 (successor of version 4.1). The KS orbitals were expanded in plane waves contained in an orthorhombic supercell of varying dimensions and using either Martins-Troullier (MT) 48 or Goedecker-Teter-Hutter (GTH) 49 pseudisation of the atomic core orbitals (the respective values for the energy cutoff are given in the results section; the supercell size for every system is available in the supporting information). Following standard practice for hybrid functionals, BLYP pseudopotentials were used for CAM-B3LYP calculations, and OLYP pseudopotentials were used for calculations with CAM-O3LYP. The density was expanded with a 4 times greater cutoff value than the one adpoted for the orbitals. The Poisson equations for the 17



isolated systems were solved using the algorithm of Martyna and Tuckerman. 50 All calculations made use of an atomic wavefunction initialization using distributed Lanczos, 51 the new distributed linear algebra algorithm by Bekas and Curioni 52 and the ‘new’ exact exchange driver by Weber et al.; 31 the cutoff in the calculation of the Fock exchange energy were not changed with respect to the standard values for orbitals and density.

5

Results and Discussion

5.1 5.1.1

Convergence of Eigenvalues HOMO-LUMO gaps in a plane wave/pseudopotential basis

When implementing the singularity correction for (screened) hybrid functionals, Broqvist et al. 33 have also assessed the convergence of the KS-HOMO-LUMO gaps with respect to both the energy cutoff and the size of the supercell. In the following, we shall present a similar assessment on an isolated water molecule. Considering that the LUMO is very diffuse in this specific case, the HOMO-LUMO gap appears to be a sensitive measure of convergence. When assessing the convergence with respect to the size of the simulation supercell, two scenarios have to be distinguished: In a periodic setup, the gap for an isolated system can only be reproduced when the molecule at the centre of the cell is sufficiently far apart from its periodic images. When the Poisson equations of the replicas are decoupled (and the requirements of the used Poisson solver appropriately met 50 ), the gap will converge with respect to the lowest G-vector components, which corresponds to increasingly longerrange components in real space as the supercell size increases. This may be especially important if the LUMO is very diffuse (we note that an unbound continuous state will only be appropriately described if the length of the simulation cell l = ∞). Hence, a change in cutoff value enhances the accuracy of the description by adding more rapidly oscillating, short-range components; the maximum ‘diffuseness’ allowed is essentially governed by the

18


Page 18 of 41

Page 19 of 41 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


choice of l for the simulation cell. Table 1 shows the gap obtained using the CAM-B3LYP xc functional for an isolated system contained in varying sizes of the simulation supercell using hard GTH pseudopotentials. Values for a fully periodic system are also given. The corresponding values for softer MT pseudopotentials are tabulated in the SI. Table 1: HOMO-LUMO gap of a water molecule calculated using CAM-B3LYP and GTH pseudopotentials at increasing energy cutoffs and varying cubic simulation cell lengths l using periodic or isolated system boundary conditions for solving the Poisson equation. Ecut [Ry] 70 80 100 120 150 180

∆i [eV], isolated 10 Å 20 Å 10.080 10.361 10.073 10.420 10.150 10.490 10.259 10.531 10.345 10.560 10.383 10.574

system 30 Å 10.399 10.447 10.510 10.548 10.577 10.591

∆i [eV], periodic system 10 Å 20 Å 30 Å 10.105 10.362 10.399 10.100 10.420 10.448 10.177 10.491 10.509 10.284 10.532 10.548 10.367 10.561 10.577 10.404 10.575 10.591

The gaps show convergence at 150Ry or all systems. For a small cubic simulation supercell (l = 10Å), chosing a lower cutoff value of 100Ry introduces a substantial error of 0.23eV. A notable error is still present at 100Ry even for the two larger simulation cells, but it becomes less relevant for practical purposes, since the maximal deviation of < 0.1eV lies below the typical accuracy of the functional itself. Whereas the values in the smallest of the supercells still have an error of about 0.1eV at a cutoff of 120Ry, the corresponding values have converged in the 20Å and 30Å supercells, with errors being lower than 0.05eV, and full convergence is reached at 150Ry for all of the three supercells considered. The convergence behaviour is analogous to the one observed for GGA or standard hybrid functionals once the simulation cell is of sufficient size: Changes in the gap are still substantial when increasing the length l of the cubic simulation cell from 10 Å to 20 Å, with changes in the converged gap of 0.2eV. The change in gap is within the usual numerical tolerance (< 0.05eV) for a further extension to 30 Å, emphasising again the importance of an appropriately large cell for the correct description of the system’s LUMO. 19



The results for the fully periodic system (where the Poisson equations are not decoupled) show the same trend, with the only relevant difference with respect to the isolated system occuring in the 10 Å box. Still, these differences are lower than those observed when enlarging the supercell. Given the trends observed for the isolated system, this is most likely attributed to spurious interactions between periodic images at this intermolecular distance. The influence of these interactions on the gap supports the view that the requests on the Tuckerman-Martyna Poisson solver are not yet met either, since the simulation cell must span at least twice the spatial extent of the charge density. The strong changes in gaps when increasing the cutoff within the small simulation cell is hence due to an insufficient cell size for both the isolated and periodic system, resulting in an incorrect description of the electron density. The same considerations hold for the gaps obtained with the softer MT pseudopotentials (cf. the SI), albeit convergence is achieved at lower cutoff values. The maximum deviation in the converged gaps with respect to GTH pseudisation is very small, ∆∆ = 0.02eV, illustrating that the influence of the pseudisation of the cusp on the gaps is negligibly low. These results confirm that the use of the CAM in a plane wave/pseudopotential formalism does not introduce any convergence issues or additional restrictions, given that the size of the simulation supercell is chosen sufficiently large in order to accomodate the whole spectrum of G-vectors that are required to span the range of the µ-dependent switching function.

5.1.2

Comparison with atom-centred basis sets

Plane waves inherently contain diffuse components even at a low cutoff energy, and the character of the diffuse functions is restricted only by the size of the supercell, and convergence is reached with respect to high-frequency components needed to describe the region around the pseudised core. In turn, the number of short-range components and the description of the orbitals around the nuclei can be systematically improved by using pseudopotentials of increasing hardness (which requires a simultaneous increase in the cutoff value). The com-

20


Page 20 of 41

Page 21 of 41 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


plete basis set limit can therefore be systematically reached (once completeness holds with respect to the adopted pseudopotential). The situation presents itself substantially different in atom-centred basis sets. In a Gaussian basis, the orbitals around the nuclei are well described, but for an accurate description of most molecular properties, the compact basis usually has to be enhanced by augmentation with long-range functions. When the KS orbitals are expanded in Gaussians, a single set of diffuse functions is often sufficient for routine applications. (A more realistic description of the cusp and a correct decay of the basis at long range can be obtained by resorting to a Slater-type basis.) A comparison between plane wave and atom-centred basis set calculations can therefore reveal the influence of the longest-range components (described well within plane waves) and the explicit description of the orbitals around the nuclear cusp (reproduced well using Gaussian functions). In the following, we will compare the HOMO-LUMO gaps of the preceeding section with the corresponding results obtained from various Gaussian basis sets of increasing accuracy. Additional tests were performed using a Slater-type basis in order to obtain a systematic analysis with respect to the cusp condition. The results are shown in table 2. When increasing the the size of the basis from double to hextuple zeta, changes in the gap are considerable for the non-augmented basis sets, spanning a range of a total of 0.9eV. Once a single diffuse function is included, the gaps become much more uniform, with a difference between aug-cc-pVDZ and aug-cc-pVTZ of only 25 meV. When increasing zeta to ζ = 6, the gap fluctuates within a negligible range of 3 meV. Adding more diffuse functions never changes the gap by more than 5 meV. The aug-cc-pVDZ basis can therefore be considered sufficiently accurate when calculating HOMO-LUMO gaps with CAM-B3LYP. The trend is very similar for a Slater-type basis, where we have only included single-zeta values for comparison. As for the Gaussian basis set, the omission of diffuse functions leads to an insufficient description of the gap. However, as soon as a single set of diffuse functions

21



Page 22 of 41

is included, the gap again converges rapidly. The change from augmented double to triple zeta is only about 8meV. Even an augmented single zeta basis seems to yield a suprisingly accurate gap, with the difference being only 0.08eV with respect to the augmented triple zeta basis. While the difference between the converged gaps in plane waves and the corresponding values in a Gaussian basis is about 0.1eV, it is only 0.03eV when compared to a Slater type basis (augmented triple zeta vs. GTH/180 Ry/30 Å cell). Even though the orbitals are pseudised around the core, the plane wave/pseudopotential approach yields results which are virtually indistinguishable from all-electron calculations with an atom-centred basis. The slightly larger deviation with respect to the Gaussian basis may be attributed to differences in the long-range decay and the description of the cusp, but they still lie well within what is usually deemed chemical accuracy. Table 2: HOMO-LUMO gap of a water molecule calculated using atom-centred basis sets augmented with a varying number of diffuse functions. Augmentation fcts. cc-pVDZ cc-pVTZ cc-pVQZ cc-pV5Z cc-pV6Z SZ DZ TZ

5.2

0 1 Gaussian basis 9.832 10.679 10.387 10.702 10.576 10.705 10.691 10.705 10.714 10.704 Slater basis 17.300 10.711 12.318 10.631 11.310 10.639

2

3

10.667 10.696 10.670 10.700 10.700

10.666 10.695 10.700 10.700 10.700

Excitation energies

The most frequent use of Coulomb-attenuated functionals is the description of excited states (which is influenced by the accuracy of the KS eigenvalues examined in the previous section through the linear response equations). We therefore conclude the assessment of our implementation of the CAM in plane waves by comparing the results to excited state energies 22


Page 23 of 41 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


obtained using Gaussian bases. In line with the trends observed in the preceeding section, only the cc-pVDZ basis was considered, which was augmented with a different number of diffuse functions. The results are depicted in table 3. The excitation energies were computed for the set of molecules introduced above, which apart from DMABN all contain transitions with charge-transfer or Rydberg character. Since the plane wave/pseudopotential implementation of LR-TDDFT in the CPMD code is limited to the Tamm-Dancoff approximation (TDA), the electronic spectra in a Gaussian basis were obtained within the same approximation for the sake of comparison. While it has been found that high-overlap singlet transitions can be substantially affected by the use of the TDA, it has also been reported that the results obtained from the TDA compare more favourably to high-level reference values than those obtained from full TDDFT. 53 This effect can be understood in terms of the triplet stability measure, the corresponding values for a subset of the molecules considered here have been reported and discussed in reference 53. For the dipeptide, DMABN and 7-azaindole, the results of the singly-augmented Gaussian basis sets already exhibit a negligible difference to those obtained with the plane wave/pseudopotential approach. Neither the character nor the energetics of the transitions do change when using a doubly augmented Gaussian basis. Both Gaussians and plane waves yield the same ordering of states and differences in the excitation energies, which are smaller than 0.05eV, i.e. they lie within a range that we have previously considered as converged. The small deviations may be attributed to both the pseudisation of the orbitals in plane waves, as well as the limited spatial extent of a localised atom-centred basis, along with its predefined decay properties. The situation is different in naphthalene and p-nitroaniline. In naphthalene, the S1 state dominated by a HOMO-1 → LUMO transition is only predicted using plane waves or the singly augmented atom-centred Gaussian basis. However, the energetics of the S2 state (mainly HOMO → LUMO) is consistent between the singly and doubly augmented Gaussian basis. The excitation energies agree well between our plane wave implementation

23



Page 24 of 41

Table 3: Excitation energies of the first two to four excited singlet states obtained for various systems and basis sets using LR-TDDFT and the CAM-B3LYP xc functional within the Tamm-Dancoff approximation. Values for a customised CAM-O3LYP functional are indicated, too. Literature reference values from high-level wavefunction calculations are given where available. In p-nitroaniline, the first forbidden transition is not predicted in a plane-wave basis when using CAM-B3LYP; its oscillator strength is f = 0 in the other calculations. The character of the transition is indicated whenever it is not consistent within different methods. Reference values for naphthalene, DMABN and the dipeptide correspond to those adopted by Peach et al. 9 State cc-pVDZ Naphthalene 54 S1 – S2 4.73 DMABN 55 S1 4.91 S2 5.29 7−azaindole 56 S1 5.19 1 La 5.24 1 Lb S2 S3 5.36 26 Dipeptide S1 5.67 S2 5.92 S3 7.09 S4 7.33 p-nitroaniline 57 S1 4.02 S2 4.57 S3 – S4 4.79

CAM-B3LYP aug-cc-pVDZ d-aug-cc-pVDZ

MT/80 Ry

CAM-O3LYP MT/80 Ry

Ref.

4.66 4.86

– 4.86

4.65 4.79

4.63 4.88

4.46 4.88

4.76 4.98

4.76 4.96

4.77 4.95

4.68 4.81

4.25 4.56

5.06 1 La 5.08 1 Lb 5.32

5.05 1 La 5.08 1 Lb 5.32

5.03 1 La 5.07 1 Lb 5.33

5.05 1 Lb 5.09 1 La 5.35

4.22 1 Lb 4.49 1 La 5.27

5.67 5.89 6.24 6.55

5.67 5.88 6.24 6.51

5.72 5.93 6.20 6.44

5.71 5.91 6.16 6.38

5.62 5.79 7.18 8.07

4.00 4.53 4.57 4.87

4.00 4.53 – 4.87

– 4.54 – 4.87

4.19 4.59 4.69 4.92

4.12 4.66 4.68 4.75

24


Page 25 of 41 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 aug-cc-pVDZ. The basis-set dependency issue becomes more involved in the case of p-nitroaniline. Whereas the S2 π → π ∗ transition is predicted in all cases, the S3 n → π ∗ state at 4.68eV only appears when one single diffuse function is used. Once a further set of diffuse functions is included, the S3 n → π ∗ transition disappears again. It can therefore be concluded that in the case of p-nitroaniline, the description of the S3 state using CAM-B3LYP is particularly sensitive with respect to the choice of basis. Since the S3 transition disappears once the basis is enlarged, this indicates that the (quantitatively correct) prediction using aug-cc-pVDZ is a mere artefact, and that CAM-B3LYP is not properly able to describe the excitations in the limit of a complete basis. This view is further supported by the results obtained in plane waves, where the S3 state is absent: The ordering of the states obtained at a cutoff of 120Ry coincides with the one obtained using a doubly augmented Gaussian basis. It has to be noted that the S1 state has an oscillator strength f = 0, and that it is not predicted when using plane waves. Overall, the results obtained using a plane wave pseudopotential approach are in excellent agreement with the ones obtained with all-electron calculations in a Gaussian basis set. If any of the orbitals included in the transitions of interest is highly diffuse, plane waves fare better than a singly augmented Gaussian basis. There is no indication that the presence (or absence) of pseudisation has any (relevant) influence on the spectra, with the remaining differences between the doubly augmented Gaussian basis and the plane waves being vanishingly small. Although the excitation energies usually do not change considerably when approaching the diffuse limit, the character and number of states may, which can be important for applications such as excited state dynamics. The observations made in p-nitroaniline further stress the importance of a sufficiently large basis, since seemingly correct predictions may be an artefact due to an incomplete basis. Only the use of a very diffuse Gaussian basis set or plane waves reveal that CAM-B3LYP does not properly predict one of the transitions. Given that most standard applications of Coulomb-attenuated functionals use only a single set of

25



augmentation functions, plane waves hence offer the advantage of converging much more rapidly towards the basis set limit, thanks to their inherently diffuse character. This is further illustrated by the haphazard description of the S1 state in naphthalene, which is easily recovered in plane waves.

5.3

Performance of new customised ‘à-la-carte’ Coulomb-attenuated functionals

Our generalised implementation of the Coulomb-attenuation method renders adjustments to the established xc functionals straightforward, for instance by using different rangeseparation parameters α, β, or by changing µ (commonly referred to as γ-tuning). The flexible form of Kσx also opens the distinct possibility of assembling new long-range corrected functionals based on simple physico-chemical considerations. CAM-B3LYP has become a valuable tool in the calculation of excited states, as illustrated by the substantial improvements of the excitation energies with respect to a simple GGA published in the aformentioned studies. Still, the values in table 3 reveal that not all transitions may be accurately captured quantiatively or qualitatively. The basis-set sensitivity of the excitations in p-nitroaniline are more of a practical issue, since reliable values may still be recovered in an accurately large basis, although at a substantial computational cost. However, the wrong ordering of states in 7-azaindole is one example that can only be overcome by resorting to generally less widely applicable xc functionals such as LC-BLYP, with a similar situation occuring in the β-dipeptide. The orbital localization problems in the AT base pair represent a further challenge. Such problems are often attributed to stem from an imbalance between short- and long-range exchange. The usual remedy to this weakness is the use of LC-BLYP, which fares worse for the systems for which CAM-B3LYP excels, but offers a more reasonable description of compounds where the ordering of CT or Rydberg states proves to be inaccurate. In the following, we have attempted to construct a functional that is sufficiently accurate for both systems where CAM-B3LYP is of satisfactory accuracy, 26


Page 26 of 41

Page 27 of 41 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


as well as for the few charge-transfer systems where the functional has its known weaknesses. A modified version of CAM-B3LYP with 80% exact exchange at long range was reported to yield eigenvalue differences closer to LC-BLYP, but performed worse for other properties where the ‘standard’ CAM-B3LYP yields accurate results. 11 LC-BLYP includes no exact exchange at shortest range, whereas CAM-B3LYP still includes 19%, a value close to the 20% used in standard B3LYP. A suitable compromise between LC-BLYP and CAM-B3LYP could therefore lie in attenuating an existing hybrid functional with more GGA exchange at short range, and less at long range. Following the CAM-B3LYP approach, the short-range contribution of the exact exchange should lie close to the value used in the conventional hybrid to ensure proper balance at short-range. Handy’s OPTX 58 functional in conjecture with Lee-Yang-Parr correlation 59 has on several occasions been shown to be superior to Becke’s 1988 exchange functional, 60–63 and hybrids including OPTX such as O3LYP 58 include a lower percentage of exact exchange than the famous B3LYP while reatining comparable accuracy. We therefore assumed a Coulombattenuated version of O3LYP with 80% exact exchange at long range (as in the CAM-B3LYP declination in reference 11) and only 11% at short range (as in the O3LYP-hybrid) to offer the same benefits as LC-BLYP or CAM-B3LYP with 80% exact exchange at long range, but with an improved description of the short-range region due to the inclusion of the more accurate OPTX. The performance of the CAM-O3LYP functional on the test set used in the previous chapter, including the problematic 7-azaindole, is summarised in table 3. In the dipeptide, much alike the retinal protonated Schiff base, CAM-O3LYP yields virtually indistinguishalbe results from CAM-B3LYP, with a maximum deviation in S4 of 0.06eV. The S1 and S2 states are therefore accurately described, whereas the deviations for S3 and S4 remain too large for practical applications in both of the functionals. We should note that the good agreement between CAM-B3LYP and reference values reported by Peach et al. were based on values obtained in a basis without diffuse functions and without the TDA; the higher energies of the

27



CT states we report for the cc-pVDZ basis are in line with their results. With CAM-O3LYP, the state ordering in 7-azaindole is now correctly reproduced, although the absolute errors are still considerable, with the energetic difference between the 1 La and 1 Lb states being too low. (This is also observed for CAM-B3LYP and may be attributed to the use of the TDA, since in general, the energy difference between the two states becomes larger if the TDA is not employed.) In p-nitroaniline, the spectrum substantially improves with the use of CAM-O3LYP, due to an improved description of the S3 state, which is now predicted even when approaching the basis set limit. For systems in which a classical hybrid functional is preferrable to CAM-B3LYP, such as DMABN, the error due to CAM-O3LYP is comparable to the one of CAM-B3LYP, even slightly smaller in the case of the S2 state. Table 4 shows the excitation energies for a structure of the retinal protonated Schiff base from a molecular dynamics snapshot calculated with both CAM-B3LYP and CAM-O3LYP. Since CAM-B3LYP has been shown to yield very accurate excitation energies in this system, this can serve as an additional benchmark for the accuracy of our new functional. Indeed, CAM-O3LYP yields virtually indistinguishable excitation energies. Table 4: Comparison of excitation energies between CAM-B3LYP and CAM-O3LYP for a structure of the retinal protonated Schiff base obtained from a molecular dynamics snapshot. CAM-B3LYP is known to give good agreement with respect to high-level wavefunction methods for this molecule. 17,38 State CAM-B3LYP CAM-O3LYP retinal protonated Schiff base S1 2.66 2.63 S2 4.13 4.25 S3 4.97 4.95 S4 5.12 5.08 The remaining ‘problematic cases’, the AT base pair and the β-dipeptide are presented in table 5, where excitation energies obtained from CAM-O3LYP and CAM-B3LYP are compared to reference values. For one of the two notorious cases, CAM-O3LYP outperforms CAM-B3LYP on a qualitative level: The spectrum of the AT base pair is qualitatively correctly reproduced, with the 28


Page 28 of 41

Page 29 of 41 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


HOMO lying on adenine rather than thymine, in contrast to CAM-B3LYP results. Overall, both functionals capture the energetics of all but the S1 state accurately, but CAM-O3LYP gives a correct theoretical description of the orbital localization. Only in the β-dipeptide there is no improvement with respect to CAM-B3LYP, but the qualitative and quantitative behaviour is again very similar. Table 5: Comparison of excitation energies between CAM-B3LYP, CAM-O3LYP and reference values for some typical usage cases of range-separated functionals known to be problematic when treated with CAM-B3LYP. The character of the transition is given where it differs with respect to the reference. The reference values for the β-dipeptide correspond to those adopted by Peach et al.; 9 remaining values for molecules included therein and not explicitly mentioned so far are tabluated in the Supporting Information. State CAM-B3LYP AT base pair 64 S1 5.25 5.26 S2 S3 5.36 S4 5.38 26 β-dipeptide S1 5.69 n1 → π1∗ S2 5.76 n2 → π2∗ S3 6.06 π1 → π2∗ 6.14 n1 → π2∗ S4

CAM-O3LYP

Ref.

5.24 5.29 5.40 5.45

4.94 5.21 5.40 5.47

5.66 5.74 6.00 6.01

n1 n2 π1 n1

→ π1∗ → π2∗ → π2∗ → π2∗

5.10 5.40 7.99 9.13

n2 n1 π1 n1

→ π2∗ → π1∗ → π2∗ → π2∗

Overall, CAM-O3LYP shows identical or superior performance to CAM-B3LYP for all of the excitations; it appears to be more versatile in the description of systems that require a larger percentage of exact exchange at long range, where it is able to remedy some of the pitfalls encountered with CAM-B3LYP. These are systems where LC-BLYP has typically been used so far. While the mean absolute error (MAE) over all the molecules considered here is about MAE = 0.65eV for CAM-B3LYP, it improves for CAM-O3LYP, where MAE = 0.55eV. This demonstrates that based on simple considerations, a new Coulomb-attenuated functional can be constructed. CAM-O3LYP predicts qualitatively correct state ordering where CAM-B3LYP fails, and in these systems yields excitation energies similar to the latter. With CAM-O3LYP, it is therefore possible to cover a larger range of systems than with 29



CAM-B3LYP with slightly improved accuracy. This is further reflected by the MAE of both functionals, which is about 0.1eV lower for the CAM-O3LYP presented here. For other delicate systems and specific problems, our implementation offers the possibilty of adapting existing xc functionals at hand, or to assemble entirely new Coulomb-attenuated functionals.

6

Conclusions

We have presented a new, efficient and fully flexible implementation of Coulomb-attenuated functionals in the plane wave/pseudopotential code CPMD which allows for a customised composition of exchange-correlation functionals. On the base of a comprehensive test suite, we could demonstrate that the results obtained within the plane wave/pseudopotential framework do not significantly deviate with respect to results obtained in all-electron calculations with Gaussian bases. The results indicate that the complete basis set limit is more easily reached in plane waves, and that the pseudisation of the nuclear cusp is of no relevant influence on HOMO-LUMO gaps or excitation energies. We have also shown that based on the same considerations that led to the construction of CAM-B3LYP, a new xc functional CAMO3LYP can be constructed, which shows improved performance over CAM-B3LYP in systems where the latter fails, and yields comparable accuracy for systems where CAM-B3LYP typically fares well. This demonstrates that the flexibility of ‘à-la-carte’ combinations of xc functionals can help in obtaining excitations of higher accuracy over a larger range of systems.

Acknowledgement The authors thank J. Magnus Haugaard Olsen for providing a molecular dynamics snapshot of the retinal protonated Schiff base, Valery Weber for comments on the implementation of the new xc driver, and they thank the CADMOS project for providing computing time on the BG/Q. UR gratefully acknowledges financial support from the Swiss National Science 30


Page 30 of 41

Page 31 of 41 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


Foundation Grant No. 200020-146645 and the NCCR MUST.

Supporting Information Available • Convergence of the HOMO-LUMO gaps in water using MT pseudopotentials • Simulation supercell dimensions for all systems • Additional excitation energies for the molecules of ref. 9 that have not been tabulated in the main text

References (1) Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 1964, 136, B864– B871. (2) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133–A1138. (3) Runge, E.; Gross, E. K. U. Density-Functional Theory for Time-Dependent Systems. Phys. Rev. Lett. 1984, 52, 997–1000. (4) Casida, M.; Ipatov, A.; Cordova, F. In Time-Dependent Density Functional Theory; Marques, M. A., Ullrich, C. A., Nogueira, F., Rubio, A., Burke, K., Gross, E. K. U., Eds.; Springer Berlin Heidelberg: Berlin, Heidelberg, 2006; pp 243–257. (5) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. A long-range-corrected time-dependent density functional theory. The Journal of Chemical Physics 2004, 120, 8425–8433. (6) Langreth, D. C.; Perdew, J. P. Theory of nonuniform electronic systems. I. Analysis of

31



the gradient approximation and a generalization that works. Phys. Rev. B 1980, 21, 5469–5493. (7) Tozer, D. J.; Amos, R. D.; Handy, N. C.; Roos, B. O.; Serrano-ANDRES, L. Does density functional theory contribute to the understanding of excited states of unsaturated organic compounds? Molecular Physics 1999, 97, 859–868. (8) Becke, A. D. A new mixing of Hartree-Fock and local density-functional theories. The Journal of Chemical Physics 1993, 98, 1372–1377. (9) Peach, M. J. G.; Benfield, P.; Helgaker, T.; Tozer, D. J. Excitation energies in density functional theory: An evaluation and a diagnostic test. The Journal of Chemical Physics 2008, 128, 044118. (10) Iikura, H.; Tsuneda, T.; Yanai, T.; Hirao, K. A long-range correction scheme for generalized-gradient-approximation exchange functionals. The Journal of Chemical Physics 2001, 115, 3540–3544. (11) Yanai, T.; Tew, D. P.; Handy, N. C. A new hybrid exchange-correlation functional using the Coulomb-attenuating method (CAM-B3LYP). Chemical Physics Letters 2004, 393, 51 – 57. (12) Peverati, R.; Truhlar, D. G. Quest for a universal density functional: the accuracy of density functionals across a broad spectrum of databases in chemistry and physics. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 2014, 372 . (13) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. The Journal of Physical Chemistry 1994, 98, 11623–11627.

32


Page 32 of 41

Page 33 of 41 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


(14) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (15) Jensen, L.; Govind, N. Excited States of DNA Base Pairs Using Long-Range Corrected Time-Dependent Density Functional Theory. The Journal of Physical Chemistry A 2009, 113, 9761–9765. (16) Arulmozhiraja, S.; Coote, M. L. 1 La and 1 Lb States of Indole and Azaindole: Is Density Functional Theory Inadequate? Journal of Chemical Theory and Computation 2012, 8, 575–584. (17) Valsson, O.; Campomanes, P.; Tavernelli, I.; Rothlisberger, U.; Filippi, C. Rhodopsin Absorption from First Principles: Bypassing Common Pitfalls. Journal of Chemical Theory and Computation 2013, 9, 2441–2454. (18) Laktionov, A.; Rothlisberger, U. 2014, unpublished. (19) Perdew, J. P.; Ernzerhof, M.; Burke, K. Rationale for mixing exact exchange with density functional approximations. The Journal of Chemical Physics 1996, 105, 9982– 9985. (20) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. The Journal of Chemical Physics 1999, 110, 6158–6170. (21) Peach, M. J. G.; Cohen, A. J.; Tozer, D. J. Influence of Coulomb-attenuation on exchange-correlation functional quality. Phys. Chem. Chem. Phys. 2006, 8, 4543–4549. (22) Rohrdanz, M. A.; Herbert, J. M. Simultaneous benchmarking of ground- and excitedstate properties with long-range-corrected density functional theory. The Journal of Chemical Physics 2008, 129, 034107. (23) Riffet, V.; Jacquemin, D.; Cauët, E.; Frison, G. Benchmarking DFT and TD-DFT

33



Functionals for the Ground and Excited States of Hydrogen-Rich Peptide Radicals. Journal of Chemical Theory and Computation 2014, 10, 3308–3318. (24) Eriksen, J. J.; Sauer, S. P.; Mikkelsen, K. V.; Christiansen, O.; Jensen, H. J. A.; Kongsted, J. Failures of TDDFT in describing the lowest intramolecular charge-transfer excitation in para-nitroaniline. Molecular Physics 2013, 111, 1235–1248. (25) Sevilla, M. D. Comment on ‘Excited States of DNA Base Pairs Using Long-Range Corrected Time-Dependent Density Functional Theory’. The Journal of Physical Chemistry A 2009, 113, 11093–11094. (26) Serrano-Andrés, L.; Fülscher, M. P. Theoretical Study of the Electronic Spectroscopy of Peptides. III. Charge-Transfer Transitions in Polypeptides. Journal of the American Chemical Society 1998, 120, 10912–10920. (27) Odelius, M.; Laikov, D.; Hutter, J. Excited state geometries within time-dependent and restricted open-shell density functional theories. Journal of Molecular Structure: THEOCHEM 2003, 630, 163 – 175, WATOC ’02 Special Issue. (28) Kronik, L.; Stein, T.; Refaely-Abramson, S.; Baer, R. Excitation Gaps of Finite-Sized Systems from Optimally Tuned Range-Separated Hybrid Functionals. Journal of Chemical Theory and Computation 2012, 8, 1515–1531. (29) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. The Journal of Chemical Physics 2003, 118, 8207–8215. (30) www.cpmd.org, CPMD. Copyright IBM Corp 1990-2008, Copyright MPI für Festkörperforschung Stuttgart 1997-2001. (31) Weber, V.; Bekas, C.; Laino, T.; Curioni, A.; Bertsch, A.; Futral, S. Shedding Light on Lithium/Air Batteries Using Millions of Threads on the BG/Q Supercomputer.

34


Page 34 of 41

Page 35 of 41 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


2014 IEEE 28th International Parallel and Distributed Processing Symposium. 2014; pp 735–744. (32) Gygi, F.; Baldereschi, A. Self-consistent Hartree-Fock and screened-exchange calculations in solids: Application to silicon. Phys. Rev. B 1986, 34, 4405–4408. (33) Broqvist, P.; Alkauskas, A.; Pasquarello, A. Hybrid-functional calculations with planewave basis sets: Effect of singularity correction on total energies, energy eigenvalues, and defect energy levels. Phys. Rev. B 2009, 80, 085114. (34) Xu, X.; Goddard, W. A. The X3LYP extended density functional for accurate descriptions of nonbond interactions, spin states, and thermochemical properties. Proceedings of the National Academy of Sciences of the United States of America 2004, 101, 2673– 2677. (35) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 1988, 38, 3098–3100. (36) Perdew, J. P. Electronic structure of solids; Akademie Verlag, Berlin, 1991; Vol. 11. (37) Rostov, I. V.; Amos, R. D.; Kobayashi, R.; Scalmani, G.; Frisch, M. J. Studies of the Ground and Excited-State Surfaces of the Retinal Chromophore using CAM-B3LYP. The Journal of Physical Chemistry B 2010, 114, 5547–5555. (38) Sneskov, K.; Olsen, J. M. H.; Schwabe, T.; Hattig, C.; Christiansen, O.; Kongsted, J. Computational screening of one- and two-photon spectrally tuned channelrhodopsin mutants. Phys. Chem. Chem. Phys. 2013, 15, 7567–7576. (39) Palczewska, G.; Vinberg, F.; Stremplewski, P.; Bircher, M. P.; Salom, D.; Komar, K.; Zhang, J.; Cascella, M.; Wojtkowski, M.; Kefalov, V. J.; Palczewski, K. Human infrared vision is triggered by two-photon chromophore isomerization. Proceedings of the National Academy of Sciences 2014, 111, E5445–E5454. 35



(40) Aidas, K.; Angeli, C.; Bak, K. L.; Bakken, V.; Bast, R.; Boman, L.; Christiansen, O.; Cimiraglia, R.; Coriani, S.; Dahle, P.; Dalskov, E. K.; Ekström, U.; Enevoldsen, T.; Eriksen, J. J.; Ettenhuber, P.; Fernández, B.; Ferrighi, L.; Fliegl, H.; Frediani, L.; Hald, K.; Halkier, A.; Hättig, C.; Heiberg, H.; Helgaker, T.; Hennum, A. C.; Hettema, H.; Hjertenæs, E.; Høst, S.; Høyvik, I.-M.; Iozzi, M. F.; Jansík, B.; Jensen, H. J. Aa.; Jonsson, D.; Jørgensen, P.; Kauczor, J.; Kirpekar, S.; Kjærgaard, T.; Klopper, W.; Knecht, S.; Kobayashi, R.; Koch, H.; Kongsted, J.; Krapp, A.; Kristensen, K.; Ligabue, A.; Lutnæs, O. B.; Melo, J. I.; Mikkelsen, K. V.; Myhre, R. H.; Neiss, C.; Nielsen, C. B.; Norman, P.; Olsen, J.; Olsen, J. M. H.; Osted, A.; Packer, M. J.; Pawlowski, F.; Pedersen, T. B.; Provasi, P. F.; Reine, S.; Rinkevicius, Z.; Ruden, T. A.; Ruud, K.; Rybkin, V. V.; Sałek, P.; Samson, C. C. M.; de Merás, A. S.; Saue, T.; Sauer, S. P. A.; Schimmelpfennig, B.; Sneskov, K.; Steindal, A. H.; SylvesterHvid, K. O.; Taylor, P. R.; Teale, A. M.; Tellgren, E. I.; Tew, D. P.; Thorvaldsen, A. J.; Thøgersen, L.; Vahtras, O.; Watson, M. A.; Wilson, D. J. D.; Ziolkowski, M.; Ågren, H. The Dalton quantum chemistry program system. WIREs Comput. Mol. Sci. 2014, 4, 269–284. (41) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannen36


Page 36 of 41

Page 37 of 41 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


berg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ã.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian09 Revision D.01. Gaussian Inc. Wallingford CT 2009. (42) Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. The Journal of Chemical Physics 1989, 90, 1007–1023. (43) Van Lenthe, E.; Baerends, E. J. Optimized Slater-type basis sets for the elements 1-118. Journal of Computational Chemistry 2003, 24, 1142–1156. (44) te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca Guerra, C.; van Gisbergen, S. J. A.; Snijders, J. G.; Ziegler, T. Chemistry with ADF. J. Comput. Chem. 2001, 22, 931–967. (45) Baerends, E. J.; Ziegler, T.; Atkins, A. J.; Autschbach, J.; Bashford, D.; Baseggio, O.; Bérces, A.; Bickelhaupt, F. M.; Bo, C.; Boerritger, P. M.; Cavallo, L.; Daul, C.; Chong, D. P.; Chulhai, D. V.; Deng, L.; Dickson, R. M.; Dieterich, J. M.; Ellis, D. E.; van Faassen, M.; Ghysels, A.; Giammona, A.; van Gisbergen, S. J. A.; Goez, A.; Götz, A. W.; Gusarov, S.; Harris, F. E.; van den Hoek, P.; Hu, Z.; Jacob, C. R.; Jacobsen, H.; Jensen, L.; Joubert, L.; Kaminski, J. W.; van Kessel, G.; König, C.; Kootstra, F.; Kovalenko, A.; Krykunov, M.; van Lenthe, E.; McCormack, D. A.; Michalak, A.; Mitoraj, M.; Morton, S. M.; Neugebauer, J.; Nicu, V. P.; Noodleman, L.; Osinga, V. P.; Patchkovskii, S.; Pavanello, M.; Peeples, C. A.; Philipsen, P. H. T.; Post, D.; Pye, C. C.; Ramanantoanina, H.; Ramos, P.; Ravenek, W.; Rodríguez, J. I.; Ros, P.; Rüger, R.; Schipper, P. R. T.; Schlüns, D.; van Schoot, H.; Schreckenbach, G.; Seldenthuis, J. S.; Seth, M.; Snijders, J. G.; Solà, M.; M., S.; Swart, M.; Swerhone, D.; te Velde, G.; Tognetti, V.; Vernooijs, P.; Versluis, L.; Visscher, L.; Visser, O.; Wang, F.; Wesolowski, T. A.; van Wezenbeek, E. M.; Wiesenekker, G.; Wolff, S. K.; Woo, T. K.;

37



Yakovlev, A. L. ADF2017, SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam, The Netherlands, https://www.scm.com. (46) Fonseca Guerra, C.; Snijders, J. G.; te Velde, G.; Baerends, E. J. Towards an order-N DFT method. Theoretical Chemistry Accounts 1998, 99, 391–403. (47) Hirata, S.; Head-Gordon, M. Time-dependent density functional theory within the Tamm-Dancoff approximation. Chemical Physics Letters 1999, 314, 291 – 299. (48) Troullier, N.; Martins, J. L. Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B 1991, 43, 1993–2006. (49) Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials. Phys. Rev. B 1996, 54, 1703–1710. (50) Martyna, G. J.; Tuckerman, M. E. A reciprocal space based method for treating long range interactions in ab initio and force-field-based calculations in clusters. The Journal of Chemical Physics 1999, 110, 2810–2821. (51) Bekas, C.; Curioni, A.; Andreoni, W. Atomic wavefunction initialization in ab initio molecular dynamics using distributed Lanczos. Parallel Computing 2008, 34, 441 – 450, Parallel Matrix Algorithms and Applications. (52) Bekas, C.; Curioni, A. Very large scale wavefunction orthogonalization in Density Functional Theory electronic structure calculations. Computer Physics Communications 2010, 181, 1057 – 1068. (53) Peach, M. J. G.; Tozer, D. J. Overcoming Low Orbital Overlap and Triplet Instability Problems in TDDFT. The Journal of Physical Chemistry A 2012, 116, 9783–9789. (54) Grimme, S.; Parac, M. Substantial Errors from Time-Dependent Density Functional Theory for the Calculation of Excited States of Large π Systems. ChemPhysChem 2003, 4, 292–295. 38


Page 38 of 41

Page 39 of 41 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


(55) Bulliard, C.; Allan, M.; Wirtz, G.; Haselbach, E.; Zachariasse, K. A.; Detzer, N.; Grimme, S. Electron Energy Loss and DFT/SCI Study of the Singlet and Triplet Excited States of Aminobenzonitriles and Benzoquinuclidines:âĂĽ Role of the Amino Group Twist Angle. The Journal of Physical Chemistry A 1999, 103, 7766–7772. (56) Svartsov, Y. N.; Schmitt, M. Electronically excited states of water clusters of 7azaindole: Structures, relative energies, and electronic nature of the excited states. The Journal of Chemical Physics 2008, 128, 214310. (57) Kosenkov, D.; Slipchenko, L. V. Solvent Effects on the Electronic Transitions of pNitroaniline: A QM/EFP Study. The Journal of Physical Chemistry A 2011, 115, 392–401. (58) Handy, N. C.; Cohen, A. J. Left-right correlation energy. Molecular Physics 2001, 99, 403–412. (59) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785–789. (60) Hoe, W.-M.; Cohen, A. J.; Handy, N. C. Assessment of a new local exchange functional OPTX. Chemical Physics Letters 2001, 341, 319 – 328. (61) Baker, J.; Pulay, P. Assessment of the Handy-Cohen optimized exchange density functional for organic reactions. The Journal of Chemical Physics 2002, 117, 1441–1449. (62) Baker, J.; Pulay, P. Assessment of the OLYP and O3LYP density functionals for firstrow transition metals. Journal of Computational Chemistry 2003, 24, 1184–1191. (63) Xu, X.; Goddard, W. A. Assessment of Handy-Cohen Optimized Exchange Density Functional (OPTX). The Journal of Physical Chemistry A 2004, 108, 8495–8504. (64) Lange, A. W.; Herbert, J. M. Both Intra- and Interstrand Charge-Transfer Excited States in Aqueous B-DNA Are Present at Energies Comparable To, or Just Above, the 39



1ππ ∗ Excitonic Bright States. Journal of the American Chemical Society 2009, 131, 3913–3922.

40


Page 40 of 41

Page 41 of 41 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


Figure 2: For Table of Contents Only

41


Plane-Wave Implementation and Performance of Ã -la-Carte Coulomb

Recommend Documents