Química universitaria — tercer año · Bachelor Year 3

3Química computacional: Hartree–Fock y DFT

Se cree que el ion hidruro de helio, $\ce{HeH+}$ , fue la primera molécula que se formó al enfriarse el Universo joven: un átomo de helio que retiene un protón. Se obtuvo en un tubo de descarga de laboratorio en 1925, pero se buscó en el espacio durante décadas sin éxito, hasta que en 2019 se detectó su raya de rotación en una nebulosa planetaria. Para entonces, su longitud de enlace, sus frecuencias de vibración y las longitudes de onda de sus rayas se conocían desde hacía mucho —por cálculo. Un programa informático que resuelve, de forma aproximada, la ecuación de Schrödinger de los electrones puede predecir la forma, la energía y el espectro de una molécula que nadie ha metido en un frasco. Este capítulo explica cómo funcionan esos programas, desde la separación de núcleos y electrones hasta la teoría de Hartree–Fock y la teoría del funcional de la densidad, y para qué se puede confiar en sus números y para qué no.

Lo que ya sabes

El Capítulo 1 aportó el hamiltoniano, la ecuación de Schrödinger y el átomo de hidrógeno; el Capítulo 2, el determinante de Slater y la integral de intercambio. El volumen del segundo año construyó orbitales moleculares como combinaciones lineales de orbitales atómicos (CLOA), con integrales de solapamiento, de Coulomb y de resonancia y un determinante secular, y resolvió el método de Hückel; también usó los espectros fotoelectrónicos para ver las energías orbitales.

Un clúster de cálculo: filas de servidores en los que se ejecutan hoy la mayoría de los cálculos de química cuántica, a menudo durante horas o días por molécula.

3.1 Born–Oppenheimer y la superficie de energía potencial

El hamiltoniano de una molécula contiene las energías cinéticas de sus núcleos y electrones y todas sus interacciones coulombianas. Los núcleos son al menos 1836 veces más pesados que los electrones: para las mismas fuerzas se mueven mucho más despacio, y los electrones se ajustan casi instantáneamente a cada posición de los núcleos.

Definición 3.1 (Aproximación de Born–Oppenheimer)

En la aproximación de Born–Oppenheimer, la ecuación de Schrödinger electrónica se resuelve con los núcleos fijos en las posiciones $\mathbf R$ ,

\hat H_{\mathrm{el}}\,\psi_{\mathrm{el}}(\mathbf r;\mathbf R) = E_{\mathrm{el}}(\mathbf R)\,\psi_{\mathrm{el}}(\mathbf r;\mathbf R),

y los núcleos se mueven después en el potencial $E_{\mathrm{el}}(\mathbf R) + V_{\mathrm{nn}}(\mathbf R)$ , la energía electrónica más la repulsión de los núcleos.

Definición 3.2 (Superficie de energía potencial, punto estacionario, optimización de geometría)

La superficie de energía potencial (SEP) de una molécula es la función $U(\mathbf R) = E_{\mathrm{el}}(\mathbf R) + V_{\mathrm{nn}}(\mathbf R)$ de sus coordenadas nucleares. Un punto estacionario es un punto en el que se anulan todas las derivadas primeras de $U$ ; un mínimo es una geometría de equilibrio. Buscar uno a partir de una estructura de partida, siguiendo las fuerzas $-\nabla U$ , es una optimización de geometría.

Para una molécula diatómica, la SEP es una curva $U(R)$ ; para una triatómica, una función de tres coordenadas; para $N$ átomos, de $3N - 6$ . La masa de los núcleos no interviene en $\hat H_{\mathrm{el}}$ : $\ce{H2}$ , $\ce{HD}$ y $\ce{D2}$ comparten la misma SEP, la misma longitud de enlace y la misma constante de fuerza, y solo difieren en cómo se mueven los núcleos sobre ella (energía de punto cero, Capítulo 1).

Proposición 3.3 (Frecuencias a partir de la hessiana)

Cerca de un mínimo, $U \approx U_0 + \frac12\sum_{ij}H_{ij}\,\delta q_i\,\delta q_j$ , con la hessiana $H_{ij} = \partial^2U/\partial q_i\partial q_j$ . En coordenadas ponderadas por la masa, $\delta q_i\sqrt{m_i}$ , los valores propios $\lambda_k$ de la hessiana ponderada por la masa dan las frecuencias angulares armónicas $\omega_k = \sqrt{\lambda_k}$ de los modos normales. En un mínimo, todos los $\lambda_k$ son positivos (seis, o cinco para una molécula lineal, son nulos: traslaciones y rotaciones); en un punto de silla de primer orden, exactamente uno es negativo.

Demostración parcial. Para una diatómica, $U \approx U_0 + \frac12k(R - R_e)^2$ con $k = U''(R_e)$ , y el oscilador del Capítulo 1 da $\omega = \sqrt{k/\mu}$ . En general, las ecuaciones clásicas $m_i\ddot q_i = -\sum_jH_{ij}q_j$ se convierten, en coordenadas ponderadas por la masa $x_i = \sqrt{m_i}q_i$ , en $\ddot{\mathbf x} = -\tilde H\mathbf x$ con la matriz simétrica $\tilde H_{ij} = H_{ij}/\sqrt{m_im_j}$ ; sus vectores propios oscilan independientemente a $\sqrt{\lambda_k}$ . Un $\lambda_k$ negativo da una frecuencia imaginaria: en esa dirección la energía baja; es un máximo. Los modos normales del Capítulo 5 son estos vectores propios. ∎

Una superficie de energía potencial vista desde arriba, como un mapa de curvas de nivel sobre dos coordenadas nucleares (esquema). Dos mínimos están unidos, a través de un punto de silla (cruz), por el camino de menor energía (trazo discontinuo): el camino de energía mínima estudiado en el . — Una superficie de energía potencial vista desde arriba, como un mapa de curvas de nivel sobre dos coordenadas nucleares (esquema). Dos mínimos están unidos, a través de un punto de silla (cruz), por el camino de menor energía (trazo discontinuo): el camino de energía mínima estudiado en el Capítulo 12.

