Protein–Ligand Interaction Energies with Dispersion Corrected

ARTICLE pubs.acs.org/JPCA

Protein
Ligand Interaction Energies with Dispersion Corrected Density Functional Theory and High-Level Wave Function Based Methods Jens Antony,† Stefan Grimme,*,† Dimitrios G. Liakos,‡ and Frank Neese‡ † ‡

Organisch-Chemisches Institut, Universit€at M€unster, Corrensstrasse 40, D-48149 M€unster, Germany Lehrstuhl f€ur Theoretische Chemie, Universit€at Bonn, Wegelerstr. 12, D-53115 Bonn, Germany ABSTRACT: With dispersion-corrected density functional theory (DFT-D3) intermolecular interaction energies for a diverse set of noncovalently bound protein
ligand complexes from the Protein Data Bank are calculated. The focus is on major contacts occurring between the drug molecule and the binding site. Generalized gradient approximation (GGA), metaGGA, and hybrid functionals are used. DFT-D3 interaction energies are benchmarked against the best available wave function based results that are provided by the estimated complete basis set (CBS) limit of the local pair natural orbital coupled-electron pair approximation (LPNO-CEPA/1) and compared to MP2 and semiempirical data. The size of the complexes and their interaction energies (ΔEPL) varies between 50 and 300 atoms and from
1 to
65 kcal/mol, respectively. Basis set effects are considered by applying extended sets of triple- to quadruple-ζ quality. Computed total ΔEPL values show a good correlation with the dispersion contribution despite the fact that the protein
ligand complexes contain many hydrogen bonds. It is concluded that an adequate, for example, asymptotically correct, treatment of dispersion interactions is necessary for the realistic modeling of protein
ligand binding. Inclusion of the dispersion correction drastically reduces the dependence of the computed interaction energies on the density functional compared to uncorrected DFT results. DFT-D3 methods provide results that are consistent with LPNO-CEPA/1 and MP2, the differences of about 1
2 kcal/mol on average ( X), the two cardinal numbers of the two is the correlation energy calculated basis sets used, E(LPNO;X) corr with LPNO-CEPA/1 within the smaller basis set and EYSCF the ,E(MP2;Y) SCF energy calculated with the larger basis set. E(MP2;X) corr corr are the MP2 correlation energies calculated with the two basis the two-point extrapolated MP2 correlation sets and E(MP2;∞) corr energy estimated as follows: ¼ EðMP2;∞Þ corr

Þ X 3 EðMP2;XÞ
Y 3 EðMP2;Y corr corr X3
Y 3

The basis set superposition error (BSSE) has been treated by counterpoise (CP) correction79 when indicated. Density fitting techniques, also called resolution-of-the-identity (RI) approximation, were used throughout to speed up MP280 and DFT GGA9,10 computations roughly by a factor of 10
20 at insignificant extra inaccuracy. The default auxiliary basis sets for Coulomb fitting and correlation calculations were taken. All results have been obtained with a slightly modified version of the Turbomole81 and the ORCA82 program packages. For singlepoint energy calculations and the geometry optimization of the hydrogen atom positions, tight SCF (scfconv = 7 in Turbomole, TightSCF in ORCA) and optimization (TightOpt in ORCA) convergence critera are used. For the quadrature of the exchangecorrelation terms, the multiple grid (m4 in Turbomole, Grid4 in ORCA) has been specified. Semiempirical calculations were done with MOPAC83 and the DFTB84 program, M0685 results were obtained with Gaussian.86

III. RESULTS AND DISCUSSION A. Protein
Ligand Model Complexes. In case of protein
ligand interactions, the natural source of three-dimensional structure information are experimentally (by X-ray crystallography or NMR spectroscopy) determined data from the Protein Data Bank (PDB). Because more than 10000 structures of protein
ligand complexes are deposited, the question of a representative but still realistically sized subset for which quantum calculations are feasible has to be addressed. For the purpose of scoring function development and evaluation, databases as AffinDB,87 BindingDB,88 Binding MOAD,89 the protein-small-molecule database,90 or 11212

dx.doi.org/10.1021/jp203963f |J. Phys. Chem. A 2011, 115, 11210–11220

The Journal of Physical Chemistry A

ARTICLE

Figure 1. Three examples of model structures used for protein
ligand interaction energy calculation. Top left: mitochondrial matrix adenylate kinase and its substrate AMP (2AK3). Top right: the peptide CF3-Leu-Ala-NHC6H4-CF3 with porcine pancreatic elastase (7EST). Bottom: C2 symmetrybased diol inhibitor of HIV-1 protease (1HVI). The full line encloses the drug molecule.

PDBbind91 are available, supplying affinity data for structurally resolved protein
ligand complexes from the Protein Data Bank (PDB). We started from a compilation of Raha and Merz, which has already been used for linear scaling quantum chemical calculations at the semiempirical level.16 Our focus is on modeling major contacts occurring between the drug molecule and the binding site. Heteroatoms not in the active site and water molecules are excluded. All amino acids within 3.2 Å of the ligand and the ligand were extracted from the coordinate file. Depending on the size of the drug molecule and its binding geometry, the cutoff is reduced to 2.7 Å to avoid too large complexes or increased to 3.6 Å when no part or not enough of the protein was included for smaller values. Hydrogen atoms were added on the basis of standard geometries. Ionizable residues were protonated assuming neutral side chains. Histidines are protonated at Nε. N and C termini of the broken backbone are kept neutral. The positions of hydrogens were optimized with the ORCA program82 using the BP86 functional68,71,92 with the old (-D2) dispersion correction67 and a valence double-ζ Gaussian basis set (SVP).76 Heavy atoms were fixed at experimental positions with the ORCA option optimizehydrogens. The resulting rms deviations of the heavy atoms after optimization from the experimental positions were below 0.01 Å except for four cases were the rms deviation increased up to 0.08 Å. Figure 1 shows three examples of the resulting complexes.

B. Protein
Ligand Interaction Energies by MP2 and LPNOCEPA/1. For benchmark purposes, MP2/CBS and LPNO-CEPA/

