Benchmark Interaction Energies for Biologically Relevant Noncovalent

Dec 19, 2011 - ... Daniel M. Chipman , Christopher J. Cramer , William A. Goddard , Mark S. Gordon , Warren J. Hehre , Andreas Klamt , Henry F. Schaef...
1 downloads 13 Views 880KB Size
ARTICLE pubs.acs.org/JPCA

Benchmark Interaction Energies for Biologically Relevant Noncovalent Complexes Containing Divalent Sulfur Benjamin J. Mintz† and Jerry M. Parks*,‡ †

Oak Ridge Leadership Computing Facility, Oak Ridge National Laboratory, 1 Bethel Valley Road, Oak Ridge, Tennessee 37831-6164, United States ‡ UT/ORNL Center for Molecular Biophysics, Oak Ridge National Laboratory, P.O. Box 2008, Oak Ridge, Tennessee 37831-6309, United States

bS Supporting Information ABSTRACT: Molecules containing divalent sulfur can participate in significant noncovalent interactions. Computing accurate noncovalent interaction energies using ab initio quantum chemical methods requires a proper description of electron correlation effects. Coupled-cluster theory with single and double substitutions and perturbative triple substitutions, CCSD(T), using extrapolation to the complete basis set (CBS) limit has become the method of choice for computing accurate interaction energies of noncovalently bound complexes. Here, interaction energies are computed for several biologically relevant hydrogen-bonded and dispersion-bound complexes that contain divalent sulfur. Eight-point estimated CCSD(T)/CBS dissociation curves along the noncovalent interaction vector are computed for each complex. As a comparison of high-accuracy ab initio methods, interaction energies are also calculated for each complex using the correlation-consistent Composite Approach (ccCA). We find that, on average, the two methods yield energies within 0.1 kcal mol1 of each other. The interaction energies provided here should be useful for developing and assessing the accuracy of more approximate ab initio, density functional theory, semiempirical, and classical force field approaches.

’ INTRODUCTION Noncovalent interactions play a critical role in defining the structure and function of biological systems, but describing these interactions accurately using quantum chemistry can be quite challenging. To aid in developing and assessing the accuracy of various theoretical methods for computing interaction energies, high-accuracy benchmark data sets are essential. The S22 set of noncovalent interaction energies developed by Hobza and co-workers1 was created as a benchmark database for developing and assessing the accuracy of more approximate ab initio, density functional theory, semiempirical, and classical force field methods. In the original work, interaction energies were computed for 22 complexes of biological relevance using coupled-cluster theory including single, double, and perturbative triple substitutions, CCSD(T),2 using extrapolations to the complete, one-particle basis set (CBS) limit. The S22 interaction energies were computed using the following formulation: þ ΔCCSDðTÞ Eint ¼ EHF=CBS þ ΔEMP2=CBS corr

ð1Þ

where EHF/CBS and ΔEMP2/CBS are the CBS-extrapolated corr HartreeFock (HF) reference and MP2 correlation energies, respectively, and ΔCCSD(T) is the difference between the CCSD(T) and MP2 correlation energies computed using a relatively small basis set. The S22 set was later extended by r 2011 American Chemical Society

Riley and Hobza to include four additional hydrogen-bonded complexes, which they named S26-07.3 Grafova et al. computed five-point CCSD(T)/CBS dissociation curves along the noncovalent interaction coordinate for each complex in the S22 set.4 In their approach, the interaction coordinate was scaled to generate complexes with one shortened and three elongated intermonomer distances. Merz and co-workers computed 14-point dissociation curves and derived analytic Morse curves for 20 of the 22 complexes in the S22 set.5 The approach for computing interaction energies using eq 1 is based on the assumption that the MP2 energies contain most of the correlation effects, and that the ΔCCSD(T) correction is small and converges faster with respect to basis set size than .6 Whereas this approximation worked well for the ΔEMP2/CBS corr smaller complexes in the S22 set, where relatively large basis sets (i.e., up to cc-pVQZ) were used to compute ΔCCSD(T), it broke down for the larger complexes, which necessitated the use of smaller basis sets (i.e., as small as cc-pVDZ). Podeszwa et al. improved the S22 interaction energy estimates by using basis sets as large as aug-cc-pVTZ with an additional set of 3s3p2d2f midbond functions when computing the ΔCCSD(T) term.6 Received: October 3, 2011 Revised: December 15, 2011 Published: December 19, 2011 1086