3.2 El principio variacional

La ecuación de Schrödinger electrónica no puede resolverse exactamente para ninguna molécula con más de un electrón. Las aproximaciones se juzgan con un solo teorema.

Teorema 3.4 (El principio variacional)

Para toda función de prueba normalizable $\phi$ que cumpla las condiciones de contorno del problema,

E[\phi] = \frac{\langle\phi|\hat H|\phi\rangle}{\langle\phi|\phi\rangle} \ge E_0,

donde $E_0$ es el valor propio más bajo de $\hat H$ ; la igualdad solo se da si $\phi$ es una función propia del estado fundamental. Es el principio variacional.

Cuanto menor es la energía de una función de prueba, mejor es; todo parámetro ajustable se fija minimizando $E$ .

Proposición 3.5 (Helio con una carga apantallada)

Para el helio, la función de prueba $\phi = \eu^{-\zeta(r_1 + r_2)}$ (en unidades atómicas: longitudes en $a_0$ , energías en $E_h$ ) da $E(\zeta) = \zeta^2 - \frac{27}{8}\zeta$ , mínima para $\zeta = \frac{27}{16}$ , con $E = -(\frac{27}{16})^2 E_h = -2.8477\,E_\mathrm{h}$ .

Demostración. Para un electrón en $\eu^{-\zeta r}$ , $\langle T\rangle = \zeta^2/2$ y $\langle 1/r\rangle = \zeta$ ; la repulsión de dos electrones en la misma función $1s$ es $\frac58\zeta$ (una integral clásica, que se admite). Así, $E = 2\cdot\frac12 \zeta^2 - 2\cdot2\zeta + \frac58\zeta = \zeta^2 - \frac{27}{8}\zeta$ , y $\dd E/\dd \zeta = 0$ da $\zeta = 27/16 = 1.6875$ . ∎

Ejemplo 3.6 (¿Cuán buena es?)

La energía total del helio es menos la suma de sus dos energías de ionización, $-(24.5874 + 54.4178)\,\mathrm{eV} = -2.9034\,E_\mathrm{h}$ . La función de un parámetro queda $0.056\,E_\mathrm{h}$ ( $1.5\,\mathrm{eV}$ ) demasiado alta: cada electrón ve un núcleo apantallado hasta $Z_{\mathrm{eff}} = 1.69$ por el otro, como sugieren las reglas de Slater, pero los electrones también se esquivan de un modo que ningún producto de dos funciones puede describir.

Proposición 3.7 (Variación lineal)

Para una función de prueba $\phi = \sum_{i=1}^n c_i\chi_i$ construida sobre $n$ funciones fijas, con $H_{ij} = \langle\chi_i|\hat H|\chi_j\rangle$ y $S_{ij} = \langle\chi_i|\chi_j \rangle$ , los valores estacionarios de $E$ son las raíces de $\det(H - ES) = 0$ , y la raíz más baja es una cota superior de $E_0$ .

Demostración. $E\sum_{ij}c_ic_jS_{ij} = \sum_{ij}c_ic_jH_{ij}$ (coeficientes reales). Se deriva respecto de $c_k$ y se impone $\partial E/\partial c_k = 0$ : $\sum_j(H_{kj} - ES_{kj})c_j = 0$ para todo $k$ . Una solución no nula exige que se anule el determinante. La raíz más baja es $E[\phi]$ para su vector propio, luego $\ge E_0$ . ∎

El método de Hückel del volumen del segundo año es esta proposición con orbitales $p$ y $H_{ij}$ empíricos. Todos los métodos siguientes son la misma idea con mejores funciones y mejores elementos de matriz.

3.3 Conjuntos de base

Definición 3.8 (Conjunto de base, base mínima, orbitales de tipo Slater y de tipo gaussiano)

Un conjunto de base es el conjunto de funciones fijas $\chi_i$ , centradas en los átomos, a partir de las cuales se construyen los orbitales moleculares. Una base mínima tiene una función por cada orbital atómico ocupado (una $1s$ para H, cinco para C). Un orbital de tipo Slater (STO) tiene la forma radial $r^{n-1}\eu^{-\zeta r}$ ; un orbital de tipo gaussiano (GTO) tiene la forma $\eu^{-\alpha r^2}$ multiplicada por un polinomio en $x$ , $y$ , $z$ . Una función gaussiana contraída es una combinación fija $\sum_kd_k\,\eu^{-\alpha_kr^2}$ de gaussianas primitivas.

Las funciones de Slater tienen la forma correcta (un punto anguloso en el núcleo y una cola exponencial), pero sus integrales bielectrónicas sobre varios centros son muy costosas; las funciones gaussianas tienen la forma equivocada, pero el producto de dos gaussianas centradas en puntos distintos es una gaussiana centrada en un punto entre ambos, de modo que toda integral se reduce a una fórmula cerrada. Las contracciones obtienen lo mejor de ambas: en la base STO-3G, cada función de Slater se sustituye por una suma fija de tres gaussianas ajustadas a ella.

Ejemplo 3.9 (La función STO-3G del hidrógeno)

La función $1s$ del hidrógeno en STO-3G es $0.15433\,g(3.42525) + 0.53533\, g(0.62391) + 0.44463\,g(0.16886)$ , donde $g(\alpha)$ es la gaussiana normalizada $(2\alpha/\pi)^{3/4}\eu^{-\alpha r^2}$ (exponentes en $a_0^{-2}$ ). Imita una función de Slater de exponente $\zeta = 1.24$ —un átomo de hidrógeno ligeramente comprimido, como lo está en las moléculas. Su solapamiento con la función $1s$ exacta del mismo exponente es $0.9998$ ; solo se pierden el punto anguloso en el núcleo y la cola lejana (figura siguiente).

