Implementation of Analytical Energy Gradient of Spin-Dependent

Subscriber access provided by Purdue University Libraries

Article

Implementation of analytical energy gradient of spin-dependent general Hartree-Fock method based on infinite-order Douglas– Kroll–Hess relativistic Hamiltonian with local unitary transformation Yuya Nakajima, Junji Seino, and Hiromi Nakai J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.5b00928 • Publication Date (Web): 05 Apr 2016 Downloaded from http://pubs.acs.org on April 9, 2016

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

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

Page 1 of 35

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

Implementation of analytical energy gradient of spin-dependent general Hartree-Fock method based on the infinite-order Douglas-Kroll-Hess relativistic Hamiltonian with local unitary transformation Yuya Nakajima,† Junji Seino,‡ and Hiromi Nakai∗,†,‡,¶,§ †Department of Chemistry and Biochemistry, School of Advanced Science and Engineering, Waseda University, Tokyo 169-8555, Japan ‡Research Institute for Science and Engineering, Waseda University, Tokyo 169-8555, Japan ¶CREST, Japan Science and Technology Agency, 4-1-8 Honcho, Kawaguchi, Saitama 332-0012, Japan §Elements Strategy Initiative for Catalysts and Batteries (ESICB), Kyoto University, Katsura, Kyoto 615-8520, Japan E-mail: [email protected] Phone: +81 (0)3 5286 3452. Fax: +81 (0)3 3205 2504

1 ACS Paragon Plus Environment


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

Abstract An analytical energy gradient for the spin-dependent general Hartree-Fock method based on the infinite-order Douglas–Kroll–Hess (IODKH) method was developed. To treat realistic systems, the local unitary transformation (LUT) scheme was employed both in energy and energy gradient calculations. The present energy gradient method was numerically assessed to investigate the accuracy in several diatomic molecules containing fifth- and sixth-period elements and to examine the efficiency in one-, two, and three-dimensional silver clusters. To arrive at a practical calculation, we also determined the geometrical parameters of fac-Ir(ppy)3 and investigated the efficiency. The numerical results confirmed that the present method describes a highly accurate relativistic effect with high efficiency. The present method can be a powerful scheme to determine geometries of large molecules, including heavy-element atoms.


Page 2 of 35

Page 3 of 35

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


1

Introduction

A practical computational method that includes high-level relativistic effects has become one of the attracting goals for theoretical chemists. A fundamental treatment is based on a four-component (4c) Dirac theory, which adopts a one-particle Dirac operator. For many electron systems, the 4c theory normally adopts the non-relativistic (NR) Coulomb interaction 1 or a lower order of quantum electrodynamic effects such as the Breit 2 or Gaunt 3 interaction. Although the 4c relativistic calculations succeeded in accurately predicting molecular geometries 4 and electronic/magnetic properties 5,6 for small molecules, including at most 50 atoms, their application to more realistic systems is difficult. Variational collapse arising from the existence of small-component Dirac spinors and/or positronic-state spinors has been a computational problem in the 4c calculations. The kinetic balance condition, which is a technique to resolve the variational problem, requires a large number of basis functions for the small-component spinors in comparison to those for the large-component spinors. To address the expensive computational cost in the 4c calculations, the basis spinor scheme has been proposed by Nakajima et al., 7 which uses common expansion coefficients in each spinor. Recently, Kelley and Shiozaki have improved the algorithm using density fitting technique combined with a highly parallelized implementation for a large-scale calculation and have extended the algorithm to an analytical nuclear gradients. 8,9 An alternative approach is a two-component (2c) relativistic theory, which treats only electronic-state (or large-component) information. Recent development in 2c treatment for many-electron systems has exhibited high accuracy, even when compared to the 4c theory. In particular, the Dirac-exact approaches, which give equivalent results to the 4c treatments, have been proposed for the one-electron Dirac Hamiltonian such as the normalized elimination of small components (NESC), 10 the exact two-component (X2C), 11 and the infinite-order Douglas-Kroll-Hess (IODKH) 12 methods, and the many-electron DiracCoulomb(-Gaunt) ones. 13–18 Expensive computations for the decoupling between electronic and positronic states were also accelerated by introducing the approximate local transfor3 ACS Paragon Plus Environment


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

mations. 19–26 Among the approximations, the local unitary transformation (LUT) by our group and local approximation to the unitary decoupling transformation (DLU) by Reiher’s group can accurately and efficiently reproduce the results of the original 2c treatments since they include not only the diagonal relativistic components but also off-diagonal ones between adjacent atoms for the one-electron Dirac Hamiltonian. The LUT scheme, furthermore, was successfully extended to the many-electron Dirac-Coulomb Hamiltonian. 22,23 Analytical energy gradient schemes for relativistic theories are essential for the geometry optimization of heavy-element systems. The relativistic effects are roughly categorized into spin-free (SF) and spin-dependent (SD) effects. The SF effect directly contributes to bond contractions/elongations due to the relativistic shrinking/expansion of molecular orbitals. 4,27–29 The SD effect, which indirectly contributes to molecular structures due to energy-level splittings of degenerate orbitals, is non-negligible in heavy-element chemistry. For example, the bond lengths are shortened/elongated by approximately 0.1 ˚ A in At2 30 and 0.2 ˚ A in both (113)H and (115)H 31 by including the SD relativistic effect. A contribution of a two-electron spin-orbit coupling is significant to describe molecular properties such as magnetic shielding constants, 32,33 excitation energies, 34–36 and spectroscopic constants. 37 The analytical energy gradient schemes for the SF relativistic effect have been made available in Dirac-exact 2c approaches such as NESC, 38–40 X2C, 41,42 and IODKH, 43 as well as in the approximate 2c approaches such as the regular approximation (RA) schemes, 44–46 the relativistic ESC (RESC) method, 47 and the finite-order DKH methods. 48–54 Note that these three Dirac-exact 2c approaches give identical results although they need different procedures to obtain energies and properties. Recently, the analytical energy gradient for the 2c NESC method with a screened nucleus potential using effective nuclear charges has been proposed by Zou et al. 55 The scheme has shown that a contribution of a spin-orbit coupling to bond lengths is small in general, but relatively large for van der Waals complexes. The progress of the analytical energy derivative in the relativistic quantum chemistry has been well described in recent review papers. 40,56,57


Page 4 of 35

Page 5 of 35

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


From the computational point of view, the SF effect is more easily taken into account compared to the SD effect because modifications from the non-relativistic (NR) quantum chemical program are limited. In fact, the SF relativistic theory for the one-electron Dirac Hamiltonian requires only the implementation related to the one-electron integrals. On the contrary, the implementation of the SD relativistic theory needs a considerable change of the NR program. For example, many variables such as orbital coefficients, one- and two-electron integrals for the SD operators, and Fock matrices, should be treated as complex numbers. In consequence, the quantum chemical programs that can deal with the analytical energy gradient of the SD relativistic theories are very limited. In the present study, we implemented a computational program for the analytical energy gradient, as well as the energy at the SD relativistic level. Our previous study 43 derived the formulation of the analytical energy gradient for the infinite-order Douglas–Kroll–Hess (IODKH) method with/without the LUT treatment, but we implemented only the SF computation. The local approximation schemes are effective to describe chemical bonds in both light- and heavy-element systems, especially for comparatively strong bonds such as covalent bond, metallic bond, and ionic bond because the schemes reproduce the total energy obtained by corresponding conventional methods will less than 1 millihartree deviations. Note that an extension of the approximate local transformations to other Dirac-exact methods are possible, whose derivation and implementation might be easy. The IODKH method denoted here is equivalent to the infinite-order two-component (IOTC) 12 method or the Barysz–Sadlej–Snijders (BSS) 58,59 method. The present paper is organized as follows. Theory and implementation are explained in Sec. 2. Numerical assessments in terms of both accuracy and efficiency are presented in Sec. 3. Finally, concluding remarks are given in Sec. 4.



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