dx.doi.org/10.1021/jp209536e | J. Phys. Chem. A 2012, 116, 1086–1092

The Journal of Physical Chemistry A

ARTICLE

Figure 1. Representations of the MP2/TZ-CP optimized geometries of 10 hydrogen bonded complexes containing sulfur: (a) water dimer analogues, (b) cysteine analogues, and (c) methionine analogues. Complexes do not necessarily correspond to global minima. Cartesian coordinates for all complexes are provided as Supporting Information. Figure rendered using XYZViewer.47

A second issue of the original S22 study was the lack of consistency in the chosen basis sets. For example, the ΔEMP2/CBS corr term for H2O dimer was computed by extrapolating energies obtained using cc-pVQZ and cc-pV5Z basis sets, but cc-pVTZ and cc-pVQZ basis sets were used for uracil dimer. Podeszwa et al. also used basis sets of varying size to compute the ΔCCSD(T) terms.6 Sherrill and co-workers applied larger basis sets consistently across the entire set,7 and also modified eq 1 by extrapolating the MP2 and CCSD(T) correlation energies in the ΔCCSD(T) term using aug-cc-pVDZ and aug-cc-pVTZ basis sets.7 A systematic study of the convergence behavior of the ΔCCSD(T) term subsequently revealed that this term does not converge monotonically when smaller basis sets are used, and extrapolation using aug-cc-pVTZ and aug-cc-pVQZ basis sets or better may be required to achieve convergence.8  ezac et al. extended the S22 set to include 66 benchmark R noncovalent interaction energies (the S66 set) and eight-point dissociation curves along the interaction coordinate of each complex (the S66  8 set).9 A major emphasis of these sets is on the

coverage of the most important types of noncovalent interactions in biomolecular systems. Complexes dominated by hydrogen bonding, dispersion, and mixed interactions are included. Consistent basis sets were used throughout, and aug-cc-pVDZ basis sets were used to compute the ΔCCSD(T) terms. The S66 set was subsequently extended10 to include interaction energies in which the ΔCCSD(T) term was extrapolated as in ref 7, and additional angular-displaced, nonequilibrium geometries were also included. The S22 and S66 test suites are chemically diverse, but neither set includes molecules containing sulfur. As compared to their oxygen and nitrogen counterparts, relatively few theoretical studies have been carried out for sulfur-containing hydrogen-bonded systems. Nevertheless, noncovalent sulfur interactions have continued to gain attention in theoretical studies,6,7,1124 as they are known to be important in biological systems. For example, statistical analyses of structures in crystallographic databases have shown that sulfur-containing residues in proteins participate in significant hydrogen-bonding and dispersion interactions.13,17,2529 1087

dx.doi.org/10.1021/jp209536e |J. Phys. Chem. A 2012, 116, 1086–1092

