19.4: Más detalles sobre la implementación de métodos multiconfiguracionales

El método MCSCF

La optimización simultánea de los coeficientes LCAO-MO y CI realizada dentro de un cálculo MCSCF es una tarea bastante formidable. La energía variacional funcional es una función cuadrática de los coeficientes CI, por lo que se pueden expresar las condiciones estacionarias para estas variables en la forma secular:

\[ \sum\limits_J\text{H}_{I,J}C_J = \text{E C}_I . \]

Sin embargo, E es una función cuártica de los\(_{\nu ,i}\) coeficientes C porque cada elemento de la matriz\( \Phi_I \big| \text{ H } \big| \Phi_I \rangle \) involucra integrales de uno y dos electrones sobre el mos\(\phi_i\), y las integrales de dos electrones dependen cuartalmente de los\(C_{\nu ,i}\) coeficientes. Las condiciones estacionarias respecto a estos\(C_{\nu ,i}\) parámetros deben resolverse iterativamente debido a esta dependencia cuártica.

Es bien sabido que la minimización de una función (E) de varios parámetros no lineales (the\(C_{\nu ,i}\)) es una tarea difícil que puede sufrir una convergencia deficiente y puede localizar mínimos locales en lugar de globales. En una función de onda MCSCF que contiene muchos CSFs, la energía solo depende débilmente de los orbitales que están débilmente ocupados (es decir, los que aparecen en los MCS con\(C_I\) valores pequeños); en contraste, E es fuertemente dependiente de los\(C_{\nu ,i}\) coeficientes de aquellos orbitales que aparecen en los CSFs con mayores \(C_I\)valores. Por lo tanto, uno se enfrenta a minimizar una función de muchas variables (puede haber tantas\(C_{\nu ,i}\) como el cuadrado del número de funciones de base orbital) que depende fuertemente de varias de las variables y débilmente de muchas otras. Este es un trabajo muy difícil.

Por estas razones, en el método MCSCF, el número de CSFs generalmente se mantiene en un número pequeño a moderado (por ejemplo, de unos pocos a varios cientos) elegidos para describir correlaciones esenciales (es decir, cruces de configuración, disociación adecuada) e importantes correlaciones dinámicas (esas correlaciones de pares de electrones de angular , radial, izquierda-derecha, etc. naturaleza que surgen cuando están presentes orbitales 'virtuales' bajos). En una función de onda tan compacta, solo aparecen orbitales giratorios con ocupaciones razonablemente grandes (por ejemplo, caracterizadas por los elementos diagonales de la matriz de densidad de una partícula\(\gamma_{i,}\) j). Como resultado, la energía funcional se expresa en términos de variables de las que es fuertemente dependiente, en cuyo caso el proceso de optimización no lineal es menos probable que sea patológico.

Esta función de onda compacta de MCSCF está diseñada para proporcionar una buena descripción del conjunto de orbitales espín fuertemente ocupados y de las amplitudes de CI para los CSF en los que solo aparecen estos orbitales de espín. Por supuesto, no proporciona información sobre los orbitales de giro que no se utilizan para formar los CSFs en los que se basa el cálculo de MCSCF. Como resultado, la energía MCSCF es invariante a una transformación unitaria entre estos orbitales 'virtuales'.

Además de las referencias mencionadas anteriormente en ACP y MTC, los siguientes trabajos describen varios de los avances que se han hecho en el método MCSCF, especialmente con respecto a la potenciación de su tasa y rango de convergencia: E. Dalgaard y P. Jørgensen, J. Chem. Phys. 69, 3833 (1978); H. J. Aa. Jensen, P. Jørgensen, y H. ÅÅgren, J. Chem. Phys. 87, 457 (1987); B. H. Lengsfield, III y B. Liu, J. Chem. Phys. 75, 478 (1981).

El método de interacción de configuración

En el método de Interacción de Configuración (IC), se suele intentar realizar un tratamiento de alto nivel de correlación electrónica. Un conjunto de orbitales moleculares ortonormales se obtienen primero de un cálculo de SCF o MCSCF (generalmente involucra una lista pequeña a moderada de CSFs). Los coeficientes LCAO-MO de estos orbitales ya no se consideran como parámetros variacionales en el posterior cálculo de IC; solo los\(C_I\) coeficientes deben optimizarse aún más.

La función de onda CI

\[ \psi = \sum\limits_J C_J \Phi_J \]

se construye más comúnmente a partir de CSF\(\Phi_J\) que incluyen:

1. Todos los CSF en el SCF (en cuyo caso solo se incluye un solo CSF) o función de onda MCSCF que se utilizó para generar los orbitales moleculares\(\phi_i\). Este conjunto de CSF se conoce como que abarca el 'espacio de referencia' del cálculo de CI posterior, y la combinación particular de estos CSF utilizados en esta optimización orbital (es decir, la función de onda SCF o MCSCF) se llama la función de referencia.
2. Los CSFs que se generan mediante la realización de 'excitaciones' de nivel simples, dobles, triples, etc. (es decir, reemplazos orbitales) con respecto a los CFs de referencia. Las funciones de onda CI limitadas a incluir contribuciones a través de varios niveles de excitación (por ejemplo, simple, doble, etc.) se denotan S (excitada individualmente), D (doblemente), SD (individual y doblemente), SDT (individual, doblemente y triple), y así sucesivamente.

