Importance of the hybrid orbital operator derivative term for the energy gradient in the fragment molecular orbital method

Chemical Physics Letters 492 (2010) 302–308

Contents lists available at ScienceDirect

Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett

Importance of the hybrid orbital operator derivative term for the energy gradient in the fragment molecular orbital method Takeshi Nagata a,*, Dmitri G. Fedorov a, Kazuo Kitaura a,b a b

RICS, National Institute of Advanced Industrial Science and Technology (AIST), 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan Graduate School of Pharmaceutical Sciences, Kyoto University, 46-29 Yoshidashimo Adachi, Sakyo-ku, Kyoto 606-8501, Japan

a r t i c l e

i n f o

Article history: Received 24 March 2010 In final form 16 April 2010 Available online 21 April 2010

a b s t r a c t The hybrid orbital operator is crucial in the fragment molecular orbital (FMO) method for the fragmentation across covalent bonds, however, its gradients have not been properly derived. We show that these very substantial contributions are the cause of the major part of the gradient error in FMO, impeding geometry optimizations and molecular dynamics. Capped alanine decamer ðALAÞ10 and chignolin (PDB: 1UAO) solvated by 157 water molecules are used to assess the accuracy of the energy gradients, and the errors are reduced by approximately one order of magnitude. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction Many computational methods have been developed recently in order to facilitate large scale quantum-mechanical calculations [1– 10], which use fragmentation to reduce the steep scaling of ab initio methods. The fragment molecular orbital (FMO) method [11–14] has been shown to accurately reproduce ab initio properties such as the total energy and its gradient for large systems by performing molecular orbital (MO) calculations for fragments and their pairs under the embedding electrostatic potentials (ESP) due to the densities of the remaining fragments. FMO has been extended to various wave functions: second order Møller– Plesset (MP2) perturbation theory [15,16], coupled-cluster (CC) theory [17], density functional theory (DFT) [18,19], multiconfiguration self-consistent field (MCSCF) [20], configuration interaction (CI) [21,22], time-dependent DFT (TDDFT) [23–25], nuclear-electronic orbitals (NEO) [26] and open shells [27]. The accuracy of the total energy has been established in comparison to conventional MO methods at various levels of theory. The energy gradient, however, is known to be not fully analytic [28,29], in particular, the response term in the FMO energy gradients is one of the important missing terms in the FMO energy gradients, that is neglected at present because the time consuming coupled-perturbed Hartree–Fock (CPHF) equations are hard to solve for large molecules. In this study, we present a very important term missing in the energy gradient, which is the derivative of the hybrid orbital projection (HOP) term. This term is introduced to project out the * Corresponding author. Fax: +81 29 851 5426. E-mail address: [email protected] (T. Nagata). 0009-2614/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2010.04.043

redundant orbitals in fragment MOs around the detached covalent bond, which is important for FMO. In the original formulation of the FMO energy gradient [30], the HOP derivative contribution has been thought to be very small because the HOP term itself has a very small value for a converged wave function, which has been numerically verified only for the minimal basis set; we, however, show that it is responsible for a major portion of the gradient error in FMO for fragments connected by covalent bonds (for molecular clusters the HOP term does not arise). A considerable deviation of the FMO gradient from being fully analytic has caused unnecessary long geometry optimizations [31], and lead to the loss of the energy conservation in molecular dynamics [28,32]; because of this problem, the latter has been applied mainly to molecular clusters and some initial attempts to apply it to polypeptides have had to use large fragments to reduce the gradient error, leading to a very substantial increase in the simulation costs [28]. One can also expect that using gradients for further numeric derivatives (e.g., in the Hessian calculations) is also bound to produce large errors. The meaning of the HOP term (see [33] for details) is as follows. A detached covalent bond on the fragment border is formed by a bond detached atom (BDA) and a bond attached atom (BAA), divided between two fragments. In order to describe the detached bond assigned to the fragment which has the BAA, BDA from the other fragment is included in it in the ghost fashion. Thus, BDA is present in two fragments, and it is necessary to divide its atomic orbitals properly, which is accomplished by HOP. It is constructed with hybridized orbitals (see below), and pushes the orbitals to be projected out very high in energy, so that the variational space has a proper division of the BDA orbitals. The hybrid orbitals are constructed for a given type of hybridization in advance, e.g., for

