Analytical energy gradient based on spin-free infinite-order Douglas-Kroll-Hess method with local unitary transformation Yuya Nakajima, Junji Seino, and Hiromi Nakai Citation: The Journal of Chemical Physics 139, 244107 (2013); doi: 10.1063/1.4850638 View online: http://dx.doi.org/10.1063/1.4850638 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/139/24?ver=pdfcov Published by the AIP Publishing

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: 134.99.128.41 On: Tue, 31 Dec 2013 15:59:03

THE JOURNAL OF CHEMICAL PHYSICS 139, 244107 (2013)

Analytical energy gradient based on spin-free infinite-order Douglas-Kroll-Hess method with local unitary transformation Yuya Nakajima,1 Junji Seino,1 and Hiromi Nakai1,2,3,4,a) 1

Department of Chemistry and Biochemistry, School of Advanced Science and Engineering, Waseda University, Tokyo 169-8555, Japan 2 Research Institute for Science and Engineering, Waseda University, Tokyo 169-8555, Japan 3 CREST, Japan Science and Technology Agency, 4-1-8 Honcho, Kawaguchi, Saitama 332-0012, Japan 4 Elements Strategy Initiative for Catalysts and Batteries (ESICB), Kyoto University, Katsura, Kyoto 615-8520, Japan

(Received 17 September 2013; accepted 4 December 2013; published online 27 December 2013) In this study, the analytical energy gradient for the spin-free infinite-order Douglas-Kroll-Hess (IODKH) method at the levels of the Hartree-Fock (HF), density functional theory (DFT), and second-order Møller-Plesset perturbation theory (MP2) is developed. Furthermore, adopting the local unitary transformation (LUT) scheme for the IODKH method improves the efficiency in computation of the analytical energy gradient. Numerical assessments of the present gradient method are performed at the HF, DFT, and MP2 levels for the IODKH with and without the LUT scheme. The accuracies are examined for diatomic molecules such as hydrogen halides, halogen dimers, coinage metal (Cu, Ag, and Au) halides, and coinage metal dimers, and 20 metal complexes, including the fourth–sixth row transition metals. In addition, the efficiencies are investigated for one-, two-, and three-dimensional silver clusters. The numerical results confirm the accuracy and efficiency of the present method. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4850638] I. INTRODUCTION

Reliable predictions of geometries and properties for molecules including heavy elements are some of the most important tasks in quantum chemical calculations. Notably, in case of calculations for heavy (and super heavy) elements and their compounds, the inclusion of relativistic effect is essential because electrons, especially core electrons, have considerably high speeds. For example, bond contractions/elongations due to relativistic shrinking/expansion of molecular orbitals1–4 were observed in several diatomic molecules such as halogen dimers and PtAu.5, 6 The evaluations of first- and higher-order derivatives of the molecular energy enables computations of optimized geometries, various molecular properties such as molecular vibrations observed in infrared and Raman spectra and chemical shifts observed in nuclear magnetic resonance (NMR) spectra. Accordingly, the development of analytical energy derivatives with the relativistic Hamiltonian is desired in order to extend the applicability of the quantum chemical method. A highly accurate treatment in relativistic quantum chemistry is based on the four-component (4c) relativistic theory, which adopts a one-particle Dirac operator that satisfies Lorentz invariance only in terms of the motion of each electron. Numerical applications of the analytical energy gradient method with the 4c treatments confirmed the high accuracy in determining molecular geometries,1 and electric and magnetic properties.7, 8 Such applications, however, were suitable only in the case small- or medium-sized molecules because of the practical problems involving computational costs a) Electronic mail: [email protected]

0021-9606/2013/139(24)/244107/13/$30.00

and/or variational conditions9 due to the existence of smallcomponent Dirac spinors or positronic-state spinors. There have been several attempts to resolve these problems.10–12 An alternative approach is based on the two-component (2c) relativistic theory, which treats only electronic-state (or large-component) information. Recent developments in the 2c treatment of many-electron molecular systems have exhibited a higher degree of accuracy, in comparison with 4c treatment.13–18 The energy gradient methods for accurate 2c relativistic schemes are under development. The 2c relativistic theories are categorized into two types. The first approach uses an equation for the elimination of small components (ESC), such as the regular approximation (RA),19–21 the relativistic ESC (RESC),22 and the normalized elimination of the small component (NESC).23 The analytical energy gradient schemes have been made available for the zeroth-order RA and its scaling schemes,24–26 RESC,27 and NESC.28, 29 The second approach uses the decoupling matrix transformation of 4c Dirac Hamiltonian and includes approaches such as the exact quasi-relativistic method.30 In this type of method we separate the electronic and positronic Hamiltonian from 4c Dirac Hamiltonian or Fock matrix using unitary transformation. The widely used schemes are the FoldyWouthuysen (FW) transformation31 and Douglas-Kroll-Hess (DKH) transformation.32 These were extended to higherorder DKH33–36 and infinite-order DKH (IODKH)37 methods. The molecular property calculation for this type of scheme is based on the second-order38–40 and extended35, 41, 42 DKH Hamiltonian. The authors’ group has worked on the development of an accurate and efficient 2c relativistic scheme, termed local unitary transformation (LUT), at the one- and two-particle

139, 244107-1

© 2013 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: 134.99.128.41 On: Tue, 31 Dec 2013 15:59:03

244107-2

Nakajima, Seino, and Nakai

J. Chem. Phys. 139, 244107 (2013)

spin-free IODKH levels.43, 44 The LUT method utilizes the locality of the relativistic effect. Accordingly, the computational costs become drastically lower than those of conventional methods and scale linearly with respect to the system size with a small prefactor. In addition, the LUT scheme clarifies the primary relativistic effect in molecular systems, because the deviations in total energy are less than millihartrees. Recently, the LUT scheme has been combined with the linear-scaling divide-and-conquer (DC) based on HartreeFock (HF)45–47 and electron correlation48–53 methods, such as the second-order Møller-Plesset (MP2) and the coupled cluster theories with single and double excitations (CCSD).54 The combination of LUT and DC techniques leads to the first approach that achieves overall (quasi-)linear scaling with a small prefactor for relativistic electron correlation calculations. In this study, we derive the analytical energy gradient of the one-particle IODKH Hamiltonian for the accurate calculations. Furthermore, the LUT scheme is introduced for efficient calculations for large-sized molecules. This article is organized as follows. Theoretical aspects are explained in Sec. II. In Sec. III, the present scheme is numerically assessed from the viewpoints in accuracy and efficiency. Finally, the concluding remarks are provided in Sec. IV. II. THEORETICAL ASPECTS A. IODKH method

This subsection briefly summarizes the theory for the energy calculation with the two-component relativistic IODKH scheme. Many-electron two-component relativistic Hamiltonian is expressed as   h+ g++ (1) H2 = 2 (i) + 2 (i, j ), i

i>j

++ are one- and two-particle Hamiltonians, where h+ 2 and g2 respectively. As the one-particle Hamiltonian h+ 2 , we adopted the IODKH Hamiltonian given by † h+ 2 =  G,

 = (12 + Y† Y)−1/2 ,

(3)

and the first-step transformed Hamiltonian, G, is expressed as   p2 12 Y G = p2 b12 + Y† 1 − ep + M† VM + † (σ · pVσ · p), M = K(12 + bpY),  = K(αb12 − p−1 Y), K=

(4)

(10)

Here, c is the speed of light; σ , the vector of three Pauli spin matrices; p, the momentum operator; and V, the sum of the nucleus-electron attraction operators with respect to atom A,  i.e., V = A VA ; and 1n , the n × n identity matrix. Furthermore, Y operator is obtained by the decoupling condition as ep Y + Yep = α 3 (pKbVK − p−1 Kσ · pVσ · pbK) + α 2 (p−1 Kσ · pVσ · pp−1 KY − YKVK) + α 4 (pKbVbKpY − YKbσ · pVσ · pbK) + α 3 Y(Kbσ · pVσ · pKp−1 − KVKbp)Y. (11) The present study focuses on the relativistic effect on the one-particle interaction. Thus, the two-particle Hamiltonian g++ 2 is replaced as the non-relativistic (NR) Coulomb interaction: 1 g(i, j ) = . (12) rij The present treatment is denoted by IODKH/C. Note that we can extend relativistic two-particle interactions such as the Breit-Pauli approximation,55 Foldy-Wouthuysen,31 and IODKH transformation.37 In a practical calculation of Eq. (2), direct evaluation of the matrix elements of χμ |h+ 2 |χν  with respect to the atomic orbitals (AOs) {χ } is remarkably difficult. Thus, a matrix transformation method and the resolution-of-identity (RI) technique proposed by Hess are adopted for the calculation.56 The first step in Hess’ technique is the generation of the eigenvectors of p2 , {k}, which are expressed by the linear combination of the primitive functions {φ},  |k = |φr  Crk , (13) r