Definición 3.10 (Funciones de valencia desdoblada, de polarización y difusas)

Una base de valencia desdoblada describe cada orbital de valencia mediante dos (o más) funciones de tamaños distintos, de modo que el orbital pueda crecer o encogerse en una molécula. Una función de polarización tiene un $l$ mayor que los orbitales ocupados ( $p$ en H, $d$ en C) y permite que la densidad se desplace fuera del átomo; una función difusa tiene un exponente pequeño y describe aniones e interacciones débiles lejos de los núcleos.

Un nombre como 6-31G(d) dice: seis primitivas para cada función de núcleo, la valencia desdoblada en una contracción de tres y una sola primitiva, y funciones de polarización $d$ en los átomos pesados. Los resultados convergen hacia el límite de base completa a medida que crece la base; un cálculo solo vale lo que vale su base.

Izquierda: la contracción STO-3G (trazo discontinuo) reproduce de cerca la función de Slater 1s salvo en el núcleo; una sola gaussiana (punteado) no lo hace. Derecha: la energía de HeH+ a 1.4632\,a_0 durante las iteraciones de campo autoconsistente (STO-3G); converge por arriba hacia -2.8418\,E_ h en pocos ciclos, y cada iteración obedece al principio variacional. — Izquierda: la contracción STO-3G (trazo discontinuo) reproduce de cerca la función de Slater $1s$ salvo en el núcleo; una sola gaussiana (punteado) no lo hace. Derecha: la energía de $\ce{HeH+}$ a $1.4632\,a_0$ durante las iteraciones de campo autoconsistente (STO-3G); converge por arriba hacia $-2.8418\,E_\mathrm{h}$ en pocos ciclos, y cada iteración obedece al principio variacional.

3.4 Teoría de Hartree–Fock

La función de onda antisimétrica más sencilla de $N$ electrones es un único determinante de Slater. La teoría de Hartree–Fock encuentra el mejor.

Definición 3.11 (Método de Hartree–Fock, operador de Fock, campo autoconsistente)

El método de Hartree–Fock aproxima el estado fundamental por el determinante de Slater único de menor energía. Sus orbitales son funciones propias del operador de Fock

\hat f(1) = \hat h(1) + \sum_b\big[2\hat J_b(1) - \hat K_b(1)\big],

en el caso de capa cerrada: $\hat h$ es la energía cinética y la atracción hacia los núcleos de un electrón, $\hat J_b$ la repulsión de la nube de carga del orbital $b$ , y $\hat K_b$ el operador de intercambio, cuyos valores esperados son las integrales de intercambio del Capítulo 2. Como $\hat f$ depende de los orbitales que determina, las ecuaciones se resuelven por iteración hasta que los orbitales dejan de cambiar: un campo autoconsistente (SCF).

Teorema 3.12 (Ecuaciones de Roothaan–Hall)

Escribiendo cada orbital en una base, $\phi_a = \sum_\mu C_{\mu a}\chi_\mu$ , las ecuaciones de Hartree–Fock de capa cerrada se convierten en el problema matricial de valores propios

\mathbf F\mathbf C = \mathbf S\mathbf C\boldsymbol\varepsilon, \qquad F_{\mu\nu} = h_{\mu\nu} + \sum_{\lambda\sigma}P_{\lambda\sigma}\big[(\mu\nu|\sigma\lambda) - \tfrac12(\mu\lambda|\sigma\nu)\big],

con la matriz de densidad $P_{\lambda\sigma} = 2\sum_a^{\mathrm{occ}}C_{\lambda a}C_{\sigma a}$ y las integrales bielectrónicas $(\mu\nu|\lambda\sigma) = \iint\chi_\mu(1)\chi_\nu(1) r_{12}^{-1}\chi_\lambda(2)\chi_\sigma(2)$ (unidades atómicas).

Demostración parcial. La energía del determinante es $E = \sum_{\mu\nu}P_{\mu\nu}h_{\mu\nu} + \frac12 \sum P_{\mu\nu}P_{\lambda\sigma}[(\mu\nu|\sigma\lambda) - \frac12(\mu\lambda|\sigma\nu)]$ , una función cuadrática de los coeficientes. Se minimiza con las restricciones de que los orbitales sigan siendo ortonormales, $\mathbf C^{\mathsf T}\mathbf S\mathbf C = \mathbf 1$ , con un multiplicador de Lagrange para cada restricción: las condiciones de estacionariedad son $\mathbf F\mathbf C = \mathbf S\mathbf C\boldsymbol\varepsilon$ , y los multiplicadores forman una matriz que puede diagonalizarse mediante una rotación de los orbitales ocupados entre sí. El álgebra se admite. ∎

Método 3.13 (El procedimiento SCF)

Calcular una sola vez las integrales $S_{\mu\nu}$ , $h_{\mu\nu}$ y $(\mu\nu|\lambda\sigma)$ .
Proponer una matriz de densidad de partida (por ejemplo, a partir de los vectores propios de $\mathbf h$ sola, la estimación “de núcleo”).
Construir $\mathbf F$ a partir de $\mathbf P$ ; resolver $\mathbf F\mathbf C = \mathbf S\mathbf C \boldsymbol\varepsilon$ (ortogonalizar con $\mathbf S^{-1/2}$ y después diagonalizar).
Llenar los orbitales más bajos, con dos electrones cada uno; formar la nueva $\mathbf P$ y la energía.
Repetir desde el paso 3 hasta que la energía y $\mathbf P$ cambien menos que un umbral.

El bucle de campo autoconsistente de un cálculo de Hartree–Fock.

Ejemplo 3.14 ( $\ce{H2}$ en la base mínima)