303

T. Nagata et al. / Chemical Physics Letters 492 (2010) 302–308

a protein usually only single C–C bonds are detached, so one does a single methane calculation to obtain sp3 orbitals of C, which are used for all carbons serving as BDA. In our previous study [29], we compared the analytic and the numeric energy gradients, using double and triple-f basis sets with polarization, and it was found that the largest errors arise for the BDA and BAA coordinates. We show in this Letter that the density derivative term in the HOP derivative is already included in the existing FMO gradients implicitly, but other terms are neglected, which is the reason for the large errors. We introduce the other terms in the HOP derivative for proper balance and show that the accuracy is much improved. We demonstrate the efficiency of our method by optimizing chignolin (PDB: 1UAO) in the presence of solvent.

2. Theory 2.1. FMO energy gradient We briefly summarize the main FMO equations. The two-body FMO (FMO2) total energy expression is

E¼

N X

E0I þ

I

N N X X ðE0IJ  E0I  E0J Þ þ TrðDDIJ VIJ Þ; I>J

h

X X X 1 X  X @DXlm X 1 X @ Dlm Dkr  2 Dlk Dmr  vl vm jvk vr hl m þ 2 lmkr2X @a @a lm2X

þ

X @DXlm @a

lm2X

PXlm ¼ 4

X

X

X

F 0X mi ¼ hmi þ

lm2X

þ

X

X

X

Dlm Plm þ

and U a;X mi are unknowns arising from the derivative of the MO coefficient C Xli , found by solving CPHF equations [34]. PXmi is the MObased matrix element of PXlm (see below). The derivative of the ESP energy in Eq. (3) is discussed in [29]. The Fock term associated with U a;X mi is transformed using the X X relation F Xmi ¼ 0 ðm–iÞ, where F Xmi ¼ F 0X mi þ V mi and V mi are MObased matrix elements of the embedding potential [29]. With these relations, Eq. (5) becomes

4

occ occþvir X X

0X U a;X mi F mi ¼ 4

m2X

i2X

i2X

ð3Þ

lm2X

The differentiation of the internal energy E0X with respect to nuclear coordinate a leads to

X X @hXlm @E0X Dlm ¼ @a @a lm2X     @ v v jv v l m k r 1 X 1 DX DX  DX DX þ 2 lmkr2X lm kr 2 lk mr @a X

DXlm

@PXlm

þ

X @DXlm

occ occþvir X X

occ X

i2X 0X Sa;X ij F ij  4

ij2X

Sa;X ij ¼

X U a;X mi V mi

m2X

occ X vir X

X U a;X mi V mi ;

ð7Þ

i2X m2X

X

ð8Þ

X

C X li

@Slm @a

C Xmj :

ð9Þ

To derive Eq. (7), we used the relation in Eq. (8) arising from the orthogonality of MOs, SXij ¼ dij [34]. The last term in Eq. (7) cancels out when ESPs are either computed exactly or completely replaced by point charges; otherwise the cancellation holds approximately (see Eq. (15) of Ref. [29]). 0X At present, only the Sa;X ij F ij term (Eq. (7)) implicitly includes a HOP contribution term PXmi in Eq. (6), but another term P X X lm2X Dlm @P lm =@a in Eq. (4) is simply neglected. This imbalance is found to cause large errors as shown below. 2.2. Derivatives of the HOP term

bX ¼ P

X

Bk jhk ihhk j;

ð10Þ

k2X

and its contribution to the energy of X in Eq. (2) is X

hlm

@a lm2X @a h i X X X 1 X  1 X @ Dlm Dkr  2 Dlk Dmr  þ vl vm jvk vr 2 lmkr2X @a X @DXlm X @ENR P lm þ X : @a lm2X @a

X Sa;X ii F ii  4

b X is given by The HOP operator P

PX ¼

lm2X

þ

¼ 2

lm2X

where D and h are the density and one-electron Hamiltonian matrices for X, respectively. v are atomic orbitals, numbered by l; m; k and r. The HOP matrix elements PXlm are defined later. ENR X is the nuclear repulsion energy. The contribution due to the embedding electrostatic potential V [12] is given by

þ

m2X

a;X a;X U a;X ij þ U ji ¼ Sij

ð2Þ

