Mixed ExplicitâImplicit Solvation Approach for Modeling of Alkane

Subscriber access provided by Kaohsiung Medical University

Article

A Mixed Explicit-Implicit Solvation Approach for Modeling of Alkane Complexation in Water-Soluble Self-Assembled Capsules Henrik Daver, Andrés G. Algarra, Julius Rebek, Jr., Jeremy N. Harvey, and Fahmi Himo J. Am. Chem. Soc., Just Accepted Manuscript • DOI: 10.1021/jacs.8b06984 • Publication Date (Web): 06 Sep 2018 Downloaded from http://pubs.acs.org on September 6, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of the American Chemical Society

A Mixed Explicit-Implicit Solvation Approach for Modeling of Alkane Complexation in Water-Soluble Self-Assembled Capsules

Henrik Davera, Andrés G. Algarrab, Julius Rebek, Jr.c,d, Jeremy N. Harveye*, Fahmi Himoa*

a. Department of Organic Chemistry, Arrhenius Laboratory, Stockholm University, SE-106 91 Stockholm, Sweden b. Departamento de Ciencia de los Materiales e Ingeniería Metalúrgica y Química Inorgánica, Facultad de Ciencias, Instituto de Biomoléculas, Universidad de Cádiz, Puerto Real, Cádiz, 11510, Spain c. The Skaggs Institute for Chemical Biology and Department of Chemistry, The Scripps Research Institute, 10550 North Torrey Pines Road, La Jolla, California, 92037 d. Center for Supramolecular Chemistry and Catalysis, Shanghai University, Shanghai 200444, China e. Department of Chemistry, KU Leuven, Celestijnenlaan 200F, B-3001 Heverlee, Belgium

* Corresponding authors: [email protected], [email protected] ORCID: Andrés G. Algarra: 0000-0002-5062-2858 Julius Rebek, Jr.: 0000-0002-2768-0945 Jeremy N. Harvey: 0000-0002-1728-1596 Fahmi Himo: 0000-0002-1012-5611

1 ACS Paragon Plus Environment

Journal of the American Chemical Society 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract The host-guest binding properties of a water-soluble resorcinarene-based cavitand are examined using density functional theory methodology. Experimentally, the cavitand has been observed to self-assemble in aqueous solution into both 1:1 and 2:1 host:guest complexes with hydrophobic guests such as n-alkanes. For n-decane, equilibrium was observed between the 1:1 and 2:1 complexes, while 1:1 complexes are formed with shorter nalkanes and 2:1 complexes with longer ones. These findings are used to assess the standard quantum chemical methodology. It is first shown that a rather advanced computational protocol (B3LYP-D3(BJ)/6-311+G(2d,2p) with COSMO-RS and quasi-rigid rotor harmonic oscillator) gives very large errors. Systematic examination of the various elements of the methodology shows that the error stems from the implicit solvation model. A mixed explicitimplicit solvation protocol is developed that involves a parameterization of the hydration free energy of water such that water cluster formation in water is predicted to be thermoneutral. This new approach is demonstrated to lead to a major improvement in the calculated binding free energies of n-alkanes, reproducing very well the 1:1 vs. 2:1 host:guest binding trends.


Page 2 of 29

Page 3 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


I. Introduction Reversible encapsulation is a process that temporarily confines guest molecules in container host structures. The guest is isolated from other solutes and the solvent, as the container walls act as more-or-less rigid solvation shells. The containers are usually assemblies that emerge through hydrogen bonding in organic media,[1-9] while metal/ligand interactions,[10,11] ionic forces[12] or even purely hydrophobic effects can stabilize capsules in water.[13] A number of contortions accompany the confinement of flexible guest molecules in the hosts including coiling, bending, folding and other distortions that effect reactivity.[14,15] In the framework of homogeneous chemistry, the containers allow the study of species or events that cannot usually be observed in solution. Potential applications include selective binding, recognition of unstable intermediates and encapsulation for catalysis purposes, to mention just a few.[16] The hosts can also be considered as mimics of biomolecules, which usually have binding pockets that selectively bind guests with high affinities.[17] Several self-assembling complexes have been synthesized, constructed from non-covalent linkages, such as hydrogen bonds, allowing for flexible binding and release of molecules to the solution.[16,18] Of particular interest are molecular containers that can self-assemble in aqueous medium. Water is ubiquitous in biology, and, thus, in order to properly mimic biological events the hosts must be active in this medium. The aqueous environment implies competition for all polar interactions, including hydrogen bonds, making self-assembly more difficult. On the other hand, the hydrophobic effect might favor the process of self-assembly.[17] A number of hosts have been designed for self-assembly in water, and supramolecular hosts can make guest binding and reactivity possible that are usually not viable in aqueous medium.[17-21] One specific example of a water-soluble self-assembling molecule is cavitand 1 (Figure 1), developed by Rebek and co-workers,[22] building upon the design of de Mendoza.[23] This host has shown interesting binding trends for a variety of alkane guests in water.[22,24,25] Complexes of 1 host n-alkanes up to n-C18H38 (here referred to as C18 for brevity).[25] The stoichiometries of the complexes depend on the alkane length. Alkanes up to n-nonane (C9) result in 1:1 host:guest complexes,[24] whereas C11 and longer n-alkanes form 2:1 complexes in which the alkane is incarcerated inside the 12 capsule in an extended conformation.[22,24] For C10, equilibrium is observed between the 1:1 and 2:1 complexes, 3 ACS Paragon Plus Environment


