Quantum mechanical/molecular mechanical/continuum style solvation model: Second order Møller-Plesset perturbation theory Nandun M. Thellamurege, Dejun Si, Fengchao Cui, and Hui Li Citation: The Journal of Chemical Physics 140, 174115 (2014); doi: 10.1063/1.4873344 View online: http://dx.doi.org/10.1063/1.4873344 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/17?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Quantum mechanical/molecular mechanical/continuum style solvation model: Time-dependent density functional theory J. Chem. Phys. 139, 084106 (2013); 10.1063/1.4819139 Analytic gradient for second order Møller-Plesset perturbation theory with the polarizable continuum model based on the fragment molecular orbital method J. Chem. Phys. 136, 204112 (2012); 10.1063/1.4714601 Quadratically convergent algorithm for orbital optimization in the orbital-optimized coupled-cluster doubles method and in orbital-optimized second-order Møller-Plesset perturbation theory J. Chem. Phys. 135, 104103 (2011); 10.1063/1.3631129 Polarizable atomic multipole solutes in a Poisson-Boltzmann continuum J. Chem. Phys. 126, 124114 (2007); 10.1063/1.2714528 Intermolecular potentials of the methane dimer calculated with Møller-Plesset perturbation theory and density functional theory J. Chem. Phys. 125, 094312 (2006); 10.1063/1.2345198

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: 136.167.3.36 On: Sun, 02 Nov 2014 02:47:53

THE JOURNAL OF CHEMICAL PHYSICS 140, 174115 (2014)

Quantum mechanical/molecular mechanical/continuum style solvation model: Second order Møller-Plesset perturbation theory Nandun M. Thellamurege, Dejun Si, Fengchao Cui, and Hui Lia) Department of Chemistry, University of Nebraska-Lincoln, Lincoln, Nebraska 68588, USA

(Received 10 February 2014; accepted 7 April 2014; published online 7 May 2014) A combined quantum mechanical/molecular mechanical/continuum (QM/MM/C) style second order Møller-Plesset perturbation theory (MP2) method that incorporates induced dipole polarizable force field and induced surface charge continuum solvation model is established. The Z-vector method is modified to include induced dipoles and induced surface charges to determine the MP2 response density matrix, which can be used to evaluate MP2 properties. In particular, analytic nuclear gradient is derived and implemented for this method. Using the Assisted Model Building with Energy Refinement induced dipole polarizable protein force field, the QM/MM/C style MP2 method is used to study the hydrogen bonding distances and strengths of the photoactive yellow protein chromopore in the wild type and the Glu46Gln mutant. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4873344] I. INTRODUCTION

Combined quantum mechanical/molecular mechanical/ continuum (QM/MM/C) methods provide an option for modeling solvent effects.1–10 In these methods, a large molecule is divided into two regions, a relatively small QM region and a relatively large MM region, and the bulk solvent is simplified as a dielectric continuum. We have established a rigorous variational treatment of induced dipoles11–15 and induced surface charges for QM/MMpol/C style density functional theory (DFT) and Hartree-Fock (HF) self-consistent field methods,8 as well as time-dependent density functional theory (TDDFT) methods.16 It is important and advantageous to incorporate solvation models into HF and DFT methods in a variational manner because of their fundamental importance in quantum chemistry: analytic gradients and some other properties can be evaluated efficiently, and post-HF and post-DFT methods can be formulated. Second order Møller-Plesset perturbation theory method (MP2), including the spin-restricted closed shell RMP2, spinunrestricted open shell UMP2, and the Z-averaged spinrestricted open shell ZAPT2 methods, is a widely used and relatively accurate post-HF method.17–23 Pople group was the first to derive and implement the MP2 analytic gradient.24 Handy and Schaefer formulated the Z-vector method so MP2 gradient can be evaluated more efficiently.25 The ZAPT2 analytic gradient has been implemented by Fletcher et al.26 and Aikens et al.27 Recently, we implemented analytic gradients for QM/MMpol style and QM/Continuum style RMP2, UMP2, and ZAPT2 methods.28, 29 In this work, we established QM/MMpol/C style RMP2, UMP2, and ZAPT2 meth-

a) Author to whom correspondence should be addressed. Electronic mail:

[email protected]

0021-9606/2014/140(17)/174115/9/$30.00

ods, with an emphasis on deriving and implementing their analytic nuclear gradients. Photoactive yellow protein (PYP) is a 125-residue photoreceptor found in the bacterium Ectothiprhodospira (or halorhodospira) halophila.30 After blue light absorption (maximum 446 nm or 2.78 eV30 ), the chromophore in PYP undergoes an ultrafast trans-to-cis isomerization in the time scale of a few ps.31 Due to its relatively small size and thermal stability, PYP has been a model system for extensive experimental and theoretical studies. In this work we use QM/MMpol/C style MP2 method to study the ground state hydrogen bonding of the PYP chromophore in the wild type protein and its Glu46Gln mutant. According to experimental and computational studies,32–41 in the ground state the PYP chromophore is a deprotonated p-coumaric acid (pCA) in hydrogen bonding with Tyr42 and Glu46 (Figure 1 and Figure 2). When Glu46 is replaced with Gln46, pCA is expected to form a weaker hydrogen bond with Gln46 and stronger hydrogen bond with Tyr42. However, quantitative data about the hydrogen bonding strengths have not been reported. In general, accurate calculations of intermolecular interactions require electron correlation methods such as MP2. When using QM methods to compute the intermolecular interaction and hydrogen bonding energies, geometry optimization should be performed to remove close contacts that can cause strong repulsion, and to find out the best geometry these groups coexist with each other and their environment. In the study of protein active sites, the use of a continuum solvation model for the whole protein to screen the strong charge-charge interaction between ionized surface residues is of critical importance. The QM/MMpol/C style MP2 method developed in this work can be used to optimize the geometry of a protein active site at the MP2 level of theory while taking the protein and bulk solvent interactions into consideration, thus a suitable method for studying intermolecular interactions in proteins.

140, 174115-1

© 2014 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 136.167.3.36 On: Sun, 02 Nov 2014 02:47:53

174115-2

Thellamurege et al.

J. Chem. Phys. 140, 174115 (2014)

Here EHF is the HF electronic energy (including nuclear repulsion energy) of the QM region calculated using the gas phase Hamiltonian but with the force field and solvent perturbed single-determinant wavefunction; EMM is the molecular mechanics force field energy of the MM region (not including the induced dipole polarization energy). Eels in Eq. (1) is the electrostatic interaction energy between QM electrons/nuclei and MM electrostatic potentials. Usually the MM electrostatic potentials are modeled with atomic point charges  Pν Vchg,μν + EN,chg . (2) Eels = μν