h i  X DIJ  ðDI  DJ Þ VIJ ¼ DDIJlm V IJlm :

occ X

  X X U a;X mi F mi  V mi

i2X

X

TrðDDIJ VIJ Þ ¼ Tr

occ occþvir X X

where

lm2X X

ð6Þ

j2X

¼ 2

ENR X ;

ð5Þ

occ h    i X þ PXmi ; 2 /Xm /Xi j/Xj /Xj  /Xm /Xj j/Xj /Xi

ð1Þ

  1 X 1 DXlm DXkr  DXlk DXmr ðvl vm jvk vr Þ 2 lmkr2X 2

0X U a;X mi F mi ;

m2X

where F 0X mi is the internal Fock matrix element in the MO basis and m and i number MOs:

I>J

DXlm hlm þ

occ occþvir X X i2X

where X ¼ I and IJ for fragment monomers and dimers (pairs), respectively. The internal energy E0X is

E0X ¼

i

occ D  E X X bX  2 i P DXlm P Xlm ; i ¼ i2X

ð4Þ

The atomic orbital (AO) derivative terms and the Hellmann– Feynman term on the right-hand side of Eq. (4) can be calculated straightforwardly. The set of the density derivative terms is transformed into

ð11Þ

lm2X

b X jmi. where PXlm ¼ hlj P X b P acts on the MOs and projects out the orbitals, which are described in another fragment; Bk is set to a large constant value 106 a:u: (independent of k at present). k runs over all orbitals to be projected out of the variational space of X. For a typical fragment in a protein, there are two detached bonds, and k runs over hybrid orbitals on the BDAs of these two bonds. The hybrid orbitals jhk i are expanded over atomic orbitals of X; jqi, centered on the BDA for orbital k

304

PX ¼

T. Nagata et al. / Chemical Physics Letters 492 (2010) 302–308 occ X X X i2X

¼

occ occþvir X X X @SXlq X X @PX X e X SX D 4U a;X P þ Dlm ¼ mi mi @a @a qr rm m2X lm2X qr i2X

e X hrjmiC X 2C Xli hljqi D qr mi

lm2X qr

X X

e X SX : DXlm SXlq D qr rm

ð12Þ

lm2X qr

þ

X X

eX DXlm SXlq D qr

lm2X qr

e X is formally similar to the weighted density of The definition of D qr the hybrid orbitals in X (they are not occupied though), where q and r run over all AOs in X for convenience, and it is a sparse matrix with non-zero elements for AOs on BDAs, see Appendix for the defeX . inition of D qr In order to understand what contributions arise from the projection operator, we differentiate the HOP term in Eq. (12) with respect to a

eX @SXrm X X X X @ D qr X Dlm Slq þ Srm ; @a @a lm2X qr

ð13Þ

and by combining the second and third terms, X occ occþvir XX X X X X @PX X e X @Srm 4U a;X Dlm Slq D ¼ qr mi P mi þ 2 @a @a lm qr m2X i2X

þ

XX lm

qr

DXlm SXlq

eX @D qr @a

SXrm :

ð14Þ

Fig. 1. (a) a-ðALAÞ10 capped with acetyl and -NHCH3 groups. (b) Hydrated chignolin optimized at the FMO-RHF/EFP level with 6-31G (d) (colored by chemical elements as light grey (H), dark grey (C), blue (N), and red (O)). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this Letter.)

305

T. Nagata et al. / Chemical Physics Letters 492 (2010) 302–308

The matrix elements of the projection operator PXmi for the MOs optimized in SCF have very small values [35], and the total gradient in Eq. (14) is small as well, but each term in Eq. (14) is not. The FMO gradient so far included only the first term in Eq. (14) implicitly, as a part of the Fock matrix in Eq. (6). There are two different ways to alleviate this imbalance from the energy gradients: one is to remove the response term (the first term in Eq. (14)) from the terms involving the Fock matrix, i.e., subtract P Xmi from F Xmi in Eq. (7) and do not add any other correction from Eq. (14), assuming that @PX =@a is approximately zero; and another is to implement the second and third terms on the righthand side of Eq. (14) while retaining the Fock matrix in Eq. (7) as it is shown. Among the two approaches, the former is an approximation leading to the modification of the terms with U a;I mi in Eq. (7) arising from the differentiation of the electrostatic potential, and we take the latter instead and consider all terms in Eq. (14), which is exact. The second term in Eq. (14) is trivial to implement using the derivative overlap integrals readily available. The last term in Eq. (14) includes the derivatives of the expansion coefficients of hybrid orbitals. In practice, these coefficients are generated for the necessary hybridization state of BDA in a prototypical molecule (e.g., methane), and for this purpose the corresponding bond (i.e., any of the C–H bonds in methane) is oriented parallel to z-axis. However, in a system to which FMO is applied, the bond orientation can be arbitrary, and thus the hybrid orbital expansion coefficients precomputed for the z-axis have to be rotated for the actual bond direction. These rotated coefficients depend upon the atomic coordinates of BDA and BAA, which determine the bond direction. The derivative of this last term is given in Appendix.

