Introducción a la química computacional

Introducción a la química computacional
Participantes en el Congreso Solvay de 1927. De arriba a abajo y de izquierda a derecha: Piccard, Henriot, Ehrenfest, Herzen, De Donder, Schrödinger, Verschaffelt, Pauli, Heisenberg, Fowler, Brillouin; Debye, Knudsen, Bragg, Kramers, Dirac, Compton, de Broglie, Born, Bohr; Langmuir, Planck, Curie, Lorentz, Einstein, Langevin, Guye, Wilson, Richardson.
En esta foto aparecen algunas personas que hicieron grandes contribuciones a la física moderna en la que está basada la química computacional.[1]
La química computacional consiste en la predicción de las propiedades de sistemas atómicos y moleculares mediante cálculos complejos con ordenador basados en las leyes fundamentales de la mecánica cuántica.[2] El origen de estas leyes es el hecho de que el mundo microscópico se rige por normas radicalmente diferentes de aquellas que controlan el mundo macroscópico. La ecuación central de la mecánica cuántica es la ecuación de Schrödinger, que indica como calcular la función de onda de un sistema y permite determinar todas sus propiedades como su energía o su distribución en el espacio. Esta ecuación fundamental se formula de la siguiente forma: $$\left[-\frac{h^2}{8\pi^2 m}\nabla^2 + V \right] \Psi = \frac{ih}{2\pi}\frac{\partial \Psi}{\partial t}$$ donde \(\Psi\) es la función de onda del sistema, \(m\) es su masa, \(h\) es la constante de Plank y \(V\) es el campo de potencial que describe las interacciones dentro del sistema.

Las soluciones de la ecuación de Schrödinger son diferentes funciones de onda que corresponden, cada una, a un estado estacionario del sistema.[3] El producto de estas funciones de onda por su conjugado matemático (\(|\Psi^*\Psi|\)) da la distribución de probabilidades de encontrar el sistema en una cierta zona del espacio. Si el campo de potencial, \(V\), no varía en el tiempo, la ecuación de Schrödinger se puede simplificar por separación de variables. Escribiendo la función de onda como el producto de la parte de la función que depende de las coordenadas espaciales, indicadas genéricamente por como (\(x\), \(y\), \(z\)) por la parte que depende del tiempo: $$ \Psi(x, y, z, t)=\psi(x, y, z)\tau(t)$$ la ecuación de Schrödinger se puede separar en dos fórmulas diferentes. Una que depende exclusivamente del tiempo y otra que solo es función de la posición de la partícula: $$\hat{H} \psi(x,y,z) = E \psi(x,y,z)$$ donde \(E\) es la energía de la partícula y \(H\) es el denominado operador Hamiltoniano, que engloba la energías cinética y potencial de la partícula: $$\hat{H} = \hat{T} + \hat{V}$$ En un sistema molecular, \(\psi\) es función de las posiciones de todos los electrones y núcleos de la molécula. La energía cinética total es la suma de las energías cinéticas de todas las \(k\) partículas del sistema: $$\hat{T} = -\frac{h^2}{8\pi^2} \sum_{k} \frac{1}{m_k} \left( \frac{\partial^2}{\partial x_k^2}+\frac{\partial^2}{\partial y_k^2}+\frac{\partial^2}{\partial z_k^2} \right)$$ y la energía potencial es el resultado de la repulsión y atracción eléctrica entre los electrones con carga \(-e\) y los núcleos con carga \(Ze\) (donde \(Z\) es el número atómico de cada núcleo), separados por una distancia de \(r\) o \(R\): $$\hat{V} = \frac{1}{4\pi \varepsilon_0}\left[-\sum_{i}\sum_{I}\left(\frac{Z_I e^2}{\Delta r_{iI}}\right)+ \sum_{i}\sum_{j < i}\left(\frac{e^2}{\Delta r_{ij}}\right)+ \sum_{I}\sum_{J < I}\left(\frac{Z_I Z_J e^2}{\Delta R_{IJ}}\right)\right]$$ Se puede ver como, en la ecuación anterior, el primer sumando corresponde a la atracción entre núcleos y electrones, el segundo a la repulsion entre electrones y el último término a la repulsión entre núcleos.