Page 4 of 29

as expressed in Equilibrium (1-C10).[24]

2 × C10@1 (aq) ⇌ C10@12 (aq) + C10(l)

Equilibrium (1-C10)

Water has been suggested to make specific interactions with both cavitand 1 and capsule 12.[22,24,25] For cavitand 1, with its vase-like structure, water molecules are suggested to bridge between the benzimidazolone groups of the vase’s rim, forming hydrogen bonds with the urea groups.[22,24] Studies on water-soluble cyclodextrin hosts have shown that release of host-bound water molecules that are energetically ‘frustrated’, i.e. that have insufficient numbers of possible hydrogen bonds, is a driving force in the formation of inclusion complexes of other guests than water.[26,27] A similar effect has been proposed for 1.[22] It has also been suggested that water molecules can play a role in the widening of the capsule by binding along the seam of the dimer in order to accommodate folded forms of n-alkanes longer than C14.[25]

Figure 1. Cavitand 1 and capsule 12. In the picture of the capsule, the R groups are replaced with methyl groups for clarity.

Quantum chemical methodology has recently become an important tool in the study of supramolecular host-guest complexes.[28-52] The various components of the methodology have gradually improved and more accurate results are possible today. Advances in density functional theory are of course of great importance here, but also the development of various 4 ACS Paragon Plus Environment

Page 5 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


corrections, such as dispersion and free energy corrections.[33,53] A significant number of theoretical studies of supramolecular complexes in water have applied implicit solvation models.[28,30,31,33,37,38,47,50,51,53-55] However, none of these studies treated dimerization processes. In a molecular dynamics study of cyclodextrin dimerization, it was found that explicit solvent models performed better than generalized Born models.[56] Furthermore, in the SAMPL5 challenge of calculating host-guest binding free energies, in which a wide variety of theoretical methods were applied, it was found that methods involving explicit water molecules in general performed better than implicit solvent models.[57] Aqueous solvation clearly creates a challenge for standard quantum chemical protocols, and a combination of implicit and explicit solvation has been suggested to provide a general improvement for calculations of pKa, solvation free energies and reaction mechanisms in case of strong interactions between the solute and the solvent.[58-63] In the present work, we study the host-guest binding properties of 1 and 12 using DFT methodology. We will use Equilibrium (1-C10) as a benchmark, because it is well-defined and quantifiable in terms of free energy. Experimentally, the C10@1 and C10@12 complexes were observed in a roughly 1:1 ratio.[24] Hence, out of three molecules of 1, two are engaged in the C10@12 complex and one in the C10@1 complex. The concentration of cavitand 1 in solution was 0.5 mM, which means that the concentrations of both C10@1 and C10@12 were 0.167 mM. Assuming alkane C10 to be present in excess as a separate phase and its aqueous concentration to be approximately constant, the equilibrium constant K' for Equilibrium (1C10) becomes:

=

C10@

0.167 × 10 mol ∙ l = = = 6 × 10 l ∙ mol 0.167 × 10 mol ∙ l C10@

C10

