El método de Hartree-Fock (HF) es una forma aproximada de las ecuaciones de mecánica cuántica para fermiones, utilizada en física y química (donde también se conoce como método de campo autoconsistente). Esto se debe a que sus ecuaciones, basadas en orbitales de una partícula, son más accesibles computacionalmente que los métodos basados en funciones de onda de muchas partículas.

La aproximación de Hartree-Fock es el equivalente, en física computacional, a la aproximación de orbitales moleculares, de enorme utilidad conceptual para los físicos. Este esquema de cálculo es un procedimiento iterativo para calcular la mejor solución monodeterminantal a la ecuación de Schrödinger independiente del tiempo, para moléculas aisladas, tanto en su estado fundamental como en estado excitados. La interacción de un único electrón en un problema de muchos cuerpos con el resto de los electrones del sistema se aproxima promediándolo como una interacción entre dos cuerpos (tras aplicar la aproximación de Born-Oppenheimer). De esta forma, se puede obtener una aproximación a la energía total de la molécula. Como consecuencia, calcula la energía de intercambio de forma exacta, pero no tiene en absoluto en cuenta el efecto de la correlación electrónica.

La base del método de Hartree-Fock es suponer que la función de onda de muchos cuerpos es un determinante de Slater de orbitales de una partícula. Esto garantiza la antisimetría de la función de onda y considera la energía de intercambio. Sin embargo, no considera efectos de correlación que no son despreciables. A partir de esta suposición, se puede aplicar el principio variacional de mecánica cuántica, se encuentra una ecuación de autovalores para los orbitales de una partícula.

El punto de partida para el cálculo Hartree-Fock es un conjunto de orbitales aproximados. Para un cálculo atómico, estos son típicamente los orbitales de un átomo hidrogenoide (un átomo con una carga nuclear cualquiera pero con un solo electrón). Para cálculos moleculares o cristalinos, las funciones de ondas iniciales son típicamente una combinación lineal de orbitales atómicos. Esto da una colección de orbitales monoelectrónicos, que por la naturaleza fermiónica de los electrones, debe ser antisimétrica, lo que se consigue mediante el uso del determinante de Slater. El procedimiento básico fue diseñado por Hartree, y Fock añadió el antisimetrizado.

Una vez se ha construido una función de ondas inicial, se elige un electrón. Se resume el efecto de todos los demás electrones, que se usa para generar un potencial. (Por este motivo, se llama a veces a este método un procedimiento de campo promedio). Esto da un electrón en un campo definido, para el que se puede resolver la ecuación de Schrödinger, dando una función de ondas ligeramente diferente para este electrón. Entonces, el procedimiento se repite para cada uno de los otros electrones, hasta completar un paso del procedimiento. De esta forma, con la nueva distribución electrónica se tiene un nuevo potencial eléctrico. El procedimiento se repite, hasta alcanzar la convergencia (hasta que el cambio entre un paso y el siguiente es lo suficientemente pequeño).

Las ecuaciones de Hartree-Fock en su forma canónica se asemejan al problema de valores propios:

Ecuaciones de Hartree-Fock Canónicas ${\displaystyle {\hat {f}}\psi _{i}=\varepsilon _{i}\psi _{i}}$

donde ${\displaystyle {\hat {f}}}$ es el operador de Fock, que actúa sobre un electrón en el espín-orbital ${\displaystyle \psi _{i}}$ y depende de la forma de los otros ${\displaystyle N-1}$ espín-orbitales, y ${\displaystyle \varepsilon _{i}}$ es la energía orbital del espín-orbital ${\displaystyle \psi _{i}}$. Las ecuaciones de Hartree-Fock pierden la linealidad de la ecuación de Schrödinger independiente del tiempo y son ecuaciones integrodiferenciales acopladas, por lo que para encontrar el estado basal de un sistema de ${\displaystyle N}$ electrones es necesario utilizar procedimientos iterativos.

El método de Hartree-Fock para un sistema de ${\displaystyle N}$ electrones considera como función de onda un determinante de Slater construido a partir de los espín-orbitales del sistema. Debido a esto, sólo se considera la correlación estadística de los electrones con el mismo espín, lo que se refleja en que la probabilidad de encontrar a dos electrones del mismo espín en el mismo punto del espacio sea nula, pero, en cambio, si los electrones tienen espines distintos, entonces la probabilidad de que estén en el mismo punto del espacio es distinta de cero.^[1]​

Con este tipo de función de onda, el valor esperado del Hamiltoniano está dado por

${\displaystyle \langle {\hat {H}}\rangle =\sum _{i=1}^{N}h_{i}+{\frac {1}{2}}\sum _{i=1}^{N}\sum _{j\neq i}^{N}J_{ij}-K_{ij}}$,

donde ${\displaystyle h_{i}}$ es la integral monoelectrónica, que contiene la energía cinética y la atracción nuclear del electrón ${\displaystyle i}$, ${\displaystyle J_{ij}}$ es la integral de Coulomb y ${\displaystyle K_{ij}}$ es la integral de intercambio.

Siguiendo el principio variacional, el valor esperado del Hamiltoniano se minimiza bajo la condición de ortonormalidad de los espín-orbitales. Construyendo un funcional del conjunto de espín-orbitales

${\displaystyle {\mathcal {L}}[\{\psi _{i}\}]=\langle {\hat {H}}\rangle -\sum _{i=1}^{N}\sum _{j=1}^{N}\varepsilon _{ij}\left(\langle \psi _{i}|\psi _{j}\rangle -\delta _{ij}\right)}$,

donde ${\displaystyle \varepsilon _{ij}}$ son los multiplicadores de Lagrange, se puede minimizar ${\displaystyle {\mathcal {L}}}$ respecto a variaciones en la forma de los espín-orbitales. Una manera de realizar el proceso es calculando la derivada funcional de ${\displaystyle {\mathcal {L}}}$ respecto a un espín-orbital en específico, siendo esto:

${\displaystyle {\frac {\delta {\mathcal {L}}[\{\psi _{i}\}]}{\delta \psi _{k}^{*}}}=\left[{\hat {h}}+\sum _{j=1}^{N}\left(\int {\frac {|\psi _{j}(\mathbf {x} _{j})|^{2}}{r_{kj}}}d\mathbf {x} _{j}-\int {\frac {\psi _{j}^{*}(\mathbf {x} _{j}){\hat {P}}_{kj}\psi _{j}(\mathbf {x} _{j})}{r_{kj}}}d\mathbf {x} _{j}\right)\right]\psi _{k}(\mathbf {x} _{k})-\sum _{j=1}^{N}\varepsilon _{kj}\psi _{j}(\mathbf {x} _{k})=0}$,

donde ${\displaystyle {\hat {h}}}$ es el operador monoelectrónico, que incluye el operador de energía cinética y el de atracción nuclear, y ${\displaystyle {\hat {P}}_{kj}}$ es el operador de permutación entre los electrones ${\displaystyle k}$ y ${\displaystyle j}$. La relación se satisface en los puntos estacionarios del funcional, de donde surgen las ecuaciones de Hartree-Fock:

${\displaystyle \left({\hat {h}}+\sum _{j=1}^{N}\left[\int {\frac {|\psi _{j}(\mathbf {x} _{j})|^{2}}{r_{kj}}}d\mathbf {x} _{j}-\int {\frac {\psi _{j}^{*}(\mathbf {x} _{j}){\hat {P}}_{kj}\psi _{j}(\mathbf {x} _{j})}{r_{kj}}}d\mathbf {x} _{j}\right]\right)\psi _{k}(\mathbf {x} _{k})=\sum _{j=1}^{N}\varepsilon _{kj}\psi _{j}(\mathbf {x} _{k})}$

Definiendo el operador de Coulomb para el electrón ${\displaystyle k}$ como

${\displaystyle {\hat {J}}_{k}[\{\psi _{j}(\mathbf {x} _{j})\}]=\int {\frac {|\psi _{j}(\mathbf {x} _{j})|^{2}}{r_{kj}}}d\mathbf {x} _{j}}$,

el operador de intercambio para el electrón ${\displaystyle k}$:

${\displaystyle {\hat {K}}_{k}[\{\psi _{j}(\mathbf {x} _{j})\}]=\int {\frac {\psi _{j}^{*}(\mathbf {x} _{j}){\hat {P}}_{kj}\psi _{j}(\mathbf {x} _{j})}{r_{kj}}}d\mathbf {x} _{j}}$,

y el operador de Fock para el electrón ${\displaystyle k}$ como

${\displaystyle {\hat {f}}_{k}[\{\psi _{j}(\mathbf {x} _{j})\}]={\hat {h}}_{k}+\sum _{j=1}^{N}{\hat {J}}_{k}[\{\psi _{j}(\mathbf {x} _{j})\}]-{\hat {K}}_{k}[\{\psi _{j}(\mathbf {x} _{j})\}]}$,

las ecuaciones de Hartree-Fock para el electrón ${\displaystyle k}$ se pueden escribir de forma compacta como

${\displaystyle {\hat {f}}_{k}[\{\psi _{j}(\mathbf {x} _{j})\}]\psi _{k}(\mathbf {x} _{k})=\sum _{j=1}^{N}\varepsilon _{kj}\psi _{j}(\mathbf {x} _{k})}$

Esta combinación lineal sugiere que las soluciones a las ecuaciones de Hartree-Fock no son únicas. En general, una transformación unitaria de los elementos de un determinante no cambian el valor del mismo. Esta propiedad se puede aplicar a los espín-orbitales y está permitida por la forma de un único determinante de Slater de la función de onda Hartree-Fock. De acuerdo a esto, se puede demostrar que el operador de Fock y los multiplicadores de Lagrange son invariantes ante transformaciones unitarias de los espín-orbitales. Entonces, identificando a los multiplicadores de Lagrange como los elementos de matriz del operador de Fock,

${\displaystyle \langle \psi _{m}|{\hat {f}}_{k}|\psi _{n}\rangle =\sum _{n}^{N}\varepsilon _{kn}\langle \psi _{m}|\psi _{n}\rangle =\varepsilon _{mn}}$,

se puede fijar de manera única a los espín-orbitales mediante una transformación unitaria que diagonaliza la matriz de multiplicadores de Lagrange (ya que es hermítica). El conjunto de espín-orbitales que cumple con este propósito se conoce como los espín-orbitales canónicos, que es un conjunto único. Con esto, se obtienen las ecuaciones de Hartree-Fock canónicas que tienen una forma similar al problema de valores propios:

${\displaystyle {\hat {f}}_{k}[\{\psi _{j}(\mathbf {x} _{j})\}]\psi _{k}(\mathbf {x} _{k})=\varepsilon _{k}\psi _{k}(\mathbf {x} _{k})}$

Hay que destacar que el operador de Coulomb tanto como el operador de intercambio (y por consecuencia, el operador de Fock), dependen de los espín-orbitales en que se encuentran los ${\displaystyle N-1}$ electrones. Por lo tanto, la linealidad de la ecuación de Schrödinger se pierde y las ecuaciones de Hartree-Fock son ecuaciones integrodiferenciales acopladas, ya que el operador de Fock es un operador monoelectrónico que depende de los espín-orbitales de los otros electrones, cada uno con su propio operador de Fock. En otras palabras, para resolver las ecuaciones de HF para el electrón 1, se debe conocer con anticipación los espín-orbitales de los otros electrones, para luego resolver las ecuaciones de HF para los demás electrones con los espín-orbitales mejorados, hasta llegar a un límite de convergencia de la energía Hartree-Fock.