The Journal of Physical Chemistry A Ab initio calculations have shown that SHO hydrogen bonds are nearly as strong as their OHO analogues. For example, Wennmohs et al. computed CCSD(T) interaction energies and potential energy surfaces for the dimethylsulfide methanol and dimethyl ethermethanol complexes using basis sets as large as aug-cc-pVQZ.15 They found that the interaction energy of the dimethylsulfidemethanol complex was 5.46 kcal mol1, only 0.51 kcal mol1 weaker than the oxygen analogue dimethylethermethanol (5.97 kcal mol1), and that dispersion effects contribute significantly to the strongly attractive interactions. Sherrill and co-workers have used CCSD(T) calculations to study S/π interactions in the benzeneH2S complex.16,17,21,30 In the energetically most favorable orientation of the complex, which was computed to have a CBS-extrapolated interaction energy of 2.81 kcal mol1,16 the sulfur atom of H2S is located above the center of the aromatic ring, and the two hydrogens are oriented toward the ring. On the basis of symmetry-adapted perturbation theory31 (SAPT) energy decomposition analysis, they attributed the intermolecular attraction to the electrostatic interaction between the partial positive charges on the hydrogens in H2S and the electron-rich π-cloud in benzene. Later, dissociation curves for benzeneH2S were computed using various CCSD(T)/CBS estimates, with a best estimate of 2.83 kcal mol1 for the equilibrium interaction energy.30 In another work,17 additional configurations of the benzeneH2S complex were considered, the most energetically favorable of which were found to correlate with geometries of cysteine and the aromatic side chains of phenylalanine, tyrosine, and tryptophan observed in protein X-ray crystal structures. Here, following the S66  8 procedure of Hobza and coworkers,9 and revisions suggested by Sherrill and co-workers,8 benchmark-quality CCSD(T)/CBS-estimated noncovalent interaction energies and dissociation curves are computed for 10 hydrogen-bonded complexes (Figure 1) and three dispersion or mixed-character complexes (Figure 2). For the hydrogenbonded complexes, we adopt the convention that the donor is always listed first. The hydrogen-bonded complexes are subdivided into three subsets: water dimer analogues, cysteine analogues, and methionine analogues. Water dimer analogues have the same topology as the Cs-symmetric water dimer but contain one or two H2S molecules in place of H2O. These include H2OH2S, H2SH2O, and H2S dimer. Cysteine analogues, in which CH3SH is used to represent the side chain of the amino acid cysteine, include CH3SHH2O, CH3 SHNH3, CH3SHH2CO, CH3SH dimer, and formamideCH3SH. For the formamideCH3SH complex, the hydrogen from the NH bond syn to the carbonyl in formamide is the hydrogen-bond donor. Methionine analogues, in which CH3SCH3 represents the thioether group in the side chain of the amino acid methionine, include H2OCH3SCH3 and CH3OHCH3SCH3. Two S/π complexes included in the set are CH3SHbenzene, which may be classified as a mixed complex that includes both SH/π hydrogen bonding and dispersive character, and CH3SCH3 benzene, which is dispersion-dominated. CH3SHbenzene is representative of cysteine interacting with either phenylalanine, tyrosine, or tryptophan, and CH3SCH3benzene is a model of methionine interacting with these aromatic residues. An additional dispersion complex, H2SCH4, was also included in the set. For consistency and compatibility with the S66  8 set, we  ezac et al.9 for report results obtained using the procedures of R

ARTICLE

Figure 2. Representations of the MP2/TZ-CP optimized geometries of H2SCH4 and two S/π complexes containing sulfur. Complexes do not necessarily correspond to global minima. Cartesian coordinates for the complexes are provided as Supporting Information. Figure rendered using XYZViewer.47

estimating the CCSD(T)/CBS interaction energies along eightpoint dissociation curves. We also adopt the strategy developed by Takatani et al.7 for extrapolating the ΔCCSD(T) correction, which is discussed further in the Methods. However, it has been shown that the ΔCCSD(T) correction term may not begin to converge monotonically until basis sets of at least aug-cc-pVTZ quality or better are used, particularly for hydrogen-bonded complexes.8 Thus, we examine the convergence behavior of the ΔCCSD(T) term for the three smallest complexes. As a comparison of high-accuracy ab initio techniques, we also use the correlation-consistent Composite Approach3234 (ccCA) to compute interaction energies for each complex. The ccCA method is a thoroughly benchmarked composite approach that has been shown to provide results within chemical accuracy, 1 kcal mol1 of reliable experimental data.3537

’ METHODS The correlation-consistent basis sets (cc-pVnZ and aug-ccpVnZ, where n = D, T, Q) developed by Dunning and co-workers were used throughout this study.38,39 For sulfur, the cc-pV(n+d)Z and aug-cc-pV(n+d)Z basis sets40 were used, which include a tight-d function to account for defects in the original correlation-consistent basis sets for sulfur and are essential for describing sulfur-containing species correctly.41,42 Basis sets are abbreviated as nZ = cc-pVnZ or anZ = aug-cc-pVnZ, where n = D, T, or Q, and extrapolations are abbreviated as a((n  1)(n))Z. Geometry optimizations were performed at the counterpoisecorrected43 MP2/TZ level, which we abbreviate as MP2/TZ-CP. Deformation energy corrections for the monomers were neglected. In some cases, symmetry was used to expedite the calculations. Estimates of the CCSD(T)/CBS interaction energies were computed as follows. The HartreeFock (HF) self-consistent field (SCF) reference energies were computed using aQZ basis sets without extrapolation, as it was noted that the HF/aQZ results were essentially converged to the CBS limit for the S22A set.7 MP2/CBS correlation energies were approximated using a two-point X3 extrapolation44 of the aTZ and aQZ correlation energies, which we denote as MP2/CBS(a(TQ)Z). The ΔCCSD(T) correction term was calculated by taking the difference between the CCSD(T) and MP2 correlation energies using either aDZ, aTZ, or extrapolated a(DT)Z energies. This approach gives rise to three different CCSD(T)/CBS estimates, which are abbreviated as CCSD(T)/CBS(ΔaDZ), CCSD(T)/ CBS(ΔaTZ), and CCSD(T)/CBS(Δa(DT)Z), respectively, and differ only in the construction of the ΔCCSD(T) correction terms. All energies were corrected for basis set superposition 1088