Con una función STO-3G por átomo, el orbital enlazante queda fijado por la simetría, $\sigma_g \propto \chi_A + \chi_B$ , y el SCF necesita un solo paso. A $R = 1.4\,a_0$ , el solapamiento es $S_{AB} = 0.6593$ y la energía $-1.1167\,E_\mathrm{h}$ ; el mínimo está a $1.346\,a_0$ ( $71.2\,\mathrm{pm}$ , frente a los $74.1\,\mathrm{pm}$ medidos) y $-1.1175\,E_\mathrm{h}$ . La figura siguiente compara la curva con la verdadera, construida a partir de datos espectroscópicos medidos y de la energía de disociación (una curva de Morse, Capítulo 6).

La energía de H2 en función de la distancia internuclear: Hartree–Fock restringido en la base STO-3G (trazo continuo) y la curva verdadera (trazo discontinuo, forma de Morse construida a partir de los valores medidos de D_0, _e, _ex_e, r_e). Cerca del mínimo, RHF es bueno; a R grande sube muy por encima de dos átomos de hidrógeno. — La energía de $\ce{H2}$ en función de la distancia internuclear: Hartree–Fock restringido en la base STO-3G (trazo continuo) y la curva verdadera (trazo discontinuo, forma de Morse construida a partir de los valores medidos de $D_0$ , $\omega_e$ , $\omega_ex_e$ , $r_e$ ). Cerca del mínimo, RHF es bueno; a $R$ grande sube muy por encima de dos átomos de hidrógeno.

Proposición 3.15 (Por qué el Hartree–Fock restringido falla en la disociación)

En la base mínima, la función de onda RHF de $\ce{H2}$ es, salvo normalización, $\sigma_g(1)\sigma_g(2) \propto \chi_A(1)\chi_B(2) + \chi_B(1)\chi_A(2) + \chi_A(1)\chi_A(2) + \chi_B(1)\chi_B(2)$ : a cualquier distancia, mitad covalente y mitad iónica ( $\ce{H+ H-}$ ). A $R$ grande, su energía tiende a la media de dos átomos neutros y de un par de iones, no a la de dos átomos.

Demostración. Al desarrollar $(\chi_A + \chi_B)(1)(\chi_A + \chi_B)(2)$ , los dos términos cruzados ponen un electrón en cada átomo y los dos términos cuadrados ambos en el mismo átomo, con pesos iguales sea cual sea $R$ . El determinante no puede cambiar esos pesos. ∎

Definición 3.16 (Correlación electrónica, energía de correlación)

La correlación electrónica es la parte de la evitación mutua de los electrones que no describe un único determinante (en el que los electrones de espines opuestos se mueven independientemente, cada uno en el campo medio del otro). La energía de correlación es la energía exacta menos la energía de Hartree–Fock en una base completa; siempre es negativa.

La energía de correlación es de unos $-0.04\,E_\mathrm{h}$ por par de electrones, un uno por ciento de la energía total pero comparable a la energía de una reacción: los métodos que la añaden (interacción de configuraciones, teoría de perturbaciones, clústeres acoplados) cuestan mucho más que Hartree–Fock, cuyo coste ya crece aproximadamente como la cuarta potencia del número de funciones de base.

Teorema 3.17 (Teorema de Koopmans)

En la teoría de Hartree–Fock, la energía necesaria para arrancar un electrón del orbital ocupado $a$ , manteniendo congelados todos los demás orbitales, es $-\varepsilon_a$ . Las energías de ionización medidas son, por tanto, aproximadamente las opuestas de las energías orbitales: es el teorema de Koopmans.

Demostración. La energía de un determinante de capa cerrada es $\sum_a2h_{aa} + \sum_{ab}(2J_{ab} - K_{ab})$ , y $\varepsilon_a = h_{aa} + \sum_b(2J_{ab} - K_{ab})$ . Se arranca un electrón de $a$ sin cambiar los orbitales: los términos perdidos son $h_{aa}$ , las interacciones de ese electrón con todos los demás, $\sum_b(2J_{ab} - K_{ab}) - J_{aa}$ , más $J_{aa}$ con su antiguo compañero —en total, exactamente $\varepsilon_a$ . Así, $E^+ - E = -\varepsilon_a$ . ∎

Ejemplo 3.18 (Ionización de $\ce{H2}$ )

A $1.4\,a_0$ , la energía orbital STO-3G de $\sigma_g$ es $-0.5782\,E_\mathrm{h}$ : Koopmans predice $15.73\,\mathrm{eV}$ , frente a los $15.43\,\mathrm{eV}$ medidos. Los orbitales congelados sobrestiman la energía del ion (que se relajaría), y la correlación ausente actúa en sentido contrario; en las ionizaciones de valencia los dos errores se compensan en gran parte, y por eso los espectros fotoelectrónicos del volumen del segundo año pueden leerse con diagramas de orbitales.

3.5 Teoría del funcional de la densidad, y lo que puede decir un cálculo

La función de onda de $N$ electrones depende de $3N$ coordenadas; la densidad electrónica depende de tres.

Definición 3.19 (Densidad electrónica, teoría del funcional de la densidad)

La densidad electrónica $\rho(\mathbf r)$ es el número de electrones por unidad de volumen en $\mathbf r$ , sumado sobre todos los electrones: $\int\rho\,\dd\tau = N$ . La teoría del funcional de la densidad (DFT) calcula la energía del estado fundamental como un funcional de $\rho$ . En la práctica, $\rho = \sum_a|\phi_a|^2$ se construye a partir de orbitales de Kohn–Sham, los orbitales de unos electrones independientes ficticios que tienen la misma densidad que los reales, y todos los efectos de intercambio y de correlación se reúnen en un funcional de intercambio–correlación $E_{\mathrm{xc}}[\rho]$ .

Teorema 3.20 (Hohenberg–Kohn)