La forma del operador de Fock trae consigo la cancelación de los efectos de autointeracción, es decir, las integrales de Coulomb e intercambio para la interacción de un electrón consigo mismo (${\displaystyle J_{ii}}$ y ${\displaystyle K_{ii}}$, respectivamente) se anulan de manera exacta.

${\displaystyle ({\hat {J}}_{i}-{\hat {K}}_{i})\psi _{i}=0}$

Además, ${\displaystyle {\hat {f}}}$ actúa sobre un electrón a la vez y también se puede definir como la suma entre el operador Hamiltoniano monoelectrónico, ${\displaystyle {\hat {h}}_{i}}$, y un operador de potencial efectivo monoelectrónico llamado potencial de Hartree-Fock:

${\displaystyle {\hat {v}}_{i}^{HF}[\psi _{j}]=\sum _{j}^{N}{\hat {J}}_{i}[\psi _{j}]-{\hat {K}}_{i}[\psi _{j}]}$

Al calcular el valor esperado del potencial de Hartree-Fock para el electrón ${\displaystyle i}$, hay que considerar que la forma del potencial de Coulomb representa una interacción promedio del electrón ${\displaystyle i}$ con los demás electrones. Como ${\displaystyle |\psi _{i}(\mathbf {x} _{i})|^{2}|\psi _{j}(\mathbf {x} _{j})|^{2}d\mathbf {x} _{i}d\mathbf {x} _{j}}$ es la probabilidad de que el electrón ${\displaystyle i}$ en el espín-orbital ${\displaystyle \psi _{i}}$ ocupe el elemento de volumen ${\displaystyle d\mathbf {x} _{i}}$ y el electrón ${\displaystyle j}$ en el espín-orbital ${\displaystyle \psi _{j}}$ ocupe el elemento de volumen ${\displaystyle d\mathbf {x} _{j}}$, entonces la integral sobre todo el espacio es el promedio de la interacción de Coulomb ${\displaystyle r_{ij}^{-1}}$.

${\displaystyle \langle ij|ij\rangle =\langle \psi _{i}\psi _{j}|r_{ij}^{-1}|\psi _{i}\psi _{j}\rangle =\int \left(|\psi _{i}(\mathbf {x} _{i})|^{2}|\psi _{j}(\mathbf {x} _{j})|^{2}\right){\frac {1}{r_{ij}}}d\mathbf {x} _{i}d\mathbf {x} _{j}}$

Por lo tanto, se dice que el operador de Coulomb representa el potencial de Coulomb promedio local en ${\displaystyle \mathbf {x} _{i}}$, que resulta del electrón en ${\displaystyle \psi _{j}}$. Caso contrario es lo que sucede con el operador de intercambio. Como no existe un potencial simple ${\displaystyle K_{i}(\mathbf {x} _{j})}$ únicamente definido en un punto en el espacio, se dice que el operador de intercambio es un operador no-local. En otras palabras, el resultado de aplicar ${\displaystyle {\hat {K}}_{i}}$ sobre ${\displaystyle \psi _{i}(\mathbf {x} _{i})}$ depende del valor de ${\displaystyle \psi _{i}}$ en todo el espacio y no únicamente en el punto ${\displaystyle \mathbf {x} _{i}}$.

Como se ha mencionado en las secciones anteriores, el operador de Fock depende en forma funcional de los espín-orbitales ocupados. Sin embargo, luego de terminar el proceso autoconsistente, es decir, cuando la forma final de los espín-orbitales ya es conocida (dentro de un cierto margen), el operador de Fock se vuelve un operador hermítico bien definido con un número infinito de funciones propias.^[1]​

${\displaystyle {\hat {f}}|\psi _{j}\rangle =\varepsilon _{j}|\psi _{j}\rangle ,\quad j=1,2,\dots ,\infty }$

Cada una de las soluciones ${\displaystyle |\psi _{j}\rangle }$ tiene un valor propio ${\displaystyle \varepsilon _{j}}$ asociado. La interpretación de los multiplicadores de Lagrange es que representan las energías orbitales de los espín-orbitales, que ahora se pueden separar en dos tipos: los espín-orbitales ocupados y los espín-orbitales virtuales. Para un sistema de ${\displaystyle N}$ electrones, los ${\displaystyle N}$ espín-orbitales con energías más bajas pasan a llamarse espín-orbitales ocupados, mientras que el resto de los espín-orbitales que pertenece al conjunto solución se les conoce como espín-orbitales virtuales. Estos últimos son importantes cuando se requiere hacer un cálculo post-HF, cuando el determinante de referencia es la solución a las ecuaciones de Hartree-Fock.

La representación matricial del operador de Fock en la base de los espín-orbitales canónicos es diagonal, siendo sus elementos las energías orbitales:

${\displaystyle \langle \psi _{i}|{\hat {f}}|\psi _{j}\rangle =\varepsilon _{j}\langle \psi _{i}|\psi _{j}\rangle =\varepsilon _{j}\delta _{ij}}$

Por lo tanto, las energías de los espín-orbitales ocupados (etiquetados como ${\displaystyle a,b,\dots }$) están dadas por

${\displaystyle \varepsilon _{a}=h_{a}+\sum _{b\neq a}^{N}J_{ab}-K_{ab}}$,

mientras que para los espín-orbitales virtuales (etiquetados como ${\displaystyle r,s,\dots }$) se tiene

${\displaystyle \varepsilon _{r}=h_{r}+\sum _{b=1}^{N}J_{rb}-K_{rb}}$

Es importante destacar que los espín-orbitales ocupados tienen una interacción de Coulomb e intercambio con los ${\displaystyle N-1}$ orbitales restantes (lo que explica que el índice ${\displaystyle b}$ sea distinto de ${\displaystyle a}$), mientras que la energía de los espín-orbitales virtuales incluye la interacción con los ${\displaystyle N}$ electrones del estado basal Hartree-Fock. Además, la energía total en Hartree-Fock se puede escribir en términos de las energías orbitales como

