Analytic second derivative of the energy for density functional theory based on the three-body fragment molecular orbital method Hiroya Nakata, Dmitri G. Fedorov, Federico Zahariev, Michael W. Schmidt, Kazuo Kitaura, Mark S. Gordon, and Shinichiro Nakamura Citation: The Journal of Chemical Physics 142, 124101 (2015); doi: 10.1063/1.4915068 View online: http://dx.doi.org/10.1063/1.4915068 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/142/12?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Analytic second derivatives of the energy in the fragment molecular orbital method J. Chem. Phys. 138, 164103 (2013); 10.1063/1.4800990 Unrestricted Hartree-Fock based on the fragment molecular orbital method: Energy and its analytic gradient J. Chem. Phys. 137, 044110 (2012); 10.1063/1.4737860 Time-dependent density functional theory based upon the fragment molecular orbital method J. Chem. Phys. 127, 104108 (2007); 10.1063/1.2772850 Improved pseudobonds for combined ab initio quantum mechanical/molecular mechanical methods J. Chem. Phys. 122, 024114 (2005); 10.1063/1.1834899 Ab initio molecular orbital and density functional characterization of the potential energy surface of the N 2 O+Br reaction J. Chem. Phys. 109, 9410 (1998); 10.1063/1.477602

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: 137.30.242.61 On: Wed, 10 Jun 2015 20:13:34

THE JOURNAL OF CHEMICAL PHYSICS 142, 124101 (2015)

Analytic second derivative of the energy for density functional theory based on the three-body fragment molecular orbital method Hiroya Nakata ( 中田浩弥 ),1,2,3,a) Dmitri G. Fedorov,4,b) Federico Zahariev,5 Michael W. Schmidt ( 守美都舞空 ),5 Kazuo Kitaura ( 北浦和夫 ),6 Mark S. Gordon ( 強首領真空 ),5 and Shinichiro Nakamura ( 中村振一郎 )2 1

Center for Biological Resources and Informatics, Tokyo Institute of Technology, B-62 4259 Nagatsuta-cho, Midori-ku, Yokohama 226-8501, Japan 2 RIKEN, Research Cluster for Innovation, Nakamura Lab, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan 3 Japan Society for the Promotion of Science, Kojimachi Business Center Building, 5-3-1 Kojimachi, Chiyoda-ku, Tokyo 102-0083, Japan 4 NRI, National Institute of Advanced Industrial Science and Technology (AIST), 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan 5 Department of Chemistry and Ames Laboratory, US-DOE, Iowa State University, Ames, Iowa 50011, USA 6 Graduate School of System Informatics, Kobe University, 1-1 Rokkodai-cho, Nada-ku, Kobe 657-8501, Japan

(Received 26 January 2015; accepted 4 March 2015; published online 23 March 2015) Analytic second derivatives of the energy with respect to nuclear coordinates have been developed for spin restricted density functional theory (DFT) based on the fragment molecular orbital method (FMO). The derivations were carried out for the three-body expansion (FMO3), and the two-body expressions can be obtained by neglecting the three-body corrections. Also, the restricted HartreeFock (RHF) Hessian for FMO3 can be obtained by neglecting the density-functional related terms. In both the FMO-RHF and FMO-DFT Hessians, certain terms with small magnitudes are neglected for computational efficiency. The accuracy of the FMO-DFT Hessian in terms of the Gibbs free energy is evaluated for a set of polypeptides and water clusters and found to be within 1 kcal/mol of the corresponding full (non-fragmented) ab initio calculation. The FMO-DFT method is also applied to transition states in SN2 reactions and for the computation of the IR and Raman spectra of a small Trp-cage protein (PDB: 1L2Y). Some computational timing analysis is also presented. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4915068]

I. INTRODUCTION

The first and second derivatives of the energy with respect to nuclear coordinates are indispensable for quantum-chemical (QM) modelling in chemistry; for instance, for geometry optimizations, transition state searches, and the prediction of vibrational frequencies. The first derivative (gradient) is used in molecular dynamics (MD) simulations, because the forces on the atoms are given by the gradient. The second derivative (Hessian) is used to calculate harmonic frequencies, and many properties related to them, such as the zero point energy (ZPE), the vibrational contribution to the free energy, infrared (IR) intensities, and Raman activities. The Hessian of the QM energy has been introduced by pioneers such as Gerratt and Mills,1 and Pulay.2 It is usually evaluated by solving coupled perturbed Hartree-Fock (CPHF) equations,3 that are also formulated for other approaches such as the multi-configurational self-consistent field,4 coupled cluster theory with both closed and open shells,5–7 as well as for second order Møller-Plesset perturbation theory (MP2),8 configuration interaction,9 and density functional theory (DFT).10–12

a)Electronic mail: [email protected]. Fax: +81-48-467-8503. b)Electronic mail: [email protected]

0021-9606/2015/142(12)/124101/10/$30.00

The cost of computing the Hessian is very high even for single reference methods, thus, the calculations are limited to small and medium sized systems. Among single reference approaches, one of the most widely used is DFT,13 which typically involves the numeric integration of exchangecorrelation functionals on a grid.14 For computing the Hessian in large molecular systems, there are several options: combined QM and molecular mechanics (MM),15 n-layered integrated molecular orbital and molecular mechanics (ONIOM),16 and density-functional tight-binding (DFTB).17 To accelerate the Hessian evaluation for a subset of atoms, usually, in a numeric fashion, a number of partial schemes is available.18–21 The linearscaling QM algorithms22–25 have primarily been developed for energy and gradient evaluations; the efficiency of Hessian computations has also been improved.26–28 Some fragment based approaches29–40 feature analytic second derivatives.41–46 The fragment molecular orbital (FMO) method47–51 is a fragment-based approach. In the FMO method, the system is divided into a set of fragments (also called monomers). The electronic state of a fragment is embedded into an electrostatic potential (ESP) that is determined by the nuclei and electron densities of the other fragments. The monomer calculations are thus coupled and iterated to self-consistency. Next, dimer and trimer calculations are performed in the ESP of monomers (fixed at this point). The fragmentation

142, 124101-1

© 2015 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: 137.30.242.61 On: Wed, 10 Jun 2015 20:13:34

124101-2

Nakata et al.

J. Chem. Phys. 142, 124101 (2015)

