Computational study of the capability of the side chains of amino acids for setting up cation···π interactions relevant in protein stability and structure

Computational Study of the Capability of the Side Chains of Amino Acids for Setting up Cation···πInteractions Relevant in Protein Stability and Structure Ana Angustias Rodríguez Sanz Departamento de Química Física Facultade de Ciencias Lugo, xuño 2015 TESIS DE DOCTORADO UNIVERSIDADE DE SANTIAGO DE COMPOSTELA FACULTADE DE CIENCIAS DEPARTAMENTO DE QUÍMICA-FÍSICA Ph. D. Thesis Computational study of the capability of the side chains of amino acids for setting up cation···π interactions relevant in protein stability and structure Ana Angustias Rodríguez Sanz Lugo, xuño 2015 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: “Computational study of the capability of the side chains of amino acids for setting up cation···π interactions relevant in protein stability and structure”. Presentada por Dna. Ana Angustias Rodríguez Sanz 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 22 xuño de 2015 Asdo: Enrique M. Cabaleiro Lago Asdo: Jesús Rodríguez Otero Asdo: Ana Angustias Rodríguez Sanz Para mi abuela, para mis padres, para mis hermanos Agradecimientos En primer lugar me gustaría mostrar mi agradecimiento a mis directores de tesis, Jesús y Quique. Todavía recuerdo aquella conversación vía Skype que me dio tan buena impresión del proyecto y de vosotros como personas (aunque os empeñarais en quitarle importancia…). Gracias por hacer esto posible y sobre todo, gracias Quique por estar ahí siempre, con tu paciencia infinita a todos mis despistes y por todo lo que me enseñaste, por tu ayuda incondicional e inestimable. Por supuesto doy las gracias también al Ministerio de Economía y Competitividad por la beca FPI, así como la beca de Estancias Breves, en Dublín, donde viví unos meses intensos y conocí gente a la que admiro. En primer lugar quiero dar las gracias a mi jefa Isabel, que sin conocerme me dio todo su apoyo y confianza en mi trabajo desde que entré en su grupo hasta que salí el último día. Mis compañeros de grupo también fueron un gran ejemplo a seguir. Además de ser unas personas maravillosas, su pasión por el trabajo y la ciencia me animaron a continuar. Gracias Elena, Viola, P atrick, Julian, Michela…y todas aquellas personas con las que me crucé en aquella estancia y que me apoyaron y ayudaron tanto, dentro y fuera del laboratorio. Gracias también a mi compañera de despacho, Alba, a la que tanto le debo, y que sin ella todo habría sido mucho más complicado, ¡y por esos encargos de baberos! Este año, apareció una nueva compañera andaluza, Belén, que también hace mis horas frente al ordenador más amenas. Ortigueira, cómo triunfaste con tu mujer, ¡vaya pareja más chula que hacéis! También quiero darle las gracias a Stu y Ronan, que aunque nos separen miles de kilómetros, os voy a tener ahí, no sólo para acribillaros a preguntas de inglés, si no que sé que tengo dos buenos amigos. Por supuesto tengo que agradecer el apoyo constante de mi familia, mis padres, mis hermanos y mi abuela, que siempre estarán ahí para lo que necesite, igual que yo para ellos aunque estemos en la distancia. ¡Qué sería yo sin vosotros! También tengo que agradecer el apoyo de mis amigos de Málaga, Bea, Paco, Rocío, Inma, Samantha… que aunque no estemos continuamente en contacto, cuando nos vemos es como si no hubieran pasado los días. ¡Por esas rutas en bici por la playa y esos viajes por el norte! El día que llegué a Lugo sabía que empezaba una nueva etapa lejos de mi familia y amigos de siempre. Sin embargo, tuve la suerte de encontrar gente increíble desde el primer momento. Empezando por todas mis compañeras de piso que siempre han estado ahí para escuchar mis aventuras y desventuras con la tesis, aguantar mis cambios de ánimo y aún así darme todo su apoyo. Este último año he estado de maravilla, gracias Cris, Esther, Yoli e Iria. Y de entre todas mis compañeras, sé que me llevo unas amigas de verdad para toda la vida, así que gracias Zoraida, Marina, Iria… me he sentido como en casa con vosotras; luego está Alfonso, que empezó dándome consejos sobre plantas carnívoras y prestándome su linterna cuando salía tarde de trabajar, y ahora vivo enganchada a sus historias con las que no puedo parar de reír; ¡y también están cómo no las Marquesitas!: Carmen, mi compi de defensa y rutas en bici por el río; María mi compi de aventuras por tierras hostiles; Patri, mi compi de viaje y mercadillos; y Elba, la consejera amorosa por excelencia. Es una suerte haberos conocido y que ahora forméis parte de mi vida, sólo espero que nunca salgáis de ella. Ya sabéis que allá donde vaya seréis bienvenidos. También quería agradecer a las renombradas Morconas Morunas esas clases de baile y alguna cena que otra tan divertida que te hacen olvidar por unas horas todos los problemas. Y a Mar, por esos ratos de risas en la conserjería junto con el de las historias de antes. Y bueno, por qué no, también le doy las gracias a esos bichitos peludos de Nala, Bimba y Trufa, que siempre estaban ahí esperando detrás de la puerta cuando llegaba a casa para hacerme compañía. Incluso algunos aportando sus ideas en modelización… En definitiva, agradecer a todas las personas que aún sin saberlo, han aportado su granito de arena con su apoyo para ayudarme a no abandonar y seguir adelante de una forma u otra. ¡Así! Resumen iii características de las interacciones no covalentes entre las especies que se encuentren próximas, y en particular las de la interacción catión···π. Además, también sería posible que la estructura de las especies interaccionantes se modificara dependiendo del grado de exposición al disolvente. A medida que se van introduciendo moléculas de disolvente, las distintas moléculas de disolvente compiten con la nube aromática por el catión, de modo que la interacción es cada vez más débil. Estas evidencias constituyen el punto de partida de este proyecto de tesis, en el que se pretende, entre otras cosas, estudiar el efecto que el entorno podía provocar sobre las interacciones catión···π. Es por esto que uno de los aspectos que más ampliamente se ha estudiado en esta tesis ha sido el efecto de la introducción progresiva de moléculas de agua al sistema en el que aparece la interacción catión···π. El estudio de estos procesos se puede plantear de diversas maneras: una de ellas consiste en generar progresivamente nuevos sistemas con una molécula de disolvente más, a partir de las estructuras más estables obtenidas para un determinado nivel de microhidratación. La nueva molécula de agua se añade en diferentes posiciones, siguiendo la intuición química; otra alternativa consiste en hacer una exploración del espacio conformacional de cada uno de los sistemas hidratados para así obtener una muestra representativa de las estructuras más favorables en cada caso, que se emplearán como punto de partida de cálculos posteriores. El conjunto de estructuras obtenido se criba por medio de la energía para reducir la enorme cantidad de estructuras posibles y quedarnos con los sistemas más estables y por lo tanto los más probables. Todas las estructuras seleccionadas de esta forma son sometidas a cálculos mecano-cuánticos de optimización por los que se obtienen estructuras y energías con una mayor calidad. En los siguientes apartados se hará un breve resumen de cada uno de los trabajos realizados en esta tesis, que incluirá también las principales conclusiones a las que se llegó en cada uno de los casos. Efecto de la microhidratación sobre las interacciones metilamonio···fenol y amonio···fenol Este estudio, correspondiente al capítulo 4, se basa en la microhidratación progresiva de complejos formados por un catión y una molécula de fenol, de modo que se van introduciendo hasta tres moléculas de agua para así poder valorar el efecto de la presencia de este número reducido de moléculas de disolvente sobre las características de la interacción catión···π. Existen algunos trabajos previos que tratan a cerca de de la interacción de cationes simples, principalmente alcalinos, con benceno. La interacción con benceno es relativamente simple ya que sólo ofrece una localización para la interacción favorable con cationes. El fenol, perteneciente a la cadena lateral de la tirosina, tiene por el contrario la peculiaridad de que además del sistema aromático presenta un grupo hidroxilo que puede actuar como segundo punto de anclaje para la interacción con el catión, aumentando así el número de posibles complejos iv Resumen catión···π formados. Los cationes elegidos para este trabajo han sido los cationes amonio y metilamonio, en un intento de modelizar las características de la interacción con los residuos catiónicos presentes en la cadena lateral de la lisina. La utilización del catión metilamonio, con posibilidad de interaccionar tanto a través del grupo amonio como a través del grupo metilo, incrementa la complejidad del sistema estudiado. Por lo tanto, en este primer trabajo se trata de determinar las características de la interacción amonio···fenol como una primera aproximación al comportamiento de los contactos entre cadenas laterales de lisina y tirosina. El procedimiento seguido en el estudio ha sido la microhidratación progresiva de los complejos. Como se ha indicado anteriormente, se obtienen las estructuras de los complejos más estables entre cationes amonio y fenol en ausencia de moléculas de agua; empleando estas estructuras, se introduce una primera molécula de agua en las posiciones más favorables y se obtienen las estructuras optimizadas más favorables para los complejos monohidratados. El proceso se continúa introduciendo de modo análogo las siguientes moléculas de agua. La introducción de sucesivas moléculas de agua al sistema aumenta significativamente el número de posibles complejos, ya que se incrementa en gran medida la complejidad de la superficie de energía potencial. Sin embargo, los mínimos encontrados para ambos complejos, tanto con amonio como metilamonio, presentan características parecidas. Los complejos más estables localizados con ambas especies muestran un patrón cíclico de interacción muy similar, en el que suelen estar involucrados el sistema aromático, el grupo hidroxilo del fenol, el catión y al menos una molécula de agua. Este patrón cíclico se puede observar desde la introducción de la primera molécula de agua y permanece con la inclusión del resto de moléculas de agua. Por tanto, los resultados parecen indicar que al menos la primera molécula de agua interacciona de forma específica, formando un ciclo de enlaces de hidrógeno especialmente estable y que sobrevive tras la adición de más moléculas de agua. A medida que se van introduciendo moléculas de agua, las diferencias de estabilidad entre complejos con diferentes estructuras van decreciendo. Así, la estabilización es mayor para la primera molécula de agua, mientras que la inclusión del resto de moléculas de agua provoca una estabilización cada vez menos intensa. Esto se debe a que la primera molécula de agua es capaz de interaccionar simultáneamente tanto con el catión como con la molécula de fenol. Después de este paso, la carga del catión se ha reducido en parte por lo que una segunda molécula de agua interaccionará con éste con menor intensidad. Sin embargo, el grupo hidroxilo también es capaz de interaccionar con nuevas moléculas de agua, involucrando al grupo OH en un nuevo enlace de hidrógeno con el oxígeno del agua de tipo φ-OH···O. Por lo tanto, ya a partir de la segunda molécula de agua, cada vez va siendo menos relevante el lugar que ocupa una nueva molécula de agua, igualándose cada vez más las estabilidades obtenidas por interacción directa con el catión o Resumen v mediante enlaces de hidrógeno con otras moléculas de agua o el fenol. La participación del grupo hidroxilo del fenol en la red de enlaces de hidrógeno formada es significativa. Efecto de la microhidratación sobre la interacción guanidinio···π Al igual que en el trabajo anterior, en este estudio, correspondiente al capítulo 5, se ha llevado a cabo un proceso de microhidratación gradual, pero en este caso en complejos formados por el catión guanidinio con diferentes sistemas aromáticos. Como se mencionó anteriormente, los aminoácidos catiónicos más básicos son la lisina, la prolina y la arginina. Una vez realizado el estudio anterior con los cationes amonio y metilamonio pertenecientes a la lisina, se ha considerado hacer un trabajo similar en el que el catión interaccionante fuera el guanidinio, localizado en la cadena lateral del aminoácido arginina. Sin embargo, a diferencia del trabajo anterior, se incluyeron en el estudio las moléculas de benceno, fenol e indol (pertenecientes a las cadenas laterales de fenilalanina, tirosina y triptófano respectivamente). En un estudio previo de nuestro grupo de investigación, se observa que la incorporación de moléculas de agua a un complejo formado por catión guanidinio y benceno provoca un descenso apreciable en la intensidad de la interacción catión···π, pero promueve además un cambio en la estructura de los complejos, de modo que se vuelven más favorables estructuras en las que el guanidinio se sitúa en paralelo al anillo aromático. Por otra parte, tanto fenol como indol ofrecen una segunda localización favorable para la interacción con el catión, que en el caso del indol se trata de un segundo anillo aromático de cinco miembros en el que uno de los átomos es nitrógeno. Las estructuras más estables en ausencia de moléculas de agua presentan en todos los casos al catión guanidinio situado perpendicularmente a la especie aromática e interaccionando con la misma por medio de sus grupos NH2, estableciendo diferentes enlaces de hidrógeno. Los complejos con benceno sólo presentan un mínimo, mientras que con fenol e indol surgen más posibilidades debido a la participación del grupo hidroxilo del fenol y en el anillo pirrólico del indol. En todo caso, los resultados muestran que la interacción se va haciendo más intensa a medida que pasamos de benceno, a fenol y a indol. El análisis de la naturaleza de la interacción mediante cálculos SAPT(DFT) revela los motivos de estas diferentes estabilidades. Como era de esperar, la contribución más significativa a la estabilidad de los complejos proviene del término electrostático, pero se registran importantes contribuciones de la inducción y la dispersión a la estabilidad de los complejos. La contribución electrostática aumenta al pasar de benceno a fenol o indol debido a que se trata de especies polares. Asimismo, la inducción aumenta ya que se trata de especies más polarizables. Por último, en complejos con indol se registra un aumento significativo de la dispersión asociado a la mayor extensión de la nube aromática. Por lo tanto, el incremento de estabilidad observado se debe a incrementos en todas las contribuciones estabilizantes. La diferencia de estabilidades entre vi Resumen los complejos con diferentes especies aromáticas no varía significativamente incluso con la inclusión de moléculas de agua. Sin embargo, estructuralmente adoptan patrones de interacción muy distintos, aumentando también drásticamente el número de estructuras posibles a medida que vamos introduciendo moléculas de agua. La estabilidad que se gana al introducir una nueva molécula de agua es aproximadamente la misma con las tres especies aromáticas, así que la diferencia de estabilidad de los complejos se debe principalmente a las diferencias de intensidad de la interacción guanidinio···π. Interacción catión···π y efectos de la microhidratación en complejos formados por el catión pirrolidinio y especies aromáticas de las cadenas laterales de aminoácidos aromáticos En el capítulo 6 se detalla este tercer y último estudio de microhidratación, en el que al igual que en los dos anteriores, se introducen gradualmente moléculas de agua, hasta un máximo de tres, a los sistemas formados por el catión pirrolidinio y diferentes especies aromáticas para evaluar de nuevo los efectos provocados sobre la interacción catión···π. El catión pirrolidinio se escogió siguiendo el mismo criterio que para los trabajos anteriores: es una especie fácilmente encontrada en los sistemas biológicos y en particular en el aminoácido prolina. Así, con este estudio se cerraría la sección dedicada a la microhidratación en la que se han tenido en cuenta algunos de los cationes más básicos pertenecientes a las cadenas laterales de los aminoácidos esenciales. De nuevo se usaron benceno, fenol e indol como especies aromáticas para caracterizar las interacciones catión···π. Al igual que se observó para el caso del catión guanidinio, la interacción catión···π es más intensa con el indol, seguida de la interacción con el fenol. Así, la interacción benceno···catión es la menos favorecida de las tres ya que tan sólo puede establecer un punto de interacción a través de su nube electrónica. Por el contrario, los dos puntos de anclaje disponibles en el fenol e indol hacen posible la formación de una gran variedad de estructuras diferentes. Estos resultados se obtienen con independencia del método de cálculo usado. Dicho comportamiento también viene reflejado por los valores obtenidos mediante SAPT(DFT), que indican que aunque la interacción para las tres especies aromáticas es principalmente electrostática, dispersión e inducción también tienen una contribución no despreciable, principalmente para los complejos con indol. La interacción con fenol ya posee un carácter electrostático y de inducción más acusado que con el benceno debido a su grupo hidroxilo y su momento dipolar permanente. En el caso del indol, aunque también hay una importante contribución de la inducción, la contribución de la dispersión se hace visiblemente notable gracias a su sistema aromático extendido. Las estructuras más estables observadas suelen presentar la interacción catión···π entre el grupo NH del pirrolidinio y las nubes electrónicas de los sistemas aromáticos. Sin embargo, hay Resumen vii algunas excepciones en los complejos con fenol, en el que varias de las estructuras más estables presentan interacciones NH···O con el hidroxilo del sistema. Además, también se observan interacciones secundarias CH···π en bastantes casos en los complejos con fenol e indol, adoptando el catión en estos casos disposiciones paralelas respecto del anillo aromático. La introducción sucesiva de moléculas de agua en cada uno de los sistemas provoca una gran variedad de mínimos y en muchos casos con energías de coordinación muy similares, lo que indica que diferentes disposiciones dan lugar a una estabilidad similar. Además, aunque para el caso del benceno el patrón de interacción permanece más o menos simple, para los complejos con indol y fenol se hace más complejo. Sin embargo, sí que se puede observar que en la mayoría de los casos se conserva un patrón cíclico en el que están involucrados ambos puntos de anclaje para fenol e indol, con los que el catión normalmente se encuentra interaccionando, salvo alguna excepción. La inclusión de moléculas de agua también afecta a la estabilidad relativa entre los complejos de fenol e indol. Mientras que en ausencia de moléculas de agua los complejos con indol son claramente más estables que los de fenol, ya con dos moléculas de agua, las diferencias de estabilidad son prácticamente inexistentes. Finalmente, es interesante mencionar también que entre los diferentes métodos de cálculo usados, los que proporcionan una mejor relación precisión/coste computacional son el M06-2X/6-31+G* y el SCSN-MP2/AVDZ, siendo el primero especialmente indicado para los complejos no hidratados. Interacción entre el cation guanidinio y aminoácidos aromáticos Una vez concluido el estudio de microhidratación se propuso aumentar el tamaño del sistema con la inclusión de aminoácidos completos (Trp, Phe y Tyr), de modo que se pudieran evaluar las características de la interacción catión···π en un sistema más complejo. Se han tenido en cuenta tanto las estructuras zwitteriónicas como las que tienen los grupos amino y ácido en forma neutra. Para este trabajo se volvió a escoger el guanidinio como catión debido a su relevancia y a su peculiar estructura plana, que abre la posibilidad para la formación de diferentes tipos de interacciones. Los resultados indican que las interacciones entre el guanidinio y los aminoácidos Tyr y Phe, tanto en su forma zwitteriónica como neutra, son muy similares en estabilidad, mientras que la interacción con el Trp es algo más intensa. En todos los casos, los complejos más estables presentan interacciones cation···π, si bien se han localizado complejos con la forma zwitteriónica del amino ácido que son prácticamente tan estables como los formados con la forma neutra. El análisis de los distintos componentes de la energía de interacción muestra que los complejos con el amino ácido plegado de modo que permita el contacto catión···π son los más favorecidos debido a que el contacto catión···π aporta mayor es tabilización en forma de dispersión e viii Resumen inducción. Este aporte extra de estabilidad debido al contacto catión···π es capaz de superar el coste energético necesario para plegar el amino ácido, que es menor en otras estructuras más extendidas. Por lo tanto, para describir correctamente este tipo de interacciones también es necesaria una evaluación apropiada de las distintas contribuciones a la estabilidad, ya que la estabilidad final resulta del balance de contribuciones que a menudo operan en sentidos opuestos. En este trabajo se ha probado una batería de métodos de cálculo de diversa naturaleza, proporcionando todos ellos resultados bastante apropiados. En todo caso, los mejores resultados se obtienen con el método MP2.X, que es el único capaz de reproducir los valores de estabilidad relativa entre formas zwitteriónicas y neutras de los complejos. Respecto a otros métodos menos costosos, tanto M06-2X, B3LYP-D como SCSN-MP2/aug-cc-pVDZ son las mejores elecciones si tenemos en cuenta la relación precisión/coste computacional. Interacción entre el cation imidazolio y aminoácidos aromáticos En este trabajo se ha estudiado la interacción entre el catión imidazolio y la serie de aminoácidos aromáticos, incluyendo la histidina. El catión imidazolio es un anillo pentagonal que contiene dos grupos N-H. Se trata por tanto de un catión plano que puede interaccionar favorablemente mediante stacking con los diferentes anillos aromáticos de los aminoácidos. El procedimiento seguido es similar al planteado para el estudio de los complejos con catión guanidinio, empleando los métodos que ya han mostrado ser capaces de proporcionar una descripción adecuada de este tipo de sistemas. Los resultados indican que los complejos formados por el catión imidazolio con Phe, Tyr y Trp presentan características muy similares entre sí, y a la vez coincidentes con las obtenidas para los complejos con catión guanidinio. Las estructuras más estables presentan contactos catión···π, aunque complejos con los aminoácidos en forma zwitteriónica presentan estabilidades muy similares. Nuevamente, la estabilidad de los complejos aumenta al aumentar el tamaño del anillo aromático implicado en la interacción debido a que proporciona mayores contribuciones de inducción y dispersión, que compensan el coste necesario para plegar el aminoácido. Los complejos con His presentan un comportamiento totalmente diferente. No se encuentran estructuras con contactos cation···π entre los mínimos más estables, que corresponden a estructuras con el aminoácido en forma zwitteriónica y con el catión interaccionando con el aminoácido sólo mediante enlaces de hidrógeno. Los resultados indican que, tanto en fase gas como en presencia de disolvente, es posible la formación de estructuras en las que se establecen contactos catión···π con Phe, Tyr y Trp, mientras que ese tipo de estructuras no es posible para los complejos con His. Resumen ix En conjunto, podemos decir que los trabajos presentados en esta tesis ayudan a profundizar en el conocimiento de las características de los contactos catión···π que aparecen de forma habitual en sistemas biológicos. De modo general, los resultados indican que la presencia de un pequeño número de moléculas de disolvente que puede interaccionar con el contacto catión···π puede tener un impacto considerable en sus características. Así, la interacción se debilita, pero más relevantes son los cambios en la estructura de los complejos que se deben a la necesidad de acomodar las moléculas de disolvente próximas. Por otra parte, los contactos catión···π considerados están dominados por una combinación de efectos electrostáticos y de inducción, si bien el papel de la dispersión también puede ser relevante, especialmente en los cationes más estructurados. La preferencia de los aminoácidos por estructuras con contactos catión···π se debe a mayores contribuciones de la inducción y especialmente dispersión a la estabilidad, capaces de vencer la resistencia del aminoácido a adoptar conformaciones más plegadas. 1 Introduction 1.2. Non covalent interactions 9 mechanics origin, arising from quantum-induced instantaneous polarized multipoles in molecules. In symmetrical molecules like hydrogen there does not seem to be any distortion in the charge distribution to produce dipole moments. However, that is only true on average. The movement of the electrons in a molecule makes the charge distribution fluctuate continuously, and at any one instant, the electronic distribution might be not symmetrical, arising then an electric moment. This instant dipole can induce another dipole in a surrounding molecule, and thus, they can interact with each other. This synchronized movement of electrons, the so-called electron correlation, can occur even over a large number of molecules as long as the molecules are close together, spreading thus the dispersion effect. Indeed, London forces are the only long-range interactions present in all molecular interactions, being the most relevant contribution except in small polar molecules. The correlation effect favors lower-energy configurations, decreasing then the energy of the system. Also, the resulting effect is attractive since correlation becomes stronger when molecules are close each other. The charge fluctuation in a molecule that causes correlation is related to the polarizability. In this way, multipoles, either instantaneous or induced, are directly related to molecular polarizabilities and consequently dispersion interactions are too. The size and molecular shape also affect the dispersion interaction, since polarizability depends on them, dispersion forces usually increasing with the size of the molecules. The last effects included as part of the long-range forces are resonance and magnetic effects. On the one hand, resonance interaction occurs when at least one of the molecules is in a degenerate state, or if the system is formed by identical molecules, when one of them is in an excited state. Therefore, it cannot be observed in closed-shell molecules in their ground states. On the other hand, magnetic effects occur when there are unpaired spins. Anyway, both of them are negligible for the systems in this study. Short-range effects, contrary to long-range effects, need a more complex description, being the most important contribution the so-called exchange-repulsion effect.21 If the distance between nuclei is short enough, the overlap of the charge clouds becomes appreciable and large distortions are required by the Pauli Exclusion Principle. Thus, it is not possible to treat the system as composed of completely separated molecules in their respective unperturbed states, as for long-range forces. To satisfy this principle, a correction arising from the antisymmetry of the wavefunction is included in the wavefunction, the exchange term. Although the origin of some distortion is due to the Coulombic repulsion, this seems to be a secondary effect. For the systems treated in this work, with closed electronic shells, the electronic clouds of the interacting atoms tend to avoid each other. The decrease on charge density between the interacting atoms and therefore the reduction of the screening of the nuclear charges by the electrons, results in a 10 Chap. 1. Introduction repulsion effect between the nuclei and in the last instance, between molecules. However, although molecules at short distances mostly suffer this repulsive effect, this repulsion is also combined with an attractive force coming from the fact that the electrons, instead of moving around only one atom, can increase its free movement over both atoms (or molecules). At short distances, charge transfer and charge penetration have to be considered too. The former effect was firstly introduced by Mulliken, and consists on an attractive force between an acceptor (molecule deficient on electrons) and the electronic cloud of a donor (molecule rich on electrons). This effect may take place in the region between both molecules and greatly increases the interaction energy, though it can also be considered as a part of induction forces acting at short-range.13, 17 The latter effect arises from the overlap of electron clouds, and corresponds to the difference between the electrostatic and coulombic energies, which start to be different at short distances, since the multipole expansion does not provide a satisfactory description at these distances, resulting in an error. Thus, charge penetration can be seen as the attraction between the semi-shielded nuclei of one molecule and the electronic cloud associated to another molecule, making this combination unviable for description by using the usual multipolar expansion. Also, the magnitude of penetration, although attractive, decreases exponentially with increasing distance. Thus, different contributions can be defined at short distances depending on the method applied, though the overall effect is repulsive. Taking into account the cation···π interactions studied in this work, it can be said that the main contributions will be electrostatics, induction and, in larger systems, also dispersion forces. On the one hand, the interaction between the charge of the cation present in the system and the multipole distribution of the other species will be mainly electrostatic. On the other hand, the polarizable aromatic cloud of the aromatic species will be affected by the polarizing ion, leading to a high induction effect. 1.3. Interactions with aromatic systems When two species of the same or different nature are in contact, interactions between them appear and new bonds can be formed or not. These latter interactions, depending on the involved species, could belong to one or more kinds and could be of different strength. Thus, interaction between ions is purely electrostatic, while for instance in the interactions between polar species there are both electrostatic and induction contributions. In Table 1.1 different kinds of interactions as well as their approximate energy value in a simple contact between two species are shown. 1.3. Interaction with aromatic systems 11 Table 1.1. Nature and strength of different non covalent interactions. Interaction Binding Energy (kcal mol-1) Nature of the Interaction Ion-Ion 25 – 85 Electrostatic Ion-Dipole 10 – 50 Electrostatic and induction Dipole-Dipole 1 – 10 Electrostatic and induction Hydrogen Bond 1 – 30 Electrostatic and induction (dipole-dipole) XH···π 1-5 Electrostatic and induction (weak HB) Cation···π 1 – 20 Electrostatic and induction π···π Stacking 1 – 10 Weak electrostatic and dispersion Anion···π 5 – 10 Electrostatic and induction Lone pair···π 1 – 5 Electrostatic Halogen Bond 1 – 45 Weak electrostatic van der Waals Forces < 1 Dispersion However, it is important mentioning that frequently not only isolated interactions are established, but several of them jointly contribute to the stability of the system, greatly affecting their properties. A remarkable case is that one related with the DNA chain structure, stabilized by many stacking interactions together with hydrogen bonds between nitrogen bases. Nevertheless, the present work is focused on interactions involving aromatic systems, especially cation···π interactions. The presence of the delocalized electronic cloud of the aromatic units, together with their planar structure, confers to these systems some characteristic properties favoring them to be present in lots of structures and key processes of biochemistry, as protein-ligand interactions, neurotransmitters or ionic channels.22, 23 For that reason, the design of new drugs based on this kind of interactions has also gained importance.22-24 As already mentioned in previous sections, amino acids are fundamental structures whose importance relies on the fact that they are the constituent units of proteins. Some of them present aromatic units in their side chains, so it is also possible to establish interactions through them. The way in which these aromatic units interacts and their strength will be determined by the interacting species as well as by the environment in which the interaction is established. 12 Chap. 1. Introduction Therefore, each aromatic amino acid has a characteristic aromatic unit which could be used in different aspects, as shown in Figure 1.3. These amino acids can interact with other species by means of the aromatic rings in their side chain, so it is crucial to understand how aromatic species interact. Figure 1.3. Aromatic amino acids and their side chains. Commonly, an aromatic ring can establish interactions with different species, which in most cases can be classified as one of the following three types: π···π interactions, XH···π interactions and ion···π interactions. These three types of non covalent interactions will be briefly described below.25 Stacking or π···π interactions In aromatic systems, the electron density is located into π orbitals associated with the carbon framework.26 On the other hand, the σ electronic cloud, corresponding to the bonds with the peripheral hydrogens of the aromatic system, is positively polarized. Thus, two different regions coexist in the aromatic systems: an electron poor region (σ cloud) and a rich electron region, associated with the π system. This model was firstly proposed by Hunter and Sanders (1990) based on the combination of the van der Waals and electrostatic forces with the aim of explaining the different configurations observed and to predict their interaction energies.27 Each one of the previous regions is able to interact with different kind of species. For instance, if another aromatic system is in the neighborhood, different situations may occur: a perpendicular interaction between σ and π electronic clouds (T-shaped or Edge to face); or parallel interactions. In this latter situation it is possible to observe two different orientations: one with parallel displaced systems and another one with aromatic rings totally overlapped (face-toHis (imidazole) Phe (benzene) Tyr (phenol) Trp (indole) 1.3. Interaction with aromatic systems 13 face or sandwich), although only few examples exist for the sandwich arrangement.28, 29 The most common situations are the T-shaped or the parallel displaced ones, where molecules, although lie above each other, are offset enabling the complementary electron poor and rich regions to match up. Figure 1.4 shows the different types of interactions between two aromatic systems, benzene in this example, and how parallel arrangement competes with other dispositions.30, 31 Figure 1.4. Types of π···π interactions between benzene molecules. Different factors contribute to the strength of these interactions, such as the presence of substituents in the rings like electron-donors or electron acceptors, the ring being or not a heterocycle, or the number of condensed aromatic rings involved in the structure. Thus, the strength of the interaction tends to increase when electron-attractor groups are present on the ring, since the π···π repulsion decreases because this substituent attracts the electronic cloud, diminishing its density.32, 33 Also it increases with the higher number of rings involved34 or when the rings are heterocycles.35 In the latter case, as in the example of the electron-attractor substituents, the heteroatom, more electronegative than the carbon atom, causes a decrease of the electronic π density. This is a common effect when a nitrogen atom is substituting a carbon atom of the six-membered ring, as in the case of pyridine.35 Anyway, a recent review of this interpretation has been done by Wheeler pointing towards direct through space effects to explain substituent effects.36, 37 Although stacking interactions were considered to be energetically much less significant than H-bonding (0-10 kcal mol-1), accurate calculations have revealed that when they are associated, Offset face to face or Paralleldisplaced Edge to face or T-shaped Faceto face or Sandwich 14 Chap. 1. Introduction like in DNA, surprisingly large stabilization energies are observed comparable with those of strong H-bonding. It is widely known that stacking interactions together with hydrogen bond play a key role in DNA stabilization as well as in its characteristic helical structure.22, 23 This fact makes stacking interactions of great importance in biology and supramolecular chemistry, being not surprising the large number of drugs which has been developed based on them.38-40 For instance, in cancer therapy a new target has emerged, the so-called c-Met or hepatocyte growth factor receptor (HGFR). This receptor tyrosine kinase has an abnormal activation in many cancer cells,41, 42 and selective c-Met inhibitors have been developed in which π···π stacking interactions have an important role in the binding to the target. Furthermore, with the discovery of fullerenes and related carbon-rich materials, a novel type of stacking involving curved carbon networks has been introduced. This kind of materials with accessible concave surfaces appear to be good candidates for the formation of π···π stacked supramolecular assemblies with fullerenes. The buckycatcher synthesized by Sygula el al,43 which acts as a receptor of fullerene C60 molecule, is an outstanding example of this kind of interactions and its applicability. Also, since carbon nanotubes (CNTs) have been widely utilized as novel drug carriers, the π···π stacking interaction has been used by some authors as a way to bind the drug to the CNTs structure.44 XH···π interactions In this interaction a π electron cloud of an aromatic ring and an X-H group are involved, usually C-H, N-H or O-H groups. CH···π interactions were described in 1952 by Tamres, when he observed that benzene and analogous compounds dissolve exothermically in chloroform.45 In 1957, Reeves and Schneider proposed this interaction as a type of H-bond, taking into account their NMR results. Since then, CH···π interactions have been described in a vast number of small molecular systems, being Nishio’s book an excellent compilation of these observations.46, 47 Interactions between π systems and NH or OH bonds have been also described with similar characteristics, although quantitatively the strength of the different interactions follows the order OH···π > NH···π > CH···π.24 Spectroscopic measurements of benzene···water and benzene···ammonia complexes reported by Suzuki et al.48 and Rodham et al.49 showed that benzene acts as proton acceptor and that water and ammonia are positioned above the benzene plane. The NH···π interaction between a NH bond and a π-system was first reported in 1959 by Oki and Imamura using the measurements of the IR spectra of N-benzylaniline and its derivatives,50 although it was also observed in other systems,51-56 being a remarkable example the presence of this interaction in the hemoglobin-drug complex, reported by Perutz et al.57 The preferred location of the amino groups with respect to the aromatic ring was reported by Burley and Petsko, indicating that these groups tend to locate above and below the aromatic system and 1.3. Interaction with aromatic systems 15 close to the center.58 On the other hand, the OH···π (or OH/π) interaction plays an important role since it allows the interaction between non polar groups as aromatic systems and the polar water molecules of the solvent.59 This kind of interaction is the responsible of the fact that benzene, a typical nonpolar organic molecule, is soluble in water, although also the delocalization of the π cloud greatly contributes to the affinity towards the polar solvent.60 In XH···π contacts, the main contributions to the strength are the charge transfer and the electrostatic effects, which mainly control the directionality of the interactions, giving rise to a linear interaction between the CH bond and the p orbitals of the ring. However, dispersion effects also contribute in a non despicable way, suggesting that XH···π contacts have different nature than prototypical hydrogen bonds.61-64 Recent experimental and theoretical studies in the gas phase65, 66 show that the main contribution to this interaction is in many cases due to dispersion effects instead of electrostatics ones, as in the typical hydrogen bonds. The electrostatic contribution only becomes dominant when the hydrogen involved in the X-H···π interaction is very acid, for example, in cases as chloroform or acetylene. This fact makes the X-H···π interaction stable in both polar and non-polar solvents. Besides, the different dependence against the directionality of the interaction, corroborates their different nature.67 While typical hydrogen bonds have a high directionality, the X-H···π interactions are not that strict with respect to the orientation of the species involved in it. In fact, it is not rare to find non linear arrangements on X-H···π contacts.24, 68 The overall protein stability is in many cases not more than a few kcal mol-1, so that the overall stabilization energy of about 0.5-1.0 kcal mol-1 per interaction can be an important contributor.67 In fact, the formation of many complexes of proteins with special ligands or cofactors in which X-H···π interactions have a special participation have been described by different authors.46, 69-73 An example of applicability of this interaction was the design of serine protease inhibitors.74, 75 It is worth mentioning that substituent effects on XH···π interactions can impact deeply on their properties, but their magnitude as well as their physical origin depend strongly on the nature of the X group. Particularly, the higher polarity of the XH bond, the higher the effects caused by the substituents will be.76 XH···π interactions will be of great importance in this work, since the aromatic units of the aromatic side chains can interact by means of this kind of interaction with solvent molecules or the XH residues present in the molecules (intra or intermolecular interactions depending on the case). 16 Chap. 1. Introduction Cation···π interactions Aromatic systems such as benzene or phenol are commonly found in different biological or nonbiological contexts. In the biological context this kind of systems can be observed for instance as a part of aromatic amino acids as tyrosine, tryptophan or phenylalanine.26 The existence of charged species surrounding them, even forming part of the same protein (like cationic side chains of arginine or lysine) allows for the possibility of interactions between the charged groups and the electronic cloud of a π system. A particular and usual case is that one in which a cation interacts with the delocalized electron density that lies above and below the plane of an aromatic ring. In fact, an important role of this kind of interactions is that one related with the building of the correct tertiary structures in proteins. In Figure 1.5 examples of different aromatic amino acids interacting with a guanidinium cation, present in the side chain of arginine, are shown. Cation···π interactions can be appreciated acting together with hydrogen bonds. Figure 1.5. Examples of cation···π interaction. The cation···π interaction is usually explained in terms of electrostatic effects.77-80 Accordingly, this kind of interaction arises from the electrostatic interaction between the negative first nonzero multipole moment of the arene (quadrupole moment in the case benzene), and the positively charged cation. Other molecules like ethylene or acetylene and its derivates also have a quadrupole moment, which can interact in the same way by means of cation···π interactions. Although this is the larger contribution, other effects are also important, as the polarization of the π-electron system by the cation.81 However, the dispersion contribution is typically small.82 The cation···π interaction usually has a large magnitude, more appreciable in the gas phase, as one might expect. For instance, the first evidence of this cation···π interaction, published in 1981 Gnd-Phe Gnd-Tyr Gnd-Trp 1.3. Interaction with aromatic systems 17 by Kebarle,83 reported that K+ ions are better stabilized in the gas phase by benzene (ΔHº = -19 kcal mol-1) than by water (ΔHº = -18 kcal mol-1). However, contrarily to other strong non covalent forces like ion pairs or hydrogen bonds, cation···π interactions retain the strength and the specificity in a better way in water (and even across a range of solvents), and they can be even strong in this environment.84 For instance, since benzene is a hydrophobic nonpolar molecule, the first step of desolvation is not as big a problem. Therefore, the global stability will be dominated by the desolvation energy of the cation and also by the tendency of the cation to bind the benzene molecule. This fact is the responsible of the differences in the cation···π interaction tendencies when they are measured in gas or aqueous phases. Although Li+ cation has the strongest cation···π interaction in the gas phase among all the alkaline cations, when the measurement is done in water, it is the K+ cation the one presenting the largest tendency to interact with benzene. K+ cation has much lower desolvation energy than Li+ while still remains as a good π binder. Consequences of these properties are shown in reports which indicate that most cation···π interactions in proteins are located at the surface of proteins, where salt effects would be the weakest, pointing again to the importance of the cation···π interactions on protein structure.85 Due to their remarkable properties, the interest on cation···π interactions rapidly increased, giving rise to a wide amount of experimental and computational works focused on their study. Thus, only a few years after Kebarle’s report, Meot-Ner demonstrated that the more complex organic cations like ammonium or tetramethylammonium (TMA) are also good πbinders in the gas phase.86 Besides, ab initio studies of this kind of interaction have been done to obtain the theoretical minimum energy structures.87 Cation···π interactions are also present and take part on numerous structures and biological processes. For instance, cation···π interactions are common in proteins or their complexes with different ligands or even with DNA.78, 80, 87, 88 Also, cation···π contacts usually participate in the bond between neurotransmitters and receptors. A remarkable example of this is the nicotinic receptor. This protein usually binds the acetylcholine (Ach) neurotransmitter, triggering a process in which the voluntary movements of muscles are involved. This binding is established between the quaternary ammonium of the Ach and one specific tryptophan of the receptor by means of a cation···π interaction. However, the nicotine molecule, protonated in the biologic system, shows a similar structure as the neurotransmitter, with a quaternary ammonium being able to bind the receptor also by means of a cation···π interaction, resulting in a physiological response.89 Another example, postulated by Dougherty in the early 1990’s,90 was the possibility of the participation of cation···π interactions in highly selective K+ channels. As commented above, among the alkalines, K+ shows the larger cation···π interaction in water, allowing the channels to be more selective for K+ than for Na+. Once proved the efficiency of the interaction between an 18 Chap. 1. Introduction aromatic ring and a cation due to its abundance in many kinds of systems, this interaction has also been used to increase the π-face selectivity in asymmetric catalysis87, 91 or as binding places for alkaline metals. In the latter case, molecules similar to crown ethers but with aromatic systems located in the structure have been used.87, 92 Regarding to substituent effects in cation···π interactions, it is worth mentioning that not only electrostatic forces are involved, but also polarization effects on the aromatic systems play an important role.93 Thus, both effects are important for an appropriate description of the cation···π interaction, either for heterocyclic arenes or not. The size of the aromatic systems greatly affects the polarizability, leading to the increase of the induction effects. In fact, the induction effects would be enough to stabilize cation···π complexes even in cases in which the electrostatic term is repulsive.94 Substituents also influence the polarizability of the arenes and in general, it can be assumed that while electron-donor (π-electron rich) groups strengthen the cation···π interactions, the electron-acceptor (π-electron poor) groups weaken it, but not only from the local C–X dipoles, but also due to the resonance effects over the aromatic ring. However, these substituent effects can also be seen from another point of view, considering that these effects arise from interactions of the cation with the local dipole (or higher-order multipoles) associated with the substituent.95 Although it is not of relevance in this work, it is worth mentioning that interactions between anions and aromatic rings are also possible. Due to the great capability of substituent effects, when an aromatic system is completely substituted with electron-acceptors (for instance halogen, nitro or cyano groups), the π cloud becomes electron-deficient so that attractive interactions with anions is possible.96, 97 1.4. References (1) A. V. Finkelstein, O. B. Ptitsyn, Protein Physics: A Course of Lectures, Academic Press, London, 2002. (2) G. E. Schulz, R. H. Schirmer, Principles of Protein Structure, Springer, Berlin, 1979. (3) M. P. Callahan, K. E. Smith, H. J. Cleaves, J. Ruzicka, J. C. Stern, D. P. Glavin, C. H. House, J. P. Dworkin. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 13995-13998. (4) D. P. Glavin, J. L. Bada, K. L. F. Brinton, G. D. McDonald. Proc. Natl. Acad. Sci. U. S. A. 1999, 96, 8835-8838. (5) J. R. Cronin, S. Pizzarello. Science 1997, 275, 951-955. (6) J. E. Elsila, D. P. Glavin, J. P. Dworkin, Z. Martins, J. L. Bada. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, E3288. (7) J. L. Bada. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, E85. Chap. 2. Objectives 25 The interaction between cations and aromatic species in the side chains of some amino acids is one of the key factors that control the stability of proteins. Different studies suggest that this type of interaction is present in a variety of systems indicating its relevance in many processes of chemical and biological recognition, such as nerve transmission or transport through the membrane. Considering the relevance of these interactions, the main objective of this thesis is to gain a deeper comprehension of the characteristics of systems that establish interactions between aromatic and cationic species by applying computational chemistry methods. It is already know that the interaction between cations and aromatic systems is strong in the gas phase as a consequence of large contributions from electrostatics. The electrostatic contribution mainly comes from charge-dipole and charge-quadrupole contributions, and it usually dominates the interaction in model systems employing small cations and aromatic species. However, the aromatic electron cloud is easily polarizable, so large contributions from induction have been also observed, especially when large aromatic systems are involved. Therefore there is a general agreement as to consider the cation···π interaction as controlled b y a combination of electrostatic and induction. However, the data found in the literature show a significant variability with regard to the intensity of such interactions and thus the importance they actually have for protein stability. From some of these results, the interaction between the cation and one of the aromatic amino acids is intense and can make a decisive contribution to the stability of the system. On the other hand, other studies indicate that the role played by such interactions is marginal, contributing only a few tenths of kcal mol-1 to the stability of the system, although such estimates are often purely qualitative. Overall, everything seems to indicate that the environment of the interacting species may affect the characteristics of the interaction, either causing changes in intensity or changes in the geometrical disposition of the interacting fragments. This thesis proposes conducting a thorough study of this kind of interactions in systems of interest, trying to clarify essentially which are the main features of such interactions. Therefore, it also aims to quantify the role played by the environment on these interactions, and whether it can act as a modulator of them, either weakening or causing significant structural changes. Since the study of cation···π interactions is a wide field for study, this thesis will be focused on interactions involving aromatic and charged species that can be on the side chain of amino acids. Since the problems considered in the thesis are quite complex for their treatment using quantum mechanics methods, a gradual approach to the problem has been chosen. Thus, a series of initial studies has been carried out using simplified models for the amino acids where only the aromatic ring is considered in the interaction. Using these simplified models, studies have been conducted in order to ascertain the role played by the solvent molecules closer to the cation···π 26 Chap. 2. Objectives contact upon its characteristics. These studies using simplified models for the aromatic amino acid are the subject of chapters 4 to 6. Chapters 7 and 8 of this thesis show the results for a more complex model using the whole amino acid structure in the calculations. The goals of the different chapters are summarized below. Chapter 4 is devoted to the simplest of the systems considered, corresponding to complexes formed by phenol with ammonium and methylammonium cations. Therefore, the systems mimic the interaction between a phenol unit in the side chain of tyrosine and the cationic end of lysine. Most studies in literature on cation···π interactions deal with benzene complexes employing simple cations as alkali metals. The hydration pattern of an alkali cation···benzene complex is pretty simple, the water molecules coordinating to the cation while a coordination location is occupied by the benzene unit. However, the situation can be already rather different when a slightly different aromatic species such as phenol is employed. The presence of the hydroxyl group makes another position available for hydration thus leading to different hydration patterns. A microhydration study has been carried out, so a small number of water molecules has been added stepwise to the cation···π in order to ascertain their impact on the structure and energetics of the cation···π moiety. The results will help to determine the main characteristics of the competition between the aromatic cloud and the hydroxyl group in order to interact with both the cation and the water molecules. Also, the role played by the water molecules closest to the cation···π contact, and especially their impact in structure and stability, will be revealed. Phenol Methylammonium Ammonium CationsAromatic ring 4 Chap. 2. Objectives 27 Chapter 5 deals with complexes formed by guanidinium cation and benzene, phenol and indol, in an attempt to recover information on the interaction between the cationic end of arginine and the aromatic amino acids phenylalanine, tyrosine and tryptophan. The procedure to follow is similar to the one employed in chapter 4. Thus, a small number of water molecules is included in an stepwise manner in the complexes formed by guanidinium cation and the three aromatic species. The results obtained will allow a direct comparison of the interaction of the three different aromatic molecules with guanidinium, giving hints about its strength. Also, guanidinium can interact with the aromatic rings adopting both a perpendicular and stacked orientation. While the perpendicular complex is the most stable in the gas phase, it has been indicated that in benzene complexes the presence of water molecules can change the order of stability favoring the parallel arrangement. The results will give more information about this competition between different orientations for the three aromatic units considered, as well as the possible role of water molecules in changing the stability order. Chapter 6 further complicates the cationic unit employed, which in this case corresponds to pyrrolidinium cation, the protonated end of proline. Pyrrolidinium cation shows an ammonium group within a saturated five-membered ring, exhibiting larger flexibility than any of the cations considered in previous chapters. Besides, due to the ring structure of the cation it could be expected that parallel stacked structures will be favored. Reaching this point a question arises as to whether the methods employed in previous chapters are capable of properly describe the systems studied. In principle, the cation···π interaction, controlled by electrostatics and induction, could be described with common DFT methods. However, the possibility of forming stacked structures with larger contributions from dispersion suggests that more appropriate methods should be employed. Therefore, complexation energies will be obtained at the CCSD(T)/CBS Benzene Phenol Indole Guanidinium CationAromatic rings 5 28 Chap. 2. Objectives level, being used as benchmark in order to check the performance of cheaper methods. Therefore, in this chapter the goals are twofold: determining the characteristics of the interaction and the role of water molecules, as in previous chapters, and checking the performance of different methods in this kind of complexes. In chapter 7, whole aromatic amino acids are considered, focusing on their interaction with guanidinium cation. Considering the whole amino acid introduces the problem of the large number of conformers which can form complexes with guanidinium cation and similar stability. Therefore, an automatic procedure for searching the conformational space of the amino acids and their complexes with guanidinium is employed based on an empirical force field. The results are further refined by using different quantum chemistry methods. As in the previous chapter, benchmark-quality results have been obtained for the stabilities in order to check the performance of cheaper methods. The results obtained will give information on the interaction and, specially, of the role of the amino acid group trying to interact with the guanidinium cation, thus competing with the cation···π contact. Benzene Phenol Indole Pyrrolidinium CationAromatic rings 6 TryptophanPhenylalanine Tyrosine Guanidinium CationAromatic amino acids 7 Chap. 2. Objectives 29 Finally, in chapter 8 a similar study employing the whole amino acids is performed, but in this case considering imidazolium as cation, trying to model the protonated side chain of histidine. Complexes of imidazolium cation with phenylalanine, tyrosine and tryptophan as in the previous chapter, as well as those formed with histidine are considered. An exhaustive exploration is performed for these complexes, employing the methods that have already proved to properly reproduce the interaction in chapter 7. Also, the effect of the solvent as modeled by a polarizable continuum model upon the characteristics of the complexes is evaluated. Histidine Tryptophan Phenylalanine Tyrosine Imidazolium CationAromatic amino acids 8 3 Methodology 3.1. Methods based on the wavefunction 33 3.1. Methods based on the wavefunction Macroscopic systems can be easily described by Newtonian mechanics, where the energy can apparently vary continuously. However, when a system has to be described at a microscopic level, a different mechanics is required, in which quantization is taken account: this is quantum mechanics.1 The fundamental postulate of quantum mechanics states that a so-called wavefunction exists, Ψ, for any state which describes the system, and that any observable property a of the system can be obtained by applying an appropriate operator upon that wavefunction. Thus, the goal of quantum mechanics is finding the system wavefunction. This postulate can be expressed as the following equation:  aA ˆ (eq. 3.1) This equation is called the Schrödinger equation when the operator is the Hamiltonian, and the property obtained is the energy of the system:  EH ˆ (eq. 3.2) Usually, the Hamiltonian operator without relativistic effects is taken into account, constructed from five contributions to the total energy of a system, being those contributions the kinetic energies of the electrons and nuclei; the attraction of the electrons to the nuclei; and the interelectronic and internuclear repulsions. In the process of obtaining the best wavefunction which describes our system, its quality can be judged by evaluation of the energy eigenvalues associated to it. As stated by the variational principle, the wavefunction with the lowest energy will be the most accurate and presumably the best one. This trial wavefunction is usually represented by a combination of appropriate functions (LCAO approach), the basis set, that are commonly atomic orbitals (AOs) in chemistry. A large number of methods based on the wavefunction have been developed. In the next sections some of those will be briefly exposed. 3.1.1. Ab initio methods The term ab initio comes from the Latin expression “from the beginning”. This kind of methods corresponds to those computational methods which draw upon theoretical principles, without relying in any experimental data to obtain the results. Quantum mechanics needs approximate methods to solve multiple electron systems, since the equations can be solved exactly only for 34 Chap. 3. Methodology one electron systems. Thus, ab initio methods, as other many ones, were developed for this purpose using, in most cases, mathematical approximations.2 In this section, the most popular ab initio methods and those ones used for obtaining the results in this present work are shown. 3.1.1.1. Hartree-Fock approximation The Hartree-Fock approximation (HF) is a central method to quantum chemistry, and though it does not include electronic correlation, it has an important role as starting point for more accurate approximations.3 This approximation takes into account the Born-Oppenheimer approach and also the non-relativistic Hamiltonian. Thus, the method can be described as the process of obtaining the best Slater determinant3 set up by the molecular orbitals constituting the N-electron wavefunction of the system, taking into account the proper antisymmetry which satisfies the Pauli exclusion principle. All of these molecular orbitals are formed by linear combination of atomic orbitals (one-electron functions), the coefficients obtained as a result of the Hartree-Fock method. According to Hartree-Fock theory, it is possible to approximate the wavefunction of the system through just a single Slater determinant formed by an specific set of N spinorbitals {χi} that results in the best approximation to the ground state of the N-electron system described by this Hamiltonian, those which minimize the electronic energy as the variational principle states. Thus, the simplest antisymmetric wavefunction, which can be used to describe the ground state of N-electron systems, is a single Slater determinant, NbaHF  ...... 21  . (eq. 3.3) Particularly, the closed-shell Hartree-Fock electronic energy can be obtained from eq. 3.2 by using the electronic Hamiltonian, resulting in a three-term expression. According to the variational principle, the best wavefunction is the one which gives the lowest possible energy HFHFHF HE   , (eq. 3.4) where H  is the full electronic Hamiltonian. By minimizing EHF with respect to the choice of spinorbitals, one can derive an equation, called the Hartree-Fock equation, which determines the optimal spinorbitals, and which is an eigenvalue equation of the form: aa f   ˆ , (eq. 3.5) 3.1. Methods based on the wavefunction 41 If all cluster operators up to were included in , the coupled cluster wavefunction would be equivalent to full CI because all possible excited determinants would be generated. But, as in CI, a truncation of the expansion is necessary. The advantage of coupled cluster over full CI appears when this truncation of the expansion is done. Truncation in CI at the i level implies that only states up to i excitations are taken into account, making the method non-size consistent. However, when the truncation in CC is done in a specific level of excitation, higher levels are taken into account due to the own definition of the cluster operator. For instance, if only the double-excitation operator is consider (CCD), = 2 the Taylor expansion would be: HF 3 2 2 2 2HF ˆψ... !3 ˆ 2! ˆ ˆ 1ψe        TT T T CCD . (eq. 3.22) Thus, as 2 is the double-excited operator, its products give the multiple of two excitations. The third term into the parenthesis corresponds to the quadruple excitations, the fourth corresponds to the hextuple excitations, etc, ensuring size consistency.5 If also the single excitations 1 are included, the method is called CCSD. The expression for the CCSD wavefunction is then described by the following expression:       HF 3 21 2 2121 HF ) 2 ˆ 1 ˆ ( Sψ ... ˆˆ !3 1 ˆˆ !2 1 ˆˆ 1 ψe                 TT TTTT TT DCC . (eq.3.23) Thus, the exponential operator may be arranged by groups with the same excited states, ... ˆ 24 1 ˆˆ 2 1 ˆ 2 1 ˆˆˆ ˆ 6 1 ˆˆˆˆ 2 1 ˆˆ 1e 4 1 2 12 2 2314 3 1213 2 121 ˆ                     TTTTTTT TTTTTTT T (eq. 3.24) In this way, it can be appreciated that the first term corresponds to the reference HF wavefunction and the second to all singly excited determinants. In the first parentheses, with the doubly excited determinants, two kinds of terms appear: connected ( 2) and disconnected ( 12), and in the rest of parenthesis it can be observed the same, with “true” ( i) or “product” excitations which result in an i excitation. Physically, a connected type such as 4, corresponds to 42 Chap. 3. Methodology four electrons interacting simultaneous, while a disconnected term such as 22, corresponds to two non-interacting pairs of interacting electrons. This CCSD method is the most common choice since it includes up to double excitation operators, and has a computational cost not too high.1, 2, 20 However, there are other orders of CC expansion as CCSDT21 and so on, in which higher excitation operators are included. The drawback is the huge computational cost they imply. Another option to estimate the effects of the connected triples is to use CCSD(T) calculations,22 in which the triple excitations are not included by adding the triple-excited operator but by means of a perturbation approach. Indeed, CCSD(T) calculations constitute a very robust method and one of the most used as standard.23 3.2. Improving the results Any calculation performed using a wavefunction-based method is obviously affected by errors. Apparently, these errors come from two different sources: 1) from the one-electron basis functions, the building blocks of the N-electron wavefunction. 2) from the method itself, and its difficulties for describing the correlation properly. In the first case, the error is related to the fact that an incomplete basis set is employed. The larger the basis set employed, the smaller would be the associated error. In the second case, methods approaching to FCI solution as MPn or CC need to be truncated to make the process viable. This truncation causes loss of information, carrying then an associated error compared to FCI results. By including higher orders of MPn or CC, the error would decrease, but the computational costs would increase considerably. Therefore, taking into account these sources of errors, the best way to reduce them is the choice of a good correlation method combined with a large basis set, always according to our computing resources. In practice, when accurate results are desired, it is usual to employ MP2 or CCSD(T), depending on the size of the system and the available resources. 3.2.1. Extrapolation to basis limit As commented in the previous section, the error associated to the one-electron basis set could be reduced by using larger basis sets. However, the continuous shift towards larger basis sets does not produce a convergence as quick as it would be desirable. Conversely, it converges very slowly as the size of the basis set is increased, being impractical or computationally expensive the solution of using very large basis set. 3.2. Improving the results 43 Besides this procedure, explicitly correlated methods can be employed, which tend more quickly to the limit than the standard ones.24, 25 An intermediate solution consists on extrapolating the results obtained with moderate sized basis sets, obtaining then an estimation of the limiting value.26-32 However, this last solution requires a smooth behavior of the energy (smoothly varying as the basis set increases) to construct a fitting function, and then approaching the limiting value. For example, the correlation consistent cc-pVXZ family of basis sets proposed by Dunning is a good choice for this purpose.33 Once the basis set is chosen, it has to be taken into account that the behavior of the energy can be split into two different behaviors, that one due to the HF energy and that one due to the correlation contributions to the energy. Firstly, regarding to the HF energies, it can be assumed an approximately exponential behavior in molecules as in the case of atoms.24, 34, 35 Thus, the following expression for the behavior of the HF energy is frequently employed: BXHF CBS HF XAeEE   , (eq. 3.25) where X is the cardinal number of the basis set. Thus, 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 B HF X HF X HF CBS EE EE eb b bEE E 21 11 ; 1          . (eq. 3.26) Although this method gives a good approximation, it is still computational expensive since it requires three calculations of increasing size. Other two-point extrapolation schemes have been developed, as the one proposed by Karton and Martin;36, 37 )9exp()1( XXAEE HF CBS HF X . (eq. 3.27) These errors associated to HF energies are not the main problem, since HF energies converge quite quickly with the basis set, especially considering energy differences, as interaction energies. Also, since errors due to correlation are usually larger, basis sets as cc-pVTZ seem to be enough for a proper description of the HF contribution. However, regarding to the correlation energies, it has to be said that the dependency on the basis set size is stronger and also that the convergence with the basis set is even slower. For that reason, larger basis set would be needed. The correlation energy usually behaves as: 44 Chap. 3. Methodology 3  AXEE XCBS . (eq. 3.28) The correlation energy can be obtained then by performing two calculations (since it contains only two unknowns): 3  AXEE XCBS , (eq. 3.29) 3  AYEE YCBS . (eq. 3.30) So, 33 33 YX EYEX EYX exact    . (eq. 3.31) The advantage of this scheme relies on the fact that employing lower order basis sets as cc-pVTZ and cc-pVQZ, better results are obtained than using directly a huge basis set as cc-pV6Z, reducing considerably the computational costs. 3.2.2. Benchmark values Usually, MP2 estimations are used instead of better ones as CCSD(T) to obtain values at the CBS limit for the correlation energy. For these procedures, CCSD(T) method requires computational resources that in cases of systems of moderate size are unviable. However, although MP2 is less demanding, also entails an error due to the N-electron model employed. Therefore, in order to estimate the CCSD(T)/CBS values but with a reasonable computational cost, some approaches have been devised, many of them by Hobza.38, 39 These methods proposed by Hobza and collaborators, are based on the fact that the differences between the correlation energy contributions to the interaction energy of a dimer using MP2 and CCSD(T) are pretty independent of the basis set size.18, 38, 39 Thus, it is possible to perform an approximation of the CCSD(T)/CBS correlation energy by using a correction with a small basis set (term in parentheses), since the difference between CCSD(T) and MP2 is assumed constant:  smallbasiscorr MP smallbasiscorr TCCSD CBScorr MP CBScorr TCCSD EEEE , 2 , )( , 2 , )(  , (eq. 3.32) and where the MP2 contribution to the correlation energy is estimated to basis limit with the extrapolation procedures described above. Although this procedure is usually employed to obtain 3.3. Density Functional Theory methods 45 benchmark values for the interaction energies of complexes, it has to be said that still entails errors associated to the small basis calculation.40-42 The following step within this approach is to avoid the calculation of the CCSD(T)/smallbasis value, still demanding for moderate size systems. Procedures for estimating the CCSD(T) corrections have been proposed by Hobza and Rezac (MP2.543 and the following improvements MP2.X44). The MP2.X approach consists of substituting the CCSD(T) calculations by cheaper MP3 ones, and also uses and empirical scaling coefficient. This coefficient is obtained by fitting to the CCSD(T)/CBS estimates of a set of complexes of different nature. The limit is estimated as:  smallbasiscorr MP smallbasiscorr MP CBScorr MP CBScorr TCCSD EECEE , 2 , 3 , 2 , )(  , (eq. 3.33) given a proper C coefficient in every case, it has been seen that the accuracy of the MP2.X results is almost independent of the small basis set employed. Besides, even with a basis set as small as 6-31G* the results are similar to those obtained from CCSD(T) calculations.45, 46 Thus, the computational costs are considerably reduced by using both a smaller basis set and a lower-level method. 3.3. Density Functional Theory methods The density functional methods arose from the idea of an interpretable or more intuitive method than those derived from the wavefunction.1 Thus, instead of studying the wavefunction of the system, these methods work with a useful physical observable in determining the energy of the system. This physical observable is the electron density ρ because it lets the construction of the Hamiltonian operator depending only on the positions and atomic numbers of the nuclei and the total number of electrons; integrated over all space, ρ also gives the number of electrons N of the system. The electron density is a good choice since as the nuclei are effectively point charges, the local maxima in the electron density correspond to these positions. Besides, the advantage of the electron density over the wavefunction for the calculations relies on the fact that, since the electron density is the square of the wavefunction, the electron density has the same number of variables when the number of electrons increase, so it is independent of the system size, making the electron density a useful tool for that purpose. The problem is then to find the exact functional which gives the relation between the density and the ground state energy of the system. DFT methods have been developed with this aim.47-52 There were some previous approximations to the DFT methods in which only the electron density was used as a variable to try to evaluate the molecular energy, knowing that it is 46 Chap. 3. Methodology separable into kinetic and potential components. Hohenberg and Kohn in 1964 were the ones who proved two critical theorems for establishing DFT as a legitimate quantum chemical methodology, and basically as a proof that the electron density ρ determined completely the ground state electronic energy. The first Hohenberg–Kohn theorem states that knowing the electron density of a stationary non-degenerate ground state of any system, any observable existing in this state can be calculated from it.5, 53 This is because observables like the energy are functionals of the electron density of the ground state of the system.1 Thus, the electron density determines the Hamiltonian (except for an additive constant) so there is a direct relationship between the density and the wavefunction through the external potential: ][][][][  eeNe VVTE  . (eq. 3.34) T[ρ] and Vee[ρ] do not depend on the external potential and can be included within the Hohenberg-Kohn functional FHK[ρ], with which ][)()(][  HKvFdvE  rrr , (eq. 3.35) where Ev[ρ] indicates that, for a specific external potential v(r), the energy is a functional of the density. On the other hand, the second Hohenberg–Kohn theorem is a variational theorem, which states that determining the density that minimizes the energy of the ground state it is possible to theoretically obtain, in an exact way, the electron density of a non-degenerate ground state. 3.3.1. The Kohn and Sham method With those two theorems in hand, the energy may be minimized to determine the density of the ground state.20 The main problem is the unknown expression of the relation between FHK and the density, in particular the T[ρ] form. Trying to solve this problem, Kohn and Sham54 proposed an ingenious method, similar in structure to the Hartree-Fock method, to calculate the energy from ρ, using a non-interacting electron system as a reference which, after being applied an external potential vs(r), provides a wavefunction, ψs, which has the same density as the real system. As commented above, Hartree-Fock and density functional theories have a very similar mathematical development, and thus, both need the Hamiltonian of the system to get the energy.5 However, in DFT methods an approximate Hamiltonian is employed instead of the exact one, since an ideal system is considered, without electron-electron interactions, providing a potentially 3.3. Density Functional Theory methods 47 exact result with the same density than the real interacting electron system. The Hamiltonian of this system only contains the single-electron terms, meanwhile the exact wavefunction is the Slater determinant, constructed from the so-called Kohn-Sham orbitals:    N is N i N i siviihH 11 2 1 )( ˆ )( 2 1 )( ˆ ˆ , (eq. 3.36) )()...3()2()1( ! 1 321 N NNs   . (eq. 3.37) These unknown Kohn-Sham orbitals, which minimize the energy, may be determined by numerical methods, or expanded in a set of basis functions, analogously to the HF method: iiis v         )( 2 12r ; ijji   . (eq. 3.38) Also, the exact electron density for this system is expressed as a linear combination of these orbitals, and it is used to calculate the energy,   occ N ii 1 2 )()( rr  . (eq. 3.39) Following this method, the energy functional of the real system may be divided into three parts (implicitly including correlation energy in all the terms): kinetic energy, T[ρ], and potential energy split into the attraction between the nuclei and electrons, Ene[ρ], and electron–electron repulsion, Eee[ρ] (the nuclear–nuclear repulsion is a constant within the Born-Oppenheimer approximation). Also, the Eee[ρ] term may be divided into Coulomb (J[ρ]) and exchange (K[ρ]) parts. On the other hand, this Kohn-Sham procedure proposes the division of the kinetic energy functional into two terms, one which can be calculated exactly, and a small correction term due to the difference between the exact kinetic energy and that calculated by assuming non-interacting orbitals, it is taken account into the exchange-correlation term. Thus, the general DFT energy expression would be: ][][][][][ xcnesDFT  EJETE  , (eq. 3.40) 48 Chap. 3. Methodology where Ts[ρ] corresponds to the kinetic energy functional described by the Slater molecular orbitals, and the rest expressions of the components of the eq. 3.40 are below:    occ N iiis T 1 2 2 1 ][  , (eq. 3.41)   nuclei || )()( ][ N aa aa ne d Z Er rR rR   , (eq. 3.42) ' |'| )'()( 2 1 ][ rr rr rr ddJ      . (eq. 3.43) Due to the fact that the electron interaction exists, the kinetic energy provided by Ts[ρ] is not the total kinetic energy (although HF theory yields the 99% approximately) and it is necessary to take it into account the difference, although small, between the exact kinetic energy of the real system and that calculated by assuming non-interacting electrons, the correlation kinetic energy. This difference is included in eq. 3.44 , in which it is also included the difference between Eee and J potential energy terms:     ][][][][][ eesxc  JETTE  . (eq. 3.44) The first parentheses corresponds to the kinetic correlation energy and in the second one coexist both exchange and potential correlation energy, although the exchange energy is by far the largest contributor to Exc. All the contributions to the energy without a simple expression as a function of the electron density are included within the Exc[ρ]. However, this exact Exc[ρ] functional is not known, and thus, as commented above, the aim of DFT methods is to design functionals connecting the electron density with the energy.47-52 Compared with HF theory, DFT provides much better results. Also, if the exact Exc[ρ] were known DFT would provide the exact total energy including the correlation. However, there is an important problem with density functional methods that is their inability to systematically improve the results and the known failure to describe certain important features, such as van der Waals interactions. Anyway, different functionals used for DFT methods have been developed assuming different approximations being then suitable for different kind of systems with acceptable accuracy. The method to construct the different functionals vary from quantum mechanics to methods in which parametric functions are used to best reproduce experimental results.2 Which functional is the best one will have to be settled by comparing the performance with experiments or high-level wave mechanics calculations. 3.3. Density Functional Theory methods 49 3.3.2. Functional types In early attempts on the DFT development all the energy components were expressed as a functional of the electron density; however, it was seen that this led to no good performance. Thus, instead of this, the modern methods use the electron density, represented by means of an auxiliary set of orbitals to calculate the electron kinetic energy (Kohn and Sham, 1965),54 and use approximations to calculate the exchange–correlation energy, which is the only unknown functional. There are many stages of development for the functionals. The so-called Local Density Approximation, LDA, is the simplest approximation, and it is based only on the electron density. In this approximation, the electron density can be treated as a uniform electron gas or, equivalently, as a slowly varying function, with the advantage that the calculations of the correlation energy of a uniform gas have been already determined by Monte Carlo methods for a number of different densities. For that reason, LDA approximation works very well in systems in which the density is maintained approximately constant. In 1980, Vosko, Wilk and Nusair (VWN)55 constructed a suitable analytic interpolation formula in order to use these results in DFT calculations. In the LDA approximation εxc[ρ] is a functional that depends exclusively on the density. Contributions to correlation and to exchange are usually treated separately as: ][][][  LDA c LDA x LDA xc  , (eq. 3.45) being the exchange energy per particle 3 1 3 1 8 3 4 9 ][           LDA x . (eq. 3.46) Thus, the correlation energy per particle can take different values. The simplest method is known as Xα, takes εc[ρ]=0 and α=2/3, and thereby this method includes electron exchange but not correlation. When open-shell systems are used, the method is called the local spin density approximation (LSDA). Both LDA and LSDA in general underestimate the exchange energy by ~10%, creating larger errors than the whole correlation energy.5 Also, electron correlation is overestimated, often by a factor close to 2, implying at the same time the overestimation of the bond strengths. However, LSDA methods often provide results with similar accuracy to that obtained by wave mechanics HF methods, and work very well in systems in which the density is maintained 50 Chap. 3. Methodology approximately constant.5, 52 The main problem is that they still fail to describe van der Waals complexes. Different improvements have been developed regarding to the accuracy of the method or the necessity to consider, in many cases, a non-uniform electron gas.1, 2 Besides, the LDA methods are also not suitable for systems with weak bonds or for making reliable thermochemical predictions, overestimating in general the bond energy by approximately 30%.20 The simplest improvement is obtained by means of the so-called gradient-corrected methods, in which the electron density and its gradient are used to construct the set of exchange-correlation functionals. These methods are also known as Gradient Corrected or Generalized Gradient Approximation (GGA) methods, being also sometimes referred to as non-local methods. To develop these methods, Becke proposed to include an additional term starting from the LDA functional, and this is the usual way to proceed. There are many other non-local corrections for both the exchange part and the correlation part, different from Becke’s correction, although the most often used for exchange is Becke’s. Among the corrections developed for the exchange energy there are for example PW86,56 developed by Perdew and Wang, or B88 proposed by Becke,57 which include a β parameter determined by fitting to known atomic data. Among the corrections aimed to improve the correlation energy, one popular functional is LYP,58, 59 proposed by Lee, Yang and Parr, where parameters are determined by fitting to data for helium. The widely used B-LYP functional comes from the combination of these exchange and correlation functionals.57-59 Another way to correct the functionals is by using hybrid methods. These methods combine functionals from other methods with part of HF calculations, usually the exchange integrals. This approximation can be developed when a uniform gas of non-interacting electrons is considered since there is no correlation energy but only exchange energy; indeed the exact exchange energy is given by HF method (the wavefunction for these systems corresponds to the single Slater determinant). Thus, the exchange and correlation energies could be expressed as the sum of the LSDA and a gradient correction term of exchange and correlation respectively, and including also the exact exchange for the exchange energy, as in the following expression:   GGA c LSDA c GGA x exact x LSDA x method hybrid xc c baa1 EE EEEΕ   . (eq. 3.47) The a, b and c parameters depend on the chosen forms for GGA x E and GGA c E , and are determined experimentally. Becke 3 parameter functional (B3) method is an example of such hybrid methods, in which the popular GGA exchange functional B88 is used as a correction to the 3.4. Interaction energy 57 molecule, taking as reference the isolated constituting atoms. However, it has been observed that intramolecular BSSE for small molecules gives negligible errors. Different behavior has the BSSE for large molecules in which there exists interaction among different areas or sections within the same molecule. In these situations, the best choice to reduce BSSE is to employ the largest basis set as possible.111, 112 3.4.2. Many-body effects When more than a pair of molecules is involved in a system, the calculation of the interaction energy becomes a harder process due to interactions among all the fragments.113 The total energy of the system could be expressed as the sum over all fragments, but also including a correction from all pairs, trios, etc., of interacting fragments; ... , ...      i ij jik ijk i ij ij iiijk EEEE . (eq. 3.61) On the other hand, the interaction energy can be expressed similarly to eq. 3.57, in which the interaction energy is the energy of the whole complex minus the energy corresponding to the constituting molecules. In this case, instead a dimer, it is considered the trimer ABC as example. Applying then the counterpoise method, in which also every fragment has the same geometry and basis set than the whole complex, the interaction energy can be expressed as: )()()()( ABCEABCEABCEABCEE trimer C trimer B trimer A trimer ABC int ABC  . (eq. 3.62) However, another more useful expression of the interaction energy of the cluster, in this case the trimer, can be used. In a first approach, since the monomer terms do not affect to the interaction energy (only to deformation), they are not considered. Instead of that, this interaction energy can be expressed as a pairwise sum of interaction between all pairs of molecules, plus a correction involving trios in this case: body int BC int AC int AB int ABC EEEEE   3 , (eq.3.63) where the E3-body corresponds to the correction due to the trio ABC, and the rest of the terms are the interaction energies of each pair of molecules, and can be expressed as: )()()( ABCEABCEABCEE trimer B trimer A trimer AB int AB  , (eq. 3.64) 58 Chap. 3. Methodology )()()( ABCEABCEABCEE trimer C trimer A trimer AC int AC  , (eq. 3.65) )()()( ABCEABCEABCEE trimer C trimer B trimer BC int BC  . (eq. 3.66) Similarly, the interaction energy for larger systems can be obtained, just increasing the number of correction terms (E4-body, E5-body, etc) within the eq. 3.63 until the order of the system.99, 104, 114, 115 In the case of a trimer, the three-body contribution to the interaction energy of the trimer would be: int BC int AC int AB int ABCbody EEEEE  3 . (eq. 3.67) This quantity is usually small in trimers, but there are systems in which that term becomes significant, as in those systems with important polarizing effects or hydrogen bonds. However, it becomes clear that when the number of constituting molecules of the system increases, also these many-body effects become more relevant. Finally, it is important to indicate that also BSSE seems to be related to the number of fragments forming the cluster, since an increase of the BSSE has been observed when more molecules were included, even when the value of BSSE for a given dimer was small. For that reason, it is mandatory to perform a proper treatment of BSSE in large clusters. 3.5. Symmetry-Adapted Perturbation Theory (SAPT) Although the magnitude of the interaction energy for a given geometry is easily obtained applying the supermolecule method with the wavefunction-based or DFT methods, detailed information about the nature, origins or characteristics of the interaction itself is not provided. Thus, methods which give this kind of information would be desirable.94, 116 Then, focusing on the intermolecular forces, it would be interesting to employ methods which provide us with values for the different contributions, i.e. electrostatics, repulsion, polarization and dispersion effects. Interaction energies, compared with those energies due to the chemical bond or with the intrinsic electronic energies of atoms or molecules, are more than an order of magnitude smaller in the best case. Thus, as it is a small contribution to the energy of the whole system, it can be treated as a perturbation, considering the most important contributions to the energy of the system are those from isolated monomers. Particularly, Symmetry-Adapted Perturbation Theory 3.5. Symmetry-Adapted Perturbation Theory (SAPT) 59 (SAPT) is a well-known method among those based on perturbation theory which, applying a perturbational scheme (Rayleigh-Schrödinger) to a complex, provides expressions for the different contributions: electrostatics, induction and dispersion effects. The simplest partitioning of the Hamiltonian for a pair of molecules A and B corresponds to the following expression: VHVHHH BA ˆˆˆˆˆˆ 0 , (eq. 3.68) where the operators A and B are the Hamiltonians of the isolated reference systems, those from the non-interacting molecules, 0 corresponding to the unperturbed Hamiltonian of the system and the operator is the perturbation, in our case, the interaction energy. The different contributions to the interaction energy are obtained as different first order, second order, … corrections of the energy as the result of applying the Rayleigh-Schrödinger perturbation theory under the previous assumptions. Thus, the first order correction energy corresponds to the electrostatic energy, since it shows the coulombic interaction between the electron density of both monomers, BABA el VE 0000 )1( ˆ   . (eq. 3.69) On the other hand, induction and dispersion contributions appear for the correction to second order, depending on single and double excitations. For the interacting monomers A and B, the induction contribution is computed as the single excitations within one monomer induced by the close presence of the other monomer, obtaining two equivalent expressions, one for each monomer. For example, the induction of single excitations within monomer A is expressed by the following expression:    00 2 000 ˆ mAA m BA m BA A ind EE V E  . (eq. 3.70) Finally, dispersion contribution is considered as the term within second order correction which takes into account double excitations:    0;0 00 2 00 ˆ nm BAB n A m B n A m BA disp EEEE V E  . (eq. 3.71) 60 Chap. 3. Methodology Thus, by means of these contributions together with the multipole expansion, intermolecular interactions can be described with a better qualitative and quantitative information, being possible describing them as functions of molecular properties as multipoles or polarizabilities.94, 116 However, this Rayleigh-Schrödinger perturbation theory is only successful for molecules a long distance apart, being completely inappropriate at short range.94 There are some reasons which lead to method failure and although one of them is due to the breaking down of the multipole expansion, the most fundamental failure comes from the fact that the repulsion between molecules at short range is not taken into account. This is an important consequence of ignoring the exchange contribution, and thus, the overlapping of wavefunctions when molecules are close enough is missed. The problem of this loss lies in the wrong symmetry of the reference wavefunction when the Rayleigh-Schrödinger perturbation theory is used to describe intermolecular interactions.94, 116 This reference wavefunction BA 000   is not antisymmetric upon exchange of electrons between A and B, so it cannot satisfy the Pauli principle. Thus, the description of repulsion forces when molecules are close is defective. Symmetrized Rayleigh-Schrödinger theory, or nowadays commonly called Symmetry Adapted Perturbation Theory (SAPT), is the most used method to solve the antisymmetry problem. Thus, the exchange-repulsion energy is achieved by forcing the antisymmetry in the energy expressions, modifying the electron density in the proper way which causes a repulsive force on the nuclei (exchange-repulsion energy).117-119 As in Rayleigh-Schrödinger perturbation theory at large distances, the SAPT procedure consists on a series of contributions which can be associated to physical effects when the order is low. However, the forced antisymmetrization produces new exchange-repulsion terms which now accompany to the old long-range terms. Thus, when the wavefunction of closed-shell molecules overlaps significantly, a strong repulsion appears overcoming the problem. In this way, the interaction energy by this procedure can be expressed as: ... (2) dispexch (2) disp (2) indexch (2) ind (1) exch (1) elint   EEEEEEE . (eq. 3.72) The monomers’ wavefunctions are not computed exactly, but separating HF and correlation contributions, including the intramolecular electron correlation. This intramolecular contribution is an important term that cannot be avoided since molecular properties required as SAPT inputs are significantly affected by electron correlation.18, 117-119 To solve this, the Hamiltonian can be expressed as a double perturbation in both the interaction and the intramonomer correlation, 3.5. Symmetry-Adapted Perturbation Theory (SAPT) 61 VWFWFH BBAA ˆˆˆˆˆˆ  . (eq. 3.73) The Hamiltonian is written as a sum of monomer Fock operators, F ˆ , the potential of each monomer, W ˆ , and the interaction potential, V ˆ . Therefore, now it is possible to write the interaction energy as:    0;1 int ji ij exch ij pol EEE , (eq. 3.74) where ij is the order of V ˆ and W ˆ . The accuracy of the interaction’s description will increase with higher orders in both expansions. Commonly, the expansion is truncated at second-order in V ˆ , since the thirdand higher order terms are usually negligible, but different truncations can be done giving rise to different models. Indeed, induction effects are often associated to higher terms, so that in polar systems the truncations would be done at higher orders. Another way to take into account such effects is by means of a correction term: )20()20()10()10( int indexchindexchel HF HF EEEEE    , (eq. 3.75) where HF Eint is the Hartree-Fock interaction energy calculated using the supermolecular method. Although at HF level dispersion is not included, by applying perturbation theory based on HF wavefunctions this term can be obtained by SAPT. Thus, the HF+D method can be obtained by correcting the HF calculation with the dispersion contribution, so the interaction energy using HF wavefunction up to second order can be expressed as the following sum of SAPT contributions. )20()20()20()20()10()10( int dispexchdispHFindexchindexchel DHF EEEEEEE    . (eq. 3.76) For a more accurate description of the interactions, intramonomer correlation effects could be included. For example, the so-called SAPT2 level, including similar corrections to MP2, corresponds to: )22()22()12()11()12( intint indexchindexchexchel DHFSAPT EEEEEEE   . (eq. 3.77) The next step would be including contributions to third order both in the intermolecular perturbation and the intramonomer correlation perturbation. However, this procedure requires relatively significant computer resources. Other approaches have been developed, as SAPT(DFT) 62 Chap. 3. Methodology that will be commented in the following section. Also, it is worth mentioning that Hohenstein and Sherrill have developed a density fitting based approach in order to reduce the cost of SAPT calculations.18, 120, 121 This approach has been included into the PSI4 program.122 3.5.1. SAPT(DFT) As commented above, to include intramonomer correlation effects within SAPT, together with the triple perturbation theory used in the SAPT method, would require high computational costs would be required, and when larger molecular systems are employed it becomes prohibitive. A solution to this difficulty could be a SAPT approach utilizing a density functional theory (DFT) description of monomers, since if a correlated description of monomers is used, the costly correlation terms could be avoided.123 Such a method, now called SAPT(DFT), has been first proposed by Williams and Chabalowski (2001) and later developed by others.123-125 This procedure consists on changing the HF orbitals and energies by their Kohn-Sham counterparts. Since the intramonomer correlation effects are already taken care of in DFT calculations, only the perturbation on the intermolecular interaction will suffice on the SAPT calculation.123 This improvement allows SAPT to be performed on much larger systems than previously did, although the initial results with this procedure were disappointing.94 The incorrect asymptotic behavior of the exchange-correlation functionals seems to be the main reason of these inaccurate results, with the consequence of one too diffuse electron density that leads to errors in both the first-order terms. In addition, the second-order terms were calculated using uncoupled perturbation theory and were also poor. To circumvent this problem, two very similar proposals were soon developed independently by Hesselmann and Jansen126, 127 and Szalewicz,128, 129 focused on the correction for the asymptotic behavior. In Jansen’s proposal the functional employed is combined with the LB94 functional, which has a correct asymptotic behavior at long-range, but fails at short-range,127 arising a new problem consisting on inconsistency at intermediate distances between the behavior of the two functionals.130, 131 In order to solve this, a gradient-regulated connection method developed by Grüning et al. was employed.130 This is called Adiabatic Correction AC and requires the sum IP + εHOMO as the input parameter, which vanishes for 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. 3.6. Electron density analysis 63 After solving these accuracy problems, the SAPT(DFT) interaction energy can be expressed as: HFdispexchdisprespindexch respindexchel DFTSAPT EEE EEEE      )20()20()20( , )20( , )10()10( int , (eq. 3.78) where a δHF correction term is also included which takes into account higher-order contributions. With this implementation of DFT within SAPT theory, the computational costs decrease considerably for the description of intramonomer correlation, allowing the use of this method to larger systems than those allowed for SAPT based in MBPT or CC. To reduce even more the computational costs, the density fitting approach can be employed. Furthermore, the performance of these DFT-based methods seems to be better than the more conventional SAPT using triple perturbation theory. Also, simplified extrapolation schemes can be used to reach the complete basis set limit in the framework of SAPT(DFT) calculations.132 3.6. Electron density analysis Molecular properties often can be expressed in terms of a property density.133 In these cases, it is possible to obtain the contribution of an atom to that molecular property by integrating this density over the bounded volume of that atom in the molecule. Thus, a set of physical properties characterize each atom in a molecule or crystal, corresponding then to molecular properties, and being added up to those properties of the whole molecular system. In the next sections methods based on this principle will be exposed. 3.6.1. QTAIM One of the aims that led the development of the Atoms in Molecules Theory (AIM) was that most methods are only focused on obtaining approximate solutions to Schrödinger’s equation. As Bader indicated, these approximate solutions, although provide a very useful tool to predict and understand the properties of some systems of interest (energy, geometry and other many properties), are not couched in the conceptual language of chemistry, as for instance the Lewis model of the electron pair or the molecular structure hypothesis.134 Both of them treat the systems from chemical observations, and quantum mechanics lacked the incorporation of these observations into it. On the one hand, since Lewis model of the electron pair, it is known that the localized or delocalized nature of the electrons is related with aspects as resonance, reactivity, geometry or bonding. On the other hand, the molecular structure hypothesis relates the molecular structure of a system with its properties. 64 Chap. 3. Methodology The molecular structure hypothesis comes from Dalton’s atomic hypothesis formulated in 1807, in which the atoms retain their identity even in chemical combination with other atoms. This hypothesis evolved incorporating the geometrical aspects, providing the concept of molecular structure, according to which the structure of a molecule is imparted by the bonds, which form a network linking the different atoms in a molecule. This hypothesis came from the observation that atoms or groups of atoms have characteristic sets of properties which in general encompass a narrow range, so that those atoms or groups of atoms can be identified or a given behavior can be predicted. It was found that, if the distribution of charge of a given atom or group of atoms remained constant, their properties remained also constant, even their contribution to the total energy of the system. Also, if the distribution of charge changes the properties change by the same amount. According to Bader, this is because there exists a direct relationship between the properties of an atom and its spatial form even when the atom belongs to different systems, although obviously more evident when the atom or the group of atoms are transferrable. Thus, the properties of atoms in molecules can be experimentally measured, and with more precision those which are transferrable. With that premise arose the idea that the development of an AIM method was possible, with the purpose to relate molecular properties to those of its constituent atoms. The development of this method was made possible by the knowledge that the electron density distribution, and its topology in particular, plays an essential role in experimental chemistry, providing a good explanation and understanding of the observations like chemical structure, functional groups, transferability, chemical reactivity or chemical bonding.134 Also, the real three dimensional space was at last introduced into theoretical chemistry by this theory, since electron density provides a description of charge throughout real space, and thus relates the concepts of state functions in Hilbert space and the physical model of matter in space, according to Bader’s words. Bader and coworkers started to develop the foundations of AIM theory in the seventies, on the basis of their own works on molecular electron density distributions made a decade before. This theory is also-called quantum theory of atoms in molecules (QTAIM) in more recent literature due to its rigorous basis in quantum mechanics. 3.6.1.1. The topology of the electron density The nucleus of each atom plays a role as a charge attractor, being the principal force involved in the final topology of the electron density, and therefore, providing it of a maximum of electron density at the position of each nucleus. In this way, the final density distribution in a molecule is described basically by the balance of the forces that the neighboring nuclei exert on the electrons. To know how the topology of the electron density, ρ(r), is, it should be necessary to indicate how 3.6. Electron density analysis 65 it changes, and for that purpose the gradient of the electron density, noted as ∇ρ(r), is employed. As a definition of a gradient of any scalar at a point in space, ∇ρ(r) denotes 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. In this way, the gradient is going to be used to define the “critical points” (CP) whose study will determine the topology of the electron density of each system. Thus, a “critical point” in the electron density is defined as the point in space where the gradient of the electron density (the first derivatives respect to the coordinates) vanish. The results about the topology of the electron density are provided by the study of the critical points associated to the electron density of the system. At the same time, critical points are classified according to two properties called rank (ω) and signature (σ), and denoting each critical point as the couple of their values by (ω, σ). The rank is defined as the number of non-zero curvatures of ρ at the critical point. The signature is defined as the algebraic sum of the signs of the curvatures, contributing each curvature with ±1 depending on whether it is a positive or negative curvature. The most common value of ω is 3, since values of ω < 3 are mathematically unstable and tend to vanish or bifurcate, and are not usually found in equilibrium charge distributions. In this way, there are four types of stable critical points with ω = 3, corresponding with different elements of chemical structure: a) Nuclear critical point (NCP): (3, –3) Three negative curvatures: ρ is a local maximum. b) Bond critical point (BCP): (3, –1) Two negative curvatures: ρ is a maximum in the plane defined by the two associated axes, and a minimum along the axis perpendicular to this plane. c) Ring critical point (RCP): (3, +1) Two positive curvatures: ρ is a minimum in the plane defined by the two associated axes, and a maximum along the axis perpendicular to this plane. d) Cage critical point (CCP): (3, +3) Three positive curvatures: ρ is a local minimum. In a molecule or crystal there exists a relationship between the number and the type of critical points: , where n denotes the number of critical points of each type, and the set { } for a given system is called as the “characteristic set”. A useful concept is the bond path (BP), which is described as the single line of maximum electron density linking the nuclei of two chemically bonded atoms. This line can be curved to a 66 Chap. 3. Methodology greater or lesser degree depending on the ring strain of the molecule. Also, if some bond paths arrange in a plane forming a ring, a ring critical point will appear among them in the same plane. Raising the dimensionality, if some of those rings arrange forming a polyhedron, a cage critical point will appear among them into the volume enclosed by the rings. The bond path gives hints on the type of chemical bonding, i.e. weak, strong, closed-shell, and open shell interactions. Other concept commonly used is the bond critical point (BCP), which describes the point into the bond path in which the electron density of the two atoms linked has the same value. In other words, the bond critical point is the point where ∇ ρ(r) = 0, along the bond path, and also the point with the lowest value of electron density. Thus, there exists a zero-flux surface defined by all of ∇ ρ(r) trajectories verifying that ∇ ρ(r) = 0 and which intersects the bond path in the BCP. That occurs between each pair of linked atoms. With only those two concepts, bond path and bond critical point, the molecular graph can be constructed with the aim to visualize the molecular structure. This molecular graph constitutes the representation of the bond paths linking the nuclei of bonded atoms in equilibrium geometry with the associated critical points. 3.6.2. NCI index The non covalent interaction index (NCI) was developed with the aim of having a method to view and analyze non covalent interactions.135, 136 This index provides the most relevant features of the interaction in a quick and simple way, without unnecessary extra data. The background of the method relies on the fact that the so-called reduced density gradient, s, takes different values depending on the kind of interaction. This variable comes from the electron density, ρ, and its first derivative, and has the following expression:  3431 2 32 1     s . (eq. 3.79) In regions far from the molecules with very small values of density, the reduced density gradient will have very large positive values. Conversely, within regions in which there are interactions (covalent or non covalent), the reduced density gradient will assume very small values, nearly zero. To identify non covalent interactions among all small reduced gradients, s is plotted versus ρ, revealing that these interactions exhibit one or more spikes in the low-density, low-gradient region, a signature of non covalent interactions. Apparently, this feature comes from the 3.8. References 73 3.8. References (1) C. J. Cramer, Essentials of Computational Chemistry: Theories and Models, 2nd Edition, John Wiley & Sons, Chichester, 2004. (2) D. C. Young, Computational Chemistry: A Practical Guide for Applying Techniques to Real World Problems, John Wiley & Sons, Inc., New York, 2001. (3) S. Attila, S. O. Neil, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory, Dover Publications, New York, 2012. (4) R. J. Bartlett. Annu. Rev. Phys. Chem. 1981, 32, 359-401. (5) F. Jensen, Introduction to computational chemistry, John Wiley and Sons, Chichester, 2001. (6) S. Grimme. J. Chem. Phys. 2003, 118, 9095-9102. (7) S. Grimme, L. Goerigk, R. F. Fink. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2012, 2, 886906. (8) J. G. Hill, J. A. Platts. J. Chem. Theory Comput. 2006, 3, 80-85. (9) M. Gerenkamp, S. Grimme. Chem. Phys. Lett. 2004, 392, 229-235. (10) T. P. M. Goumans, A. W. Ehlers, K. Lammertsma, E.-U. Wuerthwein, S. Grimme. Chem. - Eur. J. 2004, 10, 6468-6475. (11) S. Grimme. J. Phys. Chem. A 2005, 109, 3067-3077. (12) J. Antony, S. Grimme. J. Phys. Chem. A 2007, 111, 4862-4868. (13) K. E. Riley, J. A. Platts, J. Řezáč, P. Hobza, J. G. Hill. J. Phys. Chem. A J Phys Chem A 2012, 116, 4159-4169. (14) S. Grimme. J. Comput. Chem. 2003, 24, 1529-1537. (15) I. Hyla-Kryspin, S. Grimme. Organometallics 2004, 23, 5581-5592. (16) M. Pitonak, J. Rezac, P. Hobza. Phys. Chem. Chem. Phys. 2010, 12, 9611-9614. (17) T. Takatani, E. G. Hohenstein, C. D. Sherrill. J. Chem. Phys. 2008, 128, 124111/1-124111/7. (18) E. G. Hohenstein, C. D. Sherrill. WIREs Comput. Mol. Sci. 2012, 2, 304-326. (19) A. Hellweg, S. A. Grun, C. Hattig. Phys. Chem. Chem. Phys. 2008, 10, 4119-4127. (20) J. A. Bort, J. B. Rusca, Theoretical and Computational Chemistry: Foundations, Methods and Techniques, Publicacions de la Universitat Jaume I, Castelló de la Plana, 2007. (21) T. Van Voorhis, M. Head-Gordon. J. Chem. Phys. 2000, 113, 8873-8879. (22) K. Raghavachari, G. W. Trucks, J. A. Pople, M. Head-Gordon. Chem. Phys. Lett. 1989, 157, 479-483. (23) T. Helgaker, W. Klopper, A. Halkier, K. Bak, P. Jørgensen, J. Olsen, Highly Accurate Ab Initio Computation of Thermochemical Data, in Quantum-Mechanical Prediction of Thermochemical Data, Vol. 22 (Ed. J. Cioslowski), Springer Netherlands, New York, 2001, pp.1-30. (24) C. Hättig, W. Klopper, A. Köhn, D. P. Tew. Chem. Rev. 2011, 112, 4-74. (25) L. Kong, F. A. Bischoff, E. F. Valeev. Chem. Rev. 2011, 112, 75-107. (26) T. Helgaker, W. Klopper, H. Koch, J. Noga. J. Chem. Phys. 1997, 106, 9639-9646. 74 Chap. 3. Methodology (27) A. Halkier, T. Helgaker, P. Jørgensen, W. Klopper, H. Koch, J. Olsen, A. K. Wilson. Chem. Phys. Lett. 1998, 286, 243-252. (28) D. G. Truhlar. Chem. Phys. Lett. 1998, 294, 45-48. (29) A. J. C. Varandas. J. Chem. Phys. 2007, 126, 244105/1-244105/15. (30) A. Halkier, T. Helgaker, P. Jørgensen, W. Klopper, J. Olsen. Chem. Phys. Lett. 1999, 302, 437-446. (31) D. G. Liakos, F. Neese. J. Phys. Chem. A J Phys Chem A 2012, 116, 4801-4816. (32) F. Neese, E. F. Valeev. J. Chem. Theory Comput. 2010, 7, 33-43. (33) R. A. Kendall, T. H. Dunning, R. J. Harrison. J. Chem. Phys. 1992, 96, 6796-6806. (34) T. Helgaker, P. Jørgensen, J. Olsen, Molecular electronic-structure theory John Wiley & Sons, Chichester, 2000. (35) W. Klopper, W. Kutzelnigg. J. Mol. Struct.: THEOCHEM 1986, 135, 339-356. (36) A. Karton, J. L. Martin. Theor. Chem. Acc. 2006, 115, 330-333. (37) F. Jensen. Theor. Chem. Acc. 2005, 113, 267-273. (38) K. E. Riley, P. Hobza. WIREs Comput. Mol. Sci. 2011, 1, 3-17. (39) P. Jurecka, J. Sponer, J. Cerny, P. Hobza. Phys. Chem. Chem. Phys. 2006, 8, 1985-1993. (40) J. Řezáč, K. E. Riley, P. Hobza. J. Chem. Theory Comput. 2014, 10, 1359-1360. (41) J. Řezáč, K. E. Riley, P. Hobza. J. Chem. Theory Comput. 2011, 7, 2427-2438. (42) P. Hobza, J. Šponer. J. Am. Chem. Soc. 2002, 124, 11802-11808. (43) M. Pitonak, P. Neogrady, J. Cerny, S. Grimme, P. Hobza. ChemPhysChem 2009, 10, 282289. (44) K. E. Riley, J. Rezac, P. Hobza. Phys. Chem. Chem. Phys. 2011, 13, 21121-21125. (45) R. Sedlak, K. E. Riley, J. Řezáč, M. Pito ňák, P. Hobza. ChemPhysChem 2013, 14, 698-707. (46) K. E. Riley, J. Rezac, P. Hobza. Phys. Chem. Chem. Phys. 2012, 14, 13187-13193. (47) R. G. Parr, W. Yang, Density-Functional Theory of Atoms and Molecules, Oxford University Press, USA, Oxford, 1989. (48) L. J. Bartolotti, K. Flurchick, An Introduction to Density Functional Theory, in Rev. Comput. Chem., Vol., John Wiley & Sons, Inc., Hoboken, 2007, pp.187-216. (49) A. St-Amant, Density Functional Methods in Biomolecular Modeling, in Rev. Comput. Chem., Vol., John Wiley & Sons, Inc., Hoboken, 2007, pp.217-259. (50) T. Ziegler. Chem. Rev. 1991, 91, 651-667. (51) E. J. Baerends, O. V. Gritsenko. J. Phys. Chem. A J Phys Chem A 1997, 101, 5383-5403. (52) W. Koch, M. C. Holthausen, A Chemist's guide to density functional theory, Wiley-VCH, Weinhim, 2000. (53) P. Hohenberg, W. Kohn. Phys. 1964, 136, B864-B871. (54) W. Kohn, L. J. Sham. Phys. 1965, 140, A1133-A1138. (55) S. H. Vosko, L. Wilk, M. Nusair. Can. J. Phys. 1980, 58, 1200-1211. (56) J. P. Perdew. Phys. Rev. B 1986, 33, 8822-8824. (57) A. D. Becke. Phys. Rev. A 1988, 38, 3098-3100. 3.8. References 75 (58) C. Lee, W. Yang, R. G. Parr. Phys. Rev. B 1988, 37, 785-789. (59) B. Miehlich, A. Savin, H. Stoll, H. Preuss. Chem. Phys. Lett. 1989, 157, 200-206. (60) A. D. Becke. J. Chem. Phys. 1993, 98, 1372-1377. (61) Y. Zhao, D. G. Truhlar. J. Chem. Phys. 2006, 125, 194101/1-194101/18. (62) Y. Zhao, D. Truhlar. Theor. Chem. Acc. 2008, 120, 215-241. (63) Y. Zhao, D. G. Truhlar. J. Chem. Theory Comput. 2008, 4, 1849-1868. (64) R. Peverati, D. G. Truhlar. The Journal of Physical Chemistry Letters 2011, 3, 117-124. (65) R. Peverati, D. G. Truhlar. The Journal of Physical Chemistry Letters 2011, 2, 2810-2817. (66) R. Peverati, D. G. Truhlar. Phys. Chem. Chem. Phys. 2012, 14, 16187-16191. (67) R. Peverati, D. G. Truhlar. Phys. Chem. Chem. Phys. 2012, 14, 13171-13174. (68) E. G. Hohenstein, S. T. Chill, C. D. Sherrill. J. Chem. Theory Comput. 2008, 4, 1996-2000. (69) J. Cerny, P. Hobza. Phys. Chem. Chem. Phys. 2005, 7, 1624-1626. (70) S. Tsuzuki, H. P. Lüthi. J. Chem. Phys. 2001, 114, 3949-3957. (71) M. J. Allen, D. J. Tozer. J. Chem. Phys. 2002, 117, 11113-11120. (72) P. Hobza, J. šponer, T. Reschel. J. Comput. Chem. 1995, 16, 1315-1325. (73) S. Kristyán, P. Pulay. Chem. Phys. Lett. 1994, 229, 175-180. (74) N. Kurita, H. Sekino. Int. J. Quantum Chem 2003, 91, 355-362. (75) E. R. Johnson, R. A. Wolkow, G. A. DiLabio. Chem. Phys. Lett. 2004, 394, 334-338. (76) Y. Zhao, N. E. Schultz, D. G. Truhlar. J. Chem. Phys. 2005, 123, 161103/1-161103/4. (77) Y. Zhao, N. E. Schultz, D. G. Truhlar. J. Chem. Theory Comput. 2006, 2, 364-382. (78) Y. Zhao, B. J. Lynch, D. G. Truhlar. J. Phys. Chem. A J Phys Chem A 2004, 108, 4786-4791. (79) S. Grimme. J. Chem. Phys. 2006, 124, 034108/1-034108/16. (80) S. Grimme. J. Comput. Chem. 2004, 25, 1463–1473. (81) U. Zimmerli, M. Parrinello, P. Koumoutsakos. J. Chem. Phys. 2004, 120, 2693-2699. (82) S. Grimme. J. Comput. Chem. 2006, 27, 1787–1799. (83) S. Grimme, J. Antony, S. Ehrlich, H. Krieg. J. Chem. Phys. 2010, 132, 154104/1154104/19. (84) S. Grimme, S. Ehrlich, L. Goerigk. J. Comput. Chem. 2011, 32, 1456-1465. (85) S. Grimme. WIREs Comput. Mol. Sci. 2011, 1, 211-228. (86) A. Krishtal, C. Van Alsenoy, P. Geerlings. J. Chem. Phys. 2014, 140, 184105/1-184105/14. (87) S. N. Steinmann, C. Corminboeuf. J. Chem. Theory Comput. 2010, 6, 1990-2001. (88) S. N. Steinmann, C. Corminboeuf. J. Chem. Phys. 2011, 134, 044117/1-044117/5. (89) L. A. Burns, Á. V. Mayagoitia, B. G. Sumpter, C. D. Sherrill. J. Chem. Phys. 2011, 134, 084107/1-084107/25. (90) Q. Wu, W. Yang. J. Chem. Phys. 2002, 116, 515-524. (91) O. A. Vydrov, T. Van Voorhis. J. Chem. Phys. 2010, 133, 244103/1-244103/9. (92) M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, B. I. Lundqvist. Phys. Rev. Lett. 2004, 92, 246401. (93) F. D. John, G. Tim. J. Phys.: Condens. Matter 2012, 24, 073201. 76 Chap. 3. Methodology (94) A. J. Stone, The theory of intermolecular forces, Oxford University Press, Oxford, 2013. (95) G. Chałasiński, M. M. Szczȩśniak. Chem. Rev. 2000, 100, 4227-4252. (96) K. Szalewicz, B. Jeziorski. J. Chem. Phys. 1998, 109, 1198-1200. (97) S. S. Xantheas. J. Chem. Phys. 1996, 104, 8821-8824. (98) E. Cabaleiro-Lago, J. Rodríguez-Otero, Á. Peña-Gallego. Theor. Chem. Acc. 2011, 128, 531539. (99) E. M. Cabaleiro-Lago, M. A. R os. J. Chem. Phys. 2000, 112, 2155-2163. (100) M. v. Hopffgarten, G. Frenking. WIREs Comput. Mol. Sci. 2012, 2, 43-62. (101) A. Campo-Cacharrón, E. M. Cabaleiro-Lago, J. Rodríguez-Otero. ChemPhysChem 2012, 13, 570-577. (102) E. M. Cabaleiro-Lago, Á. Peña-Gallego, J. Rodríguez-Otero. J. Chem. Phys. 2008, 128, 194311/1-194311/8. (103) G. Chalasinski, M. M. Szczesniak. Chem. Rev. 1994, 94, 1723-1765. (104) C. D. Sherrill, Computations of Non covalent π Interactions, in Rev. Comput. Chem., Vol., John Wiley & Sons, Inc., Hoboken, 2009, pp.1-38. (105) F. B. van Duijneveldt, J. G. C. M. van Duijneveldt-van de Rijdt, J. H. van Lenthe. Chem. Rev. 1994, 94, 1873-1885. (106) H. B. Jansen, P. Ros. Chem. Phys. Lett. 1969, 3, 140-143. (107) S. F. Boys, F. Bernardi. Mol. Phys. 1970, 19, 553-566. (108) J. Alvarez-Idaboy, A. Galano. Theor. Chem. Acc. 2010, 126, 75-85. (109) Ł. M. Mentel, E. J. Baerends. J. Chem. Theory Comput. 2013, 10, 252-267. (110) T. Helgaker, W. Klopper, D. P. Tew. Mol. Phys. 2008, 106, 2107-2143. (111) R. M. Balabin. J. Chem. Phys. 2010, 132, 211103/1-211103/4. (112) D. Asturiol, M. Duran, P. Salvador. J. Chem. Theory Comput. 2009, 5, 2574-2581. (113) M. J. Elrod, R. J. Saykally. Chem. Rev. 1994, 94, 1975-1997. (114) K. Walczak, J. Friedrich, M. Dolg. J. Chem. Phys. 2011, 135, 134118/1-134118/11. (115) J. F. Ouyang, M. W. Cvitkovic, R. P. A. Bettens. J. Chem. Theory Comput. 2014. (116) I. Kaplan, Intermolecular Interactions : Physical Picture , Computational Methods, John Wiley & Sons, Chichester, 2006. (117) B. Jeziorski, R. Moszynski, K. Szalewicz. Chem. Rev. 1994, 94, 1887-1930. (118) R. Moszynski. Mol. Phys. 1996, 88, 741-758. (119) R. Moszynski, B. Jeziorski, A. Ratkiewicz, S. a. Rybak. J. Chem. Phys. 1993, 99, 8856-8869. (120) E. G. Hohenstein, C. D. Sherrill. J. Chem. Phys. 2010, 133, 014101/1-014101/12. (121) E. G. Hohenstein, C. D. Sherrill. J. Chem. Phys. 2010, 132, 184111/1-184111/10. (122) 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. (123) H. L. Williams, C. F. Chabalowski. J. Phys. Chem. A J Phys Chem A 2001, 105, 646-659. 3.8. References 77 (124) A. J. Misquitta, K. Szalewicz. Chem. Phys. Lett. 2002, 357, 301-306. (125) A. J. Misquitta, B. Jeziorski, K. Szalewicz. Phys. Rev. Lett. 2003, 91, 033201. (126) A. Heßelmann, G. Jansen. Chem. Phys. Lett. 2002, 357, 464-470. (127) G. Jansen. WIREs Comput. Mol. Sci. 2014, 4, 127-144. (128) A. J. Misquitta, R. Podeszwa, B. Jeziorski, K. Szalewicz. J. Chem. Phys. 2005, 123, 214103/1-214103/14. (129) K. Szalewicz. WIREs Comput. Mol. Sci. 2012, 2, 254-272. (130) M. Grüning, O. V. Gritsenko, S. J. A. van Gisbergen, E. J. Baerends. J. Chem. Phys. 2001, 114, 652-660. (131) D. J. Tozer, N. C. Handy. Mol. Phys. 2003, 101, 2669-2675. (132) J. Řezáč, P. Hobza. J. Chem. Theory Comput. 2011, 7, 685-689. (133) Cherif F. Matta (Editor), Russell J. Boyd (Editor), The Quantum Theory of Atoms in Molecules: From Solid State to DNA and Drug Design, WILEY-VCH, Weinheim, 2007. (134) R. F. W. Bader, Atoms in Molecules: A Quantum Theory, Clarendon Press, Oxford 1990. (135) 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. (136) 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. (137) R. F. W. Bader, H. Essén. J. Chem. Phys. 1984, 80, 1943-1960. (138) R. F. W. Bader. J. Phys. Chem. A J Phys Chem A 1998, 102, 7314-7323. (139) M. Orozco, F. J. Luque. Chem. Rev. 2000, 100, 4187-4226. (140) Q. Liu, J. W. Brady. J. Phys. Chem. B 1997, 101, 1317-1321. (141) M. A. Vincent, I. J. Palmer, I. H. Hillier. J. Mol. Struct.: THEOCHEM 1997, 394, 1-9. (142) P. Manivet, M. Masella. Chem. Phys. Lett. 1998, 288, 642-646. (143) S. Koneshan, J. C. Rasaiah, R. M. Lynden-Bell, S. H. Lee. J. Phys. Chem. B 1998, 102, 4193-4204. (144) H.-P. Cheng. J. Phys. Chem. A J Phys Chem A 1998, 102, 6201-6204. (145) M. Diraison, P. Millie, S. Pommeret, T. Gustavsson, J. C. Mialocq. Chem. Phys. Lett. 1998, 282, 152-158. (146) C. Alemán, S. E. Galembeck. Chem. Phys. 1998, 232, 151-159. (147) R. Biswas, S. Bhattacharyya, B. Bagchi. J. Phys. Chem. B 1998, 102, 3252-3256. (148) R. Gratias, H. Kessler. J. Phys. Chem. B 1998, 102, 2027-2031. (149) N. U. Zhanpeisov, J. Leszczynski. J. Phys. Chem. A J Phys Chem A 1999, 103, 8317-8327. (150) N. U. Zhanpeisov, J. Leszczynski. J. Phys. Chem. A J Phys Chem A 1998, 102, 6167-6172. (151) M. L. S. Mendoza, M. A. Aguilar, F. J. O. del Valle. J. Mol. Struct.: THEOCHEM 1998, 426, 181-190. (152) K. A. Swiss, R. A. Firestone. J. Phys. Chem. A J Phys Chem A 1999, 103, 5369-5372. (153) R. Takasu, K. Hashimoto, R. Okuda, K. Fuke. J. Phys. Chem. A J Phys Chem A 1999, 103, 349-354. 78 Chap. 3. Methodology (154) M. Masamura. J. Mol. Struct.: THEOCHEM 1999, 466, 85-93. (155) V. Tran, B. J. Schwartz. J. Phys. Chem. B 1999, 103, 5570-5580. (156) K. Coutinho, N. Saavedra, S. Canuto. J. Mol. Struct.: THEOCHEM 1999, 466, 69-75. (157) T. van Mourik, S. L. Price, D. C. Clary. J. Phys. Chem. A J Phys Chem A 1999, 103, 16111618. (158) T.-M. Chang, L. X. Dang. J. Phys. Chem. B 1997, 101, 10518-10526. (159) J. M. Martínez, R. R. Pappalardo, E. S. Marcos. J. Phys. Chem. A J Phys Chem A 1997, 101, 4444-4448. (160) G. D. Scholes, R. D. Harcourt, I. R. Gould, D. Phillips. J. Phys. Chem. A J Phys Chem A 1997, 101, 678-684. (161) A. Tongraar, K. R. Liedl, B. M. Rode. Chem. Phys. Lett. 1998, 286, 56-64. (162) H. Sato, F. Hirata, A. B. Myers. J. Phys. Chem. A J Phys Chem A 1998, 102, 2065-2071. (163) K. Hashimoto, T. Kamimoto. J. Am. Chem. Soc. 1998, 120, 3560-3570. (164) J. M. Martínez, R. R. Pappalardo, E. S. Marcos, K. Refson, S. Díaz-Moreno, A. Muñoz-Páez. J. Phys. Chem. B 1998, 102, 3272-3282. (165) T.-M. Chang, L. X. Dang. J. Phys. Chem. B 1999, 103, 4714-4720. (166) A. Tongraar, B. M. Rode. J. Phys. Chem. A J Phys Chem A 1999, 103, 8524-8527. (167) H. D. Pranowo, B. M. Rode. J. Phys. Chem. A J Phys Chem A 1999, 103, 4298-4302. (168) S. W. Rick, B. J. Berne. J. Phys. Chem. B 1997, 101, 10488-10493. (169) P. E. Smith. J. Phys. Chem. B 1999, 103, 525-534. (170) B. Madan, K. Sharp. Biophys. Chem. 1999, 78, 33-41. (171) Y.-K. Cheng, P. J. Rossky. Biopolymers 1999, 50, 742-750. (172) Y.-P. Pang, J. L. Miller, P. A. Kollman. J. Am. Chem. Soc. 1999, 121, 1717-1725. (173) D. Jacquemin, C. Michaux, E. A. Perpète, G. Frison. J. Phys. Chem. B 2011, 115, 3604-3613. (174) J. P. Gallivan, D. A. Dougherty. J. Am. Chem. Soc. 2000, 122, 870-874. (175) P. E. Mason, C. E. Dempsey, G. W. Neilson, S. R. Kline, J. W. Brady. J. Am. Chem. Soc. 2009, 131, 16689-16696. (176) C. D. Tatko, M. L. Waters. Protein Sci. 2003, 12, 2443-2452. (177) A. Rodríguez-Sanz, J. Carrazana-García, E. Cabaleiro-Lago, J. Rodríguez-Otero. J. Mol. Model. 2013, 19, 1985-1994. (178) C. Adamo, G. Berthier, R. Savinelli. Theor. Chem. Acc. 2004, 111, 176-181. (179) Y. Xu, J. Shen, W. Zhu, X. Luo, K. Chen, H. Jiang. J. Phys. Chem. B 2005, 109, 5945-5949. (180) A. S. Reddy, H. Zipse, G. N. Sastry. J. Phys. Chem. B 2007, 111, 11546-11553. (181) A. Campo-Cacharrón, E. Cabaleiro-Lago, J. Rodríguez-Otero. Theor. Chem. Acc. 2012, 131, 1-13. (182) A. S. Mahadevi, G. N. Sastry. Chem. Rev. 2013, 113, 2100-2138. (183) E. M. Cabaleiro-Lago, J. Rodríguez-Otero, Á. Peña-Gallego. J. Chem. Phys. 2011, 135, 214301/1-214301/9. 3.8. References 79 (184) C. Michaux, J. Wouters, D. Jacquemin, E. A. Perpète. Chem. Phys. Lett. 2007, 445, 57-61. (185) C. Michaux, J. Wouters, E. A. Perpète, D. Jacquemin. J. Phys. Chem. B 2008, 112, 24302438. (186) C. Michaux, J. Wouters, E. A. Perpète, D. Jacquemin. J. Am. Soc. Mass. Spectrom. 2009, 20, 632-638. (187) C. Michaux, J. Wouters, E. A. Perpète, D. Jacquemin. J. Phys. Chem. B 2008, 112, 98969902. (188) C. Michaux, J. Wouters, E. A. Perpète, D. Jacquemin. J. Phys. Chem. B 2008, 112, 77027705. (189) J. Tomasi, M. Persico. Chem. Rev. 1994, 94, 2027-2094. (190) J. Tomasi, B. Mennucci, R. Cammi. Chem. Rev. 2005, 105, 2999-3094. (191) C. J. F. Böttcher, O. C. van Belle, P. Bordewijk, A. Rip, Theory of electric polarization, Elsevier Scientific Pub. Co., Amsterdam, 1978. (192) A. Klamt, G. Schuurmann. J. Chem. Soc. Perkin Trans. 2 1993, 799-805. 4 Microhydration study: complexes between NH4+ and CH3NH3+ with phenol 4.3. Results 89 However, the second most stable structure PA1-2 is only 0.5 kcal mol-1 less stable. After these two structures, there is an energy gap of about 2 kcal mol-1 to the next stable minimum. It is worth noting that inclusion of the first water molecule makes the most stable minimum to be that with the cation over the phenyl ring, contrary to the behavior observed in the cluster without water molecules. This is a consequence of the secondary interaction established by the water molecule and the phenol moiety. The water molecule can establish a hydrogen bond with the hydroxyl oxygen, which is stronger than the hydrogen bond formed with the aromatic ring (typical values for water···water hydrogen bonds are around 4-5 kcal mol-1 whereas for O-H··· amount to 2-3 kcal mol-1).50 Therefore, even though the interaction between the cation and the phenol molecule is weaker over the ring, this loss is compensated with the additional strength of the O-H···O hydrogen bond formed. This effect can be clearly seen in structures PA1-3 and PA1-4; since there are no hydrogen bonds between water and phenol, the most stable structure is PA1-3, with the cation over the hydroxyl group. Also, the complexation energy difference between PA1-1 and PA1-4 allows an estimation of the contribution of the hydrogen bond to phenol of about 3.5 kcal mol-1 (1.9 kcal mol-1 for the water contact with the aromatic ring). Methylamonium complexes behave in a similar manner, showing complexation energies around 2.5-3 kcal mol-1 less negative than the corresponding ammonium minima. As before, no changes are observed after introduction of ZPE or enthalpy corrections. Table 4.2. Complexation energies (kcal mol-1) obtained for the most stable minima of the clusters containing one water molecule (see Figure 4.2) as obtained at the MP2/6-31+G(2d,p) level. Ammonium Methylammonium ∆Ecomplex ∆EZPE ∆H298 ∆Ehyd ∆Ecomplex ∆EZPE ∆H298 ∆Ehyd PA1-1 -37.77 -34.37 -35.17 -18.56 PM1-1 -35.26 -29.00 -29.35 -18.03 PA1-2 -37.30 -34.30 -34.91 -18.03 PM1-2 -34.17 -28.20 -28.38 -16.95 PA1-3 -35.44 -32.81 -33.02 -16.07 PM1-3 -32.70 -27.12 -26.88 -15.36 PA1-4 -34.25 -31.70 -31.76 -14.94 PM1-4 -32.38 -26.94 -26.59 -15.04 The energy differences between analogous structures of ammonium and methylamonium is larger than in complexes without water, suggesting that ammonium cation is able to polarize more efficiently the water molecule, giving an extra stabilization to the complexes compared to methylammonium. Another way of quantifying the effect of the new water molecule is focusing on the energy change observed when a water molecule is added to the complexes without water. 90 Chap. 4. Microhydration study: complexes between NH4+ and CH3NH3+ with phenol Therefore, the formation of PA1-1 from PA0-3 implies an energy gain of -19.8 kcal mol-1, whereas in forming PA1-2 from PA0-1 only changes by -17.4 kcal mol-1. These 2.5 kcal mol-1 reflect the commented above about the different strength of O-H··· and O-H···O contacts. On the other hand, forming PA1-4 from PA0-3 changes the energy by -15.5 kcal mol-1, so the formation of the water···phenol hydrogen bond gives an extra stabilization of 2-3.5 kcal mol-1. The same is observed for methylammonium complexes, with changes of -17 kcal mol-1 and -15 kcal mol-1 in the presence and absence of the water···phenol hydrogen bond, respectively. The energy gain is smaller than in ammonium complexes, because in the formation of the water···phenol hydrogen bond, a contact between the methyl group and phenol has to be broken. Considering the values obtained for ∆Ehyd corresponding to the interaction between phenol and the rest of the complex considered as a single unit, it can be appreciated that the values registered are similar to those obtained in the absence of water (of course are less negative than complexation energies obtained for complexes with one water molecule since the cation···water interaction is not included). Therefore, contrary to the usual trends observed in other systems where the presence of one water molecule decreases the strength of the interaction because it competes with the aromatic molecule for interacting with the cation,26-30 in minima PA1-1 and PM1-1 ∆Ehyd amounts to -18.6 and -18.0 kcal mol-1, respectively. Though the cation··· interaction is weakened as shown by the values obtained for PA1-3 and PM1-3 which exhibit decreases in the interaction strength of more than 3 kcal mol-1, the formation of the O-H···O hydrogen bond in PA1-1 and PM1-1 compensates for this effect. In benzene complexes, the decrease in strength amounts to around 1.7 kcal mol-1 when the first water molecule is included.29 The inclusion of the second water molecule increases the complexity of the potential energy surface with lots of minima obtained. The six more stable ones for ammonium and methylammonium complexes are shown in Figure 4.3, with their complexation energies listed in Table 4.3. It can be appreciated that most of the minima shown in Figure 4.3 present a cyclic pattern of hydrogen bonds similar to those observed for the complexes containing one water molecule. Therefore, two different patterns arise: one presenting a ···H-N-H···O-H···O- hydrogen bond network, and the other with a series of -O···H-N-H···O-H··· contacts. Most of the structures in Figure 4.3 present these patterns (exceptions are PA2-5 and PM2-4), with the second water molecule occupying one of the free N-H units of ammonium cations, or hydrogen bonded as acceptor to the hydroxyl group of phenol. Therefore, as the second water molecule is included, the hydroxyl group starts participating in the hydrogen bond network of several of the most stable structures. This is an indication that the energy differences between coordinating the water molecule to the cation or to the hydroxyl group have diminished, being competitive with each other. 4.3. Results 91 Figure 4.3. Selected most stable minima for the complexes formed by ammonium and methylammonium with phenol in the presence of two water molecules as obtained at the MP2/6-31+G(2d,p) level of calculation. Table 4.3 lists the values obtained for the complexation energies of the minima shown in Figure 4.3. Considering ammonium complexes, it can be observed that the two most stable structures differ in stability by only 0.5 kcal mol-1, so including more water molecules decreases the difference in stability between structures presenting O-H···O and O-H··· contacts. The next structure in order of stability PA2-3 already presents a -OH···O hydrogen bond, and is 1.6 kcal mol-1 less stable than the most stable minimum found, and isoenergetic with the analogous structure bearing a O-H··· contact PA2-4. Minimum PA2-5 is unique since no contact between water molecules and the phenol moiety is established, with a drop in stability of 2.5 kcal mol-1 with respect to the most stable minimum. Finally, in PA2-6 the second water PA2-1 PA2-2 PM2-1 PA2-3 PA2-4 PA2-5 PA2-6 PM2-2 PM2-3 PM2-5 PM2-4 PM2-6 92 Chap. 4. Microhydration study: complexes between NH4+ and CH3NH3+ with phenol molecule is hydrogen bonded to the previous one. The complexation energy for this structure is 2.9 kcal mol-1 less negative than that of the most stable minimum. Therefore, when a second water molecule is included in the cluster, there are three main possibilities: bonding to a N-H group, bonding to the hydroxyl oxygen and bonding to the previous water molecule. As deduced from the comparison of PA2-1, PA2-3 and PA2-6 the first option is the most favorable, the second being 1.6 kcal mol-1 less stable, whereas the third one decreases complexation energy by 2.9 kcal mol-1 with respect to bonding to the N-H group, or 1.3 kcal mol-1 with respect to hydroxyl oxygen. Table 4.3. Complexation energies (kcal mol-1) obtained for the most stable minima of the clusters containing two water molecules (see Figure 4.3) as obtained at the MP2/6-31+G(2d,p) level. Ammonium Methylammonium ∆Ecomplex ∆EZPE ∆H298 ∆Ehyd ∆Ecomplex ∆EZPE ∆H298 ∆Ehyd PA2-1 -51.18 -46.38 -47.08 -15.83 PM2-1 -48.14 -40.50 -40.72 -16.00 PA2-2 -50.57 -46.01 -46.61 -15.16 PM2-2 -47.20 -38.81 -39.58 -14.55 PA2-3 -49.59 -44.06 -45.31 -13.81 PM2-3 -46.97 -39.70 -39.67 -14.80 PA2-4 -49.43 -44.45 -45.44 -13.56 PM2-4 -46.48 -38.56 -39.12 -14.56 PA2-5 -48.76 -44.59 -44.80 -13.31 PM2-5 -46.29 -38.37 -38.88 -13.61 PA2-6 -48.31 -43.13 -44.28 -13.25 PM2-6 -45.72 -37.58 -38.28 -13.86 In the case of methylammonium clusters, the behavior is pretty similar, though some differences arise. The most stable structure PM2-1 is analogous to the most stable minimum of ammonium complexes. However, the second most stable minimum, with a complexation energy only 0.8 kcal mol-1 less negative already presents a -OH···O contact, with methylamonium over the phenyl ring. It is worth noting that whereas in PM2-1 all N-H groups of the cation are occupied, in PM2-2 there is a free one, with the small impact in energies indicating that complexation with the cation or phenol already gives similar stabilization to the complex. In fact there is an inversion in the order of stability between structures PM2-2 and PM2-3 with respect to that observed in ammonium complexes, though the difference between these two structures amounts to only 0.3 kcal mol-1. Structure PM2-4 is also very close in energy, being 1.6 kcal mol-1 less stable than the global minimum. In this structure both water molecules are bound to methylammonium N-H groups, establishing both O-H···O and O-H··· contacts with phenol. Therefore in this structure the interaction takes place between a hydrated methylammonium 4.3. Results 93 cation and a phenol molecule, with no direct interaction between aromatic molecule and cation, in a similar way to that observed in benzene clusters with three water molecules.29 Finally, in PM2-6 the new water molecule is hydrogen bonded to the previous one, with a stability loss of 2.4 kcal mol-1 with respect to the most stable minimum. So, even when the behavior is similar with both cations, in methylammonium complexes there are smaller differences for water to be coordinated to any of the favorable locations within the cluster, especially between N-H and phenol O-H group. In complexes with two water molecules, inclusion of ZPE does alter the order of stability already discussed, though when this happens it is a consequence of the structures being almost isoenergetic. ∆Ehyd values shown in Table 4.3 are significantly less negative than those obtained in complexes with one water molecule. Since the second water molecule does not directly interact with phenol, there is no compensating effect for the decrease in the cation··· interaction as a consequence of the cation charge being shared among all neutral species in the complex. In fact, the charges obtained from a natural bond orbital (NBO) analysis indicate that the charge of the ammonium cation amounts to 0.93 in the ammonium···phenol complex, decreasing to 0.91, 0.89 and 0.89 when including up to three water molecules (values for methylammonium are 0.94 a.u. in the complex with phenol and 0.92, 0.89 and 0.89 a.u. as water is included). Therefore, ∆Ehyd drops by about 2 kcal mol-1 with respect to the values observed in complexes with one water molecule. The inclusion of a third molecule extremingly complicates the search for minima, since an overwhelming amount of minima can be located, the most stable of which are shown in Figure 4.4. The corresponding complexation energies are shown in Table 4.4. In the case of ammonium complexes, most structures show the cyclic hydrogen bond pattern observed in the smaller clusters, but other possibilities arise presenting a variety of hydrogen bond networks, with water molecules interacting among themselves (PA3-2 and PA3-8, for example), or participation of the hydroxyl group. The energy differences among structures are very small, with up to 7 minima within a complexation energy interval of only 1 kcal mol-1. In any case the most stable minimum exhibits the same structure observed in smaller clusters with the third water molecule occupying the last free N-H group of ammonium cation. A similar pattern is observed in PA3-2, only 0.2 kcal mol-1 less stable. However, PA3-3, with a complexation energy only 0.4 kcal mol-1 less negative than the most stable structure, already shows a -OH···O hydrogen bond, revealing a further decrease on the difference between coordinating to the cation or to the hydroxyl group. Therefore, already with three water molecules, there are plenty of structures with similar stability because the stability differences of the favorable regions for water coordination have diminished to be almost negligible. Due to the simultaneous participation of phenol hydroxyl group and aromatic ring in the hydrogen bond network, the structural patterns are similar to those found in ammonium···water clusters with n+2 water molecules.51 94 Chap. 4. Microhydration study: complexes between NH4+ and CH3NH3+ with phenol Figure 4.4. Selected most stable minima for the complexes formed by ammonium and methylammonium with phenol in the presence of three water molecules as obtained at the MP2/6-31+G(2d,p) level of calculation. The results obtained for methylammonium complexes are similar, though only four structures are within a 1 kcal mol-1 range from the most stable minimum. The main difference with ammonium complexes arises because in methylammonium cluster already with two water molecules and phenol there are no N-H free groups for the third water to be coordinated, so it must be incorporated to the hydrogen bond network. Therefore, the most stable structure already shows a -OH···O hydrogen bond, as also does PM3-3, the third most stable minimum. The rest of stable structures present O-H···O hydrogen bonds among water molecules or with water acting as hydrogen bond donor to the hydroxyl group. It becomes clear that including more water molecules already implies the formation of hydrogen bonds between water molecules on a second solvation shell or necessarily introduces -OH···O contacts. Almost any position PM3-1 PM3-11 PM3-12 PM3-5 PM3-6 PM3-7 PM3-8 PM3-9 PM3-10 PM3-3 PM3-4 PM3-2 PA3-1 PA3-11 PA3-12 PA3-5 PA3-6 PA3-7 PA3-8 PA3-9 PA3-10 PA3-3 PA3-4 PA3-2 4.3. Results 95 occupied by the water molecule will lead to minima with complexation energies of similar magnitude. In fact, analyzing the energy changes when including water molecules to the phenol···cation complexes, it becomes clear that the stabilization drops significantly when water molecules are included in the complex. Thus, when the first water molecule is included, the interaction changes by -17.9 kcal mol-1 in ammonium complexes and by -17.6 kcal mol-1 in methylammonium ones, as a consequence of both a new N-H···O contact but also due to the presence of a new O-H···O hydrogen bond to the hydroxyl group. The inclusion of the second water molecule stabilizes the complex in a significantly smaller amount, reaching -13.4 and -12.4 kcal mol-1 for ammonium and methylamonium complexes, respectively. So, the stabilization drops by 4.5 kcal mol-1 in ammonium and by 5.2 kcal mol-1 in methylammonium. This happens because the inclusion of the second water molecule only introduces a new N-H···O contact. Table 4.4. Complexation energies (kcal mol-1) obtained for the most stable minima of the clusters containing three water molecules (see Figure 4.4) as obtained at the MP2/6-31+G(2d,p) level. Ammonium Methylammonium ∆Ecomplex ∆EZPE ∆H298 ∆Ehyd ∆Ecomplex ∆EZPE ∆H298 ∆Ehyd PA3-1 -63.06 -56.79 -57.42 -13.31 PM3-1 -59.18 -49.47 -50.03 -13.54 PA3-2 -62.91 -55.65 -56.94 -13.59 PM3-2 -58.82 -49.07 -49.75 -14.43 PA3-3 -62.73 -56.56 -57.88 -13.60 PM3-3 -58.27 -48.17 -49.05 -13.14 PA3-4 -62.62 -53.84 -56.25 -13.67 PM3-4 -58.14 -46.35 -48.29 -13.30 PA3-5 -62.55 -56.43 -57.00 -13.04 PM3-5 -58.02 -48.36 -48.97 -13.10 PA3-6 -62.30 -55.83 -56.75 -12.98 PM3-6 -57.88 -48.68 -48.93 -12.19 PA3-7 -62.10 -55.53 -56.49 -12.10 PM3-7 -57.86 -48.72 -48.99 -12.99 PA3-8 -61.90 -54.10 -55.82 -12.14 PM3-8 -57.86 -47.90 -48.70 -13.07 PA3-9 -61.71 -54.69 -55.96 -12.46 PM3-9 -57.70 -48.52 -48.81 -13.28 PA3-10 -61.19 -53.59 -55.26 -11.91 PM3-10 -57.28 -47.41 -48.26 -11.82 PA3-11 -61.07 -53.92 -55.43 -11.31 PM3-11 -57.14 -47.69 -48.16 -12.22 PA3-12 -61.05 -54.25 -55.37 -11.83 PM3-12 -56.93 -47.52 -47.97 -11.98 96 Chap. 4. Microhydration study: complexes between NH4+ and CH3NH3+ with phenol Finally, when the third water molecule is included, the energy change amounts to -11.9 kcal mol-1 for ammonium and -11.5 kcal mol-1 for methylammonium, with an extra drop of 1.5 and 0.9 kcal mol-1, respectively. These values indicate that the first water molecule is tightly bound within the complex whereas as more water molecules are included they are more loosely held in the cluster. The same trends are observed in ∆Ehyd values, which become less negative as the more water molecules are included. The values for ∆Ehyd of the most stable complexes of ammonium amount to -19.9, -18.6, -15.8 and -13.3 kcal mol-1 for complexes from 0 to 3 water molecules. It becomes clear that the first water molecule hardly affects the interaction with phenol due to the compensation of the loss in cation···phenol interaction by the O-H···O hydrogen bond. However, as more water molecules are included the interaction strength changes by larger quantities, exceeding 2 kcal mol -1. In the case of methylammonium complexes (-18.2, -18.0, -16.0 and -13.5 kcal mol-1) the effect is similar, though in this case the first water molecule is able of totally recovering the loss of strength in the interaction. The inclusion of more water molecules produces a decrease in the interaction strength of more than 2 kcal mol-1. Therefore, as observed in other systems,26-30 the inclusion of water weakens the cation··· interaction. Nevertheless, the participation of the hydroxyl group interacting as hydrogen acceptor to one water molecule makes the weakening only evident in complexes with at least two water units. 4.4. Conclusions Clusters formed by one phenol molecule and an ammonium or methylamonium cation in the presence of up to three water molecules have been computationally studied at the MP2/6-31+G(2d,p) level of calculation. Both ammonium and methylammonium form complexes interacting with the aromatic ring and the hydroxyl group of phenol with similar stabilities. However, in methylammonium complexes, secondary interactions are established between the methyl group and phenol. The presence of water molecules greatly increases the complexity of the potential energy surfaces of the clusters though the minima located show similar characteristics for both cations. In any case, as one water molecule is incorporated to the system, the most stable minima present cyclic patterns with the water molecule bound to the cation and simultaneously establishing a hydrogen bond with the phenol molecule via the hydroxyl group or the aromatic ring. The inclusion of more water molecules does not break this pattern which is observed in the most stable structures of all clusters studied. 4.5. Refrences 97 As more water molecules are included, the energy differences of the favorable interaction sites allowed to the new water molecule (contact with N-H of the cation, O-H of phenol or another water molecule) become smaller, so in clusters with two and three water molecules, several of the most stable minima present a -OH···O hydrogen bond between the phenol hydroxyl group and a water molecule. The energy change upon formation of a complex with n water molecules from the n-1 one decreases as more water molecules are included. The stabilization is especially significant for the first water molecule since it interacts simultaneously with the cation and the phenol molecule. The incorporation of the second and third molecules is accompanied by significantly smaller changes. Therefore, though the presence of water weakens the phenol···cation interaction, at least two water molecules are needed to produce a noticeable effect. The first water molecule partially recovers the loss in phenol···cation interaction by means of the hydrogen bond to phenol oxygen. The results obtained in the present study can help understanding the interaction between ammonium cations and the side chain of tyrosine, especially in environments where the amino acid is only partially exposed to the solvent. 4.5. References (1) L. M. Salonen, M. Ellermann, F. Diederich. Angew. Chem. Int. Ed. 2011, 50, 4808-4842. (2) P. Hobza, R. Zaradnik, Intermolecular complexes: the role of van der Waals systems in physical chemistry and the biodisciplines, Elsevier, Amsterdam, 1988. (3) E. A. Meyer, R. K. Castellano, F. Diederich. Angew. Chem. Int. Ed. 2003, 42, 1210-1250. (4) J.-M. Lehn, Supramolecular Chemistry: concepts and perspectives, VCH, Weinheim, 1995. (5) F. Voegtle, Editor, Supramolecular Chemistry: An Introduction, Maruzen Co., Ltd., New York, 1995. (6) F. Voegtle, Editor, Comprehensive Supramolecular Chemistry, Volume 2: Molecular Recognition: Receptors for Molecular Guests, Pergamon, Oxford, 1996. (7) J. C. Ma, D. A. Dougherty. Chem. Rev. 1997, 97, 1303-1324. (8) N. S. Scrutton, A. R. Raine. Biochem. J. 1996, 319, 1-8. (9) J. P. Gallivan, D. A. Dougherty. Proc. Nart. Acad. Sci. USA 1999, 96, 9459-9464. (10) D. A. Dougherty. J. Nutr. 2007, 137, 1504S-1508S. (11) M. L. Waters. Biopolymers (Peptide Science) 2004, 76, 435-445. (12) J. P. Gallivan, D. A. Dougherty. J. Am. Chem. Soc. 2000, 122, 870-874. (13) M. A. Anderson, B. Ogbay, R. Arimoto, W. Sha, O. G. Kisselev, D. P. Cistola, G. R. Marshall. J. Am. Chem. Soc. 2006, 128, 7531-7541. 98 Chap. 4. Microhydration study: complexes between NH4+ and CH3NH3+ with phenol (14) B. W. Berry, M. M. Elvekrog, C. Tommos. J. Am. Chem. Soc. 2007, 129, 5308-5309. (15) R. M. Hughes, M. L. Benshoff, M. L. Waters. Chem. Eur. J. 2007, 13, 5753-5764. (16) R. M. Hughes, M. L. Waters. J. Am. Chem. Soc. 2005, 127, 6518-6519. (17) R. M. Hughes, M. L. Waters. J. Am. Chem. Soc. 2006, 128, 13586-13591. (18) R. M. Hughes, M. L. Waters. J. Am. Chem. Soc. 2006, 128, 12735-12742. (19) H. Khandelia, Y. N. Kaznessis. J. Phys. Chem. B 2007, 111, 242-250. (20) P. E. Mason, C. E. Dempsey, G. W. Neilson, S. R. Kline, J. W. Brady. J. Am. Chem. Soc. 2009, 131, 16689-16696. (21) A. J. Riemen, M. L. Waters. Biochemistry 2009, 48, 1525-1531. (22) Z. Shi, C. A. Olson, N. R. Kallenbach. J. Am. Chem. Soc. 2002, 124, 3284-3291. (23) C. D. Tatko, M. L. Waters. Protein Sci. 2003, 12, 2443-2452. (24) K. E. Riley, P. Hobza. WIREs Comput. Mol. Sci. 2011, 1, 3-17. (25) C. D. Sherrill, Computations of Non covalent π Interactions, in Rev. Comput. Chem., Vol., John Wiley & Sons, Inc., 2009, pp.1-38. (26) C. Adamo, G. Berthier, R. Savinelli. Theor. Chem. Acc. 2004, 111, 176-181. (27) A. S. Reddy, H. Zipse, G. N. Sastry. J. Phys. Chem. B 2007, 111, 11546-11553. (28) N. J. Singh, S. K. Min, D. Y. Kim, K. S. Kim. J. Chem. Theory Comput. 2009, 5, 515-529. (29) Y. Xu, J. Shen, W. Zhu, X. Luo, K. Chen, H. Jiang. J. Phys. Chem. B 2005, 109, 59455949. (30) E. M. Cabaleiro-Lago, J. Rodríguez-Otero, Á. Peña-Gallego. J. Chem. Phys. 2011, 135, 214301/1-214301/9. (31) D. Feller, M. W. Feyereisen. J. Comput. Chem. 1993, 14, 1027-1035. (32) H. Watanabe, S. Iwata. J. Chem. Phys. 1996, 105, 420-431. (33) M. Gerhards, K. Kleinermanns. J. Chem. Phys. 1995, 103, 7392-7400. (34) R. Wu, B. Brutschy. Chem. Phys. Lett. 2004, 390, 272-278. (35) D. M. Benoit, D. C. Clary. J. Phys. Chem. A 2000, 104, 5590-5599. (36) T. Ebata, A. Fujii, N. Mikami. Int. J. Mass Spectrom. Ion Processes 1996, 159, 111-124. (37) C. Janzen, D. Spangenberg, W. Roth, K. Kleinermanns. J. Chem. Phys. 1999, 110, 98989907. (38) W. Roth, M. Schmitt, C. Jacoby, D. Spangenberg, C. Janzen, K. Kleinermanns. Chem. Phys. 1998, 239, 1-9. (39) M. S. Marshall, R. P. Steele, K. S. Thanthiriwatte, C. D. Sherrill. J. Phys. Chem. A 2009, 113, 13628-13632. (40) H. M. Lee, P. Tarakeshwar, J. Park, M. R. Kołaski, Y. J. Yoon, H.-B. Yi, W. Y. Kim, K. S. Kim. J. Phys. Chem. A 2004, 108, 2949-2958. (41) G. Chałasiński, M. M. Szczȩśniak. Chem. Rev. 2000, 100, 4227-4252. (42) S. F. Boys, F. Bernardi. Mol. Phys. 1970, 19, 553-566. 5.2. Computational details 105 Complexation energies for the minima thus located have been obtained at the M06-2X/6-31+G* and the MP2/aug-cc-pVDZ levels of calculation. Taking into account that MP2 is known to overbind clusters where aromatic units are present,31-33 especially when stacked, empirical corrected MP2 variants have been considered. These so-called Spin Component Scaled MP2 methods are based on an empirical scaling of the contributions of parallel and antiparallel electron pairs to the correlation energy.34 Therefore, the scaled MP2 energy is obtained as: )(2,)(2,RHF E ps Ept + E = X-SCS    , (eq. 5.1) where pt = 0.33 ; ps = 1.20 for SCS-MP234 and pt = 1.76 ; ps = 0 for SCSN-MP235. Of course for the original MP2 both parameters are 1. In order to avoid the basis set superposition error (BSSE),36 the counterpoise method has been used to obtain the complexation energy of all the minima located.36, 37 Therefore, interaction energies are obtained employing the full basis set of the complex as:   i complex i complex ijkEijkEE ...)(...)( int . (eq. 5.2) Terms in parentheses indicate the basis set while superscripts refer to the geometry employed in the calculations. As the geometry of the molecules changes when the cluster is formed, an additional contribution describing this effect must be included, obtained as the energy difference between the molecules in the cluster geometry and in isolation.    i isolated i complex idef iEiEE )()( . (eq. 5.3) Finally, the complexation energy is the combination of these two quantities: defcomplex EEE  int . (eq. 5.4) In order to obtain more information about the nature of the interaction with the different aromatic systems, Symmetry Adapted Perturbation Theory38, 39 calculations based on DFT (SAPT(DFT)) have been carried out for the complexes without water.40, 41 Since many body effects are not easily treated with this method, the analysis has been restricted to non-hydrated systems. These calculations have been performed at the LPBE0AC/aug-cc-pVDZ level employing Molpro.42 This involves a correction by adding a shift to the asymptotic part of the potential. This shift is obtained as the sum of the ionization potential and the energy of the 106 Chap. 5. Microhydration study: complexes between guanidinium cation and aromatic species highest occupied molecular orbital, as obtained at the PBE0/aug-cc-pVDZ level of calculation. The DFT-SAPT calculations were performed with the aug-cc-pVDZ basis set, employing the ccpVTZ/JKFIT for Hartree–Fock and aug-cc-pVDZ/MP2FIT for the second-order dispersion terms. The rest of the calculations have been performed by using Gaussian09.43 5.3. Results After optimization a variety of minima have been obtained, their number quickly growing as more water molecules are included. Therefore, results will be presented for a selection of the most stable minima found, numbered and sorted by their stability at the MP2/aug-cc-pVDZ level. Figure 5.1 shows the minima found for the complexes formed by the guanidinium cation and the aromatic systems employed in this work: benzene, indole and phenol, whereas Table 5.1 list the values obtained for complexation energies. Figure 5.1. Structures of the minima obtained for complexes containing guanidinium cation and the aromatic molecules. Selected distances are shown in Å. BG0-1 IG0-1 IG0-2 IG0-3 PG0-1 PG0-2 PG0-3 PG0-4 2.297 2.264 2.208 2.226 2.657 2.799 1.883 2.329 1.976 1.977 2.310 2.858 2.310 2.743 5.3. Results 107 One minimum has been found for guanidinium-benzene complex, corresponding to a T-shaped structure where guanidinium cation interacts with benzene by means of a double contact with the NH2 groups, in agreement with previously published results.22, 30 In this minimum BG0-1 (the nomenclature reflects the benzene B, indole I, or phenol P unit, guanidinium G cation, the number of water molecules present in the system (0 in this case), and a numerical identifier for each structure) guanidinium interacts with the C-C midbonds, establishing a N-H···π contact at about 2.30 Å to the center of the ring. The complexation energy amounts to around -14 kcal mol-1. In the case of indole complexes three different structures have been found. IG0-1 corresponds to a minimum with guanidinium in a T-shaped structure interacting simultaneously with both rings of indole; IG0-2 is similar, but the cation is rotated and interacts only with the phenyl ring of indole; IG0-3 corresponds to a parallel-stacked structure with guanidinium interacting with both rings of indole. In the latter case, the distances from guanidinium to the ring center are longer, since there are no N-H···π hydrogen bonds (see Figure 5.1). As regards complexation energies IG0-1, as expected, is the most stable minimum, reaching around -21 kcal mol-1 with the different methods. That is, the presence of a larger aromatic molecule increases the strength of the interaction by about -7 kcal mol-1 with respect to the structure with benzene. IG0-2 is about 2-3 kcal mol-1 less stable, whereas IG0-3 is the least stable of the three minima, with a complexation energy of -16 kcal mol-1. Table 5.1. Complexation energies (kcal mol-1) for the complexes formed by guanidinium cation and the aromatic molecules considered in this work as obtained at the M06-2X/6-31+G* and different versions of MP2/aug-cc-pVDZ. ∆EM06-2X ∆EMP2 ∆ESCS-MP2 ∆ESCS(N)-MP2 BG0-1 -13.42 -14.19 -12.39 -13.83 IG0-1 -20.51 -21.26 -18.70 -20.70 IG0-2 -17.93 -18.81 -16.58 -18.24 IG0-3 -17.06 -16.40 -13.72 -15.06 PG0-1 -17.20 -17.42 -15.56 -17.07 PG0-2 -17.22 -16.13 -14.85 -15.96 PG0-3 -14.04 -14.85 -12.97 -14.55 PG0-4 -13.68 -12.34 -10.09 -11.23 108 Chap. 5. Microhydration study: complexes between guanidinium cation and aromatic species Phenol complexes with guanidinium show the four minima presented in Figure 5.1. It can be observed that PG0-1, PG0-3 and PG0-4 approximately correspond to the minima found for indole complexes, with the hydroxyl group of phenol playing the role of the second aromatic ring in indole. Besides, there is a minimum, PG0-2, showing guanidinium interacting only with the hydroxyl oxygen by means of two equivalent N-H···O hydrogen bonds. The complexation energies in Table 5.1 show that the most stable minimum is PG0-1, reaching -17 kcal mol-1, but closely followed by PG0-2. In fact, in the case of the M06-2X/6-31+G* results, PG0-2 is as stable as PG0-1, but it can be observed that, compared to MP2 values, M06-2X seems to overestimate the interaction in structures showing N-H···O hydrogen bonds. As in indole complexes, PG0-3 and the parallel stacked PG0-4 are less stable by a couple of kcal mol-1. Overall it can be observed that the strength of the interaction increases from benzene to phenol to indole, as a consequence of the double simultaneous interaction with the ring and the larger size of indole. As regards the method employed it can be confirmed that the results obtained with any of the methods are pretty similar with the exception of SCS-MP2/aug-cc-pVDZ, which tends to underestimate the strength of the interaction compared with the others. Also, there are almost no differences between the regular MP2 and SCSN-MP2, and both agree quite well with the M06-2X/6-31+G* results. In order to understand why the interaction increases from benzene to indole, SAPT(DFT) calculations have been performed for the minima shown in Figure 5.1. The results are displayed in Figure 5.2. First, in benzene complex it becomes clear that the main contribution to the interaction energy is electrostatic (-10.0 kcal mol-1). This reflects the N-H···π hydrogen bonds formed by the charged guanidinium cation and the benzene ring. Also, as a consequence of a large electrostatic interaction, a significant induction contribution is also observed (-8.4 kcal mol-1), as a result of the deformation of the polarizable π cloud by the charged species nearby. However, the contribution of dispersion is significant, reaching –6.3 kcal mol-1, despite the T-shaped orientation of guanidinium relative to benzene. Finally, a large repulsion contribution of a magnitude similar to the electrostatic one is found. These results agree with those previously published,22, 30 and confirm the important role of dispersion in this kind of complexes. As comparison, in benzene-sodium complexes, the dispersion contribution barely reaches -1.5 kcal mol-1.10, 44 5.3. Results 109 Figure 5.2. SAPT(DFT) energy decomposition for the complexes shown in Figure 5.1. Indole complexes behave quite differently. As observed in Figure 5.2, the contribution from electrostatics is larger than that observed in benzene for any of the minima shown in Figure 5.1. That is, the presence of a permanent dipole in indole produces a stronger electrostatic contribution that reaches -15.6 kcal mol-1 in the most stable complex. Quite surprisingly, the electrostatic contribution is almost as large in the parallel-stacked minimum IG0-3, of about -15.0 kcal mol-1, whereas decreases when guanidinium cation interacts only with the phenyl ring in indole (-13.9 kcal mol-1), though in all cases is still larger than in benzene complex. As expected, 110 Chap. 5. Microhydration study: complexes between guanidinium cation and aromatic species taking into account the larger size of the polarizable aromatic cloud, induction contributions are large, amounting to –11.6 kcal mol-1 in the most stable complex, and decreasing to -10.0 kcal mol-1 in IG0-2. However, in the case of the parallel stacked minimum IG0-3 there is an important decrease of polarization, which does not reach -7 kcal mol-1 and therefore is similar to that in benzene. Dispersion is also larger than in benzene complex, amounting to –8.9 and -7.4 kcal mol-1 in minima IG0-1 and IG0-2, respectively. However, in minimum IG0-3, the parallel orientation of guanidinium relative to indole produces a larger contribution of dispersion, which is the largest among the minima considered, reaching -10.2 kcal mol-1, and partially compensating the smaller contribution of polarization. Repulsion is huge in all structures though it favors IG0-2. Phenol complexes display similar characteristics. The electrostatic contribution is large, reaching -15.9 kcal mol-1 in the most stable complex, but amounting to -17.9 kcal mol-1 in PG0-2, the largest electrostatic contribution among minima in Figure 5.1. This is a consequence of the double N-H···O hydrogen bond formed in this minimum. The smaller electrostatic contribution, as in indole, occurs in PG0-3 where guanidinium interacts only with the phenyl ring. In accordance with the large electrostatic contribution, induction terms are also important, amounting to around -10 kcal mol-1 in the most stable structures. As in indole complexes, the parallel orientation of guanidinium produces a decrease in polarization contribution, which only amounts to -5.7 kcal mol-1. Dispersion is also large, being midway between the values observed for benzene and indole, and reaching -7.1 kcal mol-1 in the most stable complex. As in indole, the largest contribution is observed in the parallel PG0-4, amounting to -8.8 kcal mol-1 and compensating the loss in induction. Overall, it becomes clear that in the interaction of guanidinium with aromatic molecules, the electrostatic contribution is usually the leading term, but with important contributions from induction and dispersion. Moving from benzene to phenol to indole, all contributions to the interaction energy are increased, especially electrostatic and dispersion ones. This latter contribution is more significant in parallel minima, so it is crucial in order to explain the formation of stacked structures. Figure 5.3 shows the selected most stable minima found for the complexes formed by guanidinium and each of the aromatic species studied in the presence of one water molecule. Guanidinium can easily accommodate up to three water molecules by means of double NH···O hydrogen bonds, and at most only one of these positions is blocked in the minima in Figure 5.1. Thus, it can be expected the water molecule to be bonded to guanidinium cation whereas the cation itself interacts with the aromatic unit as in Figure 5.1. However, as already commented before, in benzene complexes a new possibility arises, with water intercalated between the cation 5.3. Results 111 and the aromatic ring. The most stable structure (Table 5.2) simply corresponds to BG0-1 with one water molecule coordinated to guanidinium, reaching -29 kcal mol-1. The other structure BG1-2 with water hydrogen bonded to benzene is less stable by about 4 kcal mol-1. Complexes with indole are more difficult to analyze since more structures are possible, so only a selection of the most stable ones are shown in Figure 5.3. The most stable structure with MP2/aug-cc-pVDZ corresponds, as in benzene, to IG0-1 plus one water molecule coordinated to guanidinium (IG1-1), whereas the cation···π contact remains almost unperturbed. In the case of IG1-2 and IG1-3, both come from IG0-3, with the cation oriented in parallel to the aromatic ring and one water molecule bonded to it. However, since the aromatic cloud of indole is more extended than that of benzene, the water molecule also interacts simultaneously with the other ring of indole (with the pyrrol ring in IG1-2 and the phenyl one in IG1-3). Figure 5.3. Structures of selected minima obtained for complexes containing one water molecule. PG1-1 PG1-2 PG1-3 PG1-4 PG1-5 IG1-1 IG1-2 IG1-3 IG1-4 IG1-5 BG1-1 BG1-2 112 Chap. 5. Microhydration study: complexes between guanidinium cation and aromatic species Table 5.2. Complexation energies (kcal mol-1) for selected complexes containing one water molecule as obtained at the M06-2X/6-31+G* and different versions of MP2/aug-cc-pVDZ. ∆EM06-2X ∆EMP2 ∆ESCS-MP2 ∆ESCS(N)-MP2 BG1-1 -30.79 -29.05 -26.42 -29.22 BG1-2 -28.38 -25.23 -22.39 -25.08 IG1-1 -37.01 -35.29 -31.93 -35.21 IG1-2 -35.26 -32.54 -28.59 -31.90 IG1-3 -34.25 -31.40 -27.44 -30.78 IG1-4 -34.20 -31.10 -27.43 -30.28 IG1-5 -31.12 -30.82 -27.35 -30.52 PG1-1 -34.04 -31.76 -29.00 -31.89 PG1-2 -34.41 -30.80 -28.64 -31.08 PG1-3 -33.90 -30.11 -26.64 -29.71 PG1-4 -30.89 -29.14 -26.14 -29.25 PG1-5 -31.91 -28.20 -25.32 -28.16 With any of the methods employed, coordination of guanidinium to the phenyl ring is favored. The same happens in phenol complexes with one water molecule. The most stable minima PG1-1 and PG1-2 correspond to PG0-1 and PG0-2 with the water molecule coordinated to one free position in guanidinium. About 2 kcal mol-1 less stable there is a structure like PG1-5 with guanidinium interacting with phenol by only one NH2 group. Also, it is worth noting that when water coordinates to the phenol hydroxyl group by means of a hydrogen bond (PG1-4), the stability is similar, so this kind of structures start being competitive with others where the water molecule is directly coordinated to the cation. Comparing PG1-1 and PG1-4 it can be observed that the energy difference between coordinating to the cation or the hydroxyl group only introduces an energy difference of 3 kcal mol-1 (as already observed in other cation-phenol complexes23, 27). In any case, in phenol complexes, the three most stable structures are virtually isoenergetic even though the interaction pattern is totally different, with guanidinium interacting with the aromatic ring, with the hydroxyl group or in a parallel orientation with water hydrogen bonded to phenol. 5.3. Results 113 Figure 5.4. Structures of selected minima obtained for complexes containing two water molecules. Including the second water molecule increases the number of minima with similar complexation energies. A selection of these is shown in Figure 5.4, with complexation energies listed in Table 5.3. In benzene complex, the most stable minimum found (-43 kcal mol-1) corresponds to the T-shaped BG0-1 minimum with the two water molecules occupying the free coordination sites of guanidinium. It would be expected that most stable structures with two water molecules will present similar characteristics, with water molecules interacting directly with the cation. However, a minimum like BG2-2 is less stable than BG2-1 by 2 kcal mol-1 but with similar stability as other structures as BG2-3 and BG2-4. The peculiarity is that in BG2-2, one of the water molecules forms a hydrogen bond to the other water unit without interacting directly with the cation. That is, already with two water molecules, hydrogen bonding between water molecules starts being competitive with the cation···water interaction. A similar behavior is observed in complexes with indole. In this case, the most stable structure corresponds to IG2-1 with the two water molecules coordinated to guanidinium and complexation energy of -48 kcal mol-1. Other possibilities with water interacting with the ring (IG2-2) or with water···water hydrogen bond (IG2-4) are around 2 kcal mol-1 less stable. In phenol complexes three structures present the same stability (around -45 kcal mol-1), corresponding to contacts of guanidinium with the phenyl ring and hydroxyl group, double contact with hydroxyl, and finally BG2-1 BG2-3 BG2-4 BG2-5BG2-2 IG2-1 IG2-3 IG2-4 IG2-5 IG2-6IG2-2 PG2-1 PG2-3 PG2-4 PG2-5 PG2-6PG2-2 114 Chap. 5. Microhydration study: complexes between guanidinium cation and aromatic species a parallel structure with a OH···O hydrogen bond to the hydroxyl group. These minima correspond to the most stable monohydrated clusters, with the second water molecule in the free NH2 groups of guanidinium cation. Table 5.3. Complexation energies (kcal mol-1) for selected complexes containing two water molecules as obtained at the M06-2X/6-31+G* and different versions of MP2/aug-cc-pVDZ. ∆EM06-2X ∆EMP2 ∆ESCS-MP2 ∆ESCS(N)-MP2 BG2-1 -46.54 -42.52 -39.12 -43.06 BG2-2 -43.52 -39.77 -36.39 -40.45 BG2-3 -44.88 -39.34 -35.70 -39.61 BG2-4 -44.97 -39.17 -36.10 -39.96 BG2-5 -40.20 -35.45 -31.84 -35.70 IG2-1 -52.05 -48.06 -43.96 -48.34 IG2-2 -51.44 -46.43 -41.76 -46.27 IG2-3 -49.74 -46.14 -42.43 -46.55 IG2-4 -49.38 -45.70 -41.61 -46.12 IG2-5 -50.22 -45.35 -40.75 -45.08 IG2-6 -47.16 -44.40 -40.15 -44.56 PG2-1 -49.44 -44.84 -41.36 -45.36 PG2-2 -50.06 -44.08 -41.08 -44.72 PG2-3 -49.86 -44.05 -39.93 -44.21 PG2-4 -47.20 -43.05 -39.24 -43.63 PG2-5 -49.84 -42.75 -38.25 -42.82 PG2-6 -46.57 -42.32 -38.79 -42.94 Guanidinium cation can establish three simultaneous contacts by means of its NH2 groups. In the most stable minima discussed so far one of these positions is occupied by the aromatic molecules, whereas the other two are available for accommodating up to two water molecules. The inclusion of the third water molecule changes this behavior because now at least one of the water molecules or the aromatic unit cannot interact directly with the NH2 groups of guanidinium. Figure 5.5 shows the most stable minima found with three water molecules, and 5.5. References 121 (22) E. M. Cabaleiro-Lago, J. Rodríguez-Otero, Á. Peña-Gallego. J. Chem. Phys. 2011, 135, 214301/1-214301/9. (23) A. Rodríguez-Sanz, J. Carrazana-García, E. Cabaleiro-Lago, J. Rodríguez-Otero. J. Mol. Model. 2013, 19, 1985-1994. (24) C. Adamo, G. Berthier, R. Savinelli. Theor. Chem. Acc. 2004, 111, 176-181. (25) Y. Xu, J. Shen, W. Zhu, X. Luo, K. Chen, H. Jiang. J. Phys. Chem. B 2005, 109, 59455949. (26) A. S. Reddy, H. Zipse, G. N. Sastry. J. Phys. Chem. B 2007, 111, 11546-11553. (27) A. Campo-Cacharrón, E. Cabaleiro-Lago, J. Rodríguez-Otero. Theor. Chem. Acc. 2012, 131, 1-13. (28) T. D. Vaden, J. M. Lisy. J. Chem. Phys. 2004, 120, 721-730. (29) Y. Zhao, D. Truhlar. Theor. Chem. Acc. 2008, 120, 215-241. (30) N. J. Singh, S. K. Min, D. Y. Kim, K. S. Kim. J. Chem. Theory Comput. 2009, 5, 515529. (31) S. Tsuzuki, T. Uchimaru. Curr. Org. Chem. 2006, 10, 745-762. (32) E. G. Hohenstein, C. D. Sherrill. WIREs Comput. Mol. Sci. 2012, 2, 304-326. (33) C. D. Sherrill, Computations of Non covalent π Interactions, in Rev. Comput. Chem., Vol., John Wiley & Sons, Inc., 2009, pp.1-38. (34) S. Grimme. J. Chem. Phys. 2003, 118, 9095-9102. (35) J. G. Hill, J. A. Platts. J. Chem. Theory Comput. 2007, 3, 80-85. (36) S. F. Boys, F. Bernardi. Mol. Phys. 1970, 19, 553-566. (37) G. Chałasiński, M. M. Szczȩśniak. Chem. Rev. 2000, 100, 4227-4252. (38) B. Jeziorski, R. Moszynski, K. Szalewicz. Chem. Rev. 1994, 94, 1887-1930. (39) C. D. Sherrill. Acc. Chem. Res. 2013, 46, 1020-1028. (40) A. J. Misquitta, K. Szalewicz. J. Chem. Phys. 2005, 122, 214109/1-214109/19. (41) A. Heßelmann, G. Jansen, M. Schütz. J. Chem. Phys. 2005, 122, 014103/1-014103/17. (42) H.-J. W. and, P. J. K. and, F. R. M. and, M. S. u. t. and, P. C. and, G. K. and, T. K. and, R. L. and, A. M. and, G. R. and, T. B. A. and, R. D. A. and, A. B. and, A. B. and, D. L. C. and, M. J. O. D. and, A. J. D. and, F. E. and, E. G. and, C. H. and, A. H. and, G. H. and, T. H. and, G. J. and, C. K. o. and, Y. L. and, A. W. L. and, R. A. M. and, A. J. M. and, S. J. M. and, W. M. and, M. E. M. and, A. N. and, P. P. and, K. P. u. and, R. P. and, M. R. and, T. S. and, H. S. and, A. J. S. and, R. T. and, T. T. and, M. W. and, A. Wolf, Molpro version 2010.1, a package of ab initio programs, see http://www.molpro.net (2010). (43) 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, 122 Chap. 5. Microhydration study: complexes between guanidinium cation and aromatic species M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. M. J. A, 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, N. J. 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 B.01, Gaussian, Inc., Wallingford, CT 2009. (44) D. Kim, S. Hu, P. Tarakeshwar, K. S. Kim, J. M. Lisy. J. Phys. Chem. A 2003, 107, 1228-1238. 6 Microhydration study: complexes between pyrrolidinium cation and aromatic species 6.1. Introduction 125 6.1. Introduction Non covalent interactions are crucial in a great variety of phenomena in physics, chemistry and biochemistry.1 The interaction among molecules not bound by chemical bonds makes possible the formation of supramolecular units which are the basis of many aspects of the promising fields of supramolecular chemistry and nanotechnology.2, 3 In biological systems, non covalent interactions are responsible of phenomena as molecular recognition and play a crucial role in enzymatic catalysis and protein folding, among other processes as ion transport or molecular recognition.1, 4-6 The ubiquity of these non covalent interactions points out the need for a deep knowledge and understanding of the characteristics of the interaction at a microscopic level.7, 8 The application of computational chemistry techniques is especially suited for this purpose, since they allow the microscopic description of the system with a detail that is often not affordable from experiment. Most especially, methods based in quantum chemistry or its approximations allow the isolation of the specific effects produced by a given interaction, and a detailed dissection of its nature.9-11 Among the different kinds of non covalent contacts present in biological systems, a special place is occupied by interactions involving aromatic units.4-6, 12 When aromatic molecules participate in a non covalent interaction, it usually corresponds to one of the following three types: ··· interaction, XH··· hydrogen bond or ion··· interaction.4, 13 In the case of proteins, the presence of aromatic species is almost guaranteed, since there are several amino acids which bear an aromatic moiety in their side chains. Therefore, the side chains of phenylalanine, tyrosine and tryptophan are usually found to be interacting with participation or their aromatic groups (benzene, phenol and indole, respectively). Quite commonly, these side chain groups interact with cationic species present in the environment which allow the formation of a stabilizing cation··· interaction. In fact, the presence of other cationic groups in the side chains of amino acids as arginine of lysine makes a quite frequent event the presence of cation··· interactions within a protein.6, 14 After the seminal work of Dougherty et al.,14-18 the cation··· interaction is nowadays recognized as one of the main structural stabilizing motifs in the structure of proteins, with an important participation in a variety of recognition processes.19 Though the cation··· interaction is strong in the gas phase as corresponds to the interaction of a bare cation with a polarizable and electron rich aromatic cloud,13, 19 the presence of solvent molecules can significantly alter its characteristics.19 Several studies have shown that coordination of electron-rich solvent molecules to the cation decreases the intensity of the interaction with the cation, as expected taking into account the decrease of the effective charge carried by the cation.20-26 Also, solvent molecules 126 Chap. 6. Microhydration study: complexes between pyrrolidinium and aromatic species can promote significant structural changes as recently shown for the guanidinium···benzene interaction.23 In any case, there is some controversy about the real role played by cation··· interaction in the stability of proteins, since different works have pointed out a negligible effect whereas others indicate that the contribution to the stability is very important.27-33 These differences have been attributed to environmental effects, mostly due to the presence of the solvent.19, 31 However, it is known that the degree of exposure to the solvent of a given cation··· interaction can vary in a wide range of values, from contacts fully exposed to the solvent where the cation··· contact behaves as immersed in bulk solvent, to situations with the cation··· buried in a hydrophobic cavity with no access for solvent molecules.22 In this respect, the application of computational chemistry methods in systems which are gradually microhydrated can reveal hints about the effects caused by the presence of the most nearby solvent molecules and how their presence affects the characteristics of the cation··· interaction. Most studies to date have been carried out employing benzene as model, and mostly with alkali cations.19, 24 To our knowledge, few studies have dealt with the interaction with more complex cations as methylammonium, tetramethylammonium or guanidinium.20-23, 26, 34, 35 The presence of a more complex cation introduces some significant differences, such as a larger structural complexity resulting in a larger variety of minima, but also on the nature of the interaction, which becomes more dispersive in nature, as shown in the case of guanidinium interaction with benzene.23 In the present work, a study is performed about the interaction of a structured cation with the aromatic residues of the aromatic amino acids: benzene, phenol and indole. The cation chosen for this purpose is the pyrrolidinium cation. Pyrrolidinium is an ammonium cation within a fivemembered saturated ring. The pyrrolidine ring is present in a variety of natural products, most significantly as part of the structure of the amino acid proline. Also, it is part of the structure of nicotine, and in fact it is believed that nicotine addiction begins with the activation of the receptor by means of a cation··· contact with the pyrrolidinium ring of nicotine.16, 18 Pyrrolidinium is also frequently used as constituent in a variety of ionic liquids.36 Taking into account the ring nature of this cation it can be expected that the tendency to form stacked structures with aromatic molecules, and therefore establish contacts with a more dispersive nature, will be larger than in simpler cations. The effect of introducing several water molecules on the interaction has been analyzed in order to obtain some clues about the role water molecules play in modulating the cation··· interactions. The results obtained will help understanding the main features of these relevant interactions and the role of the closer environment upon their characteristics. 6.2. Computational details 127 6.2. Computational details Complexes formed by pyrrolidinium cation and one aromatic unit among benzene, phenol and indole, have been computationally studied by means of ab initio and density functional theory methods. The geometry of each complex has been optimized at the M06-2X/6-31+G* level of calculation.37 Once a stationary point has been reached, a frequency analysis has been carried out to ensure that the structure corresponds to a minimum in the potential energy surface of the complex. The structures found for the complexes formed by the cation and each aromatic species has been chosen as starting point for introducing water molecules in a stepwise procedure. Water molecules have been located in the most favorable regions of the cation···aromatic complex; that is, directly bonded to the N-H groups or expanding the hydrogen bond network by interacting with previous water molecules. It is worth noting that in the case of phenol and indole, the O-H and N-H groups can also participate in the hydrogen bond network, so structures with water hydrogen-bonded to phenol or indole have also been considered. Once a minimum has been found, the complexation energy has been obtained by applying the supermolecule method, using the counterpoise procedure to avoid any basis set superposition error.38, 39 Therefore, interaction energies are obtained as:  i complexcomplex ijkEijkEE ...)(...)( int (eq. 6.1) where terms in parentheses indicate the basis set and superscripts the geometry employed in the calculations. As the geometry of the molecules changes when the cluster is formed, an additional contribution describing this effect must be included, obtained as the energy difference between the molecules in the cluster geometry and in isolation.    i isolated i complex idef iEiEE )()( (eq. 6.2) defcomplex EEE  int (eq. 6.3) Since no reference values for the complexation energies are available for these systems, highlevel calculations have been carried out for non-hydrated complexes. Thus, energies have been estimated at the CCSD(T) complete basis set limit (CBS) employing an extrapolation of the MP2 correlation energy obtained with the aug-cc-pVDZ and aug-cc-pVTZ basis sets. Following a 128 Chap. 6. Microhydration study: complexes between pyrrolidinium and aromatic species widely used extrapolation procedure,10, 40 the contribution of the correlation energy to the complexation energy is obtained as: 3; )1( )1( )1( )1( 2, 33 3 2, 33 3 2,        XE XX X E XX X E ZXAV MPcorr AVXZ MPcorr CBS MPcorr . (eq. 6.4) and assuming HF values already converged to the basis limit, the MP2 complexation energy to basis limit is: CBS MPcorr AVTZ HF CBS MP EEE 2,2  . (eq. 6.5) Finally, the deficiencies on the MP2 method are corrected by using CCSD(T)/AVDZ calculations:  AVDZ MP AVDZ TCCSD CBS MP CBS TCCSD EEEE 2)(2)(  . (eq. 6.6) The results obtained with this extrapolation procedure have been employed to assess the performance of cheaper methods that will be employed in the hydrated systems. Besides the M06-2X/6-31+G* already employed in geometry optimizations, several empirically scaled variants of MP2 have been also used for obtaining complexation energies; namely SCS-MP2 and SCSN-MP2.41 In these methods the contributions to correlation energies from same-spin and different-spin electron pairs are empirically scaled in order to improve the performance of the native MP2. In SCS-MP2 the scaling factors are 1.20 and 0.33 for oppositespin and same-spin components,41 whereas in SCSN-MP2 the factors are 0.00 and 1.76.42 Finally, a new variant proposed by Hobza in order to obtain high quality results avoiding the CCSD(T) costly calculations has also been considered.43 In this MP2.X method the extrapolation to basis limit is performed as commented above, but the CCSD(T) correction is substituted by a MP3 correction adequately scaled in order to reproduce the CCSD(T) values. In the present work MP3/6-31G* level has been employed for this purpose with an scaling factor of 0.86; that is:   *316 2 *316 32.2 86.0 G MP G MP CBS MP CBS XMP EEEE   . (eq. 6.7) 6.3. Results 129 Also, in order to obtain more information about the characteristics of the interaction, the interaction energies for non-hydrated complexes have been decomposed by applying the SAPT method with intramonomer correlation effects described at the DFT level (SAPT(DFT)).44-46 So, using the geometries optimized at the M06-2X/6-31+G* level, density fitted DFT-SAPT calculations were carried out, providing information on the individual physical components of the interaction energy. For these calculations the LPBE0AC exchange-correlation was used, involving a shift parameter obtained as the sum of the ionization potential and the energy of the highest occupied molecular orbital. Orbital energies and ionization potentials have been obtained by using the PBE0 functional with the aug-cc-pVDZ basis set. The DFT-SAPT calculations were performed with the aug-cc-pVDZ basis set, employing the cc-pVTZ/JKFIT for Hartree–Fock,and aug-cc-pVDZ/MP2FIT for the second-order dispersion terms. All SAPT and CCSD(T) calculations have been performed with Molpro.47 M06-2X and MP3 calculations have been done with Gaussian09,48 whereas MP2 calculations have been performed with Turbomole 6.3.49 In order to save computational time in MP2 calculations, the resolution of the identity approach has been employed both for the HF and correlation energies. That is; RI-JK-MP2 calculations have been performed using the aug-cc-pVXZ auxiliary basis set for correlation and the def2-TZVPP auxiliary basis set for both coulomb and exchange in the calculation of HF energies.50, 51 6.3. Results 6.3.1. Pyrrolidinium···π complexes In this section the results obtained for complexes formed by pyrrolidinium and benzene, phenol or indole in the absence of water molecules will be presented. Figure 6.1 shows the minimum energy structures found for the complexes formed by pyrrolidinium and the aromatic species considered in this study, whereas Table 6.1 lists the values obtained for the complexation energies of these structures as obtained with different computational methods. Two very similar structures have been found for pyrrolidinium···benzene complex PB (the nomenclature will indicate P for pyrrolidinium, B, P or I for benzene, phenol or indole, and a numerical identifier indicating the structure). As expected, in both minima the interaction takes place by means of hydrogen bond contacts of the N-H groups of the cation with the aromatic cloud of benzene. However, whereas in PB-1 there is only one NH··· contact at about 2.0 Å, in PB-2 the pyrrolidinium cation adopts a more vertical orientation with respect to benzene, with two N-H··· contacts, though one of them is pretty long (2.8 Å). 130 Chap. 6. Microhydration study: complexes between pyrrolidinium and aromatic species In the case of phenol complexes the behavior is much more complicated since phenol presents two favorable positions for interacting with the cation, the hydroxyl oxygen and the aromatic ring. Besides, pyrrolidinium can interact with one of these positions by means of the ammonium group while simultaneously establishes favorable contacts trough the positively charged C-H groups of the carbon skeleton. Therefore, up to six different minima have been located for pyrrolidinium···phenol (PP) complexes, as shown in Figure 6.1. These six minima are grouped by pairs of similar structures which differ on whether N-H groups interact with the oxygen or the phenyl ring. Thus, in minima PP-5 and PP-6 only one contact is registered by means of a N-H···X interaction, which is oriented towards the hydroxyl oxygen in PP-5 and towards the phenyl ring in PP-6. The second N-H unit is oriented in order to interact favorably with the free OH or phenyl ring, but the distance is quite large, about 3.3-3.6 Å. In the rest of minima two simultaneous contacts with phenol are observed, one with the N-H group and the other by means of one of the C-H groups of the cation. In minima PP-1 and PP-2 the C-H group is next to nitrogen, whereas in PP-3 and PP-4 is two bonds apart. On the other hand, in PP-1 and PP-3 there are N-H···O contacts whereas in PP-2 and PP-4 there are N-H···π ones. The distances of the N-H···O hydrogen bonds are shorter than 1.9 Å, reaching values over 2 Å for N-H···π ones. The distances are usually larger for contacts involving C-H units, as expected. The behavior is similar in indole complexes. Indole also presents two favorable positions for interacting with cations corresponding to the two rings. Therefore, minima have been found showing N-H···π contacts to the phenyl or pyrrol rings of indole. Structures showing only one NH···X contact as in phenol have been tried, but they collapse to PI-5, with pyrrolidinium located roughly over the bridge between the two rings and establishing two N-H···π contacts with both rings. The distance to the phenyl ring if much shorter (2.1 Å) than to the pyrrol ring (2.7 Å). Table 6.1 summarizes the complexation energies for the complexes shown in Figure 6.1 as obtained with a variety of methods. Skipping by the moment method performance and considering only CCSD(T)/CBS values, it can be observed that both minima located for PB present complexation energies of -16 kcal mol-1 since the energy surface is very shallow regarding rotations of the N-H groups over the phenyl ring. The interaction is slightly stronger in phenol complexes, reaching -17.5 kcal mol-1 in the most stable structure corresponding to PP-3. In this structure there is a N-H···O hydrogen bond whereas a C-H group two bonds apart participates in a contact to the π cloud. Rotating the pyrrolidinium cation as to establish a NH···π plus a C-H···O contact is only penalized by 0.6 kcal mol-1 (PP-4), showing that both the hydroxyl groups and the aromatic ring are equally capable of interacting favorably with cations, as observed in previous work.25, 26 Appendix 4 233 Edef Etot Eele Eexch Eind Eexch-ind Edisp Eexch-dis δHF Phe-A 7.5 -29.3 -38.0 32.9 -22.8 10.8 -13.9 2.2 -5.8 Phe-B 6.3 -28.8 -34.8 25.1 -17.8 7.1 -10.9 1.6 -3.7 Phe-C 1.0 -28.0 -30.4 21.0 -15.0 6.4 -7.3 1.2 -3.6 Phe-D 1.9 -27.1 -30.4 20.3 -14.6 6.1 -7.0 1.1 -3.4 Phe-E 5.7 -27.0 -36.4 30.2 -21.4 10.2 -9.8 1.8 -5.7 Phe-F 5.6 -26.7 -36.0 32.2 -23.1 11.0 -10.4 1.9 -6.2 Phe-G 5.9 -26.4 -36.5 31.3 -22.2 10.6 -9.7 1.8 -6.0 Phe-H 7.3 -24.3 -34.6 29.4 -21.0 9.7 -9.6 1.7 -5.6 Phe-I 10.5 -23.1 -34.9 23.6 -17.0 7.1 -8.6 1.4 -3.7 Phe-zw-A 18.5 -27.8 -51.5 35.2 -25.1 12.0 -10.2 2.0 -7.1 Phe-zw-B 21.0 -26.8 -52.9 35.4 -25.6 12.1 -10.1 2.0 -7.1 Table A4.4. SAPT(DFT) partitioning (kcal mol-1) for complexes with phenylalanine. Edef Etot Eele Eexch Eind Eexch-ind Edisp Eexch-dis HF Tyr-A 10.3 -29.5 -40.6 34.3 -23.6 11.8 -17.0 2.7 -4.6 Tyr-B 7.8 -29.6 -38.8 30.9 -20.8 10.1 -15.3 2.4 -3.4 Tyr-C 7.1 -29.5 -35.7 24.1 -17.6 7.2 -11.0 1.6 -3.3 Tyr-D 1.4 -28.1 -30.8 21.1 -15.1 6.4 -7.4 1.2 -3.6 Tyr-E 8.3 -27.2 -36.5 30.4 -20.6 10.6 -15.6 2.4 -3.6 Tyr-F 2.2 -27.2 -30.7 20.4 -14.7 6.2 -7.0 1.2 -3.5 Tyr-G 3.5 -26.6 -29.0 20.7 -16.7 7.6 -9.5 1.4 -3.0 Tyr-H 9.2 -26.8 -35.8 28.4 -19.7 10.0 -15.6 2.4 -3.0 Tyr-O 6.1 -26.7 -36.3 32.2 -23.2 10.9 -10.4 1.9 -6.3 Tyr-Q 7.5 -25.6 -31.8 25.6 -18.5 8.6 -13.1 1.9 -3.6 Tyr-P 6.8 -25.9 -31.6 23.7 -17.8 8.5 -12.1 1.8 -3.2 Tyr-zw-A 18.9 -28.3 -52.3 35.4 -25.3 12.0 -10.3 2.0 -7.2 Tyr-zw-B 21.4 -27.2 -53.8 35.8 -25.9 12.3 -10.2 2.1 -7.2 Tyr-zw-E 25.9 -21.3 -47.9 30.6 -22.5 9.6 -12.2 2.0 -4.8 Tyr-I 6.5 -29.4 -34.7 23.4 -17.2 6.8 -10.8 1.6 -3.3 Tyr-J 10.0 -28.7 -39.3 33.6 -23.1 11.6 -16.7 2.6 -4.6 Tyr-K 1.1 -27.5 -30.1 21.2 -15.1 6.4 -7.3 1.2 -3.7 Tyr-L 8.2 -27.7 -36.0 28.3 -19.3 9.1 -14.6 2.2 -3.1 Tyr-M 2.2 -27.2 -30.7 20.4 -14.7 6.1 -7.0 1.1 -3.5 Tyr-N 6.7 -26.4 -33.0 28.1 -19.4 9.6 -14.2 2.1 -4.0 Tyr-zw-C 18.5 -28.1 -51.8 35.7 -25.5 12.1 -10.3 2.0 -7.2 Tyr-zw-D 21.4 -27.2 -53.7 35.7 -25.8 12.2 -10.2 2.1 -7.2 Tyr-zw-F 29.1 -20.4 -52.5 37.1 -26.2 13.2 -16.8 2.9 -4.5 Table A4.5. SAPT(DFT) partitioning (kcal mol-1) for complexes with tyrosine. 234 Appendix 4 Edef Etot Eele Eexch Eind Eexch-ind Edisp Eexch-dis HF Trp-A 8.6 -31.7 -39.9 32.9 -23.5 10.7 -14.6 2.3 -5.9 Trp-B 5.9 -31.6 -36.6 29.3 -21.3 10.3 -14.8 2.3 -4.3 Trp-C 6.7 -31.7 -36.7 26.6 -19.1 7.4 -12.4 1.8 -4.0 Trp-D 6.8 -31.7 -37.0 27.5 -19.6 8.0 -13.1 1.9 -4.0 Trp-E 7.6 -31.5 -38.0 28.7 -21.1 9.8 -14.5 2.2 -3.8 Trp-F 5.1 -30.1 -32.9 24.7 -18.5 8.0 -12.0 1.7 -4.2 Trp-G 10.9 -29.9 -39.6 34.4 -23.8 11.9 -18.2 2.8 -5.4 Trp-H 9.4 -28.9 -37.8 33.2 -23.5 11.8 -16.3 2.5 -5.5 Trp-I 1.3 -28.1 -31.1 21.9 -15.7 6.7 -7.4 1.2 -3.9 Trp-J 8.7 -28.0 -35.3 29.7 -20.8 9.7 -15.0 2.2 -4.6 Trp-K 9.5 -28.3 -40.8 35.1 -25.5 12.3 -12.6 2.2 -6.4 Trp-L 2.5 -27.8 -31.6 20.7 -15.0 6.3 -7.1 1.2 -3.5 Trp-M 7.0 -27.8 -38.7 31.6 -22.5 10.7 -10.1 1.8 -6.1 Trp-N 3.1 -27.5 -32.0 21.0 -15.2 6.4 -7.2 1.2 -3.6 Trp-zw-A 18.1 -29.9 -53.5 37.0 -26.5 12.7 -10.5 2.1 -7.6 Trp-zw-B 21.3 -29.1 -55.9 36.8 -26.6 12.7 -10.4 2.1 -7.5 Trp-zw-C 21.1 -28.7 -55.1 36.5 -26.4 12.5 -10.4 2.1 -7.4 Trp-zw-D 29.0 -24.8 -55.7 39.7 -28.8 14.5 -16.8 2.9 -6.9 Table A4.6. SAPT(DFT) partitioning (kcal mol-1) for complexes with tryptophan. Appendix 5 235 Appendix 5 “Cation···π interactions between imidazolium and aromatic amino acids” Table A5.1. Complexation energies (kcal mol-1) obtained for the most stable minima of the imidazolium···phenylalanine complex. M06-2X/ 6-31+G* MP2/CBS MP2.X SAPT Edef,MP2/CBS Phe-A -30.50 -31.90 -29.24 -29.48 8.23 Phe-B -27.07 -28.38 -28.31 -26.47 20.07 Phe-C -29.51 -31.14 -28.29 -28.57 8.29 Phe-D -27.54 -28.48 -28.13 -27.92 1.68 Phe-E -26.54 -27.47 -27.45 -25.68 22.55 Phe-F -25.94 -27.76 -27.37 -25.83 20.74 Phe-G -26.78 -27.61 -27.25 -27.21 2.40 Phe-H -26.41 -27.42 -26.98 -27.06 2.56 Phe-I -27.44 -27.96 -26.68 -27.29 7.05 Phe-J -25.50 -26.88 -26.55 -25.10 23.53 Phe-K -27.41 -27.45 -26.54 -26.90 7.40 Phe-L -26.56 -27.64 -26.44 -26.22 5.49 Phe-M -26.83 -27.22 -26.43 -26.31 5.79 Phe-N -27.08 -26.88 -25.85 -26.28 7.78 Phe-O -24.82 -26.14 -25.83 -23.81 20.06 Phe-P -25.99 -26.70 -25.67 -25.63 6.49 Phe-Q -25.30 -25.11 -25.22 -22.96 23.22 Phe-R -25.92 -26.42 -25.17 -25.65 7.19 Phe-S -24.04 -25.04 -24.76 -22.72 22.25 Phe-T -24.72 -25.32 -24.45 -24.08 3.85 236 Appendix 5 Table A5.2. Complexation energies (kcal mol-1) obtained for the most stable minima of the imidazolium···tyrosine complex. M06-2X/ 6-31+G* MP2/CBS MP2.X SAPT Edef,MP2/CBS Tyr-A -31.38 -32.90 -29.77 -29.95 9.31 Tyr-B -31.28 -32.68 -29.50 -29.52 9.78 Tyr-C -30.23 -30.64 -29.38 -29.04 7.02 Tyr-D -29.85 -30.14 -28.94 -28.58 6.99 Tyr-E -30.66 -32.26 -28.77 -28.92 9.98 Tyr-F -27.26 -28.64 -28.46 -26.90 20.54 Tyr-G -27.69 -28.73 -28.27 -27.93 2.10 Tyr-H -27.10 -28.46 -28.27 -26.71 20.20 Tyr-I -29.32 -30.24 -28.24 -27.84 7.11 Tyr-J -27.21 -28.42 -27.86 -27.64 2.25 Tyr-K -27.15 -28.18 -27.72 -27.38 1.81 Tyr-L -26.81 -27.71 -27.68 -25.70 19.81 Tyr-M -26.57 -27.64 -27.52 -26.08 23.05 Tyr-N -27.94 -28.46 -27.48 -27.49 5.72 Tyr-O -26.51 -27.56 -27.45 -26.00 23.04 Tyr-P -25.97 -27.80 -27.29 -26.06 21.07 Tyr-Q -28.40 -28.91 -27.19 -26.77 8.00 Tyr-R -26.67 -27.62 -27.18 -27.24 2.78 Tyr-S -28.29 -29.24 -27.17 -26.79 7.35 Tyr-T -28.47 -28.70 -27.07 -26.62 8.62 Appendix 5 237 Table A5.3. Complexation energies (kcal mol-1) obtained for the most stable minima of the imidazolium···tryptophan complex. M06-2X/ 6-31+G* MP2/CBS MP2.X SAPT Edef,MP2/CBS Trp-A -33.03 -35.28 -31.42 -31.60 10.65 Trp-B -31.36 -33.03 -31.27 -30.98 5.66 Trp-C -29.02 -30.44 -30.30 -28.49 20.18 Trp-D -29.61 -31.71 -30.27 -30.15 7.28 Trp-E -31.53 -33.57 -30.10 -30.01 9.26 Trp-F -29.45 -31.68 -29.82 -29.62 5.74 Trp-G -28.74 -29.90 -29.60 -27.92 20.98 Trp-H -29.37 -31.95 -29.48 -29.28 5.92 Trp-I -28.29 -29.49 -29.36 -27.86 23.59 Trp-J -27.80 -29.72 -29.25 -27.80 21.26 Trp-K -30.51 -33.21 -28.95 -29.05 10.91 Trp-L -30.03 -32.39 -28.60 -28.85 11.95 Trp-M -29.18 -30.30 -28.53 -28.99 9.57 Trp-N -27.87 -28.85 -28.47 -27.94 2.03 Trp-O -26.55 -29.11 -28.47 -27.20 22.06 Trp-P -27.20 -28.82 -28.39 -27.20 24.64 Trp-Q -28.01 -30.16 -28.25 -27.92 5.31 Trp-R -27.46 -28.13 -28.15 -26.10 22.04 Trp-S -27.01 -28.54 -28.08 -26.91 24.69 Trp-T -28.06 -30.04 -28.06 -27.75 8.37 238 Appendix 5 Table A5.4a. Complexation energies (kcal mol-1) obtained for the most stable minima of the imidazolium···histidine complex. M06-2X/ 6-31+G* MP2/CBS MP2.X SAPT Edef,MP2/CBS His-A -34.35 -34.90 -34.78 -33.64 19.36 His-B -34.52 -34.46 -34.32 -33.22 21.18 His-C -33.85 -34.06 -34.04 -32.49 18.48 His-D -32.34 -32.87 -33.99 -44.46 79.67 His-E -33.02 -34.07 -33.60 -32.88 21.11 His-F -33.99 -33.56 -33.53 -32.02 20.00 His-G -33.30 -33.65 -33.19 -32.47 22.92 His-H -31.34 -31.72 -32.88 -43.42 82.80 His-I -31.54 -31.78 -32.68 -43.06 78.81 His-J -32.05 -32.80 -32.38 -31.09 19.52 His-K -31.07 -30.20 -32.16 -36.87 89.81 His-L -30.19 -29.26 -31.14 -35.12 87.34 His-M -29.20 -29.62 -30.70 -41.06 80.89 His-N -28.49 -29.62 -30.57 -40.87 74.55 His-O -29.71 -28.36 -30.39 -34.48 95.91 His-P -30.52 -30.57 -30.21 -30.18 6.05 His-Q -28.62 -28.44 -30.10 -34.26 79.17 His-R -30.69 -31.21 -30.00 -30.14 11.01 Appendix 5 239 Table A5.4b. Complexation energies (kcal mol-1) obtained for the most stable minima of the imidazolium···histidine complex. M06-2X/ 6-31+G* MP2/CBS MP2.X SAPT Edef,MP2/CBS His-S -31.22 -30.82 -29.93 -30.38 8.45 His-T -30.77 -30.45 -29.89 -30.03 7.12 His-U -30.07 -30.32 -29.83 -29.99 6.38 His-V -27.04 -26.69 -29.81 -43.20 122.61 His-W -30.52 -30.82 -29.70 -29.95 8.41 His-X -30.45 -29.99 -29.62 -29.85 6.19 His-Y -27.91 -27.59 -29.08 -33.28 77.58 His-Z -28.95 -29.32 -28.94 -29.04 9.80 His-AA -29.56 -29.81 -28.89 -29.31 7.33 His-AB -28.11 -26.63 -28.71 -33.32 100.73 His-AC -26.75 -27.64 -28.49 -40.57 88.08 His-AD -29.73 -29.41 -28.47 -28.54 14.47 His-AE -28.27 -28.44 -28.26 -27.54 37.41 His-AF -29.41 -28.32 -28.20 -26.54 5.26 His-AG -28.77 -27.92 -27.71 -25.54 6.95 His-AH -25.78 -26.96 -27.43 -24.54 75.71 His-AI -25.66 -26.97 -27.42 -23.54 80.72 240 Appendix 5 Figure A5.1. Structures of the most stable minima of Imz···Phe complex as obtained at the M06-2X/6-31+G* level. A B C D E F G H I J K L M N O P Q R S T Appendix 5 241 Figure A5.2. Structures of the most stable minima of Imz···Tyr complex as obtained at the M06-2X/6-31+G* level. A B C D E F G H I J K L M N O P Q R S T 242 Appendix 5 Figure A5.3. Structures of the most stable minima of Imz···Trp complex as obtained at the M062X/6-31+G* level. A B CD E F GH I J K L M N O P QRST Q Appendix 5 249 Table A5.5. SAPT(DFT) contributions in mEh to the interaction energy in Imz···Phe complexes as obtained employing PBE0 with the aug-cc-pVDZ basis set. Eele Erep Eind Eexch-ind Edis Eexch-dis HF Phe-A -59.84 54.14 -38.58 20.10 -25.04 3.97 -10.78 Phe-B -81.96 60.91 -45.64 22.80 -16.50 3.29 -14.52 Phe-C -58.59 53.17 -38.02 19.91 -24.72 3.91 -10.39 Phe-D -50.13 38.09 -28.64 13.32 -11.86 2.01 -8.04 Phe-E -84.56 61.66 -46.80 23.21 -16.38 3.31 -14.78 Phe-F -81.52 63.80 -47.47 23.20 -16.61 3.29 -16.33 Phe-G -50.22 37.42 -28.27 13.09 -11.44 1.96 -7.89 Phe-H -50.08 37.83 -28.54 13.11 -11.27 1.93 -8.37 Phe-I -59.46 53.71 -40.19 19.78 -17.02 3.02 -11.85 Phe-J -84.77 65.65 -49.43 23.95 -16.58 3.33 -17.08 Phe-K -59.12 48.94 -36.62 17.76 -15.41 2.76 -10.52 Phe-L -48.77 39.54 -29.64 14.26 -18.46 2.74 -7.15 Phe-M -47.91 33.80 -25.56 10.98 -17.30 2.46 -4.77 Phe-N -59.60 51.50 -38.21 18.65 -15.66 2.84 -11.32 Phe-O -75.86 55.07 -41.54 20.51 -15.69 3.05 -13.02 Phe-P -47.12 34.70 -26.95 11.79 -16.93 2.37 -6.23 Phe-Q -80.21 54.38 -41.92 21.06 -15.64 3.10 -11.95 Phe-R -55.89 46.07 -34.70 16.79 -14.88 2.64 -10.00 Phe-S -77.41 54.64 -41.81 20.59 -15.53 3.05 -12.79 Phe-T -41.85 29.86 -23.29 10.83 -15.23 2.13 -4.43 250 Appendix 5 Table A5.6. SAPT(DFT) contributions in mEh to the interaction energy in Imz···Tyr complexes as obtained employing PBE0 with the aug-cc-pVDZ basis set. Eele Erep Eind Eexch-ind Edis Eexch-dis HF Tyr-A -62.22 56.64 -39.64 21.22 -27.64 4.36 -10.81 Tyr-B -63.62 57.91 -40.89 22.49 -28.29 4.55 -10.19 Tyr-C -56.34 42.03 -29.87 13.43 -21.02 3.04 -5.27 Tyr-D -54.70 41.28 -28.96 12.36 -20.83 2.93 -5.31 Tyr-E -63.46 58.31 -40.64 22.65 -28.80 4.63 -10.01 Tyr-F -83.29 61.43 -46.08 22.87 -16.56 3.30 -14.71 Tyr-G -50.68 38.12 -28.65 13.25 -11.92 2.01 -8.07 Tyr-H -82.68 62.28 -46.63 23.18 -16.67 3.33 -15.00 Tyr-I -55.99 46.04 -33.38 17.89 -23.59 3.58 -6.38 Tyr-J -50.19 38.41 -28.91 13.25 -11.69 1.96 -8.58 Tyr-K -49.49 38.39 -28.79 13.37 -11.96 2.02 -8.14 Tyr-L -78.81 53.84 -41.12 20.57 -15.72 3.07 -11.91 Tyr-M -86.08 62.52 -47.40 23.44 -16.51 3.34 -15.06 Tyr-N -49.09 31.64 -23.76 9.72 -16.90 2.25 -3.99 Tyr-O -85.92 62.44 -47.34 23.40 -16.50 3.34 -15.04 Tyr-P -82.48 65.58 -48.82 23.69 -16.78 3.34 -17.04 Tyr-Q -54.03 43.72 -31.31 15.44 -22.71 3.33 -6.11 Tyr-R -50.84 37.72 -28.49 13.15 -11.53 1.97 -7.99 Tyr-S -53.76 45.48 -32.40 16.99 -23.72 3.53 -6.62 Tyr-T -52.81 40.76 -30.40 14.51 -22.31 3.16 -5.37 Appendix 5 251 Table A5.7. SAPT(DFT) contributions in mEh to the interaction energy in Imz···Trp complexes as obtained employing PBE0 with the aug-cc-pVDZ basis set. Eele Erep Eind Eexch-ind Edis Eexch-dis HF Trp-A -65.96 60.99 -42.83 22.69 -30.50 4.77 -11.53 Trp-B -55.41 43.61 -33.60 16.49 -21.22 3.22 -7.99 Trp-C -85.98 65.88 -49.50 24.63 -17.23 3.50 -16.21 Trp-D -56.26 42.20 -32.95 15.10 -19.48 2.92 -7.99 Trp-E -58.57 50.86 -37.45 20.12 -28.88 4.39 -8.32 Trp-F -53.89 42.89 -31.90 15.23 -19.57 2.93 -8.85 Trp-G -48.75 35.62 -27.58 14.10 -20.93 3.03 -4.99 Trp-H -54.37 45.32 -33.57 17.50 -22.70 3.47 -8.03 Trp-I -90.41 66.76 -50.48 24.99 -17.10 3.50 -16.63 Trp-J -86.14 70.15 -52.26 25.37 -17.51 3.54 -18.62 Trp-K -62.57 59.84 -41.42 22.19 -30.44 4.76 -11.09 Trp-L -66.39 62.52 -44.15 24.02 -29.96 4.77 -10.98 Trp-M -64.96 59.39 -43.98 21.41 -20.46 3.45 -13.02 Trp-N -50.88 39.45 -29.77 13.86 -12.15 2.08 -8.45 Trp-O -85.15 68.24 -51.15 24.69 -17.88 3.47 -17.95 Trp-P -90.59 70.81 -53.09 25.67 -17.34 3.54 -18.95 Trp-Q -48.75 35.62 -27.58 14.10 -20.93 3.03 -4.99 Trp-R -83.52 55.94 -43.18 21.73 -15.97 3.20 -12.43 Trp-S -90.13 71.04 -53.27 25.68 -17.35 3.53 -19.06 Trp-T -52.94 43.29 -32.57 15.50 -22.27 3.18 -8.04 252 Appendix 5 Table A5.8. SAPT(DFT) contributions in mEh to the interaction energy in Imz···His complexes as obtained employing PBE0 with the aug-cc-pVDZ basis set. Eele Erep Eind Eexch-ind Edis Eexch-dis HF His-A -94.54 70.49 -52.35 26.51 -17.89 3.74 -17.68 His-B -96.76 71.01 -52.98 26.70 -17.77 3.74 -17.92 His-C -90.39 63.42 -47.60 24.42 -17.10 3.53 -14.88 His-D -109.76 170.50 -171.33 49.54 -30.58 4.14 -105.23 His-E -95.60 77.40 -57.00 27.94 -18.50 3.85 -21.32 His-F -91.79 62.35 -47.33 24.28 -16.88 3.50 -14.44 His-G -97.86 77.60 -57.34 28.01 -18.35 3.84 -21.37 His-H -112.07 170.47 -171.78 49.51 -30.60 4.15 -105.73 His-I -108.00 170.22 -170.23 49.42 -30.46 4.12 -104.18 His-J -89.23 66.78 -49.58 24.79 -17.31 3.55 -16.99 His-K -148.93 176.99 -169.25 63.47 -32.33 6.53 -93.37 His-L -147.24 174.92 -164.47 62.44 -31.82 6.39 -90.46 His-M -109.32 169.96 -168.75 48.99 -30.39 4.14 -103.90 His-N -103.94 169.41 -165.84 48.88 -30.20 4.10 -101.29 His-O -156.38 176.03 -167.64 62.88 -32.01 6.49 -92.25 His-P -62.66 46.03 -34.38 16.54 -13.24 2.42 -10.33 His-Q -139.31 170.00 -156.05 59.49 -30.41 6.00 -85.78 His-R -70.09 55.48 -40.37 19.76 -18.92 3.24 -11.64 His-S -66.02 47.35 -35.53 17.32 -15.75 2.82 -9.58 His-T -61.29 40.42 -30.35 14.33 -15.10 2.49 -7.26 His-U -62.70 46.89 -34.91 16.64 -13.18 2.41 -11.03 His-V -130.42 191.47 -218.49 66.19 -41.17 5.70 -130.67 His-W -60.92 39.60 -30.31 13.27 -15.12 2.44 -7.66 His-X -61.70 43.45 -32.30 15.35 -12.83 2.28 -9.65 His-Y -135.25 168.14 -153.80 58.23 -29.93 5.84 -85.27 His-Z -67.20 47.31 -35.25 16.99 -13.28 2.46 -10.84 His-AA -60.69 42.75 -31.40 14.77 -15.11 2.48 -8.75 His-AB -152.48 174.68 -173.13 63.39 -31.78 6.29 -95.65 His-AC -111.02 174.52 -176.11 51.93 -35.23 4.78 -108.01 His-AD -71.13 50.91 -37.23 16.55 -16.88 2.86 -10.91 His-AE -113.49 106.76 -79.49 37.06 -22.30 4.73 -34.13 His-AF -54.73 31.80 -25.31 12.20 -11.73 2.01 -5.19 His-AG -58.13 34.80 -25.34 11.47 -12.97 2.16 -5.61 His-AH -113.56 172.19 -157.49 48.64 -34.19 4.63 -94.28 His-AI -114.69 174.22 -163.96 50.43 -34.92 4.74 -98.80