On the interaction between ions and complex aromatic systems: amino acid side chains, ion channels, and buckybowls

TESE DE DOUTORAMENTO Programa Ciencia e tecnoloxía química DEPARTAMENTO DE QUÍMICA-FÍSICA FACULTADE DE CIENCIAS LUGO SETEMBRO 2014 Alba Campo Cacharrón ON THE INTERACTION BETWEEN IONS AND COMPLEX AROMATIC SYSTEMS: amino acid side chains, ion channels, and buckybowls On the interaction between ions and complex aromatic systems: amino acid side chains, ion channels, and buckybowls Alba Campo Cacharrón Departamento de Química-Física Facultade de Ciencias LUGO, SETEMBRO 2014 Autorización dos directores da tese D. Enrique Manuel Cabaleiro Lago Profesor/a do Departamento: Química-Física (Lugo) D. Jesús Rodríguez Otero Profesor/a do Departamento: Química-Física (CIQUS) Como Directores da Tese de Doutoramento titulada: “On the interactions between ions and complex aromatic systems: amino acid side chains, ion channels and buckybowls”. Presentada por Dna. Alba Campo Cacharrón Alumna do Programa de Doutoramento en Ciencia e Tecnoloxía Química (D1121) Autorizan a presentación da tese indicada, considerando que reúne os requisitos esixidos no artigo 34 do regulamento de Estudos de Doutoramento, e que como Director da mesma non incurre nas causas de abstención establecidas na lei 30/1992. Lugo a 20 Setembro de 2014 Asdo: Enrique M. Cabaleiro Lago Asdo: Jesús Rodríguez Otero Agradecementos Quixera facer uns breves agradecementos posto que me gustan máis as demostracións con abrazos que con palabras, e máis, se estas hai que ir buscalas a libros tan pouco apetecibles para ler como este. En primeiro lugar aos directores de tese, Quique e Jesús. En especial a Quique, que máis que un director é un amigo. Gracias pola túa paciencia, os teus consellos e pola túa preocupación. Fas que todo pareza fácil e sinxelo e sempre poñendo un toque de humor. A Ortigueira, porque a tese sen esos cafés e crucigramas diarios ben seguro non houbera saído adiante. A meus pais, polo seu apoio constante, os seus ánimos e sobretodo o seu sacrificio para que eu puidera continuar e acabar esta tese aínda sen contrato. Son o meu exemplo. A meu irmán, que aínda que non falemos tan a miúdo como deberamos sei que conto co seu apoio incondicional ante calquera das miñas decisións. A meus tios e primos, porque son a mellor familia que un pode ter, sempre pensando en comer e en xogar o tenis. As miñas amigas, porque sen elas non sería quen son. Tamén a Raquel, parte fundamental destos anos en Lugo estudiando, porque sei que levo unha amiga para a vida. A Tamara e a Maritere, por ese duro pero fantástico ano traballando á vez que facía a tese. Polos nosos cotilleos ata altas horas da madrugada, polos chocolates e pola vosa preocupación pola miña vida amorosa. Gracias, porque sei que levo dúas amigas para sempre. Por último a Miguel, que aínda que foi o último en entrar na miña vida convertiuse nun dos máis importantes, sempre alentándome a estudiar, a mellorar, a emprender, sempre co seu apoio e cariño infinito. INDEX RESUMEN ................................................................................................................................. I 1. INTRODUCTION ................................................................................................................... 1 1.1. SUPRAMOLECULAR CHEMISTRY ................................................................................................... 3 1.2. INTERMOLECULAR INTERACTIONS ................................................................................................ 4 1.3. INTERACTION INVOLVING AROMATIC SYSTEMS ............................................................................... 9 1.4. INTERACTIONS WITH BUCKYBOWLS ............................................................................................ 17 1.5. REFERENCES.......................................................................................................................... 21 2. OBJECTIVES ........................................................................................................................ 25 3. METHODOLOGY ................................................................................................................. 31 3.1. INTERACTION ENERGY ............................................................................................................. 33 3.1.1. Basis Set Superposition Error ..................................................................................... 34 3.1.2. Many-body effects ..................................................................................................... 36 3.2. WAVEFUNCTION-BASED METHODS ........................................................................................... 37 3.2.1. The Hartree-Fock method .......................................................................................... 37 3.2.2. Many Body Perturbation Theory (MPn)..................................................................... 39 3.2.3. Coupled Cluster .......................................................................................................... 41 3.2.4. Errors in Wavefunction-based Methods .................................................................... 43 3.2.4.1. Extrapolating to basis limit ................................................................................................ 45 3.2.4.2. Obtaining benchmarking values........................................................................................ 47 3.3. DENSITY FUNCTIONAL THEORY METHODS ................................................................................... 49 3.3.1. Kohn-Sham procedure ............................................................................................... 50 3.3.2. Functional types ........................................................................................................ 52 3.3.3. Dispersion-corrected DFT methods ............................................................................ 54 3.3.3.1. DFT-D2 .............................................................................................................................. 55 3.3.3.2. DFT-D3 .............................................................................................................................. 56 3.4. REDUCING COMPUTATIONAL COST ............................................................................................. 57 3.5. INTERACTION ENERGY PARTITIONING ......................................................................................... 59 3.5.1. Energy Decomposition Analysis (EDA) Methods ........................................................ 59 3.5.2. Symmetry-Adapted Perturbation Theory Methods (SAPT) ........................................ 61 3.5.2.1. SAPT(DFT) ......................................................................................................................... 65 3.6. ELECTRON DENSITY ANALYSIS .................................................................................................... 66 3.6.1. NCI index .................................................................................................................... 67 3.7. SOLVATION EFFECTS ............................................................................................................... 71 3.7.1. Continuum models ..................................................................................................... 72 3.8. REFERENCES.......................................................................................................................... 75 iv menos estables, y de hecho no se han localizado en complejos con fenol debido a la tendencia de esta especie aromática a interaccionar mediante enlaces de hidrógeno. En resumen, los trímeros están principalmente condicionados por la intensidad de los contactos catión···π, pero interacciones secundarias entre especies aromáticas puede modular su comportamiento, especialmente cuando se forman enlaces de hidrógeno en complejos con fenol e indol. A pesar de esto, el balance de las diferentes contribuciones a la energía (π···π, X-H···π y M+···π) es muy delicado dependiendo de la naturaleza y de la orientación relativa de los fragmentos. El segundo bloque de esta tesis está dedicado al estudio de las interacciones no covalentes donde intervienen especies aromáticas y aniones. Así como las interacciones catión···π han sido estudiadas desde hace tiempo, el campo de las interacciones entre moléculas aromáticas y aniones es de reciente desarrollo. Esto es probablemente debido a que en una primera instancia la interacción entre un anión y una molécula aromática no parece posible debido a la capacidad de donar electrones de ambas especies. Sin embargo, pronto se hizo evidente que este carácter dador podía ser modulable en las moléculas aromáticas a través de la sustitución del anillo con grupos electroatrayentes. Así, se pueden definir las interacciones anión···π como las interacciones favorables entre aniones y sistemas aromáticos deficientes en electrones. Las interacciones anión···π han cobrado gran relevancia debido a su posible participación en importantes áreas como por ejemplo en bioquímica, ya que el ADN es un polianión y muchos cofactores y sustratos de las enzimas son aniónicos. Además, se han propuesto nuevas aplicaciones empleando este tipo de interacciones para actuar en receptores aniónicos. Por tanto, cuatro capítulos en esta tesis se corresponden con el estudio de la interacción anión···π: uno de ellos dedicado a las características de la interacción en un canal aniónico sintético y otros tres centrados en las características de las interacciones entre especies aromáticas curvas y aniones, principalmente para determinar su posible uso como receptores aniónicos. El capítulo 6 está dedicado al primer canal aniónico sintético basado en interacciones anión···π, que ha sido propuesto recientemente. Las posibilidades de este canal y su novedosa estructura compuesta por un motivo geométrico con anillos aromáticos deficientes en electrones repetido a lo largo del canal lo hacen perfecto para el estudio de las interacciones anión···π. Así, se han estudiado complejos formados por modelos simplificados del canal iónico interaccionando con cuatro aniones diferentes (Br-, Cl-, F-, y OH-) a los que se les han añadido hasta tres moléculas de agua de forma explícita. Los resultados obtenidos muestran que la interacción de las unidades que forman el canal con los aniones es intensa en fase gas, dando lugar a complejos muy estables. Sin embargo, la presencia de un pequeño número de moléculas de agua que pueda acompañar al ion dentro del canal altera de forma significativa las características de los complejos, especialmente los más estables, formados con los aniones más polarizantes. v Los resultados indican que el papel de las moléculas de agua más próximas al anión puede ser relevante en el funcionamiento del canal iónico, disminuyendo los costes asociados a la deshidratación del ión para entrar en el canal. Además, si algunas moléculas entran con el ion en el canal pueden contribuir a facilitar el proceso estableciendo interacciones atractivas con el propio canal iónico. Una última sección, formada por los tres capítulos finales, está dedicada a las interacciones anión···π en las que participan los sistemas aromáticos curvos llamados buckybowls. Los buckybowls son hidrocarburos aromáticos policíclicos formados por una serie de anillos fusionados de cinco y seis miembros que dan lugar a una estructura curvada en forma de cuenco. La curvatura tiene su origen en la propia estructura de estas moléculas, ya que no es posible formar una estructura plana combinando pentágonos y hexágonos. Este tipo de especies también son denominadas fragmentos de fulerenos, ya que su esqueleto carbonado se corresponde con partes de las estructuras de los fulerenos. Los buckybowls más sencillos que se pueden proyectar sobre el fulereno C60 son el coranuleno C20H10 (seis anillos hexagonales rodeando un anillo pentagonal central) y el sumaneno C21H12 (tres anillos hexagonales y tres pentagonales alternos en torno a un anillo hexagonal central). En esta tesis se considerarán derivados sustituidos de ambas especies. Los buckybowls muestran diferentes propiedades dependiendo de la cara cóncava o convexa en la que tiene lugar la interacción. Tanto sumaneno como coranuleno forman complejos con cationes y metales de transición, fundamentalmente por la cara convexa del bowl. Al igual que otras especies aromáticas, podría ser posible modular las características de estos bowls mediante sustituyentes apropiados, de modo que pudieran coordinar aniones preferentemente por la cara cóncava y actuar como receptores aniónicos. Este aspecto es uno de los objetivos a tratar en la presente tesis. En el capítulo 7 se intenta determinar cómo la sustitución de grupos en el borde de los buckybowls afecta a sus propiedades. Concretamente, se trata de determinar si se produce la inversión del potencial electrostático molecular (MEP) de los bowls, pudiendo ser así adecuados para un contacto favorable con aniones. Además, mediante este estudio se ha comprobado qué grupos sustituyentes son los más adecuados para favorecer la interacción buckybowl···anión. Se ha determinado el potencial electrostático del coranuleno sustituido con 5 o 10 grupos fluoruro, cloruro y nitrilo. La curvatura del coranuleno apenas sufre cambios al ser sustituido, excepto en los derivados decasustituidos con cloruro y nitrilo, en lo que se aprecia una pérdida de curvatura del bowl. Tras comprobar que se consigue la inversión del potencial electrostático (especialmente con CN), se optimizan los complejos formados por dichos buckybowls sustituidos y tres aniones diferentes, Cl-, Br- y BF4-. En dichos complejos se prueban diferentes regiones para la interacción con los aniones, comprobando así las diferencias energéticas entre complejos formados por las caras cóncava y convexa. Los contactos con el anión situado en la parte cóncava del vi bowl son los más estables con energías de coordinación que siguen aproximadamente los valores obtenidos para los MEPs, siendo los más estables aquellos formados por el coranuleno sustituido completamente con 10 grupos CN y el anión cloruro por la cara cóncava. Además, se ha estimado el efecto de varios disolventes empleando modelos de continuo, que muestran importantes pérdidas de estabilidad de los complejos en disolución. Aún así, los resultados animan a seguir investigando este tipo de sistemas, ya que parecen apuntar que los buckybowls sustituidos podrían actuar como receptores de aniones. El capítulo 8 está más centrado en aspectos metodológicos. Se ha detectado que la interacción con buckybowls es bastante compleja de describir, y que los distintos métodos de cálculo proporcionan resultados dispares. Como no existen valores de referencia para este tipo de sistemas, en este capítulo se han realizado cálculos con una batería de métodos de diferente nivel, para poder determinar valores de referencia así como para estimar qué métodos de menor coste computacional pueden ofrecer una descripción aceptable de estos sistemas. En este estudio se han considerado complejos con aniones y cationes, para poder establecer una comparación directa de sus características usando los mismos métodos con los mismos buckybowls. Se han considerado complejos formados por bowls parcialmente sustituidos con CH3, F y CN basados en coranuleno y sumaneno, e iones Cl- y Na+. En el caso del coranuleno se han sustituido cinco átomos de hidrógeno alternos en el borde del bowl. Con sumaneno hay dos posibles sustituciones, una que afecta a los átomos de hidrógeno de los grupos CH de los anillos hexagonales, y otra en la que se sustituyen los átomos de hidrógeno de los grupos CH2 de los anillos pentagonales. Se han obtenido las energías de interacción de los complejos a lo largo de una línea que pasa por el eje central del buckybowl variando las distancias a las que se sitúan los iones. Los resultados obtenidos con los diferentes métodos probados muestran un grado de correspondencia variable. Comparando con los valores obtenidos con el método de referencia empleado (MP2.X) se observa que el método SCS-MP2 extrapolado a base completa es el que mejores resultados ofrece a un coste de cálculo razonable, aunque otras opciones de menor coste como el funcional M06-2X también proporcionan una descripción aceptable. En general, el sumaneno interacciona más fuertemente que el coranuleno con ambos iones. Además, en el sumaneno el efecto de la sustitución es más pronunciado en los grupos CH que en los CH2. De forma general las interacciones siguen el mismo patrón que el observado para los MEPs, aunque pueden apreciarse ligeras desviaciones. Los resultados indican que la naturaleza de la interacción es diferente en complejos catiónicos y aniónicos. Así, la interacción con aniones está controlada por la interacción electrostática con contribuciones significativas de la dispersión, mientras que en los complejos con cationes son las contribuciones de inducción las que desempeñan un papel predominante. vii Por último en el capítulo 9 se estudia el efecto que la naturaleza del anión y la presencia de disolvente ejercen sobre las características de los complejos con buckybowls. Se han considerado una serie de variables y efectos a estudiar tales como: tres tipos de buckybowls sustituidos con grupos nitrilo, uno basado en coranuleno y dos basados en sumaneno; cuatro líneas de aproximación del anión con respecto al bowl; distintas orientaciones de los iones con respecto al bowl; seis aniones diferentes agrupados en pares de distinto tipo: monoatómico (Cl-, Br-), trigonal plano (NO3-, CO2H-) y tetraédrico (BF4-, ClO4-). Las energías de interacción se han obtenido empleando los métodos SCS-MP2 y M06-2X, tal como sugerían los resultados del capítulo 8. Los resultados indican que en todos los casos los complejos se forman por la cara cóncava del bowl, siguiendo su eje central. Igualmente, las distintas orientaciones del anión no parecen tener demasiado efecto sobre la interacción, aunque se favorecen orientaciones con el mayor número de regiones negativas orientadas hacia la pared interna del bowl. Los complejos más estables en fase gas son los formados con CO2H-, seguidos por los monoatómicos, NO3-, y finalmente los tetraédricos como los menos estables. Este comportamiento cambia radicalmente al incorporar el efecto del disolvente. Todos los complejos sufren importantes pérdidas de estabilidad, pero los más afectados son los formados con CO2H- que pasan a estar entre los menos estables en disolventes con constantes dieléctricas elevadas. Los complejos más estables en esas condiciones son los formados con Br-. Los resultados conjuntos de estos tres trabajos indican que sería plausible que los aniones fueran atrapados por los buckybowls, dando lugar a complejos que podrían ser estables incluso en disolución, abriendo nuevas posibilidades para el diseño de receptores aniónicos basados en este tipo de especies tan atractivas. 1. Introduction 1. Introduction 3 1.1. Supramolecular Chemistry Supramolecular chemistry is a field of science that deals with the relationship between molecular structure and function; it is the chemistry of the noncovalent bond, which forms the basis of highly specific recognition, transport, and regulation events that act in biological processes. Supramolecular chemistry arises from the natural result of the curiosity of humans, who have imitated phenomena of the nature during centuries. This subject was initiated with the pioneering works of Jean-Marie Lehn in 1969 about the idea of molecular recognition and led him to obtain the Nobel Prize in 1987 together with Donald J. Cram and Charles J. Pedersen.[1, 2] A tentative definition of Supramolecular Chemistry provided by Lehn is the following:[1,3] “Supramolecular chemistry directs to the development of big complex chemical systems from the components that interact between them by means of intermolecular noncovalent forces” and has opened a new field where the interest of various disciplines of Chemistry converge, including Organic and Inorganic Chemistry, Physical Chemistry, Biochemistry, Materials Science and Nanotechnology. As mentioned above, the objective of Supramolecular Chemistry is to form an aggregate or molecular entity constituted by species that are not linked covalently among them, but associated by their geometric or/and electronic affinity; that is, they are molecularly recognized. Three concepts are essential to understand these processes: fixation, recognition and coordination. These concepts were established at the end of XIX century by Paul Ehrlich (1898), Emil Fischer (1894) and Alfred Werner (1893), respectively.[1, 2] Ehrlich recognized that a molecule cannot act if it does not bind to the neighborhood of others. On the other hand, Fischer introduced the key-lock concept to explain the performance of enzimes. According to this concept, the enzyme selectively recognizes the substrate because it presents a specific geometry in its active site, the same way a key fits with its lock. The better they fit together the more efficient the complexation will be. This recognition is not only geometric; but it also implicates a chemical interaction and a conformational rearrangement thus explaining the notable specificity of enzyme catalysis. This kind of noncovalent interaction can be related to the idea of coordination as introduced by Werner. Following the same line, Cram was the first one to introduce the terms host-guest.[4, 5] Thus, one supramolecule is obtained from the process in which one molecule acts as receptor (host) and another one as substrate (guest), binding to the first one to obtain a receptor-substrate complex. Normally, the receptor is a big molecule or aggregate, as an Alba Campo Cacharrón 10 Figure 1.2. Aromatic groups in amino acid side chains. When aromatic species are present in complex systems, it is usual for the structure to show π···π contacts, XH···π interactions and interactions between ions and the π cloud, being any of these contacts relevant for determining the structural characteristics and energy of the system. Following there is a brief description of these interactions.[13] X-H···π interactions This interaction is based on the attraction between X-H groups, normally C-H, N-H or O-H groups, and the π electron cloud of an aromatic ring.[13, 17-21] Normally, the range of distances between the hydrogen atom and the centroid of the ring is 2.4-3.2 Å, larger than the distance in a typical hydrogen bond (around 2 Å). The C-H···π interaction was first postulated by Tamres in 1952, who noted that dissolving benzene in chloroform was exotermic.[22] Qualitatively, the strength of the C- H···π interaction arises mainly from charge transfer, giving rise to interaction geometries where the CH bond lies directly in line with a p-orbital on the ring. O-H···π and N-H···π interactions show similar characteristics and have been frequently recognized as a motif contributing to the stability in different systems. Besides, despite the weakness of these individual contacts, their effects can be additive and in macromolecular systems their influence can be pronounced. This has been show in many cases by Nishio, who has compiled an extensive literature database on C-H··· π interactions and their effects on energetics, reactivity and conformational changes of different species.[17, 23] Theoretical studies and experimental results in the gas phase seem to indicate that the nature of C-H···π interaction is totally different to conventional hydrogen bonds. While the hydrogen bonds are mainly due to electrostatic interactions, the electrostatic component in the C-H···π interactions is minimum, being dispersion interactions the main cause of formation. Other fundamental difference is the directionality that both interactions exhibit. While the fundamental hydrogen bond feature is its high Phenylalanine (benzene) Hystidine (imidazole) Tryptophan (indole) Tyrosine (phenol) 1. Introduction 11 directionality due to the electrostatic contribution, the C-H···π interaction does not depend so dramatically on the orientation of the interacting group. Thus, it is frequent to find X-H···π contacts clearly departing from the typical linear arrangement of hydrogen bonds.[17, 19] π···π interactions (stacking) This type of interaction takes place between the π electron clouds of stacked aromatic systems[11-13, 24, 25] Normally, one of the rings is electron-rich, while the other is electron poor, although cases have been described where both rings have the same electron wealth. This corresponds in principle to a weak interaction (bonding energy around 0-10 kcal/mol) but with a global effect of great importance from a biological and supramolecular points of view, as well as in crystallography. π···π interactions, like in benzene dimer, are usually governed by dispersion effects, since these are the unique forces acting in the case of nonpolar molecules (the quadrupole-quadrupole interaction between benzene molecules is weaker). π···π interactions play an important role in the stabilization of DNA, together with hydrogen bonds, resulting in the stacking of bases pairs and generating its characteristics helicoidal structure.[11, 12] Based on this, a number of intercalating drugs have been designed exploiting π···π stacking interactions. On the other hand, π···π stacking interactions have a lot of applications in Supramolecular Chemistry, especially for host-guest systems. A remarkable case is the one described by Sygula et al., who have synthesized a buckycatcher based on a multitude of conjugated aromatic rings that adopt a concave conformation matching perfectly with a fullerene C60 molecule, acting as receptor of this molecule through π···π stacking interactions.[26] There has been controversy concerning the physical nature of this kind of interaction. In 1990, Hunter and Sanders proposed a simple model based on the competition between the electrostatic and the van der Waals forces to explain the variety of geometries observed for these interactions and to quantitatively predict their interaction energies.[27] These authors postulate the existence of an attractive overlap in van der Waals interactions, which are proportional to the surface contact area between both π systems. This overlap is due to an attraction between the π electronic cloud negatively charged of one of the rings and the  electronic cloud positively charged from the other. The relative orientation of both rings is determined by electronic repulsions between both π systems negatively charged. For this reason, when an aromatic system is stacked in a parallel way, it is normally observed that the rings are not totally aligned, with one slightly displaced with respect to the other, minimizing the π···π repulsion and maximizing the π··· attraction. In fact, few examples exist where the aromatic rings are arranged totally overlapping.[28, 29] The most common arrangements are the paralleldisplaced and T-shaped dispositions, because π··· attractions predominate in both. The three dispositions have been usually employed in this kind of complexes to study the Alba Campo Cacharrón 12 effect of the arrangement on the interaction energy and stability of the complexes. These arrangements are shown in Figure 1.4 with a benzene dimer as example.[24, 25] It can be observed that in the case of benzene dimer (the most studied π···π interaction), the stacking interactions are competitive with other possible arrangements as the C-H···π interaction in T-shaped structures.[24, 25] This can be observed in Figure 1.4, which shows how both structures are virtually isoenergetic. It is worth mentioning that the strength of these interactions is highly affected by different factors, such as the presence of electron-donor or electron-acceptor substituents in the rings, the existence of heteroatoms forming part of the ring and the degree of annelation of the rings involved in the stacking. Thereby, the higher the number of fused rings, the more favorable the stacking.[13] Regarding substitution effects, it has been shown that the presence of electron-attractor substituents in the ring increases the strength of this kind of interactions since the electronic density of the π cloud of the ring decreases, minimizing the π···π repulsions between the rings. In fact, studies by Sherrill and coworkers have shown that substitution always leads to stronger interaction unregardingly of the donor or acceptor character of the substituents.[30, 31] Figure 1.4. Types of π···π contacts for benzene dimer with their corresponding interaction energies. When the stacking is produced between aromatic heterocycles the strength of these interactions increase. The substitution of a carbon atoms of the ring by a nitrogen atom in a six-membered rings (like in pyridine for example) causes a decrease of the electronic π density in the carbon atoms of the ring, which leads to a stabilization of the system and decreases the π···π repulsion forces as mentioned above.[32] In any case, this interpretation of substituent effect has been recently reviewed by Wheeler.[33, 34] 1. Introduction 13 Cation···π interactions A cation···π interaction arises from the electrostatic interaction of a cation with a face of a π system, as could be a benzene ring and its derivatives or π systems as ethylene. The first evidence of this type of interaction in the gas phase came out from the work by Kebarle in 1981.[35] In a systematic study of ion solvation by various solvents, Kebarle observed that benzene stabilizes K+ ions better than water in the gas phase. Kebarle proposed that this was a result of the ion interacting with the quadrupole moment of benzene. A molecule with a dipole moment, such as water, experiments a favorable electrostatic interaction with an ion if the ion is positioned near the appropriate end of the dipole. Benzene, of course, has no dipole moment, but it does have a substantial permanent quadrupole moment.[36] A quadrupole can be thought of as two dipoles aligned in such a way so that there is no net dipole. Thus, there is a permanent nonspherical charge distribution in benzene, with regions of relative negative and positive charges. Just as an ion can be attracted to the appropriate end of a dipole, so can an ion experience a favorable interaction with appropriate regions of a quadrupole. This is an electrostatic interaction and does not requires adjustment of the electronic distribution around the ion or the molecule. Importantly, there is no a priori reason to expect that such interactions will be inherently weaker when the molecule contributes a quadrupole rather than a dipole as demonstrated by Reisse and Williams.[37, 38] Most neurotransmitters have a cationic group that permits them a selective anchorage to their receptors by cation···π interactions.[39-41] The cation···π interaction has important applications in the field of Supramolecular Chemistry. [41-49] Numerous studies have reported the occurrence of cation···π interactions in protein structures and in proteinligand and protein-DNA complexes.[40, 41, 50, 51] These analyses have revealed the preferential location of amino groups in the area of aromatic rings.[39, 52] This interaction is calculated to be even more stabilizing than an analogous salt bridge, and it is not so strongly attenuated in water.[53] The side-chains of the aromatic amino acid residues, Phe, Tyr and Trp, provide a surface of negative electrostatic potential than can bind to a wide range of cations through a predominantly electrostatic interaction. It has been found that 50% of the Arg residues are in contact with an average of two aromatic side chains. Of particular interest is the interaction of the cationic Arg residue with aromatic side chains. Two limiting geometries are possible, a perpendicular arrangement in which the NH of the Arg points into the face of the aromatic unit, and a parallel or stacked arrangement of the planar guanidinium of Arg and the aromatic moiety. The stacked arrangement is more frequently found, but there seems that this is related to environment effects. Also, series of cation···π interactions involving both Arg and Lys appear in different structures containing several aromatic and cationic side chains from different strands of the protein, intercalated to form an extended array of cation···π interactions.[39, 51, 52] Another remarkable case is the acetylcholine nicotinamide receptor whose mechanism of molecular recognition to their substrate is based only on cation···π Alba Campo Cacharrón 14 interactions.[40, 54] Furthermore, systems with molecular structures similar to crown ethers, with π systems strategically placed, have shown to be very effective binding places for alkali cations.[41, 55] Additionally, cation···π interactions have also been used to increase the π-face selectively in catalysis in asymmetric catalysis.[41, 56] The strength of the cation···π interaction has been rationalized on the basis of the strong electrostatic interaction between the positively charged cation and the ring negative molecular electrostatic potential (MEP). Many studies, especially of gas-phase complexes, established that electrostatic interactions play a prominent and usually dominant role in prototypical cation···π interactions.[41, 57, 58] Electrostatic reasoning can also explain variations due to changes in the aromatic ring. The more negative the maximum in electrostatic potential over the center of the aromatic molecule, the stronger the cation···π interaction. In fact, good correlation has been found in many cases between the MEP value and the strength of the cation···π interaction. [41, 59, 60] A clear indication that electrostatics play an important role in cation···π interactions comes from a comparison of simple alkali metals binding to benzene. The observed trend in stability is Li+ > Na+ > K+ > Rb+, the classical electrostatic sequence, and exactly what would be seen if benzene were replaced by Cl- or if one was comparing hydration energies.[41, 50, 61] However, this analysis neglects the role played by polarization interaction. The electrostatic contribution does not represent the 100% of the energy of cation···π interactions; in fact, the fraction of the total binding energy that comes from electrostatics varies considerably depending on the aromatic molecule. The “nonelectrostatic” component of the cation···π interaction, sometimes the major component, reflects a combination of effects mostly related to the polarizability of the aromatic unit. Probably the most important of these for simple systems is the interaction of the ion with the induced dipole in the aromatic molecule.[62, 63] As the aromatic species becomes larger it is expected that induction contribution also becomes larger reflecting the increase in polarizability. In fact, studies show that there can be cation···π complexes held by induction contributions, whereas the electrostatic term is repulsive.[64] Induction contributions also explain the presence of off-plane cation···π stabilizing interactions for example in benzene···Na+ complexes. The location of the cation in the ring plane is destabilized electrostatically, but also shows a stabilizing induction contribution leading to a global attractive (though weak) interaction.[65] Different studies have been devoted to determine the origins of substitution effects upon the cation···π interaction. Recent studies propose that the origin of the changes has nothing to do with changes in the π cloud but depend on trough-space interaction between the cation and the substituents.[59] These results have been slightly corrected by Quiñonero et al. indicating that induction effects are also responsible of such changes.[66] 1. Introduction 15 In summary, the cation···π interaction is an intense interaction in the gas phase that can frequently appear in systems of interest and basically controlled by electrostatic and induction contributions. However, the presence of solvent or other units near the cation···π contact can modulate significantly the strength of the interactions. It is usually observed a decrease on the interaction strength, but this can also be accompanied by significant structural changes in the geometry of the complex. Solvent molecules compete with the aromatic component for interacting with the cation and, as a consequence, a decreasing of the intensity of the cation···π interaction in the presence of the solvent is normally observed. This led to some controversy regarding the role of cation···π interactions in protein stability. While some studies suggest a relevant contribution another studies estimate that the contribution is insignificant.[15, 41, 53, 58, 61, 67] These differences are normally attributed to the different degree of exposure of the cation···π contact to the solvent. If the contact is buried in a hydrophobic region, it can present an appreciable intensity, while if the contact is exposed to the solvent, it even cannot take place. Anion···π interactions While coordination of cations has been object of study one century ago, the coordination of anions has received little attention until short time ago. The advent of synthetic molecules able to coordinate cations and anions was almost simultaneous: in 1967, C. J. Pedersen prepared the first synthetic ligand able to coordinate cations,[68] and only one year after C. H. Park and H. E. Simmons synthesized the first system suitable to coordinate anions which was called “kapatinato” (from kapatinosis, which means swallow in greek).[69] Despite of this almost simultaneous discovery, the area of anion coordination was relatively unexplored in contrast to the cation coordination. The first discovery of the concept anion···π interaction in literature is dated from the year 2000 by Schneider.[2] From that moment on, various studies were carried out to understand the nature of the interaction between anions and aromatic species.[70-81] The results of those studies demonstrated that, contrary as it could be expected, the intensities of interactions between cations and anions with aromatic species can be of similar magnitude as shown in Table 1.5. A priori the idea of an interaction between an aromatic ring and an anion will be nonviable due to the electron donor capacity of both molecules, but the electron donor capacity of the aromatic ring can be modulated so, if the aromatic ring is completely substituted with groups which attract electrons, it becomes an electro-deficient aromatic cloud and then an attractive interaction with an anion is possible. Thus, anion···π interactions are described as the favorable noncovalent interactions between electron deficient aromatic systems and an anion.[70, 72] Alba Campo Cacharrón 16 Table 1.5. Interaction energies of different types of contacts with the π cloud obtained at diverse levels of calculation in kcal/mol. ΔE (kcal/mol) C6H6···C6H6 -2,5 C10H8··· C10H8 -5,7 C6H6···H2O -3,0 C6H5OH···CH3CN -6,6 C6H6···Na+ -21,3 C6H6···K+ -17,0 C6F6···Cl- -14,0 C3N3H3···Cl- -6,9 Studies revealed that anion···π interactions are, in general, dominated by the electrostatic and induction contributions.[75, 78] The electrostatic component of the interaction is directly related with the permanent quadrupole moment of the electrondeficient aromatic ring. As indicated in Figure 1.6, the quadrupole of a benzene ring is positive but it can be modulated through the substitution of groups on the rim of the ring. For instance, the quadrupole moment for a benzene ring is -8.45 B while the quadrupole moment for the benzene ring substituted with six fluoride groups is positive by 9.50 B.[70, 72] Therefore, the presence of electron-withdrawing substituents (halogen, nitro or cyano groups) or nitrogen atoms in the ring (pyridine, triazine, tetrazine, among others) favors the formation of anion···π complexes. Figure 1.6. Schematic representation of the quadrupole moments of hexafluorobenzene and benzene. Quadrupole moments (Qzz) in Buckingham and molecular polarizabilities parallel to the main symmetry axis (||) in atomic units are given. 1. Introduction 17 A topological analysis of the electrostatic potential in anion···π interactions has demonstrated that there exists a correlation between the MEP value of the aromatic ring and the electrostatic contribution to the anion···π interaction.[70, 72] Thus, the systems with a very positive quadrupole moment give more favorable interactions. Other studies show that for molecules with very positive quadrupole moments the anion···π interaction is dominated basically by the electrostatic term, while for molecules with small quadrupole moments the polarization contribution induced by the anion can be the dominating one. Finally, considering the larger polarizability of anions it has also been found that dispersion contributions tend to be larger in anion···π complexes than in similar cation···π ones.[70, 72] The properties of the anion are also important for applications of anion···π interactions in Supramolecular Chemistry.[70, 72, 82] Both the electrostatic and polarization contributions to the total interaction energy depend strongly upon the ion-arene distance. Small anions are more polarizing and show short equilibrium distances and, consequently, give rise to stronger interactions. In addition, planar and linear anions such as NO3- or N3- can also interact with the aromatic ring via π···π stacking.[70, 72, 79] Finally, a comment must be said on the interplay of the different interactions commented above, all of them involving aromatic species. Different π interactions are very important and omnipresent in a great variety of biological systems, so the study of the mutual influence of different interactions is crucial. The interplay between ion···π and π···π interactions, which can lead to strong cooperativity effects, has been shown. The cooperativity effects can be favorable or unfavorable depending on the nature of the aromatic ring and the sign of the ion.[41, 83-87] The theoretical results on ion···π···π complexes have also been used to explain an unexpected experimental finding regarding the parallel stacking of pentafluorphenyl groups in substituted ferrocenes.[88] Most often, the interactions are not isolated but inserted in long structures where multiple interactions of different type are possible. Thus, it is quite typical to find π···π interactions at the same time as cation···π or X-H···π interactions in many proteins and amino acids. Sometimes, these interactions are combined forming a larger complex where two or more different kinds of interactions are established with simultaneous participation of the aromatic units.[11, 12, 41] 1.4. Interactions with buckybowls Polycyclic aromatic hydrocarbons (PAH) are a family of hydrocarbon molecules that typically possess a structure formed by a series of fused benzene rings creating a planar structure. In graphene, material with potentialities still to discover, this structural motive extends creating large bidimensional layers. Some of these PAHs constitute important atmospheric contaminants with some implications for human health but, on the other hand, some PAHs have been located in the interstellar medium and postulated as species that could act as basis for the more primitive ways of life.[89] Alba Campo Cacharrón 18 In 1966 Barth and Lawton presented the first synthesis of the C20H10 PAH, called corannulene, that presents a non-planar structure.[90] Contrary to other PAH as helicenes the non-planarity of corannulene is not a consequence of the presence of bulky groups that force the location of part of the fused rings out of the plane. The origin of the nonplanarity of corannulene is the proper structure of the carbon rings: corannulene is formed by five benzene rings grouped around a pentagonal central ring. The presence of the five carbon central ring introduces the curvature in the corannulene structure because it is not possible to build a planar layer fusing pentagons and hexagons. According to this, the structure is deformed creating species with bowl shape that, contrary to planar PAH, present two faces, concave and convex, which can exhibit different properties. This is the motivation why these species are denominated molecular bowls, π-bowls or buckybowls.[91] The relatively recent discovery of buckminsterfullerene and other fullerenes has incremented the interest in buckybowls because they are considered as fullerene fragments and they could be a key piece in the synthesis and design of new species related to fullerenes and carbon nanotubes.[92, 93] There is a complete family of buckybowls that has been identified or synthesized in the laboratory, and it has been observed that though these species frequently share as structural motive their bowl shape, their properties change appreciably with size, curvature and substituents.[91] Figure 1.7. Buckyminsterfullerene and the smallest buckybowls related to it. Corannulene C20H10 Sumanene C21H12 Buckyminsterfullerene C60 1. Introduction 19 In any case, this area of study is recently emerging: while corannulene has been synthesized in 1966, the other smallest buckybowl related to C60, denominated sumanene (C21H12), has been synthesized so recently as in 2003.[94] Sumanene presents a structure that is formed by a hexagon as the central ring and alternating hexagons and pentagons around it. Due to its recent discovery, relatively few studies which describe properties and characteristics of sumanene are published.[94-99] These carbon aromatic nanosystems show subtle dependences between their structure, dynamic and properties, so they are a set of materials with great potentialities in a variety of areas. One of the most interesting aspects is the possibility of creating intermolecular complexes of different nature employing one or both faces of the buckybowls, creating supramolecular structures bonded by noncovalent interactions. In this sense, there is a great variety of possibilities in which buckybowls could act as receptors interacting with other buckybowls or fullerenes, transition metals, ions of different nature, etc.[91] An interesting application of these species is related to a concave-convex interaction between two buckybowls by means of π···π contacts, so they that can be employed as tweezers for catching species of interest, for instance fullerenes. In this case, buckybowls can be used as one of the most efficient buckycatchers, as proposed by Sygula.[26] Whereas stacking with planar hydrocarbons is limited to the region nearby the bottom of the bowl, the use of curved surfaces allows extending the contact area increasing the interaction. Different modifications of coranulenne and sumanene have been studied as possible fullerene receptors based on π···π interactions, though other effects as C-H···π can also help stabilizing the complex as already shown for sumanene.[100-106] Figure 1.8. Tweezers developed by Sygula et al. to catch fullerenes. 2. Objectives 27 The main objective of this thesis is to gain more insight into the characteristics of systems that establish interactions between aromatic and charged species by applying computational chemistry methods. The interest on ion···π interactions has grown in the last decades with new evidences of their importance in many fields covering from biochemistry to materials science. More specifically, the upcoming of anion····π interactions has produced a renewed interest about the characteristics and nature of the interactions with aromatic systems, as well as their relationship with other types of interactions. Computational chemistry has also recently evolved with new methods and functionals specifically designed for studying intermolecular interactions, especially regarding the treatment of dispersion interactions. These new approaches, together with techniques allowing significant reductions of the computational cost, have made possible the application of appropriate methods to systems of increasing size, thus avoiding being restricted to the simplest cases as benzene complexes with alkaline or halogen ions. The study of ion···π interactions constitutes a wide field for research, including very different systems and phenomena, so in this thesis the goal will be focused on a series of specific aspects of the cation···π (chapters 4 and 5) and anion···π (chapters 6 to 9) interactions. More precisely, the goals of the different chapters are summarized in the following. Chapter 4 is devoted to phenol···cation (K+, Na+, Li+ and Mg2+) complexes. Though the interaction of simple aromatic molecules, mostly benzene, with simple cations as alkaline ones has already been studied in detail, there are several aspects of the interaction that deserve more attention. Cation···π interactions are strong in the gas phase, but the presence of solvent molecules or other donating groups nearby can modulate the strength and characteristics of the interaction. Microhydration of complexes opens a route for isolating the effects due to specific interactions with the solvent. Clusters formed by phenol and simple cations will be subjected to stepwise microhydration in order to determine the most stable arrangements of the hydrated complexes. That way, the balance among the different attractive contacts in the clusters can be analyzed and their impact on cluster properties determined. Moreover, systems containing a small number of molecules can be monitored by infrared spectroscopy, so the vibrational spectra will be predicted for the complexes studied. The results will be compared to the experimental ones obtained by Vaden and Lisy (J. Chem. Phys. 2004, 120, 721-730) for phenol complexes with K+ and Na+ and four water molecules. Chapter 5 deals with the analysis of another possible modulating effect in cation···π complexes. Besides solvent molecules, other electron donating groups can also interact with the cation in a cation···π contact. In this chapter, the effect of a second aromatic Alba Campo Cacharrón 28 unit interacting with a cation···π contact will be analyzed. Ternary systems formed by two aromatic molecules and one cation will be studied in order to quantify the balance between the different cation···π, X-H···π and π···π interactions and their mutual influence. The systems are formed by a cation and aromatic rings selected as to represent possible contacts between amino acid side chains. Benzene, phenol and indole are employed as aromatic units, whereas guanidinium cation is selected as a model of the cationic side chain of arginine. The most stable structures of these ternary complexes will be determined and their characteristics analyzed. Consequently, possible specific interactions, as hydrogen bonds in phenol and indole complexes, or the role played by dispersion interactions would be revealed. Also, a pair energy decomposition will allow determining the magnitude of possible cooperative effects in this kind of complexes. Chapter 6 is the first one devoted to anion···π interactions. Recently, Matile et al. (J. Am. Chem. Soc. 2006, 128, 14788-14789) synthesized the first anionic channel which is supposed to be based on anion···π interactions. This anion channel is constructed by a series of oligomers of naphthalendiimides (NDIs), which arrange themselves as to create a pore through which anion transportation has been verified. However, several questions remain unanswered about how this channel works. Specifically, it has been suggested that solvent molecules can play an active role in anion transport, maybe going into the channel and facilitating the process. Two simple models will be employed for studying the interaction of anions with NDI, and the effect of water molecules will be estimated by explicitly incorporating a small number of water molecules into the model and also by using a continuous representation of the solvent. The results will help understanding the role of the solvent molecules closest to the anion in order to favor the interaction with the channel units. Chapters 7 to 9 are dedicated to an emergent field in chemistry related to the properties and interactions of curved aromatic systems (buckybowls). In these chapters, the possible use of these buckybowls as ion receptors will be analyzed by introducing different modifications in the structures of the bowls. Substituted derivatives of corannulene and sumanene will be employed in these studies. Chapter 7 is our first approach to this subject. Corannulene shows negative molecular electrostatic potential (MEP) by both faces of the bowl, but proper substitution with electron-withdrawing groups can invert its molecular electrostatic potential and provide a favorable interaction with anions. Substitution of corannulene with five and ten F, Cl and CN groups will be considered, and its effect upon the properties of the bowls will be analyzed, specially regarding their MEPs, seeking a stronger interaction with anions. The effect of these substitutions upon the interaction with anions will be studied in 2. Objectives 29 complexes formed with Cl-, Br- and BF4- anions, obtaining their optimal geometries and interaction strengths. It is intended to form remarkably stable complexes, looking for inclusion structures with the anion by the concave side of the bowl, as opposed to the convex complexation already observed in cations. Chapter 8 is more oriented to the performance of different computational methods. It has been observed that in this kind of systems a reliable method is difficult to find, the errors being quite large in many cases. Since there are no reference values, the method employed is somewhat blindly chosen. In order to alleviate this problem, a thorough study will be carried out in complexes of Na+ and Clwith buckybowls. Two typical buckybowls, corannulene and sumanene, will be considered and substituted by CH3, F and CN groups to promote changes in the molecular electrostatic potentials. Consequently, the interaction will be modulated favoring cations or anions. Reference values for the interaction will be obtained by employing high-level calculations, being used afterwards for checking the performance of a variety of more affordable methods. A reliable account on how substitution affects the interaction in both anion and cation complexes will be obtained, allowing the direct comparison of cation···π and anion···π interactions in the same set of systems. The detailed analysis of the characteristics of the interaction will hopefully reveal the intrinsic different nature of the interaction with cations and anions in these extended π systems. Chapter 9 completes the studies on the interaction between anions and buckybowls. After selecting an appropriate method in chapter 8, now this method will be applied in order to gain insight about how the shape, size and disposition of the anions affect the properties of the complexes. CN-substituted corannulene and sumanene will be employed to study the interaction with six different anions, going from monoatomic Cl- and Br-, to planar trigonal NO3- and CO2H-, and to tetrahedral BF4- and ClO4-. Several orientations of the polyatomic anions as well as different approaching lines to the bowl will be considered. The characteristics of these complexes will be analyzed in detail in order to find the key factors controlling the interaction in each case and the possible differences associated to the different anion’s nature. Also, solvent effects as modeled by a continuum model will be evaluated to check whether the proposed complexes would be formed in solution. A series of solvents with different dielectric constant will be used to quantify how the gas-phase interaction is altered and to which extent anion’s nature favors or hinders the formation of complexes in solution. 3. Methodology 3. Methodology 33 3.1. Interaction Energy Though interaction energies are several orders of magnitude smaller than electronic energies, and usually significantly smaller than bond energies, the usual approximate methods employed for solving the Schrödinger equation in chemical systems can be employed for studying noncovalent interactions.[1-3] The concept of interaction energy appears naturally within the Born-Oppenheimer approximation when a system formed by say, two molecules, atoms or ions is considered.[3] Quantum chemistry treats the system as a whole so after calculating with a given method, the energy of the complete system (supermolecule) is obtained. As indicated in Figure 3.1, a system defined by a set of nuclei at given positions together with the corresponding electrons has to be divided (arbitrarily or not) in a set of separable atoms, molecules or ions. Figure 3.1. A phenol···water dimer. Left: the set of atoms forming the system. Right: the division employed for studying the interaction between phenol and water. Once this separation is done, the energy of the system can be expressed as:[3-5] int ABBAAB EEEE  . (eq. 3.1) The energy of the system is then expressed as a sum of the energies of the isolated fragments plus a contribution coming from their interaction. Following this line, the interaction energy could then be obtained within the supermolecule approach as: )()()()( int  BAABAB EEEE , (eq. 3.2) where it has been highlighted that the same set of nuclear positions has been employed in the calculation of the three quantities.[3, 6-8] However, if we are interested in studying the process of complex formation, a new term has to be included, since the geometries of the species forming the complex can change as a consequence of the interaction and formation of the complex.[6-8] Therefore, a new term has to be included describing this effect, Alba Campo Cacharrón 34 )()()()()( 00  BABAdef EEEEE . (eq. 3.3) Along this work this term will be called the deformation energy,[9, 10] though in literature other names can be found for this quantity, such as relaxation energy (in the sense that if we think of dissociation, the geometry has to relax to that of the isolated fragments)[6-8] or preparation energy (in the sense that this is the energy needed to prepare the fragments in order to form the complex at the final geometry).[4] The sum of these two quantities defines the complexation energy, the energy change observed when a complex is formed from the isolated, relaxed molecules that form it. The negative of this quantity is often called the binding energy, related to the energy needed for separating the complex into isolated fragments. Other effects such as zero point energy corrections or thermal effect are added as usual to each of the energies needed for obtaining the magnitudes described above. In summary, the complexation process can be formally divided into two steps: deformation plus interaction, so the complexation energy is obtained as: )()()( int  defAb complex AB EEE (eq. 3.4) or )()()()( )()()()( 00   BABA BAAB complex AB EEEE EEEE (eq. 3.5) and: )()()()( 00  BAAB complex AB EEEE (eq. 3.6) So, the complexation energy could be obtained without making any reference to the interaction and deformation energies. However, it could be interesting to separate the complexation process into these two contributions in order to obtain a better interpretation of the system, especially when large deformation effects are present as consequence of important geometry changes. In such situations, a very large interaction could be hidden by opposing large deformation effects.[7, 8, 11, 12] Also, the practical calculation of complexation energies faces some problems, which make advisable to separate them into these two contributions as it will be commented below. 3.1.1. Basis Set Superposition Error Any researcher dealing with the quantum chemistry calculation of interaction energies has faced the problem of Basis Set Superposition Error (BSSE).[3, 5, 6, 13] The problems appear when using (eq. 3.2) for obtaining the interaction energy. With this expression, the interaction energy of a dimer is obtained as a difference between three quantities. 3. Methodology 35 As commented above, interaction energies are several orders of magnitude smaller than electronic energies, so the interaction energy, a small quantity, is obtained from the difference among large quantities. Therefore, it has to be ensured that all electronic energies are obtained coherently, and that any errors coming from the method employed will be properly cancelled out. However, a practical quantum chemistry calculation is performed with a finite basis set employed for constructing the wavefunction, and this is the origin of the problem. Consider a dimer AB. When monomer A approaches monomer B in order to form the dimer AB, the energy of the dimer is artificially lowered because monomer A can employ the basis set centered onto the atoms of B in order to improve the description of its electron distribution, and the same applies for monomer B using basis functions centered in A. This energy lowering is not possible in the calculations in isolated monomers where only the basis set of each monomer is present. As indicated by van Duijneveldt,[14] this energy lowering as extra basis functions are employed is not an error in itself. The error comes from an inconsistent treatment of the monomers, which cannot take benefit from the basis set of the other as the distance becomes larger. This inconsistent treatment of the monomers as the distances change is the origin of the BSSE. Therefore, BSSE depends on the geometry of the systems, and a different value is obtained for each different arrangement of nuclei. It has to be taken into account that even if BSSE is removed completely, there still remain other errors associated to the method employed and the finite basis set used. The usual procedure for correcting BSSE is the counterpoise method proposed by Boys and Bernardi, explained below.[15, 16] The uncorrected interaction energy or a dimer AB would be obtained as in eq 3.2, or: )()()()( dimerdimerdimerint BEAEABEABEBAABAB  , (eq. 3.7) where now the terms in parentheses indicate the basis set employed in the calculations and superscripts the geometry employed. These quantities are obtained in three separate calculations for the dimer, monomer A with its basis functions and monomer B with its basis functions. As commented above, the different basis sets employed in the calculations introduce BSSE, so the counterpoise correction consists on obtaining the three energies using the same basis set and energy in all calculations. )()()()( dimerdimerdimerint ABEABEABEABEBAABAB  . (eq. 3.8) Now the energy of monomer A is obtained in a calculation with the basis set of B in the same positions as in the dimer, and the same applies to monomer B. In that way, the same basis set is employed in all calculations and BSSE is corrected. Therefore, BSSE can be defined as: )()()()( dimerdimerdimerdimer ABEABEBEAEBSSE BABA  (eq. 3.9) Alba Campo Cacharrón 42 000 ˆ   . (eq. 3.30) 0 ˆ  is the wave operator, that in CI is a linear operator, but in CC is an exponential one: ... ˆˆˆ 1 ˆ3210  CCC CI (eq. 3.31) ... ˆ 6 1 ˆ 2 1 ˆ 1) ˆ exp( ˆ32 0 TTTT CC (eq. 3.32) ... ˆˆˆˆ 321  TTTT . (eq. 3.33) More explicitly, the expansion of the exponential will be (only singles and doubles included): ... ˆˆ !2 1 ˆˆˆ !2 1 ˆ !4 1 ˆ !3 1 ˆ !2 1 ˆˆ 1 ˆ2 1221 2 2 4 1 3 1 2 1210  TTTTTTTTTT . (eq. 3.34) The C and T operators are excitation operators that change one, two, …. occupied spinorbitals by virtual ones. For example:     sr ba rs ab rs ab tT  02 ˆ , (eq 2. 35) where rs ab t are the coefficients (or amplitudes) to be determined in CC calculations. The great advantage of CC methods over CI ones is that the exponential antsazt employed in CC methods ensures that the method will be size consistent, thus avoiding the major problem faced by CI calculations. Usually, CC methods are employed truncated up to a given excitation order, typically including single and double excitations, CCSD. Though amplitudes are only calculated for singles and doubles, the exponential expansion introduces estimations of higher-order excitations expressed as products of singles and doubles. For example, quadruple excitations are not included, but their effect is estimated by products of doubles. Therefore, to a given truncation level CC behaves better than CI since these so-called disconnected clusters include the effect of higher excitations. CI methods, on the other hand, only include excitations up to the truncation level. The most common choice in present calculations is to employ CCSD method combined with a suitably flexible basis set. Including triple excitations implies a huge computational effort so CCSDT is only applicable to small systems. Most often, after the CCSD solution has been reached, a perturbative estimation of the triples is added in the so-called CCSD(T) method. That way the effect of triple excitations (often very important) is included without needing to solve the full CCSDT equations. In any case, 3. Methodology 43 inclusion of triples is also a demanding task, so CCSD(T) can only be used for systems with moderate size. Nowadays, results obtained at the CCSD(T) level with a large basis set are considered as the golden standard of quantum chemistry calculations, especially in the case of intermolecular interactions dominated by dispersion effects. Very recently, results obtained with CCSD(T) have been compared with those obtained with higher order CC calculations such as CCSDT[Q] for a set of complexes, confirming the very good behavior of CCSD(T).[54-56] 3.2.4. Errors in Wavefunction-based Methods From previous sections it is clear that a typical calculation with a wavefunction-based method is affected by two sources of error as depicted in Figure 3.2. When comparing the result of an actual calculation employing a post-HF method with the exact solution (within Born-Oppenheimer approximation and non-relativistic Hamiltonian) an apparent error is observed coming from two different sources.[21] Figure 3.2. Representation of the errors associated to electronic structure calculations. Alba Campo Cacharrón 44 Considering that FCI solution is the exact solution for a given basis set, an error is introduced due to the deficiencies of our method for describing correlation. MPn and CC methods are different approaches for approximating to the FCI solution, but in actual calculations both MPn and CC have to be truncated somehow, thus introducing an error as compared with the FCI result (the N-electron error in Figure 3.2). In any case, the methods currently employed for including electron correlation already show a very good performance, so the N-electron error can be reduced including higher orders of MPn or CC, of course at the cost of much more expensive calculations. The other source of error is related with the set of one-electron basis functions employed for constructing the N-electron wavefunction. That is, for a given N-electron model (MP2, MP4, CCSD, …) there is an error associated with the basis of one electron functions employed, the basis set. When the result of an actual calculation is compared with the same calculation employing a complete basis set, a difference is observed corresponding to the basis set error. In order to reduce this error, larger basis sets have to be employed, ideally reaching the limit of complete basis set. When this limit is reached, there still remains an error related with the deficiencies of our N-electron model as commented above. Therefore, in order to obtain accurate results from a wavefunction based method, one has to take care in order to reduce both sources of errors, thus trying to employ an adequate representation of correlation by means of a good N-electron model, combined with a one-electron basis set large enough in order to reduce the errors associated to incompleteness of the basis set. That is, one has to use the obvious choice of a good correlation method combined with a large basis set. The choice of N-electron model and basis set is of course conditioned by the computing resources since an increment in the quality of any of them has a deep impact on the computation time. The currently used N-electron models perform quite a good job approaching the FCI results, including correlation energy to a great extent. In practical calculations, however, one is normally reduced to the use of MP2 if the size of the system is moderate, or to CCSD(T) if the size of the system and resources allow it. As commented above, the use of CCSD(T) combined with a flexible basis set is nowadays considered as the gold standard in quantum chemistry calculations. As indicated in a recent paper, CCSD(T) to the basis limit is able of reproducing the interaction energies of a set of complexes within 1.5% error.[54] In summary, as for the N-electron model, we are limited in the practice to choose between the practical and cheapest MP2 and the more accurate and demanding CCSD(T). Other intermediate options could be of interest in particular problems. The error associated to the one-electron basis set deserves a little more attention, so it will be analyzed in the following section. 3. Methodology 45 3.2.4.1. Extrapolating to basis limit The problem of the error associated to the one-electron basis set is that it converges very slowly as the size of the basis set is increased, so very large basis sets are needed to approach the limiting value of complete basis set (CBS). Figure 3.3 represents the convergence of the correlation energy of a water molecule as the size of the basis set is increased. It can be observed that the convergence is very slow, and even with the enormous cc-pV6Z basis set the result is still far from the limit. This bad convergence with respect to the one-electron basis set is related with the inability of the usual oneelectron basis function to reproduce the electron-electron cusp.[21, 57, 58] In any case, in practical calculations the largest basis sets to be employed as routine are cc-pVQZ or aug-cc-pVTZ or similar, so it can be observed from Figure 3.3 that the results will be far from the CBS limit. It has to be taken into account that when computing energy differences as for example in obtaining interaction energies part of these errors can be cancelled and the results could be nearer the limiting value. Figure 3.3. Variation of the correlation energy in water molecule as the size of the basis set is increased. X 0 2 4 6 8 10 Ecorr (a.u.) -0.40 -0.38 -0.36 -0.34 -0.32 -0.30 -0.28 -0.26 -0.24 -0.22 -0.20 water molecule cc-pXDZ Alba Campo Cacharrón 46 The solutions to this problem depend on using very large basis sets, which is impractical or, more recently, by using explicitly correlated methods, which tend more quickly to the limit than the standard ones.[57, 58] However, there is an intermediate solution used quite often which passes by extrapolating the results obtained with moderate-sized basis set in order to obtain an estimation of the limiting value.[59-65] These extrapolation schemes are only applicable if the behavior of the energy as the basis set grows is smooth so it can be fitted to a function and then obtain the limiting value. Thus, this kind of extrapolation schemes are mostly limited to well-balanced basis set as the correlation consistent cc-pVXZ family proposed by Dunning.[66] In this respect, one should distinguish between the behavior of HF energies and the contributions coming from correlation. In the case of HF energies, it has been found that in atoms the energies approximately follow an exponential behavior.[20, 57, 67] Assuming that this exponential decay can be also employed in molecules, an often-employed extrapolation scheme assumes that the HF energy behaves as: )exp( BXAEE HF CBS HF X . (eq. 3.36) So, if three calculations are performed with correlation consistent basis set of increasing X, the limiting value could be estimated as: HF X HF X HF X HF X HF X HF X HF CBS EE EE Bb b bEE E 21 11 )exp(; 1        . (eq. 3.37) However, the need for three calculation of increasing size is uncomfortable, so other two-point extrapolation schemes have been devised. For example, Karton and Martin proposed the following two-point procedure:[68, 69] )9exp()1( XXAEE HF CBS HF X . (eq. 3.38) In any case, HF energies converge quite quickly with the basis set, and especially when energy differences are considered, as interaction energies. Thus, it is usually found that obtaining the HF interaction energies with a basis set of cc-pVTZ quality is often enough since the errors associated to correlation are normally larger. When dealing with correlation energies, the dependency on the basis set size is still more important, so a good estimation of correlation energies often demands very large basis sets. Furthermore, the convergence of correlation energy with the basis set size is even slower than in the HF case, so larger basis sets would be needed, as represented in Figure 3.3. Fortunately, there are extrapolation schemes that seem to perform reasonably well allowing a good estimation of the correlation energy at a moderate cost. In He atom it has been observed that the error in correlation energy varies as 3. Methodology 47 3  cNEN , being N the principal quantum number.[57] That is, there is a cubic decay of the correlation energy as larger basis sets are included. Assuming a similar behavior in molecules and identifying N with the X ordinal in Dunning basis sets,[59, 60, 63] it can be assumed that the correlation energy changes as: 3  AXEE XCBS . (eq. 3.39) This expression contains only two unknowns and therefore can be solved if two calculations are performed for obtaining the correlation energy, 3  AXEE XCBS , (eq. 3.40) 3  AYEE YCBS . (eq. 3.41) So: 33 33 YX EYEX EYX exact    . (eq. 3.42) Of course the result will depend on the values used for X and Y, but it has been found that a T-Q extrapolation employing cc-pVXZ basis sets already gives results better than the direct calculation employing a cc-pV6Z basis set. In practical calculations in systems of moderate size, one is usually limited to perform T-Q extrapolations with the cc-pVXZ basis set or D-T ones with the aug-cc-pVXZ basis sets. 3.2.4.2. Obtaining benchmarking values In the previous section it has been exposed how applying an extrapolation scheme to correlation energies there is an affordable route for obtaining correlation values at the CBS limit. However, this kind of extrapolation is still quite demanding so it is usually performed with MP2 estimations of the correlation energy. Therefore, even when the MP2 values are obtained at the basis limit, there still remains an error associated to the low quality of the N-electron model employed. The straightforward solution would be employing a better model, say CCSD(T), in order to obtain better estimations of the correlation energy, and then apply an extrapolation scheme. This option, though formally appropriate, is very demanding, and for systems of moderate size cannot be applied in routine calculations due to the high cost of the CCSD(T) calculations with the larger basis set. Therefore, other approaches have been devised in order to estimate the CCSD(T)/CBS values but with a reduced computational cost, many of them proposed by Hobza. Hobza and collaborators observed that even though the correlation energy contributions to the interaction energy of a dimer are very different with MP2 and CCSD(T), the differences between both methods were pretty independent of the basis set size.[52, 70, 71] Therefore, the Alba Campo Cacharrón 48 CCSD(T)/CBS limiting value for the correlation contribution to the interaction energy could be obtained as:   smallbasiscorr MP smallbasiscorr YTCCSD CBScorr MP CBScorr TCCSD EEEE , 2 ,)( , 2 ,)(  . (eq. 3.43) In this expression, the MP2 contribution to the correlation energy is estimated to basis limit with the extrapolation procedures explained above, and the result is corrected from the inefficiencies of the N-electron model, by using a CCSD(T) calculation with a small basis set, and assuming that the difference between CCSD(T) and MP2 is fairly constant. An alternative way of understanding eq. 3.43 is by considering that the correlation contribution is obtained at the CCSD(T) level with a small basis set, and the basis set incompleteness error is estimated at the MP2 level, eq. 3.44.   smallbasiscorr MP CBScorr MP smallbasiscorr YTCCSD CBScorr TCCSD EEEE , 2 , 2 ,)( ,)(  (eq. 3.44) In any case, this procedure is commonly used in order to obtain benchmark values for the interaction energies of complexes. It should not be forgotten however that errors are still present, depending on the size of the small basis set employed for the CCSD(T) calculation. Applying this approach, several dataset have been constructed which are widely used as reference for developing approximate methods.[72-75] Still within this approach, for moderate-sized systems the bottleneck of the calculations is the CCSD(T) calculation that, even with a basis set of aug-cc-pVDZ quality can be very demanding. Thus, Hobza and Rezak proposed an alternative way of estimating the correction to the MP2/CBS value avoiding the expensive CCSD(T) calculation. The first of these proposals (MP2.5)[76] and the subsequent refinement (MP2.X)[77] substitute the CCSD(T) calculation by cheaper MP3 ones. Thus in the MP2.X approach, the correlation contribution to the interaction energy is obtained as:   smallbasiscorr MP smallbasiscorr MP CBScorr MP CBScorr TCCSD EECEE , 2 , 3 , 2 ,)(  . (eq. 3.45) Hobza and Rezak introduced an empirical scaling coefficient obtained by fitting to the CCSD(T)/CBS estimates of a set of complexes of different nature. With this approach it has been observed that the accuracy of the MP2.X results is almost independent of the basis set employed for the MP3 calculation given a proper C coefficient. Thus, modest basis sets as 6-31G* can be employed giving results pretty similar to those obtained by performing an actual CCSD(T) calculation.[78, 79] As a consequence, the MP2.X procedure allows saving computational resources by using the cheaper MP3 method, but also allowing the use of smaller basis set in order to estimate the N-electron correction to the MP2 limit. 3. Methodology 49 3.3. Density Functional Theory Methods Density Functional Theory (DFT) is an alternative to methods based on the calculation of the wavefunction as the goal for the description of the system and its properties. In the last decades DFT has been developed intensely because it becomes very useful for a lot of different kind of systems due to its advantageous relation between the quality of the results and computational cost. Like post-HF methods, DFT includes the electronic correlation term but with a cost similar to HF calculations.[32, 80] The foundation of DFT is based on the idea that the information that can be extracted from the wavefunction can also be obtained from the electronic density. Whereas the energy in wavefunction-based method is a functional of the wavefunction (N spatial coordinates + N spin coordinates in a N-electron system), which is dependent on the number of atoms, in DFT the energy is a functional of the electron density, so it only depends on the three spatial coordinates   )r(fE    , (eq. 3.46) In order to work with the density functional theory it is necessary to apply the first theorem of Hohenberg and Kohn (1964).[32, 81] “Any observable of a stationary non-degenerate ground state can be calculated, exactly in theory, from the electron density of the ground state. In other words, any observable can be written as a functional of the electron density of the ground state”. In a system with M nuclei of charge Za located in Ra, the interaction of N electrons of the system with the M nuclei can be described through V operator. rdrrrdrr Rr Z Rr Z Vi i ai a i a ai a     )()()( ||||           , (eq. 3.47) where )(r   is the effective potential of an electron and  ii rrr )()(   is the operator for the electronic density. The rest of the operators in the Hamiltonian depend exclusively on the coordinates of the electrons, so their expressions are equal in all systems and only the number of electrons changes. Therefore, the total electronic energy depends only on the total number of electrons (N) and external potential )(r   . In addition to this, other essential characteristic of the first Hohenberg-Kohn theorem is that the relationship between energy and density is univocal. Therefore, if there exist two different external potentials that generate the same electronic density, then both potentials must be the same. On the other hand, the second Hohenberg-Kohn theorem provides a variational principle:[32, 81] Alba Campo Cacharrón 50 “The electron density of a non-degenerate ground state can be calculated, exactly in theory, determining the density that minimizes the energy of the ground state”. The second Hohenberg-Kohn theorem establishes the existence of a variational principle for the energy, verifying that for a test electronic density )(r   , the energy for such system is larger than or equal to the energy for the real ground state of the system. So, the energy reaches a minimum value for the exact ground state.     )()( 00 rErE    . (eq. 3.48) Therefore, the differential equation   0 )( )(     r rE   (eq. 3.49) Is fulfilled, in which the Lagrange multiplier ensures that the electronic density is normalized for N electrons. Therefore, the density can be variationally optimized in order to approach the real ground state density. Starting from these two fundamental theorems, the main goal of the DFT methods consists on designing functionals that connect energy and electronic density but, unfortunately, Hohenberg-Kohn theorems do not establish how the exact connection between both magnitudes is. In order to solve this problem, Kohn and Sham developed a practical application of this theory through a method with a formulation similar in structure to the Hartree-Fock method. 3.3.1. Kohn-Sham procedure A general expression for the energy taking into account the Born Oppenheimer approximation could be the next one: [32, 80]          eene EETE  . (eq. 3.50) The energy is divided into three parts: kinetic energy    T , attraction between the nuclei and electrons    ne E and electron-electron repulsion    ee E (the nuclear-nuclear repulsion is a constant within the Born-Oppenheimer approximation). Furthermore, similarly as in Hartree-Fock theory, the    ee E term may be divided into a Coulomb and Exchange part    J and    K , though implicitly including correlation energy in all terms. The    ne E and    J functionals are given by their classical expressions, where the ½ factor in    J allows the integration to run over all space for both variables.   aa a ne rd rR rZ E    || )( ][   , (eq. 3.51) 3. Methodology 51 ' |'| )'()( 2 1 ][ rdrd rr rr J       . (eq. 3.52) The basic idea in the Kohn-Sham formalism is splitting the kinetic energy functional into two parts, one which can be calculated exactly and a small correction term. In the real molecular system, the electronic density and the kinetic energy will be written as: ii ii T  2 2 1 ][   , (eq. 3.53)   ii 2  . (eq. 3.54) Kohn-Sham formalism establishes then the calculation of the kinetic energy under the assumption of a reference system of non-interacting electrons S (in the same sense as HF orbitals in wave mechanics describe non-interacting electrons) but under an external potential such that the density is the same as in the real system. The solution, within the Born-Oppenheimer approximation, is provided by Schrödinger equation, and the system could be described by a Slater determinant of molecular orbitals  i which have the following exact kinetic energy: ii iiS T  2 2 1 ][   , (eq. 3.55)   iiS 2  . (eq. 3.56) Of course, the electrons interact among themselves and the previous equations do not provide the correct total kinetic energy. However, just as HF theory provides ~99% of the correct answer, the difference between the exact kinetic energy and that calculated by assuming non-interacting electrons is small. The remaining kinetic energy not included in ][  S T is absorbed into an exchange-correlation term, and a general DFT expression for the energy can be written as:            xcneSDFT EJETE  . (eq. 3.57) If ][  xc E is expanded up, it becomes easier to understand which are the contributions to this exchange-correlation energy.                JETTE eeSxc  . (eq. 3.58) The problem is similar to that encountered in wave mechanics HF theory: determining a set of orthogonal orbitals that minimize the energy. The electron density is expressed Alba Campo Cacharrón 58 As it can be observed the integral can be expressed as a product of two generalized densities. Then, these densities can be approximated by using a linear expansion employing an adequate auxiliary basis set: )()( rr  Nbas PP pq Ppq d  . (eq. 3.70) There are different methods for obtaining the expansion coefficients, but one of the most common options leads to the following fitting coefficients:   PPQ pq PPpqd])[|( 1 J (eq. 3.71) with )( 1 )()()|( 2 12 1121 rrrrr Pqp r ddPpq     (eq. 3.72) and )( 1 )( 2 12 121 rrrr QPPQ r ddJ     . (eq. 3.73) Thus, making the proper substitutions, the four-index integral can be expressed as: )|(])[|()|()|( 1rsQPpqrsQdrspq PQQ pq Q   J . (eq. 3.74) The important aspect in this expression is that the four-index integral has been factorized into an expression running through three indexes. One advantage of this procedure is that the storage requirements are greatly reduced when only three-index integrals are needed at most. Also, three- and two-index integrals are more easily evaluated that the corresponding four-index ones, thus saving computational time. This kind of approximation expressing four-index integrals in two or three-index ones, can be applied to different methods, the computational saving being different depending on each case. This technique has been applied very successfully in order to speed up DFT calculations employing pure functionals, where the time for Coulomb contribution can be greatly reduced. Also, MP2 calculations benefit from this approach since the correlation contribution comes from two-electron integrals involving two occupied and two virtual orbitals. Resolution of the Identity can also be applied to the exchange contribution (RI-JK), though in this case the speedup is not as advantageous as for the Coulomb part. 3. Methodology 59 The central point in the resolution of the identity approach is the design of proper auxiliary basis sets which are capable of reproducing the density introducing negligible errors. Different sets of auxiliary basis functions are available in literature specifically designed for applying the RI approach for the Coulomb,[112] exchange[113] or correlation[114] energy calculations. 3.5. Interaction Energy Partitioning Applying the above-described supermolecule method with any of the wavefunctionbased or DFT methods provides a magnitude for the interaction energy for a given geometry. However, it would be desirable to obtain more information about the nature, origins and characteristics of the interaction itself, but this kind of information is not provided by the supermolecule approach. Considering the classical and historical description of intermolecular forces, the interaction is usually rationalized in terms of contributions from electrostatics, repulsion, polarization and dispersion. Thus, a method providing such a kind of physical partitioning would be desirable, and it will be the subject of this section.[2, 3] 3.5.1. Energy Decomposition Analysis (EDA) Methods With the common denomination of Energy Decomposition Analysis (EDA) there is a variety of methods devised for partitioning the interaction energy obtained from variational supermolecule HF and DFT methods.[4, 115] Most EDA methods are variants of the original partitioning scheme proposed by Kitaura and Morokuma.[116] In the original formulation the interaction energy of a dimer was decomposed in electrostatic, repulsion, polarization and charge transfer contributions as obtained at the HF level. Therefore, no dispersion contribution was obtained. EDA methods have been later extended to many-body systems by Chen and Gordon, thus allowing the partitioning in trimers, tetramers or larger clusters.[117] In any case, the original EDA partitioning posed some problems, and the attempts of solving them give way to a variety of EDA methods, such as the natural energy decomposition analysis (NEDA),[118-120] the reduced variational space analysis (RVS),[117, 121] and the generalized Kohn-Sham EDA (GKS-EDA)[122] Following there is a general description of how most EDA methods work.[4] The formation of the complex is divided in a series of steps, each one associated with a given physical contribution. Consider a dimer AB formed by two separate units A and B already in the final geometry they have in the complex. 1. In the first step of dimer formation according to EDA, fragments A and B with frozen charge distribution are taken from infinite separation and brought together to the position in the dimer. The interaction between the frozen charge densities of A and B gives the electrostatic interaction: Alba Campo Cacharrón 60       12 21 )()( )()()()( r rr drdrrrdrVrrdrV R ZZ EBA A B BAABelec       . (eq. 3.75) 2. In the second step of EDA the product wavefunction employed before, normalized but not antisymmetrical, is antisymmetrized. As a consequence of the antisymmetrization the energy changes, leading to a term called Exchange or Pauli repulsion. 3. The wavefunction is allowed to relax to give the final state of the dimer AB with energy EAB. The energy lowering associated with this orbital relaxation is the orbital relaxation contribution. Therefore, the interaction energy is split into three contributions as indicated in eq. 3.76, orbitalPauliticelectrostaAB EEEE  int . (eq. 3.76) This is the global framework, but different EDA approaches further split the interaction energy by dividing Pauli and orbital relaxation terms. Depending on the method employed, Pauli repulsion can be split as Exchange + Repulsion, and orbital relaxation is further divided in A polarization, B polarization and charge transfer terms. EDA partitioning has been employed in one chapter of this thesis. The LMO-EDA method[123] used can be considered an extension and modification of the methods developed by Kitaura and Morokuma[116], Ziegler and Rauk,[124] and Hayes and Stone.[125] The main features of this partitioning are listed as follows: 1. The electrostatic, exchange and repulsion terms are obtained from the Heitler- London interaction energy as proposed by Hayes and Stone from an antisymmetric product of the monomer HF spinorbitals. 2. Polarization energy is defined as orbital relaxation energy going from monomer orbitals to dimer orbitals. Thus, no charge transfer term is defined. 3. A dispersion contribution can be computed via a supermolecule calculation with a post-HF or DFT method. Dispersion is defined as the difference between the sum of all contributions and the interaction energy obtained with the post-HF method. Though this has been common practice, it has to be taken into account that using HF partitioning, the dispersion contribution really corresponds to correlation energy contribution, containing other effects over electrostatic, induction and repulsion terms. 4. The partitioning can be applied similarly to DFT methods. 3. Methodology 61 In summary, LMO-EDA allows partitioning the global interaction energy into contributions that can give hints on the effects controlling the interaction. Within the Hartree-Fock framework the interaction energy is expressed as: polrepexchelec HF AB EEEEE  int, . (eq. 3.77) Applying a post-HF method as, say, CCSD(T), a dispersion contribution is defined: disppolrepexchelec TCCSD AB EEEEEE  )(int, , (eq. 3.78) with dispersion (correlation contribution) defined as: HF AB TCCSD ABdisp EEE  )( . (eq. 3.79) In the case of employing a DFT method a similar partition is employed, so disppolrepexchelec DFT AB EEEEEE  int, . (eq. 3.80) Dispersion is obtained as: )( int, polrepexchelec DFT ABdisp EEEEEE  . (eq. 3.81) 3.5.2. Symmetry-Adapted Perturbation Theory Methods (SAPT) Symmetry Adapted Perturbation theory methods employ a totally different approach in order to directly obtain the interaction energy.[126, 127] SAPT is based on perturbation theory, where the interaction itself is treated as the perturbation and its magnitude is directly computed. Applying a perturbational scheme (Rayleigh-Schrodinger) to a complex naturally provides expressions which can be identified with contributions from electrostatic, induction and dispersion. The most obvious partitioning of the Hamiltonian operator of an interacting pair of molecules A and B is the following: VHVHHH BA ˆˆˆˆˆˆ 0 , (eq. 3.82) where 0 ˆ H is the solution for the unperturbed system and V ˆ is the operator for the intermolecular interaction. The reference system consists on the isolated noninteracting molecules and the perturbation is the interaction. AAAA EH   and equivalently for the monomer B BBBB EH   . As a consequence, the reference wavefunction is the product of the isolated molecules’ wavefunctions BA 000   . Applying Rayleigh-Schrödinger perturbation theory under these assumptions, the typical expression can be obtained for the different corrections to first, second, … order.[3] That way, the first order correction will be: Alba Campo Cacharrón 62 BABA el VE 0000 )1(   . (eq. 3.83) This expression corresponds to the coulombic interaction between the electron density of both monomers, and is therefore associated to the electrostatic energy. To second order, terms appear depending on single and double excitations which can be associated with induction and dispersion contributions. Thus:    00 2 000 mAA m BA m BA A ind EE V E  (eq. 3.84) corresponds to single excitations within monomer A due to the presence of the nearby monomer B, and therefore is associated to the induction contribution of monomer A. An equivalent expression is also found for single excitations involving B    00 2 000 nBB n B n ABA B ind EE V E  . (eq. 3.85) Finally, to second order another term remains involving double excitations which is assigned as dispersion contribution.    0;0 00 2 00 nm BAB n A m B n A m BA disp EEEE V E  . (eq. 3.86) It is worth noting that these expressions, together with the multipole expansion, provide much of our qualitative and quantitative discussion on the role of noncovalent bonding forces, allowing to describe intermolecular interactions as functions of molecular properties such as multipoles or polarizabilities.[2, 3] The perturbation theory for well-separated molecules described above (the long-range approximation or polarization approximation) is very successful if the molecules are a long distance apart, but at short range fails completely.[3] Part of the reason for the failure of the theory as usually formulated is that the multipole expansion breaks down. A more fundamental failure is that the repulsion between molecules that occurs at short-range is completely missed from the long-range theory. This failure arises from the fact that if the molecules are close enough for their wavefunctions to overlap, exchange cannot be ignored. The source of the difficulties of Rayleigh-Schrödinger perturbation theory for describing intermolecular interactions at short range is the wrong symmetry of the reference wavefunction.[2, 3] The reference wavefuntion BA 000   is antisymmetric upon electron exchange within A or within B, but it is not upon exchange of electrons 3. Methodology 63 between A and B. Therefore, the reference wavefunction is physically unsound, thus leading to the lack of proper repulsion forces as the molecules come together. Therefore, the reference wavefunctions should be properly antisymmetrized in order to satisfy the Pauli principle. There have been different approaches to overcome the antisymmetry problem, but the prevailing one is the symmetrized Rayleigh-Schrödinger theory, a nowadays synonym of Symmetry Adapted Perturbation Theory (SAPT). Within SAPT, antisymmetry is forced in the energy expressions, modifying the electron density in such a way as to cause a repulsive force on the nuclei. This is the force corresponding to the exchangerepulsion energy.[126, 127] The outcome of SAPT procedure is a series of contributions that to low order can be associated to physical effects as in polarization theory. The main difference, however, is that each of the polarization terms is now accompanied by an exchange-repulsion term arising from proper antisymmetry. Therefore, the interaction energy is expressed in SAPT as: ... )2()2()2()2()1()1( int   dispexchdispindexchindexchel EEEEEEE . (eq. 3.87) Thus, it can be seen that antisymmetrization produces new terms that are completely missed when a simple product wavefunction is used. The most important fact is that there is a strong repulsion between closed-shell molecules when their wavefunctions overlap significantly. The description of SAPT presented above assumed that one knows the exact wavefunctions for monomers. In practice, the wavefunctions are computed separating the HF and the correlation contributions. Since electron correlation significantly affects molecular properties required as input to SAPT it is mandatory to account for intramolecular electron correlation. [52, 126, 127] If the MP decomposition of the Hamiltonians is used for monomers, SAPT becomes a double-perturbation theory according to the following splitting of the Hamiltonian, WVFH  (eq. 3.88) so the sum of monomer’s Fock operators BA FFF  is now the unperturbed Hamiltonian, whereas V and the sum of MP potentials of monomers BA WWW  are the two perturbation operators. Accordingly, the interaction energy can now be expressed as      0;1 int ji ij exch ij pol EEE , (eq. 3.89) Alba Campo Cacharrón 64 where ij is the order in V-W. A large number of terms in this expansion have been developed up to i=3 and j from 0 to 4, depending on i and the physical component. Thus, SAPT allows describing the interaction with increasing accuracy by including terms to higher orders in both expansions. Depending on the contributions included, different models can be defined, though the series expansion is typically truncated at secondorder in V, resulting in a complete neglect of third- and higher order terms. It has been observed that in polar systems higher order effects can be important, often associated with induction effects. A correction is recommended in order to at least introduce an estimation of such effects. The correction is defined as )20()20()10()10( int indexchindexchel HF HF EEEEE    , (eq. 3.90) where HF Eint is the Hartree-Fock interaction energy calculated using the supermolecular method. Therefore, SAPT can partition the interaction energy at the HF level as: HFindexchindexchel HF EEEEE    )20()20()10()10( int . (eq. 3.91) Of course, at the HF level, the interaction energy does not include dispersion, but applying perturbation theory based on HF wavefunctions the dispersion term can be obtained in SAPT. Thus, the Hartree-Fock calculation can be corrected with the dispersion contribution obtaining the HF+D method (in a similar way as the already commented DFT-D methods). )20()20()20()20()10()10( int dispexchdispHFindexchindexchel DHF EEEEEEE    . (eq. 3.92) The last expression contains all SAPT contributions to order 2 obtained employing HF wavefunctions. As commented above, higher-order terms in induction are taken care of by the HF  contribution. The next step would be including intramonomer correlation effects, leading to a more accurate description of the interactions. For example, the so-called SAPT2 level corresponds to: )22()22()12()11()12( intint indexchindexchexchel DHFSAPT EEEEEEE   , (eq. 3.93) including corrections similar to MP2. More accurate models would include contributions to third order both in the intermolecular perturbation and the intramonomer correlation perturbation. The problem of using such expression is its high computational demand, which has led to other approaches to be commented in the following section. In any case, it is worth noting that, recently, Hohenstein and Sherrill have developed a SAPT scheme exploiting the resolution of the identity approach in order to reduce the 3. Methodology 65 computational cost, which could make SAPT calculations more feasible.[52, 128, 129] This approach has been coded into the PSI4 program.[130] 3.5.2.1. SAPT(DFT) Despite the successes of SAPT, the calculation of the intramonomer correlation terms makes it computationally prohibitive for larger molecular systems. Williams and Chabalowski suggested that if a correlated description of the monomers was used, the costly correlation terms can be avoided.[131] Due to computational considerations, Williams and Chabalowski suggested that a DFT description of the monomers would be best suited. This fact allows SAPT to be performed on much larger systems than previously allowed although the initial results were rather poor. Thus, the idea was simply changing the HF orbitals and energies by their Kohn-Sham counterparts. This way, a SAPT calculation including only the perturbation on the intermolecular interaction will suffice, because intramonomer correlation effects were already taken care of in the DFT calculations. [131] The results obtained following this procedure were disappointing. It was observed that one of the reasons for this inaccuracy was the incorrect asymptotic behavior of the exchange-correlation functions for obtaining the Kohn-Sham energies. While the exchange-correlation potential should decay as 1/r for a neutral system, the standard local or/and gradient-corrected DFT exchange-correlation potentials decay too quickly. The method was soon improved by Hesselmann and Jansen[132-136] and Szalewicz[137-140] independently, leading to essentially identical methods. In both proposals, a correction for the asymptotic behavior is assumed. In Jansen’s proposal the functional employed is combined with the LB94 functional, which has a correct asymptotic behavior at longrange (but fails at short-range).[136] This fact creates a new problem that it is the breaking at intermediate distances between the behavior of the two functionals.[141, 142] In order to arrange this, a gradient-regulated connection method as developed by Grüning et al. was employed.[141] This scheme (Adiabatic Correction AC) requires the sum HOMO IP   as the input parameter (which vanishes in the case of exact KS DFT). The value of HOMO  is obtained from a calculation with the uncorrected xc functional, whereas the ionization potentials can either be taken from experiment or calculated from the difference of KS DFT calculations of the neutral and the ionized systems, respectively. Once these problems affecting accuracy are solved, the equation obtained for SAPT(DFT) is the following: HFdispexchdisprespindexchrespindexchel DFTSAPT EEEEEEE    )20()20()20(, )20(, )10()10( int . (eq. 3.94) The same correction term HF  is included in the calculation in order to take into account higher-order contributions. Despite of the difficulties in the implementation of Alba Campo Cacharrón 66 DFT in a complicated theory like SAPT, SAPT(DFT) presents many advantages in comparison to SAPT. The main reason, and the origin of this method, is the decreased computational cost for the description of intramonomer correlation, which allows the application of SAPT partitioning to larger systems that could not be afforded with SAPT based in MBPT or CC. In addition to this, it is worth noting that density fitting approach can be employed, greatly reducing the computational cost. Also, simplified extrapolation schemes have been proposed in order to reach the complete basis set limit in the framework of SAPT(DFT) calculations.[143] 3.6. Electron density analysis Atoms in Molecules quantum theory (QTAIM), that was developed by Bader, is a theory based on the topological analysis of the electron density, which brought quantum mechanics into applicability to an atom within a molecule.[144] When a molecular property can be expressed in terms of a property density, the contribution of a given atom to that molecular property can be obtained by integrating this density over the volume of the atom in the molecule. Thus, the theory relates concepts as bonding, functional groups or chemical reactivity to the topology of the underlying electron density, though the aspect of the molecular density is always the same: large cups at the nucleus and an exponential decreasing behavior in all directions. The characterization of critical points is based on the behavior of the gradient vector in the zone nearby. So, a critical point in the electron density can be defined as a point in space at which the first derivates of the density vanish, kji dz d k dy d j dx d izyx     , (eq. 3.95) meaning that each individual derivative in the gradient operator,  , is zero and not just their sum. The gradient of a scalar function such as )(r  at a point in space is a vector pointing in the direction in which )(r  undergoes the greatest rate of increase and having a magnitude equal to the rate of increase in that direction. Considering the second derivatives, the elements of the tensor   , one can discriminate between a local minimum, a local maximum, or a saddle point. There are nine second derivatives of )(r  that can be arranged in the so-called Hessian matrix, which when diagonalized at a critical point rc gives: 3. Methodology 67                                      3 2 1 ' 2 2 2 2 2 2 00 00 00 00 00 00       c rr z y x . (eq. 3.96) 1  , 2  and 3  are the eigenvalues of the Hessian matrix (ordered as 321   ) and represent the curvatures of the density with respect to the three principal axes x’, y’ and z’. There are four types of stable critical points having three non-zero eingenvalues:  1  , 2  , 3  <0 ; (3,-3): all curvatures are negative. ρ is a local maximum. This is called a Nuclear Critical Point.  1  , 2  <0, 3  >0 ; (3,-1): two negative curvatures. ρ is a maximum in the plane defined by the corresponding eigenvectors and a minimum along the third axis which is perpendicular to this plane. This is called a Bond Critical Point.  1  <0, 2  , 3  >0 ; (3,+1): two positive curvatures. ρ is a minimum in the plane defined by the corresponding eigenvectors and a maximum along the third axis which is perpendicular to this plane. This is called a Ring Critical Point.  1  , 2  , 3  >0 ; (3,+3): Three curvatures are positive. ρ is a local minimum and corresponds to a Cage Critical Point. The collection of paths linking the nuclei of bonded atoms in an equilibrium geometry with the associated critical points is known as the molecular graph, which provides an unambiguous definition of the “molecular structure” and can thus be used to locate changes in structure along a reaction path. Also, the strength of a chemical bond is reflected on the electron density at the bond critical point )( b  . b  is greater than 0.20 a.u. in shared (covalent) bonding and less than 0.10 a.u. in a closed-shell interaction (ionic, vdW, hydrogen bond, etc.). b  has been shown to be strongly correlated with the binding energy for several types of bonding interaction. In conclusion, the analysis of the electron density and its derivatives through its main topological features provides much information about the characteristics of the system and can be employed to analyze the most relevant noncovalent interactions present in the systems under study. 3.6.1. NCI index As commented in the preceding section, Atoms in Molecules Theory can be employed to analyze in detail the characteristics of the interaction in intermolecular systems. However, in large systems, QTAIM provides an overwhelming amount of information: a large number of critical points, their characteristics, a large set of bond paths, etc. Alba Campo Cacharrón 74 x f      1 )( , (eq. 3.104) where x is a fitted parameter, originally set to 0.5, though other values have been proposed. The choice of this parameter has little impact in solvents with large dielectric constant whereas the results are more affected if the solvent shows low dielectric constant. Figure 3.6. Cavity construction for the PCM calculations as performed with Gaussian. Therefore, the solvent is treated as a dielectric, with apparent charges defining the reaction field obtained applying conductor boundary conditions, and properly scaled as to represent a specific value of the dielectric constant. COSMO is especially robust with respect to artifacts that result from the small percentage of the solute’s charge reaching the outer side of the cavity, the so-called outlying charge. In any unconstrained electronic structure calculation of the solute there is a tail that extends beyond the cavity surface. Therefore, there is a small percentage of the electron density in the dielectric.[149] A proper formulation for this problem should include not only the surface polarization but also a volume polarization outside the cavity, introducing apparent volume charges at points outside the cavity. COSMO includes an approximation to volume polarization making a final calculation extending somewhat the limits of the cavity and comparing the values obtained with the original and the extended surfaces. It should be taken into account that outlying charge errors can easily 3. Methodology 75 reach 20% for neutral species, and be even larger for anions, so a proper correction for these effects should be considered. This has been a problem in D-PCM methods which led to incorporate COSMO into the PCM suite (C-PCM), though the latest developments in PCM models such as the integral equation formalism variant PCM (IEFPCM) include an approximate correction to this problem.[147] 3.8. References [1] P. Hobza, R. Zaradnik, Intermolecular complexes : the role of van der Waals systems in physical chemistry and the biodisciplines, Elsevier, Amsterdam, 1988. [2] I. Kaplan, Intermolecular Interactions : Physical Picture , Computational Methods, John Wiley & Sons, Chichester, 2006. [3] A. J. Stone, The theory of intermolecular forces, Oxford University Press, Oxford, 2013. [4] M. v. Hopffgarten, G. Frenking. WIREs Comput. Mol. Sci. 2012, 2, 43-62. [5] G. Chalasinski, M. M. Szczesniak. Chem. Rev. 1994, 94, 1723-1765. [6] G. Cha asi ski, M. M. Szcz niak. Chem. Rev. 2000, 100, 4227-4252. [7] K. Szalewicz, B. Jeziorski. J. Chem. Phys. 1998, 109, 1198-1200. [8] S. S. Xantheas. J. Chem. Phys. 1996, 104, 8821-8824. [9] E. Cabaleiro-Lago, J. Rodríguez-Otero, Á. Peña-Gallego. Theor. Chem. Acc. 2011, 128, 531-539. [10] E. M. Cabaleiro-Lago, M. A. R os. J. Chem. Phys. 2000, 112, 2155-2163. [11] A. Campo-Cacharrón, E. M. Cabaleiro-Lago, J. Rodríguez-Otero. ChemPhysChem 2012, 13, 570-577. [12] E. M. Cabaleiro-Lago, Á. Peña-Gallego, J. Rodríguez-Otero. J. Chem. Phys. 2008, 128, 194311/1-194311/8. [13] C. D. Sherrill Computations of Noncovalent π Interactions, in Reviews in Computational Chemistry Vol. 26, John Wiley & Sons, New York, 2009. [14] F. B. van Duijneveldt, J. G. C. M. van Duijneveldt-van de Rijdt, J. H. van Lenthe. Chem. Rev. 1994, 94, 1873-1885. [15] H. B. Jansen, P. Ros. Chem. Phys. Lett. 1969, 3, 140-143. [16] S. F. Boys, F. Bernardi. Mol. Phys. 1970, 19, 553-566. [17] J. Alvarez-Idaboy, A. Galano. Theor. Chem. Acc. 2010, 126, 75-85. [18] A. J. C. Varandas. J. Phys. Chem. A 2009, 114, 8505-8516. [19] Ł. M. Mentel, E. J. Baerends. J. Chem. Theory Comput. 2013, 10, 252-267. [20] T. Helgaker, P. Jørgensen, J. Olsen, Molecular electronic-structure theory John Wiley & Sons, Chichester, 2000. [21] T. Helgaker, W. Klopper, D. P. Tew. Mol. Phys. 2008, 106, 2107-2143. [22] R. M. Balabin. J. Chem. Phys. 2010, 132, 211103/1-211103/4. [23] D. Asturiol, M. Duran, P. Salvador. J. Chem. Phys. 2008, 128, 144108/2-144108/5. [24] H. Valdés, V. Klusák, M. Pitoňák, O. Exner, I. Starý, P. Hobza, L. Rul šek. J. Comput. Chem. 2008, 29, 861-870. [25] D. Asturiol, M. Duran, P. Salvador. J. Chem. Theory Comput. 2009, 5, 2574-2581. [26] M. J. Elrod, R. J. Saykally. Chem. Rev. 1994, 94, 1975-1997. [27] E. E. Dahlke, D. G. Truhlar. J. Phys. Chem. B 2006, 110, 10595-10601. Alba Campo Cacharrón 76 [28] K. Walczak, J. Friedrich, M. Dolg. J. Chem. Phys. 2011, 135, 134118/1-134118/11. [29] J. F. Ouyang, M. W. Cvitkovic, R. P. A. Bettens. J. Chem. Theory Comput. 2014, 10, 3699-3707. [30] A. Szabo, N. S. Ostlund. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory, Dover Publications, 2012. [31] C. J. Cramer, Essentials of Computational Chemistry: Theories and Models, 2nd Edition, John Wiley & Sons, Chichester, 2004. [32] F. Jensen, Introduction to computational chemistry, John Wiley and Sons, Chichester, 2001. [33] D. P. Tew, W. Klopper, T. Helgaker. J. Comput. Chem. 2007, 28, 1307-1320. [34] C. D. Sherrill, T. Takatani, E. G. Hohenstein. J. Phys. Chem. A 2009, 113, 10146- 10159. [35] S. Tsuzuki, K. Honda, T. Uchimaru, M. Mikami. J. Chem. Phys. 2004, 120, 647-659. [36] S. Tsuzuki, K. Honda, T. Uchimaru, M. Mikami, K. Tanabe. J. Am. Chem. Soc. 2001, 124, 104-112. [37] M. O. Sinnokrot, E. F. Valeev, C. D. Sherrill. J. Am. Chem. Soc. 2002, 124, 10887- 10893. [38] S. Grimme. J. Chem. Phys. 2003, 118, 9095-9102. [39] S. Grimme, L. Goerigk, R. F. Fink.WIRES Comput. Mol. Sci. 2012, 2, 886-906. [40] M. Gerenkamp, S. Grimme. Chem. Phys. Lett. 2004, 392, 229-235. [41] T. P. M. Goumans, A. W. Ehlers, K. Lammertsma, E.-U. Wuerthwein, S. Grimme. Chem. Eur. J. 2004, 10, 6468-6475. [42] S. Grimme. J. Phys. Chem. A 2005, 109, 3067-3077. [43] J. Antony, S. Grimme. J. Phys. Chem. A 2007, 111, 4862-4868. [44] F. Neese, T. Schwabe, S. Kossmann, B. Schirmer, S. Grimme. J. Chem. Theory Comput. 2009, 5, 3060-3073. [45] K. E. Riley, J. A. Platts, J Řezáč, P. Hobza, J. G. Hill. J. Phys. Chem. A 2012, 116, 4159-4169. [46] J. G. Hill, J. A. Platts. J. Chem. Theory Comput. 2006, 3, 80-85. [47] R. A. Distasio Jr, M. Head-Gordon. Mol. Phys. 2007, 105, 1073-1083. [48] S. Grimme. J. Comput. Chem. 2003, 24, 1529-1537. [49] I. Hyla-Kryspin, S. Grimme. Organometallics 2004, 23, 5581-5592. [50] M. Pitonak, J. Řezáč, P. Hobza. Phys. Chem. Chem. Phys. 2010, 12, 9611-9614. [51] T. Takatani, E. G. Hohenstein, C. D. Sherrill. J. Chem. Phys. 2008, 128, 124111/1- 124111/8. [52] E. G. Hohenstein, C. D. Sherrill. WIREs Comput. Mol. Sci. 2012, 2, 304-326. [53] A. Hellweg, S. A. Grün, C. Hättig. Phys. Chem. Chem. Phys. 2008, 10, 4119–4127. [54] J. Řezáč, P. Hobza. J. Chem. Theory Comput. 2013, 9, 2151-2155. [55] J. Řezáč, L. Šimová, P. Hobza. J. Chem. Theory Comput. 2012, 9, 364-369. [56] D. G. A. Smith, P. Jankowski, M. Slawik, H. A. Witek, K. Patkowski. J. Chem. Theory Comput. 2014, 10, 3140-3150. [57] C. Hättig, W. Klopper, A. Köhn, D. P. Tew. Chem. Rev. 2011, 112, 4-74. [58] L. Kong, F. A. Bischoff, E. F. Valeev. Chem. Rev. 2011, 112, 75-107. [59] T. Helgaker, W. Klopper, H. Koch, J. Noga. J. Chem. Phys. 1997, 106, 9639-9646. [60] A. Halkier, T. Helgaker, P. Jørgensen, W. Klopper, H. Koch, J. Olsen, A. K. Wilson. Chem. Phys. Lett. 1998, 286, 243-252. [61] D. G. Truhlar. Chem. Phys. Lett. 1998, 294, 45-48. 3. Methodology 77 [62] A. J. C. Varandas. J. Chem. Phys. 2007, 126, 244105-244120. [63] A. Halkier, T. Helgaker, P. Jørgensen, W. Klopper, J. Olsen. Chem. Phys. Lett. 1999, 302, 437-446. [64] D. G. Liakos, F. Neese. J. Phys. Chem. A 2012, 116, 4801-4816. [65] F. Neese, E. F. Valeev. J. Chem. Theory Comput. 2010, 7, 33-43. [66] R. A. Kendall, T. H. Dunning, R. J. Harrison. J. Chem. Phys. 1992, 96, 6796-6806. [67] W. Klopper, W. Kutzelnigg. J. Mol. Struct. (THEOCHEM) 1986, 135, 339-356. [68] A. Karton, J. L. Martin. Theor. Chem. Acc. 2006, 115, 330-333. [69] F. Jensen. Theor. Chem. Acc. 2005, 113, 267-273. [70] K. E. Riley, P. Hobza. WIREs Comput. Mol. Sci. 2011, 1, 3-17. [71] P. Jurecka, J. Sponer, J. Cerny, P. Hobza. Phys. Chem. Chem. Phys. 2006, 8, 1985- 1993. [72] J. Řezáč, K. E. Riley, P. Hobza. J. Chem. Theory Comput. 2014, 10, 1359-1360. [73] J. Řezáč, K. E. Riley, P. Hobza. J. Chem. Theory Comput. 2011, 7, 2427-2438. [74] P. Hobza, J. Šponer. J. Am. Chem. Soc. 2002, 124, 11802-11808. [75] P. Hobza in http://www.begdb.com. Last accesed: 15/09/14. [76] M. Pitoňák, P. Neogrády, J. Černý, S. Grimme, P. Hobza. ChemPhysChem 2009, 10, 282-289. [77] K. E. Riley, J. Rezac, P. Hobza. Phys. Chem. Chem. Phys. 2011, 13, 21121-21125. [78] R. Sedlak, K. E. Riley, J. Řezáč, M. Pitoňák, P. Hobza. ChemPhysChem 2013, 14, 698-707. [79] K. E. Riley, J. Rezac, P. Hobza. Phys. Chem. Chem. Phys. 2012, 14, 13187-13193. [80] W. Koch, M. C. Holthausen, A Chemist's guide to density functional theory, Wiley- VCH, Weinhem, 2000. [81] P. Hohenberg, W. Kohn. Phys. Rev. 1964, 136, B864-B871. [82] A. D. Becke. Phys. Rev. A 1988, 38, 3098-3100. [83] C. Lee, W. Yang, R. G. Parr. Phys. Rev. B 1988, 37, 785-789. [84] B. Miehlich, A. Savin, H. Stoll, H. Preuss. Chem. Phys. Lett. 1989, 157, 200-206. [85] A. D. Becke. J. Chem. Phys. 1993, 98, 1372-1377. [86] Y. Zhao, D. G. Truhlar. J. Chem. Phys. 2006, 125, 194101-194119. [87] Y. Zhao, D. G. Truhlar. Theor. Chem. Acc. 2007, 120, 215-241. [88] Y. Zhao, D. Truhlar. Theor. Chem. Acc. 2008, 120, 215-241. [89] Y. Zhao, D. G. Truhlar. J. Chem. Theory Comput. 2008, 4, 1849-1868. [90] R. Peverati, D. G. Truhlar. J. Phys Chem. Lett. 2011, 3, 117-124. [91] R. Peverati, D. G. Truhlar. J. Phys Chem. Lett. 2011, 2, 2810-2817. [92] R. Peverati, D. G. Truhlar. J. Chem. Theory Comput. 2012, 8, 2310-2319. [93] R. Peverati, D. G. Truhlar. Phys. Chem. Chem. Phys. 2012, 14, 16187-16191. [94] R. Peverati, D. G. Truhlar. Phys. Chem. Chem. Phys. 2012, 14, 13171-13174. [95] S. Grimme. J. Chem. Phys. 2006, 124, 0341081-0341097. [96] Y. Zhao, B. J. Lynch, D. G. Truhlar. J. Phys. Chem. A 2004, 108, 4786-4791. [97] S. Grimme. J. Comput. Chem. 2004, 25, 1463–1473. [98] S. Grimme. J. Comput. Chem. 2006, 27, 1787–1799. [99] S. Grimme, J. Antony, S. Ehrlich, H. Krieg. J. Chem. Phys. 2010, 132, 214301/1- 214320. [100] S. Grimme, S. Ehrlich, L. Goerigk. J. Comput. Chem. 2011, 32, 1456-1465. [101] S. Grimme. WIREs Comput. Mol. Sci. 2011, 1, 211-228. Alba Campo Cacharrón 78 [102] A. Krishtal, C. Van Alsenoy, P. Geerlings. J. Chem. Phys. 2014, 140, 184105/1- 184105/14. [103] S. N. Steinmann, C. Corminboeuf. J. Chem. Theory Comput. 2010, 6, 1990-2001. [104] S. N. Steinmann, C. Corminboeuf. J. Chem. Phys. 2011, 134, 044117/1-044117/5. [105] O. A. Vydrov, T. Van Voorhis. J. Chem. Phys. 2010, 133, 244103/1-244103/9. [106] M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, B. I. Lundqvist. Phys. Rev. Lett. 2004, 92, 246401. [107] F. D. John, G. Tim. J. Phys.: Condens. Matter 2012, 24, 073201. [108] L. A. Burns, Á. V. Mayagoitia, B. G. Sumpter, C. D. Sherrill. J. Chem. Phys. 2011, 134, 084107/2-084107/25. [109] E. R. Johnson, A. D. Becke. J. Chem. Phys. 2005, 123, 024101/1-024101/7. [110] C. Van Alsenoy. J. Comput. Chem. 1988, 9, 620-626. [111] H.-J. Werner, F. R. Manby, P. J. Knowles. J. Chem. Phys. 2003, 118, 8149-8160. [112] K. Eichkorn, O. Treutler, H. Öhm, M. Häser, R. Ahlrichs. Chem. Phys. Lett. 1995, 240, 283-290. [113] F. Weigend. J. Comput. Chem. 2008, 29, 167-175. [114] C. Hattig. Phys. Chem. Chem. Phys. 2005, 7, 59-66. [115] M. S. Gordon, D. G. Fedorov, S. R. Pruitt, L. V. Slipchenko. Chem. Rev. 2011, 112, 632-672. [116] K. Morokuma. J. Chem. Phys. 1971, 55, 1236-1244. [117] W. Chen, M. S. Gordon. J. Phys. Chem. 1996, 100, 14316-14328. [118] E. D. Glendening, A. Streitwieser. J. Chem. Phys. 1994, 100, 2900-2909. [119] G. K. Schenter, E. D. Glendening. J. Phys. Chem. 1996, 100, 17152-17156. [120] E. D. Glendening. J. Phys. Chem. A 2005, 109, 11936-11940. [121] W. J. Stevens, W. H. Fink. Chem. Phys. Lett. 1987, 139, 15-22. [122] P. Su, Z. Jiang, Z. Chen, W. Wu. J. Phys. Chem. A 2014, 118, 2531-2542. [123] P. Su, H. Li. J. Chem. Phys. 2009, 131, 014102/1-014102/15. [124] T. Ziegler, A. Rauk. Theor. Chim. Acta 1977, 46, 1-10. [125] I. C. Hayes, A. J. Stone. Mol. Phys. 1984, 53, 83-105. [126] B. Jeziorski, R. Moszynski, K. Szalewicz. Chem. Rev. 1994, 94, 1887-1930. [127] R. Moszynski. Mol. Phys. 1996, 88, 741-758. [128] E. G. Hohenstein, C. D. Sherrill. J. Chem. Phys. 2010, 133, 014101/1-014101/12. [129] E. G. Hohenstein, C. D. Sherrill. J. Chem. Phys. 2010, 132, 184111/1-184111/10. [130] J. M. Turney, A. C. Simmonett, R. M. Parrish, E. G. Hohenstein, F. A. Evangelista, J. T. Fermann, B. J. Mintz, L. A. Burns, J. J. Wilke, M. L. Abrams, N. J. Russ, M. L. Leininger, C. L. Janssen, E. T. Seidl, W. D. Allen, H. F. Schaefer, R. A. King, E. F. Valeev, C. D. Sherrill, T. D. Crawford. WIREs Comput. Mol. Sci. 2012, 2, 556-565. [131] H. L. Williams, C. F. Chabalowski. J. Phys. Chem. A 2000, 105, 646-659. [132] A. Heßelmann, G. Jansen. Chem. Phys. Lett. 2002, 357, 464-470. [133] A. Heßelmann, G. Jansen. Chem. Phys. Lett. 2002, 362, 319-325. [134] A. Heßelmann, G. Jansen. Chem. Phys. Lett. 2003, 367, 778-784. [135] A. Heßelmann, G. Jansen, M. Schütz. J. Chem. Phys. 2005, 122, 014103/1- 014103/17. [136] G. Jansen. WIREs Comput. Mol. Sci. 2014, 4, 127-144. [137] A. J. Misquitta, R. Podeszwa, B. Jeziorski, K. Szalewicz. J. Chem. Phys. 2005, 123, 214103/1-214103/14. [138] R. Bukowski, R. Podeszwa, K. Szalewicz. Chem. Phys. Lett. 2005, 414, 111-116. 3. Methodology 79 [139] R. Podeszwa, R. Bukowski, K. Szalewicz. J. Chem. Theory Comput. 2006, 2, 400- 412. [140] K. Szalewicz. WIREs Comput. Mol. Sci. 2012, 2, 254-272. [141] M. Grüning, O. V. Gritsenko, S. J. A. van Gisbergen, E. J. Baerends. J. Chem. Phys. 2001, 114, 652-660. [142] D. J. Tozer, N. C. Handy. Mol. Phys. 2003, 101, 2669-2675. [143] J. Řezáč, P. Hobza. J. Chem. Theory Comput. 2011, 7, 685-689. [144] R. F. W. Bader, Atoms in Molecules: A Quantum Theory, Clarendon Press, Oxford 1990. [145] J. Contreras-Garc a, E. R. Johnson, S. Keinan, R. Chaudret, J.-P. Piquemal, D. N. Beratan, W. Yang. J. Chem. Theory Comput. 2011, 7, 625-632. [146] E. R. Johnson, S. Keinan, P. Mori-Sánchez, J. Contreras-García, A. J. Cohen, W. Yang. J. Am. Chem. Soc. 2010, 132, 6498-6506. [147] J. Tomasi, B. Mennucci, R. Cammi. Chem. Rev. 2005, 105, 2999-3094. [148] A. Klamt, G. Schuurmann. J. Chem. Soc. Perkin Trans. 2 1993, 799-805. [149] M. J. Vilkas, C.-G. Zhan. J. Chem. Phys. 2008, 129, 194109/1-194109/7. 4. Effects of microhydration on the characteristics of cation-phenol complexes 4. Phenol microhydration 83 4.1. Introduction Intermolecular interactions are known to play a key role in many aspects of chemistry and biology, being crucial to phenomena as protein···ligand interaction and molecular recognition.[1, 2] Among the different intermolecular forces, those involving aromatic moieties present especial character due to the presence of conjugated electrons and a planar geometry.[3-5] Also, as regards biochemistry, interactions involving aromatic units are crucial in protein structure. It is believed that aromatic groups interact in a different manner than aliphatic units in the side chains of amino acids, and can provide specificity in protein folding.[4-8] Cation··· interactions are strong interactions in the gas phase and have been recognized as one of the structural motifs conditioning the structure in proteins.[5, 9] Since the initial works from Dougherty et al. the cation··· interaction is now regarded as one key factor in determining the characteristics of proteins, together with hydrogen bonds, stacking interactions and salt bridges.[5-7, 9, 10] The importance of these cation··· interactions in proteins is easily understood taking into account that some amino acids as phenylalanine, tyrosine, tryptophan and histidine bear an aromatic unit in their side chains, whereas other amino acids as arginine, lysine and histidine possess cationic groups depending on the pH. Thus, interactions between side chains of these amino acids are often observed in protein structure suggesting their relevance as a stabilizing motif.[5, 9, 11] Though cation··· interactions are known to be strong interactions in the gas phase, this has not to be true in solution. Different studies give contradictory results ranging from an important contribution to protein stabilization to an almost negligible effect.[11- 22] These discrepancies are generally addressed to solvent effects, depending on the degree of exposure of the cation··· contact to the solvent. Several studies exist in literature dealing with the interaction of alkali cations with benzene, showing the usual trend of stronger interaction as the size of the cation decreases.[23-26] Also, several authors have studied the effect of water coordinated to the cation on the strength of the cation···benzene interaction.[27-32] However, most studies have been carried out with benzene as a model for aromatic interaction, with few works devoted to other aromatic units.[33-36] In this work, a study of the interaction between a K+, Na+, Li+ or Mg2+ cation with phenol in the presence of a small number of water molecules has been performed by employing ab initio and density functional theory methods. Phenol is not only a common chemical but also the chromophore of the aromatic amino acid tyrosine. Contrary to benzene, phenol possesses two different regions where a cation can establish a stabilizing interaction: the aromatic cloud and the lone pairs of the hydroxyl oxygen.[37, 38] Also, the hydroxyl group can act both as donor or acceptor in hydrogen bonds with the water molecules included in the cluster. Therefore, even when phenol is Alba Campo Cacharrón 90 optimum cation···water interaction, even at the cost of loosing part of the stabilization gained by the formation of OH···O or O-H··· contacts. With the smaller cations, there are larger differences between structures with water coordinated to the cation with respect to other structures with water hydrogen bonded to phenol. In fact, it can be observed from Table 4.2 that the energy difference between structures Phe-X1-1H2O and Phe-X2-1H2O decreases as the size of the cation increases. That is, for K+ complexes the hydroxyl group of phenol can compete with the cation for interacting with water, whereas for smaller cations the direct coordination to the cation is clearly favored. Table 4.2. Complexation energies (kcal/mol) of the most stable monohydrated cation···phenol complexes as obtained at the MP2/6-31+G(2d,p)//MP2/6-31+G(d) level. Phe-X1- 1H2O Phe-X2- 1H2O Phe-X3- 1H2O Phe-X4- 1H2O Phe-X5- 1H2O Phe-X6- 1H2O K+ -31.87 -30.24 -33.30 -33.07 -29.94 -35.72 Na+ -41.80 -36.13 -42.41 -42.35 -35.56 -43.70 Li+ -61.01 -48.99 - -62.23 -50.02 - Mg+2 -171.40 -140.45 - -174.91 -140.61 - 4.3.3. Dihydrated cation···phenol complexes With the inclusion of the second water molecule, a greater variety of structures has been located, so only the most stable among the minima located will be discussed for each cation. Generally speaking, the larger the cation, the larger the number of different minima within a given energy interval, so the behavior is simpler for Li+ and Mg2+ clusters, whereas for K+ the most complex behavior arises. Figure 4.3 shows the most stable minima located for complexes formed by phenol with K+ and Na+ cations including two water molecules. It can be observed that the first four structures have similar characteristics for both cations, whereas differences appear for the other two presented in Figure 4.3. Considering K+ complexes and the values for complexation energies listed in Table 4.3, it can be observed that structures with the K+ cation over the aromatic ring of phenol are more stable than the analogous structure with K+ interacting with the oxygen atom. The three structures named as Phe-X2-2H2O, Phe-X4-2H2O and Phe-X6-2H2O present similar complexation energies, showing that in the case of K+ complexes there is no significant energy difference when water molecule coordinates to K+ directly of interacts with the OH group of phenol, and the same can be observed in the rest of structures with K+ interacting with the hydroxyl group. 4. Phenol microhydration 91 Figure 4.2. Most stable minima of the monohydrated cation···phenol complexes. The numbers correspond to distances in Å as obtained at the MP2/6-31+G(d) level of calculation. Phe-X1-1H2O 2.685 2.298 1.870 2.002 K+2.676 Na+2.276 Li+1.882 Mg+2 1.987 Phe-X2-1H2O 2.625 2.252 1.842 1.955 1.761 1.743 1.697 1.468 Phe-X3-1H2O 2.709 2.300 2.645 2.273 2.460 2.274 2.669 2.277 1.876 1.986 2.913 2.503 1.967 2.004 Phe-X4-1H2O 2.848 2.392 1.911 1.944 1.794 1.780 1.760 1.609 Phe-X5-1H2O Phe-X6-1H2O 2.936 2.513 2.626 2.276 1.963 2.036 Alba Campo Cacharrón 92 Figure 4.3. Most stable minima located at the MP2/6-31+G(d) level for the dihydrated complexes of K+ and Na+ with phenol. Phe-X1-2H2O Phe-X2-2H2O Phe-X4-2H2OPhe-X3-2H2O Phe-K5-2H2O Phe-K6-2H2O Phe-Na5-2H2O Phe-Na6-2H2O 4. Phenol microhydration 93 The most stable minima corresponds to a trigonal arrangement of the two water molecules and the phenyl ring around the K+ cation, though Phe-X6-2H2O shows similar stability, with both water molecules interacting directly with K+, simultaneously establishing a hydrogen bond between them, keeping an O···K···O angle of only 65 degrees. However, as indicated above, several of these minima show stability differences within the error of the methods employed, so the behavior could slightly change if more sophisticated (and costly) methods are employed. Table 4.3. Complexation energies (kcal/mol) of the most stable dihydrated complexes as obtained at the MP2/6-31+G(2d,p)//MP2/6-31+G(d) level of calculation. Phe-X1-2H2O Phe-X2-2H2O Phe-X3-2H2O Phe-X4-2H2O Phe-X5-2H2O Phe-X6-2H2O K+ -46.77 -48.30 -45.30 -47.07 -46.42 -47.92 Na+ -59.46 -59.89 -55.15 -55.76 -54.36 -54.00 Li+ -83.52 -81.95 -74.81 -75.14 -75.15 -75.84 Mg+2 -218.63 -218.96 -196.73 -196.97 -198.46 -201.75 In the case of Na+ complexes, though the minima are similar, the energy ordering depends more on the arrangement of water molecules around the cation, so the two most stable structures correspond to a similar trigonal arrangement, the difference being the location of the cation over the ring or close to the hydroxyl group. The rest of the structures are less stable with relative energies more than 4 kcal/mol above the most stable minimum. It is worth noting that no stable minimum was found similar to those observed for K+ cation, with water molecules bound to the cation and interacting between themselves. In the case of Na+, the interaction with water is stronger, so the distorted hydrogen bond between water molecules cannot overcome the energy loss when departing from the trigonal arrangement. For K+ complexes, minima similar to those obtained for Na+ can be found, with complexation energies amounting to -44.9 and -43.5 kcal/mol. Li+ and Mg2+ complexes shown in Figure 4.4 present similar characteristics to those found for Na+. For both cations the same minima are obtained, with energies listed in Table 4.3. The most stable structures correspond to a trigonal arrangement around the cation, with no large differences regarding whether the cation is on the ring or over the hydroxyl group. There is a significant energy gap with the rest of structures, reaching 6 kcal/mol and 17 kcal/mol for Li+ and Mg2+, respectively. It is also worth noting that in Mg2+ complexes, coordination with a water molecule in the second shell is energetically favored over coordination to the hydroxyl group of phenol, contrary to the behavior observed with the other cations, where the hydroxyl group can be competitive for coordinating water molecules compared with other water units already present in the complex. Alba Campo Cacharrón 94 Figure 4.4. Most stable minima located at the MP2/6-31+G(d) level for the dihydrated complexes of Li+ and Mg2+ with phenol. Phe-X1-2H2O Phe-X2-2H2O Phe-X5-2H2O Phe-X6-2H2O Phe-X3-2H2O Phe-X4-2H2O 4. Phenol microhydration 95 Figure 4.5. Most stable minima located for the trihydrated complexes of K+ with phenol. Numbers correspond to complexation energies in kcal/mol at the MP2/6- 31+G(2d,p)//MP2/6-31+G(d) level. Phe-K2-3H2O Phe-K1-3H2O -59.22 -59.08 Phe-K4-3H2O Phe-K3-3H2O -59.75 -58.83 Phe-K6-3H2OPhe-K5-3H2O -60.04 -60.28 Alba Campo Cacharrón 96 Figure 4.6. Most stable minima located for the trihydrated complexes of Na+, Li+ and Mg2+ with phenol. Numbers correspond to complexation energies in kcal/mol at the MP2/6-31+G(2d,p)//MP2/6-31+G(d) level. Phe-Li1-3H2OPhe-Na1-3H2O -100.54-73.63 Phe-Mg1-3H2O -257.98 Phe-Li2-3H2OPhe-Na2-3H2O -99.28-73.76 Phe-Mg2-3H2O -254.73 Phe-Li3-3H2OPhe-Na3-3H2O -98.11-72.02 Phe-Mg3-3H2O -245.94 Phe-Li4-3H2OPhe-Na4-3H2O -97.64-72.17 Phe-Mg4-3H2O -244.70 4. Phenol microhydration 97 4.3.4. Trihydrated complexes The incorporation of a third water molecule increases the complexity of the potential energy surface of the complexes giving rise to a great variety of minima, usually with similar stabilities, especially for the larger cations. Therefore, K+ complexes present the most complex behavior with the largest number of minima within a small energy interval. The most stable minima found are shown in Figure 4.5, differing in stability by less than 1 kcal/mol. Therefore, as water molecules are included in the complex, there are smaller differences between direct coordination to the K+ cation, to another water molecule or to the hydroxyl group, thus allowing for a greater variety of structures. The most stable complexes formed with Na+, Li+ and Mg2+ are shown in Figure 4.6. In Na+ and Li+ clusters, the two most stable structures show very different structural arrangements. Whereas Phe-X1-3H2O corresponds to the tetrahedral arrangement of water and phenol around the cation, in Phe-X2-3H2O the cation is surrounded by the three water molecules and does not interact directly with phenol, so the interaction takes place between a hydrated cation and the phenol moiety. Finally, in Mg2+ complexes the tetrahedral arrangements are the most stable, the rest of the minima found being less stable by more than 10 kcal/mol. 4.3.5. Tetrahydrated clusters The inclusion of a fourth water molecule complicates even more the exploration of the potential energy surface of the clusters. Therefore, a partial exploration starting from the most stable complexes found for trihydrated clusters was performed. The fourth most stable minima for each of the cations found in this work are shown in Figure 4.7 together with their complexation energies. These minima are similar to those already found by Vaden and Lisy [37] though since the levels of calculation are different, differences arise. For example, the minima in Figure 4.7 for K+ are within 1.5 kcal/mol, whereas those reported by Vaden and Lisy span over 5 kcal/mol. Though complexation energies become more negative as the polarizing power of the cation increases as before, the energy differences among structures for a given cation are almost negligible. No stable structure with the five oxygen atoms surrounding the cation was located, so most complexes present a tetrahedral arrangement of oxygen atoms around the cation, whereas the fifth oxygen atom interacts in a second hydration shell. Therefore, all tetrahydrated complexes exhibit interactions among water molecules or between water and phenol. In summary, it can be observed that as more water molecules are included, there are fewer differences between the different kinds of union: water···water, water···phenol or water···cation. For the smallest cations this implies that only a couple of structures present similar stability whereas the rest are significantly less stable. However, for Na+ complexes and, especially for K+ complexes, Alba Campo Cacharrón 98 there can be more stable structures within a small energy interval so there can be contributions for several conformers to the properties of the system, as already suggested by the infrared spectra to be discussed in the following section. If the variations on complexation energy upon incorporation of a water molecule are considered (Table 4.4), the following behaviors can be observed. For K+ complexes, the addition of the first water molecule stabilizes the complex in a larger quantity than the complexation of K+ to phenol. This is due to the strong K+···H2O interaction, but also to the formation of the O-H···O hydrogen bond. In Na+ complexes the stability gain associated to this first water molecule is slightly smaller than the cation··· interaction. For Li+ and Mg2+, even though the interaction with water is stronger, the complex is stabilized by a considerably smaller quantity than the original cation··· interaction. As more water molecules are included in the complex the stability gain decreases steadily, though for K+ complexes the energy changes upon inclusion of the second and third water molecules are equivalent, suggesting the greater capacity of potassium for accommodating water molecules and also the extra stabilization due to the formation of hydrogen bonds. For all cations, the inclusion of the fourth water molecule is accompanied by a stability gain roughly half of that obtained with the first water unit. Table 4.4. Complexation energy changes (kcal/mol) for the incorporation of one water molecule to a complex as obtained at the MP2/6-31+G(2d,p)//MP2/6-31+G(d) level of calculation. Water molecules n K+ Na+ Li+ Mg2+ 0 -17.19 -22.61 -35.99 -116.4 1 -18.56 -21.09 -26.21 -58.51 2 -12.55 -16.19 -21.32 -43.72 3 -11.98 -13.87 -17.02 -39.35 4 -10.29 -12.24 -13.36 -23.97 4. Phenol microhydration 99 Figure 4.7. Most stable minima located for tetrahydrated complexes. Numbers correspond to complexation energies in kcal/mol at the MP2/6-31+G(2d,p)//MP2/6- 31+G(d) level. Phe-Mg1-4H2O Phe-Mg2-4H2O Phe-Mg3-4H2O Phe-Mg4-4H2O Phe-K3-4H2O Phe-K4-4H2O Phe-Li1-4H2O Phe-Li2-4H2O Phe-Li3-4H2O Phe-Li4-4H2O Phe-Na1-4H2O Phe-Na2-4H2O Phe-Na3-4H2O Phe-Na4-4H2O -71.83 -71.18 -70.55 -86.00 -85.37 -85.79 -85.97 -112.55 -112.40 -113.53 -113.90 -278.95 -281.12 -281.48 -281.95 Phe-K2-4H2O Phe-K1-4H2O -71.54 Alba Campo Cacharrón 106 [58] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. J. A. Montgomery, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, D. J. Fox, Gaussian 09, revision A.02, 2009. [59] M. S. Marshall, R. P. Steele, K. S. Thanthiriwatte, C. D. Sherrill. J. Phys. Chem. A 2009, 113, 13628-13632. 5. Interaction of aromatic units of amino acids with guanidinium cation. The interplay of π···π, X- H···π and M+···π contacts 5. Guanidinium aromatic trimers 109 5.1. Introduction Non-covalent interactions play a key role in many areas of modern chemistry, especially in the field of supramolecular chemistry and molecular recognition, as well as in biochemistry.[1, 2] Among the different kinds of non-covalent contacts relevant in biological systems, a special place is occupied by interactions with participation of aromatic units.[3-6] It is believed that aromatic groups interact in a different manner than aliphatic units in the side chains of amino acids, and can provide specificity in protein folding. More specifically, three kinds of non-covalent interactions involving π systems are usually considered: ion···π, XH···π and π···π contacts, which are attractive interactions that can affect significantly to the behavior of a given system.[3-5] π···π interactions, like in benzene dimer, are usually governed by dispersion effects and play an essential role in the folding of proteins and in the structure of DNA as well as in its interactions with small molecules.[3, 5] Hydrogen bonding interactions to the aromatic cloud (XH···π) are also important as to determine the characteristics of many systems.[4, 7, 8] These hydrogen bonds are usually weaker than the typical OH···O ones, and usually exhibit a larger dispersive nature. These two kinds of π interactions are normally weak, though the combined effect of many of them can have a deep impact on the characteristics of the system. On the other hand, ion···π interactions are strong interactions in the gas phase,[4, 9-11] usually dominated by electrostatic and polarization terms as a consequence of the presence of the ion and a polarizable π cloud.[12, 13] Cation···π interactions are recognized to be an important factor in ion selectivity in potassium channels, and their importance has been demonstrated in neurotransmitter receptors.[4, 11, 14, 15] The importance of the cation···π interaction in proteins is easily understood taking into account that some amino acids such as phenylalanine, tyrosine, tryptophan, and histidine bear an aromatic unit in their side chains, whereas other amino acids as arginine, lysine, and histidine possess cationic groups depending on the pH.[6, 16] Thus, interactions between side chains of these amino acids are often observed in protein structure suggesting their relevance as a stabilizing motif. Though the cation··· interaction is usually stronger than π···π or XH···π contacts in the gas phase, as corresponds to the interaction of a bare cation with a polarizable and electron-rich aromatic cloud, the environment can significantly alter its characteristics.[11, 17] Several studies have shown how the coordination of electron-rich solvent molecules to the cation decreases the intensity of the interaction with the cation, as it would be expected taking into account the decrease of the effective charge carried by the cation due to the solvent molecules.[18-22] In any case, the behavior and stability of a system, even a simple one, will be usually the result of the interplay of different non-covalent interactions. If aromatic molecules are present it is probable for all π···π, XH···π and ion···π to play a role. Alba Campo Cacharrón 110 In the present work, the interplay among these intermolecular contacts is analysed in ternary systems containing one cationic unit and two equal or different aromatic units. Trying to represent possible contacts among amino acid side chains the aromatic molecules considered have been benzene, phenol and indole, as models for phenylalanine, tyrosine and tryptophan, respectively. Guanidinium C(NH2)3+ has been chosen as cation since it is part of the side chain or arginine. Guanidinium is widely used as a denaturant since it is believed to interact with the side chains of proteins, as well as a basic unit for constructing anion selective receptors or ionic liquids.[23, 24] Besides, due to its special planar structure guanidinium cation is more prone than simpler cations to present parallel-stacked structures which can be of importance in these systems.[20, 25] In fact, stacking has been observed between solvated guanidinium cations in water solution.[26, 27] An interesting phenomenon when dealing with intermolecular interactions is the possibility of cooperative or anticooperative effects in systems with more than two species.[28, 29] These phenomena are usually weak though they can be of importance as to characterize the behavior of the system, as indicated by several studies recently devoted to the task of evaluating these effects.[12, 30-34] Usually, cooperative effects are associated to polarization, so the combination of a cation and a polarizable aromatic cloud, like in the trimers object of this study, is indicative that cooperativity, a priori, could be significant in determining the characteristics of these systems.[12] In summary, the present work intends to shed light on the characteristics on the interaction between guanidinium and aromatic moieties of amino acids. The results obtained would help to rationalize the results observed in multiple cation···π interactions in proteins and their mutual influence. 5.2. Methods As commented above, systems consisting on a guanidinium cation and two aromatic molecules among benzene, phenol and indole have been considered in this study. These complexes have been fully optimized at the M06-2X/6-31+G* level of calculation[35] starting from a variety of different structures. Starting structures were constructed from prototypes observed in the guanidinium-aromatic and aromatic-aromatic dimers. Therefore, the parallel, displaced-parallel and T-shaped structures for benzene dimer have been employed, and analogous structures have also been used for phenol and indole.[5, 36] T-shaped and parallel orientation of the guanidinium cation with respect to the aromatic clouds have been considered,[20, 37] with the cation located between the aromatic molecules or bounded to only one of them. Other possible structures corresponding to hydrogen-bonded clusters have also been tested following chemical knowledge (i.e. O-H···O in phenol dimer, N-H···π in indole-containing clusters, etc.). After 5. Guanidinium aromatic trimers 111 one stationary point is located a frequency calculation has been carried out in order to ensure that the structure corresponds to a minimum. For the minima the complexation energy has been obtained by applying the counterpoise procedure with a variety of methods.[38, 39] Therefore, the complexation energy of the complexes has been obtained as:     i complex i complex i i isolated i complex complex iEijkEiEijkEE )(...)()(...)( (eq. 5.1) where superscripts refer to the geometry employed, subscripts to the fragment considered and terms in parentheses to the basis set employed in the calculation. Applying this procedure, M06-2X/6-31+G* complexation energies have been obtained. Also, complexation energies have been obtained with the MP2 method. Single point calculations have been carried out at the M06-2X/6-31+G* optimized geometries employing MP2 with the aug-cc-pVDZ and aug-cc-pVTZ basis sets. With these results, a two-point extrapolation of the correlation energy to the basis set limit has been performed applying[40] 3; )1( )1( )1( )1( 2, 33 3 2, 33 3 2,       XE XX X E XX X EZXpV MPcorr pVXZ MPcorr CBS MPcorr . (eq. 5.2) The MP2 complexation energy to basis limit is the estimated as CBS MPcorr pVTZ HF CBS MP EEE 2,2  . (eq. 5.3) However, it is well know that MP2 tends to overestimate the magnitude of the interaction when aromatic molecules are involved, especially when oriented in parallel.[5, 41] A variety of empirical procedures have been proposed to correct for this deficiency, mostly based in a different scaling of the contributions of parallel and antiparallel electrons to correlation. The first proposal of this kind is the Spin Component Scaled MP2 (SCS-MP2) proposed by Grimme, where opposite-spin and same-spin contributions are scaled by 1.20 and 0.33, respectively.[42] Therefore SCS-MP2 complexation energies have been obtained with the aug-cc-pVDZ and aug-cc-pVTZ basis sets, and were also extrapolated to basis set limit (SCS-MP2/CBS). In order to analyze the balance of the interactions between fragments in the trimer, pair energy contributions have also been calculated. Therefore, the interaction energy is expressed as   ij bodyij EEE 3 (eq. 5.4) where the ∆Eij are the interaction energies for each pair formed in the trimer as computed employing the whole basis set and the optimized geometry of the trimer.[28, Alba Campo Cacharrón 112 29] The difference between the summation of pair energies and the interaction energy of the trimer is the contribution from 3-body effects. This partitioning has been performed with the methods commented above. Finally, in order to have more insight into the nature of the interaction and the balance of different contributions, an energy partitioning scheme has been applied. Therefore, the interaction energies have been decomposed in electrostatic, repulsion (exchange + repulsion), polarization and dispersion components by applying the Local Molecular Orbital-Energy Decomposition Analysis (LMO-EDA).[43] The partitioning has been performed with both the M06-2X and MP2 methods. Optimizations and frequency calculations have been performed with Gaussian09;[44] LMO-EDA calculations were performed with GAMESS,[45, 46] and Turbomole has been employed in the MP2 calculations.[47] In order to save computation time, the resolution of the identity has been applied in MP2 calculations, both to the correlation calculation (RI-MP2) and to the HF one (RI-JK), employing suitable auxiliary basis sets as provided in Turbomole. Therefore the corresponding aug-cc-pVXZ fitting basis sets have been employed in the correlation part whereas def2-TZVPP has been used for RI-JK calculations.[47-49] 5.3. Results Minimum energy structures and complexation energies will be presented first for complexes formed by guanidinium and two equal aromatic units, followed by results obtained for complexes formed by guanidinium cation and two different aromatic molecules (mixed complexes). 5.3.1. Complexes with the same aromatic molecules Figure 5.1 shows the structures of the most stable minima located for the complexes formed by guanidinium cation and two benzene molecules, whereas Table 5.1 lists the values obtained for their complexation energies with different levels of calculation. Four different minima have been found for these complexes. Bz-Bz-1 corresponds to a structure with stacked aromatic rings and the guanidinium cation interacting with one of them, adopting a perpendicular arrangement which has been shown to be the most stable minimum in guanidinium···Bz complexes.[20, 37] Two N-H···π contacts are observed: one at 2.3 Å and a longer one at 2.9 Å. Bz-Bz-2 and Bz-Bz-3 correspond to structures where benzene rings display a perpendicular arrangement (T-shaped) with C-H···π contacts at around 2.5-2.6 Å to the center of the other aromatic ring. These arrangements roughly correspond to the typical structures found for benzene dimer (parallel displaced and C-H···).[5] Finally, Bz-Bz-4 corresponds to a doubly T-shaped structure, where the aromatic rings are far apart, both of them coordinating the 5. Guanidinium aromatic trimers 113 guanidinium cation in perpendicular arrangements, forming hydrogen bonds at 2.3 Å of the center of the ring. Table 5.1 shows the complexation energies of Bz-Bz complexes as obtained with a variety of methods. Taking into account the values in Table 5.1 the results obtained at the M06-2X/6-31+G* and at the SCS-MP2/CBS levels of calculation are pretty similar, whereas MP2 gives in all cases overestimated complexation energies, especially with the larger basis set. Therefore, in the following we will mainly discuss values obtained with these two levels of calculation. Table 5.1. Complexation energies (kcal/mol) obtained for complexes formed by guanidinium and two equal aromatic molecules as obtained with different methods, all employing the optimized M06-2X/6-31+G* geometry. 631+G* aug-cc-pVDZ CBS M062X MP2 SCS-MP2 MP2 SCS-MP2 Bz-Bz-1 -18.55 -21.03 -16.97 -22.79 -18.54 Bz-Bz-2 -21.03 -22.44 -18.58 -24.30 -20.27 Bz-Bz-3 -20.58 -21.41 -17.20 -23.25 -18.85 Bz-Bz-4 -25.62 -27.12 -23.52 -29.25 -25.48 Ph-Ph-1 -33.19 -31.93 -27.83 -34.05 -29.71 Ph-Ph-2 -31.21 -30.47 -25.76 -32.80 -27.85 Ph-Ph-3 -29.47 -29.93 -26.34 -31.88 -28.10 Ph-Ph-4 -28.73 -27.31 -22.36 -29.58 -24.36 Ph-Ph-5 -32.20 -31.05 -26.76 -33.20 -28.67 Ph-Ph-6 -32.21 -32.71 -29.05 -34.67 -30.81 Ph-Ph-7 -32.81 -30.63 -27.90 -32.17 -29.35 In-In-1 -29.61 -34.47 -27.43 -36.87 -29.53 In-In-2 -32.72 -36.01 -30.25 -38.37 -32.36 In-In-3 -31.54 -34.70 -28.94 -37.04 -31.04 In-In-4 -33.14 -36.39 -30.55 -38.79 -32.70 In-In-5 -38.16 -40.04 -34.70 -42.69 -37.09 Alba Campo Cacharrón 114 Figure 5.1. Most stable minima found for complexes of guanidinium with two benzene molecules as obtained at the M06-2X/6-31+G* level. Selected distances in Å. Bz-Bz-2Bz-Bz-1 Bz-Bz-3 Bz-Bz-4 2.262 2.911 2.172 3.107 2.331 2.632 3.652 2.293 2.292 2.930 2.962 5. Guanidinium aromatic trimers 115 As expected, Bz-Bz-4, with a doubly T-shaped structure is the most stable complex with a complexation energy of -25.5 kcal/mol because guanidinium can coordinate at the same time the two benzene rings, establishing two simultaneous cation··· interactions. Bz-Bz-2, with one T-shaped contact, is about 5 kcal/mol less stable (-20.3 kcal/mol). Here, one of the benzene rings coordinates to the guanidinium cation through a cation··· interaction and by means of a CH··· contact to the other benzene ring. Bz- Bz-1 corresponds to a parallel-displaced benzene dimer coordinated to guanidinium by one of the phenyl rings. The lack of direct contact between guanidinium and one of the benzene moieties makes this structure less stable, reaching only -18.5 kcal/mol. Finally, in Bz-Bz-3 the benzene ring that does not interact directly with the guanidinium cation moves away from the perpendicular arrangement. This structure could be considered as a combination of two of the minima of guanidinium···benzene complexes; one T-shaped and one parallel displaced with guanidinium parallel to the ring.[20, 37] This is the only structure which shows significant differences between M06-2X (-20.6 kcal/mol) and SCSMP2/CBS (-18.9 kcal/mol) results. Contrary to benzene, phenol possesses two different regions where a cation can establish a stabilizing interaction: the aromatic cloud and the lone pairs of the hydroxyl oxygen. Therefore, even when phenol is very similar to benzene, the presence of the hydroxyl group introduces a greater complexity on the potential energy surface of the clusters, allowing for a greater variety of stable structures. Figure 5.2 shows the structures of the most stable minima located for the complexes formed by guanidinium cation and two phenol molecules. It can be observed that phenol complexes cannot be classified into prototypical structures as easily as benzene ones, since due to the presence of the hydroxyl group there is a tendency to form O-H···O hydrogen bond contacts. Ph-Ph-6 and Ph-Ph-7 correspond to doubly T-shaped structures. In Ph-Ph-6 guanidinium cation coordinates the phenyl ring and the hydroxyl group whereas in Ph- Ph-7 it is coordinated only to the hydroxyl groups. It is worth noting that similar structures but with different disposition of the rings have been located, though only the most stable ones are included in Figure 5.2. In most cases the structures show O-H···O short hydrogen bonds at 1.7-1.9 Å, whereas guanidinium is also hydrogen-bonded to the hydroxyl group by means of N-H···O hydrogen bonds at around 1.8-1.9 Å. The contacts between guanidinium and the phenyl ring show similar behavior as that observed in benzene complexes. No typical stacked structure has been found in complexes with phenol because the hydroxyl group tends to interact establishing OH···π or OH···O hydrogen bonds as it can be observed in all minima in Figure 5.2 except the doubly T-shaped ones. As regards T- shaped structures, only Ph-Ph-3 could be considered within this arrangement because it is the only one that presents both aromatic rings interacting in a perpendicular disposition, even though the interaction takes place by means of a O-H···O contact. Alba Campo Cacharrón 122 Figure 5.4. Selected minima for complexes formed by guanidinium, benzene and phenol as obtained at the M06-2X/6-31+G* level. Selected Distances in Å. Bz-Ph-DT 2.243 2.950 1.982 1.964 2.191 2.067 2.321 1.912 1.918 2.334 2.270 Ph-Bz-PBz-Ph-P Bz-Ph-T Ph-Bz-T 5. Guanidinium aromatic trimers 123 Figure 5.5. Selected minima for complexes formed by guanidinium, benzene and indole as obtained at the M06-2X/6-31+G* level. Selected Distances in Å. Bz-In-DT 2.264 2.894 2.253 2.203 2.173 2.518 2.241 2.192 2.241 2.965 3.056 2.249 2.241 In-Bz-PBz-In-P Bz-In-T In-Bz-T Alba Campo Cacharrón 124 Figure 5.6. Selected minima for complexes formed by guanidinium, phenol and indole as obtained at the M06-2X/6-31+G* level. Selected Distances in Å. 1.962 2.698 2.408 1.956 2.207 2.205 2.252 1.849 2.239 2.436 2.285 2.230 2.306 2.405 2.314 2.292 Ph-In-PPh-In-P Ph-In-T Ph-In-T Ph-In-DT 5. Guanidinium aromatic trimers 125 Finally, for complexes containing the cation plus indole and phenol (Figure 5.6) the behavior is similar overall. For the doubly T-shaped minimum the complexation energy is halfway between those of complexes with only one kind of aromatic molecule. However, it can be observed that phenol coordination is preferred over indole for stacked and T-shaped structures. This is because in the presence of phenol, In-Ph structures correspond to our prototypes, but for Ph-In complexes there are clear deviations. Therefore whereas In-Ph-P is a typical stacked structure with a complexation energy of -27 kcal/mol, Ph-In-P minimum is not stacked, and guanidinium interacts with both aromatic units leading to a complexation energy of -31.7 kcal/mol. In fact this minimum corresponds more closely to a guanidinium coordinated to the hydroxyl group which establishes a O-H···π hydrogen bond to the pyrrol ring of indole. T-shaped minima both present hydrogen bonds. Ph-In-T forms a O-H···π hydrogen bond to the pyrrol ring of indole, whereas a N-H···π hydrogen bond is observed in In-Ph-T, this latter structure being disfavored by about 2 kcal/mol. So, these mixed complexes behave halfway the complexes formed with only one kind of aromatic molecule. This is especially evident in the case of doubly T-shaped structures, which always show complexation energies almost midway the complexes with only one kind of aromatic molecule. These doubly T-shaped structures are always the most stable found among the clusters studied. A similar behavior is observed in stacked clusters of benzene and indole, but the presence of phenol and the tendency of its hydroxyl group to form hydrogen bonds introduce other possibilities for the interaction. In the absence of extra interactions with the guanidinium cation or hydrogen bonds, stacked structures are the least stable. T-shaped minima are usually the second most stable, but again in phenol complexes the formation of hydrogen bonds can alter the order of stability, favoring coordination of phenol by the guanidinium cation. In the absence of these effects guanidinium coordinates preferentially to indole over phenol and over benzene. A question arises about whether the kind of structures already discussed is present in proteins. It is known that stacking interactions between aromatic side chains are frequent, as it is the case with the cation···aromatic interactions.[4, 11] It is expected, however that the motifs considered in this work, simultaneously involving three different side chains will be less frequent. However, searching into the Protein Data Bank Europe it is possible to find the kind of structures considered in this work.[52] As example, searching for a pattern formed by a nitrogen atom of guanidinium in arginine contacting with indole ring in tryptophan, which in turn is stacked to phenyl group of phenylalanine (like in In-Bz-P) results in more than forty coincidences. Among these structures there are several of them that resemble the patterns shown in the manuscript. Of course there are geometrical differences that are mainly associated to the more complex environment in the protein, and to the rest of the side chains hindering the free orientation of the aromatic rings and the cationic fragment to interact Alba Campo Cacharrón 126 in an optimal way. In any case, many of the structures found resemble the minima considered in this work, so it can be expected that their characteristics would be helpful in order to understand these kind of contact in proteins. Also, gas phase results would help to isolate other factors like solvent effects of other groups nearby which can alter the mutual arrangement of aromatic and cationic side chains. 5.3.3. Pair energies In order to understand more closely the behavior of these ternary systems, a decomposition on pair contributions has been performed for the clusters discussed above. Table 5.3 shows the results obtained for the pair energies of complexes as obtained at the M06-2X/6-31+G* level of calculation, together with the contributions from three-body effects. It can be observed that in all complexes the interaction of guanidinium with benzene, phenol and indole amounts to -14 kcal/mol, -18 kcal/mol and -21 kcal/mol, respectively. The strength of these interactions is almost the same in all minima considered, thus indicating that the presence of a second aromatic unit hardly affects the guanidinium···π interaction. For example guanidinium···benzene interaction varies between -12.8 and -13.9 kcal/mol, and similar variations are observed for the other aromatic species. Part of these changes can be attributed to changes in the geometry of the cation···π interaction depending on the complex considered. In any case, the effect of the second aromatic unit on the guanidinium···π interaction is described by the 3-body contribution to the interaction energy, which in most cases is not large. In doubly T-shaped structures there is a similar pattern for the pair interactions in all cases. As expected, there are two strong interactions due to guanidinium···π interactions, plus one weak repulsive interaction due to the interaction of the two aromatic molecules located far apart. For the doubly T-shaped minima the contribution of the 3-body term is destabilizing, amounting to between 1.3 kcal/mol for Bz-Bz complex to 2.7 kcal/mol in In-In complex. This is a consequence of the charge of the cation being shared by the two aromatic molecules. Mixed complexes show a similar behavior, but the two cation···π contacts are not equivalent showing a pair of values typical for the two aromatic units considered. T-shaped minima typically show one strong interaction due to the cation interacting with the closer aromatic molecule, plus two similar weak attractive interactions between the cation and the other aromatic molecule, and between both aromatic molecules. Therefore, in Bz-Bz-2 there is a guanidinium···benzene contact amounting to around -13 kcal/mol, plus an interaction of -6 kcal/mol between guanidinium and the other benzene molecule far apart. Finally, both benzene molecules contribute with -2 kcal/mol to the interaction due to a C-H···π contact. As in the case of doubly T-shaped minima, 3-body effects are repulsive and small. 5. Guanidinium aromatic trimers 127 Alba Campo Cacharrón 128 Figure 5.7. Pair energy contributions in complexes with benzene and indole as obtained at the M06-2X/6-31+G* level. C is always guanidinium cation. In the X-Y pair, A corresponds to X and B to Y. A similar behavior is observed in mixed complexes where benzene is coordinated to the cation. In the case of indole and phenol complexes there is also a strong cation···π interaction together with two weaker contributions from the other molecule pairs. However, in these complexes the interaction between the two aromatic molecules occurs by means of a hydrogen bond, so the stabilization increases with respect to benzene complexes. Thus, in the T-shaped minima containing two indole molecules, the N-H···π hydrogen bond contributes with around -8 kcal/mol to the stability of the complex. Most complexes with phenol or indole coordinated to guanidinium show similar values (though somewhat weaker, around -4 to -5 kcal/mol). Comparing the values obtained for the mixed Ph-In complexes, it can be appreciated how the O-H···π contact (-7.3 kcal/mol) is stronger than the N-H···π one (-5.7 kcal/mol). Also, in these structures, contrary to benzene ones, 3-body effects are attractive, amounting between -1.5 to -2.7 kcal/mol, these values being a consequence of the presence of the hydrogen bonds. Parallel minima show similar energy decomposition as T-shaped structures, with a strong cation···π contact plus two weak contacts from the other pairs. When phenol and indole are involved, the interaction between aromatic molecules can reach around - 5 kcal/mol, being only slightly weaker than the interaction in hydrogen-bonded structures present in T-shaped minima. Three body effects are mostly attractive but weaker than those observed in the other structural arrangements. The different pair contributions can be easily seen in Figure 5.7, where results are shown for benzene and indole-containing complexes, showing the patterns discussed above. -45 -35 -25 -15 -5 5 Bz-Bz Bz-In In-Bz In-In Bz-Bz Bz-In In-Bz In-In Bz-Bz Bz-In In-In Energy (kcal mol-1) AB AC BC 3-body Parallel T-Shaped DoublyT-Shaped Energy (kcal/mol) 5. Guanidinium aromatic trimers 129 5.3.4. LMO-EDA analysis LMO-EDA decomposition has been performed with both M06-2X and MP2 methods. It has to be taken into account that in the case of post-HF methods, the so-called dispersion contribution corresponds to the contribution of the correlation to the interaction energy. With DFT methods, on the other hand, dispersion (obtained as a difference between a sum of terms and the interaction energy from a supermolecule calculation) also recovers deviations from other contributions. Along this section only results obtained with the MP2/aug-cc-pVDZ partitioning will be considered, the M06-2X ones presented in Appendix B. The largest differences between both methods have been observed in the repulsion term, which is larger with M06-2X, thus leading as compensation to too large dispersion contributions. Figures 5.8 and 5.9 show the results obtained for this decomposition for the structures presented in Figures 5.1 to 5.6. It can be observed that in Bz-Bz complexes the repulsion is almost constant and provides a kind of background against which the attractive contributions of electrostatics, polarization and dispersion bring close together the fragments in the complex. In the doubly T-shaped minimum Bz-Bz-4 the main contribution to the stability of the complex is electrostatic, amounting to -21 kcal/mol, as a consequence of the interaction of the cation with both aromatic clouds of the benzene molecules. On the other hand, dispersion contributes with almost -13 kcal/mol to the stability of the complex. The polarization due to the cation also contributes significantly to the stability, reaching -15 kcal/mol. In the parallel structure Bz-Bz-1 the electrostatic contribution drops to -17 kcal/mol, and consequently, the contribution of polarization also decreases to -11 kcal/mol. These joined contributions introduce a difference of about 8-9 kcal/mol with respect to the doubly T-shaped one. This is partly recovered by an increase of dispersion contribution due to the stacking of the rings. As regards T-shaped structures the distribution of contributions is pretty similar to the parallel ones. Therefore, the preference of doubly T-shaped minima is a consequence of larger electrostatic and polarization contributions, whereas other structures are more favored by dispersion. In Ph-Ph clusters, the electrostatic contribution is large in all minima, even in parallel or T-shaped ones. This is related to the formation of hydrogen bonds in most of the minima. In Ph-Ph-1 and Ph-Ph-5 the electrostatic contribution is the largest, because in these minima guanidinium is coordinated to the hydroxyl group, which simultaneously forms a O-H···O (Ph-Ph-1) or O-H···π (Ph-Ph-5) hydrogen bond. The chain-like arrangement of the contacts produces large electrostatic and polarization contributions. In-In complexes resemble the behavior of benzene aggregates. The doubly T-shaped minimum presents large electrostatic, polarization and dispersion contributions, combined with the smallest repulsion term. On the other hand the parallel structure In- In-1 shows much smaller electrostatic and polarization terms, while dispersion is larger, partly cancelled out with the increase in repulsion. Alba Campo Cacharrón 130 Figure 5.8. LMO-EDA decomposition of complexes with the same aromatic molecules at the MP2/aug-cc-pVDZ level. -16.6 -17.0 -17.5 -21.2 22.7 21.4 23.3 22.2 -11.0 -12.1 -10.8 -15.4 -16.1 -14.5 -16.6 -12.6 -80 -60 -40 -20 0 20 40 Bz-Bz-1 Bz-Bz-2 Bz-Bz-3 Bz-Bz-4 Energy (kcal mol-1) Electrostatic Repulsion Polarization Dispersion -37.4 -35.5 -34.0 -30.7 -35.8 -31.7 -34.1 35.9 37.5 34.6 32.9 34.6 28.3 38.9 -17.3 -17.4 -20.4 -14.0 -15.8 -17.8 -16.2 -14.1 -16.6 -11.5 -18.1 -15.6 -12.2 -21.1 -80 -60 -40 -20 0 20 40 Ph-Ph-1 Ph-Ph-2 Ph-Ph-3 Ph-Ph-4 Ph-Ph-5 Ph-Ph-6 Ph-Ph-7 Energy (kcal mol-1) Electrostatic Repulsion Polarization Dispersion -25.8 -27.4 -27.8 -28.0 -30.9 35.4 31.2 31.4 31.3 30.1 -15.3 -17.1 -18.5 -18.3 -19.6 -28.9 -21.9 -21.0 -21.4 -19.4 -80 -60 -40 -20 0 20 40 In-In-1 In-In-2 In-In-3 In-In-4 In-In-5 Energy (kcal mol-1) Electrostatic Repulsion Polarization Dispersion Energy(kcal/mol) Energy(kcal/mol)Energy(kcal/mol) 5. Guanidinium aromatic trimers 131 Figure 5.9. LMO-EDA decomposition of complexes with different aromatic molecules at the MP2/aug-cc-pVDZ level. -19.0 -22.2 -20.3 -23.0 -38.8 -24.2 26.3 25.7 29.3 29.6 38.6 30.7 -11.7 -11.3 -12.6 -13.9 -17.2 -14.5 -18.5 -14.5 -22.1 -23.0 -19.4 -23.0 -80 -60 -40 -20 0 20 40 Bz-Ph-P Ph-Bz-P Bz-In-P In-Bz-P Ph-In-P In-Ph-P Energy (kcal mol-1) Electrostatic Repulsion Polarization Dispersion -29.3 -25.7 -20.5 -23.8 -30.1 -25.8 29.1 28.0 24.8 26.6 33.6 29.4 -14.3 -16.5 -13.3 -16.8 -18.9 -15.8 -16.2 -13.9 -18.5 -17.3 -18.4 -19.6 -80 -60 -40 -20 0 20 40 Bz-Ph-T Ph-Bz-T Bz-In-T In-Bz-T Ph-In-T In-Ph-T Energy (kcal mol-1) Electrostatic Repulsion Polarization Dispersion -26.5 -25.4 -30.6 25.3 26.9 28.0 -16.7 -16.5 -18.5 -12.5 -17.7 -15.3 -80 -60 -40 -20 0 20 40 Bz-Ph-DT Bz-In-DT Ph-In-DT Energy (kcal mol-1) Electrostatic Repulsion Polarization Dispersion Energy(kcal/mol) Energy(kcal/mol) Energy(kcal/mol) Alba Campo Cacharrón 234 Figure 9.6. NCI plots for the complexes formed by SumaCN and the anions studied. Only complexes by the I line and the most stable orientation are shown. The product of the density times the sign of the second eigenvalue of its hessian is mapped onto an isosurface of reduced density gradient with value 0.5 a.u. The color scale goes from - 0.015 a.u. (blue) to 0.015 a.u. (red). Cl- BF4- NO3- Br- ClO4- CO2H- 9. Anion’s Nature and solvent effects 235 9.3.4. NCI Analysis Figure 9.6 shows the non-covalent interaction NCI plots for complexes formed with SumaCN. NCI maps the product of the sign of the second eigenvalue of the hessian of the density times the density onto an isosurface of reduced density gradient, allowing a graphical display of the most relevant stabilizing (blue) and destabilizing (red) interaction in a system.[49, 50] Considering the plots in Figure 9.6, a region of stabilizing interactions can be observed for Clcomplexes, corresponding to the six interactions of the anion with the central carbons of the sumanene derivative. This attractive region extends over the six Ca atoms (see Figure 9.1), corresponding to six equivalent bond critical points with density of 0.0106 a.u. according to Atoms in Molecules theory.[59] Bromide complex shows a similar plot, though it can be already appreciated in Figure 9.6 that the interaction is slightly weaker (lighter blue), in accordance with the smaller values of the critical points (0.0098 a.u.). The plot is more complex for trigonal anions which, in addition, show two different orientations. NO3- complex shows three attractive regions roughly corresponding to contacts with the three phenyl rings in SumaCN. These three regions are associated to a set of six bond critical points with density of 0.0111 a.u. The different orientation of the most stable complex with CO2H- anion (B) gives rise to a different pattern, with two pairs of bond critical points connecting the anion with the Cb atoms of the pentagonal (0.0124 a.u.) and hexagonal (0.0125 a.u.) rings. An attractive interaction between the anion and the hydrogen atom of the CH2 group can also be observed in the NCI plot. Finally, for BF4- and ClO4- complexes the behavior is similar. In the most stable orientation three oxygen atoms point to three different phenyl rings, so three attractive regions are displayed in the NCI plot. These attractive regions are related to a set of six bond critical points connected to the Cc atoms, with densities of 0.0114 a.u. for BF4- and 0.0113 a.u. for ClO4-. CoraCN and SumaCN2 complexes show similar NCI plots (Appendix D). The largest qualitative difference is observed for SumaCN2 complexes which also display areas of weak interaction which extend towards the CN groups. These regions are already observed in Clcomplex, becoming larger as the size of the anion increases, and they are probably associated to the steric hindrance experienced by the anions due to the disposition of the CN groups, as commented above. Therefore, these plots allow identifying the main attractive interactions established and qualitatively obtaining information about the intensity of the interaction. 9.3.5. SAPT(DFT) Energy Analysis Figure 9.7 summarizes the results of the Symmetry Adapted Perturbation Theory calculations for the most stable complexes formed by each bowl with the anions considered in this study. In agreement with previous work,[23] in the complex formed by Clanion and CoraCN the largest contribution to the interaction energy is the electrostatic one, reaching -41 kcal/mol, with contributions from induction and Alba Campo Cacharrón 236 dispersion reaching -14 and -13 kcal/mol, respectively. Due to the large cancellation of electrostatic and repulsion contribution, the stability of the complex mainly comes from the combined effect of dispersion and induction. When bromide is coordinated to CoraCN, there is an increase in repulsion relative to Clcomplex as a consequence of the larger size of the anion. Dispersion slightly increases due to the same reason, whereas induction decreases by a similar amount (this could be related with the slightly larger equilibrium distance in Brcomplex). On the other hand, the electrostatic contribution increases with respect to that observed in chloride complex. A priori, the opposite effect would be expected, due to the longer equilibrium distance but, as indicated in previous work,[23] the electrostatic contribution comes to a significant extent from penetration effects and this could be the reason of the larger electrostatic contribution in bromide complex. Chloride and bromide complexes with SumaCN and SumaCN2 show the same behavior already described for CoraCN. Therefore, in all cases there are increments in the electrostatic, induction and dispersion contributions, though partially cancelled out by the increment in repulsion. Nitrate complexes with CoraCN and SumaCN show a large electrostatic contribution together with significant repulsion, though it must be taken into account that the balance of these two contributions is less favorable than for both Cl- and Brcomplexes. Induction is of similar magnitude as in complexes with monatomic anions, whereas dispersion exhibits a significant increase as a consequence of the larger size of nitrate anion and the close proximity to the walls of the bowl in orientation A. In the complex with SumaCN2 the repulsion contribution undergoes a large increase together with a small decrease in electrostatic contribution, so the combination of these first-order terms gives an overall balance of +7 kcal/mol (-6 kcal/mol in SumaCN). Dispersion and induction contributions partially cancel this effect showing increments of around 1-2 kcal/mol each with respect to SumaCN complex. The large repulsion contribution has its origins in the unfavorable interaction with the CN groups close to the oxygen atoms of the nitrate anion in orientation A. CO2H- anion forms the most stable complexes among all anions studied as a consequence of the large electrostatic contribution observed, just a bit smaller than in Brcomplexes. This, together with a significant increase in the induction contribution and also large dispersion, overcomes the larger repulsion consequence of the short distances between the oxygen atoms of the anion in orientation B and the carbon atoms of the bowl. Tetrahedral anion complexes exhibit almost the same characteristics for a given buckybowl. Thus, the energy contributions are almost the same for both anions, though dispersion is always larger for ClO4-, being the reason of the larger stability of its complexes as compared with BF4- ones. Electrostatic contributions are the smallest ones observed among the different complexes, as also happens to induction. 9. Anion’s Nature and solvent effects 237 Figure 9.7. SAPT(DFT) decomposition of complexes formed employing the I approaching line and the most stable orientation. Cl Br NO3 CO2H BF4 ClO4 -50 -40 -30 -20 -10 0 10 20 30 40 50 -40.8 -42.7 -34.8 -41.3 -31.1 -31.0 31.9 35.3 31.9 36.4 29.1 30.0 -14.4 -13.0 -13.1 -17.0 -12.0 -11.3 -13.0 -14.1 -18.0 -15.9 -14.7 -17.9 -36.3 -34.5 -34.0 -37.6 -28.7 -30.1 CoraCN Elec Rep Ind Disp Total E (kcal mol -1 ) Cl Br NO3 CO2H BF4 ClO4 -50 -40 -30 -20 -10 0 10 20 30 40 50 -50.0 -51.7 -42.1 -49.4 -36.1 -35.7 37.0 40.6 36.2 40.4 31.5 32.6 -16.1 -14.8 -14.5 -18.7 -12.9 -12.4 -14.8 -16.0 -19.7 -17.6 -15.8 -19.4 -43.9 -42.0 -40.1 -45.4 -33.3 -34.8 SumaCN Elec Rep Ind Disp Total E (kcal mol -1 ) Cl Br NO3 CO2H BF4 ClO4 -50 -40 -30 -20 -10 0 10 20 30 40 50 -43.0 -45.5 -33.5 -43.0 -26.9 -27.2 40.6 45.9 40.7 47.4 32.7 35.1 -16.9 -15.8 -15.6 -19.8 -13.4 -12.8 -16.4 -17.9 -21.2 -19.8 -16.5 -20.5 -35.6 -33.4 -29.6 -35.2 -24.0 -25.4 SumaCN2 Elec Rep Ind Disp Total E (kcal mol -1 ) CoraCN SumaCN SumaCN2 Cl-Br-CO2H-ClO4- BF4- NO3- Cl-Br-CO2H-ClO4- BF4- NO3- Cl-Br-CO2H-ClO4- BF4- NO3- Energy(kcal/mol) Energy(kcal/mol) Energy(kcal/mol) Alba Campo Cacharrón 238 Figure 9.8. SAPT(DFT) decomposition of complexes formed with Br-, NO3-, ClO4- and SumaCN along different approaching lines employing the most stable orientation. I O H R -50 -40 -30 -20 -10 0 10 20 30 40 50 -51.7 -37.0 -28.8 -20.9 40.6 29.4 21.4 27.4 -14.8 -14.4 -10.9 -15.4 -16.0 -9.9 -8.2 -8.5 -42.0 -31.9 -26.5 -17.4 Br Elec Rep Ind Disp Total E (kcal mol -1) I O H R -50 -40 -30 -20 -10 0 10 20 30 40 50 -42.1 -31.7 -26.1 -18.3 36.2 25.6 20.5 19.8 -14.5 -14.4 -11.6 -11.8 -19.7 -10.1 -11.1 -7.8 -40.1 -30.7 -28.3 -18.2 NO3 Elec Rep Ind Disp Total E (kcal mol -1) I O H R -50 -40 -30 -20 -10 0 10 20 30 40 50 -35.7 -26.2 -22.9 -12.3 32.6 19.0 17.5 16.0 -12.4 -10.9 -9.7 -10.2 -19.4 -10.5 -10.1 -8.3 -34.8 -28.6 -25.2 -14.8 ClO4 Elec Rep Ind Disp Total E (kcal mol -1) Br- NO3- ClO4- I O RH I O RH I O RH Energy(kcal/mol) Energy(kcal/mol)Energy(kcal/mol) 9. Anion’s Nature and solvent effects 239 On the other hand, for these bulky anions dispersion is large, only nitrate complexes showing larger dispersion contributions. Even though large repulsion contributions could be expected, these anions are located further away than smaller ones so repulsion is not as large as for other anions. Thus, repulsion contribution is the smallest among complexes studied, but its relative weight is the largest. Figure 9.8 explains the preference for I orientation over other possibilities. When the anion approaches the bowl by the central ring of the convex face of the bowl (O), there is a sharp decrease on electrostatic contribution relative to I, since the MEP by the convex face of the bowls is less positive. Repulsion also decreases since there are few atoms close by, but the combination of both contributions already favors by 3 kcal/mol the formation of the Brcomplex by the concave side. However, as the size of the anion grows, the combination of electrostatic plus repulsion becomes more favorable to the O complex (see results for ClO4-). Induction shows similar contributions by both faces of the bowl (I and O), while dispersion decreases by 6-10 kcal/mol for the O structures. Complexes by the H line are even less stable since all contributions diminish, especially the electrostatic (and induction) one. When the anion interacts directly with a C-H group (R) forming a hydrogen bond to the rim of the bowl, the electrostatic contribution is the smallest one observed for the different attacking lines. Dispersion also decreases, as expected taking into account the larger distances to other atoms. This, combined with a repulsion larger than the electrostatic contribution makes these complexes the least stable ones among those considered in this work. Figure 9.9. SAPT(DFT) decomposition of complexes formed with NO3-, ClO4- and SumaCN employing different anion orientation and the I approaching line. NO3 – A NO3 – B ClO4 – A ClO4 – B ClO4 – C -50 -40 -30 -20 -10 0 10 20 30 40 50 -42.1 -38.4 -28.8 -35.7 -33.5 36.2 30.8 21.7 32.6 28.8 -14.5 -14.2 -10.4 -12.4 -11.7 -19.7 -16.1 -14.6 -19.4 -17.4 -40.1 -37.9 -32.1 -34.8 -33.9 SAPT - SumaCN Elec Rep Ind Disp Total E (kcal mol -1) NO3--A NO3--A ClO4--A ClO4--CClO4--B Energy(kcal/mol) Alba Campo Cacharrón 240 The source of the orientation dependence of the interaction energy can be analyzed in Figure 9.9, which shows the effect of the orientation of the anion on the interaction when approaching SumaCN by the I line. In the case of nitrate complexes, orientation A introduces more repulsion since the three oxygen atoms are closer to the bowl (only two in B). The larger repulsion in A is compensated by larger dispersion and electrostatic contributions, the result being that minimum A is favored. In fact, it is the larger contribution from dispersion which makes the formation of the A complex more favorable (in the case of CO2H- induction is what makes B more stable). In ClO4- complexes orientation B is the most stable one because, despite the larger repulsion contribution, it maximizes the stabilizing contributions, reaching the largest values for electrostatic, dispersion and induction. 9.3.6. Solvent Effect The interaction energies for the complexes studied have been obtained in different solvents by applying the C-PCM method at the M062X/6-31+G* level. The presence of a dielectric medium has a deep impact onto the energetics of the complexes, but in all cases the most stable structures still correspond to complexes formed by the I line and, overall, the stability order obtained for the different attacking lines is the same as observed in the gas phase. Solvent effect is larger in complexes formed with the concave face of the bowl, where the anion has to desolvate to a larger extent since part of its surface is contacting with the bowl. On the other hand, R complexes are the least penalized since almost the whole anion surface is still exposed to the solvent. Despite this, the stability differences observed in the gas phase are so large that no changes in the stability order due to the solvent depending on the attacking line have been observed. Therefore, Table 9.4 only shows I complexes in a series of solvents with varying dielectric constant for the most stable orientation of the anion. Complexes formed with CoraCN in toluene already exhibit significant changes with respect to the values obtained in the gas phase despite the low dielectric constant of the solvent. Thus, complexation energies decrease by around one half in all cases with respect to the gas phase, so the most stable complexes barely reach -19 kcal/mol. Further decreases are observed as the dielectric constant of the solvent increases, and already in chloroform the interaction energies are around one third of the gas phase values. In dimethylchloride and solvents with larger dielectric constants the interaction energies stabilize in values around 20-15 % of the original gas phase values. Thus, in water, complexation energies hardly reach -7 kcal/mol for any of the complexes with CoraCN. SumaCN complexes also show similar changes as the solvent becomes more polar. However, as in the gas phase, complexes formed with SumaCN are more stable than those formed with CoraCN. In chloroform, the most favorable complexes exhibit complexation energies around -15 to -16 kcal/mol, which reduce even more and hardly reach -9 kcal/mol at most in water. Complexes of Br- and Clwith 9. Anion’s Nature and solvent effects 241 SumaCN2 show similar values to those obtained in SumaCN as the dielectric constant of the solvent increases. Polyatomic anions exhibit smaller stabilities as a consequence of the reduced size of the cavity due to the presence of the CN groups in the concave side of the bowl. Table 9.4. M06-2X/6-31+G* interaction energies, in kcal/mol, for complexes formed employing I approaching line, the most stable orientations and different solvents modeled with C-PCM, as obtained at the minima of the M06-2X/6-31+G* gas-phase curves. Gasphase Toluene CHCl3 THF Cl2C2H4 Ethanol DMSO H2O  1.00 2.37 4.71 7.43 10.13 24.85 46.82 78.36 CoraCN Cl- -36.97 -18.16 -11.46 -9.00 -7.87 -6.03 -5.44 -5.17 Br- -35.39 -18.42 -12.37 -10.15 -9.13 -7.47 -6.93 -6.69 NO3- -35.15 -17.90 -11.68 -9.38 -8.32 -6.58 -6.03 -5.77 CO2H- -39.47 -18.62 -10.88 -7.98 -6.62 -4.41 -3.70 -3.37 BF4- -30.77 -15.18 -9.55 -7.45 -6.49 -4.90 -4.39 -4.16 ClO4- -30.86 -16.22 -10.99 -9.06 -8.16 -6.71 -6.24 -6.03 SumaCN Cl- -45.29 -22.96 -14.91 -11.93 -10.55 -8.31 -7.59 -7.26 Br- -43.46 -23.21 -15.92 -13.21 -11.97 -9.94 -9.28 -8.99 NO3- -41.66 -21.62 -14.34 -11.63 -10.38 -8.34 -7.68 -7.38 CO2H- -48.02 -24.16 -15.24 -11.88 -10.31 -7.74 -6.91 -6.53 BF4- -35.59 -17.92 -11.51 -9.12 -8.02 -6.21 -5.63 -5.37 ClO4- -35.57 -18.90 -12.96 -10.75 -9.74 -8.08 -7.55 -7.30 SumaCN2 Cl- -37.33 -18.83 -12.60 -10.36 -9.35 -7.72 -7.20 -6.96 Br- -35.40 -19.06 -13.57 -11.61 -10.72 -9.30 -8.85 -8.64 NO3- -31.38 -14.94 -9.24 -7.17 -6.22 -4.69 -4.20 -3.98 CO2H- -38.40 -17.92 -10.53 -7.78 -6.52 -4.45 -3.79 -3.48 BF4- -25.85 -11.77 -6.96 -5.21 -4.41 -3.12 -2.71 -2.52 ClO4- -26.18 -13.02 -8.59 -7.00 -6.27 -5.11 -4.74 -4.57 Alba Campo Cacharrón 242 Figure 9.10. Variations in interaction energies relative to Clcomplexes as the dielectric constant of the solvent is changed. C-PCM results obtained at the M06-2X/6-31+G* level. dielectric constant 010 20 30 40 50 60 70 80 Eint(kcal mol-1) -2 0 2 4 6 dielectric constant 010 20 30 40 50 60 70 80 Eint(kcal mol-1) -2 0 2 4 6 dielectric constant 010 20 30 40 50 60 70 80 Eint(kcal mol-1) -2 0 2 4 6 Col 1 vs Col 2 Col 1 vs Col 3 Col 1 vs Col 4 Col 1 vs Col 5 Col 1 vs Col 6 Col 7 vs Col 8 CoraCN Col 1 vs Col 2 Col 1 vs Col 3 Col 1 vs Col 4 Col 1 vs Col 5 Col 1 vs Col 6 Col 7 vs Col 8 Col 1 vs Col 2 Col 1 vs Col 3 Col 1 vs Col 4 Col 1 vs Col 5 Col 1 vs Col 6 Col 7 vs Col 8 Cl- Br- NO3- CO2H- BF4- ClO4- Col 1 vs Col 2 Col 1 vs Col 3 Col 1 vs Col 4 Col 1 vs Col 5 Col 1 vs Col 6 Col 7 vs Col 8 Col 1 vs Col 2 Col 1 vs Col 3 Col 1 vs Col 4 Col 1 vs Col 5 Col 1 vs Col 6 Col 7 vs Col 8 Cl- Br- NO3- CO2H- BF4- ClO4- Col 1 vs Col 2 Col 1 vs Col 3 Col 1 vs Col 4 Col 1 vs Col 5 Col 1 vs Col 6 Col 7 vs Col 8 Col 1 vs Col 2 Col 1 vs Col 3 Col 1 vs Col 4 Col 1 vs Col 5 Col 1 vs Col 6 Col 7 vs Col 8 Cl- Br- NO3- CO2H- BF4- ClO4- Col 1 vs Col 2 Col 1 vs Col 3 Col 1 vs Col 4 Col 1 vs Col 5 Col 1 vs Col 6 Col 7 vs Col 8 Col 1 vs Col 2 Col 1 vs Col 3 Col 1 vs Col 4 Col 1 vs Col 5 Col 1 vs Col 6 Col 7 vs Col 8 Cl- Br- NO3- CO2H- BF4- ClO4- SumaCN2 SumaCN Energy(kcal/mol) Energy(kcal/mol)Energy(kcal/mol) 9. Anion’s Nature and solvent effects 243 However, the effect of the solvent upon the stability of a given complex is remarkably different for each anion. While bromide complex becomes the most stable complex with CoraCN as the dielectric constant grows, CO2H- complex is penalized to a larger extent, becoming one of the least stable complexes obtained. Chloride and NO3- complexes exhibit similar interaction energies, whereas ClO4- complexes are relatively favored, reaching interaction energies similar to those of NO3- complexes. These trends are more easily seen in Figure 9.10, which shows the changes in interaction energy of the complexes as the dielectric constant changes relative to chloride complexes. As commented above, stability differences in the gas phase are quite large, spanning an energy range of around 8-12 kcal/mol. However, the incorporation of the solvent quickly decreases this range and already in chloroform it only amounts to around 4 kcal/mol. It can be observed that all anion complexes with CoraCN gain stability relative to chloride ones, except in the case of CO2H-, which is the most destabilized by the solvent as commented above. For dielectric constants around 10 or larger, only BF4- and CO2H- complexes are less stable than Clones; NO3- and ClO4- exhibit similar stabilities, whereas Brcomplexes are clearly the most stable ones. The same behavior is observed for SumaCN complexes, even though in this case Clcomplexes are equally stable as NO3- and ClO4-. SumaCN2 complexes show a different behavior, and in this case only Brcomplexes are more stable than Clones, independently of the dielectric constant of the solvent. Considering the decreases in interaction energies relative to gas phase, it becomes clear that the most destabilized complexes are CO2H- ones, followed by Cl- and BF4- ones. On the other hand, Br- and ClO4- complexes are the least affected by the presence of the solvent, while NO3- occupies an intermediate position. These trends are the consequence of two main factors. First, the larger the intensity of the interaction between bowl and anion in the gas phase, the larger would be the stability in solvent, as already indicated by the fact that the stability order depending on the attacking line is preserved from the gas phase. Second, the desolvation cost of the anion has the larger impact on the final energy differences. In this case, the larger the interaction in the gas phase the larger the interaction with the solvent molecules and the larger the desolvation cost. Thus, polarizing anions such as Cl- and CO2H- are clearly disfavored by desolvation, whereas bulkier anions as Br- and ClO4- are the ones taking the larger benefit. Overall, it can be observed that complexes formed with Branion are the most stable ones, whereas the least stable are the BF4- ones. As regards the rest of the anions, there are changes in the order of stability depending on the solvent and the bowl considered. Overall, the presence of the solvent clearly favors complexes with Br-, while CO2H- undergoes the largest destabilization becoming one of the least stable complexes Alba Campo Cacharrón 251 Complexes involving ion···π interactions have been computationally studied in order to gain insight about the characteristics and factors controlling the interactions. A variety of systems have been considered ranging from cation···π contacts in simple aromatic moieties to ion···π interactions in more complex structures as buckybowls. Below, the main conclusions reached in these studies are listed.  The incorporation of water molecules to phenol···cation complexes leads to a great variety of similarly stable structures differing on the topology of the hydrogen bond network. The behavior is more complex than in benzene analogues because the hydroxyl group of phenol starts participating in hydrogen bonds when three and four water molecules are included.  Li+ and Mg2+ complexes show structures with water and phenol surrounding the ions without forming hydrogen bonds among themselves. Otherwise Na+ and K+ interact more weakly with water and phenol so hydrogen bonds start being competitive, leading to minima with hydrogen bonds between water molecules or between water and phenol.  Vibrational spectra in the region of the O-H stretching mode are quite simple in Li+ and Mg2+ complexes, with shifts associated to hydrogen bonds being only observed when the fourth water molecule is located in the second solvation shell. Na+ and K+ complexes show more complex spectra, with important red shifts associated to O-H···O hydrogen bonds between water molecules but also with participation of phenol in φ-OH···O hydrogen bonds.  Ternary complexes formed between guanidinium cation and aromatic units from amino acids show a variety of minima which can be roughly grouped as parallel stacked (two parallel rings), T-shaped (perpendicular rings) and doubly T-shaped (cation between both rings). Most stable structures are doubly T-shaped ones, followed by T- shaped and parallel stacked minima  The interaction is mostly controlled by the intensity of the cation···aromatic contacts in the cluster. Guanidinium prefers interacting with indole than with phenol and benzene. The formation of hydrogen bonds in complexes containing phenol or indole introduces extra stability in T-shaped structures.  Three-body effects are only cooperative in T-shaped minima containing indole or phenol when hydrogen bonds are formed. Though the interaction is dominated by electrostatics, contributions from dispersion and induction are also relevant. 10. Conclusions 252  The study of the interaction of anions in simplified models for a synthetic anion channel based in naphthalendiimide (NDI) units shows that anion···π interactions are highly favorable in the gas phase, with minima connecting a path that allows the anions to move along the aromatic surface.  The presence of water molecules strongly affects the interaction of fluoride and hydroxide anions with NDI, whereas the effect is less remarkable in the case of bromide and chloride complexes. As a consequence, the intensity of the interaction becomes similar among the different anions, with differences below 4 kcal/mol (in gas phase these differences reach more than 20 kcal/mol).  The results suggest that a limited number of water molecules attached to the anion could be crucial to overcome dehydration costs, also contributing to the stabilization of the complex by means of favorable water···NDI contacts.  Substitution of corannulene and sumanene is an effective method to promote great changes in the electric properties of the bowls, hardly changing their geometry and keeping their characteristic bowl shape.  Substitution allows modulating the molecular electrostatic potential (MEP) of the bowls so they can be tuned for favorable interaction with anions or cations. Electronwithdrawing groups as nitrile create positive MEP regions on both faces of the bowl; electron-donating groups as methyl and unsubstituted bowls exhibit negative MEP regions on both faces; other groups as fluoride produce almost zero MEP allowing favorable interaction with both anions and cations.  The interaction in these systems is hard to describe with accuracy. Testing a variety of methods it has been found that the best performance is given by the SCS-MP2 method extrapolated to basis limit, which matches almost perfectly the reference values obtained at the MP2.X/6-31G(0.25) level of calculation. M06-2X/6-31+G* also seems to be able to provide a quite balanced description of the concave/convex energy differences despite the errors introduced in the interaction energies.  Sumanene and its derivatives interact more favorably than corannulene with anions and cations. Substitution effects in sumanene are stronger if applied in the CH aromatic groups than in the CH2 groups of pentagonal rings. Cation complexes are dominated by a large induction component that can be even larger than the electrostatic contribution, whereas in anion complexes electrostatics and dispersion are the main stabilizing components. Alba Campo Cacharrón 253  Complexation energies roughly follow the series of MEPs, being the most stable complexes those formed with corannulene substituted with CN groups. However, there are also changes affecting to the induction and dispersion contributions which may introduce deviations over the behavior expected from a purely electrostatics point of view.  Anions are oriented so the largest negative regions point towards the positively charged surfaces of the bowl. Tetrahedral anions orient themselves with three bonds towards the bottom of the bowl. Nitrate complexes are more stable with the anion parallel to the bottom of the bowl, the oxygen atoms oriented towards the positively charged six-carbon rings. Formate complexes are more stable with the anion perpendicular to the bottom of the bowl, avoiding the interaction of the hydrogen atom with the walls of the bowl.  The most stable complexes correspond to the anion located on the symmetry axis of the bowl by the concave face. The interaction of the anions with the convex face of the bowl or with the hydrogen atoms on the rim of the bowl is significantly less favorable. The order of stability of the complexes formed with CN-substituted bowls in the gas phase as obtained for the different anions is the following: CO2H- > Cl- > Br- > NO3- >> ClO4- > BF4-.  The presence of solvent has a deep impact on the properties of the complexes, and already in solvents with low dielectric constant the interaction strength decreases dramatically. Solvent affects the different anions to distinct extents, so the more polarizing ones are the most affected whereas the bulkier ones register smaller decreases. Therefore, in solvents of moderate to large dielectric constant, the most stable complexes are those formed with Brwhereas CO2H- is among the least favorable complexes formed.  Overall, the results obtained for ion complexes formed with properly substituted buckybowls encourage to follow this topic, since this strategy promises to be suitable for designing selective anion receptors. Appendices Appendix A 257 Appendix A Table A.1. Complexation energies (kcal/mol) for monohydrated cation···phenol clusters as obtained with different levels of calculation. Phe-X1- 1H2O Phe-X2- 1H2O Phe-X3- 1H2O Phe-X4- 1H2O Phe-X5- 1H2O Phe-X6- 1H2O B3LYP(a) K+ -31.03 -28.16 -30.90 -30.77 -27.61 -32.35 MP2-A(b) -31.12 -29.39 -33.18 -32.38 -29.14 -34.98 MP2-B(c) -31.87 -30.24 -33.30 -33.07 -29.94 -35.72 MP2-C(d) -31.91 -29.29 -33.19 -32.18 -29.10 -34.87 MP2-D(e) -32.33 -30.01 -33.94 -32.99 -30.20 -36.13 MP2-E(f) -32.28 -30.26 -34.08 -33.30 -30.78 -36.62 B3LYP(a) Na+ -44.89 -37.57 -44.89 -45.01 -37.40 -44.72 MP2-A(b) -41.45 -35.23 -41.46 -41.52 -34.80 -42.94 MP2-B(c) -41.80 -36.13 -42.41 -42.35 -35.56 -43.70 MP2-C(d) -41.94 -35.65 -42.55 -41.88 -35.09 -43.23 MP2-D(e) -41.25 -35.32 -42.10 -40.95 -34.48 -42.45 MP2-E(f) -41.05 -35.35 -41.99 -40.95 -34.73 -42.63 B3LYP(a) Li+ -64.60 -51.59 -65.71 -53.44 MP2-A(b) -61.00 -48.67 -62.14 -49.74 MP2-B(c) -61.01 -48.99 -62.23 -50.02 MP2-C(d) -61.17 -48.98 -62.27 -50.02 MP2-D(e) -60.73 -49.10 -61.82 -49.93 MP2-E(f) -60.63 -49.04 -61.77 -49.99 B3LYP(a) Mg+2 -183.04 -153.00 -185.71 -151.48 MP2-A(b) -171.18 -139.90 -174.86 -140.31 MP2-B(c) -171.40 -140.45 -174.91 -140.61 MP2-C(d) -171.55 -140.50 -175.00 -140.62 MP2-D(e) -170.32 -139.71 -173.15 -139.05 MP2-E(f) -170.98 -140.37 -174.03 -139.85 a) B3LYP/6-31+G(2d,p); b) MP2/6-31+G(2d,p)//B3LYP/6-31+G(2d,p); c) MP2/6-31+G(2d,p)//MP2/6-31+G(d); d) MP2/631+G(2d,p); e) MP2/6-311++G(2d,2p)//MP2/6-31+G(d); f) MP2/6-311++G(3d,2p) //MP2/6-31+G(d). Alba Campo Cacharrón 258 Table A.2. OH stretching frequencies for selected complexes with sodium as obtained at the MP2/6-31+G(2d,p) level. Values obtained with the harmonic approximation and corrected from anharmonicity are presented. Harmonic Anharmonic Factor(*) 0.9797 1.0321 Phenol 3657.0 3657.0 Water 3672.9 3688.3 Phe-Na-O 3646.9 3649.2 Phe-Na- 3631.2 3631.8 Phe-Na1-1H2O 3758.3 3816.0 Phe-Na1-1H2O 3653.9 3685.2 Phe-Na1-1H2O 3648.9 3650.5 Phe-Na2-1H2O 3772.6 3778.8 Phe-Na2-1H2O 3650.7 3661.7 Phe-Na2-1H2O 3326.5 3251.6 Phe-Na3-1H2O 3750.0 3744.0 Phe-Na3-1H2O 3646.7 3647.0 Phe-Na3-1H2O 3629.3 3635.7 Phe-Na4-1H2O 3758.3 3829.1 Phe-Na4-1H2O 3654.0 3688.4 Phe-Na4-1H2O 3634.5 3635.7 Phe-Na5-1H2O 3766.0 3763.4 Phe-Na5-1H2O 3644.8 3654.3 Phe-Na5-1H2O 3372.9 3327.6 Phe-Na6-1H2O 3739.4 3736.8 Phe-Na6-1H2O 3626.3 3627.1 Phe-Na6-1H2O 3571.4 3570.7 (*) As commented in the text, a correction factor is applied in order to reproduce the experimental OH stretching frequency of phenol. Appendix B 259 Appendix B Table B.1. LMO-EDA partition (kcal/mol) of the interaction energy for complexes in the manuscript with equal aromatic units, as obtained at the M06-2X/aug-cc-pVDZ level. Electrostatic Exchange Repulsion Polarization Dispersion Bz-Bz-1 -15.05 -12.19 41.04 -9.69 -23.66 Bz-Bz-2 -15.89 -10.85 39.41 -10.19 -24.52 Bz-Bz-3 -16.09 -11.76 42.76 -9.13 -27.31 Bz-Bz-4 -19.61 -10.63 38.75 -12.78 -22.71 Ph-Ph-1 -35.71 -20.37 68.35 -16.88 -28.96 Ph-Ph-2 -33.56 -19.94 69.29 -16.11 -32.16 Ph-Ph-3 -32.52 -20.05 64.61 -19.78 -23.11 Ph-Ph-4 -29.29 -17.47 62.49 -12.44 -34.64 Ph-Ph-5 -34.16 -19.41 65.84 -15.35 -29.98 Ph-Ph-6 -30.16 -14.78 51.85 -15.85 -24.49 Ph-Ph-7 -34.13 -16.37 55.26 -16.20 -21.12 In-In-1 -24.12 -20.04 66.30 -14.00 -38.78 In-In-2 -26.13 -17.82 59.87 -14.72 -33.97 In-In-3 -26.70 -17.60 58.99 -15.93 -32.34 In-In-4 -26.89 -17.48 58.98 -15.72 -32.99 In-In-5 -29.66 -15.16 54.88 -16.34 -32.72 Alba Campo Cacharrón 266 . Figure C.2. Interaction energies obtained for complexes formed by unsubstituted, CH3- substituted and CN-substituted corannulene and sumanene with sodium cation by the concave (left) and convex (right) faces of the bowls. R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -34 -32 -30 -28 -26 -24 -22 -20 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -34 -32 -30 -28 -26 -24 -22 -20 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -34 -32 -30 -28 -26 -24 -22 -20 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -34 -32 -30 -28 -26 -24 -22 -20 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS MP2.X/6-31G* MP2.X/6-31+G* BLYP M062X SAPT/scaledMP2.X/6-31G(0.25) Cora-in Cora-out Suma-in Suma-out Energy (kcal/mol) Energy (kcal/mol)Energy (kcal/mol) Energy (kcal/mol) R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -18 -16 -14 -12 -10 -8 -6 -4 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -20 -18 -16 -14 -12 -10 -8 -6 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -18 -16 -14 -12 -10 -8 -6 -4 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -18 -16 -14 -12 -10 -8 -6 -4 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -18 -16 -14 -12 -10 -8 -6 -4 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -18 -16 -14 -12 -10 -8 -6 -4 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS MP2.X/6-31G* MP2.X/6-31+G* BLYP M062X SAPT/scaledMP2.X/6-31G(0.25) CoraMet-in CoraMet-out SumaMet-in SumaMet-out SumaMet2-in SumaMet2-out Energy (kcal/mol) Energy (kcal/mol)Energy (kcal/mol) Energy (kcal/mol)Energy (kcal/mol) Energy (kcal/mol) R (Å) 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 E (kcal mol-1) -4 -2 0 2 4 6 8 10 R (Å) 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 E (kcal mol-1) -4 -2 0 2 4 6 8 10 R (Å) 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 E (kcal mol-1) 2 4 6 8 10 12 14 16 R (Å) 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 E (kcal mol-1) 2 4 6 8 10 12 14 16 R (Å) 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 E (kcal mol-1) -6 -4 -2 0 2 4 6 8 R (Å) 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 E (kcal mol-1) -6 -4 -2 0 2 4 6 8 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS MP2.X/6-31G* MP2.X/6-31+G* BLYP M062X SAPT/scaledMP2.X/6-31G(0.25) CoraCN-in CoraCN-out SumaCN-in SumaCN-out SumaCN2-in SumaCN2-out Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol) Appendix C 267 Figure C.3. Relative energy (in-out) of complexes formed by chloride anion with the different bowls. Figure C.4. Relative energy (in-out) of complexes formed by sodium cation with the different bowls. R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -12 -10 -8 -6 -4 -2 0 2 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS MP2.X/6-31G* MP2.X/6-31+G* BLYP-D3 M062X SAPT/scaledMP2.X/6-31G(0.25) CoraF SumaF SumaF2 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS BLYP-D3 M062X SAPT/scaledMP2.X/6-31G(0.25) Cora Suma R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS BLYP-D3 M062X SAPT/scaledMP2.X/6-31G(0.25) Cora Suma R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS BLYP-D3 M062X SAPT/scaledMP2.X/6-31G(0.25) Cora Suma R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS BLYP-D3 M062X SAPT/scaledMP2.X/6-31G(0.25) Cora Suma Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol)Energy (kcal/mol) R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -12 -10 -8 -6 -4 -2 0 2 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -12 -10 -8 -6 -4 -2 0 2 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS MP2.X/6-31G* MP2.X/6-31+G* BLYP-D3 M062X SAPT/scaledMP2.X/6-31G(0.25) CoraMet SumaMet SumaMet2 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -12 -10 -8 -6 -4 -2 0 2 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -12 -10 -8 -6 -4 -2 0 2 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -14 -12 -10 -8 -6 -4 -2 0 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS MP2.X/6-31G* MP2.X/6-31+G* BLYP-D3 M062X SAPT/scaledMP2.X/6-31G(0.25) CoraCN SumaCN SumaCN2 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS BLYP-D3 M062X SAPT/scaledMP2.X/6-31G(0.25) Cora Suma R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS BLYP-D3 M062X SAPT/scaledMP2.X/6-31G(0.25) Cora Suma R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS BLYP-D3 M062X SAPT/scaledMP2.X/6-31G(0.25) Cora Suma Energy (kcal/mol) Energy (kcal/mol)Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol)Energy (kcal/mol) R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -4 -2 0 2 4 6 8 10 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -4 -2 0 2 4 6 8 10 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -4 -2 0 2 4 6 8 10 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS MP2.X/6-31G* MP2.X/6-31+G* BLYP M062X SAPT/scaledMP2.X/6-31G(0.25) R (Å) 1,8 2,0 2,2 2,4 2,6 2,8 3,0 3,2 E (kcal mol-1) -4 -2 0 2 4 6 8 10 R (Å) 1,8 2,0 2,2 2,4 2,6 2,8 3,0 3,2 E (kcal mol-1) -4 -2 0 2 4 6 8 10 R (Å) 1,8 2,0 2,2 2,4 2,6 2,8 3,0 3,2 E (kcal mol-1) -4 -2 0 2 4 6 8 10 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS MP2.X/6-31G* MP2.X/6-31+G* BLYP M062X SAPT/scaledMP2.X/6-31G(0.25) R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS BLYP-D3 M062X SAPT/scaledMP2.X/6-31G(0.25) Cora Suma R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS BLYP-D3 M062X SAPT/scaledMP2.X/6-31G(0.25) Cora Suma R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS BLYP-D3 M062X SAPT/scaledMP2.X/6-31G(0.25) Cora Suma CoraF SumaF Cora Suma SumaF2 Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol)Energy (kcal/mol) R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -4 -2 0 2 4 6 8 10 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -4 -2 0 2 4 6 8 10 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -6 -4 -2 0 2 4 6 8 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS MP2.X/6-31G* MP2.X/6-31+G* BLYP-D3BJ M062X SAPT/scaledMP2.X/6-31G(0.25) R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -4 -2 0 2 4 6 8 10 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -4 -2 0 2 4 6 8 10 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -6 -4 -2 0 2 4 6 8 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS MP2.X/6-31G* MP2.X/6-31+G* BLYP-D3BJ M062X SAPT/scaledMP2.X/6-31G(0.25) R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS BLYP-D3 M062X SAPT/scaledMP2.X/6-31G(0.25) Cora Suma R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS BLYP-D3 M062X SAPT/scaledMP2.X/6-31G(0.25) Cora Suma R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 R (Å) 2,4 2,6 2,8 3,0 3,2 3,4 3,6 3,8 E (kcal mol-1) -10 -8 -6 -4 -2 0 2 4 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS BLYP-D3 M062X SAPT/scaledMP2.X/6-31G(0.25) Cora Suma CoraCN SumaCN SumaCN2 CoraMet SumaMet SumaMet2 Energy (kcal/mol) Energy (kcal/mol)Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol)Energy (kcal/mol) Alba Campo Cacharrón 268 Figure C.5. SAPT(DFT) energy decomposition for complexes formed by chloride anion and unsubstituted, CH3-substituted and CN-substituted bowls. The vertical line indicates the position of the minimum obtained at the MP2.X level of calculation. R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 Eele -Erep Eind Edis Etot Eele -Erep Eind Edis Etot Cora-in Cora-out Suma-in Suma-out Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol) R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 Eele -Erep Eind Edis Etot CoraMet-in CoraMet-out SumaMet-in SumaMet-out SumaMet2-in SumaMet2-out Energy (kcal/mol) Energy (kcal/mol)Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol)Energy (kcal/mol) R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 Eele -Erep Eind Edis Etot CoraCN-in CoraCN-out SumaCN-in SumaCN-out SumaCN2-in SumaCN2-out Energy (kcal/mol) Energy (kcal/mol)Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol)Energy (kcal/mol) Appendix C 269 Figure C.6. SAPT(DFT) energy decomposition for complexes formed by sodium cation and unsubstituted, CH3-substituted and CN-substituted bowls. The vertical line indicates the position of the minimum obtained at the MP2.X level of calculation. R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS MP2.X/6-31G* MP2.X/6-31+G* BLYP-D3BJ M062X SAPT/scaledMP2.X/6-31G(0.25) Eele -Erep Eind Edis Etot Cora-in Cora-out Suma-in Suma-out Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol) R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS MP2.X/6-31G* MP2.X/6-31+G* BLYP M062X SAPT/scaledMP2.X/6-31G(0.25) Eele -Erep Eind Edis Etot CoraMet-in CoraMet-out SumaMet-in SumaMet-out SumaMet2-in SumaMet2-out Energy (kcal/mol) Energy (kcal/mol)Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol) R (Å) 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 R (Å) 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 E (kcal mol-1) -50 -40 -30 -20 -10 0 10 MP2/CBS SCS-MP2/CBS SCSN-MP2/CBS MP2.X/6-31G* MP2.X/6-31+G* BLYP-D3BJ M062X SAPT/scaledMP2.X/6-31G(0.25) Eele -Erep Eind Edis Etot CoraCN-in CoraCN-out SumaCN-in SumaCN-out SumaCN2-in SumaCN2-out Energy (kcal/mol) Energy (kcal/mol)Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol) Energy (kcal/mol) Alba Campo Cacharrón 270 Table C.1. SAPT(DFT) values (kcal/mol) for complexes formed by the bowls and chloride anion. The values are obtained by interpolation of the SAPT(DFT) curves at the equilibrium geometry obtained at the MP2.X/6-31G(0.25) level (Table 8.1). Dispersion and dispersion-exchange contributions are scaled by 1.193 as indicated in the text. Eele Erep Eind Edis Etot CONCAVE/IN cora -6.01 21.32 -10.24 -12.76 -7.68 coraMet -2.87 19.49 -11.17 -12.75 -7.31 CoraF -20.62 27.80 -11.40 -14.37 -18.60 coraCN -40.95 33.26 -14.66 -15.81 -38.16 suma -8.33 22.41 -10.84 -13.71 -10.47 sumaMet -5.46 20.21 -11.98 -13.82 -11.05 sumaF -25.63 30.12 -12.17 -15.68 -23.38 sumaCN -49.61 37.59 -16.23 -17.66 -45.91 suma -8.33 22.41 -10.84 -13.71 -10.47 sumaMet2 -10.55 24.96 -14.25 -16.29 -16.14 sumaF2 -25.46 34.77 -13.04 -16.94 -20.68 sumaCN2 -42.29 40.98 -16.86 -19.44 -37.61 CONVEX/OUT cora 0.38 16.16 -10.48 -8.26 -2.20 coraMet 3.50 15.06 -11.00 -8.14 -0.57 CoraF -11.34 21.72 -12.11 -9.54 -11.27 coraCN -34.04 33.83 -17.61 -12.08 -29.90 suma -0.64 16.52 -10.77 -8.59 -3.47 sumaMet 3.86 14.86 -11.28 -8.29 -0.85 sumaF -12.76 21.13 -11.99 -9.66 -13.28 sumaCN -39.72 33.96 -17.79 -12.43 -35.98 suma -0.64 16.52 -10.77 -8.59 -3.47 sumaMet2 -0.92 16.64 -11.72 -9.08 -5.07 sumaF2 -16.73 26.37 -13.73 -10.83 -14.92 sumaCN2 -30.67 33.23 -16.77 -12.46 -26.67 Appendix C 271 Table C.2. SAPT(DFT) values (kcal/mol) for complexes formed by the bowls and sodium cation. The values are obtained by interpolation of the SAPT(DFT) curves at the equilibrium geometry obtained at the MP2.X/6-31G(0.25) level (Table 8.2). Dispersion and dispersion-exchange contributions are scaled by 1.193 as indicated in the text. Eele Erep Eind Edis Etot CONCAVE/IN cora -8.90 8.13 -24.56 -2.66 -27.99 coraMet -12.89 9.02 -27.03 -2.83 -33.72 CoraF 4.62 6.59 -23.36 -2.41 -14.56 coraCN 24.26 4.40 -23.53 -1.99 3.15 suma -8.06 8.24 -25.62 -2.80 -28.24 sumaMet -12.36 9.33 -28.67 -3.01 -34.71 sumaF 7.53 7.04 -24.72 -2.62 -12.77 sumaCN 30.65 4.65 -25.08 -2.15 8.05 suma -8.06 8.24 -25.62 -2.80 -28.24 sumaMet2 -7.74 8.71 -29.05 -3.02 -31.10 sumaF2 7.32 5.08 -23.20 -2.31 -13.12 sumaCN2 20.53 2.69 -22.36 -1.78 -0.92 CONVEX/OUT cora -12.55 8.81 -23.66 -2.05 -29.46 coraMet -16.52 9.51 -25.52 -2.15 -34.67 CoraF -1.44 7.37 -22.65 -1.86 -18.58 coraCN 17.29 5.13 -22.61 -1.54 -1.73 suma -10.54 7.73 -23.92 -2.06 -28.78 sumaMet -15.88 8.56 -26.16 -2.18 -35.67 sumaF 1.31 6.39 -22.79 -1.85 -16.93 sumaCN 23.89 4.00 -22.70 -1.45 3.74 suma -10.54 7.73 -23.92 -2.06 -28.78 sumaMet2 -10.81 8.01 -25.42 -2.13 -30.35 sumaF2 4.16 5.21 -21.60 -1.67 -13.89 sumaCN2 15.16 3.37 -20.56 -1.34 -3.36 Appendix D 273 Appendix D Figure D.1. NCI plots for the complexes formed by CoraCN (left) and SumaCN2 (right) and the anions studied. Only complexes by the I line and the most stable orientation are shown. The product of the density times the sign of the second eigenvalue of its hessian is mapped onto an isosurface of reduced density gradient with value 0.5 a.u. The color scale goes from -0.015 a.u. (blue) to 0.015 a.u. (red). Cl- BF4- NO3- Br- ClO4- CO2H- Cl- BF4- NO3- Br- ClO4- CO2H-