Formulation and Implementation of the Spin-Restricted Ensemble

Subscriber access provided by Lancaster University Library

Quantum Electronic Structure

Formulation and Implementation of the Spin-Restricted Ensemble-Referenced Kohn-Sham Method in the Context of Density Functional Tight Binding Approach In Seong Lee, Michael Filatov, and Seung Kyu Min J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.9b00132 • Publication Date (Web): 10 Apr 2019 Downloaded from http://pubs.acs.org on April 12, 2019

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

Formulation and Implementation of the Spin-Restricted Ensemble-Referenced Kohn-Sham Method in the Context of Density Functional Tight Binding Approach In Seong Lee,† Michael Filatov,‡ and Seung Kyu Min∗,† †Department of Chemistry, School of Natural Science, Ulsan National Institute of Science and Technology (UNIST), 50 UNIST-gil, Ulsan 44919, Republic of Korea ‡Department of Chemistry, Kyungpook National University, Daegu 702-701, Republic of Korea E-mail: [email protected] Abstract The spin-restricted ensemble-referenced Kohn-Sham (REKS) method and its stateinteraction state-averaged variant (SI-SA-REKS, or SSR) provide computational platform for seamless inclusion of multi-reference effects into the density functional calculations. The SSR method enables an accurate calculation of the vertical excitation energies for the molecules with multi-reference ground states and describes conical intersections between the ground and excited states with the accuracy matching the most sophisticated ab initio multi-reference wavefunction methods. In this work, the SSR method is formulated and implemented in the context of the long-range corrected density functional tight binding (LC-DFTB) approach. The new LC-DFTB/SSR method enables calculation of the excited electronic states and the S1 /S0 conical intersections

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

of very large molecules. The LC-DFTB/SSR method is benchmarked against vertical excitation energies and conical intersection energies and geometries of several organic molecules with π/π ∗ and n/π ∗ transitions. It is demonstrated that the LC-DFTB/SSR method describes these molecules with reasonable accuracy, which can be considerably improved by a slight modification of the LC-DFTB spin polarization parameters.

Graphical TOC Entry geometry

2


Page 2 of 37

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


Introduction Phenomena occurring in the electronically excited states of molecules are central for design and application of molecular devices such as, e.g., light-driven molecular rotary motors 1–3 and photosynthetic/photovoltaic materials. 4–7 Accurate theoretical description of the shape of the excited state potential energy surfaces (PESs) and real crossings between them (conical intersections, CIs) is at the heart of understanding of these phenomena. 8–11 In spite of recent fast progress in the computer technology, computation of the excited states of large and very large molecular systems still remains largely beyond the reach of most computational methods, be it the multi-reference methods of ab initio wavefunction theory, such as the multi-reference configuration interaction 12 or the active space methods, 13–19 or the linear response time dependent density functional theory 20 (TD-DFT) methods. The DFTB methodology furnishes a semi-empirical approximate yet accurate framework for obtaining energies and geometries of very large molecular systems. 21–23 The implementation of the time-dependent DFTB (TD-DFTB) method opened up possibility to access excited electronic states of these systems. 24,25 A considerable improvement of the description of the excited states, especially the states with charge transfer characteristic, was achieved with the long-range corrected DFTB approach. 26–28 However, even in its most accurate variant, the TD-DFTB methodology is not able to describe CI between the ground and excited electronic states, as it is based on the spin-conserving linear response formalism. 29–31 CIs are of paramount importance for the mechanism of the excited state reactions, 32–34 as they provide funnels for the ultrafast population transfer between the excited and the ground electronic states. The accurate description of the double cone topology of CIs requires the use of multi-reference methods and the topography of the PESs in the CIs vicinity is sensitive to the inclusion of the dynamic electron correlation. 35,36 The inclusion of the dynamic electron correlation in the multi-reference ab initio wavefunction methods imposes enormous computational burden and restricts application of these methods to relatively small molecules. Within the DFT framework, the (spin-conserving) linear response TD-DFT fails 3