dx.doi.org/10.1021/jp209536e |J. Phys. Chem. A 2012, 116, 1086–1092

The Journal of Physical Chemistry A

ARTICLE

Table 1. Estimated CCSD(T)/CBS and ccCA Interaction Energies in kcal mol1 for Noncovalent Complexes Containing Divalent Sulfur complex (symmetry)

CBS

CBS

CBS

(ΔaDZ)a (ΔaTZ)b (Δa(DT)Z)c ccCAd

H2O dimer (Cs)e

4.89

4.96

4.98

4.85

H2OH2S (Cs)

2.85

2.90

2.92

2.81

H2SH2O (Cs)

2.60

2.66

2.69

2.59

H2S dimer (Cs) CH3SHH2O (Cs)

1.64 2.31

1.67 2.36

1.68 2.38

1.89 2.28

CH3SHNH3 (Cs)

2.89

2.95

2.98

2.88

CH3SHH2CO (C1)

3.07

3.11

3.13

3.03

CH3SH dimer (C1)

3.12

3.15

3.16

3.08

formamideCH3SH (C1)

6.24

6.33

6.37

6.21

H2OCH3SCH3 (Cs)

5.28

5.35

5.39

5.22

CH3OHCH3SCH3 (Cs)

5.62

5.69

5.72

5.57

H2SCH4 (C1) CH3SCH3benzene (C2v)

0.84 1.17

0.85 1.12

0.85 1.10

0.82 1.17

CH3SHbenzene (Cs)

3.53

3.50

3.48

3.48

ΔCCSD(T) is the difference between the MP2/aDZ and CCSD(T)/ aDZ correlation energies. b ΔCCSD(T) is the difference between the MP2/aTZ and CCSD(T)/aTZ correlation energies. c ΔCCSD(T) is the difference between the extrapolated MP2/a(DT)Z and CCSD(T)/ a(DT)Z correlation energies. d Excluding zero-point energy corrections. e Shown for comparison to water dimer analogues containing H2S. a

error using the counterpoise method.43 The frozen core approximation was used for all MP2 and CCSD(T) calculations. All calculations were performed using NWChem.45  ezac et al.,9 eight-point Following the S66  8 procedure of R CCSD(T)/CBS(ΔaTZ) dissociation curves were generated in which the noncovalent interaction coordinate in the MP2/TZCP optimized complexes was scaled by a factor of 0.9, 0.95, 1.05, 1.1, 1.25, 1.5, and 2.0. For the hydrogen-bonded complexes, the interaction coordinate is defined as the vector between the hydrogen atom of the donor hydrogen and the acceptor heavy atom. For the dispersion and mixed complexes, the definition of the interaction coordinate is system-dependent. For H2SCH4, the vector defined by the carbon atom of CH4 and the nearest hydrogen atom of H2S was used. For CH3SHbenzene, the interaction coordinate is defined as the vector between the sulfhydryl hydrogen and the center of the aromatic ring. For CH3SCH3benzene, the sulfur atom and the aromatic ring center define the interaction vector. Geometries were not reoptimized for the displaced complexes. The ccCA method32 was also used to compute interaction energies for each complex at the MP2/TZ-CP minima. The ccCA interaction energies were computed as follows. The HF/ CBS limit was obtained using a two-point a(TQ)Z extrapolation. To this HF reference energy was added a ΔMP2/CBS correction, which was computed as a three-point extrapolation of the MP2/ aDZ, aTZ, and aQZ correlation energies. A ΔCCSD(T) correction was computed as the difference between CCSD(T)/TZ and MP2/TZ correlation energies. Additionally, ccCA includes corrections to account for corevalence interactions, scalar relativity, spinorbit effects (where applicable), and zero-point vibrational energies, but zero-point corrections were not included here. All of the extrapolations and correction terms in ccCA are described in detail in the literature.3234

Figure 3. Convergence of the ΔCCSD(T) term for (top) H2OH2S, (middle) H2SH2O, and (bottom) H2S dimer using aXD (X = D,T,Q) basis sets and two-point a(DT)Z and a(TQ)Z extrapolations.

