Métodos de cálculo computacional mecánico-cuántico para el estudio de estructuras electrónicas Alberto Jesús Uribe Jiménez

1.

Introducción

Se parte, en general, de un sistema compuesto por un gran número de átomos que interactúan entre si, en cualquiera de sus formas y afectados o no por un campo eléctrico externo. Puede escribirse el Hamiltoniano como sigue: ˆ = H

∇2i i=1 2m

−~2

PN

−

PP

− ~2

∇2i I=1 2MI

PP

+

e2 2

PN

1 i,j=1;i6=j |ri −rj |

(1) I=1

ZI e2 i=1 |ri −RI |

PN

+

e2 2

PP

ZI ZJ I,J=1;J6=I |RI −RJ |

Donde R es el conjunto de las coordenadas nucleares y P el número de estas; r es el conjunto de las coordenadas electrónicas y N el número de las mismas; Z las cargas nucleares y M las masas nucleares [1, 2, 5]. Todas las propiedades están descritas por la función de onda Ψi (r, R, t), resolviendo la ecuación de Schrödinger: ˆ i (r, R, t) = Ei Ψi (r, R, t) HΨ

(2)

La resolución de esta ecuación es en la práctica una tarea muy compleja y es aquí donde los siguientes métodos nos ofrecen alternativas viables a este cálculo [1, 2, 5].

2.

Teoría de funcionales de la densidad (DFT)

En 1964 Pierre Hohenberg y Walter Kohn desarrollaron la base teórica de la teoría de funcionales de la densidad (DFT, Density Functional Theory) y demostraron que la energía está determinada por la densidad electrónica del estado fundamental. El método desarrollado, DFT, calculaba la primera a partir de la segunda: E0 = E0 [ρ0 ].[3, 4] Posteriormente Kohn y Sham demostraron que es posible extraer una ecuación para los orbitales, a partir de los cuales puede obtenerse la densidad anteriormente mencionada. Es a partir de este momento cuando el método cobra mayor relevancia.[4] 1

Hoy en día es un método fundamental y ampliamente usado en el campo de la física y la química cuántica.

2.1.

Teoremas de Hohenberg y Kohn

Hohenberg y Kohn demostraron los dos teoremas siguientes: Teorema I: Existe una correspondencia uno a uno entre la densidad del estado base de un sistema de muchos electrones y el potencial externo que se genera al suponer la aproximación de Born - Oppenheimer.[2] Es decir, el potencial externo está determinado por la densidad electrónica de forma unívoca, excepto por una constante aditiva. Corolario: Dado que la densidad electrónica determina el potencial, entonces, también determina la función de onda del estado base. Como consecuencia inmediata se tiene que el valor esperado del estado base de cualquier observable, es un único funcional de la densidad electrónica exacta del estado base.[2] Existe una única función de onda de un sistema con un gran número de electrones, definida por la ecuación de Schrödinger, y a partir de la cual se obtiene una densidad electrónica con la misma información que la mencionada función de onda.[2] Teorema II: − la densidad electronica de prueba ρ0 (→ r ) no negativa y normalizada a ´Dada − → − 0 → N, ρ ( r ) d r = N , entonces: − E [ρ0 (→ r )] ≥ E0

(3)

y ˆ 0

0

E [ρ ] = F [ρ ] +

− − − ρ0 (→ r )v(→ r )d→ r

(4)

donde D E ˆ Ψ [ρ0 ] F [ρ0 ] = Ψ [ρ0 ] Tˆ + U

(5)

siendo Ψ [ρ0 ] el estado fundamental con ρ0 (r) como densidad de estado fundaˆ la energía de la interacción electrón-electrón. mental, Tˆ la energía cinética y U [2] La densidad electrónica del estado fundamental es aquella que minimiza el − funcional de energía E [ρ(→ r )].[2, 4] Además, el funcional de energía puede ser dividido en tres partes:

2

− − − − E [ρ (→ r )] = T [ρ (→ r )] + U [ρ (→ r )] + Eext [ρ (→ r )] (6) − − donde T [ρ (→ r )] es la energía cinética, U [ρ (→ r )] la energía de interacción → − entre electrones y E [ρ ( r )] la energía externo debido a la carga del núcleo. ext

[2, 3]

2.2.

Ecuación de Kohn-Sham

