CHEMPHYSCHEM ARTICLES DOI: 10.1002/cphc.201402379

Coupled-Cluster Interaction Energies for 200-Atom Host– Guest Systems Milica Andrejic´,[b] Ulf Ryde,[c] Ricardo A. Mata,[b] and Pr Sçderhjelm*[a] We have developed a method to calculate interaction energies of large systems (such as host–guest or even protein–ligand systems) at the local coupled-cluster with singles, doubles, and perturbative triples level, and with extrapolation to the limit of a complete basis set. The method is based on the polarizable multipole interactions with supermolecular pairs molecular fractionation approach, which combines a pairwise quantummechanical evaluation of the short-range interactions with a polarizable multipole treatment of many-body effects. The

method is tested for nine guest molecules binding to an octaacid host (in total 198–207 atoms), as part of the SAMPL4 blind challenge. From the test calculations, the accuracy of the approach is found to be 10 kJ mol1 or better. Comparison with dispersion-corrected density functional theory reveals that the latter underestimates the dispersion contribution for this type of system, which leads to a difference in the ranking of the ligands.

1. Introduction An accurate prediction of binding free energies between proteins and ligands is often hindered by inadequate sampling of all relevant geometries, as well as deficiencies in the evaluation of the interaction energy at a given geometry. These problems are interconnected in the sense that, if much sampling is required, one normally has to resort to simple molecular mechanics (MM) force fields, whereas if an accurate energy evaluation by using quantum mechanical (QM) methods is required, one can normally study only a limited number of conformations. For systematic method development, it is therefore of great interest to investigate artificial, smaller model systems, usually called host-guest complexes. These systems involve the same principles and types of interactions as protein–ligand systems, but both problems are easier; the number of conformations is smaller owing to the rigidity of the host, and the system itself is smaller, reducing the computational cost of accurate energy evaluations. However, not even for host–guest systems with around 200 atoms can the interaction energy be evaluated at the coupledcluster singles, doubles, and perturbative triples [CCSD(T)] level, which is currently considered currently to be the gold [a] Dr. P. Sçderhjelm Department of Biophysical Chemistry, Lund University Chemical Center, P.O. Box 124 SE-221 00 Lund (Sweden) E-mail: [email protected] [b] M. Andrejic´, Prof. R. A. Mata Institut fr Physikalische Chemie, Universitt Gçttingen Tammannstrasse 6, 37077 Gçttingen (Germany) [c] Prof. U. Ryde Department of Theoretical Chemistry, Lund University Chemical Center, P.O. Box. 124 SE-221 00 Lund (Sweden) Supporting Information for this article is available on the WWW under http://dx.doi.org/10.1002/cphc.201402379.

 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

standard of QM calculations with an a priori reliability and accuracy in the order of 4 kJ mol1 for many interacting systems.[1] This is unfortunate, because even if there exist QM methods that can handle these system sizes and are often almost as accurate as coupled-cluster methods [primarily dispersion-corrected density functional theory (DFT)],[2] there is no guarantee that they will be sufficiently accurate for previously untested systems. This is because they typically contain an empirical component that has been fitted to the systems of interest. Much effort has been spent on reducing the high scaling of accurate quantum-chemical methods with respect to system size. For the post-Hartree–Fock (HF) family of methods, such as second-order Møller–Plesset perturbation theory (MP2) and coupled-cluster theory, a promising alternative is to use local correlation treatments. Owing to its local nature, dynamical correlation is best described through localized orbital spaces. The steep scaling of conventional wavefunction methods is linked strongly to the use of canonical orbitals, which are not amenable to approximations of this kind. Several alternatives have been suggested over the years, differing mainly in the way the virtual space is defined.[3–7] The use of projected atomic orbitals (PAOs) has been a common choice. However, the use of pair natural orbitals has become increasingly popular and can lead to pronounced savings, as has been shown recently.[8] Nevertheless, routine calculations with such methods are still not common, and limited usually to systems below 100 atoms. This is due not only to the late onset of the linearscaling regime of some of the methods, but also to the underlying HF calculations, which can become quite demanding for more than 3000 basis functions. In fact, it has been shown to dominate the computational cost already for relatively small systems.[8–10]

ChemPhysChem 2014, 15, 3270 – 3281

3270

CHEMPHYSCHEM ARTICLES Another solution to the scaling problem is given by molecular fractionation methods.[11, 12] These work by dividing the system into fragments and performing QM calculations on only a subset of the fragments at a time. Applied to the interaction energy between a central molecule and a set of surrounding fragments, the simplest approach would be to approximate the interaction by a sum of interaction energies with the individual fragments (for example, the pairwise additive method in Ref. [13]). An extension of this idea, to the case in which the surrounding fragments are connected by covalent bonds, is the molecular fractionation with conjugate caps (MFCC) method,[14] which uses an alternating sum of positive and negative terms to avoid double-counting any interactions. The disadvantage of these methods is that they assume pairwise additivity of the various fragment interactions or in other words, they neglect the many-body effects of the system. The inclusion of many-body effects in fractionation methods has been addressed in many different ways.[11] One of the most used methods is the fragment molecular orbital (FMO) method,[15] which uses an embedding in the exact electron density of the surrounding fragments. Other methods, such as the electrostatically embedded pairwise additive model,[13] the binary interaction method,[16] the generalized energy-based fragmentation (GEBF) method,[17] and the electrostatic fieldadapted (EFA) MFCC method[18] rely on point charges to account for global many-body effects. A fragmentation method may also involve several layers of QM methods,[19] for example, a combination of CCSD(T) and MP2.[20] The polarizable multipole interactions with supermolecular pairs (PMISP) method[21] was introduced to address the manybody effects in a simple, yet accurate way, through the use of classical expressions for the interactions between multipoles and polarizabilities. By keeping the treatment of many-body effects separate from the QM calculations, it is trivial to adapt the method to new QM methods and quantum chemistry packages. Conversely, the polarization model can also be replaced, as was done in the so-called effective fragment molecular orbital (EFMO) approach.[22] The PMISP method has been combined with a QM/MM treatment of long range interactions[23] and used to calculate binding free energies through an MM-Poisson–Boltzmann surface-area- (PBSA) like approach.[24] An advantage of the method is the connection to polarizable force fields; the method basically represents the limit of how accurate a force field can become if it relies on a given polarization model for treating many-body effects (for example, a nonpolarizable force field can never become more accurate than MFCC).[21] Thus, benchmarking a polarization model by using many-body effects (as will be done in the current paper) not only guarantees a good performance of PMISP, but also lays the basis for a useful polarizable force field. However, there might be other aspects of a polarizable force field that are also important,[25] such as the ability to adapt the electrostatic properties of a molecule depending on its conformation through intramolecular polarization.[26–29] In a recent study,[30] the sampling and force field problems were addressed simultaneously in predicting free energies for nine guest molecules binding to the octa-acid host. The study  2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

