Story Transcript
DEPARTAMENTO DE QUÍMICA FÍSICA FACULTAD DE CIENCIAS UNIVERSIDAD DE CÁDIZ
CÁLCULO TEÓRICO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTÁNDAR
DAVID ZORRILLA CUENCA
PUERTO REAL, JULIO DE 2001
DEPARTAMENTO DE QUÍMICA FÍSICA FACULTAD DE CIENCIAS UNIVERSIDAD DE CÁDIZ
CÁLCULO TEÓRICO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTÁNDAR
MEMORIA que presenta el Ldo.: D. DAVID ZORRILLA CUENCA para optar al grado de DOCTOR EN CIENCIAS QUÍMICAS dirigida por los doctores: D. MANUEL FERNÁNDEZ NÚÑEZ DÑA. ROSA RODRÍGUEZ HUERTAS
PUERTO REAL, JULIO DE 2001
Dña. ROSA RODRÍGUEZ HUERTAS, Catedrática de Escuela Universitaria del departamento de Estadística e Investigación Operativa de la Universidad de Cádiz y D. MANUEL FERNÁNDEZ NÚÑEZ, Catedrático del departamento de Química-Física de la Universidad de Cádiz, como sus directores
HACEN CONSTAR:
Que esta Memoria, titulada “Cálculo teórico de propiedades moleculares mediante bases no estándar”, presentada por D. David Zorrilla Cuenca, resume su trabajo de Tesis y, considerando que reúne todos los requisitos legales, autorizan su presentación y defensa para optar al grado de Doctor en Ciencias Químicas por la Universidad de Cádiz
Cádiz, Julio de 2001
Dra. Rosa Rodríguez Huertas
Dr. Manuel Fernández Núñez
AGRADECIMIENTOS
AGRADECIMIENTOS
AGRADECIMIENTOS La Tesis que se presenta a continuación es el fruto de seis años de trabajo, con sus alegrías y sus penas. Pienso que en estos seis años he sido influenciado en la forma de pensar por un gran número de personas. Muy a mi pesar, entre los nombres propios que a continuación debería citar, seguramente alguno caerá en el olvido. Por ello, antes de empezar desearía expresar mi mayor agradecimiento a todos, a los que apareceréis y a los que no. En primer lugar quiero expresar mi más profundo agradecimiento a mi Director de Tesis, D. Manuel Fernández Núñez, que no sólo ha ejercido de tutor, sino que en estos años hemos entablado una gran amistad y que en momentos familiares difíciles estuvo ahí para apoyarme y ayudarme incluso los “preciados” fines de semana. Recuerdo un 19 de Marzo, día del Padre, que estando trabajando en el programa informático de esta tesis me dijo que ya que no podía estar con sus hijos por estar estudiando fuera de Cádiz, estaba dichoso de poder estar con sus segundos hijos. Por todo esto y más, gracias. A Dña. Rosa Rodríguez Huertas, Codirectora de Tesis, que con sus sabios conocimientos matemáticos y estadísticos han influido notablemente en esta Tesis. A D. Jaime Fernández Rico y D. Rafael López por su aportación desinteresada de varias subrutinas utilizadas en el programa UCA-CMC. A D. Joaquín Martín Calleja, actual director del Departamento de Química-Física de la UCA, por sus esfuerzos por mejorar mi precaria situación laboral.
AGRADECIMIENTOS-i
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
A
mi
grupo
de
investigación
“Cálculo
Teórico
de
Propiedades
Moleculares” y a los doctorandos del grupo por sus consejos sobre la programación y diseño de los programas informáticos realizados para esta tesis. A la memoria de D. Julio Moral de Dios, por los buenos ratos que pasamos y por los buenos chistes que nos contaba en el laboratorio de Química Cuántica. A los compañeros de la Escuela Politécnica Superior de Algeciras, en especial a D. Juan Antonio Poce Fatou, por sus consejos y experiencia docente, y por los buenos ratos que hemos pasado en la Alameda y en los Carnavales. A los miembros del grupo de investigación “Síntesis, Caracterización y Evolución de Materiales” y a todo el Departamento de Química Física de la UCA por su apoyo y amistad en todos estos años. A mis Padres, por confiar en mí y dejarme hacer siempre lo que he considerado oportuno en materia de estudios. A mi hermano, familia y amigos, por apoyarme y querer siempre lo mejor para mí. A Yolanda, por su paciencia, perseverancia, ayuda y cariño que ha demostrado en todos estos años. También por sus buenos quehaceres culinarios, los cuales me han permitido no desfallecer en ningún momento. Te quiero.
GRACIAS A TODOS ...
AGRADECIMIENTOS-ii
ÍNDICE
ÍNDICE
ÍNDICE
AGRADECIMIENTOS .............................................................................I ÍNDICE .............................................................................................III INTRODUCCIÓN ............................................................................... VII 1. 2. 3.
Métodos Montecarlo y Cuasimontecarlo ....................................... vii Perspectivas del desarrollo informático ......................................... ix Descripción del contenido de este trabajo ......................................x
CAPÍTULO I: ANTECEDENTES .............................................................. 1 1. 2. 3. 4. 5. 6. 7. 8.
Introducción.............................................................................. 1 La Ecuación de Schrödinger......................................................... 2 Método Hartree-Fock .................................................................. 7 El método de campo Autoconsistente ..........................................15 Métodos post Hartree-Fock ........................................................26 Propiedades moleculares............................................................32 Integrales de los cálculos atómicos y moleculares .........................36 Integración mediante métodos MonteCarlo y Cuasimontecarlo........40
CAPÍTULO II: INTEGRACIÓN CON FUNCIONES DE BASE DE SLATER . 45 1. 2. 3. 4. 5.
Introducción.............................................................................45 Armónicos Esféricos ..................................................................47 Integrales Monocéntricas ...........................................................49 Integrales Policéntricas..............................................................56 Resultados Comparativos ...........................................................59
CAPÍTULO III: INTEGRACIÓN CON FUNCIONES DE BASE GAUSSIANAS ......................................................................................................... 65 1. 2.
Introducción.............................................................................65 Teorema del Producto (de Gaussianas) ........................................66
íNDICE-iii
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
3. 4. 5. 6. 7. 8. 9.
Integrales de Solapamiento........................................................ 70 Integrales de Energía Cinética .................................................... 72 Integrales de Momento.............................................................. 73 Integrales de Energía Potencial................................................... 75 Integrales Bielectrónicas............................................................ 80 Armónicos Esféricos y Armónicos Sólidos ..................................... 80 Resultados ............................................................................... 82
CAPÍTULO IV: INTEGRACIÓN MONTECARLO Y CUASIMONTECARLO CON FUNCIONES DE BASE STO (I) .................................................... 87 1. 2. 3. 4. 5. 6.
Introducción............................................................................. 87 Muestreo ................................................................................. 88 Integrales de Solapamiento........................................................ 96 Integrales Monoelectrónicas de Energías Cinética y Potencial ....... 102 Integrales Bielectrónicas.......................................................... 102 Uso de las Integraciones MC y CMC en Cálculos Químico Cuánticos104
CAPÍTULO V: INTEGRACIÓN MONTECARLO Y CUASIMONTECARLO CON FUNCIONES DE BASE STO (II) ........................................................ 105 1. 2. 3. 4.
Introducción........................................................................... 105 Cálculos Atómicos ................................................................... 106 Cálculos Moleculares ............................................................... 113 Problemas de Convergencia ..................................................... 116
CAPÍTULO VI: INTEGRACIÓN MONTECARLO Y CUASIMONTECARLO CON FUNCIONES DE BASE STO (III)............................................... 119 1. 2. 3. 4.
Introducción........................................................................... 119 Análisis de la Varianza: Fundamento Teórico .............................. 120 El modelo del análisis de la varianza con varios factores .............. 127 Estudio Estadístico de los errores inducidos por la INTEGRACIÓN aproximada............................................................................ 130
CAPÍTULO VII: DENSIDADES POLIGONALES................................... 149 1. 2. 3. 4. 5.
Introducción........................................................................... 149 Construcción de la función de peso “Poligonal+Exponencial” ........ 150 Muestreo de números aleatorios con función de peso Pn(r) ........... 151 Poligonales Angulares.............................................................. 153 Resultados ............................................................................. 154
CAPÍTULO VIII: CÁLCULOS CON FUNCIONES DE BASE ARBITRARIAS ....................................................................................................... 157 1.
Funciones de Base Arbitrarias................................................... 157
íNDICE-iv
ÍNDICE
2. 3. 4.
Funciones Espacialmente Restringidas. Muestreo de Integrales con Funciones de Base Recortadas .................................................. 161 Cálculos de Propiedades Atómicas en Cajas Esféricas con Funciones “STO-Recortadas” ................................................................... 168 Cálculo de Propiedades Moleculares con Funciones de Base Espacialmente Restringidas ...................................................... 174
CAPÍTULO IX: PROGRAMA UCA-CMC ............................................... 177 1. 2.
Introducción........................................................................... 177 UCA-CMC ............................................................................... 178
CAPÍTULO X: RESUMEN, CONCLUSIONES Y PERSPECTIVAS ............ 211 1. 2. 3.
Principales resultados alcanzados con esta Tesis ......................... 211 Conclusiones más destacables .................................................. 212 Perspectivas........................................................................... 213
REFERENCIAS ................................................................................. 215 1. 2. 3. 4. 5. 6.
Referencias citadas en este trabajo ........................................... 215 Referencias generales sobre química Cuántica ............................ 223 Artículos de Interés general en el tema de las integrales montecarlo, cuasimontecarlo e integrales moleculares .................................. 225 Referencias sobre el método montecarlo (Aspectos más generales)227 Referencias sobre la Integración Montecarlo ............................... 227 Recursos en la Red.................................................................. 228
TABLAS ........................................................................................... 235 Tabla 1: Armónicos esféricos ........................................................... 236 Tabla 2: Relación entre las constantes empleadas por Stewart et al y por Fernández Rico et al ................................................................ 236 Tabla 3: Armónicos reales-Armónicos sólidos .................................... 237 Tabla 4: Agunas bases estándar del programa GAUSSIAN .................. 239 Tabla 5: Molécula de H2. R(H-H)=0.7414 u.a. (Sobre eje OZ) ............. 240 Tabla 6: Molécula de H2. R(H-H)=1.401 u.a. (Sobre eje OZ) ............... 241 Tabla 7: Molécula de H2. R(H-H)=3.0 u.a. (Sobre eje OZ) .................. 242 Tabla 8: Molécula de LiH. R(Li-H)=2.515 u.a. (Sobre eje OX) .............. 243 Tabla 9: Molécula de LiH. R(Li-H)=3.015 u.a. (Sobre eje OX) .............. 244 Tabla 10: Molécula de LiH. R(Li-H)=3.515 u.a. (Sobre eje OX) ............ 245 Tabla 11: Molécula de H2O R(O-H)=1.63 u.a. ANGULO(HOH)=104.5º . 246 Tabla 12: Molécula de H2O R(O-H)=1.81 u.a. ANGULO(HOH)=104.5º . 247 Tabla 13: Molécula de H2O R(O-H)=2.17 u.a. ANGULO(HOH)=104.5º . 248 Tabla 14: Molécula de CH4. R(C-H)=1.9054 u.a................................. 249 Tabla 15: Molécula de CH4. R(C-H)=2.1089 u.a................................. 250 Tabla 16: Molécula de CH4. R(C-H)=2.5981 u.a................................. 251
íNDICE-v
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
íNDICE-vi
INTRODUCCIÓN
INTRODUCCIÓN
INTRODUCCIÓN Esta tesis continúa el trabajo iniciado en el Laboratorio de Química Cuántica de la Universidad de Cádiz por la Dra. Rosa Rodríguez con su Tesis “Integración Montecarlo en la determinación de Orbitales Atómicos y Moleculares”. Ambas tienen como objetivo principal el desarrollo de métodos que permitan el cálculo de funciones de onda de sistemas polielectrónicos con funciones de base arbitrarias, y el cálculo de funciones de onda de sistemas de partículas que interaccionan con potenciales no electrostáticos. En ambos casos el problema básico se centra en la resolución de las integrales de tipo: (pq | rs) =
∫∫ χ (1) χ (1) V(1,2) χ (2) χ (2) dτ p
q
r
s
1
dτ 2
El enfoque elegido en nuestro laboratorio para este problema ha sido intentar resolverlas por métodos numéricos, lo que a primera vista puede parecer un tanto utópico, por tratarse de integrales de seis dimensiones. Sin embargo se justifica por la sinergia entre dos factores cuyo desarrollo en los últimos años juega claramente a favor de los métodos numéricos. En primer lugar se encuentra el extraordinario perfeccionamiento que están experimentando los métodos de integración Montecarlo. En segundo lugar el aún más extraordinario incremento que viene produciéndose en la capacidad de cálculo de los ordenadores. Explicamos a continuación, un poco más detalladamente, lo que queremos decir.
1. MÉTODOS MONTECARLO Y CUASIMONTECARLO La esencia del método Montecarlo es la experimentación con números aleatorios. El procedimiento usado consiste en diseñar juegos de azar con estos números, esperando obtener de su observación conclusiones útiles para la resolución del problema que se esté estudiando. Aunque se han publicado algunos trabajos relacionados con el método de Montecarlo que no han
INTRODUCCIÓN-vii
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
precisado el uso de ordenadores, lo cierto es que la utilidad del método de Montecarlo se ha visto enormemente incrementada con el uso de las modernas computadoras. En particular, sin el auxilio de éstas, el trabajo presentado en esta memoria hubiera sido absolutamente impensable. Resulta difícil de creer que basándose en el azar puedan obtenerse conclusiones de validez general y, de hecho, quedan investigadores que desconfían todavía de las estimaciones que se consiguen con procedimientos tipo Montecarlo, a pesar de sus múltiples éxitos en el campo de la Investigación Operativa, de la Física y de otras ramas de las Ciencias, como la Biología, la Química, e incluso la Medicina. Los métodos de Montecarlo probabilistas y deterministas.
suelen
clasificarse
en
dos
tipos:
En el Montecarlo probabilista se simulan, mediante los números aleatorios, fenómenos que son aleatorios en la realidad. Los números se eligen de tal forma que reproduzcan la distribución de probabilidad de la población estudiada y, de su observación, se deducen características de ésta. Por ejemplo, la Física Nuclear suministra las funciones que rigen el movimiento de los neutrones. Reproduciendo estas leyes con números aleatorios se puede simular un reactor nuclear y experimentar con él, evitando los problemas de dinero, tiempo y seguridad que implicaría experimentar con un reactor nuclear verdadero. En el Montecarlo determinista, en cambio, se resuelven problemas que no son aleatorios en la realidad, asociándolos con algún experimento aleatorio diseñado expresamente para que se relacionen con el problema estudiado. Un ejemplo de este tipo es una parte de este trabajo: el cálculo numérico de integrales definidas. Los métodos Cuasimontecarlo se usan también para la integración numérica, con un esquema de aplicación totalmente similar al del Montecarlo. Se distinguen solamente en la forma en que se generan los puntos en que hay que evaluar el integrando. En los Cuasimontecarlo la generación de tales puntos no pretende imitar distribuciones aleatorias, sino que se usan números deterministas, diseñados especialmente para que estén “bien distribuidos”. En la actualidad estos métodos se están utilizando en una gran variedad de problemas que no siempre son estrictamente científicos, como por ejemplo en la encriptación de datos, la generación de texturas animadas, etc.
INTRODUCCIÓN-viii
INTRODUCCIÓN
2. PERSPECTIVAS DEL DESARROLLO INFORMÁTICO Cuando aparecieron los primeros ordenadores era difícil imaginar hasta que punto dependeríamos de ellos en la actualidad. Estos primeros ordenadores fueron concebidos para la investigación pero eran inalcanzables económicamente por usuarios con un presupuesto reducido. En muchas ramas de la Química los ordenadores tienen poca utilidad, excepto como herramienta ofimática, pero no así en la Química Cuántica donde se observa un crecimiento vertiginoso de su rama computacional. En los últimos años se ha disparado el uso de ordenadores personales en las universidades, tanto para el cálculo (en las distintas disciplinas) como para otras utilidades afines. Este aumento se ha debido principalmente a tres factores: -
La caída de los precios de los ordenadores
-
El aumento de su potencia
-
La simplificación de los lenguajes y técnicas de programación
La combinación de estos tres factores ha hecho que actualmente se puedan realizar cálculos químico-cuánticos sumamente complejos en ordenadores personales, con un presupuesto realmente reducido. Desde el comienzo de esta tesis, la informática ha progresado radicalmente (ley de Moore) y esperamos que así continúe. En septiembre de 1995 contábamos en el laboratorio con un par de ordenadores, de tipo 486DX, y el programa utilizado para la compilación de los programas era el FORTRAN 77 versión 5 de Microsoft, que actuaba bajo MsDos. La programación de gráficos en este sistema operativo era una ardua tarea. Además los tiempos necesarios para los cálculos mecanocuánticos del tipo que ahora se presentan eran totalmente inaceptables. En aquel momento aparecieron los Pentium y los sistemas operativos gráficos como el Windows 3.XX. Aparecieron nuevos compiladores de FORTRAN como el Powerstation de Microsoft, pero sus subrutinas y funciones graficas no estaban todavía preparadas para el sistema operativo Windows, sino que eran las mismas subrutinas de MsDos. Con la aparición de los nuevos Windows 95 y posteriormente 98, y la aparición de la programación orientada a objeto se necesitaba de una modificación profunda del FORTRAN 77, que vio la luz en el FORTRAN 90 y su actual modificación FORTRAN 95, que se parece al C, e impone una nueva filosofía de programar. Cuando aparece el Visual FORTRAN (entorno visual de Microsoft y compilador FORTRAN de Digital) es cuando
INTRODUCCIÓN-ix
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
nuestro programa toma su forma actual. Se cambian subrutinas y se amplían los métodos de cálculo de integrales y los métodos para lograr la convergencia. En 1999 aprovechamos el Visual Basic para la programación de la interfase programa-usuario. Era fácil de utilizar, y relativamente sencillo de fusionar con el FORTRAN. De ahí que se decidiéramos hacer un generador de datos de fácil uso e intuitivo, como complemento del programa FORTRAN de cálculo de propiedades moleculares. Actualmente el programa tiene un generador de datos hecho en Visual Basic (Gendat) y una parte específica para el cálculo hecho en FORTRAN 95 (OMCM), dando lugar al programa UCA_CMC que se presenta en esta Tesis. En resumen, simplemente en el tiempo que ha durado la elaboración de esta tesis, la potencia de cálculo, incluso de los laboratorios más modestos se ha multiplicado por diez o más, y la capacidad de realizar programas complejos en un tiempo razonable, también. La combinación de este hecho con el desarrollo de los procedimientos para mejorar la eficiencia de las integraciones tipo Montecarlo, es lo que ha hecho posibles los cálculos que presentaremos, y promete la posibilidad de realizar otros mucho más precisos en un próximo futuro.
3. DESCRIPCIÓN DEL CONTENIDO DE ESTE TRABAJO La presente Tesis consta de tres partes claramente diferenciadas. En primer lugar una parte matemática basada en la investigación de métodos de generación de números aleatorios de baja discrepancia (cuasialeatorios) para el cálculo de las integrales que aparecen en los cálculos atómicos y moleculares. En segundo lugar una parte informática donde se aplican los modelos teóricos y las técnicas de programación investigadas, cuyo fruto es la realización de los programas OMCM y UCA_CMC. Por último, y principal, la parte de Química Teórica donde se ha investigado el rendimiento de los métodos desarrollados en la primera parte y del programa desarrollado en la segunda parte. Deseamos resaltar, finalmente, que esta Tesis puede ser el preámbulo de futuras investigaciones dentro del campo de la Química Teórica, ya que los métodos y programas desarrollados constituyen herramientas muy versátiles. Algunas de tales investigaciones ya han comenzado a desarrollarse en nuestro laboratorio. Por ejemplo, aparte de lo presentado al final de esta memoria, se está empleando el programa UCA_CMC para realizar cálculos de orbitales moleculares con bases que cumplen con exactitud la “aproximación” ZDO, con objeto de emplearlos como guía en la parametrización de métodos semiempíricos con base extendida y correlación.
INTRODUCCIÓN-x
Capítulo I: ANTECEDENTES
ANTECEDENTES
Capítulo I: ANTECEDENTES 1. INTRODUCCIÓN 2. LA ECUACIÓN DE SCHRÖDINGER 3. MÉTODO HARTREE-FOCK 4. EL MÉTODO DE CAMPO AUTOCONSISTENTE 5. MÉTODOS POST HARTREE-FOCK 6. PROPIEDADES MOLECULARES 7. INTEGRALES DE LOS CÁLCULOS ATÓMICOS Y MOLECULARES 8. INTEGRACIÓN MEDIANTE MÉTODOS MONTECARLO Y CUASIMONTECARLO
1. INTRODUCCIÓN Nos ha parecido interesante comenzar la presentación de la Tesis que hemos realizado con un resumen, necesariamente breve, pero lo más personal posible, de los temas de carácter general que se encuentran más directamente relacionados con el tema concreto que hemos estudiado. La tesis que presentamos se engloba dentro del área de la QuímicaFísica, y más concretamente en el campo de la Química Teórica. Dentro de esta, la tesis se inscribe tanto en el campo de la Química Cuántica como en el de la Química Computacional. El campo que cubre la Química Cuántica es bien conocido. El concepto de Química Computacional es algo más reciente y se refiere a las aplicaciones del cálculo numérico en cualquier campo de la Química. En relación con nuestro trabajo, la Química Computacional nos permite calcular un conjunto básico de propiedades, cuyos exponentes mas
Capítulo I-1
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
destacados son la energía de un determinado conjunto de átomos (molécula), la geometría óptima de un sistema molecular (es decir, la conformación geométrica de los núcleos de tal forma que tengan la energía mas baja con el método de cálculo empleado) y otras propiedades moleculares como el momento dipolar, las constantes de fuerza, las frecuencias vibracionales, etc. Los métodos de Cálculo Teórico de las Propiedades Moleculares se pueden dividir en tres grupos principales (ab-initio, funcionales de la densidad y semiempíricos) dependiendo de que no usen más que las constantes físicas fundamentales o empleen datos experimentales para parametrizar parte de los cálculos. Esta tesis se refiere a los métodos ab-initio para el cálculo de las propiedades moleculares. Dentro del cálculo ab-initio uno de los problemas más arduos que hay que resolver para llegar a resultados útiles, es el cálculo de las integrales moleculares implicadas. Nuestra tesis ha profundizado en la aplicación de determinados métodos numéricos (Montecarlo y Cuasimontecarlo) a la resolución de tales integrales, y ha desembocado en la elaboración de un programa de cálculo ab-initio de propiedades moleculares que puede operar, en principio, con cualquier tipo de funciones de base.
2.
LA ECUACIÓN DE SCHRÖDINGER
2.1.
La ecuación de Schrödinger
En 1926 [Ref. 1] el físico-matemático austriaco Erwin Schrödinger obtuvo la primera ecuación de las ondas materiales que incluía el potencial de una manera satisfactoria:
r ∂Ψ(ri, t) r ˆΨ(ri, t) = ih H ∂t
(1)
Esta ecuación constituye el fundamento teórico de casi todo lo que sabemos
ˆ es el actualmente sobre el comportamiento químico de los átomos. En ella, H operador Hamiltoniano: ) H=
n
h2
∑ − 2m i =1
r ∇ 2i + V(ri , t)
(2)
i
h2 2 ∇ i es un operador relacionado con la energía cinética de cada 2m r electrón del sistema y V(ri , t) es un operador relacionado con la energía en el que −
Capítulo I-2
ANTECEDENTES
potencial del conjunto de electrones. La ecuación (1), que es la ecuación de Schrödinger original, incluye el tiempo. Cuando el potencial no depende de éste, es posible simplificar la función de onda, escribiéndola como un producto de una función de las coordenadas de las partículas por otra que depende exclusivamente del tiempo: r r Ψ(ri , t) = ψ(ri ) ⋅ ϕ(t)
(3)
Introduciendo la expresión (3) en la (1) se obtiene dos ecuaciones independientes (una para la función espacial y otra para la función temporal). De ellas, la de mayor trascendencia para nuestro trabajo es la ecuación de Schrödinger independiente del tiempo:
r r ˆ ψ(ri ) = E ψ(ri) H
(4)
ˆ es el operador Hamiltoniano que aparece en la ecuación (1). El donde H parámetro E que aparece en esta ecuación representa la energía total del sistema. La ecuación (4) admite múltiples soluciones, que corresponden (salvo en caso de degeneración) a las distintas energías posibles para el sistema. La función de onda que corresponde a la energía mas baja es la que representa el estado fundamental del sistema. En nuestro trabajo no se ha tomado en cuenta ningún efecto relativista, de manera que no se tratará aquí de la ecuación de Dirac [Ref. 2] ni de las generalizaciones de ésta a los sistemas polielectrónicos [Ref. 3].
2.2.
El Hamiltoniano Molecular
Para un sistema molecular la función de onda depende de las coordenadas de todas las partículas en el sistema. Núcleos y electrones intervienen en pie de igualdad en el Hamiltoniano molecular (los núcleos se consideran, al igual que los electrones, partículas puntuales sin estructura r interna). Denominaremos ri al vector de las coordenadas del i-ésimo electrón r y Ri al vector de las coordenadas del i-ésimo núcleo. El Hamiltoniano molecular se expresa como la suma de un término de energía cinética y otro de energía potencial: ˆ=ˆ ˆ H T+V
(5)
en la que:
Capítulo I-3
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
2
h ˆ T =− 2
1
∑m ∇
2 r ˆ(r ) = e V 4πε0
2 i
i
i
(6)
ZiZ j
∑ rr − rr i> j
i
(7)
j
siendo mi las masas de los componentes de la molécula (núcleos y electrones) y Zi sus cargas eléctricas en unidades atómicas [Ref. 4]. Para una molécula, es razonable dividir el operador de energía cinética en dos sumandos, uno sobre los electrones, y el otro sobre los núcleos. Análogamente dividimos el operador de energía potencial en tres términos, uno representando las interacciones entre núcleos, entre núcleos y electrones y entre electrones (usaremos los índices µ,ν,... para denotar los electrones y A, B,... para denotar los núcleos): r r r r r r ˆ=ˆ ˆNN(R) + V ˆeN(r , R) + V ˆee (r ) H TN(R) + ˆ Te (r ) + V
(8)
Utilizando unidades atómicas: ˆ=− H
1
∑ 2M A
A
∇2A −
1
ZAZB
∑2∇ + ∑ R µ
2 µ
A >B
AB
−
ZA
∑r
A,µ Aµ
+
1
∑r
(9)
µ > ν µν
r r r r r r donde rµν = rµ − rν , rAµ = rA − rµ , R AB = rA − rB y MA es la masa del núcleo A. La ecuación (9) es conocida como Hamiltoniano “exacto” no relativista. Sin embargo es importante recordar que este Hamiltoniano no tiene en cuenta efectos que a veces son importantes. En primer lugar, no se toman en cuenta las interacciones espín-órbita debidas al movimiento orbital del electrón alrededor del núcleo. En segundo lugar, aunque la velocidad de un electrón en un átomo es normalmente pequeña, las correcciones relativistas de masa pueden llegar a ser apreciables para electrones internos en átomos pesados.
2.3.
La aproximación de Born y Oppenheimer
Núcleos y electrones participan de forma equivalente en la ecuación de Schrödinger. Sin embargo, la masa de los núcleos es mucho mayor que la masa de los electrones y ello permite utilizar la analogía clásica de un movimiento muy rápido (el de los electrones) frente a otro bastante mas lento (el de los núcleos) para intentar un desacoplamiento de ambos movimientos. Físicamente, lo que esto implica es que los electrones reaccionan rápidamente ante cualquier cambio de la configuración nuclear y que, por lo tanto, la
Capítulo I-4
ANTECEDENTES
distribución electrónica dentro de un sistema molecular esencialmente de la posición de los núcleos y no de su velocidad.
depende
La primera justificación matemática satisfactoria de la separación de los movimientos electrónicos y nucleares fue elaborada por Born y Oppenheimer en 1927. Posteriormente se han elaborado otras, por ejemplo por el mismo Born en 1951 [Refs. 5 y 6], lo que ha dado lugar a cierta confusión en el uso del término “Aproximación de Born y Oppenheimer”. Matemáticamente, lo que implica la posibilidad de desacoplar ambos movimientos es que podemos escribir la función de onda del sistema como: r r r r r ψ(r , R) = ψ(r;R) ⋅ ψ(R)
(10)
r r donde la notación (r;R) significa que la función de onda electrónica depende r paramétricamente de la posición de los núcleos, es decir, se fija el valor de R r r r para un valor determinado de “ R A ”, y se resuelve la ecuación para ψ(r , R A ) . r Haciendo lo mismo para un rango de R obtenemos una función de energía r r potencial para el movimiento de los núcleos. La función ψ(r , R ) es función propia del hamiltoniano electrónico: r r r r ˆel = ˆ ˆeN(r , R) + V ˆee (r ) H Te (r ) + V
(11)
r y ψ(R) es función propia del hamiltoniano nuclear: r r ˆN = ˆ ˆ(R) H TN(R) + V
(12)
r ˆ(R) es una serie de Taylor en torno a una posición de equilibrio del donde V r ˆNef (R) potencial V (potencial intramolecular). Definimos el potencial r ˆNef (R) como el potencial que gobierna el movimiento de los intramolecular V núcleos de una molécula en la medida en que éste pueda considerarse independiente del movimiento de los electrones, y viene dado [Ref. 7] por la suma de la energía electrónica y la repulsión internuclear para una molécula aislada: r r r ˆNef (R) = Enel(R) + V ˆNN(R) V
(13)
Las ecuaciones que se obtienen no son triviales, debido a que se debe realizar un cambio de variables que permita separar el centro de masas del
Capítulo I-5
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
sistema y realizar una serie de manipulaciones. El Hamiltoniano del sistema polielectrónico podemos escribirlo resumidamente como: r r r r r r ˆ=ˆ ˆe (r ) + V ˆNN(R) + V ˆeN(r , R) + V ˆee (r ) H TN(R) + T
(14)
r El término ˆ TN(R) se desprecia ya que su efecto sobre la función de ondas es r mucho menor que el de ˆ Te (r ) . Lo que se acostumbra hacer, y que muchas veces se confunde con la propia aproximación de Born y Oppenheimer, es la llamada aproximación de los núcleos fijos donde lo que efectivamente se realiza es identificar como nula la componente de energía cinética de los núcleos. r r r r r ˆel = T ˆe (r ) + V ˆNN(R) + V ˆeN(r , R) + V ˆee (r ) H
(15)
Hecho esto, la ecuación que entonces queda para resolver tiene la forma: r r r r r ˆ el ψ el (r ; R 0 ) = E el (R 0 )ψ el (r ; R 0 ) H
(16)
r donde R 0 se refiere ahora a una posición fija determinada de los núcleos. Nótese que la solución de esta ecuación electrónica nos da como resultado r una energía electrónica Eel(R 0 ) que es una función de la posición de los r ˆNN(R ) ya que núcleos. Frecuentemente se prescinde también del término V r para un R fijo es un valor constante, cuyo efecto sobre la ecuación (4) es un simple cambio de origen de las energías. El hamiltoniano “útil” queda entonces reducido a: r r r r ˆel = ˆ ˆeN(r , R) + V ˆee (r ) H Te(r ) + V
(17)
Más allá de la aproximación de núcleos fijos se encuentran otras, como la de Born y Oppenheimer [Ref. 5] o la aproximación adiabática [Ref. 8]. Ambas desembocan en el tratamiento del problema “núcleos + electrones” mediante dos ecuaciones separadas, una para los electrones que es distinta para cada posición de los núcleos: r r ˆelψ el(r;R) = − 1 H 2
ZA
∑∇ − ∑R 2 µ
µ
A,µ
Aµ
+
r r r r r 1 ψ el(r;R) = Eel(R)ψ el(r;R) µν
∑r µ>ν
(18)
y otra para los núcleos, que es distinta para cada estado electrónico obtenido en (18):
Capítulo I-6
ANTECEDENTES
r ˆNψN(R) = − H
1
∑ 2M A
A
r ∇2A + Enel(R) +
r r ZAZB ψN(R) = Etot ψN(R) AB A >B
∑R
(19)
en la que el “potencial efectivo” que actúa sobre los núcleos es: r r r ˆNef (R) = Enel(R) + V ˆNN(R) V
(20)
Obsérvese que el potencial efectivo es la suma de una energía electrónica r r ˆNN(R ) , totalmente Enel(R) (valor propio de (18)) y una repulsión internuclear V clásica. Para algunos sistemas este balance es una función con al menos un mínimo, que si resulta suficientemente profundo justifica la aparición de una molécula estable. No obstante la aparición de tal mínimo no es un hecho necesario. En la mayoría de los sistemas sólo se produce para alguno de sus estados electrónicos. A pesar del enorme número de moléculas que conocemos, estas no son mas que una pequeña fracción de un número mucho mayor de estados electrónicos, de los que sólo una pequeña parte poseen el mínimo que, si es suficientemente profundo, caracterizará a una molécula. Una discusión detallada de este tema puede encontrarse en [Ref. 7]. r r Notemos que la función de onda electrónica ψ el(r;R) tiene que cumplir dos condiciones. Por una parte, dado que el cuadrado de su módulo representa la densidad de probabilidad de encontrar las partículas en puntos determinados del espacio, entonces, la integral sobre todo el espacio de esa densidad debe valer 1. Para ello, usualmente se multiplica la función de onda electrónica por una constante de normalización N, tal que se cumpla que:
∫ψ
r r
el(r; R )
2
r dr = 1
(21)
Por otra parte, dado que los electrones son fermiones, es necesario que la función de onda electrónica sea antisimétrica al intercambio interelectrónico, lo que matemáticamente se expresa:
r r r r r r r r ψ el(r1,..., rµ , rν ,..., rn ) = −ψ el(r1,..., rν , rµ ,..., rn )
(22)
3. MÉTODO HARTREE-FOCK Las funciones estacionarias de los átomos hidrogenoides suelen denominarse orbitales en recuerdo de las órbitas del primitivo modelo atómico de Bohr. Sin embargo la palabra orbital tiene una acepción mucho mas amplia
Capítulo I-7
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
que la de “solución de la ecuación de ondas de un átomo hidrogenoide”: También es una función monoelectrónica empleada para describir, en combinación con otras, la función de ondas de un sistema polielectrónico. Si se desprecia la repulsión entre los electrones de un átomo, su hamiltoniano queda reducido a una suma de operadores monoelectrónicos análogos a los de un átomo hidrogenoide. En una aproximación de electrones completamente independientes, la función de onda electrónica puede suponerse igual a un producto de orbitales hidrogenoides (diferentes, según el Principio de Exclusión). Si se quieren tener en cuenta los efectos de la degeneración de canje de manera mas cuidadosa, puede emplearse un producto antisimetrizado, es decir un determinante de Slater.
3.1.
Orbitales de Hartree. Métodos autoconsistentes
El modelo de electrones totalmente independientes constituye una aproximación muy grosera de la realidad, aunque permite entender el origen de algunas de las propiedades de los átomos (por ejemplo, sus propiedades periódicas). El primer método que tuvo éxito en la obtención de unos orbitales atómicos netamente mejores que los hidrogenoides fue elaborado por D. R. Hartree en 1928 [Ref. 9], a partir de unas ideas intuitivas publicadas por Bohr en 1922 [Ref. 10]. Hartree observó que si la función de onda electrónica de r un átomo es un producto de orbitales φ µ (rµ ) : n
ψ(1,2,K ,N) = las
funciones
r φ µ (rµ )
pueden
∏ φ (µ)
(23)
µ
µ =1
determinarse
resolviendo
la
ecuación
de
Schrödinger “efectiva”: ∇2µ Z ˆ ˆ hµ φµ = − − + Jv (µ) φµ = rµ v ≠ µ 2
∑
εµ φµ
(24)
en la que ˆ Jν (µ) =
∫
φ*ν (ν)φν (ν) r d rν rµν
(25)
r es el operador de Coulomb, que consiste en multiplicar por la función de rµ obtenida al realizar la integración respecto a ν. El operador
∑ ˆJ
ν
depende
entonces de cual sea el conjunto de orbitales utilizados por el átomo
Capítulo I-8
ANTECEDENTES
considerado y se refiere a las coordenadas de un solo electrón, pues las de los demás desaparecen al integrar. Existe una ecuación del tipo (25) para cada orbital de un átomo y el conjunto se denomina ecuaciones de Hartree. Hartree mismo ideó un procedimiento para resolverlas, en el que se considera que los orbitales φµ son funciones propias de un potencial central: φµ (µ) ≈ R(rµ) ⋅ Ylm(θµ , ϕµ )
(26)
y se introduce el método del campo autoconsistente (SCF=Self Consistent Field). Este método consiste en utilizar un conjunto de orbitales iniciales φ(µ0) para construir los operadores ˆ Jv de (24), y mediante un proceso iterativo en la resolución de esta ecuación llegar a los llamados orbitales autoconsistentes de Hartree: φ(µn) ≈ φ(µn+1) ≈ φSCF µ
(27)
Aunque hay casos en los que el procedimiento de Hartree no converge, en la mayoría de ellos basta recomenzar el cálculo partiendo de una aproximación de orden cero diferente, para lograr la convergencia. Aún así, algunas veces es imposible resolver las ecuaciones de Hartree, en unos casos por problemas numéricos y en otros casos porque la configuración supuesta no es válida (por ejemplo por ser de energía demasiado alta). En la actualidad es posible resolver las ecuaciones de Hartree sin necesidad de considerar que los orbitales son funciones propias de un potencial central. Se puede aplicar, por ejemplo, el método Roothaan que veremos más adelante, en el que la resolución de las ecuaciones sigue siendo iterativa. En la práctica se emplean otras ecuaciones publicadas posteriormente, las de Hartree-Fock, que conducen a funciones de onda atómicas y moleculares algo mejores que las de Hartree, con un esfuerzo de cálculo similar. En principio, las ecuaciones de Hartree sirven exclusivamente para generar los orbitales autoconsistentes. La energía del sistema debe calcularse a posteriori, empleando la función total obtenida y el hamiltoniano verdadero, y
nunca
sumando
los
valores
propios
εSCF µ
pues
tal
procedimiento
sobrestimaría fuertemente la repulsión interelectrónica. Cuando se estudian los átomos polielectrónicos es frecuente usar, en vez de los orbitales autoconsistentes de Hartree, orbitales hidrogenoides con
Capítulo I-9
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
carga nuclear efectiva, en los cuales la carga Z* depende del orbital y átomo considerado. Cada orbital hidrogenoide con carga efectiva es función propia de un hamiltoniano hidrogenoide en el que el potencial es (en unidades atómicas): Vi(r) = −
Zi* r
(28)
con un valor de Z i* distinto para cada orbital del átomo considerado. Esta variación respecto a la carga nuclear verdadera se puede interpretar como el resultado de un “apantallamiento” del núcleo, producido por el resto de los electrones, y se puede determinar aplicando el método variacional o aplicando las reglas de Slater, el cual propuso no solo unas reglas sencillas para obtener la Z* de (28) sino también una simplificación de los orbitales hidrogenoides que convierte, cuando se usan en átomos polielectrónicos, en los orbitales de Slater (STO) [Ref. 4].
3.2.
Orbitales de Hartree-Fock
Las ecuaciones de Hartree-Fock son las que determinan los orbitales más adecuados, desde el punto de vista variacional, para construir una función de onda polielectrónica monodeterminantal de electrones apareados: ψ = φ1φ1φ2 φ2...φnφn
(29)
El principio variacional permite deducir las ecuaciones monoelectrónicas que deben cumplir las funciones ϕi del determinante de Slater (29), para que este sea la mejor solución de una ecuación de Schrödinger electrónica dada, teniendo en cuenta la estructura admitida para la función de onda ~ ψ. Dado que los electrones son fermiones, tienen espín semientero que usualmente se representa como alfa (espín +1/2) o beta (espín -1/2). Los espín-orbitales moleculares pueden escribirse en general (dentro de las aproximaciones que estamos empleando) como un producto de una función espacial y una función de espín: r r φi(ri ) = ϕi(ri )si(σ)
(30)
Las dos funciones de espín, α y β tienen la propiedad de ser ortonormales
∫ α dσ = ∫ β dσ = 1 2
2
Capítulo I-10
∫ αβdσ = 0 .
ANTECEDENTES
Considerando entonces las funciones de espín, vemos que dos espínorbitales moleculares pueden ser distintos, ya sea porque su parte espacial (el orbital molecular) es diferente o porque su función de espín sea diferente (o, por supuesto, ambas). Ello lleva entonces a que cada orbital molecular puede usarse dos veces (una con espín alfa y otra con espín beta). Esto, que se conoce con el nombre de principio de exclusión de Pauli (dos fermiones no pueden ocupar un estado de exactamente la misma energía, por lo cual los espínorbitales deben diferir en alguno de sus componentes) conduce al principio de aufbau o de construcción. Puesto en forma sencilla, lo que dice es que los orbitales moleculares se van ocupando en forma progresiva desde aquél con energía orbital mas baja hacia arriba. El caso en que tenemos n/2 orbitales moleculares ocupados por n electrones, la mitad de ellos con espín alfa y la otra mitad con espín beta, es lo que se conoce con el nombre de configuración electrónica de capas cerradas o de electrones apareados, y es muy importante en química, ya que generalmente constituye el estado fundamental de la mayoría de las moléculas, especialmente en el caso de las orgánicas (si bien moléculas muy comunes, como el O2 por ejemplo, no tienen una configuración de capa cerrada como estado fundamental del sistema). Todas las configuraciones electrónicas que no son de capas cerradas se identifican como de “capas abiertas” o “electrones desapareados”. La multiplicidad del sistema (que se relaciona con el valor propio del operador ˆ 2 ) coincide con el número de electrones desapareados más una de spín S unidad, por lo cual, la multiplicidad de un sistema de capa cerrada será 1 y se dice que estamos en presencia de un singlete. Todos los sistemas de capa cerrada son singletes, pero no todos los singletes son de capa cerrada ni todos los estados fundamentales son singlete. Al haber dos electrones desapareados con espín alfa S=1 la multiplicidad es 3, por lo cual nos encontramos en presencia de un triplete. Es posible tener sistemas con multiplicidades mas altas (cuadrupletes, quintupletes, etc). También es posible que aún cuando dos electrones estén descritos mediante dos orbitales espacialmente diferentes ( px y py por ejemplo), no tengan el mismo espín sino espín opuesto, con lo que también nos encontraremos en presencia de un singlete, si bien este no es de capa cerrada. Finalmente, por ejemplo, si ionizamos el átomo de C para obtener el catión C+ nos encontramos con que uno de los dos electrones desapareados ahora ha desaparecido (por convención, el que tiene espín beta). Por lo tanto, ahora S=+1/2 y consecuentemente 2S+1 = 2, con lo cual estamos en presencia de un doblete. La mayoría de los radicales libres de la química orgánica son dobletes.
Capítulo I-11
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
Finalmente, mencionaremos un detalle referido a los cálculos propiamente: Si bien el hecho de que cada electrón ocupe un espín-orbital diferente nos permite utilizar dos veces cada orbital molecular para construir el determinante de Slater (una vez combinado con la función de espín adecuada), podemos levantar esa restricción y establecer que independientemente de las funciones de espín, las partes espaciales serán todas ellas diferentes. Los métodos en que se mantiene la restricción de igual parte espacial para las dos funciones de espín se conocen con el nombre de restringidos. Por el contrario los métodos que no aplican tal restricción se llaman no restringidos. Nótese que los cálculos restringidos y no restringidos pueden realizarse en sistemas de capa cerrada o capa abierta, dado que son conceptos diferentes. r Cuando los espinorbitales ϕi(r )si(σ) con que se construye (29) son ortonormales, y se encuentran apareados en la forma ...ϕiϕi... , la energía ~ aproximada E puede escribirse en función de integrales sobre coordenadas únicamente espaciales. Resulta, para sistemas que tengan sus electrones apareados: m ~ E = 2 Hii +
∑ i=1
∑∑ (2J
ij
i
− Kij )
(31)
j
en donde : Hii =
∫
∇2µ ϕi∗(µ) − − 2
∫ ∫ ϕ (µ)ϕ (µ) r
Kij =
∫ ∫ ϕ (µ)ϕ (µ) r
∗ i
∗ i
(32)
r r ϕ∗j (ν)ϕ j(ν)d rµd rν
(33)
r r ϕ∗j (ν)ϕi(ν)d rµd rν
(34)
1
Jij =
i
µν
1
j
r ZA ϕi(µ)d rµ Aµ
∑r
µν
A
En adelante se llamará orbitales a las funciones ϕi, y espinorbitales φi a las que incluyen la coordenada de espín. El valor de la energía (31) depende, para cada molécula, de cuales sean los orbitales ϕi con los que se construya su función polielectrónica (29). La energía es, por consiguiente, un funcional de los orbitales. Las condiciones ~ que deben cumplir las funciones ϕi para que el funcional E [ϕi] sea un máximo o mínimo se estudia mediante el cálculo variacional. Sin embargo, para
Capítulo I-12
ANTECEDENTES
encontrar las ecuaciones que deben cumplir las funciones ϕi no basta con aplicar directamente las reglas del cálculo variacional a la energía (31). Es necesario tener en cuenta que, para que la expresión empleada para la energía sea válida, los orbitales ϕi deben ser ortonormales: Sij =
r
∫ ϕ ϕ dr = δ ∗ i j
(35)
ij
~ Esto obliga, para obtener los extremos del funcional E [ϕi] restringidos por las condiciones (35), a aplicar un método análogo al de los multiplicadores de Lagrange, que consiste en hallar los extremos libres del funcional: ~ G[ϕi ] = E [ϕi ] +
∑ ∑ λ (S ij
i
ij
− δij )
(36)
j
Las reglas del cálculo variacional permiten obtener directamente las ecuaciones que deben cumplir las funciones ϕi para que (36) sea un extremo. Estas, que son tantas como pares de electrones, quedan en función de los parámetros desconocidos λij: hµ + ˆ
∑ [2ˆJ (µ) − Kˆ (µ)]ϕ (µ) = ∑ (−λ j
j
i
j
ij
/ 2) ⋅ ϕ j (µ)
(37)
j
siendo: 2
∇µ ˆ hµ = − − 2 ˆ Jj(µ) =
ZA
∑r A
1
∫ ϕ (ν) r ∗ j
µν
(38)
Aµ
r ϕ j (ν)d rν
(39)
ˆ j por la propiedad: y definiéndose el operador de canje K r ˆ j(µ)ϕ (µ) = ϕ (µ) ϕi∗(ν) 1 ϕ (ν)d rν K i j j rµν
∫
(40)
Las ecuaciones integrodiferenciales acopladas (37) son de tal complejidad que difícilmente resultarían resolubles, si no fuera por una afortunada coincidencia: Debido a que los determinantes son invariantes frente a las transformaciones unitarias, el conjunto de orbitales ϕi que suministran una función determinantal ~ ψ dada no es único. Por ejemplo, si se sustituyen los
Capítulo I-13
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
orbitales ϕi que determinan las columnas del determinante (29) por combinaciones lineales: ϕ'i =
∑T ϕ
(41)
ij j
j
en las que los coeficientes Tij formen una matriz unitaria (T-1=T+ ya que son reales), la función polielectrónica ~ ψ no se altera. Esta propiedad permite seleccionar entre todos los posibles conjuntos de funciones ϕi que hagan mínimo G[ϕi], uno para el que los multiplicadores λij valgan cero si i ≠ j:
λ ij = −2ε i δ ij
(42)
con lo que las m ecuaciones (37) se reducen a una sola, igual para todos los electrones: hµ + ˆ
ˆ ˆ (µ)]ϕ (µ) = ε ϕ (µ) µ − [ 2 J ( ) K ∑ j
j
j
i
i i
(43)
o abreviadamente: ˆϕi = ε i ϕi F
(44)
ˆ es el operador de Fock, común a todos los electrones. Como los donde F orbitales ϕi de una función determinantal han de ser diferentes, se deben hallar tantas soluciones de (44) como orbitales se necesiten para construir ~ ψ. Cuando lo que se está determinando es la función de ondas electrónica del estado fundamental, normalmente se eligen las m soluciones que tengan ~ menor εi, ya que aunque la energía aproximada E no sea la suma de los valores propios: ~ E ≠
∑ε
(45)
i
i
habitualmente,
~ E será mínima cuando lo sea
∑ε
i
.
i
La búsqueda de soluciones de las ecuaciones de Hartree-Fock (44) es necesariamente iterativa, pues las propias soluciones ϕi se encuentran ˆ . Su resolución es análoga al incluidas en la definición de los operadores ˆJ y K método del campo autoconsistente de Hartree y da lugar a los orbitales autoconsistentes de Hartree-Fock. Conviene señalar finalmente, que el
Capítulo I-14
ANTECEDENTES
procedimiento descrito no siempre es convergente, aunque lo será siempre que se use una aproximación inicial conveniente. Por este motivo, para resolver los problemas de convergencia suele bastar con reiniciar el cálculo fallido con unos orbitales ϕ(i0) más adecuados. Otro problema que se presenta con frecuencia, que también puede resolverse cambiando de orbitales de partida, es una convergencia hacia soluciones no deseadas. Por ejemplo, a veces se obtienen orbitales que serían adecuados para describir alguno de los estados electrónicos excitados del sistema que se esté estudiando, en vez del estado fundamental.
4. EL MÉTODO DE CAMPO AUTOCONSISTENTE A pesar de resultar un tanto elemental, hemos decidido detallar aquí la obtención y proceso de resolución de las ecuaciones de Roothaan-Hall en la forma en que luego las hemos programado. Ello facilitará, probablemente, la labor de perfeccionamiento de nuestro programa que se irá encomendando a futuros colaboradores del laboratorio. Los programas que resuelven las ecuaciones de Roothaan y sus generalizaciones a electrones desapareados son hoy día extraordinariamente numerosos. Nosotros hemos preferido programarlas una vez más, a utilizar uno de los programas existentes, con objeto de controlar al máximo el proceso en el que se introducirán las modificaciones que constituyen el objeto de esta tesis.
4.1.
Capas cerradas (RHF). Ecuaciones de Roothaan-Hall
Aunque las ecuaciones de Hartree-Fock son en principio resolubles, en la práctica sólo pueden tratarse con ellas problemas atómicos, debido a la
ˆ ⋅ ϕi = εi ⋅ ϕi que hay que complicación de la ecuación de valores propios F resolver en cada iteración. Los problemas moleculares no pudieron tratarse hasta que, en 1951, C. J. Roothaan en EE.UU., e independientemente G. G. Hall en Gran Bretaña, combinaron las ideas de Hartree y Fock con la hipótesis de que los orbitales óptimos ϕi se pudiesen expresar como combinaciones lineales de las funciones de base adecuadas. Si cada orbital ϕi del determinante de Slater (29) se expresa mediante una combinación lineal de las funciones de base χp:
Capítulo I-15
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
k
ϕ i = ∑ C pi χ p
(46)
p =1
la energía (31) puede escribirse en función de los coeficientes Cµi: ~ E =2
∑∑∑ C
c ∗ piCqiHpq
i
p
+
q
∑∑∑∑∑∑ C
∗ ∗ piCqiCrjCsj
i
en donde, para abreviar monoelectrónicas:
j
p
la
q
r
notación
∫
c Hpq =
[2(pq | rs) − (ps | rq)] (47)
s
se
han
definido
las
ˆ χ (µ)d rr χp∗ (µ) ⋅ h i q µ
integrales
(48)
y bielectrónicas: 1
∫∫ χ (µ)χ (µ) r ∗ p
(pq | rs) =
q
µν
r r χr∗(ν)χs (ν)d rµd rν
(49)
Las condiciones de ortonormalidad (35) equivalen a las condiciones:
∑ ∑ (C
∗ piCqjSpq
p
)
− δij = 0
q
∀i, j
(50)
las que: Spq =
r
∫ χ χ dr ∗ p q
(51)
son las integrales de solapamiento. Como la condición (50) es de “obligado cumplimiento” la condición de energía mínima se debe obtener por el método de los multiplicadores de Lagrange, igualando a cero las derivadas respecto a los coeficientes Cpi de la función auxiliar: ~ G(Cpi ) = E(Cpi ) +
λij
∑∑ ∑∑ i
j
p
q
CpiCqjSpq − δij
(52)
en la que λij son los multiplicadores de Lagrange. Nótese que ahora se manejan funciones, en vez de funcionales. Sustituyendo E(Cpi) por su valor (47), y eligiendo λij=-2εiδij, es decir eligiendo las soluciones que simplifican el problema, exactamente igual a como se hace al deducir las ecuaciones de Hartree-Fock:
Capítulo I-16
ANTECEDENTES
∑ ∑∑
G(Cpi ) = 2
i
p
q
∑∑ ∑∑∑∑ i
j
p
q
r
∗ c Cpi Cqi Hpq − εiSpq − 1 +
(
)
[2(pq | rs) − (ps | rq)]
∗ ∗ Cpi CqiCrj Csj
(53)
s
Derivando (53) respecto a Cqi e igualando a cero se obtiene una ecuación homogénea. Repitiendo para todos los valores posibles del índice q se obtiene un sistema de ecuaciones homogéneas:
∑ (F k
pq
p =1
− εiSpq )Cpi = 0
(q = 1,2,..., k)
(54)
en el que: c Fpq = Hpq +
∗ rj sj (pq | rs)
∑ ∑ ∑ 2C C j
r
s
−
1 (ps | rq) 2
(55)
El sistema (54) es lo que se conoce como ecuaciones de Roothaan, y tiene tantas soluciones como funciones de base hayamos empleado. Seleccionando tantas como pares de electrones tenga el sistema podemos construir con ellas el determinante de Slater que represente un estado electrónico de la molécula tratada. Normalmente se escogen las soluciones que conduzcan a la menor energía. Los orbitales obtenidos al resolver las ecuaciones y no usados son los “orbitales virtuales” que sirven para obtener aproximaciones a estados electrónicos excitados por sustitución de algún orbital “ocupado” en el estado fundamental por alguno de los virtuales. También tienen utilidad, como se verá más adelante, para construir funciones de onda menores que la HartreeFock. Aunque el sistema (54) no es lineal - pues los propios elementos Fpq contienen las incógnitas Crj - la propia forma de las ecuaciones ya sugiere una manera de resolverlas, que es la siguiente: a) Se elige un conjunto de orbitales moleculares aproximados ϕ(i0) , definidos por unos coeficientes de orden cero C(pi0) . (0) , con los coeficientes b) Se construye una matriz de Fock de orden cero Fpq
C(pi0) .
Capítulo I-17
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
(0) c) Se resuelve el sistema de ecuaciones (54) con la matriz Fpq obtenida con
los coeficientes C(pi0) , lo que conduce a unos nuevos valores C(pi1) :
∑ (F
(0) pq
− εiSpq )C(pi1) = 0
(56)
d) Con los coeficientes obtenidos C(pi1) se construye una matriz de Fock (1) mejorada Fpq , con la que se plantea un sistema de ecuaciones (54) más
próximo al exacto que el (56):
∑ (F
(1) pq
− εiSpq )C(pi2) = 0
(57)
e) Se repite el ciclo (c!d!c) hasta que los coeficientes obtenidos no se diferencien significativamente de los del ciclo precedente,
C(pin) ≈ C(pin−1)
(58)
Se dice, entonces, que "se ha alcanzado la autoconsistencia en coeficientes". A menudo, en lugar de (58) se emplea como criterio de autoconsistencia la repetición de los valores propios ε(in) ≈ ε(in+1) , que es una condición "menos exigente" (se cumple con mas facilidad), que la de autoconsistencia en coeficientes. Aunque en cada iteración las ecuaciones de Roothaan pueden resolverse como se quiera, por lo común se hace matricialmente (pues es un procedimiento muy adecuado para la programación). Las ecuaciones del tipo (57) son equivalentes a la ecuación entre matrices:
FC=SCΛ en la que (51),
C
(59)
F tiene por elementos los Fpq dados en (55), S tiene los Spq dados en
es la matriz formada por los coeficientes que se buscan y Λ es una
matriz diagonal formada por los valores propios εi. Calculando la matriz
S-1/2 y
definiendo:
D = S1/2 C
(60)
se ve que (59) es equivalente a la ecuación: (S-1/2
Capítulo I-18
F S-1/2) D = D Λ
(61)
ANTECEDENTES
que puede resolverse, como tantos otros problemas, por diagonalización. En efecto, el producto (S-1/2F S-1/2) resulta ser siempre una matriz simétrica, de manera que
D
y Λ pueden obtenerse por el método de Jacobi. Por último, los
coeficientes Cpi se obtienen mediante la relación inversa de (60):
C = S1/2 D
(62)
Es importante subrayar que, si bien hay un sistema de ecuaciones para cada orbital ϕi, resulta que todos los sistemas tienen una forma idéntica, por lo cual basta obtener las soluciones de un solo sistema. Con las m primeras, cambiando simplemente el índice de electrón, se obtienen todos los orbitales que se necesitan para construir el determinante de Slater de la molécula. Evidentemente, para poder plantear (54) es preciso usar un número de funciones de base mayor que el de orbitales a determinar, k>m, para que el sistema de ecuaciones suministre un número suficiente de soluciones. Asimismo, conviene no olvidar que la energía total del sistema polielectrónico no va a ser la suma de los valores propios εi, sino que debe calcularse explícitamente, sustituyendo los valores finales autoconsistentes de los coeficientes Cpi en la expresión (47). Esta se reescribe a menudo como: ~ E =
∑ ∑ ∑C C ∗ pi
i
p
qi
c ( Fpq + Hpq )
(63)
q
debido a que es más fácil de evaluar, y también: 1 ~ 1 c E = ∑ ∑ Ppq( Fpq + Hpq ) = Tr P(F + Hc ) 2 p q 2
[
]
(64)
en donde Tr significa "traza" (=suma de elementos de la diagonal de una matriz) y P representa a la matriz densidad correspondiente al estado polielectrónico considerado. Esta matriz se construye a partir de los coeficientes de los orbitales moleculares ocupados (es decir, los que realmente participan en la función de onda total, no los virtuales), según la fórmula: n
Ppq =
∑nC
∗ i piCqi
(65)
i=1
que, para el estado fundamental de un sistema de electrones apareados equivale a: m
Ppq =
∑ 2C
∗ piCqi
(66)
i=1
Capítulo I-19
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
Si las funciones χp empleadas para desarrollar φ i =
∑C
ip χ p
formasen
i
una base completa, es decir si el desarrollo pudiera considerarse exacto, las ecuaciones de Roothaan (55) y las de Hartree-Fock (44) serían totalmente equivalentes. En la práctica esto no ocurre, de manera que los resultados de Roothaan solo se acercan a los de Hartree-Fock con ciertas bases, que dependen de cual sea la molécula o el átomo tratado. La búsqueda de bases adecuadas para aplicar el método de Roothaan ha sido uno de los temas de investigación de la Química Cuántica durante bastante tiempo, pues tropezaba con la dificultad de que las integrales Hpq, y sobre todo las (pq|rs), que aparecen en las ecuaciones de Roothaan (55), no pueden calcularse con cualquier clase de funciones χp.
4.2.
Tipos de funciones de base
Uno de los objetivos de esta tesis es ampliar el tipo de funciones de base que puedan usarse en los cálculos con el método de Roothaan-Hall, para los que, en la actualidad sólo pueden emplearse funciones de muy pocas clases, debido a los problemas de integración. Entre las bases que se usan hoy día cabe hacer la siguiente clasificación: a) Atendiendo a la clase de funciones utilizadas, la base puede estar formada por funciones exponenciales (ETO) o puede estar formada por Orbitales Gaussianos (GTO). b) Atendiendo al número de funciones utilizadas, la base puede ser mínima, múltiple-zeta, extendida, polarizada, o flotante. Todas las funciones ETO incluyen un factor exp(-αr) ydifieren en que multiplican este factor por distintos polinomios de r. Entre las funciones exponenciales podemos distinguir los siguientes tipos: -
Funciones de Slater (STO): Las funciones STO son de la forma: χ = rne− αr Ylm(θ, ϕ)
(67)
El valor de α puede ser directamente el que corresponda al átomo aislado, ser un valor adaptado al hecho de formar parte de una molécula (pero fijo para cada clase de átomo), o ser un valor obtenido mediante una minimización de la energía independiente para cada función de la base empleada.
Capítulo I-20
ANTECEDENTES
En la expresión (67) de los STO se pueden reemplazar los armónicos esféricos por armónicos cúbicos [Ref. 11] de manera que otra forma de expresarlos es: n
χ = xnx y y znz e− αr
(68)
Las funciones de Slater son las ETO más empleadas debido a la sencillez de su expresión. -
Funciones B: χ=
π − αr e 2
n−1
(n − 1 + i)!
∑ (n − 1 − i)!i!2 (αr) i
i =0
n −i −1
Ylm(θ, ϕ)
(69)
donde n = N - L. Las funciones B fueron propuestas inicialmente por Shavitt y Karplus [Refs. 12, 13 y 14], y posteriormente por Steinborn y otros [Refs. 15 y 16]. Su principal ventaja parece ser que poseen transformadas de Fourier y Laplace especialmente simples. -
Funciones Λ (Slater-Laguerre): χ = e −αr L2nL−+12(2αr) Ylm(θ, ϕ)
(70)
donde n = N – L, y Lqp (x) son los polinomios generalizados de Laguerre definidos por [Ref. 17]: Lqp (x)
p
= (p + q)!
(−1)i xi
∑ (p − i)!(q + i)!i!
(71)
i =0
-
Funciones ψ (Sturmianas): χ = e −αr L2nL−+11(2αr) Ylm(θ, ϕ)
(72)
que tienen la particularidad de ser funciones propias del operador de energía cinética [Refs. 18 y 19]. -
Funciones Ψ: χ = e −αr L2nL−+11(2αr) r −1 Ylm(θ, ϕ)
(73)
que forman secuencia biortogonal con las Sturmianas [Ref. 20]. -
Funciones Hidrogenóides:
Capítulo I-21
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
χ = e −αr / N L2nL−+11(2αr / N) Ylm(θ, ϕ)
(74)
Las funciones GTO poseen, en lugar del factor exp(-αr) que caracterizan a los ETO, un factor exp(-αr2) que facilita decisivamente el cálculo de las integrales (pq|rs): 2
χ = rne −αr Ylm(θ, ϕ)
(75)
Las funciones GTO no suelen emplearse de esta forma, sino en forma cartesiana más o menos equivalente, en la que se utilizan los armónicos cúbicos: n
χ = xnx y y znz e − αr
2
(76)
Las funciones gaussianas que se emplean como funciones de base reciben usualmente el nombre de primitivas y a partir de ellas pueden formarse combinaciones lineales fijas que se llaman funciones contraídas (CGTO). Las funciones contraídas son las que se emplean habitualmente como conjuntos de funciones de base para los cálculos moleculares: χ=
∑C G µ
µ
(77)
µ
donde Gµ representa a funciones del tipo (76). Las bases mínimas contienen únicamente a los orbitales atómicos de las subcapas ocupadas (total o parcialmente) en el estado fundamental del átomo considerado. Por ejemplo, en un átomo de carbono - o uno de oxígeno, nitrógeno, etc. - la base mínima está formada por cinco funciones (de tipo 1s, 2s, 2px, 2py, 2pz). Las bases doble-zeta contienen doble número de funciones que las mínimas, pues emplean dos orbitales, con distinto valor del exponente α, por cada orbital de una base mínima. Las bases extendidas son las que incluyen orbitales de capas mas "altas" que las ocupadas en el estado fundamental. Por ejemplo, una base extendida en el átomo de carbono podría incluir orbitales 3s y 3p. Las bases con polarización incluyen funciones con armónicos esféricos más "altos" que los correspondientes a las subcapas ocupadas. Por ejemplo, una base que incluya orbitales de tipo p sobre los átomos de hidrógeno, o de tipo d sobre los átomos de carbono, constituye un ejemplo de base polarizada. Por último, base flotante es aquella en que las funciones no están (necesariamente) centradas sobre los núcleos de los átomos.
Capítulo I-22
ANTECEDENTES
Nombre
Explicación
Bases mínimas (BM)
Contiene solamente el número necesario de funciones para representar los orbitales atómicos de las subcapas ocupadas (total o parcialmente) de los átomos en las moléculas (1 para el H, 5 para C, O, etc)
Contiene dos (DZ) o tres (TZ) bases mínimas que se Bases Doble diferencian en sus exponentes (por cada átomo de H, por Zeta (DZ) y ejemplo, una base DZ contiene 2 funciones 1s con distintos Triple Zeta (TZ) exponentes) Es una BM para los electrones de las capas interiores y DZ (DZV) ó TZ (TZV) para los electrones de la capa de valencia Base "SplitValence" (DZV) (p.ej., para el carbono tendrá una única función contraída para representar el orbital 1s y dos ó tres contraídas para y (TZV) representar cada uno de los orbitales 2s y 2p)
Bases con Polarización
Incluye funciones de polarización, es decir, funciones que no están ocupadas en el estado fundamental del átomo considerado, pero que permiten la polarización de la densidad electrónica hacia las zonas del enlace en la molécula. Estas funciones tienen momento angular superior en uno a la de la última función de la capa de valencia (funciones 2p para el H, 3d para los elementos de la primera fila, etc)
Bases con funciones difusas
Estas añaden funciones con exponentes mucho mas pequeños que los necesarios para representar el átomo neutro. Se llaman difusas porque se extienden a zonas lejanas al núcleo. Sirven para representar la densidad electrónica en aniones, complejos débilmente enlazados y estados excitados de moléculas.
Tabla 1: Tipo de funciones de base atendiendo al número de funciones utilizadas
4.3.
Capas abiertas (UHF). Ecuaciones de Pople y Nesbet.
En general, la construcción de la función de ondas de un sistema con electrones desapareados conviene hacerla a partir de varios determinantes de Slater. Sin embargo, una función monodeterminantal que en muchos casos da buenos resultados, es la obtenida asignando espín α a unos orbitales y espín β
Capítulo I-23
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
a otros sin obligar a los electrones a estar apareados. Esto es, sin obligar a que la parte espacial de cada espinorbital con espín α esté repetida en otro con espín β. Por ejemplo, para un sistema con n=nα+nβ electrones, se puede construir un determinante: ψ = det φ1α (1)φ2α (2)...φnα (nα )φ1β (nα + 1)...φnβ (nα + nβ ) α
β
(78)
Si no se supone φiα = φiβ , es decir si no "apareamos" ningún electrón, se dice que ψ es un determinante no restringido. La obtención autoconsistente de los orbitales mas adecuados para construir un determinante de Slater "no restringido" fue propuesta por Pople y Nesbet (en EE.UU. y, simultáneamente, por G. Berthier en Francia), siguiendo unos pasos paralelos a los que conducen a las ecuaciones de Roothaan en los sistemas de electrones apareados. La obtención de la energía media correspondiente a la función de ondas (78), que puede hacerse de forma parecida a la del caso apareado, conduce a la expresión: ~ E =
1 2 i=1
n
∑
n
n
∑∑
Hii +
i =1
Jij −
j =1
nα
nα
∑∑
Kijα −
i=1 j=1
nβ
Kijβ j=1 nβ
∑∑ i=1
(79)
en la que suma “de 1 a n” debe hacerse sobre todos los orbitales, las sumas de “1 a nα” sólo sobre los orbitales con espín α, y las sumas de “1 a nβ” sólo sobre los de espín β. Las integrales K ijα son sobre orbitales de espín α Kijα =
∫∫ φ
α∗ α i (1)φ j (2)
r r 1 α∗ φ j (1)φiα (2)d r1d r2 r12
(80)
y, análogamente, las integrales K ijβ son sobre orbitales de espín β: Kijβ =
∫∫ φ
β∗ β i (1)φ j (2)
r r 1 β∗ φ j (1)φiβ (2)d r1d r2 r12
(81)
Definiendo una energía monoelectrónica para los orbitales φiα y otra para los φiβ : εiα
=
Hiiα
+
nα
∑ (J
ij
j=1
Capítulo I-24
−
Kijα )
nβ
+
∑J
ij
j=1
(82)
ANTECEDENTES
εiβ
=
Hiiβ
nβ
∑
+
(Jij − Kijβ ) +
j=1
nα
∑J
(83)
ij
j=1
y escribiendo:
φ iα = ∑ C αpi χ p
(84)
φ βi = ∑ C βpi χ p
(85)
p
p
se llega, sin demasiadas dificultades, al sistema de ecuaciones acopladas (a través de los términos Fpq):
∑C
α qi
(Fpqα − ε iα Spq ) = 0
(86)
p
∑ C (F β qi
β pq
− εiβSpq ) = 0
(87)
p
llamadas Ecuaciones de Pople-Nesbet. En ellas: c α Fpq = Hpq +
∑ ∑ [P
− Prsα (ps | rq)
∑ ∑ [P
β − Prs (ps | rq)
rs (pq | rs)
r
c β Fpq = Hpq +
s
rs (pq | rs)
r
s
α Ppq =
nα
∑n C
α α α pjCqj
]
(88)
]
(89)
(90)
j=1
β Ppq =
nβ
∑n C
β β β pjCqj
(91)
j=1
α β Ppq = Ppq + Ppq
(92)
La resolución de las ecuaciones de Pople-Nesbet es, como la de las ecuaciones de Roothaan, iterativa. Se parte de dos conjuntos de coeficientes β 0 α 0 α β { Cpi } y { Cpi } razonables. Se construye con ellos Fpq y Fpq y se resuelven,
independientemente, los dos sistemas de ecuaciones (86) y (87). El primero α (1) } mejorados, y análogamente, da un nuevo conjunto de coeficientes { Cpi β (1) α β (87) suministra otro conjunto { Cpi } . Repitiendo el cálculo de Fpq y Fpq con
Capítulo I-25
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
los nuevos coeficientes se entra en un ciclo que puede repetirse hasta alcanzar el grado de autoconsistencia que se desee. La solución de las ecuaciones (86)+(87) depende de cual sea el apareamiento impuesto a los electrones. En general los sistemas con número par de electrones tienen una energía mucho mas baja si todos están apareados que si no lo están. Pero no siempre ocurre así: una excepción muy conocida es la de la molécula de O2, que a distancias próximas a la de equilibrio posee menos energía con dos electrones desapareados que si se aparean todos los electrones.
5. MÉTODOS POST HARTREE-FOCK En la actualidad es posible obtener, mediante el empleo de bases adecuadas, funciones de onda polielectrónicas muy próximas al límite Hartree-Fock. Conviene, sin embargo, no olvidar que este límite no representa la solución exacta de la ecuación de Schrödinger considerada, sino una aproximación basada en la asignación a priori de una estructura monodeterminantal a la función de onda: ψ(1,2,..., N) = | φ1φ2...φN |
(93)
La aplicación del método variacional con esta estructura para la función de onda demuestra que los orbitales mas adecuados para construir el determinante (93) son los de Hartree-Fock. El modelo de Hartree-Fock tiene, entonces, en cuenta la repulsión interelectrónica a través del potencial ˆ ˆ j(µ) . Ahora bien, el efecto de la repulsión interelectrónica Jj(µ) − K efectivo
∑[ j
]
no queda completamente representado por este potencial efectivo, común para todos los electrones. En un modelo mas preciso, el comportamiento individual de éstos debería verse afectado en alguna medida por la posición instantánea de cada electrón del sistema, y a esto corresponde el significado de la palabra correlación. Las consecuencias de que la función de onda exacta no sea la de Hartree-Fock, se llaman efectos de correlación y la mas importante es, en términos generales, la energía de correlación: Ecorrelación = Eexacta - EHartree-Fock
Capítulo I-26
(94)
ANTECEDENTES
en la que se entiende que Eexacta (normalmente desconocida) no incluye efectos relativistas ni del movimiento nuclear. Por consiguiente no coincide por completo con la energía experimental, aunque normalmente la diferencia sea del todo despreciable. Como hemos dicho, en la aproximación de Hartree-Fock cada electrón interactúa con los demás sólo a través del potencial promedio causado por los n-1 restantes electrones, en lugar de a través de la repulsión electrónica instantánea que implica el potencial de Coulomb. Existen varios métodos desarrollados para tratar de calcular la energía de correlación, que pueden clasificarse de la siguiente forma: •
Interacción de configuraciones (CI);
•
Métodos perturbacionales (MPPT/MBPT);
•
Coupled clusters (CC);
•
Métodos multiconfiguracionales autoconsistentes (MCSCF), incluyendo los casos de SCF de espacio activo completo (CASSCF), o el método de enlace de valencia generalizado (GVB).
5.1.
Interacción de configuraciones
Al resolver las ecuaciones de Hartree-Fock de un sistema con N electrones se obtienen siempre más orbitales φi de los que se necesitan para construir la función de onda del estado fundamental: ψ 0 = | φ1φ2φ3...φN |
(95)
Si se emplean K funciones de base “sobran” (K-N) orbitales, pues el número de soluciones independientes de las ecuaciones de Roothaan es igual al de funciones de base empleadas, mientras que el número de orbitales ocupados es el de electrones del sistema. Los orbitales “sobrantes” reciben el nombre de orbitales virtuales y pueden utilizarse, entre otras cosas, para mejorar la función de onda Hartree-Fock por el método de interacción de configuraciones. Se basa en ir sustituyendo orbitales ocupados en el determinante de Slater por combinaciones de orbitales virtuales donde sus constantes se calculan por método variacional. Cada uno de las posibles estados de los N electrones en los M orbitales moleculares (M>N) recibe el nombre de configuración electrónica y está representado por un determinante de Slater. Si promovemos un único
Capítulo I-27
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
electrón desde un orbital ocupado a uno desocupado, obtenemos, la configuración electrónica de un nuevo estado. Lo que hemos hecho para conseguir esta nueva configuración electrónica es lo que se conoce como una excitación simple. Evidentemente, para un conjunto de N electrones y M orbitales existe un número finito de tales posibles excitaciones simples. Eventualmente podemos excitar un segundo electrón del sistema para obtener una configuración doblemente excitada, o un tercer electrón con lo cual tendremos una excitación triple, etc. El método de interacción de configuraciones completa (full CI) se basa en expresar la función de onda del sistema como una combinación lineal de todas las configuraciones electrónicas posibles del sistema, donde los coeficientes de la combinación lineal se optimizan empleando nuevamente el teorema variacional. La función de onda obtenida da la solución exacta de la ecuación de Schrödinger electrónica en la base considerada, pero tiene el gran inconveniente de que sólo es practicable para sistemas pequeños (esto es, con pocos electrones y/o pocos orbitales virtuales). Ello se debe al crecimiento combinatorio de la cantidad de términos con el nivel de excitación. Una forma usual de disminuir el tamaño de la CI es considerar únicamente algunas de todas las posibles excitaciones. Esto es lo que se conoce como CI truncada o limitada, cuyo ejemplo mas frecuente es la CI que incluye únicamente excitaciones simples y dobles (CISD). El truncamiento del desarrollo que no incluya todas las posibles excitaciones (i.e., que no sea fullCI) conduce a un error que se conoce como “error de inconsistencia con el tamaño” (size consistency). Básicamente, lo que sucede es que la energía CISD (o cualquier CI truncada) de un sistema de N subsistemas idénticos no es igual a la suma de las N energías CISD de cada uno de los subsistemas. Existen métodos que permiten corregir este problema, entre los cuales el mas usado es el método CI cuadrático (QCISD). Los métodos CI en general no son los mas utilizados en el cálculo rutinario de la energía de correlación por varias razones. Fundamentalmente, estas razones son: •
El error debido a la inconsistencia con el tamaño.
•
La lentitud de la convergencia del desarrollo CI.
•
El costo del método, debido a que para el cálculo de la energía necesita el cálculo de las integrales moleculares, que a su vez están relacionadas con las integrales sobre funciones de base
Referencias de este método las podemos encontrar en [Ref. 21].
Capítulo I-28
ANTECEDENTES
5.2.
El método de M/oller-Plesset
El método de interacción de configuraciones no es el único que puede emplearse para tener en cuenta los efectos de correlación. Existen muchos otros, aunque la propia coexistencia de tantos métodos para tratar la correlación pone de manifiesto que ninguno es completamente satisfactorio. El método de Mφller-Plesset, basado en la teoría de perturbaciones de Rayleigh-Schrödinger, resulta atractivo porque conduce a resultados bastante buenos con un esfuerzo de cálculo moderado. Por este motivo, es el que hemos elegido para incluir la correlación en nuestros programas. Para aplicar el método de perturbaciones al cálculo de la energía de correlación se divide el operador hamiltoniano del sistema polielectrónico en dos partes:
ˆ =H ˆHF + H ˆcorr H
(96)
en las que la primera es el hamiltoniano correspondiente al modelo de Hartree-Fock que corresponde al hamiltoniano no perturbado y el segundo término
ˆ corr = H ˆ −H ˆHF H
es
la
perturbación.
Dependiendo
del
orden
de
perturbación que se elija se obtienen los distintos niveles del método de Mφller-Plesset: MP2, MP3, MP4, .... Los métodos perturbacionales presentan la ventaja sobre los de interacción de configuraciones de ser consistentes respecto al tamaño del sistema (size consistency), por lo que resultan particularmente útiles para estudiar la energía de correlación en complejos moleculares respecto a sus componentes. Tiene, sin embargo, dos tipos de defectos: •
No son variacionales, por lo cual no convergen de forma monótona hacia la energía real del sistema, sino que pueden dar resultados por encima o por debajo de la misma.
•
La convergencia de la serie de perturbaciones puede variar mucho de sistema a sistema, incluso en casos en que los sistemas estén íntimamente relacionados (isómeros, por ejemplo). Esto hace que el porcentaje de la energía de correlación recuperado a un nivel determinado no sea el mismo en distintos sistemas, de manera que puede llevar a notorios errores en la predicción de energías relativas.
Finalmente, nótese que el método de Mφller-Plesset puede aplicarse tanto a hamiltonianos RHF (en el caso de capa cerrada) como UHF (para capa
Capítulo I-29
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
cerrada o capa abierta por igual). En el caso de cálculos restringidos de capa abierta (ROHF) no existe una elección única del hamiltoniano de referencia
ˆHF ) y por lo tanto existen varias técnicas alternativas para aplicar la teoría (H de perturbaciones a partir de ellos. Una explicación detallada encontrarse en [Ref. 21].
5.3.
del
método
de
Mφller-Plesset
puede
Coupled Clusters
Los métodos de clusters acoplados coupled clusters (CC) tienen una filosofía similar a la de los métodos perturbacionales, pero aquí el desarrollo es exponencial en lugar de lineal. Así, en lugar de expresar la función de onda como en el método de perturbaciones: Ψ = Ψ(0) + λΨ(1) + λ2Ψ(2) + λ3Ψ(3) + ...
(97)
1 2 1 ˆ3 ˆ+ ˆ Ψ = exp(ˆ T) Ψ0 = (1 T+ ˆ T + T + ...) Ψ0 2! 3!
(98)
se hace de la forma:
donde el operador ˆ T es una suma de operadores de cluster. Cada operador de cluster aplicado a la función de onda de referencia provoca todas las excitaciones de orden i, por lo cual, por ejemplo, la aplicación del operador ˆ T 2 genera todas las configuraciones doblemente excitadas, etc. De la misma forma que sucedía con la CI, realizar un tratamiento completo del problema resulta imposible en cuanto el sistema deja de ser pequeño. Por lo que, la función de onda CC puede truncarse en cualquier punto (es decir, a un cierto orden máximo de excitación). Por ejemplo, la función de onda CC mas frecuente es la que incluye sólo el operador de dobles excitaciones (CCD), en la forma : 1 2 1ˆ3 Ψ = (ˆ 1+ ˆ T2 + ˆ T2 + T2 + ...) Ψ0 2! 3!
(99)
Nótese que debido a los productos de operadores, la función de onda CCD contiene los mismos términos que la CI del mismo orden, pero también términos adicionales: Ψ = Ψ0 +
∑t
i, j,a,b
Capítulo I-30
ab ab ij Ψij
+
∑t
abcd abcd ijkl Ψijkl i, j,k,l,a,b,c,d
+ ...
(100)
ANTECEDENTES
donde los dos primeros términos son los mismos que surgen de un tratamiento CI, pero los términos siguientes están presentes sólo en CCD. Estos términos hacen que CCD, a diferencia de CI, sea consistente con el tamaño del sistema, con lo que elimina uno de los problemas de aquélla. Por otra parte, la fórmula (100) difiere también de lo que se obtendría empleando teoría de perturbaciones, pues en CCD se considera la suma de las dobles excitaciones a orden infinito, incluyendo de hecho excitaciones cuádruples, etc. Por tanto, la energía CCD incluye mucha más energía de correlación que la MP2, por ejemplo, y converge mucho más rápidamente que esta serie. Su único defecto grave respecto a Mφller-Plesset es que resulta mucho mas costosa de calcular, por lo cual las optimizaciones de geometría usando CCD, son mucho menos frecuentes que las obtenidas usando MP2. Los detalles sobre este método lo podemos encontrar también en [Ref. 21].
5.4.
Métodos multiconfiguracionales autoconsistentes
Los métodos multiconfiguracionales autoconsistentes implican la optimización simultánea de los orbitales moleculares dentro de las configuraciones electrónicas consideradas en un desarrollo del tipo CI, y de los coeficientes del propio desarrollo. El caso más conocido de método multiconfiguracional es el llamado SCF multiconfiguracional (MCSCF). En este método se seleccionan ciertas configuraciones electrónicas (por ejemplo, determinantes de Slater) formadas a partir de una cierta configuración de referencia. A menudo estos determinantes se seleccionan en base a una cierta "intuición química" ya que, obviamente, el tamaño del cálculo a realizar lo hace impracticable mas allá de ciertos límites. A continuación, se optimizan, utilizando el teorema variacional, los coeficientes de la combinación lineal de determinantes, reoptimizando simultáneamente en cada ciclo los coeficientes de los orbitales moleculares individuales (en lugar de fijarlos en sus valores HF como se hace en la CI). La función de onda MCSCF resulta útil para tratar los problemas que surgen de la correlación no dinámica, por ejemplo, la disociación homolítica de enlaces. Una forma común de elegir los determinantes que participan en la MCSCF es realizar lo que se conoce como SCF completo en el espacio activo (complete-active-space SCF, CASSCF). En este caso los determinantes se eligen de forma que se incluyan todas las excitaciones posibles dentro de un subespacio del número total de orbitales, que se conoce con el nombre de "espacio-activo". La elección de los orbitales que se incluyen dentro de este
Capítulo I-31
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
espacio activo es responsabilidad del investigador. Este método ha sido desarrollado por Roos y colaboradores [Ref. 22]. Otra forma de MCSCF es la conocida con el nombre de enlace de valencia generalizado (generalized valence bond, GVB) desarrollado por Goddard y colaboradores [Ref. 23]. En este caso el MCSCF es realizado incluyendo un número muy pequeño de configuraciones electrónicas, con el propósito de describir bien un estado electrónico que es inherentemente multiconfiguracional (como un singlete de capa abierta, por ejemplo) o la disociación homolítica de un enlace. El nombre surge porque usualmente las configuraciones electrónicas a incluir se seleccionan mediante un análisis de enlace de valencia de las distribuciones electrónicas. Finalmente debemos mencionar brevemente la existencia de métodos que permiten aplicar una CI no a una única configuración de referencia sino a un conjunto de configuraciones de referencia seleccionadas de tal manera que representen correctamente el proceso químico de interés. Este tipo de procedimiento se conoce con el nombre de CI multireferencia y, usualmente, se concreta en la realización de un cálculos CASSCF o GVB para determinar las configuraciones de referencia importantes en el problema que se estudia, seguido de una CISD incluyendo las configuraciones electrónicas obtenidas de excitar electrones a partir de todas las configuraciones de referencia. Este procedimiento es extremadamente costoso en comparación con los mencionados anteriormente. Referencias de este método las podemos encontrar en [Ref. 21].
6. PROPIEDADES MOLECULARES Una de las aplicaciones químicas más importantes de la Mecánica Cuántica es que permite calcular por vía teórica las propiedades moleculares, y permite comprender tanto la estructura como la reactividad de las moléculas. Revisaremos a continuación, muy someramente, algunas de las propiedades moleculares más adecuadas para su cálculo teórico, que nos han servido como “banco de pruebas” de los métodos de cálculo desarrollados en esta tesis. Dentro de las propiedades moleculares susceptibles de ser calculadas teóricamente cabe distinguir dos clases, según estén basadas en cálculos de las energías de las moléculas, o estén basadas en el cálculo de su distribución electrónica.
Capítulo I-32
ANTECEDENTES
Entre las propiedades “de raíz energética” cabe resaltar: -
La geometría molecular (distancias y ángulos)
-
Las energías de atomización y disociación
-
Las energías de ionización y electroafinidades
-
Las energías de excitación electrónica
-
Las energías de excitación vibracional y rotacional
-
Etc...
Entre las propiedades relacionadas con la distribución de la carga eléctrica de las moléculas cabe resaltar: -
Sus momentos dipolares y multipolares
-
Sus polarizabilidades
-
Sus susceptibilidades magnéticas
-
Determinados efectos espectroscópicos de la aplicación de campos eléctricos y magnéticos
-
Etc...
A continuación revisamos brevemente cual es el enfoque teórico del cálculo de las propiedades citadas.
6.1.
Propiedades energéticas
El cálculo de las distancias y ángulos de enlace se lleva a cabo minimizando la energía total calculada resolviendo la ecuación de Schrödinger electrónica. Como es bien sabido una molécula puede “romperse” siguiendo distintos canales. Uno de ellos es separarse en los átomos neutros que la componen, llamándose energía de atomización a la energía necesaria para esta ruptura. La energía de atomización puede calcularse como diferencia entre la energía total de la molécula y la energía de cada uno de los átomos que la componen. Otro canal es la formación de iones, o de átomos en estados excitados. A esta energía, que depende de cuales sean los productos, se la denomina energía de disociación y podemos calcularla mediante la
Capítulo I-33
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
diferencia de la energía total de la molécula y las energías de los productos de la disociación. La energías de ionización (positivas o negativas, esta última llamada electroafinidad) podemos determinarlas restando las energías totales calculadas para el sistema neutro y para el ionizado, o bien aplicando la aproximación de Koopmans. En este último caso, la energía de ionización positiva coincide con la energía del orbital ocupado más alto (HOMO). La energía de ionización negativa o electroafinidad coincide, según la aproximación de Koopmans con la energía del primer orbital vacío (LUMO).
6.2.
Propiedades eléctricas y magnéticas
El momento dipolar puede calcularse directamente conociendo las posiciones de los núcleos y la función de onda electrónica. En la aproximación Hartree-Fock-Roothaan: µ = 2.5416 donde las integrales
µ|r | ν
∑Z r − ∑P
µν
aa
µν
a
µ|r | ν
(101)
expresan el valor medio del operador de
posición respecto a las funciones de base empleadas y los elementos Pµν son los elementos de la matriz densidad, obtenidos a partir de los coeficientes de los orbitales moleculares: Pµν =
∑n C i
iµCiν
(102)
i
El momento dipolar puede considerarse como el primer término de un desarrollo del campo eléctrico generado por la molécula de la forma r a a a V(R) = 1 + 22 + 33 + ... R R R
(103)
r en el que R es el punto en que se mide el potencial creado por la molécula y r R el módulo de R (que da idea de la distancia del punto considerado a la molécula). En moléculas neutras: a1 = 0
y
r r a2 = µ x R
(104)
r siendo µ el momento dipolar. Este da una idea, por tanto, de la magnitud de las interacciones que esta molécula realizará con sus vecinas. El término que
Capítulo I-34
ANTECEDENTES
le sigue, representado por a3, es el momento cuadrupolar (cuya consideración es importante para aquéllas moléculas que por ser muy simétricas tienen momento dipolar nulo). Los programas usuales de cálculo permiten también obtener momentos multipolares de mayor orden (octupolo, hexadecapolo, etc) pero debe tenerse cuidado porque los momentos dipolares de orden alto pues, a diferencia del momento dipolar, son dependientes del origen de coordenadas. En realidad, la propiedad eléctrica clave de las moléculas es el potencial electrostático que crean a su alrededor, esto es la energía de interacción entre el sistema de cargas de la molécula y una carga puntual unitaria localizada en el punto considerado: r V(R) =
∑
r r −1 Za R a − R −
a
∑P
µν
µν
1 µ| r r |ν r −R
(105) r r
El cálculo del potential electrostático molecular sobre una superficie que rodea a la molécula (definida, por ejemplo, por una cierta densidad electrónica constante) da idea de como la molécula puede interaccionar con moléculas polares o especies cargadas. Un concepto importante en el estudio de moléculas es el de carga electrónica de los átomos de una molécula. No existe una definición unívoca de como asignar cargas electrónicas a los átomos dentro de una molécula, dado que toda partición de la densidad electrónica entre los distintos núcleos es, hasta cierto punto, arbitraria. El modelo usado desde más antiguo para particionar la densidad electrónica es el Análisis de Población de Mulliken. En este esquema, la densidad electrónica correspondiente a la interacción entre dos orbitales dados se divide simplemente por igual entre ellos, de manera que la "carga" asignada a un cierto orbital es: qµ = Pµµ +
∑P
µν
µ ν
(106)
µ≠ν
La suma de las cargas de todos los orbitales asociados a un cierto átomo dan la población electrónica de ese átomo en la molécula considerada. El análisis de población de Mulliken tiene inconvenientes importantes. Por ejemplo, es sumamente dependiente del conjunto de funciones de base empleado y, con bases extendidas, puede generar situaciones irreales, tales como cargas mayores de 2e en un orbital determinado. Esto ocurre cuando las funciones de base centradas sobre un cierto átomo están siendo usadas en realidad para describir la densidad electrónica en las cercanías de otro núcleo.
Capítulo I-35
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
Una solución es realizar el análisis empleando los orbitales naturales, que se obtienen diagonalizando la matriz densidad. Otra forma de determinar las cargas, que es menos dependiente de las condiciones del cálculo que el análisis de poblaciónes de Mulliken, consiste en ajustar las poblaciones de manera que reproduzcan el potencial electrostático molecular calculado en una cierta zona alrededor de la molécula. Este es un procedimiento más costoso que el de Mulliken, pero más adecuado si se quieren determinar cargas atómicas para usar en métodos de Mecánica Molecular. La importancia de las propiedades magnéticas en el campo de la Química se debe en gran medida a que constituyen el fundamento de las espectroscopias de RMN, que tanta utilidad poseen en la elucidación de la estructura de las moléculas, sobre todo orgánicas y biológicas. El cálculo teórico de estas propiedades se realiza a través de la evaluación de los apantallamientos y de las constantes de acoplamiento por procedimientos que pueden estudiarse, por ejemplo, en [Ref. 24].
7. INTEGRALES DE LOS CÁLCULOS ATÓMICOS Y MOLECULARES 7.1.
Clasificación
Existen varias formas de clasificar las integrales implicadas en los cálculos atómicos y moleculares. Si tenemos en cuenta el número de variables, podemos clasificarlas en monoelectrónicas y bielectrónicas, según incluyan las coordenadas de uno o de dos electrones respectivamente. En el caso de considerar el tipo de operador matemático incluido, tendremos integrales cinéticas si contiene el operador
∇ i2 , de atracción nuclear si es el
ZA 1 , y de repulsión electrónica si es el . Otro criterio de clasificación, rij rAi importante al tratar moléculas, es el número de centros implicados. Según este criterio podemos encontrar integrales monocéntricas, bicéntricas, tricéntricas y tetracéntricas. También se pueden distinguir entre integrales canónicas que son aquellas que tienen los ejes de las funciones de base alineados de forma que los ejes OZ coincidan y los OX y OY sean paralelos, y no canónicas cuando están rotados con respecto a la disposición canónica. Resumiendo los criterios indicados al principio, la clasificación de las integrales moleculares puede resumirse en una tabla como la Tabla 2.
Capítulo I-36
ANTECEDENTES
de solapamiento monocéntricas bicéntricas monocéntricas Monoelectrónicasde energía cinética bicéntricas monocéntricas de atracción nuclear bicéntricas tricéntricas monocétricas bicéntricas Bielectrónicas tricéntricas tetracéntricas Tabla 2: Clases de integrales moleculares
7.2.
Integrales monoelectrónicas
Las integrales monoelectrónicas que aparecen el los cálculos atómicos y moleculares son de la forma:
) I pq = ∫ χ ∗p ⋅ Ω ⋅ χ q dτ
(107 )
)
)
donde Ω es un operador. Cuando el operador Ω es la unidad las integrales
−
∇2 son integrales de 2
Z rAµ
son integrales “de
se llaman de solapamiento Spq, cuando es el operador
−
energía cinética Tpq, y cuando es el operador atracción nuclear” Vpq(A). Spq =
r
∫ χ (µ) ⋅ χ (µ)dr ∗ p
q
(108)
µ
Tpq = −
1 2
∫ χ (µ) ⋅ ∇ χ (µ)dr
(109)
VA,pq =
∫
χp∗ (µ) ⋅
r ZA χq(µ)d rµ rAµ
(110)
∗ p
2 µ q
r
µ
Capítulo I-37
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
Las transformar,
integrales
de
expresando
energía el
cinética
∇
operador
2 µ
monocéntricas en
coordenadas
se
pueden
polares
y
suponiendo que la función de base la podemos expresar como una parte radial R(r) por un armónico esférico, como:
Tpq
(
)
∞ lp lp + 1 1 = − δmpmq δl p l q Rp* r 2R 'q' + rR 'q dr + 2 2 0
∫ (
De forma análoga se calculan monocéntricas, resultando:
)
las
integrales
de
atracción
(111)
nuclear
∞
∫
VM,pq = −ZA δmpmq δl p l q Rp*R qrdr
(112)
0
7.3.
Integrales bielectrónicas Son de la forma:
(pq rs)=∫ ∫ χpχq(1) r1
12
r r χr χs(2) d r1d r2
(113)
Para resolverlas cuando son monocéntricas y las funciones de base son de la forma R(r) Ylm(θ, ϕ) se puede sustituir el factor r12-1 por su desarrollo: 1 = r12
∞
4π
r
l =0
∗ l ,m
(1)Υl,m(2)
(114)
lo que da lugar a la expresión:
(pq rs) = ∫ ∫ χ χ p
q
1 r r χr χsd r1d r2 = 4πNpNqNrNs ⋅ r12
l 1 r R pR q −l l = 0 2l + 1 ∞
l
(115)
La primera integral corresponde a la parte radial y las segunda y tercera a la parte angular. La integral radial se descompone en dos sumandos con el fin de separar las variables. El primero para r2< r1 y el segundo para r1< r2
Capítulo I-38
ANTECEDENTES
∞∞
∫∫
RpR q (1)
00
r< r>
∞ r1
=
l
l +1
RrR s (2)r12r22dr1dr =
∞ ∞ l r RpR q l +1 RrR sr12r22dr1dr + RpR q 1l +1 r1 r2 0 r2 =0 0 r2 =r1
r2
∫∫
l
∫ ∫
(116) RrR sr12r22dr1dr
La primera parte angular se calcula, expresando Ap y Aq como combinación lineal de armónicos esféricos:
∫
IA1 = apaq Υl
∫
mq
+ aqbp (−1)
Υl
pmp
Υl
p −mp
qmq
Υl
∫ dΩ +b b ∫ (−1)
mq
Υl∗ m dΩ + apbq (−1)
qmq
Υl∗ m
Υl
pmp
mp +mq
p q
Υl
Υl
q −mq
p −mp
Υl∗ m dΩ +
Υl
q −mq
Υl∗ m dΩ
(117)
Cada una de estas cuatro integrales se obtienen utilizando las expresiones siguientes [Ref. 25]:
∫∫
1/2
Υl
Υ Υl∗mdΩ 1m1 l 2m2
(2l1 + 1)(2l 2 + 1) = 4π(2l + 1)
C(l1l 2l;m1m2m)C(l1l 2l;000) ⋅ δm,m1 +m2 (118)
donde C(l1l2l; m1m2m) son los coeficientes de Clebs-Gordan, de Wigner [Ref. 26]. La segunda integral angular es similar a la primera. Cuando las funciones de base son GTO el desarrollo (114) no resulta adecuado pues conduce a integrales (116) que no tienen resolución analítica, ni permiten una resolución numérica eficiente. Cuando se emplean GTO hay que recurrir a una transformada de Laplace de r12-1: 1 2 = r12 π
∞
∫ exp(−r 0
u2 ) du
2 12
(119)
que introduce una nueva variable de integración, pero conduce a integraciones analíticas (en monocéntricas) o numéricas eficientes (en policéntricas). El desarrollo completo de las integrales con GTO lo veremos en el Capítulo III de esta tesis, que sigue muy de cerca el dado en [Ref. 27].
Capítulo I-39
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
8. INTEGRACIÓN MEDIANTE MÉTODOS MONTECARLO Y CUASIMONTECARLO 8.1.
Introducción.
La esencia del método Montecarlo es la experimentación con números aleatorios. El procedimiento usado consiste en diseñar juegos de azar con estos números, esperando obtener de su observación conclusiones útiles para la resolución del problema que se esté estudiando. Dentro de los métodos Montecarlo podemos distinguir dos tipos, uno probabilista donde se simulan con números aleatorios fenómenos que son realmente aleatorios, y otro determinista donde se resuelven problemas que no son aleatorios en la realidad. El cálculo numérico de integrales definidas que se realiza en esta tesis es un ejemplo de Montecarlo determinista. La integración Montecarlo no ha sido demasiado usada hasta el momento en problemas de Química Cuántica debido a la dificultad de obtener resultados suficientemente precisos con tiempos de cálculo aceptables. A pesar de ello señalaremos tres ventajas importantes de estos métodos: -
Su precisión es independiente del número de dimensiones de la integral considerada.
-
Su uso no depende del tipo de integrando, por lo que se puede usar con funciones de base no tratables con métodos analíticos.
-
Quedan favorecidos por la evolución de los medios de cálculo con los que se aplican.
Existen dos procedimientos para calcular las integrales definidas por el método de Montecarlo. El primero, llamado método de éxito-fracaso, está basado en interpretar la integral como un hipervolumen n-dimensional donde su valor viene dado por la probabilidad de encontrar un punto aleatorio dentro o fuera de dicha área. El segundo, llamado método de la media muestral, está basado en la definición del valor medio de una variable aleatoria continua [Ref. 28]. Este segundo método es el que tendremos en cuenta en el desarrollo de esta Tesis. En trabajos anteriores [Ref. 29] se ha podido comprobar que los problemas de eficiencia de este método limitaban muy drásticamente el tamaño de las bases utilizables. Por una parte el tiempo de cálculo consumido por las integraciones Montecarlo es muy largo, pues depende del cuadrado de
Capítulo I-40
ANTECEDENTES
la precisión deseada, la cual no puede ser pequeña si se pretende calcular orbitales atómicos o moleculares. Por todo ello se ha decidido investigar nuevos conjuntos de números que, aunque no sean aleatorios, proporcionan buenas estimaciones para las integrales. Estos métodos pseudoaleatorios, llamados métodos Cuasimontecarlo, se puede usar para la integración numérica. El esquema de aplicación es totalmente similar al de Montecarlo, distinguiéndose solamente en la forma en que se generan los puntos en que hay que evaluar el integrando. En este caso la generación no pretende imitar distribuciones aleatorias, sino que se usan números deterministas, diseñados especialmente para que estén “bien distribuidos”. Los metodos Cuasimontecarlo han contribuido también a mejorar las estimaciones obtenidas bajo el esquema de Montecarlo, no solo en el campo de la integración numérica, sino también en los problemas de simulación y optimización. El nombre del método fue usado por primera vez por Richtmyer en 1951. Posteriormente tenemos trabajos de Zaremba, Sugihara, Murota, Neiderreiter... En el apartado de los “puntos bien distribuidos” , que han sido usados en las últimas fases de este trabajo, hemos contado especialmente con las contribuciones de Korobov, Halwka, Zaremba y Neiderreiter.
8.2.
Integración Montecarlo
La evaluación de las integrales definidas por el Método de Montecarlo se basa en el siguiente teorema [Ref. 30]: Sean x1, x2 ,..., xn... variables aleatorias independientes, idénticamente distribuidas, con función de densidad f(x). Si gi son funciones de xi, entonces: G=
N
1 N
∑ g (x ) i
(120)
i
i =1
es una variable aleatoria que verifica: E(G) =
var(G) =
1 N
N
∑ E(g (x )) i
(121)
i
i =1
1 N2
N
∑ var[g (x )] i
i
(122)
i=1
En particular, cuando todas las g(xi) son idénticas, e iguales a g(x), se tiene que:
Capítulo I-41
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
E(G) =
1 N
1
N
∑ E(g (x )) = N NE(g(x)) = E(g(x))
(123)
1 1 N var g(x) = var [g(x)] 2 N N
(124)
i
i
i =1
var(G) =
En virtud de la definición de esperanza matemática de g(x), la primera de estas expresiones, puede escribirse en la forma: 1 E(G) = E N
N
∑ 1
g(x i ) =
∞
∫ g(x)f (x)dx = E(g(x))
(125)
−∞
Este resultado es el fundamento de la integración Montecarlo con función densidad, que consiste en muestrear una serie de números aleatorios xi con función de densidad f(x) y usarlos para calcular la media de los valores obtenidos para g(x). Esa media es una estimación de la integral, con varianza que decrece a medida que aumenta el número de puntos empleados, según se deduce de la expresión (124) para var(G). El error estándar ε de la estimación Montecarlo de una integral es la raíz cuadrada de la varianza, y toma el valor [Ref. 31]: ε=
var g σ = N N
(126)
lo que indica que para obtener una integral con un error estándar menor que una cierta cota ε prefijada se deberá elegir un N que cumpla la condición:
N≥
σ2 ε2
(127)
var g ≤c Nc
(128)
La desigualdad de Tchebycheff: P G − G ≥
suministra además una cota para la probabilidad de obtener un error mayor que el propuesto en la estimación del valor de la integral, pudiéndose siempre disminuir este error sin más que aumentar el valor de N. Pero esta forma de disminuir el error es lenta. El camino para aumentar la eficiencia del método MC en el cálculo de las integrales pasa entonces por la elección de una descomposición del integrando en un producto que sea adecuado para que la
Capítulo I-42
ANTECEDENTES
constante σ resulte lo menor posible (Integración MonteCarlo con función densidad). Si se quiere evaluar la integral n-dimensional:
∫
G=
Ω
g(X)dX
(129)
se introduce una función f(X) que sea positiva en el dominio Ω y cuya integral valga 1, para asegurar que f(X) sea una función de densidad. Entonces: G=
En principio
∫
Ω
g(X) f(X)dX f(X)
(130)
g(X) debería ser finita , pero el método funciona también en f(X)
casos donde haya una cantidad numerable de puntos con
g(X) infinito. f(X)
En este caso se evalua la función g(X)/f(X) en N valores de X que se muestrean con función de probabilidad f(X). La varianza es en este caso es: var(G) =
∫
g(X)2 f(X)dx − G2 f(X)2
(131)
La función f que minimiza la varianza es precisamente proporcional a la función g. Sin embargo emplear una función densidad así no tiene sentido, ya que para ello se requiere conocer la integral de g, y ese problema es justamente el que estamos tratando de resolver por Montecarlo. En los casos prácticos, por tanto, se elige una función que, siendo fácil de muestrear, sea tan parecida a g como sea posible.
8.3.
Métodos Cuasimontecarlo (ideas básica)
El método de integración Cuasimontecarlo tiene en común con el de Montecarlo que la integral se estima como el valor medio del integrando en un cierto conjunto de puntos, y se diferencia de aquel por la forma en que son elegidos estos puntos: -
En el método Montecarlo son elementos de una muestra con distribución aleatoria en el dominio de integración, obtenida a partir de alguna rutina de generación de números pseudoaleatorios.
Capítulo I-43
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
-
En el caso de los métodos Cuasimontecarlo los conjuntos de puntos no pretenden ser aleatorios, sino que son construidos explícitamente con el objeto de obtener buenas integraciones desde un punto de vista general.
Las sucesiones de puntos que hemos ensayado para la resolución de las integrales moleculares, una vez hecho los cambios de variable necesarios para transformar todos los intervalos de integración en (0,1) y emplear las funciones de peso necesarias, han sido las sucesiones de Halton [Ref. 32], Korobov [Refs. 33 y 34] y conjuntos de puntos de Hammersley [Refs. 33 y 35]. En estudios previos [Refs. 28, 37, 38 y 39] se ha comprobado la superioridad de las técnicas de integración Cuasimontecarlo sobre las Montecarlo ordinarias evaluando diversas clases de integrales, y obteniendo las energías de un conjunto de átomos y alguna molécula. De hecho para esta clase de cálculos, el método Cuasimontecarlo (en cualquiera de las versiones ensayadas) es entre diez y cien veces más eficiente que el de Montecarlo.
Capítulo I-44
Capítulo II: INTEGRACIÓN CON FUNCIONES DE BASE DE SLATER
INTEGRACIÓN CON FUNCIONES DE BASE DE SLATER
Capítulo II: INTEGRACIÓN CON FUNCIONES DE BASE DE SLATER 1. INTRODUCCIÓN 2. ARMÓNICOS ESFÉRICOS 3. INTEGRALES MONOCÉNTRICAS 4. INTEGRALES POLICÉNTRICAS 5. RESULTADOS COMPARATIVOS
1. INTRODUCCIÓN Las integrales atómicas y moleculares pueden calcularse analíticamente sólo para unas pocas clases de funciones de base. El objetivo primario de esta tesis ha sido poner a punto métodos numéricos que permitan emplear cualquier tipo de funciones de base. Evidentemente, para nosotros era de gran interés poder comparar nuestros cálculos con otros analíticos, en algún caso en que éstos fuesen posibles. Para que la comparación sea significativa es muy conveniente, además, llevarla a cabo con programas que solo difieran en el cálculo de las integrales. Por este motivo hemos realizado el estudio y la programación de las integrales analíticas que se presentan en este capítulo, aunque sea fácil encontrar programas estándar para su cálculo. El cálculo analítico de las integrales Spq, Hpq y (pq|rs) monocétricas con funciones de base de Slater es factible, y el objeto de este capítulo es resumir sus fundamentos y detallar un programa capaz de llevarlos a cabo y usarlos en la obtención de orbitales de átomos polielectrónicos. Las funciones de Slater reales pueden ser de tres clases: χi = Nir n e − α r Pl0 i
i i
i
(1)
Capítulo II-45
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
~ cos miϕ
(2)
~ m ~ χi = Nir n e − α r Pl senmiϕ
(3)
χi = Nir n e − α r Pl i
i i
i
~ mi i
i
i i
i
donde Pl
m
son las funciones asociadas de Legendre con
cos θ como variable,
y Ni es la constante de normalización completa de la función, que puede descomponerse como producto de tres constantes de normalización parciales:
N = Nr Nθ Nϕ
(4)
correspondientes a las variables r, θ y ϕ . Se ha utilizado el índice
~ para m
designar las funciones reales, con el convenio: ~ mi > 0 ~ mi < 0
( ) ( )
~ χi incluye cos mϕ ~ χi incluye sen mϕ
⇒ ⇒
(5)
Con objeto de simplificar en lo posible las expresiones, usaremos en lo que sigue las funciones de Legendre normalizadas con la condición :
∫ (N π
0
θ
(2l + 1) (l − Nθ = 2 (l +
)
~ m l
2
P (cos θ) sen θ dθ = 1 ⇒
~ )! 1 / 2 m ~ )! m
(6)
~
y las designaremos como Plm siendo: ~
~
Plm = Nθ Plm
(7)
N = Nr Nϕ
(8)
R (r ) = r n e − αr
(9)
Denotaremos por:
~
Plm ~ ~ ϕ A l , m~ (θ, ϕ) = Plm sen m i m~ ~ Pl cos m i ϕ Por tanto, las funciones de Slater antedicha, de la forma:
Capítulo II-46
~ =0 si m ~ >0 si m
(10)
~
l
l +1
R r R r r dr1 dr + 2 2 s 1 2
R r R s (2)r12r22 dr1 dr = ∞
∞
∫ ∫
0 r2 =r1
R pR q
r1 r2
l
l +1
(47) R r R r r dr1 dr 2 2 s 1 2
Ambas se resuelven aplicando reiteradamente la expresión, válida si n ∈Z+: n ax ∫ x e dx =
(−1)nn! eax n nxn−1 n(n − 1)xn−2 x − + − ... an a a2 a
Se obtiene, para la primera el valor:
Capítulo II-54
(48)
INTEGRACIÓN CON FUNCIONES DE BASE DE SLATER
∞
γ R0 = ∫ r1
np + nq − l +1
[
exp − (α p + α q )r1
r1
]∫r
nr + n s + l + 2
exp[− (α r + α s )r2 ] dr2 dr1
(49)
nr +ns − l +1
exp[− (α r + α s ) r2 ]dr2 dr1
(50)
2
r2 = 0
0
y para la segunda: ∞
γ R∞ = ∫ r1
np +nq + l +2
[
exp − (α p + α q )r1
∞
] ∫r
2
r2 =r1
0
Si llamamos n1 = (np + nq - l) y n2 = (nr + ns + l), se obtiene como fórmula general para cada una de las integrales radiales: γ R0 =
(n2 + 2)!(n1 + 1)! n + 3 (n2 + 2)!(n1 + n2 + 4 − k)! −∑ n + n + 5 +k k αn2 + 3 ⋅ α1n + 2 k =1 (n2 + 3 − k )! α 2 (α1 + α 2 ) 2
2
1
1
2
(51) γ R∞ =
n1 + 2
∑ k =1
3.4.2.
(n1 + 1)!(n1 + n2 + 4 − k)! αk2 (α1 + α2 )n +n +5 +k 1
2
Parte angular La primera parte de la integral angular se calcula análogamente a
como se hizo para las integrales de momento:
γ A (p, q ) = a p a q ∫ Υl p m~ p Υl q m~ q Υl∗ m~ dΩ + a p b q ∫ (−1) + a q b p ∫ (−1)
~ m q
Υl p − m~ p Υl q m~ q Υl∗ m~ dΩ + b p b q ∫ (−1)
~ m q
~ +m ~ m p q
Υl p m~ p Υl q − m~ q Υl∗ m~ dΩ + Υl p − m~ p Υl q − m~ q Υl∗ m~ dΩ
(52)
y resolviendo estas cuatro integrales igual que se hizo para las de momento. La segunda integral angular ( γ A (r , s )) es, como se aprecia fácilmente, similar a la anterior, y se reduce a ésta sin más que transformar conjugado, utilizando nuevamente la expresión
Υl , m en su armónico ~
Υl∗, m~ = (−1) m Υl , − m~ .
Una vez realizado lo anterior, basta aplicar la relación de ortonormalidad (17) de los armónicos esféricos para transformar la suma infinita de la expresión (46), donde l varía desde 0 a infinito, en una suma finita, ya que solo los valores de l que no cumplan: max ( l p + l q , l r + l s ) < l < min( | l p - l q |, | l r - l s | )
(53)
conducirán a términos nulos por ortonormalidad.
Capítulo II-55
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
4. INTEGRALES POLICÉNTRICAS Aunque hace tiempo que es posible calcular directamente las integrales bicéntricas con funciones de base STO [Ref. 40], y recientemente se han desarrollado procedimientos para calcular también las integrales tri y tetracéntricas [Ref. 41], nosotros hemos optado por calcular todas las integrales policéntricas vía desarrollos STO-nG con valores altos de n [Refs. 42 y 43]. El motivo de esta elección no es único, pero principalmente se debe a que nuestro interés último no es el manejo de las funciones STO, sino el de funciones generalizadas. Estas pueden integrarse por cuasimontecarlo (CMC) como hacemos en esta tesis, pero también podrían tratarse por desarrollo en n gaussianas y mediante procedimientos híbridos de desarrollo nG y CMC. Eligiendo los desarrollos nG para el tratamiento de las integrales con base STO, dejamos nuestro programa preparado para implementar más adelante en él los procedimientos nG-CMC que acabamos de sugerir.
4.1.
Desarrollos STO-nG
La alternativa de tipo habitual más próxima a las integraciones CMC que estudiamos en esta tesis son los desarrollos gaussianos, ampliamente utilizados en el manejo aproximado de las funciones de base de Slater. Los métodos STO-nG fueron propuestos por primera vez por Stewart [Refs. 44 y 45] y ampliamente utilizados por J. A. Pople y sus colaboradores [Ref. 46]. Estos desarrollan bases STO-nG con n comprendido entre 2 y 6 gaussianas, que luego fueron ampliados a valores de n mucho mayores [Refs. 42 y 43]. Los desarrollos STO-nG con n suficientemente grande (del orden de 18, o más) permiten calcular las integrales sobre STO con una precisión inmejorable (en el sentido de que los errores de integración llegan a tener menos importancia que los demás errores, como los originados en el proceso de diagonalización, etc.). Nosotros hemos empleado los desarrollos nG en tres contextos: a) Como método de contraste con el resto de nuestros cálculos. b) Para incluir cálculos directos con bases STO en nuestro programa. c) Para desarrollar mediante gaussianas las funciones de base recortadas que hemos estudiado en una parte de esta Tesis. Los desarrollos STO-nG se basan en la idea de representar las funciones Φ a desarrollar en la forma:
Capítulo II-56
INTEGRACIÓN CON FUNCIONES DE BASE DE SLATER
Φ≈
∑C G µ
µ
~ ≡Φ
(54)
µ
en la que Gµ representa funciones de base dependientes de algún parámetro α µ . Los valores óptimos de los parámetros α µ y de los coeficientes Cµ pueden determinarse, tal como hizo Stewart ya en 1969 [Refs. 44 y 45], aplicando la condición de mínimos cuadrados:
∫ (Φ − Φ ) dr = MÍNIMO ~2 r
(55)
A partir de esta idea básica se puede continuar por caminos ligeramente diferentes. Stewart desarrolló cada STO 1s, 2s, 3s, 4s, 2p, 3p, etc. de forma independiente, empleando como base funciones gaussianas normalizadas: 1
(2 / ξ )n +1 2 n−1 Φ= r exp(−ξr) Ylm(θ, ϕ) (2n)! 1
2 12 2 22l + 3 nG 2l +3 π l ~ Φ = r Y ( θ , ϕ ) ⋅ dk α 4 exp(−αkr 2 ) ∑ l m (2l + 1)!! k =1
(56)
En nuestro programa hemos seguido el procedimiento STO-nG propuesto por Fernández Rico et al en [Refs. 42 y 43]. En este caso se han utilizado como base las funciones gaussianas sin normalizar: 1
(2 / ξ )n+1 2 2E n−2l −1 l exp(−r) si n − l impar r Φ= r Ylm(θ, ϕ) r exp(−r) si n − l par (2n)!
(57)
nG
~ Φ = Nr r l Ylm(θ, ϕ) ∑ Ck exp(−αkr 2 ) k =1
donde E(x) es la parte entera de x, y Nr es la constante de normalización ~ radial para el desarrollo de la función Φ , dada por:
Nr
2
∫
∞
0
2E n− l −1 r 2 r l exp(−r) si n − l impar r 2dr = 1 r exp(−r) si n − l par
(58)
Con objeto de comparar los cálculos realizados y poder utilizar desarrollos superiores a STO-6G con el programa Gaussian, es necesario buscar una relación entre las constantes que utilizan Stewart et al y las que utilizan Fernández Rico et al. Esta viene determinada por la fórmula:
Capítulo II-57
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
1 2 2 22l +3 π (2l + 1)!!
1
2
dk α
2 l +3
4
= NrCk
(59)
El manejo de las funciones de base gaussianas con armónicos esféricos es algo engorroso, de manera que éstas han sido progresivamente sustituidas por las gaussianas cartesianas: ny (µ) nz (µ) − α µr 2
Gµ = N xnx (µ)y
z
e
(60)
La única ventaja de las gaussianas cartesianas es que facilitan el cálculo de las integrales finales, pero se trata de una ventaja lo suficientemente importante como para hacer que los programas de cálculo más habituales hagan uso de ellas.
4.2.
Integrales policéntricas
Las integrales monoelectrónicas correspondientes a las funciones de base gaussianas son del tipo: I=
∫
r r ˆ (r )χ d r = N N χp*O 1 q p q
r
r
∑ ∑ C C ∫ G Oˆ (r)G dr p µ
µ
q ν
p µ
1
q ν
(61)
ν
r ˆ (r ) es un operador matemático monocéntrico. Obsérvese que las donde O 1 funciones gaussianas pueden estar centradas en el mismo centro (integrales monocéntricas) o bien en centros distintos (integrales bicéntricas). En el caso r ˆ (r ) esté centrado en otro lugar distinto a los de las en el que el operador O 1
gaussianas tendremos integrales tricéntricas. Las integrales bielectrónicas son del tipo: * r p 1
r
r r
r
r
r r
∫ χ (r )χ (r ) Oˆ (r , r ) χ (r )χ (r )dr dr = r r r r r r r r ˆ (r , r ) G (r )G (r ) d r d r = N N N N ∑ ∑ ∑ ∑ C C C C ∫ G (r )G (r ) O
I=
p q r
q 1
2
1
* r
2
p µ
s
µ
ν
ρ
2
s
q r ν ρ
s σ
2
1
p µ
1
2
q ν 1
2
1
2
r ρ
2
s σ
2
1
(62)
2
σ
r r ˆ (r , r ) es un operador que puede depender de las coordenadas de donde O 2 1 2 ambos electrones. Estas integrales pueden ser monocéntricas, bicéntricas, tricéntricas o tetracéntricas, dependiendo de que el origen de coordenadas tenga cada función de base.
Capítulo II-58
INTEGRACIÓN CON FUNCIONES DE BASE DE SLATER
Cuando se recurre a desarrollos STO-nG estas integrales quedan descompuestas en la combinación lineal de integrales sobre gaussianas:
(µ ν |ρ σ) = ∫ G G (1) 1 µ
ν
r12
r r GρGσ (2) d r1d r2
(63)
que solo pueden ser monocéntricas o bicéntricas, debido a que el producto de gaussianas de distinto centro es siempre otra gaussiana centrada en un punto intermedio, tal como se comprobará en el Capítulo III.
5. RESULTADOS COMPARATIVOS 5.1.
Integrales STO monocéntricas
Con objeto de comprobar si las fórmulas de integración implementadas en nuestro programa eran correctas, se han recalculado las energías atómicas obtenidas por Clementi y Raimondi [Ref. 47] y por Clementi y Roetti [Ref. 48] en cierto número de átomos (Tabla 2). La coincidencia de resultados que se observa en esta tabla garantiza la corrección de nuestros procedimientos, muchos de los cuales son comunes para estos cálculos analíticos y para los posteriores cálculos numéricos, objeto último de nuestro estudio. La discrepancia entre resultados que se observa en algunos casos con la base doble-Z es debida a que Clementi y Roetti realizaron cálculos ROHF, mientras que nosotros los realizamos UHF. Por ello, nuestros resultados tienen energías un poco más bajas que las de referencia, pero solo en los átomos con electrones desapareados.
5.2.
Integrales STO-nG monocéntricas
En la Tabla 3 hemos comparado los métodos de integración con desarrollos gaussianos realizados con nuestro programa, con los resultados analíticos de Clementi-Raimondi [Ref. 47] con bases SZ. Hemos ensayado varios desarrollos con distinto número de gaussianas utilizadas en el mismo. Como podemos observar los desarrollos STO-12G ya dan un resultado muy aceptable, y los STO-18G son casi exactos (pero requieren un tiempo de calculo mucho mayor, por ejemplo, el Mg con 12G tarda 1’6’’ y con 18G 5’20’’).
Capítulo II-59
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
En la Tabla 4 comparamos los métodos de integración analíticos con los resultados analíticos de Clementi-Roetti [Ref. 48] con bases DZ. Se alcanza un grado de precisión parecido al de la Tabla 3, aunque se utiliza un número doble de funciones de base. Las discrepancias en las energías marcadas con un asterisco se deben, como se ha indicado anteriormente, a que los cálculos de Clementi-Roetti eran del tipo ROHF y los nuestros UHF.
SIMPLE-Z ATOMOS
ESTADO
DOBLE-Z
CLEMENTI RAIMONDI
NUESTRO PROGRAMA
CLEMENTI ROETTI
NUESTRO PROGRAMA
He
1
-2.8476563
-2.8476563
-2.8616726
-2.8616726
Li
2
-7.4184820
-7.4184820
-7.4327213
-7.4327442
Be
1
-14.556740
-14.556740
-14.572369
-14.572369
B
2
-24.498369
-24.498369
-24.527920
-24.528149
C
3
P
-37.622389
-37.622389
-37.686749
-37.687967
N
4
S
-54.268900
-54.268900
-54.397951
-54.401061
O
3
-74.540363
-74.540365
-74.804323
-74.806192
F
2
P
-98.942113
-98.942113
-99.401309
-99.401983
Ne
1
S
-127.81218
-127.81218
-128.53511
-128.53511
Na
2
-161.12392
-161.12392
-161.84999
-161.84999
Mg
1
-198.85779
-198.85779
-199.60701
-199.60701
Al
2
-241.15376
-241.15377
-241.87307
-241.87316
Si
3
P
-288.08996
-288.08998
-288.85116
-288.85130
P
4
S
-339.90988
-339.90988
-340.71595
-340.71613
S
3
-396.62762
-396.62762
-397.50229
-397.50229
Cl
2
P
-458.52369
-458.52369
-459.47960
-459.47047
Ar
1
S
-525.76525
-525.76525
-526.81511
-526.81511
S S S P
P
S S P
P
Tabla 2: Valores de la Energía total de diversos átomos comparados con los de Clementi-Raimondi [Ref. 47] y Clementi-Roetti [Ref. 48].
Capítulo II-60
INTEGRACIÓN CON FUNCIONES DE BASE DE SLATER
Átomos
Estado
Método ClementiRaimondi
Analítico OMCM
STO-6G OMCM
STO-12G OMCM
STO-18G OMCM
He
1
-2.8476563
-2.8476563
-2.8463014
-2.8476492
-2.8476561
Li
2
-7.4184820
-7.4184820
-7.4153780
-7.4184659
-7.4184817
Be
1
-14.556740
-14.556740 -14.5510412 -14.5567102 -14.5567394
B
2
-24.498369
-24.498369
-24.4893219
-24.4983220 -24.4983679
C
3
P
-37.622389
-37.622389
-37.6091833
-37.6223213 -37.6223876
N
4
S
-54.268900
-54.268900
-54.2506972
-54.2688079 -54.2689000
O
3
-74.540363
-74.540365
-74.5162941
-74.5402446 -74.5403629
F
2
P
-98.942113
-98.942113
-98.9112827
-98.9419611 -98.9421110
Ne
1
S
-127.81218
-127.81218
-127.773672
-127.811993 -127.812178
Na
2
-161.12392
-161.12392
-161.076990
-161.123693 -161.123919
Mg
1
-198.85779
-198.85779 -198.801526 -198.857512 -198.857783
Al
2
-241.15376
-241.15377
Si
3
P
-288.08996
-288.08998 -288.012724 -288.089601 -288.089973
P
4
S
-339.90988
-339.90988
-339.820416
S
3
-396.62762
-396.62762
-396.524986 -396.627118 -396.627609
Cl
2
P
-458.52369
-458.52369
-458.406947
-458.523125 -458.523682
Ar
S1
-525.76525
-525.76525
-525.633445
-525.764617 -525.765244
S S S P
P
S S P
P
-241.087871
-241.153446 -241.153765
-339.909448 -339.909877
Tabla 3: Valores de la Energía total de diversos átomos con base SZ, comparados con los de Clementi-Raimondi [Ref. 47], valor analítico y desarrollos STO-nG (n = 6, 12 y 18). Las cifras en negrita son las que coinciden con el valor analítico.
Capítulo II-61
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
Método Átomos
Estado
ClementiRoetti
Analítico OMCM
STO-6G OMCM
STO-12G OMCM
STO-18G OMCM
He
1
-2.8616726
-2.8616726
-2.8609261
-2.8616686
-2.8616724
Li
2
-7.4327213*
-7.4327442
-7.4307146
-7.4327334
-7.4327440
Be
1
-14.572369
-14.572369
-14.5681106
B
2
-24.527920*
-24.528149
-24.5210913 -24.5281113 -24.528149
C
3
P
-37.686749*
-37.687967
-37.6774081
-37.6879113 -37.687966
N
4
S
-54.397951*
-54.401061
-54.3862489
-54.4009831 -54.401060
O
3
-74.804323*
-74.806192
-74.7862617
-74.8060887 -74.806191
F
2
P
-99.401309*
-99.401983
-99.3761483
-99.4018495 -99.401981
Ne
1
S
-128.53511
-128.53511
-128.50255
-128.53494
-128.53510
Na
2
-161.84999*
-161.84999
-161.80989
-161.84979
-161.84999
Mg
1
-199.60701
-199.60701
-199.55844
-199.60676
-199.60700
Al
2
-241.87307*
-241.87316
-241.81533
-241.87286
-241.87315
Si
3
P
-288.85116*
-288.85130
-288.78309
-288.85095
-288.85130
P
4
S
-340.71595*
-340.71613
-340.63656
-340.71572
-340.71612
S
3
-397.50229*
-397.50303
-397.41108
-397.50255
-397.50302
Cl
2
P
-459.47960*
-459.48004
-459.37471
-459.47949
-459.48003
Ar
1
S
-526.81511
-526.81511
-526.69530
-526.81448
-526.81510
S S S P
P
S S P
P
-14.5723459 -14.572369
Tabla 4 : Valores de la Energía total de diversos átomos comparados con los de Clementi-Roetti [Ref. 48], valor analítico y desarrollos STOnG (n = 6, 12 y 18). Los datos bibliográficos con * han sido calculado con ROHF en vez de UHF. Las cifras en negrita son las que coinciden con el valor analítico.
Capítulo II-62
INTEGRACIÓN CON FUNCIONES DE BASE DE SLATER
5.3.
Integrales STO-nG policéntricas
En la Tabla 5 hemos comparado nuestros resultados con desarrollos STO-nG con los realizados con el programa G94W. En dicha tabla se muestran la energía y el momento dipolar calculados en la geometría experimental para algunas moléculas. Se ha preferido realizar una comparación entre energías, y entre momentos dipolares, mejor que una comparación directa entre valores de integrales. Así, se puede comprobar de forma simultánea la exactitud de todas las integrales implicadas en cada cálculo, ya que resultaría extraordinariamente improbable obtener un buen valor de E ó µ con integrales imprecisas.
Método
Molécula
H2 Re=1.40112 ua
CO Re=2.13218 ua
6G Slater (OMCM)
12G Slater (OMCM)
18G Slater (OMCM)
STO-18G (Gaussian)
E (u.a.)
µ (u.a.)
E (u.a.)
µ (u.a.)
E (u.a.)
µ (u.a.)
E (u.a.)
µ (u.a.)
-1.1252920
0.00000
-1.1255776
0.00000
-1.1255791
0.00000
-1.1255789
0.00000
-112.3384313
0.03910
-112.3386164
0.03910
-112.3386109
0.03976
-75.7014120
0.69260
-75.7015322
0.69260
-75.701526
0.69259
-113.4754432
0.61956
-113.4756283
0.61956
-113.4756562
0.61961
-112.3014485 0.03961
H2O RHO=1.80887 ua αHOH=104.522º
-75.6772411
0.69282
H2CO RCH=2.08248 ua RCO=2.27523 ua αHCO=121º41’
-113.4384329 0.61921
Tabla 5: Cálculo de la Energía (E) y el momento dipolar (µ) a la distancia y ángulos experimentales para varias moléculas.
Capítulo II-63
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
Capítulo II-64
Capítulo III: INTEGRACIÓN CON FUNCIONES DE BASE GAUSSIANAS
INTEGRACIÓN CON FUNCIONES DE BASE GAUSSIANAS
Capítulo III: INTEGRACIÓN CON FUNCIONES DE BASE GAUSSIANAS 1. INTRODUCCIÓN 2. TEOREMA DEL PRODUCTO (DE GAUSSIANAS) 3. INTEGRALES DE SOLAPAMIENTO 4. INTEGRALES DE ENERGÍA CINÉTICA 5. INTEGRALES DE MOMENTO 6. INTEGRALES DE ENERGÍA POTENCIAL 7. INTEGRALES BIELECTRÓNICAS 8. ARMÓNICOS ESFÉRICOS Y ARMÓNICOS SÓLIDOS 9. RESULTADOS
1. INTRODUCCIÓN En el capítulo precedente hemos visto cómo se realizan los desarrollos STO-nG. En este capítulo expondremos los métodos de cálculo que se han implementado en nuestro programa para llevar a cabo la integración con funciones gaussianas. Las funciones gaussianas GPµ , que utilizaremos en nuestros cálculos, son gaussianas cartesianas sin normalizar y tienen la forma: G Pµ = x nA x y
ny A
z nA z e
− α r A2
(1)
n
donde la parte preexponencial xnAx y Ay znAz se considera centrada en el mismo punto A que la parte exponencial e − α r . 2 A
Capítulo III-65
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
Las integrales con funciones de base gaussianas contienen como máximo dos funciones del tipo (1) por cada electrón, debido al hecho de que al multiplicar dos gaussianas con distintos centros se obtiene una nueva gaussiana centrada en un punto intermedio. Ello permite reducir el número de centros de las integrales y constituye el motivo más importante de la generalización del uso de funciones de base gaussianas.
2. TEOREMA DEL PRODUCTO (DE GAUSSIANAS) El teorema del producto establece que el producto de dos gaussianas con distinta parte preexponencial y distinto centro es otra gaussiana. r r Sea A el vector de posición del centro de la gaussiana y sea r el vector de r posición del electrón respecto a unos ejes principales. Por tanto, r A viene dado por (Ver Figura 1):
Figura 1: Vectores implicados en la definición de funciones de base
x A = x − A x r r r rA = r − A ⇒ y A = y − A y z A = z − A z
(2)
El producto de dos gaussianas centradas en los puntos A y B viene dado por:
Capítulo III-66
INTEGRACIÓN CON FUNCIONES DE BASE GAUSSIANAS
G A G B = x nAxA y AyA z nAzA x nBxB y ByB z nBzB exp(−αrA2 ) exp(−β rB2 ) n
n
(3)
Nótese que las gaussianas que utilizamos en nuestros desarrollos no están normalizadas. De la Figura 2, podemos ver que: r r r rA = r − A r r r rB = r − B r r r rP = r − P r r r αA + β B P= α+β
(4)
Desarrollando el exponente de las gaussianas obtenemos: r r r r − αrA2 − β rB2 = −α(r − A)2 − β(r − B)2 = r r 2 r r 2 r r r r = −α (r − P) + (P − A) − β (r − P) + (P − B) = → → r r = −α(rP + PA)2 − β(rP + PB)2 =
[
] [
__ 2
]
(5)
__ 2
r → r → = −α (rP2 + PA + 2rP ⋅ PA) − β (rP2 + PB + 2rP ⋅ PB) = = −(α +
β) rP2
__ 2
__ 2
r → r → − α PA − β PB − 2(α rP ⋅ PA + β rP ⋅ PB)
Como: r r → r r αA + β B r β PA = P − A = −A = AB α+β α+β r r → r r αA + β B r α → −B = − PB = P − B = AB α+β α+β
(6)
→ → r r αβ → αβ → rP (α PA + β PB) = rP AB− AB = 0 α+β α +β
(7)
→
por tanto:
Capítulo III-67
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
Figura 2: Vectores implicados en el teorema del producto de gaussianas.
Llamando γ = α + β , nos queda:
β 2 2 2 − αrA − βrB = − γ rP − α α + β
2 __ 2 AB − β
2 2 2 αβ __ __ 2 AB (8) AB = − γ rP − α+β α + β α
Para el producto de gaussianas tenemos que:
G A G B = x nAxA y AyA z nAzA x nBxB y ByB z nBzB K exp(− γrP2 ) n
n
(9)
donde: γ = α+β αβ ___ 2 K = exp − AB γ
(10)
Luego el producto de dos gaussianas sin parte preexponencial es otra gaussiana (también sin parte preexponencial), multiplicada por el factor constante K. La nueva gaussiana tiene exponente
γ = α + β y está centrada
en un punto del segmento AB, dado por la relación: →
→
α OA + β OB OP = α+β →
Capítulo III-68
(11)
INTEGRACIÓN CON FUNCIONES DE BASE GAUSSIANAS
Veamos ahora que ocurre con el producto de los factores preexponenciales de las gaussianas GA y GB. Expresando todas las variables
r
respecto al centro P del producto de las partes exponenciales: ___
x A = x − A x = (x − Px ) + (Px − A x ) = xP + PA x
(12)
y utilizando el desarrollo del binomio de Newton:
x
___
= ( x p + PA x )
n xA A
n xA
n xA ___ n xA − i
= ∑ PA x i =0
n xA i x p i
(13)
por tanto el producto de potencias con distintos centros viene dado por: n xA n xB ___ n xA − i ___ n xB − j
n xA n xB i + j x p = i = 0 j= 0 i j n xA + n xB n xA n xB ___ n xA − i ___ n xB − j n xA n xB k x p = = ∑ ∑∑ PA x PB x i j k =0 i = 0 j= 0 k =i + j x
=
n xA A
x
n xB B
= ∑∑ PA x
n xA + n xB
∑ k =0
___
PB x
(14)
___
f k (n xA , n xB , PA x , PB x ) x kp
Definiendo: ___
n xA n xB ___ n xA − i ___ n xB − j
___
f k (n xA , n xB , PA x , PB x ) ≡ ∑∑ PA x
PB x
i = 0 j= 0 k =i + j
n xA n xB i j
(15)
Por el teorema del producto:
G AG B = K
n xA + n xB n yA + n yB n zA + n zB
∑ ∑ ∑f
k x =0
k y =0
k z =0
f f x kp x y p y z kp z exp(− γrp2 ) = G (px )G (py )G (pz ) (16) k
kx ky kz
donde:
G (ph ) = K h
n hA + n hB
∑f
k h =0
kh
h kp h exp(− γh 2p )
(h = x , y, z)
(17)
Capítulo III-69
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
αβ ___ 2 el centro P está definido en (4), y Kh = exp − ABh . Atención: no confundir γ el índice (minúsculas) kh con la constante (mayúsculas) Kh. Nuestro programa calcula las constantes fk, que aparecen en (14) mediante la función:
FKH(L1, L2, L, PAH, PBH, PH, K)
(18)
que admite los valores L1 = nxA, nyA ó nzA; L2 = nxB, nyB ó nzB; PAH = PAx, PAy ó PAz; PBH = PBx, PBy ó PBz; K = kx, ky ó kz. Las variables L y PH se utilizan cuando hay que calcular integrales de momento. Como hemos podido observar, el producto de gaussianas sin parte preexponencial es una nueva gaussiana centrada en P. En cambio, el producto de gaussianas con parte preexponencial es una combinación lineal de gaussianas centradas en P.
3. INTEGRALES DE SOLAPAMIENTO La integral de solapamiento entre dos gaussianas GA y GB centradas en r r los puntos A y B respectivamente, se reduce al producto de tres integrales monodimensionales: +∞
SAB =
∫
r GAGBd r =
−∞
+∞
∫
+∞
+∞
Gp(x )dx Gp(y )dy Gp(z )dz
−∞
∫
∫
−∞
(19)
−∞
cada una de las cuales se pueden relacionar con la integral: +∞
Fk (γ) =
∫ exp(−γt ) t dt 2
k
(20)
−∞
que para valores de k impares es nula, y para valores de k pares vale:
γ
−(k +1) / 2 (k − 1)!! π k Fk (γ) = 2 2 0 donde (k-1)!! = (k-1)⋅(k-3)⋅(k-5)...3⋅1.
Capítulo III-70
k par k impar
(21)
INTEGRACIÓN CON FUNCIONES DE BASE GAUSSIANAS
El programa OMCM calcula el valor de Fk (γ) , mediante una función llamada (22)
FKGAMMA(GAMMA, K)
en la que GAMMA=γ y K=k. Esta función utiliza a su vez otra llamada FACT1(NUM) que devuelve el valor del semifactorial (NUM)!!. Aplicando la (20) a las integrales que aparecen en (19) obtenemos: +∞
∫G
(h ) p
d h = Kh
−∞
= Kh
n hA + n hB 2
∑f
k h =0
2k h
+∞n hA + n hB
∫ ∑f
−∞ k h = 0
kh
h exp(− γh )dh = K h kh p
(2k h − 1)!! π 2kh
2 p
γ
+∞
n hA + n hB
∑
k h =0
fk h
∫h
kh p
exp(− γh 2p ) dh =
−∞
(23)
−( 2 k h +1) / 2
Nótese que se ha hecho el cambio fk h (kh = 0, 1, 2,...) → f2k h (2kh = 0, 2, 4,...) , que tiene en cuenta que la función fk h es nula para valores impares de kh. La integral de solapamiento se puede obtener, por tanto, mediante una fórmula cerrada:
S AB = K x K y K z
n yA + n yB
n xA + n xB 2
n zA + n zB 2
∑ξ ∑ξ ∑ξ ⋅ γ
k x =0
2
kx
k y =0
ky
k z =0
−( 2 k x + 2 k y + 2 k z +3) / 2
kz
π
3
2
(24)
en la que: ξk h = f2k h y
(2kh − 1)!! 2k h
con h = x, y, z
(25)
f 2 k h está dada por (15).
El programa OMCM contiene una función para calcular la integral (19), llamada:
INT_GPH (IL1, IL2, IL, EXP1, EXP2, AH, BH)
(26)
Los valores de las variables que toman esta función son: IL1 = nxA, nyA ó nzA; IL2 = nxB, nyB ó nzB; EXP1=α; EXP2=β; AH = Ax, Ay ó Az; BH = Bx, By ó Bz. La integral de solapamiento SAB la calcula la función SG(NB1,NB2,NG1,NG2) que toma como variable las funciones de base que vamos a integrar y el número de gaussianas del desarrollo. Las variables IL y PH se utilizan cuando hay que calcular integrales de momento.
Capítulo III-71
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
4. INTEGRALES DE ENERGÍA CINÉTICA Las integrales de energía cinética con funciones de base gaussianas son del tipo: TAB = −
∂2 r ∂2 1 ∂2 GA 2 + 2 + 2 GBd r = ∂x 2 ∂z ∂y
∫
1 ∂2GB ∂2GB ∂2GB 1 dz = − [T1 + T2 + T3 ] dy + GA dx + GA = − GA 2 2 2 2 2 ∂x ∂y ∂z
∫
∫
(27)
∫
Cada una de las integrales T1, T2 y T3 se descompone en un producto de tres integrales: T1 =
∫
GxAGyAGzA
∂2 ∂x
2
(G G G )drr = ∫ G x B
y B
z B
x A
∂2GBx ∂x
2
∫
∫
dx GyAGBydy GzAGBz dz
(28)
Conviene definir: Th = < GhA | Sh = <
GhA
∂2 ∂h2
| GhB
| GhB > = >=
∫
∫
GhA
GhAGhBdh
∂2GhB ∂h2
T1 = NANB TxSySz dh ⇒ T2 = NANBSx TySz T3 = NANBSxSy Tz
(29)
con lo que la integral a calcular tiene la forma: TAB = −
[
1 TxSySz + Sx TySz + SxSy Tz 2
]
(30)
Las integrales del tipo Sh son integrales de solapamiento, ya estudiadas, y las Th se pueden calcular por partes:
Th = < GhA |
∂
2 2
∂h +∞
| GhB > =
∫G
h A
∂
2
GhB 2
∂h
∂GhA h dh u = GA → du = ∂h dh = = 2 h h ∂ ∂ G G B B ∂h2 dh = dv → v = ∂h
(31)
+∞ ∂GhA ∂GhB ∂GhB dh = GhA − ∂h ∂h ∂h −∞ −∞
∫
Para resolverla veamos antes cual es el valor de
G hB = h nBhB exp(−β h 2B ) :
Capítulo III-72
∂GhB , sabiendo que ∂h
INTEGRACIÓN CON FUNCIONES DE BASE GAUSSIANAS
n −1 n +1 2 2 ∂G hB ∂G hB ∂h B ∂G hB n hB h BhB exp(−β h B ) − 2βh BhB exp(−β h B ) = = = ∂h ∂h B ∂h ∂h B − 2β h nBhB +1 exp(−β h 2B )
n hB ≥ 1 n hB = 0
(32)
que como puede verse es una función impar y por tanto: +∞
h ∂GhB GA =0 ∂h −∞
(33)
La integral Th queda reducida a: +∞
∂G hA ∂G hB dh = ∂h ∂h −∞
Th = − ∫ +∞
(
)
(
)
= − ∫ n hA h nAhA −1 − 2αh nAhA +1 exp(−αh 2A ) n hB h nBhB −1 − 2β h nBhB +1 exp(−β h 2B ) dh = −∞
= − n hA n hB ∫ h nAhA −1h nBhB −1 exp(−αh 2A ) exp(−β h 2B ) dh +
(34)
+ 2β n hA ∫ h nAhA −1h nBhB +1 exp(−αh 2A ) exp(−β h 2B ) dh + + 2αn hB ∫ h nAhA +1 h nBhB −1 exp(−αh 2A ) exp(−β h 2B ) dh − − 4αβ ∫ h nAhA +1h nBhB +1 exp(−αh 2A ) exp(−βh 2B ) dh Las cuatro integrales obtenidas se calculan en nuestro programa con la función INT_GPH(IL1,IL2,IL,EXP1,EXP2,AH,BH), y el valor total de la integral de energía cinética se calcula en la función TG(NB1,NB2,NG1,NG2).
5. INTEGRALES DE MOMENTO Para el cálculo de algunas propiedades moleculares como el momento dipolar, momento cuadrupolar, etc., es necesario evaluar también integrales del tipo:
M AB = < G A | x n xC y
n yC
z n zC | G B > = ∫ G A x n xC y
n yC
r z n zC G B d r
(35)
En las cuales aparecen desarrollos del tipo:
Capítulo III-73
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
n hA n hB n hC ___ n hA − i ___ n hB − j
h nAhA h nBhB h n hC = ∑ ∑ ∑ PA h
PB h
i = 0 j= 0 ii = 0
=
n hA + n hB + n hC
∑
k h =0
=
n hA + n hB + n hC
∑
k h =0
n n n Phn hC − ii hA hB hC h ip+ j+ ii = i j ii
n n n n −i n −j hA hB hC ___ hA ___ hB n hC − ii n hA n hB n hC i + j+ ii k h h p x p = (36) PB h Ph ∑ ∑ ∑ PA h i j ii i = 0 j= 0 ii = 0 k h = i + j+ ii ___
___
f k h ( n hA , n hB , n hC , PA h , PB h , Ph ) h kp h
donde hemos definido: ___
___
f k h (n hA , n hB , n hC , PA h , PB h , Ph ) = n hA n hB n hC ___ n hA −i ___ n hB − j
= ∑∑∑ PA h
PB h
n hC −ii h
P
i =0 j=0 ii = 0
(37)
n hA n hB n hC i+ j+ii h p i j ii
Definiendo las funciones:
G hp = K h
n hA + n hB + n hC
∑f
k h =0
kh
h kp h exp(− γh 2p )
(38)
podemos calcular las integrales: +∞
h ∫ Gpd h = Kh
−∞
= Kh
+∞n hA + n hB + n hC
∫
−∞
n hA + n hB + n hC 2
∑f
k h =0
2k h
∑ f k h h kp h exp(− γh 2p )dh = K h
k h =0
γ −( 2 k h +1) / 2
+∞
n hA + n hB + n hC
∑
k h =0
fk h
∫h
kh p
exp(− γh 2p ) dh =
−∞
(39)
(2k h − 1)!! π 2kh
Utilizando la definición (38), la integral de momento nos quedaría como: n yC
MAB = < GA | xn y zn | GB > = xC
zC
+∞
+∞
+∞
−∞
−∞
−∞
x y z ∫ Gpdx ∫ Gpdy ∫ Gpdz
(40)
y usando (39):
M AB = K x K y K z donde:
Capítulo III-74
n xA + n xB + n xC n yA + n yB + n yC n zA + n zB + n zC 2 2 2
∑ξ
k x =0
kx
∑ξ
k y =0
ky
∑ξ ⋅ γ
k z =0
kz
−( 2 k x + 2 k y + 2 k z + 3) / 2
π
3
2
(41)
INTEGRACIÓN CON FUNCIONES DE BASE GAUSSIANAS
ξk h = f2k h
(2kh − 1)!! 2k h
con h = x, y, z
(42)
Nótese que ahora las funciones Gp(h) llevan una suma en nhA+nhB+nhC. La función:
FKH (L1, L2, L, PAH, PBH, PH, K)
(43)
INT_GPH (IL1, IL2, IL, EXP1, EXP2, AH, BH)
(44)
y la función:
aceptan los valores de nxC, nyC y nzC en las variables L e IL y el valor de Px, Py o Pz en la variable PH. En caso de no ser integrales de momento, estas variables deben ser cero. El valor de la integral de momento se obtiene de la función MG(NB1,NB2,NG1,NG2,MOMENTO).
6. INTEGRALES DE ENERGÍA POTENCIAL 6.1.
Transformada de Laplace
Las integrales monoelectrónicas de energía potencial y las bielectrónicas dependen de términos que contienen 1/r12. La transformada de Laplace permite sustituir esta función por una integral de manejo menos complicado que el del término 1/r12. La transformada (t) de la función g(ρ) es: ∞
(t)= ∫ exp(−tρ)g(ρ)dρ
(45)
0
En nuestro caso interesa usar: ∞
(r 2 ) = ∫ exp(−r 2ρ)ρλ / 2 −1dρ
(46)
0
Haciendo el cambio u=ρr 2 obtenemos: 1 (r 2 ) = r
λ
∫
∞
0
( ) ( ) λ
exp(−u)u λ / 2 −1du = 1 Γ λ2 r
(47)
Despejando:
Capítulo III-75
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
(1 / r ) = (r ) (λ / 2) = λ
Cuando
2
1 Γ λ2
(
)∫
∞
0
exp(−r 2ρ)ρλ / 2 −1dρ
(48)
λ =1 y por tanto: 1 1 = r π
∫
∞
0
exp(−r 2ρ)ρ−1 / 2dρ
(49)
y eligiendo ρ=u2 y r=r12: 2 2 1 2 ∞ r12 = π ∫0 exp(−r12u )du
6.2.
(50)
Integrales de energía potencial
Las integrales de energía potencial pueden ser mono, bi, ó tricéntricas. Estas últimas aparecen debido a que el potencial puede estar referido a un punto distinto de los centros de las funciones de base implicadas en la integral.
Figura 3: Vectores implicados en el cálculo de las integrales de energía potencial
Dado que el potencial de los cálculos moleculares es:
Capítulo III-76
INTEGRACIÓN CON FUNCIONES DE BASE GAUSSIANAS
−ZA
r V (r ) =
∑ rr − Rr
(51) A
su valor medio con funciones de base gaussianas es de la forma:
ˆ | GB >= − < GA | V
todos los núcleos
∑
Zc < G A |
c
1 | G B > = −∑ Z c Vc rc c
(52)
donde rc2 = (x − Cx )2 + (y − Cy )2 + (z − Cz )2 , y Zc es la carga nuclear del núcleo C.
Introduciendo
Vc = < G A |
Vc = =
1 π
1 π ∞
transformada
de
Laplace
1 , rc
de
la
integral
1 | G B > nos quedaría: rc
∞
∫ dρ ρ
−1 2
0
−1 2
+∞
+∞
∫ dxdydz exp(−ρr ) G G 2 c
A
B
=
−∞
∫ dρ ρ ∫ dx ⋅ x 0
la
n xA A
x
n xB B
−∞
+∞
+∞
−∞
−∞
(53)
exp(−ρx ) exp(−αx ) exp(−βx ) ∫ dy ⋅ ... ∫ dz ⋅ ... 2 c
2 A
2 B
Llamando Uh(ρ) a las integrales respecto a h (= x, y, z) tenemos: Vc =
1
∞
∫ dρ ρ π
−1 2
Ux (ρ) Uy (ρ) Uz (ρ)
(54)
0
con: +∞
U h (ρ) = ∫ dh h nAhA h nBhB exp(−ρh c2 − αh 2A − βh 2B )
(55)
−∞
Haciendo los cambios: γ = α+β
θ = γ+ρ
αβ Kh = exp − ABh γ αAh + βBh Ph = γ
γρ ___ 2 PC h Λh = exp − θ γPh + ρCh Qh = θ
hP = h − Ph
QP = h − Qh
___ 2
(56)
Capítulo III-77
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
+∞
Uh(ρ) = ΛhKh
∫ dh h
nh n'h A hB
exp(−θh2Q ) = ΛhJh
(57)
−∞
puede expresarse la integral Vc en la forma: Vc =
∞
1
dρΛ(ρ) ρ π ∫
−1
2
Jx (ρ)Jy (ρ)Jz (ρ)
(58)
0
Λ (ρ) = Λ x (ρ) ⋅ Λ y (ρ) ⋅ Λ z (ρ) . La integral obtenida (58) se puede
en la que
calcular haciendo el cambio de variable: γ t2
ρ=
2
⇒ t2 =
1−t γ dt = dρ 3 1 2θ 2ρ 2
ρ θ
(59)
con el que se llega: 3
1
γ Vc = dtΛ(t ) 2 γ π 0 1 − t 2
∫
2
2
Jx (ρ)Jy (ρ)Jz (ρ)
(60)
Para calcular cada integral Jh puede operarse como en el caso de las integrales de solapamiento, cuando desarrollamos el producto
h
n hA A
h
n hB B
=
n hA + n hB
∑
k h =0
___
h nAhA h nBhB :
___
f k h (n hA , n hB , QA h , QB h ) h Qk h
___
n hA n hB ___ n hA −i ___ n hB − j
___
f k h (n hA , n hB , QA h , QB h ) = ∑∑ QA h i = 0 j= 0
i + j= k h
QB h
n hA n hB i j
(61)
Sustituyendo en la integral Jh:
Jh = Kh
n hA + n hB
∑ f kh
k h =0
+∞
∫ dh
Q
h Qk h exp(−θh Q2 )
(62)
−∞
Utilizando (21) para valores de kh pares:
Jh = Kh
Capítulo III-78
n hA + n hB 2
∑
k h =0
f 2kh θ
−
2 k h +1 2
(2k h − 1)!! π 2kh
(63)
INTEGRACIÓN CON FUNCIONES DE BASE GAUSSIANAS
De (59) podemos deducir que: ___
___
___
___
QAh = PA h − t2 PC h
___
___
QBh = PB h − t2 PC h
(64)
Por lo tanto: 1
γ 2 Jh = Kh 2 1− t n hA n hB
___
n hA + n hB 2
∑
k h =0
___
⋅ ∑∑ (PA h − t PC h ) 2
γ 2 1− t
n hA −i
−k h
___
(2k h − 1)!! π⋅ 2kh ___
(PB h − t PC h ) 2
i = 0 j= 0 i + j= k h
n hB − j
(65)
n hA n hB i j
y sustituyendo en Vc:
2 (π ) 2 Vc = K xK yKz γ π 3
⋅ n xA n xB
n hA + n hB n hA + n hB n hA + n hB 1 2 2 2
∑ ∑ ∑
k h =0
k h =0
k h =0
1− t2 dt Λ ( t ) ∫0 γ 2
k x + k y +k z
⋅
(2k x − 1)!! (2k y − 1)!! (2k z − 1)!! ⋅ k 2kx 2kz 2 y
___
___
⋅ ∑∑ (PA x − t PC x ) 2
n xA −i
___
___
(PB x − t PC x ) 2
i = 0 j= 0 i + j= k x
(66) n xB − j
n xA n xB n yA n yB n zA n zB ∑∑ ...∑∑ ... i j i=0 j=0 i=0 j=0
Como puede verse, lo que se obtiene contiene productos de polinomios en potencias pares de t, que luego habrá que integrar respecto de t. Las integrales que restan son todas del tipo: ___ 2 Fm γ PC =
___ 2 2m t exp PC t2 dt − γ ∫0 1
(67)
En nuestro programa hemos incluido varias subrutinas análogas a las utilizadas para las integrales de solapamiento pero que, en vez de devolver un valor numérico, devuelven los coeficientes del polinomio en potencias pares de t. La subrutina que devuelve los coeficientes del polinomio es:
FKH_POL (L1, L2, EXP1, EXP2, PAH, PBH, PCH, K, PH)
(68)
Las demás sirven para sumar, multiplicar y elevar a una potencia los polinomios. Otra subrutina llamada FREQ16(Y,FV) y obtenida de [Ref. 50] calcula la función Fm. El cálculo total de la integral de energía potencial lo realiza la función VG(NB1,NB2,NG1,NG2) que llama a todas las anteriores.
Capítulo III-79
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
7. INTEGRALES BIELECTRÓNICAS El cálculo de las integrales bielectrónicas es parecido al de las integrales de energía potencial, con la salvedad de que ahora pueden ser tetracéntricas. Su resolución detallada puede encontrarse en el capítulo 5 de [Ref. 51]. En nuestro programa hemos recurrido a una subrutina cedida por Rafael López y J. Fernández Rico [Ref. 52], ligeramente modificada para adaptarla a nuestros convenios y notaciones.
8. ARMÓNICOS ESFÉRICOS Y ARMÓNICOS SÓLIDOS Cuando se trata de utilizar gaussianas para el cálculo de las integrales moleculares surge el siguiente problema: las funciones de base de tipo Slater suelen relacionarse con los armónicos esféricos: m
χ = N rn e−αr Yl Los
armónicos
(θ, ϕ) = N rn e−αrPlm (cos θ) e±imϕ
Ylm (θ, ϕ)
esféricos
pueden
sustituirse
(y
(69) se
sustituyen
habitualmente) por sus combinaciones lineales reales, dando lugar a lo que llamaremos “Armónicos Reales”. En estos, como habíamos considerado en el tema de integración analítica, el signo de m dependía de elegir un seno o un ~ coseno en la función empleada. Llamando m a este nuevo índice:
~ ≥ 0 ⇒ χ incluye cos(m ~ ϕ) χ = Nr n e − αr P m~ (cos θ) ⋅ cos(m ~ ϕ) m l ⇒ ~ n − αr m ~ < 0 ⇒ χ incluye sen (m ~ ϕ) ~ ϕ) m χ = Nr e Pl (cos θ) ⋅ sen (m
(70)
En el caso de funciones de base gaussianas, los armónicos esféricos no siempre aparecen explícitamente, sino que a menudo se sustituyen por productos de coordenadas x, y ó z elevados a potencias enteras:
χ = N ∑ C µ G µ = N ∑ Cµ x n x y y z n z e n
µ
− αr 2
(71)
µ
Estos productos de potencias de las coordenadas cartesianas se denominan “Armónicos Sólidos”, y están relacionados con los armónicos esféricos reales. Para encontrar la relación basta usar la definición de las coordenadas esféricas y propiedades trigonométricas. Se obtienen las relaciones que aparecen en la tabla III del apéndice TABLAS.
Capítulo III-80
INTEGRACIÓN CON FUNCIONES DE BASE GAUSSIANAS
La matriz de transformación que permite pasar de los armónicos esféricos a los armónicos sólidos no es cuadrada, pues para l ≥ 2 hay más armónicos sólidos que armónicos esféricos: S Px Py Pz dx2 − y2 d dxy dxz yz d3z2 −r2 f 3 2 x −3xy fy3 −3yx2 f 2 2 (x − y )z fxyz f 2 2 fxr −5xz yr2 −5yz2 f 3 2 5z −3r z
=
1
1 0 0 0 1 0 0 0 1 1 0
−1 0 0 0 0 0 0 1 0 0
0
0
0 0 1 0
0
0
0 0 0 1
−1 −1 2 0 0 0 0
0 −3
0
0
0
0
0 −3 0 0
0 1
1 0
0 0
0 0 −1 0
0 0
0 0
0 1
0 0
0 0
0 0
0 1
0 0
0 0
1 0
0 1 −3 0
0 0
0 0 −3 2
1
0
0 0 0 −4 0 0
0 0 −4 0
1 x y z x2 y2 z2 xy xz yz x3 2 0 x y 0 x 2z 0 y3 1 xy2 0 y 2z 0 z3 0 xz2 yz2 xyz
(72)
Como puede verse en la matriz correspondiente a los armónicos 0 ≤ l ≤ 3 (de
16x20), algunos armónicos esféricos están relacionados con sólo un armónico sólido, pero otros están relacionados con una combinación lineal de ellos. En el primer caso, si ya estaban normalizados, siguen estando normalizados; en el segundo caso, hay que volver a normalizarlos. Veamos, como ejemplo, la normalización de la función fx 3 −3xy 2 que por abreviar llamaremos f:
f ≡ f x3 −3xy 2 ⇒ | f = | x 3 − 3 | xy 2 | x 3 = x 3 R (r ) | xy 2 = xy 2 R (r ) x3 x3 x 3 xy 2 f f = N ∑ ∑ Cλ Cσ λ σ ⇒ λ σ = xy 2 x 3 xy 2 xy 2 λ σ r x 3 x 3 = ∫ x 6 R 2 d r = ∫ r 6 sen 6 θ cos 6 ϕr 2 sen θR 2 drdθdϕ =
1 a = a 1 (73)
= ∫ sen 7 θdθ∫ cos 6 ϕdϕ∫ r 8 R 2 dr = 1 a=
x 3 xy 2 x3 x3
∫ sen θdθ∫ sen ϕ cos ϕdϕ∫ r R dr = 1 = 5 ∫ sen θdθ∫ cos ϕdϕ∫ r R dr 7
2
7
4
6
8
8
2
2
f f = x 3 x 3 − 3 x 3 xy 2 − 3 xy 2 x 3 + 9 xy 2 xy 2 =
44 5
Capítulo III-81
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
de manera que la constante de normalización de la función
N=
f x 3 −3 xy 2 sería:
5 44
(74)
Análogamente se puede calcular la constante de normalización con las demás funciones, obteniéndose la matriz de transformación normalizada: 1
1 0 0 0 1 0 0 0 1 3
2 0 0 0 −1 2
− 3 2 0 0 0 −1 2
0 0 0 0 1
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0 5
44
0
0
0
−3 5
0
44
0 5
44
−3 5
0
0
0
0
0
0
0
0
0
− 15
0
0
0
44
0
0
15
0
0
0 15 244
0
0
0
0
0
0
0
0
0 15 244
0
0
0
0
0
0 0
15
244
0
28
15
0 −3 5
92
244
0
0
−3 5
28
92
2 5
0 − 4 15
92
0 244
0
0
− 4 15
0
0
244
(75) 0 0 0 1 0 0 0
En el programa OMCM esta matriz se encuentra almacenada en la subrutina:
SOL_REAL(L, S_R)
(76)
que toma el valor de L y devuelve una variable dimensionada que depende de L, llamada S_R.
9. RESULTADOS En la Tabla 1 se presenta una comparación de los resultados obtenidos con nuestro programa y con Gaussian 94W, para una serie de moléculas diatómicas representativas. Los cálculos con ambos programas han sido realizados con la misma base, mínima con exponentes dados por las reglas de Slater para cada átomo y desarrollos con 18 gaussianas [Ref. 53]. La base se eligió así porque permitía una primera comparación de nuestros resultados con los de [Ref. 49], que puede verse en la Tabla 2. En ella podemos comprobar que los cálculos STO-nG con los desarrollos de Fernández Rico et al [Ref. 53] pueden
Capítulo III-82
INTEGRACIÓN CON FUNCIONES DE BASE GAUSSIANAS
considerarse equivalentes a los realizados directamente con STO, a partir de unas 12 gaussianas. M o l é c u l a
Método STO-18G Slater (Gaussian)
Intervalo ajuste (u.a.)
STO-18G Slater (OMCM)
Re (u.a.)
µe (u.a.)
Ee (u.a.)
Eajuste (u.a.)
Re (u.a.)
µe (u.a.)
Ee (u.a.)
Eajuste (u.a.)
Li2 Be2
[4.9,5.5]
5.34282
0.00000
-14.841496
-14.841489
5.34283
0.00000
-14.841496
-14.841489
[3.1,3.8]
3.41411
0.00000
-28.974107
-28.974107
*
*
*
*
C2
[2.1,2.6]
2.43504
0.00000
-75.214884
-75.214898
2.43504
0.00000
-75.214885
-75.214898
N2
[1.9,2.5]
2.14781
0.00000
-108.58006
-108.58009
2.14770
0.00000
-108.58006
-108.58010
F2
[2.2,2.8]
2.44983
0.00000
-197.87877
-197.87879
2.44983
0.00000
-197.87877
-197.87879
LiH
[2.6,3.2]
3.12532
2.55015
-7.967088
-7.9670901
3.12532
2.54976
-7.967090
-7.967090
BH
[2.0,2.6]
2.59699
0.27216
-25.069872
-25.069873
2.59699
0.27218
-25.069873
-25.069873
NH
[1.8,2.4]
2.20101
0.23144
-54.671072
-54.671072
2.20101
0.23144
-54.671072
-54.671072
HF
[1.5,2.1]
1.94461
0.23709
-99.491174
-99.491187
1.94461
0.23709
-99.491174
-99.491187
CO
[1.9,2.5]
2.18195
0.26904
-112.34557
-112.34558
2.18194
0.26757
-112.34557
-112.34558
BF
[2.1,2.7]
2.34330
0.87323
-123.61606
-123.61608
2.34330
0.87298
-123.61606
-123.61608
LiF
[2.4,3.0]
2.62803
1.05767
-106.37176
-106.37177
2.62807
1.05854
-106.37176
-106.37177
Tabla 1: Cálculo de la distancia de equilibrio (Re), el momento dipolar (µe) y de la
Energía en el equilibrio tanto la ajustada (Eajustada) como la calculada a la distancia de equilibrio (Ee) para varias moléculas. Los exponentes de los desarrollos gaussianos han sido tomados de [Ref. 49]. * El caso del Be2 es especial, y se explica en el texto.
El uso de una base STO-18G con exponentes atómicos planteó algunas dificultades a la hora de realizar comparaciones con los resultados del programa Gaussian, debido a que esta base no se encuentra entre las predefinidas en G94W. Hubo que emplear la opción GEN del programa e introducir en él la base STO-18G de [Ref. 53] necesaria para cada molécula. Para ello se utilizó el programa GENGAUSSIAN, elaborado por D. Jesus Sánchez Márquez en nuestro laboratorio, que permite introducir este tipo de bases (y muchas otras) como datos para Gaussian de una forma cómoda.
Capítulo III-83
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
Re (u.a.)
STO-1G
STO-3G
STO-6G
STO-12G
STO-18G
[Ref 49]
Li2 Be2
5.051
-12.327291
-14.663610
-14.834567
-14.840718
-14.840749
-14.84075
3.780
-24.445178
-28.734033
-29.045082
-29.056394
-29.056453
-29.05645
C2
2.3475
-63.998309
-74.452463
-75.185052
-75.211083
-75.211215
-75.21122
N2
2.068
-92.674988
-107.516241
-108.537605
-108.573439
-108.573620
-108.57362
F2
2.680
-168.661103
-195.983697
-197.795567
-197.856560
-197.856858
-197.85686
LiH
3.015
-6.640727
-7.872781
-7.963393
-7.966648
-7.966665
-7.96666
BH
2.329
-21.263241
-24.799243
-25.052923
-25.062048
-25.062094
-25.06210
NH
1.976
-46.596103
-54.118327
-54.642853
-54.661094
-54.661185
-54.66119*
HF
1.733
-84.840993
-98.531031
-99.447651
-99.478386
-99.478536
-99.47854
CO
2.132
-95.929725
-111.250296
-112.306575
-112.343386
-112.343570
-112.34357
BF
2.385
-105.561518
-122.425241
-123.575905
-123.615303
-123.615498
-123.61550
LiF
2.850
-90.612189
-105.338452
-106.331503
-106.365042
-106.365207
-106.36521
Molécula
Tabla 2: Energías moleculares obtenidas con desarrollos STO-nG de [Ref. 53] y un número creciente de gaussianas por STO. Las cifras en negrita son aquellas que coinciden con el resultado de Ransil. El dato con * está calculado con el programa gaussian 94W ya que el de la bibliografía es erróneo.
La comparación presentada en la Tabla 1 se ha realizado ajustando los datos (Energía / distancia internuclear) mediante nuestro programa UCA-VIB [Ref. 89] (Ajuste de potenciales de vibración para moléculas diatómicas), y usando los intervalos que se indican en la tabla con subintervalos de 0.1 u.a. Se puede observar que la concordancia entre los valores de Re y µe obtenidos con nuestro programa y con Gaussian 94W está dentro de los márgenes de error con los que Gaussian calcula estos parámetros. La comparación de los datos Ecalculada y Eajustada se refiere a las energías calculadas en R=Re (⇒ Ecalculada) y las estimadas con el ajuste realizado por UCA_VIB (⇒ Eajustada). La buena concordancia entre ambos datos indica que el intervalo empleado para el ajuste es el adecuado (cuando el intervalo elegido para el ajuste no está bien centrado entorno al mínimo Re, o los datos E(Ri) no están bien balanceados, aparecen discrepancias entre las Ecalculada y Eajustada) e indica también que el programa UCA_VIB, que el laboratorio necesita usar en otras líneas de investigación, funciona correctamente. El uso de la comparación de los valores Re(OMCM) / Re(G94W) y µe(OMCM) / µe(G94W) como test del correcto funcionamiento de nuestro programa viene sugerido por varios motivos. En primer lugar, que en el cálculo de Re se combinan datos procedentes de cálculos a diversas
Capítulo III-84
INTEGRACIÓN CON FUNCIONES DE BASE GAUSSIANAS
distancias. Si cualquiera de estos datos tuviera demasiado error, ello se reflejaría inmediatamente en el valor calculado de Re. En segundo lugar, nos parece más interesante comparar el rendimiento de esta clase de programas en el cálculo de propiedades moleculares medibles (como Re y µe) que con el cálculo de parámetros exclusivamente teóricos. La comparación entre GAUSSIAN y OMCM en el caso de la molécula de Be2 no resulta inmediata, debido a que ambos programas convergen a estados electrónicos distintos. Si el programa G94W se utiliza de manera estándar, da lugar a una curva de energía total con mínimo y configuración
1σ g2 1σ u2 2σ g2 1π u2 , mientras
que el nuestro da una curva sin mínimo con configuración 1σ 2g 1σ u2 2σ 2g 2σ u2 . En la Figura 4 puede verse que los resultados Gaussian corresponden a un estado electrónico excitado, pero con distancia de equilibrio, mientras que el otro estado no la tiene. No es difícil obligar al programa OMCM a converger a la configuración
1σ g2 1σ u2 2σ g2 1π u2 que da el programa Gaussian, y obtener:
Ee (1σ2g 1σu2 2σ2g 1πu2;R e = 3.4141) = −28.9741068
(−28.974107 con G94W) (77)
Para obtener este resultado basta añadir funciones Py ó Pz a la base y arrancar el cálculo con el orbital 1πu biocupado.
Figura 4: Comparación de los resultados obtenidos para la molécula de Be2 con el programa OMCM y con el programa Gaussian 94W.
Capítulo III-85
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
Capítulo III-86
Capítulo IV: INTEGRACIÓN MONTECARLO Y CUASIMONTECARLO CON FUNCIONES DE BASE STO (I)
INTEGRACIÓN MONTECARLO Y CUASIMONTECARLO CON FUNCIONES DE BASE STO (I)
Capítulo IV: INTEGRACIÓN MONTECARLO Y CUASIMONTECARLO CON FUNCIONES DE BASE STO (I) 1. INTRODUCCIÓN 2. MUESTREO 3. INTEGRALES DE SOLAPAMIENTO 4. INTEGRALES MONOELECTRÓNICAS DE ENERGÍAS CINÉTICA Y POTENCIAL 5. INTEGRALES BIELECTRÓNICAS 6. USO DE LAS INTEGRACIONES MC Y CMC EN CÁLCULOS QUÍMICO CUÁNTICOS
1. INTRODUCCIÓN En este capítulo vamos a presentar las generalidades acerca de los métodos de integración Montecarlo (MC) y Cuasimontecarlo (CMC), y los vamos a aplicar en casos con resultados exactos bien conocidos, con objeto de hacer una primera estimación de sus posibilidades. Concretamente analizaremos, primero, una serie de cálculos de integrales individuales, y a continuación el cálculo de propiedades atómicas y moleculares -realizadas con bases STO- en los que se comparan con los resultados exactos varias formas de aplicar las integraciones CMC. Con objeto de introducir los conceptos básicos de una forma práctica, revisaremos en primer lugar el cálculo por el método de MC de una integral en una dimensión. Conviene reescribirla introduciendo una función auxiliar w(x):
Capítulo IV-87
CALCULO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTANDAR
b
I=
∫
b
f(x)
∫ w(x) w(x)dx
f(x)dx =
a
(1)
a
Cuando la función w(x), cumple determinadas condiciones, se llama función de densidad o función de peso. Por motivos que enseguida se verán, esta función debe “parecerse lo más posible” a la función g(x), pero ser fácil de integrar. En el método de MC la integral I se aproxima por un sumatorio de N términos: I≈
1 N
N
f(ξ j )
∑ w(ξ ) j =1
(2)
j
en el que cuantos más puntos ξj se incluyan, más se aproximará el valor obtenido al valor exacto de la integral. Los puntos ξj han de elegirse mediante un muestreo en el que haya una densidad de probabilidad w(ξj) de obtener cada punto ξj. El muestreo del número aleatorio γ puede realizarse por MC o por CMC en cualquiera de sus versiones. Nosotros hemos comparado en nuestros cálculos los resultados de emplear puntos aleatorios y cuasialeatorios generados con los siguientes métodos: a) Montecarlo ordinario b) Método de Halton [Ref. 32] c) Conjuntos de puntos de Hammersley [Ref. 33 y 35] d) Korobov - Haselgrove [Ref. 33 y 34]
2. MUESTREO La eficiencia de la suma (2) como aproximación de una integral depende de varios factores. En primer lugar del número de puntos ξj empleados, pero sobre todo depende de la forma de elegir tales puntos. Según se ha visto en la sección 8 del capítulo I, en el caso del método de MC ordinario el error cometido es: ε =
Capítulo IV-88
σ N
(3)
INTEGRACIÓN MONTECARLO Y CUASIMONTECARLO CON FUNCIONES DE BASE STO (I)
siendo σ2 la varianza del cociente g(x) = f(x) / w(x) y N el número de puntos empleados. Esta relación pone de manifiesto que aumentar la eficiencia aumentando el número N es poco adecuado, pues para conseguir un dígito más de precisión se necesitaría emplear 100 veces más puntos ξj. La mejor opción es elegir la función de peso w(x) de manera que minimice la varianza: σ = 2
∫
f(x) f 2 (x) w(x) dx − w(x) dx 2 w (x) w(x)
∫
2
(4)
En la práctica esta condición equivale a elegir w(x) lo más parecida posible a la función f(x) que quiere integrarse. La limitación principal del empleo de las funciones de peso es que w(x) debe tener una integral resoluble. La forma habitual de muestrear con densidad w(x) requiere [Ref. 54] calcular una integral llamada marginal de la variable ξ: ξ
m arg inal (ξ) = ∫ w(x)dx
(5)
a
en la que “a” representa el límite inferior de los valores de ξ. Asimismo, se necesita un generador de números aleatorios γj comprendidos en 0