Khon y Sham elaboraron su teoría partiendo de los estudios realizados por Hartree y Fock. Se aproximó un sistema de electrones que interactuaban entre ellos a uno en que no ocurría esto, sino que cada electrón está sometido a un potencial creado por las demás partículas, manteniéndose las densidades del estado base. En estas circunstancias, el funcional de energía se expresa como:[2, 3] − − − − − EKS [ρ (→ r )] = TS [ρ (→ r )] + EH [ρ (→ r )] + Eext [ρ (→ r )] + Exc [ρ (→ r )]

(7)

− − con TS [ρ (→ r )] la energía cinética de una sola partícula, Exc [ρ (→ r )] la energía → − de correlación-intercambio y EH [ρ ( r )] la energía de Hartree, que se define por con la expresión:[5] ˆ 2 − e ρs (→ r 0) 3→ d − r0 (8) EH = − − |→ r −→ r 0| Las ecuaciones con las que se minimiza el potencial de Kohn Sham son:[2] − − ˆ KS ψi (→ H r ) = Ei ψi (→ r)

(9)

y sustituyendo se obtiene[2]   ~2 2 → − → − → − − − ∇ + VH ( r ) + Vxc ( r ) + Vext ( r ) ψi (→ r ) = Ei ψi (→ r) − 2me

(10)

donde Vxc está definida por[2] − δExc [ρ (→ r )] − Vxc [ρ (→ r )] = → − δρ ( r )

(11)

− δEH [ρ (→ r )] − VH [ρ (→ r )] = → − δρ ( r )

(12)

VH es

y ψi (r) las N soluciones de las ecuaciones de Kohn-Sham que minimizan la energía. Con ellas se puede obtener la densidad electrónica del estado fundamental como

3

− ρ (→ r)=

N X

− − ni ψi∗ (→ r ) ψi (→ r)

(13)

i=1

entendiendo ψi (r) como los orbitales, ni los números de ocupación de estos orbitales y N el número total de electrones en el sistema.[4, 6] Las soluciones a este sistema se resuelven con un proceso iterativo, con una densidad inicial que nos proporciona los potenciales Vxc y VH , a partir de los cuales se obtienen los orbitales de Kohn-Sham y con ellos se obtiene de nuevo el valor de la densidad electrónica.[2, 6] Es lo que aparece en la figura 1.

Figura 1: Algoritmo para DFT con el método de Kohn-Sham.[7]

3.

Aproximaciones al potencial de correlación-intercambio

Es fundamental en la teoría de la DFT tener una buena aproximación al potencial de correlación-intercambio, relacionándose este con la energía de correlación4

intercambio como aparece en la ecuación 11. Por otra parte, pueden obtenerse las ecuaciones ˆ ˆ → − → − → − → − − − Exc [ρ ( r )] = ρ ( r ) εxc [ρ ( r )] d r = exc [ρ (→ r )] d→ r − − − δ (ρ (→ r ) εxc [ρ (→ r )]) δexc [ρ (→ r )] − Vxc [ρ (→ r )] = = − − δρ (→ r) δρ (→ r)

(14) (15)

que relaciona la energía de correlación-intercambio y el potencial de correlación− intercambio con la energía de correlación-intercambio de cada partícula εxc [ρ (→ r )] → − y con la energía de correlación intercambio por unidad de volumen exc [ρ ( r )]. Con estas relaciones se pueden plantear varias aproximaciones. LDA y GGA son las más comunes y usadas, y se tratan en los apartados posteriores.

3.1.

Aproximación de densidad local (LDA)