which can be converted to a free energy difference of -RT ln(K') = −5 kcal/mol (including the standard free energy of C10 in its pure phase). This value will be used as a reference to assess the various computational protocols. There are a couple of other pieces of information known about this reaction, and these will be mentioned in the appropriate places of the discussion below. We show that the first attempt to model Equilibrium (1-C10), using a rather advanced methodology, yields very large errors. Systematic examination of the various elements of the



methodology is then presented, showing that the main deviation stems from the solvation model. An amended computational protocol is then proposed that significantly reduces the error.

II. Reference computational protocol The total Gibbs energy of the system is calculated according to the following equation: Gtot = Egas + δGfree + δGsol where Egas is the gas-phase energy of the system, δGfree is the free energy correction to the gas-phase energy, and δGsol is the solvation free energy. First, the calculations were performed using the following computational scheme, which will be referred to as the Reference protocol. This computational protocol has been used very recently in the study of a cycloaddition reaction inside a very similar self-assembled capsule, but in an organic solvent, and yielded very good results.[45] The calculations were performed using the B3LYP-D3(BJ) functional,[53,64-66] i.e. the B3LYP functional augmented with the DFT-D3 dispersion method, utilizing the Becke-Johnson damping function, as implemented in Gaussian 09.[67] Geometry optimizations were performed with the 631G(d,p) basis set, and on the basis of these optimized geometries, single-point calculations were conducted with the 6-311+G(2d,2p) basis set. In Gaussian 09, only the two-body part of the D3(BJ) dispersion correction is included in the DFT-D3 method. The three-body terms were therefore added manually to the single-point energies calculated with the larger basis set. The sum of these contributions constitutes Egas. Next, frequency calculations were carried out at the same level of theory as the geometry optimizations to characterize the nature of the optimized structures and to calculate the free energy correction, δGfree. The quasi-rigid-rotor-harmonic-oscillator (qRRHO) approach [33] was applied, as this method corrects for the overestimation of the entropic contributions from low-frequency vibrational modes by the standard harmonic-oscillator approach. An additional benefit of this approach is that it minimizes the impact of errors in the magnitude of lowfrequency modes due to numerical error in the Hessian. For both reasons, this approach can be important for large systems with many soft modes and it has been used quite frequently in recent years, yielding generally good results.[68-73] The free energy calculations were done using rotational symmetry numbers σr = 1. Determination of the approximate molecular point 6 ACS Paragon Plus Environment

Page 6 of 29

Page 7 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


groups and symmetry numbers for all molecules was carried out in GaussView,[74] and the term RT ln(σr) was added to the calculated free energy corrections. For a couple of structures (the empty cavitand 1 and the water dimer, see Supporting Information), small imaginary frequencies were obtained that could not be eliminated by further optimization of the geometries. In these cases, to avoid inaccuracies that would arise from comparing entropic values calculated with different numbers of vibrational modes, these imaginary frequencies were treated as real in the qRRHO calculations. Finally, solvation free energy, δGsol, was calculated with the COSMOtherm software,[75,76] at 298.15 K and 1 atm pressure, at the BP86/TZVP level of theory.[77,78] This solvation method has been used in several similar applications in recent years.[33,36,51] To correct for the change in standard state in going from the ideal gas (1 atm standard state) to a water solution (1 mol/l standard state for all compounds except water, for which 55.3 mol/l is used), a value of RT ln(24.5 l/mol × 1 mol/l) = +1.9 kcal/mol is added for all solutes, and RT ln (24.5 l/mol × 55.4 mol/l) = +4.3 kcal/mol is added for water. The solvation free energies of the n-alkanes are calculated in their own phases since they are poorly soluble in water, and the corresponding liquid phase concentrations are used in the standard state corrections. An important note here concerns the conformational search. The size of the system under study requires a thorough examination of different conformations in order to identify the energetically most stable one. We have therefore explicitly considered different hydrogen bonding patterns, guest conformations and host-guest binding modes for each complex studied, and the structures reported here are the ones with the lowest free energy. A model of 1 is used in which the alkylpyridinium substituents and the Cl− counterions are replaced by methyl groups (the effect of this replacement will be examined below). In the most stable conformer of C10@1, the alkane adopts six gauche and one anti configurations for the internal C-C single bonds (Figure 2). The coiled part of the alkane is inside the cavitand while the four linearly arranged carbons are found in the solvent-exposed end. In free C10, each gauche interaction is calculated to be associated with an energy penalty of ca 0.7 kcal/mol. This penalty is, however, compensated by stronger interactions between the cavitand and the six-gauche conformation of C10 as well as a smaller solvent-exposed area. This isomer is calculated to be 1.2 kcal/mol more stable than the one with an all-anti C10. In the capsule complex C10@12, the most stable structure of C10 was found to be the all-anti



conformation, as shown in Figure 2. Using the lowest-energy structures of C10@1 and C10@12 and the computational protocol described above, the free energy of Equilibrium (1-C10) is calculated to be −46.8 kcal/mol, which is an overestimation by more than 40 kcal/mol compared to the experimental value of −5 kcal/mol. This shows that the adopted methodology, which is normally considered as relatively accurate and has been employed successfully in previous studies, fails severely in this case. In the following sections, we will examine the possible sources of this failure by considering a number of alterations of the various elements of the computational protocol.

Figure 2. B3LYP-D3(BJ)/6-31G(d,p)-optimized geometries of the most stable conformers of C10@1 and C10@12.

III. Cavitand model The alkylpyridinium substituents are crucial for the solubility in water, but these groups are not expected to strongly affect the binding properties of the cavitand. In the reference protocol presented above, these alkylpyridinium substituents, together with the Cl− counterions, were therefore replaced by methyl groups. Here, we evaluate the effect of this truncation on the free energy of Equilibrium (1-C10).


Page 8 of 29

Page 9 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1. Energies (kcal/mol) for Equilibrium (1-C10) calculated with different computational protocols.a ∆Egas

Entry Method 1

Experimental value

2

b

∆δGfree

∆δGsol

∆Gtot -5

Reference protocol

-115.9

+7.1

+62.0

-46.8

Cavitand model 3

[Full model]4+

+118.1

-99.7

+25.5

4

[Full model]4+·4Cl–

-115.5

+58.4

-50.0

Functional 5

TPSS-D3(BJ)c

6

M06-2X

7

c

c

B3LYP (no dispersion)

-117.1

+5.7

+63.3

-48.1

-95.7

+5.5

+60.1

-30.1

-105.9

+19.5

+66.2

-20.2

Dispersion model 8

XDMd

-115.8

-46.7

9

DFT-D2

-102.3

-33.2

Basis set 10

6-311+G(2d,2p) + Counterpoise

-114.7

-45.6

11

def2-TZVP

-116.1

-47.0

12

def2-TZVP + gCP

-115.1

-46.0

Free energy protocol 13

RRHO

+8.2

-45.7

14

qRRHO (only vibrational)

+3.0

-50.9

Implicit solvation model

15

SMD

+61.1

-47.7

16

C-PCM

+44.6

-64.2

17

SMD, joined cavities

+63.4

-45.4

18

C-PCM, joined cavities

+38.9

-69.9

[a]

Unless stated otherwise, all calculations are performed on basis of the geometries optimized at the B3LYP-D3(BJ)/6-31G(d,p) level of theory. [b] Truncated cavitand model with B3LYP-D3(BJ)/6-311+G(2d,2p), COSMO-RS, and qRRHO [c] The geometries are re-optimized with the respective functional and the 6-31G(d,p) basis set. [d] Based on the electron density calculated at the B3LYP/6-311+G(2d,2p) level of theory.

Using the full alkylpyridinium substituents means that cavitand 1 has 4+ charge, and bringing together two of these units in the dimer leads consequently to a large repulsion in the



gas phase. Indeed, as seen in Table 1 (entries 2 and 3), the gas-phase energy change (∆Egas) associated with Equilibrium (1-C10) for the full model is more than 230 kcal/mol higher than for the truncated model. The solvation correction (∆δGsol) of course acts in the opposite direction, but only by 160 kcal/mol. In total, thus, the model with alkylpyridinium substituents results in an energy of +25.5 kcal/mol, to be compared to −46.8 kcal/mol obtained with the truncated model. Adding eight chloride counterions to the model, such that the full model is neutral, gives energies that are very close to the truncated model (−50.0 vs. −46.8 kcal/mol, see Table 1 entries 2 and 4). This shows that the truncated model provides a good representation of the full system and can be used for the evaluation of the methodology. In water, the halide ions are not expected to be tightly bound to the alkylpyridinium substituents, so in principle the 4+ and 8+ models should yield accurate results, but we presume that continuum solvation of the multiple charged groups leads to the error obtained when including only the alkylpyridinium groups in the model (and not the chloride ions).

IV. Gas-phase energy In order to evaluate the sensitivity of the computational results to the choice of DFT functional, we have in addition to the B3LYP-D3(BJ) method also used the TPSSD3(BJ)[53,66,80] and M06-2X[81] functionals. These are popular methods and have been used previously in supramolecular chemistry applications.[30,33,36-38,40,42-44,47,54,82] The geometries were then re-optimized with these functionals, and corresponding thermal and solvation corrections were re-calculated. The results are given in Table 1 entries 5-7. While TPSS-D3(BJ) is found to give very similar results compared to B3LYP-D3(BJ), −48.1 vs. −46.8 kcal/mol, the difference is quite large using M06-2X. With this method, the overall free energy of Equilibrium (1-C10) is calculated to be −30.1 kcal/mol, which is 17 kcal/mol smaller in magnitude than with the reference protocol. As seen in Table 1, the difference stems mainly from the gas-phase energies. Interestingly, the M06-2X value is rather insensitive to whether the D3 correction is added or not. The ∆Egas value changes only by 1.7 kcal/mol upon inclusion of this correction (from −95.7 to −97.4 kcal/mol). On the other hand, The B3LYP functional without the dispersion correction gives a large difference compared to B3LYP-D3(BJ), −20.2 vs −46.8 kcal/mol (Table 1, entries 2 and 7).


Page 10 of 29

Page 11 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Notably, both the gas phase energies and the free energy corrections differ in this case. The geometry of the cavitand optimized with the B3LYP functional changes significantly with and without the dispersion correction. The diameter of the cavitand-decane complex C10@1, measured as the distance between the oxygen atoms of diametrically-opposite carbonyl groups at the rim, is ca 11 Å when the geometry is optimized with the uncorrected B3LYP functional and ca 8 Å when optimized with B3LYP-D3(BJ). As shown in Figure 3, the latter is much closer to the distances obtained when four water molecules are included in the model forming hydrogen bonds between the cavitand units (ca 7 Å).

Figure 3. Top views of the geometry of C10@1, optimized with the functionals a) B3LYP, b) B3LYP-D3(BJ), and c) B3LYP-D3(BJ) with four water molecules present. The alkane is made transparent for clarity.

