Ab initio non-relativistic spin dynamics Feizhi Ding, Joshua J. Goings, Michael J. Frisch, and Xiaosong Li Citation: The Journal of Chemical Physics 141, 214111 (2014); doi: 10.1063/1.4902884 View online: http://dx.doi.org/10.1063/1.4902884 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/21?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Ab initio Ehrenfest dynamics J. Chem. Phys. 123, 084106 (2005); 10.1063/1.2008258 Correlations, spin dynamics, defects: the highly frustrated kagomé bilayer Low Temp. Phys. 31, 704 (2005); 10.1063/1.2008131 The ab initio model potential method with the spin-free relativistic scheme by eliminating small components Hamiltonian J. Chem. Phys. 114, 6000 (2001); 10.1063/1.1356735 Spin-free relativistic no-pair ab initio core model potentials and valence basis sets for the transition metal elements Sc to Hg. II J. Chem. Phys. 111, 10436 (1999); 10.1063/1.480398 Spin dynamics in high- T C superconductors AIP Conf. Proc. 483, 207 (1999); 10.1063/1.59619
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Sun, 03 May 2015 18:40:19
THE JOURNAL OF CHEMICAL PHYSICS 141, 214111 (2014)
Ab initio non-relativistic spin dynamics Feizhi Ding,1 Joshua J. Goings,1 Michael J. Frisch,2 and Xiaosong Li1,a) 1 2
Department of Chemistry, University of Washington, Seattle, Washington 98195, USA Gaussian, Inc., 340 Quinnipiac St, Bldg 40, Wallingford, Connecticut 06492, USA
(Received 29 August 2014; accepted 17 November 2014; published online 4 December 2014) Many magnetic materials do not conform to the (anti-)ferromagnetic paradigm where all electronic spins are aligned to a global magnetization axis. Unfortunately, most electronic structure methods cannot describe such materials with noncollinear electron spin on account of formally requiring spin alignment. To overcome this limitation, it is necessary to generalize electronic structure methods and allow each electron spin to rotate freely. Here, we report the development of an ab initio time-dependent non-relativistic two-component spinor (TDN2C), which is a generalization of the time-dependent Hartree-Fock equations. Propagating the TDN2C equations in the time domain allows for the first-principles description of spin dynamics. A numerical tool based on the Hirshfeld partitioning scheme is developed to analyze the time-dependent spin magnetization. In this work, we also introduce the coupling between electron spin and a homogenous magnetic field into the TDN2C framework to simulate the response of the electronic spin degrees of freedom to an external magnetic field. This is illustrated for several model systems, including the spin-frustrated Li3 molecule. Exact agreement is found between numerical and analytic results for Larmor precession of hydrogen and lithium atoms. The TDN2C method paves the way for the ab initio description of molecular spin transport and spintronics in the time domain. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4902884] I. INTRODUCTION
Spin-based technology takes on an ever-increasing role in material design as scientists and engineers seek to utilize and manipulate spins to maximize the quantum yield of solar cells, to process quantum information, and to resolve complex molecular structure. Fundamental to these scientific and technological applications is a detailed understanding of spindependent many-electron dynamics.1–4 To date, most theoretical investigations on molecular spintronics rely on the use of effective model Hamiltonians which depend on experimental parameterization5–8 and focus on collinear spin systems in which all the spin magnetizations are aligned along a single anisotropy axis. While these methods have their merits, they are limited in their description of the time-dependent nature of the spin-dependent many-electron wave function. This is especially true for phenomena driven by strong spin non-collinearity, spin-orbit coupling, or by interaction with intense electromagnetic fields. A proper description of spin dynamics must come from a solution of the first-principles spin-dependent Hamiltonian in the time domain. Most standard electronic structure methods, however, cannot treat systems involving non-collinear spins because they formally require electronic spins to be aligned (anti-) parallel with respect to each other. In both the unrestricted Hartree-Fock (UHF) and spin-density density functional theories (SDFT), the one-electron spin orbitals are required to be eigenfunctions of the spin operator Sˆz . Practically, this means requiring the number of spin-up and spin-down electrons to a) Email: [email protected]
0021-9606/2014/141(21)/214111/8/$30.00
be conserved in the solution of the self-consistent-field (SCF) equations. As a result, spin orbitals are a product of a spatial function with either a spin-up or spin-down component. This description is equivalent to the assumption of collinear spins, i.e., there exists a global spin quantization axis for the whole system, and the magnetization vector m(r) is of the same direction at all points in space. In order to treat non-collinear spins, one needs to retain the full vector form of the magnetization m(r) and allow each spin quantization axis to rotate. This is equivalent to writing the spin orbitals as a superposition of the spin-up and spin-down states. For Hartree-Fock, this leads to the generalized Hartree-Fock (GHF) method,9–13 which is similar in structure to the wave function used in twocomponent relativistic models. For density functional theory (DFT), a generalization of the functionals to account for noncollinear spins is required.14–22 The GHF and non-collinear DFT methods have been used to study the non-collinear spin states of dissociated molecules and geometrically frustrated systems.10–13, 15, 16, 19–21, 23–26 In general, the phenomenon of non-collinear magnetism arises from two different types of interactions: the first comes from nearest-neighbor magnetic exchange interaction and is usually the dominant contribution; the second stems from the anisotropic interaction due to spin-orbital coupling and is typically one order of magnitude weaker for light elements. The first kind of interaction, i.e., the magnetic exchange, comes from the Pauli principle which is embedded in the antisymmetrized electronic wave function, even though in the non-relativistic Hamiltonian the Coulomb repulsion does not explicitly depend on spin. Therefore, a non-collinear non-relativistic treatment is still necessary for
141, 214111-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Sun, 03 May 2015 18:40:19
214111-2
Ding et al.
a proper description of spin magnetizations, especially in cases of spin frustration. It should be noted that although the non-relativistic Hamiltonian commutes with the spin operators S2 and Sz and the exact wave function may possess all the spin symmetries, approximate wave functions, such as a single Slater determinant, need not have these symmetry constraints.9–11, 13, 27 In fact, adding symmetry constraints can only raise the energy of a variationally optimized solution. Thus, by removing some symmetry constraints, it is possible to obtain a lower energy variational solution. In such a situation, one is forced to choose between an approximate wave function that is variationally optimum and a higher energy approximate wave function that contains the proper symmetries of the Hamiltonian. This dilemma was first pointed out by Löwdin, and is known as Löwdin’s symmetry dilemma.28 From the point of view of electronic dynamics, the extension of the collinear DFT and HF to the time domain gives rise to the time-dependent density functional theory (TDDFT) and time-dependent Hartree-Fock (TDHF).29–39 In particular, the real-time formalism of TDDFT/TDHF (RT-TDDFT/RTTDHF) has proven to be a promising method for revealing the physical underpinnings of dynamical multi-electron phenomena.40–44 However, these methods are unable to capture the spin dynamics driven by an external magnetic field, or an internal magnetic field generated by the motion of other electrons within the system. Implicit to the single Slater determinant based RT-TDDFT/RT-TDHF method are the constraints that spin state of each electron cannot change during the time evolution and that they are collinear with respect to each other. These constraints prevent conventional RTTDDFT/RT-TDHF from fully accounting for spin-dependent many-electron dynamics. In order to describe non-collinear spin dynamics, we extend the generalized Hartree-Fock to the time domain, giving rise to the real-time time-dependent non-relativistic two-component spinor (RT-TDN2C) method. We present in this work the mathematical formalism and numerical implementation with case studies of non-collinear spin dynamics in response to an external magnetic field. To the best of our knowledge, this work is the first attempt to address manybody spin dynamics within the full ab initio description. II. METHODOLOGY
In the following derivations, the wave function is represented in a set of atomic orbitals (AO) {χ } with μ, ν, . . . as indices. Primed notations are used for matrices in the AO basis and unprimed notations for matrices in the orthonormal basis. In particular, F and P (F and P) are Fock and density matrices in the AO (orthonormal) basis. F and P are related to their AO basis counterparts through an orthonormal transformation. σ and τ are used as spin variables (σ ∈ {α, β}; τ ∈ {α, β}). All equations are expressed in atomic units. A. Time-dependent non-relativistic two-component spinor formalism
In the TDHF theory, the N-electron wave function is represented by a single Slater determinant ψ(t) composed of N
J. Chem. Phys. 141, 214111 (2014)
time-dependent spin orbitals, {ψk (x, t)} ψ1 (x1 , t) ψ2 (x1 , t) 1 ψ1 (x2 , t) ψ2 (x2 , t) ψ(t) = √ .. .. N ! . . ψ (x , t) ψ (x , t) 1
2
N
N
··· ··· ..
.
···
ψN (x1 , t) ψN (x2 , t) , .. . ψN (xN , t) (1)
where the variable x includes both the spatial coordinate r and the spin coordinate ω for one electron. Applying the Dirac-Frenkel time-dependent variational principle45, 46 to this approximate wave function, we obtain the equations of motion for the single-particle spin orbitals, i
∂ ψ (x, t) = fˆ(t)ψk (x, t), ∂t k
(2)
where fˆ is the Fock operator. Equivalently, Eq. (2) can be written in the Liouville-von Neumann form ∂ γˆ (t) = [fˆ(t), γˆ (t)], (3) i ∂t where γˆ (t) is the reduced single-particle density operator and is defined by γˆ (t) =
N
|ψk (t)ψk (t)|.
(4)
i
In order to allow for non-collinearity, we express the time-dependent spin orbitals as a time-dependent twocomponent spinor, α φk (r, t) ψk (x, t) = , (5) β φk (r, t) β
where the spatial functions {φkα (r, t)}, {φk (r, t)} are expanded in terms of a common set of real AO basis functions {χμ (r)}, α φkα (r, t) = Cμk (t)χμ (r), (6) μ β
φk (r, t) =
β
Cμk (t)χμ (r),
(7)
μ
where the time-dependent expansion coefficients C’s may take complex values. With the definition of the twocomponent spinor, the time-dependent single-particle density operator can be formulated with the following spin-blocked structure: αα γˆ (t) γˆ αβ (t) . (8) γˆ (t) = γˆ βα (t) γˆ ββ (t) In the AO basis, its matrix form can be written as αα αβ (t) P (t) P P (t) = , Pβα (t) Pββ (t)
(9)
where each spin block Pσ τ (t) is given by σ τ Pμν (t) =
N
σ τ∗ Cμi (t) · Cνi (t).
(10)
i
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Sun, 03 May 2015 18:40:19
Ding et al.
214111-3
J. Chem. Phys. 141, 214111 (2014)
The Fock matrix in the AO basis also takes a spin-blocked form analogous to the density matrix αα αβ (t) F (t) F F (t) = , (11) Fβα (t) Fββ (t) where each spin block Fσ τ (t) is given by στ Fσ τ (t) = hσ τ (t) + δσ τ Jαα (t) + Jββ (t) − K (t), (12) where hσ τ is the one-electron Hamiltonian matrix, whose time-dependence arises from the interaction of the spin system with external electromagnetic perturbations (to be discussed in Sec. II C). J and K are Coulomb and exchange matrices, respectively, σ σ (t) Jμν
≡ μ|Jˆσ σ |ν =
N
−1 σ φj (2)χν (1) dr1 dr2 χμ∗ (1)φjσ ∗ (2)r12
j
=
N j
=
σ Cλj (t)
·
σ∗ Cκj (t)
−1 dr1 dr2 χμ∗ (1)χκ∗ (2)r12 χλ (2)χν (1)
The final density-matrix based TDN2C Liouville-von Neumann equation in the orthonormal basis can be written as ∂ Pαα (t) Pαβ (t) i ∂t Pβα (t) Pββ (t)
Fαα (t) Fαβ (t) Pαα (t) Pαβ (t) = , . (17) Fβα (t) Fββ (t) Pβα (t) Pββ (t) Equation (17) is the working equation to describe the time evolution of the non-relativistic two-component spinor. It can be considered as a time-dependent extension of the generalized Hartree-Fock equation. As with most time-dependent electronic structure theories, Eq. (17) can be solved in a response formalism in the frequency domain or propagated explicitly in the time-domain. In Secs. II B–II D, we will present a time-propagation approach and perturbation terms arising from the interaction between a spin system and a static magnetic field, as well as a method designed to analyze spin magnetization in the time domain.
λκ σ σ Pλκ (t) · (μν|κλ) ,
(13)
B. Numerical propagation using modified midpoint and unitary transformation
λκ σ τ Kμν (t)
≡ μ|Kˆ σ τ |ν =
N
−1 χν (2)φjσ (1) dr1 dr2 χμ∗ (1)φjτ ∗ (2)r12
j
=
N j
=
−1 σ τ∗ Cλj (t) · Cκj (t) dr1 dr2 χμ∗ (1)χκ∗ (2)r12 χν (2)χλ (1)
λκ σ τ Pλκ (t) · (μλ|κν) .
(14)
In this work, we use a set of atom-centered atomic orbital basis functions which form a non-orthogonal basis set. The density and Fock matrices are then transformed from the AO basis to an orthonormal basis using the following equations: Pαα (t) Pαβ (t)
Fαα (t)
Fαβ (t)
C† (tn ) · F(tn ) · C(tn ) = (tn ),
(18)
U(tn ) = exp[−i · 2t · F(tn )]
λκ
Pβα (t) Pββ (t) V 0 Pαα (t) Pαβ (t) VT = · · 0 V Pβα (t) Pββ (t) 0
To obtain a time-evolution dynamics of a spin system, the time-dependent equation for non-relativistic two-component spinors (Eq. (17)) is integrated with a modified midpoint and unitary transformation (MMUT) algorithm.40, 42 In the MMUT method, the time-evolution operator is a unitary transformation matrix U(tn ) that is constructed from the eigenvectors C(tn ) and eigenvalues (tn ) of the generalized Fock matrix at time tn ,
0 VT
,
(15)
Fβα (t) Fββ (t) V−T 0 Fαα (t) Fαβ (t) V−1 0 = · . · 0 V−T Fβα (t) Fββ (t) 0 V−1 (16) The transformation matrix V = S1/2 in the Löwdin orthonormalization method, or may be obtained by the Cholesky decomposition S = VT V, where S is the overlap matrix of AO basis functions, Sμν = μ|ν.
= C(tn ) · exp[−i · 2t · (tn )] · C† (tn ),
(19)
where t is the time step. The density matrix is then propagated from time tn−1 to tn+1 using the time-evolution operator U(tn ), P(tn+1 ) = U(tn ) · P(tn−1 ) · U† (tn ).
(20)
Note that the Fock and density matrices in Eqs. (18)–(20) are assumed to have a spin-blocked structure (Eqs. (15) and (16)). By computing the Fock matrix at the midpoint of the step, the MMUT method accounts for linear changes in the density matrix during the time step. The MMUT algorithm can be used with a large step size, while still maintaining very good control of numerical noise and integration errors.40 C. Two-component spinor in the presence of a static and uniform magnetic field
The formalism of the non-relativistic two-component spinor allows realistic simulations of non-collinear spin dynamics within the ab initio framework. In the non-relativistic framework, the interaction of an electron spin with external electromagnetic field is described by the so-called
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Sun, 03 May 2015 18:40:19
214111-4
Ding et al.
J. Chem. Phys. 141, 214111 (2014)
Schrödinger-Pauli Hamiltonian 1 hˆ Pauli = [σ · (p + A)]2 − U, 2
(21)
where A and U are the vector potential and scalar potential of the electromagnetic field. p = −i∇ is the canonical momentum operator. Embedding the Schrödinger-Pauli Hamiltonian into the time-dependent Hartree-Fock equation gives rise to the spindependent non-relativistic two-component spinor, i.e., the time-dependent Pauli-Hartree-Fock (PHF) equation i
∂ ψ (x, t) = fˆPauli (t)ψk (x, t). ∂t k
(22)
In this work, we consider a many-electron spin system interacting with an external static, uniform magnetic field B. Upon choice of the Coulomb gauge ∇ · A = 0 and A = 12 B × r, the one electron Hamiltonian of the PHF equation can be written as 1 1 1 hˆ Pauli = hˆ 0 (r) + σ · B + B · l + (B × r)2 , 2 2 8
(23)
where hˆ 0 (r) is the field-free one-electron Hamiltonian. l = r × p is the orbital angular momentum. The last two terms represent the interaction of the orbital angular momentum with the magnetic field and the higher order term that is quadratic in the field strength. These two terms will affect the time evolution of the spatial distribution of the electrons and thus affect the total energy of the system. For weak magnetic fields (e.g.,
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Sun, 03 May 2015 18:40:19
THE JOURNAL OF CHEMICAL PHYSICS 141, 214111 (2014)
Ab initio non-relativistic spin dynamics Feizhi Ding,1 Joshua J. Goings,1 Michael J. Frisch,2 and Xiaosong Li1,a) 1 2
Department of Chemistry, University of Washington, Seattle, Washington 98195, USA Gaussian, Inc., 340 Quinnipiac St, Bldg 40, Wallingford, Connecticut 06492, USA
(Received 29 August 2014; accepted 17 November 2014; published online 4 December 2014) Many magnetic materials do not conform to the (anti-)ferromagnetic paradigm where all electronic spins are aligned to a global magnetization axis. Unfortunately, most electronic structure methods cannot describe such materials with noncollinear electron spin on account of formally requiring spin alignment. To overcome this limitation, it is necessary to generalize electronic structure methods and allow each electron spin to rotate freely. Here, we report the development of an ab initio time-dependent non-relativistic two-component spinor (TDN2C), which is a generalization of the time-dependent Hartree-Fock equations. Propagating the TDN2C equations in the time domain allows for the first-principles description of spin dynamics. A numerical tool based on the Hirshfeld partitioning scheme is developed to analyze the time-dependent spin magnetization. In this work, we also introduce the coupling between electron spin and a homogenous magnetic field into the TDN2C framework to simulate the response of the electronic spin degrees of freedom to an external magnetic field. This is illustrated for several model systems, including the spin-frustrated Li3 molecule. Exact agreement is found between numerical and analytic results for Larmor precession of hydrogen and lithium atoms. The TDN2C method paves the way for the ab initio description of molecular spin transport and spintronics in the time domain. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4902884] I. INTRODUCTION
Spin-based technology takes on an ever-increasing role in material design as scientists and engineers seek to utilize and manipulate spins to maximize the quantum yield of solar cells, to process quantum information, and to resolve complex molecular structure. Fundamental to these scientific and technological applications is a detailed understanding of spindependent many-electron dynamics.1–4 To date, most theoretical investigations on molecular spintronics rely on the use of effective model Hamiltonians which depend on experimental parameterization5–8 and focus on collinear spin systems in which all the spin magnetizations are aligned along a single anisotropy axis. While these methods have their merits, they are limited in their description of the time-dependent nature of the spin-dependent many-electron wave function. This is especially true for phenomena driven by strong spin non-collinearity, spin-orbit coupling, or by interaction with intense electromagnetic fields. A proper description of spin dynamics must come from a solution of the first-principles spin-dependent Hamiltonian in the time domain. Most standard electronic structure methods, however, cannot treat systems involving non-collinear spins because they formally require electronic spins to be aligned (anti-) parallel with respect to each other. In both the unrestricted Hartree-Fock (UHF) and spin-density density functional theories (SDFT), the one-electron spin orbitals are required to be eigenfunctions of the spin operator Sˆz . Practically, this means requiring the number of spin-up and spin-down electrons to a) Email: [email protected]
0021-9606/2014/141(21)/214111/8/$30.00
be conserved in the solution of the self-consistent-field (SCF) equations. As a result, spin orbitals are a product of a spatial function with either a spin-up or spin-down component. This description is equivalent to the assumption of collinear spins, i.e., there exists a global spin quantization axis for the whole system, and the magnetization vector m(r) is of the same direction at all points in space. In order to treat non-collinear spins, one needs to retain the full vector form of the magnetization m(r) and allow each spin quantization axis to rotate. This is equivalent to writing the spin orbitals as a superposition of the spin-up and spin-down states. For Hartree-Fock, this leads to the generalized Hartree-Fock (GHF) method,9–13 which is similar in structure to the wave function used in twocomponent relativistic models. For density functional theory (DFT), a generalization of the functionals to account for noncollinear spins is required.14–22 The GHF and non-collinear DFT methods have been used to study the non-collinear spin states of dissociated molecules and geometrically frustrated systems.10–13, 15, 16, 19–21, 23–26 In general, the phenomenon of non-collinear magnetism arises from two different types of interactions: the first comes from nearest-neighbor magnetic exchange interaction and is usually the dominant contribution; the second stems from the anisotropic interaction due to spin-orbital coupling and is typically one order of magnitude weaker for light elements. The first kind of interaction, i.e., the magnetic exchange, comes from the Pauli principle which is embedded in the antisymmetrized electronic wave function, even though in the non-relativistic Hamiltonian the Coulomb repulsion does not explicitly depend on spin. Therefore, a non-collinear non-relativistic treatment is still necessary for
141, 214111-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Sun, 03 May 2015 18:40:19
214111-2
Ding et al.
a proper description of spin magnetizations, especially in cases of spin frustration. It should be noted that although the non-relativistic Hamiltonian commutes with the spin operators S2 and Sz and the exact wave function may possess all the spin symmetries, approximate wave functions, such as a single Slater determinant, need not have these symmetry constraints.9–11, 13, 27 In fact, adding symmetry constraints can only raise the energy of a variationally optimized solution. Thus, by removing some symmetry constraints, it is possible to obtain a lower energy variational solution. In such a situation, one is forced to choose between an approximate wave function that is variationally optimum and a higher energy approximate wave function that contains the proper symmetries of the Hamiltonian. This dilemma was first pointed out by Löwdin, and is known as Löwdin’s symmetry dilemma.28 From the point of view of electronic dynamics, the extension of the collinear DFT and HF to the time domain gives rise to the time-dependent density functional theory (TDDFT) and time-dependent Hartree-Fock (TDHF).29–39 In particular, the real-time formalism of TDDFT/TDHF (RT-TDDFT/RTTDHF) has proven to be a promising method for revealing the physical underpinnings of dynamical multi-electron phenomena.40–44 However, these methods are unable to capture the spin dynamics driven by an external magnetic field, or an internal magnetic field generated by the motion of other electrons within the system. Implicit to the single Slater determinant based RT-TDDFT/RT-TDHF method are the constraints that spin state of each electron cannot change during the time evolution and that they are collinear with respect to each other. These constraints prevent conventional RTTDDFT/RT-TDHF from fully accounting for spin-dependent many-electron dynamics. In order to describe non-collinear spin dynamics, we extend the generalized Hartree-Fock to the time domain, giving rise to the real-time time-dependent non-relativistic two-component spinor (RT-TDN2C) method. We present in this work the mathematical formalism and numerical implementation with case studies of non-collinear spin dynamics in response to an external magnetic field. To the best of our knowledge, this work is the first attempt to address manybody spin dynamics within the full ab initio description. II. METHODOLOGY
In the following derivations, the wave function is represented in a set of atomic orbitals (AO) {χ } with μ, ν, . . . as indices. Primed notations are used for matrices in the AO basis and unprimed notations for matrices in the orthonormal basis. In particular, F and P (F and P) are Fock and density matrices in the AO (orthonormal) basis. F and P are related to their AO basis counterparts through an orthonormal transformation. σ and τ are used as spin variables (σ ∈ {α, β}; τ ∈ {α, β}). All equations are expressed in atomic units. A. Time-dependent non-relativistic two-component spinor formalism
In the TDHF theory, the N-electron wave function is represented by a single Slater determinant ψ(t) composed of N
J. Chem. Phys. 141, 214111 (2014)
time-dependent spin orbitals, {ψk (x, t)} ψ1 (x1 , t) ψ2 (x1 , t) 1 ψ1 (x2 , t) ψ2 (x2 , t) ψ(t) = √ .. .. N ! . . ψ (x , t) ψ (x , t) 1
2
N
N
··· ··· ..
.
···
ψN (x1 , t) ψN (x2 , t) , .. . ψN (xN , t) (1)
where the variable x includes both the spatial coordinate r and the spin coordinate ω for one electron. Applying the Dirac-Frenkel time-dependent variational principle45, 46 to this approximate wave function, we obtain the equations of motion for the single-particle spin orbitals, i
∂ ψ (x, t) = fˆ(t)ψk (x, t), ∂t k
(2)
where fˆ is the Fock operator. Equivalently, Eq. (2) can be written in the Liouville-von Neumann form ∂ γˆ (t) = [fˆ(t), γˆ (t)], (3) i ∂t where γˆ (t) is the reduced single-particle density operator and is defined by γˆ (t) =
N
|ψk (t)ψk (t)|.
(4)
i
In order to allow for non-collinearity, we express the time-dependent spin orbitals as a time-dependent twocomponent spinor, α φk (r, t) ψk (x, t) = , (5) β φk (r, t) β
where the spatial functions {φkα (r, t)}, {φk (r, t)} are expanded in terms of a common set of real AO basis functions {χμ (r)}, α φkα (r, t) = Cμk (t)χμ (r), (6) μ β
φk (r, t) =
β
Cμk (t)χμ (r),
(7)
μ
where the time-dependent expansion coefficients C’s may take complex values. With the definition of the twocomponent spinor, the time-dependent single-particle density operator can be formulated with the following spin-blocked structure: αα γˆ (t) γˆ αβ (t) . (8) γˆ (t) = γˆ βα (t) γˆ ββ (t) In the AO basis, its matrix form can be written as αα αβ (t) P (t) P P (t) = , Pβα (t) Pββ (t)
(9)
where each spin block Pσ τ (t) is given by σ τ Pμν (t) =
N
σ τ∗ Cμi (t) · Cνi (t).
(10)
i
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Sun, 03 May 2015 18:40:19
Ding et al.
214111-3
J. Chem. Phys. 141, 214111 (2014)
The Fock matrix in the AO basis also takes a spin-blocked form analogous to the density matrix αα αβ (t) F (t) F F (t) = , (11) Fβα (t) Fββ (t) where each spin block Fσ τ (t) is given by στ Fσ τ (t) = hσ τ (t) + δσ τ Jαα (t) + Jββ (t) − K (t), (12) where hσ τ is the one-electron Hamiltonian matrix, whose time-dependence arises from the interaction of the spin system with external electromagnetic perturbations (to be discussed in Sec. II C). J and K are Coulomb and exchange matrices, respectively, σ σ (t) Jμν
≡ μ|Jˆσ σ |ν =
N
−1 σ φj (2)χν (1) dr1 dr2 χμ∗ (1)φjσ ∗ (2)r12
j
=
N j
=
σ Cλj (t)
·
σ∗ Cκj (t)
−1 dr1 dr2 χμ∗ (1)χκ∗ (2)r12 χλ (2)χν (1)
The final density-matrix based TDN2C Liouville-von Neumann equation in the orthonormal basis can be written as ∂ Pαα (t) Pαβ (t) i ∂t Pβα (t) Pββ (t)
Fαα (t) Fαβ (t) Pαα (t) Pαβ (t) = , . (17) Fβα (t) Fββ (t) Pβα (t) Pββ (t) Equation (17) is the working equation to describe the time evolution of the non-relativistic two-component spinor. It can be considered as a time-dependent extension of the generalized Hartree-Fock equation. As with most time-dependent electronic structure theories, Eq. (17) can be solved in a response formalism in the frequency domain or propagated explicitly in the time-domain. In Secs. II B–II D, we will present a time-propagation approach and perturbation terms arising from the interaction between a spin system and a static magnetic field, as well as a method designed to analyze spin magnetization in the time domain.
λκ σ σ Pλκ (t) · (μν|κλ) ,
(13)
B. Numerical propagation using modified midpoint and unitary transformation
λκ σ τ Kμν (t)
≡ μ|Kˆ σ τ |ν =
N
−1 χν (2)φjσ (1) dr1 dr2 χμ∗ (1)φjτ ∗ (2)r12
j
=
N j
=
−1 σ τ∗ Cλj (t) · Cκj (t) dr1 dr2 χμ∗ (1)χκ∗ (2)r12 χν (2)χλ (1)
λκ σ τ Pλκ (t) · (μλ|κν) .
(14)
In this work, we use a set of atom-centered atomic orbital basis functions which form a non-orthogonal basis set. The density and Fock matrices are then transformed from the AO basis to an orthonormal basis using the following equations: Pαα (t) Pαβ (t)
Fαα (t)
Fαβ (t)
C† (tn ) · F(tn ) · C(tn ) = (tn ),
(18)
U(tn ) = exp[−i · 2t · F(tn )]
λκ
Pβα (t) Pββ (t) V 0 Pαα (t) Pαβ (t) VT = · · 0 V Pβα (t) Pββ (t) 0
To obtain a time-evolution dynamics of a spin system, the time-dependent equation for non-relativistic two-component spinors (Eq. (17)) is integrated with a modified midpoint and unitary transformation (MMUT) algorithm.40, 42 In the MMUT method, the time-evolution operator is a unitary transformation matrix U(tn ) that is constructed from the eigenvectors C(tn ) and eigenvalues (tn ) of the generalized Fock matrix at time tn ,
0 VT
,
(15)
Fβα (t) Fββ (t) V−T 0 Fαα (t) Fαβ (t) V−1 0 = · . · 0 V−T Fβα (t) Fββ (t) 0 V−1 (16) The transformation matrix V = S1/2 in the Löwdin orthonormalization method, or may be obtained by the Cholesky decomposition S = VT V, where S is the overlap matrix of AO basis functions, Sμν = μ|ν.
= C(tn ) · exp[−i · 2t · (tn )] · C† (tn ),
(19)
where t is the time step. The density matrix is then propagated from time tn−1 to tn+1 using the time-evolution operator U(tn ), P(tn+1 ) = U(tn ) · P(tn−1 ) · U† (tn ).
(20)
Note that the Fock and density matrices in Eqs. (18)–(20) are assumed to have a spin-blocked structure (Eqs. (15) and (16)). By computing the Fock matrix at the midpoint of the step, the MMUT method accounts for linear changes in the density matrix during the time step. The MMUT algorithm can be used with a large step size, while still maintaining very good control of numerical noise and integration errors.40 C. Two-component spinor in the presence of a static and uniform magnetic field
The formalism of the non-relativistic two-component spinor allows realistic simulations of non-collinear spin dynamics within the ab initio framework. In the non-relativistic framework, the interaction of an electron spin with external electromagnetic field is described by the so-called
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Sun, 03 May 2015 18:40:19
214111-4
Ding et al.
J. Chem. Phys. 141, 214111 (2014)
Schrödinger-Pauli Hamiltonian 1 hˆ Pauli = [σ · (p + A)]2 − U, 2
(21)
where A and U are the vector potential and scalar potential of the electromagnetic field. p = −i∇ is the canonical momentum operator. Embedding the Schrödinger-Pauli Hamiltonian into the time-dependent Hartree-Fock equation gives rise to the spindependent non-relativistic two-component spinor, i.e., the time-dependent Pauli-Hartree-Fock (PHF) equation i
∂ ψ (x, t) = fˆPauli (t)ψk (x, t). ∂t k
(22)
In this work, we consider a many-electron spin system interacting with an external static, uniform magnetic field B. Upon choice of the Coulomb gauge ∇ · A = 0 and A = 12 B × r, the one electron Hamiltonian of the PHF equation can be written as 1 1 1 hˆ Pauli = hˆ 0 (r) + σ · B + B · l + (B × r)2 , 2 2 8
(23)
where hˆ 0 (r) is the field-free one-electron Hamiltonian. l = r × p is the orbital angular momentum. The last two terms represent the interaction of the orbital angular momentum with the magnetic field and the higher order term that is quadratic in the field strength. These two terms will affect the time evolution of the spatial distribution of the electrons and thus affect the total energy of the system. For weak magnetic fields (e.g.,