1(TZVP) calculations on 15 of the first 25 protein
ligand model complexes were performed. The results are summarized in Table 1. The listed total interaction energies cover almost 2 orders of magnitude, from 0.2 to
63.3 kcal/mol. The correlation energy contribution is crucial for correctly modeling protein
ligand binding, as all complexes are seriously underbound or even unbound at the Hartree
Fock level. All levels reported in Table 1 provide the same sequence of protein
ligand interaction energies. The absolute values, however, differ according to the basis set used and the treatment of electron correlation. At the MP2/TZVP level, the correlation energy contribution is overestimated by 1.0
7.0 kcal/ mol without and by 0.9
6.0 kcal/mol with CP correction compared to the LPNO-CEPA/1(TZVP) correlation energy contribution to protein
ligand binding. With the TZVP basis set, the correlation energy contains a sizable BSSE, ranging from 2.7 to 15.4 kcal/mol, with an average of 7.7 kcal/mol at the MP2 and from 2.6 to 14.4 kcal/ mol (on average 7.3 kcal/mol) at the LPNO-CEPA/1 level. The correlation energy correction, however, is affected by the BSSE only by a few tenth of a kcal/mol toward smaller values, and the Counterpoise correction to the LPNO-CEPA/1
MP2 correlation energy difference is on average equal to
0.4 kcal/mol and does not exceed
1 kcal/mol in the protein
ligand model complexes considered. The CP correction to the MP2/CBS values varies from 1.6 to 6.3 kcal/mol (on average 3.4 kcal/mol). Its effect on 11213

dx.doi.org/10.1021/jp203963f |J. Phys. Chem. A 2011, 115, 11210–11220

The Journal of Physical Chemistry A

ARTICLE

Table 1. BSSE (Un)corrected MP2/CBS and LPNO-CEPA/1(TZVP) Results (kcal/mol) ΔEPL MP2b CBS cc-pV(D/T)Z a

no CPC

CPC

ΔEPL,corr MP2c (TZVP) no CPC

CPC

differencee

LPNO-CEPA/1d (TZVP) no CPC

CPC

LPNO-CEPA/1
MP2 no CPC

CPC

total ΔEPL

“best”

“CPC/2”f

DFTg

1


3.2


0.1


15.4


9.9


13.5


8.1

1.9

1.8

0.2

1.8

2


7.2


5.4


14.9


10.4


13.0


8.8

1.9

1.6


4.5


4.6

3 4


10.2
14.5


7.9
12.9


13.2
7.7


9.2
4.2


11.6
6.5


7.7
3.2

1.6 1.2

1.4 1.0


7.5
12.6


7.5
12.0

5


17.4


14.3


18.0


11.1


15.6


9.0

2.4

2.1


13.5


12.7

6


17.6


16.0


5.0


2.3


3.9


1.4

1.0

0.9


15.8


15.7

7


18.1


16.4


5.1


2.4


3.9


1.3

1.2

1.1


16.1


16.3

8


18.7


16.9


5.7


2.6


4.3


1.4

1.3

1.2


16.5


16.7

9


21.8


18.5


19.3


11.3


16.4


8.7

3.0

2.6


17.4


17.6

10


26.6


21.5


36.4


23.2


31.4


18.7

5.0

4.5


19.3


18.4

11 12


27.6
31.8


23.6
27.1


17.6
37.9


11.2
24.9


14.6
33.2


8.5
20.8

3.0 4.7

2.6 4.1


22.8
25.0


22.7
24.6

14


42.7


36.4


39.3


25.0


32.6


19.2

6.7

5.8


33.3


32.1

15


43.8


38.4


31.6


19.5


26.8


15.1

4.9

4.4


36.5


36.3


28.3

13

16


37.2

17


40.2

18


43.3

19 20


47.1
49.2

21


48.8
59.7

22 23


72.8


66.9


38.2


22.8


31.1


16.8

7.0

6.1


63.3


64.7
69.1

24 a

Ordered according to BLYP-D3/TZVP intermolecular interaction energies. b CBS extrapolated energies with the cc-pVDZ/cc-pVTZ basis sets. c MP2(TZVP) correlation energy. d LPNO-CEPA/1(TZVP) correlation energy. e The difference ΔEPL,corr(LPNO-CEPA/1(TZVP))
ΔEPL, corr(MP2(TZVP)) calculated both with the BSSE-corrected and the BSSE-uncorrected energies. f Half of the sum of columns 1, 2, 7, and 8 (see text, section IIIB). g B97-D3/def2-TZVP (see text, section IIIC).