3. Computational details To assess the accuracy of the FMO energy gradients and verify the importance of the HOP derivatives, we chose two molecular systems from the previous study [29] for test calculations: the a-helix of alanine decamer ðALAÞ10 capped with acetyl and -NHCH3 groups and chignolin (PDB ID: 1UAO) [36] solvated by 157 water molecules using the effective fragment potentials (EFP) to describe water [37]. The structures of ðALAÞ10 and hydrated chignolin (Fig. 1) were taken from the previous study [29], and these systems were fragmented as one residue per fragment, as shown in Fig. 2. We applied the electrostatic potential point charge (ESP-PTC) and the electrostatic dimer (ES-DIM) approximations [12] for the ESPs of separate fragments. The ESP-PTC threshold was set to 2.5; for ES-DIM we used 2.5 and 2.0 for ðALAÞ10 and hydrated chignolin, respectively (the thresholds are unitless, expressed in terms of the van-der-Waals atomic radii [12]). We used the FMO gradient code [29–31] modified with the addition of the HOP derivatives. The development and all calculations were done using GAMESS [38,39] package parallelized by the generalized distributed data interface [40]. Spherical functions are used throughout and other accuracy thresholds were set to default values. Numeric gradients were computed with double differencing, with the coordinate step of 0.005 Å. 4. Results and discussion Fig. 1a depicts the a-helix structure of ðALAÞ10 , which is expected to produce large gradient errors, because of its compact form compared to the b-sheet and the extended isomers. Note that for the HOP derivatives, the second term in Eq. (14) has values over

Fig. 2. Fragmentation of (a) ðALAÞ10 and (b) chignolin (fragment residues in FMO are shifted from conventional ones by the –CO group).

Table 1 Errors of the analytic gradients with and without the HOP derivative terms relative to the corresponding numeric gradients for CH3 COðALAÞ10 NHCH3 at the (a) RHF/6-31G (d) and (b) RHF/6-311G (d) levels, and their RMSDs (all in 105 a:u:). (a) 6-31G (d)

(b) 6-311G (d)

With HOP deriv.

C8 (BDA) C9 (BAA) C18 (BDA) C19 (BAA) C28 (BDA) C29 (BAA) C38 (BDA) C39 (BAA) C48 (BDA) C49 (BAA) C58 (BDA) C59 (BAA) C68 (BDA) C69 (BAA) C78 (BDA) C79 (BAA) C88 (BDA) C89 (BAA) RMSD

Without HOP deriv.

With HOP deriv.

Without HOP deriv.

x

y

z

x

y

z

x

y

z

x

y

z

23 23 38 23 49 7 30 2 51 60 24 23 30 2 2 22 60 32 33

16 61 27 45 41 18 2 20 11 46 19 12 62 38 42 8 69 38

10 33 20 41 30 2 37 27 5 13 3 15 7 7 55 8 19 71

181 237 61 28 330 205 131 14 156 118 159 8 50 10 162 9 124 132 147

154 67 172 16 105 110 99 162 36 67 104 146 191 231 180 144 304 257

173 93 188 227 40 43 207 135 138 20 173 139 87 82 44 88 217 37

112 35 11 26 183 29 9 2 179 108 70 27 84 45 70 10 28 32 69

5 55 79 33 28 19 32 22 46 34 0 40 139 38 74 59 92 53

35 47 144 3 0 19 155 47 60 16 119 63 73 15 9 6 32 71

611 295 407 98 875 128 615 218 406 41 738 117 166 266 742 112 629 446 418