where the coefficient Crk is determined by diagonalizing the matrix with the elements φ r |p2 |φ s . The obtained (pseudo) eigenvectors, {k}, satisfy k|p2 |k   = ωk δkk ,

k  ,k 

(5) (6)

,

(7) (8)

(14)

where ωk is a (pseudo) eigenvalue. The explicit matrix expressions of Eqs. (2)–(11) are written as  †  + h2 kk = kk Gk k k k , (15) kk = k|(12 + Y† Y)−1/2 |k  ,

1/2

1 b= , ep + 1

α = 1/c.

and

with

ep + 1 2ep

(9)

(2)

where the second-step transformation, , is defined as



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

Gkk = k|p2 b12 |k   + +

(16)

  † p2 Ykk k  | 12 |k  Yk k 1 − ep k  ,k 

 † Mkk k  |V|k  Mk k k  ,k 

† + kk k  |σ · pVσ · p|k   k k

 ,

(17)

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: 134.99.128.41 On: Tue, 31 Dec 2013 15:59:03

244107-3

Nakajima, Seino, and Nakai

J. Chem. Phys. 139, 244107 (2013)

Here, the approximate RI technique is employed, which is expressed as  |k k| ≈ 1. (21)

with Mkk = Kk (1 + bk pk Ykk ) ,

(18)

  kk = Kk αbk − pk−1 Ykk ,

(19)

k

The matrix elements of functions of p2 , f(p2 ), such as K, b, and p operators in Eqs. (7)–(9) become

and

Ykk

 1 = α 3 (pk Kk bk k|V|k  Kk ek + ek

k|f (p2 )|k   = f (ωk )δkk .

(22)

B. LUT-IODKH method

− pk−1 Kk k|σ · pVσ · p|k  bk Kk ) + α2

In this subsection, the basic concept of an effective unitary transformation scheme, the LUT method with the oneelectron IODKH Hamiltonian, is presented. The LUT scheme is designed on the basis of the locality of the relativistic effect. In the LUT scheme, the entire one-particle unitary transformation U is approximated as a block-diagonal form of the subsystem contributions

 pk−1 Kk k|Kσ · pVσ · p|k  pk−1  Kk  Yk  k  k 

− Ykk Kk k  |V|k  Kk + α4





(pk Kk bk k|V|k  bk Kk pk Yk k

UTotal ≈ UA ⊕ UB ⊕ UC · · · ,

k 

where {A, B, C, · · ·} denotes the subsystems, which are adopted as atomic partitioning in this article. Considering the analytical energy gradient, atomic partitioning is essential for simplifying the formulas as explained in Sec. II E. The relativistic effect of nuclear attraction is dominant on each atom and its interactions within a small distance. The relativistic effect of the kinetic energy is dominant on each atom. The matrix representation of the one-particle Hamiltonian with basis function {χ } is expressed as

− Ykk Kk bk k  |σ · pVσ · p|k  bk Kk ) + α3



Ykk (Kk bk k  |σ · pVσ · p|k  Kk pk−1 

k  ,k 

− Kk k  |V|k  Kk bk pk )Yk k .

(23)

(20)

⎧    B χ ⎪ χμA  T+ + V+ (A = B) + VNR ⎪ C ν A ⎪ ⎪ ⎪ C=A ⎪   B

A  +  B A  LUT  B ⎨ χ A  V+ + V+ + TNR +  (A =  B, RAB ≤ τ ) VNR μ C χν A B χμ  h2 χν ≈ χμ  h2 χν = ⎪ C=A,B ⎪ ⎪

A  NR  NR  B ⎪ ⎪  VC χν (A = B, RAB > τ ), ⎪ ⎩ χμ T +

(24)

C

where TNR and VNR C denote the standard NR kinetic and nuclear-attraction operators: TNR =

1 2 p 12 2

(25) C. Analytical energy derivative for the IODKH method

and VNR C = VC . T+ and V+ C are the relativistic operators given by     p2 T+ = † p2 b12 + Y† 12 Y  1 − ep and

 †  † † V+ C =  M VC M +  (σ · pVC σ · p)  .

information. The validity of such treatments in the energy calculations has been confirmed numerically in the previous article.43

(26)

(27)

(28)

Note that Y, , M, and  operators are determined within individual subsystems because U contains only subsystem

This subsection provides an explicit expression of the analytical energy gradient for the IODKH/C method: relativistic h+ 2 for one-particle Hamiltonian and non-relativistic g for two-particle interaction. In the Hartree-Fock method, the first derivative of the energy with respect to a coordinate RA of nucleus A is written as      ∂(h+ 1  ∂(μλ νρ) ∂E 2 )μν + = Dμν Dμν Dλρ ∂RA ∂RA 2 μ,ν,λ,ρ ∂RA μ,ν   ∂Vnuc   ∂Sμν + − εi Dμν , (29) ∂RA ∂RA 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: 134.99.128.41 On: Tue, 31 Dec 2013 15:59:03

244107-4

Nakajima, Seino, and Nakai

J. Chem. Phys. 139, 244107 (2013)

where μ, ν, λ, and ρ are the indices of {χ }; (μλ||νρ), antisymmetric electron repulsion integral for {χ }; D, the density matrix; Vnuc , the nuclear repulsion potential; ε, the orbital energy; and S, the overlap matrix. In the expression of Eq. (29), ∂(μλ||νρ)/∂RA , ∂Vnuc /∂RA , and ∂Sμν /∂RA are identical to those in the NR theory. The IODKH Hamiltonian in {χ } and that in {k} have the following relation † † + h+ 2 ({χ }) = d X h2 ({k})Xd,

(30)

X = ,

(31)

and

     ∂ek 1 ∂ek ∂Ykk Ykk − Ykk = − ∂RA ek + ek ∂RA ∂RA ∂ (pk Kk bk k|V|k  Kk ) ∂RA  ∂  −1 pk Kk k|σ · pVσ · p|k  bk Kk −α 3 ∂RA  ∂   +α 2 pk−1 Kk k|σ · pVσ · p|k  pk−1  Kk  Yk  k  ∂RA k 

+α 3

with

where d denotes the contraction coefficients, and h+ 2 ({χ }) and h+ 2 ({k}) are the IODKH Hamiltonians in {χ } and {k}, respectively.  is the transformation matrix from space {χ } to orthonormalized space, while  is the transformation matrix from orthonormalized space to {k}. The analytical derivative for h+ 2 ({χ }) is expressed as  †   +  ∂h+ ∂X † ∂h2 ({k}) 2 ({χ }) h+ X = d† X + X 2 ∂RA ∂RA ∂RA  ∂X +X† h+ d, (32) 2 ∂RA with ∂h+ 2 ({k}) = ∂RA



∂† ∂RA





∂G ∂RA



∂ . ∂RA (33) Here, we can obtain the derivatives of the matrix elements of  and G by differentiating individual matrix elements with respect to {k} in Eqs. (16)–(20),  †    ∂Y ∂Y 1 ∂kk † −3/2 † Y+Y |k   = k|− (12 + Y Y) ∂RA 2 ∂RA ∂RA †

G + 

and ∂Gkk ∂ = k|p2 b12 |k   ∂RA ∂RA   ∂  † p2   Ykk k | 12 |k Yk k + ∂RA 1 − ep k  ,k  +

 k  ,k 

+

 k  ,k 

+α 4

 ∂ (pk Kk bk k|V|k  bk Kk pk Yk k ) ∂R A k 

−α 4

 ∂ (Ykk Kk bk k  |σ · pVσ · p|k  bk Kk ) ∂R A  k

+α 3

 ∂  Ykk Kk bk k  |σ · pVσ · p|k   ∂R A k  k 

× k  |V|k  Kk bk pk Yk k ) .

(38)

Note that vectors such as K, b, and p in Eqs. (34)–(38) have information about molecular structure through the function {k}. Thus, the elements should be differentiated as  2 ∂pk ∂Kk 1 , (39) = − α 2 ek−3 Kk−1 ∂RA 8 ∂RA   ∂bk 1 2 −1 2 ∂pk2 , (40) = − α ek bk ∂RA 2 ∂RA  2 ∂pk ∂ek 1 , (41) = − α 2 ek−1 ∂RA 2 ∂RA and ∂pk−1 1 = − pk−3 ∂RA 2

∂ † (M  k  |V|k  Mk k ) ∂RA kk



∂pk2 ∂RA

 .

(42)

D. Analytical derivative for the space transformation matrices

∂ † (  k  |σ · pVσ · p|k   k k ), ∂RA kk (35)