As seen, uncorrected B3LYP gives a similar overall free energy for Equilibrium (1-C10) as M06-2X, which is closer to the experimental value than the reference B3LYP-D3(BJ) result. This is somewhat unexpected since it is generally accepted that B3LYP-D3(BJ) gives better results than uncorrected B3LYP. To further investigate this issue, one can consider the following equilibrium as an auxiliary condition on the modeling. Namely, according to the experiments Equilibrium (2-C10) should be exergonic, since complex C10@1 can be observed.[24]

1 (aq) + C10 (l) ⇌ C10@1 (aq)


Equilibrium (2-C10)


Note that in solution, a 1:10 mixture of C4v-symmetric 1 and the so-called ‘velcrand’ dimer (1-C2v)2 is observed,[22] so Equilibrium (2-C10) should include both these species on the left hand side. However, modeling of (1-C2v)2 proved to be quite challenging, since the positively charged alkylpyridinium substituents and the chloride counterions were found to interact with the 1-C2v units at the monomer-monomer interface (see Supporting Information). Interestingly, when n-decane is added to the mixture no signal from (1-C2v)2 or 1 was observed,[24] which means that C10@1 is more stable than both complexes, in particular 1, as formulated in Equilibrium (2-C10). Uncorrected B3LYP predicts Equilibrium (2-C10) to be endergonic by as much as +23.3 kcal/mol, and therefore C10 is calculated to be unbound to cavitand 1. B3LYP-D3(BJ), TPSS-D3(BJ) and M06-2X also predict endergonicity (+3.0, +5.1 and +1.7 kcal/mol, respectively), but these values are significantly lower than uncorrected B3LYP, and much closer to the experimental outcome. It can therefore be concluded that the better agreement of B3LYP for Equilibrium (1-C10) as compared to B3LYP-D3(BJ) must be fortuitous. To further assess the sensitivity of the results to the specific choice of dispersion correction, we have also used the wave function-based XDM protocol[83] and the DFT-D2 correction[84], in both cases in conjunction with the B3LYP functional (Table 1, entries 8 and 9). The calculations were performed on the basis of the geometries optimized with B3LYP-D3(BJ). While the XDM and D3(BJ) corrections give very similar results (−46.7 vs. −46.8 kcal/mol), the D2 correction results in a much lower exergonicity (−33.2 kcal/mol). The latter value is closer to the experimental result, but still far from the measured −5 kcal/mol. In addition, with XDM Equilibrium (2-C10) is calculated to have a positive free energy of +5.2 kcal/mol, while B3LYP-D2 is the only method considered here that predicts this equilibrium to be exergonic, by −7.0 kcal/mol. According to these results, it seems that the B3LYP-D2 method gives somewhat better agreement with the experimental findings than B3LYP-D3(BJ) and B3LYP-XDM. The latter ones are, however, generally found to be superior to B3LYP-D2,[85] which might indicate that the D2 benefits from a fortuitous cancellation of errors. Further tests with B3LYP-D2 will be presented in the last part of the present paper. As a final element of the evaluation of the gas-phase energy, we also investigated the sensitivity of the results to the basis set superposition error (BSSE). Since eight hydrogen bonds are formed between the cavitands in the 12 dimer, the energies might be prone to