Here μ and ν are spin-orbital basis functions; Vchg,μν are the one-electron integrals of the atomic charges; EN,chg is the electrostatic energy between QM nuclei and MM atomic charges; P is the HF density matrix. In the effective fragment potential (EFP)8, 10, 42, 43 method, electric multipole points are used for the MM atoms  Pμν Vmul,μν + EN,mul . (3) Eels =

FIG. 1. The wild type photoactive yellow protein and its active site.

μν

II. THEORY A. QM/MMpol/C style Hartree-Fock gradient

We have established a variational theory for QM/ MMpol/C style Hartree-Fock methods.8 Here the main points and equations are summarized for the purpose of establishing notations for the current paper. The total energy of a molecular system described with a QM/MMpol/C style Hartree-Fock (HF) method can be written as EH F /MMpol/C = EH F + EMM + Eels + Erep + Edisp + Epolsol .

(1)

Here Vmul,μν are the one-electron integrals of the multipole potentials; EN,mul is the electrostatic energy between QM nuclei and MM electric multipole points. Erep and Edisp in Eq. (1) are the repulsion and dispersion interaction energies, respectively, between QM atoms and MM atoms. In most QM/MM methods, they are modeled with Lennard-Jones (LJ) potentials Erep + Edisp = ELJ .

In the EFP method, repulsion-dispersion is described with effective Gaussian potentials that directly act on the QM electrons   Pμν Vrep,μν + Pμν Vdisp,μν . (5) Erep + Edisp = μν

H

H

Tyr42 CH3

H

O

H

H

H O

CH3

H H 3C S

H

H

H

O

p CA

Cys69

H

H

O

Glu46

O

H

H

Tyr42 CH3

H

O

H

H

H

H

O

O

Cys69

H

H

H

p CA

N H

O

Gln46

FIG. 2. The QM region (46 and 47 atoms, respectively) in MP2/ AMBER(pol)/FixSol calculations of wild type photoactive yellow protein and its Glu46Gln mutant.

μν

Here Vrep,μν and Vdisp,μν are the one-electron integrals of the MM repulsion and dispersion potentials. Epolsol in Eq. (1) is the polarization and solvation energy of the entire QM/MMpol/C system, including the induced dipole and induced surface charge contributions, 1 Epolsol = − UT w 2 1 = − (UN + Uels + Ue )T (wN + wels + we ) 2 1 = − (UN + Uels )T (wN + wels ) 2 1  − Pμν Pρσ (Uμν )T wρσ 2 μνρσ −

CH3

H H 3C S

H

(4)

1 Pμν [(UN + Uels )T wμν 2 μν

+ (Uμν )T (wN + wels )].

(6)

Here U is a vector of the electrostatic fields at the dipole polarizability points and the electrostatic potential at the surface tesserae created by QM nuclei (UN ), HF density (Ue ),

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: 136.167.3.36 On: Sun, 02 Nov 2014 02:47:53

174115-3

Thellamurege et al.

J. Chem. Phys. 140, 174115 (2014)

and MM electrpstatic points (Uels ); the superscript T denotes transpose here and hereafter; w is a vector collecting all induced dipoles/charges, which satisfies the following linear response equation: B · w = U.

(7)

The definition of the B matrix can be found in a previous paper.8 Due to the linearity of Eq. (7), w in Eq. (6) can be determined separately for the nuclei (wN ,), HF electron density

(we ) of the QM region and the electrostatic potential (wels ) of the MM region. In addition, the induced dipoles/charges can be determined for the electrostatic field/potential due to a product of two basis functions such as μν, for example, wμν and Uμν in Eq. (6). Once EHF/MMpol/C is minimized via a self-consistent field procedure (in which the Fock operator contains induced dipole and induced surface charge terms), Pulay’s method44 can be used to determine its first derivative (nuclear gradient) with respect to a coordinate x,

x x x x x EHx F /MMpol/Con = EHx F + EMM + Erep + Edisp + Eels + Epolsol

=−



x Wμν Sμν +

μν

+ GxMM +





1  x Pμν Pρσ  μρ| |νσ x + EN,N 2 μνρσ   x x x + Pμν Vdisp,μν + Pμν Vels,μν + EN,els

Pμν hxμν +

μν x Pμν Vrep,μν

μν

μν

μν

1 1  − [(UN + Uels )T (wN + wels )]x − Pμν Pρσ [(Uμν )T wρσ ]x 2 2 μνρσ −

1 Pμν [(UN + Uels )T wμν + (Uμν )T (wN + wels )]x . 2 μν

