Significado
Las reacciones ácido-base se encuentran entre los procesos químicos más importantes. Sin embargo, carecemos de una forma sencilla de describir esta clase de reacciones en función de las coordenadas atómicas. De hecho, una vez disueltos en agua, el H+ y su anión conjugado OH- tienen una estructura altamente fluxional difícil de precisar. Aquí resolvemos esta cuestión adoptando el punto de vista de la descripción de las reacciones ácido-base como un equilibrio entre el soluto y el conjunto del disolvente. Esto permite identificar descriptores de aplicación general. Como consecuencia, ahora es posible realizar una simulación cuantitativa de muestreo mejorado de la reacción ácido-base en el agua y en otros entornos como las cavidades de las zeolitas o en las superficies.
Abstract
Las reacciones ácido-base son omnipresentes en la naturaleza. Entender sus mecanismos es crucial en muchos campos, desde la bioquímica hasta la catálisis industrial. Desgraciadamente, los experimentos sólo proporcionan una información limitada sin que se conozca mucho el comportamiento molecular. Las simulaciones atomísticas podrían complementar los experimentos y arrojar una preciosa luz sobre los mecanismos microscópicos. Sin embargo, las grandes barreras de energía libre relacionadas con la disociación de protones hacen obligatorio el uso de métodos de muestreo mejorados. Aquí realizamos una simulación de dinámica molecular (DM) ab initio y mejoramos el muestreo con la ayuda de la metadinámica. Esto ha sido posible gracias a la introducción de descriptores o variables colectivas (CVs) que se basan en una perspectiva conceptualmente diferente de los equilibrios ácido-base. Probamos con éxito nuestro enfoque en tres soluciones acuosas diferentes de ácido acético, amoníaco y bicarbonato. Estas son representativas del comportamiento ácido, básico y anfotérico.
- ácido-base
- metadinámica
- variables colectivas
- muestreo mejorado
Las reacciones ácido-base juegan un papel clave en muchas ramas de la química. Las reacciones de complejación inorgánica, el plegado de proteínas, los procesos enzimáticos, la polimerización, las reacciones catalíticas y muchas otras transformaciones en diferentes áreas son sensibles a los cambios de pH. Comprender el papel del pH en estas reacciones implica tener control sobre su reactividad y cinética.
La importancia crucial del pH ha estimulado la recopilación de una gran cantidad de datos sobre los equilibrios ácido-base. Estos se miden normalmente en fases gaseosas y condensadas, utilizando técnicas espectroscópicas y potenciométricas. Sin embargo, la precisión de estos métodos tiene limitaciones prácticas, especialmente en las fases condensadas (1). Además, es muy difícil extraer de los datos experimentales una imagen microscópica de los procesos implicados. Por lo tanto, no es sorprendente que el equilibrio ácido-base haya sido objeto de una intensa actividad teórica (1⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓-12).
La acidez de una especie química en el agua puede expresarse en términos de pKa, el logaritmo negativo de la constante de disociación del ácido. Hay dos formas de calcular estos valores, una estática y otra dinámica.
El enfoque más estándar es el estático, en el que las energías libres en fase de solución, y en consecuencia los pKas, se obtienen cerrando un ciclo de Born-Haber compuesto por las energías libres en fase gaseosa y de solvatación (1, 3⇓⇓⇓-7). Aunque es extremadamente exitoso en muchos casos, el enfoque estático tiene algunas limitaciones. Es necesario elegir un modelo de solvatación y los modelos de solvente continuo tienen una precisión limitada. Esto es particularmente cierto en sistemas como las zeolitas o las proteínas caracterizadas por cavidades irregulares en las que una descripción implícita del disolvente es un reto. Obviamente, con este enfoque no se puede obtener información dinámica. Además, puede haber reacciones competitivas que no pueden tenerse en cuenta a menos que se incluyan explícitamente en el modelo.
En principio, estas limitaciones podrían eliminarse en un enfoque más dinámico basado en simulaciones de dinámica molecular (MD) en las que las moléculas del disolvente se traten explícitamente. Si se dispusiera de un tiempo de computación ilimitado, tales simulaciones explorarían todas las vías posibles y asignarían el peso estadístico relativo a los diferentes estados. Por desgracia, la presencia de cuellos de botella cinéticos frustra esta posibilidad al atrapar el sistema en estados metaestables, ya que los diferentes estados de protonación están separados por grandes barreras. Además, en las reacciones ácido-base se rompen y forman enlaces químicos. Esto requiere el uso de MD ab initio en el que las fuerzas interatómicas se calculan sobre la marcha a partir de teorías de estructura electrónica. Esto encarece el cálculo y reduce aún más la escala de tiempo que puede explorarse.
Para superar esta dificultad, se hace obligatorio el uso de métodos de muestreo mejorado (13) que aceleran la exploración del espacio configuracional. Una clase muy popular de métodos de muestreo mejorado se basa en la identificación de los grados de libertad que están involucrados en la reacción lenta de interés. Estos grados de libertad suelen denominarse variables colectivas (CV) y se expresan como funciones explícitas de las coordenadas atómicas R. El muestreo se mejora entonces añadiendo un sesgo que es función de las CV elegidas (14⇓-16). Además, el diseño de un conjunto adecuado de buenas CVs tiene también un significado más profundo. Los CVs exitosos capturan de forma condensada la física del problema, identifican sus grados de libertad lentos y conducen a una descripción modelística útil del proceso.
En las reacciones químicas estándar, esto es relativamente sencillo ya que se pueden asignar estructuras bien definidas a los reactivos y productos (17⇓-19). Este no es el caso de las reacciones ácido-base en las que se añade o sustrae un protón del soluto. Una vez que este proceso ha tenido lugar, los iones de agua (H+ o/y OH-) se disuelven y su estructura se vuelve esquiva. De hecho, los iones de agua pueden difundirse rápidamente en el medio a través de un mecanismo de Grotthuss (20). Se vuelven altamente fluxionales y la identidad de los átomos que participan en su estructura cambia continuamente. La naturaleza de estas especies es, por tanto, difícil de capturar en una función analítica explícita de R. Sin embargo, dada la relevancia de las reacciones ácido-base, se han hecho muchos intentos de definir estas entidades (8⇓⇓⇓-12). Desgraciadamente, estos CVs tienen una naturaleza ad hoc y, aunque tienen éxito en tal o cual caso, no pueden aplicarse de forma general.
Para construir CVs generales y útiles damos dos pasos conceptuales. Uno es considerar el proceso ácido-base como una reacción en la que intervienen sólo unos pocos elementos, a saber, todo el disolvente y los residuos que reaccionan en la molécula disuelta. Por ejemplo, cuando sólo hay un tipo de residuo que se disocia, pensamos en el equilibrio ácido-base como una reacción del tipoA+H2NON⇌Bq0+H2N+q1ONq1,donde N es el número de moléculas de agua, A y B son una molécula genérica ácido-base en disolución y su especie conjugada, respectivamente, q0 y q1 son enteros que pueden asumir valores +1 y -1 según el comportamiento ácido-base de la especie, y q1+q0=0.
Esto implica que no consideramos el disolvente como un conjunto de moléculas que compiten por reaccionar con la especie ácido-base. Más bien consideramos el disolvente en su totalidad como uno de los dos aductos. Adoptar este punto de vista es especialmente relevante en disolventes polares como el agua, que se caracterizan por tener redes muy estructuradas. En este caso la presencia de un exceso o una deficiencia de protones cambia localmente la estructura de la red y esta distorsión se propaga a lo largo de toda la red.
Desde los primeros días de Wicke y Eigen (21) y Zundel y Metzger (22), los investigadores han luchado con el número de moléculas que deben incluirse en la definición de la perturbación (23⇓-25). Dada la ausencia de parámetros físicos capaces de dar una respuesta clara e inequívoca a esta cuestión, la idea de considerar el disolvente como un todo elude este problema. Así, el disolvente no es sólo un medio con un papel pasivo, sino que se considera como un conjunto de moléculas que contribuyen colectivamente a la formación del par ácido-base conjugado. Este punto de vista es mucho más cercano al original propuesto por Brønsted y Lowry en el que la reacción puede verse como un simple intercambio de un catión hidrógeno entre un par ácido-base.
Para que la reacción tenga lugar el centro de la perturbación tiene que alejarse del soluto. Por lo tanto, el segundo paso importante es controlar el centro de la perturbación. Debido a los mecanismos de Grotthuss, la perturbación se mueve a lo largo de la red. Esto puede llevar a diferentes definiciones del centro del defecto. Sin embargo, si teselamos todo el espacio utilizando poliedros de Voronoi centrados en los átomos de oxígeno del agua, podemos asignar inequívocamente cada átomo de hidrógeno a uno y sólo uno de estos poliedros. El sitio cuyo poliedro de Voronoi contiene un número anómalo de protones se toma como centro de la perturbación (Fig. 1).
Dos ejemplos de partición del espacio. (Izquierda) Mostramos una aproximación convectiva en la que se utiliza la distancia del átomo de oxígeno para definir su entorno. Pueden verse claramente superposiciones artificiales. (Derecha) La teselación de Voronoi no sufre estos defectos.
Este punto de vista confiere al método un carácter muy general, haciéndolo aplicable a todo sistema ácido-base, sin necesidad de fijar de antemano los pares reaccionantes. Así, es posible explorar todos los estados de protonación relevantes incluso en sistemas compuestos por más de un par ácido-base.
Este enfoque general permite definir CVs sin tener que imponer estructuras específicas o seleccionar la identidad de los átomos involucrados. Probamos nuestro método realizando simulaciones metadinámicas en un caso de ácido débil (ácido acético), en una base débil (amoníaco) y en una especie anfótera (bicarbonato) elegidas como puntos de referencia debido a su fuerza comparable, pero a su diferente comportamiento ácido-base.
Métodos
Como se ha comentado anteriormente, introducimos dos CVs, uno relacionado con el estado de protonación y otro que localiza los defectos de carga y mide su distancia relativa. Ambos CVs necesitan una definición robusta para asignar los átomos de hidrógeno al respectivo sitio ácido-base. Para conseguir este resultado, dividimos todo el espacio en poliedros de Voronoi centrados en los sitios ácido-base i situados en Ri. Los sitios incluyen todos los átomos capaces de romper y formar enlaces con un protón ácido. La partición estándar del espacio de Voronoi se describe mediante un conjunto de funciones de índice wi(r) centradas en los diferentes Ris, de forma que wi(r)=1 si el átomo i es el más cercano a r y wi(r) = 0 en caso contrario. Para su uso en métodos de muestreo mejorados, los CV deben ser diferenciables. Para ello, introducimos una versión suave de las funciones de índice, wis(r). Éstas se definen utilizando funciones softmaxwis(r)=e-λ|Ri-r|∑me-λ|Rm-r|, donde i y m recorren todos los sitios ácido-base y λ controla la inclinación con la que las curvas decaen a 0, es decir, la selectividad de la función. Con una elección adecuada de λ, esta definición consigue el resultado deseado, como se muestra en la Fig. 2. De este modo, un átomo de hidrógeno en una posición Rj se asigna al poliedro centrado en el sitio i con el peso wi(Rj). Entonces, el número total de átomos de hidrógeno asignados al sitio ácido-base i esWi=∑j∈Hwis(Rj),donde la suma en j corre sobre todos los átomos de hidrógeno.
Teselación suave de un espacio 2D con celdas centradas en los tres átomos de oxígeno de la molécula de agua. Las regiones azules planas representan la porción de espacio en la que la función asume un valor de 1 y las amarillas representan las fronteras entre celdas. Esta superficie se ha obtenido con un valor de λ=4.
Se puede asociar a cada sitio ácido-base un valor de referencia Wi0 que cuenta el número de átomos de hidrógeno enlazados en el estado neutro. La diferencia entre el valor instantáneo de los átomos de hidrógeno y el de referencia esδi=Wi-Wi0.Cuando es diferente de cero, δi señalará si el ith sitio ha ganado o perdido un protón. En el caso de los átomos de oxígeno del agua, un ion hidronio tiene un δi=+1 mientras que un ion hidróxido tiene δi=-1.
A continuación, agrupamos los sitios ácido-base en especies. Por ejemplo, en el caso del aminoácido más simple, la glicina, en solución acuosa, el número de especies Ns será igual a 3. Todos los átomos de oxígeno del agua pertenecen a una especie; luego se cuentan en otra especie los dos átomos de oxígeno carboxílicos, y finalmente se considera como tercera especie el átomo de nitrógeno del grupo amino.
En el espíritu de este trabajo contamos el exceso o los defectos totales de protones asociados a cada especie:qk=∑i∈kδi.Esto implica que no nos interesa la identidad específica del sitio reaccionado, sino si la kª especie en su totalidad ha aumentado (qk=+1), ha disminuido (qk=-1), o no ha cambiado su número de protones. Si consideramos un soluto con una sola fracción reactiva, entonces cada estado posible del sistema puede ser descrito por uno de los tres vectores 2D (0,0), (-1,1), o (1,-1).
En el caso general cada estado de protonación puede ser descrito por un vector q→=(q0,q1,…qNs-1) con dimensión igual al número de sitios reactivos no equivalentes, Ns. En el Apéndice del SI se ofrece una explicación más exhaustiva.
Para su uso en el muestreo mejorado, estos vectores deben expresarse como una función escalar f=f(q→) tal que, para cada q→ físicamente relevante, f alcance valores capaces de distinguir los diferentes estados de protonación globales. Hay infinitas formas de construir un escalar a partir de un vector. Posiblemente la opción más sencilla sea escribir f(q→)=X→⋅q→ y, para distinguir los diferentes estados de protonación, elegir X→=(20,21,22,…2Ns-1).
Esto conduce a la siguiente definición para el CV que se utiliza para describir el estado de protonación del sistema,sp=∑k=0Ns-12k⋅qk,donde k son los índices utilizados para etiquetar los respectivos grupos de sitios reactivos. En el Apéndice del SI se ha elaborado un ejemplo en detalle. Por supuesto, la CV se hace continua por el uso de Wi en el cálculo de la δi necesaria para evaluar qk en la Ec. 5.
La segunda CV es una suma de las distancias entre todos los sitios ácido-base multiplicada por su carga parcial δi,sd=∑i,m>i-rim⋅δi⋅δm,donde los índices i y m recorren todos los sitios ácido-base pertenecientes a diferentes grupos k, y rim es la distancia entre los dos átomos. De este modo, sólo el par ácido-base que ha intercambiado un protón da una contribución diferente de cero. La ecuación 7 sólo es válida cuando hay un solo par ácido-base conjugado. Sin embargo, debido a la acción del sesgo, puede ocurrir ocasionalmente que se formen varios pares ácido-base. Para evitar el muestreo de estos eventos tan improbables, aplicamos una restricción en el número de pares. En el Apéndice del SI se ofrecen más detalles.
Resultados
Hemos aplicado nuestro método a tres soluciones acuosas de ácido acético, amoníaco y bicarbonato como representaciones de un ácido débil, una base débil y un compuesto anfótero, respectivamente. Las configuraciones de las tres simulaciones son idénticas, excepto por la identidad de las moléculas disueltas. Esto asegura que el resultado refleja la diferente química de estos tres sistemas y que no hay ningún sesgo debido a la condición inicial.
Cada simulación de los sistemas se realizó con simulaciones MD de Born-Oppenheimer combinadas con metadinámica bien temperada (14, 26) utilizando el paquete CP2K (27) parcheado con PLUMED 2 (28) y funcional fuertemente restringido y apropiadamente normado (29) para la energía xc, Exc. Véase el Apéndice SI para más detalles.
En la Fig. 3 representamos las superficies de energía libre (FESs) en función de sp y sd. Estas FESs reproducen vívidamente el comportamiento esperado. Todas tienen un mínimo en sp=0 que corresponde al estado en el que no hay cargas en el disolvente. En el FES del ácido acético (Fig. 3A) un segundo mínimo cerca de sp=-1 refleja su comportamiento ácido. Por el contrario, el FES del amoníaco (Fig. 3B) muestra un segundo mínimo cercano a sp=1. Las formas de los FESs de amoníaco y ácido acético están aproximadamente relacionadas por una simetría de espejo que refleja su comportamiento contrastado. De forma similar, el FES simétrico del bicarbonato (Fig. 3C) refleja su carácter anfótero.
(A-C) Superficies de energía libre a lo largo de sp y sd de ácido acético (A), amoníaco (B) y bicarbonato (C) en solución acuosa. Las barras de color indican la energía libre expresada en unidades kJ⋅mol-1. El CV sd se expresa en angstroms.
A medida que se forma el par conjugado sd comienza a asumir valores positivos correspondientes a la separación y difusión del par conjugado. En comparación con el estado no disociado en el que sólo se permite sd=0, los estados en los que está presente un par conjugado muestran una forma alargada de las cuencas a lo largo de esta variable. Esto se debe al comportamiento difusivo de los iones hidronio e hidróxido en solución que hace accesible un rango continuo de distancias. Además, a lo largo de este CV se observa una barrera alrededor de 1,5 que corresponde a la ruptura del enlace covalente entre el átomo de hidrógeno y el sitio ácido-base.
Conclusiones
La aplicabilidad general de este método a sistemas de distinta naturaleza es un paso importante dado en su comprensión y descripción. El esquema puede ampliarse para incluir efectos nucleares cuánticos con el uso de la dinámica molecular integral de trayectoria (30). Esto tendría una importancia cuantitativa ya que, por ejemplo, los valores de pKa se ven afectados por la deuteración. Además, la ausencia de suposiciones o imposiciones sobre los candidatos reactivos o las trayectorias de reacción permite ampliar este método a sistemas de complejidad creciente que no pueden abordarse con los métodos tradicionales. Ejemplos de preguntas a las que ahora se puede dar respuesta son los equilibrios tautoméricos en procesos bioquímicos y el comportamiento de los ácidos en las zeolitas y en la superficie de los óxidos expuestos al agua.
Agradecimientos
Los cálculos se llevaron a cabo en el clúster ETH Euler y en el clúster Mönch del Centro Nacional de Supercomputación de Suiza. Esta investigación fue apoyada por la subvención de la Unión Europea ERC-2014-AdG-670227/VARMET.
Publicado bajo la licencia PNAS.