www.chemphyschem.org was part of the SAMPL4 challenge[31] and was set up as a blind test of several procedures combining quantum-chemical calculations with either advanced sampling methods or energy minimization. Among the tested methods was a new way of calculating interaction energies at the local CCSD(T0) [LCCSD(T0)] level extrapolated to the complete basis set (CBS) limit, for structures geometry-optimized at the DFT level. To our knowledge, this was the first time interaction energies for systems of this size were computed at such a high level of theory. The purpose of the current article is to describe this approach in detail and to test thoroughly its theoretical grounds. In particular, we introduce and test an improved treatment of intramolecular polarization and investigate the reliability of the PMISP method with respect to changes in basis set. Having established the accuracy of the method, we use the results to discuss the accuracy of other QM methods for this type of system. Finally, the interaction energies are combined with various correction terms to allow a comparison with experimental binding free energies.

Methods The PMISP Fractionation Method The PMISP method for calculating the interaction energy between two molecules A and B is based on a spatial division of one or both molecules into several fragments and calculating each fragment–fragment interaction energy with an accurate quantumchemical method. The method can be described as a combination of the MFCC method, which assumes pairwise additivity of the fragment–fragment interaction energies, and a many-body term, which takes into account the deviation from pairwise additivity and is estimated by a polarizable multipole interaction (PMI) model [Eq. (1)]: PMISP MFCC MB EAB ¼ EAB þ EAB

ð1Þ

in which [Eq. (2)]: MFCC EAB ¼

n X

ci EAsup iB

ð2Þ

i¼1

and [Eq. (3)]: MB AB

E

¼

E

PMI AB



n X

! PMI i Ai B

ð3Þ

cE

i¼1

In these equations, we assumed that only molecule A is fragmented, to simplify the notation. The coefficients ci are either 1 for a normal fragment or 1 for a conjugate cap (concap) fragment, the contribution of which should be subtracted to avoid double counting. The concap fragments are normally the overlapping parts between two adjacent normal fragments; however, in general, any choice of fragments and ci coefficients is allowed, as long as P the relation nk ðAÞ ¼ i ci nk ðAi Þ holds for every possible atom k, in which nk(X) equals 1 if atom k occurred in molecule X, and 0 otherwise. All energies EXY are rigid interaction energies between molecules X and Y, that is, assuming the same monomer geometries as sup , are in the complex. All supermolecular interaction energies, EXY ChemPhysChem 2014, 15, 3270 – 3281

3271

CHEMPHYSCHEM ARTICLES

www.chemphyschem.org

corrected for the basis-set superposition error with the counterpoise method.[32] The PMI model was constructed specifically for the actual geometry by performing several quantum-chemical calculations on each fragment (unperturbed and perturbed by an electric field in three orthogonal directions). The LoProp method[33] was used to create a classical description of each fragment in terms of distributed multipoles (charges, dipoles, quadrupoles, and octupoles) and anisotropic polarizabilities centered in the atomic nuclei and the midpoints of covalent bonds. The multipoles and polarizabilities for the whole molecule A were then obtained through a summation analogous to that of the energies [Eq. (4)]: Pk ðAÞ ¼

X

ci Pk ðAi Þ

ð4Þ

i

in which Pk is a multipole or polarizability component at center k. These “raw” multipoles need to be polarized by the surrounding fragments, because the multipoles Pk(Ai ) of each fragment are obtained in the absence of all other fragments. For conceptual reaPMI into electrostatic and insons (to enable a natural division of EAB duction contributions),[21] this step was done before introducing the molecule B, that is, as a prepolarization, and thus acted as a correction applied at the dipole level [Eq. (5)]: ind mfinal k ðAÞ ¼ mk ðAÞ þ mk ðAÞ

ð5Þ

in which mk(A) is the raw dipole moment from Equation (4), and mind k ðAÞ is the self-consistent induced dipole from an intramolecular polarization model, which will be presented in detail in the next section. Owing to the prepolarization step, the electric field used PMI does not include confor polarization in the final calculation of EAB tributions from multipoles within the same molecule. In contrast, the coupling between polarizabilities is always present (as will also be explained in the next section).

Treatment of Intramolecular Polarization From Equations (3) and (4), it can be seen that the contribution from any pair of raw multipoles, one in A and one in B, cancels exactly, so that only contributions from the polarization survive. Thus, hypothetically, if the PMI model was replaced by a pairwise addiMB PMISP MFCC ¼ 0 by construction, and EAB ¼ EAB , tive force field, then EAB [14] that is, the original MFCC method would be recovered. Even if the many-body effect is caused by polarization, it can be divided into an electrostatic contribution, which comes from the modified electrostatic properties of molecule A in the prepolarization step [Eq. (5)], and an induction contribution, which is caused by the nonadditivity of the induction energy itself [Eq. (6)]: MB AB

E

MB;ele AB

¼E

MB;ind AB

þE

¼

ele AB

E



n X

! ele i Ai B

cE

þ

ind AB

E

i¼1



n X

! ind i Ai B

cE

i¼1

ð6Þ The electrostatic contribution is often the most important and long-ranged one.[25] Moreover, this contribution depends on the treatment of polarization interactions within molecule A (intramolecular polarization), which are necessarily modeled more approximately than the polarization between different molecules. In particular, to avoid polarization catastrophes and obtain a physical behavior of the model, a polarizable force field normally includes a rule for which interactions are neglected or damped (typically be 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

tween close-lying centers). If one aims at high accuracy, special care is required on devising this rule. In the original PMISP method,[21] we derived the following rule for the static polarization (polarization caused by the static multipoles): A polarizable site m should not feel the electric field from a multipole site k if they have belonged to the same fragment in at least one LoProp calculation, that is, if nk(Ai) = nm(Ai ) = 1 for at least one i. However, the motivation of this rule (see Ref. [21]) depends on an important assumption, namely that if the pair (k,m) occurs in two fragments with ci = 1 and ci = 1, the effect on m is cancelled. This assumption holds to a good approximation if the truncation of bonds (and satisfying the dangling bonds with H atoms) occurs only at CC bonds, as in previous applications of PMISP. For the systems in the current study, however, truncations of CO bonds are unavoidable. To solve this problem, we implemented a new, more precise rule for treating the intramolecular polarization. Formally, for every polarizable site m, a set of multipoles A(m) representing the difference in environment between the full molecule and the fragment is constructed as follows [Eq. (7)]: X

nm ðAi Þ¼0

Pk ðAðmÞ Þ ¼

ci Pk ðAi Þ

ð7Þ

i

in which the summation is restricted to the fragments that do not include m. In the A(m) set, there might be nonzero multipoles positioned close to m. They tend to be small in magnitude if the fractionation scheme is chemically reasonable, but they can still cause numerical problems and unphysical polarization. Thus, all multipoles within a distance of R1 from site m are deleted from A(m) and the vector sum of the removed multipoles is spread out evenly on all sites for which the distance from m lies in the interval (R1,R2). Thus, rather than omitting the charges and multipoles they are translocated to neighboring multipole sites. We tested various values of R1 and R2, as shown in Table S2 in the Supporting Information, and found the results to be insensitive to this choice. Thus, we used the thresholds R1 = 0.53  and R2 = 1.59  throughout this study (i.e. 1 and 3 Bohr). After possible multipole translocations, the electric field from the set of multipoles A(m) at site m is computed and the procedure is repeated for every polarizable site m. Finally, this field is used as the static field in the prepolarization, that is, a self-consistent calculation of all induced dipoles mind k ðAÞ in molecule A, to be used in Equation (5). This procedure is illustrated with the example in Figure 1, which shows the assembly of multipoles from fragment calculations and the new treatment of intramolecular polarization for a simple molecule (1,3-diphenoxybenzene) bearing some resemblance to the octa-acid host. For simplicity, only charges (zeroth order multipoles) are shown, in hundredths of the elementary charge e. This avoids consideration of the bond centres (because LoProp charges in the bonds are zero by construction), but the same principle can be applied for higher-order multipoles and bond centres. The charges were taken from HF/6-31G* calculations, but to avoid cluttering, the charge of each “normal” hydrogen atom has been absorbed into its bound carbon atom. However, the link hydrogen atoms that replace a truncated part of the molecule are explicitly shown. The left part of Figure 1 shows the atomic charges of the three fragments and the “raw” charges resulting from combining these charges according to Equation (4). This is a rather inaccurate description of the target molecule, because the charges of each fragChemPhysChem 2014, 15, 3270 – 3281