La densidad del estado fundamental de un sistema de electrones determina el potencial externo (los núcleos) salvo una constante, y por tanto el hamiltoniano y todas las propiedades del estado fundamental; y el funcional de la energía de la densidad es mínimo en la densidad verdadera.

Demostración. Admitido a este nivel. ∎

La demostración corresponde a cursos más avanzados. Garantiza que existe un funcional exacto, no cuál es: todo funcional práctico es una aproximación, elegida dentro de una jerarquía —densidad local, correcciones de gradiente, híbridos que incorporan una parte de intercambio de Hartree–Fock— y calibrada sobre datos de referencia. Las ecuaciones de Kohn–Sham se resuelven con el mismo bucle SCF que Hartree–Fock, con un coste similar, pero incluyen la correlación; por eso la DFT se ha convertido en el caballo de batalla de la química.

Método 3.21 (Leer un cálculo)

Comprobar que el SCF ha convergido y que la optimización de geometría ha terminado con fuerzas por debajo del umbral.
Hacer un cálculo de frecuencias al mismo nivel: todas las frecuencias reales indican un mínimo; una frecuencia imaginaria, una estructura de transición, cuyo modo normal es el movimiento a través de la barrera.
Comparar lo comparable: diferencias de energía entre estructuras calculadas a un mismo nivel, nunca energías absolutas de niveles distintos.
Comparar con el experimento conociendo los errores típicos del nivel: las frecuencias armónicas suelen escalarse a la baja por 0.9 a 0.97; las longitudes de enlace son buenas a 1 pm aproximadamente con un funcional híbrido y una base polarizada; las energías de reacción, a unos pocos $\mathrm{kJ}/\mathrm{mol}$ en el mejor de los casos.

Método 3.22 (Elegir un nivel de teoría)

Geometrías y frecuencias de moléculas orgánicas: un funcional híbrido con una base de valencia desdoblada polarizada. Interacciones débiles: añadir una corrección de dispersión y funciones difusas. Barreras de reacción y energética precisa: métodos de función de onda correlacionados sobre geometrías DFT, hasta donde lo permita el presupuesto. Sistemas grandes (proteínas, superficies): una pequeña región tratada cuánticamente dentro de un campo de fuerzas clásico. Probar siempre el nivel elegido en una molécula emparentada cuya respuesta se conozca.

En el laboratorio — Un experimento computacional

Un cálculo se lleva a cabo como un experimento y se registra en el cuaderno: el programa y su versión, el método, la base, los umbrales de convergencia, la estructura de partida, y los ficheros de salida, conservados como datos brutos. Una secuencia típica es: construir la molécula, optimizar su geometría a un nivel modesto, calcular las frecuencias para confirmar un mínimo y después calcular la energía a un nivel más alto sobre esa geometría.

Historia — Una molécula conocida primero por el cálculo

$\ce{HeH+}$ se obtuvo en el laboratorio en 1925. Sus propiedades se calcularon con detalle a partir de la década de 1930, y se predijo que sería abundante en algunos gases astrofísicos; pero su primera raya de rotación está en el infrarrojo lejano, absorbida por la atmósfera. En 2019, un observatorio aerotransportado la detectó en la nebulosa planetaria NGC 7027, a la frecuencia que los cálculos y los espectros de laboratorio habían predicho.

3.6 Ejercicios

Ejercicio 3.1 ★

Para el átomo de hidrógeno, se toma la función de prueba $\eu^{-\alpha r^2}$ . En unidades atómicas, $E(\alpha) = \frac32\alpha - 2\sqrt{2\alpha/\pi}$ . Encuéntrense el mejor $\alpha$ y la energía correspondiente; compárese con $-\frac12E_h$ .

Solución

Solución de Ejercicio 3.1.

$\dd E/\dd\alpha = \frac32 - \sqrt{2/\pi\alpha} = 0$ da $\alpha = 8/9\pi = 0.283\, a_0^{-2}$ y $E = \frac32\alpha - 2\sqrt{2\alpha/\pi} = -4/3\pi = -0.4244\,E_\mathrm{h}$ , un 15 % por encima del valor exacto $-0.5$ : una sola gaussiana no tiene ni el punto anguloso ni la cola de la función $1s$ .

Ejercicio 3.2 ★

¿Cuántas funciones de base usa un cálculo sobre el benceno en la base STO-3G? ¿Y en 6-31G(d), con seis funciones $d$ cartesianas por carbono y dos funciones $s$ por hidrógeno?

Solución

Solución de Ejercicio 3.2.

STO-3G: 5 funciones por C ( $1s$ , $2s$ , tres $2p$ ) y 1 por H: $6 \times 5 + 6 = 36$ . 6-31G(d): por C, 1 (núcleo) + $2 \times 4$ (valencia desdoblada) + 6 ( $d$ ) $= 15$ ; por H, 2: $6 \times 15 + 6 \times 2 = 102$ .

Ejercicio 3.3 ★

Explíquese por qué $\ce{H2}$ y $\ce{D2}$ tienen la misma longitud de enlace y la misma constante de fuerza pero números de onda de vibración distintos. ¿Qué razón cabe esperar?

Solución

Solución de Ejercicio 3.3.

El hamiltoniano electrónico no contiene las masas nucleares: dentro de la aproximación de Born–Oppenheimer, la SEP, y por tanto $r_e$ y $k = U''(r_e)$ , es la misma. $\tilde\omega = \sqrt{k/\mu}/2\pi c$ depende de la masa reducida, que se duplica: razón $\sqrt2 = 1.414$ (medida: 1.413).

Ejercicio 3.4 ★

Un cálculo de frecuencias sobre una estructura de fórmula $\ce{C2H5F}$ da 17 frecuencias reales y una frecuencia imaginaria, $487\iu$ $\mathrm{cm}^{-1}$ . ¿Qué tipo de punto estacionario es? ¿Cuántas frecuencias cabía esperar en total?