the MP2/CBS total interaction energy estimate is toward less negative values and thus opposite to that of the CP correction to the LPNO-CEPA/1
MP2 correlation energy difference, overcompensating for the latter and resulting in an average CPC of 3.1 kcal/mol in the total energy estimate. In the total interaction energies, only half the CPC both for MP2/CBS and LPNO-CEPA/1
MP2 correlation energy correction is taken, because the counterpoise correction is known to overestimate the BSSE.93
96 We illustrate this here for the S22 benchmark set of intermolecular interaction energies:97 the CP uncorrected MP2/CBS from the cc-pVDZ and cc-pVTZ basis sets overestimates the MP2/CBS obtained by extrapolation with the aug-cc-pVTZ and aug-cc-pVQZ basis sets (CP corrected) with a mean absolute deviation (MAD) of 0.75 kcal/mol, while the MP2/CBS from the cc-pVDZ and cc-pVTZ basis sets with CPC underestimates the same MP2/CBS reference with a MAD of 0.34 kcal/mol. With only half of the CPC, the MP2/CBS from the cc-pVDZ and cc-pVTZ basis sets again overestimates the MP2/CBS(aug-cc-pV(T/Q)Z) with CPC, but now the MAD is only 0.21 kcal/mol. The same argument is valid if the MP2/CBS is obtained by extrapolation with the aug-cc-pVTZ and aug-ccpVQZ basis sets without CPC (the MADs amount to 0.69 kcal/ mol for MP2/CBS(cc-pV(D/T)Z) without, 0.39 kcal/mol for MP2/CBS(cc-pV(D/T)Z) with, and 0.19 kcal/mol for

MP2/CBS(cc-pV(D/T)Z) with half CPC) or for MP2/CBS(cc-pV(Q/5)Z). For the S22 set, we have applied the same scheme for obtaining LPNO-CEPA/1(CBS) values and found a small underbinding tendency on average (the MAD(MD) with respect to CCSD(T)/CBS98 is 0.5 (
0.5) kcal/mol). This behavior is more pronounced for the hydrogen-bonded systems in S22 (the typical deviation is 6
8%) than for the dispersion-dominated or mixed systems for which LPNO-CEPA/1(CBS) is very accurate (typical deviation 10 kcal/mol) systems (about 2
3 kcal/mol otherwise). The reference values from column 9 in Table 1 are abbreviated as LPNO-CEPA/1(CBS) in the following. It is of course important to have an estimate of how accurate the reference values are expected to be. Recently, extensive benchmark calculations of weak interaction energies on the basis of the LPNOCEPA/1 method have been performed. The results indicate that relative to the estimated one- and many-particle limits, the LPNOCEPA/1 method is accurate to within 0.3
0.4 kcal/mol for interaction energies (note that this estimate is based on a comparison 11214

dx.doi.org/10.1021/jp203963f |J. Phys. Chem. A 2011, 115, 11210–11220

The Journal of Physical Chemistry A

ARTICLE

Figure 2. BLYP-D3/TZVP (PP = TZVPP, QZ = QZVP) intermolecular interaction energy for 24 (21, 14) protein
ligand model complexes.

of electronic energies only).44,99 The basis sets used here are probably too small to reach a comparably high level of accuracy. However, even if the uncertaintities in the reference values would be larger by a factor of 2 compared to the benchmark results, the conclusion would still remain that the extrapolated LPNO-CEPA/1 values are certainly accurate enough for the purposes of this work. C. Protein
Ligand Interaction Energies by Dispersion Corrected Density Functional Theory and MP2. Figure 2 shows intermolecular interaction energies in the supermolecule approach without counterpoise correction for basis set superposition error (BSSE) for 24108 of the first 25 of the 165 complexes of ref 16, obtained with the BLYP functional,71,72 the D3-dispersion correction,7 and valence triple-ζ (TZVP and TZVPP)77 and quadruple-ζ (QZVP)78 Gaussian basis sets, ordered according to the interaction strength. Remaining BSSE are on the order of only 5% of ΔEPL with triple-ζ (TZ) basis sets for large complexes and are at least partially absorbed into the DFT-D3 parameters. The two empirical parameters for each DF (sr,6 and s8) were determined by a least-squares fit to (relative) reference energies, most of them based on estimated CCSD(T)/ CBS results.7 The DFT data are derived from calculations with the very large def2-QZVP AO basis sets. Thus, the resulting parameters are essentially free of basis set superposition error (BSSE) and other signifcant incompleteness effects. In the Supporting Information of ref 7, corresponding values optimized using a triple-ζ AO basis that would be more appropriate for an efficient treatment of very large molecular systems and that accounts implicitly for the then more relevant BSSE effects are given. However, in the present work we resort to the standard DFT-D3 parameters throughout. The interaction energies cover a range from about 0 kcal/mol (7EST model complex) to 80 kcal/mol (1HIH). The dispersion correction (or the error from plain DFT) increases with increasing absolute interaction energy. The contribution of the D3-correction dominates the model protein
ligand binding energies. Uncorrected DFT seriously underbinds and, in some instances, the complexes are even unbound. As is already well-known, the basis set convergence of DFT-GGA interaction energies is very fast, and we observe here (Figure 2) almost no difference between TZ and QZ results. This also demonstrates that the remaining BSSE are indeed negligible for the investigated systems.

Figure 3. PBE-D3, B97-D3, TPSS-D3, BLYP-D3, B3LYP-D3, M06, and MP2 intermolecular interaction energy for protein
ligand model complexes ((def2-)TZVP basis) compared to LPNO-CEPA/1(CBS). The plot in the bottom shows shifted energies to allow a more convenient comparison of the data.

In Figure 3, results obtained with different GGA-type density functionals (PBE,69,70 B97-D,67 and BLYP71,72), the meta-GGA TPSS,75 the hybrid density functional B3LYP,73,74 all with the D3dispersion correction7 and a valence triple-ζ Gaussian basis set (TZVP), and MP2(TZVP) with half CP correction are compared to our LPNO-CEPA/1(CBS) interaction energies. The agreement of the DFT-D3 protein
ligand interaction energies is excellent: with few exceptions do the dispersion-corrected functionals deliver the same energetic ordering between the complexes. The D3-correction strongly reduces the dependence on the density functional (compared to pure DFT). The shifted (relative) DFT-D3 interaction energies cover a range of up to 4 kcal/mol for the PBE, BLYP, TPSS, and B3LYP functionals. This is particularly clearly visible in the bottom part of Figure 3, where a linearly increasing value is added to the interaction energy (3  complex No.
3 kcal/mol), allowing a larger scale for the representation of the energy values. The largest outliers were obtained for the density functional with the weakest performance for noncovalent interactions (BP86,68 which has an MAD for the S22 benchmark set97 of 0.62 kcal/mol compared to 0.62 kcal/mol for PBE, 0.46 kcal/mol for TPSS and B3LYP, and 0.23 kcal/mol for BLYP7) and are not shown in Figure 3. With the BP86 results included, the DFT-D3 relative interaction energy range is increased from 4 to 7 kcal/mol. The uncorrected DFT relative interaction energies, in contrast, cover a range of up to more than 11215

dx.doi.org/10.1021/jp203963f |J. Phys. Chem. A 2011, 115, 11210–11220

The Journal of Physical Chemistry A

ARTICLE

Table 2. Statistics of the Deviation between DFT-D3, M06, as well as MP2 and Reference ΔEPL Values (kcal mol
1) with TZVP (TZVPP/QZVP) Basis for 15 Protein
Ligand Model Complexes (14 for QZVP)a MP2b,c

M06b

MD

1.1 (
3.6/5.9)

MADg rmsdh

1.6 (3.6/5.9) 1.5 (4.4/6.7)

MAX
MINi

3.5 (8.6/11.6)

6.3

f

B3-LYPd

TPSS

B97-De

0.5


3.6


2.6 (
1.8/
0.8)

1.7 2.1

3.6 4.1

2.6 (1.8/1.0) 3.1 (2.2/1.2)

B-LYP

PBE

BP86

0.3 (0.3)


3.2 (
2.3/
0.9)


3.5


4.7

0.5 (0.6) 0.7 (0.8)

3.2 (2.3/0.9) 3.7 (2.6/1.2)

3.5 3.9

4.7 5.8

7.7

5.9 (4.0/3.5)

2.9 (2.9)

6.9 (4.6/2.8)

7.3

10.5

MDf


2.9


1.7 (
0.8/0.1)


0.3 (
0.2)


2.1 (
1.2/0.0)


2.7


3.1

MADg

2.9

1.7 (1.0/0.9)

0.9 (0.8)

2.1 (1.3/0.6)

2.7

3.3

rmsdh

3.5

2.2 (1.5/1.1)

1.2 (1.2)

2.6 (1.6/0.8)

3.3

4.1

MAX
MINi

6.4

4.9 (3.8/3.7)

5.1 (5.3)

5.1 (3.5/3.0)

5.9

9.4

a

Top: Zero-damping7 (except for MP2 and M06); bottom: BJ-damping.101 b No dispersion correction. c TZVP basis with half (no/full) counterpoise correction. d 11 complexes. e def2-TZVP (def2-TZVPP) basis. f Mean deviation. g Mean absolute deviation. h Root mean square deviation. i Error spread (largest positive minus largest negative deviation).