3272

CHEMPHYSCHEM ARTICLES

www.chemphyschem.org

Figure 1. An example showing the assembly of multipoles for 1,3-diphenoxybenzene from fragment calculations and the new treatment of intramolecular polarization. Charges are given in hundredths of the elementary charge e. The left part of the Figure shows the atomic charges of the three fragments and the “raw” charges resulting from combining these charges according to Equation (4). The right part of the Figure illustrates how intramolecular polarization is applied to account for the influence of neighboring fragments on the multipoles.

ment have been calculated without the presence of the other fragments. The right part of Figure 1 shows how intramolecular polarization was applied to remedy this problem. In particular, we illustrate how the static electric field at any of the polarizable sites of the leftmost phenoxy moiety (marked with a) was computed in the current study (new rule) and how it would have been computed in previous PMISP studies (old rule). With the new rule, an artificial set of charges A(a) = A2A3 was created by subtracting the charges of the third fragment from those of the second fragment. Note that this set does not correspond to a real molecule (for example it contains a subtracted “antihydrogen” atom marked by X in Figure 1) but simply represents the set of charges required to polarize the a sites as to account for the difference in environments between the full molecule and the first fragment. Note also that the H-marked charge closest to the polarizable sites was small in magnitude (this holds for any reasonable fractionation scheme), but it was nevertheless unphysically close to the oxygen atom and thus needed to be translocated when calculating the field at the oxygen atom (but not at the other sites marked with a). With the old rule, the set A(a) did not include the atoms of the central benzene ring because they belonged to a common fragment A1 (or, in other words, perfect cancellation between A2 and A3 was assumed for these atoms). Thus, no translocation was needed, but the net charge of 0.21 e gave rise to a systematic overestimation of the magnitude of the electric field at the a sites and thereby led to an inaccurate description of the whole molecule.

the new rule is more reliable. There might exist simpler rules that are sufficiently accurate, but our major concern was to apply a rule that accounted rigorously for the difference in environments between the full molecule and the fragment. The computational cost of the classical calculations was at any rate negligible. A similar discussion in the context of flexible molecules can be found in Ref. [26].

For “well-behaved” truncations (e.g. breaking aliphatic CC bonds), the results do not depend significantly on which rule is applied, but, as illustrated by the example, for truncations of CO bonds,

Tij ¼ dij

 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

The polarizability coupling can be treated by applying the same exclusion rules as for the static polarization, as is done for example, in the Amber 2002 force field.[34] Alternatively, it can be treated by a more general rule including coupling between all atom pairs but damping the short-distance interactions, as is done for example, in the atomic multipole optimized energetics for biomolecular simulation (AMOEBA) force field.[35] In the original PMISP model, we followed the former approach and used the rule stating that polarizabilities that had been in at least one common fragment calculation do not couple with each other. However, with the new rule for static polarization, which is even more specialized for getting the static polarization correct, it was more appropriate to treat polarizability coupling in a more straightforward way based on damping. We thus applied a Thole damping expression,[36] in which the dipole interaction tensor Tij between two isotropic polarizabilities a1 and a2 at a distance r was modified into [Eq. (8)]:

4u3  3u4 3xi xj u  5 r3 r

ChemPhysChem 2014, 15, 3270 – 3281

ð8Þ

3273

CHEMPHYSCHEM ARTICLES

www.chemphyschem.org

in which atomic units are employed and [Eq. (9)]:  u ¼ min 1;

r 2:3268ða1 a2 Þ1=6

 ð9Þ

This equation was developed for isotropic polarizabilities, but we used the same expression also for anisotropic polarizabilities, with the modification that both ai were replaced by the largest eigenvalue of each polarizability tensor. Using instead the average eigenvalue led to polarization catastrophe in some cases, because the LoProp polarizabilities can be strongly anisotropic and the interaction could occur in the most polarizable direction.

Computational Details The supermolecular calculations were performed by using MOLPRO.[37] The orbital basis sets used were the Dunning correlation-consistent polarized valence basis sets of double and triple zeta quality (cc-pVDZ and cc-pVTZ),[38, 39] the former for benchmark runs. The basis set dependence was corrected based on extrapolated results by using augmented aug-cc-pVTZ and aug-cc-pVQZ basis sets. All wave function calculations were carried out with density fitting approximations. The default Coulumb/exchange (JKFIT) and correlation (MP2FIT) fitting basis sets were used.[40] In the local calculations, we used occupied Pipek–Mezey local orbitals,[41] with the domains selected according to the natural population analysis (NPA) criteria (TNPA = 0.03).[42] In the LCCSD(T0) calculations, only a number of (close-lying) pairs were included fully in the coupled cluster treatment, with a large majority of the orbitals treated at the LMP2 level. This was controlled by distance criteria,[3, 43] which determined the pair classes. To converge the coupled cluster results, we used Rc = 3 and Rw = 5 Bohr, close and weak pair distances, respectively. These values have been applied previously in weakly interacting systems with good results.[44] A common issue in PAO-based local calculations is the stability of the domain definition. This is particularly important in fragment calculations, in which errors can add up as one builds a total energy out of the sum of many smaller calculations. To keep the procedure as automated as possible, and given that there are several extended p-bonds in the host system, all domains sharing more than one center were merged. This somewhat increased the computational cost of each individual calculation, but was a straightforward procedure to decrease the impact of small fluctuations in the orbital charge distributions. The complete basis-set limit of the supermolecular interaction energies was obtained through n3 extrapolation at the MP2 level (without local approximations) by using two calculations with diffuse basis sets, aug-cc-pVTZ and aug-cc-pVQZ, for each fragment– guest pair. This basis-set correction was added to the LCCSD(T0) results obtained with the cc-pVTZ basis set (LCC/TZ) [Eq. (10)]: EAsup ðLCC=CBSÞ ¼ EAsup ðLCC=TZÞ  EAsup ðLMP2=TZÞ þ EAsup ðMP2=CBSÞ iB iB iB iB ð10Þ By construction, the same relation holds for the corresponding MFCC and PMISP estimates of the full EAB interaction energy. Multipoles and polarizabilities at the HF level were computed by using MOLCAS[45] through the LoProp method.[33] Multipoles and polarizabilities at the LMP2 level were computed by using MOLPRO.[37] The response of the wave function was taken into account by using the so-called response density (or variational density),[46] which was then fed into MOLCAS for LoProp analysis. All  2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