Así, para calcular la energía de una molécula con \(k\) partículas, incluyendo tanto los electrones como los núcleos atómicos, es necesario resolver la siguiente ecuación diferencial: $$\left[-\frac{h^2}{8\pi^2} \sum_{k} \frac{1}{m_k} \nabla^2 +\frac{1}{4\pi \varepsilon_0}\left[-\sum_{i}\sum_{I}\left(\frac{Z_I e^2}{\Delta r_{iI}}\right)+ \sum_{i}\sum_{j < i}\left(\frac{e^2}{\Delta r_{ij}}\right)+ \sum_{I}\sum_{J < I}\left(\frac{Z_I Z_J e^2}{\Delta R_{IJ}}\right)\right] \right]$$ $$ \psi(x,y,z)= E \psi(x,y,z)$$ Esta ecuación sólo tiene solución exacta para unos pocos sistemas sencillos, como el átomo de hidrógeno. Para estructuras más complejas, la ecuación de Schrödinger se complica tanto que es imposible de resolver. Por lo tanto, se tienen que asumir algunas aproximaciones para calcular la energía de estos sistemas.

La primera aproximación es la de Born-Oppenheimer que asume que como los electrones se mueven mucho más rápido que los núcleos atómicos debido a su diferencia de masa, se puede considerar que los núcleos se encuentran estacionarios. De esta forma, la ecuación de Schrödinger total se divide en dos ecuaciones, una de las cuales (ecuación electrónica de Schrödinger) describe el movimiento electrónico para cada posible geometría molecular. Los autovalores de esta ecuación, una función matemática de las coordenadas nucleares, actúan como campo de fuerzas para el movimiento de los núcleos en la ecuación nuclear de Schrödinger.

Método de Hartree-Fock
Las siguientes aproximaciones entran dentro de lo que se conoce como el método de Hartree-Fock. Este método sustituye la ecuación de onda por una combinación lineal de orbitales moleculares \(\phi\). Esta combinación debe cumplir las mismas condiciones que satisface cualquier función de onda, estas son: normalización (la probabilidad de encontrar la partícula en todo el espacio debe ser 1) y antisimetría (cuando dos electrones se intercambian, la función cambia de signo). Para ello, de partida, se definen los orbitales de tal forma que estén normalizados y como combinación lineal de los mismos se elige la función antisimétrica más simple que existe: el determinante.
Este determinante se compone teniendo en cuenta el principio de exclusión de Pauli, que establece que cada orbital puede alojar dos electrones siempre que tenga diferente espín. De esta forma, se define la función de onda para un sistema de \(n\) electrones como la combinación lineal de orbitales moleculares multiplicados por una función de espín, que se denomina \(\alpha\) cuando el espín vale \(+1/2\) y \(\beta\) cuando es \(-1/2\): \begin{equation} \psi = \frac{1}{\sqrt{n!}} \begin{vmatrix} \phi_1 (1) \alpha & \phi_1 (1) \beta & \phi_2 (1) \alpha & \phi_2 (1) \beta & \cdots & \phi_{\frac{n}{2}} (1) \alpha & \phi_{\frac{n}{2}} (1) \beta\\ \phi_1 (2) \alpha & \phi_1 (2) \beta & \phi_2 (2) \alpha & \phi_2 (2) \beta & \cdots & \phi_{\frac{n}{2}} (2) \alpha & \phi_{\frac{n}{2}} (2) \beta\\ \vdots & \vdots & \vdots & \vdots & & \vdots & \vdots\\ \phi_1 (n) \alpha & \phi_1 (n) \beta & \phi_2 (n) \alpha & \phi_2 (n) \beta & \cdots & \phi_{\frac{n}{2}} (n) \alpha & \phi_{\frac{n}{2}} (n) \beta\\ \end{vmatrix} \end{equation} Este determinante se denomina determinante de Slater y cada fila representa todas las posibles combinaciones para cada electron con los \(n/2\) orbitales disponibles. El factor inicial \(1/\sqrt{n!}\) permite mantener la normalización.