to produce the correct double cone topology of the PESs at the intersection point. 29–31 The spin-flip TD-DFT, 37–40 which is based on the use of open-shell triplet reference state, can recover the correct topology; 31,50 however, at the price of substantial spin-contamination of the response states, which renders automatic identification of proper singlet states impossible. Another alternative to the traditional wavefunction multi-reference methods and density functional methods is offered by methods based on ensemble density functional theory (eDFT), 41 e.g., the REKS and the SSR methods. 42–47 At a modest computational cost, the latter method enables very accurate description of the vertical excitation energies and the S1 /S0 CIs matching the results of the most sophisticated multi-reference wavefunction calculations, e.g., MRSDCI. 30,48 In this paper, the SSR method is formulated and implemented in the context of DFTB methodology, in particular the LC-DFTB method. The equations for the energy and the analytic energy gradient of the ground and excited electronic states are derived for the LC-DFTB/SSR method. The new method is implemented in the development version of the DFTB+ code 49 and benchmarked in the calculation of the geometries, branching plane vectors, and relative energies of the S1 /S0 CIs of several organic molecules, including the popular molecular model of the retinal chromophore. 50,51

Computational methods LC-DFTB/SSR methodology The SSR method has been previously described in the literature, 42–47 and only a brief account of its most important features will be given here. The SSR methodology employs ground state eDFT to describe multi-reference ground states of molecules and eDFT for ensembles of ground and excited state to obtain excitation energies from a variational time-independent formalism. 41 The use of an ensemble representation for the density and the energy of a strongly correlated molecular system leads to the occurrence of the fractional occupation 4


Page 4 of 37

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


numbers (FONs) of several frontier Kohn-Sham (KS) orbitals. Currently, the SSR method is formulated for active spaces including two electrons in two fractionally occupied orbitals (i.e., SSR(2,2)) and for four electrons in four fractionally occupied orbitals (SSR(4,4)). 52 However, in the following, the derivations will be restricted to the SSR(2,2) method. In the SSR(2,2) method, the energy of multi-reference state is expanded in terms of six microstates with the fixed and integer occupations of the active orbitals φa and φb ; (nLaα nLaβ nLbα nLbβ ) for the Lth microstate are (1100), (0011), (1001), (0110), (1010), and (0101) for L = 1, 2, · · · , 6, respectively. The energies of the ground and excited states are obtained from variational minimization (with respect to the KS orbitals and their FONs) of an ensemble of the REKS(2,2) perfectly spin-paired singlet (PPS) electronic configuration and an open-shell singlet (OSS) configuration obtained within the same (2,2) active space 53 followed by solving 2×2 secular equation PPS

 E  ∆SA

SA

∆

E OSS











  a00 a01   =  a10 a11

E0SSR

0

0

E1SSR



  a00 a01    a10 a11

(1)

to include possible coupling between the PPS and the OSS electronic configurations. The energies of PPS and OSS states, E P P S and E OSS , respectively, are expressed in the form of an ensemble average over the six microstates as

E

PPS

=

E OSS =

6 X L=1 6 X

CLP P S EL

(2)

CLOSS EL

(3)

L=3

where C1P P S = na /2, C2P P S = nb /2, C3P P S = C4P P S = −C5P P S = −C6P P S = −f (na /2)/2, and C3OSS = C4OSS = −2C5OSS = −2C6OSS = 1 with na and nb as FONs of the active KS orbitals 1− 21

φa and φb and f (na /2) = (na nb )

na nb +δ 1+δ

with δ = 0.4 subject to na + nb = 2. Within

SA-REKS(2,2) orbital optimization, the SA-REKS(2,2) energy, E SA = w0 E P P S + w1 E OSS

5



Page 6 of 37

with w0 + w1 = 1, is minimized with respect to KS orbitals φp as well as FONs na and nb under the constraint of orbital orthonormality. Typically, the weights w0 = w1 = 1/2 are chosen. As a result of orbital optimization, a single particle equation is obtained,

