Adsorcion-de-H2O-y-H2-en-superficies-de-carbono

Vista previa del material en texto

INSTITUTO DE FÍSICA, UNAM
Departamento de Estado Sólido
“ADSORCIÓN DE H2O y H2 EN
SUPERFICIES DE CARBONO”
Tesis que presenta
Eduardo Rangel Cortes
como requisito parcial
para obtener el t́ıtulo de Doctor
en F́ısica
Director: Dr. Fernando Magaña Soĺıs. México D.F. 2010
 
UNAM – Dirección General de Bibliotecas 
Tesis Digitales 
Restricciones de uso 
 
DERECHOS RESERVADOS © 
PROHIBIDA SU REPRODUCCIÓN TOTAL O PARCIAL 
 
Todo el material contenido en esta tesis esta protegido por la Ley Federal 
del Derecho de Autor (LFDA) de los Estados Unidos Mexicanos (México). 
El uso de imágenes, fragmentos de videos, y demás material que sea 
objeto de protección de los derechos de autor, será exclusivamente para 
fines educativos e informativos y deberá citar la fuente donde la obtuvo 
mencionando el autor o autores. Cualquier uso distinto como el lucro, 
reproducción, edición o modificación, será perseguido y sancionado por el 
respectivo titular de los Derechos de Autor. 
 
 
 