Todos los orbitales que tienen la misma energía constituyen una capa. Cuando existe el mismo número de electrones \(\alpha\) que electrones \(\beta\) en la función de onda, se dice que es un sistema de capa cerrada y la función de onda que hemos definido en la ecuación anterior se denomina función de onda Hartree-Fock restringida (RHF). En esta formulación se obliga a los electrones apareados a ocupar el mismo orbital espacial.
Existen dos aproximaciones en el método de Hartree-Fock para sistemas de capa abierta tales como moléculas excitadas o con electrones desapareados. Una es el método Hartree-Fock de capa abierta restringido (ROHF) que permite electrones individuales en un orbital molecular, manteniendo los electrones apareados en el mismo orbital. La otra aproximación es el método Hartree-Fock sin restringir (UHF) en el que se permite a electrones apareados alojarse en distintas funciones espaciales. Este método tiene en cuenta la repulsión de Pauli por la que los electrones del mismo espín tienden a mantenerse lejos unos de otros.

El método de Hartree-Fock también sustituye el operador de Hamilton por el operador de Fock, donde se considera la interacción promedio entre los distintos electrones, despreciándose las interacciones instantáneas.
El operador de Fock, \(\hat{F}\), es un operador de un solo electrón y nos da la energía de un electrón en un orbital individual \(j\) (ecuación de Hartree-Fock): \begin{equation} \label{hartreemotherfocker} \hat{F}(i)\phi_j (i) = \varepsilon_j \phi_j (i) \end{equation} Consta de un término denominado Hamiltoniano de core de un electrón \begin{equation} \hat{H}^\text{core}(i) = - \frac{h^2}{8 \pi m_i}\nabla^2 - \sum_\alpha \frac{Z_\alpha e^2}{r_{i\alpha}} \end{equation} que representa la energía de un electrón \(i\) en el campo de todos los núcleos y dos factores de corrección, el operador de coulomb \(\hat{J}\) y el operador de intercambio \(\hat{K}\) que, respectivamente, cuantifican la energía potencial de interacción entre el electrón \(i\) y la nube de electrones que lo rodean y aporta la antisimetría necesaria que la función de onda requiere. De esta forma, el operador de Fock es: \begin{equation} \hat{F}(i)=\hat{H}^\text{core}(i)+\sum_{j=1}^{n/2}\left[2 \hat{J}_j(i)-\hat{K}_j(i)\right] \end{equation} El factor de 2 en el operador de Coulomb se debe a que hay dos electrones en cada orbital. Como los operadores dependen de las propiedades de otros electrones, es decir, de la solución, el método requiere de un proceso iterativo hasta que se llega a la convergencia. Por eso, el método de Hartree-Fock se denomina también del campo autoconsistente (SCF).

Ante la complejidad que había alcanzado el método de Hartree-Fock, el físico holandés Clemens Roothaan propuso expresar los orbitales moleculares como combinaciones lineales de funciones monoelectrónicas sencillas a las que llamó funciones de base. Un orbital molecular se define como: \begin{equation} \label{omphi} \phi_j = \sum_{\mu = 1}^m c_{\mu j}\chi_\mu \end{equation} donde \(j\) es el indicador del orbital, \(\mu\) es el número de función de base y \(c\) es el coeficiente de expansión del orbital. El número de funciones de base para representar exactamente el orbital molecular debería ser infinito. En la práctica se elige un número finito, \(m\), de funciones de base. Si \(m\) es lo suficientemente grande y las funciones \(\chi\) se eligen bien, se pueden representar los orbitales moleculares con un error despreciable.
Con esta aproximación, la ecuación de Hartree-Fock queda así: \begin{equation} \sum_\mu c_{\mu j} \hat{F} \chi_\mu = \varepsilon_i \sum_\mu c_{\mu j} \chi_\mu \end{equation} Multiplicando por el conjugado de las funciones de base empleadas \(\chi_\nu^*\) donde \(\nu\) también va de 1 a \(m\) e integrando sobre todo el espacio ambos términos de la ecuación, se llega a la ecuación de Roothaan: \begin{equation} \textbf{FC}=\textbf{SC}\varepsilon \end{equation} o, equivalentemente: \begin{equation} \label{secomp} \sum_{\mu = 1}^m (F_{\mu \nu} - \varepsilon_j S_{\mu \nu}) c_{\mu j} = 0 \quad j=1,\ldots,n \quad \nu=1,\ldots, m \end{equation} donde \(F_{\mu \nu}\) es: \begin{equation} F_{\mu \nu}\equiv\int \chi_\nu^* \hat{F} \chi_\mu dr \end{equation} y \(S_{\mu \nu}\) es: \begin{equation} S_{\mu \nu} \equiv \int \chi_\nu^* \chi_\mu dr \end{equation} siendo \(r\) son las coordenadas espaciales electrónicas.
De esta forma, las ecuaciones de Roothaan forman un sistema de \(n\times m\) ecuaciones homogéneas lineales simultáneas para las \(n\times m\) incógnitas \(c_{\mu j}\) donde \(\mu\) es 1, 2, ..., \(m\) y los coeficientes \(c_{\mu j}\) describen al orbital molecular \(\phi_j\) definido por Roothaan. Para que la solución no sea trivial el primer factor del producto debe ser: \begin{equation} \det(F_{\mu \nu} - \varepsilon_j S_{\mu \nu})=0 \end{equation} Las raíces de esta ecuación son las energías \(\varepsilon_i\) del orbital. A partir de ellas se obtiene la energía total de la molécula. El cálculo de la ecuación de Roothaan también se realiza por iteración debido a que el operador de Fock depende de los orbitales y estos se calculan con el operador de Fock.

