CÁLCULO TEÓRICO DE PROPIEDADES MOLECULARES MEDIANTE BASES NO ESTÁNDAR

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

0 downloads 50 Views 2MB Size

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,µ



+

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)



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





∑∑

Kijα −

i=1 j=1



 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α

+



∑ (J

ij

j=1

Capítulo I-24



Kijα )



+

∑J

ij

j=1

(82)

ANTECEDENTES

εiβ

=

Hiiβ





+

(Jij − Kijβ ) +

j=1



∑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 C

α α α pjCqj

]

(88)

]

(89)

(90)

j=1

β Ppq =



∑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ónicasde 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





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

Get in touch

Social

© Copyright 2013 - 2024 MYDOKUMENT.COM - All rights reserved.