664 10 659 77 107 355 59 263 518 193 217 403 118 428 50 133 548 505

658 252 547 292 539 57 645 37 524 243 256 41 799 39 196 357 634 160

306

T. Nagata et al. / Chemical Physics Letters 492 (2010) 302–308

all atoms in X, while the third term has non-zero values only with respect to BDAs and BAAs. Nevertheless, the HOP derivatives mainly contribute to the BDA and BAA gradient elements, and we thus focus on comparing the analytic and numeric gradient elements for them in the following discussion. For ðALAÞ10 , we calculated the errors of the analytic gradients with and without the HOP derivatives relative to the numeric gradients at the RHF/6-31G (d) level. We denote the former and the latter as ‘new gradients’ and ‘old gradients’ in this Letter. Table 1a lists the gradient errors at the RHF/6-31G (d) level, where the two consecutive carbon atoms, e.g., C8 and C9, form detached covalent bonds. Fig. 3a plots all values shown in Table 1a for 54 elements (corresponding to having 9 detached bonds, two atoms in each, and three coordinates). One can readily see that in Table 1a the errors are improved by about an order of magnitude and as seen in Fig. 3a, the zigzag shape (square) of the gradient errors is flattened to small values (diamond) by introducing the HOP derivatives. We also calculated the root mean square deviations (RMSD) for both sets of the gradient errors, which are shown in the lowest row of Table 1a. The new gradients are found to be significantly improved by a factor of 4, but still have the errors on the order of 104 a.u., attributed to the CPHF-related terms [29].

It is important to assess how an increase in the basis set quality [41] affects the accuracy of gradients. Therefore, we used 6-311G (d) and repeated the accuracy tests, whose results are shown in Table 1b and Fig. 3b. As expected, the gradient errors in Table 1b and the amplitude of the zigzag shape (square) in Fig. 3b are much larger for the larger basis set. In comparison to the old values, the new gradients are improved by a factor of 6. Very large errors such as 8:75  103 a:u: are considerably reduced (to 1:83  103 a:u:). The gradient errors for the remaining atoms are much smaller, thus we expect that one can now optimize structures with FMO using large basis sets to the typical threshold of 104 a:u:, which was hardly possible without the HOP derivatives. The residual error in the gradients of the BDA (arising from the neglected terms) shows some increase with the basis set. Chignolin is a small protein consisting of 10 amino acids [36], which has the b-hairpin structure as shown in Fig. 1b. Chignolin was solvated by 157 water molecules treated with EFP (a water layer of 5 Å) and optimized at the FMO-RHF/6-31G (d) level. Table 2 shows that for the initial geometry taken from the earlier study [29] (it was optimized at the same level but without HOP derivatives) Table 2 Errors of the analytic gradients with and without the HOP derivative terms relative to the corresponding numeric gradients for hydrated chignolin at the RHF/6-31G (d) level, and their RMSDs (all in 105 a:u:). 6-31G (d)

C2 (BDA) C3 (BAA) C6 (BDA) C7 (BAA) C18 (BDA) C19 (BAA) C26 (BDA) C27 (BAA) C33 (BDA) C34 (BAA) C42 (BDA) C43 (BAA) C49 (BDA) C50 (BAA) C53 (BDA) C54 (BAA) C60 (BDA) C61 (BAA) RMSD

Fig. 3. Errors of the analytic gradient elements relative to the numeric gradient elements, plotted for BAAs and BDAs only. The system is ðALAÞ10 at (a) the FMORHF/6-31G (d) level and (b) the FMO-RHF/6-311G (d) level. Diamond: analytic gradient with the HOP derivatives. Square: analytic gradient without the HOP derivatives.

With HOP deriv.

Without HOP deriv.

x

y

z

x

y

z

4 22 53 50 34 104 37 3 10 15 20 57 19 24 25 17 1 1 39

0 4 22 3 74 29 43 10 11 8 16 0 5 129 7 71 31 10

4 9 31 17 6 87 7 9 22 4 39 83 22 25 1 71 39 25

88 179 131 77 169 159 18 149 382 184 101 204 418 102 58 185 102 34 176

27 25 9 58 94 8 156 97 115 12 25 276 180 120 43 157 126 73