− Se considera la energía de correlación-intercambio de cada partícula εxc [ρ (→ r )] como un funcional dependiente exclusivamente de la densidad electrónica. Las − [ρ (→ r )]) e intercambio contribuciones de esta energía a la correlación (εLDA c → − [ρ ( r )]) se corresponden con cada uno de los miembros de la siguien(εLDA x te ecuación[8] − − − εLDA [ρ (→ r )] = εLDA [ρ (→ r )] + εLDA [ρ (→ r )] xc x c

(16)

Se toma para la parte de la energía correspondiente al intercambio, el modelo del gas de electrones de densidad homogénea y una densidad ρ constante εLDA x

9 − [ρ (→ r )] = − α 4



3 8π

 31

1

ρ3

(17)

con una constante α = 2/3.[8] Con este resultado y la ecuación 14 se obtiene la siguiente expresión 9 − Exc [ρ (→ r )] = − α 4



3 8π

 13 ˆ

4 − ρ 3 d→ r

(18)

y haciendo uso de la ecuación 11, se consigue el siguiente resultado para el potencial [8] VxLDA

  31 − − 1 3 3 δ (ρ (→ r ) εx [ρ (→ r )]) → − =− α ρ3 (ρ ( r )) = → − 2 π δρ ( r )

(19)

Para el caso de la energía de correlación se pueden usar varios métodos, como el de Vosko, Wilk y Nusair, aunque el más simple es el denominado Xa , − consistente en considerar εc [ρ (→ r )] = 0 y α = 2/3.[8]

5

Características de LDA: Las principales características para esta aproximación son:[6] 1. Es válida para sistemas homogéneos. 2. Predice moléculas y sólidos más ligados. 3. La tendencia química suele ser correcta. 4. Para sistemas con enlaces covalentes, iónicos y metálicos, las geometrías se reproducen bien. Las propiedades dieléctricas son sobreestimadas en, aproximadamente, el 10 %. 5. Para sistemas débilmente ligados, la longitud de los enlaces es demasiado corta. Limitaciones En sistemas finitos, no decae en el vacío con −e2 /r. Esta es una de las limitaciones más importantes. LDA funciona especialmente bien para sistemas en que la densidad electrónica es aproximadamente constante. Sin embargo falla en aquellos sistemas en que existen variaciones importantes de la densidad electrónica o de la autointeracción, así como aquellos que tienen frecuentes inhomogeidades, como sistemas con enlaces moleculares débiles.[6]

3.2.

Aproximación de gradientes generalizados (GGA)

Dadas las limitaciones de LDA para los materiales con inhomogeidades, se plantea la aproximación GGA, en la cual se introducen gradientes de la densidad electrónica, considerando así, además de la densidad en un punto concreto, como varía la misma alrededor del punto. La expresión quedará entonces:[6, 8] GGA Exc

− [ρ (→ r )] =

ˆ

− − − ρ (→ r ) εxc [ρ (→ r )] d→ r +

ˆ

− − Fxc [ρ (→ r ) , ∇ρ (→ r )] dr

(20)

Se pretende corregir el comportamiento asintótico y las propiedades de escalamiento con respecto a LDA. Para ello, comparando las siguientes ecuaciones[8] ˆ → − − − Ex [ρ ( r )] = ρ (→ r ) εx [ρ] d→ r (21) 1 − Ex [ρ (→ r )] = 2

ˆ

− − ρ (→ r 1 ) d→ r1

ˆ

− − ρxc (→ r 1, → r 2) → d− r2 → − → − | r 1 − r 2|

(22)

se obtiene que[8] 1 − εx [ρ (→ r )] = 2

ˆ

− − ρxc (→ r 1, → r 2) → d− r2 → − → − | r 1 − r 2| 6

(23)

Para ver el comportamiento asintótico se realiza el límite cuando r tiende a infinito[8, 9] 1 − l´ım εx [ρ (→ r )] = − 2r

r→∞

(24)

y por tanto 1 − l´ım Ex [ρ (→ r )] = − r→∞ 2

ˆ

− ρ (→ r) − d→ r r

1 − l´ım Vx [ρ (→ r )] = − r Con respecto a LDA, considerando que el límite de la densidad es [8] r→∞

(25) (26)

− l´ım ρ (→ r ) = e−ar

r→∞

Al hace el limite de εLDA [ρ] se muestra su comportamiento asintótico[8] x − l´ım εLDA [ρ (→ r )] = 0 x

r→∞

Este no es el comportamiento adecuado y por ello se trata de corregir en GGA. De entre toda las soluciones existentes, la corrección propuesta por Becke en 1988 es la más utilizada. En ella se introduce un término a la expresión de la energía de intercambio de la LDA, esto es[8] − − − ExGGA [ρ (→ r )] = ExLDA [ρ (→ r )] + ExN LDA [ρ (→ r )]

(27)

siendo el nuevo término ExN LDA

− [ρ (→ r )] =

ˆ

4 − − r ) f N LDA (x) d→ r ρ 3 (→

(28)

donde x=

|∇ρ|

(29)

4

ρ3

Por tanto, con las ecuaciones 14 y 28, se obtiene que [8] − − LDA εN [ρ (→ r )] = ρ1/3 (→ r ) f N LDA (x) x f

N LDA

(30)

(x) puede tomar los valores siguientes[8] x 6 ln x

(31)

βx2 1 + 6βx ln x

(32)

βx2 1 + 6βx sinh−1 x

(33)

f N LDA (x) = − f N LDA (x) = − f N LDA (x) = −

7

Becke demostró que, de estos tres valores de f N LDA (x), el que proporciona una mejor aproximación es correspondiente a la ecuación 33 con β = 0,0042.[8]

3.3.

Perdew-Burke-Ernzerhof (PBE)

La expresión para la expresión de la energía correlación intercambio según este método es[10, 12] ˆ − − − − − P BE BE Exc [ρ (→ r )] = d3 rρ (→ r ) εP (rs (→ r ) , s (→ r ) , ζ (→ r )) (34) xc − donde ζ (→ r ) es la polarización de spin definida como (ρ ↑ −ρ ↓) − ζ (→ r)= ρ − rs (→ r ) es el radio de Wigner-Seitz con la expresión rs =

4πρ 3

− y s (→ r )tiene la forma de |∇ρ| − s (→ r)= 2kF ρ con kF = 3π 2 ρ


1/3

.[10]

Referencias [1] M. E. G. Méndez. Estudio teórico de la adsorción de metales de los grupos IV y V sobre superficies semiconductoras y su efecto en el crecimiento epitaxial de Si sobre Ge (001). Ensenada, Baja California, Centro de Investigación Científica y de Educación Superior de Ensenada. Centro de Ciencias de la Materia Condensada, 2000.

[2] W. L. Gonzalez Olaya. Director: Ph. D. R. Cardona Cardona. Uso de la teoría del funcional densidad (DFT) en la caracterización estructural y electrónica de la perovskita triple Sr3CoSb2O9 sintetizada en el laboratorio. Universidad Nacional de Colombia, Facultad de Ciencias, Departamento de Física. Bogota, Colombia 2013. [3] M. I. Nicolás Vazquez, E. Marín Chiñas, F. M. Castro Martínez, R. Miranda Ruvalcaba. Algunos aspectos básicos de la química computacional. Primera edición, abril 2006. Universidad Nacional Autónoma de México. ISBN 970-32-3307-4. [4] J. P. Loperena, A. R. Ballesteros, C. A. R. González, L. D. S. Vázquez. Teoría de funciones de la densidad. Principios de estructura de la materia. Universidad Nacional Autónoma de México. [5] J. Rodríguez Ruiz. Tutores: B. Biel Ruiz y F. J. García Ruiz. Estudio del efecto del strain en láminas de MOS2 mediante métodos atomísticos. Universidad de Granada. [6] Iván Scivetti. Director: Dr. A. Caro. Hidrógeno: más allá de la aproximación clásica. Instituto Balseiro, Comisión Nacional de Energía Atómica, Universidad Nacional de Cuyo, Junio 2003.

8

[7] A. F. OliveiraI, G. Seifert, T. Heine, H. A. Duarte. Density-functional based tightbinding: an approximate DFT method. J. Braz. Chem. Soc. vol.20 no.7 São Paulo 2009. [8] Juan Andrés, Juan Beltrán [eds], Química Teórica y Computacional, Castellón de La Plana, Universidad Jaime I, 2000. [9] X. Hua, X. Chen, W. A. Goddard. Generalized generalized gradient approximation: An improved density-functional theory for accurate orbital eigenvalues. Phys. Rev. B 55, 16103. June 1997. [10] M. Ernzerhof, G. E. Scuseria. Assessment of the Perdew–Burke–Ernzerhof exchangecorrelation functional. J. Chem. Phys., Vol. 110, No. 11, 15 March 1999. [11] X. Xu. The extended Perdew-Burke-Ernzerhof functional with improved accuracy for thermodynamic and electronic properties of molecular systems. J. Chem. Phys., Vol. 121, No. 9, 1 September 2004. [12] E. S. Fabian. Cálculos Computacionales de Estructuras Moleculares. Universidad de Alicante 2016-01-13.

9