20 kcal/mol with BLYP and PBE at the extremes. B97-D3/def2TZVP results also increase the range of DFT-D3 interaction energies, but in contrast to BP86 toward smaller values where the LPNO-CEPA/1(CBS) references are located. The MP2(TZVP) results are consistent with the DFT-D3 interaction energies, reassuring that the true interaction energy for the protein
ligand model complexes will not be totally off the range covered by the density functionals. Note that the performance of MP2 for noncovalent interactions of the systems in the S22 set with an MAD of 0.79 kcal/mol in the CBS limit and of 0.75 (1.10/0.70) kcal/mol for the TZVPP atomic orbital basis set with half (no/ full) counterpoise correction100 is even worse than that of the dispersion corrected density functionals above. The BP86, PBE, B-LYP, TPSS, and B3-LYP DFT-D3 methods overestimate the absolute LPNO-CEPA/1(CBS) intermolecular interaction energies, as can be seen also from the mean deviations (MDs) given in Table 2, which are all negative and of the same size as the mean absolute deviations (MADs). In contrast to this, the B97-D67 GGA, which is the best performing DFT-D3 method with an MAD of 0.5 kcal/mol, has a positive MD of 0.3 kcal/mol. The worst performing of the density functionals tested in this work is BP86, all with respect to MAD, rmsd, and error spread, which is equal to minus of the largest negative deviation, as no positive deviation occurs. CP uncorrected MP2 only beats BP86 and has a less regular behavior than all DFT-D3 results in relation to the LPNO-CEPA/1(CBS) reference (Table 2, first column, top), whereas the MADs for half the CP correction and of the hybrid meta-GGA M0685 are the same as that of TPSS-D3 for the benchmark set of this study (Table 2, second column, top). Very recently, replacement of the zero-damping function that merges the classical and density based descriptions of electron correlation by finite damping, as introduced by Becke and Johnson (BJ-damping), was proposed.101 For typical noncovalent interactions, it was shown that DFT-D3(BJ) is on average for 12 standard density functionals slightly better than “zero-damping” DFT-D3. For the 15 protein
ligand model complexes of this work where LPNO-CEPA/1(CBS) results are available, a reduction of the MAD by 0.7 kcal/mol (B3-LYP) up to 1.4 kcal/mol (BP86) is obtained upon the replacement of zero-damping by BJ-damping, with the exception of B97-D, whose MAD increases by 0.4 to 0.9 kcal/mol, which is still among the best performing DFT-D3 methods examined. Thus, in most cases, BJ-damping results in somewhat better interaction energies than zero-damping. So far, MADs of DFT-D3 results obtained with the TZVP basis have been discussed. Because the DFT-D3 parameters have

been obtained close to the basis set limit,7,101 the TZVP basis might be too small, leading to inconsistencies between the DFT and the dispersion correction description of the protein
ligand interaction complexes. As a first step to estimate basis set effects on the results, the properly polarized TZVPP set is used instead of underpolarized TZVP. We see from Table 2 that the MADs improve for both TPSS and B-LYP by 0.8 and 0.9 kcal/mol, respectively, for zero-damping and by 0.7 and 0.8 kcal/mol, respectively, for BJ-damping. This trend is continued when passing to QZVP. The MADs of B-LYP from LPNO-CEPA/ 1(CBS) intermolecular interaction energies for 14 protein
ligand model complexes are 0.9 and 0.6 kcal/mol for zeroand BJ-damping, respectively, and 1.0 and 0.9 kcal/mol for these variants in case of TPSS. Thus, with a large basis set, dispersioncorrected BLYP and TPSS perform equally well as B97-D/def2TZVP within the accuracy of the reference. The good performance of BLYP-D3/QZVP is in line with the usual benchmark studies on intermolecular interactions.7,101 It also demonstrates the reliability of the estimated LPNO-CEPA/1(CBS) reference of the current work. The basis set effect on B97-D from def2TZVP to def2-TZVPP is small due to the similarity of the two basis sets. If for TZ to QZ the basis set effect on B97-D would be of the same size as that of B-LYP, the MADs for def2-QZVP were expected to be 2.0 and 1.5 kcal/mol for zero- and BJ-damping, respectively. Thus, although the B97-D3/def2-TZVP is by far the best DFT level compared to the reference, there might be some error compensation between electronic and BSSE effects involved. However, for practical reasons, we will stick to the B97D3/def2-TZVP results for further benchmarking in the following, in particular, because the QZ results are currently available for not more structures than the LPNO-CEPA/1(CBS) interaction energies. This DFT-D3 level is also recommended for future applications to other protein
ligand interactions. D. Semiempirical Protein
Ligand Interaction Energies. In Figure 4, AM1 and PM3 semiempirical intermolecular interaction energies are compared to B97-D3/def2-TZVP results. Because pure AM1 and PM3 intermolecular interaction energies are throughout underbound in comparison to DFT-D3, we added the dispersion correction in the zero-damping variant7 with the two parameters (sr,6 and s8) preliminarily adjusted only to the interaction energies of the S22 set97 (for the description of the full set of reference data for determination of sr,6 and s8, see ref 7). Values of sr,6 = 0.498 and s8 = 0.541 amount to a minimal MAD of 1.49 kcal/mol for the dispersion-corrected AM1-D3 method, starting from an MAD of 6.81 kcal/mol for pure AM1. 11216

dx.doi.org/10.1021/jp203963f |J. Phys. Chem. A 2011, 115, 11210–11220

The Journal of Physical Chemistry A

ARTICLE

Table 3. Statistics of the Deviation between Semiempirical and Reference ΔEPL Values (kcal mol
1)a AM1-D3b PM3-D3b PM6 PM6-DH2 DFTBc DFTB-D3b,c MDd


5.6


7.6

8.1


4.9

2.3

0.3

MADe rmsdf

7.7 9.8

7.7 10.6

8.1 10.5

5.1 6.5

3.5 4.8

4.0 4.8

max
ming

33.1

22.1

28.1

16.7

14.3

17.1

MDd


9.3


9.6

11.5


7.1

2.6

0.0

MADe

10.8

10.2

11.5

7.4

3.9

3.8

rmsdf

13.7

12.8

14.0

9.6

5.0

4.7

max
ming

42.4

29.0

29.5

21.5

16.6

19.4

a