20 45 361 128 211 1 235 125 82 111 278 102 623 178 117 101 127 12

Fig. 4. Errors of the analytic gradient elements relative to the numeric gradient elements, plotted for BAAs and BDAs only. The system is hydrated chignolin at the FMO-RHF/EFP/6-31G (d) level. Diamond: analytic gradient with the HOP derivative. Square: analytic gradient without the HOP derivative.

307

T. Nagata et al. / Chemical Physics Letters 492 (2010) 302–308

most of the gradient errors are below 104 a:u: and the RMSD is comparable to that for ðALAÞ10 with the same basis set; the errors are not increased by the presence of water molecules. In Fig. 4, one can see that the accuracy of the new gradients is obviously improved relative to the old values. Next we compare the results of the two geometry optimizations for the hydrated chignolin, with and without the HOP derivatives. We note that we reoptimized both structures starting from the same initial structure to make sure that the same level of approximations and accuracy thresholds is used [37,42]. The minimum energies obtained with the old and new gradients are 3802.827646 a.u. and 3802.839738 a.u., respectively, and the latter is more stabilized by 31.7 kJ/mol. The new gradients are expected to refine the structures of proteins with a little extra computational cost.

eX ¼ D lm

X

Bk ~cllk ~clmk ¼

k

X

Bk

X

l

l

dlq dmr clqk clrk :

ð17Þ

q;r

k

l

Consequently, one has to take derivative of d because clmk and Bk are l constant in FMO calculation. Thus, one needs @dlq =@a, where a is a coordinate of either BDA or BAA. Let us define b as the vector pointing from BDA to BAA, b ¼ RBAA  RBDA and z as the z-axis unit vector (0, 0, 1). The counterclockwise rotation matrix A along the unit vector n by angle x is given by [45]:

0

nx nx nx ny A ¼ cos x1 þ ð1  cos xÞ@ ny nx ny ny nz nx nz ny 0 1 0 nz ny 0 nx A; þ sin x@ nz ny nx 0

1 nx nz ny nz A nz nz ð18Þ

where 5. Conclusions

zb ; jz  bj zb : cos x ¼ b

n¼

We have shown that the partial treatment of the hybrid orbital projection derivatives causes large errors in the FMO energy gradients. By introducing the missing terms, we have improved the gradients by about one order of magnitude, especially for the atoms forming detached bonds, as shown for the typical polypeptide molecules, ðALAÞ10 and hydrated chignolin. It has been found that the optimized structures have a considerably lower energy, while calculations of additional terms require a negligible time. We conclude that HOP derivatives should always be added to the FMO energy gradient calculation. The technique of the differentiation of the rotated molecular orbitals may be useful in some other method development. The derived and implemented HOP derivative terms are expected to improve convergence of geometry optimizations and molecular dynamics[32,43,44], although the latter may need the CPHF-related terms which we are going to implement in future. Acknowledgment This work has been supported by the Next Generation Super Computing Project, Nanoscience Program (MEXT, Japan). Appendix A

l

The matrices d are generated from the rotation matrix A, i.e., l l 0 1 d ðAÞ ¼ d ðn; xÞ. For example, d ¼ 1 for s-orbitals and d ¼ Ay for 2 p-orbitals. For d-orbitals, d of 6  6 (if Cartesian functions are used) is obtained from the explicit transformation of xx, xy, etc. according 1 to that of x, y, z in d . For complex-valued spherical functions the matrices d are irreducible representations of the rotation operators in the SOð3Þ group in the basis of spherical harmonics, known as the Wigner functions, and a closed compact form of them can be used [45]. l For Cartesian products the matrices d form reducible representations for l > 1. These matrices d are used for various point-group related operations in quantum-chemistry codes, for example, for the assignment of irreducible representation labels to CI states in non-Abelian groups [46]. Taking advantage of the explicit form of z ¼ ð0; 0; 1Þ and the a-dependent vector b, one can simplify these equations further. l The necessary derivative of @dlq =@a are thus reduced to the derivatives of the matrix elements of A or their products, all of which can be derived analytically and which we omit due to their large volume and the trivial work in derivations. To give one specific example, 1

d11 ¼

Hybrid orbitals are expanded in terms of AOs as

jhk i ¼

XX flg

vll ~cllk ;

~l

clk ¼

X

l2l

l

l

