Química Computacional: Fundamentos y Métodos
Química Computacional: Fundamentos y Métodos
Moleculares
E. San Fabián
15 de junio de 2023
2
Índice general
I
II ÍNDICE GENERAL
2. Mecánica Molecular. 67
2.1. Introducción. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
2.2. Potenciales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
2.2.1. Ecuación de la energı́a . . . . . . . . . . . . . . . . . . . . . . . 68
2.2.2. Otros ... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
2.3. Amber . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
2.3.1. Introducción. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
2.3.2. Descripción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
2.4. Métodos QM/MM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
Bibliografı́a. 121
ÍNDICE GENERAL III
Fundamentos de la Quı́mica
Computacional
Introducción
Se denomina Quı́mica Computacional (o Teórica) a la obtención de información
estructural de sistemas quı́micos por medio de cálculos matemáticos basados en leyes
fundamentales de la fı́sica.
Ver la página. web: Introduction to Computational Chemistry, en la que David
Young[2] nos muestra, entre otros, los siguientes contenidos:
1
2 Cálculos Computacionales de Estructuras Moleculares
lecular).
Para todo sistema aislado existe una función matemática, tal que en dicha función
se contiene toda la información significativa del sistema. Se la suele denominar
función de estado (funciones de onda) del sistema (Ψ).
Es una función de las coordenadas de las partı́culas que componen el sistema,
y del tiempo: Ψ(q, t), y debe cumplir ciertas condiciones, como ser continua,
simple-evaluada y de cuadrado integrable.
ÂΨ = aΨ (1.1)
El valor de la propiedad α es a.
⟨Ψ | B̂ | Ψ⟩
β = ⟨B̂⟩ = (1.2)
⟨Ψ | Ψ⟩
Por último hay un postulado muy importante que liga la función de onda con el
tiempo:
∂Ψ ∂Ψ 1
iℏ = iℏΨ̇ = ĤΨ o = ĤΨ (1.3)
∂t ∂t iℏ
donde H es el operador Hamiltoniano, que es el asociado con la energı́a del
sistema.
Existen ciertos sistemas en los que Ĥ no depende del tiempo, estos sistemas se
llaman estacionarios. En estos casos, ya que la energı́a cinética no depende del t, es en
los que el potencial depende tan sólo de las coordenadas, y la función de onda Ψ(q, t)
se puede desdoblar como el producto de una función dependiente de las coordenadas
y otra que dependa del tiempo.
Fundamentos de la Quı́mica Computacional. 3
iℏ ∂φ(t)
=W =⇒ φ(t) = Ce−iW t/ℏ
φ(t) ∂t
1
ĤΨ0 (q) = W =⇒ ĤΨ0 (q) = W Ψ0 (q) = EΨ0 (q)
Ψ0 (q)
Esta es la ecuación de autovalores para estados estacionarios.
El de incertidumbre
1 1
∆a · ∆b ≥ |⟨[Â, B̂]⟩| ó ∆x · ∆px ≥ ℏ
2 2
El de correspondencia de Bohr
Ci =< Ψi |Φ >
Alguna nota más sobre las funciones de onda:, cuando decı́amos que depende de las
coordenadas ( y del tiempo), nos referı́amos a todas las coordenadas que pueden tener
las partı́culas que componen el sistemas, las espaciales y las de espı́n.
Y deben de cumplir el principio de antisimetrı́a de Pauli (Deben ser antisimétricas
respecto al intercambio de dos electrones cualesquiera).
4 Cálculos Computacionales de Estructuras Moleculares
Principio variacional.
El Principio variacional dice que ”dado un sistema cuyo hamiltoniano es Ĥ,
si Ψ es cualquier función normalizada que se comporta bien y satisface las
condiciones de contorno, o condiciones lı́mite del problema, entonces:
⟨Ψ Ĥ Ψ⟩ ≥ E0
⟨Φ Ĥ Φ⟩
W (λ1 , λ2 , . . . , λn ) =
⟨Φ | Φ⟩
Ĥ = Ĥ 0 + Ĥ 1
siendo Ĥ 1 la perturbación que introducimos o debemos introducir.
Lo que vamos a hacer es intentar obtener un relación entre las funciones propias del
sistema perturbado y las del sistema sin perturbar. La forma de hacerlo es simulando
que la perturbación se va haciendo paulatinamente, lo cual en lenguaje matemático
equivale a multiplicar la perturbación por un parámetro λ y sumarla al hamiltoniano
sin perturbar:
Ĥ = Ĥ 0 + λĤ 1
como es lógico, para λ = 0, el sistema será el no-perturbado, y después podemos
aumentar el valor de λ hasta λ = 1, en que tendremos el Ĥ del sistema perturbado.
Claro está, si el Ĥ depende de λ , el lógico pensar que la función de onda Ψn y los
valores propios En dependerán también de λ, con lo que los puedo desarrollar en una
serie de Taylor de potencias de λ:
6 Cálculos Computacionales de Estructuras Moleculares
∂Ψn ∂ 2 Ψn λ2
Ψn = Ψn |λ=0 + λ+ + ...
∂λ λ=0 ∂λ2 λ=0 2!
∂En ∂ 2 En λ2
En = En |λ=0 + λ+ + ...
∂λ λ=0 ∂λ2 λ=0 2!
Por hipótesis, si λ = 0
Ψn = Ψ0n En = En0
y si tomamos el criterio de denominar :
∂ k Ψn 1 ∂ k En 1
Ψ(k)
n = En(k) =
∂λk λ=0 k! ∂λk λ=0 k!
entonces :
ĤΨn = EΨn
(En(0) Ψ(0) (0) (1) (1) (0) (1) (1) 2 (2) (0) 2
n + En Ψn λ + . . . + En Ψn λ + En Ψn λ + . . . + En Ψn λ + . . .)
(Ĥ 0 Ψ(0) (0) (0) 0 (1) 1 (0) (0) (1) (1) (0)
n − En Ψn ) + λ(Ĥ Ψn + Ĥ Ψn − En Ψn − En Ψn )
Como λ es un parámetro arbitrario, los coeficientes de λ deben ser igual a cero por
separado y nos queda :
Fundamentos de la Quı́mica Computacional. 7
La primera ecuación es la ecuación secular del sistema sin perturbar, cuya solución
conocemos.
En la segunda ecuación ya aparecen las correcciones de primer orden a la función
y a la energı́a, tratemos de ver como son.
Tenemos que Ĥ 0 es hermı́tico, y el conjunto de funciones propias del sistema sin
perturbar formarán un conjunto completo ortonormal, por lo que podemos desarrollar
(1)
Ψn como una combinación de dichas funciones:
X (0)
Ψ(1)
n = aj Ψj
j
X (0)
aj (Ej − En(0) )δmj = En(1) δmn − ⟨Ψ(0) 1 (0)
m | Ĥ | Ψn ⟩ ⇒
j
(0)
am (Em − En(0) ) = En(1) δmn − ⟨Ψ(0) 1 (0)
m | Ĥ | Ψn ⟩
En ∼
= En(0) + En(1) = En(0) + Hnn
1
(0)
am (Em − En(0) ) = −Hmn
1
Ĥ = T̂ + V̂ = T̂ n + T̂ e + V̂ en + V̂ ee + V̂ nn (1.6)
N
n
X 1 2
T̂ = − ∇A (1.7)
A
2m A
n
X
e 1 2
T̂ = − ∇µ (1.8)
µ
2
n X
N
en
X ZA
V̂ = − (1.9)
µ A
rµA
X′ 1 n n
′ 1
X X
ee
V̂ = = ĥµ + (1.10)
µ<ν
rµν µ µ<ν
rµν
N
′ ZA ZB
X
nn
V̂ = (1.11)
A<B
rAB
y la función de onda dependerı́a de las coordenadas de los núcleo y de las de los
electrones:
Fundamentos de la Quı́mica Computacional. 9
Ψ(RQ , xq ) (1.12)
donde en xq hemos incluido las coordenadas cartesianas y de espı́n del electrón q.
La resolución de ĤΨ = EΨ es inviable, y hemos de hacer algunas aproximaciones,
basadas en modelos fı́sicos que tengan sentido.
El primero de ellos es considerar la gran diferencia entre las masas de los núcleos
y la de los electrones (para el hidrógeno, unas 1836 veces) con lo que es de prever un
movimiento mucho más lento para los núcleos que para los electrones. Esto es lo que
considera la aproximación Born-Oppenheimer, según la cual podemos considerar los
núcleos fijos en unas posiciones determinadas y resolver el problema para los electrones.
Con esta aproximación, el V nn es constante, la T n será nula y se puede escribir la
función de onda como:
Ψ(RQ , xq ) = Φ(RQ )ΨeR (xq ) (1.13)
de forma que ahora puedo resolver:
siendo
ĤRe = T̂Re + V̂Ree + V̂ren (1.15)
Lógicamente la energı́a total será:
anteriormente con este nombre) son con frecuencia, mezclas de dos o mas estructuras
electrónicas de orbitales moleculares o EV simples los cuales tienen funciones de energı́a
que se cruzan, y que son las que usualmente calculamos con un determinante.
La representación adiabática de SEP es la que se obtiene directamente del trata-
miento exacto del movimiento nuclear.
Un ejemplo de esto es la superficie X + H2 (X halógeno), cuando X se aproxima en
la dirección que forma 60 grados con el eje H2 (Las siguientes figuras están tomadas de la ref.
[6]):
ADIABÁTICAS:
(R − Re )2
E(R) = E(Re ) + E ′ (Re )(R − Re ) + E ′′ (Re ) + ... (1.18)
2!
y lo equiparamos a un oscilador armónico, como primera aproximación,
1
U = U (Re ) + Ke (R − Re )2 (1.19)
2
vemos inmediatamente que
Ke = E ′′ (Re ) (1.20)
1/2
1 E ′′ (Re )
νe = (1.21)
2π µ
Hay otros sistemas de coordenadas que tienen un empleo más especı́fico, y que
mencionaremos de pasada:
Hay otros muchas formas de definir las posiciones de los núcleos, por ejemplo in-
dicando sólo las distancias entre núcleos (harı́an falta N (N − 1)/2), o los sistemas
de coordenadas hiper-esféricos (distancias a un origen y ángulos) de aplicación en
estudios estadı́sticos (Monte Carlo) y dinámicos, etc.
n n
1 X Zk
ĥ0i = − ∇2i − (1.28)
2 k
rik
0
Pnque0no es el real del sistema, pero tiene la ventaja respecto al no perturbado Ĥ =
i=1 ĥi , de que considera de algún modo las repulsiones inter-electrónicas.
Pero los electrones son partı́culas indistinguibles y sus funciones tienen que ser
antisimétricas. Consideremos la función de onda de Slater, que cumple con estas con-
diciones.
1
Ψ = √ | ϕ1 ϕ2 . . . ϕn | (1.29)
n!
siendo los ϕi , funciones espı́n-orbital.
La expresión de la energı́a para un sistema de capa cerrada, definido por la función:
1
Ψ = √ | ϕ1 (1)ϕ1 (2)ϕ2 (3)ϕ2 (4) . . . ϕn/2 (n − 1)ϕn/2 (n) | (1.30)
n!
14 Cálculos Computacionales de Estructuras Moleculares
siendo
1
Jij =< ϕi (µ)ϕi (µ) ϕj (ν)ϕj (ν) > (1.33)
rµν
1
Kij =< ϕi (µ)ϕj (µ) ϕj (ν)ϕi (ν) > (1.34)
rµν
además, se tienen que :
{F̂i ϕi = ϵi ϕi } (1.36)
siendo F̂i el operador de FOCK :
n/2 h i
X
F̂i = ĥ0i + Jˆj (i) − K̂j (i) (1.37)
j=1
tal que
n n
1 X Zk
ĥ0i
= − ∇2i − (1.38)
2 k
rik
Z
ϕj (1)ϕj (1)
Jˆj (i)ϕi (µ) = ϕi (µ) dr1 (1.39)
rµ1
Z
ϕj (1)ϕi (1)
K̂j (i)ϕi (µ) = ϕj (µ) dr1 (1.40)
rµ1
y ϵi es la energı́a HF del orbital iesimo, que vale:
n/2
X
ϵi =< ϕi F̂ ϕi >= ϵ0i + [2Jij − Kij ] (1.41)
j=1
Nα Nβ Nα
X X X
F̂α = ĥ + Jˆiα + Jˆiβ − K̂iα
iα iβ iα
Nα Nβ Nβ
X X X
F̂β = ĥ + Jˆiα + Jˆiβ − K̂iβ
iα iβ iβ
ΨRHF = |ϕ1 ϕ1 ϕ2 ϕ2 ϕ3 ϕ3 | ΨUHF = |ϕα1 ϕβ1 ϕα2 ϕβ2 ϕα3 ϕα4 | ΨROHF = |ϕ1 ϕ1 ϕ2 ϕ2 ϕ3 ϕ4 |
Restricted Hartree-Fock Unrestricted Hartree-Fock Restricted Open-Shell Hartree-Fock
(n+1)
ϕi ∼ (n)
= ϕi y
(n+1)
ϵi
(n)
= ϵi (1.46)
con lo que se ha logrado la auto-consistencia, y tendremos los resultados finales :
{ϕSCF
i }{ϵSCF
i } SELF CONSISTENT FIELD (1.47)
Claro que para que todo esto converja, hemos de buscar un buen punto de partida,
(0)
ϕi , es decir, ese punto de partida debe estar razonablemente próximo al exacto, es
ası́ como surge la hipótesis adicional de Hartree, que es la de suponer que el potencial
efectivo tienen simetrı́a esférica, por lo que solo depende de r (distancia al núcleo), y
ası́ empleamos funciones del tipo
Donde α es un parámetro que se puede determinar bien por las reglas de Sla-
ter, bien de forma variacional. Estas funciones presentan una convergencia muy
rápida, pero tienen la contrapartida de que no son ortogonales.
Fueron las primeras utilizadas por Roothaan y Bagus :
m
X
φiλα = χpλα Ciλp
p=1
Las hizo muy completas para átomos, Serafı́n Fraga, y después Clementi y Roetti
para los sistemas Litio-Kripton. McLean la ha ampliado hasta el Radon.
Fundamentos de la Quı́mica Computacional. 17
2
rl e−αr (1.50)
1.4.1. Gaussianas
Dentro de las GTO se trabaja con dos tipos, las GTO esféricas y las GTO’s carte-
sianas, que se definen respectivamente por:
ESFÉRICAS:
2
χpλα (r, θ, φ) = N (npλ , αpλ )rnpλ −1 e−αpλ r Yλα (θ, φ)
1 2npλ +1
N (npλ , αpλ ) = 2npλ +1 [(2npλ − 1)!!](2π)− 4 (αpλ ) 4
CARTESIANAS:2
2
χplmn (x, y, z) = N (l, αp )N (m, αp )N (n, αp )xl y m z n e−αp r
1 2 1 2k+1
N (k, α) = [(2k − 1)!!]− 2 ( ) 4 α 4
π
En las cartesianas se habla de GTO s, p, d,... según el valor l + m + n = 0, 1, 2, ...
respectivamente. Indicar que las GTO cartesianas del tipo d tienen 6 funciones, que
son equivalentes a las 5 GTO’s esféricas y una GTO esférica del tipo 3s.
La contracción de funciones consiste en generar nuevas funciones de base(N ) a partir
de combinaciones lineales adecuadas de un conjunto de funciones de bases primitivas
previamente generado(M ). Con esto reducimos un conjunto de M funciones a N, y
ahora los procesos de construir el operador de Fock y diagonalizarlo, sólo dependerán
de (N 4 ) y (N 3 ); además, cálculos post-SCF, dependerán de potencias de N mayores de
cuatro. Hay dos esquemas de contracción (Ver página 197 de ref. [10]):
Existen muchas clases de funciones GTO’s, veamos algunos de los conjuntos más em-
pleados:
Globalmente, y atendiendo a su utilización tendrı́amos :
Bases de Jensen:
Bases de Ahlrichs:
Pseudo-potenciales:
Ası́ mismo existen funciones base especiales para cálculos especiales, como por ejem-
plo cálculos de interacción espı́n-espı́n.
Según el programa que se utilice, se escriben de una u otra forma
Hay programas con unas interfaces muy agradables, como Ecce-NwChem: ECCE
O 0
S 6 1.00
8588.50000000 0.00189515
1297.23000000 0.01438590
299.29600000 0.07073200
87.37710000 0.24000100
25.67890000 0.59479700
3.74004000 0.28080200
SP 3 1.00
42.11750000 0.11388900 0.03651140
9.62837000 0.92081100 0.23715300
2.85332000 -0.00327447 0.81970200
SP 1 1.00
0.90566100 1.00000000 1.00000000
SP 1 1.00
0.25561100 1.00000000 1.00000000
D 1 1.00
5.16000000 1.00000000
D 1 1.00
1.29200000 1.00000000
D 1 1.00
0.32250000 1.00000000
F 1 1.00
1.40000000 1.00000000
SP 1 1.00
0.08450000 1.00000000 1.00000000
****
Tabla 1.1: Cálculos atómicos con distintos conjuntos de funciones de base, para el
Oxı́geno.
STO-3G 3-21G 6-311G cc-pVTZ aug-cc-pVTZ
Oxı́geno:
Func. de Base 5 9 13 30 46
Func. Primitivas 15 15 26 55 75
Integr. bi. calcul. 33 213 775 775+11153 2073+60346
E(UHF) -73.804150 -74.393657 -74.802496 -74.811756 -74.812982
E(UMP2) -73.804150 -74.443340 -74.861050 -74.954902 -74.959294
E(UMP4) -73.804150 -74.449435 -74.869170 -74.973102 -74.977907
a
cpu time en s. 5.4 5.5 8.3 14.3
a
En un PC AMD Athlon XP 2000+ (cpu MHz=1659),
bajo linux Red Hat, Release 9 (kernel 2.4.20-31.9)
1.4.3. Pseudo-potenciales
Sólo consideran los electrones de valencia moviéndose en el potencial generado por
el núcleo y los electrones del core:
22 Cálculos Computacionales de Estructuras Moleculares
Tabla 1.2: Cálculos atómicos con distintos conjuntos de funciones de base, para el
Oxı́geno. UHF, MP4(SDTQ) y B3LYP. (t) es el tiempo de cpu, en segundos, en 12
procesadores Intel(R) Xeon(R) CPU E5-2670 0 2.60GHz.
Oxı́geno:
Base NºFunc. E(UHF) E(MP4)(t) EM P4
corr E(B3LYP)(t) EB3LY
corr
P
al de pseudo-potenciales:
nv X nv
ps
X 1 2 ps 1
Ĥ = − ∇i + Vi +
i
2 r
i<j ij
donde Vl (r) es una función de r y Pl representa el proyector sobre los armónicos esféricos
de simetrı́a l.
o, la no-local:
Z X
V ps = − + Cpq |fp >< fq |
r p,q
LANL2DZ ECP
24 Cálculos Computacionales de Estructuras Moleculares
CRENBL ECP
Estos ECPs se llaman también ”consistentes en forma”, porque mantienen la forma
de los orbitales atómicos en la región de valencia.
CRENBS ECP
(b) 3d10 4s2 4p3 5d10 6s2 6p3 3d10 4s2 4p4 5s2 5p6 5d9 6s1
(*) Un sólo ciclo.
(Func. Base Primitivas time )
Fundamentos de la Quı́mica Computacional. 25
Con las bases atómicas {χA } y {χB }, para el cálculo de AB, utilizaremos una base
{χA +χB }, por lo que el sistema AB estará mejor descrito que los sistemas individuales
A y B.
La forma más usual de corregir parcialmente este error, o por lo menos de conside-
rarlos es el método de Boys y Bernardy counterpoise method, que consiste en realizar
los cálculos de los átomos con la base completa de la molécula.
0 1 0 1 0 1
H(Fragment=1) 0.00000000 0.00000000 0.58022808
Br(Fragment=1) 0.00000000 0.00000000 -0.83659185
F(Fragment=2) 0.00000000 0.00000000 2.77788358
H(Fragment=2) 0.00000000 0.00000000 3.69953441
Tabla 1.4: Error de superposición de bases, para cálculos MP4 con diferentes conjuntos
de funciones de base, del sistema OH-H, con Re =0.958 Å, α(HOH)= 104.4776 0
{⟨Tijk |Ĥ|Sia ⟩}
abc abc ab abc abc abc
{0} {⟨Tijk |Ĥ|Dij ⟩} {⟨Tijk |Ĥ|Tijk ⟩} {⟨Tijk |Ĥ|Q⟩} · · · · · ·
{0} {0} ··· ··· ··· ··· ···
El truncamiento del desarrollo Ψf ull−CI en la forma ΨCISD o cualquier otra forma que
no incluya todas las posibles excitaciones (i.e., que no sea full-CI) conduce a un error
que, en inglés, se conoce como size consistency3 y que podrı́amos llamar en español
3
Size extensivity: Se refiere al correcto escalado lineal de un método con el número de electrones.
HF, MBPT, CC y full-CI son size-extensive, CI truncados no. Si es size-extensive, es size-consistent,
pero la autoconsistencia tiene el requisito de la correcta fragmentación
Funciones de onda multideterminantales. 29
total
√ monomero
Ecorr (CISD) = N Ecorr (CISD)
en lugar de la relación lineal que deberı́a tener si el cálculo fuera consistente con
el tamaño del sistema. Existen métodos que permiten corregir este problema, entre
los cuales el más usado es el método CI cuadrático (QCISD), que se deriva del CISD,
incluyendo los términos de mayor orden que incluye el CCSD, aunque no todos. Esto
hace que sea de coste similar al CCSD y los resultados también son análogos, por lo
que casi es más coherente utilizar CCSD que QCISD.
Los métodos CI en general no son los más utilizados en el cálculo rutinario de la
energı́a de correlación por varias razones. Fundamentalmente, estas razones son:
La falta de consistencia de tamaño.
La convergencia del desarrollo CI es muy lenta, por lo cual la cantidad de energı́a
de correlación que se recupera con los tratamientos más simples es pequeña.
Es un método costoso, porque para el cálculo de la energı́a (que en el caso CISD se
representa como
XX
ECISD = EHF + Cijab [(ij|ab) − (ia|jb)] (1.53)
i<j a<b
(donde i,j son orbitales ocupados y a,b orbitales virtuales) necesita el cálculo de las
integrales moleculares, que a su vez están relacionadas con las integrales sobre funciones
de base en la forma
X
(ij|ab) = Ciµ Cjν Caλ Cbσ (µν|λσ)) (1.54)
µνλσ
Tabla 1.5: Energı́as de una y dos moléculas de agua a 25 Å, con la geometrı́a: ángulo
HOH= 104.5 y Re =0.957 Å. Base Def2-SVP.
Método H2 O 2 H2 O a 25 Å Diferencia
HF -75.961025 -151.922045 0.00000
MP2 -76.162108 -152.324211 0.00000
MP3 -76.169076 -152.338148 0.00000
MP4 -76.174329 -152.348654 0.00000
CISD -76.163652 -152.309181 0.01812
CCSD -76.171715 -152.343426 0.00000
CCSD(T) -76.174693 -152.349382 0.00000
∆EDavidson = (1 − C02 )∆ECISD = −0.005061 (H2 O)
∆EDavidson = (1 − C02 )∆ECISD = −0.016835 (2H2 O)
CISDTQ = -152.32 ≈ −152.309181 − 0.016835 = −152.326016
Podemos hacernos una idea de los resultados CI, analizando los resultados que presen-
tan Helgaker et al. (pag 183 de ref. [12]), para el H2 O a dos distintas distancias OH:
Tabla 1.6: Diferencias de energı́a respecto a la FCI, calculadas con la base cc-pVDZ,
para el H2 O, a dos geometrı́as, con ángulo HOH= 110.565 y Rref = 1.84345 a0 . W es
el peso de las funciones CI truncadas en la FCI.
R=Rref R=2·Rref
E-EF CI W E-EF CI W
RHF 0.217822 0.941050 0.363954 0.589664
CISD 0.012024 0.998047 0.072015 0.948757
CISDT 0.009043 0.998548 0.056094 0.959086
CISDTQ 0.000327 0.999964 0.005817 0.998756
CISDTQ5 0.000139 0.999985 0.002234 0.999553
CISDTQ56 0.000003 1.000000 0.000074 0.999993
ERHF -76.024038 -75.587711
EF CI -76.241860 -75.951667
Funciones de onda multideterminantales. 31
H0 Ψ(0) = E (0) Ψ(0) E (0) =< Ψ(0) |H0 |Ψ(0) > (1.56)
mientras que el operador H1 , multiplicado por un cierto parámetro λ que finalmente
haremos igual a 1, es la perturbación, que suponemos pequeña respecto al operador
de orden cero (y que será, en el caso que consideramos, la correlación electrónica). A
continuación entonces podemos expresar tanto la energı́a como la función de onda del
Hamiltoniano exacto (i.e., el que incluye la correlación) como una serie en potencias
del parámetro λ
Introduciendo estas expresiones en la ecuación de Schrödinger e igualando término
a término los coeficientes de las distintas potencias obtenemos las ecuaciones que nos
permiten obtener las correcciones perturbativas pertinentes.
Esto sirve para cualquier perturbación. Cuando H0 es el operador de Fock, la con-
creción de la MBPT planteada en las ecuaciones anteriores recibe el nombre de teorı́a
de perturbaciones de Møller-Plesset (MPPT) y es el método más comúnmente emplea-
do para calcular la energı́a de correlación. En particular, la MPPT de segundo orden,
que recibe el nombre de MP2 está programada en forma tal en la mayor parte de los
programas de cálculo que su evaluación es hoy en dı́a muy rápida. En el caso de MP2,
la energı́a de correlación toma la forma
Tabla 1.7: Energı́as MPn, calculadas con Gaussian09 (FC, 5d,7f), para el H2 O, a dos
geometrı́as, con ángulo HOH= 110.565 y Rref = 1.84345 a0 .
cc-pVDZ aug-cc-pVDZ
R=Rref R=2·Rref R1 =Rref , R2 3·Rref R=Rref R=2·Rref
RHF -76.0240386 -75.5877113 -75.6808813 -76.0389405 -75.6002312
MP2 -76.2264254 -75.8950057 -76.0037391 -76.2595807 -75.9314538
MP3 -76.2333241 -75.8808395 -76.0459271 -76.2641138 -75.9068279
MP4 -76.2387068 -75.9338987 -76.0941739 -76.2732235 -75.9688048
MP5 -76.2392642 -75.9332675 -76.0957760 -76.2719741 -75.9524374
FCIa -76.2397611 -75.9499394 -76.0604790 -76.2733209 -75.9750710
a
Con Psi4, FC y esféricas (unos 400 s)
e−T̂ ĤeT̂ |ΨCC >= e−T̂ ECC |ΨCC >= ECC e−T̂ |ΨCC >= ECC |ΨHF >
ECC =< ΨHF |e−T̂ HeT̂ |ΨHF > =< ΨHF |HeT̂ |ΨHF >
Donde ahora se trunca la exponencial
1
(1 + T̂ + T̂ 2 + · · · )
2
y también
T̂ = T̂1 + T̂2 + · · ·
De la misma forma que sucedı́a con la CI, realizar un tratamiento completo del
problema es imposible cuando el sistema no es muy pequeño. De la misma forma que
la CI, 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 más frecuente es la
que incluye sólo el operador de dobles excitaciones, en la forma
1X 1 X ab † †
T̂2 = t̂ij = t â â âi âj (1.60)
2 ij 4 ijab ij a b
donde X
t̂i ≡ tαi â†α âi (1.61)
α
X
† †
t̂ij ≡ tab
ij âa âb âj âi (1.62)
a>b
con
â†p |ϕq · · · ϕs >= |ϕp ϕq · · · ϕs > (1.63)
âp |ϕp ϕq · · · ϕs >= |ϕq · · · ϕs > (1.64)
(En general:
2 X
1 † †
T̂n = tab...
ij... â b̂ ...âi âj ... (1.65)
n! ij...ab...
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.
En efecto, de la ecuación (1.58) obtenemos
1 2
ΨCCD = Ψ0 1 + T̂2 + T̂2 + .... (1.66)
2
34 Cálculos Computacionales de Estructuras Moleculares
X X
ΨCCD = Ψ0 + tab ab
ij Ψij + tabcd abcd
ijkl Ψijkl + ... (1.67)
i,j,a,b i,j,k,l,a,b,c,d
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 expresión anterior también difiere de la 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 cuádruple excitaciones,
etc. Consecuentemente, la energı́a CCD recupera mucha más energı́a de correlación
que la MP2, por ejemplo, y converge mucho más rápidamente que esta serie.
En este caso, las ecuaciones de acoplamiento son:
0 = ⟨Φai |Ĥ 1 + T̂2 |ΦHF ⟩ (1.68)
1 2
0 = ⟨Φab ij |Ĥ 1 + T̂2 + T̂2 |ΦHF ⟩ (1.69)
2
El único defecto grave de CC respecto a MPn es que resulta mucho más costosa de
calcular, por lo cual las optimizaciones de geometrı́a usando CCD, por ejemplo, son
mucho menos frecuentes en la literatura que las obtenidas usando MP2.
Ver: [Link] .
Para un cálculo CCSD tendremos sólo el acoplamiento con simples y dobles, y el
resto de los términos se anularán, quedándonos dos conjuntos de ecuaciones acopladas
y no lineales:
a 1 2 1 3
0 = ⟨Φi |Ĥ 1 + T̂1 + T̂2 + T̂1 + T̂1 T̂2 + T̂1 |ΦHF ⟩ (1.70)
2 3!
1
0 = ⟨Φij |Ĥ 1 + T̂1 + T̂2 + T̂12 + T̂1 T̂2 +
ab
2
1 3 1 2 1 2 1 4
T̂ + T̂ + T̂ T̂2 + T̂1 |ΦHF ⟩ (1.71)
3! 1 2 2 2 1 4!
Debido al aumento de su complejidad, a partir del nivel CCSD, surgen una gran
cantidad de tratamientos aproximados, entre los que destacan los efectuados con las
triples excitaciones
CCSDT, CCSD[T], CCSD(T), CCSDT-n Ası́, el método CCSDT es muy costoso,
por lo que suelen utilizar los aproximados:
CCSD(T), como la energı́a CCSD más una estimación de las triples:
Tabla 1.8: Energı́as CC, calculadas con la base cc-pVDZ, con Cfour-MRCC, full elec-
trons, para el H2 O, a dos geometrı́as, HOH= 110.565 y Rref = 1.84345 a0 .
Método Energı́a tiempo
R=Rref R=2·Rref
RHF -76.0240385 -75.5877113
CC2 -76.2296059 -75.9145119 0.5 s
CCSD -76.2381164 -75.9296328 1. s
CC3 -76.2412740 -75.9528090 7.2 s
CCSD(T) -76.2412018 -75.9554852 1. s
CCSDT -76.2413670 -75.9530698 12 s
CC4 -76.2418754 -75.9532675 122 s
CCSDT(Q) -76.2418729 -75.9533325 33 s
CCSDTQ -76.2418413 -75.9516351 145 s
CCSDTQP -76.2418575 -75.9516455 5631 s
FCIa -76.2418601 -75.9516669 12062 s
a
Con Psi4, Full y esféricas
Tabla 1.9: Energı́as del Etileno, Base 6-311G**. Geometrı́a optimizada al nivel CCSD.
Cálculos efectuados con G09, Cfour-MRCC, full y sobre Intel(R) Xeon(R) CPU E5-
2670 0 @ 2.60GHz
Método Energı́a (CCn) Tiempo (CCn)
HF -78.053914 10”
MP2 -78.381568 35”
MP3 -78.406276 60”
MP4 -78.421513 82”
MP5 -78.421628 70’
CCSD -78.411101 (-78.384761) 2.2’[13.2’] (6”[16”])
CCSD[T] -78.423379
CCSD(T) -78.423006 2.8’
CCSDT -78.423471 (-78.423406) 186’ (25’)
CCSDT[Q] -78.423649
CCSDT(Q) -78.424104 604’
CCSDTQ -78.424072 (-78.42408 ) 1500’ [8493’] (6000’)
Funciones de onda multideterminantales. 37
Vamos, que es carı́simo y se puede optar por otros métodos más intuitivos, más
elaborados, y menos costosos.
Los métodos multi-configuracionales 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 multi-configuracional es el llamado SCF multiconfi-
guracional (MCSCF). En este método se seleccionan ciertas configuraciones electrónicas
(i.e., determinantes de Slater) formados en base a una cierta configuración de referencia.
38 Cálculos Computacionales de Estructuras Moleculares
(En la página 119 de ref.[3] podéis ver una variante del CASSCF, el denominado
RASSCF, que es otra forma de generar las configuraciones en una ventana de OMs
activos.)
Lo usual es resolver el problema MCSCF según el esquema de Roos en “dos pasos”
(Roos BO (1983) In: Diercksen GHF, Wilson S. eds. Methods in computational mole-
cular physics. Reidel, Dordrecht, Holland, pp 161-187). Es decir, se hace por separado
la optimización de los coeficientes CI y la de los orbitales.
Los coeficientes CI se obtienen por un CI convencional. Se considera que
X
Ψ(n) = Ca(n) Φa
a
con nti son los números de ocupación (con un valor de 0 ó 1) de los espı́n-orbitales i y
uat son los coeficientes de acoplamiento por simetrı́a, del estado representado por Ψ(n) .
pero este Fi es distinto del operador de Fock que obtenı́amos cuando utilizábamos
una función monodeterminantal (o mono-configuracional), ahora nos aparecen unas
expresiones más complejas, que vamos a obviar, y finalmente, se puede mostrar que
1X
Etot = < i|ni ĥ + F̂i |i >
2 i
donde X
ni = nia Ca2
a
4. Mejora de los orbitales moleculares (formación del gradiente y del hessiano de los
orbitales y resolución de las ecuaciones de Newton-Raphson).
40 Cálculos Computacionales de Estructuras Moleculares
Tabla 1.12: Tiempo de cálculo en un ordenador con procesador Intel(R) Xeon(R) Gold
6150 CPU @ 2.70GHz, (sólo un procesador para el Gaussian09 y Psi4). H2 O base
cc-pVDZ. Geometrı́a, la de siempre.
N N N N
1 X 2 XX 1 X
Ĥ = − ∇ + + v(ri ) (1.72)
2 i=1 i=1 j>i
|ri − rj | i=1
1.7.1. Funcionales
Un funcional es similar a una función, que relaciona una variable x con un valor y,
pero donde la cantidad y, en vez de depender de la coordenada x, depende de una (o,
en general, más de una) función ϕ(x). Más coloquialmente, un funcional es una función
cuya variable es otra función. El conjunto de funciones admisibles F constituye el
dominio del funcional F [ϕ].
De acuerdo con lo anterior, la energı́a de un sistema es un funcional de la función
de onda,
Z
E = E[Ψ] = Ψ(x1 , . . . , xN )ĤΨ∗ (x1 , . . . , xN ) dx1 · · · dxN (1.73)
La minimización del funcional E[Ψ] sobre todo el espacio de Hilbert permite obtener
la energı́a del estado fundamental E0 y su función de onda Ψ0 ,
δE
= 0. (1.78)
δϕi (r)
Z
Γ(x′1 , x′2 ; x1 , x2 ) = N (N − 1) ΓN (x′1 , x′2 , x3 , . . . , xN ; x1 , x2 , x3 , . . . , xN ) dx3 · · · dxN
(1.80)
N (N −1)
Existen varios criterios de normalización (por ejemplo, a utilizado por 2
,
Löwdin, a N (N − 1), utilizado por McWeeny) y hay que saber cual se utiliza.
La matriz densidad reducida de primer orden se puede obtener por reducción de la
de segundo orden,
Z
′ 1
γ(x ; x) = Γ(x′1 , x2 ; x1 , x2 ) dx2 (1.81)
N −1
Las densidades electrónica de espı́n α y β, y la densidad electrónica total del sistema
se pueden definir respectivamente como
Especı́ficamente, la densidad ρ asociada con una función de onda Ψ viene dada por
Z
ρ(r) = Ψ(x1 , x2 , . . . , xN )Ψ∗ (x1 , x2 , . . . , xN ) dσ1 dx2 · · · dxN (1.85)
N
X Z
⟨Ψ| v(ri )|Ψ⟩ = ρ(r)v(r) dr (1.86)
i=1
Z Z Z
1 2 ′
1
Γ(x1 , x2 ; x1 , x2 )
E = E[Γ] = − ∇ γ(x ; x) x′ =x dx + v(r)ρ(r) dr + dx1 dx2
2 |r2 − r1 | 2
(1.87)
Para el cálculo de la energı́a cinética se ha adoptado la convención usual de que,
antes de la integración, un operador actúa sólo sobre x, pero no sobre x′ , y, una vez
que ∇2 ha actuado sobre γ, se iguala x′ a x y se integra.
46 Cálculos Computacionales de Estructuras Moleculares
El principio variacional nos dice que, para la función de onda del estado fundamental
Ψ′0 de cualquier sistema con Hamiltoniano Ĥ ′ , se cumple que
Z
⟨Ψ′0 |Ĥ ′ |Ψ′0 ⟩ − [v ′ (r) − v(r)] ρ′ (r) dr ≥ ⟨Ψ0 |Ĥ|Ψ0 ⟩, (1.91)
N N
1X 2 X
Ĥs = − ∇ + v(ri ) (1.94)
2 i=1 i=1
El Funcional Ts [ρ]
48 Cálculos Computacionales de Estructuras Moleculares
Es muy sencillo escribir una expresión exacta para Ts [ρ] en términos de un con-
junto de N orbitales {ϕi } restringidos a
• Ser ortonormales, Z
ϕi (r)ϕj (r) dr = δij (1.97)
1 ρ(r)ρ(r′ )
Z
EXC [ρ] = Q[ρ] − ′
drdr′ + T [ρ] − Ts [ρ]. (1.99)
2 ∥r−r ∥
El término que sigue a Q[ρ] (El hipotético funcional Universal de la densidad) es la
repulsión Culombiana clásica. Como tanto Q[ρ], como la repulsión Culombiana clásica,
como Ts [ρ] son funcionales de la densidad, también lo es EXC [ρ].
Queremos hacer algunas observaciones sobre el funcional de intercambio y correla-
ción:
En la definición de EXC [ρ] intervienen dos sistemas, un sistema interactuante para
el cálculo de Q[ρ], y un sistema no interactuante para el cálculo de Ts [ρ]. Ambos
sistemas tienen la misma densidad ρ
Como el funcional Q[ρ] es desconocido, también lo es EXC [ρ]. Hay que recurrir a
aproximaciones.
Si se conociera la forma explı́cita de EXC [ρ], se podrı́a obtener la energı́a y den-
sidad de cualquier sistema interactuante por minimización del funcional
1 ρ(r)ρ(r′ )
Z Z
E[ρ] ≡ v(r)ρ(r) dr + drdr′ + Ts [ρ] + EXC [ρ]. (1.100)
2 ∥ r − r′ ∥
A continuación veremos el procedimiento para minimizar esta última expresión, que
da lugar a las conocidas ecuaciones de Kohn y Sham.
La ecuación de Euler-Lagrange nos dice que para obtener la ecuación que debe
cumplir el orbital ϕi hay que igualar la derivada funcional de E[{ϕi }, {ϕ̄i }] con respecto
a ϕi a cero. De esta igualación surgen las ecuaciones de Kohn y Sham,
1 2
− ∇ + v(r) + φ(r) + µXC (r) ψi (r) = ϵi ψi (r) (1.106)
2
para orbitales α, y
1 2
− ∇ + v(r) + φ(r) + µ̄XC (r) ψ̄i (r) = ϵ̄i ψ̄i (r), (1.107)
2
para orbitales β. En esta ecuaciones4 , φ(r) es el potencial de Coulomb,
ρ(r′ )
Z
φ(r) = ′
dr′ (1.108)
∥r−r ∥
Finalmente, vemos que hay un término que depende del espı́n. Es el potencial de
intercambio y correlación, que se define, para electrones α, como la derivada funcional
de Exc [ρα , ρβ ] con respecto a ρα ,
δEXC
µXC (r) = , (1.109)
δρα (r)
4
Obsérvese que no aparecen multiplicadores de Lagrange del tipo ϵij y ϵ̄ij , como en (1.105). Ello es
debido a que las ecuaciones de Kohn y Sham (como ocurre en el caso Hartree-Fock) son invariantes bajo
una transformación unitaria de los orbitales, y siempre es posible poner la matriz de multiplicadores
de Lagrange en forma diagonal.
Métodos del Funcional de la Densidad. 51
donde la función ϵXC depende paramétricamente de ρα , ρβ , ρα,x , ρβ,x , ρα,y , ρβ,y , ρα,z
y ρβ,z . Las ultimas seis cantidades son derivadas parciales de la densidad. Ası́,
∂ρα
. ρα,x ≡ (1.112)
∂x
Las otras cinco cantidades se definen análogamente.
El potencial de intercambio y correlación correspondiente a este funcional aproxi-
mado viene dado por
∂(ρϵXC ) ∂ ∂(ρϵXC ) ∂ ∂(ρϵXC ) ∂ ∂(ρϵXC )
µXC (r) = − − − , (1.113)
∂ρα ∂x ∂ρα,x ∂y ∂ρα,y ∂z ∂ρα,z
con una expresión análoga para µ̄XC (r).
El procedimiento habitual es la separación del potencial de intercambio-correlación
en dos, uno de intercambio y otro de correlación, ası́ han surgido innumerables fun-
cionales, tanto de intercambio, como de correlación, unos ligados entre sı́, y a menudo
utilizados indistintamente. En el apéndice A se muestra una lista, siempre sin actuali-
zar, de un conjunto de funcionales, tanto de intercambio, como de correlación.
δE
(con una expresión análoga para el componente β). Al término δραC se le denomina
potencial de correlación y se le representa por µC (r). Resulta que la derivada funcional
de EX [ρα , ρβ ] con respecto a ρα no es otra cosa que el operador de intercambio Kˆα de la
teorı́a Hartree-Fock (la demostración es inmediata como consecuencia de la definición
dada a EX [ρα , ρβ ]). Esto nos permite escribir las ecuaciones de Kohn y Sham (1.106)
como
1 2 ˆ
− ∇ + v(r) + φ(r) + Kα + µC (r) ψi (r) = ϵi ψi (r), (1.116)
2
y recordando la definición del operador de Fock Fˆα obtenemos:
h i
ˆ
Fα + µC (r) ψi (r) = ϵi ψi (r), (1.117)
o, con otras palabras, cuando se utiliza el funcional de intercambio exacto, las ecua-
ciones de Kohn y Sham toman la misma forma que las ecuaciones de Hartree-Fock,
pero con un término perturbativo debido a la energı́a de correlación (para electrones
β, las ecuaciones son análogas). Estas ecuaciones se denominan las ecuaciones de Kohn
y Sham con intercambio exacto 5 . También llamadas por Parr ecuaciones de Hartree-
Fock-Kohn-Sham.
que los funcionales utilizados son bastante insensibles a cambios en la densidad (con
mejora del conjunto de base, o el uso de la densidad exacta o de una densidad de mayor
calidad en vez de la densidad Hartree-Fock).
Otra caracterı́stica del método Kohn y Sham con intercambio exacto es que si trata-
mos el intercambio exactamente, podemos estudiar la calidad de las diversas aproxima-
ciones existentes para la correlación en forma “pura”, sin contaminación procedente de
un intercambio aproximado. Aproximar a la vez el intercambio y la correlación puede
dar una impresión errónea sobre sus respectivas cualidades, porque, por ejemplo, en
la aproximación local, el intercambio está infravalorado, mientras que la correlación es
sobrevalorada, lo que resulta en una cancelación accidental de los errores de ambos.
GGA: Funcionales que además, hacen uso de del gradiente de la densidad (∇ρα , ∇ρβ )
(Generalized Gradient Approximation).
Con Dispersión: Que incluyen términos de dispersión de largo rango, tipo van der
Waals. (GD2, D3, GD3BJ)
[Link]
Métodos del Funcional de la Densidad. 55
Penúltimo Penúltimo
Corrección de rango-separado: LC-
Para que hacerse una idea de la fecha y sus constructores:
Los de intercambio:
S: Slater, 74.
XA: Xα (Slater)
B : Becke, 86.
O: Handy, 01.
Los de correlación:
Funcionales hı́bridos:
B3LYP : Becke Three Parameter Hybrid Functional, la han liado con poner mal
LYP, es decir sumar una contribución local y la LYP , que la transforman en
LYP-local, 93.
M06 , M06HF, M062X, M11: Hı́bridos de Truhlar and Zhao, 06, 08, 11.
Métodos del Funcional de la Densidad. 57
Otros funcionales :
M11 ,M11L M12SX , MN12SX : Nuevos hı́bridos del grupo de Truhlar, 11, 12.
Tabla 1.14: Energı́as, calculadas con la base cc-pVDZ, para el H2 O, a dos geometrı́as,
con ángulo HOH= 110.565 y Rref = 1.84345 a0 .
R=Rref
RHF -76.024039
CISD -76.230237 MP2 -76.226425 CC2 -76.229606 PBEPBE -76.332323
CISDT -76.232818 MP3 -76.233324 CCSD -76.238116 B3LYP -76.418951
CISDTQ -76.240534 MP4 -76.238707 CC3 -76.241274 M062X -76.387050
CISDTQ5 -76.241722 MP5 -76.239264 CCSD(T) -76.241202 HCTH407 -76.409733
CISDTQ56 -76.241858 CCSDT -76.241367
CC4 -76.241875 CASSCF(4,4) -76.076027
CCSDTQ -76.241841 CASSCF(6,5) -76.076568
CCSDTQP -76.241858
MC(4,4)-PDFT -76.305074
FullCI -76.241861 MC(6,5)-PDFT -76.306737
R=2·Rref
RHF -75.587711
CISD -75.879650 MP2 -75.895006 CC2 -75.914512 PBEPBE -76.018813
CISDT -75.895571 MP3 -75.880840 CCSD -75.929633 B3LYP -76.080167
CISDTQ -75.945848 MP4 -75.933899 CC3 -75.952809 M062X -76.013207
CISDTQ5 -75.949431 MP5 -75.933268 CCSD(T) -75.955485 HCTH407 -76.087846
CISDTQ56 -75.951591 CCSDT -75.953070
CC4 -75.953268 CASSCF(4,4) -75.816825
CCSDTQ -75.951635 CASSCF(6,5) -75.817066
CCSDTQP -75.951646
MC(4,4)-PDFT -76.007219
FullCI -75.951665 MC(6,5)-PDFT -76.011964
60 Cálculos Computacionales de Estructuras Moleculares
Tabla 1.15: Diferencia de energı́as entre dos geometrı́as del H2O (2Rref y Rref ), calcu-
ladas con la base cc-pVDZ, con HOH= 110.565 y Rref = 1.84345 a0 .
RHF 0.436328
CISD 0.350587 MP2 0.331419 CC2 0.315094 PBEPBE 0.313510
CISDT 0.337247 MP3 0.352484 CCSD 0.308483 B3LYP 0.338784
CISDTQ 0.294686 MP4 0.304808 CC3 0.288465 M062X 0.373843
CISDTQ5 0.292291 MP5 0.305996 CCSD(T) 0.285717 HCTH407 0.321887
CISDTQ56 0.290267 CCSDT 0.288297
CC4 0.288607 CASSCF(4,4) 0.259202
CCSDTQ 0.290206 CASSCF(6,5) 0.259502
FulCI 0.290196 CCSDTQP 0.290212 MC(6,5)-PDFT 0.294773
Métodos semiempiricos. 61
< sA sA |sB sB >=< sA sA |pB pB >=< pA pA |pB pB >=< AA|BB > (1.121)
En MINDO/3:
ZB < µµ|BB >= ZB < AA|BB > (1.124)
MINDO/3:
En MNDO es similar:
Tabla 1.16: Errores promedios sin el signo. Aplicado a unos 7600 sistemas. Tomado de:
Accuracy of PM6
Propiedad AM1 PM3 PM6
δHf (kcal/mol) 22.86 18.20 8.01
Distancias enlace (Å) 0.130 0.104 0.091
Ángulos 8.77 8.50 7.86
µ (D) 0.67 0.72 0.85
P.I. (eV) 0.63 0.68 0.50
Barreras de activación.
Tabla 1.18: Cálculos realizados para la molécula H2 O con distintos métodos. Error (ε).
Base aug-cc-pVTZ.
Método (tiempo) Energy RO−H (Å) ε αHOH () ε µ (Debye) ε IE (eV) ε EA (eV)
MINDO3 (1) -0.085457 0.949 -0.009 103.72 -0.756 0.831 -1.025 12.375 -0.246 5.317
AM1 (1.3) -0.094384 0.961 0.003 103.54 -0.933 1.861 0.006 12.010 -0.611 5.977
PM6 (1.3) -0.086406 0.949 -0.009 107.61 3.128 2.068 0.213 11.661 -0.960 5.986
HF (4.5) -76.061203 0.941 -0.017 106.37 1.895 1.939 0.084 11.046 -1.575 0.799
B3LYP (5.1) -76.466197 0.962 0.004 105.10 0.621 1.847 -0.008 12.783 0.162 0.446
M062X (5.6) -76.430092 0.959 0.001 105.32 0.847 1.895 0.040 12.789 0.168 0.637
MP2 (8.2) -76.328992 0.961 0.003 104.10 -0.373 1.860 0.005 12.859 0.238 0.639
MP4 (19.2) -76.343678 0.963 0.005 103.98 -0.495 12.727 0.106 0.616
CISD (15.4) -76.322396 0.955 -0.003 104.64 0.164 1.877 0.022 12.442 -0.179 0.671
CCSD (22.7) -76.333670 0.959 0.001 104.50 0.025 1.734 -0.121 12.576 -0.045 0.634
CCSD(T) (31.9) -76.342326 0.962 0.004 104.13 -0.352 12.668 0.047 0.608
Exp 0.958 104.48 1.855 12.621
Capı́tulo 2
Mecánica Molecular.
2.1. Introducción.
Los cálculos de mecánica molecular o campo de fuerza están basados en modelos
simples de mecánica clásica de estructura molecular. La mecánica molecular nunca
puede ser considerada como una aproximación exacta a la estructura quı́mica de una
molécula.
La mecánica molecular trata las moléculas como una matriz de átomos goberna-
dos por un grupo de funciones de potencial de mecánica clásica. (Ver la página web:
Molecular Mechanics Methods)
Se basan fundamentalmente en la consideración del oscilador, armónico o anarmóni-
co, para la descripción del enlace molecular, aplicándolo a la distancia de enlace, al
ángulo de enlace e incluso a los ángulos diedros. Esto se puede complementar con
algún tratamiento de interacciones electrostáticas o de dipolo entre las cargas puntales
de los átomos que forman el enlace.
(Ver Cap.2 de ref.[3] y Cap. 2 de ref.[11])
Avogadro, Amber y CHARMM, son programas que utilizan los métodos de mecáni-
ca molecular.
Los problemas que plantea son los propios de un método cuyo modelo es muy
aproximado, y depende de muchos parámetros, no sabiendo a priori si son los correctos
para nuestro problema. Otro problema es el ajustar las interacciones electrostáticas a los
de Mecánica Molecular, (sólo el AMBER lo hace más o menos bien), y la consideración
de la anarmonicidad es muy importante.
2.2. Potenciales
Entre los más conocidos y que mejor funcionan están los MM2, MM3 y MMX.
Todos ellos usan términos anarmónicos para la distancia y el ángulo de enlace. Otro
método que proporciona resultados similares o mejores es el MMFF94, (Merch Mo-
lecular Force Field), que está parametrizado a partir de cálculos ab initio, y aunque
es análogo al MM3, está más dirigido al cálculo de proteı́nas y sistemas de interés
biológico. El AMBER es tiene su propio método de mecánica molecular, pero sólo con
términos armónicos y diagonales, es decir sin términos cruzados ángulo-distancia, que
sı́ tienen los MM anteriores. Sin embargo presenta la mejora de los parámetros de car-
67
68 Cálculos Computacionales de Estructuras Moleculares
gas atómicas. Este programa se diseñó para el cálculo de biomoléculas y modela bien
los enlaces de hidrógeno. Por último, el CHARMM de Karplus, está parametrizado con
datos experimentales.
(min.)
=0
z }| { 1 1
U (R) = U (Re ) + U ′ (Re ) (R − Re ) + U ′′ (Re )(R − Re )2 + U ′′′ (Re )(R − Re )3 + · · ·
2 3!
El primer y segundo término son cero, uno por elegirlo ası́ , el otro por ser un
mı́nimo, y si se trunca al primer orden no nulo, nos queda la muy simple expresión de
la energı́a potencial vibracional entre los átomos A y B:
1
U (RAB ) = KAB (RAB − RABe )2
2
Esta ecuación sólo es útil para distancias próximas a las distancias de equilibrio,
si se quiere utilizar a distancias mayores, se emplean truncamientos a tercer orden o
incluso cuarto orden:
1 (3) (4)
U (RAB ) = [KAB + KAB (RAB − RABe ) + KAB (RAB − RABe )2 ](RAB − RABe )2
2
También se puede considerar la expresión de la ecuación de Moorse:
2
U (RAB ) = DAB 1 − eαAB (RAB −RABe )
Métodos de Mecánica Molecular 69
7 4
U (RAB ) = DAB αAB − αAB (RAB − RABe ) + αAB (RAB − RABe ) (RAB − RABe )2
2 3 2
12
1 (3) (4)
U (θABC ) = [KABC + KABC (θABC − θABCe ) + KABC (θABC − θABCe )2 ](θABC − θABCe )2
2
donde:
θABC = ángulo entre dos enlaces adyacentes AB y BC.
KABC = constante de fuerza del ángulo de enlace ABC.
N angles
X
Eang = U (θi )
i=1
N
X tors
Etor = U (ωi )
i=1
con:
ωi = ángulo de torsión
{j} es el conjunto de periodicidades.
ϕi es ángulo de fase.
70 Cálculos Computacionales de Estructuras Moleculares
N oop
X kδ
i
Eoop = ∗ d2i
i=1
2
siendo:
kiδ = constante de enlace fuera del plano para el átomo i del tipo δ ((kcal/mol2 .)
Para enlaces fuera del plano de átomos trigonales (ej. átomos sp2 y aromáticos), di
es la altura del átomo central sobre el plano de los sustituyentes.
Energı́a electrostática
qA q B
UAB =
εAB RAB
donde qA es la carga neta sobre el átomo A, y εAB es un parámetro relacionado con la
constante dieléctrica.
NX
atoms NX
atoms
Eele = Uij
i=1 j=i+1
(también de las cargas netas) , junto con sus constantes asociadas, para el conjunto de
átomos que forman nuestro sistema.
Hay muchos campos de fuerza, y cada programa utiliza el suyo propio, incluso a
veces no están muy bien definidos, lo que puede ser un problema para su reproducibi-
lidad.
Algunos de los más conocidos son AMBER, MM2, MM3, MM4, MMFF, DREI-
DING, CHARMM. En la tabla 2.1 de ref.[11], podéis ver una amplia gama de los
existentes, ası́ como de algunos programas que los soportan.
2.3. Amber
2.3.1. Introducción.
Aquı́ disponemos del Amber 16 (2016). Ver [Link] para los cambios
más importantes (Pero ya esta el Amber22 y el AmberTools22).
2.3.2. Descripción
Programas preparatorios:
LEaP es el programa para crear y/o modificar un sistema en Amber.
ParmEd ayuda a obtener información sobre los parámetros definidos en un fi-
chero propio de dichos parámetros.
antechamber es el programa principal para desarrollar campos de fuerza para
moléculas tipo drogas o aminoácidos modificados, usando el GAFF (Campo de
Fuerza del Amber gneral).
Métodos de Mecánica Molecular 73
Programas de simulación:
pmemd (parte del Amber) es una versión del sander optimizada para cálculos
más rápidos y en paralelo.
NAB (Nucleic Acid Builder)es un lenguaje que puede utilizarse para escribir
programas que realizan simulaciones no-periódicas.
Programas de análisis:
# ONIOM(HF/6-31G(d):PM6:UFF)
# ONIOM(B3LYP/6-31G(d,p):UFF) Opt
0 1 0 1 0 1
F -1.041506214819 0.000000000000 -2.126109488809 M
F -2.033681935634 -1.142892069126 -0.412218766901 M
F -2.033681935634 1.142892069126 -0.412218766901 M
C -1.299038105677 0.000000000000 -0.750000000000 M H 5
C 0.000000000000 0.000000000000 0.000000000000 H
H 0.000000000000 0.000000000000 1.100000000000 H
O 1.125833024920 0.000000000000 -0.650000000000 H
Tenemos las capas Alta (H) y Media (M), la H la constituyen los tres últimos
átomos, y el resto forman la M, y se define un átomo de unión, que es el primer átomo
de carbono.
Como decı́amos, hay dos capas en este cálculo, usando el B3LYP para la capa H y
el campo de fuerza UFF del Amber, para la capa M.
Capı́tulo 3
Aplicaciones de la Quı́mica
Computacional
De entre todos ellos, me gusta hacer hincapié en los siguientes programas de cálculo:
GAUSSIAN
Es un paquete de programas que surgió en los 70, (Gaussian70) y se encuentra
en la versión del 2016 (Gaussian16). Aquı́ trabajaremos y utilizaremos la versión
Gaussian09 y/o Gaussian16.
Contiene métodos muy variados:
Mecánica molecular (Amber, UFF, Dreiding )
Métodos Semi-empı́ricos: AM1, PM3, PM6, CNDO, INDO, MINDO/3, MNDO
Métodos monodeterminantales ab initio: RHF, UHF, ROHF
Métodos monodeterminantales del funcional de la densidad: LSD, GGA, MetaG-
GA e hı́bridos.
Métodos tipo capas, con hasta tres capas : ONIOM
Métodos Multiconfiguracionales: CAS-SCF
Métodos de Coupled Clusters: CCD, CCSD, CCSD(T)
Métodos de Interacción de configuraciones: CIS, CISD
Métodos perturbativos: MP2,.. MP5
Métodos que combinan varios cálculos para obtener una buena precisión en los
resultados: G1,G2, G3, CBS..., W1.
Estados excitados con TD para HF y DFT, y EOM en CC
Correcciones relativistas...NO
75
76 Cálculos Computacionales de Estructuras Moleculares
ORCA
ORCA es un flexible, eficiente y fácil de usar herramienta general para quı́mica
cuántica con énfasis especı́fico en propiedades espectroscópicas de moléculas de
capa abierta.
Está constituido por un conjunto de programas controlados desde un driver prin-
cipal.
Es muy interesante para el cálculo de estados excitados y el uso de funciones
multireferenciales.
Es muy útil leerse sus ”consejos y trucos” (Some Tips and Tricks) de su manual.
PSI4
Es un programa ab initio de quı́mica cuántica, diseñado para ser muy eficiente y
proporcionar propiedades moleculares. Tiene la opción de utilizar una interface
Python
Este paquete incluye la posibilidad de cálculos Hartree-Fock, perturbativos (MP2,
... MP(n), coupled cluster (CC2,... CC3), interacción de configuraciones (CISD,
QCISD(T), CI(n), FCI), complete-active-space self-consistent-field, y multi-reference
Aplicaciones de la Quı́mica Computacional. 77
NWChem
Está pensado para trabajar con cualquier tipo de sistema, desde moléculas en fase
gas, a sólidos o sistemas biológicos. Puede utilizar tanto funciones gaussianas,
como ondas planas. Y está muy implementado para utilizarlo con GPUs
Ver sus capacidades: [Link]
CFOUR-MRCC
Programas muy rápidos, y que se usan para realizar cálculos CC.
ABINIT.
Como los anteriores, calcula estructuras electrónicas de sistemas moleculares y de
sólidos dentro de la teorı́a del funcional de la densidad, usando peudopotenciales,
y ondas planas.
78 Cálculos Computacionales de Estructuras Moleculares
Ghemical
Programas de MM y semi-empı́ricos:
Avogadro
Para construir y optimizar con MM estructuras moleculares.
Hulis
Es muy útil para realizar cálculos preliminares de sistemas π-electrónicos
GABEDIT
Interfaz gráfica para los paquetes ed quı́mica computacional: Firefly, Gamess-US,
Gauss, MOLCAS, Molpro, OpenMopac, Orca, MPQC, NWChem y Q-Chem.
4. Cálculo SCF.
Cálculos Post-SCF
Búsqueda de puntos estacionarios en la configuración geométrica.
6. Cálculo de propiedades.
Suelen tener sus propios conjuntos de funciones de base, que, en general, son
comunes, y permiten la utilización de cualquier tipo de función de base, con
las limitaciones en el número cuántico angular propias de cada programa (hoy
en dı́a es posible llegar a utilizar funciones h, aunque su uso no sea habitual).
Dentro de esto se encuentra la utilización de pseudo-potenciales, pudiendo en
algunos utilizar, además de las comúnmente utilizadas Funciones Gaussianas, las
funciones de Slater o funciones de ondas planas para sistemas periódicos.
Algunos incorporan de forma automática los cálculos de counterpoise.
Hay programas que te permiten utilizar, junto a los cuánticos, métodos de Mecáni-
ca Molecular: Amber, Dreiden, UFF.
Cálculos SCF
• El esquema RHF es único, suelen tener el UHF de Pople et al, pero respecto
al ROHF, hay diversidad de posibilidades de escribir la matrices de Fock, ver
por ejemplo en GAMESS: Guest and Saunders, Roothaan, Davidson/1988,
Binkley, Pople and Dobosh, McWeeny and Diercksen, Faegri and Manne. Y
esto implica un proceso de convergencia diferente, aunque si convergen, la
solución final sea la misma.
80 Cálculos Computacionales de Estructuras Moleculares
• Otro cálculo SCF muy común es el DFT. Dado el Babel que existe (Habrı́a
que dejar de hablar de la ”Escalera de Jacob” para empezar a hablar de ”la
Torre de Babel’), ası́, se están incorporando continuamente métodos nuevos,
lo que te exige utilizar la última versión de cada programa si los necesitas o
quieres probarlos. Pero incluso con los funcionales más utilizados y comunes
puede haber diferencias, sobre todo en los procedimientos estándares para
su integración numérica, e incluso en la definición de alguno de los más
populares funcionales (B3LYP).
Otra diferencia es la implementación de los funcionales más modernos.
• MR-SCF: Para acabar con los procesos SCF, la programación de los méto-
dos multireferenciales y sus esquemas, pueden variar de unos programas a
otros, en general utilizan determinantes, pero otros pueden utilizar funcio-
nes de configuración de estado (CSF) que ya son funciones propias de S 2 .
Ası́ mismo, pueden permitir o no utilizar distintas particiones de los OMs,
(CASSCF, RASSCF).
• Junto a estos está el muy utilizado GVB, existiendo un programa (VB2000),
especı́fico de este tipo de cálculos.
Post-SCF
Entramos en el apartado más amplio y donde cada programa es un mundo, los
procesos post-SCF.
• MPn: Comenzaremos por los más sencillos, los perturbativos, los habituales
son los perturbativos tipo Møller Pesset, comenzando con el MP2 y llegando
al MP4, MP5 en los programas habituales, al MP(n) que parece que tiene
el Psi. También hay diferencias en cuanto a que los tengan programados
para distintas funciones SCF (RHF, UHF, ROHF, CASSCF-MP2, MR-MP2,
CCSD-MP2...)
• CI: Respecto a los CI, comienzan en el CIS, llegando en algunos casos al
FullCI (Dalton, y aunque el CCSD sea común a la mayorı́a, otros como
QCISD(TQ) son más propios de ciertos programas.
• CC: Otro gran grupo son los CCs, de nuevo son comunes los CCSD, pero hay
programas más especı́ficos que llegan a CCSDTQP, y la serie CC2,...,CC5,
CC(n), como el CFOUR-MRCC. (Unos usan la descomposición de Cholesky,
..) (Ver J Chem. Theor. Comput.-9-2687)
• El estudio de estados excitados se puede llevar a cabo con los métodos MR-
SCF, los CI y con EOM-CC, pero además son muy actuales los TD-DFT,
• Indicar que se pueden realizar cálculos de los denominados ”Procedimientos
para la obtención de energı́a con gran precisión”, como G1, G2, G3, G4 ,
..., CBS-4, CBS-q, CBS-QB3, ROCBS-QB3, CBS-Q, CBS-APNO, y W1U,
W1BD, W1RO (todos ellos en el Gussian).
• Solvatación con diferentes métodos (PCM, EFP, SCRF,....)
• QM/MM, con dos o tres capas. (ONION)
• Correcciones relativistas (ZORA)
Aplicaciones de la Quı́mica Computacional. 81
Propiedades moleculares
• Mössbauer-parameters
• Polarizabilidades e hiperpolarizabilidades
Como se ve hay muchos programas que permiten hacer cálculos moleculares mecano-
cuánticos, aunque a lo largo de este curso usaremos con más asiduidad sólo algunos de
ellos: Gaussian[22, 23, 24, 25, 26], Gamess[27] NWChem[28] y ORCA[29]
Todos ellos nos permiten hacer cálculos de estructuras moleculares con métodos ab
inito y semi-empı́ricos. En general nos referiremos al Gaussian, aunque en general las
ideas son validas para el resto:
En todos ellos podéis ver sus capacidades, por ejemplo, el Gamess, tiene este cuadro:
82 Cálculos Computacionales de Estructuras Moleculares
Dentro de los semi-empı́ricos tenemos los modelos CNDO, INDO, MINDO3, MN-
DO, ZINDO/S AM1 y PM3. Los cuatro últimos usan el conjunto modificado de paráme-
tros del MOPAC. Con MINDO, MNDO, AM1, y PM3 se pueden hacer CIs limitadas.
Otro método clásico es el Hückel y extended-Hückel con diversos parámetros.
También dispone de cálculos de Mecánica Molecular, con varios métodos: AMBER,
DREIDING y UFF.
Métodos post-SCF:
MP2, MP3, MP4, MP5. Cálculos perturbativos tipo Moller-Plesset, con per-
turbaciones de segundo (MP2), tercero (MP3), cuarto (MP4) y quinto or-
den(MP5). El MP4 se puede realizar considerando las SDQ o las SDTQ
excitaciones.
CI Método de Interacción de Configuraciones, con SD , que puede ampliarse a
SDT o SDTQ excitaciones.
CCD Método de ”coupled clusters”, que puede ser con dobles excitaciones, o
considerar las simples también (CCSD). (CCSD(T))
84 Cálculos Computacionales de Estructuras Moleculares
LYP Lee, Yang y Parr (C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785
(1988), ver también B. Miehlich, A. Savin, H. Stoll, and H. Preuss, Chem.
Phys. Lett. 157, 200 (1989)).
P86 Correcciones a su funcional local, por Perdew 1986 (J. P. Perdew, Phys.
Rev. B 33, 8822 (1986).
PW91 Perdew 1991 (Ver el de intercambio)
B95 Definido como parte de su funcional hı́brido parametrizado , con correccio-
nes ’time-dependent’ (A. D. Becke, J. Chem. Phys. 104, 1040 (1996)).
PBE El funcional con corrección de gradientes de Perdew, Burke and Ernzerhof.
Aplicaciones de la Quı́mica Computacional. 85
Funcionales hı́bridos:
Ee = Ecombinada + ∆(HLC)
E corr,X = E corrCBS + cX −3
La de Peterson y Woon[34, 35], con tres parámetros lineales:
2
E X = E CBS + Ae−(X−1) + Be−(X−1)
Métodos Semi-empı́ricos
Los métodos que permite calcular son:
HUCKEL
ZINDO/S
AM1
PM3
PM6
Efectos de solvente:
Propiedades Moleculares
Frecuencias vibracionales y modos normales
Polarizabilidades e hyperpolarizabilidades
Análisis de población
Densidades de espı́n
Espectro Raman
Termoquı́mica en Gaussian
90 Cálculos Computacionales de Estructuras Moleculares
3.3. Ejemplos:
Veamos un ejemplo concreto, un cálculo de optimización de geometrı́a de la molécula
de agua, con una base cc-PVDZ, mediante un cálculo hı́brido DFT como el B3LYP:
3.3.1. Gaussian09-16:
%mem=1Gb
%nproc=3
%Chk=[Link]
#P rB3LYP/cc-pvdz 6D 10F Opt=() SCF=() iop(5/33=1) Pop=(Full)
h2o-G03
3.3.2. Gamess:
123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789
3.3.3. NWChem:
scratch_dir /home/scr
permanent_dir /home/sanfa/h2o-NWC
Title "h2o-NWC"
Start h2o-NWC
echo
charge 0
ecce_print /home/sanfa/h2o-NWC/[Link]
dft
mult 1
XC b3lyp
grid fine
mulliken
end
driver
default
end
3.3.4. ORCA:
! B3LYP cc-pvDZ
% pal nprocs 8
end
* xyz 0 1
O 0.00000 0.00000 0.00000
H 0.922641 0.652406 0.00000
H -0.922641 0.652406 0.00000
*
%mem=1200MW
%nproc=3
# B3lyp/aug-cc-pVDZ opt pop=regular gfprint
0 1
F 0.0 0.0 0.0
F 0.0 0.0 0.7
end_file
scratch_dir /home/scr
permanent_dir /home/sanfa/pru-pro
Title "Molecula F2 optimizacion"
start f2_opt
charge 0
#ecce_print /home/sanfa/pru-prog/[Link]
dft
mult 1
# iterations 50
# print kinetic_energy
xc b3lyp
# mulliken
# decomp
end
memory 750 mb
molecule F2 {
0 1
F
F 1 0.7
}
optimize (’b3lyp’)
#
# Molecula de Fluor calculada con Orca
#
%pal nprocs 8 # any number (positive integer)
end
R=0.70000
Aplicaciones de la Quı́mica Computacional. 95
*ACES2(CALC=CCSD(T),BASIS=AUG-PVDZ
GEO_CONV=10
CC_CONV=12
LINEQ_CONV=12
SCF_CONV=12
MEMORY=600000000)
# Molecula F2
basis=aug-cc-pVDZ
calc=CC4
mem=100GB
gopt=full
core=frozen
optalg=simplex
geom
F
F 1 R
R = 0.7
BASIS
aug-cc-pVDZ
F2 molecule. Basis: aug-cc-pVDZ.
Geometria equilibrio
Atomtypes=1 Generators=1 Y X Angstrom
Charge=9.0 Atoms=2
F 0.000000 0.000000 0.000000
F 0.000000 0.000000 0.700000
**DALTON INPUT
.OPTIMIZE
**WAVE FUNCTION
.HF
*END OF INPUT
&control
calculation = ’scf’ ,
tstress = .false.
/
&system
ibrav = 1,
celldm(1) = 12.0,
nat = 2,
ntyp = 1,
ecutwfc = 80,
ecutfock=160,
input_dft = ’B3LYP’
nspin = 1
exxdiv_treatment = ’gygi-baldereschi’
x_gamma_extrapolation = .TRUE.
/
&electrons
conv_thr = 1.0d-7
/
ATOMIC_SPECIES
F 9.00d0 [Link]
ATOMIC_POSITIONS (bohr)
F 0.0000 0.0000 0.0000
F 0.0000 0.0000 1.33823
K_POINTS gamma
......................
******************************************
Gaussian 16: ES64L-G16RevC.01 3-Jul-2019
7-May-2021
******************************************
%nproc=3
-------------------------------------------
# B3lyp/aug-cc-pVDZ opt pop=regular gfprint
-------------------------------------------
1/18=20,19=15,26=3,38=1/1,3;
Aplicaciones de la Quı́mica Computacional. 97
2/9=110,12=2,17=6,18=5,40=1/2;
3/5=16,7=10,11=2,24=100,25=1,30=1,71=1,74=-5/1,2,3;
4//1;
5/5=2,38=5/2;
6/28=1/1;
7//1,2,3,16;
1/18=20,19=15,26=3/3(2);
2/9=110/2;
99//99;
2/9=110/2;
3/5=16,7=10,11=2,25=1,30=1,71=1,74=-5/1,2,3;
4/5=5,16=3,69=1/1;
5/5=2,38=5/2;
7//1,2,3,16;
1/18=20,19=15,26=3/3(-5);
2/9=110/2;
6/19=2,28=1/1;
99/9=1/99;
--------------------------------
title "Molecula F2 optimizacion"
--------------------------------
Symbolic Z-matrix:
Charge = 0 Multiplicity = 1
F 0. 0. 0.
F 0. 0. 0.7
GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad
Berny optimization.
Initialization pass.
----------------------------
! Initial Parameters !
! (Angstroms and Degrees) !
--------------------- ----------------------
! Name Definition Value Derivative Info. !
-----------------------------------------------------------------------
! R1 R(1,2) 0.7 estimate D2E/DX2 !
-----------------------------------------------------------------------
Input orientation:
---------------------------------------------------------------------
Center Atomic Atomic Coordinates (Angstroms)
Number Number Type X Y Z
98 Cálculos Computacionales de Estructuras Moleculares
---------------------------------------------------------------------
1 9 0 0.000000 0.000000 0.000000
2 9 0 0.000000 0.000000 0.700000
---------------------------------------------------------------------
Stoichiometry F2
Framework group D*H[C*(F.F)]
Deg. of freedom 1
Full point group D*H NOp 8
Largest Abelian subgroup D2H NOp 8
Largest concise Abelian subgroup C2 NOp 2
Standard orientation:
---------------------------------------------------------------------
Center Atomic Atomic Coordinates (Angstroms)
Number Number Type X Y Z
---------------------------------------------------------------------
1 9 0 0.000000 0.000000 0.350000
2 9 0 0.000000 0.000000 -0.350000
---------------------------------------------------------------------
Rotational constants (GHZ): 0.0000000 108.5760439 108.5760439
Standard basis: Aug-CC-pVDZ (5D, 7F)
AO basis set (Overlap normalization):
Atom F1 Shell 1 S 7 bf 1 - 1 0.0000000 0.0000000 0.661404143598
0.1471000000D+05 0.7229535153D-03
0.2207000000D+04 0.5569055564D-02
0.5028000000D+03 0.2834429748D-01
0.1426000000D+03 0.1067956983D+00
0.4647000000D+02 0.2878097307D+00
0.1670000000D+02 0.4517054881D+00
0.6356000000D+01 0.2668829077D+00
Atom F1 Shell 2 S 7 bf 2 - 2 0.0000000 0.0000000 0.661404143598
0.1471000000D+05 0.9329717475D-05
0.5028000000D+03 0.3153039638D-03
0.1426000000D+03 -0.3125687006D-02
0.4647000000D+02 -0.1184270573D-01
0.1670000000D+02 -0.1257376908D+00
0.6356000000D+01 -0.9650219096D-01
0.1316000000D+01 0.1094036315D+01
Atom F1 Shell 3 S 1 bf 3 - 3 0.0000000 0.0000000 0.661404143598
0.3897000000D+00 0.1000000000D+01
Atom F1 Shell 4 S 1 bf 4 - 4 0.0000000 0.0000000 0.661404143598
0.9863000000D-01 0.1000000000D+01
Atom F1 Shell 5 P 3 bf 5 - 7 0.0000000 0.0000000 0.661404143598
0.2267000000D+02 0.6483402149D-01
0.4977000000D+01 0.3405353598D+00
0.1347000000D+01 0.7346464068D+00
Atom F1 Shell 6 P 1 bf 8 - 10 0.0000000 0.0000000 0.661404143598
0.3471000000D+00 0.1000000000D+01
Aplicaciones de la Quı́mica Computacional. 99
...................................
Harris functional with IExCor= 402 and IRadAn= 5 diagonalized for initial gues
**********************************************************************
**********************************************************************
Orbital symmetries:
Occupied (SGG) (SGU) (SGG) (PIU) (PIU) (SGU) (SGG) (PIG)
(PIG)
Virtual (SGU) (SGG) (PIU) (PIU) (SGG) (PIG) (PIG) (SGU)
(DLTG) (DLTG) (SGG) (PIU) (PIU) (SGU) (PIG) (PIG)
(SGG) (SGU) (DLTU) (DLTU) (PIU) (PIU) (PIG) (PIG)
(SGU) (SGG) (SGU) (DLTG) (DLTG) (PIU) (PIU) (DLTU)
(DLTU) (SGG) (PIG) (PIG) (SGU)
The electronic state is 1-SGG.
Alpha occ. eigenvalues -- -24.93431 -24.91790 -2.42203 -1.11794 -1.11794
Alpha occ. eigenvalues -- -0.91677 -0.85952 -0.12829 -0.12829
Alpha virt. eigenvalues -- 0.07746 0.08506 0.09337 0.09337 0.19954
Alpha virt. eigenvalues -- 0.21531 0.21531 0.39058 0.83710 0.83710
Alpha virt. eigenvalues -- 0.94045 0.99334 0.99334 1.05717 1.21137
Alpha virt. eigenvalues -- 1.21137 1.22955 1.39192 1.43247 1.43247
Alpha virt. eigenvalues -- 1.49810 1.49810 2.22510 2.22510 2.54377
100 Cálculos Computacionales de Estructuras Moleculares
..................................
Density Matrix:
-------------------------------------------------------------------
Cartesian Forces: Max 7.731867899 RMS 4.463996013
GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad
Berny optimization.
FormGI is forming the generalized inverse of G from B-inverse, IUseBI=4.
Internal Forces: Max 7.731867899 RMS 7.731867899
Search for a local minimum.
Step number 1 out of a maximum of 20
All quantities printed in internal units (Hartrees-Bohrs-Radians)
Mixed Optimization -- RFO/linear search
Second derivative matrix not updated -- first step.
The second derivative matrix:
R1
R1 17.34000
ITU= 0
Eigenvalues --- 17.34000
RFO step: Lambda=-2.94682750D+00 EMin= 1.73400000D+01
Linear search not attempted -- first point.
Maximum step size ( 0.300) exceeded in Quadratic search.
-- Step size scaled by 0.787
Iteration 1 RMS(Cart)= 0.14142136 RMS(Int)= 0.10000000
Iteration 2 RMS(Cart)= 0.07071068 RMS(Int)= 0.00000000
Iteration 3 RMS(Cart)= 0.00000000 RMS(Int)= 0.00000000
ClnCor: largest displacement from symmetrization is 5.63D-17 for atom 1.
Variable Old X -DE/DX Delta X Delta X Delta X New X
(Linear) (Quad) (Total)
R1 1.32281 7.73187 0.00000 0.30000 0.30000 1.62281
Item Value Threshold Converged?
Maximum Force 7.731868 0.000450 NO
RMS Force 7.731868 0.000300 NO
Maximum Displacement 0.150000 0.001800 NO
RMS Displacement 0.212132 0.001200 NO
Predicted change in Energy=-1.539260D+00
GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad
..............................................
-------------------------------------------------------------------
Center Atomic Forces (Hartrees/Bohr)
Number Number X Y Z
-------------------------------------------------------------------
1 9 0.000000000 0.000000000 -0.000008754
2 9 -0.000000000 -0.000000000 0.000008754
-------------------------------------------------------------------
Cartesian Forces: Max 0.000008754 RMS 0.000005054
GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad
Berny optimization.
Using GEDIIS/GDIIS optimizer.
Internal Forces: Max 0.000008754 RMS 0.000008754
Search for a local minimum.
Step number 6 out of a maximum of 20
All quantities printed in internal units (Hartrees-Bohrs-Radians)
Mixed Optimization -- En-DIIS/RFO-DIIS
Swapping is turned off.
Update second derivatives using D2CorX and points 5 6
DE= -3.12D-07 DEPred=-4.34D-07 R= 7.18D-01
Trust test= 7.18D-01 RLast= 1.31D-03 DXMaxT set to 4.24D-01
The second derivative matrix:
R1
R1 0.37695
ITU= 0 1
Eigenvalues --- 0.37695
En-DIIS/RFO-DIIS/Sim-DIIS IScMMF= -3 using points: 6 5
RFO step: Lambda=-2.03310341D-10.
DidBck=F Rises=F RFO-DIIS coefs: 0.98225 0.01775
Iteration 1 RMS(Cart)= 0.00001642 RMS(Int)= 0.00000000
Iteration 2 RMS(Cart)= 0.00000000 RMS(Int)= 0.00000000
ClnCor: largest displacement from symmetrization is 7.11D-22 for atom 2.
Variable Old X -DE/DX Delta X Delta X Delta X New X
(Linear) (Quad) (Total)
R1 2.65087 0.00001 0.00002 0.00000 0.00002 2.65090
Item Value Threshold Converged?
Maximum Force 0.000009 0.000450 YES
RMS Force 0.000009 0.000300 YES
Maximum Displacement 0.000012 0.001800 YES
RMS Displacement 0.000016 0.001200 YES
Predicted change in Energy=-1.016552D-10
Optimization completed.
-- Stationary point found.
----------------------------
! Optimized Parameters !
! (Angstroms and Degrees) !
-------------------------- --------------------------
104 Cálculos Computacionales de Estructuras Moleculares
Input orientation:
---------------------------------------------------------------------
Center Atomic Atomic Coordinates (Angstroms)
Number Number Type X Y Z
---------------------------------------------------------------------
1 9 0 0.000000 0.000000 -0.351391
2 9 0 0.000000 0.000000 1.051391
---------------------------------------------------------------------
Stoichiometry F2
Framework group D*H[C*(F.F)]
Deg. of freedom 1
Full point group D*H NOp 8
Largest Abelian subgroup D2H NOp 8
Largest concise Abelian subgroup C2 NOp 2
Standard orientation:
---------------------------------------------------------------------
Center Atomic Atomic Coordinates (Angstroms)
Number Number Type X Y Z
---------------------------------------------------------------------
1 9 0 0.000000 0.000000 0.701391
2 9 0 0.000000 0.000000 -0.701391
---------------------------------------------------------------------
Rotational constants (GHZ): 0.0000000 27.0364328 27.03643
**********************************************************************
.................................................
1\1\GINC-CUANT5\FOpt\RB3LYP\Aug-CC-pVDZ\F2\SANFA\07-May-2021\0\\# B3ly
p/aug-cc-pVDZ opt pop=regular gfprint\\title "Molecula F2 optimizacion
"\\0,1\F,0.,0.,-0.3513912699\F,0.,0.,1.0513912699\\Version=ES64L-G16Re
vC.01\State=1-SGG\HF=-199.5353652\RMSD=4.143e-10\RMSF=5.054e-06\Dipole
=0.,0.,0.\Quadrupole=-0.2348554,-0.2348554,0.4697109,0.,0.,0.\PG=D*H [
C*(F1.F1)]\\@
The archive entry for this job was punched.
Herramientas gráficas
No son fundamentales para el análisis de los resultados, pero lo hace más agradable,
y además estamos en el tiempo de las imágenes.
Ayudan a ver la geometrı́a del sistema, su evolución en una optimización o un
camino de reacción, la forma de su densidad electrónica, su volumen, los orbitales
moleculares y, si se han calculado, los potenciales electrostáticos y modos vibracionales.
En las páginas web
Molecular Graphics Applications in the CSB y Universal Molecular Modeling Soft-
ware List, tenemos una lista bastante amplia del software gráfico (y otros) que existe.
Otro ejemplo de página es: [Link] ( Department of Chemistry del
Imperial College London.)
A veces interesa hacer una ”pelı́cula”, puedes tomar diversas imágenes y juntarlas
con utilidades como mencoder:
En Linux con un simple grupo de comandos desde la consola puedes crear archivos
de vı́deo a partir de imágenes jpg que se encuentren en un determinado directorio.
107
108 Cálculos Computacionales de Estructuras Moleculares
La mayorı́a tiene versiones para Windows y Linux, y para algunos hay que comprar
la licencia de uso.
También se pueden incluir aquı́ otros paquetes ya vistos, con interfaces gráficas,
como Ecce.
4.1.1. Molden
4.1.2. Jmol
Como su página web dice: Un visor Java de código abierto para estructuras quı́micas
en tres dimensiones. con prestaciones para compuestos quı́micos, cristales, materiales
y biomoléculas.
Lee fácilmente de Gaussian, GAMESS, NWChem, Amber y VASP. Las densidades
y Orbitales Moleculares se representan fácilmente
4.1.3. Mercury
The Cambridge Crystallographic Data Centre (CCDC), muy agradable y potente,
sobre todo para estructuras cristalinas.
110 Cálculos Computacionales de Estructuras Moleculares
4.1.4. Molekel
Paquete gráfico molecular para la visualización de resultados de estructura atómica
y molecular obtenidos con GAUSSIAN, GAMESS-US, ADF y otros programas.
Muy sencillo de instalar en Linux, y con versión para Windows.
Con el Gaussian03, tenı́a un problema, que se resuelve simplemente con editar el
fichero de salida del Gaussian y cambiar Gaussian 03 por Gaussian 98.
Con manual-tutorial incluido.
4.1.5. RasMol
Este visor de estructuras moleculares lee ficheros tipo ”pdb” y ”xyz”.
Y toda la información sobre él la podéis encontrar en :
RasMol Home Page Contents
Manual
Es el origen de Jmol y Chime.
Herramientas Gráficas. 111
4.1.6. MGLTools
Este paquete contiene varios programas, por ejemplo:
vision
4.2.1. Avogadro
4.2.2. Gabedit
Es un paquete gráfico, que actualmente viene incluido y mantenido por la comu-
nidad Linux, y actúa como interface para ejecutar los programas Gaussian, Gamess,
Orca, Molcas, Molpro, Mopac, MPQC, PCGamess, Tinker.
Ver el manual.
Herramientas Gráficas. 113
4.2.3. Ghemical
4.2.4. Ecce
114 Cálculos Computacionales de Estructuras Moleculares
4.2.5. Hulis
Es muy útil para realizar cálculos preliminares de sistemas π-electrónicos.
4.3.2. gnuplot
Programa para representaciones bi- y tri-dimensionales.
Suele venir con las distribuciones de linux.
Y tiene buena documentación.
Ejercicios prácticos:
rg16
g16 [Link] & (El resultado en [Link])
%mem=1GB
%nproc=3
#P B3LYP/3-21G SCF=Tight
gfinput iop(5/33=1,6/7=3)
Carga Multiplicidad
F
--link1--
Otro calculo
117
118 Cálculos Computacionales de Estructuras Moleculares
CBSTcc−pV
P
XZ
aug-cc-pVDZ
aug-cc-pVTZ
aug-cc-pVQZ
CBSEE
aug−cc−pV XZ
CBSTaug−cc−pV
P
XZ [Link] [Link] [Link]
Tabla 5.1: Energı́as (E) en hartrees para el átomo de Fluor, con distintos conjuntos
de funciones de base y extrapolados: Extrapolación exponencial (EE) y Tres pa-
rametros con exponentes enteros (TP)
Cálculos:
#P method=(FULL)/cc-pVDZ
Cuestiones:
%mem=1GB
%nproc=3
# b3lyp/cc-pVTZ scan guess=(mix)
0 1
F
F 1 r
variables:
r=0.9 41 -0.1
end_file
# Titulo y comentarios
x1 y1 y’1 y’’1 ...
x2 y2 y’2 y’’2 ...
Al final de un cálculo con ”scan”, podemos encontrar los datos que nos
ayuden.
1
Realizar un barrido desde 0.9 Å hasta 5.0 Å, para los cuatro tipos de cálculo, y al revés, desde
5.0 hasta 0.9, para los cálculos ”unrestricted”.
Ejercicios. 121
Prácticas avanzadas:
a) Con ayuda del programa ”Gabedit” y utilizado el programa Gaussian[25,
26] o el Orca[29] (porque esto es de un tutorial suyo) , realizar la siguiente
práctica:
Obtener la configuración nuclear óptima de una de las siguientes molécu-
las:
C2 H4 (Etileno) , C2 H4 O (Acetaldehido), C2 H4 O2 (Acid. acético) o C2 H2 Cl2
(Etileno 1,2 dicloro - cis - trans).
Hacerlos a nivel CCSD/6-31+G* y M06-2X/6-31+G*2 .
Para dicha configuración de equilibrio:
• Obtener el espectro IR y compararlo con el experimental (gas de
mayor resolución)
• Representar gráficamente los OMs HOMO y LUMO)3 .
• Compara los resultados con los obtenidos en el curso de DFT.
NIST-chemistry:[Link]
b) Para la molécula elegida previamente, calcular el potencial de ionización
vertical.
Para ello utilizar diversos métodos: MPn (n=2,,..5)4 . Para el catión realizar
los cálculos ’unrestricted’ y utilizar el programa Gaussian.
c) Formar un dı́mero del sistema elegido, y analizar la energı́a de interacción,
considerando el error de superposición de bases, y utilizando los métodos:
RHF, MP2, CISD, CCSD.
Realizarlo con el programa NwChem[39]. Y si funciona, con Ecce.
Se puede hacer un análisis de las interacciones débiles tipo NCI (Noncovalent
Interactions). (Ver Tutorial nci-analysis).
d ) Seguir un tutorial del Amber.
Por ejemplo: An Introduction to Molecular Dynamics Simulations using AM-
BER, o Simulating a DNA polyA-polyT Decamer:.
e) Práctica con algún sistema de interés para el alumno:
De un aminoácido, como la alanina(CH3 − CH(N H2 )COOH) a nivel
RHF, con la base 3-21G y con los métodos semiempı́ricos AM1 y PM6.
La geometrı́a inicial se puede obtener con ”ghemical”
Para el sistema cuya geometrı́a inicial se ha construido con ”ghemical”
o ”gabedit”. Primero con un método semi-empı́rico (AM1, PM3, PM6),
y después a nivel HF y B3LYP con la base 3-21G.
Cálculo de algún sistema en que estéis interesados, a un nivel aceptable
para su tamaño.
2
Asegurase de que considera la simetrı́a del sistema
3
Recordad mandar escribir la funciones de base y orbitales en el .log, usando, por ejemplo: gfinput
pop=full iop(6/6=3)
4
Basta con realizar el cálculo MP5
122 Cálculos Computacionales de Estructuras Moleculares
Bibliografı́a
[1] G. Li Manni et al., Journal of Chemical Theory and Computation 10, 3669
(2014), pMID: 26588512.
[3] F. Jensen, Introduction to Computational Chemistry (John Wiley & Sons, West
Sussex, 2002).
[4] A. Szabo y N. S. Ostlund, Modern Quantum Chemistry (Dover, New York, 1996).
[5] R. M. Martin, Electronic Structure: Basic Theory and Practical Methods (Cam-
bridge University Press, Cambridge, 2005).
[13] G. Herzberg, Molecular Spectra and Molecular Structure. III. Electronic Spectra
and Electronic Structure of Polyatomic Molecules (D. Van Nostrand, Princeton,
1966).
123
124 Cálculos Computacionales de Estructuras Moleculares
[22] M. J. Frisch et al., Gaussian94, Revision D.4 (Gaussian Inc., Pittsburgh PA,
1995).
[23] M. J. Frisch et al., Gaussian98, Revision A.7 (Gaussian Inc., Pittsburgh PA,
1998).
[24] M. J. Frisch et al., Gaussian 03, Revision C.02, Gaussian, Inc., Wallingford, CT,
2004.
[26] M. J. Frisch et al., Gaussian 16 Revision C.01, 2016, gaussian Inc. Wallingford
CT.
[27] G. M. J. Barca et al., The Journal of Chemical Physics 152, 154102 (2020).
[35] D. E. Woon y T. H. Dunning, The Journal of Chemical Physics 101, 8877 (1994).
[50] P. S. Svendsen y U. von Barth, Physical Review B (Condensed Matter) 54, 17402
(1996).
[54] J. P. Perdew, K. Burke, y M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
[55] J. P. Perdew, K. Burke, y M. Ernzerhof, Phys. Rev. Lett. 78, 1396 (1997).
[59] A. D. Boese y N. C. Handy, The Journal of Chemical Physics 114, 5497 (2001).
[63] W. J. Carr, Jr, R. A. Coldwell-Horsfall, y A. E. Fein, Phys. Rev. 124, 747 (1961).
[65] P. Gombas, Die Statistiche Theorie des Atoms und ihre Anwendungen (Springer
Verlag, Viena, 1949).
[76] H. Stoll, C. M. E. Pavlidou, y H. Preuss, Theor. Chim. Acta 49, 143 (1978).
[79] A. Savin, H. Stoll, y H. Preuss, Theor. Chim. Acta 70, 407 (1986).
[109] B. Miehlich, A. Savin, H. Stoll, y H. Preuss, Chem. Phys. Lett. 157, 200 (1989).
[117] A. Savin, Int. J. Quantum Chem.: Quantum Chem. Symp. 22, 59 (1988).
[130] R. Colle, F. Moscardo, P. Riani, y O. Salvetti, Theor. Chim. Acta 44, 1 (1977).
[132] F. Moscardó, M. Paniagua, y E. San-Fabián, Theor. Chim. Acta 53, 377 (1979).
[136] R. McWeeny, en The New World of Quantum Chemistry, editado por B. Pullman
y R. G. Parr (Reidel, Dordrecht, 1977), pp. 3–31.
[147] C. Adamo y V. Barone, The Journal of Chemical Physics 110, 6158 (1999).
[150] J.-D. Chai y M. Head-Gordon, Phys. Chem. Chem. Phys. 10, 6615 (2008).
[151] J.-D. Chai y M. Head-Gordon, The Journal of Chemical Physics 128, 084106
(2008).
Funcionales DFT
El problema de la teorı́a del Funcional densidad es que para obtener las energı́as
y densidades del estado fundamental, precisamos conocer el funcional de intercambio-
correlación exacto.
Dado que no se conoce, se han de introducir aproximaciones, y se formulan funciona-
les aproximados, que, históricamente, se han desglosado en funcionales de intercambio
junto a otros funcionales de correlación.
EX [ρα , ρβ ] EC [ρα , ρβ ]
Entre los funcionales de intercambio vamos a diferenciar los locales y aquellos otros
que incluyen correcciones del gradiente de la densidad.
Para los funcionales de energı́a de correlación, distinguiremos a su vez los locales, con
sin corrección de auto-correlación, aquellos que incluyen corrección de gradiente, y por
último los funcionales que además de depender de la densidad, dependen de las matrices
densidad reducidas de primer o de segundo orden1 . Se dice que una aproximación al
funcional de energı́a de intercambio-correlación es local cuando viene dado por una
expresión de la forma
Z
EXC [ρα , ρβ ] = ρ(r)ϵXC (ρα , ρβ ) dr, (A.1)
∂(ρϵXC )
µXC (r) = , (A.2)
∂ρα
para electrones α, y
∂(ρϵXC )
µ̄XC (r) = , (A.3)
∂ρβ
1
No todos los funcionales pueden clasificarse dentro de este esquema. Como ejemplo, tenemos los
funcionales deducidos a partir del concepto de hueco de intercambio y correlación, como la “average-
density approximation” de Gunnarsson, Jonson y Lundqvist [40], el “weighted-density method” de
Gunnarsson, Jonson y Lundqvist [41], y de Alonso y Girifalco [42, 43], o el “modified weighted-density
scheme” de Gunnarsson y Jones [44]. Para una revisión sobre el tema ver [45, sección 8.5].
131
132 Cálculos Computacionales de Estructuras Moleculares
para electrones β.
Un caso particular muy importante de aproximación local es la aproximación local
gas de electrones [17, 46], caso en que la función ϵXC (ρα , ρβ ) es la energı́a de intercambio-
correlación por electrón de un gas de electrones uniforme con densidades de espı́n ρα
y ρβ . Es la popular aproximación Local Spin Density o LSD.
La aproximación local gas de electrones (tanto para el intercambio como para la
correlación) fue, históricamente, la primera que se utilizó [47, 17, 48, 46], y continua
siendo una de las más empleadas actualmente.
que es la distancia media entre un par de electrones, dice que el potencial es propor-
cional a r1s (Toma el modelo de una esfera de radio rs cargada uniformemente), o lo
que es lo mismo, proporcional a ρ1/3 (r) y el intercambio se describe como:
Z
Ex [ρ(r)] ≃ Cx ρ(r)ρ1/3 (r(dr
LDA
XZ x2σ
EX = EX − β ρσ4/3 −1 x
dr (A.7)
σ
1 + 6βx σ sinh σ
134 Cálculos Computacionales de Estructuras Moleculares
1/3
3 3
CX = (A.9)
2 4π
|∇ρσ |
xσ = 4/3
(A.10)
ρσ
y la constante β se ha ajustado empı́ricamente a un valor de 0.0042 u.a.
P W 91 1 X P W 91
EX [ρα , ρβ ] = E [ρσ ] (A.11)
2 σ X
donde Z
P W 91
EX [ρ] = ρεx (rs , 0)F (s)dr (A.12)
3 |∇ρ|
εx (rs , 0) = − kF s= (A.13)
4π 2kF ρ
1/3
kF = (3π 2 ρ) ,y
2
1 + 0.19645s · sinh−1 (7.7956s) + (0.2743 − 0.1508e−100s )s2
F (s) = (A.14)
1 + 0.19645s · sinh−1 (7.7956s) + 0.004s4
0.804
F (s) = 1 + 0.804 − . (A.16)
1 + 0.21951
0.804
s2
1/3
donde s = |∇ρ| / (2ρkF ), con kF = (3π 2 ρ) .
Ejercicios. 135
2
1 + 0.21516s arcsinh(7.7956s) + 0.30042 − 0.17696e−100s s2
F mPW (s) = . (A.18)
1 + 0.21516s arcsinh(7.7956s) + 0.00228s3.73
2
s2 s2
F (s) = 1 + 0.21951 − 0.015 . (A.20)
1 + 0.157s2 1 + 0.157s2
[58]
Z Z
ExTPSS [ρ] = fxTPSS 3
(ρ(r), |∇ρ(r)| , t(r)) d r = −Cx ρ4/3 (r)F TPSS (p(r), z(r)) d3 r,
(A.21)
1/3
where Cx = (3/4) (3/π) and
0.804
F TPSS (p, z) = 1 + 0.804 − x . (A.22)
1 + 0.804
136 Cálculos Computacionales de Estructuras Moleculares
In Eq. (A.22),
s
2
2
1 10 z 146 2 73 1 3 p2
x = √ 2 + c p + q̃ − q̃ b z +
(1 + ep) 81 (1 + z 2 )2 2025 b 405 2 5 2
2 2 #
1 10 √ 10 3
+ p2 + 2 e z + 0.21951ep3 , (A.23)
0.804 81 81 5
PN
where c = 1.59096, e = 1.537, z = tW /t with t = (1/2) i=1 ∇ψi∗ · ∇ψi and tW =
(1/8) |∇ρ|2 /ρ,
|∇ρ|2
p= , (A.24)
4 (3π 2 )2/3 ρ8/3
and
9 α−1 2
q̃b = p + p, (A.25)
20 1 + 0.4α (α − 1) 3
2/3
where α = t − tW /tTF with tTF = (3/10) (3π 2 ) ρ5/3 .
donde
VWN P αC (x) F P αC (x) 4
ϵC (ρα , ρβ ) = ϵC (x) + ′′ f (ξ) + ϵC (x) − ϵC (x) − ′′ ξ f (ξ), (A.27)
f (0) f (0)
con
ρα − ρβ
ξ = (A.28)
ρα + ρβ
(1 + ξ)4/3 + (1 − ξ)4/3 − 2
f (ξ) = (A.29)
2(21/3 − 1)
4
f ′′ (0) = (A.30)
9(21/3 − 1)
s
1/3
√ 1 3
x = rs = , (A.31)
a0 4π(ρα + ρβ )
siendo a0 el radio de Bohr.
(Notar que : a30 rs3 4π
3
= ρ1 )
P F
Las funciones ϵC (x), ϵC (x) y αC (x), están todas ellas expresadas por la siguiente
fórmula
x2 (x − x0 )2 2b(c − x20 )
bx0 Q
A ln − ln + arctan , (A.32)
X(x) X(x0 ) X(x) X(x0 )Q 2x + b
con
X(x) = x2 + bx + c (A.33)
√
Q = 4c − b2 . (A.34)
138 Cálculos Computacionales de Estructuras Moleculares
Los parámetros A, x0 , b y c tienen valores diferentes para cada una de las tres funciones
P F
ϵC (x), ϵC (x) y αC (x), como se especifica en la tabla A.1
A x0 b c
P
ϵC (x) 0.0310907 -0.10498 3.72744 12.9352
F
ϵC (x) 0.0155453 -0.32500 7.06042 18.0578
1
αC (x) − 6π -0.0047584 1.13107 13.0045
P F
Tabla A.1: Parámetros para las funciones ϵC y ϵC que representan respectivamente
las energı́as de correlación por electrón para los casos paramagnético y ferromagnético,
mientras que αC es la espı́n stiffness. El resultado obtenido con estos parámetros está
en unidades atómicas.
Con el funcional de Vosko, Wilk y Nusair se obtendrán los mejores resultados para
gases de electrones, pero no necesariamente que se obtendrán también los mejores
resultados para átomos y moléculas. Se pueden obtener resultados más parecidos a
los experimentales utilizando aproximaciones locales más inexactas (inexactas en el
sentido de dar lugar a resultados peores para gases de electrones). Por ejemplo, Wilk
y Vosko [73] han encontrado que, para átomos, cuando se emplea en conjunción con
una aproximación local para el intercambio, este funcional da lugar a energı́as totales
peores que las obtenidas con el funcional de Gunnarsson y Lundqvist [69]. En cambio,
si se utiliza el intercambio exacto, con el funcional de Vosko, Wilk y Nusair se obtienen
las mejores energı́as totales (en este caso concreto de átomos).
Para rs ≥ 1,
γ
√ . (A.35)
1 + β1 rs + β2 rs
Para rs < 1,
A ln rs + B + Crs ln rs + Drs . (A.36)
Cada función ϵPC (rs ) y ϵFC (rs ) viene determinada por un conjunto diferente de paráme-
tros, como se especifica en la tabla A.2. La elección de los parámetros es tal que cada
función y su primera derivada son continuas para rs = 1.
Utilizando las funciones ϵPC (rs ) y ϵFC (rs ), Perdew y Zunger aproximaron la energı́a
de correlación por electrón de un gas de electrones con densidades de espı́n ρα y ρβ
como
PZ
ϵC (ρα , ρβ ) = ϵPC (rs ) + ϵFC (x) − ϵPC (x) f (ξ),
(A.37)
Ejercicios. 139
γ β1 β2 A B C D
P
ϵC (rs ) -0.1423 1.0529 0.3334 0.0311 -0.048 0.0020 -0.0116
ϵFC (rs ) -0.0843 1.3981 0.2611 0.01555 -0.0269 0.0007 -0.0048
Tabla A.2: Parámetros para las funciones ϵPC y ϵFC que representan respectivamente las
energı́as de correlación por electrón para los casos paramagnético y ferromagnético. El
resultado obtenido con estos parámetros está en unidades atómicas.
donde
ρα − ρβ
ξ = (A.38)
ρα + ρβ
(1 + ξ)4/3 + (1 − ξ)4/3 − 2
f (ξ) = (A.39)
2(21/3 − 1)
1/3
1 3
rs = , (A.40)
a0 4π(ρα + ρβ )
siendo a0 el radio de Bohr. Perdew y Zunger reconocieron que esta interpolación para
polarizaciones de espı́n intermedias es menos realista que la utilizada por Vosko, Wilk
y Nusair. A pesar de esto, la diferencia entre ambas interpolaciones nunca es mayor
que un 1.3 %.
Con esta última expresión, el funcional local gas de electrones de Perdew y Zunger
queda como4 .
Z
EC [ρα , ρβ ] = ρ(r)ϵPZ
PZ
C
(ρα , ρβ ) dr. (A.41)
Hay un esquema para corregir esta autocorrelación, desarrollado por Stoll, Pavlidou
y Preuss [76], y que esta basado en la razonable suposición de que, aunque la correlación
entre electrones con el mismo espı́n es importante en un gas de electrones, ésta debe ser
bastante pequeña en sistemas finitos como átomos y moléculas5 . Por ello, para tener un
funcional adecuado para estos últimos sistemas, uno debe sustraer del funcional local,
ECLSD [ρα , ρβ ], la contribución debida a la correlación entre electrones del mismo espı́n.
Es decir,
ECSIC [ρα , ρβ ] = ECLSD [ρα , ρβ ] − ECLSD [ρα , 0] − ECLSD [0, ρβ ]. (A.42)
Nótese que este funcional es local en el sentido de que sólo depende de las densidades
de espı́n ρα y ρβ . Obviamente, este funcional falla cuando se aplica a gases de electro-
nes uniformes, pero da una energı́a de correlación nula en sistemas mono-electrónicos.
Además, es un funcional que tiene consistencia de tamaño [78]. Para átomos y molécu-
las, donde la suposición realizada anteriormente es razonable, se obtienen resultados
bastante buenos para la energı́a de correlación [76, 78, 79].
Ma y Brueckner
Unos de los primeros[87], desarrollado por series perturbativas. Presenta problemas
de raı́ces imaginarias.
ϵ0C (ρ)
Z
MB 1
EC = (A.45)
2 2 0.32
1 − B(ρ)|∇ρ|
0.32ϵ0 (ρ) C
con
donde a0 es el radio de Bohr. Este funcional debe ser usado junto con la aproximación
de Perdew y Zunger [72] para ECLSD [ρα , ρβ ]. El parámetro f˜ de la ecuación (A.49) se
Ejercicios. 143
ha ajustado a 0.11 para reproducir la energı́a de correlación exacta del átomo de neón.
Señalaremos también que Perdew ha dado [74] la expresión del potencial de correlación
correspondiente a este funcional.
El término dependiente del gradiente se anula para una densidad uniforme, de ma-
nera que este funcional recupera la forma local para gases de electrones uniformes,
por lo que es exacto para estos sistemas. El funcional se ha probado con sistemas
mono-electrónicos, dando una correlación casi nula [74], de manera que puede consi-
derarse como un funcional sin auto-correlación. También proporciona buenas energı́as
de correlación para átomos y moléculas [74, 79, 102].
El Funcional de Becke
La derivación del funcional de Becke [77] está basada en el concepto de hueco, un
concepto familiar en Quı́mica Cuántica (ver, por ejemplo, sección 4.8 de [103], o, sección
2.4 de [45]).
Dentro de la teorı́a del funcional de la densidad también son útiles los huecos (de
intercambio-correlación, de intercambio, y de correlación), aunque la definición es dife-
rente. Esta definición requiere la llamada conexión adiabática [69, 90]: un Hamiltoniano
Ĥλ dependiente de un parámetro λ en la forma
N N N N
1 X 2 XX λ X
Ĥλ = − ∇ + + vλ (ri ), (A.54)
2 i=1 i=1 j>i
∥ r i − rj ∥ i=1
y
′
ρσ (r1 )hσσ (r1 , r2 )
Z Z
1X C
EC = dr1 dr2 . (A.59)
2 σ,σ′ ∥ r1 − r2 ∥
donde γα y γβ son las matrices densidad de primer orden construidas con los
orbitales de Kohn y Sham 7 .
Es bien sabido [69] que sólo el promedio esférico de un hueco alrededor del electrón
′
de referencia situado en r1 contribuye a la energı́a. Denotaremos por hσσ C,λ
(r1 , s)
σσ ′
al promedio esférico del hueco hC,λ (r1 , r2 ) cuando el electrón de referencia está en
r1 , y el segundo electrón se encuentra situado en una esfera de radio s centrada
en r1 . Obsérvese que la notación empleada distingue entre huecos promediados
esféricamente y huecos sin promediar únicamente por los argumentos, (r1 , s) y
(r1 , r2 ) respectivamente. Análogas definiciones se aplican a los restantes huecos.
7
Obsérvese que el empleo de los orbitales de Kohn y Sham en vez de los orbitales Hartree-Fock para
construir el hueco de intercambio implica que se está utilizando una definición de intercambio diferente
de la definición de intercambio Hartree-Fock. La diferencia entre las dos definiciones es esta última
utiliza los orbitales que dan lugar al determinante con menor energı́a total (los orbitales Hartree-Fock),
mientras que la primera emplea los orbitales que dan lugar a la mı́nima energı́a cinética. Para más
detalles ver [104, 92, 105, 106].
Ejercicios. 145
1ρ σσ ′
σ (r1 )hC,λ (r1 , s)
Z Z Z
1X
EC = dr1 ds dλ. (A.72)
2 σ,σ′ 0 s
a.u.)
Z
BEC 2ln(1 + zαβ )
EC [ρα , ρβ ] = −0.8 1−
ρα ρβ zαβ dr
zαβ
Z
4 2 zαα
− 0.01 ρα Dα zαα 1 − ln 1 + dr
zαα 2
Z
4 2 zββ
− 0.01 ρβ Dβ zββ 1 − ln 1 + dr, (A.73)
zββ 2
con
X (∇ρα )2
Dα = |∇ϕiα |2 − (A.74)
iα
4ρα
X (∇ρβ )2
Dβ = |∇ϕiβ |2 − , (A.75)
iβ
4ρβ
donde ϕiα y ϕiβ representan respectivamente orbitales de Kohn y Sham con espı́n α y
β 9.
Finalmente, falta determinar las longitudes de correlación zαβ , zαα y zββ . Becke
supuso que la longitud zσσ′ es proporcional a la suma de los radios del hueco de inter-
cambio (o de Fermi) para electrones de espı́n σ y σ ′ , es decir,
′
zσσ′ = cσσ′ (RFσ + RFσ ) (A.76)
A continuación se da una definición precisa de radio del hueco de Fermi. De las
propiedades resumidas anteriormente, se deduce que el intercambio puede escribirse
como Z Z
1 1
EX = − −1
ρα ⟨sα ⟩ dr − ρβ ⟨s−1
β ⟩ dr, (A.77)
2 2
donde ⟨s−1 σσ
σ ⟩ es la media del inverso del radio del hueco de Fermi |hX | en el punto de
referencia r, Z
−1 1 σσ
⟨sσ ⟩ = |h (r, s)| ds. (A.78)
s X
Esta interpretación sugiere la siguiente definición para el radio de Fermi RFσ ,
1
RFσ = . (A.79)
⟨s−1
σ ⟩
podemos, tras comparar con (A.77), expresar los radios de Fermi de la siguiente forma,
1
RFα = − (A.81)
2ϵαX
1
RFβ = − β. (A.82)
2ϵX
9
En la práctica, se suelen emplear los orbitales Hartree-Fock en vez de los de Kohn y Sham, cambio
que introduce un error despreciable [109].
Ejercicios. 147
Las energı́as de correlación obtenidas con este funcional son similares a las obtenidas
con el funcional de Perdew [74]. Para aplicaciones del funcional de Becke en átomos y
moléculas, ver el trabajo de Miehlich y colaboradores [109]. Para el comportamiento
del funcional de Becke ante el escalado de coordenadas, ver el trabajo de Wilson y Levy
[101] y Levy [100]
ρ = ρα + ρβ (A.88)
1/3
1 3
rs = (A.89)
a0 4π(ρα + ρβ )
ρα − ρβ
ξ = , (A.90)
ρα + ρβ
Z
1 −2/3 5/3 1 1 2 −cρ−1/3
Ec = −a ρ + bρ CF ρ − 2tw + tw + ∇ ρ e dr,
1 + dρ−1/3 9 18
(A.91)
Z
γ(r) n
−5/3
h
8/3
Ec = −a ρ + 2bρ 22/3 CF ρα8/3 + 22/3 CF ρβ − ρtw +
1 + dρ−1/3
1 α β 1 2 2 −cρ−1/3
(ρα tw + ρβ tw ) + (ρα ∇ ρα + ρbeta ∇ ρbeta ) e dr, (A.92)
9 18
donde a = 0.04918, b = 0.132, c = 0.2533 y d = 0.349
1 |∇ρ|2 1 2
tw = − ∇ρ (A.93)
8 ρ 8
3
cF = (3π 2 )2/3 = 2.871234 (A.94)
10
ρ2α (r) + ρ2β (r)
γ(r) = 2 1 − (A.95)
ρ2 (r)
1/3
kF = (3π 2 ρ) , and
εPW91
c (rs , ζ, t) = εPW92
c (rs , ζ) + H0 (rs , ζ, t) + H1 (rs , ζ, t) . (A.97)
g3β 2 t2 + At4
2α
H0 (rs , ζ, t) = ln 1 + , (A.98)
2α β 1 + At2 + A2 t4
1/3
where α = 0.09, β = νCc (0) with ν = (16/π) (3π 2 ) and Cc (0) = 0.004235, and
2α 1
A= . (A.99)
β e−2αεc /(g3 β 2 ) − 1
PW92
H1 (rs , ζ, t) is given by
3
H1 (rs , ζ, t) = ν Cc (rs ) − Cc (0) − cx g 3 t2 e−100g (ks /kF )t ,
4 2 2 2
(A.100)
7
where cx = −0.001667 and
0.002568 + 0.023266rs + 7.389 · 10−6 rs2
Cc (rs ) = + 0.001667. (A.101)
1 + 8.723rs + 0.472rs2 + 0.07389rs3
Ejercicios. 149
1/3
kF = (3π 2 ρ) , and
εPBE
c (rs , ζ, t) = εPW92
c (rs , ζ) + H (rs , ζ, t) . (A.103)
t2 + At4
3 β
H (rs , ζ, t) = γg ln 1 + , (A.104)
γ 1 + At2 + A2 t4
Funcional (PL)
Funcional (P86)
Funcional (PW91)
Funcional (B95)
Funcional (PBE)
where 3 !
tW
fcTPSS = ρεrevPKZB
c 1 + 2.8 εrevPKZB
c . (A.107)
t
150 Cálculos Computacionales de Estructuras Moleculares
PN
In Eq. (A.107), t = (1/2) i=1 ∇ψi∗ · ∇ψi , tW = (1/8) |∇ρ|2 /ρ, and
2 ! 2 X
tW tW
ρσ
εrevPKZB
c = εPBE
c 1 + C (ζ, ξ) − (1 + C (ζ, ξ)) ε̃c,σ , (A.108)
t t σ=↑,↓
ρ
where εPBE
c is Eq. (A.103),
donde Γ(R) es la parte diagonal de la matriz reducida de segundo orden cuando ambos
electrones se encuentran en la misma posición R. De estas relaciones podemos despejar
las cantidades ρα y ρβ en función de ρ(R) y Γ(R) para obtener
p
ρ(R) + ρ(R)2 − 2Γ(R)
ρα (R) = (A.113)
p 2
ρ(R) − ρ(R)2 − 2Γ(R)
ρβ (R) = . (A.114)
2
Estas últimas relaciones sólo son ciertas para el caso Hartree-Fock. Pues bien, el pro-
cedimiento de Moscardó y SanFabián [114] para evitar la doble cuenta de la energı́a de
correlación consiste simplemente en reemplazar, en cualquier funcional dependiente de
las densidades de espı́n, los valores de ρα y ρβ por las ecuaciones (A.113) y (A.114).
El procedimiento mejora considerablemente el comportamiento de estos funcionales
cuando se emplea con densidades CI.
A lo largo de los años, se han propuesto cierto número de funcionales de energı́a
de correlación desarrollados de acuerdo con la definición de energı́a de correlación CI.
Por ejemplo, los trabajos de Lie-Clementi [115, 116], Savin [117, 118], Colle-Salvetti
[119, 120], Gritsenko y colaboradores [121, 122], o Moscardó-SanFabián [123]. Un resu-
men interesante de Savin ha sido publicado recientemente [124]. Aquı́ vamos a revisar
brevemente algunos funcionales, concretamente los de Lie-Clementi, Colle-Salvetti y
Moscardó-SanFabián.
ρ1/3
Z
LC 1/3
EC [ρ] = a1 + b1 ln(1 + b2 ρ ) ρ dr , (A.115)
a2 + ρ1/3
y reemplazaron la densidad ρ por una densidad modificada ρm que depende de la
densidades de los orbitales naturales ρi y de los correspondientes números de ocupación
ni ,
2
X
ρm = ni e−(2−ni ) /2 ρi (A.116)
i
152 Cálculos Computacionales de Estructuras Moleculares
e−0.58/β
Z
CS ΓCI (R, R; R, R)
EC [ρ, ΓCI , ΓHF ] = −0.37588π 1 + 0.173K dR, (A.117)
(β + 0.8)β 2 β2
donde K y β se definen como
▽r ΓCI (R − 2r , R + 2r ; R − 2r , R + 2r )
2
K= (A.118)
ΓCI (R − 2r , R + 2r ; R − 2r , R + 2r ) r=0
y
(▽21 + ▽22 )ΓCI (r1 , r2 ; r,1 , r,2 )
ν
β = q 1+
2 ΓCI (r1 , r2 ; r,1 , r,2 )
#
2 2 , , 2
(▽ + ▽2 )ΓHF (r1 , r2 ; r1 , r2 ) r1 =r1 ,
− 1 , , ,
r2 =r2 ρ(R)1/3 , (A.119)
ΓHF (r1 , r2 ; r1 , r2 ) r=0
Lee, Yang y Parr [111] la han modificado para obtener un verdadero funcional de la
densidad que depende únicamente de las densidades de spin.
Z
MSF
EC [ρ, ΓCI ] = (N − 1) ϵC [ρ, t, β] dr (A.120)
154 Cálculos Computacionales de Estructuras Moleculares
2ρ2
2
αβπ 1/2 α2
2 β
ϵC [ρ, t, β] = ϕ + + 2
π 1/2 β 4 t3/2 2a 2a3/2 2a
1/2
1 1 2 απ 1 1
+ϕ − β + − β
b a 2 b3/2 a3/2
1 1
− − β2 , (A.121)
b 2a
2 2/3
−1 (N − 1)ρ
t(r) = π , (A.122)
ΓCI (r)
p
C22 + 4C1 C3 − C2
ϕ = , (A.123)
2C1
π 1/2 2 2α 3α2 π 1/2
C1 = β + β + , (A.124)
2a3/2 a2 4a
5/2
1/2 1 1 2 1 1
C2 = π − β + 2α 2 − 2 β, (A.125)
b3/2 a3/2 b a
1 1
C3 = π 1/2 3/2 − 3/2 β 2 , (A.126)
b 2a
2/3
ρ
a = 2+ 2, (A.127)
tβ
ρ2/3
b = 1+ 2, (A.128)
tβ
1
α = , (A.129)
2
y β se toma como en la ecuación (A.119), pero con q = 1.8 y ν = 1.5 (hay otras
opciones para β [123]). La forma de la expresión de ϵC [ρ, t, β] es tal que cuando ΨCI es
la función de onda exacta (β = ∞), la correlación obtenida por el funcional (A.120)
también es cero.
Al igual que para el funcional de Colle y Salvetti, en este caso también se ha
tomado como función de onda Hartree-Fock al determinante con mayor coeficiente (en
valor absoluto) de ΨCI . Como consecuencia, cuando se utiliza con funciones de onda
GVB o MCSCF, este método tampoco tiene consistencia de tamaño. Por otra parte, su
dependencia explı́cita en N da lugar, incluso cuando ΨCI es la función Hartree-Fock, a
una correlación para un sistema disociado que es diferente de la suma de las energı́as
de correlación de los componentes aislados. Afortunadamente, ambos efectos tienden a
cancelarse mutuamente.
Una excelente revisión del funcional de Moscardó y San-Fabián, donde se incluye
un estudio compativo con otros funcionales, ha sido publicado recientemente [144].
N N
ψi∗ (r)ψj (r)ψj∗ (r′ )ψi (r′ ) 3 3 ′
Z Z
1 XX
ExHF [{ψi }] =− δσ σ d rd r . (A.130)
2 i=1 j=1 i j |r − r′ |
Ex (LSDA)+Ec (LSDA)+a0 ∗(Ex (HF )−Ex (LSDA))+ax ∗∆−Ex (B88)+ac ∗∆−Ec (P 91)
Pero el funcional utilizado en el Gaussian tiene la forma:
(1−a0 )∗Ex (LSDA)+a0 ∗Ex (HF )+ax ∗∆−Ex (B88)+ac ∗Ec (LY P )+(1−ac )∗Ec (V W N )
Donde tenemos el intercambio de Slater, el intercambio Hartree-Fock y la corrección
del funcional de intercambio de Becke-88, el funcional de correlación de Lee, Yang y
Parr y el funcional local de correlación de Vosko, Wilk y Nusairi[71], respectivamente.
Notar que el funcional e correlación de LYP incluye el VWN como local.
Los valores de los coeficientes determinados por Becke son:
xc HFexch 0.20 slater 0.80 becke88 nonlocal 0.72 lyp 0.81 vwn_5 0.19
Pero esto no sigue ası́ en los restantes Gaussian, ni coincide con los propuestos por
Becke en 1993, hay que ver las opciones del Gaussian que utilizas en cada momento.
Ahora, para el G92/DFT, un funcional se puede expresar como:
P2 ∗Ex (HF )+P1 ∗(P4 ∗Ex (local)+P3 ∗Ex (non−local))+P6 ∗Ec (local)+P5 ∗Ec (non−local)
Por ejemplo, la ruta para generar el funcional de Becke tipo ”half-and-half” seria:
PBE0 PBE 1
ExHF [{ψi }] − ExPBE [ρ] ,
Exc [{ψi }] = Exc [ρ] + (A.131)
4
donde ExPBE [ρ] y EcPBE [ρ] vienen dadas por las Eqs. (A.16) y (A.102), respectivamente.
B2PLYP[148, 149]:
(2006) Considera la energı́a de correlación perturbativa a orden 2 y los funcionales
de Becke y LYP, con tan sólo dos parámetros:
B2PLYP-D:
Añade al anterior una corrección empı́rica de dispersión (D) (R−6
ij ), que depende de
las posiciones de los núcleos.
TPSSh:
Es el hı́brido de TPSS:
On-top density functional: La densidad de espı́n m(r) = ρα (r) − ρβ (r), para un sólo
determinane, se puede escribir como
4Π(r)
R(r) =
ρ(r)2
Con R ≤ 1 en todos los puntos del espacio.
El problema surge con las funciones multidererminantales, donde R puede se mayor
que 1.
Ellos hacen:
Johnson B.G. (1995), In: Modern Density Functional Theory – A Tool for Chemistry,
Theoretical and Computational Chemistry, Volume 2, Seminario J. and Politzer
P., Editors, pp. 169 – 219, Elsevier, Amsterdam, 1995.
Johnson B.G., Gill P.M.W, Pople J.A. (1993), J. Chem. Phys., 98, 5612.
Slater J.C. (1974), Quantum Theory of Molecules and Solids, Volume 4, The self-
consistent field for molecules and solids, McGraw-Hill, New York.
159
160 Cálculos Computacionales de Estructuras Moleculares
cd cambiar de directorio
cd Pruebas
cd ./UNIX
cd ../vi
cd /u/ls
cd
mv renombrar
mv resultado [Link]
cp copiar
cp /u/ls/[Link] [Link]
cp /u/ls/[Link] .
i
162 Cálculos Computacionales de Estructuras Moleculares
Comando Descripción
• apropos palabra Ver comandos relacionados con palabra. Ver también threadsafe
which comando Ver la ruta completa de comando
time comando Medir cuanto tarda comando
• time cat Iniciar cronómetro. Ctrl-d para detenerlo. Ver también sw
• nice info Lanzar comando con prioridad baja (info en este ejemplo)
• renice 19 -p $$ Darle prioridad baja al shell (guión). Usar para tareas no interactivas
dir navegación
• cd - Volver al directorio anterior
• cd Ir al directorio personal (home)
(cd dir && comando) Ir a dir, ejecutar comando y volver al directorio inicial
• pushd . Guardar el directorio actual en la pila para luego, poder hacer popd y volver al mismo
búsquedas de archivo
• alias l=’ls -l - -color=auto’ listado de directorio rápido
• ls -lrt Listar archivos por fecha. Ver también newest
• ls /usr/bin | pr -T9 -W$COLUMNS Imprimir 9 columnas en ancho de la terminal
find -name ’*.[ch]’ | xargs grep -E ’expre’ Buscar ’expre’ en este directorio y subdirectorios. Ver también findrepo
find -type f -print0 | xargs -r0 grep -F ’ejemplo’ Buscar ’ejemplo’ en todos los archivos regulares en este directorio y subdirectorios
find -maxdepth 1 -type f | xargs grep -F ’ejemplo’ Buscar ’ejemplo’ en todos los archivos regulares de este directorio
find -maxdepth 1 -type d | while read dir; do echo $dir; echo cmd2; done Procesar cada elemento con muchos comandos (con un bucle while)
• find -type f ! -perm -444 Hallar archivos sin permiso general de lectura (util para sedes web)
• find -type d ! -perm -111 Hallar directorios sin permiso general de acceso (util para sedes web)
• locate -r ’file[/̂]*\.txt’ Buscar nombres en indice en cache. Este re es igual a glob *file*.txt
• look referencia Búsqueda rápida (ordenada) de prefijo en diccionario
• grep - -color referencia /usr/share/dict/palabras Resaltar ocurrencias de expresión regular en diccionario
archivos
gpg -c file Encriptar archivo
gpg [Link] Desencriptar archivo
tar -c dir/ | bzip2 > [Link].bz2 Crear archivo compacto de dir/
bzip2 -dc [Link].bz2 | tar -x Extraer archivo compacto (usar gzip en vez de bzip2 para archivos [Link] )
tar -c dir/ | gzip | gpg -c | ssh user@remoto ’dd of=[Link]’ Crear compactado encriptado de dir/ en equipo remoto
find dir/ -name ’*.txt’ | tar -c - -files-from=- | bzip2 > dir [Link].bz2 Crear compactado de subconjunto de dir/ y subdirectorios
find dir/ -name ’*.txt’ | xargs cp -a - -target-directory=dir txt/ - -parents Copiar subconjunto de dir/ y subdirectorios
( tar -c /dire/de/copiame ) | ( cd /este/dir/ && tar -x -p ) Copiar (con permisos) directorio copiame/ a directorio /este/dir/
( cd /dire/de/copiame && tar -c . ) | ( cd /este/dir/ && tar -x -p ) Copiar (con permisos) contenido del directorio copiame/ a directorio /este/dir/
( tar -c /dire/de/copiame ) | ssh -C user@remoto ’cd /este/dir/ && tar -x -p’ Copiar (con permisos) directorio copiame/ a directorio remoto /este/dir/
dd bs=1M if=/dev/hda | gzip | ssh user@remoto ’dd of=[Link]’ Respaldo de disco duro en equipo remoto
rsync (Usar la opción - -dry-run para probarlo)
rsync -P rsync://[Link]/ruta/a/archivo archivo Obtenerr solo diffs. Repetir muchas veces para descargas conflictivas
rsync - -bwlimit=1000 desde archivo al archivo Copia local con taza lı́mite. Parecido a nice para E/S (I/O)
rsync -az -e ssh - -delete ∼/public html/ [Link]:’∼/public html’ Espejo de sede web (usando compresión y encriptado)
rsync -auz -e ssh remote:/dir/ . && rsync -auz -e ssh . remote:/dir/ Sincronizando directorio actual con uno remoto
wget (herramienta de descargas multiuso)
• cd cmdline && wget -nd -pHEKk [Link] Guardar en directorio actual una versión navegable de una página web
wget -c [Link] Retomar descarga de un archivo parcialmente descargado
wget -r -nd -np -l1 -A ’*.jpg’ [Link] Descargar una serie de archivos en el directorio actual
wget [Link] FTP permite globalizaciones directas
• wget -q -O- [Link] | grep ’a href’ | head Procesando directamente la salida
echo ’wget url’ | at 01:00 Descargar la url a 1AM al directorio en que esté
wget - -limit-rate=20k url Hacer descargas de baja prioridad (en este caso, no exceder los 20KB/s)
wget -nv - -spider - -force-html -i [Link] Revisando los enlaces de una página
wget - -mirror [Link] Actualizar eficientemente una copia local de una página web (útil si usamos cron)
redes (Nota los comandos ifconfig, route, mii-tool, nslookup son obsoletos)
ethtool interface Listar estado de interfase
• ip link show Listar interfases
ip link set dev eth0 name wan Renombrar eth0 a wan
ip addr add [Link]/24 brd + dev eth0 Agregar ip y máscara ([Link])
ip link set dev interface up Subir (o bajar) la interfase
ip route add default via [Link] Establecer [Link] como valor por omisión para la puerta de enlace.
• tc qdisc add dev lo root handle 1:0 netem delay 20msec Agregarle 20ms de espera al dispositivo de retorno (para hacer pruebas)
• tc qdisc del dev lo root Quitar la espera agregada antes.
• host [Link] Obtener la dirección ip para el dominio o al revés
• hostname -i Obtener la dirección ip local (equivale al anfitrión ‘hostname‘)
• netstat -tupl Listar los servicios de internet de un sistema
• netstat -tup Listar las conexiones activas de/hacia un sistema
windows (nota samba es el paquete que permite todos estos comandos de redes de windows )
• smbtree Hallar equipos windows. Ver también findsmb
nmblookup -A [Link] Hallar el nombre (netbios) de windows asociado con la dirección ip
smbclient -L windows box Listar archivos compartidos en equipos windows o servidor samba
mount -t smbfs -o fmask=666,guest //windows box/share /mnt/share Montar un directorio compartido
echo ’mensaje’ | smbclient -M windows box Enviar mensaje emergente al equipo windows (desactivado por omisión en XP sp2)
math
• echo ’(1 + sqrt(5))/2’ | bc -l Cuentas rápidas (Calcular ¿). Ver también bc
• echo ’obase=16; ibase=10; 64206’ | bc Conversiones de base (decimal a hexadecimal)
• echo $((0x2dec)) Conversiones de base (hex a dec) ((expansión aritmética del shell))
• echo ’pad=20; min=64; (100*106̂)/((pad+min)*8)’ | bc Mas complejo (int) [Link]. Ejemplo: tasa máxima de paquetes FastE
• echo ’pad=20; min=64; print (100E6)/((pad+min)*8)’ | python Python maneja notación cientı́fica
• echo ’pad=20; plot [64:1518] (100*10**6)/((pad+x)*8)’ | gnuplot -persist Graficar tasa de paquetes FastE vs. tamaño de paquetes
• seq 100 | (tr ’\n’ +; echo 0) | bc Agregar una columna de números. Ver también add y funcpy
Ejercicios. 163
manejo de textos (nota: como sed usa stdin y stdout, para editar archivos, agregar... <viejoarchivo >nuevoarchivo)
sed ’s/cadena1/cadena2/g’ Remplaza cadena1 por cadena2
sed ’s/\(.*\)1/\12/g’ Modificar cualquiercadena1 con cualquiercadena2
sed ’/ *#/d; / ˆ *$/d’ Quitar comentarios y lineas en blanco
sed ’:a; /\\$/N; s/\\\n//; ta’ Concatenar lineas con al final
sed ’s/[ \t]*$//’ Quitar blancos finales de las lineas
sed ’s/\([\\‘\\”$\\\\]\)/\\\1/g’ Escapar metacaracteres activos del shell dentro de comillas dobles
sed -n ’1000p;1000q’ Listar la linea 1000
sed -n ’10,20p;20q’ Listar de la linea 10 a la 20
sed -n ’s/.*<title>\(.*\)<\/title¿.*/\1/ip;T;q’ Extraer titulo de página web en HTML
sort -t. -k1,1n -k2,2n -k3,3n -k4,4n Sort de direcciones ip de tipo IPV4
• echo ’Test’ | tr ’[:lower:]’ ’[:upper:]’ Conversión de cajas
• tr -dc ’[:print:]’ < /dev/urandom Filtrando caracteres no imprimibles
• grep ’processor’ /proc/cpuinfo | wc -l Contar lineas
definir operaciones (Nota export LANG=C es para acelerar, aquı́ también se supone que no hay lı́neas duplicadas en los archivos)
sort archivo1 archivo2 | uniq Union de archivos sin ordenar
sort archivo1 archivo2 | uniq -d Intersección de archivos sin ordenar
sort archivo1 archivo1 archivo2 | uniq -u Diferencia de archivos sin ordenar
sort archivo1 archivo2 | uniq -u Diferencia Simétrica de archivos sin ordenar
comm archivo1 archivo2 | sed ’s/ ˆ \t*//’ Unión de archivos ordenados
comm -12 archivo1 archivo2 Intersección de archivos ordenados
comm -13 archivo1 archivo2 Diferencia de archivos ordenados
comm -3 archivo1 archivo2 | sed ’s/ ˆ \t*//’ Diferencia Simétrica de archivos ordenados
calendario
• cal -3 Mostrar calendario
• cal 9 1752 Mostrar calendario para mes y año determinado
• date -d fri Que dı́a cae este viernes. Ver también day
• date - -date=’25 Dec’ + %A ¿En que dı́a cae la Navidad, este año?
• date - -date ’1970-01-01 UTC 1234567890 seconds’ Convertir total de segundos desde la época a una fecha
• TZ=’:America/Los Angeles’ date ¿Que hora es en la Costa Oeste de EEUU (usar tzselect para hallar TZ)
echo ”mail -s ’tomar el tren’ P@[Link] < /dev/null”| at 17:45 Recordatorio por email
• echo ”DISPLAY=$DISPLAY xmessage cooker”| at ”NOW + 30 minutes” Recordatorio emergente
locales
• printf ” %’d\n”1234 Imprimir numero agrupado por miles de acuerdo a su locale
• BLOCK SIZE=1́ ls -l pedir que ls agrupe por miles de acuerdo a su locale
• echo ”Yo vivo en ‘locale territory‘” Extraer información de la base de datos del locale
• LANG=en IE.utf8 locale int prefix Buscar información de locale para determinado paı́s. Ver también ccodes
• locale | cut -d= -f1 | xargs locale -kc | less Listar campos en base de datos del locale
recode (obsoletos: iconv, dos2unix, unix2dos)
• recode -l | less Ver conversiones disponibles (aliases en cada lı́nea)
recode windows-1252.. archivo a [Link] .ansi”de Windows a tabla de caracteres locales (auto hace conversión CRLF)
recode utf-8/CRLF.. archivo a [Link] utf8 de Windows a tabla de caracteres locales
recode iso-8859-15..utf8 archivo a [Link] Latin9 (Europa oriental) a utf8
recode ../b64 < [Link] ¿archivo.b64 Codificado Base64
recode /qp.. < [Link] ¿[Link] Decodificado de citas imprimibles (qp)
recode ..HTML < [Link] ¿[Link] Texto a HTML
• recode -lf windows-1252 | grep euro Buscar tabla de caracteres
• echo -n 0x80 | recode latin-9/x1..dump Mostrar representación de un código en tabla de caracteres latin-9
• echo -n 0x20AC | recode ucs-2/x2..latin-9/x Ver codificado latin-9
• echo -n 0x20AC | recode ucs-2/x2..utf-8/x Ver codificado utf-8
CDs
gzip < /dev/cdrom ¿[Link] Guardar una copia de los datos de cdrom
mkisofs -V NOMBRE -r dir | gzip ¿[Link] Crear imagen de cdrom con el contenido de dir
mount -o loop [Link] /mnt/dir Montar la imagen cdrom en /mnt/dir (solo lectura)
cdrecord -v dev=/dev/cdrom blank=fast Limpiar un CDRW
gzip -dc [Link] | cdrecord -v dev=/dev/cdrom - Grabar un cdrom con imagen (usar dev=ATAPI -scanbus para confirmar ruta dev)
cdparanoia -B Extraer pistas de audio desde un CD a archivos wav en directorio actual
cdrecord -v dev=/dev/cdrom -audio *.wav Armar un CD de audio con todos los wavs en directorio actual (ver también cdrdao)
oggenc - -tracknum=’pista’ [Link] -o ’[Link]’ Crear un archivo ogg con un archivo wav
espacio de disco (Ver también FSlint)
• ls -lSr Mostrar archivos, de menor a mayor
• du -s * | sort -k1,1rn | head Mostrar usuarios de disco principales en el directorio actual. Ver también dutop
• df -h Mostrar espacio libre de disco
• df -i Mostrar inodos libres
• fdisk -l Mostrar tamaños y tipos de particiones de disco (pedir como root)
• rpm -q -a - -qf ’ %10SIZE\t %NAME\n’ | sort -k1,1n Listar todos los paquetes por tamaño instalado (Bytes) de distribuciones RPMs
• dpkg-query -W -f=’$Installed-Size;10\t$Package\n’ | sort -k1,1n Listar todos los paquetes por tamaño instalado (Kbytes) de distribuciones deb
• dd bs=1 seek=2TB if=/dev/null of=[Link] Crear un gran archivo de prueba (sin ocupar espacio). Ver también truncate
monitoreo/rastreo
• strace -c ls ¿/dev/null Resumir/perfil de llamadas al sistema hechas con comando
• strace -f -e open ls ¿/dev/null Listar llamadas al sistema hechas con comando
• ltrace -f -e getenv ls ¿/dev/null Listar llamadas a librerı́as hechas con comando
• lsof -p $$ Listar las rutas que abrió el id de proceso
• lsof ∼ Listar procesos que solicitaron apertura de rutas
• tcpdump not port 22 Ver tráfico de redes excepto ssh. Ver también tcpdump not me
• ps -e -o pid,args - -forest Listar procesos de una jerarquı́a
• ps -e -o pcpu,cpu,nice,state,cputime,args - -sort pcpu | sed ’/ ˆ 0.0 /d’ Listar procesos por % de uso de cpu
• ps -e -orss=,args= | sort -b -k1,1n | pr -TW$COLUMNS Listar procesos por uso de memoria. Ver también ps [Link]
• ps -C firefox-bin -L -o pid,tid,pcpu,state Listar todos los hilos de un proceso determinado
• ps -p 1,2 Listar información de un ID determinado
• last reboot Ver historia de reencendido del sistema
• free -m Ver cantidad de RAM (que queda) (-m muestra en MB)
• watch -n.1 ’cat /proc/interrupts’ Observar continuamente los datos que van cambiando
164 Cálculos Computacionales de Estructuras Moleculares
1
gvim se ha popularizado últimamente, quizás de uso más amigable al contar con un entorno de
tipo gráfico.
165
166 Cálculos Computacionales de Estructuras Moleculares
+ Ir a la linea siguiente
- Ir a la linea anterior
:+8 Ir a la linea que esta 8 puestos más abajo
:-9 Ir a la linea que esta 9 puestos más arriba
:6 Ir a la linea numero 6
:P,U s/texto viejo/texto nuevo/ Substituir texto desde las lineas P a U; solo la
primera vez que aparezca en cada linea. Ejemplos:
:1,$ s/hola/adiós/ substituir el primer ”hola”de
cada linea del fichero por ’adiós’
168 Cálculos Computacionales de Estructuras Moleculares
NOTA: Cada vez que se borra texto, el texto borrado pasa a un ALMACÉN tem-
poral, de donde elimina lo que estuviese almacenado previamente.
Apéndice D
D.1. Introducción:
La mayorı́a de programas de tratamiento de textos son de composición visual o
WYSIWYG (What You See Is What You Get), lo que el autor escribe y ve en la
pantalla del ordenador es, aproximadamente, lo que aparecerá impreso.
Son muy fáciles de aprender a usar, pero no permiten indicar el estilo en el que
se quiere componer cada una de las partes de un trabajo cientı́fico, ni dan apoyo a la
indicación de la estructura del trabajo: división del trabajo en apartados, la situación
del resumen y de la lista de referencias, etc.
La escritura de textos cientı́ficos necesita no solo caracteres de texto y mecanismos
para indicar el estilo de cada una de las partes de un escrito, sino también muchos
caracteres especiales, mecanismos para escribir fórmulas matemáticas complejas, para
componer tablas, para incluir referencias bibliográficas, etc. Ası́ llegamos al desarrollo
de los nombrados sistemas de composición lógica.
TeX es un sistema de composición lógica, creado por D. Knuth, que se ha convertido
en lo más habitual para el tratamiento de textos académicos y cientı́ficos, sobretodo
de textos que incluyen muchos elementos de notación matemática, como tesinas y tesis
doctorales, apuntes de asignaturas, artı́culos de revistas cientı́ficas, comunicaciones a
congresos, libros técnicos y cientı́ficos, etc., con una calidad tipográfica comparable a
la de las mejores editoriales cientı́ficas.
LaTex es una versión de Tex creada por L . Lamport, que incluye definiciones de la
estructura de los principales tipos de trabajos cientı́ficos, artı́culos, libros, informes de
investigación, transparencias, etc., además de facilidades para programar otros tipos
de estructuras.
D.2. Contenidos:
1. Introducción
2. Aspectos básicos
169
170 Cálculos Computacionales de Estructuras Moleculares
Efectos especiales
Referencias cruzadas
Composición de textos matemáticos
Tratamiento de la bibliografı́a
Dibujos, figuras y tablas
3. Aspectos avanzados
Matemáticas superiores
Gráficos y colores
Uso de BibTex
Creación de un ı́ndice alfabético
Modificación de clases y estilos de documentos
Manual-pdf
Manual-html
Manuales de CervanTEX
Un ejemplo:
\documentclass[a4paper,12pt]{article}
\usepackage[spanish]{babel}
%\usepackage[latin1]{inputenc}
\usepackage[utf8x]{inputenc} % recode -t UTF-8..ISO8859-1 file
\usepackage{graphicx}
\usepackage{color}
\textwidth 15.5cm
\textheight 25.0cm
\topmargin -1.0cm
% define el tı́tulo
\title{{\Huge \bf CÁLCULOS COMPUTACIONALES DE ESTRUCTURA DE LA MATERIA} \\
\mbox{ }\\
\mbox{\Large Prácticas}\\
\mbox{ }\\}
\author{ \mbox{ }\\
\mbox{ }\\
{\Large \bf Yo con mi apellido} \\
\mbox{ }\\}
\date{ \hspace{4.5cm} Alicante, \today}
Ejercicios. 171
\begin{document}
% genera el titulo
\maketitle
Texto introductorio. Los párrafos se separan por lineas en blanco. Hay secciones,
\section{Introducción}
Hola
soy un documento de prueba
Más texto
% Ejemplo 1
\ldots cuando Einstein introdujo su formula
\begin{equation}
e = m \cdot c^2 \; ,
\end{equation}
la cual es al mismo tiempo la mas conocida y la
menos bien entendida de la formulas fı́sicas.
% Ejemplo 2
\ldots from which follows Kirchoff s current law:
\begin{flushleft}
Este texto esta\\ alineado a la izquierda.
\LaTeX{} no está tratando de dejar cada lı́nea del mismo largo.
\begin{equation}
\hat{H} \cdot \Psi = \epsilon \cdot \Psi
\label{hamiltonian}
\end{equation}
$$\sum_{k=1}^{n} I_k i= 0 \; .$$
Kirchhoff s voltage law can be derived
\end{flushleft}
\begin{figure}[!hbp]
172 Cálculos Computacionales de Estructuras Moleculares
\begin{center}
\includegraphics[width=0.4\textwidth]{test}
\end{center}
\caption{Todo muy bonito aqui arriba.}
\label{Fig:hola}
\end{figure}
\begin{flushright}
Este texto está \\alineado a la derecha.
\LaTeX{} no está tratando de dejar cada lı́nea del mismo largo.
\end{flushright}
\begin{center}
En el centro\\de la Tierra
\end{center}
\begin{equation}
i \hbar \frac{\partial\psi(r,t)}{\partial t} = H \Psi(r,t)
\label{EQ:HAMIL}
\end{equation}
\item Cálculos
\end{itemize}
Con Tablas:
\begin{table}[h]
\hspace{3cm}\begin{tabular}{|l|rl|r| }
\hline
\hline
SD-CI & 183. & 508. & 1353. \\
CAS-CI & 2. & 5. & 4. \\
Ejercicios. 173
\section{Sección de figuras}
\clearpage
\begin{figure}[!hbp]
\begin{center}
\includegraphics[angle=90, width=0.3\textwidth]{test}
\end{center}
\caption{Una gráfica hecha con xmgrace.}
\label{Fig:test}
\end{figure}
\begin{figure}[!hbp]
\includegraphics[scale=0.30]{cep}
\caption{Una figura de no sé que, en formato eps o png}
\end{figure}
\LARGE
\ldots{} Y aquı́ termina.
\normalsize
\end{document}
174 Cálculos Computacionales de Estructuras Moleculares