Page 12 of 29

Page 13 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


BSSE, where the basis functions on one monomer extend to the other. This would result in an artificial stabilization of the supramolecular complex with respect to the monomers, which would favor the right hand side of Equilibrium (1-C10). To examine this effect, counterpoise calculations [86,87] were performed with the 6-311+G(2d,2p) basis set. As seen in Table 1 (entries 2 and 10), the difference is negligible, about 1 kcal/mol. In addition, calculations were also performed with the def2-TZVP basis set[88], and the geometrical counterpoise correction (gCP)[89] (entries 11 and 12). The calculated difference compared to the reference protocol is less than 1 kcal/mol, showing thus that the free energy of Equilibrium (1-C10) is insensitive to the BSSE.

V. Free energy correction For the inclusion complexes studied here, typically around ten vibrational modes feature eigenfrequencies below 50 cm−1, which contribute significantly to the free energy computed using the harmonic-oscillator model of vibrational modes. In the reference computational protocol discussed above, the free energy correction is calculated using the quasi-rigid-rotorharmonic-oscillator (qRRHO) method.[33] To evaluate the sensitivity of the results to this choice, we also calculated the correction using the standard RRHO method. As seen in Table 1 entries 2 and 13, the two methods give results within 1 kcal/mol from each other. It has been argued that the translational and rotational degrees of freedom are significantly suppressed in the solution phase,[90-92] and therefore the free energy corrections should be somewhere between the RRHO value and the one calculated neglecting these degrees of freedom.[81] In principle, this effect should be accounted for by the continuum model and the inclusion of dispersion corrections.[93] Nevertheless, to quantify it, we have also considered the free energy correction using the qRRHO excluding the translational and rotational contributions (Table 1 entry 14). Interestingly, the energy changes by 4.1 kcal/mol in the wrong direction compared to the experimental value.

VI. Implicit solvation model The final part of the methodology to be evaluated concerns the choice of implicit solvation model. In the reference protocol above, the COSMO-RS approach, as implemented in the 13 ACS Paragon Plus Environment