’ RESULTS AND DISCUSSION Geometries. MP2/TZ-CP gradient optimizations were performed for each complex. The optimized geometries for the hydrogen-bonded complexes are shown in Figure 1, and the dispersion and S/π complexes are shown in Figure 2. These geometries were used to generate dissociation curves as a function of the noncovalent interaction coordinate, with CCSD(T)/CBS(ΔaTZ) energies being computed for each point. Interaction Energies. Three different estimates of the CCSD(T)/CBS interaction energies for each complex, which differ only in the way the ΔCCSD(T) term is computed, are shown in Table 1. For consistency with the S66 set, we computed CBS(ΔaDZ) interaction energies. We also computed CBS(ΔaTZ) and CBS(Δa(DT)Z) energies, the latter of which are consistent with the approach of Takatani et al.7 and also with the revised S66 set.10 CBS(ΔaTZ) interaction energies for the dissociation curves are provided in the Supporting Information. The differences between the CBS(ΔaDZ) and CBS(ΔaTZ) interaction energies range from 0.01 to 0.09 kcal mol1 (calculated as CBS(ΔaTZ)  CBS(ΔaDZ)), with the largest 1089

dx.doi.org/10.1021/jp209536e |J. Phys. Chem. A 2012, 116, 1086–1092

The Journal of Physical Chemistry A

ARTICLE

Figure 4. CCSD(T)/CBS(ΔaTZ) eight-point dissociation curves for (a) water dimer and water dimer analogues containing H2S, (b) cysteine analogues, (c) methionine analogues, and (d) H2SCH4 and Sπ complexes. Nonequilibrium geometries were obtained by scaling the noncovalent interaction vector in the MP2/TZ-CP optimized complexes by factors of 0.9, 0.95, 1.05, 1.1, 1.25, 1.5, and 2.0. Energies and geometries are provided as Supporting Information.

deviation coming from formamideCH3SH. Therefore, CCSD(T)/CBS(ΔaDZ) appears to be a reasonable approximation for estimating CCSD(T) interaction energies near the CBS limit for sulfur-containing complexes, while avoiding costly CCSD(T)/ aTZ computations. Convergence of the ΔCCSD(T) Term. On the basis of an analysis of the convergence behavior of the ΔCCSD(T) term in several small noncovalent complexes,8 a similar analysis was performed here in which ΔCCSD(T) was computed using basis sets up to aQZ and two-point extrapolations up to a(TQ)Z for the three smallest sulfur-containing complexes considered in this work, H2OH2S, H2SH2O, and H2S dimer. Consistent with previous results, we found that the a(DT)Z extrapolations overcorrect relative to the a(TQ)Z values (Figure 3). It was also noted previously that using double-ζ basis sets lowers the accuracy significantly in extrapolations of the form used here and should not be included.46 We therefore consider CCSD(T)/ CBS(ΔaTZ) energies to be the most accurate that can be computed feasibly and dependably, although the Δa(DT)Z extrapolation also performs well for the cases considered here. All energies discussed in the remainder of this work are CCSD(T)/ CBS(ΔaTZ) values unless otherwise noted. Water Dimer Analogues. For comparison with the three H2S-containing water dimer analogs, we also computed CCSD(T)/CBS energies for the Cs-symmetric H2O dimer at its MP2/ TZ-CP minimum. As expected, none of the three complexes containing H2S interacts as strongly as H2O dimer, which has an interaction energy of 4.96 kcal mol1 (Table 1). Nevertheless, the computed interaction energies of the sulfur-containing systems are significant and non-negligible. For H2OH2S, the interaction energy is 2.90 kcal mol1. For H2SH2O, an