with ∂Mkk ∂ ∂Kk = (1 + bk pk Ykk ) + Kk (bk pk Ykk ), ∂RA ∂RA ∂RA

 ∂ (Ykk Kk k  |V|k  Kk ) ∂R A  k

 ∂  3 (Ykk Kk ×Kk pk−1  Yk  k  − α ∂R A   k k

 + † G

(34)

−α 2

(36)

  ∂ kk ∂  ∂Kk  αbk − pk−1 Ykk + Kk αbk − pk−1 Ykk , = ∂RA ∂RA ∂RA (37)

This subsection presents the explicit derivative of space transformation matrix X and ∂pk2 /∂RA . This is the same procedure as that in Ref. 38. First, the procedure for the derivative of  in Eq. (31) is explained.  is determined by the condition  T S = 1,

(43)

with the overlap matrix S. In more details, S is first diagonalized as WT SW = s.

(44)

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: 134.99.128.41 On: Tue, 31 Dec 2013 15:59:03

244107-5

Nakajima, Seino, and Nakai

J. Chem. Phys. 139, 244107 (2013)

 is obtained using unitary transformation W and eigenvalues of s,  = Ws−1/2 .

(45)

Here, the derivative of T({k}) is expressed as   ∂T({φ}) ∂T({k})  − [PV , T({k})]− , = † ∂RA ∂RA

The explicit expression of the derivative of  is as follows:     ∂s ∂W −1/2 1 ∂ −3/2 s . (46) = − Ws ∂RA ∂RA 2 ∂RA

with

The derivative of s is given by       ∂s ∂S ∂WT ∂W T T SW + W W+W S = ∂RA ∂RA ∂RA ∂RA

Note that ∂T({k})/∂RA has a relationship with ∂p/∂RA :

 =W

∂S ∂RA

 W − [PW , s]− ,



∂ ∂RA

 .

∂T ∂p =2 . ∂RA ∂RA

(47) T

PV = †

(48)

(55)

(56)

(57)

To obtain ∂T({φ})/∂RA , Eq. (53) is differentiated with respect to RA :       ∂TW −1/2 1 −1 ∂s ∂T({φ}) s , T({χ }) , s = s−1/2 − ∂RA ∂RA 2 ∂RA +

with

(58) ∂W PW = WT . ∂RA

(49)

TW = WT TNR W.

Here, square bracket represents the anti-commuter/commuter of matrices as [A, B]± = AB ± BA.

(50)

Note that the offdiagonal elements of ∂s/∂RA and the diagonal elements of PW are zero because s is a diagonal matrix and PW is an antisymmetric matrix, respectively. From Eq. (48) and these conditions, the matrix expression of PW is obtained as (PW )kk

{WT (∂S/∂RA )W}kk = (1 − δkk ). sk − sk 

Equation (58) is rewritten using the derivative of TW :  NR  ∂TW T ∂T W − [PW , TW ]− . =W ∂RA ∂RA

∂W = WPW , ∂RA

(52)

∂/∂RA is determined from Eq. (46) using ∂s/∂RA and ∂W/∂RA . Let us explain the procedure for the derivative of  in Eq. (31). The first step is a transformation of TNR from {φ} to orthonormalized space as T({φ}) =  † TNR .

(53)

By diagonalizing T({φ}),  is obtained: T({k}) = † T({φ}).

(54)

(PV )kk =

{† (∂TW /∂RA ) }kk (1 − δkk ). T({k})k − T({k})k

(60)

(61)

Finally, the derivative of  becomes ∂ = PV . ∂RA

(62)

Using the derivatives of  and , the derivative of space transformation matrix X can be obtained:     ∂ ∂ ∂X + . (63) = ∂RA ∂RA ∂RA E. Analytical derivative for the LUT-IODKH method

This subsection explains an analytical derivative of the LUT Hamiltonian for an efficient gradient calculation. The derivative expression of the LUT-IODKH matrix of Eq. (24) is given by

  NR  B χμA  VC χν



(A = B)

C=A

  + B χμA  V+ A + VB χν   B ∂ A  NR  + χμ T + VNR C χν ∂RA C=A,B

A  NR  NR  B χμ  T + VC χν

(59)

Let us think about the properties of Eq. (55). Offdiagonal elements on the left-hand side in Eq. (55) are zero because T({k}) is diagonal. The diagonal elements of PV are zero due to the antisymmetry of PV . These properties provide

(51)

Thus, ∂s/∂RA is obtained by substituting Eq. (51) to Eq. (48). Since ∂W/∂RA is also given by

⎧ ∂ ⎪ ⎪ ⎪ ⎪ ∂RA ⎪ ⎪ ⎪ ⎪ ⎪ ∂ ⎪ ⎪ ∂ A  LUT  B ⎨ ∂RA χν = χ h ⎪ ∂RA μ 2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ∂ ⎪ ⎪ ⎪ ⎩ ∂R A

where

(A = B, RAB ≤ τ )

(64)

(A = B, RAB > τ ),

C

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: 134.99.128.41 On: Tue, 31 Dec 2013 15:59:03

244107-6

where ∂V+ C ∂RA

Nakajima, Seino, and Nakai

J. Chem. Phys. 139, 244107 (2013)

    ∂ † † ∂VC † = M M+ σ · pVC σ · p  . ∂RA ∂RA (65)

Note that the IODKH transformation is only applied to nucleus-electron interaction for different atoms with a small distance. If A = B, the derivative of the IODKH transformed matrices such as T+ and V+ C become zero because the transformations of LUT-IODKH have closure forms within individual atoms. If A = B with RAB ≤ τ , the nuclear attraction elements V+ C for the LUT-IODKH transformation are obtained by differentiation with respect to RA as well as the NR kinetic terms TNR and nuclear attraction ones VNR . F. Implementation

This subsection presents the scheme for evaluating the energy gradient with the (LUT-) IODKH Hamiltonian. Table I summarizes the procedures for (a) NR/C, (b) IODKH/C, and (c) LUT-IODKH/C. For the NR/C case, the density matrix D and orbital energy ε obtained from energy calculation are preserved in the energy gradient calculation. Then the derivatives of NR one-electron integrals TNR and VNR are calculated in contracted formalism. For the conventional IODKH case, integrals and momentum factors in uncontracted formalism such as S, TNR , VNR , , Y, G, h, M, , K, b, ep , and p−1 as well as the density matrix D and orbital energy ε obtained from energy calculation are preserved in energy gradient calculation. Furthermore, the derivatives of NR and relativistic integral, S, TNR , VNR , and σ · pVσ · p, are calculated in uncontracted formalism. Using S and TNR , the derivatives of transformation matrices  and  are calculated using Eqs. (43)–(62). At the same time, the derivatives of kinematic factors of Eqs. (39)–(42)

can be calculated. Using kinematic factors and its derivatives, the derivatives of , G, and Y in {k} are calculated using Eqs. (34), (35), and (38), respectively. Finally, the derivative of relativistic one-particle Hamiltonian in Eq. (33) is obtained by back transformation into {χ } using the derivatives of  and . For the LUT-IODKH case, the procedure becomes simpler than the procedure of the conventional IODKH. First, the transformation factors , Y, M, , K, b, ep , and p−1 in each subsystem as well as D and ε are preserved in energy calculation. Then, the NR integrals for TNR and VNR are calculated in contracted formalism. In order to obtain the transformed integrals for interatomic interactions with RAB ≤ τ in Eq. (64), the derivatives of VNR and σ · pVσ · p are calculated in uncontracted formalism. After that, the transformation is performed using Eq. (65). Thus, the derivative of LUT-IODKH Hamiltonian matrix can be obtained. Note that this complicated procedure involving Eqs. (43)–(62) is unnecessary in the case of the LUT-IODKH scheme. In addition, the range of transformation is extremely small. Consequently, the LUT-IODKH scheme is suitable for the efficient calculation for analytical derivative of one-electron integral.

III. NUMERICAL ASSESSMENTS A. Computational details

This section presents the assessment of the performance of the analytical energy gradient for the IODKH and LUTIODKH schemes. We adopted the spin-free formalism of IODKH/C: the one-electron Dirac Hamiltonian and twoelectron non-relativistic Coulomb interaction. For LUT, the threshold for cutoff of relativistic interaction, τ , was set to 3.5 Å. For comparison, the spin-free 4c Dirac-Coulomb and non-relativistic (NR/C) methods were also calculated.

TABLE I. Procedure for computing the analytical energy gradient of (a) NR/C, (b) IODKH/C, and (c) LUT-IODKH/C.

(I)

From energy calculation

(II)

Derivatives of NR one-electron integrals

(III)

Derivatives of NR one-electron integrals for relativistic transformation (uncontracted form) Derivatives of transformation matrices

(IV) (V)

(VI) (VII)

(a) NR/C

(b) IODKH/C

(c) LUT-IODKH/C

D, ε

D, ε S, TNR , VNR , , Y, G, h, M, , K, b, ep , p−1 (for whole subsystem) ...

D, ε , Y, M, , K, b, ep , p−1 (for each system) ∂TNR ∂VNR ∂RA , ∂RA (except for interatomic interaction with RAB ≤ τ ) ∂VNR ∂(σ · pVσ · p) ∂RA , ∂RA (for interatomic interaction with RAB ≤ τ ) ...

∂TNR ∂VNR ∂RA , ∂RA

Derivatives of relativistic one-electron integrals

Derivatives of NR two-electron integral ∂(μλ||νρ) ∂RA Energy gradient calculation from Eq. (29).

...

... ...

∂S ∂TNR ∂VNR ∂(σ · pVσ · p) ∂RA , ∂RA , ∂RA , ∂RA

∂ ∂ ∂RA , ∂RA [from Eqs. (46) and ∂ep ∂p−1 ∂K ∂b ∂RA , ∂RA , ∂RA , ∂RA

(62)]

[from Eqs. (39)–(42)] ∂Y ∂RA [from Eq. (38)] ∂M ∂ ∂RA , ∂RA [from Eqs. (36) and (37)] ∂ ∂G ∂RA , ∂RA [from Eqs. (34) and (35)]

∂V+ ∂RA

[from Eq. (65)] (for interatomic interaction with RAB ≤ τ )

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: 134.99.128.41 On: Tue, 31 Dec 2013 15:59:03

244107-7

Nakajima, Seino, and Nakai

J. Chem. Phys. 139, 244107 (2013) TABLE II. Basis sets for diatomic molecules and 20 metal complexes used at the levels of HF, MP2, and B3LYP in NR/C, IODKH/C, and large component of 4c. Period 1 2 3 4

5 FIG. 1. Geometry of (a) one-dimensional silver cluster Ag10 , (b) twodimensional silver cluster Ag13 (2 × 2), and (c) three-dimensional silver cluster Ag63 (2 × 2 × 2).

The calculations of parts of small-component integrals, i.e., (SS|SS)-type ones, were skipped in the 4c method.57 Numerical tests were performed for various closed-shell diatomic molecules: hydrogen halides HX, halogen dimmers X2 (X = F, Cl, Br, I, and At), coinage metal hydrides MH, and coinage metal dimers M2 (M = Cu, Ag, and Au). In addition, one-, two-, and three-dimensional silver clusters Agn were calculated to investigate the efficiency of the LUT scheme. Figure 1 illustrates the structures of Agn clusters. For onedimensional silver clusters with Agn (n = 1, 2, · · ·, 10), the Ag–Ag bond distance used here was 2.706 Å. For two- and three-dimensional silver clusters, a face-centered cubic structure with the experimental lattice constant, 4.086 Å,58 was adopted. The unit cells for two- and three-dimensional clusters contain 5 and 14 atoms, respectively. We adopted several sizes of Agn clusters; the total numbers of atoms in the twodimensional clusters were 5, 13, 25, 41, 61, 85, 113, 145, and 181, which correspond to 1 × 1, 2 × 2, . . . , and 9 × 9 unit cells, respectively. The total numbers of atoms in the threedimensional clusters were 14, 23, 32, 38, 53, 63, 74, 88, 123, and 172, which correspond to 1 × 1 × 1, 2 × 1 × 1, 3 × 1 × 1, 2 × 2 × 1, 3 × 2 × 1, 2 × 2 × 2, 3 × 3 × 1, 3 × 2 × 2, 3 × 3 × 2, and 3 × 3 × 3. For more practical systems, 20 metal complexes including the fourth–sixth row-transition metals were examined. The calculations were performed at the HF, MP2, and density functional theory (DFT) with B3LYP functional levels.59 The basis sets for all-electron calculation in diatomic molecules and 20 metal complexes are summarized in Table II. All basis sets were used in uncontracted form. For the 4c method, the kinetically balanced small-component basis sets were adopted. MCP/NOSeC-VTZP60–70 and Stuttgart/Dresden (SDD71–74 ) type effective core potential are used in RECP. In SDD calculation, the first–third row atoms were described with cc-pVTZ.75, 76 The IODKH energy gradient method with/without the LUT scheme was implemented in the GAMESS program.77 The 4c calculations were performed in the DIRAC12 program.78 The Gaussian09 program was used in SDD calculations.79

6

Element

Description

Basis sets

H C, N, O, F P, Cl Ti, V Cr Co, Ni Cu Zn Br Zr Nb Mo Ru Ag Cd I Ta, W, Re Os, Pt, Au, Hg At

(6s3p2d) (12s8p3d2f) (16s11p4d2f) (19s13p7d3f) (20s13p8d3f) (20s12p9d3f) (20s12p8d3f) (20s13p9d3f) (20s15p10d2f) (22s17p13d3f) (22s17p12d3f) (22s16p11d3f) (21s16p12d3f) (22s16p12d3f) (22s16p13d3f) (22s18p13d2f) (26s22p16d12f) (26s21p16d12f) (26s23p17d13f)

Ref. 83 Refs. 61 and 84 Refs. 62 and 84 Refs. 64 and 84 Refs. 64 and 84 Refs. 64 and 84 Refs. 64 and 84 Refs. 64 and 84 Refs. 65 and 84 Refs. 66 and 84 Refs. 66 and 84 Refs. 66 and 84 Refs. 66 and 84 Refs. 66 and 84 Refs. 66 and 84 Refs. 65 and 84 Refs. 68 and 84 Refs. 68 and 84 Refs. 67 and 84

B. Numerical gradient values

This subsection investigates the accuracy of the analytical gradient values for IODKH/C and LUT-IODKH/C Hamiltonians by comparing with the numerical ones in HAt, At2 , AuH, and Au2 molecules. Table III compares between the analytical and numerical energy gradient values calculated by the NR/C, IODKH/C, LUT-IODKH/C, and 4c at the HF level. We examined the gradient values in 1.0, 1.5, 2.0, and 3.0 times of the experimental bond length, re , which were 1.742, 3.046, 1.5324, and 2.4719 Å for HAt,80 At2 ,5 AuH,81 and Au2 ,82 respectively. The numerical gradient calculations were performed using a central difference method with 0.01 Å stepsize. The differences between the analytical g ana and numerical g num energy gradient values, that is, g ana , are shown in parentheses, which represent the correctness of the analytical scheme. The deviations from 4c, g 4c , are given in square brackets, which represent the two-electron relativistic effect of contribution of (SS|SS)-type integrals. In IODKH/C, the absolute deviations from the numerical results, i.e., g ana , are less than 0.0004 hartree/bohr. Such small deviations confirm the correctness of the present implementation for the analytical energy gradient of IODKH/C. The deviations from 4c, i.e., | g 4c |, are comparatively large in At2 and Au2 . This is due to the lack of two-electron relativistic effect in 2c or (SS|SS)-type integral in 4c. In LUTIODKH/C, | g ana | are also small. The largest deviation is 0.0001 hartree/bohr. In addition, the absolute values of differences between IODKH/C and LUT-IODKH/C, which are given in | g LUT |, are less than 0.0003 hartree/bohr. This means that the LUT treatment evaluates the relativistic energy gradient without loss of accuracy.

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: 134.99.128.41 On: Tue, 31 Dec 2013 15:59:03

244107-8

Nakajima, Seino, and Nakai

J. Chem. Phys. 139, 244107 (2013)

TABLE III. Analytical and/or numerical gradient values (hartree/bohr) of HAt, At2 , AuH, and Au2 calculated by NR/C, IODKH/C, LUT-IODKH/C, and 4c at HF level. NR/C Molecules HAt

At2

AuH

Au2

IODKH/C

LUT-IODKH/C

4c

Lengtha

g ana

g ana

g num

g anab

g 4cc

g ana

g num

g anab

g LUTd

g ana

1.0re 1.5re 2re 3re 1.0re 1.5re 2re 3re 1.0re 1.5re 2re 3re 1.0re 1.5re 2re 3re

0.0113 0.0675 0.0379 0.0098 0.0256 0.0273 0.0082 0.0019 0.0823 0.0336 0.0313 0.0149 0.0798 0.0170 0.0127 0.0045

0.0192 0.0657 0.0366 0.0096 0.0304 0.0259 0.0079 0.0019 0.0135 0.0540 0.0365 0.0138 0.0305 0.0270 0.0134 0.0036

0.0195 0.0654 0.0365 0.0097 0.0305 0.0255 0.0078 0.0015 0.0135 0.0540 0.0365 0.0138 0.0304 0.0270 0.0133 0.0037

( − 0.0003) (0.0003) (0.0001) ( − 0.0001) ( − 0.0001) (0.0004) (0.0001) (0.0004) (0.0000) (0.0000) (0.0000) (0.0000) (0.0001) (0.0000) (0.0001) ( − 0.0001)

[ − 0.0002] [0.0001] [ − 0.0001] ... [ − 0.0110] [ − 0.0047] [ − 0.0027] [ − 0.0012] [0.0011] [0.0000] [ − 0.0001] [ − 0.0002] [0.0117] [ − 0.0053] [ − 0.0031] [ − 0.0013]

0.0194 0.0655 0.0365 0.0096 0.0304 0.0256 0.0078 0.0019 0.0135 0.0540 0.0365 0.0138 0.0305 0.0268 0.0133 0.0036

0.0193 0.0655 0.0366 0.0096 0.0305 0.0256 0.0078 0.0020 0.0135 0.0541 0.0365 0.0138 0.0305 0.0269 0.0133 0.0035

(0.0001) (0.0000) ( − 0.0001) (0.0000) ( − 0.0001) (0.0000) (0.0000) ( − 0.0001) (0.0000) ( − 0.0001) (0.0000) (0.0000) (0.0000) ( − 0.0001) (0.0000) (0.0001)

[0.0002] [ − 0.0002] [ − 0.0001] [0.0000] [0.0000] [ − 0.0003] [ − 0.0001] [0.0000] [0.0000] [0.0000] [0.0000] [0.0000] [0.0000] [ − 0.0002] [ − 0.0001] [0.0000]

0.0194 0.0656 0.0367 N. C.e 0.0414 0.0306 0.0106 0.0031 0.0124 0.0540 0.0367 0.0140 0.0188 0.0323 0.0165 0.0048

a

re denotes the experimental bond length taken from Refs. 5 and 80–82. Difference in analytical gradient value from the numerical one is given in parentheses. c Difference in analytical gradient value from 4c is given in square brackets. d Difference between IODKH/C and LUT-IODKH/C in analytical gradient value is given in square brackets. e SCF calculations were not converged. b

C. Accuracies of IODKH/C and LUT-IODKH/C methods

This subsection examines the accuracy of the analytical energy gradient scheme for IODKH/C Hamiltonian. Table IV shows the bond lengths (Å) of diatomic molecules calculated by the NR/C, IODKH/C, and 4c methods at the HF level. The difference r rel of IODKH/C and 4c methods from NR/C are given in square brackets, representing the relativistic effect accounted by IODKH/C and 4c. The deviations r 4c of TABLE IV. Bond length (Å) of HX, X2 (X = F, Cl, Br, I, and At), MH, and M2 (M = Cu, Ag, and Au) molecules obtained by NR/C, IODKH/C, and 4c methods at HF level. NR

IODKH/C

4c

Mole.

r NR

r IODKH

r rela

r 4cb

r 4c

r rela

HF HCl HBr HI HAt F2 Cl2 Br2 I2 At2 CuH AgH AuH Cu2 Ag2 Au2

0.897 1.264 1.470 1.608 1.711 1.328 1.975 2.274 2.675 2.890 1.568 1.778 1.828 2.440 2.818 2.924

0.896 1.263 1.404 1.600 1.684 1.328 1.976 2.272 2.663 2.854 1.540 1.700 1.568 2.401 2.706 2.604

[ − 0.001] [ − 0.001] [ − 0.066] [ − 0.008] [ − 0.027] [0.000] [0.001] [ − 0.002] [ − 0.012] [ − 0.036] [ − 0.028] [ − 0.078] [ − 0.260] [ − 0.039] [ − 0.112] [ − 0.320]

( − 0.001) ( − 0.001) ( − 0.001) (0.000) (0.001) (0.000) (0.000) ( − 0.002) ( − 0.002) (0.010) (0.000) (0.001) ( − 0.001) ( − 0.001) (0.000) (0.002)

0.897 1.264 1.405 1.600 1.683 1.328 1.976 2.274 2.665 2.844 1.540 1.699 1.569 2.402 2.706 2.602

[0.000] [0.000] [ − 0.065] [ − 0.008] [ − 0.028] [0.000] [0.001] [0.000] [ − 0.010] [ − 0.046] [ − 0.028] [ − 0.079] [ − 0.259] [ − 0.038] [ − 0.112] [ − 0.322]

a b

Relativistic correction is in square brackets. Difference in bond length from that obtained using 4c is in parentheses.

IODKH/C from 4c are shown in parentheses, indicating the accuracy of the IODKH/C treatment. The relativistic effect on bond lengths, r rel , is negative except for F2 and Cl2 . The absolute values increase as heavier elements are added to the systems. For example, the values in halogen dimers are 0.000, 0.001, −0.002, −0.012, and −0.036 Å for F2 , Cl2 , Br2 , I2 , and At2 , respectively. This might be the common tendency in the relativistic effect. Furthermore, the deviations in halogen molecules are smaller than those in coinage metal molecules. This is because the coinage metal molecules have σ -bonding formed by each s valence orbital, which is largely shrunken by the relativistic effect in comparison with σ -bonding formed by p valence orbitals seen in halogen molecules. The deviations of IODKH/C from 4c, r 4c , are comparatively small. The largest deviation is 0.010 Å for At2 . This indicates that the spin-free two-electron relativistic corrections reasonably demonstrate the bond lengths of elements, at least, up to the sixth row of the periodic table. Table V shows the bond lengths (Å) of diatomic molecules calculated by the NR/C and IODKH/C at the HF, MP2, and B3LYP levels. The differences r corr of NR/C and IODKH/C methods with respect to the corresponding HF results are given in square brackets, representing the correlation effect accounted by NR/C and IODKH/C. The deviations r exp of NR/C and IODKH/C methods from reference are shown in parentheses, indicating the accuracy of the NR/C and IODKH/C treatment. In additions, the mean absolute error (MAE) and the mean error (ME) are given for r exp . In the MP2 results, the correlation effect on bond lengths, r corr , of IODKH/C shows the same tendency as NR/C except for HBr and HAt. Hydrogen halides except for HAt and

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: 134.99.128.41 On: Tue, 31 Dec 2013 15:59:03

244107-9

Nakajima, Seino, and Nakai

J. Chem. Phys. 139, 244107 (2013)

TABLE V. Bond length (Å) of HX, X2 (X = F, Cl, Br, I, and At), MH, and M2 (M = Cu, Ag, and Au) molecules calculated by the NR/C and IODKH/C at HF, MP2, and B3LYP level. HF

MP2

NR/C IODKH/C Molecules HF HCl HBr HI HAt F2 Cl2 Br2 I2 At2 CuH AgH AuH Cu2 Ag2 Au2

r NR 0.897 1.264 1.470 1.608 1.711 1.328 1.975 2.274 2.675 2.890 1.568 1.778 1.828 2.440 2.818 2.924

r IODKH 0.896 1.263 1.404 1.600 1.684 1.328 1.976 2.272 2.663 2.854 1.540 1.700 1.568 2.401 2.706 2.604

NR/C r NR 0.917 1.270 1.412 1.611 1.714 1.400 1.985 2.284 2.677 2.892 1.451 1.662 1.709 2.240 2.591 2.710

r corra [0.020] [0.006] [−0.058] [0.003] [0.003] [0.072] [0.010] [0.010] [0.002] [0.002] [−0.117] [−0.116] [−0.119] [−0.200] [−0.227] [−0.214]

MAE ME a b

B3LYP IODKH/C

r expb (−0.001) (−0.005) (−0.003) (0.001) ... (−0.012) (−0.003) (0.003) (0.011) ... (−0.015) (0.044) (0.177) (0.021) (0.061) (0.238)

r IODKH 0.917 1.270 1.410 1.603 1.682 1.400 1.984 2.281 2.669 2.860 1.424 1.584 1.495 2.203 2.505 2.468

r corra [0.021] [0.007] [0.006] [0.003] [−0.002] [0.072] [0.008] [0.009] [0.006] [0.006] [−0.116] [−0.116] [−0.073] [−0.198] [−0.201] [−0.136]

0.042 0.037

NR/C

r expb (−0.001) (−0.005) (−0.005) (−0.007) ... (−0.012) (−0.004) (0.000) (0.003) ... (−0.042) (−0.034) (−0.037) (−0.016) (−0.025) (−0.004)

r NR 0.922 1.281 1.425 1.624 1.725 1.398 2.012 2.316 2.708 2.919 1.482 1.697 1.746 2.281 2.681 2.804

0.014 − 0.013

IODKH/C

r corra r expb [0.025] (0.004) [0.017] (0.006) [−0.045] (0.010) [0.016] (0.014) [0.014] ... [0.070] (−0.014) [0.037] (0.024) [0.042] (0.035) [0.033] (0.042) [0.029] ... [−0.086] (0.016) [−0.081] (0.079) [−0.082] (0.214) [−0.159] (0.062) [−0.137] (0.151) [−0.120] (0.332)

r IODKH 0.922 1.281 1.423 1.618 1.696 1.398 2.012 2.313 2.702 2.883 1.459 1.625 1.538 2.248 2.595 2.547

r corra r expb [0.026] (0.004) [0.018] (0.006) [0.019] (0.008) [0.018] (0.008) [0.012] ... [0.070] (−0.014) [0.036] (0.024) [0.041] (0.032) [0.039] (0.036) [0.029] ... [−0.081] (−0.007) [−0.075] (0.007) [−0.030] (0.006) [−0.153] (0.029) [−0.111] (0.065) [−0.057] (0.075)

0.072 0.070

Expt. Reference 0.9176 82 1.2749 82 1.4146 82 1.6099 82 ... ... 1.4119 82 1.9879 82 2.2811 82 2.6663 82 ... ... 1.4658 81 1.6179 81 1.5324 81 2.2192 85 2.5303 86 2.4719 82

0.023 0.020

Difference in bond length from the corresponding HF result is given in square brackets. Difference in bond length from the experimental result is given in parentheses.

halogen dimers have positive values, whereas coinage metal halides and coinage metal dimers have negative values. The deviations of NR/C from reference, r exp , increase as the heavier elements are present in the molecules. On the other hand, the deviations of IODKH/C from reference remain approximately constant. For example, the values in coinage metal halides are −0.042, −0.034, and −0.037 Å for CuH,

AgH, and AuH, respectively. Furthermore, MAE and ME of IODKH/C are relatively smaller than those of NR/C. This indicates that both correlation and relativistic effects are essential for the molecules containing heavier elements. In the B3LYP results, the correlation effects on bond lengths, r corr , of IODKH/C and NR/C show similar trends except for HBr. The deviations of NR/C from reference,

TABLE VI. Bond length (Å) of HX, X2 (X = F, Cl, Br, I, and At), MH, and M2 (M = Cu, Ag, and Au) molecules calculated by the IODKH/C and LUTIODKH/C at HF, MP2, and B3LYP level. HF Molecules HF HCl HBr HI HAt F2 Cl2 Br2 I2 At2 CuH AgH AuH Cu2 Ag2 Au2 MAE ME a

MP2

B3LYP

r IODKH

r LUT

r LUTa

r IODKH

r LUT

r LUTa

r IODKH

r LUT

r LUTa

0.896 1.263 1.404 1.600 1.684 1.328 1.976 2.272 2.663 2.854 1.540 1.700 1.568 2.401 2.706 2.604

0.896 1.264 1.404 1.601 1.682 1.328 1.976 2.273 2.664 2.840 1.541 1.700 1.570 2.401 2.705 2.600

(0.000) (0.001) (0.000) (0.001) ( − 0.002) (0.000) (0.000) (0.001) (0.001) ( − 0.014) (0.001) (0.000) (0.002) (0.000) ( − 0.001) ( − 0.004)

0.917 1.270 1.410 1.603 1.682 1.400 1.984 2.281 2.669 2.860 1.424 1.584 1.495 2.203 2.505 2.468

0.917 1.270 1.410 1.604 1.689 1.400 1.985 2.281 2.667 2.847 1.425 1.584 1.491 2.203 2.506 2.462

(0.000) (0.000) (0.000) (0.001) (0.007) (0.000) (0.001) (0.000) ( − 0.002) ( − 0.013) (0.001) (0.000) ( − 0.004) (0.000) (0.001) ( − 0.006)

0.922 1.281 1.423 1.618 1.696 1.398 2.012 2.313 2.702 2.883 1.459 1.625 1.538 2.248 2.595 2.547

0.922 1.281 1.423 1.618 1.705 1.398 2.013 2.313 2.700 2.879 1.460 1.626 1.536 2.247 2.595 2.540

(0.000) (0.000) (0.000) (0.000) (0.009) (0.000) (0.001) (0.000) ( − 0.002) ( − 0.004) (0.001) (0.001) ( − 0.002) ( − 0.001) (0.000) ( − 0.007)

0.002 − 0.001

0.003 − 0.001

0.002 0.000

Difference in bond length from that obtained using the IODKH/C method.

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: 134.99.128.41 On: Tue, 31 Dec 2013 15:59:03

244107-10

Nakajima, Seino, and Nakai

J. Chem. Phys. 139, 244107 (2013)

D. Computational cost of the LUT scheme

O n2.52

O n1.31 n FIG. 2. System-size dependence of CPU time (s) of the one-electron analytical energy gradient in Agn (n = 1, 2, . . . , 10) obtained with conventional IODKH and LUT-IODKH methods. A Xeon X5470/3.33GHz (Quad Core) processor was used on a single core.

r exp , are positive except for F2 . The values increase as the elements in the system are heavier. For example, the values in coinage dimers are 0.062, 0.151, and 0.332 Å for Cu2 , Ag2 , and Au2 , respectively. The deviations of IODKH/C are positive except for F2 and CuH and slightly increase in the case of heavier elements. MAE and ME of IODKH/C are smaller than those of NR/C because both correlation and relativistic effects are considered. Table VI shows the bond length (Å) of diatomic molecules calculated by the IODKH/C and LUT-IODKH/C at the HF, MP2, and B3LYP levels. The deviation r LUT of the LUT-IODKH/C method from IODKH/C is shown in parentheses, indicating the accuracy of the LUT-IODKH/C treatment. The deviations are both positive and negative. The largest absolute deviations in HF, MP2, and B3LYP are 0.014 Å for At2 , 0.013 Å for At2 , and 0.009 Å for HAt, respectively. Since MAE and ME are reasonably small, the LUT-IODKH/C treatment can sufficiently describe the relativistic effect for estimating the bond lengths.

On

In this subsection, the computational costs of the IODKH/C and LUT-IODKH/C methods are investigated on one-, two-, and three-dimensional silver clusters. Figure 2 shows the system-size dependence of the central processing unit (CPU) time for the one-electron analytical energy gradient calculation in the one-electron IODKH and LUT-IODKH transformations of one-dimensional silver clusters Agn (n = 1, 2, · · ·, 10). The vertical axis represents the CPU time in seconds. The horizontal axis shows the system size n. The CPU times were measured with the Xeon X5470/3.33 GHz (Quad Core) on a single core. The computational costs of the IODKH transformation are comparatively large: the scaling is O(n2.52 ). On the other hand, the LUT technique drastically reduces the computational cost compared with the conventional IODKH method; the scaling is O(n1.31 ) with a small prefactor. These scalings reasonably agree well with the theoretical ones, namely, O(n3 ) and O(n1 ), respectively. Figure 3 shows the system-size dependence of the CPU times in the LUT-IODKH transformations of two- and threedimensional silver clusters. Some of linear- (or low-) scaling techniques in quantum chemistry work only for onedimensional systems. However, the LUT treatments for the energy and gradient calculations theoretically scale linear not only one- but also two- or three-dimensional systems. In fact, Figure 3 shows O(n1.21 ) for both two- and three-dimensional cases. Thus, the present LUT-IODKH/C method is confirmed as a highly efficient scheme to evaluate the analytical gradient for arbitrary dimensional systems. Table VII shows the CPU times for the seven steps in the analytical energy gradient calculations at the HF level in an Ag10 cluster: the one-electron integral (OEI), the unitary transformation of the OEI, an initial guess obtained with the Hückel calculation, the two-electron integral (TEI), the SCF procedure, the derivative of OEI (dOEI), and the derivative of TEI (dTEI). The Hamiltonians used here were NR/C, IODKH/C, and LUT-IODKH/C. The percentage of the CPU time required for each step is given in parentheses. In the energy calculations, which include OEI, UT of OEI, Guess, TEI, and SCF, the CPU time for unitary transformation is

O n1.21

1.21

n

FIG. 3. System-size dependence of CPU time (s) of the one-electron analytical energy gradient in multi-dimensional silver clusters Agn obtained with LUTIODKH methods. A Xeon X5470/3.33GHz (Quad Core) processor was used on a single core. (a) Two-dimensional cluster. (b) Three-dimensional cluster.

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: 134.99.128.41 On: Tue, 31 Dec 2013 15:59:03

244107-11

Nakajima, Seino, and Nakai

J. Chem. Phys. 139, 244107 (2013)

TABLE VII. CPU time (s) for seven steps in Hartree-Fock calculation obtained with NR/C, IODKH/C, and LUT-IODKH/C in Ag10 linear cluster: one-electron integral (OEI), unitary transformation (UT) of OEI, initial guess (Guess) calculated with Hückel scheme, two-electron integral (TEI), SCF procedure (SCF), derivative of OEI (dOEI), and derivative of TEI (dTEI). A Quad core Xeon/3.3 GHz processor was used on a single core. NR Time

IODKH/C %

Time

LUT-IODKH/C %

Time

LUT-IODKH/C method can dramatically reduce CPU times. Note that the entire CPU time for the LUT-IODKH/C is close to those of the NR/C method. These results indicate that the LUT-IODKH/C scheme gives accurate results that are close to that of the conventional IODKH/C ones as well as 4c ones with the computational costs that are comparable with that of the NR method.

%

OEI 0.22 (0.01) 0.22 (0.00) 0.23 (0.01) UT of OEI 0.00 (0.00) 41 902.32 (41.56) 12.54 (0.31) Guess 254.07 (6.52) 262.55 (0.26) 267.26 (6.50) TEI 886.03 (22.74) 928.95 (0.92) 886.19 (21.57) SCF 1242.38 (31.89) 2446.54 (2.43) 1328.72 (32.34) dOEI 2.31 (0.06) 53 755.54 (53.31) 15.91 (0.39) dTEI 1510.97 (38.78) 1536.95 (1.52) 1598.19 (38.89) Total 3895.98 (100.00) 100 833.07 (100.00) 4109.04 (100.00)

drastically reduced in LUT-IODKH/C compared with IODKH/C. In the conventional IODKH/C for analytical gradient calculation, the CPU time for the unitary transformation of the dOEI is one of the bottlenecks in the entire procedure, accounting for 53.31% of the CPU time. However, the

E. Metal complexes

This subsection provides optimized geometries for 20 transition metal complexes. Table VIII shows the bond lengths for the 20 transition metal complexes calculated by NR, RECP, MCP, and LUT-IODKH at the B3LYP level. A-B indicates the bond length between A and B. The differences r rel of RECP, MCP, and LUT-IODKH methods from NR/C are given in brackets, representing the relativistic effect accounted by RECP, MCP, and LUT-IODKH. Figure 4 shows the difference in bond length from that obtained using the NR/C method. The bond lengths in the sixth row metal complexes calculated by relativistic methods are shorter than that of NR. Especially in the twelfth group, Zn, Cd, and Hg, the bond lengths are shorter than those in the other groups.

TABLE VIII. Bond lengths (Å) for the 20 transition-metal complexes calculated by NR, RECP, MCP, and LUT-IODKH at the B3LYP level. NR Period 4

5

6

a

RECP

MCP

LUT-IODKH

Group

Complex

Bond

r NR

r RECP

r rela

r MCP

r rela

r IODKH

r rela

4

TiOCl2

5

TiF4 VO(OEt)3

6 9

Cr(CO)6 Co(CO)4 H

10

5

Ni(PF3 )4 Ni(CO)4 ZnMe2 ZrF4 ZrBr4 NbCl5

8 12 5

MoBr4 RuO4 CdEt2 TaCl5

6

WOCl4

7

Re(CO)5 Br

8 10

OsO4 cis-Platin

12

HgMe2

Ti–O Ti–Cl Ti–F V–O V–OEt Cr–C Co–H Co–C(ax.) Co–C(eq.) Ni–P Ni–C Zn–C Zr–F Zr–Br Nb–Cl(ax.) Nb–Cl(eq.) Mo–Br Ru–O Cd–C Ta–Cl(ax.) Ta–Cl(eq.) W–O W–Cl(ax.) W–Cl(eq.) Re–C(ax.) Re–C(eq.) Re–Br Os–O Pt–N Pt–Cl Hg–C

1.598 2.239 1.754 1.577 1.768 1.926 1.472 1.822 1.809 2.131 1.847 1.958 1.906 2.495 2.342 2.234 2.400 1.693 2.201 2.361 2.312 1.721 2.431 2.299 1.994 2.077 2.703 1.737 2.195 2.370 2.247

1.593 2.232 1.749 1.576 1.763 1.920 1.475 1.802 1.816 2.128 1.832 1.945 1.906 2.495 2.339 2.289 2.406 1.682 2.153 2.337 2.288 1.696 2.404 2.266 1.944 2.023 2.656 1.697 2.106 2.307 2.123

( − 0.005) ( − 0.007) ( − 0.005) ( − 0.001) ( − 0.005) ( − 0.006) (0.003) ( − 0.020) (0.007) ( − 0.003) ( − 0.015) ( − 0.013) (0.000) (0.000) ( − 0.003) (0.055) (0.006) ( − 0.011) ( − 0.048) ( − 0.024) ( − 0.024) ( − 0.025) ( − 0.027) ( − 0.033) ( − 0.050) ( − 0.054) ( − 0.047) ( − 0.040) ( − 0.089) ( − 0.063) ( − 0.124)

1.601 2.237 1.755 1.587 1.768 1.931 1.478 1.804 1.824 2.119 1.820 1.933 1.907 2.494 2.342 2.298 2.402 1.687 2.171 2.304 2.264 1.689 2.365 2.241 1.931 2.006 2.662 1.684 2.057 2.294 2.111

(0.003) ( − 0.002) (0.001) (0.010) (0.000) (0.005) (0.006) ( − 0.018) (0.015) ( − 0.012) ( − 0.027) ( − 0.025) (0.001) ( − 0.001) (0.000) (0.064) (0.002) ( − 0.006) ( − 0.030) ( − 0.057) ( − 0.048) ( − 0.032) ( − 0.066) ( − 0.058) ( − 0.063) ( − 0.071) ( − 0.041) ( − 0.053) ( − 0.138) ( − 0.076) ( − 0.136)

1.595 2.236 1.752 1.576 1.764 1.920 1.468 1.811 1.798 2.115 1.830 1.942 1.905 2.488 2.336 2.288 2.390 1.681 2.158 2.306 2.260 1.700 2.405 2.264 1.944 2.022 2.652 1.700 2.086 2.303 2.118

( − 0.003) ( − 0.003) ( − 0.002) ( − 0.001) ( − 0.004) ( − 0.006) ( − 0.004) ( − 0.011) ( − 0.011) ( − 0.016) ( − 0.017) ( − 0.016) ( − 0.001) ( − 0.007) ( − 0.006) (0.054) ( − 0.010) ( − 0.012) ( − 0.043) ( − 0.055) ( − 0.052) ( − 0.021) ( − 0.026) ( − 0.035) ( − 0.050) ( − 0.055) ( − 0.051) ( − 0.037) ( − 0.109) ( − 0.067) ( − 0.129)

12 4

Relativistic correction in brackets.

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: 134.99.128.41 On: Tue, 31 Dec 2013 15:59:03

244107-12

Nakajima, Seino, and Nakai

J. Chem. Phys. 139, 244107 (2013) 1 P.

Pyykkö, Chem. Rev. 88, 563 (1988). Dyall and K. Fægri, Introduction to Relativistic Quantum Chemistry (Oxford University Press, 2007). 3 M. Reiher and A. Wolf, Relativistic Quantum Chemistry (Wiley-VCH, Germany, 2009). 4 J. Autschbach, J. Chem. Phys. 136, 150902 (2012). 5 L. Visscher and K. G. Dyall, J. Chem. Phys. 104, 9040 (1996). 6 M. Abe, S. Mori, T. Nakajima, and K. Hirao, Chem. Phys. 311, 129 (2005). 7 P. Salek, T. Helgaker, and T. Saue, Chem. Phys. 311, 187 (2005). 8 S. Komorovský, M. Repsiský, O. L. Malkina, V. G. Malkin, I. Malkinaí, and M. Kaupp, J. Chem. Phys. 128, 104101 (2008). 9 W. H. E. Schwarz and E. Wechsei-Trakowski, Chem. Phys. Lett. 85, 94 (1982). 10 T. Nakajima, T. Yanai, and K. Hirao, J. Comput. Chem. 23, 847 (2002). 11 W. Liu and D. Peng, J. Chem. Phys. 125, 044102 (2006). 12 M. S. Kelley and T. Shiozaki, J. Chem. Phys. 138, 204113 (2013). 13 R. Samzow, B. A. Hess, and G. Jansen, J. Chem. Phys. 96, 1227 (1992). 14 T. Nakajima and K. Hirao, J. Chem. Phys. 119, 4105 (2003). 15 C. van Wüllen and C. Michauk, J. Chem. Phys. 123, 204113 (2005). 16 J. Seino and M. Hada, Chem. Phys. Lett. 442, 134 (2007). 17 J. Seino and M. Hada, Chem. Phys. Lett. 461, 327 (2008). 18 J. Sikkema, L. Visscher, T. Saue, and M. lliaš, J. Chem. Phys. 131, 124116 (2009). 19 C. Chang, M. Pelissier, and P. Durand, Phys. Scr. 34, 394 (1986). 20 E. van Lenthe, E. J. Baerends, and J. G. Snijders, J. Chem. Phys. 99, 4597 (1993). 21 E. van Lenthe, E. J. Baerends, and J. G. Snijders, J. Chem. Phys. 101, 9783 (1994). 22 T. Nakajima and K. Hirao, Chem. Phys. Lett. 302, 383 (1999). 23 K. G. Dyall, J. Chem. Phys. 106, 9618 (1997). 24 E. van Lenthe, A. Ehlers, and E.-J. Baerends, J. Chem. Phys. 110, 8943 (1999). 25 J. H. van Lenthe, S. Faas, and J. G. Snijders, Chem. Phys. Lett. 328, 107 (2000). 26 M. Filatov and D. Cremer, J. Chem. Phys. 118, 6741 (2003). 27 D. G. Fedorov, T. Nakajima, and K. Hirao, Chem. Phys. Lett. 335, 183 (2001). 28 M. Filatov and D. Cremer, Chem. Phys. Lett. 370, 647 (2003). 29 W. Zou, M. Filatov, and D. Cremer, J. Chem. Phys. 134, 244117 (2011). 30 W. Kutzelnigg and W. Liu, J. Chem. Phys. 123, 241102 (2005). 31 L. L. Foldy and S. A. Wouthuysen, Phys. Rev. 78, 29 (1950). 32 M. Douglas and N. Kroll, Ann. Phys. 82, 89 (1974). 33 T. Nakajima and K. Hirao, J. Chem. Phys. 113, 7786 (2000). 34 M. Barysz, J. Chem. Phys. 114, 9315 (2001). 35 W. A. de Jong, R. J. Harrison, and D. A. Dixon, J. Chem. Phys. 114, 48 (2001). 36 M. Reiher and A. Wolf, J. Chem. Phys. 121, 10945 (2004). 37 M. Barysz and A. J. Sadlej, J. Chem. Phys. 116, 2696 (2002). 38 V. A. Nasluzov and N. Rösch, Chem. Phys. 210, 413 (1996). 39 I. Malkina, O. L. Malkina, and V. G. Malkin, Chem. Phys. Lett. 361, 231 (2002). 40 A. V. Matveev, V. A. Nasluzov, and N. Rösch, Int. J. Quantum Chem. 107, 3236 (2007). 41 F. Neese, A. Wolf, T. Fleig, M. Reiher, and B. A. Hess, J. Chem. Phys. 122, 204107 (2005). 42 R. Mastalerz, G. Barone, R. Lindh, and M. Reiher, J. Chem. Phys. 127, 074105 (2007). 43 J. Seino and H. Nakai, J. Chem. Phys. 136, 244102 (2012). 44 J. Seino and H. Nakai, J. Chem. Phys. 137, 144101 (2012). 45 W. Yang, Phys. Rev. Lett. 66, 1438 (1991). 46 W. Yang and T.-S. Lee, J. Chem. Phys. 103, 5674 (1995). 47 T. Akama, M. Kobayashi, and H. Nakai, J. Comput. Chem. 28, 2003 (2007). 48 M. Kobayashi, T. Akama, and H. Nakai, J. Chem. Phys. 125, 204106 (2006). 49 M. Kobayashi, Y. Imamura, and H. Nakai, J. Chem. Phys. 127, 074103 (2007). 50 M. Kobayashi and H. Nakai, J. Chem. Phys. 129, 044103 (2008). 51 M. Kobayashi and H. Nakai, Int. J. Quantum Chem. 109, 2227 (2009). 52 M. Kobayashi and H. Nakai, J. Chem. Phys. 131, 114108 (2009). 53 M. Kobayashi and H. Nakai, Phys. Chem. Chem. Phys. 14, 7629 (2012). 54 J. Seino and H. Nakai, J. Chem. Phys. 139, 034109 (2013). 55 H. A. Bethe and E. E. Salpeter, Quantum Mechanics of One- and TwoElectron Atoms (Springer, Berlin, 1957).

r

2 K.

FIG. 4. Difference in bond length calculated by the RECP, MCP, and LUTIODKH/C from that obtained using the NR/C method.

IV. CONCLUSION

In this study, we have developed the analytical energy gradient for the IODKH method and proposed the efficient scheme by adopting the LUT technique. The matrix evaluation technique proposed by Hess was applied to analytical energy gradient. The accuracy of the present method was numerically confirmed for the diatomic molecules, halogen halides HX, halogen dimers X2 (X = F, Cl, Br, I, and At), coinage metal halides MH, and coinage metal dimers M2 (M = Cu, Ag, and Au). The bond lengths in diatomic molecules calculated by IODKH method agree very well with the results obtained by using the 4c method. Especially for coinage dimers, the spin-free effect directly affect bond lengths. In addition, the efficiency in the calculation of analytical energy gradient for relativistic transformation was examined with the comparison of the CPU time using the one, two-, and three-dimensional silver clusters. Applying the LUT scheme to the gradient resulted in a significant reduction in the computational cost for unitary transformation, while the accuracy was the same as that achieved using the conventional IODKH method. The present analytical energy gradient method for LUT-IODKH was found to be applicable to more practical systems, namely, metal complexes.

ACKNOWLEDGMENTS

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), the Ministry of Education, Culture, Sports, Science and Technology (MEXT), and the Computational Materials Science Initiative (CMSI), Japan; by 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 the Japan Society for the Promotion of Science (JSPS) for a Research Fellowship.

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: 134.99.128.41 On: Tue, 31 Dec 2013 15:59:03

