Combined Computational and Experimental NMR Study of Calix [4

Article pubs.acs.org/JPCC

Combined Computational and Experimental NMR Study of Calix[4]arene Derivatives Vincenzo Verdolino,*,† Laura Baldini,‡ Ferruccio Palazzesi,† Federico Giberti,† and Michele Parrinello† †

Department of Chemistry and Applied Biosciences, ETH Zurich, and Facoltà di Informatica, Istituto di Scienze Computazionali, Università della Svizzera Italiana, Via G. Buffi 13, 6900 Lugano, Switzerland ‡ Dipartimento di Chimica Organica e Industriale, Università di Parma, Parco Area delle Scienze 17/A, I-43124 Parma, Italy S Supporting Information *

ABSTRACT: A combined computational and experimental study of a complex supramolecular system constituted by calix[4]arene derivatives that dimerize upon CO2 binding is presented. The theoretical investigation includes ab initio density functional theory, molecular dynamics, and metadynamics analysis of both monomers and dimers. The ab initio calculation of the dimerization energy demonstrates the exergonic character of the process, due to the formation of a strong hydrogen bond network between ammonium and carbamate groups. The dimerization is driven by −31.1 kcal/ mol in the case of the fully outward orientation of the carbamic hydrogens, while it results in a weaker process when different carbamic orientations are considered. The molecular dynamics simulations show the critical conformational degrees of freedom driving monomers and dimers toward common structures. These conformations show tilted orientations of the carbamic groups highlighting the fundamental role of dynamics in evaluating the most stable configurations. Metadynamics simulations describe, in agreement with the other computational tools, the conformational free energy surface of these calix[4]arenes defining three stable conformational families. ROESY and variable temperature 1H NMR experiments are in agreement with our simulations. The presented approach aims to be the reference for investigating complex supramolecular systems.

I. INTRODUCTION Modern computational techniques can rationalize ambiguous experimental evidence and, in some cases, suggest different interpretations. In particular, consolidated experimental routines and the extensive experience gained on a field can, sometimes, suggest apparently straightforward conclusions and interpretations. In most cases such conclusions are in agreement with experimental data and the general knowledge that scientists have on that particular subject. However, when the chemical system is particularly complex, experimental data can lead to partially incorrect conclusions. In this work we highlight this point focusing on a selected case which can be easily exported to many others in chemistry. One of the areas in which our laboratories are interested is supramolecular chemistry and, in particular, the harnessing of weak and reversible noncovalent interactions such as hydrogen bonds, hydrophobic interactions, metal coordination, and electrostatic interactions. Exhaustive analysis and reviews of these subjects are reported in literature.1 Also in this field, a significant contribution in understanding formation mechanisms, prevision of chemical properties, and rationalizations comes from theoretical and computational approach. In relatively few cases combined experimental and computational techniques have been employed in supramolecular chemistry. Among the examples we cite here the work of Mafra and co-workers2 who investigated the supramolecular © 2012 American Chemical Society

assembly of two different polymorphs employing NMR, XRD, and DFT calculations. Molecular dynamics combined with NMR, UV, and ESMS spectroscopy have been successfully used for post-self-assembly engineering of metal−organic supramolecular species,3 and quantitative evaluations of π−π interaction at high ab initio level of theory (CCSD(T) and DFT) have been employed in order to rationalize the effect of aryl substituent on stacking interactions.4 Conformational sampling through Monte Carlo (MC) algorithm has been employed by Gorrea and co-workers in order to investigate the ordered self-assembly aggregation of cyclobutane-containing dipeptides,5 whereas Ilot and co-workers investigated both solid state structures and dynamics of fluoronaphtalene derivates at an atomistic level.6 Here we extend the computational approach combining DFT and MD with enhanced sampling methods in order to demonstrate the critical role played by conformational degrees of freedom. An example where a complementary approach demonstrates how detailed and accurate the proposed procedure can be is in the case of CO2 molecular trapping by calix[4]arene derivatives. Tetraaminocalix[4]arene 1 (Scheme 1) rapidly and efficiently absorbs 2 equiv of CO2 thanks to the Received: August 8, 2012 Revised: October 9, 2012 Published: October 9, 2012 23441

dx.doi.org/10.1021/jp307895x | J. Phys. Chem. C 2012, 116, 23441−23452

The Journal of Physical Chemistry C

Article

Scheme 1. Reaction of 1 with CO2 and Following Dimerization of 2 → 3a

a

Two possible orientations of the carbamic hydrogens (inward “a” and outward “b”) are shown.

with the quantitative knowledge of the dimerization energy, the dynamics of both monomers and dimers in solution, and the exploration of the conformational free energy surface (FES), we gain a solid overview of these CO2 trapping scaffolds. Moreover, the similarity of the investigated systems with several other “capsulated” supramolecular settings8−11 extends our computational procedure to a wide scope of applications. This paper demonstrates how complementary techniques can address unsolved or apparently solved questions and define a solid computational approach in investigating supramolecular systems.

reaction of its four amino groups that form two ammonium and two carbamate groups (2) as recently reported by us.7 The reaction product, studied in solution by NMR spectroscopy (through a combination of 1H and 13C mono- and twodimensional techniques, including DOSY), is a hydrogenbonded dimer 3 of two bis(ammonium) bis(carbamate) intramolecular salts that is highly stable. All the NH ammonium and the CO carbonyl groups are involved in a seam of 12 strong intra- and intermolecular hydrogen bonds. Interestingly, we never observed a monofunctionalized 1-CO2 adduct. The addition of less than 2 equiv of CO2 resulted in the presence in solution of a mixture of free 1 and dimer 3. In ref 7 we attempted to assign the molecular structure speculating on the conformational orientation of the carbamic hydrogens in the dimer. Combining 2D-NMR data and a static model of the calixarene structure, we came to the conclusion that the four carbamic N−Hα (Scheme 1 top and bottom calixarene in species) 3 once oriented inward would represent the most stable conformation. However, some details of our experimental data were not fully accounted for by this model. Prompted by this we performed a combined computational and experimental investigation of the principal species involved. Unexpectedly, our results show that the most stable conformation is not the one previously proposed, but the one with the N−Hα facing outward. We also show that the system has a far greater flexibility than one could expect on the basis of standard DFT optimization. Our result is strengthened by the fact that the NOESY peaks can be better explained if one takes into account this conformational flexibility. The orientation of the NH group might appear as a minor detail. However, if we consider the role of the amino terminations in the dimerization process we can easily imagine how their flexibility may influence both stability and CO2 affinity. For this reason, future rational design of more efficient CO2 trapping molecules has to consider the complexity and the crucial role of conformational freedom for supramolecular systems. Moreover, a complete description of hydrogen-bonded systems, very common in supramolecular chemistry, must include a combination of NMR data and computational techniques as standardly done in biochemistry. Nonetheless,

II. COMPUTATIONAL DETAILS Ab Initio. Molecular orbital calculations were carried out using the G09 version of the GAUSSIAN series of programs.12 Geometries were optimized in the gas phase at the B3LYP density functional level of theory13−15 using the LanL2DZ basis set (D95 V on the first row of elements).16 Optimized geometries have been tested by means of the Hessian calculation to confirm effective convergence criteria at the minimum and in order to exclude any imaginary frequencies associated with saddle points. In order to take into account the solvent effects, single point energy calculations have been carried out on the gas phase optimized geometries by means of the polarizable continuum model (PCM)17−19 in chloroform. Starting geometries were taken from semiempirical optimizations and standard conformational Monte Carlo sampling.20 Molecular Dynamics. MD simulations were performed with AMBER 10.021 on both the monomers and the dimers. We employed the generic amber force field (GAFF)22 for all the simulations, and the atomic charges were calculated using the restrained electrostatic method (RESP)23 implemented in the antechamber module24 available in the AMBER 10.0 package.21 RESP charges have been calculated by Gaussian 09 on the previously optimized geometries taking into account the solvent effect as previously described for single point calculations. All the dynamics were run in explicit solvent setting up a truncated dodecahedron box of chloroform with a 10 Å buffer. In this way, at the start of the simulations, each molecule is separated by at least 20 Å from its periodic image. 23442



Article

Scheme 2. Front and Top View of Some of the Possible Relative N−H Bond Orientationsa

a (a) From the front view, bonds and atoms are reported in bold with respect to the transparent behind the plane σ; (b) from the top the arrows point inward/outward of the octagonal indicating IN or OUT configurations; oblique arrows indicate tilted displacements from ideal configurations. Different colors are used to identify different couples of hydrogens.

opens, we set a repulsive potential wall which forces the two distances to vary in the range of 4−13 Å. We ran MD-MTD simulations for 280 ns with 2 fs time steps keeping the temperature constant at 300 K using the Bussi et al.31 thermostat and a γ factor of 2 ps−1. We added Gaussianshaped bias to the system with hills of 0.25 kcal/mol (∼0.4 kBT) height and 0.02 Å widths every 1 ps. For the metadynamics calculations we employed the GROMACS4.5.1 package32−35 and the plug in PLUMED-1.2.2.36 The generalized amber force field GAFF has been converted in the GROMACS format through the available acpype.pl script.37,38

For monomers 2a and 2b (Scheme 1), the number of solvent molecules was resulting in 242 and 245 translating in a total number of atoms equal to 832 and 841 and fitted in cubes of sizes equal to 29.78 × 29.78 × 29.78 Å and 34.26 × 34.26 × 34.26 Å, respectively. In the case of the two dimers 3a and 3b, we added 662 and 695 solvent molecules leading to a total amount of atoms equal to 2198 and 2297 and, again, fitted in boxes of 41.68 × 41.68 × 41.68 Å and 42.70 × 42.70 × 42.70 Å dimensions. The slight difference in the number of solvent molecules did not affect the final result, and both systems ended toward common structures due to their flexibility. All the systems have been minimized employing the “sander” module in AMBER 10.0, and the long-range Coulomb interaction has been taken into account using a particle mesh Ewald25 with a cutoff for nonbonded atoms of 18 Å. In performing an initial relaxation, we constrained the solvent and then performed a full local minimization. We used an MD time step of 2 fs and constrained the hydrogen positions using the SHAKE algorithm.26 We used an NVT ensemble with the temperature controlled by a Langevin thermostat.27 The temperature was raised from 0 to 300 K in 50 ps using a restraint on the calixarene atomic positions of force constant equal to 10 kcal/mol·Å2. The last step in the equilibration procedure consists of removing the constraint and further running the MD for 100 ps. Finally we recorded the MD trajectories running the simulations for 100 ns in an NPT ensemble. Metadynamics. Free energy surfaces relative to the carbamic hydrogen orientations were computed by means of well-tempered metadynamics.28−30 We executed long MD metadynamics (MD-MTD) runs in order to reach convergence and to find the conformational minima. Collective variables (CVs) were selected in order to investigate the relative orientations of the carbamic hydrogens. Among all the possibilities investigated, the most suitable CV resulted to be the relative distance of the two facing carbamic hydrogens. Moreover, in order to avoid the sampling of extreme configurations in which for instance the system collapses or

III. EXPERIMENTAL DETAILS NMR Spectroscopy. 1H NMR spectra (400 MHz,) were recorded on a Bruker AV400 spectrometer using partially deuterated solvents as internal standards. Variable temperature 1 H NMR experiments were recorded in the range 193−298 K in CD2Cl2. Sample temperatures were controlled with the variable temperature unit of the instrument. Synthesis and characterization of 2 and 3 are elsewhere reported.7 IV. RESULTS AND DISCUSSION In the notation adopted in Scheme 2, the structures 3a and 3b (Scheme 1) are labeled as (IN/IN)−(IN/IN) and (OUT/ OUT)−(OUT/OUT) indicating that both of the two NH carbamates belonging to different calixarenes point inward and outward, respectively. These two are limiting cases of intermediate configurations having the four hydrogens oriented (IN/OUT)−(IN/OUT), (OUT/OUT)−(IN/IN), and (OUT/OUT)−(IN/OUT) as in Scheme 2. As we shall see below the flexibility leads to strong deviation from these idealized structures and allows for tilted intermediate configurations. In order to investigate the chemical properties of these systems, we perform several analyses: (a) we calculate the dimerization energy involved in the process reported in Scheme 1 at the DFT level; (b) we perform MD simulations for both monomers and dimers in explicit solvent (chloroform as 23443



Article

Scheme 3. Reaction of the Modified Calix[4]arene 4 with CO2 and Potential Dimerization of 5 → 6

Figure 1. DFT-optimized geometries of monomers 2a−2b (dft2a, dft2b) and dimers 3a−3b (dft3a, dft3b). For simplicity, we draw the halfstructure for the dimers including the interacting residues from the other half.

metadynamics procedure.28 As a way of testing our system we also considered the intermediates 5 and 6 reported in Scheme 3 where the calixarene has been modified allowing one CO2

in the experiments); (c) we measure VT 1H NMR and ROESY spectra of the dimer generated by 2; (d) we calculate the FES following an optimal collective variable (CV) by means of the 23444



Article

monomers. Comparing the two dimers, the most significant difference is associated with the longer intermolecular C O···H distances which are 1.80 Å (dft3a) and 1.61 Å (dft3b). This structural analysis demonstrates that both the inward and outward orientations represent local minima and that the most remarkable differences are related to the alkylic chain parameters and the monomer/monomer interface when they dimerize. Dimerization and Stabilization Energy. We focus now on the energy contribution involved in the dimerization process of each single monomer in chloroform calculated at the DFT/ PCM level of theory. In Table 2 we report the dimerization energies in kcal/mol of the monomers 2a and 2b and 5 (structure and energy details related to 5 and 6 are reported in the Supporting Information and labeled as dft5 and dft6). The data reported also include the energies calculated with no thermal contributions and zero point energy (ZPE) corrections (values in parentheses). Both the 2a → 3a and 2b → 3b processes are extremely exergonic due to the formation of several hydrogen bonds. In other words the process that is entropically disfavored due to the dimerization in one single unit is ruled by the enthalpy. In order to verify the capability of describing the hydrogen bond interaction at the presented level of theory, we analyze a modified model. We consider the same scaffold shown for 2 and 3 but reducing the number of possible hydrogen bond interactions (Scheme 3 structure 5). This compound has been already considered in our previous publication7 where we demonstrated experimentally that reducing the number of possible hydrogen bonds the dimerization process is not observed. In Table 2 we report the energy stabilization in gas phase. The data show that when decreasing the number of hydrogen bonds the stabilization energy decreases by more than 50% with respect to the 2b → 3b dimerization. Further calculations in chloroform do not add significantly to the analysis, and we are not going to discuss this model anymore. (Structure and details are reported in the Supporting Information, f1e and f1f.) The Gibbs free energy gained in the dimerization 2b → 3b is roughly 15 kcal/mol higher than the equivalent 2a → 3a (−31.1 vs −16.20 kcal/mol calculated in solution). The molecules involved in these two processes are differently stabilized by the solvent. The formation energy of 3a is favored by the presence of the solvent by 3.1 kcal/mol compared to the formation of 2a (i.e., −23.2 + 26.3 kcal/mol, see column ΔGsolv, Table 2). Similarly, for the dimer 3b, the stabilization is 3.9 kcal/mol higher if compared to 2b. Such a significant difference can be connected to the variation of the dipole moment and more in general to the electron reorganization. Indeed, in both processes the molecular dipole moment varies from 7.41−7.74 D to 0.04−0.11 D with consequent stabilization in a poorly polar solvent such as chloroform. As expected, the solvent effect in the dimerization process is small (∼3−4 kcal/mol) in comparison with the ΔGsdim (∼−16, −31 kcal/mol), and it does not justify the huge difference between the two dimerizations. The thermal contributions which include enthalpy, entropy, and zero point energy correction (ZPE) are 19.3 and 17.2 kcal/ mol for the 2a → 3a and 2b → 3b reactions respectively. In conclusion all together solvent effects and thermal contributions stabilize the 2b → 3b process over the 2a → 3a by ∼3 kcal/mol. Therefore, the difference of ∼15 kcal/mol that favors the 2b → 3b dimerization has to be explained in terms of other effects such as electronic structures, charge displacements,

reacting site only. It is to be expected that in this case the dimerization is weaker as indeed it turned out to be the case. We report below the description of the molecular structures of the DFT-optimized monomers 2 and dimers 3 (details for structures 5 and 6 in the Supporting Information). This is followed by the discussion on their relative stability and their dimerization energy and by the MD and MD-MTD analysis in comparison with the ROESY and VT- 1H NMR experiments. DFT Geometry Optimization. Although our interest is mainly focused on properties and behavior of the dimers in solution, it is natural to investigate the monomers first. Moreover, it is reasonable to assume that the dimerization process can take place only if the electronic properties and the proper orientation of the monomers allow the specific hydrogen bond interactions to take place. In Figure 1 we report the optimized geometries of the monomers and dimers. The two cases dft2a and dft2b represent two extreme conformations generated by rotation around the aliphatic carbon and the carbamic nitrogen bond. Both structures present a similar configuration of the four intramolecular hydrogen bonds. The flexibility of the alkylic chains allows two stable flattened cone conformations where the carbamic hydrogens can point either inward or outward in relation to the core of the cage (Scheme 2). This has significant consequences on the structural organization of the cage effecting the IN/IN which is more compressed and with shorter alkylic arms. Moreover, dft2a is flatter than dft2b as demonstrated by the shorter H1−H3 (labels as in Scheme 1 species 2) distance (4.32 Å in dft2a and 5.33 Å in dft2b, Figure 1). Another significant structural difference between the inward and the outward orientation is the distance between N1−Hα hydrogen and Oδ. In the structure dft2b the amidic hydrogen cannot be oriented toward the oxygen Oδ (distance equal to 5.12 Å) excluding the formation of the H bond. On the contrary this distance becomes shorter in structure dft2a (2.44 Å). This observation mainly drove our original conclusions reported in the previous publication7 suggesting a preferential structure due to the H··O stabilization. Finally the Oδ−Cα− Cγ−N1 dihedral angle switches from 11.2° to 120.9° converting the IN/IN in the OUT/OUT configuration denoting a significant modification of the alkylic chain itself. The other structural parameters are substantially unchanged for both monomers as shown in Table 1. The dimerization process does not substantially affect on the calixarene structure itself, as shown in Table 1. The only significant modification is associated with the expected stretching of the ammonium N−H and CO bonds which were not involved in the hydrogen bond network of the Table 1. Most Relevant Structural Parameters of the DFTOptimized Monomers and Dimers 2a, 2b, 3a, and 3ba

a

parameter

dft2a

dft2b

dft3a

dft3b

Oδ−O′δ H2−H4 H−N1 CO CO(h-bond) NH2−H(o-bond) NH2−H

6.41 11.20 1.01 1.26 1.35 1.08 1.02

5.90 10.97 1.01 1.26 1.35 1.08 1.02

6.42 11.09 1.02 1.29 1.33 1.06 1.05

5.99 11.01 1.01 1.29 1.32 1.05 1.07

Values are reported in Å. 23445



Article

Table 2. Dimerization Energy of Monomers 2a, 2b, and 5a ΔEgdim

b

structure 2a → 3a 2b → 3b 5→6

−36.2 (−55.6) −50.4 (−67.6) −20.2 (−34.9)

c

ΔGsdim

d

−16.20 (−35.5) −31.1 (−48.2)

ΔGsolv

e

−23.2 (2a); −26.3 (3a) −23.3 (2b); −27.2(3b)

ΔEtherm dim 19.3 17.2

f

ΔEgel

−4.7(2b−2a); −21.4(3b−3a)

ΔEsel

g

−4.8; −22.3

a

Variation energies are calculated at the B3LYP/LANL2DZ level of theory both in the gas phase and in solution. Geometries are optimized in the gas phase, and energies are corrected with PCM where applicable. Data expressed in kcal/mol. (Energies without thermal contributions reported in parentheses.) bElectronic energy of dimer − 2*electronic energy of monomer. cFree energy of dimer − 2*free energy of monomer in solution. dFree energy of “x” in solution − electronic energy of “x” in gas. eThermal contribution to the dimerization energy. fElectronic energy of “(2−3)b” − electronic energy of “(2−3)a” in gas (both monomers and dimers). gAs ΔEgel but in solution.

evidence is in agreement with what is known in terms of charge splitting and higher energy configuration (i.e., zwitterions) particularly in nonpolar solvents. Ammonium and carbamic nitrogens are also affected by the choice of the configuration: the NH3 nitrogens change from −0.5 (2a, 3a) to 0.0 (2b, 3b) of their partial charge, and following an opposite trend the N1 nitrogens experience electron enrichment from −0.7 (2a), −0.5 (3a) to −0.9 (2b) and −0.9 (3b). This evidence along with the shorter O···H distances in 3b (1.61 Å vs 1.80 Å Figure 1 dft3a and dft3b) support the conclusion that in the OUT/OUT configuration the hydrogen bond network is stronger. Other analysis including the electrostatic potential surface (EPS) mapped on the MK fitted charges supports this conclusion. (See in the Supporting Information the EPS for monomers and dimers.) Molecular Dynamics. In this section we report the analysis of the dynamics in solution of the monomers and dimers considered so far with the exception of structures 5 and 6. All the simulations start from the preoptimized DFT geometry described in the previous section and equilibrated as in the computational details. We analyze several parameters in order to identify predominant conformations and, eventually, transitions among conformations. Once all the monomers and relative dimers are equilibrated in solution we run 100 ns long MD simulation. (a). Monomers IN/IN (2a) and OUT/OUT (2b) Runs. During the equilibration process the structure of 2a underwent significant conformational changes displacing the two carbamic hydrogens in a configuration intermediate between the OUT/ OUT and the IN/OUT. A very similar structure was also reached in the OUT/OUT run without large conformational changes since the final structures of the IN/IN and OUT/OUT are much closer to the DFT-optimized OUT/OUT structure. The dynamical instability of the 2a structure is fully coherent with the DFT calculations which predict a difference in energy of 4.8 kcal/mol between the IN/IN (2a) and the OUT/OUT (2b) structures. After the initial equilibration both the IN/IN and the OUT/OUT runs follow similar dynamics as described in the Supporting Information. The most interesting motions of these compounds are associated with the flexibility of the alkylic arms, and it expresses itself mainly in the rotation of the terminal groups. Since in this motion the carbamic chains are more hindered, the ammonium and the carbamic arms follow different dynamics. The overlap of selected snapshots during the dynamics is useful to highlight the major structural changes encountered during the simulations (see Figure 2). It is possible to notice that there are two distinct groups of structures which differ mainly in the alkylic conformation of the carbamic arms. The dominant set of structures is stretched on the long direction which allows the orientation of the carbamic

conformational constraints, or hindering. Another interesting observation is about the electronic energy difference between dimers and monomers as follows: 4.8 kcal/mol (4.7 kcal/mol in gas) for 2a smaller than the analogous 2b, whereas 3a is 22.3 kcal/mol (21.4 kcal/mol in gas) less stable than the dimer 3b. This evidence suggests that a major role is played by the difference in the formation energy of the two dimers. Moreover, the electronic energy seems to be dominant, and the combination of a different charge distribution and some structural hindering can be the origin of the different stabilization. All these data coherently enforce the idea that the inward configuration, both in the monomer and in the dimer, would not be the preferred as thought until now. This energy analysis has been carried out at a higher level of theory treating the dispersion interactions with the empirical Grimme39,40 scheme and including corrections to the basis set superposition error (BSSE) confirming the same results and consolidating the aforementioned conclusions. (See the Supporting Information and Table 2bis in the Supporting Information.) An important contribution in understanding this system is given by a combined structural and electronic charge analysis at the interface of two interacting monomers. Species 3b, the most stable DFT structure, has the shortest intermolecular oxygen− hydrogen distance (1.61 Å, dft3b in Figure 1) denoting a stronger hydrogen bond network. This is supported by the Merz−Kollman (MK) charge analysis, and in Table 3 we report Table 3. Merz−Kollman Charges Calculated for Structures 2a, 2b, 3a, and 3ba

a

atom

2a

2b

3a

3b

N1 C(OO) O(C) O−h-bond(C) N(H3) H(bridging)

−0.7 1.2 −0.8 −1.2 −0.5 0.4

−0.9 0.9 −0.8 −0.4 0.0 0.0

−0.5 0.9 −0.8 −0.9 −0.5 0.4

−0.9 1.0 −0.9 −0.5 0.0 0.0

Charges in atomic units for the residues out of parentheses.

the most significant MK atomic partial charges (other representations and analysis are reported in the Supporting Information). The charge differences between the two functional groups −(OO)···N(H3) are 1.7 (2a), 0.9 (2b), 1.4 (3a), and 1.0 (3b) au. These data show that the electronic charge distribution of these two functional groups is more significantly affected by the relative carbamic IN or OUT orientation rather than the dimerization process itself. In other words the OUT/OUT configuration, both in the monomer and in the dimer, shows delocalized charges in contrast with the marked separation observed in the IN/IN structure. This 23446



Article

be accomplished without the need to disrupt the internal structure of the monomer. In conclusion the trajectories show that the hydrogen bond network is stable and persistent during the entire dynamics allowing an optimal orientation of the alkylic arms in order to eventually establish dimerization. The intramolecular hydrogen bond network does not freeze out the carboxylic functionalities leaving a considerable freedom. (b). Dimers (IN/IN)-(IN/IN) (3a) and (OUT/OUT)-(OUT/OUT) (3b) Runs. If for the monomers the calculated DFT energy gap between 2a and 2b suggested a conformational conversion 2a → 2b, the data obtained for the dimers make this suggestion even stronger. As we reported in the previous section (see Table 1), the energy difference between 3a and 3b is 22.3 kcal/ mol. Unless the two conformations were separated by inaccessible energy barrier at room temperature, the large difference between 3a and 3b has to be observed in the MD simulations. As noticed for the monomer, we can observe two major structures characterized by two different alkylic chain conformations (see “A” and “B” in Figure 4). Also in the case of the dimer the alkylic chain has enough flexibility to allow conformations similar to those discussed above (see Figure 3). One of the parameters that reflects the relative carbamic orientations is the distance between the N1−H hydrogens. After the equilibration of the initial (IN/IN)−(IN/IN) structure (3a), the distances of the two N1−H are close to the one optimized at the DFT level (∼5.6 Å). This conformation turns quickly into an IN/OUT (∼8.0 Å). In about 10−15 ns the structure of the dimer switches between the IN/OUT and the dominant OUT/OUT (∼10.0 Å). This last configuration is finally the one reached at the equilibrium. The dynamics followed by the starting OUT/OUT dimer (3b) is much more stable and quickly equilibrates converging into the same configuration reached by 3a (see Figure 5a, distance d2 in black). Interestingly we found that the motion of the two carbamic hydrogens of 3b and the distance between the two calixarenes (see caption in Figure 5) are correlated as shown in Figure 5a (comparison between d1 and d2). This evidence is supported by the simultaneous transition of the distance d1 and d2 by ∼1 and ∼2 Å respectively and associated with the (OUT/OUT)−(OUT/OUT) → (IN/OUT) transition. As observed in the monomer the alkylic chains act as a spring playing on the calixarene−calixarene distance and allowing the carbamic arm conformation to switch. Further supports to this thesis arise analyzing the Oδ−Cγ−Cβ−Cα

Figure 2. Overlap of the structures of monomer 2 sampled during 100 ns long MD starting from the DFT-optimized structure IN/IN. It is possible to recognize two principal situations having the alkylic chain either in a “compressed” (A) or in a “stretched” (B) conformation. In green bonds Cγ−Cβ in order to highlight the alkylic conformational change. (Picture made using VMDVisual Molecular Dynamics.41)

hydrogens outward. Instead, those structures that are less relevant and appearing in the first 5 ns are shorter in the alkylic chain and characterized by an irregular carbamic hydrogen orientation (tilted IN/OUT and rarely IN/IN). Following this motion, the calixarene structure acts as a spring allowing elongation and compression of the entire structure within ∼1.4 Å (see Figure 3). Together with this motion we also highlight a significant distortion of the calixarene macrocycle where the typical flattened shape continuously changes. This combined motion of the carbamic chain and of the aromatic macrocycle allows the terminal carboxy groups to vary their orientation (follow blue and red arrows in Figure 3) while at the same time maintaining energetically favorable hydrogen bond interaction. This property is crucial in the formation of the dimer which can

Figure 3. Representation of the two limiting conformations of the carbamic arms studied in this paper. The two orientations are permitted by previous alkylic chain reorganization inducing the 180° rotation along the Cα−Cβ bond and consecutive alkylic chain elongation. 23447



Article

and one (OUT/OUT)−(IN/OUT) (3e). In addition, since the (IN/IN)−(OUT/OUT) configuration was never encountered during the MD runs, we also considered the structure 3f which mostly approximates this setting (see Supporting Information for details on all the considered structures). In Table 4 we report the single point energies corrected in the meaning of thermal contributions and referred to the structure 3a and 3b. As shown, in all cases the (IN/IN)−(IN/IN) orientation is thermodynamically less stable than all the other conformations. The (OUT/OUT)−(OUT/OUT) structure is still the most stable setting. However, it is possible to notice that all the intermediate configurations extracted from the MD trajectories are reasonably close in energy to the (OUT/OUT)−(OUT/ OUT) structure, and this is coherent with the dynamics observed. Combined NMR Experiments and MD Analysis. Intrigued by the apparent mismatch between our MD simulations and what was previously believed to be the structure, we perform a combined ROESY NMR and MD analysis of 3. If the NH group pointed inward, all three sets of alkylic hydrogens would interact significantly with the carbamic hydrogen (DFT distances equal to 4.2, 2.4, and 4.2 Å from δ, ε, and ζ; Figure 6, right panel). An NMR spectrum has been recorded in CDCl3, and it shows correlation peaks between the carbamate NH protons and the ε and ζ methylene protons but not with those in δ (Figure 6, left panel). This is the proof that this group does not necessarily point inward as expected by the DFT-optimized geometry. However, an ideal (OUT/OUT)− (OUT/OUT) configuration is not possible as well because the relative distances would be longer than 5 Å. These results support that the system can be in several intermediate orientations spanning around the IN/OUT or OUT/OUT and many of them detected by NMR two-dimensional techniques. The question is if there is any fast conversion between IN → OUT orientations that cannot be detected by the NMR time scale. This is intrinsically related to the rotational degrees of freedom of the carbamic arms. In order to shed light on this point, we perform variable temperature (VT) 1H NMR experiments of 3 (Figure 7). The resulting data show that by lowering the temperature the signal of the ammonium NH3+ protons passes from a fast exchange to a slow exchange regime. In fact, below the coalescence temperature (∼200 K), the two protons involved in the intramolecular hydrogen bonds resonate at higher fields than the proton involved in the intermolecular interaction giving rise to two signals of 2:1 ratio (top left magnification). On the contrary, the carbamate NH signal does not undergo any modification throughout the whole temperature range. This can be interpreted in two ways that exclude each other: (a) the carbamate has one predominant configuration (probably an intermediate IN/OUT based on the previous ROESY analysis) generating an unchanged signal throughout the VT experiment; (b) the carbamate rotation (flipping) is still fast at 193 K. However, we believe that this last case can be excluded. Indeed, the two principal effects acting against the free rotation of the −NH3 and the carbamic group are due to the hydrogen bonds established by the functional group and the inertial momentum. Both functional groups are involved in three hydrogen−oxygen interactions of comparable energy, and the carbamic structure is much heavier and rotates more slowly around the CH2−N bond. If decreasing the temperature down to 193 K the rotation of the −NH3 is frozen,

Figure 4. Overlap of the structures of dimer 3 sampled each 120 frames during 100 ns long MD starting from the DFT-optimized structure (IN/IN)−(IN/IN). It is possible to recognize the two major situations having the alkylic chain either in a “compressed” (A) or in a “stretched” (B) conformation as reported for the monomer. Colors as in Figure 2. (Picture made using VMDVisual Molecular Dynamics.41)

dihedral angle as shown in Figure 5b and 5c where it is clearly shown that the (OUT/OUT)−(OUT/OUT) conformation can occur only when the alkylic arm is elongated and the Oδ− Cγ−Cβ−Cα dihedral is ∼60°. In conclusion the trajectories show that the dimer is stable during the dynamics, and we do not observe any dissociation. The high stability of the dimer is supported by an extensive hydrogen bond network which involves both inter- and intramolecular interactions between two calixarenes. Moreover, a behavior similar to that found in the monomer is observed also in the dimer which shows elongation/compression of the structure along the axis due to alkylic arm conformation changes. Finally the carbamic hydrogen distances, the two calixarene distances, and the alkylic arm conformation are well correlated with each other. (c). Ab Initio Analysis on Selected Configurations. Having performed static DFT optimizations which showed that there is a big energy gap between the two extreme conformations, we want to reconcile this result with the observed conformational flexibility. This is in order to make sure that the conformational picture is not an artifacts consequence of the force field used in the MD simulation. We selected from the dynamics two representative (IN/OUT)−(IN/OUT) structures (3c and 3d) 23448



Article

Figure 5. Plot of 100 ns dynamics related to (a) the relative N1−H···N1−H′ distances regarding the dimer 3b (black line d2) and the distances of the center of masses of the two opposite sets of four etheric oxygens (red line d1); (b) comparison with the Oδ−Cγ−Cβ−Cα dihedral angle; (c) correlation plot between the distance d2 and the dihedral angle Oδ−Cγ−Cβ−Cα. All the data are considered in their absolute value due to periodic conditions.

Cζ of the ammonium arms. Based on a maximum threshold of NOE’s effect equal to 5 Å, we determine the plot reported in Figure 8. The trend in Figure 8 is in overall agreement with the experimental data being substantially below or within the 5 Å threshold for both hydrogens in ε and ζ. Metadynamics. In Figure 9 we plot the distribution of relative orientations of the carbamic N−H residues obtained by reweighting algorithm.42 Such orientations are analyzed in function of two angles α1 and α2 defined as the dihedral angles made by Oδ−Oδ′−N1−Hα (Scheme 2) of two alkylic chains. These two angles define an OUT/OUT conformation when α1 = α2 = 180°, an IN/OUT when α1 = 180° and α2 = 0° or α1 = 0° and α2 = 180°, and finally an IN/IN when α1 = α2 = 0°. It is seen that the dominant conformation is one of the (OUT/ OUT)−(OUT/OUT) type. Note, however, the strong anharmonicity which would correspond to the appearance of many tilted conformations. The other minima correspond in order to IN/OUT arrangement and less likely to IN/IN structure. While these configurations exhibit large fluctuations from the ideal arrangement, it is clear that three different states of the calixarene can be identified. Note that for symmetry reasons the IN/OUT and the OUT/IN configurations are equivalent, and they represent the same physical minimum.

Table 4. Single Point Relative Energy in kcal/mol of Intermediate Structures Referred with Respect to Either dft3a or dft3b structure

vs dft3a

vs dft3b

3c 3d 3e 3f

−16.4 −18.9 −18.2 −5.9

5.2 2.6 3.4 15.7

we would expect a similar behavior from the carbamate function. It is for this reason that we believe that the dominant configuration is an intermediate IN/OUT and OUT/OUT in which the NH hydrogens are oriented toward the opposite alkylic chains albeit with some degrees of freedom. In order to prove that the MD simulations showed so far are in agreement with these new and previously published data,7 we can extract the information needed to predict the average N1− H···H−N1 interactions. In particular we demonstrate that these simulations agree with the NMR measurements also confuting the absolute assumption that only the (IN/IN)−(IN/IN) configuration can be matched. Here we analyze the distance between the carbamic N1−H and the hydrogens in the Cε and 23449



Article

Figure 6. Portion between 1.5 and 5.5 ppm of the ROESY NMR (400 MHz) spectrum of 3 recorded in CDCl3. Highlighted in color are the offdiagonal correlation peaks. On the right side is reported a schematic representation of 3 including the nomenclature used in this section.

Figure 7. 1H NMR experiments of 3 recorded in CD2Cl2 at variable temperature in the range between 298 and 193 K. Lowering the temperature below 240 K, the broadening of the −NH3 signal above 9.0 ppm is in contrast with the −NH signal between 5.0 and 4.5 ppm which is not modified. Magnification of signals between 11 and 9 ppm reported on the top left of the plot for the spectrum recorded at 193 K.

MD simulations point out the necessity of including dynamical effects. The static equilibrium DFT structures provide only a single conformation which is not a very realistic description of these molecules in solution. If we set the zero energy with the most stable (OUT/OUT)−(OUT/OUT) configuration, then the other minima are at 1.8 and 4.1 kcal/ mol for the IN/OUT and the (IN/IN)−(IN/IN) structures, respectively. The transition from the (IN/IN)−(IN/IN) into one of the IN/OUT setting is substantially barrierless, the requested energy being lower than 0.4 kcal/mol which is lower than kbT at room temperature (∼0.6 kcal/mol). The structure (OUT/OUT)−(OUT/OUT) is obtained providing almost 1 kcal/mol. This analysis very well explains the reason why MD simulations starting from an IN/IN configuration turn quickly into an IN/OUT or into an OUT/OUT conformation. Analyzing the process backward any (OUT/OUT)−(OUT/ OUT) structure needs 1.8 kcal/mol in order to switch only one orientation of the carbamidic hydrogens and a further 3.8 kcal/ mol in order to set its configuration in the (IN/IN)−(IN/IN) orientation. (See FES in the Supporting Information.) The

Figure 8. Plot (100 ns) of the distance between the N1−H and Cξ (a) and Cε (b) hydrogens. The horizontal line marks the upper threshold defining the longest H···H distance in order to record a correlation signal during the ROESY NMR experiment.

(OUT/OUT)−(OUT/OUT) and the possible IN/OUT structures are at room temperature in equilibrium allowing great flexibility as suggested by the marked anharmonicity of each minima. This implies that experimentally we would expect to observe an average situation of these fast motions among the two deepest wells. This conclusion is in agreement with our 23450



Article

10. This analysis also explains the reason why the carbamate flipping observed in our simulations never happens in one step but following a multistep progression via intermediate configurations.

V. CONCLUSIONS In this paper we report the results obtained by ab initio DFT calculations on the dimerization energy of calix[4]arene derivatives upon CO2 binding. These supramolecular dimers showed large stabilization energy as high as −31.1 kcal/mol in agreement with the experimental observations. The static DFT approach and the decomposition analysis in several contributions suggest that the OUT/OUT conformation is always favored over the IN/IN. This conclusion is verified for both the monomer and the dimer. Employing a different scaffold with reduced CO2 attaching sites, we were able to confirm that the hydrogen bond network is the major trading force in the dimerization process. Moreover, the electronic properties and the electron charge distribution play a critical role in the inward/outward configuration as demonstrated by the MK population analysis. MD simulations confirm the fragile stability of the IN/IN configuration which undergoes conversion into tilted structures just after the equilibration process converging into a common configuration when starting from the OUT/OUT structure. Similar behavior has been observed simulating the dimer which, after relatively short time, fluctuates between (OUT/OUT)− (OUT/OUT) and tilted configurations. None of our MD simulations showed recrossing between (IN/IN)−(IN/IN) and (OUT/OUT)−(OUT/OUT) structures. Such evidence has been supported by VT 1H NMR and ROESY 1H NMR in good agreement. Indeed, VT 1H NMR experiments enforce the hypothesis that (at least within the NMR time scale) the carbamic arms occupy a preferential configuration instead of flipping quickly between inward and outward. Moreover, 2DNMR draws a clear picture which is in full agreement with MD simulations placing the four carbamic arms of the dimers in a tilted fashion. We also employed, for the first time, the metadynamics on a supramolecular system showing how powerful and also qualitatively useful it can be in order to investigate a relatively large conformational space. Moreover, quantitative previsions

Figure 9. Free energy surface of the O−O−N−H dihedral angles α1 and α2 reweighted on the distance NH···HN CVs. Isoenergy lines have been plotted every 0.5 kcal/mol. (For further details see the Supporting Information.)

MD simulations, the DFT analysis, and our NMR experiments. Our simulations also allow to investigate the correlation between the configurations of the two calix[4]arenes in the dimer. In Figure 10 we plot two sections of the FES discussed above representing the configurations seen by one of the two calixarenes when the other is populating the state (a) (IN/IN) or the state (b) (OUT/OUT). (Populating the IN/OUT configuration leads to an equivalent FES of the OUT/OUT reported in case “b”.) When the configuration of one calixarene is (IN/IN) the configuration of the second calixarene is (IN/ OUT), which is roughly 2 kcal/mol more stable than the (OUT/OUT) (Figure 10a). On the contrary when one calixarene is in the (OUT/OUT) configuration (Figure 10b), the second calixarene has almost the same probability to be in the (IN/OUT) or (OUT/OUT) state due to a very small energy gap (less than 1 kcal/mol). Such correlation can be translated in a configuration induction due to the fact that when one calixarene populates a given state the other is forced to change its configuration following the FES reported in Figure

Figure 10. Free energy surface projection of the O−O−N−H dihedral angles (α1 for simplicity) reweighted on the distance NH···HN CVs. When one of the two calixarene is set in the (IN/IN) configuration (panel a), the other one prefers to populate the configuration (IN/OUT) on the (OUT/OUT) by roughly 2 kcal/mol. On the contrary (panel b) when one calixarene is set as OUT/OUT, the other one has the same probability to be in the (IN/OUT) or (OUT/OUT) configuration. 23451



Article

(9) Conn, M. M.; Rebek, J. Chem. Rev. 1997, 97, 1647−1668. (10) Hof, F.; Iovine, P. M.; Rebek, J. Abstr. Pap., Am. Chem. Soc. 2001, 221, U122−U122. (11) Ajami, D.; Rebek, J. Angew. Chem., Int. Ed. 2008, 47, 6059− 6061. (12) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09; Gaussian, Inc.: Wallingford, CT, 2009. (13) Becke, A. D. Phys. Rev. A 1988, 38, 3098−3100. (14) Becke, A. D. J. Chem. Phys. 1993, 98, 5648−5652. (15) Lee, C. T.; Yang, W. T.; Parr, R. G. Phys. Rev. B 1988, 37, 785− 789. (16) Dunnning, T. H., Jr.; Hay, P. J. Modern Theoretical Chemistry; Schaefer. H. F., III, Ed.; Plenum: New York, 1976; Vol. 3, pp 1−28. (17) Cances, E.; Mennucci, B.; Tomasi, J. J. Chem. Phys. 1997, 107, 3032−3041. (18) Cossi, M.; Barone, V.; Mennucci, B.; Tomasi, J. Chem. Phys. Lett. 1998, 286, 253−260. (19) Mennucci, B.; Tomasi, J. J. Chem. Phys. 1997, 106, 5151−5158. (20) Hastings, W. K. Biometrika 1970, 57, 97−109. (21) Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. J. Comput. Chem. 2005, 26, 1668−1688. (22) Wang, J. M.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157−1174. (23) Cieplak, P.; Cornell, W. D.; Bayly, C.; Kollman, P. A. J. Comput. Chem. 1995, 16, 1357−1377. (24) Wang, J. M.; Wang, W.; Kollman, P. A.; Case, D. A. J. Mol. Graphics Modell. 2006, 25, 247−260. (25) Essman, U.; Perela, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577−8592. (26) Gunsteren, W. F. v.; Berendsen, H. J. C. Mol. Phys. 1977, 34, 1311−1327. (27) Pastor, R. W.; Brooks, B. R.; Szabo, A. Mol. Phys. 1988, 65, 1409−1419. (28) Barducci, A.; Bussi, G.; Parrinello, M. Phys. Rev. Lett. 2008, 100. (29) Laio, A.; Parrinello, M. P. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562−12566. (30) Barducci, A.; Bonomi, M.; Parrinello, M. Wires Comput. Mol. Sci. 2011, 1, 826−843. (31) Bussi, G.; Donadio, D.; Parrinello, M. J. Chem. Phys. 2007, 126. (32) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435−447. (33) Van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701−1718. (34) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306−317. (35) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Commun. 1995, 91, 43−56. (36) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; et al. Comput. Phys. Commun. 2009, 180, 1961−1972. (37) Batista, P. R.; Wilter, A.; Durham, E. H. A. B.; Pascutti, P. G. Cell. Biochem. Biophys. 2006, 44, 395−404. (38) Sousa Da Silva, A. W.; Vranken, W. F.; Laue, E. D. Manuscript to be submitted. (39) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. J. Chem. Phys. 2010, 132, 154104−154123. (40) Grimme, S.; Ehrlich, S.; Goerigk, L. J. Comput. Chem. 2011, 32, 1456−1465. (41) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics Modell. 1996, 14, 33−38. (42) Bonomi, M.; Barducci, A.; Parrinello, M. J. Comput. Chem. 2009, 30, 1615−1621.

have been proposed on the FES which agree both with accurate ab initio techniques and force field based MD. In conclusion we presented an efficient employment of several well-established techniques but that rarely have been combined in order to explore sensitive supramolecular properties and parameters such as the electronic structure and the conformational space. The case presented is an interesting example which shows how combined techniques can shed light on ambiguous trusting or unclear experimental evidence. Indeed, we believe that the same procedure can be applied on a much broader scope in supramolecular chemistry allowing more trustable experiment rationalization, molecular engineering, and chemical functionalization.

■

ASSOCIATED CONTENT

S Supporting Information *

DFT-optimized structures, xyz atomic coordinates, energy contributions, and Hessian analysis; population analysis and RESP representation; most significant parameter analysis on the MD trajectories along with structure overlap analysis and RMSD; FES of selected CVs. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Author Contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Alessandro Barducci for fruitful discussions. Calculations have been carried out on a BRUTUS cluster at ETH Zurich and are supported by a grant from the Swiss National Supercomputing CenterCSCS under project ID s312 and s309.

■

ABBREVIATIONS DFT, density functional theory; MD, molecular dynamics; MDMTD, molecular dynamics-metadynamics; CVs, collective variables; FES, free energy surface; GAFF, generic Amber force field; RESP, restrained electrostatic potential; PCM, polarizable continuum model; MK, Merz−Kollman.

■

REFERENCES

(1) Supramolecular Chemistry: From Molecules to Nanomaterials; Gale, P. A.SteedJ. W., Eds.; Wiley & Sons: Chichester, U.K., 2012. (2) Mafra, L.; Santos, S. M.; Siegel, R.; Alves, I.; Paz, F. A. A.; Dudenko, D.; Spiess, H. W. J. Am. Chem. Soc. 2012, 134, 71−74. (3) Zheng, Y. R.; Lan, W. J.; Wang, M.; Cook, T. R.; Stang, P. J. J. Am. Chem. Soc. 2011, 133, 17045−17055. (4) Wheeler, S. E. J. Am. Chem. Soc. 2011, 133, 10262−10274. (5) Gorrea, E.; Nolis, P.; Torres, E.; Da Silva, E.; Amabilino, D. B.; Branchadell, V.; Ortuno, R. M. Chem. Eur. J. 2011, 17, 4588−4597. (6) Ilott, A. J.; Palucha, S.; Batsanov, A. S.; Wilson, M. R.; Hodgkinson, P. J. Am. Chem. Soc. 2010, 132, 5179−5185. (7) Casnati, A.; Baldini, L.; Melegari, M.; Bagnacani, V.; Dalcanale, E.; Sansone, F.; Ungaro, R. J. Org. Chem. 2011, 76, 3720−3732. (8) Tokunaga, Y.; Rudkevich, D. M.; Rebek, J. Angew. Chem., Int. Ed. 1997, 36, 2656−2659. 23452


Combined Computational and Experimental NMR Study of Calix [4

Recommend Documents