reduces the computational timing significantly. Hence, the FMO method has been applied to proteins,52,53 DNA,54 and inorganic systems.55–60 The discussion in this work focuses on the three-body FMO method (FMO3)61–63 combined with DFT.64,65 The corresponding two-body FMO expressions can be easily obtained by ignoring the three-body correction terms. A comparison of the accuracy of two and three-body expansions of DFT for the energy and gradient is found elsewhere,65 and the cost scaling of FMO with system size has already been discussed in general.48 The FMO analytic energy gradient has been developed for closed and open-shell Hartree-Fock (HF)66–71 and DFT.72 The analytic second derivative, which can be used to evaluate Raman activities,73 has been developed only for HF46,74 at the two-body level (FMO2). The FMO-RHF Hessian is not fully analytic, because some terms (e.g., contributions from second order responses) that are expected to be small are neglected.46,74 All of these terms are also neglected in FMODFT Hessians developed here, and no further approximations are introduced. In Secs. II–V, the FMO3 Hessian is first formulated for the general case. Then, the detailed derivations, specifically for DFT, are presented. All of the DFT-related exchangecorrelation terms introduced in this work are fully analytic. The accuracy of the FMO Hessian is evaluated by comparing some properties to those obtained with full (non-fragmented) DFT. Simulations of the IR and Raman spectra of the Trp cage protein (1L2Y)75 and a study of SN2 reactions in solution are also presented. The computational efficiency of FMO/DFT Hessians is illustrated for water clusters.

∂ 2 E FMO3  ∂ 2 E I′  ∂ 2∆E I J = + ∂a∂b ∂a∂b I >J ∂a∂b I  ∂2 [∆E I J K − ∆E I J − ∆EJ K − ∆EK I ] . + ∂a∂b I >J >K (6) To obtain the analytic second derivative of the FMO3 total energy, one needs to evaluate the X-mer quantities ∂ 2 E X′ /∂a∂b. The derivatives of the many-body terms ∆E I J and ∆E I J K include two types of contributions: (a) differences
of E X′ values, and (b) density transfer terms Tr ∆D X V X . For X = I J, the corresponding second derivative expressions can be found in the FMO2-HF Hessian derivation.46,74 The following focuses on the new terms that arise for FMO3 and for DFT. DFT specific terms arise only in the derivatives of the internal energies E X′ , because the ESP does not contain exchange-correlation terms.

B. Analytic second derivative of the internal energy for DFT

The second derivative of the internal fragment energy is ∂ E X′ ∂a∂b occ    hiiab, X + Piiab, X + Fii′ab, X − Siiab, X ϵ iX = 2

i ∈X

+4

Sia,j X Sib,j X ϵ iX

i, j ∈X vir  occ  b, X + Umi

II. THEORY

m ∈X i ∈X

A. Analytic second derivative for FMO3

+

61,62

The FMO3 total energy is given by   E FMO3 = E I′ + ∆E I J 

[∆E I J K − ∆E I J − ∆EJ K − ∆EK I ] ,

occ 

Sib,j X

i, j ∈X

I >J

I

+

occ 

+ (1)

I >J >K

where the two- and three-body energy contributions are, respectively,
∆E I J = E I′ J − EJ′ − E I′ + Tr ∆D I J V I J , (2)
′ IJK IJK ′ ′ ′ . (3) ∆E I J K = E I J K − EK − EJ − E I + Tr ∆D V E X′ is the internal energy of monomers (X = I), dimers (X = I J), or trimers X = I JK, and V X is the corresponding electrostatic potential. ∆D I J (∆D I J K ) is the difference between the density matrices of a dimer (trimer) and monomers, defined as
∆D I J = D I J − D I ⊕ DJ , (4)
IJK IJK I J K ∆D =D − D ⊕D ⊕D . (5) The second derivative of the energy with respect to nuclear coordinates a and b is

occ  i, j ∈X

Sia,j X

  occ  4F ′a, X − 4S a, X ϵ X − 2 a, X ′ S j k A j k, mi  mi i  i m  j,k ∈X    occ  2F ′a, X − 1  a, X ′ S A  i j 2 k,l ∈X k l i j,k l      occ  2F ′b, X − 1  b, X ′ S A  i j 2 k,l ∈X k l i j,k l   

X ∂ 2 ENR ∂ 2 Exc ab, X, X + −U , (7) ∂a∂b ∂a∂b where i, j, k, l, and m denote molecular orbitals (MO) (m is used for virtual orbitals, while i, j, k, and l are used for occupied orbitals). hiiab, X and Piiab, X are the second derivatives of the one-electron core integrals and the hybrid projection operator integrals, respectively, with respect to coordinates a X and b. Sia,j X and Siab, are the first and second derivatives of the j X X overlap integrals, respectively. Fi′a, and Fi′ab, are matrices j j formed from the first and second integral derivatives of the internal Fock matrix, respectively. Fi′,j X is obtained from the conventional Fock matrix FiXj by subtracting the ESP,

+

Fi′,j X = FiXj − ViXj .

(8)

Uia,j X are the orbital response terms, related to the dependence of the MO coefficients on the nuclear coordinates,3 ϵ iX is

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: 137.30.242.61 On: Wed, 10 Jun 2015 20:13:34

124101-3

Nakata et al.

J. Chem. Phys. 142, 124101 (2015) ab, X, X

the energy of MO i in fragment X. U is the second order orbital response, whose definition is provided in Subsection II C. To extend the HF formulation to DFT, one has to modify the internal energy expression by adding the exchangecorrelation energy derivatives,   X ∂ 2 Exc ∂ 2     = w f (ρ , ρ , γ , γ , γ ) At At α β αα α β β β  , ∂a∂b ∂a∂b  A∈X t ∈ A   (9) where w At and f At are the quadrature weights and the exchange-correlation functional at the grid point t belonging to atom A. ρ and γ are the electronic density and its derivatives. α and β are spin labels. The exchange-correlation contributions are evaluated by numerical integration over all atoms A and grid points t. Consequently, the internal Fock matrix Fi′,j X and the orbital Hessian matrix Ai′ j,k l for DFT become Fi′,j X = hiXj +

occ 


∂ 2Tr ∆D I J K V I J K ∂a∂b   ∂ 2   IJK IJK I IJK  = D µν Vµν − D µν Vµν ∂a∂b  µ,ν ∈I J K µ,ν ∈I     K I JK J IJK D µνVµν  , (13) D µνVµν − −  µ,ν ∈K µ,ν ∈J where µ and ν are the atomic orbital indices. Eq. (13) is a combination of four terms, each of which is given by (for X = I JK, I, J, or K) ( ) X IJK  ∂ 2 D µν Vµν

=

+2

Viiab, I J K −

occ  occ ( 

Via,j I J K Sib,j X + Vib,j I J K Sia,j X

)

i ∈X j ∈X vir  occ 

b, X a, I J K Umi Vmi + Rab, X + U

ab, X, I J K

,

m ∈X i ∈X

k ∈X

(10)

(14)

(11)

where Rab, X is a two-electron integral response term.79 ab, X, I J K U is the orbital response term, whose definition is given below. Inserting Eqs. (7)–(14) into Eq. (6), one can calculate the FMO3 Hessian. ab, X,Y The second order orbital response U (Y , like X, denotes some monomer, dimer, or trimer) appears in the above ab, X, X expressions only in the Y = X or Y = I JK forms, U and ab, X, I J K U . The dimer terms for Y = I J are also present in the two-body expressions, which are not explicitly listed above and can be found elsewhere.46 The second order responses arise from the derivative vir occ   a, X Y of Umi Vmi with respect to nuclear coordinate b. An