interaction energy of 2.66 kcal mol1 was found. The H2S dimer interaction energy is only 1.67 kcal mol1, confirming that SHS hydrogen bonds are significantly weaker than OHS and SHO hydrogen bonds, and showing the importance of considering each type of hydrogen bond in benchmark sets. Dissociation curves for H2O dimer and its analogues are shown in Figure 4. Clearly, equilibrium hydrogen-bond distances in the H2S-containing complexes are longer than in H2O dimer, and the curves near their respective minima are broader. At their respective minima, the OHS distance in H2OH2S (2.62 Å) is ∼0.3 Å longer than the SHO distance in H2SH2O (2.25 Å), even though its interaction energy is slightly stronger (Figure 4). The dissociation curve for H2S dimer is quite shallow and broad, with its minimum at a donoracceptor distance of 2.88 Å. Cysteine Analogues. Next, we consider hydrogen-bonded systems representative of the side chain of the amino acid cysteine. Systems for which CH3SH is the hydrogen-bond donor have interaction energies similar in magnitude to the H2S complexes. CH3 SHH 2 O is the weakest in this series, with an interaction energy of 2.36 kcal mol 1 (Table 1). CH3SHNH3, CH3SHH2CO, and CH3SH dimer each have interaction energies of ∼3 kcal mol1. In the formamide CH3SH complex, CH3SH is the hydrogen-bond acceptor, and an NH bond from the amide interacts strongly with the thiol. The interaction energy for CH3SHformamide is 6.33 kcal mol1, a consequence of the strong polar interactions between the amide and the sulfhydryl group. Methionine Analogues. Interaction energies computed for systems representing the side chain of methionine were found to be quite strong. In both complexes considered here, the thioether 1090

dx.doi.org/10.1021/jp209536e |J. Phys. Chem. A 2012, 116, 1086–1092

The Journal of Physical Chemistry A

ARTICLE

Table 2. Individual Corrections to the ccCA Interaction Energies in kcal mol1 for Noncovalent Complexes Containing Divalent Sulfur ΔCCa

ΔCVb

ΔDKc

H2O dimer (Cs)d

0.08

0.03

0.01

H2OH2S (Cs)

0.16

0.02

0.00

H2SH2O (Cs)

0.22

0.03

0.01

complex (symmetry)

H2S dimer (Cs)

0.25

0.28

0.01

CH3SHH2O (Cs) CH3SHNH3 (Cs)

0.19 0.28

0.02 0.03

0.01 0.02 0.00

CH3SHH2CO (C1)

0.23

0.03

CH3SH dimer (C1)

0.48

0.06

0.00

formamideCH3SH (C1)

0.25

0.05

0.01

H2OCH3SCH3 (Cs)

0.32

0.05

0.02

CH3OHCH3SCH3 (Cs)

0.35

0.06

0.02

H2SCH4 (C1)

0.09

0.02

0.00

CH3SCH3benzene (C2v) CH3SHbenzene (Cs)

0.89 1.16

0.05 0.05

0.01 0.02

ΔCC = CCSD(T)/cc-pVTZ  MP2/cc-pVTZ. b ΔCV = MP2(AE)/ccpCVTZ  MP2(FC)/cc-pVTZ. c ΔDK = MP2(DK)/cc-pVTZ-DK  MP2/cc-pVTZ. d Shown for comparison to water dimer analogues containing H2S. a

group serves as the hydrogen-bond acceptor, and the hydroxyl group from either H2O or CH3OH is the donor. The interaction energy of H2OCH3SCH3 was found to be 5.35 kcal mol1, and the interaction energy for CH3SCH3CH3OH is 5.69 kcal mol1 (Table 1). The dissociation curves for the two methionine analogues are quite similar to each other, particularly at displacements of ∼3 Å or greater (Figure 4). Dispersion and Sπ Complexes. Interaction energies for H2SCH4, CH3SHbenzene, and CH3SCH3benzene are shown in Table 1. For H2SCH4, the interaction energy is 0.85 kcal mol1, indicative of weak dispersive character. For CH3SCH3benzene, the interaction energy is 1.12 kcal mol1, again consistent with the expected magnitude for dispersiondominated complexes. However, the interaction energy for CH3SHbenzene is 3.50 kcal mol1, indicating mixed hydrogen bonding and dispersion characteristics. The interaction between CH3SH and benzene is slightly stronger than for H2Sbenzene in the “hydrogens down” configuration,17 likely a consequence of additional dispersive interactions between the methyl group of the thiol and the aromatic ring. Dissociation curves for the three dispersion/mixed complexes are shown in Figure 4. Comparison of CCSD(T)/CBS(ΔaTZ) and ccCA Interaction Energies. As a comparison of two distinct ab initio approaches for computing interaction energies, ccCA energies were computed for each complex at the MP2/TZ-CP geometries. The ccCA energies are shown in Table 1. Encouragingly, the ccCA and CCSD(T)/CBS(Δa(DT)Z) energies agree to within ∼0.1 kcal mol1, with one exception. Using ccCA, H2S dimer is 0.22 kcal mol1 lower in interaction energy than the CCSD(T)/ CBS(ΔaTZ) value. This difference is attributed to a larger ΔCV correction as compared to the other complexes (Table 2). H2S dimer and CH3SH dimer are the only complexes in the set that contain two sulfurs, but CH3SH dimer does not exhibit the same behavior in the ΔCV term. Presumably, this effect derives from the dominant sulfursulfur interactions in H2S dimer, which are less pronounced in the larger CH3SH dimer. The ΔDK