DFT results were obtained by using Turbomole 6.4.[47, 48] Turbomole was also used for the COSMO-RS calculations[49–51] (see below), whereas MOLPRO was used for the LMP2/COSMO calculations.[52]

Systems As a test case, we used nine guest molecules binding to the octaacid host, that is, the octa-acid part of the SAMPL4 challenge.[31] The octa-acid host is shown in the top left corner of Figure 2 and the guest molecules are shown in Figure 3. The host–guest complexes were geometry-optimized at the TPSS-D3/def2-SV(P) level in COSMO(e=80) implicit solvent, as described in our earlier study.[30] The optimized structures are shown in Figure 7 b of Ref. [30] and are available in XYZ format in the Supporting Information. Earlier applications of PMISP have dealt with peptide chains and relied on a standard way of fractionation, in which each capped amino acid constitutes one fragment and all concaps are identical.[21] For the current systems, the fractionation was instead done manually. The fractionation scheme is shown schematically in Figure 4. The fractionation led to 24 “normal” fragments with ci = 1 (encircled in Figure 4), 24 concap fragments with ci = 1, 4 fragments with ci = 2, and 4 supporting fragments with ci = 1. The need for the last two categories arose from the complex pattern of overlapping fragments due to the branched structure of the octaacid molecule. A complete list of fragments, and the atoms included in each of them, is given in Table S1. For testing the PMISP method, we also defined three artificial systems (see Figure 2), which contained only a subset of the atoms of the full host as well as hydrogen atoms capping the dangling bonds.

Free Energies To the vacuum energy we added three correction terms to obtain an estimate of the binding free energy, which can be compared to experimental results [Eq. (11)]:[30] PMISP DGcalc ðLCC=CBSÞ þ DGsolv þ DGthermal þ DEstrain bind ¼ EAB

ð11Þ

in which DGsolv is the change in solvation energy between the complex and isolated monomers, estimated by using two calculations at the LMP2/cc-pVDZ level. One estimation was in vacuum and one in COSMO implicit solvent with e = 80, and with the nonelectrostatic terms taken from the COSMO-RS method.[49–51] In the corresponding DFT results, the standard COSMO-RS method was used to determine DGsolv, with COSMO calculations at the BP86/ TZVP level and e = 1, as described in Ref. [30]. The LMP2 level was chosen to provide a better account of the electronic distribution, more consistent with the LCCSD(T0) level used for the vacuum interaction energy. DGthermal is the thermal correction to the Gibbs free energy at T = 298 K and 1 atm pressure (including zero-point vibrational energy, entropy and enthalpy corrections) obtained from vibrational frequencies calculated at the MM level (but the TPSS/def2-SV(P) geometry for the rotational contributions) by using the rigid-rotor harmonic-oscillator approximation.[30] Lowlying vibrational modes were treated by the free-rotor approximation, by using the interpolation model suggested by Grimme and w0 = 100 cm1.[52] Finally, DEstrain is the difference of the energy of the isolated guest evaluated at the optimized complex geometry and at the optimized guest geometry, respectively, calculated at the TPSS/def2-QZVP’ level on geometries optimized with the def2SV(P) basis set. The exact procedure for computing all these corrections has been detailed in the previous study.[30] As explained in that reference, we tried to estimate also the strain energy of the ChemPhysChem 2014, 15, 3270 – 3281

3274

CHEMPHYSCHEM ARTICLES

www.chemphyschem.org

Figure 2. The full octa-acid host and three smaller subsystems used to assess the accuracy of the PMISP method. The molecular charge is 8 for the full host and 4 for the ring, whereas the truncated host (with the propionate chains removed and the benzoates turned into benzene rings) and the truncated rings are uncharged. Hydrogen atoms are omitted for clarity.

Figure 3. The nine guest molecules.

host, but found severe problems with the host molecule existing in different local minima in the various optimized host–guest complexes (especially with respect to the orientation of the four propionate chains as well as a breathing motion of the host). Therefore

 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

Figure 4. Fractionation of the octa-acid host. The schematic picture shows one flattened side of the host; the same pattern repeats itself all around the cyclic molecule so that each chemical group appears four times. Only normal fragments with ci = 1 are shown; the concap fragments mainly represent the overlapping regions of the shown fragments. A list of the atoms included in each fragment is given in Table S1 in the Supporting Information.

we concluded that, for this system, the use of rigid interaction energies without considering the host strain is in fact more reliable.

ChemPhysChem 2014, 15, 3270 – 3281

3275

CHEMPHYSCHEM ARTICLES 2. Results and Discussion We have computed interaction energies for the octa-acid host–guest system with nine different ligands at the LCCSD(T0) level with the PMISP approach. We first investigate the accuracy of the various approximations involved in the calculations and then compare the results to experimental data, by using the three corrections in Equation (11).

2.1. Assessment of the PMISP Method The PMISP method has previously been used mainly for protein–ligand systems, in which the fractionation follows straightforwardly from the polymeric nature of the amino acid chain. In the octa-acid host, which is a much more branched molecule, the fragments are smaller, more overlapping, and contain a greater number of bond truncations than in the previous applications. Therefore, we need to assess the accuracy of the PMISP method for these systems.

2.1.1. Hartree–Fock Level The assessment is done in several steps. To scrutinize the performance of the PMI model without getting confused by electron-correlation effects, we first compare the PMISP energy at the HF/cc-pVDZ level to reference results obtained with the whole molecule treated in a single calculation. The net interaction energy is large in vacuum (  1000 kJ mol1), because it is dominated by the Coulomb interaction between the host and guest molecules with net charges of 8 and 1 e, respectively. MB , Therefore, we concentrate on the many-body effect, EAB which is compared to the many-body effect given by reference sup MFCC  EAB . The results are given in Table 1. The calculations, EAB mean absolute error over the nine guest molecules is 2.6 kJ mol1 (27 % of the average magnitude of the many-body effect), that is, similar to those obtained in a previous study,[21] despite the increased complexity of the fractionation in the

www.chemphyschem.org current study. It is particularly satisfying that the maximum error over the nine systems is only 5.7 kJ mol1. To further validate the stability of the method, we designed three artificial hosts (see Figure 2) by truncating the octa-acid molecule in various ways, thus varying the size and the charge (0, 4, and 8). For all truncations and all guest molecules, we calculated the reference energies and compared them with the PMISP estimates. Moreover, we computed the PMI representation of the full molecule and compared the host–guest electrostatic interaction energy predicted by this representation and the representation computed fragment-wise by Equation (5). Again we concentrated on the treatment of intramolecular polarization by subtracting the pairwise contribution, Pn ele i¼1 ci EAi B . The results of both these investigations are shown in Figure 5. The excellent correlation (R2 = 0.91 and 0.92, respectively) and the small mean absolute errors (2.7 and 2.2 kJ mol1, respectively) show that the model is robust and that it gives correct results for the right reason. Evidently, an accurate model for the intramolecular polarization is crucial for deriving the electrostatic properties of a molecule in terms of the properties of smaller fragments. The old exclusion rule (see Figure 1) gave a mean absolute error of 9.6 kJ mol1, that is, more than three times higher than the new rule, and a systematic energy shift (see Figure S1). If the old rule for static polarization was combined with the new (Thole) model for polarizability coupling, the mean absolute error only decreased to 8.9 kJ mol1, indicating that the major improvement with the new rule comes from the static polarization. 2.1.2. Electron-correlated Level