Solución

Solución de Ejercicio 3.4.

Ocho átomos dan $3N - 6 = 18$ vibraciones: 17 reales y una imaginaria. Un valor propio negativo de la hessiana: un punto de silla de primer orden, una estructura de transición (aquí, por ejemplo, de una rotación o de una eliminación).

Ejercicio 3.5 ★★

Con $E(\zeta) = \zeta^2 - \frac{27}{8}\zeta$ , calcúlese la mejor energía variacional del helio en $E_\mathrm{h}$ y en eV, y su error respecto del valor exacto $-2.9034\,E_\mathrm{h}$ . ¿Qué fracción de la energía total representa el error?

Solución

Solución de Ejercicio 3.5.

$\zeta = 27/16$ , $E = -(27/16)^2 = -2.8477\,E_\mathrm{h} = -77.49\,\mathrm{eV}$ . Error: $-2.8477 + 2.9034 = 0.0557\,E_\mathrm{h} = 1.52\,\mathrm{eV}$ , es decir, un 1.9 % de la energía total —pero más que muchas energías de reacción.

Ejercicio 3.6 ★★

Dos funciones de base dan $H_{11} = -1.0\,E_\mathrm{h}$ , $H_{22} = -0.5\,E_\mathrm{h}$ , $H_{12} = -0.2\,E_\mathrm{h}$ , $S_{12} = 0.3$ . Resuélvase la ecuación secular. ¿Queda la raíz más baja por debajo de $H_{11}$ ?

Solución

Solución de Ejercicio 3.6.

$(-1 - E)(-0.5 - E) - (-0.2 - 0.3E)^2 = 0$ , es decir, $0.91E^2 + 1.38E + 0.46 = 0$ : $E = -1.0217\,E_\mathrm{h}$ y $-0.4947\,E_\mathrm{h}$ . Sí: incorporar la segunda función baja la energía por debajo de $H_{11}$ , como permite el principio variacional.

Ejercicio 3.7 ★★

Un cálculo de Hartree–Fock tarda 10 minutos para el benceno en STO-3G. Suponiendo un coste proporcional a la cuarta potencia del número de funciones de base, estímese el tiempo en 6-31G(d).

Solución

Solución de Ejercicio 3.7.

$(102/36)^4 = 64$ : unos 640 minutos, unas 11 horas.

Ejercicio 3.8 ★★

La energía orbital STO-3G de $\ce{H2}$ a $1.4\,a_0$ es $-0.5782\,E_\mathrm{h}$ . Dese la estimación de Koopmans de la energía de ionización en eV y compárese con $15.43\,\mathrm{eV}$ . Nómbrense dos razones de la diferencia y el sentido de cada una.

Solución

Solución de Ejercicio 3.8.

$0.5782 \times 27.211 = 15.73\,\mathrm{eV}$ , 0.30 eV por encima del valor medido. Orbitales congelados: el ion real se relaja y queda más bajo, de modo que el valor de Koopmans es demasiado alto. Correlación ausente: la molécula neutra (dos electrones) tiene más energía de correlación que el ion (un electrón), lo que hace mayor el valor verdadero. Aquí gana el primer error.

Ejercicio 3.9 ★★

En STO-3G, la energía de un átomo de hidrógeno es $-0.4666\,E_\mathrm{h}$ y la energía RHF de $\ce{H2}$ a $10\,a_0$ es $-0.5960\,E_\mathrm{h}$ . ¿Cuánto queda la energía RHF por encima de la de dos átomos a esa distancia, en eV? Explíquese con los términos iónicos.

Solución

Solución de Ejercicio 3.9.

$-0.5960 - 2(-0.4666) = 0.3372\,E_\mathrm{h} = 9.18\,\mathrm{eV}$ por encima de dos átomos. La función RHF conserva un 50 % de $\ce{H+ H-}$ , y separar un protón y un hidruro cuesta la diferencia entre la energía de ionización de H y la afinidad electrónica de H, que la energía incorpora en la media.

Ejercicio 3.10 ★★★

Para $\phi = c_1\chi_1 + c_2\chi_2$ con funciones reales, escríbase $E(c_1,c_2)$ y dedúzcanse las dos ecuaciones seculares imponiendo $\partial E/\partial c_1 = \partial E/\partial c_2 = 0$ .

Solución

Solución de Ejercicio 3.10.

$E(c_1^2 + 2c_1c_2S + c_2^2) = c_1^2H_{11} + 2c_1c_2H_{12} + c_2^2H_{22}$ . Derivando respecto de $c_1$ con $\partial E/\partial c_1 = 0$ : $E(2c_1 + 2c_2S) = 2c_1H_{11} + 2c_2H_{12}$ , es decir, $(H_{11} - E)c_1 + (H_{12} - ES)c_2 = 0$ ; del mismo modo, $(H_{12} - ES)c_1 + (H_{22} - E)c_2 = 0$ .

Ejercicio 3.11 ★★★

La energía RHF/STO-3G de $\ce{H2}$ vale $-1.116871$ , $-1.117501$ y $-1.116\,714\,E_\mathrm{h}$ a $R = 1.30$ , 1.35 y $1.40\,a_0$ . Estímese la constante de fuerza por diferencias finitas, conviértase a $\mathrm{N}/\mathrm{m}$ ( $1\,E_\mathrm{h}/{a_0}^{2} = 1556.9\,\mathrm{N}/\mathrm{m}$ ) y calcúlese el número de onda armónico. Compárese con $\tilde\omega_e = 4401\,\mathrm{cm}^{-1}$ .

Solución

Solución de Ejercicio 3.11.