dlm cmk ;

ð16Þ

m2l

where m 2 l restricts m to the same l-block as l, such as l and m run over x; y; z for l ¼ 1. The last term in Eq. (14) requires taking deriveX ative of D qr

RBAA  RBDA z z b RBAA  RBDA z þ 1 z b

ð15Þ

where index l denotes a block of coefficients (e.g., a p-block has three coefficients l ¼ x; y; z) and the sum flg runs over all such blocks. We show explicit l values for convenience of further derivations, where l is either the angular momentum if spherical atomic orbitals are used, or the rank of the Cartesian product, i.e., l ¼ i þ j þ k, for a product xi yj zk . vl are atomic basis functions on the BDA (for orbital k). The coefficients ~cllk for the specific bond l orientation are expanded in terms of the transformation matrix d describing the rotation of the fixed coefficients cllk of the precalculated hybrid orbitals for the detached bond aligned along z-axis to the actual bond orientation:

ð19Þ

!

ðRBAA  RBDA Þ2 y y  2  2 ; RBAA  RBDA þ RBAA  RBDA x x y y ð20Þ

which can now be differentiated with respect to a, one of the elements of RBAA or RBDA ðb ¼ jRBAA  RBDA jÞ. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]

A.P. Rahalkar, V. Ganesh, S.R. Gadre, J. Chem. Phys. 129 (2008) 234101. W. Hua, T. Fang, W. Li, J.-G. Yu, S. Li, J. Phys. Chem. A 112 (2008) 10864. P. Söderhjelm, U. Ryde, J. Phys. Chem. A 113 (2009) 617. H.R. Leverentz, D.G. Truhlar, J. Chem. Theory Comput. 5 (2009) 1573. M.S. Gordon, J.M. Mullin, S.R. Pruitt, L.B. Roskop, L.V. Slipchenko, J.A. Boatz, J. Phys. Chem. B 113 (2009) 9646. R.A. Mata, H. Stoll, B.J.C. Cabral, J. Chem. Theory Comput. 5 (2009) 1829. W. Xie, M. Orozco, D.G. Truhlar, J. Gao, J. Chem. Theory. Comput. 5 (2009) 459. A. Pomogaeva, F.L. Gu, A. Imamura, Y. Aoki, Theor. Chem. Acc. 125 (2010) 453. T. Touma, M. Kobayashi, H. Nakai, Chem. Phys. Lett. 485 (2010) 247. X. He, K.M. Merz, J. Chem. Theory Comput. 6 (2010) 405. K. Kitaura, E. Ikeo, T. Asada, T. Nakano, M. Uebayasi, Chem. Phys. Lett. 313 (3,4) (1999) 701. T. Nakano, T. Kaminuma, T. Sato, K. Fukuzawa, Y. Akiyama, M. Uebayasi, K. Kitaura, Chem. Phys. Lett. 351 (5,6) (2002) 475.

308

T. Nagata et al. / Chemical Physics Letters 492 (2010) 302–308