COSMOtherm software, is used to calculate free energies of solvation. We have also considered two other widely used implicit solvent models, SMD [94] and C-PCM[95,96], and the results are given in Table 1 entries 15 and 16. SMD gives a very similar change in solvation free energy for Equilibrium (1-C10) as COSMO-RS (+61.1 vs +62.0 kcal/mol), while C-PCM predicts a much lower effect, +44.6 kcal/mol. With C-PCM, the overall energy of the equilibrium deviates therefore even more from the experimental value. To understand the source of the discrepancy between the different solvation methods, one can consider the urea molecule, which is akin to the imidazolone units at the rim of the cavitand, and for which the hydration free energy has previously been calculated to be −13.8 ± 0.4 kcal/mol using Monte Carlo simulations.[97] Also for this molecule, COSMO-RS and SMD give almost identical results (−11.3 and −11.2 kcal/mol, respectively), while C-PCM gives the significantly lower estimate of −7.8 kcal/mol. This strongly suggests that the difference between C-PCM and the other methods stems thus from the fact that the solvation free energy of the cavitand monomer is underestimated. These results also hint to the possibility that inclusion of explicit water molecules at the rim may provide a systematic way to improve the results. This will be investigated in detail below. Before that, however, it is important to discuss an aspect of the implicit solvation methods that can be problematic in the case of host-guest systems and that is related to the way in which the molecular cavities used in these solvation calculations are generated. Specifically, for the inclusion complexes studied here, two disjoint cavities are formed: one around the host and one around the bound guest. Between them there is a void, which the implicit solvent models decorate with point charges in order to model the electric field created by the solvent. The effect is thus that the incarcerated guests are modeled as if they were directly interacting with the solvent, which of course is an artifact of the solvation model. To estimate the effect of this on the solvation free energy of Equilibrium (1-C10), the two cavities were joined by adding extra spheres (using the ‘ExtraSph’ keywords in Gaussian) between the host and the guest. Quite reassuringly, with SMD it turns out that this is a quite minor problem, as the procedure described above gives a negligible effect of +2.3 kcal/mol, (Table 1, entries 15 and 17). With C-PCM, the effect is −5.7 kcal/mol (Table 1, entries 16 and 18), i.e. in the wrong direction compared to experiments.


Page 14 of 29

Page 15 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


VII. Mixed explicit-implicit solvation model The imidazolone units at the rim of 1 are good hydrogen bonding donors and acceptors, and many hydrogen bonds can be formed with the water solution. It has been suggested that the dehydration of these moieties in the dimerization process can be a driving force for the complexation with hydrophobic guests, with a gain of entropy when the water molecules previously bound to 1 are released to the bulk.[22] Water molecules have also been suggested to play structural roles at the rim of 1 [22,24] and when the 12 capsule is widened to accommodate long, folded alkanes.[25] While such effects should in principle be captured by the continuum model, the results above concerning e.g. urea solvation by the implicit solvation models indicate that the lack of explicit solvent molecules might be the source of the large discrepancy between the reference computational protocol and the experimental value. In this section, we will therefore consider mixed explicit-implicit solvation models. In Equilibrium (1-C10), the two C10@1 monomers have a number of explicit water molecules bound at the rim, and these must be released to the bulk upon dimerization to form C10@12. It is therefore important to develop some criterion for knowing how many explicit water molecules are needed in the model, and also to have an accurate way of treating them. Goddard et. al. have proposed that comparison of species with different hydration numbers can be done by treating surplus water molecules on either side of the equilibrium as water clusters.[59] This procedure leads to hydrogen-bonded complexes of similar size on both sides of the equilibrium, which results in the cancellation of systematic errors in the calculated free energies that are due to the model size.[59] A shortcoming of this approach is that the calculated stabilities of water clusters of different sizes can differ (see Figure 4a), introducing an error that will depend on the size of the water cluster. For example, an (H2O)16 cluster is calculated to be 11.2 kcal/mol higher in energy than two (H2O)8 clusters with the reference computational protocol. It has also been pointed out by Goddard et. al. that water clusters of size (6+4N) are calculated to be more stable than those of size (4+4N).[59] In the present study, we use a mixed explicit-implicit solvation method subject to the criterion that formation of water clusters in water, i.e. Equilibrium (3), should be thermoneutral. That is, addition of one more water molecule to a water cluster – in water – should take place with no gain or loss in free energy, which renders the specific size of the water cluster irrelevant. N×H2O (aq) ⇌ (H2O)N (aq) 15 ACS Paragon Plus Environment

Equilibrium (3)


Interestingly, using the reference protocol presented above, the C-PCM solvation method satisfies this condition while COSMO-RS and SMD overshoot, as shown in Figure 4a (open symbols).[98] A more elaborate computational scheme including both the counterpoise correction and a free energy correction term RTln(N!), where N is the number of exchangeable molecules in the complex,[99] brings the calculated free energies of the cluster formation closer to thermoneutrality with COSMO-RS and SMD, but shifts the C-PCM values away from zero (Figure 4a, filled symbols).

Figure 4. a) Calculated free energies of formation of water clusters, Equilibrium (3), using the reference protocol (open symbols) and with counterpoise and exchangeability corrections (filled symbols). Three implicit solvation models are used: COSMO-RS, SMD, and C-PCM. Lines are fitted to the filled symbols. b) Calculated free energies of cluster formation after correcting the solvation free energy of individual water molecules as explained in the text.


Page 16 of 29