${\displaystyle E_{\text{HF}}=\sum _{a=1}^{N}\varepsilon _{a}-{\frac {1}{2}}\sum _{a=1}^{N}\sum _{b\neq a}^{N}J_{ab}-K_{ab}}$

Que la energía total no sea una suma simple de las energías orbitales se ha considerado como una muestra de que los electrones en la teoría HF están estadísticamente correlacionados. Si bien, los electrones se tratan como partículas independientes moviéndose en un potencial efectivo, el potencial de HF se construye a partir de la interacción promediada con los otros electrones del sistema, por lo que las partículas no se consideran completamente independientes. Sin embargo, se ha encontrado que la energía total de una molécula con todos los electrones apareados por espín, en la geometría de equilibrio, está aproximadamente dada por^[2]​^[3]​

${\displaystyle E_{\text{HF}}\approx 1.55\sum _{i=1}^{N/2}2\varepsilon _{i}}$

Existen distintos tipos de teorías Hartree-Fock que son clasificados de acuerdo al tipo de restricción a la simetría de la función de onda. La forma más general, conocida como método Hartree-Fock generalizado, consiste en aplicar el principio variacional sin restricciones como que la función de onda final sea función propia de los operadores de espín ${\displaystyle {\hat {S}}_{z}}$ y ${\displaystyle {\hat {S}}^{2}}$, de reversión temporal, ${\displaystyle {\hat {\Theta }}}$, y en general de operadores del grupo puntual del sistema.^[4]​ Tiene dos variantes principales que están relacionadas a considerar si la función de onda es compleja o real.

En la teoría GHF, los espín orbitales tienen la forma

${\displaystyle \psi _{i}(\mathbf {r} ,\sigma )=\phi _{i}^{\alpha }(\mathbf {r} )\alpha (\sigma )+\phi _{i}^{\beta }(\mathbf {r} )\beta (\sigma )}$

donde ${\displaystyle {\ce {\phi_{i}({\mathbf {r}})}}}$ es la parte espacial del espín-orbital ${\displaystyle \psi _{i}}$ y ${\displaystyle \alpha (\sigma )}$ y ${\displaystyle \beta (\sigma )}$ son funciones de las coordenadas de espín del electrón, ${\displaystyle \sigma }$. Es conveniente analizar las propiedades estadísticas de la función de onda GHF mediante la probabilidad de encontrar a un electrón en un cierto punto del espacio, con los otros electrones en cualquier parte, y la probabilidad de encontrar simultáneamente a un electrón en un punto dado y a otro electrón en otro punto del espacio, que puede ser el mismo que el anterior o no, con los otros electrones en cualquier otra parte. La densidad de probabilidad de encontrar a un electrón con las coordenadas de espacio y espín ${\displaystyle \mathbf {x} }$ para una función de onda GHF está dada por

${\displaystyle \rho _{1}^{GHF}(\mathbf {x} )=\sum _{i}^{N}|\psi _{i}(\mathbf {x} )|^{2}}$