np Fˆp φp =

X

SA pq φq ,

(4)

q

where the Fock operator Fˆp of the pth orbital φp is given by,

Fˆp =

X L

np =

P

L

CLSA

nLpα FˆαL + nLpβ FˆβL ; p ∈ core, a, b, np

(5)

CLSA (nLpα + nLpβ ) with CLSA = w0 CLP P S + w1 CLOSS , and SA pq is the Lagrangian matrix

SA SA∗ is not invariant under a element with SA pq = qp . Since the total SA-REKS(2,2) energy E

unitary transformation of the KS orbitals, Eq. (4) cannot be solved by simple diagonalization. Instead, we solve Eq. (4) with the use of the coupling operator technique. 54 The interstate coupling parameter ∆SA is calculated using the SA-REKS(2,2) Lagrangian matrix element SA ab between the active orbitals φa and φb as √ √ ∆SA = ( na − nb )SA ab .

(6)

Implementation of DFTB formalism in SSR(2,2) approach is straightforward. The stanP dard DFTB formalism is based on the linear combination of atomic orbitals as φi = µ Cµi χµ with atomic orbitals χµ where µ and other Greek letters except α and β represent atomic orbital (AO) indices. The total energy of the Lth microstate, EL , with electronic configurations {nLiσ } is obtained by long-range corrected spin-polarized DFTB framework 21–23,26,27,55,56

6


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

EL =

X µν

−

L 0 Pµν Hµν +

1 XX 1 X X 0L 0L L ∆Pµν ∆PτLγ Sµν Sτ γ Γfrµν,τ γ + P P Sµν Sτ γ Wµν,τ γ 8 µν τ γ 8 µν τ γ µν τ γ