Top: estimated LPNO-CEPA/1(CBS) reference for 15 protein
ligand model complexes. Bottom: B97-D3/def2-TZVP “best” DFT reference for 24 protein
ligand model complexes. b Zero-damping.7 c 11/20 complexes. d Mean deviation. e Mean absolute deviation. f Root mean square deviation. g Error spread (largest positive minus largest negative deviation).

Figure 4. AM1 and PM3 intermolecular interaction energies without and with dispersion correction, PM6, PM6-DH2, and DF-TB(-D3 zerodamping) results in comparison to B97-D3/def2-TZVP.

Similarly, the D3-correction with sr,6 = 0.602 and s8 = 0.551 reduces the MAD of the pure PM3 interaction energies for the S22 set from 5.92 to 1.42 kcal/mol for PM3-D3. The remaining error of the dispersion corrected AM1 and PM3 methods is mainly due to a poor description of hydrogen bonds. Owing to the small sr,6 values for AM1 and PM3, the dispersion correction is very short ranged. Thus, the basic electronic structure part of both methods is over-repulsive. Returning to the present data on protein
ligand interaction energies, in 10 out of the 24 observed cases the semiempirical interaction energies match the DFT-D3 results well within their accuracy range (Figure 4). However, the remaining dispersioncorrected AM1 and PM3 interaction energies are with one exception (AM1-D3) in part seriously overbound (up to more than 20 kcal/mol). Furthermore, the energetic sequence given by the B97-D3/def2-TZVP results is not matched by AM1-D3 and PM3-D3 in all cases. All in all, the empirical dispersion correction enables a simple and significant improvement of the AM1 and PM3 methods. Pure semiempirical results themselves are unsatisfactory. However, without a subsequent adjustment of the Hamiltonian in the presence of the dispersion correction used, the approach cannot be recommended. We also include the PM6-DH2 method and the dispersion corrected density functional based tight binding (DFTB) approach into the comparison. The PM6-DH2 method was parametrized to reproduce interaction energies for geometries obtained from high-level quantum mechanical calculations to

correct binding errors of the PM6 method.102,103 As can be seen from Figure 4, dispersion contribution and hydrogen-bonding terms correct the PM6 intermolecular interaction energies, which are underbinding compared to B97-D3/def2-TZVP, toward and beyond the DFT-D3 results. Nine of the PM6-DH2 results are lying outside the region of the density functional interaction energies (4 kcal/mol) and are overbinding. The density functional based tight binding method104 has been extended for the description of long-range dispersive interactions.105 DFTB interaction energies are mostly too small in comparison to B97-D3/def2TZVP, but are all lying within 4 kcal/mol of the DFT-D3 reference and reproduce the relative ordering of the energies in general. Substituting the Slater
Kirkwood type model implemented in DFTB105 by the D3-dispersion correction adjusted to the S22 benchmark set of intermolecular interaction energies (parameters sr,6 = 1.699 and s8 = 1.504, total MAD for S22 1.14 kcal/mol with zero-damping, and a1 = 0.5705, s8 = 1.7851, and a2 = 4.5504, total MAD for S22 1.32 kcal/mol with BJ-damping) leads to similar accuracy. This result is consistent with the above findings for AM1 and PM3 that errors from the underlying semiempirical Hamiltonian are the accuracy limiting factor and not the details of a dispersion correction. Because the D3-correction is less empirical and more general than the standard treatment of dispersion in DFTB, it is recommended as a default. DFTB results for four complexes are missing in Figure 4 because of the nonavailability of parameters for P, Cl, and F in our version of the DFTB program. Dispersion-corrected AM1 and PM3 methods and PM6 have MADs of about 8 kcal/mol from the estimated LPNO-CEPA/ 1(CBS) reference (see Table 3). The PM6-DH2 approach improves on this with an MAD of 5 kcal/mol, while DFTB has a similar MAD as dispersion corrected PBE and B3-LYP or MP2(TZVP) without CP correction (see Table 2). Compared to the ”best” DFT reference for all 24 structures, MADs of dispersion corrected AM1 and PM3 as well as of PM6 further deteriorate to 11 kcal/mol and of PM6-DH2 to 7 kcal/mol, while the DFTB results in the two versions remain constant at about 4 kcal/mol for 20 structures. To address the impact of strong H-bond interactions on outliers, the number of H-bonds between the drug and receptor model complex is determined for the 24 structures by the application of geometric (distance and angle) criteria to donor
acceptor pairs with molden106 and related to the PM6-DH2 error, 11217

dx.doi.org/10.1021/jp203963f |J. Phys. Chem. A 2011, 115, 11210–11220

The Journal of Physical Chemistry A that is, the deviation of PM6-DH2 intermolecular interaction energies from the B97-D3/def2-TZVP “best” DFT reference. The number of H-bonds between drug and receptor model complex varies from 0 to 10. On average, the PM6-DH2 error increases with an increasing number of H-bonds from 0 to
17.5 kcal/mol (linear regression: PM6-DH2 error =
0.23691
1.7594  number of H-bonds). However, a considerable spread of data is observed, reflected in the correlation coefficient
0.7930889. In particular, for a fixed number of three H-bonds, the PM6-DH2 error varies between
14.5 and 1.3 kcal/mol, and the strongest outlier of
19.6 kcal/mol is obtained for 7 instead of 10 H-bonds. These findings are equally valid if the PM6-DH2 error is based on the LPNO-CEPA/1(CBS) reference instead of the “best” DFT reference though at a somewhat smaller data basis (linear regression: PM6-DH2 error =
1.369
1.3404  number of H-bonds, correlation coefficient
0.7118532). Therefore, no simple relation between the, thus, determined number of H-bonds and the PM6DH2 error can be established, and we expect the same for similarly basic measures of aromatic
aromatic and aromatic
nonpolar interactions.

IV. CONCLUSIONS Truncated protein
ligand model complexes have been investigated with dispersion corrected density functional theory (DFT-D3) and LPNO-CEPA/1(CBS) wave function based methods as a nonempirical reference. A wide range of representative noncovalent interactions (hydrogen bonds, stacking, and so on) in more or less randomly chosen complexes is considered resulting in sizable systems from about 50 to 300 atoms. It is found that an additional dispersion term in DFT is of utmost importance even for a qualitatively correct description of the interactions. The standard DFT-D3 dispersion model has been applied in two variants (damping functions) but without any special adjustment for the investigated systems. The correction brings the interaction energies in quantitative agreement with the reference values and reduces the spread due to different density functionals (Figure 3). Overall, the consistency in the quantum mechanical results is extremely high. Because grossly different electronic structure methods have been used, this provides a high degree of confidence in the quality of the reference values. It also, once more, underlines that models as sophisticated as LPNO-CEPA/1 can be used in a black box fashion. However, while the routine application of such methods to systems of the size investigated here should probably be considered as a significant advance in electronic structure methodology, the computational effort involved in these calculations is still too high to allow for their routine use in screening studies. Hence, we concentrate below on the quality of the DFTD3 results that are obtained with very high efficiency. The quality of the DFT-D3 interaction energies suggests that the range of interactions in the full protein, which is omitted in the model complexes will be represented correctly with this method, too. No bias in the selection of the structures has been introduced, thus, the set of protein
ligand model complexes can be considered as representative for those intermolecular interactions occurring in real protein
ligand binding reactions and hence our conclusions are expected to hold generally. Concerning the intermolecular interaction energies, a consensus among reliable, first-principles quantum chemical methods is obtained. For small interaction energies the apparent DFTD3 error is in the range of 1
2 kcal/mol. For large interaction