Esta función también es conocida como la densidad electrónica cuando se integra sobre las coordenadas de espín, siendo también la diagonal de la matriz de densidad reducida sin espín de primer orden. Por otra parte, la densidad de probabilidad de encontrar simultáneamente a un electrón en las coordenadas ${\displaystyle \mathbf {x} }$, a otro electrón en ${\displaystyle \mathbf {x} '}$ y al resto en cualquier parte está dada por

${\displaystyle \rho _{2}^{GHF}(\mathbf {x} ,\mathbf {x} ')=\rho _{1}^{GHF}(\mathbf {x} )\rho _{1}^{GHF}(\mathbf {x} ')-\sum _{i}^{N}\sum _{j}^{N}\psi _{i}(\mathbf {x} )\psi _{i}^{*}(\mathbf {x} ')\psi _{j}(\mathbf {x} ')\psi _{j}^{*}(\mathbf {x} )}$

Si los electrones fueran completamente independientes, la densidad de pares, ${\displaystyle \rho _{2}(\mathbf {x} ,\mathbf {x} ')}$, sólo sería la multiplicación de las funciones de distribución de cada partícula, o sea, solamente el primer término. Sin embargo, para la función de onda GHF se obtiene un término extra que surge debido a la correlación de Fermi (intercambio) de los electrones, lo que se asocia a la forma antisimétrica de la función de onda.

Las ecuaciones de Hartree-Fock generalizadas usando los espín-orbitales canónicos está dada por

${\displaystyle {\hat {f}}_{i}^{GHF}[\psi _{j}]\psi _{i}(\mathbf {x} _{i})=\varepsilon _{i}\psi _{i}(\mathbf {x} _{i})}$,

que es una forma similar a la encontrada en la sección anterior, con la diferencia de que la forma de los espín-orbitales está identificada.

El uso de las funciones de onda GHF es limitado, ya sea en su variante de soluciones complejas o reales. Aun así, se han encontrado soluciones GHF para el átomo de berilio, en el BH, en el H₄ y en otros pocos casos.^[4]​

A pesar del nombre, este método impone algunas restricciones a la función de onda comparado con el método generalizado. La forma que tienen los espín-orbitales ya no es una mezcla entre funciones de espín ${\displaystyle \alpha }$ y ${\displaystyle \beta }$, y también posee dos variantes en que las soluciones pueden ser reales o complejas. Un espín-orbital en la teoría UHF se escribe como^[1]​^[2]​^[5]​^[6]​

${\displaystyle \psi _{i}(\mathbf {r} ,\sigma )=\phi _{i}^{s_{i}}(\mathbf {r} )s_{i}(\sigma )}$

Es decir, ahora se compone por una parte espacial ${\displaystyle \phi _{i}}$ asociada a una parte de espín, donde ${\displaystyle s_{i}}$ puede ser ${\displaystyle \alpha }$ o ${\displaystyle \beta }$. Además, la función de onda total es también función propia del operador de espín ${\displaystyle {\hat {S}}_{z}}$. La densidad electrónica para la función de onda UHF estará dada por

${\displaystyle \rho _{1}^{UHF}(\mathbf {r} )=\sum _{a=1}^{N}|\phi _{a}(\mathbf {r} )|^{2}}$

La densidad de probabilidad de encontrar simultáneamente a un electrón en ${\displaystyle \mathbf {r} }$ y a otro en ${\displaystyle \mathbf {r} '}$, con los otros en cualquier parte, será

${\displaystyle \rho _{2}^{UHF}(\mathbf {r} ,\mathbf {r} ')=\rho _{1}^{UHF}(\mathbf {r} )\rho _{1}^{UHF}(\mathbf {r} ')-\sum _{a=1}^{N}\sum _{b=1}^{N}\delta _{\sigma \sigma '}\phi _{a}(\mathbf {r} )\phi _{a}^{*}(\mathbf {r} ')\phi _{b}(\mathbf {r} ')\phi _{b}^{*}(\mathbf {r} )}$

donde ${\displaystyle \delta _{\sigma \sigma '}}$ es la delta de Kronecker para indicar que si los espines de los electrones considerados son iguales, entonces el segundo término que correlaciona el movimiento de éstos existirá (${\displaystyle \delta _{\sigma \sigma }=1}$). En caso contrario, cuando los dos electrones tienen espines opuestos, la delta es igual a cero y la densidad de probabilidad de pares se vuelve

${\displaystyle \rho _{2}^{UHF}(\mathbf {r} ,\mathbf {r} ')=\rho _{1}^{UHF}(\mathbf {r} )\rho _{1}^{UHF}(\mathbf {r} ')}$

Esto sugiere que los electrones de espín opuesto no están estadísticamente correlacionados, mientras que los electrones del mismo espín si "pueden verse" entre ellos. Además, en el caso en que los dos electrones considerados tengan el mismo espín y la misma posición en el espacio (${\displaystyle \mathbf {r} '=\mathbf {r} }$), la densidad de probabilidad de pares es 0.

Comparando al caso general, en la teoría UHF existen menos multiplicadores de Lagrange debido a la restricción a la simetría respecto al espín sobre un eje. Los espín-orbitales ahora se dividen en dos conjuntos de acuerdo a su espín, por lo que existirán operadores de Fock para los espín-orbitales ${\displaystyle \alpha }$ y operadores de Fock para los espín-orbitales ${\displaystyle \beta }$. Esto se hace evidente al aplicar el operador de Fock sobre un espín-orbital ${\displaystyle \psi (\mathbf {x} )=\phi _{i}^{\alpha }(\mathbf {r} )\alpha (\sigma )}$, donde el superíndice hace explícita la asociación entre el orbital ${\displaystyle \phi _{i}}$ a su forma cuando el electrón que lo ocupa tiene espín ${\displaystyle \alpha }$.

${\displaystyle {\hat {f}}\psi (\mathbf {x} )={\hat {f}}\phi _{i}^{\alpha }(\mathbf {r} )\alpha (\sigma )=\varepsilon _{i}^{\alpha }\phi _{i}^{\alpha }(\mathbf {r} )\alpha (\sigma )}$

Multiplicando por la izquierda por ${\displaystyle \alpha ^{*}(\sigma )}$ e integrando esta ecuación respecto a las coordenadas de espín del electrón, se obtiene

${\displaystyle \langle \alpha |{\hat {f}}|\alpha \rangle \phi _{i}^{\alpha }(\mathbf {r} )={\hat {f}}^{\alpha }[\phi _{j}]\phi _{i}^{\alpha }(\mathbf {r} )=\varepsilon _{i}^{\alpha }\phi _{i}^{\alpha }(\mathbf {r} )}$

donde se define el operador de Fock espacial, ${\displaystyle {\hat {f}}^{\alpha }[\phi _{j}]}$ para los electrones de espín ${\displaystyle \alpha }$, que depende funcionalmente de los orbitales espaciales de los otros electrones. De manera análoga, se define el operador de Fock espacial para los electrones de espín ${\displaystyle \beta }$. El operador de Fock para espín ${\displaystyle \alpha }$ contiene al Hamiltoniano monoelectrónico y al potencial efectivo que actúa sobre el electrón de espín ${\displaystyle \alpha }$, que está dado por la interacción de Coulomb e intercambio con todos los otros electrones de espín ${\displaystyle \alpha }$ y la interacción de Coulomb con los electrones de espín ${\displaystyle \beta }$, ya que la correlación de los electrones de espín opuesto no es considerada. Entonces, la forma de los operadores de Fock para ambos espín es

${\displaystyle {\hat {f}}^{\alpha }[\phi _{a}]={\hat {h}}+\sum _{a}^{N^{\alpha }}{\hat {J}}^{\alpha }[\phi _{a}^{\alpha }]-{\hat {K}}^{\alpha }[\phi _{a}^{\alpha }]+\sum _{a}^{N^{\beta }}{\hat {J}}^{\alpha }[\phi _{a}^{\beta }]}$

${\displaystyle {\hat {f}}^{\beta }[\phi _{a}]={\hat {h}}+\sum _{a}^{N^{\beta }}{\hat {J}}^{\beta }[\phi _{a}^{\beta }]-{\hat {K}}^{\beta }[\phi _{a}^{\beta }]+\sum _{a}^{N^{\alpha }}{\hat {J}}^{\alpha }[\phi _{a}^{\alpha }]}$

donde las definiciones para los operadores de Coulomb e interacambio para cada espín son análogas a las formas generales, pero usando los orbitales espaciales. Las energías orbitales del operador de Fock para cada espín son

${\displaystyle \varepsilon _{a}^{\alpha }=h_{i}^{\alpha }+\sum _{b}^{N^{\alpha }}J_{ab}^{\alpha \alpha }-K_{ab}^{\alpha \alpha }+\sum _{b}^{N^{\beta }}J_{ab}^{\alpha \beta }}$

${\displaystyle \varepsilon _{a}^{\beta }=h_{i}^{\beta }+\sum _{b}^{N^{\beta }}J_{ab}^{\beta \beta }-K_{ab}^{\beta \beta }+\sum _{b}^{N^{\alpha }}J_{ab}^{\beta \alpha }}$

Así, la energía total en el método UHF, en términos de las energías orbitales, se lee como

${\displaystyle E_{UHF}=\sum _{a}^{N^{\alpha }}\varepsilon _{a}^{\alpha }+\sum _{a}^{N^{\beta }}\varepsilon _{a}^{\beta }-{\frac {1}{2}}\sum _{a}^{N^{\alpha }}\sum _{b}^{N^{\alpha }}J_{ab}^{\alpha \alpha }-K_{ab}^{\alpha \alpha }-{\frac {1}{2}}\sum _{a}^{N^{\beta }}\sum _{b}^{N^{\beta }}J_{ab}^{\beta \beta }-K_{ab}^{\beta \beta }-\sum _{a}^{N^{\alpha }}\sum _{b}^{N^{\beta }}J_{ab}^{\alpha \beta }}$

Para los casos en que el sistema a considerar tiene un número par de electrones, con todos ellos apareados por espín, existe una formulación más sencilla del método Hartree-Fock conocida como Hartree-Fock Restringido. En este método, se requiere que los espín-orbitales sean un producto de funciones espaciales y funciones de espín que sean funciones propias del operador ${\displaystyle {\hat {S}}_{z}}$ y el número de electrones ${\displaystyle \alpha }$ debe ser igual al número de electrones ${\displaystyle \beta }$. Además, los orbitales espaciales se comparten por dos electrones de distinto espín, es decir, un orbital ${\displaystyle \phi _{i}(\mathbf {r} )}$ está multiplicado por ${\displaystyle \alpha (\sigma )}$ para el electrón de espín ${\displaystyle \alpha }$ y el mismo orbital ${\displaystyle \phi _{i}(\mathbf {r} )}$ está multiplicado por ${\displaystyle \beta (\sigma ')}$ para el electrón de espín ${\displaystyle \beta }$ que lo ocupa. Esta situación también se conoce como un sistema de capa cerrada, donde los espín-orbitales para el conjunto ${\displaystyle \alpha }$ y para el conjunto ${\displaystyle \beta }$ están dados por^[1]​^[2]​

${\displaystyle \psi _{i}(\mathbf {x} )=\phi _{i}(\mathbf {r} )s_{i}(\sigma )}$

donde ${\displaystyle s_{i}}$ puede ser ${\displaystyle \alpha }$ o ${\displaystyle \beta }$ y la parte espacial será la misma en cualquiera de los dos casos. Debido a esta restricción, la función de onda RHF es también una función propia del operador de espín ${\displaystyle {\hat {S}}^{2}}$. Al ser útil solo en sistemas de capa cerrada, la función de onda RHF está asociada al valor propio 0 para los dos operadores de espín. Al igual que en los dos casos anteriores, también posee dos variantes que dependen de si los espín-orbitales son complejos o reales.

En este caso, la densidad electrónica está dada por

${\displaystyle \rho _{1}^{RHF}(\mathbf {r} )=\sum _{a=1}^{N^{\alpha }}|\phi _{a}(\mathbf {r} )|^{2}+\sum _{a=1}^{N^{\beta }}|\phi _{a}(\mathbf {r} )|^{2}=2\sum _{a=1}^{N/2}|\phi _{a}(\mathbf {r} )|^{2}}$

donde ${\displaystyle N^{\alpha }=N^{\beta }}$ y, como la parte espacial es igual para ambos espines, se cuenta hasta ${\displaystyle N/2}$ en ambos casos. A su vez, la densidad de probabilidad de pares está dada por

${\displaystyle \rho _{2}^{RHF}(\mathbf {r} ,\mathbf {r} ')=\rho _{1}^{RHF}(\mathbf {r} )\rho _{1}^{RHF}(\mathbf {r} ')-2\sum _{i}^{N/2}\sum _{j}^{N/2}\delta _{\sigma \sigma '}\phi _{i}(\mathbf {r} )\phi _{i}^{*}(\mathbf {r} ')\phi _{j}(\mathbf {r} ')\phi _{j}^{*}(\mathbf {r} )}$

donde, al igual que en UHF, existe un término extra debido a la correlación estadística entre los electrones del mismo espín. Cuando los electrones tienen distinto espín, ${\displaystyle \sigma \neq \sigma '}$, el último término es cero y los electrones son independientes estadísticamente.

El operador de Fock pasa a ser el operador de Fock espacial como en UHF, sin embargo, ahora los conjuntos ${\displaystyle \alpha }$ y ${\displaystyle \beta }$ son iguales, por lo que sólo existirá un operador de Fock espacial que actúe sobre electrones ${\displaystyle \alpha }$ y ${\displaystyle \beta }$. Entonces, la forma que tiene el operador de Fock es

${\displaystyle {\hat {f}}[\phi _{b}]={\hat {h}}+\sum _{b}^{N/2}2{\hat {J}}[\phi _{b}]-{\hat {K}}[\phi _{b}]}$

donde los operadores de Coulomb e intercambio están definidos usando las funciones espaciales. Así, las ecuaciones del método Hartree-Fock Restringido se leen como

${\displaystyle {\hat {f}}[\phi _{b}]\phi _{a}(\mathbf {r} )=\varepsilon _{a}\phi _{a}(\mathbf {r} )}$

Las energías orbitales toman la forma

${\displaystyle \varepsilon _{a}=h_{a}+\sum _{b}^{N/2}2J_{ab}-K_{ab}}$

y la energía total está dada por

${\displaystyle {\begin{array}{lcl}E_{RHF}&=&2\sum _{a}^{N/2}h_{a}+\sum _{a}^{N/2}\sum _{b}^{N/2}2J_{ab}-K_{ab}=2\sum _{a}^{N/2}\varepsilon _{a}-\sum _{a}^{N/2}\sum _{b}^{N/2}2J_{ab}-K_{ab}\\&=&\sum _{a}^{N/2}\varepsilon _{a}+h_{a}\end{array}}}$

Una forma de resolver las ecuaciones de Hartree-Fock es mediante la introducción de un conjunto base finito, lo que se conoce como aproximación algebraica. El fundamento de este método está en la combinación lineal de orbitales atómicos introducidos por Lennard-Jones en 1929 para la descripción de moléculas diatomicas. En este sentido, para comenzar a resolver las ecuaciones de Hartree-Fock es necesario tener un conocimiento previo del sistema y así poder tener una forma inicial del operador de Fock, entonces, como los espín-orbitales no tienen una forma conocida, se expanden en términos de funciones que sí son conocidas (como los orbitales de un átomo hidrogenoide), las que pueden incluir parámetros variacionales, además de los coeficientes de expansión, que permitan encontrar el punto estacionario del funcional de energía restringido.

Para un sistema de capa cerrada, donde es necesario resolver las ecuaciones de HF restringidas sin considerar el espín, los orbitales se pueden expandir en un conjunto finito de ${\displaystyle K}$ funciones base, ${\displaystyle \{\chi _{\mu }\}}$:^[1]​

${\displaystyle \phi _{i}(\mathbf {r} )=\sum _{\mu =1}^{K}C_{\mu i}\chi _{\mu }(\mathbf {r} ),\qquad i=1,2,\dots ,{\frac {N}{2}};\quad K\geq {\frac {N}{2}}}$

donde ${\displaystyle C_{\mu i}}$ son los coeficientes de expansión. Si el conjunto base es completo, entonces esta expansión es exacta. Sin embargo, debido a las restricciones computacionales, siempre se restringe este conjunto a un cierto número de funciones base. Introduciendo esta combinación lineal en las ecuaciones de HF, se obtiene

${\displaystyle {\hat {f}}_{i}[\phi _{b}]\sum _{\nu }^{K}C_{\nu a}\chi _{\nu }(\mathbf {r} _{i})=\varepsilon _{a}\sum _{\nu }^{K}C_{\nu a}\chi _{\nu }(\mathbf {r} _{i})}$

Multiplicando desde la izquierda por ${\displaystyle \chi _{\mu }^{*}(\mathbf {r} _{i})}$ e integrando, se llega a

${\displaystyle \sum _{\nu }C_{\nu a}\int \chi _{\mu }^{*}(\mathbf {r} _{i}){\hat {f}}_{i}\chi _{\nu }(\mathbf {r} _{i})d\mathbf {r} _{i}=\varepsilon _{a}\sum _{\nu }C_{\nu a}\int \chi _{\mu }^{*}(\mathbf {r} _{i})\chi _{\nu }(\mathbf {r} _{i})d\mathbf {r} _{i}}$

Definiendo la matriz de solapamiento

${\displaystyle S_{\mu \nu }=\int \chi _{\mu }^{*}(\mathbf {r} _{i})\chi _{\nu }(\mathbf {r} _{i})d\mathbf {r} _{i}}$

y la matriz de Fock

${\displaystyle F_{\mu \nu }=\int \chi _{\mu }^{*}(\mathbf {r} _{i}){\hat {f}}_{i}\chi _{\nu }(\mathbf {r} _{i})d\mathbf {r} _{i}}$

las ecuaciones de Hartree-Fock toman la forma matricial

${\displaystyle \sum _{\nu }^{K}F_{\mu \nu }C_{\mu a}=\varepsilon _{a}\sum _{\nu }^{K}S_{\mu \nu }C_{\nu a},\qquad {\text{o bien,}}\qquad \mathbf {FC} =\mathbf {SC} {\boldsymbol {\varepsilon }}}$

donde ${\displaystyle \mathbf {C} }$ es la matriz de los coeficientes de expansión y ${\displaystyle {\boldsymbol {\varepsilon }}}$ es la matriz de las energías orbitales. Esta forma es conocida como las ecuaciones de Roothaan-Hall.^[7]​^[8]​ Hay que destacar que se asume que la base es linealmente independiente y está normalizada, sin embargo, no necesariamente es ortogonal. Por lo tanto, la magnitud de los elementos de la matriz de solapamiento están en el intervalo ${\displaystyle [0,1]}$, siendo los elementos de la diagonal iguales a 1.

La expansión de los orbitales en términos de un conjunto base es conocida como aproximación algebraica debido a que el problema se transforma de ${\displaystyle N}$ ecuaciones integrodiferenciales acopladas (${\displaystyle N/2}$ para capa cerrada) a ${\displaystyle N}$ ecuaciones algebraicas (${\displaystyle N/2}$ para capa cerrada). Las ecuaciones de Roothaan-Hall sólo son útiles en el caso en que se resuelven las ecuaciones de HF restringidas, por lo que es necesario una generalización al caso de sistemas de capa abierta.

Siguiendo un procedimiento similar a la formulación de Roothaan-Hall, los orbitales se expanden en términos de una base ${\displaystyle \{\chi _{\mu }|\mu =1,2,\dots ,K\}}$ en que los coeficientes varían dependiendo si el espín-orbital es ${\displaystyle \alpha }$ o ${\displaystyle \beta }$.

${\displaystyle \phi _{i}^{\alpha }=\sum _{\mu }^{K}C_{\mu i}^{\alpha }\chi _{\mu }}$

${\displaystyle \phi _{i}^{\beta }=\sum _{\mu }^{K}C_{\mu i}^{\beta }\chi _{\mu }}$

En este caso, las dos ecuaciones de Hartree-Fock sin restringir para cada espín garantizan que los conjuntos ${\displaystyle \{\phi _{i}^{\alpha }\}}$ y ${\displaystyle \{\phi _{i}^{\beta }\}}$ sean conjuntos ortonormales por separado, sin embargo, un miembro del conjunto ${\displaystyle \alpha }$ no tiene por qué ser ortogonal a un miembro del conjunto ${\displaystyle \beta }$. Sustituyendo esta expansión en las ecuaciones UHF y siguiendo el mismo procedimiento descrito antes, se obtienen dos ecuaciones matriciales, conocidas como las ecuaciones de Pople-Nesbet:^[9]​

${\displaystyle \mathbf {F} ^{\alpha }\mathbf {C} ^{\alpha }=\mathbf {S} \mathbf {C} ^{\alpha }{\boldsymbol {\varepsilon }}^{\alpha }}$

${\displaystyle \mathbf {F} ^{\beta }\mathbf {C} ^{\beta }=\mathbf {S} \mathbf {C} ^{\beta }{\boldsymbol {\varepsilon }}^{\beta }}$

donde las dos matrices de energías orbitales son diagonales.

1. Especificar el sistema:
   1. Conjunto de coordenadas nucleares, asociadas a los correspondientes números atómicos
   2. Número total de electrones
   3. Funciones de base. La elección de una base puede ser crítico para llegar a una convergencia adecuada, y con sentido físico, y no hay un procedimiento general con garantía de éxito. Los científicos cuánticos hablan del arte de escoger bien la base de funciones.
2. Calcular todas las integrales (interacciones) relevantes para las funciones de base: las energías cinéticas medias, la atracción electrón-núcleo, las repulsiones bielectrónicas. Como las funciones de base se mantienen a lo largo de todo el cálculo, no es necesario volver a evaluar las integrales. Dependiendo de las limitaciones técnicas del momento y de la talla del sistema, las integrales pueden o no mantenerse en la RAM. En caso de que no se mantengan, la estrategia óptima puede ser guardarlas en un disco duro o cinta, o bien recalcularlas en cada momento en que son necesarias.
3. Construir, con las integrales calculadas, la matriz de traslape S, que mide la desviación de la ortogonalidad de las funciones de la base, y, a partir de ella, la matriz de transformación X, que ortogonaliza la base.
4. Obtener una estimación de la matriz densidad P que, a partir de un conjunto de funciones de base, especifica completamente la distribución de densidad electrónica. Nuevamente, la primera estimación no es obvia, y puede precisar de inspiración artística. Un cálculo de Hückel extendido puede suponer una buena aproximación.
5. Conociendo la matriz densidad y las integrales bielectrónicas de las funciones de base, calcular el operador de interacción entre electrones, la matriz G.
6. Construir la matriz de Fock como suma del hamiltoniano "fijo" (integrales monoelectrónicas) y la matriz G
7. Transformar, con la matriz de transformación, la matriz de Fock en su expresión para la base ortonormal, F'
8. Diagonalizar F', obtener C' y e (vectores y valores propios)
9. De C' y la matriz de transformación, recuperar C, que será la expresión en las funciones de base originales
10. C define una nueva matriz densidad P
11. Si la nueva matriz densidad difiere de la anterior más que un criterio previamente fijado (no ha convergido), volver al punto 5.
12. En caso contrario, usar C, P y F para calcular los valores esperados de magnitudes observables, y otras cantidades de interés.

Si el cálculo diverge, o converge con lentitud, o llega a una solución que no es una descripción adecuada de los fenómenos que son de interés,

• o bien se corrigen los dos puntos artesanales, por ejemplo, dando más flexibilidad a las funciones de base, o, por el contrario, restringiéndolas a la parte fundamental de la física, u obtener una mejor primera estimación de la matriz densidad P
• o bien se aplican métodos que van más allá de la aproximación de Hartree-Fock

Se usa a menudo en la misma área de cálculos que la Teoría del Funcional de la Densidad, que puede dar soluciones aproximadas para las energías de intercambio y de correlación. De hecho, es común el uso de cálculos que son híbridos de los dos métodos. Adicionalmente, los cálculos a nivel Hartree-Fock se usan como punto de partida para métodos más sofisticados, como la teoría perturbacional de muchos cuerpos, o cálculos cuánticos de Monte-Carlo.

La inestabilidad numérica es un problema de este método, y hay varias vías para combatirla. Una de las más básicas y más aplicadas es la mezcla-F. Con la mezcla-F, no se usa directamente la función de ondas de un electrón conforme se ha obtenido. En lugar de esto, se usa una combinación lineal de la función obtenida con las previas, por ejemplo con la inmediatamente previa. Otro truco, empleado por Hartree, es aumentar la carga nuclear para comprimir a los electrones; tras la estabilización del sistema, se reduce gradualmente la carga hasta llegar a la carga correcta.

Desarrollos más allá del campo autoconsistente o SCF son el CASSCF y la interacción de configuraciones. Los cálculos de este tipo son relativamente económicos frente a otros de la química cuántica. De esta forma, en ordenadores personales es posible resolver moléculas pequeñas en muy poco tiempo. Las moléculas más grandes, o los desarrollos más sofisticados, para obtener resultados más exactos, siguen realizándose en superordenadores. Existen múltiples paquetes informáticos que implementan el método de campo autoconsistente, entre los que pueden destacarse Gaussian, NWChem, MOLPRO, MOLCAS y LOWDIN.

• Atkins, P.W.; Friedman, R. (2005). Molecular Quantum Mechanics (4th edición). Oxford University Press. ISBN 978-0-19-927498-7. 
• Cedillo, A. (2023). Curso de Química Cuántica. Universidad Autónoma Metropolitana. ISBN 978-607-28-2870-4.