Los orbitales de los que se eliminan los electrones y aquellos en los que se excitan los electrones pueden restringirse para centrar la atención en las correlaciones entre ciertos orbitales. Por ejemplo, si se excluyen las excitaciones fuera de los electrones centrales, se calcula una energía total que no contiene correcciones de correlación para estos orbitales centrales. A menudo es posible limitar así la naturaleza de las excitaciones orbitales para enfocarse en las cantidades energéticas de interés (por ejemplo, la ruptura del enlace CC en etano requiere correlación de la\(\sigma_{CC}\) órbita, pero los orbitales del núcleo de carbono 1s y los orbitales del enlace CH pueden tratarse de una manera no correlacionada).

Claramente, el número de MCSF incluidos en el cálculo de CI puede ser muy superior al número considerado en los cálculos típicos de MCSCF; las funciones de onda de CI que incluyen 5,000 a 50,000 MCA se utilizan rutinariamente, y las funciones con uno a varios millones de MCSF están dentro del ámbito de la practicidad (véase, para ejemplo, J. Olsen, B. Roos, Poul Jørgensen, y H. J. Aa. Jensen, J. Chem. Phys. 89, 2185 (1988) y J. Olsen, P. Jørgensen, y J. Simons, Chem. Phys. Letras 169, 463 (1990)).

La necesidad de expansiones de CSF tan grandes no debería ser una sorpresa una vez que se considera que (i) cada par de electrones requiere al menos dos CSF (digamos que requiere P de ellos, en promedio, uno dominante y otros P-1 que están doblemente excitados) para formar pares orbitales polarizados, (ii) hay de los orden de N (N-1) /2 = X pares de electrones en un átomo o molécula que contiene N electrones, y (iii) que el número de términos en la función de onda CI escala como\(P^X\). Entonces, para una molécula H\(_2\) O que contiene diez electrones, habría\(^{55}\) términos P en la expansión del CI. Esto es 3.6 x10\(^{16}\) términos si P=2 y 1.7 x10\(^{26}\) términos si P=3. Sin lugar a dudas, se trata de una sobreestimación del número de CSF necesarios para describir la correlación electrónica en H\(_2\) O, pero demuestra cuán rápido puede crecer el número de CSF con el número de electrones en el sistema.

Las\(H_{I,J}\) matrices que surgen en los cálculos de CI se evalúan en términos de integrales de uno y dos electrones sobre los orbitales moleculares utilizando el equivalente de las reglas de Slater-Condon. Para los cálculos de CI grandes, la\(H_{I,J}\) matriz completa no es realmente evaluada y almacenada en la memoria de la computadora (o en su disco); más bien, se utilizan los llamados métodos 'IC directo' (ver el artículo de Roos y Siegbahn en MTC) para calcular y sumar inmediatamente las contribuciones a la suma\( \sum\limits_J \text{H}_{I,J}C_I \) en términos de integrales, elementos de matriz de densidad y valores aproximados de las\(C_J\) amplitudes. Métodos iterativos (ver, por ejemplo, E. R. Davidson, J. Comput. Phys. 17, 87 (1975)), en los que los valores aproximados para los\(C_J\) coeficientes y la energía E se refinan a través de la aplicación secuencial de\( \sum\limits_I\text{H}_{I,J} \) a la estimación anterior del\(C_J\) vector, se emplean para resolver estos problemas de valores propios de gran matriz de CI.

El método MPPT/MBPT

En el método MPPT/MBPT, una vez que se elige el CSF de referencia y se determinan los orbitales SCF que pertenecen a este CSF, la función de onda\(\psi\) y la energía E se determinan de manera orden por orden. Esta es una de las principales fortalezas de la técnica MPPT/MBPT; no requiere que se realicen elecciones adicionales (potencialmente arbitrarias) una vez que se especifican el conjunto de bases y la configuración dominante (SCF). A diferencia de los tratamientos MCSCF e IC, no es necesario tomar decisiones de los CSFs para incluirlos o excluirlos de ellos\(\psi\). Las ecuaciones de perturbación MPPT/MBPT determinan qué CSF deben incluirse a través de cualquier orden en particular.

Por ejemplo, la corrección de función de onda de primer orden\(\psi^1 \text {(i.e., } \psi = \Phi + \psi^1\) a través de primer orden) viene dada por:

\[ \begin{align} \psi^1 &= -\sum\limits_{i<j,m>n} \dfrac{\langle \Phi_{i,j}^{m,n}\big| \text{ H - H}^0 \big| \Phi \rangle}{\epsilon_m - \epsilon_i + \epsilon_n - \epsilon_j} \bigg| \Phi_{i,j}^{m,n} \rangle \\ &= - \sum\limits_{i<j,M<n} \dfrac{\langle i,j \big| g \big| m,n \rangle - \langle i,j \big| g \big| n,m \rangle}{\epsilon_m - \epsilon_i + \epsilon_n - \epsilon_j}\bigg| \Phi_{i,j}^{m,n} \rangle \end{align} \]

donde se denotan las energías orbitales SCF\(\epsilon_k \text{ and } \Phi_{i,j}^{m,n}\) representa un CSF que está doblemente excitado en relación con\ Phi. Así, solo los CSF doblemente excitados contribuyen a la función de onda de primer orden; como resultado, la energía E se da a través de segundo orden como:

\[ \begin{align} E &= \langle \big| \text{H}^0 \big| \Phi \rangle + \langle \Phi \Big| \text{ H - H}^0 \big| \Phi \rangle + \langle \Phi \big| \text{ H - H}^0 \big| \psi^1 \rangle \\ &= \langle \Phi \big| \text{ H } \big| \Phi \rangle - \sum\limits_{i<j,m<n} \dfrac{\big| \langle \Phi_{i,j}^{m,n} \big| \text{ H - H}^0 \big| \Phi \rangle \big|^2}{\epsilon_m - \epsilon_i + \epsilon_n -\epsilon_j} \\ &= \text{E}_{\text{SCF}} - \sum\limits_{i<j,m<n} \dfrac{\big| \langle i,j \big| g \big| m,n \rangle - \langle i,j \big| g \big| n,m \rangle \big|^2}{\epsilon_m - \epsilon_i + \epsilon_n - \epsilon_j} \\ &= \text{E}^0 + \text{E}^1 + \text{E}^2. \end{align} \]

Estas contribuciones se han expresado, utilizando las reglas SC, en términos de las integrales de dos electrones que\( \langle \text{i,j} \big| g \big|\text{m,n} \rangle \) acoplan los orbitales espín excitados a los orbitales espín-orbitales de los cuales se excitaron los electrones, así como las diferencias de energía orbitales que\( [ \epsilon_m - \epsilon_i + \epsilon_n - \epsilon_j ] \) acompañan a tales excitaciones. De esta forma, queda claro que las principales contribuciones a la energía de correlación del par de orbitales ocupados\(\phi_i \phi_j\) se realizan mediante excitaciones dobles en orbitales virtuales\(\phi_m \phi_n\) que tienen acoplamiento grande (i.. e.,\( \langle i,j \big| \text{ g }\big| m,n \rangle \) integrales grandes) y pequeñas brechas de energía orbital,\( [ \epsilon_m - \epsilon_i + \epsilon_n - \epsilon_j ] \).

En orden superior, las correcciones a la función de onda y a la energía, también ingresan contribuciones de los MCS que están solos, triplicados, etc. relativos a\(\Phi\) aparecer, y contribuciones adicionales de los MCS doblemente excitados. Es relativamente común llevar cálculos MPPT/MBPT (ver las referencias dadas anteriormente en el Capítulo 19.I.3 donde se documentan las contribuciones de los grupos Pople y Bartlett al desarrollo de MPPT/MBPT) hasta el tercer orden en la energía (cuya evaluación se puede demostrar que solo requiere\(\psi^0 \text{ and } \psi^1\)). Toda la serie de programas GAUSSIAN-8X, que han sido utilizados en miles de importantes estudios químicos, calculan E a través de tercer orden de esta manera.

Además de ser extensivo de tamaño y no requerir uno para especificar entrada más allá del conjunto base y el CSF dominante, el enfoque MPPT/MBPT es capaz de incluir el efecto de todos los CSFs (que contribuyen a cualquier orden dado) sin tener que encontrar ningún valor propio de una matriz. Esta es una ventaja importante porque la determinación del valor propio de la matriz, que es necesaria en los cálculos de MCSCF y CI, requiere tiempo de computadora en proporción a la tercera potencia de la dimensión de la\(\text{H}_{I,J}\) matriz. A pesar de todas estas ventajas, es importante recordar las principales desventajas del enfoque MPPT/MBPT; su energía no es un límite superior a la energía verdadera y puede no ser capaz de tratar casos para los que dos o más CSF tienen amplitudes iguales o casi iguales porque obtiene las amplitudes de todos menos el CSF dominante a partir de fórmulas de teoría de perturbación que asumen que la perturbación es 'pequeña'.

El Método de Cluster Acoplado

La implementación del método Coupled-Cluster comienza tanto como en el caso MPPT/MBPT; se selecciona un CSF de referencia que se utiliza en el proceso SCF para generar un conjunto de espin-orbitales que se utilizarán en el posterior cálculo correlacionado. El conjunto de ecuaciones de trabajo de la técnica Coupled-Cluster se puede escribir explícitamente introduciendo la forma del llamado operador cluster T,

\[ \text{T} = \sum\limits_{i,m}t_i^m m^+ i + \sum\limits_{i,j,m,n}t_{i,j}^{m,n}m^+ n^+ \text{ j i + ...,} \]

donde la combinación de operadores\(m^+ i\) denota la creación de un electrón en orbital de espín virtual\(\phi_m\) y la eliminación de un electrón de espín-orbital ocupado\(\phi_i\) para generar una sola excitación. Por lo tanto, la operación\(m^+ n^+\) j i representa una doble excitación de\(\phi_i \phi_j\) a\(\phi_m \phi_n\). Expresar el operador de clúster T en términos de las amplitudes,\(t_i^m \text{, } t_{i,j}^{m,n} \text{, }\) etc. para CSFs individuales, doblemente, etc. excitados, y expandir los operadores exponenciales en

\[e^{-T} \text{ H } e^{T}\]

se obtiene:

\[ \langle \Phi_i^m \big| \text{ H + [H,T] + }\dfrac{1}{2} \text{[[H,T],T] + }\dfrac{1}{6}\text{[[[H,T],T],T] + }\dfrac{1}{24} \text{[[[[H,T],T],T],T] }\big| \Phi \rangle = 0; \]

\[ \langle \Phi_{i,j}^{m,n} \big| \text{ H + [H,T] + }\dfrac{1}{2} \text{[[H,T],T] + }\dfrac{1}{6}\text{[[[H,T],T],T] + }\dfrac{1}{24} \text{[[[[H,T],T],T],T] }\big| \Phi \rangle = 0; \]