1 XX L L ∆PτLγβ )Sµτ Sνγ Γlrµτ,γν + Erep . ∆PτLγα + ∆Pµνβ (∆Pµνα 8 µν τ γ

(7)

The first term in Eq. (7) represents the zeroth order energy contribution from reference density matrix P0 which is usually chosen as the superposition of neutral atomic density 0 is evaluated at P0 and the total matrices. The zeroth order Hamiltonian matrix element Hµν P P L L L L with Pµνσ = i nLiσ Cµi Cνi . 21,22 Here is defined as Pµν = σ Pµνσ density matrix element Pµν

we note that we consider real-valued molecular orbital (MO) coefficients Cµi ’s throughout this paper. The second term represents a sum of Coulomb interaction energies among Mulliken 0 L L and Sµν is an overlap matrix element. We − Pµν = Pµν atomic charge pairs where ∆Pµν fr fr fr fr fr define Γfrµν,τ γ as Γfrµν,τ γ = γµτ + γµγ + γντ + γνγ where γµν is a full range two-electron integral

which is parametrized as a function of interatomic distances. 23 The third term arises when we deal with spin-polarized electronic configurations such as microstates with L = 3, 4, 5 and L L 0L 0L while Wµν,τ γ is − Pµνβ = Pµνα are defined as Pµν 6. The spin density matrix elements Pµν

Wµν,τ γ = Wµτ +Wµγ +Wντ +Wνγ where Wµν is an atomic spin constant which is nonzero only if µ and ν are on-site elements. 55,56 The fourth term is long-range corrected energy contribution L L 0 where ∆Pµνσ = Pµνσ − Pµνσ with a reference density matrix P0σ for a given spin σ and Γlrµτ,γν lr lr lr is a sum of long-range two-electron integrals γµν as Γlrµτ,γν = γµγ + γµν + γτlrγ + γτlrν . 26,27 The fr lr explicit expressions of γµν and γµν can be found in the literature. 23,26,27 The last term of

Eq. (7) represents a repulsive energy of atomic sites which depends on the reference density P0 . L The corresponding Fock matrix element Fµνσ in terms of atomic orbitals to construct the

7



Page 8 of 37

Fock operator in Eq. (5) is given by 1X 1 X 0L ∆PτLγ Sµν Sτ γ Γfrµν,τ γ + δσ P Sµν Sτ γ Wµν,τ γ 4 τγ 4 τγ τγ 1X − ∆PτLγσ Sµτ Sνγ Γlrµτ,γν 4 τγ

L 0 Fµνσ =Hµν +

(8)

where δσ = ±1 for σ = α and σ = β, respectively. Here we note that the spin-restricted scheme yields identities; Fˆα1 = Fˆβ1 , Fˆα2 = Fˆβ2 , Fˆα3 = Fˆβ4 , Fˆα4 = Fˆβ3 , Fˆα5 = Fˆβ6 , and Fˆα6 = Fˆβ5 . The overall self-consistent procedure of LC-DFTB/SSR approach is summarized as follows: (a) From an initial guess of MO coefficients Cµi ’s, which are usually obtained from the diagonalization of the zeroth-order Hamiltonian matrix H0 , the density matrix elements L L for each microstate L are constructed according to and the Fock matrix elements Fµνσ Pµνσ

Eq. (8). (b) The microstate energy EL with given electronic configuration {nLiσ } is calculated by Eq. (7). (c) The FONs (na and nb ) are optimized by the Newton-Raphson method subject to na + nb = 2 to minimize the SA-REKS(2,2) energy, E SA = w0 E P P S + w1 E OSS . (d) From the optimized FONs or CLSA ’s, the single particle Fock operator are calculated according to Eq. (5), and the single particle equation (4) is solved by the coupling operator technique to obtain new KS orbitals. (e) The steps (a-d) are repeated until the density matrix P is converged. Then, E P P S , E OSS , and ∆SA are calculated by Eq. (2), (3), and (6), respectively. After solving the secular equation (1) with given optimized KS orbitals and FONs, the ground (E0SSR ) and the lowest excited (E1SSR ) state energies are obtained as EISSR = a20I E P P S + a21I E OSS + 2a0I a1I ∆SA ; I = 0 and 1.

(9)

Analytical gradients The analytical gradients of the individual electronic states of SA-REKS and SSR with respect to nuclear coordinates and the nonadiabatic coupling vectors between the SSR states are obtained as follows. According to the Hellmann-Feynman theorem, 57 differentiating the 8


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


energies in Eq. (9) with respect to an external parameter λ, such as an external field or an atomic coordinate, leads to ∂E P P S ∂E OSS ∂∆SA ∂EISSR = a20I + a21I + 2a0I a1I . ∂λ ∂λ ∂λ ∂λ

(10)

The derivatives of the individual energies E P P S and E OSS with respect to λ are given by 47 X ∂E X ∂ 0 EL 1 X i X j X ∂ 0 Sji X i X j X λ = CLX − + ij + ij ij − ij Uji ∂λ ∂λ 2 ∂λ ij ij L p X where p X ij for X=PPS, OSS, or SA is ij =

P

L

CLX

P

σ

nLpσ Lijσ with Lijσ =

P

µν

(11)

L Cµi Fµνσ Cνi

calculated with the real-valued MO coefficients Cµi and the overlap matrix Sij between the MOs. The notation

∂0 ∂λ

means that only the integrals over the basis functions are to be

differentiated, not the MO coefficients. With the DFTB expression Eq. (7), the derivative ∂ 0 EL ∂λ

in Eq. (11) becomes 0 ∂ 0 EL X L ∂Hµν ∂Sτ γ 1 XX ∂Sµν L L Pµν = + ∆Pµν ∆Pτ γ Sτ γ + Sµν Γfrµν,τ γ ∂λ ∂λ 8 ∂λ ∂λ µν µν τ γ X X 1 ∂Sτ γ 0L 0L ∂Sµν + Pµν Pτ γ Sτ γ + Sµν Wµν,τ γ 8 µν τ γ ∂λ ∂λ ∂S 1 XX ∂Sτ γ lr µν L L − (∆Pµνα Sτ γ + Sµν Γµτ,γν ∆PτLγα + ∆Pµνβ ∆PτLγβ ) 8 µν τ γ ∂λ ∂λ +

∂Γfrµν,τ γ 1 XX L ∆Pµν ∆PτLγ Sµν Sτ γ 8 µν τ γ ∂λ

−

∂Γlrµτ,γν ∂Erep 1 XX L L (∆Pµνα ∆PτLγα + ∆Pµνβ ∆PτLγβ )Sµν Sτ γ + . 8 µν τ γ ∂λ ∂λ

The orbital response matrix Uji0λ is defined as

∂ C ∂λ µi

=

P

j

(12)

Cµj Uji0λ , where Uji0λ can be separated

into the symmetric and anti-symmetric parts as Uji0λ = − 12

∂ 0 Sji ∂λ

+ Ujiλ ; with anti-symmetric

Ujiλ . As the SA-REKS orbitals φp are not eigenfunctions of the Fock operator of an individual state in SA-REKS(2,2) formalism, the last term in Eq. (11) does not vanish and the orbital

9



Page 10 of 37

response matrix Uλ should be calculated. The orbital response matrix is obtained from solving the coupled-perturbed (CP) equaSA∗ tions derived from differentiation of the identity SA pq = qp , which represents the varia-

tional condition for the SA-REKS orbitals. 47 Differentiating the condition hφp |np Fˆp |φq i = hφp |nq Fˆq |φq i, with respect to λ yields 47 X

δip

p SA qj

λ X Uji + δiq − q SA qj

X ∂C SA X L

L

+

λ X (p−q)(k−l) − q SA Uji + (pq|f˜Hxc |kl)Ulkλ pj

ij

ij

+

p SA pj

1X 2 t

∂λ

kl

nLpσ − nLqσ Lpqσ =

X

σ

kl

∂ 0S lk pk qk (pq|f˜Hxc |kl) − (pq|f˜Hxc |kl) ∂λ

0

p SA tq

− q SA tq

∂ Stp 1 X + ∂λ 2 t

p SA pt

− q SA pt

∂ 0 Lpqσ ∂ 0 Stq X SA X L CL npσ − nLqσ − ∂λ ∂λ σ L (13)

P SA P (i−j)(k−l) L L L L ˜σσ0 ˜pq where (ij|f˜Hxc |kl) = σσ 0 (niσ − njσ )(nkσ 0 − nlσ 0 )(ij|fHxc |kl), (ij|fHxc |kl) = L CL P SA P L L ˜σσ0 |kl) and ∂ 0 L /∂λ = P Cµp ∂ 0 F L /∂λCνq . Here (pq|f˜σσ0 |kl) = n n (ij| f C 0 0 pσ pqσ µνσ Hxc Hxc L qσ σσ µν L P P 1 σσ 0 σσ 0 µν τ γ Cµp Cνq (µν|fHxc |τ γ) + (µν|fHxc |γτ ) Cτ k Cγl is a symmetrized Hartree exchange2 0

σσ correlation kernel in MO basis with (µν|fHxc |τ γ), the Hartree exchange-correlation kernel in

AO basis, as 1 1 σσ 0 (µν|fHxc |τ γ) = Sµν Sτ γ Γfrµν,τ γ + Kσσ0 Wµν,τ γ − δσσ0 Sµτ Sνγ Γlrµτ,γν 4 4

(14)

in DFTB formalism from Eq. (8). We use identities as Kσσ0 = 1 and δσσ0 = 1 if σ = σ 0 while Kσσ0 = −1 and δσσ0 = 0 if σ 6= σ 0 . Eq. (13) can be expressed in the form of a matrix equation as 47

AUλ = Bλ or

X

Aij,kl Uklλ = Bijλ

(15)

k

Formulation and Implementation of the Spin-Restricted Ensemble

Recommend Documents