ARTICLE

energies (>10
20 kcal/mol) the error is in the range of 2
4 kcal/ mol, but which is only about 10% of the binding energy. Deviations of about the same magnitude have been observed earlier in many benchmark studies on much smaller complexes. Aside from a small and consistent overbinding of DFT-D3 compared to the LPNOCEPA/1(CBS) results, no outliers are found and in particular the energetic sequence of the protein
ligand model structures is agreed upon in general by the different methods. Although in our study the most sophisticated wave function based methods for complexes of several hundreds of atoms have been applied, the remaining LPNO-CEPA/1(CBS) error easily amounts about 5
10% in the reference, too. From benchmark studies of the popular S22 set, we know that LPNO-CEPA/1(CBS) slightly underbinds compared to CCSD(T)/CBS, which perfectly fits to the comparative DFT-D3/LPNO-CEPA/1 results obtained here for the protein
ligand model complexes. The reliability of LPNOCEPA/1(CBS) is further supported by the sequence of the DFTD3 MADs (with respect to LPNO-CEPA/1(CBS)) for different functionals found in the present study, which is the same as that for the S22 benchmark7,101 with estimated CCSD(T)/CBS reference data.97,98 For DFT calculations, TZ basis sets deliver almost converged interaction energies and we see no reason to go beyond this level in typical applications. In recent years we have gathered substantial evidence that an arbitrary noncovalent interaction in a biochemical system calculated by DFT-D3/TZ is in error by typically only 5%, maximally by 10%. Often the dispersion energy contributes 50
100% to binding in the investigated protein
ligand complexes. Except for multiply charged systems this can be expected to hold in general although one has to keep in mind that our -D3 dispersion energy also contains a short- to medium-range contribution. In any case it is clear from the present study that methods which do not account for dispersion properly cannot be expected to be even qualitatively correct. With TZVP GGA based methods plus dispersion correction, much more and also larger structures can be calculated if desired. Thus, for example, reference data for benchmarking of fast scoring functions can be provided routinely. We propose DFTD3 as the standard method that perfectly fits for this purpose. Standard quantum chemical reference methods like the “gold standard” CCSD(T) are not applicable to systems of this size routinely. Other approaches suitable for large molecules such as (approximate) local MP2 can not be recommended due to the known deficiences of MP2 in modeling intermolecular interactions,107 its strong basis set dependence and the necessary, computationally demanding Hartree
Fock step, which is required for MP2. All in all, we recommend for future benchmark purposes the 15 estimated LPNO-CEPA/1(CBS) or alternatively the 24 “best” DFT (B97-D3/def2-TZVP) values.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT Financial support by the german research foundation (DFG) by Grant AN 793/1-1 is gratefully acknowledged. F.N. and D.G.L. greatly appreciate financial support of this work by the SFB 624. We also thank Ms. Ute Becker for skillful technical assistance in the development of the parallelized LPNO-CEPA method. 11218

dx.doi.org/10.1021/jp203963f |J. Phys. Chem. A 2011, 115, 11210–11220

The Journal of Physical Chemistry A

’ REFERENCES (1) Reddy, A. S.; Pati, S. P.; Kumar, P. P.; Pradeep, H. N.; Sastry, G. N. Curr. Protein Pept. Sci. 2007, 8 (4), 329–351. (2) Cavasotto, C. N.; Orry, A. J.; Abagyan, R. A. Curr. Comput.-Aided Drug Des. 2005, 1 (4), 423–440. (3) Kellogg, G. E.; Fornabaio, M.; Spyrakis, F.; Lodola, A.; Cozzini, P.; Mozzarelli, A.; Abraham, D. J. J. Mol. Graphics Modell. 2004, 22 (6), 479–486. (4) Gilson, M. K.; Zhou, H.-X. Annu. Rev. Biophys. Biomol. Struct. 2007, 36, 21–42. (5) Wang, R.; Lu, Y.; Wang, S. J. Med. Chem. 2003, 46 (12), 2287–2303. (6) Huang, S.-Y.; Grinter, S. Z.; Zou, X. Phys. Chem. Chem. Phys. 2010, 12, 12899–12908. (7) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. J. Chem. Phys. 2010, 132, 154104. (8) Grimme, S. WIREs Comput. Mol. Sci. 2011, 1, 211–228. € (9) Eichkorn, K.; Treutler, O.; Ohm, H.; H€aser, M.; Ahlrichs, R. Chem. Phys. Lett. 1995, 240 (4), 283–289. € (10) Eichkorn, K.; Treutler, O.; Ohm, H.; H€aser, M.; Ahlrichs, R. Chem. Phys. Lett. 1995, 242 (6), 652–660. (11) Peters, M. B.; Raha, K.; Merz, K. M., Jr. Curr. Opin. Drug Discovery Dev. 2006, 9 (3), 370–379. (12) Raha, K.; Peters, M. B.; Wang, B.; Yu, N.; Wallacott, A. M.; Westerhoff, L. M.; Merz, K. M., Jr. Drug Discovery Today 2007, 12, 725–731. (13) S€oderhjelm, P.; Kongstedt, J.; Genheden, S.; Ryde, U. Interdiscip. Sci. Comput. Life Sci. 2010, 2, 21–37. (14) Gohlke, H.; Klebe, G. Angew. Chem., Int. Ed. 2002, 41, 2644–2676. (15) Raha, K.; Merz, K. M., Jr. J. Am. Chem. Soc. 2004, 126 (4), 1020–1021. (16) Raha, K.; Merz, K. M., Jr. J. Med. Chem. 2005, 48 (14), 4558–4575. (17) Hayik, S. A.; Dunbrack, R., Jr.; Merz, K. M., Jr. J. Chem. Theory Comput. 2010, 6, 3079–3091. (18) Peters, M. B.; Merz, K. M., Jr. J. Chem. Theory Comput. 2006, 2, 383–399. (19) Raha, K.; van der Vaart, A. J.; Riley, K. E.; Peters, M. B.; Westerhoff, L. M.; Kim, H.; Merz, K. M., Jr. J. Am. Chem. Soc. 2005, 127 (18), 6583–6594. (20) Merz, K. M., Jr. J. Chem. Theory Comput. 2010, 6, 1769–1776. (21) Faver, J. C.; Benson, M. L.; He, X.; Roberts, B. P.; Wang, B.; Marshall, M. S.; Kennedy, M. R.; Sherrill, C. D.; Merz, K. M., Jr. J. Chem. Theory Comput. 2011, 7, 790–797. (22) Dobes, P.; Otyepka, M.; Strnad, M.; Hobza, P. Chem.—Eur. J. 2006, 12, 4297–4304. (23) Vasilyev, V.; Bliznyuk, A. Theor. Chem. Acc. 2004, 112 (4), 313–317. (24) S€oderhjelm, P.; Kongstedt, J.; Ryde, U. J. Chem. Theory Comput. 2010, 6, 1726–1737. (25) S€oderhjelm, P.; Ryde, U. J. Phys. Chem. A 2009, 113, 617–627. (26) S€oderhjelm, P.; Aquilante, F.; Ryde, U. J. Phys Chem. B 2009, 113, 11085–11094. (27) Duan, L. L.; Tong, Y.; Mei, Y.; Zhang, Q. G.; Zhang, J. Z. H. J. Chem. Phys. 2007, 127, 145101. (28) Ishikawa, T.; Ishikura, T.; Kuwata, K. J. Comput. Chem. 2009, 30, 2594–2601. (29) Fanfrlík, J.; Bronowska, A. K.; Rezac, J.; Prenosil, O.; Konvalinka, J.; Hobza, P. J. Phys. Chem. B 2010, 114 (39), 12666–12678. (30) Beierlein, F.; Lanig, H.; Sch€urer, G.; Horn, A. H. C.; Clark, T. Mol. Phys. 2003, 101 (15), 2469–2480. (31) Gr€ater, F.; Schwarzl, S. M.; Dejaegere, A.; Fischer, S.; Smith, J. C. J. Phys. Chem. B 2005, 109 (20), 10474–10483. (32) Khandelwal, A.; Lukacova, V.; Comez, D.; Kroll, D. M.; Raha, S.; Balaz, S. J. Med. Chem. 2005, 48 (17), 5437–5447. (33) Lameira, J.; Alves, C. N.; Moliner, V.; Martí, S.; Kanaan, N.; Tunon, I. J. Phys. Chem. B 2008, 112 (45), 14260–14266.