2

Page 6 of 35

Theory and Implementation

2.1

Analytical energy gradient of general Hartree-Fock method

An implementation of the SD relativistic formalism in a self-consistent field (SCF) calculation requires considerable modifications from the standard NR program because the SD operators include Pauli spin matrices in a complex number, which mix α- and β-spin components. To treat the SD operator adequately, a general Hartree-Fock (GHF) method is adopted in the present paper. The present GHF schemes can be extended to post-HF methods. For the extension, the energy and its gradient expressions are close to those of NR, although these are evaluated in a complex number. This subsection provides the expressions of the GHF energy and its analytical gradient, as well as the difference between the SF and SD formalisms in the implementation. Since the SD relativistic formalism in the GHF method does not allow us to treat α- and β-spin separately, spinors are defined with a linear combination of α- and β-spin orbitals as

ϕi =

∑

α Cµi χµ α

µ

+

∑

{α,β} β Cµi χµ β

µ

=

∑∑ ω

ω Cµi χµ ω,

(1)

µ

where C is the complex spinor coefficient matrix, χ is the basis function, and ω denotes the spin function and indices of spins. The HF energy, E HF , in a general form is expressed as

E HF =

∑ i

[ϕi | h | ϕi ] +

1∑ ([ϕi ϕj | ϕi ϕj ] − [ϕi ϕj | ϕj ϕi ]) + Vnuc , 2 i,j

(2)

where i and j are indices of spinors. The first, second, and third terms of Eq. (2) are the one-electron Hamiltonian term, the Coulomb and exchange interactions, and the nuclear repulsion term, respectively. The one-electron relativistic Hamiltonian matrix, which contains


Page 7 of 35

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


a spin operator through the Pauli spin matrices, is given by ∑

{α,β}

[ϕi | h | ϕi ] =

∑∑ ω,ω ′

i

′

ωω Dµν [χµ ω | h | χν ω ′ ] ,

(3)

µ,ν

where D is the density matrix defined as ′

ωω Dµν =

∑

′

ω ω Ciµ Cνi .

(4)

i

ω ω Here, Ciµ . Likewise, the Coulomb and exchange terms in is a conjugate transpose of Cµi

GHF are given by

and

{α,β} {α,β} 1∑ 1 ∑ ∑ ∑ ωω′ τ τ ′ [ϕi ϕj | ϕi ϕj ] = Dµν Dλρ [χµ ωχν ω ′ | χλ τ χρ τ ′ ] 2 i,j 2 ω,ω′ τ,τ ′ µ,ν,λ,ρ

(5)

{α,β} {α,β} 1∑ 1 ∑ ∑ ∑ ωτ τ ′ ω′ D D [χµ ωχν τ ′ | χλ τ χρ ω ′ ] . [ϕi ϕj | ϕj ϕi ] = 2 i,j 2 ω,ω′ τ,τ ′ µ,ν,λ,ρ µλ νρ

(6)