Limitaciones del método de Hartree-Fock
Las principales fuentes de error del método de Hartree-Fock son tres:
La primera es aplicar de forma incorrecta la aproximación de Born-Oppenheimer, generalmente en moléculas excitadas.
La segunda se produce en moléculas con átomos pesados, cuyos electrones cercanos al núcleo sufren efectos relativistas al moverse a una velocidad cercana a la de la luz. En este caso es necesario utilizar la versión relativista de la ecuación de Schrödinger: la ecuación de Dirac.
La última fuente de error y fundamental es el tratamiento inadecuado de la correlación electrónica. El método de Hartree-Fock considera que el electrón se encuentra en un campo estacionario, homogéneo e isótropo. No obstante, ante la presencia de un electrón, el resto de electrones se alejan por las repulsiones eléctrica y de Pauli, reduciéndose su energía potencial.

Se han desarrollado varios métodos que tienen en cuenta la correlación electrónica. El método interacción de configuraciones corrige este efecto en Hartree-Fock representando la función de onda con una combinación lineal de determinantes de Slater: \begin{equation} \psi = c_0 \psi_0 + \sum_{s=1}^N b_s \psi_s \end{equation} Cuanto más determinantes se añadan, más se acercará la solución al valor real de la energía.

El método de Møller-Plesset resuelve la ecuación de Schrödinger mediante una aproximación matemática denominada método de perturbaciones. Consiste en calcular la energía de cada orbital con un operador que sí tenga solución, como el hamiltoniano de core de un electrón y posteriormente corregir su valor con la solución de aplicar al sistema un operador formado por la diferencia entre el operador problema y el conocido. La aproximación se puede aplicar a varios órdenes con el inconveniente de que, al no ser un método variacional, la energía calculada puede ser menor que la real. En cálculos mecano-cuánticos, el método de perturbaciones se aplica habitualmente a orden dos (MP2) o cuatro (MP4).

Existe un tercer método que al no estar basado en la resolución de la ecuación de Schrödinger no comete el error de la correlación electrónica. En 1964 el físico austriaco Walter Kohn demostró que la energía de un sistema mecano-cuántico viene únicamente determinada por la densidad electrónica y propuso un método al que llamó del funcional de la densidad (DFT) que facilita la elaboración de ecuaciones cuya solución es la densidad electrónica y la energía de un sistema. Es un método que se aplica frecuentemente a moléculas grandes como las de los sistemas biológicos.

Funciones de base
Todos los métodos anteriores comienzan con la elección de las funciones de base \(\chi_i\). Estas funciones suelen ser orbitales tipo Slater (STO): \begin{equation} f = N r^{n-1} e^{-\zeta r} \end{equation} donde \(N\) es la constante de normalización, \(n\) es un número natural que actúa como número cuántico principal, \(\zeta\) es la carga efectiva del núcleo calculada con las reglas de Slater y \(r\) es la distancia del electrón al núcleo; o funciones tipo gausiana (GTF) con la forma: \begin{equation} g=N x^n y^m z^l e^{-\alpha r^2} \end{equation} donde \(\alpha\) es una constante que determina la anchura de la función, \(N\) es la constante de normalización que depende de \(\alpha\), \(x\), \(y\) y \(z\) son las coordenadas cartesianas con el origen en el núcleo y \(l\), \(m\) y \(n\) son números naturales que proporcionan el mismo comportamiento angular que los orbitales moleculares. Cuando \(i + j + k = 0\), se dice que la gausiana es tipo \(s\); cuando \(i + j + k = 1\), es la gausiana tipo \(p\) y cuando \(i + j + k = 2\), la gausiana es tipo \(d\).