244107-13

Nakajima, Seino, and Nakai

J. Chem. Phys. 139, 244107 (2013)

56 B.

75 J.

57 L.

76 D.

A. Hess, Phys. Rev. A 32, 756 (1985). Visscher, Theor. Chem. Acc. 98, 68 (1997). 58 I.-K. Suh, H. Ohta, and Y. Waseda, J. Mater. Sci. 23, 757 (1988). 59 A. D. Becke, J. Chem. Phys. 98, 5648 (1993). 60 Y. Sakai, E. Miyoshi, M. Klobukowski, and S. Huzinaga, J. Chem. Phys. 106, 8084 (1997). 61 T. Noro, M. Sekiya, and T. Koga, Theor. Chem. Acc. 98, 25 (1997). 62 M. Sekiya, T. Noro, T. Koga, and H. Matsuyam, J. Mol. Struct.: THEOCHEM 451, 51 (1998). 63 E. Miyoshi, Y. Sakai, K. Tanaka, and M. Masamura, J. Mol. Struct.: THEOCHEM 451, 73 (1998). 64 T. Noro, M. Sekiya, T. Koga, and H. Matsuyama, Theor. Chem. Acc. 104, 146 (2000). 65 M. Sekiya, T. Noro, Y. Osanai, and T. Koga, Theor. Chem. Acc. 106, 297 (2001). 66 Y. Osanai, M. Sekiya, T. Noro, and T. Koga, Mol. Phys. 101, 65 (2003). 67 T. Noro, M. Sekiya, Y. Osanai, E. Miyoshi, and T. Koga, J. Chem. Phys. 119, 5142 (2003). 68 Y. Osanai, T. Noro, E. Miyoshi, M. Sekiya, and T. Koga, J. Chem. Phys. 120, 6408 (2004). 69 Y. Osanai, M. S. Mon, T. Noro, H. Mori, H. Nakashima, M. Klobukowski, and E. Miyoshi, Chem. Phys. Lett. 452, 210 (2008). 70 H. Mori, K. Ueno-Noto, Y. Osanai, T. Noro, T. Fujiwara, M. Klobukowski, and E. Miyoshi, Chem. Phys. Lett. 476, 317 (2009). 71 M. Dolg, U. Wedig, H. Stoll, and H. Preuss, J. Chem. Phys. 86, 866 (1987). 72 D. Figgen, G. Rauhut, M. Dolg, and H. Stoll, Chem. Phys. 311, 227 (2005). 73 K. A. Peterson, D. Figgen, M. Dolg, and H. Stoll, J. Chem. Phys. 126, 124101 (2007). 74 D. Figgen, K. A. Peterson, M. Dolg, and H. Stoll, J. Chem. Phys. 130, 164108 (2009).

