Implementation of Noncollinear Spin-Flip Time-Dependent Density-Functional Theory in deMon2k Miquel Huix-Rotllant and Mark E. Casida March 27, 2009

Contents 1 Overview of the Method

1

2 Keyword NONCOL

5

3 Input and the Output examples 3.1 Input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5 5 5

4 Subroutines used in the present implemenation of SF TDDFT

8

5 Known Problems, Suggested Improvements

8

Abstract This document, written by Miquel Huix-Rotllant and Mark E. Casida, describes Andrei Ipatov’s implementation of Wang and Ziegler’s noncollinear spin-flip linear response timedependent density-functional theory (spin-flip TDDFT) in deMon2k. Andrei Ipatov’s routines have been β tested by Miquel Huix-Rotllant, Bhaarathi Natarajan, and C. Muhavini Wawire both on H2 dissociation and for the ring-opening of oxirane (CH2 OCH2 ) along various pathways.

1

Overview of the Method

Careful consideration of Casida’s original formulation [1] of linear-response time-dependent densityfunctional theory (TDDFT) shows that it includes spin-α excitations and spin-β excitations (RHS of Fig. 1) but not spin-flip (SF) excitations turning spin-α electrons into spin-β electrons and vice versa (LHS of Fig. 1). This arises because only spin-independent perturbations were considered. Thus spins remain “collinear,” meaning always oriented along the z-direction. Up until Wang and Ziegler’s remarkable work [2, 3], including SF perturbations which would lead to noncollinear spins seemed to require some contribution of Hartree-Fock exchange [4]. This was not a particularly satisfying situation because the fraction of Hartree-Fock exchange needed for SF TDDFT differed considerably from what is optimal for hybrid functionals used for ground state calculations. Wang and Ziegler overcame this problem by using a “trick” from relativistic density-functional theory. In traditional density-functional theory (DFT) the exchange-correlation (xc) functional is collinear because it assumes that the spins are oriented α and β with respect to the z-axis, Exc ≡ Exc [ρα , ρβ ] . 1

(1)

Figure 1: Two-orbital model of TDDFT exctations.

Wang and Ziegler’s idea was to make this collinear functional into a noncollinear functional, Exc ≡ Exc [ρ+ , ρ− ] ,

(2)

by replacing ρα and ρβ with,

1 (ρ ± s) . (3) 2 The precise definition of the density, ρ, and the spin-polarization, s, involves recognizing that the spin-densities become a spin-density matrix in a proper noncollinear formulation involving 2component spinors, " # ρα,α ρα,β ρ= . (4) ρβ,α ρβ,β ρ± =

Then ρ = trρ = ρα,α + ρβ,β 



s2 = (ρα,α + ρβ,β )2 + 2 ρ2α,β + ρ2β,α ,

(5)

are properly rotationally invariant. We can now rederive linear-response TDDFT for a spindependent perturbation, and so arrive at a generalization of the traditional approach. In the end, however, we take the collinear limit, ρ = trρ = ρα + ρβ s = ρα − ρβ .

2

(6)

Figure 2: Two-orbital model of SF TDDFT excitations.

(Since there are no off-diagonal matrix elements in the collinear limit, the usual simpler notation is used, namely ρα = ρα,α and ρβ = ρβ,β .) The xc kernel then becomes,     

αα,αα fxc ββ,αα fxc αβ,αα fxc βα,αα fxc

αα,ββ fxc ββ,ββ fxc ββ,ββ fxc βα,ββ fxc

αα,αβ fxc ββ,αβ fxc αβ,αβ fxc βα,αβ fxc

αα,βα fxc ββ,βα fxc αβ,βα fxc βα,βα fxc

    

=

α,α α,β fxc fxc  f β,α f β,β  xc xc

    

0 0

0 0 0

0

0

α −v β vxc xc ρα −ρβ

0

0

0

α −v β vxc xc ρα −ρβ

      

(7)