Para un cálculo SCF, cada átomo de la molécula se representa utilizando como al menos un orbital tipo Slater por cada orbital atómico de capa interna y cada orbital atómico de la capa de valencia. Esta representación se denomina base mínima.

Una base extendida es cualquier base mayor que la mínima. Los cálculos SCF con bases mínimas son más fáciles que los cálculos con bases extendidas, aunque estos últimos son considerablemente más precisos. Una base extendida muy común es la base doble zeta (DZ) que reemplaza cada orbital de Slater por una combinación lineal de dos orbitales de Slater con diferentes exponentes orbitales \(\zeta\). También existen bases triple zeta y bases que, por simplicidad, sólo modifican las funciones de la capa de valencia, como la función doble zeta de valencia (VDZ) y triple zeta de valencia (VTZ). Para tener en cuenta la distorsión de los orbitales atómicos debida a la formación de la molécula, se pueden añadir bases polarizadas en las que sus números cuánticos son mayores que el número \(l\) máximo de la capa de valencia del estado fundamental del átomo. Así, existen bases doble zeta más polarización (DZP), triple zeta más polarización (TZP)...

Con las funciones tipo gausiana, en lugar de utilizar funciones individuales se emplean combinaciones lineales de varias funciones: \begin{equation} \chi_\mu = \sum_p d_{\mu p}g_p \end{equation} donde las constantes \(d_{\mu p}\) son los coeficientes de contracción. En este caso, \(\chi_\mu\) se denomina función tipo gausiana contraída (CGTF) y las \(g_p\) gausianas primitivas.

Todas las definiciones de bases con orbitales tipo Slater también se aplican a las bases CGTF. Los cálculos con gausianas requieren definir primero los coeficientes de contracción. Para ello se emplea el resultado de un cálculo preliminar con orbitales STO o GTF. Así la base STO-3G parte de un cálculo con orbitales de Slater que se perfecciona con una combinación lineal de tres gausianas.
La serie 3-21G y la serie 6-31G son funciones de base de valencia desdoblada con CGTF. En la serie 3-21G, cada orbital atómica de capa interna se representa por una combinación lineal de tres gausianas primitivas y en cada orbital de valencia se utilizan dos funciones de base, una combinación de dos primitivas gausianas y una gausiana difusa, que no es más que una primitiva con un exponente orbital muy pequeño.
La serie 6-31G utiliza seis primitivas en la capa interna y tres primitivas y una gausiana difusa en la capa de valencia.
La base 6-31G* añade a la serie 6-31G seis funciones de polarización gausianas tipo \(d\) a los átomos con número atómico superior al litio y diez funciones de polarización tipo \(f\) para los átomos desde el escandio hasta el zinc, el átomo más pesado que puede emplearse con esta base. La base 6-31G** añade a la base 6-31G* tres funciones polarizadas tipo \(p\) a los átomos de hidrógeno y helio.

Los aniones, compuestos con pares solitarios y los dímeros con enlaces de hidrógeno tienen una densidad electrónica significativa a grandes distancias del núcleo. Para mejorar la precisión para tales compuestos se forman las bases 3-21+G y 6-31+G que añaden cuatro funciones altamente difusas (tipo \(s\), \(p_x\), \(p_y\) y \(p_z\)) a las bases 3-21G y 6-31G.
Las series 3-21++G y 6-31++G añaden además una función \(s\) altamente difusa.


[1] - Fuente de la imagen: https://forrifarg.se/
[2] - Exploring chemistry with electronic structure methods, James B. Foresman, Æleen Frisch, Gaussian, Inc., 1996
[3] - Química cuántica, Ira N. Levine, Pearson Educación, 2001
Publicado el 28 de octubre de 2013