By using the same small basis set (cc-pVDZ), we were also able to compute reference interaction energies for the full system at the LMP2 level and thus evaluate the accuracy of the PMISP method if electron-correlation effects are taken into account. The results are also included in Table 1. At first sight, the mere observation that the reference many-body effect differs between the HF and LMP2 levels (by 5.5 kJ mol1 in average) seems to violate the basic assumption of PMISP, namely that Table 1. Many-body contribution obtained from supermolecular calculations with the cc-pVDZ basis set (referthe dispersion contribution is MB estimate (PMISP), and the difference between these numbers (error). Results are ence), the corresponding EAB pairwise additive, and thus given at the HF and LMP2 levels with a consistent level of theory used for computing multipoles and polarizaseems to add a substantial error bilities (the 6-31G* basis set was used for HF properties, but the difference is negligible). Finally, the error for MB to PMISP. However, the effect term at the MP2 level is also given. The mean absolute error a mixed model employing the HF-derived EAB (MAE) and maximum error over the set of nine systems are also reported. All energies are in kJ mol1. does not really come from the dispersion (by definition, disperHF level LMP2 level Mixed sion is additive at the MP2 level), Guest Reference PMISP Error Reference PMISP Error Error but from electron-correlation efBz 3.7 1.1 2.6 13.1 4.8 8.3 12.0 fects on other physical terms MeBz 11.4 14.9 3.5 16.7 10.7 6.0 1.8 such as polarization and repulEtBz 11.7 15.6 3.9 13.0 14.2 1.2 2.6 pClBz 0.5 0.5 0.0 4.3 7.8 3.6 3.8 sion. Some of those effects can mClBz 11.4 14.0 2.6 23.1 26.5 3.4 9.2 in fact be modeled by PMISP. Hx 8.4 9.0 0.6 4.2 2.6 1.7 4.7 Indeed, we see from Table 1 that MeHx 9.0 12.0 3.0 4.1 9.1 4.9 7.9 if we compute the fragment Pen 10.9 5.3 5.7 7.2 0.5 6.7 2.0 Hep 17.4 18.5 1.1 11.8 17.0 5.2 6.7 multipoles and polarizabilities inMAE 2.6 4.5 5.6 cluding electron correlation (at Max Err 5.7 8.3 12.0 the LMP2 level), the mean abso-

 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

ChemPhysChem 2014, 15, 3270 – 3281

3276

CHEMPHYSCHEM ARTICLES

www.chemphyschem.org

Figure 6. Correlation energy from the reference LMP2 calculation plotted against the correlation energy estimated by PMISP by using a mixed model (keeping HF-derived multipoles) or a consistent model (changing to LMP2derived multipoles). The line indicates perfect agreement.

Figure 5. Accuracy of the PMISP method at the HF level evaluated for the MB versus the correspondfull and truncated systems in Figure 2. Top: total EAB ing reference from a supermolecular calculation. Bottom: electrostatic part MB;ele versus the corresponding reference from multipoles computed for the EAB molecule A as a whole (multipoles for the full host were only computed for mClBz). The line indicates perfect agreement.

lute PMISP error increases, but only to 4.5 kJ mol1 (maximum 8.3 kJ mol1). In contrast, by using a mixed model (i.e. employing HF-derived multipoles throughout) gives a slightly larger mean absolute error (5.6 kJ mol1; maximum 12.0 kJ mol1). A further illustration of this effect is provided in Figure 6, in which the correlation energy from the reference calculation is plotted against the correlation energy estimated by PMISP. The use of a consistent model (i.e. changing the multipoles to those obtained at the LMP2 level) pushes many points closer to the line representing perfect agreement and the mean absolute error decreases from 5.5 to 3.5 kJ mol1. It should also be noted that the electrostatics in the reference MP2 calculation does not really correspond to MP2 calculations of both monomers,[53] so perfect agreement should not be expected; in addition, there are other approximations involved. The PMISP method builds on the assumption that the total dispersion contribution to the interaction energy can be writ 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

ten as a linear combination of fragment–guest dispersion energies, that is, that the long range many-body effects on the dispersion energy are negligible. As already mentioned, this assumption cannot be validated at the LMP2 (or MP2) level, because many-body effects on the dispersion energy appear only at a higher level of theory. Recently, it has been suggested that such effects are important for this type of system.[54] A systematic evaluation of the full systems at the coupled-cluster level is too computationally demanding. Therefore, we apply an approximate method, MBD,[55, 56] to investigate the importance of including the long range many-body dispersion. The MBD energy is computed by using the coupled Quantum Harmonic Oscillator model Hamiltonian. As we are interested mainly in the order of magnitude of the many-body effects, we employ a set of standard polarizabilities[57] instead of system-specific Hirshfeld polarizabilities recommended for accurate results. The detailed results are shown in Table S4. Interestingly, the many-body contribution is 1–2 kJ mol1 for all nine systems and always with a negative sign. Thus, compared with other errors in the method, it is negligible in both absolute and relative terms. 2.1.3. Basis-set Dependence Having established that the PMISP method accurately reproduces the reference interaction energies for all investigated systems at the HF and LMP2 levels with small basis sets, it remains to verify that the performance is not deteriorated if the size of the basis set is increased. To this end, we performed reference calculations of the truncated host (Trunc) systems at the LMP2 level by using two different basis sets, cc-pVDZ and cc-pVTZ, and studied how well the change in interaction energy was reproduced by the PMISP estimates. As shown in Table 2, the error was below 3 kJ mol1 for all except one ChemPhysChem 2014, 15, 3270 – 3281

3277

CHEMPHYSCHEM ARTICLES

www.chemphyschem.org

Table 2. Effect (kJ mol1) of changing the basis set from cc-pVDZ to ccpVTZ in the reference LMP2 calculation of the whole complex and in the corresponding PMISP estimate, respectively. The difference between PMISP and the reference is also given. Guest

Reference

PMISP

Difference

Bz MeBz EtBz pClBz mClBz Hx MeHx Pen Hep

28.8 32.6 31.6 35.4 39.2 32.6 25.2 32.2 34.5

27.1 35.2 38.3 33.1 40.4 35.6 25.9 34.8 35.3

1.7 2.6 6.7 2.4 1.2 3.0 0.7 2.6 0.8

2.2. Free Energy Prediction

system, despite the fact that the same multipoles were employed in the two sets of calculations. This indicates that the basis-set difference is mainly due to dispersion interactions and that these are well approximated by a pairwise additive model. We also investigated the magnitude of the total counterpoise correction (i.e. the PMISP sum of 56 counterpoise corrections for the individual fragments). For a sufficiently large basis set, or if extrapolating the result to the CBS limit, the correction [or equivalently, the basis set superposition error (BSSE)] is expected to go towards zero. From Table 3, we see that the