[13] D.G. Fedorov, K. Kitaura, J. Phys. Chem. C 111 (30) (2007) 6904. [14] D.G. Fedorov, K. Kitaura (Eds.), The Fragment Molecular Orbital Method: Practical Applications to Large Molecular Systems, CRC Press, Boca Raton, FL, 2009. [15] D.G. Fedorov, K. Kitaura, J. Chem. Phys. 121 (6) (2004) 2483. [16] Y. Mochizuki, T. Nakano, S. Koikegami, S. Tanimori, Y. Abe, U. Nagashima, K. Kitaura, Theor. Chem. Acc. 112 (2004) 442. [17] D.G. Fedorov, K. Kitaura, J. Chem. Phys. 123 (13) (2005) 134103/1. [18] S.-I. Sugiki, N. Kurita, Y. Sengoku, H. Sekino, Chem. Phys. Lett. 382 (2003) 611. [19] D.G. Fedorov, K. Kitaura, Chem. Phys. Lett. 389 (2004) 129. [20] D.G. Fedorov, K. Kitaura, J. Chem. Phys. 122 (2005) 054108/1. [21] Y. Mochizuki, S. Koikegami, S. Amari, K. Segawa, K. Kitaura, T. Nakano, Chem. Phys. Lett. 406 (2005) 283. [22] Y. Mochizuki, K. Tanaka, K. Yamashita, T. Ishikawa, T. Nakano, S. Amari, K. Segawa, T. Murase, H. Tokiwa, M. Sakurai, Theor. Chem. Acc. 117 (2007) 541. [23] M. Chiba, D.G. Fedorov, K. Kitaura, Chem. Phys. Lett. 444 (2007) 346. [24] M. Chiba, D.G. Fedorov, K. Kitaura, J. Chem. Phys. 127 (2007) 104108. [25] M. Chiba, D.G. Fedorov, K. Kitaura, J. Comput. Chem. 29 (2008) 2667. [26] B. Auer, M.V. Pak, S. Hammes-Schiffer, J. Phys. Chem. C 114 (2010) 5582. [27] S.R. Pruitt, D.G. Fedorov, K. Kitaura, M.S. Gordon, J. Chem. Theory Comp. 6 (2010) 1. [28] Y. Komeiji, T. Nakano, K. Fukuzawa, Y. Ueno, Y. Inadomi, T. Nemoto, M. Uebayasi, D.G. Fedorov, K. Kitaura, Chem. Phys. Lett. 372 (2003) 342. [29] T. Nagata, D.G. Fedorov, K. Kitaura, Chem. Phys. Lett. 475 (2009) 124. [30] K. Kitaura, S.I. Sugiki, T. Nakano, Y. Komeiji, M. Uebayasi, Chem. Phys. Lett. 336 (2001) 163. [31] D.G. Fedorov, T. Ishida, M. Uebayasi, K. Kitaura, J. Phys. Chem. C 111 (2007) 2722.

[32] Y. Komeiji, Y. Mochizuki, T. Nakano, D.G. Fedorov, J. Mol. Struct. (Theochem) 898 (2009) 2. [33] D.G. Fedorov, J.H. Jensen, R.C. Deka, K. Kitaura, J. Phys. Chem. A 112 (2008) 11808. [34] Y. Yamaguchi, I. Schaefer, F. Henry, Y. Osamura, J. Goddard, A New Dimension to Quantum Chemistry: Analytical Derivative Methods in Ab Initio Molecular Electronic Structure Theory, Oxford University Press, New York, 1994. [35] D.G. Fedorov, K. Kitaura, J. Comput. Chem. 28 (2007) 222. [36] S. Honda, K. Yamasaki, Y. Sawada, H. Morii, Structure 12 (2004) 1507. [37] T. Nagata, D.G. Fedorov, K. Kitaura, M.S. Gordon, J. Chem. Phys. 131 (2009) 024101. [38] M.W. Schmidt, K.K. Baldridge, J.A. Boatz, S.T. Elbert, M.S. Gordon, J.J. Jensen, S. Koseki, N. Matsunaga, K.A. Nguyen, S. Su, T.L. Windus, M. Dupuis, J.A. Montgomery, J. Comput. Chem. 14 (1993) 1347. [39] M.S. Gordon, M.W. Schmidt, Theory and Applications of Computational Chemistry, the First Forty Years, Elsevier, Amsterdam, 2005. [40] D.G. Fedorov, R.M. Olson, K. Kitaura, M.S. Gordon, S. Koseki, J. Comput. Chem. 25 (2004) 872. [41] D.G. Fedorov, K. Ishimura, T. Ishida, K. Kitaura, P. Pulay, S. Nagase, J. Comput. Chem. 28 (2007) 1476. [42] T. Nagata, D.G. Fedorov, T. Sawada, K. Kitaura, M.S. Gordon, in preparation. [43] T. Fujita, H. Watanabe, S. Tanaka, J. Phys. Soc. Jpn. 78 (2009) 104723. [44] T. Fujiwara, Y. Mochizuki, Y. Komeiji, Y. Okiyama, H. Mori, T. Nakano, E. Miyoshi, Chem. Phys. Lett. 490 (2010) 41. [45] D.A. Varshalovich, L.N. Moskalev, V.K. Kersunskii, Quantum Theory of the Angular Momentum, Nauka, Moscow, 1975 (in Russian). [46] D.G. Fedorov, M.S. Gordon, in: M.R. Hoffmann, K.G. Dyall (Eds.), Low-lying Potential Energy Surfaces, ACS Symposium Series 828, Washington, DC, 2002.