Page 17 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Interestingly, the deviation from thermoneutrality increases almost linearly with the number of water molecules. This suggests that the error stems from the description of the solvation free energy of a water molecule in water, and a simple correction term could be sufficient to adjust the trend to become thermoneutral. Indeed, by line fitting, the correction term is calculated to be +0.81, +1.93 and −0.45 kcal/mol for COSMO-RS, SMD and C-PCM, respectively, such that the hydration free energy of a water molecule becomes −5.0, −5.8, and −4.8 kcal/mol with the respective solvation model. The addition of this term leads to very small deviations from thermoneutrality, as shown in Figure 4b. An important advantage of this approach is that water molecules can be added one by one to any complex, and the total free energies calculated using the parameterized hydration free energy of water can be compared to assess whether this water molecule will interact more strongly with the studied complex or instead remain in the bulk solution. The model with the lowest total free energy would constitute a better representation of the system, and water molecules calculated to bind less well than to the bulk are accordingly not included in the model.[100] Using this revised approach for calculating the solvation free energy of water, explicit water molecules were added to C10@1 and their apparent binding free energies evaluated in order to assess whether they need to be included in the model of this complex or not, i.e. if addition of the water molecules is calculated to be exergonic or endergonic. Both the exchangeability correction and the quite expensive counterpoise correction are included in these calculations, and the SMD solvation model is chosen since it gave the best linear fit to the calculated free energies of cluster formation (Figure 4). First, a complex with four water molecules added to C10@1 was considered, called C10@1·(H2O)4, in which each water molecule resides between two benzimidazolone panels of the vase cavitand (Figure 5). The complex becomes more symmetric and compact compared to that without explicit water molecules (Figure 3). However, the directionalities of the hydrogen bonds are not ideal, and the four water molecules are in total calculated to be unbound by 0.5 kcal/mol compared to C10@1. Next, the addition of four more water molecules was considered, to form the C10@1·(H2O)8 complex in which two water molecules are bridging each pair of benzimidazolone side panels in the cavitand. In the lowest energy conformation of this structure, the water molecules were found to be arranged as shown in Figure 5. Binding of the



water molecules in this manner is calculated to be quite favorable and the complex is calculated to be 12.6 kcal/mol lower in energy compared to C10@1 without water molecules, i.e. each water is calculated to be bound by 1.6 kcal/mol. Addition of four more water molecules, such that a network of three water molecules are now bridging each pair of benzimidazolones (Figure 5), results in further stabilization of the complex by an additional 5.0 kcal/mol compared to C10@1·(H2O)8. That is, each of the four new water molecules is bound by 1.25 kcal/mol.

Figure 5. Arrangement of water molecules around the cavitand complex with C10. The corresponding free energies relative to C10@1 without water molecules bound are given.


Page 18 of 29

Page 19 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The four three-water networks can be connected by a four-water ring on top to yield C10@1·(H2O)16, as shown in Figure 5. This complex is calculated to be 2.3 kcal/mol more stable than C10@1·(H2O)12, i.e. each new water molecule is bound by 0.6 kcal/mol. Finally, another set of four water molecules can be added in the way shown in Figure 5, such that the C10@1·(H2O)20 complex is formed. Each of these water molecules is bound by only 0.3 kcal/mol, such that C10@1·(H2O)20 is calculated to be 1.3 kcal/mol lower in energy than C10@1·(H2O)16. Hence the free energies can be considered converged with 20 water molecules per cavitand complex. When constructing these increasingly larger mixed explicit-implicit solvation models, an extensive conformational search had to be performed in order to confirm that the complexes presented are the lowest in free energy. Many different arrangements were generated by manually placing the water molecules in order to yield favorable hydrogen-bonding interactions with the other water molecules and with the rim of the cavitand. From these investigations it was found that the symmetric models shown in Figure 5 were the lowest in energy. Despite this thorough search, it is still possible that the reported structures do not correspond to the global minima. For the purpose of generality, it is here nevertheless assumed that the complexes presented in Figure 5 are representative. Moreover, a benefit of using the highly ordered complexes in Figure 5 is that the hydrogen-bonding networks proved to be stable with respect to different density functional methods, which simplifies further comparisons (vide infra). Interestingly, inclusion of explicit water molecules to the dimer complex was not found to give any stabilization. The most stable complex located between C10@12 and water molecules, in which 12 explicit water molecules are interacting with the seam of the inclusion complex, was calculated to be endergonic by 4.6 kcal/mol compared to the situation in which the explicit water molecules are replaced by implicit solvent (optimized geometry shown in Supporting Information). No explicit water molecules are thus used to model this species. The explicit water molecules do not only affect the geometry of the cavitand, but also the conformation of the C10 alkane within the host-guest structure. With 8 explicit water molecules, the all-gauche conformation of C10 is calculated to be slightly lower in free energy than the conformation with one anti interaction, which is calculated to be lower in energy when no explicit water molecules are included in the model. The calculated preference for the all-gauche conformer of C10 increases with the number of water molecules



Page 20 of 29

added, because it is seen in the geometry optimizations that the portion of C10 that is exposed to solvent causes more disruption to the hydrogen-bonding network of the solvent in the case of the anti-configuration than in the all-gauche.

Back to Equilibrium (1-C10), employing this new computational protocol, the addition of 20 explicit water molecules to each monomer complex results in a new equilibrium, (4-C10).

2×C10@1·(H2O)20 (aq) ⇌ C10@12 (aq) + C10 (l) + 40×H2O (aq)

Equilibrium (4-C10)