We now see that the conventional (∆Ms = 0) and SF excitations (∆Ms = ±1) are completely decoupled from each other, as are the two SF excitations (∆Ms = +1 and ∆Ms = −1). Application to the two-orbital model of Fig. 1 indicates that the same results are found as with conventional TDDFT provided the local density approximation (LDA) is used. For generalized gradient approximations (GGAs) differences are in principle possible. Usually in SF-TDDFT only the SF part of the coupling matrix is taken. There are obviously two flavors (∆Ms = ±1), but typically only ∆Ms = −1 is used in practice. A typical example is given in Fig. 2 where ∆Ms = −1. SF excitations lead to the possibility of taking linear combinations of ground and doubly excited configurations. This leads to the correct dissociation of the ground-state of H2 [2] since, 1 √ (|σg , σ ¯g | − |σu , σ ¯u |) = 2

1 √ (|sA , s¯A | + |sA , s¯B | + |sB , s¯A | + |sB , s¯B |) 2 2 1 √ (|sA , s¯A | − |sA , s¯B | − |sB , s¯A | + |sB , s¯B |) − 2 2 1 = √ (|sA , s¯B | + |sB , s¯A |) 2 ↔ [H↑ + H↓ ↔ H↓ + H↑] . 3

(8)

Figure 3: Application of SF TDDFT to the ground and excited states of dissociating H2 .

It is also clear from Fig. 2 that SF-TDDFT leads to excitation to what should be an energetically degenerate triplet state. The deviation of this excitation energy from zero may be taken as one test of the method. Typically it is negligeable on the scale of the other excitation energies in the system. Figure 3 shows the application of SF TDDFT to the ground and excited states of dissociating H2 . The ground 11 Σg and lowest triplet 13 Σu states should dissociate to two ground state H atoms in the ground state (2E(H) ≈ −27eV). It is well-known that this happens neither in spin-restricted (i.e., same-orbitals-for-different-spin) Hartree-Fock calculations nor in DFT with approximate functionals for the ground 11 Σg , but does happen for the SF-TDDFT state (red curve.) Though not shown, the SF-TDDFT 13 Σu curve (black) is nearly degenerate with the 13 Σu SCF curve with two spin α electrons. The excited singlets 11 Σu (green curve) and 21 Σg (blue curve) should (with the basis used in this calculation) also dissociate to the same [H+ + H− ↔ H− + H+ ] ionic limit and this is indeed the case with the SF-TDDFT limit. See Ref. [5] for an application to oxirane ring opening. The main drawbacks of SF-TDDFT are the following: 1. DFT is formally only rigorous when applied to the ground state of the system of interest. The lowest energy triplet state is usually not the ground state. 2. Even if we neglect problem 1, problems may arise from inaccuracies in the xc functional, making conventional and SF TDDFT give different results. 3. Only some configurations are available from SF-TDDFT. They may not alone be enough to treat particular multiplet problems. [6] Also (vide infra) impure spin-states can arise in SF-TDDFT. 4. Detailed equations have yet to be worked out to determine SF-TDDFT excited-state spincontamination in the different-orbitals-for-different-spin (DODS) formalism.

4

On the plus side, often what makes the SF-TDDFT method attractive is that it solves problems arising when the HOMO and LUMO become nearly degenerate. In such situations the SDF for the lowest singlet state is difficult to converge while that for the lowest triplet state is much easier.

Keyword NONCOL

2

This keyword activates the Spin-Flip calculation. No options are present. Note that this keyword can only be used together with the keyword EXCIT. For EXCIT the TDA option is automatically selected in the case of a Spin-Flip calculation. The options SING, TRIP and BOTH are not active.

3

Input and the Output examples

We calculate the SF-TDDFT spectra of H2 with a triplet as a reference state and different orbitals for different spins. In this calculation a developement version of deMon2k has been used, based on the master version 2.2.[7]

3.1

Input

We have the following input:

TITLE SF-TDDFT of H2 SCFTYPE UKS MULTI 3 EXCITATION TDA NONCOL BASIS(STO-3G) GEOMETRY ZMATRIX ANGSTROM H H 1 0.7

3.2

Output

The program performs two independent SF-TDDFT calculations, one after the other. The first describes α → β transitions (so ∆Ms = −1). The secondd describes β → α transitions (so ∆Ms = +1). This is the order of results in the output. NOTES: • The number NCISAB is the number of transitions that are going to appear if a full calculation is requested. For brevety the transitions have been cut to 5. • hSˆ2 i is not calculated for SF TDDFT excited states. • The “CI coefficients” (really TDA coefficients) for the different transitions are listed as “i–>j ( xx.yyyy eV) Coeff = .zzzz” The number xx.yyyy eV is the difference of the orbital energies for the jth spin β orbital and the ith spin α orbital. NCISAB =

20

HOMO-ALPHA =

-.6970 eV 5

LUMO-BETA

=

-8.5296 eV

TAMM-DANCOFF APPROXIMATION : REFERENCE VALUE OF S(S+1):

2.0000

CALCULATED VALUE OF S**2 :

2.0000

RESULTING MULTIPLICITY:

3.0000

+-------------------------------------------------+ ENERGY (eV) d MULT. +-------------------------------------------------+ 1 -10.1114 2 -.0002 3 3.4624 4 13.2538 5 17.9097 +-----------------------------------------------+ 1 Transition energy -10.11143 eV : +-----------------------------------------------+ 1--> 2 ( 19.47640 eV) Coeff = .65132E-01 1--> 4 ( 41.85260 eV) Coeff = -.52105E-01 2--> 1 ( -7.83255 eV) Coeff = -.99594E+00 +-----------------------------------------------+ 2 Transition energy -.00023 eV : +-----------------------------------------------+ 1--> 1 ( 7.72967 eV) Coeff = .70040E+00 1--> 3 ( 30.78423 eV) Coeff = .97438E-01 2--> 2 ( 3.91418 eV) Coeff = -.70185E+00 2--> 4 ( 26.29039 eV) Coeff = .85202E-01 +-----------------------------------------------+ 3 Transition energy 3.46244 eV : +-----------------------------------------------+ 1--> 1 ( 7.72967 eV) Coeff = .68681E+00 1--> 3 ( 30.78423 eV) Coeff = .11059E+00 2--> 2 ( 3.91418 eV) Coeff = .71220E+00 2--> 4 ( 26.29039 eV) Coeff = .93150E-01 +-----------------------------------------------+ 4 Transition energy 13.25380 eV : +-----------------------------------------------+ 1--> 4 ( 41.85260 eV) Coeff = .75544E-01 2--> 3 ( 15.22201 eV) Coeff = .99610E+00 +-----------------------------------------------+ 5 Transition energy 17.90971 eV : +-----------------------------------------------+ 1--> 2 ( 19.47640 eV) Coeff = -.99097E+00 1--> 4 ( 41.85260 eV) Coeff = .91640E-01 2--> 1 ( -7.83255 eV) Coeff = -.70990E-01 6

All transitions start from an α orbital and finish in a β orbital. In the case of H2 , the triplet state contains no β electrons. Transition 1 corresponds to the decay from the triplet state to the ground state (2–> 1). 2 -|--

2 ----->

-1

1 -|--

1 -||-

The negative transition energy should not be considered as unrealistic, it is simply an artifact of the fact that we are starting from an excited state. Transition 2 is a triplet combination (1–> 1 - 2–> 2). 2 -|---> 1 -|--

1 ---sqrt(2)

2 -|-1 --|-

1 ---sqrt(2)

2 --|1 -|--

In SF this transition should always be close to zero, since it is the same as the reference state. Transition 3 corresponds to the singlet combination (1–> 1 + 2–> 2). 2 -|---> 1 -|--

1 ---sqrt(2)

2 -|-+ 1 --|-

1 ---sqrt(2)

2 --|1 -|--