Table 3. Counterpoise correction (kJ mol1) for the whole host–guest complex (through MFCC summation) at the LCCSD(T0)/cc-pVTZ (LCC/TZ), LMP2/cc-pVTZ, and MP2/CBS levels, as well as the total correction obtained for the terms that are combined as in Equation (10). Guest

LCC/TZ

LMP2/TZ

MP2/CBS

Total

Bz MeBz EtBz pClBz mClBz Hx MeHx Pen Hep

21.4 15.5 22.0 17.1 20.2 23.7 13.2 20.7 16.3

23.4 16.8 24.2 19.3 22.0 25.5 14.2 22.3 18.2

2.3 2.1 2.5 2.6 2.4 2.8 2.4 1.9 2.4

0.3 0.7 0.4 0.4 0.6 1.0 1.4 0.3 0.5

BSSE is substantial ( 20 kJ mol1) if using the ccpVTZ basis set, but it almost cancels (within 2 kJ mol1) between the LMP2 and LCCSD(T0) levels. Moreover, at the extrapolated CBS limit, the effect of employing the counterpoise correction is only 2.4 kJ mol1 on average, verifying that the extrapolation procedure works. This is satisfying, because there has been much discussion in the literature regarding to what extent the

BSSE correction should be applied. Apparently, for the systems in this study, there is no significant difference between applying the correction or not, as long as the results are extrapolated to the CBS limit (the results before extrapolation are given in Table S3).

The final vacuum interaction energies calculated at the LCCSD(T0)/CBS level (shortened to LCC in the following) are given in Table 4. Besides the many-body effect computed directly from the PMI model, the total energy also contains a correction for the PMISP error, estimated at the LMP2/cc-pVDZ level (i.e. the error in the sixth column in Table 1). Apart from approximations in the PMI model itself, this correction includes nonclassical many-body effects, such as those observed in the exchange repulsion term.[58] The investigations above showed that these effects are reasonably independent of the basis set; thus they can be estimated at the LMP2/cc-pVDZ level without significantly increasing the overall computational cost of the calculation. Notably the many-body contribution (i.e. the sum of the three contributions in Table 4) differs significantly between the various guests, from 23 kJ mol1 for mClBz to +17 kJ mol1 for MeBz. Thus, it affects the relative energies and, in some cases, even the predicted ordering of the host–guest binding strength (e.g. reversing the order of Bz and Pen). The electrostatic part is seen to be the most important part of the manybody contribution, whereas the induction part is smaller and less varying among the ligands. This is mainly because the electrostatic contribution to the total interaction energy is larger and more varying than the induction contribution, and it supports the conclusion from earlier studies[25] that inclusion of polarization into force fields is not primarily a means to improve induction interactions but to improve the electrostatic description of molecules built up from smaller fragments. The strongly repulsive interaction energy is a vacuum effect (trivial repulsion between two negative molecules), and any comparison with aqueous phase experiments requires a treatment of solvation effects. In addition, thermal effects should be added to produce proper binding free energies. Table 5

Table 4. The total PMISP energies together with their various energy contributions (kJ mol1). The MFCC contriMFCC is given at three different levels of theory, in which LCC stands for LCCSD(T0) and TZ for cc-pVTZ. bution EAB MB;ele MB;ind and EAB , The many-body contribution has been divided into electrostatic and induction contributions, EAB as well as a correction term subtracting the PMISP error at the LMP2/cc-PVDZ level.

Guest

MFCC contribution LCC/TZ LMP2/TZ

MP2/CBS

Many-body contribution Ele Ind

Corr

Total PMISP LCC/CBS

Bz MeBz EtBz pClBz mClBz Hx MeHx Pen Hep

923.4 879.0 868.0 971.7 979.8 933.5 911.3 901.5 919.3

839.4 791.5 771.8 879.4 884.4 853.4 832.5 826.3 839.7

10.0 9.4 12.1 13.9 34.0 1.6 5.9 2.6 12.1

8.3 6.0 1.2 3.6 3.4 1.7 4.9 6.7 5.2

840.5 826.3 802.2 893.2 875.8 874.0 849.8 848.9 867.3

 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

909.2 860.8 850.6 953.6 965.3 917.1 898.1 886.1 903.6

5.2 1.3 2.1 6.0 7.4 4.2 3.1 3.1 4.9

ChemPhysChem 2014, 15, 3270 – 3281

3278

CHEMPHYSCHEM ARTICLES

www.chemphyschem.org

Table 5. Contributions (kJ mol1) to the binding free energy for each of the nine guest molecules. Vacuum interaction energies at the DFT-D (TPSS-D3/ def2-QZVP’) and LCC (LCCSD(T0)/CBS/PMISP) levels are given, together with the classical PMI interaction energy. These are followed by the difference in solvation free energy evaluated with COSMO-RS at the BP86/TZVP and LMP2/cc-pVDZ levels, respectively, and two other correction terms: thermal corrections (ZPE, entropy, and enthalpy corrections), and the strain energy of the guest at the geometry of the complex (evaluated with DFT). Finally, the resulting free energies at the DFT-D and LCC levels are given together with the experimental values. The mean absolute errors (MAE) and the mean absolute errors after subtracting the median errors (MAE tr) are also given.

Guest Bz MeBz EtBz pClBz mClBz Hx MeHx Pen Hep MAE MAE tr

Vacuum energy DFT-D LCC 886.1 856.9 843.0 937.6 936.2 897.7 886.2 878.5 899.0

840.5 826.3 802.2 893.2 875.8 874.0 849.8 848.9 867.3

PMI 962.6 934.5 931.4 1031.5 1022.6 958.7 976.4 915.6 983.1

[DGsolv] DFT[a] 948.3 937.7 927.2 1021.1 1015.3 957.2 973.2 926.7 982.5

Other corrections[a] [DGthermal] [DGstrain]

LMP2 945.6 930.8 929.3 1011.3 1005.2 959.6 963.5 927.2 976.1

47.9 50.5 50.5 50.7 50.7 48.5 49.7 47.1 57.3

1.2 1.9 3.2 1.3 1.1 1.9 1.2 5.8 3.5

Free energy DFT-D 13.1 28.4 30.5 31.5 27.3 9.1 36.1 4.7 22.7 7.0 6.7

[b]

LCC

Exp

55.9 52.0 73.4 66.0 77.6 35.2 62.8 25.3 48.0 31.3 12.4

15.6 24.5 26.2 28.1 21.9 23.5 31.8 15.6 27.6

[a] The values of the corrections are identical to those used for calculating the rigid free energies in Table 6 of Ref. [30], although the individual corrections are not listed in that table (only relaxed values are given). [b] This value differs by 4 kJ mol1 from the value in Table 6 of Ref. [30] owing to a slightly different geometry for the host–Bz complex that was consistently used throughout this study.