\[ \langle \Phi_{i,j,k}^{m,n,p} \big| \text{ H + [H,T] + }\dfrac{1}{2} \text{[[H,T],T] + }\dfrac{1}{6}\text{[[[H,T],T],T] + }\dfrac{1}{24} \text{[[[[H,T],T],T],T] }\big| \Phi \rangle = 0; \]

y así sucesivamente para los CSF excitados de orden superior. Se puede demostrar, debido a la naturaleza de operador de uno y dos electrones de H, que la expansión de los operadores exponenciales trunca exactamente a la cuarta potencia; es decir términos como [[[[H, T], T], T], T], T], T], T] y los conmutadores superiores desaparecen de manera idéntica (esto se demuestra en el Capítulo 4 de la Segunda Métodos Basados en Cuantización en Química Cuántica, P. Jørgensen y J. Simons, Academic Press, Nueva York (1981).

Figura 19.3.1: Representación esquemática de la clasificación basada en expansión de clústeres. La correlación completa está compuesta por singletes, dobletes, trillizos y correlaciones de orden superior, todas ellas definidas de manera única por el enfoque de expansión de racimos. Cada esfera azul corresponde a un operador de partículas y círculos/elipses amarillos a correlaciones. El número de esferas dentro de una correlación identifica el número de conglomerados. Imagen utilizada con permiso (CC SA-BY-3.0; Christoph N. Böttge, “Phonon-assistierte Lasertätigkeit in Mikroresonatoren”).

Como resultado, las ecuaciones exactas de Cluster Acoplado son ecuaciones cuárticas para las amplitudes\(t_i^m \text{, } t_{i,j}^{m,n} \text{, }\) etc. Si bien es una tarea bastante formidable evaluar todos los elementos de la matriz del conmutador que aparecen en las ecuaciones de Cluster Acoplado anteriores, puede ser y se ha hecho (las referencias dadas anteriormente a Purvis y Bartlett son especialmente relevantes en este contexto). El resultado es expresar cada uno de esos elementos de la matriz, a través de las reglas de Slater-Condon, en términos de integrales de uno y dos electrones sobre los orbitales espín utilizados para determinar\(\Phi\), incluyendo aquellos en\(\Phi\) sí mismos y los orbitales 'virtuales' que no están en\(\Phi\).

En general, estas ecuaciones cuárticas deben resolverse de manera iterativa y son susceptibles de dificultades de convergencia similares a las que surgen en los cálculos de tipo MCSCF. En cualquier proceso iterativo de este tipo, es importante comenzar con una aproximación (a las t amplitudes, en este caso) que esté razonablemente cerca del resultado convergente final. Tal aproximación a menudo se logra, por ejemplo, al descuidar todos los términos que son no lineales en las t amplitudes (porque se supone que estas amplitudes son menores que la unidad en magnitud). Esto lleva, para las ecuaciones de trabajo de Cluster Acoplado obtenidas al proyectar sobre los CSFs doblemente excitados, a:

\[ \langle \text{ i j } \big| \text{ g } \big| \text{ m,n } \rangle ' + \left[ \epsilon_m -\epsilon_i +\epsilon_n -\epsilon_j \right] t_{i,j}^{m,n} + \sum\limits_{i',j',m',n'} \langle \Phi_{i,j}^{m,n} \big| \text{ H - H}^0 \big| \Phi_{i',j'}^{m',n'} \Phi_{i',j'}^{m',n'} \rangle t_{i',j'}^{m',n'} = 0, \]

donde\( \frac{\langle \text{ i,j }\big| \text{ g }\big| \text{ m,n }\rangle '}{\epsilon_m -\epsilon_i +\epsilon_n -\epsilon_j} \) se usa la notación para denotar la diferencia integral de dos electrones\(\rangle \text{ i,j } \big| \text{ g } \big| \text{ m,n } \rangle \text{ - } \langle \text{ i,j } \big| \text{ g } \big| \text{ n,m } \rangle\). Si, además, se ignoran los factores que acoplan diferentes CSF doblemente excitados (es decir, la suma sobre i', j', m', n'), las ecuaciones para las amplitudes t se reducen a las ecuaciones para las amplitudes de CSF de la función de onda MPPT/MBPT de primer orden:

\[ t_{i,j}^{m,n} = - \dfrac{\langle \text{ i,j }\big| \text{ g } \big| \text{ m,n } \rangle '}{\epsilon_m -\epsilon_i +\epsilon_n -\epsilon_j} \]

Como Bartlett y Pople han demostrado, existe, de hecho, una estrecha relación entre los métodos MPPT/MBPT y Cluster Acoplado cuando las ecuaciones de Cluster Acoplado se resuelven iterativamente comenzando con tal 'conjetura' inicial similar a MPPT/MBPT para estas amplitudes de doble excitación.

El método Coupled-Cluster, tal como se presenta aquí, adolece de los mismos inconvenientes que el enfoque MPPT/MBPT; su energía no es un límite superior y puede no ser capaz de describir con precisión las funciones de onda que tienen dos o más CSF con amplitud aproximadamente igual. Además, la solución de las ecuaciones de Cluster Acoplado no lineal puede ser difícil y lentamente (si acaso) convergente. Tiene las mismas ventajas que el método MPPT/MBPT; su energía es de tamaño extenso, no requiere una solución de gran valor propio de matriz, y su energía y función de onda se determinan una vez que se especifica la base y el CSF dominante.

Métodos Funcionales de Densidad

Estos enfoques proporcionan alternativas a las herramientas convencionales de la química cuántica. Los métodos CI, MCSCF, MPPT/MBPT y CC van más allá de la imagen de configuración única al agregar a la función de onda más configuraciones cuyas amplitudes determinan cada una a su manera. Esto puede llevar a un número muy grande de CSF en la función de onda correlacionada, y, como resultado, una necesidad de recursos informáticos extraordinarios.

Los enfoques funcionales de densidad son diferentes. Aquí se resuelve un conjunto de ecuaciones orbitales

\[ \left[ -\dfrac{\hbar^2}{2m_e}\nabla^2 -\sum\limits_A \dfrac{Z_Ae^2}{|\textbf{r} - \textbf{R}_A}| + \int\dfrac{\rho (\textbf{r}'e^2)}{|\textbf{r}-\textbf{r}'|}\textbf{dr}' + \text{U}(\textbf{r}) \right] \phi_i = \epsilon_i\phi_i \]

en el que los orbitales {\(\phi_i\)} 'sienten' potenciales debido a los centros nucleares (que tienen cargas Z\(_{\text{A}}\)), interacción coulómbica con la densidad electrónica total\(\rho\) (r '), y un llamado potencial de correlación de intercambio denotado U (r'). El estado electrónico particular para el que se está realizando el cálculo se especifica formando una densidad correspondiente\(\rho\) (r '). Antes de ir más allá en la descripción de cómo se llevan a cabo los cálculos de DFT, examinemos los orígenes que subyacen a esta teoría.

El llamado teorema de Hohenberg-Kohn establece que la densidad electrónica del estado fundamental\(\rho\) (r) que describe un sistema de electrones N determina de manera única el potencial V (r) en el Hamiltoniano

\[ \text{H } = \sum\limits_J \dfrac{\hbar^2}{2m_e}\nabla_j^2 + \text{V(}\textbf{r}_j\text{)} + \dfrac{e^2}{2} \sum\limits_{k \ne j} \dfrac{1}{r_{j,k}} , \]

y, debido a que H determina la energía del estado del suelo y la función de onda del sistema, la densidad del estado del suelo\(\rho (\textbf{r})\) determina las propiedades del estado del suelo del sistema. La prueba de este teorema procede de la siguiente manera:

• \(\rho (\textbf{r}) \text{ determines N because } \int \rho (\textbf{r}) d^3r = \text{N} .\)
• Supongamos que hay dos potenciales distintos (aparte de una constante aditiva que simplemente desplaza el cero de la energía total) V (r) y V' (r) que, cuando se usan en H y H', respectivamente, para resolver para un estado fundamental producen\[ \text{E}_0 \text{, } \psi \text{(}\textbf{r}\text{) and E}_0 ’ \text{, } \psi \text{’(}\textbf{r}\text{) }\] que tienen la misma densidad de un electrón:\[ \int \big| \psi \big|^2 \text{dr}_2 \text{ dr}_3 \text{ ... dr}_{\text{N}} = \rho (\textbf{r}) = \int \big|\psi '\big| \text{dr}_2 \text{ dr}_3 \text{ ... dr}_{\text{N}} \]
• Si pensamos en la función de\(\psi\) onda variacional de prueba para la H hamiltoniana, sabemos que\[ \text{E}_0 < \langle \psi ' \big| \text{H}\big| \psi '\rangle = \langle \psi '\big| \text{H'} \big| \psi ' \rangle + \int \rho(\textbf{r})\left[V(\textbf{r}) - V'(\textbf{r}) \right]\text{d}^3\text{r} = \text{E}_o ' + \int \rho (\textbf{r}) \left[ V(\text{r}) - V'(\textbf{r}) \right]\text{d}^3\textbf{r} \label{B3} \]
• De igual manera, tomando\(\psi\) como función de prueba para el H' hamiltoniano, se encuentra que\[ E_0 ' \langle E_0 + \int \rho(\textbf{r}) \left[ V'(\textbf{r}) - V(\textbf{r}) \right] \text{d}^3\textbf{r}. \label{B4}\]
• Sumando las Ecuaciones\ ref {B3} y\ ref {B4} da\[ \text{E}_0 \text{ + } \text{E}_0 \text{ + } \text{E}_0 ', \] y una clara contradicción.

Por lo tanto, no puede haber dos potenciales distintos V y V' que den el mismo estado base\(\rho (\textbf{r}\)). Entonces, la densidad del estado fundamental\(\rho (\textbf{r}\)) determina de manera única N y V, y por lo tanto H, y por lo tanto\(\psi \text{ and E}_0. \text{ Furthermore, because } \psi\) determina todas las propiedades del estado fundamental\(\rho (\textbf{r})\), entonces, en principio, determina todas esas propiedades. Esto significa que incluso la energía cinética y la energía de interacción electrón-electrón del estado fundamental están determinadas por\(\rho (\textbf{r}\)). Es fácil ver que\( \int \rho (\textbf{r}) \text{V(}\textbf{r}\text{)d}^3\textbf{r} = \text{V[}\rho ]\) da el valor promedio de la interacción nuclear del electrón (más cualquier potencial aditivo adicional de un electrón) en términos de la densidad del estado fundamental\(\rho (\textbf{r})\), pero cómo se expresan la energía cinética T [\(\rho\)] y la interacción electrón-electrón V\(_{ee}\text{[}\rho\)] energía en términos de\(\rho\)?

La principal dificultad con DFT es que el teorema de Hohenberg-Kohn muestra que los valores de estado fundamental de T\(_{ee}\), V, V, etc. son todos funcionales únicos del estado fundamental\(\rho\) (es decir, que pueden, en principio, determinarse una vez que\(\rho\) se da), pero no nos dice cuáles son estas relaciones funcionales.

Para ver cómo podría tener sentido que una propiedad como la energía cinética, cuyo operador\( -\frac{\hbar^2}{2m_e}\nabla^2 \) involucra derivados, pueda relacionarse con la densidad de electrones, considere un sistema simple de N electrones no interactuantes que se mueven en un potencial de “caja” cúbica tridimensional. Se sabe que los estados energéticos de tales electrones son

\[ E = \left( \dfrac{h^2}{2m_e L^2} \right) \left( n_x^2 + n_y^2 + n_z^2 \right), \]

donde L es la longitud de la caja a lo largo de los tres ejes, y\(n_x \text{, } n_y \text{, and } n_z\) son los números cuánticos que describen el estado. Podemos ver\(n_x^2 \text{ + } n_y^2\text{ + } n_z^2 = R^2\) como definiendo el radio cuadrado de una esfera en tres dimensiones, y nos damos cuenta de que la densidad de los estados cuánticos en este espacio es un estado por unidad de volumen en el\(n_x \text{, } n_y \text{, } n_z\) espacio. Debido a que\(n_x \text{, } n_y \text{, and } n_z\) deben ser enteros positivos, el volumen que cubre todos los estados con energía menor o igual a una energía especificada E =\(\left(\frac{h^2}{2m_eL^2}\right)R^2 \text{ is } \frac{1}{8}\) el volumen de la esfera de radio R:

\[ \Phi (\text{E}) = \dfrac{1}{8}\left( \dfrac{4\pi}{3}\right)R^3 = \left( \dfrac{\pi}{6}\right) \sqrt{9m_e E\dfrac{L^2}{h^2}}^3 . \]

Dado que hay un estado por unidad de dicho volumen,\(\Phi\) (E) es también el número de estados con energía menor o igual a E, y se denomina densidad integrada de estados. El número de estados g (E) dE con energía entre E y E+de, la densidad de estados, es la derivada de\(\Phi\):

\[ g(E) = \dfrac{d\Phi}{dE} = \dfrac{\pi}{4}\sqrt{8m_e\dfrac{L^2}{h^2}}^3 \sqrt{E}. \]

Si calculamos la energía total para N electrones, teniendo los estados energías hasta la llamada energía Fermi (es decir, la energía del HOMO orbital molecular ocupado más alto) doblemente ocupada, obtenemos la energía del estado fundamental:

\[ E_0 = 2\int\limits_0^{E_F}g(E)EdE = \dfrac{8\pi}{5}\sqrt{\dfrac{2m_e}{h^2}}^3L^3\sqrt{E_F}^5 . \]

El número total de electrones N se puede expresar como

\[ N = 2\int\limits_0^{E_F}g(E)dE = \dfrac{8\pi}{3}\sqrt{\dfrac{2m_e}{h^2}}^3 L^3 \sqrt{E_F}^3 , \]

que se puede resolver\(E_F\) en términos de N para luego expresar\(E_0\) en términos de N en lugar de\(E_F\):

\[ E_0 = \dfrac{3h^2}{10m_e}\sqrt[3]{\dfrac{3}{8\pi}}^2L^3 \sqrt[3]{\dfrac{N}{L}}^5 . \]

Esto da la energía total, que también es la energía cinética en este caso debido a que la energía potencial es cero dentro de la “caja”, en términos de la densidad electrónica\(\rho\) (x, y, z) = (\(\frac{N}{L^3}\)). Por lo tanto, puede ser plausible expresar energías cinéticas en términos de densidades de electrones\(\rho(\textbf{r}\)), pero de ninguna manera está claro cómo hacerlo para átomos y moléculas “reales” con interacciones electrón-nuclear y electrón-electrón operativas.

En uno de los primeros modelos de DFT, la teoría de Thomas-Fermi, la energía cinética de un átomo o molécula se aproxima utilizando el tipo de tratamiento anterior a nivel “local”. Es decir, para cada elemento de volumen en el espacio r, se asume que la expresión dada anteriormente es válida, y luego se integra sobre todo r para calcular la energía cinética total:

\[ T_{TF}[\rho ] = \int \dfrac{3h^2}{10m_e}\sqrt[3]{\dfrac{3}{8\pi}}^2 \sqrt[3]{\rho(\textbf{r})}^5 \text{d}^3\textbf{r} = C_F\int \sqrt[3]{\rho (\textbf{r})}^5 \text{d}^3\textbf{r} , \]

donde la última igualdad simplemente define la\(_F\) constante C (que es 2.8712 en unidades atómicas). Ignorando las contribuciones de correlación e intercambio a la energía total, esta T se combina con las energías de potencial electrón-nuclear V y Coulombic electrón-electrón para dar la energía total Thomas-Fermi:

\[ E_{0,TF}[\rho ] = C_F \int \sqrt[3]{\rho (\textbf{r})}^5d^3\textbf{r} + \int V(\textbf{r})\rho (\textbf{r}) d^3\textbf{r} + \dfrac{e^2}{2}\int \rho (\textbf{r})\dfrac{\rho (\textbf{r}')}{|\textbf{r}-\textbf{r}' |}d^3\textbf{r} d^3\textbf{r}' , \]

Esta expresión es un ejemplo de cómo\(_0\) se da E como aproximación funcional de densidad local (LDA). El término local significa que la energía se da como un funcional (es decir, una función de\(\rho\)) que depende solo de\(\rho(\textbf{r})\) en puntos en el espacio pero no de\(\rho(\textbf{r})\) en más de un punto en el espacio.

Desafortunadamente, la energía funcional Thomas-Fermi no produce resultados que sean de precisión suficientemente alta para ser de gran utilidad en química. Lo que falta en esta teoría son a. la energía de intercambio y b. la energía de correlación; además, la energía cinética se trata sólo de la manera aproximada descrita.

En el libro de Parr y Yang, se muestra cómo Dirac fue capaz de abordar la energía de intercambio para el 'gas electrónico uniforme' (N Coulomb interactuando electrones moviéndose en una carga de fondo positiva uniforme cuya magnitud equilibra la carga de los N electrones). Si la expresión exacta para la energía de intercambio del gas electrónico uniforme se aplica a nivel local, se obtiene la aproximación de densidad local Dirac comúnmente utilizada a la energía de intercambio:

\[ \text{E}_{\text{ex, Dirac}}[\rho ] = - C_x \int \sqrt[3]{\rho (\textbf{r})}^4 d^3\textbf{r} , \]

con\(C_x = \frac{3}{4}|sqrt[3]{\frac{3}{\pi}} = 0.7386\) en unidades atómicas. Al agregar esta energía de intercambio a la energía total E de Thomas-Fermi, se\(_{0,TF} [\rho]\) obtiene la llamada energía Thomas-Fermi-Dirac (TFD) funcional.

Debido a que las densidades de electrones varían bastante fuertemente espacialmente cerca de los núcleos,\(\rho\text{] and E}_{ex.Dirac}\) se necesitan correcciones a las aproximaciones anteriores a T [. Una de las llamadas aproximaciones corregidas por gradiente más comúnmente utilizadas es la inventada por Becke, y conocida como la función de intercambio Becke88:

\[ E_{ex}(Becke88) = E_{ex,Dirac}[\rho ] -\gamma \int x^2\sqrt[3]{\rho}^4\left[ \dfrac{1}{1 + 6\gamma\text{ x }csch(x)} \right] \textbf{dr,} \]

donde x =\(\frac{1}{\sqrt[3]{\rho}^4} \big|\nabla \rho \big| \text{, and } \gamma\) es un parámetro elegido para que la energía de intercambio anterior pueda reproducir mejor las energías de intercambio conocidas de estados electrónicos específicos de los átomos de gas inerte (Becke encuentra igual\(\gamma\) a 0.0042). Una corrección de gradiente común a la T [\(\rho\)] anterior se llama corrección Weizsacker y viene dada por

\[ \delta\text{T}_{\text{Weizsacker}} = \dfrac{1}{72}\dfrac{\hbar}{m_e}\int \dfrac{\big| \nabla \rho (\textbf{r}) \big|^2}{\rho (\textbf{r})} d\textbf{r}. \]

Aunque la discusión anterior sugiere cómo se podría calcular la energía del estado fundamental una vez que\(\rho (\textbf{r})\) se da la densidad del estado fundamental, todavía se necesita saber cómo obtener\(\rho\). Kohn y Sham (KS) introdujeron un conjunto de los llamados orbitales KS obedeciendo la siguiente ecuación:

\[ \left[ -\dfrac{1}{2}\nabla^2 + V(\textbf{r}) + \dfrac{e^2}{2}\int \dfrac{\rho(\textbf{r}') }{\big|\textbf{r}-\textbf{r}'\big|} \textbf{dr}' + U_{xc}(\textbf{r}) \right] \phi_j = \epsilon_j \phi_j , \]

donde el llamado potencial de correlación de intercambio\(U_{xc} (\textbf{r}) = \delta\text{E}_{xc}\dfrac{[\rho ]}{\delta\rho (\textbf{r})}\) podría obtenerse por diferenciación funcional si se\(\text{E}_{xc}[\rho ]\) conociera la energía de intercambio-correlación funcional. KS también mostró que los orbitales KS {\(\phi_j\)} podrían ser utilizados para calcular la densidad\(\rho\) simplemente sumando las densidades orbitales multiplicadas por ocupaciones orbitales n\(_j\):

\[ \rho (\textbf{r} ) = \sum\limits_j n_j \big| \phi_j (\textbf{r})\big|^2 . \]

(aquí\(n_j\) =0,1, o 2 es el número de ocupación de la órbita\(\phi_j\) en el estado en estudio) y que la energía cinética debe calcularse como

\[ \text{T} = \sum\limits_j n_j \langle \phi_j (\textbf{r}) \big| -\dfrac{1}{2} \nabla^2 \big| \phi_j (\textbf{r}) \rangle . \]

Las mismas investigaciones del idealizado 'gas electrónico uniforme' que identificaron el intercambio funcional Dirac, encontraron que la energía de correlación (por electrón) también podría escribirse exactamente como una función de la densidad electrónica\(\rho\) del sistema, pero sólo en dos casos limitantes- el alto- límite de densidad (grande\(\rho\)) y el límite de baja densidad. Todavía no existe una expresión exacta para la energía de correlación incluso para el gas electrónico uniforme que es válido a valores arbitrarios de\(\rho\). Por lo tanto, se ha dedicado mucho trabajo a crear fórmulas de interpolación eficientes y precisas que conecten las expresiones uniformes de gas electrónico de baja y alta densidad. Una de esas expresiones es

\[ \text{E}_c [\rho ] = \int \rho (\textbf{r})\epsilon_c(\rho )\text{d}\textbf{r} , \]

donde

\[ \epsilon_c(\rho ) = \dfrac{A}{2}\left[ ln\left( \dfrac{x}{X} \right) + \dfrac{2b}{Q tan^{-1}\left( \dfrac{Q}{2x + b} \right)} -\dfrac{bx_0}{X_0} \left[ ln\left( \dfrac{(x-x_0)^2}{X} + \dfrac{2(b+x_0)}{Q tan^{-1}\left( \dfrac{Q}{2x+b} \right)} \right) \right] \right] \]

es la energía de correlación por electrón. Aquí x =\(\sqrt{r_s} \text{, } X=x^2 +bx+c \text{, } X_0 = x_0^2 + bx_0 + c \text{ and } Q = \sqrt{(4c - b^2)} \text{, } A = 0.0621814\text{, } x_0\) = -0.409286, b = 13.0720, y c = 42.7198. El parámetro\(r_s \text{ is how the density } \rho \text{ enters since } \frac{4}{3} \pi r_s^3\) es igual a\(\frac{1}{\rho}; \text{ that is, } r_s\) es el radio de una esfera cuyo volumen es el volumen efectivo ocupado por un electrón. Una aproximación razonable al completo\(\text{E}_{xc}[\rho ]\) contendría el intercambio funcional Dirac (y quizás corregido por gradiente) más lo anterior\(\text{E}_C[\rho ]\), pero hay muchas aproximaciones alternativas a la función de energía de intercambio-correlación. Actualmente, muchos trabajadores están haciendo todo lo posible para “cocinar” los funcionales para la correlación e intercambiar energías, pero todavía nadie ha inventado los funcionales que son tan confiables que la mayoría de los trabajadores acuerdan usarlos.

Para resumir, al implementar cualquier DFT, se suele proceder de la siguiente manera:

1. Se elige una base orbital atómica en términos de la cual se van a expandir los orbitales KS.
2. Se hace cierta suposición inicial para los coeficientes de expansión LCAO-KS\(\text{C}_{jj,a}: \phi_j = \sum\limits_a \text{C}_{j,a}\chi_a\).
3. La densidad se calcula como\(\rho (\textbf{r}) = \sum\limits_j n_j \big|\phi_j (\textbf{r})\big|^2 \). A menudo,\(\rho (\textbf{r})\) se expande en una base orbital atómica, que no necesita ser la misma que la base utilizada para el\(\phi_j\), y los coeficientes de expansión de\(\rho\) se calculan en términos de los de la\(\phi_j\). También es común utilizar una base orbital atómica para expandirse\(\sqrt[3]{\rho (\textbf{r})} \text{ which, together with } \rho\), es necesario evaluar la contribución de la función de correlación de intercambio a E\(_0\).
4. La densidad de la iteración actual se utiliza en las ecuaciones KS para determinar el hamiltoniano\(\left[-\frac{1}{2}\nabla^2 \text{ + V}(\textbf{r}) + \frac{e^2}{2} \int \frac{\rho (\textbf{r}’)}{|\textbf{r}-\textbf{r}’|} \text{ d}\textbf{r}’ + \text{ U}_{xc}(\textbf{r}) \right]\) cuyas “nuevas” funciones propias {\(\phi_j\)} y valores propios {\(\epsilon_j\)} se encuentran resolviendo las ecuaciones KS.
5. Estas nuevas\(\phi_j\) se utilizan para calcular una nueva densidad, la cual, a su vez, se utiliza para resolver un nuevo conjunto de ecuaciones KS. Este proceso se continúa hasta que se alcanza la convergencia (es decir, hasta que los\(\phi_j\) utilizados para determinar las iteraciones actuales\(\rho\) son los mismos\(\phi_j\) que surgen como soluciones en la siguiente iteración.
6. Una vez que\(\rho (\textbf{r})\) se determina la convergente, la energía se puede calcular usando la expresión anterior
   
   \[ \text{E} [\rho ] = \sum\limits_j n_j \langle \phi_j (\textbf{r})\big| -\dfrac{1}{2} \nabla^2 \big| \phi_j (\textbf{r})\rangle + \int V(\textbf{r}) \rho (\textbf{r}) \text{d}\textbf{r} + \dfrac{e^2}{2} \int \rho (\textbf{r})\dfrac{ \rho (\textbf{r}')}{\big|\textbf{r}-\textbf{r}'\big|}\text{d}\textbf{r }\text{d}\textbf{r}' + \text{E}_{xc}[\rho ]. \]

Al cerrar esta sección, hay que recalcar una vez más que esta zona se encuentra actualmente en crecimiento explosivo y mucho escrutinio. Como resultado, es casi seguro que muchos de los funcionales específicos discutidos anteriormente serán reemplazados en un futuro próximo por versiones mejoradas y más rigurosamente justificadas. También es probable que las extensiones de DFT a estados excitados (muchos trabajadores lo persiguen activamente) se coloquen sobre un terreno más sólido y se hagan aplicables a los sistemas moleculares. Debido a que el esfuerzo computacional involucrado en estos enfoques escala mucho menos fuertemente con el tamaño del conjunto base que para los métodos convencionales (SCF, MCSCF, CI, etc.), los métodos funcionales de densidad ofrecen una gran promesa y es probable que contribuyan mucho a la química cuántica en la próxima década.

Colaboradores y Atribuciones

• Jack Simons (Henry Eyring Scientist and Professor of Chemistry, U. Utah) Telluride Schools on Theoretical Chemistry and  Jeff A. Nichols (Oak Ridge National Laboratory)