Transition 4 corresponds is an excitation which ends in a impure spin state (2–> 3). 3 ---2 -|--

3 --|-->

1 -|--

2 ---1 -|--

The problem here is that SF TDDFT is not generating the partner state. 3 ---2 -|--

3 -|---X->

2 ----

1 -|--

1 --|-

In fact, the rest of the SF transitions (other than decay to the ground state, triplet and singlet combinations and the double excitations) are not pure spin states! This is one of the main drawbacks of the SF approach. Transition 5 corresponds to a double excitation (1–> 2). 2 -|--

2 -||-->

1 -|--

-1 1 ----

Note that in fact it is a double-excitation with respect to the real ground-state of H2 but a single SF excitation with respect to the reference triplet.

7

4

Subroutines used in the present implemenation of SF TDDFT

The main subroutines used by the present implementation of SF-TDDFT are briefly reviewed here. This implementation mainly reuses previously implemented TDDFT routines. However the header for tddft is modified to include the SF dimension spaces. The parts reserved for SF calculations can be identified by the NONCOL logical variable within an IF statement. tddrv.f This is the driver which controls the initialization, calculation of kernels and integrals, the diagonalization, and output. tdini.f This routine is responsible for initializing the CIS space–i.e., calculation of the number of transitions. tdkernel.f Calculates the SF kernel. This involves multiplication by auxiliary function overlap matrices and storage in files. cisptr.f Generate the lookup table for the excitations (NOTE: This routine is a function). vsfnumxc.f Numerical calculation of the kernel and storage in files. numxcsf.f Adds up contributions from different buffers read from disk to form the SF kernel. bldmoeris.f Calculate three-center integrals needed to construct the coupling matrix (called by calcmoeris.f). sfguess.f Construct the initial guess for the SF implentation of Davison diagnolization. tdkrylov.f Initialization and solution of the TDDFT equations using the davidson algorithm. (NOTE: in this routine the SF parts of the code are identified by the character variable OPTION. The values of OPTION for SF are ALPHA-BETA for alpha to beta transitions and BETA-ALPHA for beta to alpha transitions). getsomsf.f Gets the molecular orbital energy differences for the SF part. davidson1.f Calculate the excitations with the davidson algorithm. (NOTE: The same OPTION value is used to follow the SF parts of the code. See tdkrylov.f). tdmatrix.f Initialization and full solution of the TDDFT equations. getkmo.f Calculation of the coupling matrix. kmoini.f Calculation of the coupling matrix. noncolout.f Ouput the results of the SF calculation.

5

Known Problems, Suggested Improvements • The output could be more clear. However this documentation should help considerably. • Some method should be devised to overcome the generation of impure spin states or at least identify them automatically.

8

References [1] M. E. Casida, Time-dependent density-functional response theory for molecules, in Recent Advances in Density Functional Methods, edited by D. P. Chong, page 155, Singapore, World Scientific, 1995. [2] . F. Wang and T. Ziegler, Time-dependent density functional theory based un a noncollinear formulation of the exchange-correlation potential, J. Chem. Phys. 121, 12191 (2004). [3] F. Wang and T. Ziegler, Use of noncollinear exchange-correlation potentials in multiplet resolutions by time-dependent density functional theory, Int. J. Quant. Chem. 106, 2545 (2006). [4] Y. Shao, M. Head-Gordon, and A. I. Krylov, The spin-flip approach within time-dependent density functional theory: Theory and applications to diradicals, J. Chem. Phys 118, 4807 (2003). [5] M. Huix-Rotllant, B. Natarajan, C. M. Wawire, A. Ipatov, and M. E. Casida, in progress, Unpublished, 2009. [6] . F. Wang and T. Ziegler, The performance of time-dependent density functional theory based on a noncollinear exchange-correlation potential in the calculation of excitation energies, J. Chem. Phys. 122, 074109 (2005). [7] A. M. Koester et al., deMon2k, Version 1.8, The deMon Developers, 2005.

9