$k \approx [E(1.30) - 2E(1.35) + E(1.40)]/(0.05)^2 = 0.001417/0.0025 = 0.567\,E_\mathrm{h}/{a_0}^{2} = 882\,\mathrm{N}/\mathrm{m}$ . Con $\mu = m_H/2 = 8.37 \times 10^{-28}\,\mathrm{kg}$ , $\omega = \sqrt{k/\mu} = 1.03 \times 10^{15}\,\mathrm{s}^{-1}$ y $\tilde\omega = \omega/2\pi c = 5452\,\mathrm{cm}^{-1}$ , un 24 % por encima del $\tilde\omega_e$ medido: una curva RHF en base mínima es demasiado empinada (sin correlación, con una base demasiado pequeña), y por eso las frecuencias calculadas se escalan a la baja.

Ejercicio 3.12 ★★★

Debe calcularse la diferencia de energía entre dos confórmeros de un azúcar, que se sabe que difieren en unos pocos $\mathrm{kJ}/\mathrm{mol}$ , uno de ellos con un enlace de hidrógeno intramolecular. Elíjanse y justifíquense un método, un conjunto de base y las comprobaciones que se harían.

Solución

Solución de Ejercicio 3.12.

Unas diferencias de pocos $\mathrm{kJ}/\mathrm{mol}$ en las que interviene un enlace de hidrógeno exigen correlación y dispersión: un funcional híbrido con corrección de dispersión, o mejor un método de función de onda correlacionado para las energías finales, con una base triple desdoblada polarizada que incluya funciones difusas. Se optimizan ambos confórmeros al mismo nivel, se confirman los mínimos con las frecuencias (todas reales), se añaden las energías de punto cero, se pone a prueba la base con un cálculo mayor y se compara con un sistema emparentado de energía conformacional conocida.

3.7 Problema: HeH $^+$ , la primera molécula

Problema 3.1

Problema de fin de semana — un cálculo de Hartree–Fock en base mínima del ion hidruro de helio: la base, la estimación de núcleo, el orbital convergido y la unión de un protón al helio

Se calcula $\ce{HeH+}$ a $R = 1.4632\,a_0$ en la base STO-3G (exponentes de He 6.3624, 1.1589, 0.31365; exponentes de H 3.4253, 0.62391, 0.16886; los mismos tres coeficientes). La función 1 está sobre He, la función 2 sobre H. El programa da (unidades atómicas)

\mathbf S = \begin{pmatrix}1 & 0.5368\\ 0.5368 & 1\end{pmatrix},\quad \mathbf h = \begin{pmatrix}-2.5983 & -1.4318\\ -1.4318 & -1.7318\end{pmatrix},\quad \mathbf F_{\mathrm{final}} = \begin{pmatrix}-1.5902 & -1.0610\\ -1.0610 & -0.8340 \end{pmatrix},

el orbital ocupado $\phi = 0.8766\chi_1 + 0.2025\chi_2$ , las energías orbitales $-1.6328$ y $-0.1725\,E_\mathrm{h}$ , y las energías SCF $-2.7978$ , $-2.8404$ , $-2.8418$ , $-2.8418$ . Un átomo de helio en la misma base tiene $-2.8078\,E_\mathrm{h}$ . $1\,E_\mathrm{h} = 27.211\,\mathrm{eV} = 2625.5\,\mathrm{kJ}/\mathrm{mol}$ .

Parte I — La base.

¿Qué significa “STO-3G”?
¿Por qué los exponentes del helio son mayores que los del hidrógeno? Compruébese que su razón es $(1.69/1.24)^2$ .
¿Cuántas funciones de base, electrones y orbitales ocupados hay?
Calcúlese la energía de repulsión nuclear.
Interprétese $S_{12} = 0.5368$ .
¿Por qué $h_{11}$ es más bajo que $h_{22}$ ?

Parte II — La estimación de núcleo.

Escríbase la ecuación secular $\det(\mathbf h - \varepsilon\mathbf S) = 0$ como una ecuación de segundo grado en $\varepsilon$ .
Resuélvase.
¿Por qué no es más que un punto de partida?
Dese la fórmula de la matriz de Fock en función de $\mathbf h$ , de la matriz de densidad y de las integrales bielectrónicas.
¿Qué interacciones entre electrones contiene $\mathbf F$ que le faltan a $\mathbf h$ ?
Explíquese por qué el procedimiento debe iterarse.

Parte III — El orbital convergido.

Compruébese que $\phi$ está normalizado.
Las poblaciones son los elementos diagonales de $P\mathbf S$ , $2c_1^2 + 2c_1c_2S$ sobre He y $2c_2^2 + 2c_1c_2S$ sobre H. Calcúlense.
¿Dónde está la carga positiva? ¿Se describe mejor la molécula como $\ce{He + H+}$ o como $\ce{He+ + H}$ ?
Úsese el teorema de Koopmans para estimar la energía necesaria para arrancar un electrón de $\ce{HeH+}$ .
Calcúlese la energía electrónica $E_{\mathrm{el}} = E - V_{\mathrm{nn}}$ .
¿Cuántos ciclos fueron necesarios para converger a $10^{-4}\,E_\mathrm{h}$ ? ¿Por qué disminuyó la energía en cada ciclo?
Con los exponentes del helio escalados a $\zeta = 2.0925$ en lugar de 1.69, el mismo programa da $-2.8607\,E_\mathrm{h}$ . ¿Qué base es mejor, y por qué puede afirmarse?
¿Qué representa la energía del orbital virtual, $-0.1725\,E_\mathrm{h}$ ?

Parte IV — La unión de un protón.