ARTICLE

(34) Senn, H. M.; Thiel, W. Top. Curr. Chem. 2007, 268, 173–290. (35) Wittayanarakul, K.; Aruksakunwong, O.; Saen-oon, S.; Chantratita, W.; Parasuk, V.; Sompornpisut, P.; Hannongbua, S. Biophys. J. 2005, 88 (2), 867–879. (36) Fanfrlík, J.; Brynda, J.; Rezac, J.; Hobza, P.; Lepsík, M. J. Phys. Chem. B 2008, 112 (47), 15094–15102. (37) Gleeson, M. P.; Gleeson, D. J. Chem. Inf. Model 2009, 49, 670–677. (38) Cho, A. E.; Chung, J. Y.; Kim, M.; Park, K. J. Chem. Phys. 2009, 131, 134108. (39) Cho, A. E.; Rinaldo, D. J. Comput. Chem. 2009, 30, 2609–2616. (40) Morgado, C.; Hillier, I. H.; Burton, N. A.; McDouall, J. J. W. Phys. Chem. Chem. Phys. 2008, 10, 2706–2714. (41) Riley, K. E.; Hobza, P. J. Phys. Chem. B 2008, 112 (10), 3157–3163. (42) Neese, F.; Wennmohs, F.; Hansen, A. J. Chem. Phys. 2009, 130, 114108. (43) Neese, F.; Hansen, A.; Liakos, D. G. J. Chem. Phys. 2009, 131, 064103. (44) Liakos, D. G.; Hansen, A.; Neese, F. J. Chem. Theory Comput. 2011, 7, 76–87. (45) Neese, F.; Hansen, A.; Wennmohs, F.; Grimme, S. Acc. Chem. Res. 2009, 42, 641–648. (46) London, F. Z. Physik. 1930, 60, 245
279. (47) Grimme, S.; Antony, J.; Schwabe, T.; M€uck-Lichtenfeld, C. Org. Biomol. Chem. 2007, 5, 741–758.  .; Sumpter, B. G.; Sherrill, (48) Burns, L. A.; Vazquez-Mayagoitia, A C. D. J. Chem. Phys. 2011, 134, 084107. (49) Antony, J.; Grimme, S. Phys. Chem. Chem. Phys. 2006, 8, 5287–5293. (50) Grimme, S.; M€uck-Lichtenfeld, C.; Antony, J. J. Phys. Chem. C 2007, 111 (30), 11199–11207. (51) Antony, J.; Grimme, S. Phys. Chem. Chem. Phys. 2008, 10, 2722–2729. (52) Zou, B.; Dreger, K.; M€uck-Lichtenfeld, C.; Grimme, S.; Sch€afer, H. J.; Fuchs, H.; Chi, L. Langmuir 2005, 21, 1364–1370. (53) Piacenza, M.; Grimme, S. Chem. Phys. Chem. 2005, 6, 1554–1558. (54) Piacenza, M.; Grimme, S. J. Am. Chem. Soc. 2005, 127 (42), 14841–14848. (55) Parac, M.; Etinski, M.; Peric, M.; Grimme, S. J. Chem. Theory Comput. 2005, 1, 1110–1118. (56) M€uck-Lichtenfeld, C.; Grimme, S. Mol. Phys. 2007, 105 (19), 2793–2798. (57) H€aber, T.; Seefeld, K.; Engler, G.; Grimme, S.; Kleinermanns, K. Phys. Chem. Chem. Phys. 2008, 10, 2844–2851. (58) Ducere, J.-M.; Cavallo, L. J. Phys. Chem. B 2007, 111 (45), 13124–13134. (59) Hayama, T.; Baldridge, K. K.; Wu, Y.-T.; Linden, A.; Siegel, J. S. J. Am. Chem. Soc. 2008, 130 (5), 1583–1591. (60) Jurecka, P.; Cerny, J.; Hobza, P.; Salahub, D. J. Comput. Chem. 2007, 28 (2), 555–569. (61) Morgado, C.; Vincent, M. A.; Hillier, I. H.; Shan, X. Phys. Chem. Chem. Phys. 2007, 9, 448–451. (62) Pankewitz, T.; Klopper, W. J. Phys. Chem. C 2007, 111 (51), 18917–18926. (63) Barone, V.; Casarin, M.; Forrer, D.; Pavone, M.; Sambi, M.; Vittadini, A. J. Comput. Chem. 2009, 30 (6), 934–939. (64) Kerber, T.; Sierka, M.; Sauer, J. J. Comput. Chem. 2008, 29 (13), 2088–2097. (65) Schmidt, W. G.; Seino, K.; Preuss, M.; Hermann, A.; Ortmann, F.; Bechstedt, F. Appl. Phys. A: Mater. Sci. Process. 2006, 85 (4), 387–397. (66) Grimme, S. J. Comput. Chem. 2004, 25 (12), 1463–1473. (67) Grimme, S. J. Comput. Chem. 2006, 27 (15), 1787–1799. (68) Perdew, J. P. Phys. Rev. B 1986, 33 (12), 8822–8824. (69) Perdew, J. P.; Wang, Y. Phys. Rev. B 1992, 45 (23), 13244–13249. (70) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77 (18), 3865–3868. 11219