summarizes our account of the various contributions to the free energy, some taken from our earlier study.[30] The solvation energy was reevaluated at the LMP2 level to be more consistent with the LCC vacuum interaction energies. Our final estimates of the binding free energies are also listed, as well as the experimental values, which were unknown to us at the time of the calculations. It should be noted that Ref. [30] also contains an estimate of the LCC energies, but those are less reliable owing to the use of a less accurate model for intramolecular polarization (cf. Figure S1), an inconsistent definition of domains in the localization procedure for the LCC calculations (see Methods section), and the use of DFT solvation energies for the LCC free energy. The predicted LCC free energies are plotted against the experimentally measured free energies in Figure 7. A full agreement of the absolute free energies cannot be expected, because some terms are designed to be compatible with results at a certain level of theory. It is more problematic that even the relative energies are inaccurate, and actually worse than the DFT-D estimates. After subtraction of a constant, the mean absolute errors for LCC and DFT-D are 12.4 and 6.7 kJ mol1, respectively, that is, worse than for the trivial model predicting all binding affinities to be equal (MAE 4.2 kJ mol1). At first sight, a possible reason for this failure is that the solvation free-energy contribution is evaluated at the LMP2 level and not at the LCC level, thus introducing an inconsistency in the methodology. It is clear from Table 5 that the vacuum interaction energy is to a large extent cancelled by the solvation contribution (as is also observed for protein–ligand systems), and therefore any inconsistency between these terms might cause large errors. However, based on previous experience, the difference in solvation between the LCC and LMP2 levels is expected to be smaller in magnitude than the difference between LMP2 and DFT and, as shown in Table 5, the latter difference is only 5.6 kJ mol1 on average (maximum 10 kJ mol1), that is,  2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

Figure 7. Predicted DFT-D and LCC binding free energies versus experimental results. The line indicates perfect agreement.

smaller than the discrepancy with experimental values. In fact, in terms of agreement with experimental values, it does not matter if the LCC energies are combined with DFT or LMP2 solvation energies. Thus, our results suggest that the use of inconsistent levels of theory is not the reason for the large discrepancy with experimental values. An illustrative example of the importance of the cancellation between electrostatics and solvation is provided by the mClBz and Hx guests. As can be seen in Table 5, they have almost identical LCC vacuum energies and almost identical experimental binding free energies. Thus, the big difference in DGsolv (45 kJ mol1) appears suspicious at first sight. However, the PMI contribution to the vacuum energy differs by 64 kJ mol1, sugChemPhysChem 2014, 15, 3270 – 3281

3279

CHEMPHYSCHEM ARTICLES gesting that the difference in DGsolv is actually reasonable. Together with the method insensitivity of DGsolv, this example seems to indicate that there are additional problems. One such problem could be the use of single optimized geometries; perhaps a proper sampling method such as free-energy perturbation is required to obtain reliable results. Unfortunately, our previous study showed that this is not done straightforwardly even at the DFT-D level,[30] and it becomes even more demanding at the coupled-cluster level. An advantage of having reliable interaction energies computed at a high level of theory is that it allows us to benchmark other, computationally cheaper methods. From a drug development perspective, the most interesting quantities in such comparisons are relative interaction energies of the various ligands. Therefore, the relative energies calculated at various levels of theory are displayed in Figure 8. As expected, the

Figure 8. Relative binding energies of the nine guests calculated in vacuum at four different levels of theory (DFT = [TPSS/def2-QZVP’], DFT-D = [TPSSD3/def2-QZVP’], MP2 = [MP2/CBS/PMISP], and LCC = [LCCSD(T0)/CBS/ PMISP]), free energies at the LCC level (including contributions from solvation, strain, and thermal effects), and experimental free energies. For each method, the average energy over the nine ligands has been subtracted.

DFT results without dispersion correction are not reliable for this type of system; the absolute interaction energies differ from the LCC energies by 160 kJ mol1 on average and six pairs of guests are ordered incorrectly. It is more surprising that even dispersion-corrected DFT appears to systematically underestimate the effect of dispersion, as shown by the tendency of most DFT-D points (all except Bz and MeBz) in Figure 8 to lie between the DFT and LCC results. In absolute terms, the average error of DFT-D is 38 kJ mol1, and two pairs of guests (Hep-Hx and Bz-Pen) are still ordered incorrectly. MP2, on the other hand, overestimates the dispersion energy by 16 kJ mol1 on average. This is a well known problem with MP2. However, the error is almost constant for the various ligands and hardly affects the relative energies. The change in relative energies on addition of solvation and thermal corrections to the LCC result is also shown in Figure 8, as well as the comparison to experimental values. Clearly, it is not sufficient to consider vacuum interaction energies when predicting the ranking of guests binding to the octa-acid host,  2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

www.chemphyschem.org and as discussed above, the solvation contribution used in this study does not offer any significant improvement. It can also be noted that the range of the experimental binding free energies is much narrower than for any of the theoretical estimates. This is a general tendency in free-energy calculations.

3. Conclusions We have developed a method to calculate interaction energies for large systems at the LCCSD(T0) level of theory extrapolated to the complete basis set limit. To accomplish this, we employed a combination of molecular fractionation and a multipole-polarizability expansion, denoted PMISP. Through a series of tests, we have established that the errors in the PMISP method itself and in the local treatment of electron-correlation effects are controllable and below 10 kJ mol1 for this test case. Compared with previous applications of PMISP, we were obliged to improve the treatment of intramolecular polarization to reach the desired level of accuracy. It is especially satisfying that the results were insensitive to the few parameters introduced and thus the method can be regarded as fully automated, once provided with the fractionation scheme and the usual choices in quantum chemistry, such as the basis set extrapolation scheme. The method was tested for the octa-acid host–guest system, enabling the whole complex (198–207 atoms) to be treated consistently. The same method can be applied to protein– ligand systems, provided that the ligand is fairly small (or possible to fractionate) so that the size of the interacting subsystems is small enough to allow a high level of theory and basis set. To reduce the computational cost for larger systems, the method can be combined with a purely classical treatment of long range interactions.[23] The high-level vacuum interaction energies were combined with solvation effects and thermal contributions to produce predictions of binding free energies submitted to the SAMPL4 challenge. Unfortunately, only weak correlation with experimental values was obtained (R2 = 0.2). We attribute the discrepancy to the approximation to use single optimized geometries instead of a proper sampling of all relevant geometries. An exhaustive study of the sampling issue, albeit at the semiempirical level, is currently underway.

Acknowledgements This investigation has been supported by grants from the Knut and Alice Wallenberg foundation (project 2013.0022) and the Swedish Research Council (agreement C0020401 and project 2010-5025). The computations were performed on computer resources provided by the Swedish National Infrastructure for Computing (SNIC) at Lunarc at Lund University and HPC2N at Ume University. The collaboration between the Universities of Lund and Gçttingen has been carried out within the framework of the International Research Training Group 1422 “Metal Sites in Biomolecules — Structures, Regulation, Mechanisms” and M. A. is ChemPhysChem 2014, 15, 3270 – 3281

3280