Jurado de Examen Doctoral 
U" .• m 
OFICIO: PCF/084/2010 
ASUNTO: Designación de jurado 
El Comité Académico del Posgrado en Ciencias Fís icas en su sesión del 9 de 
Febrero del presente ha designado como Jurado del estud iante EDUARDO 
RANGEL CORTES con número de cuenta 95049065, para dictaminar si el trabajo 
desarrollado como tesis titulado : ''Adsorción de H20 y H2 en superficies de 
carbono", dirigido por el Dr. Luís Fernando Magaña Salís, tiene los méritos para 
obtener el grado de Doctor en Ciencias (Física) conforme al plan de estudios 
5057. 
Propietario: 
Propietario: 
Propietario: 
Propietario : 
Propietario: 
Suplente : 
Suplente: 
Atentamente . 
Dr. Luis Fernando Magaña Salís 
Dr. Rubén Santa maría Ortíz 
Dr. Juan Salvador Arellano Peraza 
Dr. Gerardo Garcia Naumis 
Dr. Miguel Robles Pérez 
Dr. Juan Adrián Reyes Cervantes 
Dr. Pablo de la Mora y Palomar Askinasy 
"POR MI RAZA HABLARÁ EL EspíRITU" 
Ciudad Universitaria, D. F., a 12 de Febrero de 2010 
Coordinador del Posgrado en Cienc ias Fisicas 
21-~~~ -
Dr. Manuel Torres Labansat 
c.c.p.- Cada mi mbro del slnodo 
C.C.p - Interesado. 
C.C. p .- Exped lent . 
Agradecimientos
Es un verdadero placer dedicar este espacio
para expresar mis agradecimientos a las
personas e instituciones que han facilitado
las cosas para que este trabajo llegue a un
feliz término. Debo agradecer al Consejo
de Ciencia y Tecnoloǵıa (CONACYT) y al
Instituto de Ciencia y Tecnoloǵıa del
Distrito Federal (ICYT-DF) por el
patrocinio y el apoyo para la realización del
proyecto. De igual manera agradezco a la
Dirección General de Servicios de
Cómputo Académico, DGSCA y al centro
de supercómputo Kanbalam, UNAM, por
facilitarme recursos de cómputo con los
cuales se realizaron gran parte de los
cálculos de este trabajo. Debo agradecer de
manera especial al Dr. Luis Fernando
Magaña por aceptarme para hacer esta
tesis doctoral bajo su dirección. Por darme
los medios y recursos suficientes para mi
desarrollo como cient́ıfico. Agradezco a
toda mi familia, a mis padres y en especial
a mi madre por que siempre será mi fuerza
y inspiración para alcanzar mis metas.
Agradezco a mi novia Mari por todo su
amor y apoyo incondicional. A todos mis
amigos del pasado y del presente que han
estado ah́ı para apoyarme; David, Don
Beto, Benja, Pacheco, Gabriel, Alex, Alex
Valderrama, Andrea, Andrea Gutiérrez,
Marisol, Josue, Enrique, Betazo, Gregorio,
Abdellah, Vladimir, Keila, Edgar, Lupita,
Juanito, Ivan y todos los démas. Quiero
agradecerles por que todos ustedes son
participes de esta alegŕıa, por enseñarme
que todo se aprende y que al final todo el
esfuerzo conlleva a un triunfo.
Índice general
1. Introducción 1
2. Teoŕıa Funcional de la Densidad 5
2.1. Ecuación de Schrödinger . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.1. Aproximación Born-Oppenheimer . . . . . . . . . . . . . . . . . 5
2.2. Construcción de la Teoŕıa Funcional de la Densidad. . . . . . . . . . . . 6
2.2.1. Condiciones sobre la densidad electrónica ρ(r) para ser empleada
en la teoŕıa del funcional de la densidad. . . . . . . . . . . . . . 9
2.2.2. Método de Kohn-Sham para escribir la enerǵıa en función de los
orbitales moleculares. . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3. Enerǵıa de intercambio y correlación . . . . . . . . . . . . . . . . . . . 13
2.3.1. Aproximación de densidad local (LDA) . . . . . . . . . . . . . . 13
2.3.2. Aproximación de Gradiente Generalizado (GGA) . . . . . . . . 14
2.4. Densidad de estados . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.5. Análisis poblacional de Löwdin . . . . . . . . . . . . . . . . . . . . . . 15
2.6. Base de ondas planas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3. Pseudo-potenciales 17
3.0.1. Definición de pseudo-potencial . . . . . . . . . . . . . . . . . . . 17
3.0.2. Construcción de pseudo-potenciales . . . . . . . . . . . . . . . . 20
3.0.3. Pseudo-potenciales Troullier-Martins . . . . . . . . . . . . . . . 21
3.0.4. Pseudo-potenciales ultrasuaves-Vanderbilt . . . . . . . . . . . . 22
4. Dinámica Molecular 25
4.0.5. Termostatos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.0.6. Barostatos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
5. Adsorción de agua en grafeno-titanio. 29
5.1. Método . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
5.2. Discusión y resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
5.2.1. Grafeno-titanio . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
5.2.2. Naturaleza del enlace grafeno-titanio . . . . . . . . . . . . . . . 32
5.2.3. C2Ti-agua . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
i
ii ÍNDICE GENERAL
5.2.4. Titanio-agua . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.2.5. C8Ti − agua . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
6. Adsorción de H2 en un nanotubo-litio 43
6.1. Método . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
6.2. Discusión y resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
6.2.1. Estructura de los nanotubos de carbono . . . . . . . . . . . . . 45
6.2.2. Estructura electrónica de los nanotubos . . . . . . . . . . . . . . 46
6.2.3. Nanotubo decorado con litio . . . . . . . . . . . . . . . . . . . . 46
6.2.4. Adsorción de H2 en el nanotubo decorado con litio . . . . . . . 49
7. Adsorción de H2 en nanotubo-nitrógeno 57
7.1. Método . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
7.2. Discusión y resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
7.2.1. Nanotubo decorado con nitrógeno . . . . . . . . . . . . . . . . . 58
7.2.2. Adsorción de hidrógeno en el nanotubo decorado con nitrógeno . 59
7.2.3. Almacenamiento de hidrógeno en manojos de nanotubos decora-
dos con nitrógeno atómico. . . . . . . . . . . . . . . . . . . . . . 65
8. Adsorción de H2 en grafeno-litio 67
8.1. Método . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
8.2. Discusión y resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
8.2.1. Grafeno-litio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
8.2.2. Grafeno-litio-hidrógeno . . . . . . . . . . . . . . . . . . . . . . . 70
8.2.3. Grafito-litio-hidrógeno . . . . . . . . . . . . . . . . . . . . . . . 75
9. Sumario y Conclusiones 81
A. Funcional de correlación 87
B. Hibridación del átomo de carbono 89
B.0.4. Isómeros de carbono . . . . . . . . . . . . . . . . . . . . . . . . 89
C. Art́ıculos publicados 91
Bibliograf́ıa 105
Índice de figuras
1.1. Volumen de 4 kg de hidrógeno compactada en diferentes formas,
comparada con el tamaño de un automóvil [3]. . . . . . . . . . . . . 2
2.1. Comparación de la enerǵıa de enlace de la molécula de H2 con dis-
tintos funcionales. PZ (LDA), PW91 (GGA), LYP (GGA), PBE
(GGA). . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . 14
3.1. Estados atómicos para el ox́ıgeno. . . . . . . . . . . . . . . . . . . . . 18
3.2. Pseudo-funciones Troullier-Martins vs autofunciones para el átomo
de aluminio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.3. Pseudo-potenciales Troullier-Martins para el aluminio. Este pseudo-
potencial fue generado con el programa fhi98PP. . . . . . . . . . . . 22
5.1. Configuración de átomos de Ti en una hoja de grafeno. a) Celda
unidad que contiene cuatro átomos de Ti por ocho de C (C2Ti). b)
Celda unidad que contiene un átomo de Ti por ocho de C (C8Ti) . 30
5.2. a) La superficie de Ti(0 0 1) tiene un arreglo hexagonal con todos
los átomos en un solo plano. b) La configuración de C2Ti, también
tiene un arreglo hexagonal pero los átomos de Ti en color obscuro,
se encuentra en un plano inferior con respecto a los átomos de Ti
en color claro. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5.3. Isosuperficies de redistribución de densidad electrónica para la ad-
sorción de Ti en C50Ti. (a) La parte clara y obscura corresponden
a +0.14 y -0.14 e/Ȧ3 respectivamente. (b) La parte clara y obscura
corresponden a +0.034 y -0.034 e/Ȧ3 respectivamente. . . . . . . . 33
5.4. Isosuperficies ρdif−1(r). (a) La parte clara y obscura corresponden a
+0.013 y -0.13 e/Ȧ3 respectivamente. (b) La parte clara y obscura
corresponden a +0.005 y -0.005 e/Ȧ3 respectivamente. . . . . . . . 34
5.5. Isosuperficies ρdif−2(r). (a) La parte clara y obscura corresponden a
+0.013 y -0.13 e/Ȧ3 respectivamente. (b) La parte clara y obscura
corresponden a +0.005 y -0.005 e/Ȧ3 respectivamente. . . . . . . . 35
iii
iv ÍNDICE DE FIGURAS
5.6. Isosuperficies ρdif−3(r). (a) La parte clara y obscura corresponden a
+0.013 y -0.13 e/Ȧ3 respectivamente. (b) La parte clara y obscura
corresponden a +0.005 y -0.005 e/Ȧ3 respectivamente. . . . . . . . 35
5.7. Adsorción de agua en grafeno-titanio, para la configuración C2Ti. La
orientación inicial de la mólecula es con el átomo de O sobre la ĺınea
de los átomos de H. La configuración final demuestra la disociación
completa de la molécula en H y O. . . . . . . . . . . . . . . . . . . . 36
5.8. Adsorción de agua en grafeno-titanio, para la configuración C2Ti.
La orientación inicial de la mólecula es con el átomo de O por abajo
de la ĺınea de los átomos de H. La configuración final demuestra la
disociación en OH y H. . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.9. Índices de Miller correspondientes a tres planos del titanio. . . . . 38
5.10. A la izquierda se muestra la configuración inicial y a la derecha la
configuración final de la molécula de agua adsorbida en la superficie
(0 0 1). Las enerǵıas de adsorción son 2.42 eV, 5.93 eV y 5.22 eV
respectivamente. La celda unidad consta de 16 átomos de Ti. . . . 39
5.11. A la izquierda se muestra la configuración inicial y a la derecha la
configuración final de la molécula de agua adsorbida en la superficie
(1 0 0). Las enerǵıas de adsorción son 8.863 eV, 10.90 y 9.51 eV
respectivamente. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.12. Adsorción de agua en grafeno-titanio, para la configuración C8Ti.
La configuración final demuestra la adsorción de la molécula sin
disociación. La orientación inicial del agua no cambia el resultado
final. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.13. Isosuperficies de redistribución de densidad electrónica para la ad-
sorción de H2O en C8Ti. (a) La parte clara y obscura corresponden a
+0.067 y -0.067 e/Ȧ3 respectivamente. (b) La parte clara y obscura
corresponden a +0.033 y -0.033 e/Ȧ3 respectivamente. . . . . . . . 42
6.1. Imágenes de microscopia electrónica de nanotubos de carbono de
capa simple. a) Nanotubo de carbono aislado. b) Imagen de un
manojo de nanotubos con un arreglo periódico hexagonal [66]. . . . 44
6.2. Cuando conectamos los puntos O con A, y B con B́ un nanotubo
puede ser construido. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
6.3. Clasificación de los nanotubos de carbono: a) armchair, b) zig-zag,
c) quiral. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
6.4. Densidad de estados y bandas de enerǵıa a lo largo del nanotubo
(8,0). La enerǵıa de Fermi está en E = 0. El gap de enerǵıa es de
0.46 eV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
6.5. Optimización de átomos de Li en la superficie del nanotubo (8,0).
a) C64Li, b) C8Li, c) C4Li y d) C4Li-manojos. . . . . . . . . . . . . . 48
ÍNDICE DE FIGURAS v
6.6. Dinámica molecular a 400 K de un nanotubo con un cluster de ocho
átomos de Li. a) Cluster de ocho átomos de Li sobre un nanotubo
(8,0), b) Átomos dispersados de Li al final de la dinámica. . . . . . 49
6.7. a) Optimización de tres moléculas de H2 en el sistema C64Li. b)
Optimización de tres moléculas de H2 en C4Li−manojo. . . . . . . 50
6.8. Isosuperficies de redistribución de densidad electrónica para la ad-
sorción de H2 en C64Li. (a) La parte clara y obscura corresponden a
-0.007 y +0.007 e/Ȧ3 respectivamente. (b) La parte clara y obscura
corresponden a -0.034 y +0.034 e/Ȧ3 respectivamente. . . . . . . . 51
6.9. Bandas de enerǵıa a lo largo del nanotubo (8,0). Los K puntos están
a lo largo del tubo. La enerǵıa de Fermi está en E = 0. La adsor-
ción de Li en el nanotubo cambio las propiedades electrónicas. Al
desaparecer el gap de enerǵıa (figura 6.4), el sistema cambio de
semiconductor a conductor. . . . . . . . . . . . . . . . . . . . . . . . 51
6.10.Densidad de estados localizados en Li. Se observa que el nivel 2p
también participa en la interacción. La enerǵıa de Fermi para el
sistema C64Li esta en -2.127 eV. . . . . . . . . . . . . . . . . . . . . 52
6.11. Relación funcional entra la carga neta de Li y la enerǵıa de adsor-
ción de H2; al aumentar la densidad de Li, disminuye la enerǵıa de
adsorción de H2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6.12. Aqúı se muestra la evolución en el tiempo de C4Li, a 300 K y a 500 K
usando LDA y MD. a) La configuración final muestra 18 moléculas
de H2 enlazadas, con enerǵıa de adsorción promedio de 0.174 eV/H2;
que corresponde al 7.6 % en peso de H. b) La configuración final
muestra 13 moléculas de H2 enlazadas, con enerǵıa de adsorción
promedio de 0.100 eV/H2; que corresponde al 5.57 % en peso de H.
El paso de tiempo es de un femtosegundo. . . . . . . . . . . . . . . . 54
6.13. Aqúı se muestra la configuración final del sistema C4Li a 673 K;
Cinco moléculas de H2 son enlazadas. La región encerrada en ćırculo
corresponde a la coalescencia de átomos de Li y a una molécula de
H2 disociada. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
7.1. Geometŕıa optimizada de seis moléculas adsorbidas en un átomo de
N . Todas las distancias están en angstroms (Ȧ). . . . . . . . . . . . 60
7.2. Bandas de enerǵıa a lo largo del nanotubo (8,0). Los K puntos están
a lo largo del tubo. La enerǵıa de Fermi está en E = 0. La adsor-
ción de N en el nanotubo cambio las propiedades electrónicas. Al
desaparecer el gap de enerǵıa (figura 6.4), el sistema cambio de
semiconductor a conductor. . . . . . . . . . . . . . . . . . . . . . . . 60
7.3. Densidad de estados localizados para el N . La enerǵıa de Fermi
para el sistema C64N es -3.16 eV. El sistema tiene comportamiento
metálico. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
vi ÍNDICE DE FIGURAS
7.4. Enerǵıa (eV) vs pasos de tiempo (1 femtosegundo) para adsorción
y disociación de H2 en la superficie del sistema C8N . El estado ini-
cial es el sistema optimizado (a la izquierda). La configuración final
corresponde a 14 moléculas de H2 más dos moléculas de H2 disocia-
das (a la derecha). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
7.5. Enerǵıa (eV) vs pasos de tiempo (1femtosegundo) para adsorción
y disociación de H2 en la superficie del sistema C8N . Este es la
configuración final del sistema que se muestra en la figura “a” con
cuatro moléculas adicionales. Al final se tiene 17 moléculas de H2
más una molécula disociada. . . . . . . . . . . . . . . . . . . . . . . . 62
7.6. Enerǵıa (eV) vs pasos de tiempo (1 femtosegundo) para adsorción
y disociación de H2 en la superficie del sistema C8N . Para este caso
se agregan 4 moléculas de H2. Dichas moléculas quedan adsorbidas
en una segunda capa. La configuración final tiene 21 moléculas de
H2 más 6 átomos adsorbidos. . . . . . . . . . . . . . . . . . . . . . . 63
7.7. Sistema C32N4H6 − 11H2 a 300k. . . . . . . . . . . . . . . . . . . . . . 64
7.8. Densidad de carga para el sistema C32N4H6. a) Isosuperficie de carga
0.3 eV/u.a3; el color obscuro es la región en donde hay más densidad
de carga. b) Isosuperficie de carga 0.2 eV/u.a3. . . . . . . . . . . . . 65
7.9. Configuración de C4N . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
7.10. Configuración de C4N en manojos infinitos con hidrógeno almace-
nado. Al final el sistema almacena 1.79 % en peso de H. . . . . . . 66
8.1. Izquierda-Observación de un defecto topológico inducido por irra-
diación electrónica, antes (a) y después (b). c) Modelo del par
pentágono-heptágono en el grafeno. d) Simulación de HR-TEM mos-
trando una buena comparación con la imagen HR-TEM. Derecha-
Observación de vacancias formadas en grafeno. a) antes de la irra-
diación. b) Después de la irradiación. c) y d) Simulación de imágenes
HR-TEM. e) Simulación de imágenes HR-TEM comparadas con la
observada experimentalmente en (b) [85]. . . . . . . . . . . . . . . . 68
8.2. Optimización de átomos de Li en la superficie del grafeno a) C7Li,
b) C17Li . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
8.3. Dinámica a 473 K de grafeno con un cluster de 13 átomos de Li. a)
Configuración inicial. b) Después de 500 femtosegundos los átomos
de Li se dispersan. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
8.4. Comparación de la densidad de carga (e/ua3) para dos planos XY
antes y después del dopamiento. Las figuras a la izquierda corres-
ponden al grafeno solo y las de la derecha al grafeno contaminado
con Li. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
8.5. Optimización de las moléculas de H2 en a) C7Li, b) C17Li. . . . . . 73
8.6. Densidad de estados localizados para el sistema Li; Se observa que
el nivel 2p también participa en la interacción. La enerǵıa de Fermi
para el sistema C7Li está en 0.97 eV . . . . . . . . . . . . . . . . . . 73
8.7. Enerǵıa (eV) vs pasos de tiempo (1 femtosegundo) para adsorción
de H2 en la superficie del sistema C17Li. Al final se tiene 3 moléculas
de H2 adsorbidas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
8.8. Enerǵıa (eV) vs pasos de tiempo (1 femtosegundo) para adsorción
de H2 en la superficie del sistema C7Li. Al final se tiene 3 moléculas
de H2 adsorbidas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
8.9. Dos configuraciones del grafito contaminado con Li (C7Li). a) A-
A, b) A-B. En esta configuración las hoja A está desplazada con
respecto a la hoja B. . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
8.10. Enerǵıa potencial (eV) para dos hojas de C7Li. Se estudian dos
configuraciones A-A y A-B. La distancia esta en angstroms (Ȧ). . 76
8.11.Moléculas de H2 relajadas en la configuración A-A (C7Li). a) 48 GP,
b) Presión atmosférica. . . . . . . . . . . . . . . . . . . . . . . . . . . 77
8.12.Moléculas de H2 relajadas en la configuración A-B (C7Li). a) 48 GP,
b) Presión atmosférica. . . . . . . . . . . . . . . . . . . . . . . . . . . 77
8.13.Densidad de estados localizados para el sistema A-B; a) En Li se
observa que el nivel 2p también participa en la interacción. b) En
cualquier átomo de carbono 1, 2 o 3. La enerǵıa de Fermi para el
sistema A-B esta en 6.25 eV. . . . . . . . . . . . . . . . . . . . . . . 78
8.14. Isosuperficie de densidad de carga 0.15 e/u.a3 a) A-B, b) A-A. . . . 79
8.15. Isosuperficie de densidad de carga 0.10 e/u.a3 a) A-B, b) A-A. . . . 79
9.1. Caracterización de la adsorción de hidrógeno en nanotubos de car-
bono a 300 K [71]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
9.2. Imágenes TEM de a)Ti, b)Ni, c)Pd, d)Au, e)Al, f)Fe de 5nm de
espesor contaminando la superficie de una nanotubo de carbono [38]. 83
Resumen
Por un lado, estudiamos la adsorción de agua en grafeno decorado con titanio y
por otro, la capacidad de almacenamiento de hidrógeno en nanotubos de carbono y
en grafeno. Para este segundo caso, consideramos nanotubos decorados con litio, con
nitrógeno y al grafeno decorado con litio substituyendo a algún átomo de carbono.
El estudio se realizó utilizando la Teoŕıa de Funcionales de la Densidad y Dinámica
Molecular. Encontramos que el agua es adsorbida en el sistema grafeno-titanio. Cuando
se considera una cobertura de titanio en el grafeno de C2Ti, el agua se disocia al ser
adsorbida. Esta disociación ocurre de dos maneras. Una es formando H, O y H. La otra
es formando H y OH. Cuando la cobertura de Ti en el grafeno es C8Ti, la molécula
es adsorbida sin disociarse. Por otra parte, encontramos que un nanotubo de carbono
puede almacenar 7.5% en peso de hidrógeno cuando está decorado con Li y 6.0% en
peso de hidrógeno cuando está decorado con N . Estas cifras estań dentro de lo requerido
por el Departamento de Enerǵıa de los Estados Unidos para aplicaciones prácticas de los
sistemas de almacenamiento de hidrógeno. Sin embargo, cuando consideramos manojos
de nanotubos decorados estas cifras descienden a 2.7% y 1.8% respectivamente. Por
último encontramos, que el grafeno decorado con litio substitucional puede almacenar
6.2% en peso de hidrógeno.
Abstract
On one hand, we studied the water adsorption in graphene decorated with titanium
and on the other hand, we studied the capacity of hydrogen storage in carbon nanotu-
bes and in graphene. For this second case, we considered carbon nanotubes decorated
with lithium, with nitrogen and to the graphene decorated with lithium replacing so-
me carbon atom. We used Density Functional Theory and Molecular Dynamics. We
found that the water molecule is adsorbed in the system graphene-titanium. When we
considered titanium coverage of C2Ti, the water molecule is dissociated when adsor-
bed. This dissociation happens in two ways. One is forming H, O and H. The other
is forming H and OH. When the titanium coverage is C8Ti, the molecule is adsorbed
without dissociation. We also found that the carbon nanotube can store up to 7.5% in
weight of hydrogen when it is decorated with lithium and 6.0% in weight of hydrogen
when it is decorated with nitrogen. These figures are within the requirements from the
Department of Energy of the United States for practical applications of hydrogen sto-
rage systems. Nevertheless, when we considered bundles of decorated carbon nanotubes
these numbers descend respectively to 2.7% and 1.8%. Finally, we found that graphene
decorated with lithium can store up to 6.2% in weight of hydrogen.
Caṕıtulo 1
Introducción
Son diversas las tecnoloǵıas que se están desarrollando actualmente, con el fin de
sustituir a la gasolina en motores de combustión interna. El agotamiento de las reservas
de petróleo en el futuro, aśı como el uso de combustibles mas limpios y menos costosos,
dan importancia al hecho de la búsqueda de un combustible alternativo.
El hidrógeno surge como una alternativa a los combustibles fósiles; desde el punto
de vista ecológico, el hidrógeno no genera en su combustión los contaminantes de los
combustibles convencionales, tales como los gases de efecto invernadero. Por otro lado
desde el punto de vista económico la tecnoloǵıa del hidrógeno podŕıa revolucionar el
mercado del transporte yde la enerǵıa, lo que ha comenzado a conocerse como la
economı́a del hidrógeno [1].
Hay que hacer notar que el hidrógeno no es una enerǵıa primaŕıa y es preciso ob-
tenerlo mediante alguna reacción qúımica, bien de disociación1 del agua o por alguna
reacción en que el hidrógeno es un producto de la reacción. Los métodos industriales
de producción de hidrógeno se basan principalmente en los combustibles fósiles. Por
tanto se encuentra que para obtener concentraciones de hidrógeno prácticas a partir de
la disociación de agua, es preciso temperaturas superiores a los 2000 K(1727oC)[2]. A
esta temperatura la disociación del agua es solo del 1%, por lo que habŕıa de operarse
al menos a 2200 K(2073oC).
El empleo de catalizadores de platino o rutenio permite operar a temperaturas
inferiores a 1000oC [2]. La necesidad de soslayar las altas temperaturas y las bajas
concentraciones de hidrógeno de la disociación térmica directa del agua, nos lleva a
emprender la búsqueda de nuevos catalizadores, capaces de reducir las temperaturas de
operación y aumentar las concentraciones de hidrógeno.
Por otro lado, también es necesario encontrar la forma más eficiente de almacenar
hidrógeno, para que pueda ser una alternativa viable a los combustibles actuales. Hoy
en d́ıa ya existen sistemas que se han utilizado para almacenar hidrógeno, y ninguno
es satisfactorio en su totalidad [3]. Una forma es almacenar hidrógeno ĺıquido a tem-
1 Disociación en qúımica es un proceso general en el cual complejos, moléculas o sales se separan
en moléculas mas pequeñas, iones o radicales.
1
2 CAPÍTULO 1. INTRODUCCIÓN
Figura 1.1: Volumen de 4 kg de hidrógeno compactada en diferentes formas, com-
parada con el tamaño de un automóvil [3].
peratura de 21 K(−252oC) alcanzandose grandes densidades de hidrógeno molecular
(H2), sin embargo este sistema resulta ser muy costoso y requiere mucha enerǵıa (figura
1.1). Otra forma de almacenamiento de hidrógeno es como gas comprimido a 20 MP
(200 atmósferas); aunque mas barato que el hidrógeno ĺıquido, ocupa mayor volumen
y requiere contenedores pesados y de grandes dimensiones que complicaŕıan su uso en
automóviles.
Además del almacenamiento como gas comprimido o ĺıquido, existen otros métodos;
como el almacenamiento en hidruros métalicos (figura 1.1), que absorben 5 a 8% en
peso, pero liberan hidrógeno a un ritmo muy lento [3]. Otro forma de adsorción esta
en los sistemas porosos como son los nanotubos de carbono2 de pared única; hace
algunos años se descubrió que grandes cantidades de hidrógeno pod́ıan ser almacenadas
en estructuras de grafito microscópicamente pequeñas con forma de tubo [4]. Es un
campo de investigación actual, pero hasta hoy los resultados sobre su capacidad de
almacenamiento no son muy alentadores. Se necesitan avances cient́ıficos y técnicos que
ratifiquen el alto potencial de esta tecnoloǵıa. Las propiedades de almacenamiento de
hidrógeno de los nanotubos de carbono y de las hojas de grafito no están completamente
exploradas ni entendidas. Sin embargo, hay interés cient́ıfico particularmente en los
nanotubos y en las hojas de grafeno contaminadas con elementos qúımicos, ya que
podŕıan ser un medio prometedor para el almacenamiento seguro de hidrógeno3 .
La mayoŕıa de los investigadores coinciden en que se deben emprender investigacio-
nes en el desarrollo de métodos para una adecuada caracterización de estos materiales
carbonosos.
Actualmente existen algunos trabajos teóricos y experimentales de almacenamiento
de hidrógeno en materiales de carbono contaminados con elementos qúımicos. Se ha
2 Se sabe que el carbono forma diferentes estructuras. Esto se debe a la hibridación de los orbitales
atómicos (Véase apéndice B)
3 La propuesta del Departamento de Enerǵıa de E. U. A., es la búsqueda de materiales que alma-
cenen el 6 % en peso de hidrógeno [5].
3
encontrado que los metales de transición como contaminantes sobre los nanotubos,
no son buenos candidatos, ya que tienden a formar nanopart́ıculas que degradan el
almacenamiento de hidrógeno [6, 7]. Para el caso de nanotubos contaminados con Li,
existen controversias teóricas y experimentales que se tratarán en el caṕıtulo 6.
En este trabajo, se realizó un estudio a nivel teórico para encontrar materiales que
puedan almacenar el 6% en peso de hidrógeno. Estudiamos los nanotubos de carbono
y las hojas de grafeno contaminados con algunos elementos qúımicos. También se estu-
dió al titanio depositado en hojas de grafeno como catalizador para la disociación de
agua y obtención de hidrógeno.
Los objetivos de esta tesis son los siguientes:
• Estudiar el almacenamiento de hidrógeno molecular (H2) en nanotubos de car-
bono de capa única puros y contaminados con átomos de nitrógeno y litio a temperatura
ambiente y presión atmosférica.
• Estudiar el almacenamiento de hidrógeno molecular (H2) en hojas de grafeno
puras y contaminadas con aluminio, litio y azufre a temperatura ambiente y presión
atmosférica.
• Estudiar la disociación de agua en hojas de grafeno contaminadas con titanio.
• Determinar bajo qué condiciones, las hojas de grafeno contaminadas con titanio
disocian a las moléculas de agua.
Para estudiar el almacenamiento de hidrógeno y la obtención de hidrógeno a partir
de la disociación del agua, en los materiales mencionados anteriormente, será necesario
abordar el problema mediante metodoloǵıas que nos permitan describir el comporta-
miento de los sistemas a nivel cuántico. Aśı que en el caṕıtulo 2 se introduce el concepto
de Teoŕıa Funcional de la Densidad. Se empieza con el planteamiento del hamiltoniano
molecular y la aproximación Born-Oppenheimer. Se plantean las ecuaciones de Kohn-
Sham como una forma para obtener la densidad electrónica. Para conocer la reactividad
de una superficie se introduce el concepto de la densidad de estados y del análisis de
Löwdin.
En el caṕıtulo 3 se plantea el uso de pseudo-potenciales como una aproximación, para
sustituir al núcleo de los átomos y a los electrones internos por un potencial externo.
Por tanto, solo se considera los electrones de valencia para los enlaces qúımicos. El uso
de pseudo-potenciales nos permite, que el tiempo de computo sea mas corto.
Para abordar la evolución temporal de los núcleos de un sistema, en el caṕıtulo
4, se introduce el concepto de Dinámica Molecular. El movimiento de los núcleos son
gobernados por las ecuaciones de Newton. Para simular la temperatura y la presión se
plantea el uso de termostatos y barostatos.
Finalmente en los caṕıtulos 5, 6, 7 y 8 presentamos los resultados obtenidos y en el
4 CAPÍTULO 1. INTRODUCCIÓN
caṕıtulo 9 presentamos una discusión de las conclusiones.
Caṕıtulo 2
Teoŕıa Funcional de la Densidad
2.1. Ecuación de Schrödinger
La ecuación de Schrödinger fue desarrollada por el f́ısico austŕıaco Erwin Schrödinger
en 1925. El problema fundamental reside, como sabemos, en la tarea de resolver la
ecuación diferencial dada por
Ĥψ = i�
∂
∂t
ψ (2.1)
La forma del operador Hamiltoniano para un sistema de M núcleos y N electrones
es:
Ĥ =
N∑
i=1
− �
2
2me
∇2i +
M∑
I=1
− �
2
2MI
∇2I +
1
4π�0
N,M∑
i,I=1
e2oZI
| ri −RI | +
+
1
8π�0
N,M∑
i,j=1,i �=j
e2o
| ri − rj | +
1
8π�0
M,M∑
I,J=1,I �=J
e2oZIZJ
| RI −RJ | (2.2)
donde ri es la coordenada del electrón i, RI es la coodenada del núcleo I, ZI la carga
del núcleo I. Los términos de la ecuación 2.2 corresponde a la enerǵıa cinética de los
electrones, la enerǵıa cinética de los núcleos, la interacción núcleo-electrón, la interac-
ción electrón-electrón y la interacción núcleo-núcleo respectivamente. Es el término de
interacción entre electrones y iones, el que hace muy dif́ıcil la resolución de la ecuación
de Schrödinger ya que acopla las coordenadas de todos los pares de part́ıculas. La solu-
ción anaĺıtica exacta de dicha ecuación es desconocida incluso para átomoso moléculas
muy sencillas con lo que es necesario recurrir a métodos aproximados.
2.1.1. Aproximación Born-Oppenheimer
La primera aproximación realizada para la resolución de la ecuación de Schröringer
para un sistema de M núcleos y N electrones, es la aproximación de Born-Oppenheimer.
Esta aproximación se basa en que los núcleos tienen mucho mayor masa que los elec-
trones y por tanto el movimiento de los núcleos y los electrones puede considerarse en
5
6 CAPÍTULO 2. TEORÍA FUNCIONAL DE LA DENSIDAD
forma independiente. Esta aproximación simplifica los cálculos, al permitir separación
de variables fijando la posición de los núcleos y resolviendo numéricamente la ecuación
electrónica. Por tanto podemos suponer que los electrones seguirán instantáneamente
la dinámica nuclear, dado que la variación de las posiciones de los núcleos será extre-
madamente lenta en la escala de tiempo del movimiento electrónico.
Resolveremos el problema electrónico a partir de la ecuación Schröringer conside-
rando fijos los núcleos:
Ĥelec({RI})ψ({RI}) = �elec({RI})ψ({RI}) (2.3)
Donde el hamiltoniano electrónico es:
Ĥelec =
N∑
i=1
− �
2
2me
∇2i +
1
4π�0
N,M∑
i,I=1
e2oZI
| ri −RI | +
1
8π�0
N,M∑
i,j=1,i�=j
e2o
| ri − rj | (2.4)
La enerǵıa total Etot es la enerǵıa electrónica �elec más la enerǵıa de repulsión núcleo-
núcleo
Etot = �elec({RI}) + 1
8π�0
M,M∑
I,J=1,I �=J
e2oZIZJ
| RI −RJ | (2.5)
Para resolver numéricamente la ecuación 2.3, con el hamiltoniano mostrado en la
ecuacion 2.4, uno puede escoger entre una lista variada de metodoloǵıas. Elegir un
método u otro depende de la naturaleza de nuestro sistema en estudio, las propiedades
que se requieren calcular, los recursos computacionales con que se cuenten, la exactitud
requerida en los resultados, entre otras cosas. En nuestro caso particular, hemos elegido
la Teoŕıa de Funcional de la Densidad porque nos permite abordar problemas de muchos
cuerpos con relativa facilidad y con un tiempo de cálculo requerido relativamente corto
con respecto a otras metodoloǵıas y si a esto le agregamos la implementación de pseudo-
potenciales el tiempo se reduce aún más. Por otro lado, para tratar el movimiento de los
núcleos se ocupara dinámica molecular. En caṕıtulos posteriores se tratara este tema.
2.2. Construcción de la Teoŕıa Funcional de la Den-
sidad.
La Teoŕıa de Funcional de la densidad parte de la teoŕıa de Thomas-Fermi-Dirac,
la cual señala que las propiedades f́ısicas de un sistema pueden ser calculadas a partir
de la densidad electrónica [8].
Sin embargo se observa de la ecuación 2.4, que en los términos están involucrados
operadores de uno y dos electrones y por tanto parece no ser necesaria la función de onda
completa para evaluar la enerǵıa. Para calcular estos términos nos alcanza con conocer
promedios de la función de onda sobre las configuraciones de los démas electrones, como
se demostrará a continuación.
2.2. CONSTRUCCIÓN DE LA TEORÍA FUNCIONAL DE LA DENSIDAD. 7
Primero construimos la densidad electrónica a partir de la función de onda total Ψ.
La densidad electrónica en un punto es la suma de las probabilidades de encontrar a
cada uno de los electrones en ese punto independiente de la posición de los demás
ρ(r) =
∫
· · ·
∫
Ψ∗(r1, r2, ..., rn)Ψ(r1, r2, ..., rn)dr2dr3 · · · drN +∫
· · ·
∫
Ψ∗(r1, r2, ..., rn)Ψ(r1, r2, ..., rn)dr1dr3 · · · drN +
+ · · ·
∫
· · ·
∫
Ψ∗(r1, r2, ..., rn)Ψ(r1, r2, ..., rn)dr2dr3 · · · drN−1 (2.6)
En donde tenemos N términos. Podemos intercambiar las coordenadas de un par de
electrones si tenemos en cuenta que el signo cambia y convertimos a cualquiera de las
integrales de la ecuación anterior a una integral igual a la primera, por tanto:
ρ(r) = N
∫
· · ·
∫
Ψ∗(r1, r2, ..., rn)Ψ(r1, r2, ..., rn)dr2dr3 · · · drN (2.7)
Por tanto el valor esperado para un operador de un electrón para una función de
onda es:
〈Ψ|
∑
i
Oi(r)|Ψ〉 =
∑
i
∫
· · ·
∫
Oi(ri)Ψ
∗(r1, r2, ..., rn)Ψ(r1, r2, ..., rn)dr1dr2dri · · · dN =
= N
∫
· · ·
∫
Oi(ri)|Ψ(r1, r2, ..., rn)|2dr1dr2dri · · · dN =
∫
ρ(r)O(r)dr
(2.8)
Análogamente para un operador de dos electrones tenemos que:
〈
∑
i
Oi(r, ŕ)〉 =
∫ ∫
ρ(r)ρ(ŕ)O(r, ŕ)drdŕ (2.9)
De modo que para los términos de interacción electrón-núcleo y electrón-electrón
tenemos los siguientes valores de expectación respectivamente:
〈Ψ|
∑
i
Zi
| R− ri | |Ψ〉 =
∫
ρ(r)
| R− r |dr =
∫
vnuc(r)ρ(r)dr (2.10)
〈Ψ|
N,M∑
i,j=1,i�=j
e2o
| ri − rj | |Ψ〉 =
∫ ∫
ρ(r)ρ(ŕ)
| r− ŕ | drdŕ = J [ρ] (2.11)
8 CAPÍTULO 2. TEORÍA FUNCIONAL DE LA DENSIDAD
Donde vnuc y J [ρ] es el potencial del núcleo y el término clásico de interacción
electrón-electrón.
Hohenberg y Kohn sentaron las bases matemáticas de la teoŕıa del funcional de la
densidad, dándole formalismo al trabajo iniciado por Thomas, Fermi y Dirac. Sin entrar
en muchos detalles, podemos decir que las ideas fundamentales que sentaron las bases
matemáticas de la teoŕıa de la funcional de la densidad son [8]:
• Las variables que determinan las propiedades f́ısicas del sistema son el número de
part́ıculas n y el potencial externo vext(r), porque dependen de la densidad electrónica
ρ(r). Por tanto, la enerǵıa total depende de estas dos variables: E = En,vext .
• Las variables n y vext, son suficientes y necesarias para obtener la información del
estado base del sistema en estudio.
• A partir de la densidad electrónica ρ(r), podemos calcular el número de part́ıculas
en la forma
∫
ρ(r)d3r = n.
Hohenberg y Kohn legitimaron el uso de la densidad electrónica para obtener las
propiedades del estado base de un sistema, proponiendo dos teoremas, que nos limita-
remos a mencionar. Su desarrollo puede verse en el art́ıculo original de Hohenberg y
Kohn publicado en 1964 [9].
Primer teorema de Hohenberg-Kohn: Existe una relación biuńıvoca entre el
potencial externo vext(r), salvo por una constante aditiva, y la densidad electrónica ρ(r).
Este teorema nos indica que, una distribución electrónica dada está asociada a
un único potencial externo. A grandes rasgos, podemos resumir el primer teorema de
Hohenberg y Kohn de la siguiente manera: “dado un potencial externo vext, la enerǵıa
Ev del sistema es un funcional de la densidad Ev[ρ]”.
El problema es que seguimos teniendo, para cualquier sistema, dos cantidades des-
conocidas: Ev[ρ] y la propia ρ(r). Y aunque supongamos que conocemos al funcional de
enerǵıa Ev[ρ], no tenemos forma de conocer a la variable ρ(r). Sin embargo, el segundo
teorema de Hohenberg y Kohn nos proporciona un método para encontrar la ρ(r) que
nos de un valor Ev[ρ] cercano o igual al de la enerǵıa real del sistema en su estado basal.
Segundo teorema de Hohenberg-Kohn: Para una densidad de prueba ρ̃(r), tal
que ρ̃(r) ≥ 0, ∫ ρ̃(r)d3r = n, el funcional Evext [ρ̃] es una cota superior a la enerǵıa real
del estado base Ev,o, es decir:
Evext,o ≤ Evext [ρ̃] (2.12)
En busca de una breve interpretación de este teorema, tomemos como función de
onda de prueba de nuestro sistema a ψ̃, para calcular el valor esperado de la enerǵıa
Ev , dado un potencial externo vext,∫
ψ̃∗Hψ̃ d3r = Ev[ρ̃] =
∫
ρ̃(r)vext(r) d
3r + FHK [ρ̃] (2.13)
2.2. CONSTRUCCIÓN DE LA TEORÍA FUNCIONAL DE LA DENSIDAD. 9
La enerǵıa aśı calculada debe ser una cota superior a la enerǵıa real del estado base
dado el potencial externo vext(r),
Ev,o[ρ] ≤ Ev[ρ̃]
donde la variable ρ(r) es la densidad electrónica real del estado basal. Esto se demuestra
empleando el principio variacional.
Como consecuencia, el segundo teorema de Hohenberg y Kohn justifica el empleo
del cálculo variacional en la teoŕıa de Thomas-Fermi-Dirac al proponer a ETFD como
una aproximación de E[ρ].
Los teoremas de Hohenberg y Kohn demostraron la existencia y unicidad del funcio-
nal de la densidad FHK [ρ], que involucra la enerǵıa cinética del sistema y el potencial
de interacción electrón-electrón:
FHK [ρ] = T [ρ] + Vee[ρ]
2.2.1. Condiciones sobrela densidad electrónica ρ(r) para ser
empleada en la teoŕıa del funcional de la densidad.
Las dos cantidades básicas en DFT son el funcional de enerǵıa E[ρ] y la densidad
electrónica ρ(r), como se discutió con anterioridad. En el esquema de Hohenberg y
Kohn la densidad electrónica a emplear será aquella que minimice la enerǵıa del estado
basal.
Aún cuando podemos proponer una ρ(r) que satisfaga el requerimiento de minimizar
la enerǵıa E[ρ], no cualquier ρ(r) tiene sentido f́ısico. De acuerdo a los postulados de
la Mecánica Cuántica, toda densidad electrónica debe tener su origen en una función
de onda ψ para que represente a un sistema f́ısico. Por tanto, resulta indispensable
establecer las condiciones para que la ρ(r) propuesta se derive de una función de onda.
La condición de que ρ(r) deba obtenerse de una función de onda ψ, establece ciertos
requerimientos para determinar esta ρ(r).
Debe ser un número positivo o cero,
ρ(r) ≥ 0
Al integrarla sobre todo el espacio, debemos recuperar
el número total de electrones del sistema∫
ρ(r) d3r = n
y, debe ser cuadráticamente integrable∫
|∇ρ(r)1/2|2 d3r < ∞
En el marco DFT, es suficiente que la densidad cumpla los tres requerimientos arriba
citados, porque entonces la ρ(r) representa a un sistema f́ısico.
10 CAPÍTULO 2. TEORÍA FUNCIONAL DE LA DENSIDAD
2.2.2. Método de Kohn-Sham para escribir la enerǵıa en fun-
ción de los orbitales moleculares.
Dado el desconocimiento del funcional de enerǵıa cinética (desconocido aún en las
ecuaciones de Hohenberg y Kohn), Kohn y Sham [10] propusieron redefinir las ecuacio-
nes de DFT en términos de orbitales moleculares, en lugar de la densidad electrónica
ρ(r).
El valor de la enerǵıa cinética se escribe en términos de estos orbitales como:
Tnirs =
n∑
i=1
∫ ∞
0
φ∗i
(
−1
2
∇2
)
φi d
3r (2.14)
El sub́ındice nirs hace referencia a que se trata de la enerǵıa cinética de un sistema
de referencia no interactuante. La densidad electrónica también se puede escribir en
términos de orbitales moleculares (al sustituir el operador ∇2 por la identidad).
ρ(r) =
n∑
i=1
∑
s
|φi(r, s)|2 (2.15)
Las ecuaciones (2.14) y (2.15), presentan una consideración adicional, los orbitales
involucrados deben estar plenamente ocupados, esto es, no hay orbitales parcialmente
ocupados o desocupados. Estas ecuaciones implican un costo computacional menor
debido al número restringido de orbitales.
Sin olvidar que el problema principal es resolver un sistema de n part́ıculas interac-
tuantes, Kohn y Sham tuvieron la brillante idea de mapear un sistema interactuante a
un sistema no interactuante de part́ıculas.
Para un sistema no interactuante, la funcional de enerǵıa cinética está dado por la
ecuación (2.14) y por su parte la densidad electrónica está dada por la ecuación (2.15).
Para que estas ecuaciones sean empleadas para describir un sistema interactuante, Kohn
y Sham hicieron una serie de modificaciones en el potencial externo a considerar [10].
Para desarrollar claramente esta idea, partimos de la definición de la funcional universal
F [ρ]:
F [ρ] = T [ρ] + Vee (2.16)
ya que T [ρ] �= Tnirs[ρ], Kohn y Sham propusieron reescribir el funcional F [ρ] como:
F [ρ] = Tnirs[ρ] + J [ρ] + Exc[ρ] (2.17)
donde Tnirs[ρ] es la funcional de enerǵıa cinética del sistema no-interactuante, J [ρ] es
el término clásico de interacción electrón-electrón y Exc[ρ] es el el término que absorbe
la interacción no clásica entre las part́ıculas. A éste se le llama término de intercambio
y correlación y se define como:1
1 En algunos esquemas de primeros principios, la enerǵıa de correlación se define como la diferencia
entre la enerǵıa no-relativista exacta y la enerǵıa calculada mediante el método de Hartree-Fock [11].
La correlación se da por la superposición espacial de los orbitales que caracterizan a las part́ıculas,
2.2. CONSTRUCCIÓN DE LA TEORÍA FUNCIONAL DE LA DENSIDAD. 11
Exc ≡ T [ρ] − Tnirs[ρ] + Vee[ρ] − J [ρ] (2.18)
Obsérvese que este término contempla la diferencia T [ρ] − Tnirs[ρ], entre el funcio-
nal de enerǵıa cinética exacto del sistema interactuante y el funcional correspondiente
al sistema de referencia no-interactuante. Además, contempla también la diferencia
Vee[ρ]−J [ρ], que contribuye con un término no clásico de interacción electrón-electrón.
Ahora el funcional de enerǵıa Kohn-Sham lo podemos escribir como:
E[ρ] =
∫
ρ(r) vnuc(r) d
3r + Tnirs[ρ] + J [ρ] + Exc[ρ] (2.19)
donde el primer término y el tercer término del lado derecho de la ecuación correspon-
den a los términos de las ecuaciones 2.10 y 2.11 respectivamente. Para determinar la
densidad electrónica ρ(r) que minimice a la enerǵıa debemos de imponer el principio
de mı́nima enerǵıa:
δE[ρ]
δρ(r)
= 0 (2.20)
Sin embargo, como la densidad electrónica está descrita en términos de los orbitales
moleculares, la ecuación de Euler correspondiente es:[
−1
2
∇2 + veff (r)
]
φi = εi φi (2.21)
en esta última ecuación, veff (r) es el potencial efectivo de Kohn y Sham definido como
2
:
veff (r) = vnuc(r) +
δJ [ρ]
δρ(r)
+ vxc(r) (2.22)
y la contribución vxc es el potencial de intercambio y correlación que depende de la
densidad electrónia ρ(r).
vxc(r) =
δExc[ρ]
δρ(r)
(2.23)
La ecuación (2.21) es la ecuación de Schrödinger de una part́ıcula no interactuante
dentro del potencial veff . Si consideramos que veff puede jugar el papel de un poten-
cial externo para la part́ıcula en cuestión, entonces habremos alcanzado el mapeo que
nos permite describir a un sistema interactuante en términos de uno sin interacción.
En otras palabras, el mapeo de un sistema interactuante a un sistema sin interacción,
propuesto por Kohn y Sham [10], se logra modulando el potencial externo vext de tal
lo que provoca que el comportamiento de una part́ıcula se de con base al comportamiento de las
otras. Mientras que el fenómeno de intercambio es una consecuencia de la interacción esṕın-esṕın, y
f́ısicamente se refiere a la probabilidad de una part́ıcula de encontrarse en el orbital φi y en el orbital
φj [12].
2 El potencial veff es de tipo efectivo porque es el potencial producido por todas las part́ıculas
(electrones y núcleos) en conjunto. Es decir, el potencial veff no toma en cuenta a cada part́ıcula por
separado, sino el efecto conjunto de todas ellas.
12 CAPÍTULO 2. TEORÍA FUNCIONAL DE LA DENSIDAD
forma que incluya al potencial nuclear vnuc, al potencial coulombiano δJ [ρ]/δρ(r) y al
potencial de intercambio y correlación vxc(r). Bajo estas condiciones es posible usar el
formalismo “convencional” de DFT porque existe un potencial externo que gobierna al
hamiltoniano, aunque este potencial externo dependa de la densidad electrónica ρ(r) a
través de vxc(r), es decir, el esquema Kohn-Sham se fundamenta en el esquema Hohen-
berg y Kohn. También, mediante el mapeo antes mencionado hemos conseguido que la
ρnirs(r) del sistema no-interactuante corresponda precisamente a la ρis(r) del sistema
interactuante, debido a que la inequivalencia entre sistemas dada por la expĺıcita di-
ferencia de los potenciales externo vext y efectivo veff es considerada al construir los
orbitales, como puede verse en la ecuación (2.21).
Con el método Kohn-Sham no requerimos conocer la funcional de enerǵıa cinética
porque ahora lo podemos escribir de manera exacta en términos de orbitales molecu-
lares3 . Esto aumenta la precisión de los cálculos, en contraste con los funcionales en
términos de la densidad electrónica ρ(r) porque éstos se desconocen [14]. No obstante,
este método aún demanda la forma expĺıcita de la funcional de intercambio y correla-
ción Exc[ρ(r)]. El que desconozcamos la forma exacta de Exc[ρ] no es un problema grave
para la precisión de los resultados porque siempre se puede construir una funcional de
intercambio y correlación relativamente precisa.
Ahora tenemos que resolver n ecuaciones del tipo:[
−1
2
∇2 + veff (r)
]
φi = εiφi
donde:
veff (r) =
∑I
ZI
| r−RI | +
∫
ρ(ŕ)
| r− ŕ |dŕ+ vxc(r)
(una para cada orbital) en forma auto-consistente, en vez de una única ecuación para
la densidad electrónica dentro del marco estricto de DFT establecido por Hohenberg y
Kohn.
Para resolver las ecuaciones dentro del esquema Kohn-Sham, básicamente se propone
una densidad electrónica de prueba ρ̃(r) , que posteriormente empleamos para definir
el potencial efectivo dado por la ecuación (2.22). Una vez establecido el potencial se
procede a resolver las n ecuaciones de Schrödinger de una part́ıcula. Mediante este
proceso obtenemos n nuevos orbitales que usamos para construir una nueva densidad
electrónica. El proceso se repite hasta alcanzar la convergencia entre las densidades
electrónicas de prueba y la resultante. Al final la enerǵıa exacta es:
ETotal =
∑
ocupados
εi
ρ(r) =
∑
ocupados
|φi(r)|2 (2.24)
3 En 1990 March y Santamaŕıa encontraron una expresión que relaciona la funcional de la enerǵıa
cinética con la funcional de intercambio[13]
2.3. ENERGÍA DE INTERCAMBIO Y CORRELACIÓN 13
2.3. Enerǵıa de intercambio y correlación
Como se mencionó anteriormente, la expresión para Exc[ρ] es desconocida. Por tanto
el desaf́ıo consiste en construir una funcional de intercambio y correlación (Exc[ρ]). En
la actualidad existen excelentes aproximaciones para Exc[ρ] que han permitido a DFT
(Teoria Funcional de la Densidad) ocupar un lugar privilegiado dentro de los métodos
de aproximación ab-initio.
2.3.1. Aproximación de densidad local (LDA)
La primera aproximación para Exc[ρ], consiste en considerar un gas uniforme de
electrones. La expresión es la siguiente:
ELDAxc [ρ] =
∫
ρ(r)�xc[ρ]dr (2.25)
donde �xc[ρ] es la enerǵıa de intercambio y correlación por part́ıcula de un gas uniforme
de electrones. La función �xc[ρ] puede dividirse en una contribución de intercambio y
otra de correlación,
�xc(ρ) = �x(ρ) + �c(ρ) (2.26)
La parte de intercambio ya es conocida, dada por la funcional de intercambio de
Dirac [8];
�x(ρ) = −Cx�x(ρ) 13 (2.27)
donde Cx =
3
4
( 3
π
)
1
3 .
La parte de �c(ρ), se desarrollo gracias a cálculos de montecarlo cuánticos de Cerpe-
ley y Alder (1980). Estos valores fueron interpolados para proveer una forma anaĺıtica
para �c(ρ) (Vosko, Will and Nusair 1980) [8] (Para mas detalles de las funcionales ver
apéndice A.).
Las caracteŕısticas de LDA son [15, 16, 17]:
• Favorece a los sistemas más homogéneos.
• Sobre enlaza a moléculas y sólidos.
• Para puentes de hidrógeno, fuerzas de Van der Waals las distancias de enlace ob-
tenidas son cortas.
A pesar de estos inconvenientes es normalmente muy dif́ıcil mejorar los resultados
obtenidos para las enerǵıas de adsorción de hidrógeno molecular en sistemas carbonosos.
14 CAPÍTULO 2. TEORÍA FUNCIONAL DE LA DENSIDAD
Figura 2.1: Comparación de la enerǵıa de enlace de la molécula de H2 con distintos
funcionales. PZ (LDA), PW91 (GGA), LYP (GGA), PBE (GGA).
2.3.2. Aproximación de Gradiente Generalizado (GGA)
La aproximación GGA, expresa la enerǵıa de intercambio y correlación en forma
local, en términos de la densidad y del gradiente de la densidad.
Las caracteŕısticas de GGA son [8, 18] :
• Las constantes de red se incrementan comparadas con LDA.
• Las enerǵıas cohesivas decrecen comparadas con LDA.
• Subestima las enerǵıas de enlace de hidrógeno molecular en superficies de carbono.
• Mejora las enerǵıas de atomización moleculares.
Para saber cual de las funcionales de intercambio y correlación es la más adecuada
para predecir la enerǵıa de enlace en superficies de carbono, se realizarón cálculos, y se
construyeron las curvas de enerǵıa de potencial. La gráfica 2.1 muestra que el funcional
LDA predice una enerǵıa de enlace de -0.082 eV, PW91 predice enerǵıa de enlace de
-0.018 eV, PBE de -0.009 eV y LYP no predice adsorción. La adsorción es un proceso
por el cual átomos, iones o moléculas son enlazadas o retenidas en la superficie de un
material. El resultado experimental de la enerǵıa de adsorción de H2 en grafito es de
0.060 eV [19]. Por tanto el resultado que mejor se ajusta, es con LDA. Para conocer el
tipo de enlace de moléculas con superficies, aśı como la reactividad de las superficies,
2.4. DENSIDAD DE ESTADOS 15
se utilizará la densidad de estados y el análisis de Löwdin.
2.4. Densidad de estados
La densidad de estados del sistema puede calcularse a partir del espectro de autova-
lores de las ecuaciones de Kohn-Sham. Los autovalores de las ecuaciones de Kohn-Sham
no tienen significado f́ısico; sin embargo se utilizan como una aproximación a la estruc-
tura de bandas. La densidad de estados de un sistema se puede definir como:
g(�) =
∑
n
∫
BZ
1
4π3
δ(� − �n(k))dk (2.28)
donde �n(k) representa la enerǵıa de la banda n en el punto k. La suma se realiza
sobre todas las bandas. Para el calculo numérico dicha integral se transforma a una
suma sobre un número finitos de puntos k. Para estimar la densidad de estados total,
se ensancha cada uno de los valores obtenidos mediante una gaussiana de un ancho
adecuado σ:
g(�) =
∑
n
∑
k
wk
1
σ
√
2π
exp
(
� − �n(k)
2σ2
)
(2.29)
donde wk es el peso asignado para cada punto k de acuerdo al muestreo del espacio k.
La densidad de estados nos da información acerca de la reactividad de la superficies
[20].
2.5. Análisis poblacional de Löwdin
La densidad de carga electrónica en términos de los orbitales está dada por:
ρ(r) = 2
n/2∑
i=1
|φi(r)|2 (2.30)
donde n es el número de electrones y hay dos electrones por cada orbital molecular.
La condición de normalización implica:∫
ρ(r)dr = 2
n/2∑
i=1
∫
|φi(r)|2dr = 2
n/2∑
i=1
1 = N (2.31)
Si ahora remplazamos por los orbitales atómicos para saber como es la forma de la
densidad electrónica en la base utilizada en los cálculos:
ρ(r) = 2
n/2∑
i=1
∑
ν
C∗νiϕ
∗
ν(r)
∑
μ
Cμiϕμ(r) =
∑
μν
⎛
⎝2 n/2∑
i=1
C∗μiνi
⎞
⎠ϕμϕ∗ν = ∑
μν
Pμνϕμϕ
∗
ν
(2.32)
16 CAPÍTULO 2. TEORÍA FUNCIONAL DE LA DENSIDAD
integrando la ecuación 2.32, tenemos que esta integral es igual al número de átomos:
N =
∑
μν
Pμν
∫
ϕ∗ν(r)ϕμ(r)dr =
∑
μν
PμνSνμ =
∑
μ
(PS)μμ = tr(PS) (2.33)
De este modo es posible asignar cuantos electrones hay por cada orbital atómico
[21]. El elemento de la diagonal del producto (PS)μμ es el número de electrones en el
orbital atómico ϕμ.
Se puede notar, que las definiciones en el análisis poblacional depende de la elección
de las bases atómicas utilizadas. Como por ejemplo:
N =
∑
μ
(S1−αSαP)μμ =
∑
μ
(S1−αPSα)μμ (2.34)
donde se utiliza la propiedad tr(AB)=tr(BA). Por tanto α representa la libertad que
tengo en cuanto una elección de base. En particular si elegimos α = 1/2, obtenemos una
matriz de densidad de carga en la base de orbitales atómicos diagonalizados simétrica-
mente.
2.6. Base de ondas planas
Los orbitales Kohn-Sham son representados por una base de ondas planas,
φi,k(r) =
∑
G,1
2
|G+k|2≤Ecut
ci,G+ke
i((G+k)r) (2.35)
truncada hasta la enerǵıa de corte Ecut. En la práctica se utiliza una enerǵıa de
corte los suficientemente grande para lograr la convergencia en la enerǵıa total [22].
Con esta base las ecuaciones de Kohn-Sham toman la forma
Ecut∑
G
〈Ǵ+ k|FKS|G+ k〉 = �i,G+kci,G+k (2.36)
Para cualquiera de los potenciales locales la expresión para los elementos de matriz es:
Ecut∑
G
〈Ǵ+ k|V (r)|G+ k〉 = V (Ǵ−G) (2.37)
donde V (Ǵ−G) es la transformada de Fourier del potencial valuada en la diferen-
cia entre los vectores G y Ǵ; por lo tanto todos los elementos de matriz (salvo los del
pseudo-potencial no local) pueden calcularse eficientemente a partir de las transforma-
das de Fourier de los potenciales.
Caṕıtulo 3
Pseudo-potenciales
3.0.1. Definición de pseudo-potencial
Cuando los orbitales de Kohn-Sham de un sistema molecular o de una fase conden-
sada se expresan en una base, se utiliza una gran cantidad de funciones, para desarrollar
los estados del carozo atómico formado por los electrones más internosdel átomo. Estos
estados no contribuyen al enlace y permanecen prácticamente inalterados en los distin-
tos entornos qúımicos en los que se puede encontrar el elemento. Por lo tanto es útil
considerar expĺıcitamente sólo los electrones de valencia y considerar la interacción con
el núcleo y los electrones internos como parte de un potencial externo [23]. Esta aproxi-
mación se conoce como aproximación de pseudopotenciales. Los electrones internos no
intervienen en el enlace qúımico sino por medio de la interacción con los electrones de
valencia. Dado que los estados atómicos son ortogonales entre los estados de carozo y
valencia estos últimos tienen fuertes oscilaciones cerca del núcleo, como se muestra en
la figura 3.1.
La gran cantidad de oscilaciones imponen una prueba muy dura a la flexibilidad del
conjunto de funciones base utilizado para desarrollar estos orbitales; se necesitan una
gran cantidad de funciones base para reproducir las oscilaciones y esto es prácticamente
imposible cuando se utilizan ondas planas. Una de las caracteŕısticas es exigir al pseudo-
potencial que no reproduzca la interacción carozo valencia exactamente cerca del núcleo
sino que la pseudofunción de valencia (la autofunción del pseudo-potencial) sea suave
en esta región.
Planteemos las ecuaciones de Konh-Sham para el átomo:
[−1
2
d2
dr2
+
l(l + 1)
2r2
+ V [n](r)]Rnl(r) = εnlRnl(r). (3.1)
Donde V[n](r) es el potencial efectivo:
V [n](r) = −Z
r
+ VH [n](r) + Vxc[n](r), (3.2)
17
18 CAPÍTULO 3. PSEUDO-POTENCIALES
Figura 3.1: Estados atómicos para el ox́ıgeno.
que hemos supuesto central para poder separar el orbital en una parte radial y una
parte angular de la forma:
ψnlm = Rnl(r)Ylm(θ, φ)
.
Pretendemos encontrar un pseudo-potencial V PP [n](r), tal que al tomar el lugar del
potencial efectivo total reproduzca el potencial que sienten los electrones de valencia.
Por otra parte, podŕıamos requerir que las funciones propias Rnl(r) del pseudo-potencial
tengan alguna caracteŕıstica particular, por ejemplo, que no tengan nodos en la región
del carozo, de forma tal que faciliten los cálculos en situaciones más complejas. Para
que este pseudo-potencial sea útil debe ser posible aplicarlo en distintos ambientes
qúımicos, eso significa que el pseudo-potencial sea transferible1 . Al ser aplicado a un
entorno distinto al atómico la forma más general en la que podemos representar este
pseudo-potencial es la de un operador no local V (r, r′) que actua sobre la función ψ(r)
para dar Ψ(r) de la forma:
Ψ(r) =
∫
V(r, r′)ψ(r′)dr′
Estudiemos más detalladamente las caracteŕısticas de este operador. Tal como lo
1 Por ejemplo, si ya se tiene el pseudo-potencial ajustado a las propiedades experimentales atómicas
deberá ser posible utilizarlo para el cálculo de moléculas del mismo elemento o para un cristal, a esto
se le conoce como transferabilidad.
19
planteamos previamente, esperamos que este pseudo-potencial aśı como el potencial
atómico efectivo tenga simetŕıa esférica, por lo tanto V (r, r′) es función sólo de los
módulos r y r’y del ángulo θ entre ellos V (r, r′, cos(θ)). Podemos desarrollar el potencial
en cos(θ) en una serie de polinomios de Legendre2 :
V (r, r′) =
∞∑
l=0
Vl(r, r
′)Pl(cos(θ)), (3.3)
donde:
Vl(r, r
′) = Vl(r)δ(r− r′)
y por tanto la expresión completa es:
V (r, r′) =
∞∑
l=0
l∑
m=−l
Y∗lm(θ, φ)Vl(r)Ylm(θ
′, φ′). (3.4)
Esta forma del pseudo-potencial es semilocal, local en r y no local en θ y φ. Aunque
en forma general puede ser no local cuando:
Vl(r, r
′) = Vloc(r)δ(r− r′) +ΔVl(r, r′).
Esperamos que las expresiones anteriores, tanto la forma semilocal como la forma
no local, lejos del núcleo se parezcan al potencial electrostático generado por el núcleo
y los electrones del carozo; de acuerdo al teorema de Gauss este potencial será el de una
carga puntual igual a la carga de valencia fuera del carozo. Este potencial electrostático
debe de ser local y por lo tanto todos los componentes Vl(r, r
′) deben de tender asintóti-
camente a Zval
r
y se desviarán cerca del núcleo debido a la no localidad inducida por la
ortogonalidad con los estados del carozo. Resumiendo
ĺım
r→∞
Vl(r) = −Zval
r
,
para todo l.
Podemos dividir el potencial de cada término de la forma:
Vl(r) = Vloc(r) + ΔVl(r),
donde Vloc(r) es una función arbitraria que tiene el comportamiento asintótico −Zvalr
y ΔVl(r) = Vl(r) − Vloc(r) . Hemos dividido al potencial V (r, r′) en general en una
parte no local de corto alcance y una parte local arbitraria de largo alcance que tiene
el comportamiento asintótico −Zval
r
2 Los polinomios de Legendre forman un conjunto completo para desarrollar una función f(x) en el
intervalo [-1,1]
20 CAPÍTULO 3. PSEUDO-POTENCIALES
V (r, r′) = Vloc(r) +
∞∑
l=0
l∑
m=−l
Y∗lm(θ, φ)ΔVl(r)Ylm(θ
′, φ′).
Para l grande en la ecuación (3.1) el término que domina a r pequeño es el potencial
centŕıfugo l(l+1)
2r2
, por lo tanto, para algún l = lmax el potencial ΔVl(r) de corto alcance
será despreciable con respecto al potencial centŕıfugo y lo podemos considerar nulo para
l ≥ lmax. De esta forma podemos truncar la serie y la expresión para el pseudo-potencial
queda como sigue:
V (r, r′) = Vloc(r) +
lmax∑
l=0
l∑
m=−l
Y∗lm(θ, φ)ΔVl(r)Ylm(θ
′, φ′) (3.5)
En la práctica, normalmente la serie se trunca para l = 2 o l = 3 dependiendo del
átomo.
3.0.2. Construcción de pseudo-potenciales
Existe una gran libertad para la construcción de pseudo-potenciales, éste debe de
reproducir las propiedades dispersivas del núcleo y los electrones del carozo fuera de
éste: pero puede tener una forma relativamente arbitraria dentro del carozo siempre y
cuando cumpla con las siguientes condiciones:
• Las pseudo-funciones de onda de valencia generadas por el pseudo-potencial no
deben contener nodos. Esto asegura que las pseudo-funciones sean suaves y por lo tanto
requieran menos funciones base para su desarrollo.
• La pseudofunción y la función propia del potencial total deben de ser idénticas a
partir de un determinado radio de corte rcl:
RPPl (r) = Rl(r)
para r ≥ rcl
Esta condición asegura que el pseudo-potencial reproduzca la densidad electrónica
de valencia y sus propiedades fuera del carozo.
• La carga contenida dentro del radio del corte de la pseudofunción y de la auto-
función deben de ser iguales:∫ rcl
0
|RPPl (r)|2dr =
∫ rcl
0
|Rl(r)|2dr
Esto asegura el comportamiento asintótico −Zval
r
por el teorema de Gauss y asegura
ciertas propiedades dispersivas.
• Los autovalores de la pseudofunción y la autofunción deben de ser idénticos:
�PPl = �l
21
Los pseudo-potenciales que cumplen con estas condiciones se denominan “pseudo-
potenciales que conservan la norma” [24]. Un ejemplo de pseudo-potenciales de con-
servacón de norma son los de Troullier-Martins y Hamann [25] . También existen los
pseudo-potenciales ultrasuaves que no conservan la norma [26].
3.0.3. Pseudo-potenciales Troullier-Martins
La pseudofunción Troullier-Martins tiene la siguiente forma:
RPPl (r) =
⎧⎨
⎩
RAEnl (r) if r > rcl;
rlep(r) if r < rcl.
donde,
p(r) = c0 + c2r
2 + c4r
4 + c6r
6 + c8r
8 + c10r
10 + c12r
12
y RAEnl (r) es la autofuncion. Los coeficientes ci se ajustan de las condiciones de con-
tinuidad y de la curvatura cero del pseudo-potencial en el origen. En la figura 3.2 se
muestra las pseudofunciones comparadas con las autofunciones de onda para el átomo
de aluminio. Se observa que después de un cierto radio de corte (rc), las pseudofunciones
y las autofunciones son iguales.
Figura 3.2: Pseudo-funciones Troullier-Martins vs autofunciones para el átomo de
aluminio.
A las pseudo-funciones mostradas en la figura 3.2, les corresponde los pseudo-
potenciales mostrados en la figura 3.3.
22 CAPÍTULO 3. PSEUDO-POTENCIALES
Figura 3.3: Pseudo-potenciales Troullier-Martins para el aluminio. Este pseudo-
potencial fue generado con el programa fhi98PP.
Los pseudo-potenciales mostrados en la figura 3.3 son obtenidosrestando al pseu-
dopotencial apantallado (V PP,scrl (r)) el término electrostático (V
H [ρPP ; r]) y el término
de intercambio y correlación (V XC [ρPP ; r]),
V PPl = V
PP,scr
l (r)− VH [ρPP ; r]− VXC [ρPP ; r], (3.6)
con
ρPP (r) =
1
4π
∑
l=0
fl
∣∣∣∣RPPlr
∣∣∣∣
2
, donde
V PP,scrl (r) = �l −
l(l + 1)
2r2
+
1
2RPPl (r)
d2
dr2
RPPl (r)
Este potencial, se convierte en el potencial del átomo con todos sus electrones para
r > rcl. El término V
PP
l = ΔVl(r) se sustituye en la la ecuación 3.5.
3.0.4. Pseudo-potenciales ultrasuaves-Vanderbilt
Como su nombre lo sugiere, los pseudo-potenciales ultrasuaves son mucho más suaves
que los Troullier-Martins. Por tanto, se requiere mucho menos ondas planas para los
23
cálculos. En este esquema la densidad de carga de valencia n(r) esta formada por dos
contribuciones,
n(r) =
∑
n
[
|φPPl (r)|2 +
∑
ij
Qij(r)〈φPPl |βj〉〈βi|φPPl 〉
]
donde βi son funciones de proyección que dependen de las posiciones ionicas y de la
función Qij(r) que está dada por
Qij = R
∗AE
nl (r)R
AE
nl (r)− φi ∗ (r)φj(r)
, RAEnl son las autofunciones del átomo, y φj(r) son las pseudofunciones ultrasuaves
construidas sin cumplir con la condición de norma Qij = 0. Las pseudo-potenciales
ultrasuaves son obtenidos del mismo modo que los Troullier-Martins, sustituyendo la
pseudofuncion ultrasuave en la ecuación 3.6. Los pseudo-potenciales Troullier-Martins
se usarán en los caṕıtulos 5 y 8, mientras que los ultrasuaves-Vanderbilt se usarán en
los caṕıtulos 6 y 7.
Caṕıtulo 4
Dinámica Molecular
La dinámica molecular (MD) es una técnica que nos permite estudiar las propie-
dades de sistemas de muchas part́ıculas tanto en equilibrio como fuera de equilibrio.
La información que genera una corrida de MD es la posición y la velocidad de cada
part́ıcula del sistema en cada instante de tiempo, lo cual significa obtener la trayecto-
ria de un punto del espacio fase en función del tiempo. Suponiendo que el sistema es
ergódico, podemos asociar directamente el promedio temporal de cierta observable con
el promedio usual de un ensamble de la mecánica estad́ıstica. La clave en la aproxi-
mación de MD es: la separación del movimiento lento de los núcleos, del movimiento
rápido de los electrones con la aproximación de Born-Oppenheimer. El movimiento de
los núcleos son gobernados por las ecuaciones de Newton,
MI
d2
dt2
RI = − ∂
∂RI
E0({RI}) (4.1)
donde MI y RI son las masas y las coordenadas de N átomos, respectivamente, y
E0({RI}) es la enerǵıa total del estado base de la ecuación (2.5).
La enerǵıa del estado base (E0({RI})) y la densidad electrónica son obtenidas re-
solviendo autoconsistentemente las ecuaciones de Kohn-Sham como se menciono en el
capitulo 2.2. La enerǵıa total es:
E[ρ] =
∫
ρ(r) vnuc(r) d
3r + Tnirs[ρ] + J [ρ] + Exc[ρ] + E
nuc−nuc. (4.2)
Donde Tnirs[ρ] es la enerǵıa cinética de electrones no interactuantes, J [ρ] es la in-
teracción clásica electrón-electrón y Exc[ρ] la enerǵıa de intercambio y correlación. La
enerǵıa electrón-núcleo y núcleo-núcleo son:
Ee−nuc[ρ] =
∫
ρ(r) vnuc(r) d
3r (4.3)
Enuc−nuc[ρ] =
1
2
M,M∑
I,J=1,I �=J
e2oZIZJ
| RI −RJ | (4.4)
25
26 CAPÍTULO 4. DINÁMICA MOLECULAR
donde ZI y ZJ son las cargas correspondientes de los núcleos.
El sistema es representado por una supercelda periódica. Las coordenadas Rj de un
núcleo o su imagen periódica es:
RJ = τIs,Ia +R (4.5)
con J = Is, Ia,R, donde el ı́ndice Is es para las especies de un núcleo, Ia es para el
átomo mismo, τIs,Ia es un vector que va del núcleo (Is, Ia) al origen de la supercelda y
R es un vector que va del sistema de referencia al origen de la supercelda.
El efecto de los electrones del carozo y el potencial de Coulomb es remplazado por
pseudo-potenciales suaves que permite la eficiencia del uso de una base de ondas planas,
V̂e−nuc =
∑
R
∑
Is,Ia
V̂e−nucIs (r − τIs,Ia −R) (4.6)
Se emplea pseudo-potenciales construidos siguiendo los esquemas de Troullier- Mar-
tins [25] y ultrasuaves [26]. En algunos elementos, como el litio, es esencial incluir la
corrección del carozo no lineal, en la enerǵıa de intercambio y correlación (NLCV-XC).
En este caso, el efecto de la densidad electrónica del carozo en la enerǵıa de intercam-
bio y correlación, es descrita por una densidad de pseudo-carozo ñcore, el cual es una
superposición de densidades de carozo suavizadas ñcore(|r− τIs,Ia −R|), construidas
junto con los pseudo-potenciales. Los funcionales de intercambio y correlación Exc[n] y
Vxc[n] entonces son remplazados por Exc[n+ ñcore] y Vxc[n+ ñcore].
Una vez que la E0({RI}) se ha calculado, las fuerzas (FIs,Ia) que actúan entre los
núcleos, son obtenidas a partir de la ecuación 4.1, usando el algoritmo de Verlet [27]
τIs,Ia(t+ δtnuc) = τIs,Ia(t)− τIs,Ia(t − δtnuc) + δt2nuc
FIs,Ia(τIs,Ia(t))
MIs
(4.7)
donde MIs es la masa del núcleo y δtnuc es el paso de tiempo. Una buena selección para
δtnuc es
1
15
del periodo del fonon.
Las velocidades vIs,Ia de los núcleos son calculadas por
vIs,Ia(t) = 3τIs,Ia(t)− 4τIs,Ia(t − δtnuc) + τIs,Ia(t − 2δtnuc) (4.8)
la cual es aproximada a segundo orden en δt. Se nota que las velocidades no son direc-
tamente integradas por el algoritmo de Verlet, pero se obtienen de su trayectoria.
Para encontrar la estructura optimizada de un sistema, se ocupa la dinámica de
newton amortiguada. Se empieza de una configuración inicial para los núcleos que se
mueven de acuerdo al esquema iterativo
τ
(nit+1)
Is,Ia
= (1 + λIs)τIs,Ia
(nit) − λIsτ (nit−1)Is,Ia + μIsFIs,Ia(τ (nit)Is,Ia ) (4.9)
donde λIs y μIs son los parámetros de amortiguamiento y la masa reciproca respectiva-
mente [28]. El parámetro nit son los pasos de la relajación del sistema. En los materiales
ha estudiar, primero optimizaremos la estructura y posteriormente para introducir los
efectos de la temperatura y la presión, usaremos la dinámica molecular.
27
4.0.5. Termostatos
Para desarrollar simulaciones en el ensamble canónico (N V T), se emplean termos-
tatos para el control de temperatura. Los efectos del calentamiento o enfriamiento en
la simulación fueron incorporados con rescalamientos en las velocidades.
La distribución inicial de velocidades vIs,Ia son dadas de una distribución Maxwell-
Boltzman para una temperatura dada.
La enerǵıa total del sistema esta dada por:
E =
N∑
Is=1
MIsv
2
Is,Ia
2
+
N∑
i>j=1
V (rij) (4.10)
De modo que las ecuaciones de movimiento se ven modificadas de la siguiente ma-
nera:
MIs
d2
dt2
RI = FTotal (4.11)
FTotal = FIs,Ia −
MIsdvIs,Ia
dt
(4.12)
El último término de la ecuación 4.12 le introduce una fuerza extra que acelera o
desacelera los núcleos. Para obtener la temperatura deseada del sistema, cada velocidad
atómica es escalada por un factor, después de un número determinado de pasos de
tiempo de la dinámica.
v́Is,Ia = λvIs,Ia (4.13)
donde
λ =
√
Td/Ta
, Ta es la temperatura actual y Td es la temperatura aplicada. La temperatura actual
se obtiene a partir de la equipartición de la enerǵıa cinética de la ecuación 4.10,
Ta =
2EK
GKB
, donde EK = enerǵıa cinética y G = 3N −6 total de grados de libertad independientes.
4.0.6. Barostatos
Para simular un sistema a presión constante, se usan los barostatos para el control
de la presión. Parrinello y Rahman propusieron un barostato en el cual se permite el
cambio en el volumen y en la supercelda [29]. Ellos usaron como variables dinámicas
las componentes cartesianas
hij = eij
28 CAPÍTULO 4. DINÁMICA MOLECULAR
de los tres vectores aj, que definen la periodicidad de la supercelda. Aqúı ei son los tres
vectores ortonormales que definen el sistema cartesiano. Para realizar la dinámica, se
incluye en el langragiano una enerǵıa cinética ficticia de la celda
KPRcell =
W PR
2
3∑
i=1
3∑
j=1
(ḣij)
2 (4.14)
donde W PR es una masa ficticia. En el limite de N grande, el principio de equipartición
nos dice que la enerǵıa cinética de las 9 variables de la celda son pequeñascomparadas
con la enerǵıa cinética de las 3N−6 variables asociadas a las posiciones de las part́ıculas,
de modo que el método simula un ensamble isobárico-isoentalpico.
Caṕıtulo 5
Adsorción de agua en
grafeno-titanio.
La interacción de agua con superficies grafiticas ha sido estudiada minuciosamente
por diversos grupos [30, 31, 32, 33, 34]. Recientemente, el sistema de una hoja de grafeno
decorada con Ti, fue estudiado teóricamente, exponiendo dicho sistema a moléculas de
hidrógeno y aire [35]. Se encontró que el ox́ıgeno es una especie muy reactiva; oxidando a
los átomos de Ti, aún después de que las moléculas de H2 ya estaban enlazadas al átomo
de Ti. Las moléculas de nitrógeno y agua también son quimisorbidas, sin disociación.
El sistema de grafeno− Ti considerado en la referencia [35] fue estudiado con un baja
densidad de Ti, despreciando la interacción Ti − Ti. Por otro lado, las técnicas que se
han desarrollado para contaminar la superficie de grafito con titanio, también podŕıan
servir para contaminar a las hojas de grafeno [36]. En este trabajo, proponemos a las
hojas de grafeno contaminadas con titanio como catalizadores, para la obtención de
hidrógeno a partir del agua.
Sin embargo, en la práctica, podŕıa ser dif́ıcil dispersar a los átomos de Ti en las hojas
de grafeno, ya que los átomos de Ti podŕıan segregarse para formar islas, del mismo
modo que le sucede en el grafito [37]. Zhang encontró que el Al y Au depositados
en nanotubos migran y para formar nanopart́ıculas [38]. También encontróó que los
nanotubos con un pretratamiento de moléculas surfactantes, posibilitan la uniformidad
de Al y de Au depositados en los nanotubos. Este hecho sugiere que este procedimiento
podŕıa también facilitar la dispersión de átomos de Ti en la superficie de grafeno.
En esta sección se estudiará la adsorción de una molécula de agua en un sistema de
grafeno, decorado con una gran concentración de Ti1 . Se consideró dos configuraciones
estables para el sistema grafeno-titanio: la primera configuración contiene ocho átomos
de C por uno de Ti (C8Ti); la segunda contiene un átomo de Ti por dos de C (C2Ti).
Esta última configuración contiene la máxima concentración posible de Ti; hay un
átomo de Ti sobre el centro de cada hexágono de la hoja de grafeno.
1 Eduardo Rangel, Gregorio Ruiz-Chavarria, L.F. Magana.“Water molecule adsorption
on a titanium-graphene system with high metal coverage”. Carbon 47, 2009, págs: 531−533.
29
30 CAPÍTULO 5. ADSORCIÓN DE AGUA EN GRAFENO-TITANIO.
Para representar el sistema C2Ti se usa una supercelda hexagonal con a = b = 4,875
Ȧ, con condiciones periódicas, que contiene cuatro átomos de Ti por ocho de C (ver
Figura 5.1-a). Para el sistema C8Ti se usa la misma unidad de celda, que contiene un
átomo de Ti por ocho de C (ver Figura 5.1-b). En ambos casos, para estar seguros de
tratar con una hoja de grafeno aislada, el parámetro c de la celda hexagonal se tomó de
20 Ȧ; con esta distancia se anula la interacción entre las hojas de grafeno − titanio
adyacentes.
Figura 5.1: Configuración de átomos de Ti en una hoja de grafeno. a) Celda unidad
que contiene cuatro átomos de Ti por ocho de C (C2Ti). b) Celda unidad que
contiene un átomo de Ti por ocho de C (C8Ti)
5.1. Método
Se remplazó el carozo de todos los átomos por pseudopotenciales de conservación
de norma Troullier-Martins [25], en la forma completamente separable de Kleinman-
Bylander [24], los cuales fueron desarrollados con el código fhi98PP [39]. Se consideraron
los siguientes electrones de valencia: para Carbono 2s22p2 ; para Hidrógeno 1s2; para
Titanio 3s23p63d24s2; para Ox́ıgeno 2s22p4. Se ocuparon los siguientes radios de corte
(en Angstroms) para los pseudo-potenciales: Para H: rs =0.423 Ȧ, donde s es la com-
ponente local del pseudo-potencial; para C: rs =0.794 Ȧ, rp =0.815 Ȧ, y s como la
componente local; para O: rs = rp =0.661 Ȧ, con s como la componente local; y para
Ti: rs =0.807 Ȧ, rp =0.750 Ȧ, rd =1.257 Ȧ, con p como la componente local. No se
utilizó cálculos con esṕın polarizado y relativistas. La enerǵıa de corte fue de 1080 eV;
con esta enerǵıa se logra la convergencia en la enerǵıa total (sección 2.6). Se tomaron
34 puntos K en el esquema de Monkhorst-Pack [40]. El criterio de convergencia en la
enerǵıa fue de 10−4 eV. Para intercambio y correlación se utilizó la aproximación de
densidad local (LDA) con el funcional Perdew-Zunder [17].
5.2. DISCUSIÓN Y RESULTADOS 31
Usando los pseudo-potenciales de H y O se obtuvo por relajación una longitud
de enlace para H − O de 0.9578 Ȧ y para el ángulo H − O − H 104.47 (los valores
experimentales son 0.9578 y 104.50 respectivamente[41]). De la misma manera con el
pseudo-potencial de C, se obtuvo para la longitud de enlace qúımico en el grafito C−C
1.4076 Ȧ y para la distancia inter-planar 3.320 Ȧ (los valores experimentales son 1.415
y 3.350 Ȧ, respectivamente [42]). Y por último, usando el pseudo-potencial para Ti, se
obtuvieron los siguientes parámetros para la red cristalina: a =2.863 Ȧ y c =4.544 Ȧ, los
valores experimentales son 2.950 y 4.683 Ȧ, respectivamente [41]. Todos los cálculos de
DFT-MD fueron realizados con el programa Quantum-Espresso, que usa ondas planas
y pseudopotenciales [22].
5.2. Discusión y resultados
5.2.1. Grafeno-titanio
Primero se optimizan las posiciones de los átomos de Ti en la hoja de grafeno,
modelado con los sistemas: C2Ti y C8Ti, como se muestra en la figura 5.1-a y figura
5.1-b.
Se define la enerǵıa de enlace de los átomos de Ti como:
Eb = E(grafeno) + nE(Ti)− E(Grafeno+ nTi) (5.1)
Donde E(grafeno) + nE(Ti), es la enerǵıa del sistema inicial que corresponde a la
enerǵıa de la hoja de grafeno aislada y a la enerǵıa de los átomos de Ti, respectivamente.
El término E(Grafeno+nTi), es la enerǵıa para el sistema final y el valor n corresponde
al número de átomos de Ti. La enerǵıa de adsorción para C8Ti es -3.36 eV/Ti. La
distancia Ti− Ti es 4.87 Ȧ y la distancia Ti− grafeno es 1.67 Ȧ (Figura 5.1-b). Para
C2Ti, la enerǵıa de adsorción es -5.03 eV/Ti. En este caso los átomos de Ti están
contenidos en dos planos diferentes. El átomo más cercano al plano de grafeno está a
1.82 Ȧ. El otro plano de átomos de Ti, se encuentra a una distancia de 3.239 Ȧ del plano
de grafeno como se ilustra en la figura 5.1-a. De esta forma la distancia mas corta de
un átomo de Ti del plano inferior a un átomo de Ti en el plano superior es de 2.815 Ȧ.
La distancia de Ti−Ti es 4.87 y 2.44 Ȧ en el plano inferior y superior respectivamente
(Figura 5.1-a).
Es claro que la configuración C2Ti tiene una distribución hexagonal después de la
optimización; cada Ti esta sobre un hexágono de C (Figura 5.2-b). Por otro lado se
sabe que la superficie Ti (0 0 1) es hexágonal como se muestra en la figura 5.2-a, con
los átomos en un solo plano, y una distancia Ti − Ti de 2.95 Ȧ [41]. De este modo la
superficie C2Ti es parecida a la superficie Ti (0 0 1), pero no idéntica. Aśı que se espera
que la interacción de agua con C2Ti sea similar a la interacción con la superficie Ti
(0 0 1).
32 CAPÍTULO 5. ADSORCIÓN DE AGUA EN GRAFENO-TITANIO.
Figura 5.2: a) La superficie de Ti(0 0 1) tiene un arreglo hexagonal con todos los
átomos en un solo plano. b) La configuración de C2Ti, también tiene un arreglo
hexagonal pero los átomos de Ti en color obscuro, se encuentra en un plano
inferior con respecto a los átomos de Ti en color claro.
5.2.2. Naturaleza del enlace grafeno-titanio
En esta sección investigamos el tipo de enlace que se forma, de la interacción de una
hoja de grafeno con distintas densidades de titanio. Se consideró C2Ti, C8Ti, C18Ti
y C50Ti. Para las ultimas tres configuraciones, el átomo de Ti puede ser adsorbido en
tres diferentes sitios: sobre el centro del hexágono, sobre un átomo de carbono, y sobre
el enlace C −C. La posición mas estable es sobre el centro del hexágono y la distancia