where cHF represents the fraction of HF exchange in the functional: cHF = 0 for pure DFT (no HF exchange), and cHF is between 0 and 1 for hybrid functionals. (Vxc)i j and ( f xc)i j,k l are the matrix elements of the first and second functional derivatives of the exchange-correlation functional with respect to the electron density. In this study, the Becke multicenter numerical integration scheme14 was applied, but with some minor modifications according to Murray et al.76 The second derivative of the exchange-correlation energy contribution can be written as X    ∂ 2w At ∂ 2 Exc = f At (ρα , ρ β , γαα , γα β , γ β β ) ∂a∂b A∈X t ∈ A ∂a∂b ∂ 2 f At (ρα , ρ β , γαα , γα β , γ β β ) ∂a∂b ∂w At ∂ f At (ρα , ρ β , γαα , γα β , γ β β ) + ∂a ∂b  ∂w At ∂ f At (ρα , ρ β , γαα , γα β , γ β β ) + . ∂b ∂a

occ  i ∈X

{2(i j|k k) − cHF(ik| j k)} + PiXj + (Vxc)iXj ,

Ai′ j,k l = 4(i j|kl) − cHF {(ik| jl) − (il| j k)} + 4( f xc)i j,k l ,

∂a∂b

µ,ν ∈X

+ w At

(12)

All of these terms, including the second derivatives of the quadrature weights, have been implemented in the quantum chemistry program GAMESS, and the other terms are the same as our previous FMO-HF Hessian.46,74 It is important to note that the contribution of the weight derivative terms should be included even when very fine grids are used, contrary to what is suggested in other studies.77,78 Because of minor differences to other implementations, additional details of the derivations are provided for interested readers.79

m i

important finding in the previous studies46,67 is that the sum of all orbital response terms cancels out in FMO2 gradients and Hessians if the ESP point charge (PC) approximation, ESPPC, is applied to all ESPs, or if the ESP-PC approximation is not used at all;80,81 otherwise, the orbital response terms need to be evaluated.82 Now, it is shown that the orbital response ab, X,Y terms U also cancel out for FMO3, for the same two conditions. The combination of all three-body orbital response terms appearing in the last sum of Eq. (6) is as follows:  ( ab, I, I J K ab,FMO3 ab, J, I J K ab, K, I J K ∆U = −U −U −U I >J >K ab, I, I J

+U +U

+U

ab, K, J K

−U

ab, I, I

ab, J, I J

+U

−U

+U

ab, I, I K

ab, J, J

ab, J, J K

+U

−U

ab, K, I K

ab, K, K

)

.

(15)

ab, I J K, I J K

C. Analytic second derivative of the density transfer terms

The second derivative of the three-body ESP contribution is

Note that U cancels out, because it appears with opposite signs in the internal energy and the density transfer terms, see Eqs. (7) and (14). This cancellation occurs because the ESP is included in the DFT (HF) trimer calculations, and the ESP term was only formally separated into a different term

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: 137.30.242.61 On: Wed, 10 Jun 2015 20:13:34

124101-4

Nakata et al.

J. Chem. Phys. 142, 124101 (2015)

in the above expressions. By separating out the contribution I (J ) Vµν of a single external fragment J,   −Z A   ν + I (J ) J D λσ (µν|λσ), (16) Vµν = µ |r − R | A λ,σ ∈J A∈J one can identify a response contribution U individual fragments I,  ab, I, I (J ) ab, I, X U = U ,

ab, I, I (J )

for the (17)

J ,X

where U

ab, X, I (J )

= −4 +4

occ  occ  vir  i ∈X j ∈X m ∈X occ  occ  occ 

b, X X Umi (VmI (Jj ) Sia,j X + ViIj(J ) S a, jm )

Eq. (15), where they cancel out. To obtain the first order response terms Urai , one has to solve the CPHF equations. The CPHF equations in GAMESS have been extended to DFT by using an appropriate DFT form of the orbital Hessian Ai′ j,k l and Fock matrix Fi′,j X , following Eqs. (10) and (11), respectively. IR intensities and Raman activities in FMO-DFT are calculated analogously to FMO-RHF73 using DFT instead of RHF in a trivial extension. Finally, note that there is a FMO3 specific approximation,61 that neglects the contribution of those trimers in which fragments are far apart. This approximation may be used with all of the above derivations, as it simply restricts the sum of  from running over all distinct fragment triples (trimers) I >J >K

to a subset of compact trimers. b, X I (J ) a, X Smi Vi j S j m

i ∈X j ∈X m ∈X

−

occ  occ 

III. COMPUTATIONAL DETAILS

X I (J ) 2Siab, Vi j j

i ∈X j ∈X I (J ) a, X vir  occ  ∂(Vmi Umi ) . + 4 ∂b m ∈X i ∈X

(18)

The orbital response terms in Eq. (15), related to fragment I, can be expanded as U

ab, I, I J

ab, I, I K

ab, I, I J K

ab, I, I

+U −U  −U  ab, I, I (L) ab, I, I (L) = U + U L,I, J



−

U

L,I, K  ab, I, I (L)

−

L,I, J, K

U

ab, I, I (L)

= 0. (19)

L,I

Thus, the orbital response terms related to J and K likewise cancel. Therefore, the sum of all three-body terms in Eq. (15) is zero, and it is not necessary to calculate these terms explicitly (unless the ESP-PC approximation is used for some, but not all, ESPs). If two fragments I and J are far from each other, the ESDIM approximation may be used. In this case, ∆E I J obtained from a DFT dimer I J is approximated80 and the contribution ab, I, I to the orbital responses in the last term of Eq. (6), −U ab, J, J ab, I, I J ab, J, I J ab, I, I (J ) −U +U +U , is replaced by −U ab, J, J (I ) −U as derived elsewhere67 for the gradient. A similar ab, I, I ab, J, J relation holds for the Hessian, and −U −U ab, I, I J ab, J, I J ab, I, I (J ) ab, J, J (I ) +U and −U −U are the +U same, because −U −U

ab, I, I

ab, J, J

+U

ab, I, I J

+U

ab, J, I J

= −U

ab, I, I (J )

= −U

ab, J, J (I )

, .

(20)

Therefore, the use of the ES-DIM approximation does not ab,FMO3 affect the fact that ∆U is equal to zero. If the ESP-PC approximation is not used, the second order responses are not needed. The first order responses must be calculated for all monomers, dimers, and trimers (they appear in Eqs. (7) and (14)). To obtain the FMO3 Hessian, one must add terms related to the internal energy derivatives in Eq. (7), as well as the density transfer in Eq. (14). The second order response terms in each of these equations are separated in