relativistic corrections are all very small for these systems and may be safely neglected.

’ CONCLUSIONS The CCSD(T)/CBS(ΔaTZ) approach has been used to compute benchmark-quality noncovalent interaction energies for several hydrogen-bonded and dispersion-bound complexes of biological relevance that contain divalent sulfur. Eight-point dissociation curves were computed for each complex along the corresponding noncovalent interaction vector. Interaction energies were also computed using the ccCA approach and were found to be within 0.1 kcal mol1 of the CCSD(T)/CBS values, on average. Both approaches provide accurate energies and should be useful for developing and assessing the accuracy of more approximate ab initio, DFT, semiempirical, and classical force field methods, particularly for noncovalent complexes containing divalent sulfur. ’ ASSOCIATED CONTENT

bS

Supporting Information. Optimized Cartesian coordinates and energies for all complexes, and additional MP2 and ΔCCSD(T) energies. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT This research was supported by the Office of Biological and Environmental Research, U.S. Department of Energy (DOE) through the Mercury Science Focus Area (SFA) Program at Oak Ridge National Laboratory (ORNL). ORNL is managed by UTBattelle, LLC, for the U.S. Department of Energy under contract DE-AC05-00OR22725. Computer time was provided by the National Center for Computational Sciences at Oak Ridge National Laboratory and the National Science Foundation through TeraGrid/XSEDE resources provided by NCSA (Grants TG-MCA08X032, UT-TENN0004, and TG-CHE090035). Most of the calculations were performed on the Kraken supercomputer at the National Institute for Computational Sciences (http://www. nics.tennessee.edu/). ’ REFERENCES

 erny, J.; Hobza, P. Phys. Chem. Chem. (1) Jurecka, P.; Sponer, J.; C Phys. 2006, 8, 1985–1993. (2) Raghavachari, K.; Trucks, G.; Pople, J.; Head-Gordon, M. Chem. Phys. Lett. 1989, 157, 479–483. (3) Riley, K.; Hobza, P. J. Phys. Chem. A 2007, 111, 8257–8263.  ezac, J.; Hobza, P. J. Chem. Theory (4) Grafova, L.; Pito nak, M.; R Comput. 2010, 6, 2365–2376. (5) Molnar, L.; He, X.; Wang, B.; Merz, K., Jr. J. Chem. Phys. 2009, 131, 065102. (6) Podeszwa, R.; Patkowski, K.; Szalewicz, K. Phys. Chem. Chem. Phys. 2010, 12, 5974–5979. (7) Takatani, T.; Hohenstein, E. G.; Malagoli, M.; Marshall, M. S.; Sherrill, C. D. J. Chem. Phys. 2010, 132, 144104. (8) Marshall, M. S.; Burns, L. A.; Sherrill, C. D. J. Chem. Phys. 2011, 135, 194102.  ezac, J.; Riley, K. E.; Hobza, P. J. Chem. Theory Comput. 2011, (9) R 7, 2427–2438.

1091

dx.doi.org/10.1021/jp209536e |J. Phys. Chem. A 2012, 116, 1086–1092

The Journal of Physical Chemistry A

ARTICLE

 ezac, J.; Riley, K. E.; Hobza, P. J. Chem. Theory Comput. 2011, (10) R 7, 3466–3470. (11) Markham, G.; Bock, C. J. Phys. Chem. 1995, 99, 10118–10129. (12) Rablen, P.; Lockman, J.; Jorgensen, W. J. Phys. Chem. A 1998, 102, 3782–3797. (13) Duan, G.; Smith, V.; Weaver, D. Mol. Phys. 2001, 99, 1689–1699. (14) Vila, A.; Mosquera, R. Chem. Phys. 2003, 291, 73–80. (15) Wennmohs, F.; Staemmler, V.; Schindler, M. J. Chem. Phys. 2003, 119, 3208–3218. (16) Tauer, T.; Derrick, M.; Sherrill, C. D. J. Phys. Chem. A 2005, 109, 191–196. (17) Ringer, A.; Senenko, A.; Sherrill, C. D. Protein Sci. 2007, 16, 2216–2223. (18) Raub, S.; Marian, C. J. Comput. Chem. 2007, 28, 1503–1515. (19) Howard, D.; Kjaergaard, H. Phys. Chem. Chem. Phys. 2008, 10, 4113–4118. (20) Biswal, H.; Shirhatti, P.; Wategaonkar, S. J. Phys. Chem. A 2009, 113, 5633–5643. (21) Sherrill, C. D.; Sumpter, B.; Sinnokrot, M.; Marshall, M.; Hohenstein, E.; Walker, R.; Gould, I. J. Comput. Chem. 2009, 30, 2187–2193. (22) Biswal, H.; Wategaonkar, S. J. Phys. Chem. A 2009, 113, 12774–12782. (23) Solimannejad, M.; Gharabaghi, M.; Scheiner, S. J. Chem. Phys. 2011, 134, 024312. (24) Cabaleiro-Lago, E. M.; Rodríguez-Otero, J.; Pe~na-Gallego, A. J. Chem. Phys. 2011, 135, 134310. (25) Reid, K.; Lindley, P.; Thornton, J. FEBS Lett. 1985, 190, 209–213. (26) Allen, F.; Bird, C.; Rowland, R.; Raithby, P. Acta Crystallogr., Sect. B 1997, 53, 696–701. (27) Zauhar, R.; Colbert, C.; Morgan, R.; Welsh, W. Biopolymers 2000, 53, 233–248. (28) Meyer, E.; Castellano, R.; Diederich, F. Angew. Chem., Int. Ed. 2003, 42, 1210–1250. (29) Zhou, P.; Tian, F.; Lv, F.; Shang, Z. Proteins: Struct., Funct., Bioinf. 2009, 76, 151–163. (30) Sherrill, C. D.; Takatani, T.; Hohenstein, E. J. Phys. Chem. A 2009, 113, 10146–10159. (31) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Chem. Rev. 1994, 94, 1887–1930. (32) DeYonker, N.; Cundari, T.; Wilson, A. K. J. Phys. Chem. 2006, 124, 114104. (33) DeYonker, N.; Wilson, B.; Pierpont, A.; Cundari, T.; Wilson, A. K. Mol. Phys. 2009, 107, 1107. (34) Mintz, B. J.; Williams, T.; Howard, L.; Wilson, A. K. J. Chem. Phys. 2009, 130, 234104. (35) DeYonker, N.; Peterson, K.; Steyl, G.; Wilson, A. K.; Cundari, T. J. Phys. Chem. A 2007, 111, 11269–11277. (36) Williams, T.; Wilson, A. K. J. Sulfur Chem. 2008, 29, 353. (37) DeYonker, N.; Mintz, B. J.; Cundari, T.; Wilson, A. K. J. Chem. Theory Comput. 2008, 4, 328–334. (38) Dunning, T. J. Chem. Phys. 1989, 90, 1007. (39) Kendall, R.; Dunning, T., Jr.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6796. (40) Dunning, T., Jr.; Peterson, K.; Wilson, A. K. J. Chem. Phys. 2001, 114, 9244. (41) Wilson, A. K.; Dunning, T., Jr. J. Phys. Chem. A 2004, 108, 3129–3133. (42) Bell, R.; Wilson, A. K. Chem. Phys. Lett. 2004, 394, 105–109. (43) Boys, S.; Bernardi, F. Mol. Phys. 1970, 19, 553–566. (44) Halkier, A.; Klopper, W.; Helgaker, T.; Jorgensen, P.; Taylor, P. R. J. Chem. Phys. 1999, 111, 9157–9167. (45) Valiev, M.; Bylaska, E. J.; Govind, N.; Kowalski, K.; Straatsma, T. P.; van Dam, H. J. J.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T. L.; de Jong, W. A. Comput. Phys. Commun. 2010, 181, 1477–1489. (46) Halkier, A.; Helgaker, T.; Jorgensen, P.; Klopper, W.; Kock, H.; Olsen, J.; Wilson, A. K. Chem. Phys. Lett. 1998, 286, 243–252. (47) De Marothy, S. XYZViewer Version 0.97; http://www.physto. se/∼sven/. 1092

dx.doi.org/10.1021/jp209536e |J. Phys. Chem. A 2012, 116, 1086–1092