dx.doi.org/10.1021/jp203963f |J. Phys. Chem. A 2011, 115, 11210–11220

The Journal of Physical Chemistry A (71) Becke, A. D. Phys. Rev. A 1988, 38 (6), 3098–3100. (72) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37 (2), 785–789. (73) Becke, A. D. J. Chem. Phys. 1993, 98, 5648–5652. (74) Stephens, P. J.; Devlin, F. J.; Chablowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623–11627. (75) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Phys. Rev. Lett. 2003, 91, 146401. (76) Sch€afer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97 (4), 2571–2577. (77) Sch€afer, A.; Huber, C.; Ahlrichs, R. J. Chem. Phys. 1994, 100 (8), 5829–5835. (78) Weigend, F.; Furche, F.; Ahlrichs, R. J. Chem. Phys. 2003, 119 (24), 12753–12762. (79) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19 (4), 553–566. (80) Weigend, F.; H€aser, M.; Patzelt, H.; Ahlrichs, R. Chem. Phys. Lett. 1998, 294 (1
3), 143–152. (81) Turbomole, Version 6.0, Program Package for ab initio Electronic Structure Calculations; TURBOMOLE GmbH, Karlsruhe, Germany, 2009; http://www.cosmologic.de/QuantumChemistry/ main_qChemistry.html. (82) Neese, F. Orca, an ab initio, DFT and semiempirical SCF-MO package, Version 2.7-00; Chair of Theoretical Chemistry, University of Bonn, Germany, 2009. (83) Stewart, J. J. P. MOPAC2009, Molecular Orbital PACkage; Version 10.040L; Stewart Computational Chemistry, Colorado Springs, CO, 2009; http://openmopac.net. (84) Frauenheim, T. DFTB+, Density Functional based Tight Binding; Release: 1.0; Bremen Center for Computational Materials Science, University of Bremen, Germany, 2008; http://www.dftb.org/. (85) Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215–241. (86) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; , Fox, D. J. Gaussian 09, Revision A.02, Gaussian, Inc.: Wallingford, CT, 2009. (87) Block, P.; Sotriffer, C. A.; Dramburg, I.; Klebe, G. Nucleic Acids Res. 2006, 34, D522–D526. (88) Liu, T.; Lin, Y.; Wen, X.; Jorissen, R. N.; Gilson, M. K. Nucleic Acids Res. 2007, 35, D198–D201. (89) Benson, M. L.; Smith, R. D.; Khazanov, N. A.; Dimcheff, B.; Beaver, J.; Dresslar, P.; Nerothin, J.; Carlson, H. A. Nucleic Acids Res. 2008, 36, D674–D678. (90) Wallach, I.; Lilien, R. Bioinformatics 2009, 25 (5), 615–620. (91) Wang, R.; Fang, X.; Lu, Y.; Yang, C.-Y.; Wang, S. J. Med. Chem. 2005, 48 (12), 4111–4119. (92) Perdew, J. P. Phys. Rev. B 1986, 34 (10), 7406. (93) Halkier, A.; Klopper, W.; Helgaker, T.; Jørgensen, P.; Taylor, P. R. J. Chem. Phys. 1999, 111 (20), 9157–9167. (94) Tarakeshwar, P.; Choi, H. S.; Lee, S. J.; Lee, J. Y.; Kim, K. S.; Ha, T.-K.; Jang, J. H.; Lee, J. G.; Lee, H. J. Chem. Phys. 1999, 111 (13), 5838–5850. (95) Kim, K. S.; Tarakeshwar, P.; Lee, J. Y. Chem. Rev. 2000, 100, 4145–4186. (96) Valdes, H.; Sordo, J. A. J. Phys. Chem. A 2002, 106 (15), 3690–3701. (97) Jurecka, P.; Sponer, J.; Cerny, J.; Hobza, P. Phys. Chem. Chem. Phys. 2006, 8, 1985–1993.

ARTICLE

(98) Takatani, T.; Hohenstein, E. G.; Malagoli, M.; Marshall, M. S.; Sherrill, C. D. J. Chem. Phys. 2010, 132, 144104. (99) Liakos, D. G.; Neese, F. J. Phys. Chem. A 2011, submitted for publication. (100) Antony, J.; Grimme, S. J. Phys. Chem. A 2007, 111 (22), 4862–4868. (101) Grimme, S.; Ehrlich, S.; Goerigk, L. J. Comput. Chem. 2011, 32, 1456–1465. (102) Korth, M.; Pito nak, M.; Rezac, J.; Hobza, P. J. Chem. Theor. Comput. 2010, 6, 344–352. (103) Rezac, J.; Fanfrlik, J.; Salahub, D.; Hobza, P. J. Chem. Theor. Comput. 2009, 5, 1749–1760. (104) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Phys. Rev. B 1998, 58 (11), 7260–7268. (105) Elstner, M.; Hobza, P.; Frauenheim, T.; Suhai, S.; Kaxiras, E. J. Chem. Phys. 2001, 114 (12), 5149–5155. (106) Schaftenaar, G.; Noordik, J. H. J. Comput.-Aided Mol. Des. 2000, 14, 123–134. (107) Tsuzuki, S.; Honda, K.; Uchimura, T.; Mikami, M. J. Chem. Phys. 2004, 120, 647–659. (108) Structure 23 has been omitted due to an O
O distance