First, a conventional analytic DFT Hessian10–12 program was implemented in GAMESS,83,84 and parallelized similarly to the existing RHF Hessian.85 Then, the FMO-DFT analytic Hessian was implemented in GAMESS and was parallelized using the generalized distributed data interface (GDDI).61,86,87 For all calculations, the 6-31G(d) basis set with spherical harmonics was used. First, the accuracy of the FMO3-DFT Hessian is evaluated in comparison to the full DFT Hessian. For (H2O)10, three different types of exchange-correlation functionals (B3LYP, BLYP, and LC-BLYP) were used,88–90 with the default grid (96 radial points in the Euler-MacLaurin quadrature and 302 angular points in the Lebedev grid). For a test of the accuracy in systems in which fragments were connected by covalent bonds, a capped alanine decamer (ALA)10 and a styrene oligomer were also calculated (Fig. 1). Second, the FMO-DFT Hessian (with the LC-BLYP functional) was calculated to obtain the vibrational spectrum of a small Trp-cage miniprotein construct (PDB: 1L2Y), consisting of 304 atoms, whose geometry was optimized for aqueous solution using the polarizable continuum model (PCM) at the FMO3-DFT/PCM⟨1⟩ level91 until the root mean square gradient dropped below 4 × 10−4 hartree/bohrs. Third, the transition states for the SN2 reaction XCH3 + OH− → CH3 + X−, X = Br, Cl, or F, were studied in explicit solvent, namely, 41 waters (the total number of atoms is 130), using B3LYP. In the calculations of the PCM solvated Trp-cage protein and explicitly solvated SN2 reactions, an empirical dispersion (D) correction version 392 was used (denoted DFTD), with the standard grid (SG1)93 in DFT. Finally, the computational efficiency was demonstrated on a cubic box of water molecules, containing 16, 21, 54, 122, 136, and 165 water molecules using BLYP. With the exception of the styrene oligomer, all FMO calculations were performed with FMO3; in styrene, the magnitude of the three-body corrections is small, so FMO2 was used. Hessian calculations were performed on stationary points. Unless otherwise indicated, geometry optimizations and transition state searches used a convergence criterion of 0.0001 hartree/bohr, the ES-DIM approximation was used

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: 137.30.242.61 On: Wed, 10 Jun 2015 20:13:34

124101-5

Nakata et al.

J. Chem. Phys. 142, 124101 (2015)