In the case of the two-electron Coulomb operator as a two-electron interaction, summations of Eqs. (5) and (6) are rewritten as {α,β} {α,β} ) 1 ∑ ∑ ∑ ( ωω′ τ τ ′ ωτ τ ′ ω ′ Dµν Dλρ δωω′ δτ τ ′ − Dµλ Dνρ δωτ ′ δτ ω′ (χµ χν | χλ χρ ) , 2 ω,ω′ τ,τ ′ µ,ν,λ,ρ

(7)

where δ is the Kronecker delta and (χµ χν | χλ χρ ) is the two-electron integral represented in basis function space. Finally, the GHF energy can be expressed as E HF =

∑∑

′

ωω Dµν [χµ ω | h | χν ω ′ ]

ω,ω ′ µ,ν

+

) 1 ∑ ∑ ∑ ( ωω′ τ τ ′ ωτ τ ′ ω ′ Dµν Dλρ δωω′ δτ τ ′ − Dµλ Dνρ δωτ ′ δτ ω′ (χµ χν | χλ χρ ) 2 ω,ω′ τ,τ ′ µ,ν,λ,ρ

+ Vnuc .


(8)


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

Page 8 of 35

Note that the summations over the spin indices, ω and τ , are added to the energy expression in an NR or SF relativistic formalism, which show the mixing of α- and β-spin components by the SD operators. The analytical gradient of the GHF energy with respect to a coordinate of a nucleus RA is derived by directly differentiating each term in Eq. (8) as follows: HF

∂E ∂RA

(

) ( ) ωω ′ ωω ′ ∑ ∑ ∂h ∂D ′ µν µν ωω ′ = Dµν + hωω µν ∂RA ∂RA ω,ω ′ µ,ν ω,ω ′ µ,ν ) ∂ 1 ∑ ∑ ∑ ( ωω′ τ τ ′ ωτ τ ′ ω ′ Dνρ δωτ ′ δτ ω′ + Dµν Dλρ δωω′ δτ τ ′ − Dµλ (χµ χν | χλ χρ ) 2 ω,ω′ τ,τ ′ µ,ν,λ,ρ ∂RA ( ) ωω ′ 1 ∑ ∑ ∑ ∂Dµν ′ ′ ′ ωτ τ ω Dνρ δωτ ′ δτ ω′ (χµ χν | χλ χρ ) + Dτ τ δωω′ δτ τ ′ − Dµλ 2 ω,ω′ τ,τ ′ µ,ν,λ,ρ ∂RA λρ ( ) ττ′ ∂D 1 ∑∑ ∑ ′ ′ ′ λρ ωω ωτ τ ω Dµν + δωω′ δτ τ ′ − Dµλ Dνρ δωτ ′ δτ ω′ (χµ χν | χλ χρ ) 2 ω,ω′ τ,τ ′ µ,ν,λ,ρ ∂RA ( ) ωτ ∂Dµλ 1 ∑∑ ∑ ωω ′ τ τ ′ τ ′ ω′ Dµν Dλρ δωω′ δτ τ ′ − D δωτ ′ δτ ω′ (χµ χν | χλ χρ ) + 2 ω,ω′ τ,τ ′ µ,ν,λ,ρ ∂RA νρ ( ) τ ′ ω′ ∂D 1 ∑∑ ∑ ′ ′ νρ ωω ττ ωτ Dµν Dλρ δωω′ δτ τ ′ − Dµλ δωτ ′ δτ ω′ (χµ χν | χλ χρ ) + 2 ω,ω′ τ,τ ′ µ,ν,λ,ρ ∂RA ∑∑

+

(9)

∂Vnuc . ∂RA

Here, a derivative of an orthogonal condition for the spinor at i = j is written as ∑ µ,ν

(

ω ∂Ciµ ∂RA

) ωω ′

ω′

Sµν Cνi +

∑ µ,ν

( ω Ciµ

ωω ∂Sµν ∂RA

′

) ω′ Cνi

+

∑ µ,ν


ω ωω ′ Ciµ Sµν

(

′

ω ∂Cνi ∂RA

) = 0,

(10)

Page 9 of 35

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


where S is the overlap matrix. Using the relation in Eq. (10), Eq. (9) is rewritten as ) ′ ∂hωω µν ∂RA ) ∂ 1 ∑ ∑ ∑ ( ωω′ τ τ ′ ωτ τ ′ ω ′ Dµν Dλρ δωω′ δτ τ ′ − Dµλ Dνρ δωτ ′ δτ ω′ (χµ χν | χλ χρ ) (11) + 2 ω,ω′ τ,τ ′ µ,ν,λ,ρ ∂RA ( ) ωω ′ ∑∑ ∂S ∂Vnuc µν ωω ′ + Wµν δωω′ + , ∂RA ∂RA ω,ω ′ µ,ν

∂E HF ∑ ∑ ωω′ = Dµν ∂RA ω,ω ′ µ,ν

(

where W is the energy-weighted density matrix with an orbital energy εi , which is defined as ′

ωω Wµν =

∑

′

ω ω εi Ciµ Cνi .

(12)

i

As well as the energy representation, the summations over the spin indices are added in the expression of analytical energy gradient.

2.2

Analytical energy gradient for the LUT-IODKH method

This subsection gives the derivative of the 2c IODKH Hamiltonian, including the SD operator with/without the LUT scheme with respect to nucleus coordinates from the viewpoint of implementation. The details of derivation were explained in our previous paper, 43 which only includes the implementation of the SF part. In the 2c IODKH method, because the matrix element of h+ with respect to basis function {χ} is difficult to evaluate directly, a matrix transformation scheme and a resolution-ofidentity technique can be adopted. 60 The technique alters the space for the Hamiltonian from the basis function space {χ} to a p2 -space {k}. Here, p is the momentum operator. In the gradient calculation, the direct differentiation of the 2c IODKH Hamitonian is expressed



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

Page 10 of 35

as ) ∂X † h+ 2 ({k})X d ∂RA ( + ) ∂h2 ({k}) † X† d +d X ∂RA ( ) ∂X† † + + d Xh2 ({k}) d, ∂RA

h+ 2 ({χ}) =d† ∂RA

(

(13)

where d denotes the contraction coefficient matrix and X is the transformation matrix, which includes the transformation from a primitive function {φ} space to an orthonormalized {φ} space and from the orthonormalized space to the {k} space. Note that {k} is expressed by a linear combination of {φ} as |k⟩ =

∑

Crk |φr ⟩ ,

(14)

r

where the coefficient Crk is determined by a diagonalization of matrix elements ⟨φr | p2 | φs ⟩. The resulting (pseudo) eigenvectors satisfy the condition ⟨ 2 ′⟩ k p k = ωk δkk′ ,

(15)

where ωk is a (pseudo) eigenvalue. Note that the gradient of the contraction coefficient d is zero because the coefficient only contains the atomic information. In addition, the analytical gradient of the space transformation matrix ∂X/∂RA is obtained by the procedure described 2 in the previous article. 43,48 An explicit gradient expression of h+ 2 ({k}) in the p -space is given

by a straight-forward differentiation of the IODKH Hamiltonian and its factorization. } { ∂ h+ 2 ({k}) kk′ ∂RA ( ) ( ) ( ) ∑ ∂Ω† ′′ ∂Gk′′ k′′′ ∂Ωk′′′ k′ † † kk = Gk′′ k′′′ Ωk′′′ k′ + Ωkk′′ Ωk′′′ k′ + Ωkk′′ Gk′′ k′′′ , ∂R ∂R ∂R A A A ′′ ′′′ k ,k


(16)

Page 11 of 35

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


with Gkk′

( ⟨ ′⟩ ∑ ⟨ 2 † = k p b12 k + Ykk′′ k ′′ +

∑ (

k′′ ,k′′′

† Mkk ′′

′′

⟩ ) ′′′ p2 12 k Yk′′′ k′ 1 − ep

′′′

⟨k | V | k ⟩ Mk′′′ k′ +

Π†kk′′

′′

′′′

)

(17)

⟨k | σ · pVσ · p | k ⟩ Πk′′′ k′ ,

k′′ ,k′′′

⟨ ( )−1/2 ′ ⟩ Ωkk′ = k 12 + Y† Y k ,

(18)

⟩ ∂Gkk′ ∂ ⟨ 2 = k p b12 k ′ ∂RA ∂RA { ⟩ ⟨ ( 2 ) ⟩ ∑ ∂Y † ′′ ⟨ p2 ′′′ p † ′′ ′′ ∂ kk + k 12 k Yk′′′ k′ + Ykk′′ k 12 k ′′′ Yk′′′ k′ ∂RA 1 − ep ∂RA 1 − ep k′′ ,k′′′ ⟨ ⟩ } p2 ∂Yk′′′ k′ † 12 k ′′′ + Ykk k ′′ ′′ 1 − ep ∂RA { ⟩ ⟨ † ∑ ∂M ′′ ′′′ ∂V † kk k + ⟨k ′′ | V | k ′′′ ⟩ Mk′′′ k′ + Mkk′′ k ′′ Mk′′′ k′ ∂R ∂R A A k′′ ,k′′′ ∂Mk′′′ k′ ∂Π†kk′′ ′′ + ⟨k | σ · pVσ · p | k ′′′ ⟩ Πk′′′ k′ + ⟨k | V | k ⟩ ∂R ∂R A A ⟩ } ⟨ † † ′′ ∂σ · pVσ · p ′′′ ′′ ′′′ ∂Πk′′′ k′ , + Πkk′′ k Πk′′′ k′ + Πkk′′ ⟨k | σ · pVσ · p | k ⟩ k ∂RA ∂RA † Mkk ′′

′′

′′′

(19) and )−3/2 ′′ ⟩ ∂Ωkk′ 1 ∑ ⟨ ( † k 12 + Y Y =− k ∂RA 2 k′′ ,k′′′

{(

∂Yk†′′ k′′′ ∂RA

) Yk′′′ k′ + Yk†′′ k′′′

(

∂Yk′′′ k′ ∂RA

)} , (20)

where Mkk′ = Kk + αpk Kk bk Ykk′ ,

(21)

Πkk′ = αbk Kk − pk−1 Kk Ykk′ , ( )1/2 ep + 1 K= , 2ep

(22)


(23)


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

Page 12 of 35

b = (ep + 1)−1 ,

(24)

( )1/2 ep = 1 + α2 p2 ,

(25)

α = c−1 .

(26)

and

Here Y is the operator, which is determined by the condition that the Dirac Hamiltonian is completely decoupled by the unitary transformation. { ( ) 1 ′ Ykk′ = α3 pk Kk bk ⟨k | V | k ′ ⟩ Kk′ − p−1 k Kk ⟨k | σ · pVσ · p | k ⟩ bk′ Kk′ ek + ek′ ∑ ( ) − Ykk′′ α2 Kk′′ ⟨k ′′ | V | k ′ ⟩ Kk′ + α4 Kk′′ bk′′ ⟨k ′′ | σ · pVσ · p | k ′ ⟩ bk′ Kk′ k′′

+

∑(

) −1 ′′ 4 ′′ α2 p−1 k Kk ⟨k | σ · pVσ · p | k ⟩ pk′′ Kk′′ + α pk Kk bk ⟨k | V | k ⟩ bk′′ Kk′′ pk′′ Yk′′ k′

k′′

+

∑

(

)

}

′′ ′′′ α3 Ykk′′ Kk′′ bk′′ ⟨k ′′ | σ · pVσ · p | k ′′′ ⟩ Kk′′′ p−1 k′′′ − Kk′′ ⟨k | V | k ⟩ Kk′′′ bk′′′ pk′′′ Yk′′′ k′

k′′ ,k′′′

(27)


,

Page 13 of 35

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 its derivative expression is as follows: [ ∂Ykk′ 1 ∂ek′ ∂ek = −Ykk′ − Ykk′ ∂RA ek + ek ′ ∂RA ∂RA ) ∂ ( ′ + α3 pk Kk bk ⟨k | V | k ′ ⟩ Kk′ − p−1 k Kk ⟨k | σ · pVσ · p | k ⟩ bk′ Kk′ ∂RA } ∑{ ∂ ( ) −1 2 −1 ′′ 4 ′′ + α pk Kk ⟨k | σ · pVσ · p | k ⟩ pk′′ Kk′′ + α pk Kk bk ⟨k | V | k ⟩ bk′′ Kk′′ pk′′ Yk′′ k′ ∂R A ′′ k } { ∑ ) ∂ ( 2 ′′ ′ 4 ′′ ′ α Kk′′ ⟨k | V | k ⟩ Kk′ + α Kk′′ bk′′ ⟨k | σ · pVσ · p | k ⟩ bk′ Kk′ Ykk′′ − ∂R A k′′ } { ∑ ) ( ∂ −1 ′′ ′′′ ′′ ′′′ 3 +α Ykk′′ Kk′′ bk′′ ⟨k | σ · pVσ · p | k ⟩ Kk′′′ pk′′′ − Kk′′ ⟨k | V | k ⟩ Kk′′′ bk′′′ pk′′′ Yk′′′ k′ Ykk′′ ∂RA k′′ ,k′′′ +

∑ {(

−1 ′′′ 4 ′′′ α2 p−1 k Kk ⟨k | σ · pVσ · p | k ⟩ pk′′′ Kk′′′ + α pk Kk bk ⟨k | V | k ⟩ bk′′′ Kk′′′ pk′′′

)

k′′ ,k′′′

( )} ∂Yk′′′ k′ ′′ ′′′ +α3 Ykk′′ Kk′′ bk′′ ⟨k ′′ | σ · pVσ · p | k ′′′ ⟩ Kk′′′ p−1 k′′′ − Kk′′ ⟨k | V | k ⟩ Kk′′′ bk′′′ pk′′′ ∂RA ∑ ( ∂Ykk′′ ) { ( ) + − α2 Kk′′ ⟨k ′′ | V | k ′ ⟩ Kk′ + α4 Kk′′ bk′′ ⟨k ′′ | σ · pVσ · p | k ′ ⟩ bk′ Kk′ ∂RA k′′ ,k′′′ ( ) }] ′′ ′′′ ′′ ⟨k | V | k ⟩ Kk ′′′ bk ′′′ pk ′′′ Yk ′′′ k ′ +α3 Kk′′ bk′′ ⟨k ′′ | σ · pVσ · p | k ′′′ ⟩ Kk′′′ p−1 − K . ′′′ k k (28) Here, σ · pVσ · p operator can be divided into the SF and SD parts using the Dirac identity as follows: σ · pVσ · p = pV · p + iσ · (pV × p) ,

(29)

where i denotes the imaginary unit. The real and imaginary parts of Eq. (29) show the SF and SD components, respectively. The explicit expressions in the space {χ}, are written as ⟨χ | pV · p | χ′ ⟩ = ⟨χ | px V · px + py V · py + pz V · pz | χ′ ⟩


(30)


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

Page 14 of 35

and ⟩ ⟨ 0 χ (pV × p)y χ′   ′ ⟨χ | iσ · (pV × p) | χ ⟩ =  ⟨ ⟩  − χ (pV × p)y χ′ 0   ′ ′ ⟨χ | (pV × p)z | χ ⟩ ⟨χ | (pV × p)x | χ ⟩  + i  ⟨χ | (pV × p)x | χ′ ⟩ − ⟨χ | (pV × p)z | χ′ ⟩ 

with

⟩ ⟨ ⟩] ∂χ ∂χ′ ∂χ ∂χ′ ⟨χ | (pV × p)x | χ ⟩ = V − V , ∂y ∂z ∂z ∂y ⟩ [⟨ ∂χ ∂χ′ ⟩ ⟨ ∂χ ∂χ′ ⟩] ⟨ V V χ (pV × p)y χ′ = − , ∂x ∂z ∂z ∂x

(31)

[⟨

′

and

[⟨ ⟨χ | (pV × p)z | χ′ ⟩ =

⟩ ⟨ ⟩] ∂χ ∂χ′ ∂χ ∂χ′ V − V . ∂x ∂y ∂y ∂x

(32) (33)

(34)

Thus, the SD Hamiltonian demands four types of integrals and a complex formalism, although the SF one needs one type of Eq. (30) and a real formalism. On the other hand, the analytical energy gradient of the LUT-IODKH scheme is remarkably simplified as ⟩ ⟨  ∑  ∂ NR B A   χ V (A = B)  µ C χν  ∂R  A  ⟩  ⟨ C̸=A   ∑ ⟨ ⟩ ∂ ∂ NR B + + NR A A LUT B χµ VA + VB + T + VC χν (A = ̸ B, RAB ≤ τ ) χ h χν =  ∂RA ∂RA µ 2  C̸ = A,B  ⟩ ⟨    ∑  ∂ NR B A NR   χ T + V (A = ̸ B, RAB > τ )  ∂R µ C χν A

C

(35) with [ ( ) ] ∂V+ ∂V+ ∂ † ∂VA A B † † + =Ω M M+Π σ · pVA σ · p Π Ω ∂RA ∂RA ∂RA ∂RA [ ) ] ( ∂ † ∂VB † † +Ω M M+Π σ · pVB σ · p Π Ω. ∂RA ∂RA 14 ACS Paragon Plus Environment

(36)

Page 15 of 35

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


Note that the LUT-IODKH scheme does not require the analytical gradient of the space transformation matrix X because the transformation is determined only in each atom.

2.3

Implementation

This subsection presents the implementation of evaluating the analytical energy gradient of the two-component (LUT-)IODKH method. To clarify the difference in implementation among NR, SF, and SD relativistic methods, we show the evaluation processes for all the methods in detail. Table 1 summarizes the procedures of energy and energy gradient calculations and data types. In the energy calculations, the procedures are divided into three parts: one-electron integral (including TNR , VNR , σ · pVσ · p, matrix transformation, and relativistic transformation), two-electron integral, and SCF (including construction and diagonalization of Fock matrix). In the gradient calculations, the procedures are also divided into three part: one-electron integral, two-electron integral, and the energy derivative defined in Eq. (11). In detail for the NR and SF calculations, all integrals, their derivatives without relativistic integrals, and transformations are evaluated in a real form. For the SF-IODKH calculation, steps for the spin-free parts of the σ · pVσ · p integrals, matrix transformation, relativistic transformation, and these derivatives are added to the procedure in NR. On the other hand, the LUT-IODKH method simplifies the calculation using a relativistic transformation determined in each atom in the energy calculation. For the SD calculations, the spin-dependent parts of the σ ·pVσ ·p integrals in a complex form are added to the SF procedure. Furthermore, the atomic orbital size increases to twice its original size in order to describe the α- and β-spin mixing. Therefore, the relativistic transformation, SCF procedure, and these derivatives require complex calculations with twice the number of basis sets.


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

TNR VNR σ · pVσ · p Matrix transformation X Relativistic transformation 1/rij Fock Diagonalization

TNR VNR σ · pVσ · p Matrix transformation X Relativistic transformation Derivative of two-electron integral 1/rij Summation for total gradient

Energy gradient Derivative of one-electron integral

Two-electron integral SCF

Energy One-electron integral

real real – – – real real

real real – – – real real real

NR

real real real real real real real

real real real real real real real real

IODKH

real real real – real real real

real real real real (atom) real real real real

LUT-IODKH

SF

real real complex real complex real complex

real real complex real complex real complex complex

IODKH

real real complex – complex real complex

real real complex real (atom) complex real complex complex

LUT-IODKH

SD

Table 1: Procedures of energy and energy gradient calculations and data types in NR, (LUT-)SF-IODKH, and (LUT-)SDIODKH.


ACS Paragon Plus Environment

16

Page 16 of 35

Page 17 of 35

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


3

Numerical assessments

3.1

Computational details

This section examines the performance of the analytical energy gradients for the 2c IODKH and LUT-IODKH methods in terms of accuracy and efficiency. As the 2c relativistic Hamiltonian, the one-electron SD-IODKH Hamiltonian and two-electron Coulomb interaction, denoted by SD-IODKH, were adopted. In the LUT scheme, the threshold for cutoff of relativistic interaction, τ , was set to 3.5 ˚ A. For comparison, the 4c Dirac-Coulomb, (LUT-) SF-IODKH, and NR methods were also employed. Twelve closed-shell diatomic molecules containing fifth- and sixth-period atoms were calculated as numerical tests: InH, SnO, Sb2 , HI, I2 , AuH, Au2 , TlH, PbO, Bi2 , HAt, and At2 . Furthermore, multi-dimensional silver clusters were employed to investigate the efficiency in LUT-SD-IODKH. The Ag-Ag bond length in one-dimensional (1D) Agn clusters (n = 2, · · · , 10) was fixed at 2.706 ˚ A. In two- and three-dimensional (2D/3D) Agn clusters, an experiment lattice constant of 4.806 ˚ A with a face-centered cubic structure was adopted. 61 The total numbers of atoms in the 2D Agn clusters were 5, 8, 11, 13, 14, 17, 18, 23, 25, 28, 32, 39, 41, and 50, which correspond to 1 × 1, 1 × 2, 1 × 3, 2 × 2, 1 × 4, 1 × 5, 2 × 3, 2 × 4, 3 × 3, 2 × 5, 3 × 4, 3 × 5, 4 × 4, and 4 × 5 unit cells, respectively. The total numbers in the 3D Agn clusters were 14, 23, 32, 41, and 53, which correspond to 1×1×1, 1×1×2, 1×3×1, 1×1×4, and 1 × 2 × 3, respectively. For a practical calculation, fac–tris(2–phenylpyridine)iridium, fac−Ir(ppy)3 , was also examined. All the NR, SF, and SD calculations were performed at the RHF, RHF, and GHF levels, respectively. For the diatomic molecules and Agn clusters, the Sapporo-DZP-2012 62 and DKH3-Gen-TK/NOSeC-V-TZP 63–65 basis sets with uncontracted forms were employed in the H atom and other atoms, respectively. In fac−Ir(ppy)3 , the contracted DKH3-GenTK/NOSeC-V-TZP for Ir and cc-pVDZ 66 for H, C, and N atoms were adopted. Note that in the NR calculation for fac−Ir(ppy)3 , the basis set for Ir, TZP, 67 was adopted. The ana-



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

lytical energy gradient for the (LUT-)SD-IODKH method was implemented in the GAMESS program. 68 The DIRAC12 program was utilized for the 4c calculations. 69

N Ir N

N

Figure 1: Structure of fac–tris(2–phenylpyridine)iridium.

3.2

Accuracy of SD-IODKH and LUT-SD-IODKH methods

This section examines the accuracy of the analytical energy gradient method for the 2c IODKH method with/without the LUT scheme. Table 2 shows the analytical and numerical gradient values at equilibrium lengths of re obtained by NR, SD-IODKH, LUT-SD-IODKH, and 4c at the HF level. Here, the experimental bond lengths were used for re , namely, re = 1.742 ˚ A for HAt, 31 3.046 ˚ A for At2 , 70 1.5324 ˚ A for AuH, 71 and 2.4719 ˚ A for Au2 . 72 Furthermore, 1.5re and 2.0re for AuH and HAt were also calculated. The numerical gradient values were obtained using a five point numerical differentiation method at the 0.001 ˚ A step size. The inherent accuracy of this numerical differentiation is a fourth-order accuracy. The convergence threshold of the SCF calculation is 5.0 × 10−6 for the density changes. A quadruple precision routines for evaluating one-electron integrals including overlap and kinetic energy, and for diagonalizing the integral matrices were used to keep accuracy in the Hess’ numerical technique. 60 The deviation of the analytical gradient with respect to the numerical gradient, ∆g ana , shows the reliability of the present implementation. The maximum absolute deviations of NR, SD-IODKH, LUT-SD-IODKH, and 4c are 0.000002, 0.000013, 0.000008, and 0.000002, respectively. The deviations of SD-IODKH and LUT-SD-IODKH are slightly larger than that of NR. This is due to the resolution18 ACS Paragon Plus Environment

Page 18 of 35

Page 19 of 35

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-identity (RI) approximation, which is essential for Hess’ numerical calculation 60 in the IODKH method. Considering these theoretical demands, the results indicate the correct implementation of the present analytical energy gradients in SD-IODKH and LUT-SD-IODKH within the Hess approximation.


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

1.0re 1.5re 2.0re 1.0re 1.5re 2.0re 1.0re 1.0re

AuH

Au2 At2

HAt

Length

Molecules

∆g

0.082251 0.000001 0.033559 0.000000 0.031298 0.000000 0.011267 0.000002 0.067515 0.000001 0.037875 0.000000 0.079816 0.000002 0.025557 -0.000001

g

NR

0.012531 0.054765 0.036769 0.008891 0.054690 0.000126 0.028908 0.007775

g 0.000013 0.000007 0.000000 0.000009 −0.000005 0.000000 0.000006 −0.000002

∆g

SD-IODKH

0.012518 0.054768 0.036770 0.008908 0.054686 0.000126 0.028880 0.007796

g 0.000007 0.000001 0.000002 0.000001 0.000000 0.000008 0.000004 0.000000

∆g

LUT-SD-IODKH ∆g

0.011711 0.000000 0.054581 0.000000 0.036850 0.000000 0.009392 0.000000 0.055145 0.000001 0.031849 0.000000 0.029315 0.000000 0.008351 -0.000002

g

4c

Table 2: Analytical and/or numerical gradient values (hartree/bohr) of HAt, At2 , AuH, and Au2 calculated by NR, SD-IODKH, LUT-SD-IODKH, and 4c at the HF level.



20

Page 20 of 35

Page 21 of 35

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


Table 3 presents the bond lengths (˚ A) of diatomic molecules containing fifth- and sixthperiod elements calculated by the NR, SF-IODK/C, SD-IODKH, LUT-SD-IODKH, and 4c methods. The deviations of the NR, SF-IODKH, and SD-IODKH methods from 4c, ∆r4c , indicates the fully relativistic effect, the SD relativistic effect in one-electron terms, and the fully relativistic effect in two-electron terms, respectively. The deviation of SD-IODKH from SF-IODKH, ∆rSD , is shown in square brackets, representing the SD relativistic effect in one-electron terms. The deviation of LUT-SD-IODKH from SD-IODKH, ∆rLUT , indicates the accuracy of the LUT scheme. The fully relativistic effect on bond lengths, ∆r4c of NR, is mostly positive whereas I2 , Bi2 , and At2 are opposite. The absolute values are large in AuH and Au2 : 0.263 ˚ A for AuH and 0.330 ˚ A for Au2 . The large deviations are reduced by the SF relativistic effect. Actually, ∆r4c of SF-IODKH are 0.007 ˚ A for AuH and 0.008 ˚ A for Au2 . This indicates that AuH and Au2 have an s-bonding nature formed by each s valence orbital, which is largely contracted by the relativistic effect. However, the absolute values |∆r4c | of SF-IODKH exhibit large differences, especially in Bi2 and At2 : −0.095 ˚ A for Bi2 and −0.129 ˚ A for At2 . On the other hand, |∆r4c | of SD-IODKH are small in all molecules. This indicates that the SD relativistic effect in heavy-element systems is important when the bond has non-zero angular momentum such as p, d, and f orbitals. The contributions of two-electron relativistic effects can be small in bond lengths of the systems used here, although the contributions sometimes can be important. 55 Furthermore, TlH has a negative value by the SD relativistic effect, even though Bi2 , HAt, and At2 have positive values. The SD relativistic effect gives two different contributions to bond length: a splitting of degenerate orbitals and an overlap of orbitals by a second-order spin-orbit interaction. The former is when p, d, and f orbitals are split into {p1/2 , p3/2 }, {d3/2 , d5/2 }, and {f5/2 , f7/2 }, respectively. The orbitals with less (more) total angular momentum are stabilized and contracted (destabilized and expanded). The latter contributes to the stabilization of a molecule by an overlap of orbitals with the same total angular momentum and spatial symmetry, i.e., σ1/2g and π1/2g or π1/2u and σ1/2u ,



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

Page 22 of 35

which are split by a spin-orbit interaction. In the case of TlH, the former contribution is dominant because two electrons occupy the σ1/2g orbital, which is stabilized by the SD relativistic effect. In the case of Bi2 , HAt, and At2 , the latter contribution is dominant because electrons occupy bonding and/or anti-bonding orbitals. In this situation, the SD relativistic effect stabilizes a molecule largely when the bond length is long, although the opposite is true for the SF relativistic effect. As a result, the bonds are elongated by the SD relativistic effect. In all cases, ∆rLUT is considerably small. This represents that the LUT treatment can sufficiently describe both SF and SD relativistic effects for estimating the bond lengths. Note that in systems whose bond lengths are longer such as van der Waals complexes, we carefully set the cutoff radius including nearest neighbor atoms in order to keep accuracy. Details are given in Supporting Information. Table 3: Bond lengths (˚ A) of diatomic molecules containing fifth- and sixth-period elements calculated by the NR, SF-IODKH, SD-IODKH, LUT-SD-IODKH, and 4c methods at the HF level.

Molecule

InH SnO Sb2 HI I2 AuH Au2 TlH PbO Bi2 HAt At2

IODKH

NR SF

4c

SD

LUT-SD

rNR

∆r4c

rSF

∆r4c

rSD

∆r4c

∆rSD

rLUT

∆rLUT

r4c

1.856 1.798 2.460 1.608 2.675 1.828 2.924 1.932 1.890 2.647 1.711 2.890

0.012 0.008 0.018 0.005 −0.003 0.263 0.330 0.069 0.019 −0.015 0.001 −0.081

1.845 1.789 2.436 1.600 2.664 1.572 2.602 1.890 1.864 2.567 1.683 2.842

0.001 −0.001 −0.006 −0.003 −0.014 0.007 0.008 0.027 −0.007 −0.095 −0.027 −0.129

1.843 1.790 2.443 1.603 2.679 1.569 2.593 1.859 1.873 2.667 1.712 2.976

−0.001 0.000 0.001 0.000 0.001 0.004 −0.001 −0.004 0.002 0.005 0.002 0.005

−0.002 0.001 0.007 0.003 0.015 −0.003 −0.009 −0.031 0.009 0.100 0.029 0.134

1.843 1.790 2.443 1.603 2.679 1.568 2.593 1.859 1.873 2.667 1.712 2.978

0.000 0.000 0.000 0.000 0.000 −0.001 0.000 0.000 0.000 0.000 0.000 0.002

1.844 1.790 2.442 1.603 2.678 1.565 2.594 1.863 1.871 2.662 1.710 2.971


Page 23 of 35

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


3.3

Computational cost of (LUT-)SD-IODKH method

This subsection examines the computational cost of the analytical energy gradient for the SDIODKH method with/without the LUT scheme on 1D, 2D, and 3D silver clusters. Figure 2 shows the system-size dependence of the central processing unit (CPU) time for the analytical gradient of the one-electron integrals in 1D, 2D, and 3D silver clusters obtained by the SDIODKH and LUT-SD-IODKH methods. The horizontal axis indicates the total numbers of the atoms n. The vertical axis indicates the CPU time in seconds. A Xeon E5-2690/2.90 GHz (Octa core) processor on a single core was utilized to measure the CPU time. In 1D cluster, the CPU time in SD-IODKH is long; the scaling is O(n2.98 ), which corresponds to an order in a matrix multiplication. On the other hand, the LUT scheme can drastically reduce the CPU time; the scaling is quasi-linear, O(n1.33 ). The results also represent that the LUT-SD-IODKH method achieves quasi-linear scaling for 2D and 3D clusters: O(n1.35 ) for 2D and O(n1.36 ) for 3D. Table 4 shows the CPU time for an analytical energy gradient calculation of Ag10 . The NR, SF-IODKH, LUT-SF-IODKH, SD-IODKH, and LUT-SD-IODKH methods were adopted. The calculation is divided into six steps: namely, one-electron integrals including the 2c transformation (OEI), an initial guess (Guess) by H¨ uckel scheme, two-electron integrals (TEI), the SCF procedure, derivatives of OEI including the 2c transformation (dOEI), and derivatives of TEI (dTEI). The percentage of CPU time in each step is given in parentheses. In single-point energy calculations including OEI, Guess, TEI, and SCF, SF- and SDIODKH for OEI require long CPU times in comparison with NR: 0.12 s for NR, 75.72 s for SF-IODKH, and 1487.83 s for SD-IODKH. The long CPU times for OEI are reduced by the LUT schemes: 4.47 s for LUT-SF-IODKH and 24.69 s for LUT-SD-IODKH. In the SCF step, the CPU time at the SD level is ≥10 times as much as that at the SF level. This is because the GHF method for the SD effect requires a calculation with a complex form and twice as large basis sets as the RHF one. 23 ACS Paragon Plus Environment


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

In gradient calculations, including dOEI and dTEI, SF- and SD-IODKH for dOEI demand large CPU time in comparison to NR: 1.11 s for NR, 2652.28 s for SF-IODKH, and 61247.55 s for SD-IODKH. The LUT scheme reduces the CPU times: 129.36 s for LUT-SF-IODKH and 263.35 s for LUT-SD-IODKH. Furthermore, in dTEI, the CPU times at the SD level are longer than those at the SF one. This is because additional CPU times are required to multiply the density matrix with the GHF formalism to the derivative of TEI, although the integral calculations are the same. These results indicate that the LUT scheme resolves one of the bottlenecks, i.e., long CPU times in the OEI and dOEI steps.


Page 24 of 35

2

4 n

6

)

0 8

10

0

10

0

1.33

1000

30

40

20

O(n

O(n2.98)

50

60

2000

3000

4000

5000

IODKH LUT-IODKH

70

0

(b)

10

20

30 n

2D-cluster 3D-cluster

40

50

O(n1.35)

O(n1.36)

60

Figure 2: System-size dependence of CPU time (s) of the one-electron analytical energy gradient obtained with the SD-IODKH and LUT-SD-IODKH methods: (a) 1D Agn cluster and (b) 2D and 3D Agn cluster. A Xeon E5-2690/2.90 GHz (Octa Core) processor was used on a single core.

CPU time (s)

6000

(a)

CPU time (s)

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

7000

Page 25 of 35 Journal of Chemical Theory and Computation


25

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

OEI Guess TEI SCF dOEI dTEI Total

%

0.12 (0.00) 20.55 (0.75) 675.33 (24.66) 817.20 (29.84) 1.11 (0.04) 1224.57 (44.71) 2738.88 (100.00)

Time (s)

NR %

75.72 (1.35) 20.63 (0.37) 725.13 (12.91) 904.83 (16.11) 2652.28 (47.21) 1238.95 (22.06) 5617.54 (100.00)

Time (s)

SF-IODKH %

4.47 (0.15) 19.52 (0.67) 687.98 (23.65) 814.73 (28.01) 129.36 (4.45) 1253.02 (43.07) 2909.08 (100.00)

Time (s)

LUT-SF-IODKH

1487.83 20.45 690.26 10507.45 61427.55 2046.20 76179.74

Time (s) (1.95) (0.03) (0.91) (13.79) (80.64) (2.69) (100.00)

%

SD-IODKH

%

24.69 (0.18) 20.39 (0.15) 693.53 (5.04) 10687.16 (77.74) 263.35 (1.92) 2058.81 (14.98) 13747.93 (100.00)

Time (s)

LUT-SD-IODKH

Table 4: CPU time for six steps in the restricted or general Hartree-Fock calculation obtained with NR, SF-IODKH, LUT-SFIODKH, SD-IODKH, and LUT-SD-IODKH in Ag10 linear cluster: one-electron integral (OEI), initial guess (Guess) calculated with the H¨ uckel scheme, two-electron integral (TEI), SCF procedure (SCF), derivative of OEI (dOEI), and derivative of TEI (dTEI). A Xeon E5-2690/2.90 GHz processor (Octa core) was used on a single core.



26

Page 26 of 35

Page 27 of 35

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


3.4

Application in fac-Ir(ppy)3

This subsection discusses the performance of the LUT-SF-IODKH and LUT-SD-IODKH methods in fac-Ir(ppy)3 . It should be noted that this molecule cannot be treated without the LUT technique for the relativistic calculation. Table 5 shows the geometric parameters calculated at the Hartree-Fock level of theory. A-B and A-B-C denote bond lengths between atoms A and B and angles among atoms A, B, and C, respectively. ∆SF is the difference of the geometrical parameters between NR and LUT-SF-IODKH, indicating the SF relativistic effect on the geometry. ∆SD is the difference of the geometrical parameters between LUT-SFIODKH and LUT-SD-IODKH, indicating the SD relativistic effect. ∆Rel is the summation of ∆SF and ∆SD , indicating that total relativistic effect. Note that we used the nearest neighboring atom from the central metal Ir to estimate the bond lengths and angles. In the bond lengths of fac-Ir(ppy)3 , the SF relativistic effect is smaller than the SD relativistic effect. For Ir-N and Ir-C, {∆SF , ∆SD } are {−0.027, −0.049} and {0.020, −0.087} in ˚ A, respectively. On the other hand, the three bond angles, C-Ir-N, C-Ir-C, and N-Ir-N, are slightly affected by both SF and SD effects. For C-Ir-N, C-Ir-C, and N-Ir-N, {∆SF , ∆SD } are {3.2, −2.2}, {1.6, −4.4}, and {−1.2, 0.0} in degrees, respectively. Table 5: Geometric parameters of Ir(ppy)3 obtained by NR, LUT-SF-IODKH and LUT-SDIODKH at the Hartree-Fock level. Parameters Ir-N Ir-C C-Ir-N C-Ir-C N-Ir-N

NR

LUT-SF-IODKH

LUT-SD-IODKH

∆SF

∆SD

∆Rel

2.225 2.031 172.2 95.3 97.7

2.198 2.051 175.4 96.9 96.5

2.149 1.964 173.2 92.5 96.5

−0.027 0.020 3.2 1.6 −1.2

−0.049 −0.087 −2.2 −4.4 0.0

−0.076 −0.067 1.0 −2.8 −1.2



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

Conclusion

The present study develops the analytical energy gradient for the GHF method based on the 2c IODKH method with/without the efficient LUT scheme in a 2c relativistic transformation. The accuracy of the present method was assessed using several diatomic molecules containing fifth- and sixth-period elements. The optimized bond lengths are in good agreement with those of the 4c method. Furthermore, the deviations of the LUT scheme are less than 0.003 ˚ A in bond length. The efficiency of the LUT scheme was also examined using 1D, 2D, and 3D silver clusters. The CPU time for a one-electron relativistic transformation scales quasilinear in these clusters. A practical calculation using a fac-Ir(ppy)3 molecule was performed to estimate the accuracy and efficiency of the proposed method. The results show that the SD relativistic effect shortens the bond length by approximately 0.1 ˚ A. In addition, although the SD relativistic program is much more complicated than that of NR, the LUT-SD-IODKH scheme demands only about twice as much CPU time as NR. The present analytical energy gradient method for LUT-IODKH with the SD relativistic effect was found to be applicable to practical systems.

Acknowledgement Some of the present calculations were performed at the Research Center for Computational Science (RCCS), Okazaki Research Facilities, National Institutes of Natural Sciences (NINS). This study was supported in part by the Strategic Programs for Innovative Research (SPIRE), Ministry of Education Culture, Sports , Science and Technology (MEXT), and the Computational Materials Science Initiative (CMSI), Japan; by the MEXT program “Elements Strategy Initiative to Form Core Research Center” (since 2012), MEXT, Japan; and by Core Research for Evolutional Science and Technology (CREST) Program “Theoretical Design of Materials with Innovative Functions Based on Relativistic Electronic Theory” of Japan Science and Technology Agency (JST). One of the authors (J.S.) is grateful to 28 ACS Paragon Plus Environment

Page 28 of 35

Page 29 of 35

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 Japan Society for the Promotion of Science (JSPS) for a Research Fellowship. One of the authors (Y.N.) is grateful for “Early Bird” grant for young researchers at the Research Institute for Science and Engineering (RISE) of Waseda University.

Supporting Information Available Bond lengths of van der Waals complexes obtained by (LUT-)IODKH. This material is available free of charge via the Internet at http://pubs.acs.org/.

References (1) Moss, R. E. Advanced Molecular Quantum Mechanics; Chapman and Hall, London, 1973. (2) Breit, G. Phys. Rev. 1929, 34, 553–573. (3) Gaunt, J. A. Proc. R. Soc. Lond. A 1929, 124, 163–176. (4) Pyykkö, P. Chem. Rev. 1988, 88, 563–594. (5) Salek, P.; Helgaker, T.; Saue, T. Chem. Phys. 2005, 311, 187–201. (6) Komorovsk´ y, S.; Repsisk´ y, M.; Malkina, O. L.; Malkin, V. G.; Malkin, I. O.; Kaupp, M. J. Chem. Phys. 2008, 128, 104101. (7) Nakajima, T.; Yanai, M.; Hirao, K. J. Comput. Chem. 2002, 23, 847–860. (8) Kelly, M. S.; Shiozaki, T. J. Chem. Phys. 2013, 138, 204113. (9) Shiozaki, T. J. Chem. Theory Comput. 2013, 9, 4300–4303. (10) Dyall, K. G. J. Chem. Phys. 1997, 106, 9618–9626. (11) Liu, W.; Peng, D. J. Chem. Phys. 2006, 125, 044102. 29 ACS Paragon Plus Environment


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

(12) Barysz, M.; Sadlej, A. J. J. Chem. Phys. 2002, 116, 2696–2704. (13) Samzow, R.; Hess, B. A.; Jansen, G. J. Chem. Phys. 1992, 96, 1227–1231. (14) Nakajima, T.; Hirao, K. J. Chem. Phys. 2003, 119, 4105. (15) van W¨ ullen, C.; Michauk, C. J. Chem. Phys. 2005, 123, 204113. (16) Seino, J.; Hada, M. Chem. Phys. Lett. 2007, 442, 134–139. (17) Seino, J.; Hada, M. Chem. Phys. Lett. 2008, 461, 327–331. (18) Sikkema, J.; Visscher, L.; Saue, T.; lliaˇs, M. J. Chem. Phys. 2009, 131, 124116. (19) Wahlgren, U.; Schimmelpfennig, B.; Jusuf, S.; Strömsnes, H.; Gropen, O.; Maron, L. Chem. Phys. Lett. 1998, 287, 525–530. (20) Thar, J.; Kirchner, B. J. Chem. Phys. 2009, 130, 124103. (21) Seino, J.; Nakai, H. J. Chem. Phys. 2012, 136, 244102. (22) Seino, J.; Nakai, H. J. Chem. Phys. 2012, 137, 144101. (23) Seino, J.; Nakai, H. J. Chem. Phys. 2013, 139, 034109. (24) Seino, J.; Nakai, H. Int. J. Quantum Chem. 2015, 115, 253–257. (25) Seino, J.; Nakai, H. J. Comput. Chem. Jpn. 2014, 13, 1–17. (26) Peng, D.; Reiher, M. J. Chem. Phys. 2012, 136, 244108. (27) Dyall, K.; Fægri, K. Introduction to Relativistic Quantum Chemistry; Oxford University Press, 2007. (28) Reiher, M.; Wolf, A. Relativistic Quantum Chemistry; Wiley-VCH, Germany, 2009. (29) Autschbach, J. J. Chem. Phys. 2012, 136, 150902.


Page 30 of 35

Page 31 of 35

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


(30) Wang, Z.; Wang, F. Phys. Chem. Chem. Phys. 2013, 15, 17922–17928. (31) Han, Y.-K.; Bae, C.; Son, S.-K.; Lee, Y. S. J. Chem. Phys. 2000, 112, 2684–2691. (32) Vaara, J.; Ruud, K.; Vahtras, O.; ˚ Agren, H.; Jokisaari, J. J. Chem. Phys. 1998, 109, 1212–1222. (33) Malkina, O.; Schimmelpfennig, B.; Kaupp, M.; Hess, B.; Chandra, P.; Wahlgren, U.; Malkin, V. Chem. Phys. Lett. 1998, 296, 93–104. (34) Launila, O.; Schimmelpfennig, B.; Fagerli, H.; Gropen, O.; Taklif, A. G.; Wahlgren, U. J. Mol. Spectrosc. 1997, 186, 131–143. (35) Fagerli, H.; Schimmelpfennig, B.; Gropen, O.; Wahlgren, U. J. Mol. Struct. (THEOCHEM) 1998, 451, 227 – 235. (36) Rakowitz, F.; Marian, C. M.; Schimmelpfennig, B. Phys. Chem. Chem. Phys. 2000, 2, 2481–2488. (37) Malmqvist, P. ˚ A.; Roos, B. O.; Schimmelpfennig, B. Chem. Phys. Lett. 2002, 357, 230–240. (38) Filatov, M.; Cremer, D. Chem. Phys. Lett. 2003, 370, 647–653. (39) Zou, W.; Filatov, M.; Cremer, D. J. Chem. Phys. 2011, 134, 244117. (40) Filatov, M.; Zou, W.; Cremer, D. Int. J. Quantum Chem. 2014, 114, 993–1005. (41) Cheng, L.; Gauss, J. J. Chem. Phys. 2011, 135, 084114. (42) Cheng, L.; Gauss, J.; Stanton, J. F. J. Chem. Phys. 2013, 139, 054105. (43) Nakajima, Y.; Seino, J.; Nakai, H. J. Chem. Phys. 2013, 139, 244107. (44) van Lenthe, E.; Ehlers, A.; Baerends, E.-J. J. Chem. Phys. 1999, 110, 8943–8953.



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

(45) van Lenthe, J. H.; Faas, S.; Snijders, J. G. Chem. Phys. Lett. 2000, 328, 107–112. (46) Filatov, M.; Cremer, D. J. Chem. Phys. 2003, 118, 6741–6750. (47) Fedrov, D. G.; Nakajima, T.; Hirao, K. Chem. Phys. Lett. 2001, 335, 183–187. (48) Nasluzov, V. A.; Rösch, N. Chem. Phys 1996, 210, 413–425. (49) de Jong, W. A.; Harrison, R. J.; Dixon, D. A. J. Chem. Phys. 2001, 114, 48–53. (50) Malkina, I.; Malkina, O. L.; Malkin, V. G. Chem. Phys. Lett. 2002, 361, 231–236. (51) Neese, F.; Wolf, A.; Fleig, T.; Reiher, M.; Hess, B. A. J. Chem. Phys. 2005, 122, 204107. (52) Matveev, A. V.; Nasluzov, V. A.; Rösch, N. Int. J. Quantum Chem. 2007, 107, 3236– 3249. (53) Mastalerz, R.; Barone, G.; Lindh, R.; Reiher, M. J. Chem. Phys. 2007, 127, 074105. (54) Autschbach, J.; Peng, D.; Reiher, M. J. Chem. Theory Comput. 2012, 8, 4239–4248. (55) Zou, W.; Filatov, M.; Cremer, D. J. Chem. Phys. 2015, 142, 214106. (56) Cheng, L.; Stopkowicz, S.; Gauss, J. Int. J. Quantum Chem. 2014, 114, 1108–1127. (57) Cremer, D.; Zou, W.; Filatov, M. WIREs Comput. Mol. Sci. 2014, 4, 436–467. (58) Barysz, M.; Sadlej, A. J.; Snijders, J. G. Int. J. Quantum Chem. 1997, 65, 225–239. (59) Barysz, M. J. Chem. Phys. 2001, 114, 9315–9324. (60) Hess, B. A. Phys. Rev. A 1985, 32, 756–763. (61) Suh, I.-K.; Ohta, H.; Waseda, Y. J. Mater. Sci. 1988, 23, 757–760. (62) Noro, T.; Sekiya, M.; Koga, T. Theor. Chem. Acc. 2003, 109, 85–90.


Page 32 of 35

Page 33 of 35

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


(63) Sekiya, M.; Noro, T.; Osanai, Y.; Koga, T. Theor. Chem. Acc. 2001, 106, 297–300. (64) Osanai, Y.; Noro, T.; Miyoshi, E.; Sekiya, M.; Koga, T. J. Chem. Phys. 2004, 120, 6408. (65) Tatewaki, H.; Koga, T. Chem. Phys. Lett. 2000, 328, 473–482. (66) Dunning Jr., T. H. J. Chem. Phys. 1989, 90, 1007–1023. (67) Martins, L. S. C.; Jorge, F. E.; Machado, S. F. Mol. Phys. 2015, (68) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery Jr., J. A. J. Comput. Chem. 1993, 14, 1347–1363. (69) DIRAC, a relativistic ab initio electronic structure program, Release DIRAC12 (2012), written by H. J. Aa. Jensen, R. Bast, T. Saue, and L. Visscher, with contributions from V. Bakken, K. G. Dyall, S. Dubillard, U. Ekström, E. Eliav, T. Enevoldsen, T. Fleig, O. Fossgaard, A. S. P. Gomes, T. Helgaker, J. K. Lærdahl, Y. S. Lee, J. Henriksson, M. Iliaˇs, Ch. R. Jacob, S. Knecht, S. Komorovsk´ y, O. Kullie, C. V. Larsen, H. S. Nataraj, P. Norman, G. Olejniczak, J. Olsen, Y. C. Park, J. K. Pedersen, M. Pernpointner, K. Ruud, P. Salek, B. Schimmelpfennig, J. Sikkema, A. J. Thorvaldsen, J. Thyssen, J. van Stralen, S. Villaume, O. Visser, T. Winther, and S. Yamamoto (see http://www.diracprogram.org). (70) Visscher, L.; Dyall, K. G. J. Chem. Phys. 1996, 104, 9040–9046. (71) Seto, J. Y.; Morbi, Z.; Charron, F.; Lee, S. K.; Bernath, P. F.; Roy, R. J. L. J. Chem. Phys. 1999, 110, 11756–11767. (72) Huber, K. P.; Herzberg, G. Chemistry Spectra and Molecular Structure IV. Constants of Diatomic Molecules; Van Nostrand, New York, 1979.



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

Graphical TOC Entry 0.15 0.15 0.10 0.10

Diff. from 4c bond length

0.05 0 .05 in

SF-IODKH SF S

AuH A Au u2 TlH H PbO bO OB Bii2 HAt At A Att2

SD-IODKH KH

LUT-SD-IODKH T SD-IIODKH


Page 34 of 35

Page 35 of 35

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


Graphical abstract 90x35mm (300 x 300 DPI)


Implementation of Analytical Energy Gradient of Spin-Dependent

Recommend Documents