T. H. Dunning, J. Chem. Phys. 90, 1007 (1989). E. Woon and J. T. H. Dunning, J. Chem. Phys. 98, 1358 (1993). 77 M. W. Schmidt, J. A. B. Kim, K. Baldride, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis, and J. A. Montgomery, J. Comput. Chem. 14, 1347 (1993). 78 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š, Ch. R. Jacob, S. Knecht, S. Komorovský, 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. Sałek, 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. 79 M. J. Frisch, G. W. Trucks, H. B. Schlegel, et al., Gaussian 09, Revision A.2, Gaussian, Inc., Wallingford, CT, 2009. 80 Y.-K. Han, C. Bae, S.-K. Son, and Y. S. Lee, J. Chem. Phys. 112, 2684 (2000). 81 J. Y. Seto, Z. Morbi, F. Charron, S. K. Lee, P. F. Bernath, and R. J. L. Roy, J. Chem. Phys. 110, 11756 (1999). 82 K. P. Huber and G. Herzberg, Chemistry Spectra and Molecular Structure IV. Constants of Diatomic Molecules (Van Nostrand, New York, 1979). 83 T. Noro, M. Sekiya, and T. Koga, Theor. Chem. Acc. 109, 85 (2003). 84 T. Koga, H. Tatewaki, and T. Shimazaki, Chem. Phys. Lett. 328, 473 (2000). 85 R. S. Ram, C. N. Jarman, and P. F. Bernath, J. Mol. Spectrosc. 156, 468 (1992). 86 H.-G. Krämer, V. Beutel, K. Weyers, and W. Demtröder, Chem. Phys. Lett. 193, 331 (1992).

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: 134.99.128.41 On: Tue, 31 Dec 2013 15:59:03