Tema2 Teoría Del Funcional
Cargado por
scribed123scribedTema2 Teoría Del Funcional
Cargado por
scribed123scribedUNED
Máster FSC. Funcional de la densidad
UNED
UNED Tema II
UNED
EL FORMALISMO DEL
FUNCIONAL DE LA DENSIDAD
PARA EL ESTADO FUNDAMENTAL
UNED
UNED
UNED 103
Máster FSC. Funcional de la densidad
UNED
UNED
UNED
UNED
UNED
UNED
UNED
UNED
Máster FSC. Funcional de la densidad
Capítulo 4
EL GAS HOMOGÉNEO DE ELECTRONES
UNED
OBJETIVOS:
1. Estudio de las propiedades del estado fundamental del gas homogéneo de electrones en la
UNED
aproximación Hartree-Fock.
2. Obtención de la energía de interacción bajo la aproximación RPA.
3. Cálculo de propiedades de excitación en la aproximación RPA.
LECTURAS COMPLEMENTARIAS:
UNED
1. Elementary Excitations in Solids, D. Pines, capítulo 3
ORIENTACIONES AL PLAN DE TRABAJO:
Este capítulo aborda el estudio de algunas propiedades básicas de un sistema modelo de muchos
electrones: el gas homogéneo de electrones.
10 horas a este capítulo desglosadas de la siguiente
UNED
A título orientativo, se espera que dedique
manera:
• 3 horas de lectura detallada y estudio del material didáctico correspondiente (el presente
capítulo 4)
• 7 horas para el a anzamiento de los contenidos, incluyendo consultas y las lecturas
complementarias [optativas].
UNED
Ni en este capítulo ni en los siguientes se proponen problemas.
GLOSARIO:
Gas de electrones homogéneo
Modelo de jellium
UNED
Radio de Wigner
Energía de Hartree
Cristal de Wigner
Fase paramagnética del gas de electrones
Fase ferromagnética del gas de electrones
Gas de electrones libres
Momento de Fermi
Energía de Fermi
105
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
Teorema de uctuación-disipación
Conexión adiabática
UNED
Teoría de perturbaciones de muchos cuerpos
Función de Lindhard
Apantallamiento
Excitación electrón-hueco
UNED
Plasmón
UNED
UNED
UNED
UNED
106
UNED
El gas homogéneo de electrones
Máster FSC. Funcional de la densidad
4.1. INTRODUCCIÓN
1. En este tema vamos a aplicar la teoría que hemos visto hasta la fecha a un sistema muy sencillo:
UNED
el gas de electrones homogéneo. Este modelo, debido originalmente a Sommerfeld e inspirado en el
modelo atómico de Drude, es una descripción burda de un sólido metálico cristalino en la cual la red
iónica es sustituida por un fondo de carga positiva uniforme de densidad n0 , donde n0 es la carga iónica
por unidad de volumen. Así, la densidad electrónica de electrones de valencia, n (~r), se superpone al
mencionado fondo de carga positiva. Aun a pesar de su simplicidad, en este sistema (llamado jellium
en la literatura) se esconden bastantes aspectos esenciales de la Teoría de Muchos Cuerpos y todavía se
sigue utilizando en investigación activa como laboratorio a la hora de ilustrar nuevas aproximaciones
UNED
relacionadas con las propiedades electrónicas de sistemas extensos.
2. Suponiendo que el sistema es eléctricamente neutro, el número de electrones de valencia es
N = V n0 , donde V es el volumen del sólido (imponemos, como es ya habitual, condiciones de contorno
periódicas de Born - von Karman). Así, el sistema queda únicamente caracterizado por el valor de n0
o, equivalentemente, por el llamado radio de Wigner rs dado por
1/3
3
rs =
4πn0
UNED
que es preferible a n0 ya que para sólidos típicos, el valor de rs asociado a la densidad de los electrones
de valencia está en el rango 2 . rs . 6 (en unidades atómicas de Hartree).
El hamiltoniano del sistema es
N 2
X p
b i 1X 1
H
b = + vext (b
ri ) + + Ered =
2 2 |b
ri − b
rj |
i=1 i6=j
UNED
(4.1)
N
b 2i b (~r2 ) − δ (~r1 − ~r2 ) n
Z Z
X p 1 n
b (~r1 ) n b (~r1 )
= + d3~r vext (~r) n
b (~r) + d3~r1 d3~r2 + Ered
2 2 |~r1 − ~r2 |
i=1
donde Z
n0
vext (~r) = − d3~r 0
|~r − ~r 0 |
1
es el potencial creado por el fondo de carga positiva y donde hemos incluido también el término de
UNED
interacción clásica del fondo de carga consigo mismo:
n20
Z Z
1 3 3 1
Ered = d ~r1 d ~r2 =− d3~r vext (~r) n0
2 |~r1 − ~r2 | 2
En el límite termodinámico, V → ∞, las energías por unidad de volumen de la interacción
electrón-ión, electrón-electrón e ión-ión divergen por separado debido al largo alcance del potencial
de Coulomb. Sin embargo, los términos divergentes se cancelan resultando en un valor nito para
UNED
la energía total del sistema por unidad de volumen. En un sólido cristalino la eliminación de dicha
divergencia se hace mediante la conocida técnica de sumación de Ewald, pero en el jellium dicha
divergencia se puede eliminar simplemente reordenando términos:
N
n (~r1 ) − n0 ] [b
n (~r2 ) − n0 ] − δ (~r1 − ~r2 ) n
Z
X 1 1 [b b (~r1 )
H
b = b 2i
p + d3~r1 d3~r2 (4.2)
2 2 |~r1 − ~r2 |
i=1
Así, la interacción electrostática clásica total es
[n (~r1 ) − n0 ] [n (~r2 ) − n0 ]
Z
clas
D E 1
Eint = Vext + WH [n] + Ered =
b d3~r1 d3~r2
2 |~r1 − ~r2 |
1
Recordemos de nuevo el abuso del lenguaje usual al denominar de igual manera al potencial y a la energía potencial.
107
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
con WH [n] la energía de Hartree que introdujimos en el capítulo 2. Es inmediato ver que el potencial
electrostático total que actúa sobre un electrón es
UNED
n (~r 0 ) − n0
Z
vtot (~r) = vext (~r) + vH (~r) = d3~r 0
|~r − ~r 0 |
con vH (~r) el potencial de Hartree (el potencial electrostático creado por la densidad electrónica).
3. Usando los operadores densidad en la representación de momentos, es fácil ver que
N
b =1 1 X 4π nh †
X ih i o
H b 2i +
p n (~
q ) − N
bδ ~ n (~
q ) − bδ ~ − N
N b
UNED
2 2V q2
b q~,0 b q~,0
i=1 q~
donde la suma está extendida a todos los momentos compatibles con las condiciones BvK y donde
hemos usado que b(~q = ~0) = N
n b , siendo N
b el operador número de partículas. Como estamos interesados
en estados de N electrones, N = N y, así
b
N
b =1 1 X 4π n † b + lı́m 4π n0
X o
H b 2i +
p n (~
q ) n (~
q ) − N
q2 2
b b
2 2V q~→~0 q
UNED
i=1 q~6=0
Sin embargo, el último término no contribuye a la energía por unidad de volumen en el límite
2
termodinámico, por lo que en de nitiva
N
b =1 1 X 4π n †
X o
H b 2i +
p n (~
q ) n (~
q ) − N
b (4.3)
q2
b b
2 2V
i=1 q~6=0
UNED
4. Nuestro objetivo va a ser evaluar y discutir las propiedades del estado fundamental y de respuesta
de este sistema modelo. Podemos anticipar que la densidad del estado fundamental |Ψ0 i va a ser
.
homogénea: n (~
r) = n0 = n, excepto para valores extremadamente bajos de la densidad media n0 (del
orden de rs ∼ 100). Así, salvo en este límite (conocido como cristal de Wigner ), la utilización del
término gas de electrones homogéneo (HEG de sus siglas en inglés) está justi cada para denotar al
estado de los electrones de valencia en el metal de jellium.
UNED
Para el gas homogéneo podemos de nir la polarización de espín como
N+ − N− n+ − n−
ζ= = (4.4)
N n
donde Nσ y nσ son el número de electrones y la densidad electróncia con espín σ , respectivamente. Para
densidades metálicas [rs ∈ (2, 6)] también podemos anticipar que el estado fundamental exhibirá
compensación en espines (ζ = 0, fase paramagnética ). Ahora bien, para densidades menores los
UNED
efectos cuánticos pueden favorecer una alineación de espínes (ζ 6= 0) y el sistema estaría en una fase
ferromagnética. El estudio de las propiedades del estado fundamental, incluyendo estas fases asociadas
al estado de espín, será el objeto de la siguiente sección.
2
En efecto:
2/3
q ·~
Z
1 4π n0 exp (−i~ r) 3
lı́m n0 = lı́m d3 ~
r = πn0 V −1/3
V q~→~0 q 2 2V q~→~0 r 4π
igual a cero en el límite termodinámico.
108
UNED
El gas homogéneo de electrones
Máster FSC. Funcional de la densidad
4.2. PROPIEDADES DEL ESTADO FUNDAMENTAL
4.2.a. Planteamiento
UNED
1. Abordemos en primer lugar las propiedades del estado fundamental del gas de electrones usando
el método de Hartree-Fock (en principio no restringido). Así, consideremos el determinante de Slater
(+) (+) (−) (−)
|Φ0 i formado por los orbitales HF ψi (~
κ ) = φi (~r) α (σ) y ψi (~
κ ) = φi (~r) α (σ) que obedecen
las ecuaciones HF:
Z
1~2 (σ) (σ) (σ) (σ)
+ vext (~r) + vH (~r) φi (~r) + d3~r 0 Σσσ r, ~r 0 φi ~r 0 = εi φi (~r)
− ∇ X ~
2
UNED
Nσ (σ) (σ) 0 ∗
n (~r) − n0
Z Z
1~2 (σ)
X φi (~r) φi (~r ) (σ) 0 (σ) (σ)
+ d3~r 0 φi (~r) − d3~r 0
− ∇ 0 0
φi ~r = εi φi (~r)
2 |~r − ~r | |~r − ~r |
j=1
Aquí, σ = +, − y Σσσ
X es el operador de Fock para los electrones con espín σ. La energía fundamental
exacta es, entonces, ! !
(σ) (σ)
X X
clas
E0 = TS + Eint + EX + EC (4.5)
UNED
σ σ
(σ)
donde TS es la energía cinética en la aproximación HF de los electrones con espín σ, el segundo
(σ)
término es la energía electróstática clásica total, EX es la energía de intercambio de los electrones
con espín σ y EC es la energía de correlación, ausente en la aproximación HF.
Si suponemos que la densidad electrónica del estado fundamental es uniforme e igual a n0 , los
términos clásicos se cancelan y
UNED
Z
1 ~ 2 (σ) (σ) (σ) (σ)
d3~r 0 Σσσ r, ~r 0 φi ~r 0 = εi φi (~r) , σ = +, −
− ∇ φi (~r) + X ~ (4.6)
2
son las ecuaciones HF para el gas de electrones. Naturalmente, esta suposición ha de comprobarse a
posteriori. Observese que en el método HF no hay ningún acoplamiento entre electrones con diferente
espín y, en este sentido, el gas homogéneo en la aproximación HF es la superposición de dos gases
independientes, cada uno formado por electrones con igual orientación del espín.
UNED
4.2.b. El gas de electrones libres
2. Resulta ilustrativo calcular primero la función de ondas del sistema y la energía fundamental bajo
la aproximación de Hartree (ΣX = 0) Así, los orbitales obedecen las ecuaciones
1 ~ 2 (σ) (σ) (σ)
− ∇ φi (~r) = εi φi (~r)
UNED
2
y el problema se reduce al del gas de fermiones libres. Los orbitales (normalizados sobre el volumen
V) son los correspondientes a una partícula libre, esto es, las ondas planas autoestados del momento
lineal:
1
θ~k (~r) = √ exp i ~k · ~r
V
2
con energías ε~k = ~k /2 ≡ k 2 /2 y ~k un momento compatible con las condiciones BvK. Nótese que
2
θ~k (~r) = 1/V , por lo que la densidad asociada es, en efecto, uniforme.
Así, si Nσ es el número de electrones con espín σ, los orbitales ocupados serán todos aquellos para
los que se cumpla
(σ)
k ≤ kF
109
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
(σ)
siendo kF el llamado momento de Fermi para los electrones con espín σ, que está dado por
(σ)
X
Θ kF − k = Nσ
UNED
~k
donde la suma está extendida a todos los momentos permitidos y Θ es la función de Heavyside (1 si
el argumento es positivo, 0 si es negativo).
Ya que V 0, la diferencia entre los momentos permitidos es muy pequeña, por lo que la suma
sobre ~k puede transformarse en una integral:
UNED
V
X Z
−→ d3~k
~k
(2π)3
De esta manera
V
Z
(σ) (σ)
X
Θ kF − k = Nσ ⇒ d3~k Θ kF − k = Nσ
~k
(2π)3
(σ)
UNED
La integral es simplemente el volumen en el espacio de momentos de una esfera de radio kF , por
tanto
(σ) 3
V 4π kF 1/3
(σ)
= Nσ ⇒ kF = 2 × 3π 2 nσ (4.7)
(2π)3 3
La energía cinética (igual a la energía total en esta aproximación de Hartree) de los electrones con
espín σ es, entonces
UNED
(σ)
k2 k2 kF
V V k2
Z Z
(σ) (σ) (σ)
X
3~
TS = Θ kF −k = d kΘ kF −k = dk 4π k 2
~k
2 (2π)3 2 (2π)3 0 2
y la energía cinética por partícula para partículas de espín σ será
(σ) 1/3
(σ) T 3 (σ) 2 3 × 2 × 3π 2 nσ 3 (σ)
tS = H = k = = εF (4.8)
10 F
UNED
Nσ 10 5
(σ) 2
(σ)
siendo εF = kF /2 la llamada energía de Fermi para los electrones con espín σ.
Si expresamos ahora nσ en términos de la densidad total y de la polarización en espín ζ:
1±ζ
n± = n,
2
UNED
tenemos que
(±) 1/3
kF = (1 ± ζ)1/3 3π 2 n ≡ (1 ± ζ)1/3 kF
1/3
con kF = 3π 2 n el momento de Fermi (a secas) del gas de electrones libre con densidad n.
(+) (−)
Naturalmente, para la fase paramagnética (ζ = 0) se cumple que kF = kF = kF . Con esta notación,
la energía cinética por partícula es
(+) (−)
TS + TS
tS (n, ζ) =
N
y
(1 + ζ)5/3 + (1 − ζ)5/3 3kF2 (1 + ζ)5/3 + (1 − ζ)5/3 3εF
tS (n, ζ) = = (4.9)
2 10 2 5
110
UNED
El gas homogéneo de electrones
Máster FSC. Funcional de la densidad
con εF = kF2 /2 la energía de Fermi (a secas) del gas de electrones libre con densidad n. Observamos que
para cualquier densidad la fase más estable es aquella en la que no hay polarización en espín (ζ = 0),
2
consecuencia lógica del carácter no interactuante de las partículas. En este caso tS (n, ζ = 0) = 3εF /5.
UNED
Si en lugar de usar εF usamos el radio de Wigner,
2/3
(1 + ζ)5/3 + (1 − ζ)5/3 3 (1 + ζ)5/3 + (1 − ζ)5/3 1, 105
9π 1
tS (n, ζ) = ' Ha
2 10 4 rs2 2 rs2
lo que permite estimar rápidamente el orden de magnitud de las energías cinéticas por partícula para
las densidades electrónicas de valencia características en física del estado sólido.
UNED
4.2.c. La solución Hartree-Fock
3. Pasemos ahora al resultado Hartree-Fock. Como sabemos, las ecuaciones HF son autoconsistentes,
por lo que lo primero que hemos de hacer es calcular el operador de Fock para una primera solución
(guess), la cual puede ser la formada por los orbitales libres de la descripción no interactuante. Para
éstos:
UNED
θ (~r1 )∗ θ (~r2 ) exp i~k · (~r2 − ~r1 )
−1
Z
(σ) ~ ~k (σ)
X
Σσσ r1 , ~r2 ) = −
X (~ Θ kF − k k = d3~k Θ kF − k
~k
|~r2 − ~r1 | (2π)3 |~r2 − ~r1 |
Observamos que Σσσ r1 , ~r2 ) depende únicamente de ~r2 −~r1 (y, de hecho, sólo de |~r2 − ~r1 |) consecuencia
X (~
de la homogeneidad e isotropía de la solución de partículas libres. Por tanto, el operador de Fock va
a ser diagonal en la base orbital de ondas planas, por lo que éstas van a ser también solución de las
ecuaciones HF. Comprobémoslo explícitamente:
UNED
~
!∗ Z ~0
Z
ei k·~r1 d3 ~q ei~q·(~r2 −~r1 ) ei k ·~r2
~k , ~k 0 = θ~ Σ (σ)
Σσσ
X
b σσ
X θ~k 0 = − d3~r1 d3~r2 √ Θ k F − q √
k
V (2π)3 |~r2 − ~r1 | V
Para efectuar las integrales hagamos un cambio de variables de jacobiano unidad: ~ =
~r = ~r1 − ~r2 ; R
(~r2 + ~r1 ) /2. Operando
UNED
Z 3~ ei~q·~r
d3 ~q
Z Z
~k , ~k 0 = − d R ei(~k 0 −~k)·R~ d3~r e−i ~k·~r/2 e−i ~k 0 ·~r/2
(σ)
Σσσ
X Θ kF − ~
q
V (2π)3 r
La integral sobre ~
R es simplemente δ~k,~k 0 por ortonormalidad de las ondas planas. Así:
~
d3 ~q ei(q~−k)·~r
Z Z
.
~k , ~k 0 (σ) ~ (σ)
Σσσ
X = δ~k,~k 0 ΣX k = −δ~k,~k 0 Θ kF − q d3~r =
(2π)3 r
UNED
d3 ~q
Z 4π
(σ)
= −δ~k,~k 0 Θ kF − q
(2π)3 ~q − ~k
2
donde hemos usado la transformada de Fourier del potencial Coulombiano
exp (−i~q · ~r)
Z
4π
w (q) = d3~r = 2
r q
Ya sólo basta evaluar la integral. Sin pérdida de generalidad supondremos que ~k = k ~ez y, en
coordenadas esféricas (r, θ, φ)
Z k(σ) Z +1 Z 2π
~k = − 1 1
F
(σ) 2
ΣX 2
q dq dγ dφ 2 2
2π 0 −1 0 q + k − 2kqγ
111
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
con γ = cos θ. Esta integral triple es analítica:
!2
UNED
k k
(σ)
1− (σ) (σ)
+1
(σ)
k kF k
ΣX ~k = − F 1 + ln F (4.10)
π k k
2 (σ) −1
(σ)
kF kF
(σ)
(obsérvese el escalamiento del operador de Fock en k/kF ).
4. De esta manera, las autoenergías de los orbitales HF son
UNED
(σ) k2 (σ)
ε~ = + ΣX ~k
k 2
y la energía de intercambio para los electrones con espín σ es (recuérdese el epígrafe 3 de la sección
2.2.b)
(σ) 1 X D b σσ θ~ = 1
E X (σ)
EX = θ~k Σ X k
ΣX ~k
UNED
2 (σ)
2 (σ)
k≤kF k≤kF
Susituyendo y operando,
(σ)
V
Z
(σ) 1 X (σ) ~ 3~
(σ)
(σ) ~ 3kF
EX = ΣX k = d k Θ kF − k ΣX k = ... = − Nσ
2 (σ)
2 (2π)3 4π
k≤kF
UNED
por lo que la energía de intercambio por partícula para los electrones con espín σ es:
(σ) (σ)
(σ) E 3k
εX = X =− F (4.11)
Nσ 4π
Siguiendo un procedimiento similar al que vimos para tS , podemos escribir la energía de intercambio
por partícula para un gas con densidad total n y polarización ζ. Operando:
UNED
(1 + ζ)4/3 + (1 − ζ)4/3 3kF
εX (n, ζ) = − (4.12)
2 4π
o bien, en función del radio de Wigner y la polarización ζ:
1/3
(1 + ζ)4/3 + (1 − ζ)4/3 3 (1 + ζ)4/3 + (1 − ζ)4/3 0, 458
9 1
εX (n, ζ) = − =− Ha (4.13)
2 4 4π 2 rs 2 rs
UNED
5. Uniendo los dos resultados, la energía media por partícula del gas homogéneo en la aproximación
HF es
2/3 1/3
(1 + ζ)5/3 + (1 − ζ)5/3 3 (1 + ζ)4/3 + (1 − ζ)4/3 3
9π 1 9 1
εtot (n, ζ) = −
2 10 4 rs2 2 4 4π 2 rs
En particular
2/3 1/3
3 9π 1 3 9 1
εtot (n, 0) = 2
− ; (ζ = 0, fase paramagnética)
10 4 rs 4 4π 2 rs
2/3 1/3
3 9π 1 3 9 1
εtot (n, ±1) = 2
− ; (ζ = ±1, fase ferromagnética)
10 2 rs 4 2π 2 rs
112
UNED
El gas homogéneo de electrones
Máster FSC. Funcional de la densidad
UNED
εWRW P+DH
ζ
UNED
ζ
UV DX
UNED
Figura 4.1. Energía por partícula del gas homógeneo bajo la aproximación HF (trazos) comparada con los resultados
exactos Monte Carlo (línea continua). Obsérvese la transición de fase en el modelo HF a una densidad rs ' 5, 45, que en
los cálculos Monte Carlo se produce a una densidad mucho más baja.
valores que representamos en la gura 4.1. Observamos que para rs grande (densidades pequeñas)
el intercambio predomina, mientras que para rs pequeño (densidades grandes) es la energía cinética
UNED
quien predomina. A su vez la energía de intercambio es menor (más negativa) si el gas está polarizado
en espín (ζ = 1), lo que contrasta con el hecho de que la energía cinética es menor si el gas no está
polarizado (ζ = 0). Esto se traduce en el hecho de que exista una transición de fase de primer orden
entre las fases paramagnética y ferromagnética a una densidad dada por εtot (n, 0) = εtot (n, ±1). Esta
densidad crítica es
2 1 − 2−2/3
rscrit = 32/3 π 4/3 1/3 ' 5, 45
5 2 −1
UNED
de modo que si rs > rscrit la fase más estable (la de menor energía) es la ferromagnética, mientras que
si rs < crit
rs la más estable es la paramagnética. Obsérvese que este resultado es una ilustración del
principio de Hund (el intercambio tiende a favorecer la alineación de espines).
6. Para nalizar, obtengamos la función de distribución de pares Hartree-Fock. A partir de la
expresión para un determinante de Slater
Nσ 2
UNED
1 1X X (σ) (σ)
ρ2 (~r1 , ~r2 ) = n (~r1 ) n (~r2 ) − φi (~r1 )∗ φi (~r2 )
2 2 σ
i=1
tenemos que
2
d3~k (σ)
Z
1 1X
~
ρ2 (~r1 , ~r2 ) = n2 − 3 Θ k F − k eik·(~r2 −~r1 )
2 2 σ (2π)
La integral sobre ~k es, de nuevo, analítica, lo que conduce a
2
(σ) (σ) (σ)
1 1 X 2 sin kF r12 − kF r12 cos kF r12
ρ2 (~r1 , ~r2 ) = n2 1 − 2 9nσ
3
2 n σ
(σ)
kF r12
113
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
UNED
ζ
K UU
ζ
UNED
N)U
UNED
Figura 4.2. Función de correlación de pares del gas homógeneo bajo la aproximación HF.
De esta manera, para la fase paramagnética y ferromagnética
2
1 9 sin (kF r12 ) − (kF r12 ) cos (kF r12 )
UNED
ρ2 (~r1 , ~r2 ) = n2 1 − 3 ; si ζ = 0
2 2
(σ)
kF r12
!2
sin 21/3 kF r12 − 21/3 kF r12 cos 21/3 kF r12
1 2 9
ρ2 (~r1 , ~r2 ) = n 1 − ; si ζ = ±1
2 4 (kF r12 )3
UNED
Observemos, de nuevo, el escalamiento de la función de pares y, además, el hecho de que dependa
de r12 = |~r2 − ~r1 |, consecuencia de la homogeneidad e isotropía del gas homogéneo. En la gura
4.2 representamos la función de correlación de pares h (~r1 , ~r2 ) = −1 + 2ρ2 (~r1 , ~r2 ) /n2 . Nótese que la
correlación entre partículas es prácticamente nula si kF r & 4.
4.2.d. La energía de correlación RPA
UNED
7. La aproximación HF desprecia cualquier efecto de correlación coulombiana entre los electrones ya
que, recordemos, en el método HF sólo tenemos en cuenta los efectos debidos al principio de exclusión
(dos electrones no pueden estar en el mismo estado). Como segunda consecuencia, no contempla
correlación alguna entre electrones de distinto espín. La inclusión de efectos de correlación coulombiana
conllevará a que la probabilidad de encontrar dos electrones a una distancia corta entre sí sea menor,
lo que conduce a una disminución de la función de densidad de pares a cortas distancia y a un valor
esperado de la energía de interacción menor. De esta manera, la energía HF es una cota superior a la
energía real (resultado que también sabemos debido al carácter variacional del método HF).
El cálculo aproximado de la energía de correlación del gas homogéneo se remonta a los años 30
del siglo pasado, siendo Wigner el primero que estimó su valor mediante una serie de argumentos
semiheurísticos (véase el texto de Pines). Si nos ceñimos a lo que hemos visto en los temas anteriores,
la manera más sencilla de obtener una aproximación a la energía de correlación sería el método MP2
114
UNED
El gas homogéneo de electrones
Máster FSC. Funcional de la densidad
(perturbaciones a segundo orden). Para el gas de electrones, recordemos que
N
b =1 1 X 4π n †
UNED
X o
H b 2i +
p n (~
q ) n (~
q ) − b = Tb + Ŵ 0
N
q2
b b
2 2V
i=1 q~6=0
donde Ŵ 0 es la interacción entre los electrones pero sustrayendo todos los términos clásicos divergentes.
Esta igualdad implica que la formulación MP de la teoría de perturbaciones es equivalente a considerar
que Ŵ 0 es el hamiltoniano perturbador sobre el gas de electrones libre. La corrección a primer orden
es la energía de intercambio, mientras que las correcciones a orden dos y superior nos darán la energía
UNED
de correlación. Sin embargo para el gas homogéneo (de nuevo véase el texto de Pines para una
demostración detallada), la energía de correlación MP2 por partícula está dada por
εMP2
c (n, ζ) = ∞ (4.14)
Ello no quiere decir que la energía de correlación por partícula sea in nita para el gas de electrones.
Sino que la corrección a segundo orden de teoría de perturbaciones tiene un término divergente, que
ha de cancelarse con las correcciones a orden superior. De hecho, si tomamos una corrección de un
UNED
orden dado ésta será in nita: es la suma de todas las correcciones la que da el valor correcto de la
energía de correlación:
∞
X
εc (n, ζ) = εMPs
c (n, ζ)
s=2
por lo que es evidente que la teoría de perturbaciones usual no es útil para el gas de electrones.
8. Una alternativa consiste en la utilización del teorema de uctuación disipación. En principio no
UNED
cubre enteramente nuestros propósitos ya que sólo nos permite evaluar el valor esperado de la energía de
interacción, mientras que la energía de correlación contiene también correcciones a la energía cinética
HF, que no es exacta. Este problema se puede solventar mediante una transformación adiabática, que
merece la pena discutir con algo de detalle para el gas homogéneo, dejando su generalización para
cuando estudiemos la implementación de la teoría del funcional de la densidad.
Así, consideremos el gas homogéneo con Hamiltoniano Tb + Ŵ 0 , energía fundamental E0 y estado
fundamental |Ψ0 i. De namos el hamiltoniano dependiente de un parámetro λ
UNED
b (λ) = Tb + λŴ 0
H ; λ ∈ [0, 1]
E
(λ)
y sean E0 (λ) y Ψ0 las energías y estados fundamentales respectivos. Naturalmente
E E
(1) (0)
E0 (1) = E0 ; Ψ0 = |Ψ0 i y E0 (0) = TS ; Ψ0 = |Φ0 i
UNED
donde TS es la energía cinética del gas de electrones libre cuyo estado fundamental es el determinante
de Slater |Φ0 i. Si ahora usamos el teorema de Hellman-Feynman tenemos que
Z 1 Z 1 Z 1
dE0 (λ) D
(λ) ∂H
b (λ) (λ) E D
(λ) (λ)
E
E0 (1) − E0 (0) = dλ = dλ Ψ0 Ψ0 = dλ Ψ0 Ŵ 0 Ψ0
0 dλ 0 ∂λ 0
y como, por de nición para el gas homogéneo, E0 (1) = E0 = TS + EX + EC :
Z 1 hD E D Ei
(λ) (λ) (0) (0)
EC = dλ Ψ0 Ŵ 0 Ψ0 − Ψ0 Ŵ 0 Ψ0
0
donde hemos usado D E
(0) (0)
EX = hΦ0 | Ŵ 0 |Φ0 i = Ψ0 Ŵ 0 Ψ0
115
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
Ahora bien, si sustituimos Ŵ 0 por W
c dentro del corchete, los términos divergentes clásicos que aparecen
en ambos valores esperados se cancelan. Llegamos en de nitiva a que
UNED
Z 1 hD E D Ei
(λ) c Ψ(λ) − Ψ(0) W
c Ψ(0)
EC = dλ Ψ0 W 0 0 0 (4.15)
0
y si usamos el teorema de uctuación-disipación
∞ 1
χλ (~r1 , ~r2 ; ω) − χ0 (~r1 , ~r2 ; ω)
Z Z Z
dω
EC = − Im d3~r1 d3~r2 dλ (4.16)
0 2π 0 |~r1 − ~r2 |
UNED
donde χλ es la función respuesta del gas homogéneo de electrones pero que interactuan entre sí con
un potencial λ/r12 y, por tanto, χ0 es la función respuesta de partículas libres. Usando una notación
operacional simbólica, dicha corrección es
Z ∞ Z 1
dω
EC = − Im tr χλ (ω) − χ
dλ [b b 0 (ω)] w
b (4.17)
0 2π 0
donde tr es la traza, que ha de entenderse como
UNED
Z
b = d3~r A (~r, ~r)
tr A
y w
b, como siempre, es el potencial de Coulomb w (~r1 , ~r2 ) = 1/r12 .
9. Naturalmente hemos de evaluar la función respuesta para distintos valores de λ. En la
aproximación RPA que vimos en el tema anterior se tiene que
UNED
χ
b λ (ω) = χ
b 0 (ω) + χ
b 0 (ω) λwb
bχλ (ω)
(basta sustituir la interacción real w
b por la escalada λw
b). Operando,
1
χλ (ω) − χ
[b b 0 (ω)] w
b= b−χ
b 0 (ω) w
χ b 0 (ω) w
b
bI − λb
χ0 (ω) w
b
pero, en términos de operadores
UNED
Z 1
1
b = − ln bI − A
b + 1Ab2 + 1 A
b3 + ...
dλ A b =A
0 bI − λA
b 2 3
por lo que, en de nitiva, la energía de correlación en la aproximación RPA es
Z ∞ ∞ Z ∞
dω h b i X1 dω
RPA
EC = Im tr ln I − χ
b 0 (ω) w
b +χ b = − Im
b 0 (ω) w tr [(b b s]
χ0 (ω) w) (4.18)
2π s 2π
UNED
0 s=2 0
Este resultado, que hemos deducido para el gas homogéneo, tiene una validez general. Observemos
que en
RPA
EC aparecen términos que van como b2 , w
w b3 , w
b4 , etc. Así, este resultado puede entenderse
como una suma parcial de los in nitos términos de la serie perturbativa que nos da la energía de
correlación exacta. Como veremos en la próxima sección, la función respuesta χ
b 0 (ω) para el gas
homogéneo admite una expresión analítica y las integrales que aparecen en (2.15) se pueden realizar
numéricamente sin ningún esfuerzo.
10. Las energías de correlación RPA son, sin embargo, demasiado grandes (en valor absoluto)
comparadas con los valores exactos obtenidos mediante simulaciones Monte Carlo por Ceperley y
Adler alrededor de 1980 (véase la gura 4.3). Ello implica que la función de correlación de pares
que se obtendría a partir del teorema de uctuación-disipación con la respuesta RPA es negativa en
exceso para distancias relativas pequeñas. Así, aunque la RPA contemple mecanismos de correlación
116
UNED
El gas homogéneo de electrones
Máster FSC. Funcional de la densidad
UNED
40&
ε& P+DH
53$
UNED
UV DX
UNED
Figura 4.3. Energía de correlación RPA por partícula para el gas homgéneo sin polarización de espines (trazos) comparada
con los resultados exactos Quantum Monte Carlo. Obsérvese que la RPA sobrestima en valor absoluto la energía de
correlación en torno a un 50 %.
importantes, éstos no son su cientes para obtener energías enteramente satisfactorias. Recordemos
que la RPA considera que la respuesta de un sistema de electrones es clásica, en el sentido de que
UNED
responden como si fuesen partículas independientes frente al campo total (exterior + inducido). Los
efectos mecanocuánticos en la respuesta, que están ausentes en la RPA, son así los responsables de
estos errores.
Lo anterior no es óbice para considerar a la energía RPA como una mejora cualitativamente
fundamental de los resultados Hartree-Fock. A partir de una hipótesis sencilla relacionada con la
forma en la que el sistema de N electrones responde a una perturbación exterior hemos incorporado
efectos de correlación coulombiana que, sin ser completamente satisfactorios, incluyen muchos
UNED
efectos mecanocuánticos que están completamente fuera del alcance de la venerable aproximación
de Hartree-Fock. Los aparentemente malos resultados de la aproximación RPA no deben minimizar
su importancia.
UNED 117
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
4.3. EXCITACIONES EN EL GAS HOMOGÉNEO DE ELECTRONES
4.3.a. La función respuesta de partículas independientes
UNED
1. En el tema anterior vimos que la función respuesta χ (~r1 , ~r2 ; ω) era una propiedad esencial a la
hora de describir las propiedades de excitación de un sistema de N electrones cuando sobre él actuaba
bien un campo longitudinal bien una campo EM de longitud de onda larga. En la aproximación RPA,
la función respuesta es solución de la ecuación integral
χ
b (ω) = χ
b 0 (ω) + χ
b 0 (ω) wb
bχ (ω)
UNED
(seguimos usando la notación operacional).
Además de las aproximaciones implícitas a la RPA, existe una segunda di cultad ¾cuál es la mejor
forma de construir el hamiltoniano efectivo H
b0 de particulas independientes para evaluar la función
respuesta χ
b 0 (ω)? Lo inmediato sería usar el hamiltoniano HF, sin embargo por motivos que veremos
más adelante esta es una muy mala idea en sistemas extensos (como el gas de electrones). La teoría
del funcional de la densidad nos dará una pista en este sentido para el caso más general, pero en la
situación que nos ocupa (el gas de electrones) la respuesta es bien sencilla.
UNED
Para el HEG, las partículas independientes son los autoestados del momento lineal. La duda
consiste en elegir entre
PN b H
b 0 = Tb (es decir, el sistema de partículas no interactuantes) o H
b 0 = Tb +
i=1 ΣX (i) (el hamiltoniano HF). Por lo que hemos dicho (aunque sin explicación), la segunda es
una mala opción, por lo que la función respuesta χ
b 0 (ω) será la correspondiente al gas de electrones
libres. Esta elección puede justi carse usando argumentos de Teoría de Perturbaciones de Muchos
Cuerpos (MBPT de sus siglas en inglés) algo técnicos pero relacionados con el mero hecho de que
para el HEG, H c0,
b = Tb + W donde recordemos que c0
W es la interacción excluidos términos clásicos. La
UNED
idea básica es que, puesto que la respuesta interactuante está asociada a la existencia del término c0
W
PN b
(que implícitamente contiene el operador de Fock i=1 ΣX (i)), es natural escoger como respuesta no
interactuante la del hamiltoniano H
b0 = H c0.
b −W
2. Restringiéndonos al gas de electrones sin polarización de espín y recordando lo visto en la
subsección 3.4.e, la función respuesta para el gas de electrones libres está dada por
XX θ∗q~1 (~r1 ) θq~1 (~r2 ) θq~2 (~r1 ) θ∗q~1 (~r2 )
χ0 (~r1 , ~r2 ; ω) = 2 × fq~1 − fq~2
UNED
ω − εq~2 − εq~1 + iη
q~1 q~2
√
con θq~ (~r) = exp (i~q · ~r) / V , εq~ = q 2 /2 y fq~ es el factor de ocupación (1 si q ≤ kF y 0 si q > kF ). Si
sustituimos
2 XX exp [i (~q1 − ~q2 ) · (~r2 − ~r1 )]
χ0 (~r1 , ~r2 ; ω) = f q
~ − f q
~
V2 1 2
q 2 − q12
q~1 q~2 ω− 2 + iη
2
Vemos inmediatamente que χ0 (~r1 , ~r2 ; ω) sólo depende de ~r2 − ~r1 por la homogeneidad del sistema.
UNED
Esto también re eja de que no hay efectos de campo local, por lo que en la representación recíproca
χ
b 0 (ω) será diagonal. Explícitamente, haciendo el cambio ~k = ~q1 − ~q2 , ~q = ~q2
2 X −i~k·~r1 X f~k+~q − fq~ ~
χ0 (~r1 , ~r2 ; ω) = e 2 eik·~r2 =
V2
~k q~ q 2 − ~k + ~q
ω− + iη
2
f − f
1 X X ~
−ik·~
2 X ~
k+~
q q
~
0
~
= e r1
2 δ~k,~
k 0 eik ·~r2
V V
q 2 − ~k + ~q
~k ~k0
q~
ω− + iη
2
118
UNED
El gas homogéneo de electrones
Máster FSC. Funcional de la densidad
Si ahora recordamos la relación entre la representación en espacio real y espacio recíproco,
inmediatamente llegamos a que
UNED
2X f~k+~q − fq~ .
χ0 ~k, ~k 0 ; ω = δ~k,~k0 2 = δ~k,~k0 χ0 ~k; ω
V
q~ q 2 − ~k + ~q
ω− + iη
2
donde, haciendo el paso al continuo en la suma sobre los momentos ~q,
1
Z Θ kF − ~k + ~q − Θ (kF − |~q|)
UNED
χ0 ~k; ω = 2 × d3 ~q (4.19)
(2π)3 q 2 − ~k + ~q
2
ω− + iη
2
La anterior expresión es razonable, ya que recordemos que χ(~k, ω) está asociada a excitaciones que
transmiten un momento ~k al sistema. Es por ello por lo que en (4.19) aparecen todos los términos
relacionados con la excitación de un electron en un estado θq~ ocupado a otro θ~k+~q desocupado.
UNED
Sorprendentemente, la integral de (4.19) es analítica, como demostró Lindhard hace ya más de
medio siglo. El resultado (función respuesta del gas de electrones libres o, simplemente, función de
Lindhard ) es
kF ω + iη ω + iη
χ0 (k, ω) = − k− η = 0+
e e
2e
k+Ξ e
k+ +Ξ e ; (4.20)
4π 2 e
k k
e k
e
UNED
donde
1 2 z+2
Ξ (z) = 1 − z ln
4 z−2
es una función de variable compleja y e e = 2ω/kF2 = ω/εF
k = k/kF , ω son momentos y frecuencias
escalados. Nótese que debido a la discontinuidad de la función logarítmica en el plano complejo, la
inclusión del término i0+ es esencial. A su vez, la respuesta sólo depende de |~k|, re ejo de la isotropía
del gas homogéneo.
UNED
4.3.b. La aproximación RPA
3. De acuerdo con la RPA, la función respuesta del gas de electrones es
χ0 (q; ω)
χ (q; ω) =
4π
1 − 2 χ0 (q; ω)
UNED
q
donde hemos tenido en cuenta el carácter diagonal en el espacio recíproco de χ
b 0 (ω) y w
b. Como
consecuencia inmediata, la función dieléctrica y su inversa son:
4π 4π 1
(q; ω) = 1 − χ (q; ω) ; −1 (q; ω) = 1 + χ (q; ω) = (4.21)
q2 0 q2 4π
1 − 2 χ0 (q; ω)
q
4. Centrémosnos en la respuesta a campos longitudinales, para los que hay varios aspectos
interesantes que merecen destacarse. En primer lugar, si sobre el sistema actúa un potencial
perturbador φext (~q, ω), tenemos que
φtot (~q, ω) = φext (~q, ω) + φind (~q, ω) = −1 (q, ω) φext (~q, ω) (4.22)
119
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
donde φind (~q, ω) es el potencial creado por la variación de la densidad inducida por el potencial exterior
(relación que ya habíamos visto en el caso general en el epígrafe 22 de la subsección 3.4.f ). La función
−1 (q, ω) es compleja, lo que implica que va a existir un desfasaje entre el campo exterior y el total.
UNED
Ya que por otro lado la parte imaginaria de −
−1 (q, ω) es proporcional a la probabilidad de que se
produzca una excitación de momento ~q y frecuencia ω , el desfasaje entre campo exterior y campo total
3 −1 (q, ω) ' 0
está asociado a efectos de absorción de energía por el gas de electrones. A su vez, si
tenemos que el campo total es prácticamente nulo y decimos que el campo exterior ha sido apantallado
por el sistema electrónico: el campo exterior ha inducido una densidad de carga que crea a su vez un
campo que se opone al campo exterior, compensándolo. Por el contrario, si −1 (q, ω) ' 1 tenemos que
el campo exterior apenas ha afectado al sistema de electrones y, como consecuencia, el apantallamiento
UNED
es prácticamente nulo. En general, el gas de electrones es incapaz de apantallar campos de longitud
de onda muy corta (momento grande) y/o frecuencias muy altas.
Como ilustración, consideremos el gas homogéneo con rs = 3,93 (valor correspondiente a la
densidad de electrones de valencia del sodio) y analicemos la función dieléctrica inversa −1 (q, ω)
particularizada a frecuencia nula y a una frecuencia ω = εF = 0,119 Ha ( gura 4.4). Observamos que
−1 (q, 0) es real, lo cual es lógico ya que un campo estático es incapaz de excitar el sistema. A su
vez, para momentos bajos (longitudes de onda larga), −1 (q, 0) es cercano a cero lo que implica un
UNED
apantallamiento muy efectivo. Sin embargo para momentos altos (longitudes de onda corta) −1 (q, 0)
es cercano a la unidad (apantallamiento poco efectivo). Por otro lado, −1 (q, εF ) tiene una estructura
algo más complicada, con una parte imaginaria no nula dentro de un cierto rango de momentos
(posiblidad de que existan excitaciones de momento q y energía εF ), pero las tendencias principales
son las mismas. Obsérvese que en este caso la función dieléctrica es negativa para momentos bajos:
la redistribución de carga sobreapantalla al campo exterior y el campo eléctrico total tendría sentido
opuesto.
5.
UNED
En segundo lugar repitamos que − Im −1 (q, ω) es proporcional a la probabilidad de que se
produzca una excitación de momento ~q y frecuencia ω . En nuestra aproximación RPA la forma más
obvia de que se produzca dicha excitación es mediante la promoción de un electrón de un estado
ocupado con momento ~k a otro desocupado con momento ~k + ~q y energía |~k + ~q|2 /2 = ω + k 2 /2, esto
es, mediante una excitación tipo electrón-hueco.
Si q ≤ 2kF es posible realizar una excitación con momento ~q y frecuencia prácticamente nula
simplemente promocionando un electrón de la super cie de Fermi (los estados con módulo del momento
UNED
igual a kF ) a otro punto de la super cie de Fermi. Sin embargo, si q > 2kF existe una cota inferior a
las energías de excitación: la excitación de menor energía se produce promoviendo un electrón de la
super cie de Fermi cuyo momento sea opuesto a ~q. La energía de esta excitación mínima es
(q − kF )2 kF2 1 2
ω= − = q − 2qkF
2 2 2
Por otro lado, independientemente del módulo de ~q, la excitación de máxima energía se produce para
UNED
un electrón en la super cie de Fermi cuyo momento es paralelo a ~q, en cuyo caso la energía de la
excitación es
(q + kF )2 kF2 1 2
ω= − = q + 2qkF
2 2 2
En de nitiva, las energías de excitación electrón-hueco con momento ~q permitidas son
1 2
ω ∈ 0, q + 2qkF si q ≤ 2kF
2
1 2 1 2
ω ∈ q − 2qkF , q + 2qkF si q > 2kF
2 2
3
Este no es un resultado cuántico, sino uno muy general como se mani esta (por ejemplo) en el oscilador armónico clásico
forzado y amortiguado.
120
UNED
El gas homogéneo de electrones
Máster FSC. Funcional de la densidad
5H∈
UNED
∈ Tω
,P∈
ω
UNED
5H∈
∈ Tω
UNED
,P∈
ω ε )
TN)
Figura 4.4. Función dieléctrica inversa para el HEG (rs = 3,93) particularizada a frecuencias ω=0 y εF . Obsérvese que
UNED
en ambos casos el gas es incapaz de apantallar campos con momento q & 4kF y el rango de momentos para los que se
producen excitaciones de energía ω.
y para este rango de valores de la frecuencia, Im −1 (q, ω) 6= 0.
Para ilustrar este hecho, en la gura 4.5 representamos −1 (q, ω) para momentos q = 2,5kF , 1,5kF y
0,5kF . En el primer caso, las excitaciones permitidas están en el rango [1,25εF , 11,25εF ], en el segundo
UNED
están en el intervalo [0 , 5,25εF ] y, en el tercero, dentro del intervalo [0 , 1,25εF ]. Sin embargo, en este
último caso
−1 (q, ω) es singular para ω ' 2,05ε lo que se mani esta en la aparición de un pico
F
muy bien de nido (de hecho es una delta de Dirac en el gas homogéneo) en la parte imaginaria de
−1 (q, ω). Esto indica que, bajo ciertas circunstancias, es posible una excitación de momento ~q con una
energía fuera del rango permitido para excitaciones electrón-hueco individuales o, en otras palabras,
una excitación que no tiene su equivalente en el límite no interactuante. Este es un modo colectivo o
plasmón.
6. La existencia de esta excitación ya había sido anticipada en el capítulo 3 en base a un argumento
UNED
semifenomenológico basado en la función de pérdida. Podemos ahora profundizar un poco más en ella
si consideramos la relación entre campo exterior y campo total, pero escrita en términos de la función
dieléctrica:
(q, ω) φtot (~q, ω) = φext (~q, ω) ⇔ χ−1 (q, ω) δn (~q, ω) = vext (~q, ω)
Si para un cierto momento ~q existe una frecuencia ω p (q) tal que (q, ω p (q)) ' 0, es posible la existencia
de un campo nito en el sistema en ausencia de un campo exterior (o como resultado de una respuesta
a un campo exterior prácticamente nulo). Este campo total es el creado por una densidad de carga
δn (~q, ω p (q)) autosostenida que constituye el plasmón. Por tanto, el plasmón es un estado excitado del
gas de electrones igual a un modo oscilatorio de la densidad de momento ~q y frecuencia ω p (q), dada
ésta por la solución de la ecuación
(q, ω p (q)) = 0 (4.23)
121
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
5H∈
UNED
× ,P∈
T N)
UNED
5H∈
∈ Tω
,P∈
T N)
UNED
5H∈
× ,P∈
UNED
T N)
ωε )
Figura 4.5. Función dieléctrica inversa para el HEG (rs = 3,93) particularizada a momentos q/kF = 0,5 , 1,5 y 2,5.
−1
UNED
Obsérvese el plasmón para q = 0,5kF (ausente en los otros dos casos) que se mani esta en la asíntota de Re y la delta
−1
de Dirac en Im . También es interesante ver los rangos para los que las excitaciones electrón-hueco están permitidas.
−5 −1
Se ha usado un broadening η = 10 en las expresiones analíticas para .
En el caso que nos ocupa (rs = 3,93), un análisis detallado de la función dieléctrica indica que (2.20)
sólo admite una solución si q . 0,95kF . Para este momento límite, ω p (q) ' 2,8εF ' 8,8 eV. Por otra
parte, en el límite q→0 puede probarse que
UNED
ω 2p
(0, ω) = 1 − (4.24)
ω2
donde s
3 √
ωp = 3
= 4πn (4.25)
rs
es la llamada frecuencia de plasma del gas homogéneo, ya que ω p es la frecuencia del plasmón en el
límite de longitud de onda muy larga. Para rs = 3,93, se tiene que ω p = ω p (0) ' 1,88εF ' 0,216 Ha '
5,9 eV. Así, para este gas homogéneo los posibles plasmones tienen frecuencia en el rango (5,9 , 8,8)
eV y momentos q ∈ [0 , 0,95kF ].
7. Vemos así que las excitaciones en el HEG resultado de un campo longitudinal débil independiente
del espín son de dos tipos: excitaciones electrón-hueco, en las que grosso modo un electrón de momento
122
UNED
El gas homogéneo de electrones
Máster FSC. Funcional de la densidad
k ≤ kF se excita a un estado con momento ~k + ~q tal que ~k + ~q > kF , y excitaciones colectivas
o plasmones, en las que el campo exterior entra en resonancia con un modo de oscilación del gas
UNED
produciéndose una onda de carga que se propaga a lo largo del metal con momento ~q y frecuencia
ω p (q). Naturalmente estamos usando un lenguaje gurado: para el sistema interactuante real, la
correlación coulombiana hace que hablar de partículas independientes haya de entenderse siempre de
manera heurística. Sin embargo esto abre el camino a la de nición de excitación elemental en un
sistema extenso, que permite dotar de rigor formal a este tipo de a rmaciones. Puesto que volveremos
a este punto, baste de momento esta descripción semiheurísitca de las excitaciones en un metal.
8. Por último, merece la pena describir brevemente las propieades ópticas del gas homogéno (aunque
no sean muy fascinantes, la verdad). La función dieléctrica macroscópica es
UNED
ω 2p
M (ω) = (0, ω) = 1 − (4.26)
ω2
De esta manera, el coe ciente de extinción y el índice de refracción son, respectivamente,
r
ω 2p r 0 si ω < ω p
− 1 si ω < ω p
UNED
β (ω) = ω2 ; ν (ω) = ω2 (4.27)
1− p si ω > ω p
0 si ω > ω p
ω2
lo que indica que el gas homogéneo es opaco [transparente] para radiación electromagnética de
frecuencia menor [mayor] que la frecuencia de plasma. Así, un material cuya estructura electrónica
pueda describirse de manera aproximada con el modelo de Sommerfeld, re ejará la luz que incida
sobre el con frecuencia menor que ω p . Siguiendo con el ejemplo del sodio, la longitud de onda de la
radiación EM correspondiente es λcrit ' 2 × 10−7 m, que está en el rango del ultravioleta.
UNED
UNED
UNED 123
Máster FSC. Funcional de la densidad
UNED
UNED
UNED
UNED
UNED
UNED
UNED
UNED
Máster FSC. Funcional de la densidad
Capítulo 5
EL TEOREMA DE HOHENBERG-KOHN
UNED
OBJETIVOS:
1. Conocer someramente el concepto matemático de funcional.
UNED
2. De nición formal del potencial químico de un sistema electrónico.
3. Exposición del teorema básico de la Teoría del Funcional de la Densidad.
4. Breve análisis de funcionales de la energía con dependencia directa en la densidad.
LECTURAS COMPLEMENTARIAS:
UNED
1. Funcionales (concepto matemático):
a) Density Functional Theory of Atoms and Molecules (Parr y Yang): Apéndice A
2. El potencial químico:
a) Parr y Yang: Capítulo 4
UNED
3. El Teorema de Hohenberg y Kohn:
a) Parr y Yang: Capítulo 3
4. El modelo de Thomas-Fermi:
a) Parr y Yang: Capítulo 6
ORIENTACIONES AL PLAN DE TRABAJO:
UNED
En los capítulos anteriores hemos resumido algunos aspectos importantes referentes a la
descripción del estado fundamental y de los estados excitados de un sistema electrónico. En este
capítulo y en los dos siguiente nos vamos a centrar exclusivamente en el estado fundamental.
Veremos que su obtención no precisa resolver la ecuación de Schrödinger N -electrónica, sino que
se pueden calcular las propiedades de dicho estado fundamental usando una teoría alternativa:
la teoría del funcional de la densidad (DFT, de sus siglas en inglés).
UNED
A título orientativo, se espera que dedique 10 horas a este capítulo desglosadas como sigue:
• 4 horas de lectura detallada y estudio del presente material didáctico.
• 6 horas para el a anzamiento de los contenidos, incluyendo consultas y las lecturas
complementarias [optativas].
GLOSARIO:
Funcional
Derivada variacional
Funcional de Weizsäcker
125
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
Hamiltoniano convexo
Energía de ionización
UNED
A nidad electrónica
Potencial químico
Funcional de Hohenberg-Kohn
Teorema de Hohenberg-Kohn
UNED
Teoría del funcional de la densidad (DFT)
Ecuación de Euler-Lagrange
Modelo de Thomas-Fermi
Funcional cinético de Thomas-Fermi
Funcional de la energía de intercambio de Dirca
UNED
Funcional cinético de Thomas-Fermi-Weizsäcker
UNED
UNED
UNED
126
UNED
El teorema de Hohenberg-Kohn
Máster FSC. Funcional de la densidad
5.1. FUNCIONALES DE LA DENSIDAD
1. F R3
Consideremos el conjunto formado por todas las funciones de nidas en el espacio
UNED
R3 . 3
tridimensional Un funcional
es cualquier aplicación que a cada función de F R le hace
1
corresponder un escalar (en general, complejo). Así, dado un sistema de N electrones, un funcional
de la densidad electrónica g[n] es toda aplicación que hace corresponder a cada posible densidad física
n (~r) un escalar (real o complejo).
Nótese que nosotros ya conocemos dos funcionales de la densidad. El primero es el valor esperado
de la energía potencial, ya que conocido el potencial exterior, vext (~r), se tiene que
UNED
Z
.
D E
Vbext = Vext [n] = d3~r n (~r) vext (~r)
(de hecho, éste es un funcional lineal de la densidad). El segundo es la energía electrostática clásica
de Hartree: Z
1 n (~r1 ) n (~r2 )
WH [n] = d3~r1 d3~r2
2 |~r1 − ~r2 |
UNED
2. Consideremos un funcional cualquiera de la densidad g[n]. Si ésta última sufre una variación
δn (~r), el valor que toma el funcional cambiará en una cantidad
δg = g[n + δn] − g [n]
Siguiendo un procedimiento muy similar al conocido en Analisis Matemático se tiene que
δ 2 g [n]
Z Z
δg [n] 1
UNED
3 3 3
δg = d ~r δn (~r) + d ~r1 d ~r2 δn (~r1 ) δn (~r2 ) + ... (5.1)
δn (~r) 2 δn (~r1 ) δn (~r2 )
donde la función δg [n] /δn (~r) se conoce como derivada variacional del funcional g [n]. Análogamente
δ 2 g [n]
δ δg [n] δ δg [n]
= =
δn (~r1 ) δn (~r2 ) δn (~r1 ) δn (~r2 ) δn (~r2 ) δn (~r1 )
es la segunda derivada variacional del funcional g [n]. Nótese como la función δg [n] /δn (~r) depende
UNED
de la densidad. Por tanto, si jamos la posición ~r, la derivada variacional δg [n] /δn (~r) es un funcional
de la densidad en sí mismo. Al igual que sucede en el cálculo de funciones de una variable, si nos
limitamos a variaciones lineales Z
3 δg [n]
δg = d ~r δn (~r)
δn (~r)
3. A partir de la de nición (5.1) es inmediato probar que si g [n] y h [n] son dos funcionales de la
UNED
densidad y α, β dos escalares,
δ δg [n] δh [n]
(α g [n] + β h [n]) = α +β
δn (~r) δn (~r) δn (~r)
esto es, la derivada variacional es lineal. A su vez, si G (u) es una función derivable de una variable y
h [n] un funcional de la densidad
δ δh [n]
G (h [n]) = G 0 (h [n])
δn (~r) δn (~r)
donde G 0 (u) es la derivada de la función G.
1
Naturalmente la de nición se puede extender a cualquier otro conjunto de funciones.
127
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
Un caso que merece la pena analizar es el siguiente. Consideremos un funcional g [n] de la forma
Z
g [n] = ~ (~r)
d3~r F ~r, n (~r) , ∇n
UNED
donde F es una función de siete variables (las tres componentes de la posición, la densidad, y las tres
componentes del gradiente de la densidad). Bajo una variación δn (~r), aplicando análisis elemental al
integrando
Z ~ (~r)
∂ F ~r, n (~r) , ∇n 3 ∂F ~ ~ (~r)
r, n (~r) , ∇n
X ∂
δg = d3~r
δn (~
r ) + δn (~r)
UNED
∂n (~r) ∂n (~r) ∂xα
α=1 ∂
∂xα
Si ahora integramos por partes los tres últimos sumandos y aplicamos la de nición (5.1), tenemos que
~ (~r)
∂ F ~r, n (~r) , ∇n 3 ~ (~r)
∂ F ~r, n (~r) , ∇n
δg [n] X ∂
= − (5.2)
δn (~r) ∂n (~r) ∂xα ∂n (~r)
UNED
α=1 ∂
∂xα
resultado que usaremos repetidas ocasiones. En particular, a partir de (5.2) es fácil ver que si G (u) y
H (u) son funciones de una variable y f (~r) una función de la posición
Z
δ
d3~r 0 f ~r 0 G n ~r 0 = f (~r) G 0 (n (~r))
δn (~r)
UNED
Z
δ
d3~r 0 G n ~r 0 H n ~r 0 = G 0 (n (~r)) H (n (~r)) + G (n (~r)) H 0 (n (~r))
δn (~r)
4. Como ilustración de los resultados matemáticos que acabamos de estudiar, veamos unos ejemplos:
UNED
Ejemplo 5.a. Si G (u) es una función derivable, pruebe que
δG (n (~r 0 ))
= δ ~r − ~r 0 G 0 (n (~r))
δn (~r)
Solución:
Es inmediata. Si usamos la delta de Dirac tenemos que
UNED
Z
0
d3~r 00 δ ~r 00 − ~r 0 G n ~r 00
G n ~r =
y entonces
δG (n (~r 0 ))
Z
δ
d3~r 00 δ ~r 00 − ~r 0 G n ~r 00 = δ ~r − ~r 0 G 0 (n (~r))
=
δn (~r) δn (~r)
donde hemos usado Z
δ
d3~r 00 f ~r 00 G n ~r 00 = f (~r) G 0 (n (~r))
δn (~r)
128
UNED
El teorema de Hohenberg-Kohn
Máster FSC. Funcional de la densidad
Ejemplo 5.b. A partir de la de nición (5.1), obtenga la primera y segunda derivadas
variaciones de los funcionales Vext [n] y WH [n].
UNED
Solución:
Basta expandir los funcionales actuando sobre n + δn. Para la energía potencial
Z Z
3
Vext [n + δn] = d ~r vext (~r) [n (~r) + δn (~r)] = Vext [n] + d3~r vext (~r) δn (~r)
y, por tanto,
UNED
δVext [n] δ 2 Vext [n]
= vext (~r) ; =0
δn (~r) δn (~r) δn (~r 0 )
A su vez, para la energía de Hartree e intercambiando índices si es preciso
[n (~r) + δn (~r)] [n (~r 0 ) + δn (~r 0 )]
Z
1
WH [n + δn] = d3~r d3~r 0 =
2 |~r − ~r 0 |
n (~r 0 ) r) δn (~r 0 )
Z Z Z
1 3 3 0 δn (~
= WH [n] + d3~r d3~r 0
UNED
δn (~
r ) + d ~
r d ~
r
|~r − ~r 0 | 2 |~r − ~r 0 |
y entonces
n (~r 0 ) δ 2 WH [n]
Z
δWH [n] 1
d3~r 0 0
= = vH (~r) ; = = w ~
r , ~
r
δn (~r) |~r − ~r 0 | δn (~r) δn (~r 0 ) |~r − ~r 0 |
UNED
Ejemplo 5.c. Demuestre que el llamado funcional de Weizsäcker
2
Z ~ n (~r)
∇ Z 3 2
1 3 1 3 1 X ∂n (~r)
TW [n] = d ~r = d ~r
8 n (~r) 8 n (~r) ∂xα
α=1
UNED
evaluado sobre la densidad del estado fundamental n0 (~r) de un sistema no interactuante
de N = 2 partículas sometidas a un potencial exterior independiente del espín, es igual
a la energía cinética exacta de dicho estado. Usando (5.2), obtenga la primera derivada
variacional de TW [n].
Solución:
Para el sistema indicado, las dos partículas están en el mismo orbital espacial φ0 (~r) con espines
UNED
opuestos. Sin pérdida de generalidad podemos imponer que φ0 (~r) es real, por lo que la densidad
fundamental está dada por n0 (~r) = 2 × φ0 (~r)2 . Entonces:
2 2
Z ~ n (~r)
∇ Z ~ φ0 (~r)2
∇
31 1 3
TW [n] = d ~r = d ~r =
8 n (~r) 4 φ0 (~r)2
Z ~ 0 (~r) · φ0 (~r) ∇φ
φ0 (~r) ∇φ ~ 0 (~r) Z n o
= 3
d ~r ~ 0 (~r) · ∇φ
= d3~r ∇φ ~ 0 (~r)
φ0 (~r) φ0 (~r)
Si ahora integramos por partes teniendo en cuenta que las funciones de onda se anulan en el in nito:
Z
1 3 2
TW [n] = 2 × − d ~r φ0 (~r) ∇ φ0 (~r)
2
129
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
que es la suma de las energías cinéticas de las dos partículas.
La derivada variacional se puede efectuar usando las reglas generales que hemos visto en la teoría,
aunque también se puede obtener directamente expandiendo TW [n + δn]. Siguiendo este segundo
UNED
método
2
Z ~ n (~r) + ∇
∇ ~ δn (~r)
1
TW [n + δn] = d3~r =
8 n (~r) + δn (~r)
Z
1 3 1 ~ ~
2 δn (~r)
= d ~r ∇ n (~r) + ∇ δn (~r) 1− =
8 n (~r) n (~r)
UNED
Z
1 3 1 ~
2
~ ~ δn (~r)
= d ~r ∇n (~r) + 2∇n (~r) · ∇δn (~r) 1−
8 n (~r) n (~r)
donde hemos despreciado todos los términos de orden δn (~r)2 y superiores. Si ahora expandemos el
producto
Z Z
1 3 1 ~ (~r) δn (~r) + 1
2 1 ~ ~ (~r)
TW [n + δn] = TW [n] − d ~r 2 ∇n d3~r ∇n (~r) · ∇δn
8 n (~r) 4 n (~r)
UNED
e integramos por partes el último término
Z Z
1 3 1 ~ (~r) δn (~r) − 1
2
3~ · 1 ~
TW [n + δn] = TW [n] − d ~r 2 ∇n d ~r ∇ ∇n (~r) δn (~r)
8 n (~r) 4 n (~r)
De este modo
2
UNED
~ ~ 2 n (~r)
1 ∇n (~r)
δTW [n] 1 1 ~ (~r) − 1 ∇
~ ·
2 1 ~ 1∇
=− 2 ∇n ∇n (~r) = −
δn (~r) 8 n (~r) 4 n (~r) 8 n2 (~r) 4 n (~r)
5.2. EL POTENCIAL QUÍMICO
1. Consideremos un sistema electrónico en el cual el potencial exterior vext (~r) es local e independiente
del espín. Esto es, el operador Hamiltoniano es
UNED
Z
H
b = Tb + W
c+ d3~r vext (~r) n
b (~r)
con n
b (~r) el operador densidad en el punto ~r. Impongamos que dicho hamiltoniano es convexo, es decir,
(N ) 2
que si E0 es la energía fundamental del sistema con N electrones se cumple que
(N ) (N +1) (N −1) (N )
E0 − E0 ≤ E0 − E0 ; N = 1, 2, 3, ...
UNED
Nótese que si de nimos las energías de ionización y la a nidad electrónica del sistema de N electrones
en su estado fundamental como
(N −1) (N ) (N ) (N +1)
I (N ) = E0 − E0 ; A(N ) = E0 − E0 (5.3)
la condición de convexidad es equivalente a I (N ) ≥ A(N ) ; N = 1, 2, 3, ...
De acuerdo con la de nición, si el sistema de N electrones se encuentra en su estado fundamental,
I (N ) es la energía que hay que comunicar al sistema para arrancar uno de ellos de manera que el sistema
resultante de N − 1 electrones quede en su estado fundamental. Análogamente, A
(N ) es la energía que
2
Puede parecer una restricción fuerte, pero todos los operadores energía electrónicos de interés en tres dimensiones son convexos.
130
UNED
El teorema de Hohenberg-Kohn
Máster FSC. Funcional de la densidad
cedería el sistema al añadir un electrón si el sistema resultante de N + 1 electrones quedase en su
estado fundamental. A su vez, como consecuencia de la de nición I (N ) = A(N −1) y A(N ) = I (N +1) .
UNED
2. Esta condición de convexidad es importante si consideramos situaciones en las que el número de
partículas pueda variar (por ejemplo si el sistema está en contacto con un baño térmico con el cual
puede intercambiar partículas). Ya que H
b conserva el número de partículas, para un estado genérico
del espacio de Fock
∞
X E ∞
X E
|Ψi = cJ Ψ (J)
; |cJ |2 = 1 ; Ψ(J) ∈ HJ
J=0 J=0
UNED
los valores esperados de la energía y del número de partículas son
D E X∞ D E D E X∞
Ψ H Ψ =
b |cJ |2 Ψ(J) H
b Ψ(J) ; Ψ N Ψ =
b |cJ |2 J
J=0 J=0
Si ahora nos preguntamos cuál es el estado con N electrones (en el sentido general hN
bi = N) con
UNED
menor energía, ha de cumplirse que este estado sea
∞ E ∞ ∞
(J)
X X X
|Ψi = cJ Ψ0 con |cJ |2 J = N y |cJ |2 = 1
J=0 J=0 J=0
E
(J)
donde Ψ0 es el estado fundamental del sistema con un número J de electrones bien de nido, es
decir, el estado fundamental de H
b restringido al espacio HJ . Los coe cientes cJ , con las restricciones
UNED
indicadas, han de minimizar la cantidad
∞ D E
(J) (J) (J) b (J)
X
|cJ |2 E0 ; E0 = Ψ0 H Ψ0 .
J=0
La condición de convexidad garantiza que siN ∈ N, |cN | = 1 y cJ = 0 para J 6= N . Esto es, de
entre todos los estados con valor esperado del número de electrones igual a
E N (N entero), el de menor
(N )
energía es Ψ0 : el estado fundamental de H restringido a estados con un número de partículas bien
UNED
b
de nido e igual a N.
A su vez, consideremos el caso hN
b i = N + ∆, con ∆ ∈ (0, 1). La condición de convexidad garantiza
que, ahora, el estado normalizado de menor energía ha de tener la forma
√ (N )
E √ (N +1)
E
1 − ∆ Ψ0 + ∆eiφ Ψ0
con φ una fase arbitraria. Ello implica que el valor mínimo del valor esperado de la energía es
UNED
(N ) (N +1)
(1 − ∆) E0 + ∆E0
por lo que la función E0 (M ), igual a la energía mínima de un sistema electrónico con valor esperado
del número de partículas igual a M, es una función lineal a trozos, tal como se ilustra en la gura 5.1.
3. Obsérvese que no estamos haciendo otra cosa que estudiar el sistema electrónico desde el punto de
vista de la colectividad macrocanónica en el límite de temperatura nula (aunque sacri cando el rigor
que exigiría el describir los estados mediante mezclas estadísticas). Por tanto, el potencial químico a
temperatura nula para un sistema con un valor esperado M del número de partículas es
∂E0 (M )
µ (M ) = (5.4)
∂M
131
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
1
_Ψ !
UNED
1 1
1 α _Ψ !β _Ψ !
( 0
,
1
_Ψ !
1
$
UNED
1 1 1
Figura 5.1. Energia fundamental E0 (M ) para un sistema de electrones cuyo valor esperado del número de partículas es
M (no necesariamente entero).
UNED
A la vista del resultado anterior, el potencial químico puede sufrir una discontinuidad cuando M es
igual a un número entero N. En concreto, en torno a un valor entero N se cumple
−I (N )
si M = N − ∆
µ (M ) = (5.5)
−A(N ) = −I (N +1) si M = N + ∆
estando inde nido para M = N, aunque por convenio se suele tomar el valor medio en la
UNED
discontinuidad:
I (N ) + A(N ) I (N ) + I (N +1)
µ (N ) = − =− (5.6)
2 2
5.3. EL TEOREMA DE HOHENBERG Y KOHN
5.3.a. El teorema básico de la Teoría del Funcional de la Densidad
UNED
1. El teorema de Hohenberg y Kohn [Phys. Rev. 136, B864 (1964)] establece que la densidad es
la única variable estrictamente necesaria para conocer el estado fundamental de un sistema de N
electrones sometidos a la acción de un potencial exterior local vext (~r). Expongamos dicho teorema
siguiendo una formulación más moderna, debida sobre todo a Levy y Lieb.
Consideremos el conjunto formado por todos los estados cuánticos electrónicos |Ψi del espacio de
Fock F. De nimos el funcional de Hohenberg y Kohn como
UNED
D E
FHK [n] = mı́n Ψ Tb + W
c Ψ (5.7)
|Ψi→n
es decir, FHK [n] es el mínimo del valor esperado de la suma de la energía cinética y de interacción sobre
todos los estados cuánticos cuya densidad es n (~r). Nótese como FHK [n] es un funcional universal de
la densidad, en el sentido de que su de nición no depende de ningún sistema físico (esto es, de ningún
potencial exterior en concreto).
Esta de nición de FHK [n] no está restringida a densidades que normalizen a un número entero de
partículas. En este caso, los estados cuánticos|Ψi cuya densidad sea igual a n (~r) son, necesariamente,
superposición de estados |ΨJ i con diferente número de partículas. Incluso en el caso de que n (~r) fuese
una densidad que normalize a un número entero de electrones, bien pudiera ocurrir que el estado |Ψi
D E
con densidad n (~
r) que minimize Ψ Tb + W c Ψ no sea autoestado de N b.
132
UNED
El teorema de Hohenberg-Kohn
Máster FSC. Funcional de la densidad
Consideramos ahora un sistema eléctrónico de N electrones cuyo Hamiltoniano es
Z
H
b = Tb + W
c+ d3~r vext (~r) n
b (~r)
UNED
al cual únicamente impondremos la condición de convexidad. Si de nimos el funcional de la densidad
Z
Evext [n] = FHK [n] + d3~r n (~r) vext (~r) , (5.8)
se cumple el Teorema de Hohenberg y Kohn:
UNED
a) Si n0 (~r) es la densidad del estado fundamental
3
se cumple que
Z
E0 = Evext [n0 ] = FHK [n0 ] + d3~r n0 (~r) vext (~r) (5.9)
b) E0 = Evext [n0 ] es el mínimo del funcional Evext [n] sobre todas las densidades físicas que
normalizan aN electrones:
UNED
Z
E0 = mı́n Evext [n] = mı́n FHK [n] + d3~r n (~r) vext (~r) (5.10)
n(~
r)→N n(~
r)→N
La demostración es relativamente sencilla. Observemos que
D E Z D E
Evext [n0 ] = mı́n Ψ T + W Ψ + d3~r n0 (~r) vext (~r) = mı́n Ψ H
b c b Ψ
|Ψi→n0 |Ψi→n0
UNED
Evext [n0 ] es el mínimo del valor esperado de la energía sobre todas las funciones de onda cuya
esto es,
r). Entre estas funciones de onda está el propio estado fundamental |Ψ0 i por lo que
densidad es n0 (~
de acuerdo con el principio variacional, dicho mínimo ha de alcanzarse para |Ψi = |Ψ0 i. Por tanto
D E
Evext [n0 ] = Ψ0 Hb Ψ0 = E0
quedando demostrada la igualdad. Vemos que el carácter convexo del hamiltoniano garantiza que
UNED
ésta se cumpla incluso aunque contemplemos en la de nición de FHK [n] estados |Ψi sin número de
partículas bien de nido. Consecuencia inmediata es el hecho de que
D E
FHK [n0 ] = Ψ0 Tb + W
c Ψ
b0 (5.11)
Por tanto, si n0 (~r) es densidad fundamental de un sistema N -électrónico sometido a un potencial
exterior local, FHK [n0 ] es igual al valor esperado de T
b+W
c sobre el estado fundamental correspondiente.
UNED
Ahora bien, dicho funcional se construye mediante la búsqueda del mínimo de T b+Wc restringido al
conjunto de estados cuya densidad es n0 (~ r). En de nitiva, vía el funcional de Hohenberg y Kohn, la
función de ondas del estado fundamental es funcional de la densidad y, por ende, todas las propiedades
del estado fundamental, incluida la energía, son funcionales de la densidad. A la vista de todo esto,
la condición de mínimo es inmediata. Si existiese una densidad
E ñ (~r) diferente a n0D(~r) paraEla cual
Evext [ñ] ≤ E0 existiría un estado cuántico Ψ
e , cuya densidad es ñ (~r), para el que Ψe H e ≤ E0 ,
b Ψ
lo que está prohibido por el principio variacional.
2. Todo lo anterior es válido incluso si el nivel fundamental es degenerado. En este caso el
funcional Evext [n] alcanzaría su mínimo absoluto sobre cualquier miembro de la familia de densidades
3
O una densidad del nivel fundamental si éste es degenerado.
133
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
correspondientes a cualquier estado del nivel fundamental. Además podríamos considerar densidades
correspondientes a estados mezcla sin más que extender la de nición de FHK [n] a este tipo de estados
(lo que es pertinente si se desea tratar sistemas electrónicos a temperatura distinta de cero). En el
UNED
muy hipotético caso de que el hamiltoniano no fuese convexo basta restringir la de nición de FHK [n]
a una minimización sobre estados con un número N de electrones bien de nido aunque, como veremos
a continuación, perderíamos la conexión entre el teorema de Hohenberg y Kohn y el potencial químico
del sistema.
Nótese también que, si como estamos suponiendo, el potencial exterior no depende del espín, la
variable fundamental es la densidad total, no la densidad dependiente del espín. En este contexto,
UNED
los valores esperados relativos al espín son propiedades del estado fundamental. Por el contrario, si el
potencial exterior sí dependiese del espín, la variable fundamental sería la densidad dependiente del
espín n (~r, σ) ≡ nσ (~r). En esta situación basta cambiar la variable ~r por la variable ~ = (~r, σ)
κ en
toda la teoría anterior, teniéndose así de manera natural una extensión de la DFT a sistemas bajo
potenciales dependientes del espín.
Sin embargo, la teoría tal y como la hemos presentado no es aplicable al caso de un sistema sometido
a un campo electromagnético estático ya que entonces el término de interacción con el campo exterior
UNED
depende, además, de la densidad de corriente. La DFT también puede extenderse también a este tipo
de situaciones (current-density functional theory), pero el tratamiento es demasiado avanzado y no lo
trataremos en estas notas.
3. Puesto que Evext [n] alcanza un mínimo para n0 (~r):
Z
3
UNED
δ Evext [n] − µ d ~r n (~r) =0
n=n0
donde µ es un multiplicador de Lagrange cuyo signi cado físico veremos a continuación. Admitiendo
la derivabilidad variacional del funcional FHK [n] en un entorno de n0 (~r),4 se tiene que la densidad
fundamental satisface la siguiente ecuación tipo Euler-Lagrange
δEvext [n] δFHK [n]
UNED
= + vext (~r) = µ (5.12)
δn (~r) n=n0 δn (~r) n=n0
junto con la condición de normalización.
Z
d3~r n0 (~r) = N
UNED
Esto es, el problema 3N -dimensional (excluidas las variables de espín) queda formalmente reducido a
un simple problema tridimensional. Empero, el teorema de Hohenberg-Kohn expresa la existencia del
funcional universal FHK [n], pero nada dice sobre su forma y propiedades. La aplicación del teorema
básico de la Teoría del Funcional de la Densidad exige aproximar dicho funcional.
4. Tal y como hemos de nido el funcional FHK [n], el teorema de Hohenberg-Kohn y la ecuación
de Euler-Lagrange correspondiente se pueden extender a situaciones en las que el valor esperado del
(M )
número de electrones no sea entero. Sea entonces n0 (~r) la densidad fundamental cuando hNbi = M
(M +∆) (M )
y n0 (~r) = n0 (~r) + δn (~r) la densidad fundamental cuando hN i = M + ∆, con ∆ muy pequeño.
b
4
Este es un problema matemático sutil que ha dado lugar a una muy amplia literatura. En general la derivada variacional sólo
está de nida para densidades v -representables, de nidas como aquellas que son densidad fundamental de algún hamiltoniano.
134
UNED
El teorema de Hohenberg-Kohn
Máster FSC. Funcional de la densidad
Así
h i
(M )
E0 (M ) = Evext n0
UNED
h i h i Z δEvext [n]
(M ) (M )
E0 (M + ∆) = Evext n0 + δn = Evext n0 + d3~r δn (~r)
δn (~r) n=n0
(M )
y teniendo en cuenta la ecuación de Euler-Lagrange
UNED
E0 (M + ∆) = E0 (M ) + µ (M ) d3~r δn (~r)
Ahora bien, δn (~r) normaliza a ∆, ya que es la diferencia entre una densidad con M +∆ partículas y
otra con M partículas. En de nitiva
E0 (M + ∆) − E0 (M ) ∂E0 (M ) δEvext [n] δFHK [n]
µ (M ) = = = = + vext (~r)
∆ ∂M δn (~r) n=n0
(M ) δn (~r) n=n0
(M )
UNED
(5.13)
y, por tanto, el multiplicador de Lagrange que aparece es, sencillamente, el potencial químico del
sistema.
Puesto que si I (N ) 6= A(N ) , µ (M ) tiene una discontinuidad en el número entero N de partículas, la
(N )
derivada variacional de FHK [n] evaluada en la densidad n0 (~r) está de nida salvo constante aditiva
(lo que a efectos prácticos de resolución de la ecuación de Euler-Lagrange no es relevante). Esta
UNED
discontinuidad puede expresarse como
(N −∆)
δFHK [n] −I (N ) − vext (~r)
si n (~r) = n0 (~r)
= ; ∆ = 0+ (5.14)
δn (~r) −A(N ) − v (~r) si n (~r) = n(N +∆) (~r)
ext 0
5.3.b. Relación con la teoría de la respuesta lineal
UNED
5. Consideremos, de nuevo, un Hamiltoniano de N electrones bajo la acción de un potencial exterior
vext (~r). Sean |Ψ0 i, n0 (~r) y E0 el estado, la densidad y la energía fundamentales, respectivamente.
Supongamos que sobre el sistema actúa una perturbación estática δv (~r). De acuerdo con la teoría
general de perturbaciones estáticas, la variación en primer orden de la densidad y de segundo orden
en la energía están dadas por
UNED
Z
δn (~r1 ) = d3~r2 χ (~r1 , ~r2 ) δv (~r2 )
Z Z
1
δE0 = d3~r1 n0 (~r1 ) δv (~r1 ) + d3~r1 d3~r2 χ (~r1 , ~r2 ) δv (~r1 ) δv (~r2 )
2
donde
∞
X hΨ0 |b
n (~r1 )| Ψi i hΨi |b
n (~r2 )| Ψ0 i + hΨ0 |b
n (~r2 )| Ψi i hΨi |b
n (~r1 )| Ψ0 i
χ (~r1 , ~r2 ) = − = χ (~r2 , ~r1 )
Ei − E0
i=1
es la función respuesta lineal estática χ (~r1 , ~r2 ) = lı́mω→0 χ (~r1 , ~r2 ; ω) escrita en términos de las energías
Ei de los estados excitados |Ψi i del Hamiltoniano sin perturbar.
135
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
Puesto que E0 es un funcional de la densidad, también podemos escribir
δ 2 Evext [n]
Z Z
3 δEvext [n] 1
δE0 = d ~r δn (~r) + d3~r1 d3~r2 δn (~r1 ) δn (~r2 ) +
UNED
δn (~r) n=n0 2 δn (~r1 ) δn (~r2 ) n=n0
Z
+ d3~r1 [n0 (~r1 ) + δn (~r2 )] δv (~r1 )
donde los dos primeros sumandos son el cambio en el valor esperado del Hamiltoniano original, mientras
que el tercero es el valor medio de la perturbación sobre el nuevo estado. El primer sumando es nulo,
ya que Evext [n] es estacionario en torno a n0 (~r). Teniendo además en cuenta que
UNED
δ 2 Vext [n]
=0
δn (~r1 ) δn (~r2 )
que resulta del hecho de que Vext [n] es un funcional lineal de la densidad, la variación de la energía
fundamental está dada por:
δ 2 FHK [n]
Z Z
1 3 3
δE0 = d ~r1 d ~r2 δn (~r1 ) δn (~r2 ) + d3~r1 [n0 (~r1 ) + δn (~r1 )] δv (~r1 )
UNED
2 δn (~r1 ) δn (~r2 ) n=n0
Si igualamos las dos expresiones para δE0 tenemos que
δ 2 FHK [n]
Z Z
1 3 3 1
d ~r1 d ~r2 χ (~r1 , ~r2 ) δv (~r1 ) δv (~r2 ) = d3~r1 d3~r2 δn (~r1 ) δn (~r2 ) +
2 2 δn (~r1 ) δn (~r2 ) n=n0
Z
+ d3~r1 δn (~r1 ) δv (~r1 )
UNED
y escribiendo δn (~r1 ) en términos de δv y de la función respuesta lineal, llegamos a la igualdad (en
sentido matricial) K
b HK χb = −bI, donde
δ 2 FHK [n]
KHK (~r1 , ~r2 ) =
δn (~r1 ) δn (~r2 ) n=n0
UNED
En de nitiva
. δ 2 FHK [n]
KHK (~r1 , ~r2 ) = = −χ−1 (~r1 , ~r2 ) (5.15)
δn (~r1 ) δn (~r2 ) n=n0
importante expresión que relaciona la segunda derivada variacional del funcional de Hohenberg y Kohn
para la densidad de un estado fundamental con la función respuesta lineal estática de dicho estado.
5.4. EL MODELO DE THOMAS-FERMI
UNED
1. Recordemos el resultado central del teorema de Hohenberg-Kohn. Si un sistema de N electrones
está sometido a la acción de un potencial exterior vext (~r) independiente del espín, la energía
fundamental de dicho sistema es el mínimo del funcional
Z
E [n] = FHK [n] + d3~r vext (~r) n (~r)
5
sobre el conjunto de todas las densidades n (~r) asociadas a un estado cuántico de N electrones. Siendo
el resultado exacto, surgen dos di cultades obvias: determinar la forma del funcional FHK [n] y saber,
a priori, que densidades n (~r) son realmente densidades de un estado cuántico de N electrones. La
5
Por lo general, para descargar la notación se suele escribir E [n] en lugar de Evext [n]
136
UNED
El teorema de Hohenberg-Kohn
Máster FSC. Funcional de la densidad
respuesta a la segunda pregunta es relativamente sencilla: prácticamente todas las funciones continuas
no negativas que normalicen a N y para las que el funcional de Weizsäcker sea nito:
UNED
2
Z Z ~ n (~r)
∇
1
n (~r) ≥ 0 ; d3~r n (~r) = N ; TW [n] = d3~r <∞,
8 n (~r)
son densidad de un estado de N
electrones (fermiones). La respuesta a la primera no lo es tanto, y es
6
necesario formular aproximaciones al funcional universal de Hohenberg-Kohn.
2. Históricamente, aunque en un contexto diferente, el modelo estadístico de Thomas-Fermi puede
UNED
considerarse como la implementación más sencilla de la DFT. Sin discutir la derivación original de
dicha aproximación, en el lenguaje de la DFT es equivalente a suponer que
FHK [n] ' TTF [n] + WH [n] (5.16)
donde TTF [n] es el llamado funcional cinético de Thomas-Fermi dado por
2/3 Z
3 3π 2
Z
UNED
TTF [n] = d3~r n (~r) thom
s (n (~r)) = d3~r n (~r)5/3 (5.17)
10
siendo ts
hom
(n) la energía cinética de un gas de electrones libres de densidad n, mientras que WH [n] es
la energía electrostática clásica de Hartree. Observemos que en esta aproximación funcional, si n0 (~ r)
es la densidad fundamental de nuestro sistema electrónico, el valor esperado de la energía cinética está
dado por TTF [n0 ] y el de la energía de interacción electrónica por WH [n]. De esta manera, el modelo
Thomas-Fermi desprecia cualquier efecto cuántico en la energía de interacción mientras que supone
UNED
que la densidad de energía cinética en un punto ~r (es decir, la energía cinética por unidad de volumen
en dicho punto) es igual a la correspondiente a un gas de electrones libres cuya densidad es n (~r).
Ambas aproximaciones son muy crudas, pero bajo éstas la implementación de la ecuación de
Euler-Lagrange (5.12) es inmediata. Basta ver que
2/3
3π 2 n (~r 0 ) .
Z
δTTF [n] δWH [n]
= n (~r)2/3 ; = d3~r 0 = vH [n] (~r)
δn (~r) 2 δn (~r) |~r − ~r 0 |
UNED
y, entonces, la densidad fundamental veri ca la ecuación integral
2/3
3π 2 n (~r 0 )
Z
2/3
n (~r) + d3~r 0 + vext (~r) = µ (5.18)
2 |~r − ~r 0 |
sujeta a la condición de normalización
UNED
Z
d3~r n (~r) = N
La ecuación (5.17) puede resolverse numéricamente sin ningún problema (y en algunos casos, como los
átomos, admite una solución analítica). Naturalmente, el método de Thomas-Fermi no es preciso en
absoluto, pero ofrece buenas estimaciones de los órdenes de magnitud de las energías fundamentales.
3. Obsérvese que el modelo de Thomas-Fermi desprecia los efectos de intercambio-correlación en
la energía de interacción. Por ello puede considerarse como una aproximación funcional al método
de Hartree en el que la energía cinética es evaluada no en términos de las energías cinéticas de los
6
Hay un aspecto sútil que no hemos comentado. El teorema de Hohenberg-Kohn sólo es estrictamente válido para sistemas
físicos con un número nito de electrones. Sin embargo, los problemas formales que se plantearían en un sistema extenso (que es
una idealización) no son problemáticos en absoluto.
137
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
orbitales sino mediante el promedio estadístico (5.17). De esta manera tenemos dos fuentes de error:
las inherentes a la aproximación de Hartree y la debida a la aproximación funcional TTF [n].
UNED
Teniendo esto en cuenta, Dirac propuso mejorar el modelo Thomas-Fermi incluyendo un término
que contemplase los efectos de intercambio y construido de igual manera que el funcional de
Thomas-Fermi TTF [n]. Si recordamos los resultados para el gas homogéno, Dirac sugirió que la energía
de intercambio para un sistema electrónico arbitrario de densidad n (~r) puede aproximarse mediante
la siguiente expresión funcional
Z 1/3 Z
3 3
EX ' Dir [n]
EX = d3~r n (~r) εhom (n (~r)) =− d3~r n (~r)4/3 (5.19)
UNED
X
4 π
donde εhom
X (n) es la energía de intercambio por partícula para un gas homogéneo de electrones de
densidad n. Así
Dir
FHK [n] ' TTF [n] + WH [n] + EX [n] (5.20)
es una aproximación al funcional de Hohenberg y Kohn que puede verse como una simpli cación del
UNED
método Hartree-Fock y que es conocida en la literatura como el funcional de Thomas-Fermi-Dirac
(TFD). La ecuación de Euler-Lagrange correspondiente es
2/3 1/3
3π 2 r 0)
Z Z
2/3 3 0 n (~ 3 1/3
n (~r) + d ~r − n (~r) + vext (~r) = µ con d3~r n (~r) = N (5.21)
2 |~r − ~r 0 | π
que, de nuevo, no plantea problema numérico especial alguno.
UNED
4. Incluso con la adición del funcional de Dirac, esta aproximación funcional sigue sin ser
precisa. Sin embargo puede considerarse el punto de partida para un planteamiento del problema
electrónico bastante atractivo: buscar buenas aproximaciones al funcional de Hohenberg-Kohn con una
dependencia explícita en la densidad electrónica. Si esto fuese posible se reduciría de manera efectiva
(no simplemente formal) el problema de 3N variables asociado a la función de ondas multielectrónica
(ecuación de Schrödinger) a un problema de 3 variables asociado a la densidad electrónica n (~r). Esto
UNED
conduce a modelos semiclásicos o hidrodinámicos que se utilizan con bastante éxito (si tenemos en
cuenta el coste computacional asociado) en el estudio de uidos complejos en el que las interacciones
efectivas entre los iones no se modelizan, sino que se obtienen en base a cálculos de primeros principios
usando este tipo de funcionales. Ahora bien como veremos en el próximo tema, la forma habitual de
implementar la DFT no es esta, sino una bastante diferente debida a Kohn y Sham. En todo caso, y
como ejemplo de la manera de proceder en estos modelos hidrodinámicos, veamos el siguiente ejemplo
UNED
Ejemplo 5.d. Considérese el siguiente modelo funcional para FHK [n]
Dir
FHK [n] = TW [n] + TTF [n] + WH [n] + EX [n]
donde TW [n] es el funcional de Weizsäcker introducido en el ejemplo
5.c. Esta aproximación (conocida, por motivos, obvios, como funcional de
Thomas-Fermi-Dirac-Weizsäcker) introduce en el funcional de Hohenberg-Kohn un
término característico de sistemas con pocos electrones aislados (TW [n]) que intenta
complementar el funcional original TFD, más adecuado para sistemas extensos.
a) Usando el resultado del ejemplo 5.c, plantéese la ecuación de Euler-Lagrange
correspondiente para un sistema de N electrones sometidos a la acción de un potencial
vext (~r).
138
UNED
El teorema de Hohenberg-Kohn
Máster FSC. Funcional de la densidad
b)
p
Si de nimos el orbital normalizado Ψ (~
r) = n0 (~r) /N , donde n0 (~r) es la densidad
fundamental del sistema, demuéstrese que Ψ (~r) y µ veri can la ecuación tipo Schrödinger
UNED
2/3
2 !1/3
1~2
Z
N Ψ (~r 0 )2 3π 2 N Ψ (~r) 3N Ψ (~r)2
− ∇ + vext (~r) + d3~r 0 + − Ψ (~r) = µΨ (~r)
2 |~r − ~r 0 | 2 π
que puede resolverse autoconsistentemente. Así, bajo este tipo de modelos funcionales el
problema de N electrones queda reducido al problema trivial de una única partícula efectiva
cuya función de ondas es Ψ (~r).
UNED
Solución:
a) La ecuación de Euler-Lagrange es
δFHK [n]
+ vext (~r) = µ
δn (~r)
UNED
y usando las derivadas variacionales de las distintas componentes de FHK [n], ésta es igual a
2
~ 2 2/3
1 ∇n (~r) ~ n (~r) 1/3
3π 2 n (~r) n (~r 0 )
1∇
Z
3 0 3n (~r)
− + + d ~r − + vext (~r) = µ
8 n2 (~r) 4 n (~r) 2 |~r − ~r 0 | π
b) Si ahora sustituimos n (~r) por N Ψ (~r)2 , tras efectuar operaciones elementales tenemos que la
ecuación de Euler-Lagrange es
UNED
2/3
2 N Ψ (~ 2 !1/3
1
1~2
3π r ) Z
N Ψ (~r)2 3N Ψ (~r)2
− ∇ Ψ (~r) + + d3~r 0 − + vext (~r) = µ
Ψ (~r) 2 2 |~r − ~r 0 | π
y multiplicando los dos miembros por Ψ (~r) se tiene la ecuación tipo Schrödinger indicada en el
enunciado.
UNED
UNED 139
Máster FSC. Funcional de la densidad
UNED
UNED
UNED
UNED
UNED
UNED
UNED
UNED
Máster FSC. Funcional de la densidad
Capítulo 6
MÉTODO DE KOHN-SHAM I
UNED
OBJETIVOS:
1. De nición de los funcionales de la energía cinética, intercambio y correlación.
UNED
2. Tratamiento de la implementación Kohn-Sham de la DFT
3. Análisis de algunas propiedades exactas del funcional de la energía XC.
4. Discusión de las aproximaciones LDA y GGA.
UNED
LECTURAS COMPLEMENTARIAS:
1. A primer in Density Functional Theory [Fiolhais et al ( Eds.)]: Capítulo 1
2. Fiolhais et al ( Eds.): Capítulo 2
ORIENTACIONES AL PLAN DE TRABAJO:
UNED
En este capítulo nos centraremos en la implementación más popular de la Teoría del Funcional
de la Densidad, debida a Kohn y Sham. La idea básica consiste en de nir un sistema cticio
de electrones no interactuantes cuya densidad es igual a la del sistema real. De acuerdo con el
teorema fundamental de la DFT, el conocimiento del sistema cticio permite (a través de la
densidad del estado fundamental) obtener las propiedades del estado fundamental del problema
UNED
electrónico.
A título orientativo, se espera que dedique 10 horas a este capítulo desglosadas como sigue:
• 5 horas de lectura detallada y estudio del presente material didáctico.
• 5 horas para el a anzamiento de los contenidos, incluyendo consultas y las lecturas
complementarias [optativas].
UNED
En este capítulo no se proponen problemas.
GLOSARIO:
Funcional de la energía cinética de Kohn-Sham
Funcional de la energía de intercambio-correlación
Ecuaciones de Kohn-Sham
Potencial de intercambio-correlación
Sistema cticio Kohn-Sham
Conexión adiabática
141
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
Función de correlación de pares [integrada]
Hueco de intercambio-correlación
UNED
Teorema ACFDT
Teorema de compresibilidad
Aproximación local (LDA)
Aproximación local dependiente del espín (LSDA)
UNED
Aproximación de gradientes generalizada (GGA)
Funcional PBE (Perdew-Burke-Erzenhof )
UNED
UNED
UNED
UNED
142
UNED
Método de Kohn-Sham I
Máster FSC. Funcional de la densidad
6.1. EL FORMALISMO DE KOHN Y SHAM
6.1.a. Los funcionales cinético y de intercambio-correlación
UNED
1. En el capítulo anterior hemos visto que para un sistema de N electrones sometido a la acción de
un potencial exterior longitudinal vext (~r), la energía fundamental del sistema está dada por
Z
3
E0 = mı́n E [n] = mı́n FHK [n] + d ~r vext (~r) n (~r) ,
n(~
r)→N n(~
r)→N
es decir, por el mínimo del funcional de la densidad E [n] sobre todas las posibles densidades
UNED
electrónicas físicas de N electrones. Dicho mínimo se alcanza sobre la densidad del estado
fundamental n0 (~r) (o sobre una densidad fundamental si el nivel fundamental es degenerado), de
manera que
D E Z
E0 = E [n0 ] = mı́n c Ψ + d3~r vext (~r) n0 (~r) .
Ψ Tb + W
|Ψi→n0 (~
r)
Como consecuencia, si n0 (~r) es densidad fundamental de algún sistema físico de N electrones, el
el valor esperado de la energía
UNED
funcional de Hohenberg-Kohn FHK [n0 ] puede interpretarse como
cinética más la energía potencial de interacción sobre el estado fundamental de densidad
n0 (~r).
2. Existe completa libertad a la hora de descomponer el funcional de Hohenberg-Kohn FHK [n] con
la idea (razonable) de desarrollar aproximaciones adecuadas para cada parte resultante. La estrategia
más habitual (aunque no la única) es la llamada implementación de Kohn y Sham de la DFT. La
UNED
idea básica es de nir el funcional de la densidad para la energía cinética ( o, simplemente, funcional
cinético) como
D E
TS [n] = mı́n Ψ Tb Ψ (6.1)
|Ψi→n
Si consideramos sistemas de electrones no interactuantes (W
c = 0), TS [n] sería el equivalente al
funcional de Hohenberg-Kohn FHK [n]. Por tanto los resultados principales del tema anterior son
UNED
igualmente aplicables. En concreto, si n0 es densidad de un estado fundamental |ΦS i de un sistema de
electrones no interactuantes,
D E
TS [n0 ] = ΦS Tb ΦS (6.2)
por lo que TS [n0 ] la energía cinética fundamental de un sistema no
puede interpretarse como
interactuante cuya densidad fundamental es n0 (~r) .
UNED
Consecuencia inmediata de la de nición (6.1) es el hecho de que si |Ψ0 i es el estado fundamental
de un sistema electrónico interactuante de densidad n0 (~r),
D E
TS [n0 ] ≤ Ψ0 Tb Ψ0 ,
D E
y TS [n0 ] puede entenderse también como una aproximación al valor esperado Ψ0 Tb Ψ0 bajo una
prescripción de partículas independientes.
3. Siendo TS [n0 ] una aproximación a la energía cinética fundamental, el funcional de Hartree
Z
1 n (~r1 ) n (~r2 )
WH [n] = d3~r1 d3~r2 (6.3)
2 |~r2 − ~r1 |
143
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
es, igualmente, una aproximación a la energía de interacción de un estado cuya densidad es n (~r).
Por tanto, si de nimos el funcional de la densidad para la energía de intercambio-correlación (o
simplemente funcional XC) como
UNED
D E D E
EXC [n] = FHK [n] − TS [n] − WH [n] = mı́n c Ψ − mı́n Φ Tb Φ − WH [n]
Ψ Tb + W (6.4)
|Ψi→n |Φi→n
vemos que EXC [n] contiene los efectos de muchos cuerpos que están ausentes al aproximar la energía
cinética y de interacción por TS [n] y WH [n], respectivamente. Observemos que, en este contexto, el
método de Hartree es equivalente a despreciar la energía XC EXC [n].
UNED
6.1.b. Expresión exacta para el funcional cinético
4. Consideremos una densidad N -electrónica n0 (~r) correspondiente al estado fundamental |ΦS i de
un sistema de electrones no interactuantes y sea vS (~r) el potencial exterior de dicho sistema. De
acuerdo con la discusión de la subsección anterior
UNED
D E XN N
X Z
2
TS [n0 ] = ΦS T ΦS =
b φi t φi =
b εi − d~
κ |φi (~
κ )| vS (~r) (6.5)
i=1 i=1
donde φi (~
κ) es eli-ésimo orbital del Hamiltoniano monopartícula t + vbS
b con autoenergía εi . Puesto
que, insistamos, n0 (~r) es la densidad fundamental,
UNED
XX
n0 (~r) = κ )|2 .
|φi (~ (6.6)
σ i=1
Por tanto, para una densidad fundamental n0 (~r) de un sistema de partículas no interactuantes, TS [n0 ]
se puede calcular exactamente si somos capaces de obtener el potencial vS (~ r) de cuyo Hamiltoniano no
interactuante asociado n0 (~
r) es densidad fundamental. Esto último es posible numéricamente mediante
un proceso iterativo, por lo que la expresión (6.5) debe entenderse como la forma operativa exacta del
funcional cinético.
UNED
5. Si ahora consideramos el mismo sistema no interactuante pero con N − ∆ partículas, el estado
fundamental será el mismo excepto que la ocupación del último orbital, φN , es 1 − ∆. Análogamente,
si el número de partículas es N + ∆ el estado fundamental será el mismo excepto que ahora el orbital
φN +1 tiene una ocupación parcial ∆. Por tanto:
UNED
δTS [n] εN − vS (~r) si hN
bi = N − ∆
= ; ∆ = 0+ (6.7)
δn (~r) n=n0 εN +1 − vS (~r) si hN
bi = N + ∆
observando que la derivada variacional de TS [n] es discontinua cuando el valor esperado del número
de partículas es entero.
6.1.c. Las ecuaciones de Kohn y Sham
6. Volvamos a nuestro sistema de N electrones de verdad bajo el potencial vext (~r) y sea n0 (~r) la
densidad fundamental. La ecuación de Euler-Lagrange que ha de veri car n0 (~r) es
Z
δFHK [n] δTS [n] δEXC [n] n (~r)
−I = + vext (~r) = + + d3~r 0 + vext (~r)
δn (~r) δn (~r) δn (~r) |~r − ~r 0 |
144
UNED
Método de Kohn-Sham I
Máster FSC. Funcional de la densidad
donde el penúltimo término es el potencial de Hartree en el punto ~r y, a partir de ahora y mientras
no se diga lo contrario, todas las derivadas variacionales se toman considerando el límite inferior al
número de partículas (lo que nos ha permitido igualar µ al opuesto del potencial de ionización I ).
UNED
Supongamos que n0 (~r) es también densidad fundamental de un sistema de N electrones no
interactuantes bajo el potencial vS (~r). De ser esto cierto, teniendo en cuenta (6.7),
Z
δEXC [n] 3 0 n (~r)
−I − εN = + d ~r + vext (~r) − vS (~r) .
δn (~r) |~r − ~r 0 | n=n0
El miembro de la izquierda es una constante, mientras que el de la derecha es una función
UNED
r. Así, puesto que vS (~r) está de nido salvo constante aditiva,1 la igualdad anterior con rma
de ~
nuestra suposición: n0 (~r) es también densidad fundamental de un sistema cticio de partículas no
interactuantes bajo el potencial
Z
δEXC [n] n0 (~r)
vS (~r) = + d3~r 0 + vext (~r) (6.8)
δn (~r) n0 (~
r) |~r − ~r 0 |
UNED
y además, la energía del N -ésimo orbital del sistema cticio es igual (salvo signo) a la energía de
ionización exacta del sistema real:
εN = −I (6.9)
Este resultado, conocido como el Teorema de Kohn y Sham [Phys. Rev 140, A1133 (1965)]
sugiere cuál es el método de resolución de la ecuación de Euler-Lagrange bajo el tratamiento exacto
del funcional κ )} y {εi } son los orbitales y autoenergías del sistema cticio de partículas
TS [n]. Si {φi (~
no interactuantes (denominado, a partir de ahora, sistema KS), el conjunto de ecuaciones acopladas
UNED
(ecuaciones KS)
1~2
− ∇ + vXC (~r) + vH (~r) + vext (~r) φi (~
κ ) = εi φi (~
κ ) ; i = 1, ..., N
2
Z
δEXC [n] n (~r)
vXC (~r) = ; vH (~r) = d3~r 0
δn (~r) |~r − ~r 0 |
N X
X
κ )|2
|φi (~
UNED
n (~r) = (6.10)
i=1 σ
determinan, simultáneamente, la densidad fundamental del sistema n0 (~r) y el sistema cticio KS
asociado. La energía del estado fundamental es, entonces
N Z Z
X
3 1 n0 (~r1 ) n0 (~r2 )
E0 = φi b
t φi +
d ~r n0 (~r) vext (~r) + d3~r1 d3~r2 + EXC [n0 ] =
2 |~r1 − ~r2 |
i=1
UNED
N Z Z
X
3 1 n0 (~r1 ) n0 (~r2 )
= εi − d ~r n0 (~r) vXC (~r) − d3~r1 d3~r2 + EXC [n0 ] (6.11)
2 |~r1 − ~r2 |
i=1
Las ecuaciones KS son un conjunto de ecuaciones autoconsistentes, al igual que las ecuaciones
HF. La diferencia fundamental radica en que, en principio, las ecuaciones KS nos dan exactamente
la energía y la densidad fundamental del sistema siempre y cuando conozcamos el funcional de
intercambio-correlación EXC [n] (el cual, en la práctica, ha de aproximarse). A su vez, en el
hamiltoniano KS el potencial que aparece es estrictamente local, a diferencia del potencial efectivo HF
que era no local al incluir el operador de Fock.
1
Si añadimos una constante al potencial exterior únicamente estamos rede niendo el origen de energías.
145
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
7. El método Kohn-Sham puede entenderse como un procedimiento bajo el cual se transforma el
sistema real de electrones interactuantes en un sistema cticio de electrones no interactuantes pero que
comparten la misma densidad y el mismo potencial de ionización. La conexión entre ambos sistema
UNED
se hace a partir del potencial de Hartree vH (~r) y del llamado potencial XC vXC (~r) que es la derivada
variacional del funcional para la energía XC. Puesto que el sistema cticio KS no es interactuante,
todos los efectos de muchos cuerpos han de estar contenidos, de alguna forma, en vXC (~r). Esto anticipa
que la forma exacta de EXC [n] y de su derivada variacional hayan de ser altamente no triviales.
Si el sistema electrónico es degenerado, las ecuaciones KS pueden tener diferentes soluciones pero
todas ellas dando lugar a la misma energía E0 . En la práctica, debido al desconocimiento del funcional
XC, este tipo de situaciones plantean problemas serios ya que no es extraño que las diferentes soluciones
UNED
autoconsistentes a las ecuaciones KS tengan energías muy parecidas, pero diferentes.
8. Como hemos visto, para un sistema de N electrones la energía εN del sistema cticio KS (energía
del nivel HOMO -highest occupied molecular orbital-) es igual a −I , el opuesto de la energía de
ionización del sistema de N electrones. Cabría preguntarse si εN +1 (energía del nivel LUMO -lowest
unoccupied molecular orbital-) es igual a −A (el opuesto de la a nidad electrónica). La respuesta
UNED
sería a rmativa si y solo si la derivada variacional del potencial XC fuese estrictamente continua, es
decir, que TS [n] absorbiese toda la discontinuidad del funcional de Hohenberg y Kohn FHK [n] que
aparece cuando n es densidad fundamental de un sistema de N partículas (N entero). En general
podemos a rmar que no es el caso: calculos tipo CI o mecano-cuánticos muy avanzados han permitido
obtener la densidad fundamental exacta para una gran variedad de sistemas electrónicos, al igual que
su energía, su potencial de ionización y su a nidad electrónica. El conocimiento de la densidad exacta
permite, como hemos dicho anteriormente, evaluar el potencial KS asociado vS (~r). Varios autores han
comprobado que en este caso el nivel HOMO del sistema cticio KS es igual (dentro de las inevitables
UNED
barras de error numéricas) a menos la energía de ionización exacta pero, por el contrario, el nivel
LUMO del sistema cticio KS no es igual a menos la a nidad electrónica del sistema de N electrones.
Esta discusión es de especial importancia en la interpretación de los resultados KS en semiconductores
y aislantes.
9. La implementación práctica del método KS exige únicamente aproximar el funcional de la energía
XC. Como ya hemos comentado, EXC [n] ha de ser altamente no trivial, puesto que el funcional ha
UNED
de describir todos los efectos de muchos cuerpos asociados a la interacción entre los electrones. Sin
embargo, aproximaciones muy simples a EXC [n] permiten obtener excelentes resultados en una gama
muy amplia de sistemas físicos, lo que explica la popularidad de la DFT.
No conviene olvidar que el método KS transforma el problema interactuante en otro equivalente
no interactuante. La losofía original del teorema de Hohenberg y Kohn reducía el problema aun mas,
al enfatizar el carácter fundamental de la densidad (algo escondido en el formalismo de KS, ya que la
UNED
expresión exacta de TS [n] se obtiene mediante la dependencia funcional implícita del potencial efectivo
2
KS vS (~r) con la densidad ). Es por ello por lo que, a pesar de la formulación de Kohn y Sham, el
desarrollo de aproximaciones funcionales razonablemente precisas para TS [n] con dependencia explícita
en la densidad siga siendo tema activo de investigación.
6.1.d. Sistemas con dependencia en el espín
10. De acuerdo con las ecuaciones Kohn-Sham para un potencial exterior vext (~r), el potencial efectivo
KS y, por ende, el potencial de intercambio-correlación vXC (~r) no exhiben una dependencia explícita
en el espín. Esto quiere decir que en el sistema cticio KS el potencial efectivo que actúa sobre los
2
Es por ello por lo que también se suele escribir vS [n] (~
r) para enfatizar que el potencial del sistema cticio KS es un funcional
de la densidad.
146
UNED
Método de Kohn-Sham I
Máster FSC. Funcional de la densidad
orbitales con espín hacia arriba es el mismo que el que actúa sobre los orbitales con espín hacia abajo.
Entonces, si el estado fundamental posee espín total no nulo (por ejemplo, el átomo de nitrógeno)
la implementación del formalismo KS tiene una di cultad añadida: hay que intentar describir este
UNED
hecho mediante un funcional (el de intercambio-correlación) que depende únicamente de la densidad
total. En este sentido, el valor esperado del espín es una propiedad más del estado fundamental, como
puede ser el valor esperado del momento angular y, de acuerdo con el teorema de Hohenberg y Kohn,
puede expresarse como un funcional de la densidad total. Dicha dependencia funcional no es trivial
y, de hecho, uno no puede esperar a priori que el valor esperado del espín del sistema real sea igual
a dicho valor esperado en el sistema cticio KS, de igual manera que no se puede asegurar que los
valores esperados del momento angular coincidan en el sistema real y en el cticio KS. Insistamos que
UNED
lo único que podemos a rmar es que, tratando exactamente el intercambio-correlación, la densidad y
la energía de ionización de ambos sistemas son iguales.
11. Esta di cultad puede solventarse si reformulamos la DFT pero considerando que la variable
fundamental es la densidad dependiente del espín, no la densidad total. A su vez, esta generalización
incorpora de manera natural situciones en las que el potencial exterior depende del espín. Así,
rede namos los funcionales de Hohenberg-Kohn, cinético y de intercambio-correlación como
UNED
D E
FHK [n+ , n− ] = mı́n Ψ Tb + W
c Ψ
|Ψi→nσ (~
r)
D E
TS [n+ , n− ] = mı́n Φ Tb + W
c Φ (6.12)
|Φi→nσ (~
r)
EXC [n+ , n− ] = FHK [n+ , n− ] − TS [n+ , n− ] − WH [n+ + n− ]
Obsérvese que las de niciones de FHK [n+ , n− ] y TS [n+ , n− ] son ahora más restrictivas, ya que el
UNED
mínimo se obtiene sobre todos los estados cuya densidad dependiente del espín es nσ (~r), mientras
que en los funcionales originales el mínimo se de nía sobre todos los estados cuya densidad total es
n+ (~r) + n− (~r). De este modo
FHK [n+ , n− ] ≥ FHK [n+ + n− ]
TS [n+ , n− ] ≥ TS [n+ + n− ]
y, en general
UNED
EXC [n+ , n− ] 6= EXC [n+ + n− ]
Sin embargo, si nσ (~r) es densidad fundamental de un sistema real en el que el potencial exterior no
depende del espín, los dos funcionales HK dan el mismo resultado: el valor esperado de la energía
cinética y de interacción en el estado fundamental correspondiente. Igualmente, si nσ (~r) es densidad
fundamental de un sistema no interactuante con potencial independiente del espín, los dos funcionales
cinéticos coinciden y son iguales al valor esperado de la energía cinética en el estado fundamental
3
correspondiente. El único funcional que, en cualquier caso, es idéntico en ambas formulaciones es
UNED
WH [n], la energía electrostática clásica de interacción.
Ahora el potencial XC actúa de diferente manera sobre los orbitales KS dependiendo de su espín.
En concreto
(±) δEXC [n+ , n− ]
vXC (~r) = (6.13)
δn± (~r)
y las ecuaciones que obedecen los orbitales KS son
1~2 (±) (±)
− ∇ + vext (~r) + vH (~r) + vXC (~r) φ(±)
n (~r) = ε(±) (±)
n φn (~ r) (6.14)
2
3
Aquí surge un problema sobre el que apenas hemos discutido. Si r) 6= n− (~
n+ (~ r) es densidad fundamental de un sistema
no interactuante con potencial independiente del espín, necesariamente este tiene una con guración fundamental de capa abierta,
siendo degenerado Así, el estado fundamental no tiene porque ser un único determinante de Slater, sino una combinación lineal de
varios determinantes de Slater con igual energía. Esto complica algo el formalismo, pero las ideas esenciales permanecen.
147
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
(σ) (σ)
donde φn (~r) es la parte espacial del orbital n-ésimo KS con espín σ, siendo εn su autoenergía
(obsérvese cómo vH (~r) no depende del espín, puesto que es el potencial electrostático clásico
creado por la densidad total). Los orbitales KS ocupados son los N orbitales con menor energía,
UNED
independientemente de su espín. Así, si hay N+ orbitales ocupados con espín + y N− orbitales con
espín −,
N±
X 2
n± (~r) = φ(±)
n (~r) (6.15)
n=1
cerrándose el ciclo autoconsistente. Alcanzada la autoconsistencia, la energía fundamental del sistema
es
UNED
E0 = TS [n+ , n− ] + Vext [n+ , n− ] + EXC [n+ , n− ] + WH [n] = (6.16)
Nσ D E XZ
(σ)
XX
= φ(±)
n t
b φ (±)
n + d3~r vext (~r) nσ (~r) + EXC [n+ , n− ] + WH [n+ + n− ]
σ n=1 σ
12. Insistamos que si el potencial es independiente del espín, en el formalimo exacto los resultados
son los mismos bajo la prescripción original y bajo la implementación dependiente del espín. Sin
UNED
embargo ya hemos visto que los funcionales EXC [n+ , n− ] y EXC [n] son, por de nición, claramente
diferentes. En la construcción de EXC [n+ , n− ] hay más datos y resulta más sencillo de modelizar
que el funcional EXC [n], en el que la única entrada es la densidad total. Además, sí podemos a rmar
que la densidad dependiente del espín es igual en el sistema real y en el cticio KS y, como consecuencia
inmediata, el valor esperado de la tercera componente del espín es igual en ambos sistemas, el real y
el cticio KS.
UNED
6.2. PROPIEDADES EXACTAS DEL FUNCIONAL XC
1. Como hemos visto, la implementación de Kohn y Sham de la DFT permite obtener exactamente,
y de manera directa a través de las ecuaciones de KS, la densidad y la energía del estado
fundamental de un sistema de electrones. Para ello es preciso conocer el funcional de la energía
de intercambio-correlación EXC [n] y su derivada variacional vXC (~r). Ahora bien, su expresión es
desconocida, por lo que resulta imprescindible construir aproximaciones razonables adecuadas para su
UNED
aplicación a sistemas de interés.
La gran di cultad del funcional XC es el hecho de que contiene, simultáneamente, correcciones a
la energía cinética y a la energía de interacción. En el primer caso podemos decir que las correcciones
están debidas a la ausencia de correlación coulombiana al aproximar la energía cinética mediante la
correspondiente a un determinante de Slater cuya densidad es igual a la densidad fundamental del
sistema. En el segundo caso hay que tener en cuenta efectos de intercambio (correlación de Pauli)
y de correlación (de Coulomb), lo que explica el nombre dado al funcional. Es por ello por lo que
UNED
resulta muy interesante conocer algunas propiedades de este funcional, que sirvan de guía tanto a
la hora de desarrollar aproximaciones funcionales como para juzgar la pertinencia o no de algunas
de las aproximaciones que se utilizan habitualmente. Este es el objetivo de esta sección, en la que
nos limitaremos a la formulación original sin dependencia en el espín, aunque la generalización es
inmediata.
6.2.a. Relación con la teoría de la respuesta lineal
2. En el tema anterior hemos visto que hay una relación entre la segunda derivada variacional del
funcional de Hohenberg-Kohn y la función respuesta estática del sistema electrónico:
δ 2 FHK [n]
KHK (~r1 , ~r2 ) = = −χ−1 (~r1 , ~r2 ) (6.17)
δn (~r1 ) δn (~r2 ) n=n0
148
UNED
Método de Kohn-Sham I
Máster FSC. Funcional de la densidad
Se puede hacer un razonamiento idéntico pero para sistemas no interactuantes. Recordando que en
este caso FHK [n] se reduce a TS [n], tendríamos que
UNED
δ 2 TS [n]
= −χ−1
0 (~
r1 , ~r2 ) (6.18)
δn (~r1 ) δn (~r2 ) n=n0
donde χ0 (~r1 , ~r2 ) es la función respuesta estática del sistema cticio KS asociado a la densidad n0 (~r).
A partir de la expresión que vimos en su día para la respuesta lineal de un sistema no interactuante,
se tiene que
XX φ∗j (~r 1 σ) φj (~r2 σ) φi (~r1 σ) φ∗i (~r2 σ)
χ0 (~r1 , ~r2 ) = (fj − fi )
UNED
σ
εj − εi
j,i
con fi = 0 ó 1 si el orbital KS está desocupado u ocupado y donde hemos extendido el resultado ya
conocido al caso general de dependencia en el espín.
3. Puesto que FHK [n] es igual a la suma de TS , EXC y WH [n] el llamado kernel estático XC, de nido
como
δ 2 EXC [n]
UNED
KXC (~r1 , ~r2 ) =
δn (~r1 ) δn (~r2 ) n=n0
ha de cumplir
δ 2 EXC [n] 1
KXC (~r1 , ~r2 ) = = χ−1 r1 , ~r2 ) −
0 (~ − χ−1 (~r1 , ~r2 ) (6.19)
δn (~r1 ) δn (~r2 ) n=n0 |~r1 − ~r2 |
donde hemos usado el resultado δ 2 WH [n] / (δn (~r1 ) δn (~r2 )) = 1/ |~r1 − ~r2 |. Esta igualdad es muy
UNED
importante, puesto que relaciona el kérnel XC con la función respuesta del sistema interactuante.
Así, toda aproximación a EXC [n] implica, inmediatamente, una aproximación a la función respuesta
estática que incluye efectos más allá de la RPA. Además, como veremos más adelante, este tipo de
relaciones son muy útiles a la hora de diseñar aproximaciones funcionales a EXC [n] y TS [n].
6.2.b. La conexión adiabática
UNED
4. Recordemos que si n (~r) |Ψ0 i de un sistema físico, la
es la densidad del estado fundamental
energía XC se de nía exactamente como EXC [n] = hΨ0 | T + W |Ψ0 i − TS [n] − WH [n]. WH [n] es la
b c
energía de Hartree y TS [n] el funcional cinético dado por TS [n] = hΦS | T
bS |ΦS i, siendo |ΦS i el estado
fundamental de un sistema no interactuante (el sistema KS) cuya densidad fundamental es también
n (~r). De acuerdo con los teoremas de Hohenberg, Kohn y Sham, ambos estados fundamentales son,
en rigor, funcionales de la densidad n (~r).
UNED
Desde una perspectiva puramente teórica (incluyendo, por ejemplo, la aplicación de la teoría de
perturbaciones de muchos cuerpos) resulta más sencillo caracterizar los efectos mecanocuánticos en la
energía de interacción que en la energía cinética. La idea central de la llamada conexión adiabática es
determinar los efectos cuánticos en la energía cinética en función de aquellos que son relevantes en la
energía de interacción. Esta idea, propuesta originalmente por W. Pauli, no se limita a la propia DFT
y resulta una herramienta habitual en teoría de perturbaciones de muchos cuerpos.
Así, consideremos un sistema físico de electrones sometidos a la acción del potencial exterior vext (~r).
De namos un operador Hamiltoniano dependiente de un parámetro λ ∈ [0, 1] dado por
Z
b (λ) = Tb + λW
H c+ d3~r v (λ) (~r) n
b (~r)
tal que para cualquier valor de λ la densidad fundamental sea igual a la densidad fundamental n (~r)
de nuestro sistema. Observemos que, entonces, v (λ=1) (~r) es igual al potencial exterior de nuestro
149
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
E
(λ)
sistema, mientras que v (λ=0) (~r) es el potencial KS vS (~r) asociado a la densidad n (~r). Si Ψ0 es
el estado fundamental correspondiente al Hamiltoniano b (λ) ,
H tenemos que la energía fundamental
UNED
correspondiente es
D E D E Z
(λ) (λ) (λ) (λ) (λ) b (λ)
E0 = Ψ0 H
b Ψ0 = Ψ0 T + λW Ψ0
c + d3~r v (λ) (~r) n (~r)
donde la densidad fundamental, n (~r), no depende del parámetro λ ya que por construcción es la misma
cualquiera que sea el valor de λ. Si ahora aplicamos el teorema de Hellmann-Feynman,
* +
(λ) (λ)
dE0 (λ) ∂ H
b (λ)
UNED
= Ψ0 Ψ0 ,
dλ ∂λ
tenemos que
(λ)
dE0 D
(λ) (λ)
E Z ∂ v (λ) (~r)
= Ψ0 W Ψ0
c + d3~r n (~r)
dλ ∂λ
(obsérvese cómo ha desaparecido la energía cinética, ya que ésta no depende explícitamente del
parámetro λ). Integremos ahora esta expresión entre 0y 1, obteniendo
UNED
Z 1 D E Z Z
(λ=1) (λ=0) (λ) c Ψ(λ) + 3
E0 = E0 + dλ Ψ0 W 0 d ~r vext (~r) n (~r) − d3~r vS (~r) n (~r)
0
(λ=0)
Ahora bien, E0 es la energía total del sistema KS:
Z
(λ=0)
E0 = TS [n] + d3~r vS (~r) n (~r)
UNED
(λ=1)
mientras que E0 es la energía fundamental E0 de nuestro sistema físico. Eliminando el superíndice
λ = 1, que ahora es redundante, llegamos a que:
D E Z 1 D E
(λ) c (λ)
Ψ0 T + W Ψ0 = TS [n] +
b c dλ Ψ0 W Ψ0 (6.20)
0
Es decir, el valor esperado de la energía cinética se puede obtener a partir del conocimiento del
UNED
funcional cinético y de los valores esperados de la energía de interacción evaluados sobre los estados
fundamentales de una familia de sistemas cticios en los que la interacción es λWc y cuya densidad
fundamental sea igual a la del sistema real. Ya sólo falta sustituir hΨ0 | Tb + W
c |Ψ0 i en términos de los
funcionales cinético, XC y de Hartree para llegar a la siguiente expresión exacta para la energía de
intercambio-correlación
Z 1 D E
(λ) c (λ)
EXC [n] = dλ Ψ0 W Ψ0 − WH [n] (6.21)
0
UNED
de la cual ha desaparecido toda mención al operador de energía cinética.
5.
E
(λ)
Si h(λ) (~r1 , ~r2 ) es la función de correlación correspondiente al estado cuántico Ψ0 y
(λ)
= 21 n (~r1 ) n (~r2 ) 1 + h(λ) (~r1 , ~r2 )
ρ2 (~r1 , ~r2 ) la densidad de pares correspondiente, sabemos que
E Z (λ) Z
3 ρ2 (~ r1 , ~r2 ) 1 n (~r1 ) n (~r2 ) (λ)
D
(λ) (λ) 3
Ψ0 W Ψ0
c = d ~r1 d ~r2 = WH [n] + d3~r1 d3~r2 h (~r1 , ~r2 )
|~r2 − ~r1 | 2 |~r2 − ~r1 |
E
(λ)
(insistamos que la densidad es la misma para todos los estados Ψ0 . De esta manera
Z Z 1
1 n (~r1 ) n (~r2 )
EXC [n] = d3~r1 d3~r2 dλ h(λ) (~r1 , ~r2 ) (6.22)
2 |~r2 − ~r1 | 0
150
UNED
Método de Kohn-Sham I
Máster FSC. Funcional de la densidad
La función Z 1
hXC (~r1 , ~r2 ) = dλ h(λ) (~r1 , ~r2 ) (6.23)
0
UNED
es, entonces, una función de correlación de pares integrada (o efectiva) construida a partir de la
E
(λ)
integración de las funciones de correlación de pares de los estados Ψ0 de modo que
Z
1 n (~r1 ) n (~r2 ) hXC (~r1 , ~r2 )
EXC [n] = d3~r1 d3~r2 (6.24)
2 |~r2 − ~r1 |
Naturalmente podemos usar la función hueco en lugar de la función de correlación de pares. Así, si
UNED
de nimos el hueco de intercambio-correlación
E centrado en ~r como la integral de las funciones hueco
(λ)
correspondientes a los estados Ψ0 :
η XC ~r ; ~r 0 = n ~r 0 hXC ~r, ~r 0
(6.25)
tenemos otra expresión alternativa para la energía XC:
η XC (~r ; ~r 0 )
Z Z Z
1 0
EXC [n] = d3~r n (~r ) 3
d ~r ≡ d3~r n (~r ) εXC [n] (~r) (6.26)
UNED
2 |~r − ~r 0 |
donde de nimos εXC [n] (~r) como la energía de XC por partícula en el punto ~r (aunque la notación sea
algo confusa, preferimos mantener la dependencia funcional para así enfatizarla). Por tanto, la energía
de intercambio-correlación por partícula en el punto ~r es igual a un medio del potencial coulombiano
creado por el hueco de intercambio-correlación centrado en ~r.4 Naturalmente, el hueco XC hereda la
regla de suma de la función hueco y, por tanto,
Z
UNED
d3~r 0 η XC ~r ; ~r 0 = −1 ∀ ~r
(6.27)
Estas expresiones formales tienen su interés. En primer lugar sugieren que una buena forma de
modelizar el funcional XC es hacerlo a través del hueco XC (o, equivalentemente, de la función de
correlación de pares integrada). En segundo lugar nos permiten deducir algunas propiedades exactas
del funcional XC que nos pueden servir para juzgar la bondad de ciertas aproximaciones. En particular
consideremos un sistema electrónico localizado en torno a ~r = ~0. Si r→∞ el hueco XC centrado en
r 0 cercanos al origen ya que (salvo el detalle de la
UNED
~r va a permancer localizado en aquellos puntos ~
integración en el parámetro λ), el hueco XC nos informa sobre el cambio en la densidad electrónica
supuesta la existencia de un electrón en ~r, por lo que dicho cambio ha de ser relevante sólo en aquellas
regiones en las que la densidad de partida sea no nula [recuérdese que η XC (~r ; ~r 0 ) = n (~r 0 ) hXC (~r, ~r 0 )].
Por tanto
η (~r ; ~r 0 ) −1
Z
1
lı́m εXC [n] (~r) = lı́m d3~r 0 XC 0
= (6.28)
r→∞ 2 r→∞ |~r − ~r | 2r
ya que el potencial creado por una densidad localizada que integra a −1 se comporta como −1/r a
UNED
5
grandes distancias. A su vez (aunque la demostración es más técnica y la omitimos aquí)
1
lı́m vXC [n] (~r) = − (6.29)
r→∞ r
4
Recordemos que, en nuestro contexto, potencial es igual a energía potencial, por lo que el potencial al que nos estamos
re riendo es estrictamente el opuesto al potencial electrostático.
5
En todo caso, el esquema general de la demostración es el siguiente. Como
Z
1 n (~
r1 ) n (~
r2 ) hXC (~r1 , ~
r2 )
EXC [n] = d3 ~
r1 d3 ~
r2
2 |~
r2 − ~r1 |
tenemos que
r 0 ) hXC (~ r 0)
Z Z
δEXC [n] n (~ r, ~ 1 n (~
r1 ) n (~
r2 ) δ hXC (~r1 , ~
r2 )
vXC (~
r) = = ... = d3 ~
r1 d3 ~
r2 0
+ d3 ~
r1 d3 ~
r2
δn (~
r) |~r−~r | 2 |~
r2 − ~r1 | δ n (~
r)
El primer término es 2εXC [n] (~
r). El segundo (y esta es la parte técnica) decae a cero más rápidamente que r−1 cuando r 0.
151
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
6. Por otro lado, si recordamos el teorema de uctuación-disipación, tenemos que
UNED
Z +∞
Z πn (~r1 ) δ (~r2 − ~r1 ) + Im dω χ(λ) (~r1 , ~r2 ; ω)
D
(λ) c (λ)
E 1 3 3 0
Ψ0 W Ψ0 = WH [n] − d ~r1 d ~r2
2π |~r2 − ~r1 |
con χ(λ) es la función respuesta del sistema cticio de partículas que interaccionan a través de un
potencial λ/r12 y cuya densidad es igual a la densidad del sistema real. Por tanto, χ(λ=1) ≡ χ es la
función respuesta del sistema real, mientras que χ
(λ=0) ≡ χ0 es la función respuesta del sistema cticio
UNED
KS de partículas no interactuantes cuya expresión es inmediata a partir de los orbitales KS φi (~r σ) y
de sus autoenergías εi
XX φ∗j (~r1 σ) φj (~r2 σ) φi (~r1 σ) φ∗i (~r2 σ)
χ0 (~r1 , ~r2 ; ω) = (fj − fi )
σ
ω − (εi − εj ) + iη
j,i
con fi = 0 ó 1 si el orbital KS está desocupado u ocupado.
De esta manera
UNED
Z Z +∞ Z 1
1 n (~r1 ) 1 dω
EXC [n] = − d3~r1 d3~r2 δ (~r2 − ~r1 ) + Im dλχ(λ) (~r1 , ~r2 ; ω) (6.30)
2 |~r2 − ~r1 | n (~r1 ) 0 π 0
expresión conocida como ACFDT (de las siglás en inglés de conexión adiabática y teorema de
uctuación-disipación). Esto implica que la energía XC por partícula en el punto ~r1 es
Z Z +∞ Z 1
1 1 1 dω
UNED
3 (λ)
εXC [n] (~r1 ) = − d ~r2 δ (~r2 − ~r1 ) + Im dλ χ (~r1 , ~r2 ; ω) (6.31)
2 |~r2 − ~r1 | n (~r1 ) 0 π 0
y que la función de correlación de pares integrada y el hueco XC son, respectivamente,
Z +∞ Z 1
1 dω (λ)
η XC (~r1 ; ~r2 ) = − δ (~r2 − ~r1 ) + Im dλχ (~r1 , ~r2 ; ω)
n (~r1 ) 0 π 0
Z +∞ Z 1
1 1 dω (λ)
hXC (~r1 , ~r2 ) = − δ (~r2 − ~r1 ) + Im dλχ (~r1 , ~r2 ; ω)
UNED
n (~r2 ) n (~r1 ) n (~r2 ) 0 π 0
Esto sugiere que la modelización de la función respuesta (por ejemplo usando la RPA) es otra forma
para obtener la energía XC de un sistema de electrones.
6.2.c. Funcional de intercambio y funcional de correlación
7.
UNED
Aunque la discriminación de los efectos de intercambio y correlación es, en cierta medida,
arbitraria, podemos implementar dicha división en el funcional de la energía de intercambio-correlación.
Ya que el sistema cticio KS tiene en cuenta el principio de exclusión, estando ausente cualquier
mención directa a la energía de interacción, podemos de nir el funcional de intercambio EX [n] como
la contribución asociada exclusivamente al sistema cticio KS en la expresión ACFDT del funcional
EXC [n]. Así
Z Z Z +∞
3 1 3 1 1 dω
EX [n] = − d ~r1 n (~r1 ) d ~r2 δ (~r2 − ~r1 ) + Im χ (~r1 , ~r2 ; ω)
2 |~r2 − ~r1 | n (~r1 ) 0 π 0
y, por tanto, el funcional de la energía de correlación será
Z Z +∞
1 1 dω (λ)
d3~r1 d3~r2
EC [n] = − Im χ (~r1 , ~r2 ; ω) − χ0 (~r1 , ~r2 ; ω) (6.32)
2 |~r2 − ~r1 | 0 π
152
UNED
Método de Kohn-Sham I
Máster FSC. Funcional de la densidad
Aunque tedioso, es sencillo comprobar que la anterior expresión para el funcional de intercambio es
exactamente igual a
UNED
occ ∗
φj (~r1 σ) φj (~r2 σ) φi (~r1 σ) φ∗i (~r2 σ)
Z
1 XX
EX [n] = − d3~r1 d3~r2 (6.33)
2 σ
|~r2 − ~r1 |
i,j
es decir, la energía de intercambio evaluada sobre el determinante de Slater |ΨS i que es el estado
fundamental del sistema cticio KS.
Como veremos, el conocimiento de la expresión exacta para el funcional de intercambio abre el
UNED
camino a una implementación muy precisa del formalismo KS, aunque con la di cultad de la muy
complicada evaluación de la derivada variacional de EX [n] (ya que la densidad aparece de forma
implícita, no explícita).
8. Si de nimos hX (~r1 , ~r2 ) como la función de correlación de pares del sistema cticio KS
occ
1 XX
hX (~r1 , ~r2 ) = − φ∗j (~r1 σ) φj (~r2 σ) φi (~r1 σ) φ∗i (~r2 σ) (6.34)
UNED
n (~r1 ) n (~r2 ) σ
i,j
inmediatamente,
Z
1 n (~r1 ) n (~r2 ) hX (~r1 , ~r2 )
EX [n] = d3~r1 d3~r2 (6.35)
2 |~r1 − ~r 2 |
y hX (~r1 , ~r2 ) puede interpretarse como la contribución de intercambio a la función de correlación de
pares integrada hXC (~r 1 , ~r2 ).
UNED
En general, hX (~r1 , ~r2 ) tiene más alcance que hXC (~r 1 , ~r2 ), en el sentido de que hX (~r1 , ~r2 ) decrece
más lentamente que hXC (~ r 1 , ~r2 ) a medida que |~
r2 − ~r1 | → ∞. Ello es lógico puesto que en un
sistema interactuante el fenómeno de apantallamiento hace que el alcance efectivo de la interacción
coulombiana sea menor. Por tanto, aunque las energías de intercambio y correlación son negativas, sus
efectos no locales tienden a compensarse, en el sentido de que las contribuciones con|~r2 − ~r1 | 0 a
EX [n] y EC [n] en buena parte cancelan entre sí. Esta compensación no es muy importante en sistemas
atómicos, en los que la energía de intercambio predomina sobre la correlación. Sin embargo es esencial
UNED
en sistemas extensos (en los que el fenómeno de apantallamiento es importante) y también a la hora
de analizar en detalle el enlace químico, ya que las correlaciones entre los electrones pertenecientes
a átomos diferentes son cruciales.
6.3. LAS APROXIMACIONES LDA Y GGA
6.3.a. Algunas propiedades del gas homogéneo
UNED
1. El método Kohn-Sham reduce la aplicación práctica de la DFT a la necesidad de aproximar
la energía de intercambio-correlación EXC [n]. En esta sección discutiremos en profundidad la
aproximación más sencilla (y popular) al funcional XC: la llamada aproximación de densidad local
(LDA). La idea central de la LDA es modelizar la energía de intercambio y correlación de un sistema
arbitrario basándose en las propiedades del gas homogéneo. Conviene por ello recordar algunas de las
propiedades que ya hemos visto para dicho sistema. Así, consideremos un gas homogéneo de electrones
con densidad n, en principio en su fase paramagnética. El radio de Wigner y el momento de Fermi
venían dados, respectivamente, por
1/3 1/3
2
1/3 1/3 3 9π 1
kF = 3π n ' 3,094 n ; rs = =
4πn 4 kF
153
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
UNED ε;& UV +D
,QWHUFDPELR
&RUUHODFLyQ
UNED
UV
Figura 6.1. Energia de intercambio y energía de correlación por partícula para el gas homogéneo de electrones. Los valores
correspondientes a densidades metálicas están destacados.
2. En el lenguaje de la DFT, la energía cinética no interactuante del gas homogéneo por partícula
UNED
era
2 2/3
3 (9π/4)2/3 1
3 3 3π
thom
S (n) = kF2 = n2/3 =
10 10 10 rs2
mientras que la energía de intercambio por partícula era
1/3
3 3 3π 2 3 (9π/4)1/3 1
εhom
X (n) = − kF = − n1/3 = − (6.36)
4π 4π 4π rs
UNED
En este caso, la energía de correlación no admite forma analítica. Sin embargo, mediante simulaciones
numéricas Monte Carlo complementadas con técnicas teóricas más avanzadas se puede obtener la
energía total por partícula εhom
tot (n). Recordando que en el gas homogéneo la energía total de interacción
electrostática clásica es nula, la energía de correlación por partícula es
εC (n) = εhom hom
tot (n) − εX (n) − thom
S (n)
UNED
función que puede parametrizarse en la forma debida a Perdew y Wang [Phys. Rev. B 45, 13244
(1992)]:
1
εhom
C (n) = −2c0 (1 + α1 rs ) ln 1 + (6.37)
1/2 3/2
2c0 β 1 rs + β 2 rs + β 3 rs + β 4 rs2
con c0 = 0,031091, α1 = 0,21370, β 1 = 7,595 6, β 2 = 3,5875, β 3 = 1,6382, β 4 = 0,49294. Hay otras
UNED
parametrizaciones, pero esta es la más aceptada. Observemos que si la densidad es grande (rs pequeño),
εhom
C (n 0) ' 2c0 ln rs , que siempre es menor (en valor absoluto) que la energía de intercambio: a
densidades altas la correlación es despreciable frente al interecambio. Sin embargo, a densidades altas
(rs grande), εhom
C (n 0) ' −α1 (β 4 rs )−1 y la energía de correlación es del mismo orden que la de
intercambio, aunque siempre menor en valor absoluto.
3. Es igualmente interesante P
resumir las propiedades de respuesta del gas homogéneo frente a una
perturbación estática δv (~r) = q ) exp (i~q · ~r). De acuerdo con la teoría de respuesta lineal, la
q~ δv (~
componente de Fourier ~q de la densidad variará en la forma δn (~q) = χ (q) δv (~q), donde χ (q) δ q q0 es
la representación matricial de la función respuesta del gas homogéneo en el espacio recíproco. En su
momento ya vimos la expresión para χ (q, ω) tanto en una aproximación de electrones no interactuantes
como dentro del esquema RPA.
154
UNED
Método de Kohn-Sham I
Máster FSC. Funcional de la densidad
UNED
1RLQWHUDFWXDQWH
χ T DX
,QWHUDFWXDQWH 53$
UV
N)
NV
T DX
UNED
Figura 6.2. Funciones respuestas estáticas interactuante (RPA) y no interactuante para un gas homogéneo con densidad
rs = 3,83. Obsérvese el muy diferente comportamiento frente a perturbaciones de larga longitud de onda.
Si consideramos una aproximación de electrones no interactuantes, el único efecto cuántico que
estamos teniendo en cuenta es el intercambio (Principio de Exclusión). En ese caso
UNED
q 2
q
1− 1+
kF
1 2kF 2kF
χ (q) ' χ0 (q) = − 2 + ln (6.38)
π 2
q
q
4 1−
2kF 2kF
observando que
UNED
2
kF 8kF kF
χ0 (q ∼ 0) ' − 2 ; χ0 (q 0) ' − 2
π 3π q2
Si ahora abordamos el problema usando la aproximación RPA
χ0 (q)
χ (q) = (6.39)
4π
1 − 2 χ0 (q)
q
UNED
de manera que si q∼0
2 r 1/6
q2
q 4 9π 1
χ (q ∼ 0) ' − = −χ0 (q) ; ks = √ ' 1,985 × n1/6 (6.40)
4π ks π 4 rs
Este es un resultado fundamental. Al considerar la interacción entre los electrones, un potencial con
longitud de onda larga apenas tiene efecto en la densidad del sistema debido, naturalmente, al fenómeno
UNED
de apantallamiento. La redistribución de carga es muy pequeña porque el potencial coulombiano creado
por la misma apantalla la perturbación exterior. Por otro lado, si q 0
2
8kF kF
χ (q 0) ' χ0 (q 0) = − 2
3π q2
por lo que perturbaciones de corta longitud de onda apenas son apantalladas por el gas homogéneo y
éste responde como si fuese no interactuante. Las diferencias entre la función respuesta no interactuante
y la RPA pueden observarse en la gura 6.2 para un gas homogéneo con rs = 3,84.
Por último, si tenemos en cuenta la relación general
1
KXC (~r1 , ~r2 ) = χ0−1 (~r1 , ~r2 ) − − χ−1 (~r1 , ~r2 )
|~r1 − ~r2 |
155
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
con KXC (~r1 , ~r2 ) la segunda derivada variacional de EXC [n], para el gas homogéneo se veri ca que
χ0 (q)
χ (q) = (6.41)
UNED
4π hom
1− + KXC (q) χ0 (q)
q2
donde
hom (q) δ 0
KXC es la representación matricial en el espacio recíproco de KXC (~r1 , ~r2 ) para el
qq
gas homogéneo. La inclusión de los efectos más allá de la RPA cambia cuantitativamente, aunque
no cualitativamente, los resultados. De hecho, la función
hom (q)
KXC puede obtenerse usando teoría
avanzada y simulaciones numéricas, pudiendo demostrarse que
UNED
hom d2 h hom i
lı́m KXC (q) = nε XC (n) <0 (6.42)
q→0 dn2
expresión conocida como teorema de compresibilidad y que
hom
lı́m KXC (q) = cte
q→∞
Así, el comportamiento de la función χ (q) exacta es similar a la RPA tanto si q'0 como si q 0.
UNED
4. Los resultados anteriores sugieren la existencia de dos longitudes naturales en el gas homogéneo. La
primera sería λF = kF−1 , que es del orden de rs y está asociada a los efectos de intercambio. La segunda
−1 √
sería λs = ks , del orden de rs (recuérdese que estamos trabajando en unidades atómicas) y estando
asociada a los efectos de correlación. De hecho a λs se le denomina longitud de apantallamiento puesto
que toda perturbación cuya longitud de onda sea sustancialmente mayor que λs será prácticamente
apantallada por la (pequeña) redistribución de carga.
UNED
6.3.b. La aproximación local a la energía de intercambio-correlación (LDA)
5. La aproximación de densidad local (LDA de sus siglas en inglés) al funcional de
intercambio-correlación se de ne como
Z
EXC [n] ' LDA [n]
EXC = d3~r n (~r) εhom
XC (n (~
r)) (6.43)
UNED
esto es, la energía XC en el punto ~r se aproxima por el valor correspondiente a un gas homogéneo de
densidad n (~r). Así, se supone que el comportamiento local del sistema electrónico es equivalente al
del gas homogéneo, de ahí el nombre de la aproximación.
La derivada variacional de
LDA [n],
EXC igual al potencial XC en la LDA, es muy sencilla. Usando las
reglas básicas de la derivación variacional tenemos que
UNED
LDA (~ ∂ εhom
XC (n (~ r))
vXC r) = εhom
XC (n (~
r)) + n (~r) (6.44)
∂ n (~r)
La LDA puede parecer excesivamente simple a la vista de la complejidad de la energía de
intercambio-correlación real. Sin embargo ofrece resultados sorprendentemente buenos, especialmente
para sistemas extensos de interés en Física del Estado Sólido. La LDA es también una aproximación
razonable en el ámbito de la Química Cuántica (átomos, moléculas, agregados), aunque no puede
competir con otros métodos más precisos que, sin embargo, son númericamente más costosos que
el método Kohn-Sham bajo la LDA. De hecho, el éxito de la LDA es tal que, en muchos casos, la
aproximación LDA se sobrentienda al hablar del formalismo KS.
La ventaja indudable de la LDA es su sencillez. La energía y el potencial XC se reducen a meras
expresiones algebraicas en términos de la densidad y en las que se utilizan una expresión parametrizada
156
UNED
Método de Kohn-Sham I
Máster FSC. Funcional de la densidad
para la energía XC del gas homogéneo y su derivada con respecto de la densidad. De esta manera la
resolución de las ecuaciones KS es abordable incluso en sistemas muy complejos. En concreto, la LDA se
utiliza hoy en día incluso para diseñar estructuras bioquímicas muy complejas y analizar su evolución
UNED
dinámica bajo la aproximación Born-Oppenheimer (dinámica molecular de primeros principios).
6. Como hemos comentado, la idea subyacente de la LDA es suponer que el sistema físico se comporta
localmente como un gas homogéneo. Por tanto, la LDA será una buena aproximación en tanto que el
per l de densidad del sistema electrónico bajo estudio sea razonablemente homogéneo. Sin embargo
esta a rmación debe matizarse, ya que es evidente que la densidad de un átomo o una molécula dista
mucho de parecerse a una densidad homogénea. Una manera sensata de cuanti car cuán inhomogéneo
UNED
es un per l de densidad en un punto ~r es de nir la siguiente cantidad, con dimensiones de momento
lineal
~ n (~r)
∇
κ (~r) = (6.45)
2 n (~r)
por lo que cuanto mayor sea κ (~r) mayor será la inhomogeneidad del sistema. Sin embargo esto no
basta para saber si el sistema es realmente inhomogeneo, en el sentido de que su comportamiento local
UNED
sea muy diferente al del gas de electrones uniforme. Para ello conviene de nir las longitudes λs y λF
asociadas a la densidad local:
λs (n (~r)) ' 0,504 × n (~r)−1/6 ; λF (n (~r)) ' 0,323 × n (~r)−1/3 (6.46)
y evaluar las cantidades adimensionales (funciones de Perdew)
~ n (~r)
∇ ~ n (~r)
∇
t (n (~r)) = λs (~r) κ (~r) ' 0,252 × s (n (~r)) = λF (~r) κ (~r) ' 0,162 ×
UNED
; (6.47)
n (~r)7/6 n (~r)4/3
Estas dos funciones nos permiten estimar si la aproximación local es adecuada para el tratamiento del
sistema. Cuanto menores sean s (n (~r)) y t (n (~r)) mejores son, en principio, las aproximaciones locales
al intercambio y a la correlacion, respectivamente. En todo caso este criterio ha de considerarse como
una mera estimación heurística, ya que los efectos cuánticos obedecen a leyes mucho más complejas.
En un sólido cristalino, si nos restringimos a los electrones de valencia, λs ' 1 y λF ∈ (1, 2) en el
UNED
rango de densidades típicas de dichos electrones de valencia. Por otro lado, puesto que la variación de
la densidad de los electrones de valencia es generalmente pequeña, κ (~r) va a ser sustancialmente menor
que la constante de red del sólido, que es del orden de unas pocas unidades atómicas. En de nitiva,
los parámetros t y s en un sólido cristalino típico van a ser, en el peor de los casos, del orden de la
unidad. Esto justi ca el mencionado éxito de la LDA en Física del estado sólido, sobre todo cuando
tratamos sistemas en los que los electrones de valencia son s o p. En materiales formados por elementos
de transición entran en juego electrones d, el per l de densidad va a ser mas inhomogéneo, y la LDA
UNED
puede empezar a fallar, pero casi nunca de forma dramática.
En un átomo o una molécula, siempre que no nos alejemos mucho de los núcleos, la función s va
a ser menor que 3 (típicamente) y la función t va a ser bastante menor. Las regiones en las que la
densidad es apreciable son las que más contribuyen a la energía del sistema, por lo que la LDA es
adecuada a la hora de calcular energías XC de un átomo. Sin embargo, si nos alejamos de los núcleos,
las densidades decaen exponencialmente y, en este caso, t y s divergen. Por tanto la LDA no va a ser
tan buena en un sistema localizado como en un sólido y, además, no es adecuada en absoluto a la hora
de tratar las regiones con densidad muy baja en un sistema localizado, lo que es esencial al analizar
el enlace químico. Además, dentro de una tónica general buena, la bondad de la LDA va a depender
sensiblemente de la geometría del sistema. Por tanto, el cálculo de propiedades estructurales (basado
en las diferencias de energía entre geometrías distintas) puede resultar cuantitativamente erróneo en
la LDA y, en los casos en los que las diferencias energéticas entre diferentes formas alotrópicas son
pequeñas, incluso cuantitativamente. Así, aunque las distancias de enlace en moléculas y sólidos están
157
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
afectadas por un error del orden del 1 % (por defecto) en la LDA (que resulta excelente, aunque en
consecuencia tiende a favorecer estructuras más compactas que en la realidad), las frecuencias de
vibración pueden tener errores del orden del 10 % o incluso más. Aun peor, las energías de formación
UNED
pueden verse afectadas por errores de hasta el 30 %, especialmente en sólidos donde ha de compararse la
energía por átomo del solido con la del átomo individual. Por último, las barreras energéticas asociadas
a reacciones químicas evaluadas bajo la LDA son típicamente la mitad que los valores exactos. En
cualquier caso, conviene insistir que estas limitaciones son realmente pequeñas si las comparamos con
su sencillez de implementación.
7. Una crítica más fundamental a la LDA está basada en el propio carácter local de la aproximación.
UNED
En sistemas físicos poco densos (en el sentido de que los núcleos están muy alejados entre sí) o
en sistemas formados por distintos subsistemas (un átomo o una molécula interaccionando con una
super cie, por ejemplo), la densidad electrónica tiende a localizarse en regiones diferenciadas. Los
efectos de correlación (sobre todo coulombiana) entre diferentes regiones pueden resultar importantes,
pero no van a ser captados en absoluto por la LDA ya que, recordemos, la energía XC por partícula
en un punto sólo depende de la densidad en ese punto, sin tener en cuenta la in uencia ni del entorno
ni de la densidad en otras regiones alejadas. Esto tiene importancia especial en la ya mencionada
UNED
evaluación de barreras de reacción y en el cálculo de la interacción entre sistemas neutros alejados
entre sí. En este último caso la fuerza de interacción es de tipo van der Waals, y este aspecto no va a
ser contemplado en absoluto por la LDA.
A su vez, siendo la LDA una buena aproximación a la energía de intercambio-correlación, no lo es
tanto si consideramos el intercambio y la correlación por separado. Esto quiere decir que la LDA se ve
bene ciada por una cancelación sistemática de errores entre la energía de intercambio y la correlación,
UNED
cancelación que se basa en el hecho de que ambos efectos se abordan desde la misma perspectiva
teórica. Esta cancelación es especialmente afortunada en sistemas laminares, como el gra to, y sirve
para enfatizar un hecho esencial en el diseño de aproximaciones funcionales a la energía XC: no es
en absoluto una buena idea el modelizar la energía de intercambio y la de correlación por separado
usando diferentes niveles de so sticación.
Esta dependencia directa en la densidad también explica otro error típico de la LDA: posee una
UNED
energía XC de autointeracción espúrea. En un sistema con un único electrón (por ejemplo, el átomo
de hidrógeno), la energía XC ha de ser exactamente igual a menos la energía de Hartree, puesto que
la energía de interacción total es cero. En un sistema de dos electrones con espín total igual a cero,
ya hemos visto que la energía de intercambio exacta ha de ser, salvo signo, la mitad de la energía de
Hartree. Estas correcciones de autointeracción no se reproducen exactamente en la LDA, por lo que
hay un error sistemático que es mas importante a medida que el número de electrones es pequeño y,
en general, en aquellos casos en los que las funciones de onda de los orbitales KS tiendan a localizarse.
UNED
A su vez, insistamos en que εXC (~r) y vXC (~r) dependen únicamente de la densidad en n (~r). Así,
en sistemas localizados, el comportamiento asintótico de estas cantidades va a ser exponencialmente
decreciente, re ejando el decaimiento que tiene la densidad. Sin embargo hemos visto que, en realidad,
el comportamiento asintótico ha de ser proporcional a −1/r. A efectos de energías este problema
no es relevante, ya que la energía XC por unidad de volumen va a ser n (~r) εXC (~r) ∼ 0 lejos del
centro del sistema. Sin embargo el potencial XC no va a exhibir este decaimiento y el potencial
efectivo Kohn-Sham vS (~r) bajo la LDA va a decaer exponencialmente en sistemas neutros. Esto se
traduce en que las energías de los orbitales KS bajo la LDA no van a ser buenas si se comparan
con los valores exactos, sobre todo a medida que nos acercamos al nivel HOMO del sistema cticio
KS. No cabe esperar entonces que en un sistema de N electrones, la energía del nivel KS N -ésimo,
εN , sea igual a menos la energía de ionización del sistema. En el contexto de la LDA es preferible
calcular dicha energía de ionización comparando las energías fundamentales del sistema con N
electrones y del primer ion positivo con N −1 electrones. Por otro lado, la ausencia de ese término
158
UNED
Método de Kohn-Sham I
Máster FSC. Funcional de la densidad
6
−1/r hace que no exista la serie Rydberg de orbitales KS no ocupados. Esto es irrelavente en el
calculo de la energía fundamental, ya que dichos estados no contribuyen a la densidad (que, no
nos olvidemos, es la variable fundamental al calcular las propiedades del estado fundamental). Sin
UNED
embargo la función respuesta del sistema KS, χ0 (~r1 , ~r2 ; ω), va a ser muy diferente si comparamos los
resultados LDA con los que se tendrían en un tratamiento hipotéticamente exacto. En χ0 (~r1 , ~r2 ; ω)
aparecen contempladas excitaciones electrón-hueco que están ausentes en la respuesta LDA debido a la
mencionada desaparición de orbitales Rydberg desocupados. Esta limitación es relevante en cualquier
cálculo de propiedades de estados excitados basado en el cálculo de la función respuesta interactuante
a partir de una función respuesta no interactuante. Peor aun es la situación en iones negativos, en los
N , pero hay N + 1 (o incluso más) electrones. En este caso el decaimiento
que la carga nuclear neta es
UNED
del potencial KS LDA es +1/r , puesto que la carga neta del sistema es +1 y es el potencial de Hartree
el que domina a largas distancias. Sin embargo, en el tratamiento exacto el comportamiento dominante
de vXC (~r) ha de compensar al potencial de Hartree si r 0. Esto se traduce en que, en la LDA, el
potencial efectivo Kohn-Sham no es capaz de ligar N + 1 estados, por lo que en la LDA muchos iones
negativos son, erróneamente, inestables.
Por último, en la LDA la dependencia de la densidad es explícita y analítica. El funcional XC no
exhibe las discontinuidades que debería poseer el funcional XC exacto cuando consideremos pequeñas
UNED
variaciones del número total de electrones. Esto puede parecer irrelevante, pero tiene su importancia a
la hora de interpretar las energías de los orbitales KS en sistemas extensos. Además, pone de mani esto
cuan difícil es, a priori, modelizar el funcional de la energía XC que es, en cierta medida, patológico.
6.3.c. La aproximación local dependiente del espín (LSDA)
8. En la extensión de la DFT a sistemas con dependencia en el espín, la aproximación LDA al
UNED
funcional XC es
Z
LSDA [n , n ] =
EXC d3~r n (~r) εhom
+ − XC (n (~
r) , ζ (~r)) (6.48)
donde ζ (~r) = [n+ (~r) − n− (~r)] /n (~r) es la polarización local en espín, n (~r) = n+ (~r) + n− (~r) y
εhom
XC (n, ζ) es la energía XC por partícula para un gas homogéneo con densidad total n y polarización
en espines ζ . Puesto que el estado fundamental del gas homogéneo tiene una polarización en espines
hom
bien de nida, εXC (n, ζ) se de ne como
UNED
1 D b cE 1
εhom
XC (n, ζ) = T +W − tS (n, ζ) − WH [n] (6.49)
N n,ζ N
D E
donde Tb + W
c /N es el mínimo del valor esperado de Tb + W
c por partícula sobre todos los estados
n,ζ
cuánticos del gas homogéneo cuya densidad sea n y su polarización ζ, mientras que tS (n, ζ) es la
energía cinética por partícula del gas no interactuante con densidad n y polarización ζ. Al igual que
UNED
sucedía en la versión independiente del espín, la componente de intercambio es conocida, mientras que
la de correlación se puede obtener combinando teoría avanzada y simulaciones numéricas, estando esta
última convenientemente parametrizada. Es operativo comprobar que, en la LSDA
(±) ∂εhom
XC (n (~r) , ζ (~r)) ∂εhom (n (~r) , ζ (~r))
vXC (~r) = εhom
XC (n (~
r) , ζ (~r)) + n (~r) ± [1 ∓ ζ (~r)] XC (6.50)
∂ n (~r) ∂ ζ (~r)
9. Para sistemas en los que el estado fundamental no exhibe polarización en espines, la LDA y
la LSDA dan idénticos resultados. Sin embargo, en sistemas en los que el estado fundamental posee
polarización en espines hay diferencias sustanciales. Sirva como ejemplo el cálculo de la energía de
6
Esta es la serie de in nitos estados ligados por el decaimiento −1/r del potencial.
159
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
N = 2 y cuando N = 1.
ionización del helio a partir de la diferencia de energías del sistema cuando
La energía fundamental del átomo neutro es −2,83 tanto en la LDA como en la LSDA (compárese
con el valor exacto −2,90). Sin embargo, cuando N = 1 la energía LSDA es −1,93 mientras que la
UNED
energía LDA es −1,86 (el valor exacto es el resultado hidrogenoide −2). En este caso, la LSDA da
como resultado una densidad totalmente polarizada en espín, mientras que la LDA evalua la energía
XC considerando simplemente la densidad no polarizada en espin. Así, de acuerdo con la LSDA, la
energía de ionización del Helio es 0.90, según la LDA es 0.97 y el resultado exacto es 0.90. Sirva de
comentario que si hubiésemos calculado la energía de ionización a partir de la energía del nivel HOMO
en la LSDA, el resultado hubiese sido 0.57. En cualquier caso, las ventajas y limitaciones de la LSDA
son las mismas que las de la LDA salvo las inevitables diferencias cuantitativas.
UNED
6.3.d. Las aproximaciones de expansiones en gradientes (GGA)
10. La mejora más evidente de la aproximación LDA es la inclusión explícita de las funciones
dependientes de la densidad
~ n (~r)
∇ ~ n (~r)
∇
UNED
t (~r) = 0,251 93 × ; s (~r) = 0,161 62 ×
n (~r)7/6 n (~r)4/3
para así intentar corregir la energía XC bajo la LDA en aquellas regiones en las que la hipótesis
de comportamiento similar al HEG deja de ser válida. Esta es la idea central de las llamadas
aproximaciones de gradientes generalizadas (GGA). Sin entrar en excesivos detalles (sólo hay una
LDA, pero hay docenas de GGA's), la energía XC se aproxima mediante la forma funcional
UNED
GGA [n] d3~r n (~r) εhom
EXC [n] ' EXC = XC (n (~
r)) FGGA (n (~r) , t (~r) , s (~r)) (6.51)
siendo FGGA una función que hay que determinar (naturalmente se puede -y se debe, por los motivos
indicados en la sección anterior- de nir la GGA bajo el formalismo KS dependiente del espín). Así, en
la GGA la energía XC por partícula en el punto ~r es igual al valor correspondiente al gas homogéneo
con densidad n (~r) multiplicado por el factor de corrección FGGA que depende de la densidad y
de las funciones t (~r) y s (~r) a través de las cuales se pretenden modelizar efectos de correlación e
UNED
intercambio más allá de la aproximación local, respectivamente. Puesto que FGGA será una función
analítica, aunque más complicada, la derivada variacional de EXC [n] se podrá expresar explícitamente
en términos de la densidad, su gradiente y su laplaciana.
11. Las GGA's se pueden clasi car en dos grandes grupos: funcionales GGA de primeros principios,
en cuya construcción se utilizan exclusivamente propiedades generales del funcional XC y del gas
homogéneo; y funcionales GGA semiempíricos, en los que se reproducen algunas propiedades exactas
UNED
pero, sobre todo, la función FGGA se modeliza de manera que las energías XC para algunos sistemas
seleccionados (átomos y moléculas, generalmente) sean muy parecidas a las exactas, obtenidas
mediante técnicas de Química Cuántica.
El ejemplo más característico de GGA de primeros principios es el funcional de Perdew, Burke
y Ernzerhof (PBE) [Phys. Rev. Lett. 77, 3865 (1996)]. Está diseñado de manera que se cumplen
exactamente algunas propiedades de escalamiento del funcional EXC (que pueden encontrarse descritas
en el capítulo 1 del libro editado por Fiolhais et al.) y, sobre todo, del gas homogéno incluyendo
propiedades de respuesta estática en el límite de largas longitudes de onda. Los detalles del desarrollo
y la expresión explícita (complicada pero analítica) del funcional PBE pueden encontrarse en el texto
citado.
En cuanto a los funcionales GGA empíricos, el más popular es el denominado BLPY (puesto
que es igual a un funcional GGA para el intercambio desarrollado por Becke [Phys. Rev. A 38, 3098
160
UNED
Método de Kohn-Sham I
Máster FSC. Funcional de la densidad
(1988)] construido empíricamente para reproducir en buena medida las energías de intercambio exactas
de algunos átomos y moléculas; LPY viene de Lee, Parr y Yang, quienes desarrollaron una forma
funcional GGA para la correlación [Phys. Rev. B 37, 785 (1988)] a partir de resultados teóricos de
UNED
química cuántica e imponiendo la reproducción exacta de la energía de correlación para el átomo
de He). Puesto que el funcional BLPY está, por construcción, optimizado para sistemas atómicos y
moleculares, sus resultados son realmente buenos en este ámbito. En cierta medida, el funcional BLPY
es el responsable de que la DFT se halla convertido en una herramienta popular en Química Cuántica,
venciendo las reticencias tradicionales que los químicos teóricos mostraban a la DFT. Sin embargo,
estos tipos de funcionales optimizados para sistemas localizados son peores que la LDA a medida que
tratamos sistemas más extensos. De hecho, el funcional BLPY ni siquiera reproduce la energía XC
UNED
exacta del gas homogéneo (!).
12. Puede decirse que, hoy en día, las GGAs (incluyendo la dependencia en espín) son el estándar
en los cálculos DFT, sobre todo en sistemas atómicos y moleculares. Cuantitativamente son mejores
que la LDA aunque hay excepciones notables, principalmente en sistemas extensos. Por ejemplo, bajo
la GGA-PBE el gra to es inestable y, en general, todos los solidos laminares son muy mal descritos
por las GGA's. Sin embargo las GGA's no solventan las limitaciones fundamentales de la LDA . Bajo
UNED
esta perspectiva, las GGA no dejan de ser meras extensiones semilocales de la LDA. La construcción
de modelos funcionales para la energía XC más precisos requiere un cambio conceptual en su diseño,
que introduzca de manera explícita la naturaleza no-local de los efectos de intercambio-correlación.
Trataremos someramente estas aproximaciones en el siguiente capítulo.
UNED
UNED
UNED 161
Máster FSC. Funcional de la densidad
UNED
UNED
UNED
UNED
UNED
UNED
UNED
UNED
Máster FSC. Funcional de la densidad
Capítulo 7
MÉTODO DE KOHN-SHAM II
UNED
OBJETIVOS:
1. Discusión de aproximaciones más avanzadas a los funcionales de la energía XC.
UNED
2. Análisis del problema del gap
3. Descripción de la implementación del método Kohn-Sham para sistemas periódicos.
LECTURAS COMPLEMENTARIAS:
UNED
1. A primer in Density Functional Theory [Fiolhais et al ( Eds.)]: Capítulos 5 y 6
2. Fiolhais et al ( Eds.): Capítulo 2 (avanzado).
ORIENTACIONES AL PLAN DE TRABAJO:
UNED
En la primera parte de este capítulo estudiaremos brevemente otras aproximaciones a los
funcionales de la energía de intercambio y correlación que, recuerde, es la única parte del funcional
de la energía total que no está descrito a priori de manera exacta en el método de Kohn-Sham.
A continuación analizaremos brevemente el llamado problema del gap, de gran importancia en
la interpretación de los resultados Kohn-Sham. Finalmente describiremos brevemente como se
implementa el método de Kohn-Sham a sistemas periódicos, aspecto que esta relacionado con
las actividades prácticas que ilustran este tema.
UNED
A título orientativo, se espera que dedique 10 horas a este capítulo desglosadas como sigue:
• 6 horas de lectura detallada y estudio del presente material didáctico.
• 4 horas para el a anzamiento de los contenidos, incluyendo consultas y las lecturas
complementarias [optativas].
En este capítulo no se proponen problemas.
UNED
GLOSARIO:
Aproximación WDA
Función de correlación de pares
Aproximación ADA
Función de respuesta lineal
Aproximación meta-GGA
Funcionales XC híbridos
163
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
Funcionales RPA+
Implementación exacta del intercambio
UNED
Aproximación KLI al potencial de intercambio
Problema del gap
Ondas planas
Teorema de Bloch
UNED
Pseudpotenciales
UNED
UNED
UNED
UNED
164
UNED
Método de Kohn-Sham II
Máster FSC. Funcional de la densidad
7.1. APROXIMACIONES NO LOCALES AL FUNCIONAL XC
7.1.a. La aproximación no local WDA
UNED
1. Como ya sabemos, el funcional de la energía XC puede escribirse como
Z
EXC [n] = d3~r n (~r) εXC [n] (~r)
donde εXC [n] (~r) es la energía XC por partícula en el punto ~r. Las aproximaciones LDA y GGA
modelizan directamente esta cantidad a través de la expresión general
UNED
εXC [n] (~r) ' εhom
XC (n (~
r)) FGGA (n (~r) , t (~r) , s (~r))
donde, recordemos,
t (~r) = 0,251 93 × n (~r)−7/6 ∇
~ n (~r)
s (~r) = 0,16162 × n (~r)−4/3 ∇
~ n (~r)
UNED
La aproximación WDA (Weighted Density Approximation, aproximación de densidad ponderada)
abandona esta losofía, modelizando el funcional XC a partir de un nivel de teoría más profundo. El
punto de partida es la expresión exacta en términos del hueco XC centrado en ~r, η XC (~r ; ~r 0 ):
r ; ~r 0 ) r 0 ) hXC (~r, ~r 0 )
Z Z Z Z
3 1 3 0 η XC (~ 3 1 3 0 n (~
EXC [n] = d ~r n (~r) d ~r = d ~r n (~r ) d ~r
2 |~r − ~r 0 | 2 |~r − ~r 0 |
y modelizar la función de correlación de pares integrada hXC (~r, ~r 0 ). En concreto, la idea fundamental
UNED
de la WDA es suponer que la función de correlación de pares integrada de un sistema arbitrario
puede aproximarse por la del gas homogéneo convenientemente reescalada. Especí camente, si
hhom r − ~r 0 |) es la fcp integrada de un gas homogéneo con densidad n, la hipótesis WDA consiste
XC (n, |~
en aproximar la función hueco XC centrado en ~r de la siguiente manera:
η XC (~r ; ~r 0 ) = n (~r 0 ) hXC (~r, ~r 0 ) ' n (~r 0 ) hhom r) , |~r − ~r 0 |)
XC (ñ (~ (7.1)
UNED
donde ñ (~r) es una densidad ponderada (weighted density) que, para cada punto ~r, se obtiene
imponiendo la condición de normalización del hueco XC:
Z Z
3 0 0
d3~r 0 n (~r 0 ) hhom r) , |~r − ~r 0 |) = −1
d ~r η XC ~r ; ~r = −1 ⇒ XC (ñ (~ ∀ ~r (7.2)
La función hhom
XC (n, |~r − ~r 0 |) es conocida, existiendo expresiones parametrizadas y, de esta manera,
para cada punto ~r resolveríamos la ecuación anterior. Calculada la densidad ponderada, la energía XC
UNED
por partícula en el punto ~r será
n (~r 0 ) hhom r) , |~r − ~r 0 |)
Z
1 XC (ñ (~
εXC [n] (~r) ' d3~r 0 (7.3)
2 |~r − ~r 0 |
2. Técnicamente, la WDA es mucho más costosa y delicada que la LDA/GGA. En primer lugar, para
cada punto ~r tenemos que obtener la densidad ponderada ñ (~r), lo que no es difícil pero cuyo esfuerzo
computacional crece sensiblemente con el tamaño y complejidad del sistema. Por otra parte, tanto
la evaluación de εXC [n] (~r) como de vXC (~r) implican integrales sobre todo el espacio que ralentizan
considerablemente el cálculo. Ciertamente, la WDA es una alternativa válida a la LDA/GGA sólo si
la ganancia en precisión compensa el esfuerzo computacional subyacente.
165
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
3. Desde un punto de vista fundamental, dos son las ventajas de la aproximación WDA respecto
UNED
de las LDA/GGA. En primer lugar, los errores de autointeracción son muy pequeños, lo que la hace
especialmente atractiva en el estudio de metales de transición. En segundo lugar, puede demostrarse
que el potencial XC exhibe el comportamiento asintótico correcto en sistemas nitos:
1
lı́m v WDA (~r) =−
r→∞ XC r
Ello implica que, por lo general, la energía de ionización calculada a partir del nivel HOMO del sistema
UNED
cticio Kohn-Sham sea mejor que la que se obtendría a partir de los modelos LDA/GGA. A su vez, la
WDA liga iones negativos, solventándose así una de las limitaciones características de las LDA/GGA.
Por último, en sistemas muy inhomogéneos (sistemas cuasi-bidimensionales, por ejemplo) la WDA es
sustancialmente mejor que la LDA/GGA.
Sin embargo, en aquellos sistemas en los que las LDA/GGA funcionan razonablemente, la WDA
no suele exhibir mejoras sustanciales por lo que no puede considerarse un avance fundamental hacia
la búsqueda de un funcional XC con precisión química e, incluso, en sistemas de referencia sencillos
UNED
(átomos y moléculas) los resultados estructurales pueden ser algo peores que los LDA/GGA. Además,
la WDA modeliza el hueco XC a partir del correspondiente a un metal (el gas homogéneo), por lo
que su uso no estaría muy justi cado en sistemas extensos que sean semiconductores o aislantes. De
hecho, la gran limitación de la WDA (además de su coste computacional) es el hecho de que los
resultados numéricos dependen sensiblemente de la parametrización que se utilice para describir la
fcp del gas homogéneo hhom r − ~r 0 |).
XC (n, |~ En general, podemos decir que mientras que el potencial XC
bajo la WDA es mejor que el LDA/GGA, las energías XC dadas por la WDA no son sustancialmente
UNED
mejores que las LDA/GGA, por lo que no cabe esperarse una mejor descripción de las propiedades
estructurales. De hecho, el interés de los modelos WDA van dirigidos a la búsqueda de buenos
potenciales XC, hecho que es importante a la hora del cálculo de energías de excitación neutras
dentro del ámbito de la llamada teoría del funcional de la densidad dependiente del tiempo (TDDFT).
7.1.b. La aproximación no local ADA
UNED
4. Otra alternativa es la llamada aproximación de densidad promediada (Averaged Density
Approximation). En este modelo no se pretende modelizar el hueco XC del sistema físico, sino
reproducir las propiedades de respuesta lineal del gas homogéneo. Como ya hemos visto, una de
las propiedades fundamentales del funcional XC es la siguiente:
δ 2 EXC [n] 1
= χ−1 r1 , ~r2 ) −
0 (~ − χ−1 (~r1 , ~r2 )
δn (~r1 ) δn (~r2 ) n=n0 |~r1 − ~r2 |
UNED
donde χ0 y χ son las funciones respuesta del sistema cticio KS y del sistema real, respectivamente.
Partiendo de esta idea, el modelo ADA propone aproximar la energía XC por partícula en el punto ~r
como
εXC [n] (~r) ' εhom
XC (n (~
r)) (7.4)
donde n (~r) es una densidad promediada dada por la ecuación integral
Z
d3~r 0 n ~r 0 Ω n (~r) , ~r − ~r 0
n (~r) = (7.5)
La función peso Ω (n (~r) , |~r − ~r 0 |) está normalizada,
Z
d3~r 0 Ω n (~r) , ~r − ~r 0
= +1 (7.6)
166
UNED
Método de Kohn-Sham II
Máster FSC. Funcional de la densidad
y se determina imponiendo que la segunda derivada variacional del funcional XC en el límite homogéneo
reproduzca la forma exacta (conocida a partir de parametrizaciones).
UNED
5. Desde el punto de vista computacional, el modelo ADA es igual de costoso que el WDA. Tiene el
atractivo, respecto del modelo WDA, de reproducir una propiedad energética (la respuesta lineal en el
límite homogéneo) en lugar de modelizar la función de correlación de pares integrada. Los errores de
autointeracción son pequeños y el comportamiento asintótico del potencial XC en sistemas localizados
exhibe un decaimiento potencial:
c
lı́m v ADA (~r) =−
UNED
r→∞ XC r
aunque el factor c no es exactamente la unidad (en todo caso es cercano a 1). Sin embargo, quizás
por el hecho de que computacionalmente son costosos, los modelos ADA no han sido explorados
sistemáticamente en la literatura.
6. En cualquier caso, tanto los modelos ADA como los WDA estan construidos a partir de una
UNED
dependencia explícita en la densidad, por lo que no exhiben ningún tipo de comportamiento no analítico
característico del funcional XC exacto. Además, a pesar de su no localidad (la energía XC por partícula
en el punto ~r depende de la densidad en todos los puntos del sistema) la descripción que hacen de
fuerzas de largo alcance tipo van der Waals no es correcta siquiera cualitativamente.
7. Terminemos esta sección comentado que los modelos ADA y WDA, aunque no se han explorado
su cientemente a la hora de modelizar la energía XC si se han estudiado con detalle en la construcción
UNED
de aproximaciones explícitas de la densidad al funcional de Hohenberg y Kohn FHK [n].
7.2. FUNCIONALES XC CON DEPENDENCIA EN LOS ORBITALES KS
1. Los funcionales LDA, GGA, WDA y ADA para la energía XC dependen de manera explícita en
la densidad electrónica n (~r). Sin embargo, como ya hemos comentado repetidas veces, el funcional XC
UNED
exacto posee una dependencia no trivial en la densidad. Una alternativa completamente diferente es la
modelización del funcional XC incluyendo los propios orbitales KS y, eventualmente, sus autoenergías.
Formalmente es una losofía válida, puesto que el sistema cticio KS depende "funcionalmente"de la
densidad electrónica a través de las propias ecuaciones KS. Estos funcionales .avanzadosçonstituyen
un campo de investigación muy activo (véase al respecto el reciente artículo de revisión publicado por
Kümmel y Kronik [Rev. Mod. Phys. 80, 3 (2008)]), aunque en estas notas nos limitaremos a un breve
resumen, nada exhaustivo.
UNED
7.2.a. Algunos ejemplos
2. Un primer ejemplo de estos funcionales son los denominados meta-GGA, propuestos en los
últimos años por Perdew y colaboradores. Sin entrar en los detalles especí cos de su construcción,
estos funcionales XC pueden considerarse una forma más so sticada de los modelos GGA en la que la
energía XC por partícula tiene la forma
εXC [n] (~r) ' εhom
XC (n (~
r)) FmGGA (n (~r) , t (~r) , s (~r) , τ (~r)) (7.7)
donde ahora la función correctora FmGGA incluye un término, τ (~r), que es igual a la densidad de
energía cinética KS:
N
1 XX ~
τ (~r) = ∇φi (~r σ)
2 σ
i=1
167
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
De este modo, el funcional depende de los orbitales KS ocupados φi (~r). Generalmente, la idea central
que guía estos modelos es la eliminación casi completa de los errores de autointeracción. Este tipo
de funcionales dan excelentes resultados en sistemas moleculares (incluyendo macromoléculas) y
UNED
representan un excelente compromiso entre abilidad y coste computacional.
3. Una losofía diferente se basa en el hecho de que en la descomposición EXC [n] = EX [n] + EC [n]
de la energía XC, la primera contribución es conocida:
occ ∗
φj (~r 1 σ) φj (~r2 σ) φi (~r1 σ) φ∗i (~r2 σ)
Z
ex 1 3 3
XX
UNED
EX [n] =− d ~r1 d ~r2 (7.8)
2 σ
|~r2 − ~r1 |
i,j
(siguiendo una costumbre estéticamente infame, indicamos explícitamente que el funcional es el
exacto). Vemos que la dependencia funcional de
ex [n]
EX está contenida enteramente en los orbitales
ocupados del sistema cticio KS asociado a la densidad n. Así, resulta sugerente usar la expresión
exacta del funcional de intercambio y aproximar la energía de correlación:
ex app
UNED
EXC [n] ' EX [n] + EC [n] (7.9)
Bajo esta perspectiva, es el funcional de correlación el único que necesita modelizarse, y podemos usar
prescripciones tipo LDA, GGA o incluso metaGGA para este n. De este modo tendríamos funcionales
para la energía XC que podemos notar como EXX/LDA, EXX/GGA y EXX/mGGA, respectivamente.
1
Aunque la energía de intercambio exacta (7.8) es poco costosa de evaluar, no sucede lo mismo
con su derivada variacional ya que, como veremos en la siguiente subsección, su cálculo exige resolver
UNED
una ecuación integral. Sin embargo, su cálculo es factible y en la actualidad este tratamiento exacto
del intercambio en la DFT está implementado en muchos códigos de química cuántica. Más aún, dicha
derivada variacional se puede evaluar de manera aproximada y computacionalmente muy e ciente pero
sin perder mucha precisión en el cálculo, lo que aumenta el campo de aplicabilidad de este tipo de
funcionales.
Sin embargo, aunque atractiva a primera vista, la estructura funcional (7.9) es muy arriesgada:
los funcionales LDA y GGA para la energía XC se bene cian de cancelaciones de errores entre el
UNED
intercambio y la correlación. Por tanto, si evaluamos exactamente el intercambio deberíamos aproximar
la energía de correlación de manera muy so sticada, lo que no es posible bajo los modelos LDA y GGA.
Como ilustración de esta limitación, discutamos los resultados para la distancia internuclear y de la
energía de atomización en la molécula de nitrógeno, cuyos valores exactos son, respectivamente, 2.075
Bohr y 229 kcal/mol. La LDA da como resultado 2.068 Bohr y 266 kcal/mol, mientras que la GGA-PBE
resulta en 2.079 Bohr y 242 kcal/mol. Observamos que la con guración de equilibrio es excelente en
las LDA y GGA, pero la LDA sobreliga (overbinds) los átomos de nitrógeno, tendencia que sigue
UNED
la GGA aunque el error es aproximadamente la mitad que el de la LDA. Si abordamos el problema
despreciando la correlación y tratando exáctamente el intercambio, los valores que se obtienen son
2.011 Bohr y 112 kcal/mol. La geometría de equilibrio es bastante peor pero lo mas destacable es que,
aunque el valor de las energías de correlación sea pequeño en comparación con las de intercambio,
su papel a la hora de evaluar energías de atomización es muy relevante: la correlación es responsable
del 50 % de la energía de atomización. Ciertamente estos resultados pueden mejorarse incluyendo la
correlación, pero si describimos ésta usando una GGA (con una LDA los resultados son similares), la
distancia de equilibrio es 1.997 Bohr (½incluso peor que despreciando el intercambio!) y la energía de
atomización de 171 kcal/mol. En general, las propiedades estructurales de sistemas moleculares bajo
el tratamiento exacto del intercambio junto con una aproximación sencilla a la correlación (LDA o
GGA) son muy imprecisas y, de hecho, el tratamiento EXX/GGA resulta en que la molécula de uor es
1
Sin embargo hay detalles técnicos sutiles que siempre hay que tener en cuenta.
168
UNED
Método de Kohn-Sham II
Máster FSC. Funcional de la densidad
inestable (?). Algo similar ocurre en sólidos laminares, como el gra to y el nitruro de boro en su fase
alotrópica hexagonal, que son inestables bajo la EXX/LDA y la EXX/GGA (estas aproximaciones
son incapaces de predecir la existencia de una distancia de equilibrio mecánico entre las láminas).
UNED
Por el contrario, la LDA (pero no la GGA) da resultados bastante razonables para las geometrías de
equilibrio de estos sistemas.
Sin embargo, la superioridad de los modelos EXX/LDA y EXX/GGA es mani esta si calculamos
energías de ionización a partir del nivel HOMO del sistema cticio KS y si tenemos en cuenta que
describen correctamente la existencia de iones negativos. En estos dos casos, la razón fundamental
del éxito de estos funcionales es la corrección (al menos a nivel de intercambio) de los errores de
autointeracción, corrección que conduce a la reproducción del comportamiento asintótico −1/r del
UNED
potencial XC. De una forma cualitativa, puesto que los modelos LDA y GGA poseen errores de
autointeracción, el electrón del nivel HOMO del sistema cticio KS está sometido a un potencial
efectivo que contiene un término repulsivo espúreo, por lo que su energía es mayor (menos negativa)
que lo que correspondería al tratamiento KS exacto. Esto explica por qué la LDA y la GGA subestiman
las energías de ionización calculadas a partir de la energía del nivel HOMO. En el caso EXX (sin
correlación) no hay errores de autointeracción, el electrón del nivel HOMO no sufre el potencial
repulsivo espúreo y su energía está más cercana al valor exacto. Si tratásemos la correlación de
UNED
forma muy so sticada, la energía del nivel HOMO sería entonces muy similar al valor exacto. Sin
embargo, si tratamos la correlación bajo la LDA o la GGA estamos introduciendo de nuevo efectos de
autointeracción que compensan el bene cio de incluir efectos de correlación aun aproximadamente: las
energías de ionización EXX/LDA y EXX/GGA son algo peores que a nivel de intercambio solamente,
aunque mucho mejores que las LDA o GGA. Como ilustración, para la molécula de N2 tenemos que (en
eV) I EXX = 15,5, I EXX/LDA = 16,9, I EXX/GGA = 16,4, I LDA ' I GGA = 10,3. El valor exacto es 15.6
eV, prácticamente igual al que se obtiene con el tratamiento exacto del intercambio pero excluyendo
UNED
la correlación lo que, en cierta medida, con rma los problemas de este tipo de modelos apuntados en
los dos párrafos anteriores.
4. Los funcionales híbridos, que han alcanzado una enorme popularidad en Química Cuántica, poseen
la estructura
ex app app
EXC [n] ' αEX [n] + (1 − α) EX [n] + EC [n] (7.10)
UNED
donde α es un parámetro (dependiente del modelo funcional) dentro del intervalo (0, 1). Así, en estos
funcionales híbridos se introduce el intercambio exacto pero las cancelaciones entre intercambio y
correlación se pretenden modelizar quitando peso a
ex [n]. El funcional B3LYP, en el que los términos
EX
app app
EX [n] y EC [n] se describen usando funcionales GGA empíricos, resulta muy preciso en un número
muy amplio de aplicaciones de química cuántica.
UNED
5. Aumentando el nivel de so sticación,podemos intentar evaluar el funcional XC de manera exacta
en la medida de lo posible, intentando reducir al máximo las aproximaciones funcionales. Así, la energía
de intercambio se podría obtener a partir de (7.8), mientras que la energía de correlación es evaluable
bajo la RPA:
Z +∞ Z 1
dω h
(λ)
i
Z Im dλ χRPA (~r1 , ~r2 ; ω) − χ0 (~r1 , ~r2 ; ω)
RPA 1 π
EC [n] ' EC [n] = − d3~r1 d3~r2 0 0
(7.11)
2 |~r2 − ~r1 |
donde, en notación matricial,
h i−1
(λ)
χRPA (ω) = bI − λb
χ0 (ω) w
b χ
b 0 (ω)
169
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
es la función respuesta RPA del sistema cticio de partículas que interactuan a través de un potencial
λw
b. Así, si de nimos ∆EC [n] como el error cometido al evaluar la energía de correlación usando la
RPA, tendríamos que
UNED
ex RPA
EXC [n] = EX [n] + EC [n] + ∆EC [n] . (7.12)
Así, el único término desconocido es ∆EC [n] que, por ejemplo, puede modelizarse à la LDA teniendo
en cuenta que la energía RPA del gas homogéneo por partícula, εRPA
C (n), se puede calcular de manera
relativamente simple:
Z h i
∆EC [n] ' d3~r n (~r) εhom
C (n (~r)) − εRPA
C (n (~r))
UNED
Este esquema es conocido en la literatura como EXX/RPA(+) (intercambio exacto + correlación
RPA + término corrector). El primer problema es que la evaluación de
RPA [n]
EC exige, en primer
lugar, calcular la función respuesta χ0 (~r1 , ~r2 ; ω) para diferentes frecuencias y, a continuación, para
esas frecuencias y diferentes valores de λ hay que evaluar χ b (λ) (ω), lo que precisa la inversión de
una matriz. Puesto que en la función respuesta χ0 (~ r1 , ~r2 ; ω) entran todos los orbitales ocupados y
desocupados del sistema cticio KS, un parámetro de convergencia fundamental es el cuto en estados
UNED
no ocupados lo que da lugar a que la implementación técnica sea delicada y muy costosa. Sin embargo,
la principal di cultad radica en la evaluación del potencial vC (~r), muchísimo mas complicada que la
asociada al potencial de intercambio exacto. Es por ello por lo que en muchos casos, el cálculo de
la energía XC se haga de manera no autoconsistente, considerando que el potencial vXC (~r) se puede
aproximar razonablemente mediante algunos de los métodos apuntados en las anteriores secciones.
En todo caso, es interesante poder evaluar vXC (~r) a partir de las expresiones que estamos dando
aquí para mantener la consistencia funcional del método y porque es importante obtener una buena
UNED
representación del potencial XC si queremos estudiar propiedades de excitación usando TDDFT. Esto
último es especialmente relevante en sistemas localizados.
El funcional (7.12) es, hasta la fecha, la implementación más precisa del método KS pero el número
de situaciones a las que se ha aplicado es limitado. En todo caso es un esquema prometedor pero que
necesita de un trabajo importante tanto en aspectos técnicos relativos a su programación como en
aspectos teóricos referentes a la modelización del único término desconocido, ∆EC [n].
UNED
6. Otra aproximación muy relacionada con la anterior es evaluar EC [n] usando teoría de
perturbaciones tipo MP2, pero con unas ligeras modi caciones necesarias para su inclusión formal
como implementación dentro la teoría del funcional de la densidad. Esta aproximación funcional está
discutida con bastante detalle en el segundo capítulo del texto de Fiolhais et al y también (más
super cialmente) en el review de Kümmel y Kronik. Siendo interestante esta forma de proceder, uno
no puede esperar que este método funcional funcione en aquellas situaciones en las que la teoría de
UNED
perturbaciones estándar a segundo orden de la química cuántica ofrezca resultados limitados. A su
vez, trabajos recientes han demostrado que existen problemas fundamentales relativos al potencial de
correlación asociado.
7. En general, la inclusión de los orbitales KS en el funcional de intercambio-correlación permite
obtener aproximaciones mucho más precisas que la LDA o la GGA, pero el precio que hay que pagar
es el mayor coste en la evaluación del propio funcional y, sobre todo, de la derivada variacional. Así se
pierde el principal atractivo de la LDA/GGA: su simplicidad y facilidad computacional. En todo caso,
como ya hemos anticipado, este tipo de funcionales se están abriendo camino cada vez más gracias
a la mejora de los códigos de cálculo y al aumento de la capacidad computacional accesible a los
investigadores.
170
UNED
Método de Kohn-Sham II
Máster FSC. Funcional de la densidad
7.2.b. Derivadas variacionales
8. La implementación completa de los modelos descritos en la subsección anterior exige el cálculo de
UNED
las derivadas variacionales de funcionales que dependen de los orbitales KS y sus autoenergías. Veamos
en esta subsección cuál es el método general para un funcional genérico B con dicha dependencia:
B [n] = B [{φn } , {εn }]
P PN
(nótese que n (~r) = σ r σ)|2 , por lo que también estamos considerando la situación
n=1 |φn (~
habitual en la que el funcional depende explícitamente de la densidad). Para simpli car, supongamos
que N es par y que los orbitales KS exhiben compensación de espines, esto es, que son de la forma
UNED
{φ1 (~r) α (σ) , φ1 (~r) β (σ) , φ2 (~r) α (σ) , φ2 (~r) β (σ) , ...}
La idea esencial es aplicar la regla de la cadena:
Z
δB [n] δB [n] δvS (~r1 )
= d3~r1
δn (~r) δvS (~r1 ) δn (~r)
donde vS (~r) es el potencial efectivo Kohn-Sham, y evaluar δB [n] /δvS (~r1 ) de la siguiente forma:
UNED
δB δφ∗m (~r2 )
XZ X
δB [n] δB δφm (~r2 ) ∂B δεm
= d3~r2 + ∗ +
δvS (~r) m
δφm (~r2 ) δvS (~r1 ) δφm (~r2 ) δvS (~r1 ) m
∂ε m δvS (~
r1 )
Las derivadas parciales respecto de las funciones de onda y de las autoenergías son accesibles. En
cuanto al resto de derivadas variacionales, la primera es conocida de acuerdo con la teoría de la
respuesta lineal:
δvS (~r1 )
= χ−1
0 (~
r, ~r1 )
UNED
δn (~r)
Considerando variaciones en primer orden de los orbitales bajo cambios in nitesimales del potencial
exterior en un hamiltoniano monopartícula, la segunda derivada variacional, δφm (~r2 ) /δvS (~r1 ) es
δφm (~r2 )
= −φm (~r1 ) Gm (~r2 , ~r1 )
δvS (~r1 )
donde
X φ (~r2 ) φ∗ (~r1 )
i i
= G∗m (~r1 , ~r2 )
UNED
Gm (~r2 , ~r1 ) = (7.13)
εi − εm
i6=m
Por último, aplicando igualmente teoría de perturbaciones
δεm
= |φm (~r1 )|2
δvS (~r1 )
En de nitiva Z
δB [n]
d3~r1 χ−1
UNED
= 0 (~
r, ~r1 ) ΛB (~r1 ) (7.14)
δn (~r)
con
X 2 ∂B
Z
δB
3
ΛB (~r1 ) = |φm (~r1 )| − d ~r2 φm (~r1 ) Gm (~r2 , ~r1 ) + c.c. (7.15)
m
∂ εm δφm (~r2 )
Así, la evaluación de δB [n] /δn (~r) precisa obtener las funciones Gm (~r2 , ~r1 ) para todos los orbitales
2
(ocupados o no), la función respuesta no interctuante , su inversa y el vector ΛB (~ r1 ).
2
Se puede usar el cálculo de las G ya que es fácil probar que
N/2
X
φ∗j (~ r1 ) φ∗j (~
χ0 (~ r2 ) = −2 ×
r1 , ~ r1 ) φj (~
r2 ) Gj (~
r1 , ~
r2 ) + φj (~ r2 ) Gj (~
r2 , ~
r1 )
j=1
171
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
9. En el caso especí co de EXex [n], sólo tenemos que evaluar δEXex [n] /δφm (~r2 σ) con m = 1, ..., N/2
UNED
ya que para el caso con compensación de espines,
N/2 N/2 ∗
1
Z X X φj (~r) φj (~r 0) φi (~r) φ∗i (~r 0)
ex 3 3 0
EX [n] = −2 × d ~r d ~r
2 |~r 0 − ~r|
i=1 j=1
Así
UNED
occ
ex [n]
φi (~r) φ∗i (~r2 ) ∗
Z Z
δEX X N
= −2 × 3
d ~r φm (~r) = 2 × d3~r ΣX (~r, ~r2 ) φ∗m (~r) si m ≤
δφm (~r2 σ) |~r2 − ~r| 2
i
ex [n]
δEX N
= 0 si m >
δφm (~r2 σ) 2
donde hemos introducido el operador de Fock, mientras que
UNED
Z
ex [n] ∗ 2 × d3~r Σ (~r , ~r) φ (~r) si m ≤ N
δEX δEX [n] X 2 m 2
= =
δφ∗m (~r2 ) δφm (~r2 ) 0 si m > N
2
Tras algunas manipulaciones algebraicas y cambios de índices, se llega a que
N/2 ∞ ∗
b X φi φm (~r1 ) φi (~r1 ) + c.c.
UNED
X X D E
ΛX (~r1 ) = 2 × φm Σ
εm − εi
m=1 i= N +1
2
donde
D E Z
(σ)
φm ΣX φi = d3~r d3~r2 φ∗m (~r) ΣX (~r, ~r2 ) φi (~r2 )
b
es el elemento de matriz del operador de Fock entre orbitales ocupados y desocupados. El potencial
UNED
de intercambio exacto será entonces solución de la ecuación integral
Z
d3~r χ0 ~r, ~r 0 vX ~r 0 = ΛX (~r)
Un detalle técnico que es importante constatar es que
UNED
d3~r χ0 ~r, ~r 0 = 0
por lo que la ecuación integral da como resultado vX (~r ) salvo constante multiplicativa, que puede
d3~r vX (~r ) = 0 en sistemas
R
obtenerse imponiendo que lı́mr→∞ vX (~r) = 0 en sistemas localizados y que
extensos.
10. En la implementación de Görling [Phys. Rev. Lett. 83, 5459 (1999)] no se calcula directamente el
potencial vX (~r), ρX (~r) que genera un potencial de Hartree igual a vX (~r). La ventaja
sino la densidad
de esta propuesta radica en que ρX (~ r) (que normaliza a −1) está localizada en sistemas nitos y
resulta más sencilla de expresar en una expansion en funciones ortonormales que el potencial vX (~ r),
el cual recordemos decae como −1/r .
172
UNED
Método de Kohn-Sham II
Máster FSC. Funcional de la densidad
7.2.c. El método KLI
11. La di cultad evidente de la derivada variacional de un funcional dependiente de los orbitales
UNED
KS es la evaluación de la función respuesta estática χ0 (~r, ~r 0 ) y su inversión. Una simpli cación muy
interesante es debida a Krieger, Li y Iafrate (KLI), que consiste en sustituir todas las diferencias
energéticas εm − εi por un único valor ε. De esta manera, para el caso general,
N/2
1 X X
χ0 (~r1 , ~r 2 ) ' −2 × φ∗m (~r1 ) φm (~r2 ) φi (~r1 ) φ∗i (~r2 ) + c.c.
ε
m=1 i6=m
UNED
∞ Z
1 X ∂B X δB
ΛB (~r1 ) ' ε |φm (~r1 )|2 − d3~r2 φm (~r1 ) φi (~r2 ) φ∗i (~r1 ) + c.c. .
ε ∂ εm δφm (~r2 )
m=1 i6=m
r1 ) φ∗i (~r2 )
P
Si ahora tenemos en cuenta que i φi (~ = δ (~r1 − ~r2 ) ,
N/2
1 X
χ0 (~r1 , ~r2 ) ' 2 × −δ (~r1 − ~r2 ) |φm (~r1 )|2 + c.c.
ε
UNED
m=1
N/2
1 X
= 2× 2 |φm (~r1 )|2 |φm (~r2 )|2 − n (~r1 ) δ (~r1 − ~r2 )
ε
m=1
Por otro lado
∞ Z
1X 2 ∂B δB 2 3 δB
ΛB (~r1 ) ' ε |φm (~r1 )| − φm (~r1 ) − |φm (~r1 )| d ~r2 φm (~r2 ) + c.c.
UNED
ε ∂ εm δφm (~r1 ) δφm (~r2 )
m=1
Ahora bien, estamos suponiendo que εm − εi es constante, por lo que resulta razonable despreciar el
término ∂ B/∂εm y, entonces
∞ Z
1X 2 3 δB δB
ΛB (~r1 ) ' |φm (~r1 )| d ~r2 φm (~r2 ) − φm (~r1 ) + c.c.
ε δφm (~r2 ) δφm (~r1 )
m=1
UNED
De este modo, la derivada variacional vB (~r1 ) = δB [n] /δn (~r1 ) veri ca la igualdad
P∞ 2 δB δB
|φm (~r1 )| 2 × Θm hφm |vB | φm i − d3~r2 φm (~r2 )
R
m=1 + φm (~r1 )
δφm (~r2 ) δφm (~r1 )
vB (~r1 ) ' +c.c.
2 × n (~r1 )
con Θm = 1 para estados ocupados y Θm = 0 para orbitales desocupados. Esta expresión nos evita
calcular e invertir la función respuesta resolviendo la anterior igualdad iterativamente.
UNED
12. Para la expresión exacta del intercambio, la aproximación KLI se reduciría a
h D Ei Z
PN/2 2
m=1 |φm (~
r1 )| hφm |b
vX | φm i − φm b X φm
Σ + Re φm (~r1 ) d3~r Σ r, ~r1 ) φ∗m (~r)
X (~
vX (~r1 ) ' 2
n (~r1 )
ecuación que, como hemos dicho, puede resolverse iterativamente. En este caso, la aproximación KLI
es muy satisfactoria, sobre todo en sistemas localizados. A pesar de las simpli caciones realizadas
sigue exhibiendo el comportamiento asintótico exacto −1/r y la gran mayoría de los tests realizados
indican que no hay diferencias sustanciales entre los resultados KLI y los que se tienen obteniendo el
potencial de intercambio a partir de las expresiones exactas.
173
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
7.3. EL PROBLEMA DEL GAP
7.3.a. Planteamiento
UNED
1. Recordemos que dado un sistema de N electrones, la energía de ionización y la a nidad electrónica
están de nidas, respectivamente, como
(N −1) (N ) (N ) (N +1)
I = E0 − E0 ; A = E0 − E0
(M )
donde E0 es la energía fundamental de un sistema de M electrones. De acuerdo con la teoría general
UNED
(N ) (N +1)
I = −εN ; A = −εN +1
(M )
siendo εm la energía del nivel KS m-ésimo del sistema con M electrones. De nimos el gap como la
diferencia entre la energía de ionización y la a nidad electrónica:
(N −1) (N +1) (N ) (N +1) (N )
Egap = I − A = E0 + E0 − 2E0 = εN +1 − εN ≥0
Así, dentro del formalismo KS podemos calcular el gap de dos maneras. En primer lugar obteniendo las
UNED
energías fundamentales de los sistemas con N − 1, N y N + 1 electrones. En segundo lugar calculando
los niveles HOMO de los sistemas con N y N + 1 electrones. En ambos casos deberíamos obtener los
mismos resultados si usásemos el funcional XC exacto.
Obsérvese que también podemos calcular el gap de la siguiente manera:
h i h i
(N ) (N ) (N +1) (N ) KS
Egap = εN +1 − εN + εN +1 − εN +1 ≡ Egap + ∆KS
UNED
(N ) (N )
El primer término, εN +1 − εN , es el gap del sistema cticio KS con N electrones ya que para un
sistema no interactuante, el potencial de ionización y la a nidad electrónica son los opuestos de los
(N +1) (N )
niveles HUMO y LUMO, respectivamente. El segundo término, εN +1 − εN +1 , indica el cambio en
el nivel KS N + 1-ésimo cuando consideramos N +1 electrones y N electrones. Para sistemas nitos
esta forma de escribir el gap no tiene mucha utilidad, sin embargo es esencial en el caso en el que
consideremos sistemas extensos (sólidos), ya que N es, formalmente, in nito y no podemos realizar el
cálculo del gap directamente dentro del marco de la DFT. La única forma de abordarlo sería realizar
UNED
una serie de cálculos con sistemas nitos de tamaños crecientes y extrapolar el resultado en el límite
de tamaño in nito lo que, en la práctica, es tremendamente complicado.
Ahora bien, en un sistema extenso, la adición de un único electrón hace que la densidad aumente
un in nitésimo y permanezca, de facto, invariable. Por tanto el per l del potencial efectivo KS va a
ser el mismo independientemente de que consideremos N electrones o N +1 electrones. Sin embargo
no podemos olvidar la discontinuidad del potencial XC. Como añadir un electrón es, remarquemos,
añadir un in nitesimal a la densidad, tendremos que
UNED
(N +1) (N )
vS (~r) − vS (~r) = ∆XC
donde
δEXC [n] δEXC [n]
∆XC = −
δn (~r) (N +1)
n0 δn (~r) (N )
n0
es la discontinuidad del potencial XC. Por tanto, en sistemas extensos
KS
Egap = Egap + ∆XC
y sólo podremos calcular directamente el gap en la DFT si la discontinuidad ∆XC es cero, esto es, si
toda la discontinuidad del funcional de Hohenberg-Kohn está contenida en el funcional cinético.
174
UNED
Método de Kohn-Sham II
Máster FSC. Funcional de la densidad
2. En un metal, Egap = 0. No se ha encontrado hasta la fecha caso alguno de metal cuyo sistema
cticio KS asociado sea aislante. En el hipotético caso de que ocurriese,
KS .
∆XC = −Egap
UNED
Sin embargo, en semiconductores y aislantes, el gap del sistema KS bajo la aproximación
LDA/GGA es bastante menor que el valor experimental del gap real. De hecho, en sistemas como
el germanio y en fases de alta presión del silicio, el sistema KS es metálico en la LDA/GGA, mientras
que el sistema real es semiconductor. Más aún, para todos los aislantes las aproximaciones LDA/GGA
dan como resultado un gap KS que es sistemáticamente menor que el gap real. La pregunta que ha
condicionado buena parte de la investigación básica en la DFT en los últimos veinticinco años es
la siguiente: si ∆XC ' 0 para sólidos no conductores, entonces una buena aproximación funcional
UNED
nos permitirá obtener el gap del sistema, y los malos resultados LDA/GGA son consecuencia de las
aproximaciones inherentes que se realizan. Sin embargo, si ∆XC 6= 0 debemos renunciar a calcular el
gap de un semiconductor a partir de la DFT, debiendo acudir o bien al cálculo extrapolador que hemos
comentado anteriormente o a técnicas de teoría de perturbaciones de muchos cuerpos. Obsérvese que
la teoría básica de la DFT ni a rma ni niega cuál es el valor de ∆XC . Este es el problema del gap en
la DFT.
7.3.b. Resultados y discusión
UNED
3. El problema del gap pareció estar resuelto (en sentido negativo para la DFT) tras una serie
de trabajos realizados por Godby, Sham y Schlüter (GSS) alrededor de 1990. Esencialmente, estos
autores calcularon, usando funciones de onda LDA, el potencial XC para una serie de semiconductores
y aislantes típicos (Si, GaAs, AlAs y C -fase diamante-) tratando exactamente el intercambio y la
correlación bajo la RPA. Observaron que, tras este cálculo, el gap KS apenas cambiaba en comparación
con el valor LDA, por lo que la conclusión es que el gap KS exacto se puede aproximar bastante bien con
UNED
la LDA y que no es igual al gap real, para cuya obtención habría de tenerse en cuenta también ∆XC .
Muy recientemente, Gruëning, Marini y Rubio han con rmado estos resultados tras una evaluación
autoconsistente del potencial EXX/RPA y con unos parámetros de convergencia más exigentes que en
el cálculo original GSS.
Estas conclusiones parecen contradecir nuestras conclusiones para sistemas localizados. En la
LDA/GGA hay errores de autointeracción y, por tanto, el electron del nivel KS HOMO siente
UNED
un potencial XC en exceso repulsivo. Esto implica que su energía LDA/GGA ha de ser mayor que
la exacta reduciéndose así el gap. Téngase en cuenta que el nivel LUMO habría de tener la misma
tendencia, pero cuantitativamente en menor escala. Ahora bien, no debemos olvidar que en un sistema
extenso los efectos de largo alcance (incluyendo el apantallamiento) son más importantes que en
sistemas nitos o, dicho de otra manera, que en sistemas extensos (mientras no existan electrones
localizados) los efectos de la autointeracción espúrea no van a ser tan determinantes. Esto explica
por qué los gaps KS en la LDA/GGA son sólo algo menores que los quasi-exactos, ya que estas
aproximaciones tratan de forma muy razonable los efectos conjuntos de intercambio y correlación.
UNED
4. Sin embargo, Görling y colaboradores han demostrado que el gap KS bajo una aproximación
EXX/LDA o EXX/GGA (intercambio exacto, correlación bajo la LDA o la GGA) es muy similar al
gap real para toda una serie de semiconductores y aislantes. Este resultado es esperanzador desde un
punto de vista práctico, puesto que abre de nuevo la puerta al cálculo DFT de gaps en materiales
extensos sin acudir a los costosísimos cálculos de teoría de muchos cuerpos. Ahora bien, se basa en
una muy afortunada cancelación de errores que conviene explicar físicamente, para lo que hay que
atender a dos factores: el ya mencionado de la autointeracción y el apantallamiento. Como ya hemos
comentado, eliminar la autointeracción espúrea disminuye el potencial repulsivo que siente un electrón
sobre todo en las regiones en las que la densidad electrónica es mayor: baja la energía del nivel HOMO
y aumenta el gap. Por otro lado, considerar la correlación supone tener en cuenta de manera más
175
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
efectiva el apantallamiento, que se traduce en que el potencial efectivo que siente un electrón es más
uniforme que áquel que actuaría sobre el mismo si no tuviésemos en cuenta la correlación: en sistemas
extensos la correlación, tratada adecuadamente, debe tender a disminuir el gap KS ya que suaviza
UNED
el potencial efectivo total que sufre el electrón.
Ya hemos comentado que en las aproximaciónes LDAy GGA, la compensación entre intercambio
y correlación se trata de manera razonable. Si sustituyésemos el intercambio LDA o GGA por el
exacto estamos disminuyendo el apantallamiento y corrigiendo efectos de autointeracción, por lo que
EXX/LDA
vXC (~r) va a ser más abrupto aunque, sobre todo, por el primer motivo: el gap KS aumenta
bastante con respecto de la LDAo de la GGA, acercándose así al valor experimental. Sin embargo, si
ahora sustitutimos la correlación LDA por un tratamiento adecuado (por ejemplo la RPA), se vuelve a
UNED
suavizar el potencial total, recuperándose esencialmente la forma original LDA o GGA aunque el gap
KS quasi-exacto va a ser algo mayor (pero no mucho) que el LDAo GGA gracias a que ahora ya no
hay efectos de autointeracción. Interesantemente, la contribución de intercambio a la discontinuidad
∆XC (que se puede calcular exactamente) es varias veces mayor que Egap . Esto pone de mani esto
otra vez la compensación entre intercambio y correlación.
7.4. EL MÉTODO DE KS EN SISTEMAS PERIÓDICOS
UNED
7.4.a. El método de ondas planas
1. Consideremos un cierto sistema físico que exhibe una simetría periódica descrita por una red
cristalina. Cada punto de la red cristalina (CL) lo designaremos genéricamente por ~,
R mientras que
cada punto de la red recíproca (RL) será designado por ~.
G Por de nición de red periódica
UNED
~ ·R
exp iG ~ =1
Si f (~r) es una función periódica cualquiera, sabemos que ésta puede escribirse como una expansión
en ondas planas con momentos ~
G pertenecientes a la red recíproca,
CL
X
f (~r) = ~ exp iG
f G ~ · ~r ,
~
G
UNED
donde ~
f (G) es la transformada de Fourier de la función periódica f (~r) de nida mediante una integral
sobre la celda unidad (UC) cuyo volumen es Ω:
Z
~ 1
~ · ~r
f G = d3~r f (~r) exp −iG
Ω UC
A su vez, si g (~r) es una función no periódica, su expansión en ondas planas no está restringida
UNED
a componentes con momentos de la red recíproca, sino que contiene contribuciones de momento
arbitrario. Puesto que todo punto del espacio recíproco puede escribirse unívocamente como suma
de un vector ~
G
de la red recíproca y de un vector ~k en el seno de la primera zona de Brillouin (BZ)
3
del espacio recíproco, podemos escribir que
CL
d3~k X ~ ~
Z h i
g (~r) = g k + G exp i ~k + G
~ · ~r
BZ (2π)3 ~
G
donde ahora
Z h i
~ ~
g k+G = g (~r) exp −i ~k + G
~ · ~r
R3
3
Recuérdense las correspondencias BZ↔UC y RL↔CL.
176
UNED
Método de Kohn-Sham II
Máster FSC. Funcional de la densidad
2. Imaginemos ahora un sistema de N fermiones no interactuantes que están sometidos a un potencial
UNED
periódico
CL
X
vS (~r) = ~ exp(iG
vS (G) ~ · ~r).
~
G
De acuerdo con el Teorema de Bloch, los órbitales del sistema electrónico (excluida la dependencia en
espín) van a tener la forma
φ~k,n (~r) = u~k,n (~r) exp i~k · ~r
UNED
donde ~k es un vector de la primera zona de Brillouin y n un índice de banda. La función u~k,n (~r) es
periódica y, por tanto, puede expandirse en ondas planas con momentos ~
G
RL
X
u~k,n (~r) = ~ exp iG
c~k,n G ~ · ~r
~
G
UNED
A efectos prácticos, debemos escoger un convenio de normalización para las funciones u~k,n (~r). El
más habitual es
RL
X 2
~
c~k,n G =1
~
G
que equivale a Z
2
d3~r u~k,n (~r) =Ω
UNED
UC
Empero, siempre debemos comprobar este tipo de detalles (al igual que los coe cientes multiplicativos
en las de niciones de las transformaciones de Fourier) al consultar diferentes textos y/o programas.
Usando este convenio, la densidad electrónica está dada por
d3~k X
Z 2
n (~r) = 2 × Θ µ − ε~ u~k,n (~r)
BZ (2π)3 n k,n
UNED
donde Θ es la función escalón de Heavyside, µ es el potencial químico del sistema, ε~k,n es la energía
del orbital φ~k,n (~r) y el factor 2 aparece como resultado de la suma sobre los grados de libertad de
espín (la generalización a sistemas con polarización de espín es inmediata y la omitimos). Es sencillo
probar que la representación de la densidad en la red recíproca es
RL
d3~k X
Z X
~ ∗ ~ 0 ~ 0 ~
n G =2× Θ µ − ε~k,n c~k,n G c~k,n G + G
(2π)3 n
UNED
BZ ~0 G
En particular, la densidad electrónica media del sistema es n = N/Ω, siendo N el número de electrones
por celda unidad. Se cumple que
d3~k X
Z
N
~ ~
n= =2× Θ µ − ε~k,n = n G = 0 ,
Ω BZ (2π)3 n
expresión que relaciona el potencial químico del sistema con su densidad media (igual bajo el convenio
que estamos siguiendo a la componente con ~ = ~0
G de la densidad electrónica).
El carácter periódico de las funciones u~k,n (~r) sugiere resolver la ecuación de Schrödinger
monoelectrónica E E E
h φ~k,n = b
b t + vbS φ~k,n = ε~k,n φ~k,n
177
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
en espacio recíproco, esto es, en la base de ondas planas formada por los autoestados del momento con
autovalor ~k + G
~. Usando la metodología que ya vimos en la primera sesión práctica, inmediatamente
tenemos que los coe cientes ~
c~k,n (G) veri can la ecuación matricial de autovalores
UNED
RL RL
X
~
X
~ 0 1 ~ ~ 2
~ ~ 0
~ 0 = ε~ c~
~
hG,
~ G~ 0 k c~k,n G = k + G δ G,
~ G~ 0 + vS G − G c~k,n G k,n k,n
G
2
~0
G ~0
G
ecuación que, naturalmente, está restringida a vectores de la red recíproca. Naturalmente esta no es
la única opción posible, puesto que siempre podremos escoger otra base de representación diferente
UNED
para las autofunciones φ~k,n (o, equivalentemente, para las funciones u~k,n ).
Las ventajas de la representación en ondas planas son las siguientes. En primer lugar es bastante
"natural", ya que explota directamente las propiedades de una función periódica en el espacio recíproco.
En segundo lugar aprovecha que el operador de energía cinética es diagonal en esta base. En tercer
lugar permite usar de manera rutinaria transformaciones rápidas de Fourier (FFT) para pasar de la
representación en espacio recíproco a la representación en espacio real y viceversa cuando sea necesario.
Por último, los parámetros de convergencia son muy sencillos: el número de vectores de la red recíproca
UNED
que se usan (de nidos, como ya sabemos, a partir de un cuto en energías) y el muestreo (sampling)
de la primera zona de Brillouin.
Este último es general (independientemente del método de representación), ya que está relacionado
con la necesaria discretización de un continuo de estados. Afortunadamente, existe una amplísima
literatura relativa a las elecciones óptimas de los puntos ~k basadas en la evaluación precisa de las
integrales del tipo
d3~k ~
Z
g k
UNED
BZ (2π)3
usando relativamente pocos puntos ~k (teoría de puntos especiales). En general, se escogen una serie
de puntos ~
k de modo que
d3~k ~ X
Z
g k ' w s g ~ks
(2π)3
e
BZ s
donde w
es es un peso de integración asociado al s-ésimo punto ~k de nuestro muestreo de la zona de
UNED
Brillouin.
7.4.b. El método KS en la representación de ondas planas
3. Puesto que el método KS se basa en la obtención de la energía y la densidad fundamental de
un sistema electrónico mediante la construcción de un sistema cticio de partículas independientes
sometidas a un potencial vS (~r), la aplicación del método de ondas planas descrito en la subsección
UNED
anterior únicamentr requiere el conocer cómo se evalúan las diferentes componentes de vS (~r) y,
también, las distintas componentes de la energía total.
4. Así, supongamos que dentro de la celda unidad hay Nat átomos, cada uno de los cuales está
localizado en un punto ~τ i de la misma. De esta forma, el potencial exterior que actúa sobre los
electrones es
Nat
CL X
X −Zi
vext (~r) =
~ i=1
~ − ~τ i
~r − R
R
Zi es
donde
el número atómico del átomo i-ésimo de la celda unidad. Evidentemente, vext (~r) =
~
vext ~r + R si R~ ∈ CL, ya que el potencial exterior exhibe la simetría periódica del sistema. Si ahora
178
UNED
Método de Kohn-Sham II
Máster FSC. Funcional de la densidad
2
tenemos en cuenta que 4π/ ~k + G
~ es la transformada de Fourier de la función 1/r,
UNED
Nat
CL X RL
−Zi
Z
4π ~ ~ ~
ei(k+G)·(~r−R−~τ i )
X X
3~
vext (~r) = 3 d k 2
~ i=1
(2π) BZ ~ ~
k+G ~
R G
Pero
CL CL
~ ~ ~ ~ ~ 1
e−i(k+G)·R =
X X
e−ik·R = (2π)3 δ ~k
Ω
~
R ~
R
UNED
y entonces, la transformada de Fourier del potencial exterior es
RL
" N # Nat
at
X 1X −4πZi −iG·~
~ τi ~ r
iG·~ ~ = 1
X −4πZi −iG·~
~ τi
vext (~r) = 2
e e ⇒ vext G 2
e
Ω G Ω G
~
G i=1 i=1
El potencial clásico de Hartree viene dado por
UNED
Z
1
d3~r0 n ~r0
vH (~r) =
R3 |~r − ~r0 |
Si ahora sustituimos n (~r) y 1/ |~r − ~r0 | por sus transformadas de Fourier tenemos que
Z RL Z RL
~ r0 1 4π ~ ~ 0
ei(k+G)·(~r−~r )
X X
3 0 ~ 0 eiG·~
vH (~r) = d ~r n G d3~k
R3 ~0
(2π)3 BZ ~ ~k + G
~
2
G G
UNED
y operando
RL 4π 4πn G~
X
~ ~ r
iG·~ ~
vH (~r) = n G e ⇒ vH G =
G2 G2
~0
G
Esto es, la expresión en el espacio recíproco del potencial de Hartree es especialmente sencilla.
UNED
5. La energía de Hartree por celda unidad es
Z
1
WH [n] = d3~r n (~r) vH (~r)
2 UC
De nuevo, si introducimos las transformadas de Fourier correspondientes vemos que
UNED
2
Z RL RL ~
RL n G
1 X
~
~ 0 ~ G
i(G+ ~ 0 )·~
r Ω X
~
~
X
WH [n] = d3~r n G vH G e = n G vH −G = 2πΩ
2 UC 2 G2
~ G
G, ~0 ~
G ~
G
Análogamente, la energía de interacción entre los núcleos y los electrones es, por celda unidad,
Z
Vext [n] = d3~r n (~r) vext (~r)
UC
y procediendo de igual manera
RL RL Nat
X
~
~
X X
~ Zi iG·~
~
Vext [n] = Ω n G vext −G = −4π n G e τi
G2
~
G ~
G i=1
179
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
Finalmente, la energía de interacción entre los núcleos por celda unidad está dada por
Nat
Z !
1 X
d3~r
UNED
Enn =− Zi δ (~r − ~τ i ) vext (~r)
2 UC i=1
y operando
RL Nat X
Nat
2π X X Zi Zj iG·(~
~
Enn = e τ i −~τ j )
Ω G2
~ i=1 j=1
G
Observamos que hemos escrito los tres términos de la energía de interacción clásica como sumas
UNED
sobre el espacio recíproco. Ahora bien
2
~
RL n G
X 2πN 2 1
lı́m 2πΩ =
G→0 G 2 Ω G2
~
G
RL Nat 2
X X
~ Zi iG·~
~ τi
= − 4πN 1
lı́m −4π n G e
G→0 G2 Ω G2
UNED
G~ i=1
Nat X
RL X Nat 2
2π X Z Z ~ τ i −~τ j )
i j iG·(~ = 2πN 1
lı́m e
G→0 Ω G2 Ω G2
~ i=1 j=1
G
donde hemos usado la neutralidad del sistema:
P
Z i = N = n ~ = 0 . Por tanto, en las sumas sobre
G
i
elementos de la red recíproca, los términos con G = 0 divergen separadamente para WH [n], Vext [n] y
UNED
Enn , pero estas divergencias se cancelan entre sí. Por tanto, es legítimo eliminar dichas divergencias y
considerar que
2
RL
X ~
n G
WH [n] = 2πΩ
G2
~ =~0
G6
RL Nat
X X
~ Zi iG·~
~
Vext [n] = −4π n G e τi
UNED
G 2
~ =~0
G6 i=1
RL Nat X
Nat
2π X X Zi Zj iG·(~
~
Enn = e τ i −~τ j ) ≡ EEwald
Ω G2
~ =~0 i=1 j=1
G6
El último término es la llamada energía de Ewald.
UNED
6. En el contexto de la teoría de Kohn-Sham, el potencial de intercambio-correlación, vXC (~r), y
la energía XC por partícula en el punto ~r, εXC (~r), también van a ser periódicos. Por lo general, se
evalúan en espacio real y mediante FFT se obtienen sus formas de Fourier ~
vXC G y ~ .
εXC G De
esta manera, la energía XC por celda unidad es
RL
X
EXC [n] = Ω ~ εXC −G
n G ~
~
G
y el potencial efectivo KS es
Nat
!
~ = 1
4π X
~
vS G ~ −
4πn G Zi e−iG·~τ i ~
+ vXC G
G2 Ω
i=1
180
UNED
Método de Kohn-Sham II
Máster FSC. Funcional de la densidad
(nótese como, de nuevo, las divergencias en G=0 se cancelan).
UNED
7. Aunque esta forma de resolver las ecuaciones KS es especialmente adecuada para sistemas
periódicos, también sirve para resolver sistemas localizados con tal de considerar una celda unidad
su cientemente grande (supercelda) y situar una réplica del sistema en dicha supercelda. En la
situación ideal (las réplicas no interactuan entre sí), la dependencia en ~k de las autoenergías ε~k,n
y las autofunciones ~
c~k,n (G) es muy pequeña, por lo que basta resolver las ecuaciones KS para ~k = ~0
y, en consecuencia, la integración de cualquier función sobre la primera zona de Brillouin se reduce a
UNED
d3~k ~
Z
1 ~ ~
g k → g k=0
BZ (2π)3 Ω
donde hemos usado el hecho de que el volumen de la primera zona de Brillouin es (2π)3 /Ω.
7.4.c. Pseudopotenciales
UNED
8. Tal y como hemos planteado el problema, todos los electrones del sistema aparecen en el juego.
Sin embargo es lógico pensar que los electrones más cercanos al núcleo (electrones de core) van a
estar descritos por orbitales KS prácticamente independientes de la geometría considerada. Sólo los
electrones de valencia, responsables de las propiedades de enlace sico-químico van a estar descritos por
estados dependientes del sistema físico especí co. De esta manera, podemos imaginar que los orbitales
KS de los electrones del core son los mismos para el sistema que para el átomo aislado, de manera que los
UNED
electrones de valencia están sometidos a un potencial efectivo, denominado pseudopotencial, resultado
de la acción conjunta de los núcleos atómicos y los electrones de core. Naturalmente el pseudopotencial
incorpora efectos púramente cuánticos (intercambio-correlación) y ello explica su naturaleza no local.
Como se puede ver en el material complementario, existen diferentes procedimientos para generar
pseudopotenciales de primeros principios, aunque hemos de recordar que la implementación práctica
de un pseudopotencial di ere de su expresión exacta. Ello implica que el pseudopotencial es, en rigor,
una aproximación (simpli cación) del problema original.
UNED
En general, dada una especie atómica a con carga iónica Zα (es decir, Zα es la diferencia entre el
número atómico y el número de electrones del core), si el ion está localizado en el origen, la actuación
del pseudopotencial w
bσ sobre un orbital KS está dada por
Z
bα φ (~r) = wαloc (~r) φ (~r) + d3~r 0 wαn.l. ~r, ~r 0 φ ~r 0
w
wαloc (~r) es la llamada parte local del
UNED
donde pseudopotencial, que tiende a −Za /r cuando r 0,
mientras que wα (~ r 0 ) es la llamada parte
n.l. r , ~ no local, que posee menor alcance efectivo. Hemos de
mencionar que el pseudopotencial se obtiene a partir de la resolución KS de la especie α. Por tanto,
por consistencia, si queremos resolver el problema global usando una cierta aproximación a la energía
XC, esa misma aproximación ha de usarse a la hora de generar el pseudopotencial.
9. La metodolgía general es muy similar al caso all-electron que hemos visto en los epígrafes
anteriores. Si de nimos las formas de Fourier de las dos componentes del pseudopotencial asociado a
la especie α mediante las igualdades
Z Z
wαloc (~q) = d3~r wαloc (~r) e−i~q·~r ; wαn.l (~q1 , ~q2 ) = d3~r1 d3~r2 e−i~q1 ·~r1 wαn.l. (~r1 , ~r2 ) e−i~q2 ·~r2
R3 R3
181
UNED
FUNCIONAL DE LA DENSIDAD.
Máster FSC. Funcional de la densidad
se tiene que el pseudopotencial total está dado, en la base de ondas planas, por
atN
~ 1 X
~ · ~τ α
wloc G wαloc (~q) exp −iG
UNED
=
Ω
α=1
Nat
1 X
wn.l ~k + G
~ 1 , ~k + G
~2 = ~ 1 · ~τ α wn.l ~k + G
exp iG α
~ 1 , ~k + G
~ 2 exp −iG
~ 2 · ~τ α
Ω
α=1
de modo que el hamiltoniano KS está de nido por los elementos de matriz
1 2
UNED
hG,
~ G~0
~
k = ~
k + ~
G δ ~ G~0 + v H
~
G − ~
G 0 +v
XC
~
G − ~
G 0 + w loc G~ − ~
G 0 + w n.l ~k + ~
G, ~
k + ~
G 0
2 G,
A efectos de obtención de energías, debemos sustituir la energía de interacción con el potencial
exterior por la energía de interacción del sistema electrónico con el pseudopotencial. Eliminando las
divergencias cuando ~ → ~0,
G se tiene que
loc n.l
Vext [n] = Vext [n] + Vext [n]
UNED
con
RL N X Nat Z
loc
X
~ loc
~ 1 3 loc Zα
Vext [n] = Ω n G w −G + d ~r wα (~r) +
Ω Ω r
G~ α=1
d3~k X X
Z
n.l ∗ ~ ~ n.l ~ ~ ~ ~
Vext [n] = 3 Θ µ − ε~k,n c~k,n G 1 c~k,n G 2 w α k + G 1 , k + G 2
BZ (2π) n ~ ~ G 1 ,G 2
UNED
El segundo término de
loc [n]
Vext sólo depende de los iones (y, por tanto, de los electrones del core)
y es, por tanto, independiente de las posiciones de los iones dentro de la celda unidad. En muchos
programas estos dos términos se dan por separado. Evidentemente, en estas expresiones, N es el
número de electrones de valencia por celda unidad.
7.5. CONCLUSIONES Y PROBLEMAS ABIERTOS
UNED
1 En este tema hemos tratado los fundamentos básicos de la teoría del funcional de la densidad sobre
todo en la formulación de Kohn-Sham. En la actualidad, el método KS bajo las aproximaciones LDA
y GGA es un estándar en Física de Materia Condensada. El hecho de que el problema de N electrones
interactuantes se pueda reducir a un problema equivalente de N partículas independientes permite
abordar problemas con decenas de miles de electrones (naturalmente con recursos computacionales
muy avanzados - cálculo paralelo masivo -). Las aplicaciones de la DFT van, pues, desde problemas de
UNED
física atómica y molecular, biofísica, física de estado sólido y diseño de materiales en nanotecnología.
2 No nos olvidemos que la DFT es una teoría para el estado fundamental. Es decir, la solución
autoconsistente de las ecuaciones KS sólo nos da la densidad y la energía fundamentales. No debe
entonces buscarse ningún signi cado físico al potencial efectivo KS ni a los autovalores del sistema
cticio KS (con la expresión del HOMO que, en la teoría exacta, es igual a menos la energía de
ionización del sistema) tal y como hemos ilustrado en la discusión del problema del gap. Empero, como
hemos visto a lo largo de estas notas, conocida la energía del sistema electrónico podemos evaluar (en
la aproximación Born-Oppenheimer) la fuerza neta que actua sobre cada núcleo del sistema global, lo
que abre la puerta al análisis detallado de sus propiedades esctructurales. La evaluación de primeros
principios de propiedades de estados excitados ha de hacerse mediante otras herramientas que, en
muchos casos, incluyen cálculos KS como pasos intermedios.
182
UNED
Método de Kohn-Sham II
Máster FSC. Funcional de la densidad
3 Obviamente, como corresponde a un curso de posgrado, no hemos visto todos los tópicos de interés
UNED
en el campo de la DFT (que supondría un curso en exceso extenso). Algunos aspectos interesantes
que no hemos abordado son, por ejemplo:
El formalismo DFT para sistemas relativistas.
La DFT aplicada a sistemas con varios tipos de fermiones.
La extensión de la DFT a sistemas con temperatura nita y su conexión con la Mecánica
UNED
Estadística.
La relación entre aspectos formales de la DFT y otros característicos de Mecánica Estadística.
Otras implementaciones de la DFT distintas a la de Kohn y Sham.
Propiedades más formales (reglas de suma) del funcional XC y del funcional cinético.
La implementación de la DFT en códigos de dinámica molecular (el método de Car Parrinello).
UNED
4 Repitamos que la DFT es hoy en día un estándar de cálculo, y bien puede decirse que es una teoría
tremendamente robusta. Sin embargo aun hay retos abiertos (algunos de naturaleza matemática). Sin
ánimo alguno de ser exhaustivo, las líneas generales de investigación sobre la DFT (no líneas de
investigación que usan la DFT y que son tan amplias como se quieran habida cuenta la aplicabilidad
general del método) están dirigidas sobre dos vías principales:
UNED
1. Desarrollo de códigos e cientes que permitan tratar sistemas con centenares de miles de
electrones. Ello implica una investigación muy técnica orientada, fundamentalmente, al desarrollo
de algoritmos y a la representación óptima de funciones de onda.
2. Búsqueda de funcionales XC con precisión química capaces de describir procesos de interés en
quimica y bioquímica con exactitud similar a los métodos CI o DMC con un coste computacional
menor.
UNED
En cierta medida, estas dos vías son incompatibles. La DFT se puede aplicar a sistemas muy
complejos gracias a que el funcional XC se aproxima mediante la LDA o la GGA. En el momento en
el que se usen aproximaciones más so sticadas, generalmente incluyendo los propios orbitales KS, el
cálculo de la energía y el potencial XC sería incluso más costoso que la propia solución de las ecuaciones
KS. Ahora bien, no es tan grave. Al estudiar un sistema muy complicado formado por miles de átomos,
habitualmente interesa tener sólo una idea general de las estructuras y de su evolución (experimentos
en el ordenador) y en el hipotético caso de que las limitaciones de la LDA/GGA fuesen importantes,
UNED
los aspectos que se pierden se pueden incluir a mano de manera aproximada. Por ejemplo, ni la LDA
ni la GGA describen fuerzas de van der Waals, pero estas se pueden introducir de manera efectiva
en los códigos DFT orientados, por ejemplo, a problemas biofísicos. El interés en mejores funcionales
XC está sobre todo motivado por la caracterización de propiedades de sistemas más simples, para
los cuales sí resulta interesante el disponer de herramientas muy precisas con una fuerte capacidad
predictiva orientadas al conocimiento preciso de los mecanismos evolutivos (reacciones químicas, por
ejemplo) a nivel molecular.
En cualquier caso, el conocimiento completo de las propiedades electrónicas exige abordar también
y simultáneamente el problema de los estados excitados. La teoría del funcional de la densidad
dependiente del tiempo (TDDFT) es una posibilidad, que está ganando en los últimos años, combinada
con la DFT estática que hemos tratado aquí, un lugar cada vez más prominente en Química Teórica,
Bioquímica y Física de la Materia Condensada.
183
Máster FSC. Funcional de la densidad
UNED
UNED
UNED
UNED
UNED
UNED
UNED
También podría gustarte
- Gas de Fermi de Electrones LibresAún no hay calificacionesGas de Fermi de Electrones Libres23 páginas
- Teoría del Funcional de la Densidad DFTAún no hay calificacionesTeoría del Funcional de la Densidad DFT7 páginas
- Física Estadística: Límite ClásicoAún no hay calificacionesFísica Estadística: Límite Clásico16 páginas
- Funcional de Densidad Local en FísicaAún no hay calificacionesFuncional de Densidad Local en Física10 páginas
- Practico7-Electrones Libres en MetalesAún no hay calificacionesPractico7-Electrones Libres en Metales8 páginas
- Relacion7 - Gases Cuanticos Soluciones v0Aún no hay calificacionesRelacion7 - Gases Cuanticos Soluciones v016 páginas
- Funciones de Partición en Gases CuánticosAún no hay calificacionesFunciones de Partición en Gases Cuánticos21 páginas
- Tema 4. Química Cuántica: 1. Introducción A La Mecánica CuánticaAún no hay calificacionesTema 4. Química Cuántica: 1. Introducción A La Mecánica Cuántica19 páginas
- Ejercicios de Física del Estado SólidoAún no hay calificacionesEjercicios de Física del Estado Sólido4 páginas
- Materia Condensada, Marder, SolucionesAún no hay calificacionesMateria Condensada, Marder, Soluciones54 páginas
- Modelo de Electrones Libres en MetalesAún no hay calificacionesModelo de Electrones Libres en Metales18 páginas
- Gas de Electrones Libres y SuperfluidezAún no hay calificacionesGas de Electrones Libres y Superfluidez8 páginas
- Calor específico de metales a bajas temperaturasAún no hay calificacionesCalor específico de metales a bajas temperaturas16 páginas
- Física de Sólidos: Problemas ResueltosAún no hay calificacionesFísica de Sólidos: Problemas Resueltos5 páginas
- Entropía y Teoremas en Sistemas FísicosAún no hay calificacionesEntropía y Teoremas en Sistemas Físicos3 páginas
- Distribucion de Fermi Dirac y Bose EinsteinAún no hay calificacionesDistribucion de Fermi Dirac y Bose Einstein4 páginas
- 4 Problemas Drude Sommerfield PDFAún no hay calificaciones4 Problemas Drude Sommerfield PDF2 páginas
- Átomos Multielectrónicos y Mecánica CuánticaAún no hay calificacionesÁtomos Multielectrónicos y Mecánica Cuántica20 páginas
- Teoría Funcionales Densidad en CatálisisAún no hay calificacionesTeoría Funcionales Densidad en Catálisis18 páginas
- Examen Final. Solucionario. CF4F1 2023-IAún no hay calificacionesExamen Final. Solucionario. CF4F1 2023-I18 páginas
- Ordinario MecanicaEstadistica 2022a2023Aún no hay calificacionesOrdinario MecanicaEstadistica 2022a20235 páginas
- Discucion 4 Gas de Electrones de Fermi 2025 2Aún no hay calificacionesDiscucion 4 Gas de Electrones de Fermi 2025 25 páginas
- Gases Ideales de Bosones y Fermiones. Gas de Bose y Gas de FermiAún no hay calificacionesGases Ideales de Bosones y Fermiones. Gas de Bose y Gas de Fermi15 páginas
- Ecuación Van der Waals y Punto CríticoAún no hay calificacionesEcuación Van der Waals y Punto Crítico12 páginas
- Determinación del Periodo de SemidesintegraciónAún no hay calificacionesDeterminación del Periodo de Semidesintegración3 páginas
- Efecto Peltier. Celdas de PeltierAún no hay calificacionesEfecto Peltier. Celdas de Peltier6 páginas
- Efecto Fotoeléctrico - Fundamento TeóricoAún no hay calificacionesEfecto Fotoeléctrico - Fundamento Teórico7 páginas