The calculated free energies for Equilibria (1-C10) and (4-C10) differ thus by two times the binding energy of 20 water molecules to C10@1. With the SMD solvation and the counterpoise and exchangeability corrections, the free energy of Equilibrium (1-C10) is calculated to be −46.9 kcal/mol. Stabilizing the left hand side of that equation by 2×21.2 kcal/mol yields thus a free energy of Equilibrium (4-C10) of −4.5 kcal/mol. This value is in excellent agreement with the experimentally derived value of −5 kcal/mol, which shows that the approach with the explicit-implicit solvation in combination with the parameterized value for the hydration free energy of water according to the protocol presented above represents a great improvement to the reference protocol discussed above. A final important point here is that the counterpoise correction with this large number of water molecules and the 6-311+G(2d,2p) basis set is quite demanding computationally. We have found that the much cheaper gCP correction yields almost identical results for Equilibrium (4-C10). In this case, the def2-TZVP basis set has to be used since the parameterization of the method was done for it. At this level of theory, the corrected value of δGsol(H2O) becomes −6.7 kcal/mol, and the corresponding values for the complexes shown in Figure 5 then become +1.4, −12.2, −18.3, −20.1 and −21.8 kcal/mol, respectively. A complex with 24 coordinated water molecules was also considered with this method and was calculated to be somewhat less stable than the complex with 20 water molecules (see Supporting Information). Overall, the free energy of Equilibrium (4-C10) is thus calculated to be –3.7 kcal/mol, which matches very well both the experimental value (−5 kcal/mol) and the energy obtained with the explicit counterpoise correction (−4.5 kcal/mol). To summarize this section, we have demonstrated that the large error associated with the reference protocol stems mainly from the implicit solvation model. An explicit-implicit


Page 21 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


solvation model is proposed, in which the number of explicit water molecules is determined by evaluating the free energy of the added waters compared to their hydration free energy in bulk according to a parameterized value, depending on the level of theory employed. This protocol is shown to yield good results for equilibrium involving the binding of n-decane. In the following section, we will examine how this new approach performs for other energies that can be deduced from the experiments.

VIII. Further assessment of mixed explicit-implicit model In addition to Equilibrium (1-C10), a couple of analogous equilibria can be formulated concerning the monomer-dimer relationships with the C9 and C11 alkanes.

2×C9@1 (aq) ⇌ C9@12 (aq) + C9 (l)

Equilibrium (1-C9)

2×C11@1 (aq) ⇌ C11@12 (aq) + C11 (l)

Equilibrium (1-C11)

As mentioned in the introduction, it was experimentally observed that for C9, monomeric complexes C9@1 are exclusively formed, while for C11 only encapsulation complex C11@12 is observed.[24] In relation to the experimental free energy of Equilibrium (1-C10), deduced above to be −5 kcal/mol, it is thus expected that Equilibrium (1-C9) is less exergonic (or even endergonic), whereas Equilibrium (1-C11) should be more exergonic. Following the explicit-implicit approach devised above, we have calculated that 16 explicit water molecules are bound to C9@1 in the most stable structure, and 20 explicit water molecules are bound to C11@1 (see energies of the complexes with different numbers of waters in Table S1 in the Supporting Information). This leads to the following revised equilibria:

2×C9@1·(H2O)16 (aq) ⇌ C9@12 (aq) + C9 (l) + 32×H2O (aq)

Equilibrium (4-C9)

2×C11@1·(H2O)20 (aq) ⇌ C11@12 (aq) + C11 (l) + 40×H2O (aq)

Equilibrium (4-C11)

The energies are thus calculated at the B3LYP-D3(BJ)/def2-TZVP level of theory, using the SMD solvation model, the gCP correction, the exchangeability RTln(N!) term, and qRRHO for thermal corrections. For comparison, the B3LYP-D2 and M06-2X methods with 21 ACS Paragon Plus Environment


Page 22 of 29

the same basis set are evaluated. The results are presented in Table 2. Structurally, it is worth noting that for all three of C9@1, C10@1 and C11@1, the allgauche conformation of the n-alkanes is calculated to be the most stable one in the presence of explicit water molecules. Regarding the complexes with the 12 capsule, for C9 and C10 the all-anti conformation of the alkane was found to be the most stable one, while for C11@12 it was calculated that the energy is lowest when the alkane is shortened by featuring one internal gauche interaction. This is in line with experimental observations made for the capsule in nonpolar solvents, whereby gauche interactions were observed for encapsulated C11.[101] As seen in Table 2, using only the implicit solvation methods (Equilibria (1-C9), (1-C10) and (1-C11)) the three examined functionals reproduce the experimental trend in that the 2:1 cavitand:alkane complex becomes more stable relative to the corresponding 1:1 complex with increasing alkane length. However, the absolute values deviate by 25-40 kcal/mol from the experimental value. The corrected mixed explicit-implicit solvation treatment using B3LYPD3(BJ) reproduces the trend, with values in excellent agreement with the experiments. The total free energies are calculated to be +3.4, −3.7, and −10.9 kcal/mol, for the equilibria with C9, C10 and C11, respectively.

Table 2. Calculated free energies (kcal/mol) of experimentally deduced equilibria calculated with the new computational protocola and utilizing the parameterization approach for δGsol(H2O) introduced in Section VII. Method

Equilibrium 1-C9

Experimental value

4-C9

1-C10 4-C10 1-C11 4-C11 2-C10

> -5

-5

< -5

5-C10

Mixed ExplicitâImplicit Solvation Approach for Modeling of Alkane

Recommend Documents

Mixed ExplicitâImplicit Solvation Approach for Modeling of Alkane