CHEMPHYSCHEM ARTICLES supported through a Ph.D. scholarship in this international research training group. Keywords: cluster compounds · ligand binding · molecular fractionation · polarizable force fields · quantum chemistry [1] P. Jurecˇka, J. Sˇponer, J. Cˇerny´, P. Hobza, Phys. Chem. Chem. Phys. 2006, 8, 1985. [2] S. Grimme, Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1, 211 – 228. [3] C. Hampel, H. J. Werner, J. Chem. Phys. 1996, 104, 6286. [4] G. E. Scuseria, P. Y. Ayala, J. Chem. Phys. 1999, 111, 8330. [5] M. Schtz, F. Manby, Phys. Chem. Chem. Phys. 2003, 5, 3349 – 3358. [6] J. E. Subotnik, A. Sodt, M. Head-Gordon, J. Chem. Phys. 2006, 125, 074116. [7] J. Yang, G. K.-L. Chan, F. R. Manby, M. Schtz, H.-J. Werner, J. Chem. Phys. 2013, 136, 144105. [8] C. Riplinger, F. Neese, J. Chem. Phys. 2012, 138, 034106. [9] J. Kaminsky, R. A. Mata, H.-J. Werner, F. Jensen, Mol. Phys. 2008, 106, 1899. [10] J. M. Dieterich, J. C. A. Oliveira, R. A. Mata, J. Chem. Theory Comput. 2012, 8, 3053 – 3060. [11] M. Gordon, D. Fedorov, S. Pruitt, L. Slipchenko, Chem. Rev. 2012, 112, 632 – 672. [12] P. Sçderhjelm, S. Genheden, U. Ryde in Protein – Ligand Interactions (Ed.: H. Gohlke), Wiley-VCH, Weinheim, 2012. [13] E. E. Dahlke, D. G. Truhlar, J. Chem. Theory Comput. 2007, 3, 46. [14] D. W. Zhang, J. Z. H. Zhang, J. Chem. Phys. 2003, 119, 3599. [15] K. Kitaura, E. Ikeo, T. Asada, T. Nakano, M. Uebayasi, Chem. Phys. Lett. 1999, 313, 701 – 706. [16] S. Hirata, M. Valiev, M. Dupuis, S. Xantheas, S. Sugiki, Mol. Phys. 2005, 103, 2255. [17] W. Li, S. Li, Y. Jiang, J. Phys. Chem. A 2007, 111, 2193. [18] N. Jiang, J. Ma, Y. Jiang, J. Chem. Phys. 2006, 124, 114112. [19] N. Mayhall, K. Raghavachari, J. Chem. Theory Comput. 2011, 7, 1336 – 1343. [20] G. S. Tschumper, Chem. Phys. Lett. 2006, 427, 185. [21] P. Sçderhjelm, U. Ryde, J. Phys. Chem. A 2009, 113, 617. [22] C. Steinmann, D. Fedorov, J. Jensen, J. Phys. Chem. A 2010, 114, 8705 – 8712. [23] P. Sçderhjelm, F. Aquilante, U. Ryde, J. Phys. Chem. B 2009, 113, 11085. [24] P. Sçderhjelm, J. Kongsted, U. Ryde, J. Chem. Theory Comput. 2010, 6, 1726. [25] P. Sçderhjelm, Theor. Chem. Acc. 2012, 131, 1159. [26] P. Ren, J. W. Ponder, J. Comput. Chem. 2002, 23, 1497. [27] O. Engkvist, P.-O. strand, G. Karlstrçm, J. Phys. Chem. 1996, 100, 6950. [28] G. Tiraboschi, M. Fourni-Zaluski, B.-P. Roques, N. Gresh, J. Comput. Chem. 2001, 22, 1038 – 1047.

 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

www.chemphyschem.org [29] A. Holt, G. Karlstrçm, J. Comput. Chem. 2008, 29, 1084. [30] P. Mikulskis, D. Cioloboc, M. Andrejic´, S. Khare, J. Brorsson, S. Genheden, R. A. Mata, P. Sçderhjelm, U. Ryde, J. Comput.-Aided Mol. Des. 2014, 28, 375 – 400. [31] H. S. Muddana, A. T. Fenley, D. L. Mobley, M. K. Gilson, J. Comput.-Aided Mol. Des. 2014, 28, 305 – 317. [32] N. R. Kestner, J. Chem. Phys. 1968, 48, 252. [33] L. Gagliardi, R. Lindh, G. Karlstrçm, J. Chem. Phys. 2004, 121, 4494. [34] P. Cieplak, J. Caldwell, P. Kollman, J. Comput. Chem. 2001, 22, 1048. [35] J. Ponder, C. Wu, P. Ren, V. Pande, J. Chodera, M. Schnieders, I. Haque, D. Mobley, D. Lambrecht, R. DiStasio, Jr., J. Phys. Chem. B 2010, 114, 2549 – 2564. [36] B. T. Thole, Chem. Phys. 1981, 59, 341. [37] H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schtz, et al., MOLPRO, version 2012.1, a package of ab initio programs, 2012. [38] T. H. Dunning, Jr., J. Chem. Phys. 1989, 90, 1007 – 1023. [39] R. A. Kendall, T. H. Dunning, R. J. Harrison, J. Chem. Phys. 1992, 96, 6796. [40] F. Weigend, A. Kçhn, C. Httig, J. Chem. Phys. 2002, 116, 3175. [41] J. Pipek, P. G. Mezey, J. Chem. Phys. 1989, 90, 4916. [42] R. A. Mata, H.-J. Werner, Mol. Phys. 2007, 105, 2753 – 2761. [43] R. A. Mata, H.-J. Werner, S. Thiel, W. Thiel, J. Chem. Phys. 2008, 128, 025104. [44] N. Lttschwager, T. N. Wassermann, R. A. Mata, M. Suhm, Angew. Chem. Int. Ed. 2013, 52, 463 – 466; Angew. Chem. 2013, 125, 482 – 485. [45] MOLCAS 7, University of Lund, Sweden, 2007. [46] A. J. Sadlej, J. Chem. Phys. 1981, 75, 320 – 331. [47] R. Ahlrichs, M. Br, M. Hser, H. Horn, C. Kçlmel, Chem. Phys. Lett. 1989, 162, 165. [48] O. Treutler, R. Ahlrichs, Chem. Phys. 1995, 102, 346. [49] A. Klamt, J. Phys. Chem. 1995, 99, 2224. [50] F. Eckert, A. Klamt, AIChE Journal 2002, 48, 369. [51] F. Eckert, A. Klamt, COSMOtherm, version c30, release 1301. COSMOlogic GmbH & Co KG, Leverkusen, 2010. [52] S. Grimme, J. Antony, S. Ehrlich, H. Krieg, J. Chem. Phys. 2010, 132, 154104. [53] G. ChaN´asin´ski, M. M. Szczs´niak, Mol. Phys. 1988, 63, 205 – 224. [54] A. Tkatchenko, D. S. Alf, K. S. Kim, J. Chem. Theory Comput. 2012, 8, 4317. [55] A. Tkatchenko, M. Scheffler, Phys. Rev. Lett. 2009, 102, 073005. [56] A. Tkatchenko, R. A. DiStasio Jr, R. Car, M. Scheffler, Phys. Rev. Lett. 2012, 108, 236402. [57] O. A. von Lilienfeld, A. Tkatchenko, J. Chem. Phys. 2010, 132, 234109. [58] V. F. Lotrich, K. Szalewicz, J. Chem. Phys. 1997, 106, 9668.

Received: May 29, 2014 Published online on September 26, 2014

ChemPhysChem 2014, 15, 3270 – 3281

3281

Coupled-cluster interaction energies for 200-atom host-guest systems.

We have developed a method to calculate interaction energies of large systems (such as host-guest or even protein-ligand systems) at the local coupled...
1MB Sizes 0 Downloads 6 Views