applied. The vibrational (vib) contributions are given by96 1 νi , (21) ZPE = 2 i  νi Hvib = + ZPE, (22) νi exp( kT ) − 1 i  (  νi νi ) Svib = ) − k ln 1 − exp(− , νi kT T(exp( kT ) − 1) i (23) Gvib = Hvib − T Svib,

(24)

where νi is the vibrational frequency of normal mode i, k is the Boltzmann constant, and T is the temperature. Thus, the higher the frequency, the larger is the contribution to the ZPE. On the other hand, as discussed elsewhere, for the vibrational partition function and quantities derived from it (the entropy and free energies), low frequency vibrations make a large contribution.46 The vibrational contributions shown above are added to the translational and rotational terms to obtain the full free energies H and G. Therefore, comparisons between the full ab initio and FMO predictions for ZPE, H, and G are good benchmarks for the accuracy of high and low vibrational frequencies. IV. RESULTS AND DISCUSSION A. FMO accuracy

A summary of prominent IR peaks for (H2O)10 is given in Table I, and IR spectra are shown in Fig. 2(a) for the B3LYP density functional.79 The three peaks in Table I were chosen to represent high, medium, and low frequency modes. Both frequencies and intensities obtained with FMO3-DFT agree reasonably well with the full DFT for all three exchangecorrelation functionals used in this study. For example, in the case of B3LYP, the FMO frequency errors compared to the full DFT Hessian are 2, 4, and 19 cm−1 for ν1, ν2, and ν3, respectively, and the corresponding intensity errors FIG. 1. Fragmentation of (a) capped alanine decamer, (b) styrene oligomer, (c) Trp-cage protein (1L2Y).

with a unitless threshold80 of 2.0, and the default long distance trimer cutoff values were applied for all calculations. The ESPPC approximation was not used during all Hessian calculations. In the simulations of vibrational spectra, Gaussian fitting was used with a peak width of 20 cm−1. For comparison to experiment, all frequencies for the Trp-cage protein were scaled94 by 0.9915. Molecular systems were divided into 1 water molecule or 1 amino acid residue per fragment. For styrene, the system was separated into two monomer units per fragment. For SN2 reactions, the reaction center was treated as 1 fragment, and each water molecule was treated as a fragment. FragIt95 was used to generate input files. In order to calculate the ZPE, Helmholtz free energy (H), and Gibbs free energy (G), statistical thermodynamics based on the harmonic oscillator/rigid rotor approximation is

TABLE I. Three representative IR peaks ν i in (H2O)10, computed with 6-31G(d). The units for frequencies and intensities are cm−1 and D2/(u Å2), respectively. FMO3-DFT Hessian Frequency

B3LYP BLYP LC-BLYP

IR intensity

ν1

ν2

ν3

ν1

ν2

ν3

737 770 839

1763 1735 1720

3329 3187 3232

0.49 0.44 0.50

0.33 0.30 0.29

1.45 1.22 1.49

Full DFT Hessian Frequency

B3LYP BLYP LC-BLYP

IR intensity

ν1

ν2

ν3

ν1

ν2

ν3

735 746 828

1759 1730 1715

3310 3174 3213

0.44 0.39 0.46

0.32 0.27 0.28

1.15 1.22 1.59

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: 137.30.242.61 On: Wed, 10 Jun 2015 20:13:34

124101-6

Nakata et al.

J. Chem. Phys. 142, 124101 (2015) TABLE II. Three prominent IR peaks ν i , computed with LC-BLYP-D/631G(d). The units for frequencies and intensities are cm−1 and D2/(u Å2), respectively. FMO-DFT Hessian Frequency

Polyalaninea Styrene oligomerb

IR intensity

ν1

ν2

ν3

ν1

ν2

ν3

645 721

1795 1523

3523 3065

0.20 0.22

2.36 0.07

1.96 0.17

Full DFT Hessian Frequency

Polyalanine Styrene oligomer

IR intensity

ν1

ν2

ν3

ν1

ν2

ν3

653 722

1797 1523

3536 3064

0.27 0.22

2.17 0.07

1.83 0.18

a FMO3. b FMO2.

FIG. 2. IR spectra of (a) water cluster (H2O)10 computed with B3LYP/631G(d), (b) alanine decamer computed with LC-BLYP-D/6-31G(d), (c) styrene oligomer computed with LC-BLYP-D/6-31G(d). FMO-DFT and full DFT spectra are shown as red solid and blue dashed lines, respectively.

are 0.05, 0.01, and 0.30 (D2/(u Å2)). In general, the FMO intensities tend to have somewhat larger relative errors than the frequencies,46 due to the fact that intensities are more sensitive to the delocalization of normal modes. Also, the intensity of one of the peaks computed by B3LYP, ν3, features a deviation of about 20%, the largest intensity deviation. This peak corresponds to one of the asymmetric stretches in water, whose normal modes extend over multiple water molecules according to the full DFT calculation. By including trimers of water molecules in FMO3, it is possible to accurately evaluate all of the frequencies, and the error in the intensity is 20% or less (depending on the functional). In general, Fig. 2(a) shows that FMO3 reliably predicts the IR spectra of water clusters.

Results for the polyalanine and styrene oligomer are shown in Figs. 2(b) and 2(c), and Table II. Once again, the FMO-DFT results are in reasonable agreement with full DFT. The errors of the three prominent peaks in polyalanine are 8, 2, and 13 cm−1, and the errors in IR intensity are 0.07, 0.19, and 0.13 D2/(u Å2), respectively. In the styrene oligomer, the errors in vibrational frequency and IR intensity do not exceed 1 cm−1 and 0.01 D2/(u Å2), respectively. Water and the polypeptides have hydrogen bonding with significant threebody effects, whose accurate treatment necessitates FMO3, whereas in the styrene oligomer, the three-body effects are relatively small and FMO2 was used. The root mean square deviation (RMSD) between FMO-DFT and full DFT for all vibrational frequencies for the alanine decamer is 2.0 cm−1, and 0.5 cm−1 for the styrene oligomer. The maximum error among all frequencies is 10 cm−1 for the alanine decamer, and 5 cm−1 for the styrene oligomer. A comparison of the vibrational energies is given in Table III. The differences between FMO-DFT and full DFT quantities for ZPE, Gibbs free energy (G), and Helmholtz free energies (H) are small: the errors in polyalanine are 0.19, 0.47, and 0.09 kcal/mol, respectively; for the styrene oligomer, they are 0.01, 0.98, and 0.55 kcal/mol, respectively. For the water cluster, the largest error is 0.56 kcal/mol. Using FMO3 may presumably reduce the error of 0.98 kcal/mol in the FMO2 Gibbs free energy for the styrene oligomer. The fragmentation errors in the FMO3 thermodynamic quantities described here are smaller than those previously reported46 for FMO2-HF. The Gibbs free energies reported here are accurate to 1 kcal/mol, although they depend on the low frequency vibrations via the entropy, and the low frequency vibrations are often delocalized and notoriously difficult to estimate. B. IR and Raman spectra of the Trp-cage protein

The IR and Raman spectra of the Trp-cage protein were simulated with the FMO3-DFT-D Hessian, computed at the FMO3-DFT-D/PCM⟨1⟩ optimized geometry. The RMSD of the optimized structure from the NMR experiment (the first

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: 137.30.242.61 On: Wed, 10 Jun 2015 20:13:34

124101-7

Nakata et al.

J. Chem. Phys. 142, 124101 (2015)

TABLE III. ZPE and vibrational contributions to the Gibbs (Gvib) and Helmholtz (Hvib) free energies, computed with 6-31G(d) at 298.15 K, all in kcal/mol.

TABLE IV. Prominent IR peaks (in cm−1) computed with FMO3-DFT (LCBLYP-D/6-31G(d)) compared to experimental resonance Raman peaks in the Trp-cage protein (1L2Y).

FMO-DFT Hessiana

Calculateda

Experiment97

1709–1726 1563–1577 1259–1268 1217

1662 1564 1266 1210

System b

(H2O)10 (H2O)10 c (H2O)10 d Polyalanined Styrene oligomerd

ZPE

Gvib

Hvib

159.79 154.51 159.79 600.81 713.55

126.74 122.09 128.15 543.83 653.85

174.96 169.36 173.83 638.61 748.25

Full DFT Hessian System b

(H2O)10 (H2O)10 c (H2O)10 d Polyalanined Styrene oligomerd a FMO2

ZPE

Gvib

Hvib

159.71 154.19 159.40 601.00 713.56

126.90 121.65 127.95 544.33 652.87

174.84 169.15 173.47 638.70 748.80

for styrene oligomer, otherwise, FMO3.

b B3LYP. c BLYP. d LC-BLYP.

structure in the PDB)75 is 0.8343 Å. The largest contributions to the RMSD are in flexible side chains and the two termini, whereas the backbones are similar. The calculated and experimental structures are shown in Fig. 3(a). A comparison of the calculated and experimental vibrational frequencies is given in Table IV, where a range of

Amide I Amide II Amide III Tyr symmetric a Scaled94

by 0.9915.

calculated frequencies for each type of vibration is specified for normal modes of the same nature (because of the coupling of multiple groups of the same kind). Most of the calculated vibrational frequencies agree well with the experiment.97 The difference is 7 cm−1 for the Tyr symmetric mode, and the Amide II and III peaks very closely reproduce the experimental values. On the other hand, the Amide I peak has a deviation of about +50 cm−1 from the experiment. The reason for the discrepancy may be the lack of solvent in the Hessian calculations. It was previously reported that Amide I peaks are shifted by −30 cm−1 by adding explicit solvent in a small polypeptide,74 because of the coupling of the Amide I normal modes to water bending modes. Likewise, COO−1 modes are strongly affected by solvent (+59 cm−1 as reported elsewhere74), whereas Amide II were found74 to be weakly affected. The calculated IR and Raman spectra of the Trp-cage protein are shown in Figure 3(b). The vibrational modes around 1611 cm−1 are mainly localized on the ammonium cation in Arg, and the modes around 1769 cm−1 are Amide I vibrational modes that reflect inter-residue hydrogen bond coupling in the protein. The most Raman active vibrational peak is the CH2 stretching mode along the backbone main chain; these CH2 modes are not IR active modes, since the dipole derivative contributions are expected to be small. The symmetric and antisymmetric COO−1 stretching modes are observed at 1445 and 1718 cm−1, respectively, whose frequencies may be underestimated somewhat by the omission of solvent, if the trend found in another polypeptide74 holds in general.

C. Transition state analysis for SN2 reactions

FIG. 3. (a) Calculated (optimized with FMO3-DFT-D/PCM⟨1⟩ using LCBLYP-D/6-31G(d)) and experimental (NMR) structure of the Trp-cage protein (1L2Y) is shown in red and green, respectively. (b) Calculated IR and Raman spectra, shown as red solid and blue dashed lines, respectively.

The activation free energies for the SN2 reactions between XCH3 and OH− (X = Br, Cl, or F) immersed in 41 water molecules are summarized in Table V. The errors in the free energy barrier compared to full DFT are 2-3 kcal/mol. The trends within the series of halogen atoms are correctly described by the FMO3-DFT method. In order to improve the accuracy in practical applications, one could try to merge some water molecules of the first solvation shell with the reaction center into one fragment, however, this was not investigated in the present work. To understand their origin, some components of the total free energy barriers ∆G are listed in Table V. High frequencies that make large contributions to the ZPE and low frequencies that make large contributions to the vibrational free energy have some effect. Also, the error due to vibrational terms is comparable to the error in the electronic energy. This system

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: 137.30.242.61 On: Wed, 10 Jun 2015 20:13:34

124101-8

Nakata et al.

J. Chem. Phys. 142, 124101 (2015)

TABLE V. Imaginary vibrational frequencies ω (cm−1) and contributions (kcal/mol) to the SN2 reaction barriers for XCH3 and OH− (X = Br, Cl, and F) solvated in 41 water molecules: electronic energy ∆E, zero point energy ∆ ZPE, vibrational contribution ∆Gvib to the Gibbs free energy ∆G (computed with B3LYP/6-31G(d) at T = 298.15 K). Method FMO3-DFT Full-DFT FMO3-DFT Full DFT FMO3-DFT Full DFT FMO3-DFT Full DFT FMO3-DFT Full DFT

∆E ∆E ∆ ZPE ∆ ZPE ∆Gvib ∆Gvib ∆G ∆G ω ω

F

Cl

Br

24.57 23.88 0.59 0.01 3.63 1.39 28.20 25.27 486 470

14.18 15.05 0.28 1.33 3.38 4.55 17.56 19.60 419 419

7.33 5.66 3.14 3.10 4.56 3.84 11.89 9.50 324 374

is an anion in explicit solvent; due to solute-solvent charge delocalization, some FMO-DFT error is observed. The imaginary vibrational frequencies obtained with the FMO-DFT Hessian agree reasonably well with those obtained with the full DFT Hessian. The larger error found for Br is attributed to the fact that the energy barrier is small, and the potential energy surface is flat. Nevertheless, the monotonic decrease of frequency from F to Cl to Br is properly reproduced by the FMO method. The decrease in the barrier height and imaginary frequency in the series F, Cl, and Br is attributed to the larger atomic radius and the core shielding effect on the valence electrons, as well as the increase in the atomic mass. D. Computational efficiency

The timings for the computation of the FMO3 Hessians were measured on a PC cluster containing 16 2.93 GHz Xeon nodes (8 CPU cores and 12 GB memory per node). For this purpose, several water molecular clusters were used, with system sizes ranging from 16 to 165 water molecules. The results are shown in Fig. 4. The FMO3-BLYP Hessian took about 21 h for the largest system (495 atoms) in this study. For comparison, the timing of the largest cluster for B3LYP is 20.5 h, while for BLYP it is 21.4 h, as self-consistent monomer calculations for B3LYP converge faster than for BLYP. The largest water cluster in this study (495 atoms) took 21 and 18 h with FMO3-DFT and FMO3-HF, respectively, and the overall trend in Fig. 4 shows that DFT and HF Hessians for FMO are comparable in terms of computational time requirements. This is because the main cost is in the evaluation of two-electron integrals and orbital responses, while the extra cost of DFT-related grid integration is relatively minor. Finally, it is of interest to compare the time requirements for FMO3-DFT to full DFT Hessian calculations. The latter demands a significant amount of memory, so a different single node was used, with 24 cores and 128 GB of memory. Total run times for FMO3-DFT or full DFT Hessians are 50 or 187 min for (H2O)16, respectively, and 248 or 5142 min for (H2O)31. Thus, FMO3 decreases the run time by a factor of 21 for 31 waters, which is roughly the largest full Hessian run that may

FIG. 4. Computational timings of FMO3-BLYP/6-31G(d) Hessian calculations of water clusters consisting of 48, 93, 162, 249, 366, 408, and 495 atoms, measured on a PC cluster containing 16 2.93 GHz Xeons nodes (8 CPU cores and 12 GB memory per node). BLYP and HF data are shown as red solid and blue dashed lines, respectively. The number of water molecules is shown above each data point.

be done with the memory that is available on the nodes used for this study.

V. CONCLUSIONS

Analytic second derivatives of the energy with respect to nuclear coordinates have been derived for closed shell spin-restricted FMO3-DFT and implemented into GAMESS. As described in a previous paper,46,74 certain terms that are expected to be small are neglected in the FMO-RHF and FMO-DFT Hessians, so the FMO Hessian code implemented in GAMESS is not fully analytic but the magnitude of the neglected terms is typically very small. FMO2-DFT and FMO3-HF Hessian calculations are enabled as well. Full DFT Hessians have also been coded, including an implementation for long-range corrected functionals,90 such as LC-BLYP, which improve the treatment of the electrostatics. The proper account of electrostatics is thought to be important in polar systems, such as water or proteins. The new development of FMO Hessians permits simulations of IR and Raman spectra and vibrational contributions to the energetics, as well as DFT transition state searches. It has been demonstrated on a representative set of test systems (IR spectra and chemical reactions in solution) that the accuracy of FMO3-DFT properties derived from the Hessian is reasonably high. The FMO approach significantly reduces the computational timing (by a factor of 21 for 93 atoms as reported in this study), and the memory requirements, both of which preclude full DFT Hessian evaluations in large molecular systems. As a demonstrative application, a FMO3-DFT geometry optimization and Hessian evaluation for the small Trp-cage protein was performed. The vibrational frequencies obtained show reasonable agreement with experiment. The present development enables DFT Hessian calculations for large systems. As opposed to the gradient, where individual orbital responses are not needed owing to the formulation of the Z-vector3,66 equations, the Hessian requires

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: 137.30.242.61 On: Wed, 10 Jun 2015 20:13:34

124101-9

Nakata et al.

these responses, and this makes such Hessian calculations rather expensive. Nevertheless, as demonstrated in this work, the FMO3-DFT Hessian is accurate and tractable for systems containing several hundreds of atoms. ACKNOWLEDGMENTS

This work was in part supported by the Next Generation Super Computing Project, Nanoscience Program (MEXT, Japan) and Computational Materials Science Initiative (CMSI, Japan), and JSPS KAKENHI Grant No. 262235. H.N. acknowledges the use of the following computer resources: TSUBAME2.5 at the Global Scientific Information and Computing Center of Tokyo Institute of Technology, RIKEN Integrated Cluster of Clusters (RICC) at RIKEN. M.W.S., F.Z., and M.S.G. were supported by funds provided by the Department of Energy, Basic Energy Sciences, to the Ames Laboratory, administered by Iowa State University. The Ames Laboratory is operated for the U.S. Department of Energy by Iowa State University under Contract No. DE-AC0207CH11358. 1J.

Gerratt and I. M. Mills, J. Chem. Phys. 49, 1719 (1968). Pulay, Mol. Phys. 17, 197 (1969). 3Y. Yamaguchi, H. F. Schaefer III, Y. Osamura, and J. Goddard, A New Dimension to Quantum Chemistry: Analytical Derivative Methods in Ab Initio Molecular Electronic Structure Theory (Oxford University Press, New York, 1994). 4T. U. Helgaker, J. Almlöf, H. J. A. Jensen, and P. Jørgensen, J. Chem. Phys. 84, 6266 (1986). 5H. Koch, H. J. A. Jensen, P. Jørgensen, T. Helgaker, G. E. Scuseria, and H. F. Schaefer, J. Chem. Phys. 92, 4924 (1990). 6J. Gauss and J. F. Stanton, Chem. Phys. Lett. 276, 70 (1997). 7P. G. Szalay, J. Gauss, and J. F. Stanton, Theor. Chem. Acc. 100, 5 (1998). 8J. A. Pople, R. Krishnan, H. B. Schlegel, and J. S. Binkley, Int. J. Quantum Chem. 16, 225 (1979). 9R. Krishnan, H. B. Schlegel, and J. A. Pople, J. Chem. Phys. 72, 4654 (1980). 10A. Komornicki and G. Fitzgerald, J. Chem. Phys. 98, 1398 (1993). 11P. Deglmann, F. Furche, and R. Ahlrichs, Chem. Phys. Lett. 362, 511 (2002). 12B. G. Johnson and M. J. Frisch, Chem. Phys. Lett. 216, 133 (1993). 13A. D. Becke, J. Chem. Phys. 140, 18A301 (2014). 14A. D. Becke, J. Chem. Phys. 88, 2547 (1988). 15Q. Cui and M. Karplus, J. Chem. Phys. 112, 1133 (2000). 16S. Dapprich, I. Komáromi, K. S. Byun, K. Morokuma, and M. J. Frisch, J. Mol. Struct.: THEOCHEM 461-462, 1 (1999). 17H. A. Witek, S. Irle, and K. Morokuma, J. Chem. Phys. 121, 5163 (2004). 18J. D. Head, Int. J. Quantum Chem. 65, 827 (1997). 19W. J. Zheng and B. R. Brooks, Biophys. J. 89, 167 (2005). 20H. Li and J. H. Jensen, Theor. Chem. Acc. 107, 211 (2002). 21A. Ghysels, D. Van Neck, V. Van Speybroeck, T. Verstraelen, and M. Waroquier, J. Chem. Phys. 126, 224102 (2007). 22S. Goedecker and G. E. Scuseria, Comput. Sci. Eng. 5, 14 (2003). 23Linear-Scaling Techniques in Computational Chemistry and Physics, edited by P. G. Mezey and J. Leszczynski (Springer, New York, 2011). 24A. Nakata, D. R. Bowler, and T. Miyazaki, J. Chem. Theory Comput. 10, 4813 (2014). 25K. A. Wilkinson, N. D. M. Hine, and C.-K. Skylaris, J. Chem. Theory Comput. 10, 4782 (2014). 26R. A. Kendall, E. Aprà, D. E. Bernholdt, E. J. Bylaska, M. Dupuis, G. I. Fann, R. J. Harrison, J. Ju, J. A. Nichols, J. Nieplocha, T. Straatsma, T. L. Windus, and A. T. Wong, Comput. Phys. Commun. 128, 260 (2000). 27T. L. Windus, M. W. Schmidt, and M. S. Gordon, Chem. Phys. Lett. 216, 375 (1993). 28T. J. Dudley, R. M. Olson, M. W. Schmidt, and M. S. Gordon, J. Comput. Chem. 27, 352 (2006). 29M. S. Gordon, D. G. Fedorov, S. R. Pruitt, and L. V. Slipchenko, Chem. Rev. 112, 632 (2012). 30P. Otto and J. Ladik, Chem. Phys. 8, 192 (1975). 31W. Yang, Phys. Rev. Lett. 66, 1438 (1991). 2P.

J. Chem. Phys. 142, 124101 (2015) 32J.

Gao, J. Phys. Chem. B 101, 657 (1997). He and K. M. Merz, J. Chem. Theory Comput. 6, 405 (2010). 34M. Kobayashi and H. Nakai, J. Chem. Phys. 138, 044102 (2013). 35J. Friedrich, H. Yu, H. R. Leverentz, P. Bai, J. I. Siepmann, and D. G. Truhlar, J. Phys. Chem. Lett. 5, 666 (2014). 36Y. Tong, Y. Mei, J. Z. H. Zhang, L. L. Duan, and Q. G. Zhang, J. Theor. Comput. Chem. 8, 1265 (2009). 37P. Söderhjelm, J. Kongsted, and U. Ryde, J. Chem. Theory Comput. 6, 1726 (2010). 38A. Frank, H. M. Möller, and T. E. Exner, J. Chem. Theory Comput. 8, 1480 (2012). 39M. S. Gordon, Q. A. Smith, P. Xu, and L. V. Slipchenko, Annu. Rev. Phys. Chem. 64, 553 (2013). 40K. Kiewisch, C. R. Jacob, and J. Visscher, J. Chem. Theory Comput. 9, 2425 (2013). 41S. Sakai and S. Morita, J. Phys. Chem. A 109, 8424 (2005). 42W. Hua, T. Fang, W. Li, J.-G. Yu, and S. Li, J. Phys. Chem. A 112, 10864 (2008). 43M. A. Collins, J. Chem. Phys. 141, 094108 (2014). 44A. P. Rahalkar, V. Ganesh, and S. R. Gadre, J. Chem. Phys. 129, 234101 (2008). 45J. C. Howard and G. S. Tschumper, J. Chem. Phys. 139, 184113 (2013). 46H. Nakata, T. Nagata, D. G. Fedorov, S. Yokojima, K. Kitaura, and S. Nakamura, J. Chem. Phys. 138, 164103 (2013). 47K. Kitaura, E. Ikeo, T. Asada, T. Nakano, and M. Uebayasi, Chem. Phys. Lett. 313, 701 (1999). 48The Fragment Molecular Orbital Method: Practical Applications to Large Molecular Systems, edited by D. G. Fedorov and K. Kitaura (CRC press, Boca Raton, FL, 2009). 49D. G. Fedorov and K. Kitaura, J. Phys. Chem. A. 111, 6904 (2007). 50D. G. Fedorov, T. Nagata, and K. Kitaura, Phys. Chem. Chem. Phys. 14, 7562 (2012). 51S. Tanaka, Y. Mochizuki, Y. Komeiji, Y. Okiyama, and K. Fukuzawa, Phys. Chem. Chem. Phys. 16, 10310 (2014). 52M. P. Mazanetz, O. Ichihara, R. J. Law, and M. Whittaker, J. Cheminf. 3, 2 (2011). 53T. Sawada, D. G. Fedorov, and K. Kitaura, J. Am. Chem. Soc. 132, 16862 (2010). 54T. Watanabe, Y. Inadomi, K. Fukuzawa, T. Nakano, S. Tanaka, L. Nilsson, and U. Nagashima, J. Phys. Chem. B 111, 9621 (2007). 55D. G. Fedorov, P. V. Avramov, J. H. Jensen, and K. Kitaura, Chem. Phys. Lett. 477, 169 (2009). 56P. J. Carlson, S. Bose, D. W. Armstrong, T. Hawkins, M. S. Gordon, and J. W. Petrich, J. Phys. Chem. B 116, 503 (2012). 57P. V. Avramov, D. G. Fedorov, P. B. Sorokin, S. Sakai, S. Entani, M. Ohtomo, Y. Matsumoto, and H. Naramoto, J. Phys. Chem. Lett. 3, 2003 (2012). 58Y. Okiyama, T. Tsukamoto, C. Watanabe, K. Fukuzawa, S. Tanaka, and Y. Mochizuki, Chem. Phys. Lett. 566, 25 (2013). 59L. Roskop, D. G. Fedorov, and M. S. Gordon, Mol. Phys. 111, 1622 (2013). 60Y. Nishimoto, D. G. Fedorov, and S. Irle, J. Chem. Theory Comput. 10, 4801 (2014). 61D. G. Fedorov and K. Kitaura, J. Chem. Phys. 120, 6832 (2004). 62D. G. Fedorov and K. Kitaura, Chem. Phys. Lett. 433, 182 (2006). 63Y. Komeiji, Y. Mochizuki, and T. Nakano, Chem. Phys. Lett. 484, 380 (2010). 64S.-I. Sugiki, N. Kurita, Y. Sengoku, and H. Sekino, Chem. Phys. Lett. 382, 611 (2003). 65D. G. Fedorov and K. Kitaura, Chem. Phys. Lett. 389, 129 (2004). 66K. Kitaura, S. I. Sugiki, T. Nakano, Y. Komeiji, and M. Uebayasi, Chem. Phys. Lett. 336, 163 (2001). 67T. Nagata, K. Brorsen, D. G. Fedorov, K. Kitaura, and M. S. Gordon, J. Chem. Phys. 134, 124115 (2011). 68T. Nagata, D. G. Fedorov, and K. Kitaura, Theor. Chem. Acc. 131, 1136 (2012). 69H. Nakata, M. W. Schmidt, D. G. Fedorov, K. Kitaura, S. Nakamura, and M. S. Gordon, J. Phys. Chem. A 118, 9762 (2014). 70H. Nakata, D. G. Fedorov, T. Nagata, S. Yokojima, K. Ogata, K. Kitaura, and S. Nakamura, J. Chem. Phys. 137, 044110 (2012). 71T. Nagata, D. G. Fedorov, K. Ishimura, and K. Kitaura, J. Chem. Phys. 135, 044110 (2011). 72K. R. Brorsen, F. Zahariev, H. Nakata, D. G. Fedorov, and M. S. Gordon, J. Chem. Theory Comput. 10, 5297 (2014). 73H. Nakata, D. G. Fedorov, S. Yokojima, K. Kitaura, and S. Nakamura, J. Chem. Theory Comput. 10, 3689 (2014). 33X.

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: 137.30.242.61 On: Wed, 10 Jun 2015 20:13:34

124101-10

Nakata et al.

74H. Nakata, D. G. Fedorov, S. Yokojima, K. Kitaura, and S. Nakamura, Chem.

Phys. Lett. 603, 67 (2014). 75J. W. Neidigh, R. M. Fesinmeyer, and N. H. Andersen, Nat. Struct. Biol. 9, 425 (2002). 76C. W. Murray, N. C. Handy, and G. J. Laming, Mol. Phys. 78, 997 (1993). 77M. Malagoli and J. Baker, J. Chem. Phys. 119, 12763 (2003). 78J. Baker, J. Andzelm, A. Scheiner, and B. Technologies, J. Chem. Phys. 101, 8894 (1994). 79See supplementary material at http://dx.doi.org/10.1063/1.4915068 for the details of the grid related derivative terms, two-electron integral response terms, IR spectra of water using other functionals. 80T. Nakano, T. Kaminuma, T. Sato, K. Fukuzawa, Y. Akiyama, M. Uebayasi, and K. Kitaura, Chem. Phys. Lett. 351, 475 (2002). 81T. Nagata, D. G. Fedorov, and K. Kitaura, Chem. Phys. Lett. 475, 124 (2009). 82T. Nagata, D. G. Fedorov, and K. Kitaura, Chem. Phys. Lett. 544, 87 (2012). 83M. 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, and J. A. Montgomery, J. Comput. Chem. 14, 1347 (1993). 84M. S. Gordon and M. W. Schmidt, “Advances in electronic structure theory: GAMESS a decade later,” in Theory and Applications of Computational Chemistry: The First Forty Years, edited by C. E. Dykstra, G. Frenking, K. S. Kim, and G. E. Scuseria (Elsevier, Amsterdam, 2005), pp. 1167–1189.

J. Chem. Phys. 142, 124101 (2015) 85Y.

Alexeev, M. W. Schmidt, T. L. Windus, and M. S. Gordon, J. Comput. Chem. 28, 1685 (2007). 86D. G. Fedorov, R. M. Olson, K. Kitaura, M. S. Gordon, and S. Koseki, J. Comput. Chem. 25, 872 (2004). 87Y. Alexeev, M. P. Mazanetz, O. Ichihara, and D. G. Fedorov, Curr. Top. Med. Chem. 12, 2013 (2012). 88A. D. Becke, J. Chem. Phys. 98, 1372 (1993). 89A. D. Becke, J. Chem. Phys. 104, 1040 (1996). 90H. Iikura, T. Tsuneda, T. Yanai, and K. Hirao, J. Chem. Phys. 115, 3540 (2001). 91T. Nagata, D. G. Fedorov, H. Li, and K. Kitaura, J. Chem. Phys. 136, 204112 (2012). 92S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010). 93P. M. Gill, B. G. Johnson, and J. A. Pople, Chem. Phys. Lett. 209, 506 (1993). 94A. P. Scott and L. Radom, J. Phys. Chem. 100, 16502 (1996). 95C. Steinmann, M. W. Ibsen, A. S. Hansen, and J. H. Jensen, PLoS One 7, e44480 (2012). 96N. Davidson, Statistical Mechanics (McGraw-Hill, New York, 1962). 97Z. Ahmed, I. A. Beta, A. V. Mikhonin, and S. A. Asher, J. Am. Chem. Soc. 127, 10943 (2005).

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: 137.30.242.61 On: Wed, 10 Jun 2015 20:13:34