Here P and W are the HF density matrix and energy-weighted density matrix, respectively. There is a general expression for the derivative of the dot product of an electric field/potential vector UK and an induced dipole/charge vector wL , ˜ K )T Bx wL + (w ˜ K )T UxL [(UK )T wL ]x = (UxK )T wL + (w

(9)

with T  ˜ K = B−1 UK . w

(10)

So the last three terms in Eq. (8) can be written as 1 − [(UN + Uels )T (wN + wels )]x 2 1  − Pμν Pρσ [(Uμν )T wρσ ]x 2 μνρσ −

1 Pμν [(UN + Uels )T wμν + (Uμν )T (wN + wels )]x 2 μν

1 ˜ T Bx w + (w) ˜ T Ux ]. = − [(Ux )T w + (w) 2

(11)

The first and third terms in the right side of Eq. (11) represent the forces and torques between the induced dipoles/charges and the QM nuclei, QM HF density, and the MM electrostatic points; the second term represents the forces and torques among the induced dipoles and induced charges. If the B ma˜ K = wK . In trix is symmetric, we have (B−1 )T = B−1 and w the EFP methods, the polarizabilities are asymmetric, so B is asymmetric.

(8)

B. QM/MMpol/C style MP2 gradient

In the QM/MMpol/C style MP2 method established in this work, QM/MMpol/C style HF orbitals and orbital energies are used to obtain the second order energy correction E(2) . Compared to Eq. (1), the total energy of a molecular system described with such a QM/MMpol/C style MP2 method has one more term, the second order energy correction E(2) , EMP 2/MMpol/C = EH F /MMpol/C + E (2) .

(12)

The HF orbitals, MM induced dipoles and induced surface charges are variationally determined in the presence of QM nuclear and MM electrostatic, repulsion and dispersion potentials.8 The MM electrostatic, repulsion, dispersion potentials, induced dipoles, and the induced surface charges in the continuum model implicitly affect E(2) through the HF orbitals and energies. For RMP2, UMP2, and ZAPT2 methods, the first derivative (gradient) of the second order energy correction E(2) with respect to a coordinate x in the QM/MMpol/C style MP2 method can be written using density matrices and basis function integrals,45, 46   (2) x (2) x E (2),x = − Wμν Sμν + Pμν hμν μν

+



μν (2) Pμν Pρσ  μρ| |νσ x +

μνρσ



NS  μρ| |νσ x μνρσ

μνρσ

   (2) x (2) x (2) x + Pμν Vrep,μν + Pμν Vdisp,μν + Pμν Vels,μν μν

μν

μν

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: 136.167.3.36 On: Sun, 02 Nov 2014 02:47:53

Thellamurege et al.

174115-4

− −

J. Chem. Phys. 140, 174115 (2014)

 x  1  (2)  T Pμν Uμν (wN + wels ) + UTels + UTN wμν 2 μν  T x 1   (2) (2) Pμν Pρσ + Pμν Pρσ Uμν wρσ . 2 μνρσ

(13)

Here W(2) is the MP2 correction to the HF energy-weighted density matrix; P(2) is the MP2 correction to the HF density matrix, which can be determined via a modified Z-vector NS is the nonseparable two-particle density mamethod, μνρσ trix. The first four terms in Eq. (13) have the same forms as those in the regular (gas phase) MP2 energy gradient formula,25 and can be evaluated using the same computer code after W(2) and P(2) are obtained. The 5th, 6th, 7th, and 8th terms in Eq. (13) are analogs of the nuclear term in hμν . The 9th term is an analog of the Coulomb term μρ||νσ . The 8th and 9th terms in Eq. (13) contain induced dipoles/charges and can be written as 1 1 (2) T x ˜ ) B w ˜ − (w − (U(2),x )T (w + w) 2 2 1 1 ˜ T Bx w(2) − (w(2) + w ˜ (2) )T Ux . − (w) 2 2

(14)

μν



(2) −1 Pμν B Uμν ,

(16)

(2) Pμν (B−1 )T Uμν .

(17)

μν

˜ (2) = w



The QM/MMpol/C style MP2 energy defined in Eq. (12) is the QM/MMpol/C style HF energy plus the second order energy correction E(2) . Using this formulation the analytic gradient can be derived. However, the polarization-solvation energy Epolsol in Eq. (12) is obtained at the HF level instead of MP2. Following our recent work,28, 29 here we describe a simple method to obtain the solvation energy using the MP2 relaxed density. The polarization-solvation energy at the HF level is defined as [see Eq. (6)], 1 HF = − (UN + Uels + Ue )T (wN + wels + we ) . (18) Epolsol 2 MP 2 Similar to Eq. (18), a polarization-solvation energy Epolsol can be defined using the MP2 relaxed density,

1 MP 2 = − (UN + Uels + Ue + U(2) )T Epolsol 2 × (wN + wels + we + w(2) )

Here U(2), x is the derivative of the electrostatic field/potential ˜ (2) are the dipoles/charges induced created by P(2) , w(2) , and w by the electrostatic field/potential created by P(2) ,  (2) x U(2),x = Pμν Uμν , (15)

w(2) =

D. Polarization-solvation energy correction

μν

The first term in Eq. (14) represents the forces and torques between w (induced by QM nuclei, QM HF density P and MM electrostatic points) and P(2) ; the second and third terms ˜ (2) ; the represent the forces and torques between w(2) and w fourth term represents the forces and torques between w(2) and QM nuclei, QM HF electron density P and MM electrostatic points.

C. Z-vector equations

By taking into account induced dipoles and induced charges, the Z-vector method originally proposed by Handy and Schaefer25 can be used to determine the P(2) for QM/ MMpol/C style MP2 methods. The induced dipole and induced charge contributions to the Z-vector equations are similar to the induced charge contributions in the MP2/continuum methods29, 47 and the induced dipole contributions in the QM/MMpol style MP2 methods.28 The equations for RMP2, UMP2, and ZAPT2 cases are presented in the Appendix.

1 = − (UN + Uels + Ue )T (wN + wels + we ) 2 1 − (UN + Uels + Ue )T (w(2) ) 2 1 1 − (U(2) )T (wN + wels + we ) − (U(2) )T (w(2) ). 2 2 (19) Here U(2) and w(2) are the field/potential and induced HF defined in Eq. (18) dipole/charge created by P(2) . The Epolsol for HF/MMpol/C calculation is the first term in Eq. (19). The second and third terms in Eq. (19) are equal to each other and can be combined. We choose to use the third term (doubled). Therefore, the correction of the polarization-solvation energy due to P(2) is (2) MP 2 HF Epolsol = Epolsol − Epolsol

1 = −(U(2) )T (wN + wels + we ) − (U(2) )T (w(2) ) 2   1 (2) (2) T = −(U ) wN + wels + we + w 2   1 = −(U(2) )T w + w(2) . (20) 2 Here w represents the induced dipoles and induced charges in the QM/MMpol/C style HF calculation. In the sense of perturbation theory, the correction of the polarization/solvation energy in Eq. (20) is actually at the fourth perturbation order. It is difficult to derive the expression of the analytic nuclear (2) MP 2 gradient for Epolsol or Epolsol . Numerical results in our recent (2) 28, 29 and the current work show that the Epolsol is almost work a constant for a molecule in a given electronic state, and us(2) correction is similar to applying a level-shift ing the Epolsol to the potential energy surface defined in Eq. (12). Therefore, Eq. (20) is best used as a correction to the energy expression in Eq. (12).

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: 136.167.3.36 On: Sun, 02 Nov 2014 02:47:53

174115-5

Thellamurege et al.

III. IMPLEMENTATION

We implemented two different versions of the QM/ MMpol/C style MP2 methods in the General Atomic and Molecular Electronic Structure System (GAMESS48, 49 ) package: (1) One is for the conductorlike polarizable continuum model (CPCM,50, 51 a variant of COSMO52 ) and the effective fragment potential (EFP)8, 10, 42, 43 method. (2) The other is for the FixSol53 solvation model and general induced dipole polarizable force field methods in the quantum chemistry polarizable force field (QuanPol54 ) program, which is designed for general molecular mechanics (MM) and combined quantum mechanics and molecular mechanics (QM/MM) calculations. The QM methods can be Hartree-Fock (HF), generalized valence bond theory (GVB55 ), multiconfiguration selfconsistent-field (MCSCF56–59 ), density functional theory (DFT60 ), time-dependent density functional theory (TDDFT61, 62 ) and second order Møller-Plesset perturbation theory (MP217 ) methods. In both implementations, the MP2 programs include the serial RMP2 and UMP2 programs,21 the parallel RMP2 program implemented by Fletcher et al.63 the parallel RMP2 program implemented by Ishimura, Pulay and Nagase,64, 65 the parallel ZAPT2 program implemented by Fletcher et al.26 and Aikens et al.,27 and the parallel UMP2 program implemented by Aikens et al.66, 67 Various tests show that the accuracy of the analytic gradients is typically better than 10−6 hartree/bohr. It is worth noting that various continuum solvation models can be used in the QM/MMpol/Continuum style methods. The only requirement is that the solvent reaction field be described with a set of linear equations (i.e., matrix equation) because we have only considered linear cases and related formulas. The CPCM and FixSol methods implemented here are based on Poisson equation, and are good for non-ionic solutions. For ionic solutions, in principle, continuum solvation models based on the linearized Poisson-Boltzmann equation should be used. It is straightforward to extend the formulas of our QM/MMpol/Continuum style methods to any continuum models that can be formulated in matrix equations. The same is true for induced dipole polarizable force fields.

IV. COMPUTATIONAL METHODS

QM/MMpol/C style MP2 geometry optimizations were performed for wild type photoactive yellow protein (PYP) and its Glu46Gln mutant. The X-ray structure file 2PHY36 was obtained from the Protein Data Bank (PDB68 ). We note that the 2PHY file does not have the highest resolution for PYP in PDB. All water molecules in the PDB file 2PHY were deleted. Hydrogen atoms were added to the X-ray structure using the WHAT IF web interface.69 The Glu46Gln mutant was created by manually editing the 2PHY file. The total number of atoms is 1929 for the wild type PYP and 1930 for the Glu46Gln mutant, including the chromophore pCA linked to Cys69. The total charge of the pro-

J. Chem. Phys. 140, 174115 (2014)

tein (including chromophore) is −6 e. All amino acids (1908 or 1909 atoms without the chromophore) were described with the Assisted Model Building with Energy Refinement (AMBER) induced dipole polarizable force field (topology and parameter files “all_amino02.in,” “all_aminoct02.in,” “all_aminont02.in,” and “parm99.dat”).70–72 As in standard AMBER force field, the charge-charge and LJ interactions between atoms separated by three covalent bonds were scaled by 1/1.2 and 0.50, respectively. The interactions between induced dipoles and the interactions between charges and induced dipoles were excluded for 1-2 and 1-3 atom pairs, and were included with no scaling for 1-4 atom pairs and beyond. In the QM/MMpol/C calculation, the QM region had 46 or 47 atoms: the side chain of Tyr42, the partial side chain of Glu46 or Gln46, and the chromophore (Figure 2). We used the capping H atom method implemented in QuanPol to treat the covalent bonds between QM and MM atoms.54 Two QM capping H atoms for Tyr42 and Cys69 were placed at the positions of their alpha carbon atoms, while the QM capping H atom for Glu46 and Gln46 was placed at the position of its beta carbon atom. In all QM/MMpol/C calculations the 6-31++G(d,p)73 basis set was used. In the QM/MM calculation, atomic charges of the QM atoms were not used because they interact with the MM atoms using their electrons and nuclear charges; the QM atoms use LJ potentials to interact with the MM atoms. The LJ parameters of the chromophore were obtained from the general AMBER force field (GAFF)74 by using the AmberTools1275 and have been published in an earlier paper.16 Atomic charges and the covalent terms were not used in the QM/MM calculation for the QM chromophore. In general, the LJ potentials parameterized together with atomic charge potentials for pure MM calculations are not transferable to QM/MM calculations in that QM atoms use nuclear charges and electron density to interact with MM atoms. A straightforward adoption of the LJ potential may lead to very inaccurate QM-MM interaction potentials. In the current work, the hydrogen bonding groups under consideration are all treated with QM methods, so the use the AMBER GAFF LJ potentials has only a secondary effect on the QMQM interaction potentials. The aqueous solvent was described with the FixSol model and 60 initial tesserae per sphere. In the FixSol calculations, atomic radii of 0.00, 2.10, 2.00, 1.90, and 2.40 Å were used for H, C, N, O, and S atoms, respectively, to define the molecular cavity (using zero for H atoms means that they do not contribute to form the surface). For all of the 1929 or 1930 QM and MM atoms in the calculation, FixSol used ∼16 500 surface tesserae and ∼8200 Å2 total solvent accessible surface area. A dielectric constant of 78.39 was used to represent aqueous solution at 298.15 K. The induced surface charges were determined by an iterative procedure with no charge renormalization.9, 76 The computational cost of adding the AMBER polarizable force field and the FixSol continuum model to MP2 calculation depends on different cases. We have run several single point energy and gradient calculations using 8 compute cores (two Intel Xeon E5430 CPUs in one node) running at 2.66 GHz. For the wild type PYP case, closed-shell Hartree-Fock calculation with the 6-31++G(d,p)

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: 136.167.3.36 On: Sun, 02 Nov 2014 02:47:53

174115-6

Thellamurege et al.

basis set for the 46 QM atoms took 9.7 min of wall clock time; the corresponding MP2 calculation took 120 min. When the whole protein was included as ∼1900 MM force field points, the MP2/AMBER(pol) calculation took 134 min. When the protein and aqueous solvent were included as ∼1900 force field points and 16 411 surface charge points, the MP2/AMBER(pol)/FixSol calculation took 229 min. In this case, the increase of the computational time is mainly due to the large number of surface charge points, which must be iteratively solved together with the iterative solution of the Z-vector equations (the Appendix) when calculating the MP2 gradient.

J. Chem. Phys. 140, 174115 (2014) TABLE I. Comparison of hydrogen bond distances (Å) in wild type and Glu46Gln photoactive yellow protein (PYP). Wild type PYP

OpCA –OTyr HTyr –OTyr OpCA –HTyr OpCA –OGlu (NGln ) HGlu –OGlu (NGln ) OpCA –HGlu (HGln )

Glu46Gln PYP

X-raya

X-rayb

Neutronc

MP2d

2.51 NA NA 2.58 NA NA

2.52 NA NA 2.56 NA NA

2.52 0.96 1.65 2.57 1.21 1.37

2.52 1.03 1.49 2.59 1.02 1.57

X-raye X-rayf MP2g 2.48 NA NA 2.86 NA NA

2.49 NA NA 2.87 NA NA

2.50 1.05 1.46 2.93 1.02 1.91

a

X-ray structure 1OTB with resolution of 1.10 Å, Ref. 79. X-ray structure 2ZOH with resolution of 1.25 Å, Ref. 77. Neutron diffraction measured structure 2OZI (no resolution is given), Ref. 77. d MP2/AMBER(pol)/FixSol geometry optimization from the X-ray structure 2PHY. e X-ray structure 1UGU with resolution of 1.20 Å, Ref. 80. f X-ray structure 1OTA with resolution of 1.10 Å, Ref. 79. g MP2/AMBER(pol)/FixSol geometry optimization from the X-ray structure 2PHY (resolution of 1.40 Å) that was manually edited to be the Glu46Gln mutant. b c

V. RESULTS AND DISCUSSION A. Hydrogen bonding distances

We performed MP2/AMBER(pol)/FixSol geometry optimization from the PDB structure in that Tyr42 and Glu46(Gln46) were neutral, pCA was anionic. All of the 1929 or 1930 QM and MM atoms were optimized. For typical small molecules, we usually converge the maximum and root mean square gradients to 1.0 × 10−4 and 0.33 × 10−4 hartree/bohr, respectively. It is difficult to optimize the geometry of a protein with a QM/MM/C method to such a tight convergence criterion. For the wild type case, we optimized for 593 steps, after which the maximum gradient and root mean square gradient are 3.4 × 10−4 and 0.29 × 10−4 hartree/bohr, respectively. If only the 46 QM atoms are considered, the maximum gradient and root mean square gradient are 3.4 × 10−4 and 0.44 × 10−4 hartree/bohr, respectively. For the Glu46Gln case, we optimized for 432 steps, after which the maximum gradient and root mean square gradient are 6.7 × 10−4 and 0.56 × 10−4 hartree/bohr, respectively. If only the 47 QM atoms are considered, the maximum gradient and root mean square gradient are 0.54 × 10−4 and 0.18 × 10−4 hartree/bohr, respectively. After optimization, the overall geometry of the protein is very similar to the original PDB geometry. In the wild type PYP, Tyr42 remained neutral, with the HTyr –OTyr distance optimized to 1.03 Å, the OpCA –HTyr distance optimized to 1.49 Å, and the OpCA –OTyr distance optimized to 2.52 Å. Glu46 remained neutral, with the HGlu –OGlu distance optimized to 1.02 Å, the OpCA –HGlu distance optimized to 1.56 Å, and the OpCA –OGlu distance optimized to 2.59 Å (Table I). The OpCA –OTyr distance is shorter than the OpCA –OGlu distance by 0.07 Å. In the Glu46Gln mutant, Tyr42 remained neutral, with the HTyr –OTyr distance optimized to 1.05 Å, the OpCA –HTyr distance optimized to 1.46 Å, and the OpCA –OTyr distance optimized to 2.50 Å. The HGln –NGln distance optimized to 1.02 Å, the OpCA –HGln distance optimized to 1.91 Å, and the OpCA –NGln distance optimized to 2.93 Å (Table I). The OpCA –OTyr distance is shorter than the OpCA –NGln distance by 0.43 Å. This difference in the hydrogen bond length suggests that the hydrogen bonding strength between the chromophore pCA and Gln46 is much weaker than that between pCA and Glu46 and that between pCA and Tyr42.

Based on neutron diffraction study of ground state PYP, Yamaguchi et al.77 assigned a HGlu –OGlu distance of 1.21 Å and a OpCA –HGlu distance of 1.37 Å for the pCA-Glu46 hydrogen bond, and a HTyr –OTyr distance of 0.96 Å and a OpCA – HTyr distance of 1.65 Å for the pCA-Tyr42 hydrogen bond (Table I). They proposed that the ground state chromophore phenolate forms a low-barrier hydrogen bond (LBHB) with Glu46 and a short ionic hydrogen bond (SIHB) with Tyr42. Saito and Ishikita78 performed QM/MM geometry optimization of PYP, and found a HGlu –OGlu distance of 1.00 Å and a OpCA –HGlu distance of 1.58 Å for the pCA-Glu46 hydrogen bond, and a HTyr –OTyr distance of 1.01 Å and a OpCA – HTyr distance of 1.50 Å for the pCA-Tyr42 hydrogen bond. Our recent DFT/AMBER(pol)/FixSol optimized distances16 are very similar to Saito and Ishikita’s results. Here the MP2/AMBER(pol)/FixSol results (Table I) are also different from the results of Yamaguchi et al.77 Although the 2PHY file we used to start the geometry optimization does not have a high resolution, the MP2/AMBER(pol)/FixSol optimized active site geometry is in good agreement with accurate X-ray diffraction results in PDB files 1OTB,79 2ZOH,77 1UGU,80 and 1OTA,79 as indicated by the hydrogen bond lengths (Table I). B. Hydrogen bonding strength

Using the MP2/AMBER(pol)/FixSol optimized geometry of the 46-atom or 47-atom QM region, we performed localized molecular orbital energy decomposition analysis (LMO-EDA)81 for the pCA-Tyr42 pair, the pCA-Glu46 pair, and the pCA-Gln46 pair (Figure 2, pure QM with no MM, and no FixSol). The LMO-EDA calculation was performed at the MP2/aug-cc-pVTZ82 level of theory with the counterpoise method proposed by Boys and Bernardi83 for correcting the basis set superposition errors. Since we used the aug-cc-pVTZ basis set (not a very large basis set) and the counterpoise method, the derived MP2 interaction energies are typically underestimated by ∼1 kcal/mol. These interaction energies are two-body values for the two involved

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: 136.167.3.36 On: Sun, 02 Nov 2014 02:47:53

174115-7

Thellamurege et al.

J. Chem. Phys. 140, 174115 (2014)

TABLE II. Intermolecular interaction energy (kcal/mol) of the pCA-Tyr42 pair, the pCA-Glu46 pair, and the pCA-Gln46 pair in wild type photoactive yellow protein and its Glu46Gln mutant. Protein Wild type Glu46Gln Wild type Glu46Gln

Pair

Eele

Eex

pCA-Tyr42 pCA-Tyr42 pCA-Glu46 pCA-Gln46

−37.71 −39.08 −39.53 −21.66

−50.18 −53.08 −37.36 −17.23

Erep

Epol

bond is weaker than the pCA-Tyr42 hydrogen bond by ∼8 kcal/mol (Table II). The pCA-Tyr42 hydrogen bond shows similar strengths in the wild type and Glu46Gln mutant.

Edisp EMP2

+92.27 − 26.81 −4.35 −26.77 +98.16 − 29.02 −4.16 −27.17 +68.51 − 21.63 −0.96 −30.97 +29.92 − 8.61 −1.62 −19.20

groups. When other protein groups are involved, the total interaction energy tends to be smaller in magnitude than the sum of the two-body values, due to many-body effects. Nevertheless, the two-body interaction energies are still very meaningful for comparing the relative interaction strengths. These interactions contain mainly hydrogen bonding, and are general indicators of the hydrogen bonding strengths. In the wild type PYP, pCA shows a 4.20 kcal/mol stronger interaction with Glu46 (−30.97 kcal/mol) than with Tyr42 (−26.77 kcal/mol). In the Glu46Gln mutant, pCA shows a 7.97 kcal/mol weaker interaction with Gln46 (−19.20 kcal/mol) than with Tyr42 (−27.17 kcal/mol). Therefore, the pCA-Gln46 hydrogen bond is ∼12 kcal/mol weaker than the pCA-Glu46 one. The difference can be understood from the LMO-EDA analysis (Table II). The leading term is the electrostatic interaction. The Glu46 has a highly polar O-H group, which is attracted by the anionic oxygen of pCA to form a short and strong hydrogen bond, with an electrostatic energy of −39.53 kcal/mol. Due to the short bond length and significant molecular orbital overlap, the exchange, repulsion and polarization energies are all significant. The Gln46 has a less polar N-H group, which forms a long and weak hydrogen bond with pCA, with an electrostatic energy of −21.66 kcal/mol. Due to the longer bond length and less molecular orbital overlap, the exchange, repulsion and polarization energies are all smaller in magnitude than those in the Glu46 case. VI. CONCLUSION

The analytic gradient of the QM/MMpol/C style MP2 methods is derived and implemented for the Effective Fragment Potential (EFP) method and the QuanPol method. Using the analytic gradient, geometry optimization can be performed. Using MP2/6-31++G(d,p) to describe the active site (pCA, Tyr42, Glu46, or Gln46), AMBER induced dipole polarizable force field to describe the rest of the protein, and the FixSol solvation model to describe the aqueous solvent, geometry optimizations were performed for all of the atoms of the ground state of photoactive yellow protein (PYP) and its Glu46Gln mutant. The QM/MMpol/C style MP2 geometry optimization suggests that in both the wild type and the Glu46Gln mutant the chromophore pCA is anionic, forming two normal hydrogen bonds with neutral Tyr42 and neutral Glu46 (or Gln46). LMO-EDA calculations at the MP2/augcc-pVTZ level of theory reveal that the pCA-Glu46 hydrogen bond is stronger than the pCA-Tyr42 hydrogen bond by ∼4 kcal/mol (Table II), and that the pCA-Gln46 hydrogen

ACKNOWLEDGMENTS

This work was supported by the USA National Science Foundation (NSF Award No. 1010674). APPENDIX: Z-VECTOR EQUATIONS WITH INDUCED DIPOLE AND INDUCED CHARGE CONTRIBUTIONS

In the following, the doubly occupied molecular orbitals are denoted as i, j, and k; the singly occupied ones are denoted as x and y; the virtual ones are denoted as a, b, and c; general molecular orbitals (both occupied and virtual) are denoted as p and q. For RMP2, the Z-vector equation can be written as  Aaibj Pbj(2) + (εa − εi )Pai(2) = Lai . (A1) bj

Here Aaibj = ij||ab + ib||aj is the orbital Hessian matrix, εa and εi are the orbital energies, and L is the Lagrangian.25 For QM/MMpol/C style RMP2, induced dipole and induced charge terms should be added to the orbital Hessian matrix A and the Lagrangian L of the Z-vector equation, and the occupied-occupied block of the W(2) matrix, polsol

Aaibj polsol

Lai

˜ bj ), = Aaibj − UTai (wbj + w

= Lai − −

(2),polsol

Wij

1  (2) T ˜ jk) P U (wj k + w 2 j k j k ai

1  (2) T ˜ bc ), P U (wbc + w 2 bc bc ai

= Wij(2) +

(A2)

1  (2) T ˜ pq ). P U (wpq + w 4 pq pq ij

(A3)

(A4)

For UMP2, two associated Z-vector equations for the α and β orbitals should be solved. For QM/MMpol/C style UMP2, induced dipole and induced charge terms should be added to the orbital Hessian matrix A, the Lagrangian L, and the occupied-occupied block of the W(2) matrix, polsol

(A5)

polsol

(A6)

˜ bα j α ), Aa α i α bα j α = Aa α i α bα j α − UTaα i α (wbα j α + w ˜ bβ j β ), Aa α i α bβ j β = Aa α i α bβ j β − UTaα i α (wbβ j β + w polsol

La α i α

= La α i α − −





T ˜ j α kα ) Pj(2) α k α Ua α i α (wj α k α + w

j α kα T ˜ j β kβ ) Pj(2) β k β Ua α i α (wj β k β + w

j β kβ

−



T ˜ bα cα ) Pb(2) α cα Ua α i α (wbα cα + w

bα cα

−



T ˜ bβ cβ ), Pb(2) β cβ Ua α i α (wbβ cβ + w

(A7)

bβ 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: 136.167.3.36 On: Sun, 02 Nov 2014 02:47:53

Thellamurege et al.

174115-8

(2),polsol Wi α j α

=

Wi(2) αjα +

J. Chem. Phys. 140, 174115 (2014)

1  (2) T ˜ pα q α ) + P α α U α α (wpα q α + w 2 pα q α p q i j

1  (2) T ˜ pβ q β ). P β β U α α (wpβ q β + w 2 β β pq ij

= Lax − Lpolsol ax −

(A8)

For QM/MMpol/C style ZAPT2, induced dipole and induced charge terms should be added to the orbital Hessian matrix A, the Lagrangian L and the double-double, double-single and single-single blocks of the W(2) matrix, polsol

Axiyj

1 ˜ yj ), = Axiyj − UTxi (wyj + w 2

(A10)

T ˜ jk) Pj(2) k Uax (wj k + w

jk (2) T ˜ xy ) Pxy Uax (wxy + w

xy

p q

Only the α equations are given because the β equations can be obtained by switching α and β. For ZAPT2, a single Z-vector equation can be constructed using a nine-block (only six independent blocks) orbital Hessian matrix A for doubly occupied, singly occupied and virtual orbitals,27 ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ Axi,yj Axi,bj Axi,by Pyj Lxi ⎝ Aai,yj Aai,bj Aai,by ⎠ × ⎝ Pbj ⎠ = ⎝ Lai ⎠ . (A9) Aax,yj Aax,bj Aax,by Pby Lax





−



(2) T ˜ bc ), Pbc Uax (wbc + w

(A18)

bc (2),polsol

Wij

= Wij(2) +



(2) T ˜ pq ), Ppq Uij (wpq + w

(A19)

pq

= Wix(2) +

1  (2) T ˜ pq ), P U (wpq + w 2 pq pq ix

(A20)

(2),polsol (2) = Wxy + Wxy

1  (2) T ˜ pq ). P U (wpq + w 2 pq pq xy

(A21)

(2),polsol

Wix

The expressions of A, L, and W(2) for gas phase UMP2 and ZAPT2 can be found in the literature.27, 66 In adiition, corresponding modifications can be made for frozen core cases.27, 66 1 A.

polsol Aaiyj

= Aaiyj −

UTai (wyj

˜ yj ), +w

1 polsol ˜ yj ), Aaxyj = Aaxyj − UTax (wyj + w 2 polsol

(A11)

(A12)

˜ bj ), = Aaibj − 2UTai (wbj + w

(A13)

˜ bj ), Aaxbj = Aaxbj − UTax (wbj + w

(A14)

1 polsol ˜ by ), Aaxby = Aaxby − UTax (wby + w 2

(A15)

Aaibj

polsol

polsol

Lxi

= Lxi − −





T ˜ jk) Pj(2) k Uxi (wj k + w

jk (2) T ˜ xy ) Pxy Uxi (wxy + w

xy

−



(2) T ˜ bc ), Pbc Uxi (wbc + w

(A16)

bc

polsol

Lai

= Lai − 2 −2





T ˜ jk) Pj(2) k Uai (wj k + w

jk (2) T ˜ xy ) Pxy Uai (wxy + w

xy

−2

 bc

(2) T ˜ bc ), Pbc Uai (wbc + w

(A17)

H. Devries, P. T. van Duijnen, A. H. Juffer, J. A. C. Rullmann, J. P. Dijkman, H. Merenga, and B. T. Thole, J. Comput. Chem. 16, 37 (1995). 2 P. Bandyopadhyay and M. S. Gordon, J. Chem. Phys. 113, 1104 (2000). 3 A. Ohrn and G. Karlstrom, Mol. Phys. 104, 3087 (2006). 4 Q. Cui, J. Chem. Phys. 117, 4720 (2002). 5 P. Bandyopadhyay, M. S. Gordon, B. Mennucci, and J. Tomasi, J. Chem. Phys. 116, 5023 (2002). 6 S. Chalmet, D. Rinaldi, and M. F. Ruiz-López, Int. J. Quantum Chem. 84, 559 (2001). 7 A. H. Steindal, K. Ruud, L. Frediani, K. Aidas, and J. Kongsted, J. Phys. Chem. B 115, 3027 (2011). 8 H. Li, J. Chem. Phys. 131, 184103 (2009). 9 H. Li, C. S. Pomelli, and J. H. Jensen, Theor. Chem. Acc. 109, 71 (2003). 10 H. Li and M. S. Gordon, J. Chem. Phys. 126, 124112 (2007). 11 F. J. Vesely, J. Comput. Phys. 24, 361 (1977). 12 M. Neumann, F. J. Vesely, O. Steinhauser, and P. Schuster, Mol. Phys. 35, 841 (1978). 13 F. H. Stillinger and C. W. David, J. Chem. Phys. 69, 1473 (1978). 14 P. Barnes, J. L. Finney, J. D. Nicholas, and J. E. Quinn, Nature (London) 282, 459 (1979). 15 A. Warshel, J. Phys. Chem. 83, 1640 (1979). 16 N. M. Thellamurege, F. Cui, and H. Li, J. Chem. Phys. 139, 084106 (2013). 17 C. Møller and M. S. Plesset, Phys. Rev. 46, 618 (1934). 18 J. A. Pople, J. S. Binkley, and R. Seeger, Int. J. Quantum Chem. 10, 1 (1976). 19 T. J. Lee and D. Jayatilaka, Chem. Phys. Lett. 201, 1 (1993). 20 J. S. Binkley and J. A. Pople, Int. J. Quantum Chem. 9, 229 (1975). 21 M. J. Frisch, M. Head-Gordon, and J. A. Pople, Chem. Phys. Lett. 166, 275 (1990). 22 T. J. Lee, A. P. Rendell, K. G. Dyall, and D. Jayatilaka, J. Chem. Phys. 100, 7400 (1994). 23 I. M. B. Nielsen and E. T. Seidl, J. Comput. Chem. 16, 1301 (1995). 24 J. A. Pople, R. Krishnan, H. B. Schlegel, and J. S. Binkley, Int. J. Quantum Chem. 16, 225 (1979). 25 N. C. Handy and H. F. Schaefer III, J. Chem. Phys. 81, 5031 (1984). 26 G. D. Fletcher, M. S. Gordon, and R. S. Bell, Theor. Chem. Acc. 107, 57 (2002). 27 C. M. Aikens, G. D. Fletcher, M. W. Schmidt, and M. S. Gordon, J. Chem. Phys. 124, 014107 (2006). 28 H. Li, J. Phys. Chem. A 115, 11824 (2011). 29 D. J. Si and H. Li, J. Chem. Phys. 135, 144107 (2011). 30 T. E. Meyer, Biochim. Biophys. Acta 806, 175 (1985). 31 U. K. Genick, S. M. Soltis, P. Kuhn, I. L. Canestrelli, and E. D. Getzoff, Nature (London) 392, 206 (1998).

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: 136.167.3.36 On: Sun, 02 Nov 2014 02:47:53

174115-9 32 A.

Thellamurege et al.

Xie, W. D. Hoff, A. R. Kroon, and K. J. Hellingwerf, Biochemistry 35, 14671 (1996). 33 W. D. Hoff, A. Xie, I. H. M. Van Stokkum, X.-j. Tang, J. Gural, A. R. Kroon, and K. J. Hellingwerf, Biochemistry 33, 1009 (1999). 34 E. D. Getzoff, K. N. Gutwin, and U. K. Genick, Nat. Struct. Mol. Biol. 5, 568 (1998). 35 M. Baca, G. E. O. Borgstahl, M. Boissinot, P. M. Burke, D. R. Williams, K. A. Slater, and E. D. Getzoff, Biochemistry 33, 14369 (1994). 36 G. E. O. Borgstahl, D. R. Williams, and E. D. Getzoff, Biochemistry 34, 6278 (1995). 37 M. Kim, R. A. Mathies, W. D. Hoff, and K. J. Hellingwerf, Biochemistry 34, 12669 (1995). 38 E. Demchuk, U. K. Genick, T. T. Woo, E. D. Getzoff, and D. Bashford, Biochemistry 39, 1100 (2000). 39 M. Yoda, Y. Inoue, and M. Sakurai, J. Phys. Chem. B 107, 14569 (2003). 40 A. R. Kroon, W. D. Hoff, H. P. M. Fennema, J. Gijzen, G.-J. Koomen, J. W. Verhoeven, W. Crielaard, and K. J. Hellingwerf, J. Biol. Chem. 271, 31949 (1996). 41 J. J. Vanbeeumen, B. V. Devreese, S. M. Vanbun, W. D. Hoff, K. J. Hellingwerf, T. E. Meyer, D. E. McRee, and M. A. Cusanovich, Protein Sci. 2, 1114 (1993). 42 P. N. Day, J. H. Jensen, M. S. Gordon, S. P. Webb, W. J. Stevens, M. Krauss, D. Garmer, H. Basch, and D. Cohen, J. Chem. Phys. 105, 1968 (1996). 43 H. Li, H. M. Netzloff, and M. S. Gordon, J. Chem. Phys. 125, 194103 (2006). 44 P. Pulay, Mol. Phys. 17, 197 (1969). 45 J. E. Rice and R. D. Amos, Chem. Phys. Lett. 122, 585 (1985). 46 E. D. Simandiras, R. D. Amos, and N. C. Handy, Chem. Phys. 114, 9 (1987). 47 R. Cammi, B. Mennucci, and J. Tomasi, J. Phys. Chem. A 103, 9100 (1999). 48 M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. J. Su, T. L. Windus, M. Dupuis, and J. A. Montgomery, J. Comput. Chem. 14, 1347 (1993). 49 M. S. Gordon and M. W. Schmidt, in Theory and Applications of Computational Chemistry, edited by C. E. Dykstra, G. Frenking, K. S. Kim, and G. E. Scuseria, (Elsevier, Amsterdam, 2005). 50 H. Li and J. H. Jensen, J. Comput. Chem. 25, 1449 (2004). 51 V. Barone and M. Cossi, J. Phys. Chem. A 102, 1995 (1998). 52 A. Klamt and G. Schuurmann, J. Chem. Soc., Perkin Trans. 2, 799 (1993). 53 N. M. Thellamurege and H. Li, J. Chem. Phys. 137, 246101 (2012). 54 N. M. Thellamurege, D. Si, F. Cui, H. Zhu, R. Lai, and H. Li, J. Comput. Chem. 34, 2816 (2013). 55 W. A. Goddard, T. H. Dunning, W. J. Hunt, and P. J. Hay, Acc. Chem. Res. 6, 368 (1973). 56 P. E. M. Siegbahn, J. Almlof, A. Heiberg, and B. O. Roos, J. Chem. Phys. 74, 2384 (1981). 57 J. Ivanic and K. Ruedenberg, J. Comput. Chem. 24, 1250 (2003). 58 J. Ivanic and K. Ruedenberg, Theor. Chem. Acc. 106, 339 (2001).

J. Chem. Phys. 140, 174115 (2014) 59 J.

M. Rintelman, M. S. Gordon, G. D. Fletcher, and J. Ivanic, J. Chem. Phys. 124, 034303 (2006). 60 W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 61 M. E. Casida, in Recent Advances in Density Functional Methods, edited by D. P. Chong (World Scientific, Singapore, 1995), p. 155. 62 M. E. Casida, C. Jamorski, K. C. Casida, and D. R. Salahub, J. Chem. Phys. 108, 4439 (1998). 63 G. D. Fletcher, M. W. Schmidt, and M. S. Gordon, in Advances in Chemical Physics, Vol. 110 (John Wiley & Sons, Inc., 1999), p. 267. 64 K. Ishimura, P. Pulay, and S. Nagase, J. Comput. Chem. 27, 407 (2006). 65 K. Ishimura, P. Pulay, and S. Nagase, J. Comput. Chem. 28, 2034 (2007). 66 C. M. Aikens, S. P. Webb, R. L. Bell, G. D. Fletcher, M. W. Schmidt, and M. S. Gordon, Theor. Chem. Acc. 110, 233 (2003). 67 C. M. Aikens and M. S. Gordon, J. Phys. Chem. A 108, 3103 (2004). 68 H. M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T. N. Bhat, H. Weissig, I. N. Shindyalov, and P. E. Bourne, Nucl. Acids Res. 28, 235 (2000). 69 R. Rodriguez, G. Chinea, N. Lopez, T. Pons, and G. Vriend, Bioinformatics 14, 523 (1998). 70 P. Cieplak, J. Caldwell, and P. Kollman, J. Comput. Chem. 22, 1048 (2001). 71 D. A. Case, T. E. Cheatham, T. Darden, H. Gohlke, R. Luo, K. M. Merz, A. Onufriev, C. Simmerling, B. Wang, and R. J. Woods, J. Comput. Chem. 26, 1668 (2005). 72 Z. X. Wang, W. Zhang, C. Wu, H. X. Lei, P. Cieplak, and Y. Duan, J. Comput. Chem. 27, 781 (2006). 73 M. M. Francl, W. J. Pietro, W. J. Hehre, J. S. Binkley, M. S. Gordon, D. J. Defrees, and J. A. Pople, J. Chem. Phys. 77, 3654 (1982). 74 J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, and D. A. Case, J. Comput. Chem. 25, 1157 (2004). 75 D. A. Case, T. A. Darden, T. E. Cheatham III, C. L. Simmerling, J. Wang, R. E. Duke, R. Luo, R. C. Walker, W. Zhang, K. M. Merz, B. Roberts, S. Hayik, A. Roitberg, G. Seabra, J. Swails, A. W. Goetz, I. Kolossváry, K. F. Wong, F. Paesani, J. Vanicek, R. M. Wolf, J. Liu, X. Wu, S. R. Brozell, T. Steinbrecher, H. Gohlke, Q. Cai, X. Ye, J. Wang, M. J. Hsieh, G. Cui, D. R. Roe, D. H. Mathews, M. G. Seetin, R. Salomon-Ferrer, C. Sagui, V. Babin, T. Luchko, S. Gusarov, A. Kovalenko, and P. A. Kollman, 12th ed. (University of California, San Francisco, 2012). 76 C. S. Pomelli, J. Tomasi, and V. Barone, Theor. Chem. Acc. 105, 446 (2001). 77 S. Yamaguchi, H. Kamikubo, K. Kurihara, R. Kuroki, N. Niimura, N. Shimizu, Y. Yamazaki, and M. Kataoka, Proc. Natl. Acad. Sci. U.S.A. 106, 440 (2009). 78 K. Saito and H. Ishikita, Proc. Natl. Acad. Sci. U.S.A. 109, 167 (2011). 79 S. Anderson, S. Crosson, and K. Moffat, Acta Crystallogr., Sect. D 60, 1008 (2004). 80 M. Sugishima, N. Tanimoto, K. Soda, N. Hamada, F. Tokunaga, and K. Fukuyama, Acta Crystallogr., Sect. D 60, 2305 (2004). 81 P. Su and H. Li, J. Chem. Phys. 131, 014102 (2009). 82 T. H. Dunning, J. Chem. Phys. 90, 1007 (1989). 83 S. F. Boys and F. Bernardi, Mol. Phys. 19, 553 (1970).

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: 136.167.3.36 On: Sun, 02 Nov 2014 02:47:53