¿Cuál es la energía de un protón desnudo? ¿Y la de $\ce{He + H+}$ muy alejados, en esta base?
Calcúlese la energía liberada por $\ce{He + H+ -> HeH+}$ a este nivel, en $\mathrm{kJ}/\mathrm{mol}$ .
La afinidad protónica medida del helio es $177.8\,\mathrm{kJ}/\mathrm{mol}$ . Coméntese.
La energía exacta del helio es menos la suma de sus energías de ionización, 24.5874 y $54.4178\,\mathrm{eV}$ . Calcúlese en $E_\mathrm{h}$ , así como el error del átomo de helio STO-3G.
¿Por qué un error de este tamaño puede dar aun así geometrías moleculares útiles?
Enúnciese el resultado: la energía total RHF/STO-3G de $\ce{HeH+}$ a $1.4632\,a_0$ .

Solución

Solución de Problema 3.1.

1. Cada orbital de tipo Slater de una base mínima se sustituye por una contracción fija de tres gaussianas ajustadas a él. 2. El núcleo de helio atrae más a sus electrones, de modo que su función $1s$ es más compacta (exponentes mayores). $6.3624/3.4253 = 1.857 = (1.69/1.24)^2$ . 3. Dos funciones de base, dos electrones, un orbital ocupado (y uno virtual). 4. $V_{\mathrm{nn}} = Z_{\mathrm{He}}Z_{\mathrm H}/R = 2/1.4632 = 1.3669\,E_\mathrm{h}$ . 5. Las dos funciones solapan mucho (54 %): los átomos están lo bastante cerca para enlazarse. 6. $h_{11}$ es la energía de un electrón en la función de He con los dos núcleos pero sin ningún otro electrón; el núcleo de He (carga 2) lo retiene con mucha más fuerza. 7. $(h_{11} - \varepsilon)(h_{22} - \varepsilon) - (h_{12} - \varepsilon S)^2 = 0$ : $0.71185\varepsilon^2 + 2.79292\varepsilon + 2.44969 = 0$ . 8. $\varepsilon = -2.600\,E_\mathrm{h}$ y $-1.324\,E_\mathrm{h}$ . 9. $\mathbf h$ ignora la repulsión entre los dos electrones; el orbital de núcleo es demasiado contraído y demasiado bajo. 10. $F_{\mu\nu} = h_{\mu\nu} + \sum_{\lambda\sigma}P_{\lambda\sigma}[(\mu\nu|\sigma\lambda) - \frac12(\mu\lambda|\sigma\nu)]$ , con $P_{\lambda\sigma} = 2C_{\lambda1}C_{\sigma1}$ . 11. La repulsión coulombiana de cada electrón por la nube de carga del otro (y, en general, el intercambio). 12. $\mathbf F$ depende de $\mathbf P$ , que depende de los orbitales obtenidos a partir de $\mathbf F$ : las ecuaciones no son lineales y se resuelven por aproximaciones sucesivas hasta la autoconsistencia. 13. $c_1^2 + c_2^2 + 2c_1c_2S = 0.7684 + 0.0410 + 0.1906 = 1.0000$ . 14. He: $1.5369 + 0.1906 = 1.727$ ; H: $0.0820 + 0.1906 = 0.273$ electrones. 15. Cargas: He $+0.27$ , H $+0.73$ . La carga positiva está sobre todo en el hidrógeno: $\ce{HeH+}$ es un protón unido a un átomo de helio, $\ce{He + H+}$ , como cabía esperar, ya que la energía de ionización de He (24.6 eV) supera con mucho la de H (13.6 eV). 16. $-\varepsilon_1 = 1.633\,E_\mathrm{h} = 44.4\,\mathrm{eV}$ . 17. $E_{\mathrm{el}} = -2.8418 - 1.3669 = -4.2087\,E_\mathrm{h}$ . 18. Tres ciclos alcanzan $-2.8418\,E_\mathrm{h}$ . Cada ciclo da un determinante, cuya energía es una cota superior de la energía de Hartree–Fock convergida (principio variacional), y las iteraciones la mejoran. 19. La base con $\zeta = 2.0925$ da la energía más baja, luego es la mejor para esta molécula (principio variacional): la función del helio es más compacta en $\ce{HeH+}$ que en el átomo libre, cuyo $\zeta = 1.69$ copia la base estándar. 20. El orbital antienlazante vacío $\sigma^*$ ; su energía aproximaría menos la afinidad electrónica de $\ce{HeH+}$ (mal, en una base tan pequeña). 21. Cero (ningún electrón); $-2.8078\,E_\mathrm{h}$ . 22. $-2.8078 - (-2.8418) = 0.0340\,E_\mathrm{h} = 89\,\mathrm{kJ}/\mathrm{mol}$ . 23. Alrededor de la mitad del valor medido: una base mínima no puede polarizar la densidad del helio hacia el protón (no hay funciones $p$ ), de modo que el enlace es demasiado débil. (Las correcciones de punto cero y a $298\,\mathrm{K}$ son pequeñas a su lado.) 24. $-(24.5874 + 54.4178)/27.211 = -2.9034\,E_\mathrm{h}$ ; el átomo STO-3G queda $0.0956\,E_\mathrm{h}$ ( $2.6\,\mathrm{eV}$ ) demasiado alto. 25. Las geometrías dependen de cómo cambia la energía con $R$ ; la mayor parte del error, concentrada en el núcleo, es casi la misma en todas las geometrías y se cancela. 26. $E(\ce{HeH+}) = -2.8418\,E_\mathrm{h}$ en RHF/STO-3G (base estándar), a $1.4632\,a_0$ (y $-2.8607\,E_\mathrm{h}$ con el clásico $\zeta_{\mathrm{He}} = 2.0925$ ).

Términos definidos en este capítulo

Ver los 852 términos del glosario

3.1 Born–Oppenheimer y la superficie de energía potencial

3.2 El principio variacional

3.3 Conjuntos de base

3.4 Teoría de Hartree–Fock

3.5 Teoría del funcional de la densidad, y lo que puede decir un cálculo

3.6 Ejercicios

3.7 Problema: HeH+^++, la primera molécula

Términos definidos en este capítulo

3.7 Problema: HeH $^+$ , la primera molécula