Main
Main
numerici
Contents
1 Basi di Algebra lineare e metodi numerici 2
1.1 Algebra lineare in 3 dimensioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Matrici, vettori e loro proprietà . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Spazi Vettoriali, applicazioni lineari e cambi di base . . . . . . . . . . . . . . . . . . 8
1.4 Spazi vettoriali complessi, di Banach e di Hilbert . . . . . . . . . . . . . . . . . . . . 11
1.5 Teoria spettrale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.6 Norme matriciali . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.7 Singular Value Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.8 Successione di matrice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.9 Limiti dell’aritmetica nei calcolatori . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.10 Sistemi lineari perturbati . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
1.11 Metodi Diretti e Precisione Numerica . . . . . . . . . . . . . . . . . . . . . . . . . . 27
1.12 Metodi diretti di Risoluzione di sistemi lineari . . . . . . . . . . . . . . . . . . . . . . 29
1.12.1 Metodo di Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
1.12.2 Fattorizzazione LU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
1.12.3 Decomposizione di Cholesky . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
1.13 Metodi Iterativi per la Risoluzione di Sistemi Lineari . . . . . . . . . . . . . . . . . . 40
1.13.1 Fixed points methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
1.13.2 Krylov subspaces e Gradiente coniugato . . . . . . . . . . . . . . . . . . . . . 46
1.14 Calcolo di autovalori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
1.14.1 Il metodo di Davidson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
1.14.2 Problemi Numerici del metodo di Davidson . . . . . . . . . . . . . . . . . . . 64
1
Basi di Algebra lineare e metodi nu-
merici
1.1 Algebra lineare in 3 dimensioni
Un vettore in 3 dimensioni può essere rappresentato dai suoi componenti ai rispetto ad un insieme
di 3 vettori unitari mutualmente perpendicolari {êi }:
X
⃗r = ê1 a1 + ê2 a2 + ê3 a3 = êi ai ([Link])
i
L’insieme {êi } si dice formare una base ed ognuno dei 3 vettori si dice versore. La base è detta
completa nel senso che ogni vettore generico ⃗a può essere scritto come combinazione lineare dei
versori. Ovviamente non esiste una sola base per lo spazio e quindi possiamo esprimere lo stesso
vettore ⃗a come:
3
X
⃗r = ϵ̂i bi ([Link])
i
3
dove {ϵ̂i } è un’altra base di R , e {bi } è l’insieme dei coefficienti di ⃗a rispetto alla nuova [Link]
ora un modo di esprimere il vettore in una notazione più compatta; per farlo definiamo il prodotto
scalare tra due vettori: X
⃗a · ⃗b = ai bi ([Link])
i
e notiamo che: X
⃗a · ⃗a = a2i = |⃗a|2 ([Link])
i
ora sostituendo Eq. [Link] nella definizione di prodotto scalare possiamo scrivere:
XX
⃗a · ⃗b = êi · êj ai bj ([Link])
i j
l’unico modo per cui Eq. [Link] e Eq. [Link] siano uguali è che valga la relazione:
êi · êj = δij ([Link])
che dimostra che la base è ortonormale, ovvero i 3 versori sono ortogonali tra loro e unitari. Possiamo
anche calcolare la proiezione di un vettore su un versore come:
X
êi · ⃗r = êi · êj aj = δij aj = ai ([Link])
j
2
Questa relazione è detta di completezza e I⃗ è detto diadico unitario. Concludiamo enunciando
le proprietà del prodotto scalare; queste proprietà saranno quelle che da ora in poi definiranno
l’operazione di prodotto scalare:
O⃗r = ⃗b ([Link])
O(x⃗a + y⃗b) = xO⃗a + yO⃗b ([Link])
ovvero quando fatto operare su un vettore lo trasporta in un altro vettore nello spazio. Se facciamo
operare questo operatore su uno dei versori otterremmo un altro vettore generico dello spazio e
possiamo quindi scrivere: X
Oêi = ⃗o = êj Oij ([Link])
j
dove Oij è la componente del vettore êi lungo il versore êj . Siccome ogni vettore può essere espresso
come combinazione lineare dei versori per comprendere l’effetto di O su un vettore generico può
essere sensato esprimere l’operatore in funzione della base:
X
O= êi Oij ([Link])
i,j
possiamo considerare O come la rappresentazione di O rispetto ad una generica base {êi }. Per co-
modità scriviamo questa rappresentazione come una matrice, ovvero come una lista bidimensionale
di elementi:
O11 O12 O13
O = O21 O22 O23 ([Link])
O31 O32 O33
questa matrice specifica completamente come l’operatore O agisce su un generico vettore nello
spazio, visto esso può essere espresso come combinazione linare dei versori di base. Prendiamo
quindi un operatore C = AB e notiamo che:
X
Cêj = êi Cij ([Link])
i
= ABêj ([Link])
X
=A êk Bkj ([Link])
i
X
= êi Aik Bkj ([Link])
k,i
3
X
Cij = Aik Bkj ([Link])
k
Ovvero possiamo scrivere la rappresentazione matriciale di C in {êi } come:
X
C= Aik Bkj ([Link])
i,j,k
possiamo prendere questa come definizione del prodotto tra matrici. Notiamo come ci dica che ogni
elemento della matrice prodotto possa essere visto come il prodotto scalare tra la riga i-esima della
prima matrice e la colonna j-esima della seconda matrice partendo da sinistra. Il prodotto è quindi
un prodotto non commutativo.
Per via delle considerazioni precedenti e poiché non abbiamo fatto alcuna assunzione sulle di-
mensioni delle matrici (si vede dalla definizione che però il numero di righe di una deve essere uguale
al numero di colonne dell’altra e viceversa) possiamo usare la definizione di prodotto matriciale per
definire a sua volta un prodotto tra due vettori:
X
⃗a⃗b = êi êj ai bj = ([Link])
i,j
applicando la stessa definizione al prodotto tra versori e considerando la loro ortonormalità otteni-
amo: X
⃗a⃗b = a i bi ([Link])
i
che è proprio la definizione del prodotto scalare (si può infatti verificare che questo prodotto ma-
triciale soddisfa tutte le proprietà del prodotto scalare). In notazione matriciale possiamo quindi
scrivere il prodotto scalare come:
b1
⃗a · ⃗b = (a1 , a2 , a3 ) b2 ([Link])
b3
Usando queste nozioni possiamo scrivere Eq. [Link] come:
r = ea ([Link])
dove e è un vettore riga in cui ogni elemento è un versore (possiamo quindi anche vederlo come una
matrice che ha come colonne i versori) e a è un vettore riga in cui il generico elemento i-esimo è il
coefficiente ai del vettore r nella base di e. Concludiamo il paragrafo notando che:
r1 · r2 = (er1 )T (er2 ) ([Link])
dove con aT indichiamo la matrice ottenuta scambiando righe e colonne di a (nel caso di un vettore
riga esso diventa un vettore colonna e viceversa). Si vede facilmente che la relazione può essere
riscritta come:
r1 · r2 = rT1 eT er2 = rT1 Mr2 ([Link])
dove con M intendiamo la matrice metrica dello spazio:
ê1
M = eT e = ê2 (ê1 , ê2 , ê3 ) ([Link])
ê3
che nel caso di una base ortonormale è semplicemente la matrice identità I.
4
1.2 Matrici, vettori e loro proprietà
Ora che abbiamo visto come la nozione di matrice e vettore emerge naturalmente nello spazio in 3
dimensioni, passiamo ad una descrizione più generale ed assiomatica di questi oggetti:
Definition 1.1. Preso un insieme di valori in generale complessi {Aij } con pedici ordinati i =
1, 2, ...., N e j = 1, 2, ...., M , possiamo definire la matrice rettangolare A (N × M ) con N righe ed
M colonne come:
A11 A12 .. A1M
A21 A22 .. A2M
A= :
([Link])
: :
AN 1 AN 2 .. AN M
Se si verifica che N = M la matrice è detta quadrata. Il prodotto tra due generiche matrici
(N × M ) e (M × P ) è definito in Eq. [Link] ed il risultato è un’altra matrice (N × P ):
C = AB ([Link])
N X
X P N X
X P X
M
C= Cij = Aik Bkj ([Link])
i j i j k
Aa = b ([Link])
5
Notiamo subito che se la matrice A è reale allora la sua aggiunta è uguale alla trasposta; ovvero
la matrice ottenuta invertendo righe e colonne:
L’aggiunta di una matrice colonna è una matrice riga e viceversa. Consideriamo ora il prodotto
matriciale tra l’aggiunta di una matrice colonna e un’altra matrice colonna:
b1
M
† ∗ ∗ ∗ b2
X
a b = (a1 , a2 , ..., aM ) = a∗i bi ([Link])
..
i=1
bM
dove abbiamo usato la definizione di prodotto matriciale. Nel caso in cui a sia reale riotteniamo
la definizione di prodotto scalare, vedremo infatti che Eq. [Link] sarà la base per introdurre un
operazione analoga al prodotto scalare nel campo dei complessi.
† † †
Prendendo l’aggiunta di Eq. [Link] (si dimostra facilmente che (AB) = B A ):
† † †
b =a A ([Link])
M M
X ∗
∗ †
X
∗
bi = aj Aji = aj Aij ([Link])
j=1 j=1
† † †
dove b , a sono matrici colonna con M elementi ed A è una matrice (N × M ).
Elenchiamo quindi alcune proprietà e definizioni importanti per le matrici quadrate:
con Aii ̸= 0.
6
5. Una matrice è detta unitaria se la sua inversa è uguale alla sua aggiunta:
†
A−1 = A ([Link])
A−1 = AT ([Link])
si vede facilmente che una matrice ortogonale preserva la matrice metrica dello spazio vetto-
riale.
6. Una matrice è detta Hermitiana se è uguale alla sua aggiunta (una matrice Hermitiana reale
è detta simmetrica):
†
A=A ([Link])
dove Pi è un operatore di permutazione degli indici delle colonne e corre su tutte le N ! permutazioni,
mentre pi è il numero di trasposizioni necessarie per portare una permutazione generica all’ordine
naturale.
Il determinante ha le seguenti proprietà:
1. Se ogni elemento di una riga o di una colonna è nullo allora il determinante è nullo
2. Se la matrice è diagonale allora: Y
|A| = Aii ([Link])
i
3. Uno scambio tra due righe o colonne fa si che il determinante cambi di segno
†
4. |A| = (|A |)∗
5. |AB = |A||B|
6. Se due colonne o righe generiche di un determinante sono uguali, allora il determinante è nullo
7. |A−1 | = (|A|)−1 |
†
8. Se AA = I allora |A|(|A|)∗ = 1 e possiamo anche scrivere:
†
A BA = C ([Link])
|B| = |C| ([Link])
9.
Si può dimostrare inoltre che A−1 esiste se e solo se |A| ̸= 0.
7
1.3 Spazi Vettoriali, applicazioni lineari e cambi di base
Dopo aver formalizzato meglio il concetto di matrice (e quindi di vettore) e le proprietà di questo
oggetto passiamo a formalizzare l’idea di spazio vettoriale:
Definition 1.5. Definiamo uno spazio vettoriale su un campo numerico X come l’insieme degli
elementi che soddisfano le seguenti proprietà:
• ∀u, v, w ∈ V → (u + v) + w = u + (v + w)
• ∃O ∈ V t.c. : O + u = u + O = u∀u ∈ V
• ∀u ∈ V ∃(−u) ∈ V t.c. : u + (−u) = O
• ∀u, v ∈ V : u + v = v + u
• ∀u, v ∈ V e ∀c ∈ X : c(u + v) = cu + cv
• ∀a, b ∈ X e ∀v ∈ V : (a + b)v = av + bv
• ∀a, b ∈ X e ∀v ∈ V : (ab)v = a(bv)
• ∀u ∈ V ∃a ∈ X t.c. : au = u
Ogni elemento v ∈ V è detto vettore.
Possiamo quindi definire una combinazione lineare come:
n
X
x i vn ([Link])
i
con xi ∈ X∀xi ed vi ∈ V ∀vi . Si vede facilmente che l’insieme delle possibili combinazioni lineari è
un sottospazio vettoriale W di V . Questo sottospazio è definito come il sottospazio generato dai
vettori {vi } e nel caso in cui abbiamo che W = V allora diciamo che l’insieme {vi } genera V e ogni
elemento di V può essere scritto come combinazione lineare di {vi }.
Possiamo ora definire una base per il nostro spazio vettoriale:
Definition 1.6. Diciamo che {vi } è una base di V se l’insieme genera lo spazio e se vale la seguente
condizione: X
ai vi = 0 ↔ ai = 0∀ai ([Link])
i
allora xi = yi ∀i.
8
P
la dimostrazione è triviale. Abbiamo quindi dimostrato che il vettore v = i xi vi è considerabile
come il vettore delle coordinate {xi } ∈ X∀i di v rispetto ad {vi }.
Visto in modo più rigoroso la definizione di spazio vettoriale passiamo alla definizione di appli-
cazione lineare:
Definition 1.7. Definiamo un applicazione lineare come come una funzione tra due spazi vettoriali
sullo stesso campo con le seguenti proprietà:
• f (x + y) = f (x) + f (y)
• f (ax) = af (x)
con x, y ∈ V e a ∈ X.
Definita meglio l’idea di applicazione lineare passiamo a sviluppare delle relazioni fondamentali
per questo tipo di funzioni. In particolare iniziamo notando che Eq. [Link] può essere riscritta,
sfruttando le proprietà degli operatori lineari, in forma matriciale come:
dove e è la matrice riga dei versori e Oi è una matrice colonna con i coefficienti di ê′i rispetto ad e.
A questo punto si vede facilmente che:
Ôe = eO ([Link])
dove e è la matrice riga con un versore in ogni elemento e O è una matrice (N ×N ) con N dimensione
dello spazio definito da e. Ogni colonna della matrice O è una matrice colonna Oi . Prendiamo per
esempio il caso di N = 3:
La colonna Ox quindi ci dirà le componenti del vettore Ôêx = ê′x rispetto alla base e; un ragiona-
mento analogo può essere fatto per le altre due colonne. Sappiamo inoltre che, usando Eq. [Link],
possiamo scrivere:
Ô⃗r = Ôer ([Link])
E sostituendo Eq. [Link]:
Ô⃗r = Ôer = eOr = r⃗′ ([Link])
il vettore r⃗′ tuttavia può sempre essere scritto nella base e e quindi possiamo scrivere:
da cui deriva:
r′ = Or ([Link])
Questa relazione è fondamentale in quanto ci dice che le componenti di un generico vettore nella
base e′ possono essere ottenute come il prodotto delle componenti dello stesso vettore nella base e
per una matrice, detta di cambiamento di base, che esprime le componenti della base e′ rispetto
alla base e.
Per completare la trattazione sulle applicazioni lineari passiamo a dimostrare la relazione di
similarità. Consideriamo l’applicazione lineare ed invertibile P̂ tale che:
P̂ e = e′′ = eP ([Link])
9
consideriamo ora un applicazione definita come il prodotto R̂P̂ con R̂ generica applicazione lineare:
R̂P̂ e = R̂e′′ = e′′ R′′ ([Link])
dove R′′ è la matrice di rappresentazione di R′′ rispetto alla base e′′ . Usando Eq. [Link] in Eq.
[Link] la possiamo riscrivere come:
R̂P̂ e = ePR′′ ([Link])
abbiamo appena dimostrato che la matrice PR′′ è la rappresentazione di R̂P̂ rispetto alla base e.
Se consideriamo ora l’effetto su un vettore generico abbiamo:
R̂P̂ ⃗r = ePR′′ r ([Link])
Tuttavia possiamo scrivere Eq. [Link] anche come:
R̂P̂ ⃗r = R̂ePr = eRPr ([Link])
Ovvero abbiamo dimostrato che R′′ P e RP sono due rappresentazioni dello stesso operatore R̂P̂
nella stessa base e. Siccome esiste sono una matrice di rappresentazione di un operatore in una
certa base deve essere:
R′′ P = RP ([Link])
e siccome abbiamo assunto che P fosse invertibile:
R′′ = P−1 RP ([Link])
Questa è la vera e propria relazione di similarità e ci dice che le rappresentazioni di uno stesso
operatore rispetto a due basi diverse possono essere legate tra loro attraverso un’altra trasformazione
invertibile che rappresenta il cambiamento di base da e a e′′ . Inoltre in Eq. [Link] abbiamo
dimostrato che la matrice di rappresentazione di un operatore prodotto in una base è il
prodotto tra le matrici di rappresentazione nella stessa base. Si vede facilmente che due
matrici simili hanno stesso determinante e stessi autovalori.
Concludiamo il paragrafo notando che tutte le considerazioni fatte su spazi vettoriali e loro basi
possono essere estese considerando come base un insieme di funzioni:
f = (f1 (x)f2 (x), ...., fn (x)) ([Link])
Una generica funzione può essere scritta come:
X
F (x) = ci fi (x) ([Link])
i
per uno spazio del genere possiamo considerare come prodotto scalare l’integrale tra due
funzioni dello spazio. Per comodità introduciamo la notazione di Dirac, per la quale possiamo
scrivere un vettore generico come: X
|Ψ⟩ = ci |fi ⟩ ([Link])
i
il generico prodotto scalare è indicato come:
⟨Ψj ||Ψi ⟩ ([Link])
ed è facile riottenere tutte le relazioni dei capitoli precedenti con questa notazione. D’ora in poi
questa è la notazione che sarà usata per comodità.
Nel prossimo paragrafo estenderemo la nozione di spazio vettoriale ai complessi definendo in
esso anche un prodotto scalare, per poi definire un idea di metrica e distanza per spazi astratti.
10
1.4 Spazi vettoriali complessi, di Banach e di Hilbert
Abbiamo visto in Eq. [Link] che è possibile definire un prodotto simile al prodotto scalare nel caso
dei numeri complessi. Un prodotto del genere è detto prodotto Hermitiano e soddisfa:
• Sesquilinearità:
Queste proprietà fanno si che il prodotto scalare tra due vettori sia sempre un reale. Nel caso
in cui si trovi un prodotto che non soddisfi la simmetria Hermitiana esso sarà definito una forma
sesquilineare. Se il nostro prodotto soddisfa:
• ∀v ∈ V, ⟨v||v⟩ ≥ 0
• ⟨v||v⟩ = 0 ↔ v = 0
allora esso sarà detto definito positivo. Usando un prodotto definito positivo possiamo sempre
definire una norma; ovvero un applicazione V → R+ tale che:
• ||v|| ≥ 0∀v ∈ V, ||v|| = 0 ↔ v = 0
• ∀v ∈ V, ∀λ ∈ C, ||λv|| = |λ|||v||
Rr = λr ([Link])
11
stiamo fondamentelmente dicendo che il vettore r è invariante rispetto alla trasformazione indotta
da R e viene da essa solo scalato. r è detto autovettore rispetto ad R, mentre λ è detto autovalore
di r rispetto ad R. Si dimostra facilmente che λ è autovalore di R se e solo se:
det(R − λI) = 0 ([Link])
ovvero se λ è radice del polinomio caratteristico di R, si vede da Eq. [Link] che se la matrice
è diagonale gli elementi diagonali sono gli autovettori. Si dimostra sempre facilmente che due
autovettori associati ad autovalori differenti sono tra loro linearmente indipendenti, da ciò deriva
che esiste un sottospazio vettoriale di dimensione uguale al numero di autovalori non
nulli e unici. Diciamo quindi che una matrice è diagonalizzabile se esiste una base di V composta
di autovettori della matrice. Inoltre matrici simili hanno autovalori uguali e quindi questi
possono essere calcolati rispetto a qualsiasi base:
det(A − λI) = det(A − λBB−1 ) =det(B(B−1 A − λB−1 )) = det(BB−1 (B−1 AB − λI)) ([Link])
det(B−1 AB − λI) ([Link])
Passiamo ora ad enunciare il teorema spettrale, centro del paragrafo.
Theorem 1.2. Sia A ∈ Mn,n (C) una matrice hermitiana a coefficienti complessi. Esiste una base
ortonormale e di autovettori per A. La matrice sarà inoltre diagonalizzabile.
Dimostrazione: La dimostrazione va per induzione. Per il caso base n = 1 avremo:
A = (a11 ) ([Link])
che è per definizione diagonale e quindi ha un autovalore a11 . Siccome gli autovalori sono uguali
per matrici simili è banale che:
A = (a11 ) = (1)(A′ )(1) = (1)(a11 )(1) = P−1 A′ P ([Link])
ed ogni vettore sarà un autovettore. Passo induttivo: Supponiamo che il teorema valga per tutte
le matrici hermitiane di ordine n − 1 e dimostriamolo per una matrice A di ordine n.
• Passo 1 - Esistenza di un autovalore: Poiché lavoriamo in C, il polinomio caratteristico
det(A − λI) ha almeno una radice λ1 ∈ C. Quindi esiste un autovettore ⃗v1 ̸= ⃗0 tale che
A⃗v1 = λ1⃗v1 . Normalizziamo ⃗v1 in modo che ∥⃗v1 ∥ = 1.
• Passo 2 - Proprietà delle matrici hermitiane: Poiché A è hermitiana (A = A∗ ),
l’autovalore λ1 è reale. Infatti, A⃗v1 = λ1⃗v1 implica A∗⃗v1 = λ1⃗v1 = λ1⃗v1 .
• Passo 3 - Il sottospazio ortogonale è invariante: Sia W = {⃗v1 }⊥ = {w ⃗ ∈ Cn :
⟨w,
⃗ ⃗v1 ⟩ = 0} il sottospazio ortogonale a ⃗v1 . Mostriamo che A(W ) ⊆ W . Sia w
⃗ ∈ W , cioè
⟨w,
⃗ ⃗v1 ⟩ = 0. Calcoliamo:
⟨Aw, ⃗ A∗⃗v1 ⟩
⃗ ⃗v1 ⟩ = ⟨w, (per la proprietà del prodotto hermitiano) ([Link])
= ⟨w,
⃗ A⃗v1 ⟩ (poiché A è hermitiana) ([Link])
= ⟨w,
⃗ λ1⃗v1 ⟩ ([Link])
= λ1 ⟨w,
⃗ ⃗v1 ⟩ ([Link])
= λ1 · 0 = 0 ([Link])
⃗ ∈ W per ogni w
Quindi Aw ⃗ ∈ W , ossia W è invariante per A.
12
• Passo 4 - Restrizione e dimensione: Sia A|W la restrizione di A al sottospazio W .
Osserviamo che:
– dim W = n − 1 (abbiamo tolto la dimensione generata da ⃗v1 )
– A|W : W → W è ancora hermitiana, poiché per ogni ⃗u, w
⃗ ∈ W:
⟨A|W ⃗u, w⟩
⃗ = ⟨A⃗u, w⟩
⃗ = ⟨⃗u, Aw⟩
⃗ = ⟨⃗u, A|W w⟩
⃗ ([Link])
13
è detta norma della matrice A se verifica in più:
• ||AB|| ≤ ||A||||B|| ∀A ∈ Cm×n , B ∈ Cn×l
allora è detta norma matriciale submoltiplicativa. Un esempio di norma matriciale è la norma
di frobenius:
m X
X n 1/2
||A||F = |aij |2 ([Link])
i j
che si vede facilmente soddisfare tutti i requisiti. Inoltre da questa norma possiamo definire un
prodotto scalare detto di frobenius:
Concentriamoci ora sulle norme matriciali indotte (o naturali), le quali sono definite dalla
seguente quantità:
||Ax||
||A|| = maxx̸=0 ([Link])
||x||
Che è equivalente ad:
||Ax||
||A|| ≥
||x||
Da cui vediamo che:
||Ax|| ≥ ||A||||x||
Che una relazione che useremo frequentemente. Vediamo inoltre che la definizione di norma matri-
ciale indotta è equivalente ad:
||A|| = max||x||=1 ||Ax|| ([Link])
e per mostrarlo prendiamo un generico vettore:
x
u= ([Link])
||x||
||Ax|| ||Au||||x||
||A|| = maxx̸=0 = maxx̸=0 = maxx̸=0 ||Au|| ([Link])
||x|| ||x||
ma siccome u è per costruzione unitario ci basta cercare il massimo tra i vettori unitari e quindi:
come volevasi dimostrare. Vediamo quindi la norma spettrale (norma naturale 2):
e notiamo che se la eleviamo al quadrato allora questa sarà il prodotto scalare del vettore y = Ax
con se stesso:
14
dobbiamo quindi massimizzare la matrice M = A∗ A. Notiamo che questa è sempre una matrice
quadrata semidefinita positiva e vale quindi il teorema spettrale. Inoltre ricordiamo che una ma-
trice unitaria non varia la norma euclidea di un vettore e quindi non varia neanche
quella indotta di una matrice. Usando il teorema spettrale possiamo quindi riscrivere Eq.
[Link] come:
Xn
||A||22 = max||y||=1 y∗ Λy = max||y||=1 λi |yi |2 ([Link])
i
e porle uguali a 0: (
λi yi = γyi
([Link])
||y||22 = 1
che possiamo riarrangiare nella forma matriciale come:
(
Λy = γy
([Link])
||y||22 = 1
ed otteniamo che i punti stazionari della forma quadratica sono tutti i punti (x1 , x2 , ..., xn ) in cui
abbiamo che x è autovettore di M con autovalore γ e ha norma unitaria. Sostituiamo quindi questi
risultati nella forma quadratica:
la forma quadratica assume come valori nei punti stazionari proprio i valori degli autovalori; ne
segue che la norma naturale euclidea sarà:
p
||A||2 = ρ(A∗ A) ([Link])
15
dove con ρ(A∗ A) abbiamo indicato il massimo autovalore (spesso chiamato raggio spettrale).
Notiamo poi che nel caso in cui A sia hermitiana allora Eq. [Link] si riduce ad:
||A||2 = ρ(A) ([Link])
vedremo successivamente che per matrici rettangolari Eq. [Link] corrisponde proprio al valore
singolare maggiore.
Possiamo quindi definire il rango della matrice A come il numero delle righe o delle colonne
linearmente indipendenti (è facile dimostrare che i due numeri sono esattamente uguali); per
questa quantità abbiamo la seguente relazione:
rank(A) + dim(Ker(A)) = n ([Link])
se però prendiamo una matrice Hermitiana e quadrata (ovvero una per cui vale il teorema
spettrale) sappiamo che esiste una base {⃗vi } di autovettori tali che:
B⃗vi = λi⃗vi ([Link])
ma il set {λi⃗vi } è ancora una base ed è per definizione l’immagine della matrice B e la sua
dimensione (ovvero quella del rango) sarà uguale al numero di autovalori non nulli. Per una
matrice hermitiana quadrata avremo quindi:
rank(B) = #autovalori non nulli
([Link])
Ker(B) = #autovalori nulli
16
Inoltre è facile verificare che, presa la base {⃗xi } di r vettori del row space di A e la base {⃗yi }
di k vettori del kernel di A allora abbiamo che l’insieme:
AT A⃗
vi = λi v⃗i ([Link])
per ogni vi della base. Possiamo definire i valori singolari di A come i seguenti parametri:
p
σi = λi = |µi | ([Link])
che sono sempre reali visto che AT A è semidefinita positiva. Prendiamo ora una matrice di m
colonne, ognuna formata da un vettore della base {⃗vi }:
dove Σ è una matrice diagonale con i seguenti elementi: Σij = µii δij . Assumiamo ora che
rank(A) = r ≤ m; è facile verificare che rank(A) = rank(AT A) e quindi avremo che:
(
σi se i=1,2,...,r
σi = ([Link])
0 se i=r+1,...,n
1 1 σj2 σj
⟨⃗ui , ⃗ui ⟩ = ⟨A⃗vi , A⃗vj ⟩ = ⟨⃗vi , AT A⃗vj ⟩ = ⟨⃗vi , ⃗vj ⟩ = δij ([Link])
σi σj σi σj σi σj σi
E possiamo prendere i rimanenti m − r vettori (i quali formano una base per kern(AT A)) in modo
da formare una base ortonormale di Rn e formiamo la seguente matrice:
17
e calcoliamo il seguente prodotto:
UΣVT = [⃗u1 σ1 |...|⃗ur σr |0|0|...|0]VT
([Link])
= [A⃗v1 |...|A⃗vr |0|0|...|0]VT = AVVT = A
Ovvero possiamo scrivere qualsiasi matrice rettangolare come:
σ1 0 0 ... 0 ... 0 1
v1 v12 ... v1m
0 σ2 0 ... 0 ... 0 rank(A)
v21 v22 ... v2m X
σi ⃗ui⃗viT ([Link])
0
[⃗u1 , ⃗u2 , ..., ⃗um ] 0 σ3 ... 0 ... 0 =
.. .. ... ..
0 0 0 ... σr ... 0 1 2 m i
vm vm ... vm
0 0 0 ... 0 ... 0
con ⃗ui , ⃗vi rispettivamente le colonne di U e V, le quali sono matrici unitarie. Possiamo quindi
esprimere una matrice rettangolare in come prodotto di due matrici unitarie e una
matrice quasi diagonale. Possiamo quindi fare con facilità quindi praticamente tutto quello che
potremmo fare con la decomposizione spettrale di una matrice quadrata.
La SVD è utile poiché permette anche di approssimare in modo efficace e con un errore razional-
izzabile ogni matrice. Per mostrarlo prendiamo una matrice A di rango r e definiamo la seguente
matrice di rango k ≤ r:
X k
Ak = σi ui viT ([Link])
i
vale che:
minrank(X)=K ||A − Ak ||2 = σk+1 ([Link])
(m×n)
per dimstrarlo prendiamo una generica matrice B ∈ C con rank(B) = k; avremo che:
dim(Kern(B)) = n − k ([Link])
Consideriamo ora lo spazio vettoriale spannato dai k + 1 vettori vi colonne di V; l’intersezione di
tale spazio con Kern(B) è sicuramente non nulla in quanto (n − k) + k + 1 = n + 1 > n. Prendiamo
quindi un vettore z con ||z||2 = 1 appartenente a questa intersezione, avremo per questo vettore:
Bz = 0 ([Link])
e scomponendo la matrice A abbiamo:
r
X r
X
Az = σi ui viT z = σi ui αi ([Link])
i i
con αi = viT z i coefficienti della combinazione lineare di z nella base dei vi . Siccome però z è
spannato solo da k + 1 ≤ rr, allora tutti gli altri coefficienti saranno nulli:
k+1
X
Az = σi ui αi ([Link])
i
e avremo che:
k+1
X
||Az||22 = σi2 |αi |2 ([Link])
i
18
dove abbiamo sfruttato l’unitarietà dei ui . Usando poi Eq. [Link] possiamo riscrivere Eq. [Link]
come:
||Az||22 = ||(A − B)z||22 ≤ ||(A − B)||22 ||z||22 = ||(A − B)||22 ([Link])
per costruzione. Otteniamo quindi:
k+1
X
||(A − B)||22 ≥ σi2 |αi |2 ([Link])
i
possiamo ora supporre senza perdità di generalità che i valori singolari siano ordinati in ordine
decrescente:
σ1 ≥ σ2 ≥ ... ≥ σr ([Link])
σ12 ≥ σ22 ≥ ... ≥ σr2 ([Link])
Per vedere il valore minimo facciamo una minimizzazione di lagrange vincolata:
k+1
X k+1
X
L(⃗
α) = σi2 |αi |2 − λ( |αi |2 − 1) ([Link])
i i
abbiamo quindi dimostrato che possiamo approssimare ogni matrice A di rango r con una
2
matrice B di rango k < r e il minimo errore possibile è σk+1 . E’ facile dimostrare che:
||A − Ak ||2 = σk+1 ([Link])
ricordando che la norma 2 è invariante per trasformazioni unitarie, quindi:
||A − Ak ||2 = ||UT (A − Ak )V||2 = ||U(A − Ak )VT ||2 = ||UAVT − UAk VT ||2 ([Link])
= ||diag(0, ..., σk+1 , ..., σr )|| = σk+1 ([Link])
per definizione della norma 2. Quindi l’approssimazione Ak è sotto un certo punto di vista la
migliore approssimazione possibile per la matrice A e quindi la SVD consente dunque di ot-
tenere facilmente approssimazioni di rango inferiore di una data matrice, ovvero di condensare
l’informazione contenuta in A in una matrice di rango inferiore.
19
1.8 Successione di matrice
Se abbiamo una successione di matrici, possiamo definire il limite anche in questo contesto, poiché
siamo in grado di valutare la distanza tra due matrici utilizzando una norma matriciale. Analoga-
mente ai numeri reali, possiamo parlare anche di serie di matrici. Come per le serie numeriche,
una serie di matrici converge se e solo se la successione delle somme parziali converge, cioè
Sk = A1 + A2 + · · · + Ak
Questo tipo di serie è rilevante perché richiama direttamente le serie di Taylor, che permettono di
definire funzioni analitiche di matrici come l’esponenziale, il seno, il coseno, ecc. Osserviamo che,
per una matrice A, il limite
lim An = 0
n→∞
Questa condizione equivale a dire che tutte le colonne (e righe) di An tendono a zero. Quindi,
An → 0 come matrice, ossia ogni elemento di An tende a zero. Molto spesso non è pratico osservare
direttamente la matrice An per grandi n. Questa equivalenza ci dice che possiamo invece guardare
l’effetto di An su qualsiasi vettore x: se ogni volta che applichiamo A più e più volte il risultato si
avvicina a 0, allora anche la matrice stessa tende a 0.
Utilizzando la disuguaglianza triangolare e la proprietà submoltiplicativa delle norme:
||A||n2 = ρn (A),
dove ρ(A) è il raggio spettrale di A, ossia il suo autovalore di modulo massimo. Nel caso gen-
erale, ||A||2 coincide con il valore singolare massimo di A. Si ottiene quindi la condizione di
convergenza:
lim An = 0 ⇐⇒ ρ(A) < 1.
n→∞
xk+1 = Axk = Ak x0 .
Se Ak non tende a zero, l’algoritmo non converge. Per questo, la convergenza si dimostra imponendo
che il raggio spettrale della matrice di propagazione sia minore di 1.
20
Consideriamo ora una serie di potenze scalare:
∞
X
an xn ,
n=1
che converge per |x| < R dove R è il raggio di convergenza. Si può dimostrare che, se una matrice
A ha raggio spettrale (nel caso simmetrico posi- tivo definito) o valore singolare massimo (nel caso
generale) minore di R, allora anche la serie matriciale
∞
X
an An
n=0
converge. Questa proprietà è cruciale, poiché rappresenta la forma generale di una serie di Taylor.
Pertanto, se abbiamo una funzione analitica f (x) (come sin x, cos x, ex , log x, o un polinomio
razionale), possiamo definire la funzione di una matrice f (A) tramite:
∞
X
f (A) = an An ,
n=0
a condizione che la serie converga, cioè che ρ(A) sia minore del raggio di convergenza della serie. In
questo modo, possiamo calcolare funzioni di matrici utilizzando solo somme e prodotti di matrici,
senza dover ricorrere a decomposizioni ai valori singolari (che, pur essendo il metodo più accurato,
è anche più costoso computazionalmente). Esempio. Prendiamo come esempio la funzione:
∞
X An
eA =
n=0
n!
21
• n è l’esponente;
• b è la base del sistema di numerazione.
Le basi più comuni sono:
• b = 10: base decimale (umana);
• b = 2: base binaria (informatica);
• b = 16: base esadecimale (usata in programmazione e grafica).
Nel sistema binario, un numero in virgola mobile si rappresenta come:
(1.d1 d2 ...dp )2 × 2n ,
dove la parte (1.d1 d2 ...dp ) è detta mantissa. Il fatto più importante (e ”scioccante”) di questa
rappresentazione è che esistono un numero massimo e un numero minimo rappresentabili. L’insieme
dei numeri rappresentabili non è infinito. Ad esempio:
• il numero minimo ha forma 1.000...0 × 2nmin ;
• il numero massimo ha forma 1.111...1 × 2nmax .
Ma il problema non è tanto il loro valore assoluto, quanto il fatto che tali limiti esistano. Un
fatto ancora più rilevante è che esiste una differenza minima apprezzabile tra due numeri
rappresentabili. Se due numeri differiscono per meno di una certa quantità (ad esempio xmin
o 2xmin ), il calcolatore li considera uguali. In altre parole, la rappresentazione finita introduce
inevitabilmente errori di troncamento: a un certo punto si perdono cifre significative. Ad
esempio, se conosciamo solo le prime 8 o 16 cifre di π, tutto il resto ”non esiste” per il computer.
Questo comporta che, sottraendo un numero molto grande da uno molto piccolo, si possa perdere
completamente precisione, poiché le cifre del numero piccolo vengono ”cancellate”.
Le due rappresentazioni più diffuse per numeri reali nei calcolatori sono:
• 32-bit (binary32 o float), con circa 8 cifre decimali significative;
• 64-bit (binary64 o double), con circa 16 cifre significative.
Nella rappresentazione a 64 bit, disponiamo di circa 16 cifre decimali significative. Questo significa
che, se in un calcolatore sommiamo 1 e 109 , il risultato sarà ancora 109 . A prima vista potremmo
pensare che non sia un problema: dopotutto 1 + 109 ≈ 109 .
Tuttavia, disporre solo di poche cifre significative può essere molto pericoloso: anche errori piccol-
issimi possono avere effetti catastrofici, poiché gli errori numerici si propagano sempre. Un
errore iniziale dell’ordine di 10−8 può facilmente crescere, dopo alcune operazioni, fino a 10−2 , ossia
un errore dell’1%.
Nel formato binary32 (single precision) utilizziamo 32 bit:
• 1 bit per il segno;
• 8 bit per l’esponente;
• 23 bit per la mantissa (più un bit implicito).
22
L’esponente può variare approssimativamente tra −126 e 127.
Nel formato binary64 (double precision) abbiamo invece:
• 1 bit per il segno;
• 11 bit per l’esponente;
• 52 bit per la mantissa.
Ogni numero reale x viene approssimato dal suo equivalente in floating point. Poiché il numero
di cifre è finito, esiste una costante positiva, detta epsilon di macchina (ϵmach ), tale che tutti i
numeri che differiscono da x di meno di ϵmach vengono rappresentati nello stesso modo.
(
2−23 ≈ 1.19 × 10−7 , single precision,
ϵmach ≈
2−52 ≈ 2.22 × 10−16 , double precision.
Applicare una matrice A a un vettore può amplificare l’errore di un fattore proporzionale a ||A||. Se
||A|| ≈ 1012 , in single precision non ha senso effettuare il calcolo, poiché l’errore diventa maggiore
della precisione disponibile.
Sviluppando:
xϵ = (A + ϵB)−1 (b + ϵc)
xϵ = A−1 (I + ϵA−1 B)−1 (b + ϵc)
Il termine scritto sopra può essere sviluppato in serie:
Allora:
xϵ = A−1 (I + ϵA−1 B)−1 (b + ϵc) ≈ (I − ϵA−1 B)A−1 (b + ϵc) + o(ϵ)
Ne segue che:
||δx|| = ||xϵ − x|| ≤ ϵ||A−1 (c − Bx)|| ≤ ϵ||A−1 c|| + ϵ||A−1 Bx||
Anche se A è invertibile in teoria, A + ϵB potrebbe non esserlo in aritmetica floating point. Se A ha
autovalori molto piccoli, ad esempio dell’ordine di 10−12 , in double precision è ancora invertibile, ma
in single precision tali autovalori vengono ”annullati”. Il calcolatore, quindi, interpreta la matrice
come non invertibile, pur essendo teoricamente invertibile.
23
Vogliamo ora stimare la norma di ||δx||, ovvero l’errore nella soluzione del sistema lineare. Nel
primo termine, moltiplichiamo e dividiamo per ||x||:
||A−1 c||||x||
ϵ||A−1 c|| = ϵ
||x||
Nel secondo termine, moltiplichiamo e dividiamo per ||A||:
ϵ||A||||A−1 Bx||
ϵ||A−1 Bx|| =
||A||
κ(A) ≤ 108 .
Oltre questa soglia, la precisione numerica è persa: numeri anche molto diversi diventano
indistinguibili.
24
• In precisione doppia (64 bit) abbiamo circa 16 cifre significative, quindi possiamo gestire
matrici fino a κ(A) ≈ 1016 .
Possiamo quindi dire in parole semplici che:
||δx|| ||δA|| ||δb||
≈ κ(A) +
||x|| ||A|| ||b||
L’errore relativo sulla soluzione è uguale, in prima approssimazione, all’errore relativo sui dati e
sulla matrice, moltiplicato per il numero di condizionamento.
Facciamo quindi qualche osservazione sulla propagazione degli errori
• Nel calcolo di autovalori e autovettori, ad esempio nel metodo di Gram- Schmidt, anche matrici
ben condizionate possono generare numeri di condizionamento dell’ordine di 109 o 1012 . Per
questo, conoscere e controllare il condizionamento è essenziale per la stabilità numerica.
Casi particolari
L’ineguaglianza generale può essere specializzata in due casi:
1. Errore solo sui dati (δb ̸= 0, δA = 0):
||δx|| ||δb||
≤ κ(A)
||x|| ||b||
||δx|| ||δA||
≤ κ(A)
||x|| ||A||
κ(A−1 ) = κ(A)
25
Il numero di condizionamento più rilevante è quello rispetto alla norma 2:
σmax (A)
κ2 (A) = ||A||2 ||A−1 ||2 = .
σmin (A)
dove σmax (A) e σmin (A) sono rispettivamente i massimi e minimi valori singolari di A. Se A è
simmetrica, le singolarità coincidono con i suoi autovalori (in valore assoluto), quindi:
|λmax (A)|
κ2 (A) =
|λmin (A)|
||A|| = ||A−1 || = 1
da cui segue:
κ(A) = 1
Come anticipato all’inizio della lezione, le matrici ortogonali (o unitarie, nel caso complesso) non
propagano errori. Ciò non significa che li eliminino, ma semplicemente che non li amplificano.
In altre parole, moltiplicare un vettore per una matrice ortogonale non peggiora l’errore:
Nella pratica comunque , capita spesso di lavorare con matrici rettangolari. Per estendere la
definizione di condizionamento, utilizziamo la Singular Value Decomposition (SVD). Defini-
amo quindi:
κ(A) = ||A|| · ||A+ ||,
dove A+ è la pseudoinversa. Anche per matrici quadrate non invertibili, questa definizione per-
mette di stimare una soluzione di minimi quadrati e di norma minima.
Sia A simmetrica e definita positiva. Osservando il suo spettro (o i suoi valori singolari), possi-
amo dire che:
• più gli autovalori sono vicini tra loro, più la matrice è ben condizionata;
• più sono dispersi, più la matrice è mal condizionata.
Esempio: Se λmax = 10 e λmin = 1 ⇒ κ(A) = 10. Ma se λmax = 100 e λmin = 10−8 , allora
κ(A) = 1010 e la matrice è estremamente mal condizionata.
26
Esempio: basi gaussiane e dipendenza lineare
Questo fenomeno si verifica spesso in chimica quantistica, quando si usano basi gaussiane. Pren-
diamo una matrice A che per colonne le gaussiane del sistema in esame. Le basi gaussiane sono
spesso formate da funzioni linearmente indipendenti tra loro ma spesso non ortonormali. Tuttavia
può succedere che, se due gaussiane della base sono prese su centri molto vicini, venga a mancare
l’indipenza lineare. Questo è noto come problema dell’ over-completness.
Supponiamo quindi che le colonne della matrice A siano quasi linearmente dipendenti; in tal
caso la dimensione del kernel sarà ̸= 0 e quindi ci sarà un vettore non nullo tale che:
Ax = 0
AT Ax = Sx = AT 0 = 0 = λx
ma se questa deve essere valida per x ̸= 0 allora deve essere che avremo un autovalore che si avvicina
sempre più allo zero tanto più due colonne sono linearmente dipendenti, e quindi:
λmax
κ(S) = ≫1
λmin
La matrice risulta dunque mal condizionata. Questo è uno dei motivi per cui le simulazioni con
basi gaussiane possono presentare gravi problemi di convergenza.
La lezione fondamentale è che dobbiamo sempre valutare se l’errore (inevitabile) può compromettere
i risultati. Un algoritmo può essere estremamente efficiente, ma se produce una soluzione con
un errore del 2000%, non è utile. La stabilità numerica è quindi un aspetto centrale di ogni
implementazione. Garantire la stabilità significa progettare algoritmi che non amplifichino gli
errori dovuti alla rappresentazione finita dei numeri o alla propagazione degli errori di arrotonda-
mento.
27
• una matrice 1000 × 1000 richiede circa 8 megabyte di memoria;
• una matrice 10 000 × 10 000 richiede circa 800 megabyte;
• una matrice 100 000 × 100 000 richiede circa 80 gigabyte.
Il problema principale non è tanto l’entità assoluta della memoria necessaria, quanto la crescita
quadratica della stessa con la dimensione della matrice. Tale crescita rende rapidamente imprati-
cabile la memorizzazione esplicita di matrici di grandi dimensioni. In chimica computazionale, ad
esempio, nel calcolo delle equazioni di risposta lineare per proprietà come le schermature NMR a
livello DFT, non è raro dover risolvere sistemi lineari di dimensione 100 000 × 100 000. Se si lavora
a livello CASSCF, tali dimensioni possono facilmente salire a 100 milioni per 100 milioni, e ancora
di più nei metodi di tipo Coupled Cluster. In ogni caso i metodi di cui ci occuperemo in questa
parte sono detti metodi densi, poiché assumono che la matrice sia piena, cioè che la maggior
parte dei suoi elementi siano non nulli. Non facendo dunque distinzione tra elementi nulli e non
nulli, l’intera matrice viene memorizzata e trattata come tale. Questi metodi, perciò, sono
applicabili solo a problemi di dimensioni moderate.
Il costo computazionale di questi algoritmi cresce con il cubo della dimensione della matrice. In
altre parole, risolvere un sistema lineare di dimensione n × n richiede un numero di operazioni pro-
porzionale a n3 . Ad esempio, per un sistema di dimensione 1000 × 1000, sono necessarie dell’ordine
di 109 (un miliardo) di operazioni in aritmetica a virgola mobile. A questo punto, può essere utile
chiedersi: quante operazioni in virgola mobile può eseguire un comune laptop in un
secondo? La frequenza di clock di un processore, misurata in gigahertz (GHz), rappresenta la
velocità con cui esso può eseguire istruzioni. Un processore da alcuni gigahertz è dunque in grado
di gestire diversi miliardi di istruzioni al secondo. Tuttavia, ogni istruzione può operare su più dati
contemporaneamente. Ad esempio, un’operazione di somma vettoriale può sommare più elementi
di due vettori in parallelo, anziché uno alla volta. La quantità di dati che possono essere elaborati
simultaneamente dipende dalle estensioni vettoriali supportate dalla CPU. Tra queste troviamo le
famiglie di istruzioni: SSE, SSE2, SSE3, SSE4, SSE4.1, AVX, AVX2, AVX-512. Queste estensioni
indicano quante operazioni vettoriali possono essere eseguite in parallelo, ovvero quanti bit possono
essere processati simultaneamente. Ad esempio, le istruzioni AVX-512 consentono di processare
blocchi di 512 bit alla volta, permettendo quindi di sommare fino a 8 numeri con una singola
istruzione di clock. Questo tipo di parallelismo interno è ciò che consente ai moderni processori di
raggiungere elevate prestazioni nel calcolo scientifico. Naturalmente, quanto discusso in precedenza
rappresenta un caso ideale, valido solo se si riesce a scrivere un algoritmo capace di sfruttare al
100% la potenza di calcolo della CPU. Nella pratica, questo non avviene mai. Le operazioni più effi-
cienti che si possano eseguire su un computer, come ad esempio la moltiplicazione matrice-matrice,
raggiungono in genere un’efficienza dell’ordine dell’80%. Tuttavia, se si scrive un codice in proprio,
senza un’ottimizzazione particolarmente sofisticata, l’efficienza effettiva è molto più bassa:
• un programmatore medio può aspettarsi circa l’1% di efficienza:
• un buon programmatore può raggiungere il 10%;
• un programmatore eccellente può spingersi fino al 20%.
Consideriamo ora l’impatto pratico di tali valori. Un sistema lineare di dimensione 1000 × 1000
richiede dell’ordine di 109 operazioni in virgola mobile. Su un normale laptop, questo equivale a
un tempo di calcolo di circa un quarto di secondo. Poiché il costo scala come n3 un sistema di
28
dimensione 10 000 × 10 000 richiede circa 1012 operazioni, ovvero dell’ordine di 103 secondi (circa 15
minuti). Un sistema di dimensione 100 000 × 100 000 comporta circa 1015 operazioni, cioè milioni
di secondi, equivalenti a diversi giorni di calcolo.
MA = U
dove U è una matrice triangolare superiore. A quel punto, il sistema lineare iniziale
Ax = b
Ax = b
può essere risolto partendo dall’ultima equazione e risalendo verso la prima. Per l’ultima riga del
sistema, si ha una sola incognita:
ann xn = bn
da cui si ottiene immediatamente:
bn
xn =
ann
29
Salendo di una riga, si ha:
an−1,n−1 xn−1 + an−1,n xn = bn−1
e poiché xn è già noto, si può calcolare facilmente xn−1 . Procedendo in questo modo, si ottiene la
soluzione completa per sostituzione all’indietro.
Vediamo quindi quante operazioni sono necessarie. Per l’ultima equazione si eseguono una divisione
e una moltiplicazione; per la penultima, una divisione e due moltiplicazioni; poi tre, quattro, e cosı̀
via. In totale, il numero di operazioni (moltiplicazioni e somme) è approssimativamente:
n(n + 1)
2
più n divisioni. Pertanto, la complessità della risoluzione di un sistema triangolare è proporzionale
a n2 , cioè quadratica, molto inferiore rispetto alla complessità cubica dei metodi generali.
M Ax = M b
assicurandoci che ogni trasformazione intermedia sia ottenuta tramite una matrice invertibile.
Questa è la strategia fondamentale alla base dell’eliminazione di Gauss.
Partiamo dicendo che per ogni colonna k della matrice A, chiamiamo pivot l’elemento diagonale
appartenente a quella colonna, ovvero akk . Il pivot sarà il riferimento rispetto al quale elimineremo
tutti gli altri elementi della stessa colonna.
Posto A(1) = A consideriamo la prima colonna di A(1) . Indichiamo con p1 il pivot, ossia l’elemento
a11 , e distinguiamo tre casi principali:
1. Caso 1: Tutta la prima colonna è nulla. In questo caso, la matrice A ha una colonna
interamente nulla, e dunque non può essere invertibile. Non è quindi possibile procedere con
l’eliminazione; poniamo semplicemente A(2) = A(1) e passiamo oltre.
2. Caso 2: Il pivot p1 = a11 è nullo, ma nella colonna esiste almeno un elemento
non nullo. In questo caso, applichiamo una permutazione di righe per portare un elemento
non nullo in posizione (1, 1), rendendo quindi il pivot diverso da zero. Questa operazione si
effettua mediante una matrice di permutazione P (1) . In algebra lineare, permutare due righe
di una matrice A equivale a moltiplicare A a sinistra per una matrice P ottenuta a partire
30
dall’identità scambiando le righe corrispondenti. Ad esempio, se vogliamo scambiare la prima
e la terza riga, prendiamo la matrice identità e scambiamo le stesse righe:
0 0 1
P (1) = 0 1 0
1 0 0
31
3. Se akk ̸= 0, poniamo P (k) = I e proseguiamo con l’eliminazione.
La matrice di trasformazione E (k) viene costruita analogamente al caso k = 1:
k−1
I 0
E (k) =
0 M̃ (k)
dove M̃ (k) è una matrice identità modificata nei coefficienti al di sotto della diagonale nella colonna
k, in modo che:
aik
(M̃ (k) )ik = − per i > k
akk
Tutti gli altri elementi di M̃ (k) coincidono con quelli dell’identità. Dopo l’applicazione della trasfor-
mazione, otteniamo:
A(k+1) = E (k) P (k) A(k)
Dopo k iterazioni, le prime k colonne di A(k) risultano già nella forma desiderata, cioè triangolare
superiore. Ogni passo è costruito tramite moltiplicazione per matrici invertibili (E (k) e P (k) ),
dunque l’intera trasformazione globale:
è ottenuta da una matrice invertibile M = E (n−1) P (n−1) ...E (1) P (1) . Al passo n-esimo, otteniamo
una matrice A(n) che, per costruzione, risulta triangolare superiore. Abbiamo quindi costruito
l’eliminazione di Gauss.
Ora ciò che resta da verificare è se la matrice M, ottenuta come prodotto di tutte le matrici di
trasformazione e di permutazione utilizzate nei vari passaggi, sia effettivamente invertibile. Consid-
eriamo innanzitutto le matrici di permutazione. Queste si ottengono a partire dalla matrice identità,
semplicemente scambiando due righe (o due colonne). Scambiare due righe (o colonne) cambia il
segno del determinante ma non il suo valore assoluto. Poiché la matrice identità ha determinante
pari a 1, segue che ogni matrice di permutazione P (i) ha
det(P (i) ) = ±1
Passiamo ora alle matrici di trasformazione E (i) . Ciascuna di esse è una matrice triangolare inferiore
con tutti 1 sulla diagonale principale. Il determinante di una matrice di questo tipo è il prodotto
degli elementi diagonali, quindi:
det(E (i) ) = 1
Il determinante della matrice complessiva M è dunque:
Y Y
det(M ) = det(P (i) ) · det(E (j) )
i j
ossia un numero ±1 certamente diverso da zero. Ne consegue che M è invertibile. Abbiamo quindi
dimostrato che esiste una matrice invertibile M tale che:
MA = U
dove U è una matrice triangolare superiore. In altre parole, l’eliminazione di Gauss esiste ed è
sempre realizzabile.
32
Quando si applica il metodo di Gauss a un sistema lineare Ax = b, non si opera solo sulla matrice
A ma anche sul vettore dei termini noti b. Ogni volta che si applica una permutazione o una
trasformazione, si deve applicare la stessa operazione anche a b:
M Ax = M b
U x = b′
Alla fine delle operazioni si ottiene quindi un sistema triangolare superiore, che può essere risolto
facilmente per sostituzione all’indietro.
Esempio pratico
Consideriamo un esempio semplice. Partiamo dall’elemento (1, 1) e supponiamo che il pivot p1 sia
diverso da zero, quindi non è necessaria alcuna permutazione:
P (1) = I
La matrice di eliminazione E (1) serve ad annullare tutti gli elementi sotto il pivot nella prima
colonna. I valori inseriti nella prima colonna di E (1) sono − api11 . Nel nostro caso: (−2, −1). Appli-
cando questa trasformazione, otteniamo:
1 2 1
A(1) = 0 1 0
0 0 1
La matrice risulta già triangolare superiore. Anche il vettore dei termini noti deve essere trasfor-
mato:
1
b′ = E (1) b = −2
−1
Il sistema triangolare risultante è:
z = 1
y = −2
x+y+z =1
Da cui si ottiene:
x=2 y = −2 z=1
Costo computazionale
Abbiamo ora un metodo per trasformare una matrice in forma triangolare superiore e risolvere
un sistema lineare in modo efficiente. Ma quanto è costoso questo procedimento in termini di
operazioni? Per la prima colonna occorrono: (n − 1) moltiplicazioni, (n − 1) somme, (n − 1)
divisioni. Per la seconda colonna: (n − 2) operazioni di ciascun tipo, e cosı̀ via. Sommando:
n(n − 1)
1 + 2 + · · · + (n − 1) =
2
33
Poiché queste operazioni vanno effettuate per ciascuna riga, il numero totale cresce proporzional-
mente a:
n(n − 1) 1
× n ≈ n3
2 2
Il costo computazionale dell’eliminazione di Gauss è dunque dell’ordine di:
O(n3 )
1.12.2 Fattorizzazione LU
Nella pratica, l’eliminazione di Gauss non viene quasi mai utilizzata direttamente per risolvere sis-
temi lineari di grandi dimensioni, poiché esistono algoritmi più efficienti che permettono di ottenere
lo stesso risultato a un costo computazionale inferiore. In particolare, l’operazione più costosa nel
metodo di Gauss è la fattorizzazione della matrice. Alcuni algoritmi moderni riescono a eseguire
la stessa operazione con un costo pari a una frazione di quello originale: rispettivamente 31 n3 o ad-
dirittura 16 n3 operazioni. Questi due algoritmi sono oggi i più comunemente utilizzati per risolvere
sistemi lineari.
Per mostrare l’idea di questi algoritmi consideriamo il caso in cui la matrice A sia invertibile.
L’idea di base nasce proprio dall’eliminazione di Gauss: vogliamo fattorizzare la matrice A come il
prodotto di due matrici triangolari, una inferiore e una superiore:
A = LU
dove:
• L è una matrice triangolare inferiore con diagonale unitaria (Lii = 1);
• U è una matrice triangolare superiore.
Questa decomposizione prende il nome di fattorizzazione LU. Si può dimostrare che se tutti i minori
principali di A (cioè le sottomatrici ottenute eliminando le prime k righe e colonne per k = 1, 2, ..., n)
sono invertibili, allora è possibile effettuare la decomposizione senza alcun pivoting, ossia senza dover
scambiare righe durante il processo di eliminazione. In questo caso, il procedimento è stabile e il
sistema può essere risolto direttamente.
Nel caso in cui la matrice non soddisfi la condizione precedente, è comunque possibile general-
izzare la decomposizione LU introducendo una matrice di permutazione P tale che:
P A = LU
Questa variante è detta decomposizione LU con pivoting. Il pivoting non solo garantisce l’esistenza
della decomposizione anche per matrici non ben condizionate, ma migliora sensibilmente la stabilità
numerica del metodo. In pratica, anche quando non è strettamente necessario, si preferisce sem-
pre applicare il pivoting, scegliendo a ogni passo come pivot l’elemento di modulo massimo della
colonna. Questo procedimento prende il nome di pivoting parziale. Il pivoting completo (su righe e
colonne) è possibile ma raramente utilizzato, poiché comporta costi aggiuntivi nella gestione delle
permutazioni. Nel nostro caso, per semplicità, ignoreremo il pivoting e deriveremo l’algoritmo LU
senza pivoting, consapevoli che in implementazioni reali esso è sempre incluso.
34
Consideriamo ora il prodotto esplicito:
A = LU
Poiché L è triangolare inferiore, risulta Lik = 0 per k > i, mentre per U vale Ukj = 0 per k > j. Di
conseguenza, la somma effettiva si riduce a:
min(i,j)
X
Aij = Lik Ukj
k=1
Questa osservazione è fondamentale, perché implica che il calcolo di ciascun elemento della fat-
torizzazione può essere effettuato utilizzando solo elementi già noti. Questa struttura suggerisce
un ordine naturale con cui procedere nel calcolo degli elementi di L e U. L’algoritmo LU è
dunque un insieme finito di istruzioni elementari (moltiplicazioni, somme, divisioni)
che permettono di costruire progressivamente le due matrici triangolari.
Prima colonna: Consideriamo ora la prima colonna di A (cioè j = 1). Anche qui la somma va
da k = 1 a 1:
Ai1 = Li1 U11
Poiché U11 è noto e L11 = 1 si ottiene
Ai1 Ai1
Li1 = =
U11 A11
Questa relazione è valida per tutti gli elementi della prima colonna di L. Osserviamo che L coincide
(a meno del segno) con la prima colonna della matrice costruita nel processo di eliminazione di
Gauss.
poiché L11 = 1.
Abbiamo quindi determinato la prima colonna di L e la prima riga di U.
35
Seconda colonna: Consideriamo ora A22 . La somma si estende da k = 1 a 2:
Da cui si deduce:
• Se i < j (cioè nella parte superiore della diagonale), si calcola:
i−1
X
Uij = Aij − Lik Ukj
k=1
Tutti i termini Lik Ukj sono già noti dai passi precedenti, quindi non serve risolvere alcun
sistema: è una semplice somma e sottrazione.
• Se i = j, otteniamo:
i−1
X
Uii = Aii − Lik Uki
k=1
Anche qui tutti i termini del lato destro sono già noti, tranne Ujj , che è non nullo se la matrice
A è invertibile.
Si noti che, sebbene Ujj non sia mai nullo per una matrice invertibile, esso può risultare nu-
mericamente piccolo, motivo per cui in pratica si effettua sempre il pivoting, anche quando non
strettamente necessario, per garantire stabilità numerica.
36
• il processo continua col calcolo alternato di righe e colonne, fino a completare entrambe le
matrici L e U.
Ogni elemento calcolato dipende solo da elementi già determinati, grazie alla struttura triango-
lare delle due matrici. Infatti, per calcolare un elemento Aij servono solo: Li1 , Li2 , ..., Lij e
U1j , U2j , ..., Uij e tutti questi valori sono disponibili dai passi precedenti, poiché ogni moltiplicazione
oltre min(i, j) coinvolgerebbe solo zeri.
Esempio illustrativo: Consideriamo ora un esempio concreto per verificare che il procedimento
funziona effettivamente. Sia la matrice
2 ∗ ∗
A = ∗ ∗ ∗
∗ ∗ ∗
L’elemento A11 deriva dal prodotto della prima riga di L con la prima colonna di U, e i soli termini
non nulli sono quelli sulla diagonale:
A11 = L11 U11
Analogamente, A21 dipende solo da L21 e U11 , mentre A22 coinvolge anche L22 e U22 , e cosı̀
via. In altri termini, L31 risulta essere pari a 3. Passiamo ora alla seconda colonna. L’elemento
corrispondente si ottiene moltiplicando la seconda riga di L per la seconda colonna di U. In questo
caso abbiamo:
A12 = L11 U12 + 0 + 0 = 3
mentre per l’elemento successivo,
da cui ricaviamo il valore di U22 . Infine, per il termine A32 (terza riga, seconda colonna), il termine
incognito è L32 :
A32 = L31 U12 + L32 U22 + 0
L32 = 9
Terza colonna: Procediamo ora con l’ultima colonna. Per la prima riga, terza colonna, otteniamo
semplicemente
A13 = U13
Per la seconda riga, terza colonna:
Si nota dunque che, ad ogni passo, utilizziamo esclusivamente informazioni già calcolate nei passaggi
precedenti. Verificando il prodotto LU, si ritrova esattamente la matrice iniziale A.
37
Risoluzione del sistema lineare con LU Perché la decomposizione LU è cosı̀ utile? La risposta
risiede nella semplicità del metodo di risoluzione dei sistemi lineari che ne deriva. Infatti, una volta
ottenuta la fattorizzazione A = LU , si risolve il sistema lineare
Ax = b
Validità generale: La decomposizione LU esiste per ogni matrice quadrata invertibile (a meno
di pivoting), indipendentemente dalla sua struttura. Non abbiamo fatto ipotesi specifiche, se non
l’invertibilità di A. Tuttavia, esiste una classe di matrici particolarmente importante: le matrici
simmetriche definite positive (SPD).
Perché le matrici simmetriche definite positive sono cosı̀ comuni? Le matrici simmetriche
compaiono naturalmente in molti problemi fisici: per esempio, nelle rappresentazioni matriciali di
operatori in meccanica quantistica o nei problemi di ottimizzazione, dove il sistema di Newton
assume la forma:
H∆x = −∇f
con H matrice hessiana, che è simmetrica per definizione. Se ci troviamo in un punto di minimo,
H è anche definita positiva. In altri contesti, come nei metodi ai minimi quadrati, moltiplicare un
sistema rettangolare per la sua trasposta genera una matrice simmetrica e definita positiva:
AT A
A = LLT
38
ossia la matrice U coincide con la trasposta della matrice L. Questa decomposizione prende il nome
da André-Louis Cholesky, ufficiale dell’esercito francese ucciso durante la Prima Guerra Mondiale,
che ideò il metodo per risolvere sistemi lineari per scopi militari. I suoi appunti furono scoperti e
pubblicati postumi da un suo collega.
Algoritmo di Cholesky: L’algoritmo può essere descritto con un ciclo (do-loop) che scorre sugli
indici j da 1 a n. Per ciascun passo, si calcolano prima gli elementi diagonali, poi quelli al di sotto
della diagonale.
Da cui: v
u
u j−1
X
Ljj = tAjj − L2jk
k=1
Se la matrice A è definita positiva, il termine sotto radice è sempre positivo e il calcolo è ben
definito. Un test pratico di positività definita consiste proprio nel tentare una decomposizione
di Cholesky: se il valore sotto radice risulta negativo, la matrice non è definita positiva.
2. Calcolo degli elementi sotto la diagonale (i > j):
j−1
!
1 X
Lij = Aij − Lik Ljk
Ljj
k=1
Anche in questo caso, tutti i termini del lato destro sono già noti dai passi precedenti.
39
Ulteriori Considerazioni sull’Algoritmo di Cholesky e Uso Pratico
Il procedimento di decomposizione di Cholesky può essere implementato con tre cicli annidati: per
ogni j = 1, ..., n, per ogni i = j + 1, ..., n, per ogni k = 1, ..., j − 1. Questa struttura genera
la complessità O(n3 ), ma poiché i loop non percorrono tutti gli elementi della matrice (solo la
metà inferiore o superiore), il fattore complessivo risulta essere 16 rispetto alla decomposizione LU
completa.
Esercizio pratico: Si consiglia di prendere una matrice 3 × 3, moltiplicarla per la sua trasposta
in modo da ottenere una matrice simmetrica semidefinita positiva, e calcolare manualmente la
decomposizione di Cholesky. L’obiettivo è comprendere che l’algoritmo si basa su una sequenza
chiara e finita di istruzioni semplici, esattamente ciò che un computer può eseguire.
Uso di librerie ottimizzate: Nella pratica non conviene mai implementare da zero operazioni
di algebra lineare di base, come:
• somma di vettori,
• prodotto matrice-vettore,
• prodotto matrice-matrice,
40
come ordine di grandezza. Questo rappresenta una differenza enorme: ad esempio, per una matrice
106 × 106 ,
n2 = 1012 operazioni
n3 = 1018 operazioni
corrispondente a un fattore di 106 in termini di numero di operazioni. Per dare un’idea, un milione
di secondi corrisponde a molti giorni, considerando che un giorno ha 86,400 secondi. Tra i metodi
iterativi più noti troviamo Jacobi e Gauss-Seidel, che applicano lo stesso principio: approssimazioni
successive fino al raggiungimento della precisione desiderata.
Ax = b
dove A è una matrice n × n, x è il vettore delle incognite e b è il vettore dei termini noti.
Il principio fondamentale consiste nel riscrivere il sistema Ax = b nella forma di punto fisso o
forma di iterazione:
x = Px + q
dove P è la matrice di iterazione e q è un vettore noto.
Data una stima iniziale arbitraria x(0) , l’iterazione procede come segue:
x(k+1) = P x(k) + q
per k = 0, 1, 2, . . .
Per ottenere la forma di punto fisso, la matrice A viene decomposta (splitting) come:
A=M −N
dove M è una matrice non singolare e facilmente invertibile. Sostituendo nel sistema Ax = b e
moltiplicando per M −1 si ottiene:
Mx = Nx + b
xk = M −1 N xk−1 + M −1 b
da cui si identificano P e q:
P = M −1 N e q = M −1 b
Si dice che un metodo iterativo stazionario converge alla soluzione esatta x (indipendentemente
dalla scelta di x(0) ) se:
lim ek = 0
k→∞
41
• Residui: rk = Axk − b = Aek
• Incrementi: δk = xk − xk−1 = ek − ek−1
I residui vanno esattamente come l’errore e quindi da un punto di vista pratico possiamo
determinare la convergenza come il punto in cui si annullano i residui; se cerchiamo
invece il punto in cui si annullano gli incrementi (sebbene questi siano molto più facili da calcolare)
troviamo solo il punto in cui l’errore tra due iterazioni è molto simile che, in linea di principio non
è il punto di convergenza convergenza, tuttavia praticamente spesso i due punti coincidono. Una
buona regola comunque è quella di determinare la convergenza con i residui qualora
questo procedimento sia computazionalmente possibile in tempo abbastanza breve.
Si può dimostrare comunque che l’algoritmo converge solo se il raggio spettrale ρ(P ) della matrice
di iterazione P è strettamente minore di 1. Per dimostrarlo prendiamo la definizione di errore
all’iterazione k:
e(k) = x(k) − x
e sostituendo xk = M −1 b + M −1 N xk−1 otteniamo:
Ciò è equivalente a richiedere che la norma del vettore errore tenda a zero:
lim ∥e(k) ∥ = 0
k→∞
Dalla relazione e(k) = P k e(0) , per un’arbitraria norma vettoriale e la corrispondente norma matri-
ciale indotta, abbiamo:
∥e(k) ∥ = ∥P k e(0) ∥ ≤ ∥P k ∥∥e(0) ∥
La convergenza è garantita se:
lim ∥P k ∥ = 0
k→∞
che abbiamo già visto essere equivalente alla condizione che il raggio spettrale di P sia strettamente
minore di 1:
lim ∥P k ∥ = 0 ⇐⇒ ρ(P ) < 1
k→∞
42
Finora comunque abbiamo considerato come unica fonte di errore quella derivata dal metodo, ma
sappiamo che per ogni calcolo abbiamo un errore di machine precision; riscriviamo l’equazione
aggiungendo questo contributo:
xk = M −1 b + M −1 N xk−1 + ϵk
la domanda è se il sistema può ancora convergere. Per rispondere a questa domanda riscriviamo
l’errore:
ek+1 = xk+1 − x = M −1 N ek + ϵk
che in forma generale diventa:
k
X
ek+1 = P k+1 e0 + (P )i ϵk−i−1
i=0
e vogliamo che questo si annulli per k → ∞; come al solito deve essere che, rispetto ad una norma
generica, tenda a 0 l’errore:
k
X
k+1
||ek+1 ||s ≤ ||P e0 ||s + || (P )i ϵk−i−1 ||s
i=0
ma siccome siamo in uno spazio di dimensione finita ed ogni norma è equivalente allora il primo
termine è equivalente alla serie di potenza di matrice in norma 2 e si annullerà se ρ(P ) < 1.
Supponendo che questo si verifichi ci basta studiare l’andamento del secondo termine:
k
X k
X
|| (P )i ϵk−i−1 ||s ≤ || (P )i ||s ||ϵk−i−1 ||s
i=0 i=0
ma se l’errore deve essere ≤ al membro a sinistra allora dovra esserlo anche del membro a destra e
possiamo quindi scrivere:
k
X
||ek+1 ||s ≤ ||P k+1 e0 ||s + || (P )i ||s ||ϵk−i−1 ||s
i=0
a questo punto consideriamo che, siccome il numero di operazioni che facciamo è finito, il modulo
dell’errore totale derivato dalla finite machine precision deve essere uguale o minore ad una costante
per la machine precision:
||ϵk || ≤ Cϵ
inoltre notiamo che il termine nella sommatoria può essere scritto come una serie geometrica:
1 − ||P ||k+1
||ek+1 ||s ≤ C ϵ + ||P k+1 ||s ||e0 ||s
1 − ||P ||
43
otteniamo quindi che:
C
lim ||ek ||s = ϵ = C ′ϵ
k→∞ 1 − ||P ||
Che quindi ci dimostra che anche con l’errore di macchina questi metodi convergono a
qualcosa, anche se non è più la soluzione esatta. Una buona rule of thumb è di fidarsi di
una costante che sia all’incirca ”machine precision × numero di condizionamento della matrice”.
Metodo di Jacobi
Il metodo di Jacobi (o metodo delle sostituzioni simultanee) utilizza solo la componente diagonale
D per la matrice M dello splitting: M = D e N = −(E + F ).
(k+1)
Formula per Componenti Il metodo calcola tutte le nuove componenti xi in modo simul-
taneo utilizzando solo i valori del passo precedente k:
n
(k+1) 1 bi −
X (k)
xi = aij xj
aii
j=1,j̸=i
Metodo di Gauss-Seidel
Il metodo di Gauss-Seidel (o metodo delle sostituzioni successive) utilizza la parte triangolare infe-
riore di A per lo splitting: M = D + E e N = −F .
(k+1)
Formula per Componenti Il metodo sfrutta i valori delle componenti xj non appena cal-
colati nello stesso passo di iterazione. Le componenti con indice j < i sono già aggiornate:
i−1 n
(k+1) 1 bi −
X (k+1)
X (k)
xi = aij xj − aij xj
aii j=1 j=i+1
44
Condizione di Convergenza
Entrambi i metodi convergono se e solo se il raggio spettrale della rispettiva matrice di iterazione
è strettamente minore di 1: ρ(P ) < 1. La convergenza è sempre garantita se la matrice A è a
predominanza diagonale in senso stretto. Se invece A è simmetrica e definita positiva, Gauss-
Seidel converge sempre, e spesso anche Jacobi (sebbene non sempre garantito da questa singola
condizione)
Accelerazione di Anderson
Visti alcuni metodi iterativi vediamo un metodo per accelerarne la convergenza. Definiamo la
seguente quantità:
k
X
x̃k = di xi
i=0
P
dove xi è il risultato dell’ i-esima iterazione e di sono dei coefficienti tali che i di = 1. Possiamo
definire dei residui: X X X
r̃k = Ax̃k − b = A di xi − di b = di ri
i i i
facciamo una minimizzazione di lagrange e per convenienza algebrica prendiamo un vincolo λ′ = 2λ:
X k
X
L= ⟨ri , rj ⟩di dj + 2λ( di − 1)
i,j i
dove nel primo caso abbiamo sfruttato che ⟨ri , rj ⟩ = ⟨rj , ri ⟩. Possiamo riscrivere quindi il primo
caso in forma matriciale:
Gd⃗ = −⃗λ
dove ⃗λ è un vettore con solo λ. Possiamo riunire questi due vincoli in un unica equazione matriciale:
45
Il metodo di Anderson non si limita a usare l’output dell’ultima iterazione, ma sfrutta i parametri
già calcolati (i vettori xi e i residui ri della storia recente) per calcolare la k-esima iterazione
x(k) in modo che sia la migliore possibile in un senso ben definito: minimizzare il residuo. La
(k)
”migliore possibile” iterazione accelerata xaccel viene definita come la combinazione
lineare convessa delle iterazioni passate.
L’ultima cosa da chiederci prima di passare al prossimo argomento riguarda la gestione della
convergenza. In definiamo la convergenza come un valore di ”treshold” oltre il quale consideriamo
il sistema risolto:
dove il valore preso per τ ′ è una buona ”rule of thub”. Viene considerata la convergenza di tutte
e due le norme poiché offrono una prospettiva complementare: l’RMSD da una prospettiva sulla
media dei resisui, mentre la norma del massimo sul massimo residuo (errore) che commettiamo;
combinandoli siamo sicuri di fare un errore trascurabile.
46
la seguente partizione per la matrice dei coefficienti:
1
M= I
α
1
N= I−A
α
con α reale e non nullo. E possiamo scrivere ogni iterazione come:
(
x0 = x0
xk+1 = M−1 Nxk + M−1 b = xk + αrk
dove gli rk sono i residui. E’ possibile dimostrare che se la matrice è solo semidefinita (positiva o
negativa) allora il metodo non converge per nessun valore di α; mentre se la matrice è definita
(positiva o negativa) allora abbiamo convergenza per ogni:
2
0<α<
λn
con un αopt che deve essere:
λn − λ1
αopt =
λn + λ1
Il metodo si chiama del gradiente poiché prendiamo in considerazione il seguente funzionale:
n n
1 1 X X
f (x) = ⟨Ax, x⟩ − ⟨b, x⟩ = aij xi xj − bi x i
2 2 i,j=1 i=1
∇f (x) = Ax − b
A questo punto notiamo che se A è definita positiva allora non abbiamo autovalori nulli e ci basta
minimizzare ognuno dei termini quadratici (x̂i − λb̂ii )2 ; esiste quindi un unico vettore x0 che ha come
47
componenti (x̂i = λb̂ii ) che minimizza il funzionale e questo vettore è anche la soluzione unica del
sistema lineare Ax = b. Se A non è definita positiva ma solo semidefinita allora avremo almeno
un autovalore nullo; se per ogni autovalore nullo non abbiamo che b̂i = 0 allora il sistema non ha
soluzione (e quindi minimo) e il funzionale diverge a −∞. Se invece b̂i = 0 per ogni i : λi = 0 (ovvero
b appartiene all’immagine di A) allora il minimo di f (x) è ottenuto in tutti i punti x = Ã−1 b + y
dove Ã−1 è la restrizione di A alla sua immagine (ovvero prendiamo tutti i valori non nulli) e y è
lo span del Ker(A); inoltre nel caso in cui la matrice A sia semidefinita o definita negativa si può
dimostrare che il funzionale diverge comunque.
Possiamo ora dimostrare due importanti teoremi:
⟨Ax0 − b, y⟩ = 0 ∀y ∈ F
• (2): Se A è una matrice simmetrica definita positiva allora la forma funzionale f (x) =
1
2 ⟨Ax, x⟩ − ⟨b, x⟩ ha le seguenti proprietà:
Esplicitando otteniamo:
1 ∗
f (P y) = ⟨P AP y, y⟩ − ⟨P ∗ b, y⟩
2
Ora notiamo che il proiettore ortogonale al massimo annullerà qualche elemento di A ma non ne
cambierà mai il segno. La matrice P∗ AP sarà quindi al massimo semidefinita positiva; tuttavia
se è semidefinita allora potremmo avere che b non appartiene alla sua immagine e quindi avremmo
che f (P y) = −∞, tuttavia deve anche essere che:
48
arrivando ad un assurdo. Deve essere quindi che P ∗ b ∈ Im(P∗ AP) e ci sarà quindi almeno una
soluzione della forma P∗ APx = b. Dobbiamo ora solo mostrare che questa soluzione è unica e lo
faremo per assurdo; diciamo che esistano due soluzioni distinte in F :
P∗ APx1 = P ∗ b = P∗ APx2
P∗ AP(x1 − x2 ) = 0
Ma se questo deve essere nullo (e sappiamo che la matrice non è nulla) dovrà essere nullo il vettore
(x1 − x2 ) e quindi anche il seguente prodotto scalare:
Ma se come abbiamo detto prima la matrice A deve essere non nulla, allora sarà nullo:
P(x1 − x2 ) = 0
Da cui deduciamo che Px1 = Px2 , ovvero le due soluzioni sono esattamente uguali, il che è un
assurdo in base alle nostre assunzioni; deve quindi esistere un’unica soluzione all’interno del generico
sottospazio F . Moltiplicando quindi la forma normale della soluzione per un generico y otteniamo:
P∗ yAPx0 = P ∗ yb
Dove x0 è la soluzione unica in F ; ora siccome sia x0 che y sono ∈ F e usando che l’operatore di
proiezione ortogonale è idempotente otteniamo:
y(Ax0 − b) = 0
Ma di nuovo se questa deve valere per qualsiasi y ∈ F allora Ax − b deve essere un vettore nullo e
quindi il seguente prodotto scalare deve essere egualmente nullo:
⟨Ax0 − b, y⟩ = 0 ∀y ∈ F
Dimostrazione secondo teorema : Sappiamo già che f ha il suo minimo in un solo punto x
che è soluzione di Ax = b; in questo punto x abbiamo che ∇f (x) = 0. Prendiamo quindi un punto
x tale che ∇f (x) ̸= 0 e prendiamo δ = −α(Ax − b). Siccome A è simmetrica allora:
1 1
f (x + δ) = ⟨A(x + δ), (x + δ)⟩ − ⟨b, (x + δ)⟩ = f (x) + ⟨Aδ, δ⟩ + ⟨Ax − b, δ⟩
2 2
Notiamo quindi che il termine ⟨Aδ, δ⟩ è esattamente la norma di A indotta dal vettore δ; Avremo
quindi che ⟨Aδ, δ⟩ ≤ ||A||2 ||δ||2 e, siccome A è simmetrica positiva e reale, allora avremo anche che
||A|| = ρ(A), per cui possiamo scrivere:
ρ(A))
f (x + δ) ≤ ρ(A)⟨δ, δ⟩ − α⟨δ, δ⟩ = ||Ax − b||2 (α2 − α)
2
49
e quindi nel caso in cui 0 < α < 2/ρ(A) abbiamo che f (x + δ) < f (x) come volevasi dimostrare.
Da questi teoremi possiamo costruire un metodo iterativo per minimizzare f . Costruiamo
una sequenza di punti xk tali che la sequenza di f (xk ) è decrescente:
xk+1 = xk − α∇f (xk ) = xk + α(b − Ax)
Ma questo è esattamente il metodo iterativo che avevamo definito all’inizio per la risoluzione del
sistema lineare; abbiamo quindi dimostrato che risolvere un sistema lineare la cui matrice è
simmetrica e definita positiva è equivalente a trovare il minimo in Rn del funzionale
f (x) .
Riprendiamo ora la definizione di sottospazio di Krylov:
Definizione: Sia ⃗r un vettore in Rn , chiamiamo il sottospazio vettoriale di Rn spannato dai k + 1
vettori {⃗r, A⃗r, ...., Ak ⃗r} = Kk (A, ⃗r) il sottospazio di Krylov associato al vettore ⃗r e alla matrice A.
Questi sottospazi formano una sequenza crescente di sottospazi vettoriali e per ogni r0 ̸= 0 esiste
un k0 ∈ {0, 1, ..., n − 1} tale che:
(
dim(Kk ) = k + 1 se 0 ≤ k ≤ k0
dim(Kk ) = k0 + 1 se k0 ≤ k
Fondamentalmente stiamo dicendo che per ogni k la dimensione dei sottospazi aumenta fino ad
arrivare ad un valore k0 tale che Kk0 +1 = Rn . Riprendiamo ora il metodo del gradiente:
(
x0 ∈ Rn scelta iniziale
xk+1 = xk + αk (b − Axk )
E si vede facilmente per induzione che questa formula definisce esattamente lo span {r0 , Ar0 , ...., Ak r0 }
del sottospazio di Krylov relativo. E con una simile induzione si dimostra anche che xk+1 ∈ [x0 +Kk ].
Vogliamo ora dimostare che rk+1 ∈ Kk+1 . Se però xk+1 ∈ [x0 + Kk ] deve poter essere possibile
scrivere:
Xk
xk+1 = x0 + αi Ai r0
i=0
Moltiplicando A e sottraendo per b otteniamo:
k
X
rk+1 = r0 + αi Ai+1 r0
i=0
50
E quindi deve essere che rk+1 ∈ Kk+1 come volevasi dimostrare. Per ora limitiamoci a definire xk+1
nello spazio affine [x0 + Kk ] e, siccome in linea di principio ci sono infiniti modi di scegliere xk+1 ,
per determinarlo in modo univoco enunciamo due criteri:
• Principio di ortogonalizzazione: Scegliamo xk+1 ∈ [x0 + Kk ] tale che rk+1 ⊥ Kk
• Principio di minimizzazione: Scegliamo xk+1 ∈ [x0 + Kk ] che minimizzi in [x0 + Kk ] il
seguente funzionale:
1
f (x) = ⟨Ax, x⟩ − ⟨b, x⟩
2
E, se A è una matrice simmetrica definita positiva, allora per le due definizioni esiste un
vettore unico xk+1 ∈ [x0 + Kk ]. Inoltre tutte e due le definizioni corrispondono allo
stesso algoritmo nel senso che portano allo stesso valore di xk+1 e convergono in al più n
iterazioni alla soluzione vera del sistema lineare Ax = b. Minimizzare il funzionale quadratico
nella sequenza di questi sottospazi vettoriali è per definizione il metodo del gradiente
coniugato. Passiamo quindi a dimostarre che i due principi sono esattamente la stessa cosa e
definiscono in modo univovo xk+1 . Definiamo ∀y ∈ Kk il seguente funzionale:
1
g(y) = f (x0 + y) = ⟨Ay, y⟩ − ⟨r0 , y⟩ + f (x0 )
2
Vediamo che minimizzare f in [x0 + Kk ] è equivalente a minimizzare g in Kk ; abbiamo inoltre già
dimostrato che g(y) ammette un minimo unico in Kk e questo minimo è per definizione xk+1 − x0
e quindi il secondo principio definisce un unico xk+1 tale per cui vale anche che:
⟨Axk+1 − b, y⟩ = 0 ∀y ∈ Kk
Ma Axk+1 − b per definizione è uguale al residuo rk+1 e quindi abbiamo:
⟨rk+1 , y⟩ = 0
che, essendo valida ∀y ∈ Kk , è essenzialmente la condizione che rk+1 ⊥ Kk . I due algoritmi sono
quindi identici e danno lo stesso valore unico xk+1 . Facciamo quindi qualche considerazione
sulla dimensione dei sottospazi di Krylov. Se il valore di k0 è uguale ad n − 1 allora abbiamo
dim(Kk0 ) = n e il sottospazio affine [x0 + Kk0 ] coincide esattamente con Rn ; in tal caso abbiamo
dimostrato che xk0 +1 = xn è il minimo di f in Rn che soddisfa Axn = b e quindi il metodo
converge con n iterazioni. Se invece k0 < n − 1 allora abbiamo visto che ∀k ≥ k0 dim(Kk ) =
k0 + 1; in particolare Ak0 +1 r0 ∈ Kk0 e quindi
k0
X
k0 +1
A r0 = Ai r0
i=0
e il coefficiente α0 deve essere non nullo; altrimenti infatti potremmo moltiplicare per A−1 (tanto
è definita positiva e quindi non singolare) e mostrare che Ak0 +1 è una combinazione lineare dei
vettori (r0 , Ar0 , ...., Ak0 −1 ) il che implicherebbe che Kk0 = Kk0 −1 , che sarebbe una contraddizione
della definizione di dimensione critica k0 . A questo punto possiamo concludere notando che
sostituendo all’equazione precedente otteniamo:
k0
1 k0 X αi i−1
A A r0 − A r0 + x0 = b
α0 α
i=1 0
51
Ovvero la soluzione di Ax = b appartiene allo spazio affine [x0 + Kk0 ] e quindi alla (k0 + 1)-esima
iterazione il minimo di f in Rn coincide con il minimo nello spazio affine [x0 + Kk0 ] e quindi xk0 +1
è la soluzione esatta. Abbiamo quindi convergenza in k0 + 1 < n iterazioni.
Per concludere questa sezione ricapitoliamo alcune cose. Abbiamo dimostrato che cercare
minimizzare il funzionale f (x) in una successione di spazi affini è esattamente identico
a risolvere le iterazioni del metodo del gradiente; inoltre abbiamo che la soluzione
esatta ed unica xesatta del sistema sarà raggiunta alla peggio con n iterazioni.
Una sequenza di vettori coniugati rispetto ad A è a tutti gli effetti ortogonale con rispetto al
prodotto scalare definito da ⟨Ax, y⟩ (questo è esattamente un prodotto scalare visto che A è sim-
metrica e definita positiva). Questa caratteristica da il nome al metodo.
Enunciamo ora un teorema di vitale importanza nella costruzione dell’algoritmo per implementare
il metodo del gradiente coniugato.
Teorema: Sia A una matrice simmetrica definita positiva e sia (xk ) la sequenza di soluzioni
approssimate del metodo del gradiente coniugato. Sia (rk = b − Axk ) la sequenza di residui
associata. Allora esiste una sequenza (pk ) A − coniugata tale che:
xk+1 = xk + αk pk
p0 = r0 = b − Ax0 , e per 0 ≤ k ≤ k0 , rk+1 = rk − αk Apk
pk+1 = rk+1 + βk pk
Con:
||rk ||2
αk =
⟨Apk , pk ⟩
||rk+1 ||2
βk =
||rk ||2
52
Il valore di questo teorema è che da un modo pratico di calcolare la k-esima iterazione senza dover
minimizzare il funzionale f (x) oppure ortogonalizzare rk rispetto ad Kk . Infatti otteniamo per
semplice sostituzione:
(xk+1 − xk )
pk = = rk + βk−1 pk−1
αk
αk βk−1
xk+1 = xk + αk (b − Axk ) + (xk − xk−1 )
αk−1
che è un induzione su xk simil-gradiente, ma con delle correzioni che accelerano la convergenza.
Abbiamo quindi un algoritmo semplice per calcolare ogni iterazione.
i grandi vantaggi di un algoritmo del genere sono: ad ogni iterazione è calcolato effettiva-
mente un solo prodotto matrice vettore (visto che rk è calcolato con le formule di induzione e
non con la definizione b − Axk ) ed per implementare l’algoritmo non è necessario salvare la
matrice A in memoria se sappiamo come calcolare il prodotto Ay per ogni y. E’ inoltre
possibile dimostrare che per il metodo del gradiente coniugato:
p k
p k2 (A) − 1
||ek ||2 = ||xk − x||2 ≤ 2 k2 (A) p ||x0 − x||2
k2 (A) + 1
da cui vediamo che: possiamo usare il gradiente coniugato come metodo iterativo, infatti
più iterazioni facciamo più diminuisce l’errore e la velocità di convergenza dipende dalla
radice del numero di condizionamento (mentre per il metodo del gradiente semplice dipendeva
dal numero di condizionamento) e quindi il gradiente coniugato converge più velocemente
ripetto al metodo del gradiente e più vicino è k2 (A) ad 1 più veloce sarà la convergenza
(è possibile provare che la velocità di convergenza del gradiente coniugato è ottimale, ovvero non
esiste un metodo più veloce per sistemi simmetrici definiti positivi). Tuttavia, siccome il fattore
di convergenza è moltiplicato per la radice del numero di condizionamento, per matrici mal-
condizionate diventa complesso l’uso di questo metodo per via della lenta convergenza
che ne consegue. Per ovviare al problema del numero di condizionamento della matrice possiamo
eseguire un precondizionamento; ovvero cerchiamo una matrice M simmetrica e definita positiva
che è in qualche senso una buona approssimazione di A−1 . Con tale matrice potremmo risolvere:
MAx = Mb
53
con MA ≃ I e k2 (MA) << k2 (A) che convergerà più velocemente e dal quale possiamo poi
riottenere per sostituzione le soluzioni. La matrice M è detta precondizionatore e il suo ruolo
è quello di comprimere gli autovalori di A in un intervallo il più piccolo possibile (idealmente
tutti allo stesso valore numerico). Ovviamente se sapessimo esattamente A−1 allora il problema
sarebbe risolto ma siccome in generale trovare l’inversa di una matrice non è cosı̀ semplice c’é
bisogno di trovare buone approssimazioni di ques’ultima; per esempio se la matrice A è a diagonale
dominante (ovvero l’elemento diagonale della riga i è maggiore in modulo della somma dei moduli
dei restanti elementi della riga) allora l’inversa della diagonale M = (diag(A))−1 non è una cattiva
approssimazione dell’inversa di A. Se poi diciamo di poter scrivere A = D + O (con D la matrice
degli elementi diagonali e O la matrice degli elementi non diagonali) e assumiamo che rispetto a
qualche norma valga ||D||s >> ||D||s allora possiamo pensare di espandere l’inversa con una serie
di taylor:
−1
A−1 = D(I + D−1 O) ≃ D−1 − D−1 OD−1 + ...
Se invece per esempio la nostra matrice è composta da pochi elementi grandi e molti piccoli possi-
amo approssimarla come una matrice sparsa (poniamo gli elementi piccoli uguali a 0) ed esistono
algoritmi per calcolare con uno scaling lineare dei costi computazionali le inverse di questo tipo
di matrici. Da tutti questi metodi sorge però un problema, ovvero che, anche se A è simmet-
rica definita positiva, non è detto che MA sia ancora simmetrica. Per risolvere bastano i
seguenti passaggi:
Inoltre nell’implementazione pratica dell’algoritmo non dobbiamo neanche fare il gradiente coniu-
gato per à ma basta introdurre un nuovo vettore, che indicheremo con z, definito come il risultato
dell’applicazione del precondizionatore al residuo:
zk = M −1 rk
z0 = M −1 r0
o classico, e il prodotto scalare che compare nelle formule del metodγk = ⟨rk , rk ⟩, viene sostituito
da
γk = ⟨rk , zk ⟩
L’algoritmo del gradiente coniugato precondizionato (PCG) differisce dunque dal CG standard solo
per pochi dettagli: ogni volta che compare il residuo rk , esso viene sostituito dal vettore precon-
dizionato zk , e i coefficienti αk e βk vengono calcolati usando tali quantità. In forma compatta,
l’algoritmo PCG può essere descritto come segue:
1. Inizializzazione:
x0 r0 = b − Ax0 z0 = M −1 r0 p0 = z0 γ0 = ⟨r0 , z0 ⟩
54
2. Iterazione per k = 0, 1, 2, . . . fino a convergenza:
γk
αk =
⟨pk , Apk ⟩
xk+1 = xk + αk pk
rk+1 = rk − αk Apk
zk+1 = M −1 rk+1
γk+1 = ⟨rk+1 , zk+1 ⟩
γk+1
βk =
γk
pk+1 = zk+1 + βk pk
A ∈ RN ×N con neig ≪ N
Non vogliamo dunque calcolare l’intero spettro, ma solo pochi autovalori: si tratta di un’esigenza
molto comune in chimica quantistica. Ad esempio, per ottenere l’energia dello stato fondamentale
basta il primo autovalore ed anche volendo calcolare alcuni stati eccitati, tipicamente ci si limita ai
primi 10 autovalori (e già questo può essere considerato un numero elevato). Situazioni di questo tipo
emergono frequentemente nei calcoli quantistici. Non è cosı̀, tuttavia, per i metodi self-consistent
field, in cui la dimensione della base (realizzata spesso con orbitali di tipo gaussiano) rimane relati-
vamente contenuta: anche calcoli molto grandi raramente superano qualche migliaio di funzioni di
base, con valori estremi intorno a 10000. In tali casi, una diagonalizzazione completa con algebra
lineare densa è più che adeguata, e non è necessario risparmiare sulle operazioni; per esempio invece
la situazione cambia radicalmente nella fisica dello stato solido con basi d’onda piana, dove è comune
lavorare con centinaia di milioni di funzioni di base. Analogamente, nella chimica quantistica, se
55
si vogliono calcolare stati eccitati tramite teoria del funzionale della densità dipendente dal tempo
(TD-DFT), è necessario determinare i primi ∼ 50 autovalori di matrici che possono raggiungere
dimensioni di 100 000 × 100 000. Assumiamo inoltre che la matrice sia simmetrica, condizione molto
comune per motivi legati alla meccanica quantistica. Le uniche matrici non simmetriche che
si incontrano tipicamente in questo ambito sono quelle provenienti dai metodi di cou-
pled cluster per stati eccitati. In questi casi occorre prestare attenzione, poiché gli autovalori
potrebbero non essere reali; autovalori complessi rappresenterebbero un problema fisico, e possono
effettivamente emergere in certi calcoli.
⟨x, Ax⟩
R(A, x) =
⟨x, x⟩
È possibile dimostrare che per una matrice simmetrica il minimo del quoziente di Rayleigh coincide
con il più piccolo autovalore della matrice. Questo risultato è noto come teorema di Reyleigh-Ritz
e per dimostrarlo notiamo che se la matrice A è simmetrica (nel caso complesso basta considerare
una matrice hermitiana e seguire lo stesso ragionamento) allora sarà scrivibile come:
A = QΛQ−1
Dove Λ è la matrice diagonale degli autovalori e Q è una matrice ortogonale (ancora nel caso
complesso semplicemente prenderemmo una trasformazione unitaria). Supponiamo sempre che gli
autovalori siano ordinati in ordine crescente e scriviamo il quoziente di Reyleigh come:
Pn 2
xT QΛQT x cT QT QΛQT QC cT Λc ci λ i
R(A, x) = = Pn 2 T
= Pn 2 = Pi=1 n 2
xT x i=1 ci qi qi c
i=1 i i=1 ci
dove abbiamo sfruttato l’ortogonalità di Q ed il fatto che dal teorema spettrale sappiamo che le sue
colonne formano una base di autovettori di A. Concentriamoci ora sul numeratore e sottraiamogli
il vettore λ1 c2 :
Xn n
X Xn
c2i λi − λ1 c2i = c2i (λi − λ1 )
i=1 i=1 i=1
56
e siccome abbiamo supposto l’ordine crescente degli autovalori abbiamo che la somma è sempre
maggiore o uguale a 0:
Xn
c2i (λi − λ1 ) ≥ 0
i=1
n
X n
X
c2i λi ≥ c2i λ1
i=1 i=1
Pn 2
ci λ i
Pi=1
n 2 ≥ λ1
i=1 ci
e con un procedimento analogo si dimostra che:
n
X
c2i (λi − λn ) ≤ λn
i=1
n
X n
X
c2i λi ≤ λn c2i
i=1 i=1
Pn 2
ci λi
Pi=1
n 2 ≤ λn
i=1 ci
Potremmo ora chiederci come calcolare un generico autovalore λk della nostra matrice A? La
risposta sta nel principio di Courant Fischer (o min-max principle). Vediamo infatti che,
calcolato il primo autovalore e il corrispettivo autovettore, possiamo definire il seguente sottospazio
vettoriale ortogonale:
W = {w : ⟨w, q1 ⟩ = 0}
e per definizione del sottospazio possiamo scrivere:
X
w= ci qi = Qc con c1 = 0
i=1
Otteniamo quindi una formulazione analoga del quoziente di Reyleigh nel sottospazio W :
c2 λ i
P
R(A, w) = Pi=2 i 2
i=2 ci
λ2 = minw̸=0,w∈W R(A, w)
λmax = maxw̸=0,w∈W R(A, w)
calcolati λ2 , q2 possiamo calcolare il terzo autovalore e cosı̀ via finché ci è necessario. In realtà
esiste un forma più generale del principio di Courant-fischer che permette di calcolare il generico
57
autovalore λk senza il calcolo degli autovettori o dei sottospazi ortogonali ma la sua implemen-
tazione è relativamente complicata e in ogni caso in chimica quantistica non presenta una grande
utilità. Nel paragrafo precedente, discutendo il metodo del gradiente coniugato, abbiamo visto che
la minimizzazione di un problema quadratico su un sottospazio limitato possiede un unico minimo
(a meno di degenerazioni). È quindi naturale cercare un insieme di minimizzatori ortogonali
del quoziente di Rayleigh per determinare i più piccoli autovalori della matrice. Questa
è precisamente la strategia adottata nel metodo di Davidson. Vediamone quindi la formulazione:
come accennato, il metodo di Davidson è un metodo basato sui sottospazi e l’idea è quella di es-
pandere gli autovettori cercati all’interno di un sottospazio V. Indichiamo con x gli autovettori che
vogliamo approssimare:
Xm
xi = uij vj
j=1
V = span{v1 , v2 , ..., vm }
⟨vi , vj ⟩ = δij
Neig ≤ m < N
Questo è intuitivo, infatti non possiamo ricavare cinque autovettori da un sottospazio bidimension-
ale; tuttavia allo stesso tempo, desideriamo che il sottospazio sia molto più piccolo dello spazio
originale, altrimenti non vi sarebbe alcun vantaggio computazionale. L’incognita del problema
diventa, quindi, la scelta del sottospazio di espansione e una strategia efficace per am-
pliarlo quando l’approssimazione corrente non risulta soddisfacente. Dobbiamo inoltre
determinare i coefficienti di espansione uij degli autovettori all’interno di V.
AX = XΛ
58
dove Λ è una matrice diagonale contenente gli autovalori cercati:
e, poiché siamo interessati solo ai primi neig autovalori, X è una matrice rettangolare. La nostra
ansatz consiste nello scrivere gli autovettori come combinazioni lineari delle colonne di V:
X = Vu
A(Vu) = (Vu)Λ
proiettiamo ora questa equazione sul sottospazio di espansione moltiplicando a sinistra per VT :
VT AVu = uΛ
notiamo che la matrice VT AV è una matrice (m × m) molto più piccola dell’originale. Diciamo
quindi che VVT agisce come un proiettore ortogonale da Rn sul sottospazio V e la matrice:
a = VT AV
au = uΛ
Le dimensioni in gioco sono piccole: se cerchiamo, ad esempio, 100 autovalori, in R2000 potremmo
usare un sottospazio di dimensione 2000 × 100, quindi una matrice ridotta di dimensione 100 × 100.
Con le librerie LAPACK moderne, una diagonalizzazione simmetrica di questa dimensione
è pressoché istantanea. A titolo informativo, le routine LAPACK per autovalori simmetrici
hanno nomi che iniziano con la lettera che rappresenta la precisione (S, D, C, Z), seguite da due
lettere per il tipo di matrice (ad esempio SY per matrici simmetriche, GE per matrici generali) e
infine una sigla che specifica l’operazione, come EV per autovalori e autovettori.
Una volta risolto il problema proiettato e ottenuti i coefficienti (u)ij , possiamo ricostruire gli au-
tovettori approssimati X = Vu e calcolare i residui:
R = AX − XΛ
Espansione del sottospazio: Se i residui non sono sufficientemente piccoli, nasce la ”domanda
da un milione di dollari”: come ampliare il sottospazio? Il principio è lo stesso dei metodi del
gradiente: i residui contengono informazioni sulle direzioni che mancano al sottospazio corrente.
Nei metodi del gradiente coniugato, i residui sono automaticamente ortogonali tra loro; qui, in
generale, ciò non accade, ma i residui restano comunque una fonte informativa preziosa.
59
Prima di vedere come ampliare effettivamente lo spazio introduciamo dei residui precondizionati.
L’idea, mutuata dai metodi iterativi per sistemi lineari, è quella di applicare un precondizionatore
che ”metta in evidenza” l’informazione più utile. Poiché il metodo fu sviluppato originariamente
per matrici diagonalmente dominanti (come nelle matrici CI o hamiltoniane elettroniche), un ottimo
precondizionatore è spesso la diagonale della matrice. Si introduce infatti:
D = diag(A)
Y = (D − λ1 I)−1 R
dove λ1 è il più piccolo autovalore approssimato. Sarebbe tentante precondizionare ogni residuo
con il proprio autovalore stimato, ma è una cattiva idea: usare precondizionatori diversi per
vettori diversi rende più difficile interpretare la struttura dello spazio espanso. È molto
meglio usare un singolo precondizionatore per tutti, tipicamente costruito usando il più piccolo
autovalore corrente. Questo costituisce la regola pratica più efficace. In ogni caso nei sistemi lineari
sappiamo che il miglior precondizionatore possibile sarebbe A−1 ma questo non si verifica nel caso
del problema agli autovalori. Infatti se il residuo è:
Y = (A − λ1 I)−1 R
se stiamo cercando un solo autovettore (Chiedi a Lipparini), questo residuo è parallelo all’autovettore
stesso e svolgendo il prodotto otterremo:
y = (A − λ1 I)−1 (A − λ1 I)x = x
• a differenza del gradiente coniugato, qui non possiamo garantire che ciò che aggiungiamo al
sottospazio sia automaticamente ortogonale alle direzioni già presenti.
Per questi motivi, serve un precondizionatore ”ragionevolmente buono”, ma non troppo accurato.
La diagonale della matrice comunque funziona quasi sempre molto bene. Se la matrice è quasi diag-
onalmente dominante tranne per un blocco problematico (per esempio, un blocco mal condizionato
o con autovalori molto ravvicinati), allora si può utilizzare:
• la diagonale di A per la parte ”buona”;
• un precondizionatore migliore (ad esempio un’approssimazione locale alla matrice) per il
blocco problematico.
60
Tuttavia, il precondizionamento per problemi agli autovalori è concettualmente molto più difficile
rispetto ai sistemi lineari, e non esiste una ricetta generale migliore della diagonale. Un’altra pos-
sibilità è usare una decomposizione incompleta (ad esempio incomplete LU), mantenendo solo gli
elementi sopra una certa soglia e invertendo l’approssimazione ottenuta. Tuttavia, questo approccio
rimane delicato e dipende molto dal problema.
Y1 = (I − VVT )Y
Y2 = orto(Y1 )
Vnew = V ∪ Y2
abbiamo cosı̀ un nuovo sottospazio più grande, su cui ripetere nuovamente la procedura:
calcolo del problema ridotto, calcolo dei residui, controllo di convergenza, espansione del sottospazio,
e cosı̀ via.
V ∈ Rn×mmax W ∈ Rn×mmax
61
• i vettori del sottospazio;
• i prodotti AV corrispondenti.
Questi contenitori vengono dimensionati per contenere il numero massimo di vettori che vogliamo
gestire. Tipicamente:
mmax = neig · m
inoltre quando il sottospazio cresce troppo, si effettua un restart: si scarta il sottospazio, lo si
rimpiazza con gli attuali autovettori approssimati, e si ricomincia l’espansione. Questo permette di
limitare il costo computazionale e la memoria utilizzata. Prima di procedere con la procedura vera
e propria però introduciamo una comoda notazione per le matrici:
G(a : b, c : d)
con questa notazione vogliamo dire di prendere le righe di G che vanno dalla riga a alla riga b e le
colonne della stessa matrice che vanno dalla colonna c alla colonna d. Spesso se vogliamo dire di
prendere o tutte le righe o tutte le colonne usiamo semplicemente ” : ”.
Possiamo quindi iniziare la procedura vera e propria assumendo come vettori iniziali di prova gli
autovettori associati ai più piccoli autovalori della diagonale di A:
diag(A)X = ΛXguess
con Xguess ∈ Rn×neig . Inizializziamo le prime neig colonne di V con questi vettori:
k =k+1
beg(k) = (k − 1)neig + 1
end(k) = kneig
la prima iterazione ricalcolerà quindi le prime neig colonne, la seconda partira dalla colonna neig + 1
e calcolerà altre neig colonne e cosı̀ via fino ad arrivare ad mmax . Questa scelta deriva da consid-
erazioni implementative: memorizziamo non solo i vettori che costituiscono la base ortonormale
dello spazio di espansione, ma anche l’applicazione della matrice a tali vettori. Ciò significa che
all’iterazione 2 non è necessario ricalcolare Av(1) , perché è già memorizzato. Una volta disponibile
la base aggiornata, calcoliamo la matrice proiettata. In realtà potremmo evitare questo step in
quanto la matrice ridotta a[k] è simmetrica e alla k-esima iterazione il blocco
62
ha dimensione knig × knig . Poiché alla seconda iterazione il blocco in alto a sinistra è già noto,
sarebbe possibile evitare ricalcoli. Tuttavia, dato che il costo di questo passaggio è O(m2 n), ossia
lineare nella dimensione del problema, lo spreco computazionale è trascurabile rispetto al costo
dell’applicazione di A, che è O(n2 m). A questo punto risolviamo il problema agli autovalori ridotto
a[k] u[k] = u[k] Λ[k]
e costruiamo una nuova approssimazione degli autovettori completa:
X[k] = V[k] u[k]
dove u[k] contiene gli autovettori del problema proiettato. Possiamo, inoltre, calcolare:
R[k] = AX[k] − X[k] Λ[k]
Sostituendo X[k] = V[k] u[k] :
R[k] = AV[k] u[k] − X[k] Λ[k] = W[k] u[k] − X[k] Λ[k]
Poiché
AX[k] = A(V[k] u[k] ) = (AV[k] )u[k]
preservando in memoria i vettori AV[k] possiamo evitare un’ulteriore moltiplicazione matrice-
vettore. Successivamente completiamo il precondizionamento dei residui:
Y[k] = (D − λ1 I)−1 R[k]
Proiettiamo Y[k] sul sottospazio:
[k]
Y1 = (I − V[k] V[k]T )Y[k]
e ortogonalizziamo:
[k] [k]
Y2 = orto(Y1 )
A questo punto possiamo aggiornare il sottospazio:
[k]
V[k+1] (:, beg(k + 1) : end(k + 1)) = Y2
Ovviamente è necessario verificare la convergenza dei residui e gestire la procedura di restart. Se la
dimensione della base diventa troppo grande per essere memorizzata, si reimposta k = 1 e si riparte
utilizzando come nuovo guess l’approssimazione corrente degli autovettori.
Un punto essenziale da trattare riguarda l’ortonormalizzazione della base. Produrre vettori tra
loro ortonormali tramite sottrazione successiva delle proiezioni è concettualmente semplice, ma
nella pratica numerica è altamente instabile, poiché porta ad accumulo degli errori. La risposta
corretta, in generale, alla domanda: ”come ortonormalizzare un insieme di vettori?” è quasi sempre
tramite decomposizione ai valori singolari (SVD). Per una matrice rettangolare A, la sua SVD:
A = UΣV⊤
fornisce nelle prime colonne di U una base ortonormale dello spazio immagine di A. Quindi, per
ortonormalizzare le colonne di A è sufficiente sostituirle con le colonne di U:
[k] [k]
Y1 = UΣVT −→ Y2 = U
Si tratta del metodo numericamente più stabile. Tuttavia, è anche costoso: per matrici sim-
metriche si può utilizzare la routine LAPACK DSYSVD, mentre per matrici rettangolari generiche
occorre ricorrere a DGESVD.
63
1.14.2 Problemi Numerici del metodo di Davidson
Sebbene l’algoritmo appaia completo e corretto sulla carta, la sua implementazione diretta porta
a risultati deludenti. Anche per matrici ben condizionate, il metodo può non convergere ai primi
10 autovalori, oppure convergere verso valori che non appartengono allo spettro della matrice origi-
nale: questi sono i cosiddetti intruder states (in letteratura matematica il fenomeno è noto come
spectral pollution). Possiamo individuare facilmente quale sia lo step numericamente delicato
nel metodo di Davidson. Il punto chiave è che, indipendentemente dal condizionamento
della matrice originale, la matrice dei residui costruita durante l’algoritmo risulterà
inevitabilmente molto mal condizionata. È importante ricordare che una matrice numeri-
camente mal condizionata è una matrice le cui colonne diventano, dal punto di vista numerico,
quasi linearmente dipendenti, a causa della limitata precisione disponibile. Un esempio intuitivo:
se costruiamo una matrice le cui colonne hanno norma dell’ordine di 10−9 , quanto possiamo as-
pettarci che siano linearmente indipendenti in doppia precisione? Al massimo per un numero di
cifre paragonabile alle nove cifre significative delle norme; in singola precisione, colonne cosı̀ piccole
risulterebbero invece praticamente identiche. Lo step realmente problematico in Davidson è quindi
l’ortogonalizzazione, ovvero la trasformazione dei vettori residui rispetto allo spazio già costruito.
L’operazione si presenta schematicamente come:
[k] [k]
Y1 (VΣ−1 ) = Y1 T = U
l’ortogonalizzazione implica quindi una trasformazione del tipo:
[k] [k]
Y2 = Y1 T = U
Finché V è una matrice ortogonale, non introduce problemi: il numero di condizionamento di una
matrice ortogonale è pari a 1 nel norma-2. Il problema nasce da Y1 : i vettori residui tendono per
costruzione a diventare sempre più piccoli man mano che l’iterazione prosegue, e successivamente
vengono ulteriormente ridotti dal processo di proiezione sul complemento ortogonale dello spazio
già costruito. Questo li rende ancora più piccoli e più numericamente instabili. Se consideriamo i
valori singolari di Y1 , poiché le sue colonne non sono nulle, esiste almeno un vettore linearmente
indipendente, e dunque almeno un valore singolare sarà dell’ordine dell’unità. Tuttavia, tutte le
altre direzioni saranno praticamente uguali al vettore principale fino all’ordine di grandezza delle
loro norme che sono estremamente piccole per costruzione); di conseguenza gli altri valori singolari
saranno anch’essi molto piccoli, ad esempio dell’ordine di 10−9 e il numero di condizionameno di
Y1 sarà quindi:
σmax
κ(Y1 ) = ≈ 109
σmin
un valore numericamente pessimo. Questa instabilità non dipende dalla matrice originaria:
si verifica per qualunque matrice perché i residui, per definizione, devono convergere
a zero.
Consideriamo ora uno specifico vettore yi in Y1 . Sia ŷi la sua versione computata. Dopo la
proiezione, l’errore di ortogonalità rispetto allo spazio già costruito è dell’ordine della precisione di
macchina:
ϵorth = O(ϵmach )
Ma quando ortogonalizziamo ulteriormente per ottenere Y2 , l’errore viene amplificato dal numero
di condizionamento della trasformazione:
||δy2 || ≈ κ(Y1 )ϵmach
64
Se κ(Y1 ) ≈ 109 ed ϵmach ≈ 10−16 (doppia precisione), allora l’errore risultante è dell’ordine di
10−7 , che è enormemente più grande della precisione macchina e rappresenta una perdita completa
dell’ortogonalità. Un errore di ortogonalità di 10−7 è già molto grave: anche 10−10 o 10−12 sarebbe
problematico, poiché l’intero metodo si basa sull’uso di basi ortonormali. Errori cosı̀ grandi
portano l’implementazione a comportarsi di fatto come un sofisticato generatore di
numeri quasi casuali, poiché i vettori perdono ogni legame con il problema originale.
Se continuiamo ad aggiungere vettori che non sono esattamente ortogonali, quando formiamo il
problema proiettato,
V⊤ AVu = VT VuΛ
utilizziamo implicitamente l’identità V⊤ V = I. Ma se i vettori non sono esattamente ortogonali,
questa identità è falsata. In linea teorica dovremmo risolvere un problema agli autovalori general-
izzato:
(V⊤ AV)u = λ(V⊤ V)u
tenendo conto dell’errore di ortogonalità. Poiché però non lo facciamo, finiamo per risolvere
un problema che non corrisponde più all’originale, compromettendo completamente la
correttezza del metodo. Per ovviare a questo problema consideriamo ora la matrice Y2 . Le sue
colonne, fatte salve piccole imperfezioni dovute all’ortogonalizzazione, sono ortonormali tra loro;
ciò significa che ogni colonna ha norma prossima a 1. Questo è un oggetto molto più stabile di
quanto non fossero i vettori iniziali. Rimane però un problema: Y2 non è esattamente ortogonale
rispetto allo spazio di base V. La soluzione consiste semplicemente nel ripetere il processo:
Y3 = (I − VV⊤ )Y2
e successivamente:
Y4 = orto(Y3 )
Perché questo funziona? Perché mentre in precedenza la matrice da ortogonalizzare conteneva
colonne con norma piccolissima (e dunque pessime da un punto di vista numerico), ora contiene
vettori con norma ≈ 1. L’operazione
(I − VV⊤ )
elimina l’errore di ortogonalità residuo, dell’ordine di 10−4 nel caso peggiore. Sottrarre da un vettore
di norma 1 una componente di ampiezza 10−4 lascia comunque un vettore di norma ≈ 1 − 10−4 ≈ 1.
Dunque, al passo successivo l’ortogonalizzazione è stabile: i vettori sono già quasi ortonormali.
Tuttavia non tutti i problemi del metodo di Davidson scompaiono: la matrice Y2 iniziale era lontana
dall’essere ortogonale, ma almeno aveva colonne di norma dell’ordine dell’unità, non dell’ordine di
10−8 . Questo la rende un ottimo punto di partenza. Dopo la seconda ortogonalizzazione, i valori
singolari risultano tutti molto vicini a 1 e quindi il numero di condizione della trasformazione è
vicino a 1, come confermato dalla figura.
Facciamo inoltre qualche considerazione sulla proprietà associativa delle matrici; nel corso
dell’algoritmo la utilizziamo spesso implicitamente:
A(VU) = (AV)U
Questa equivalenza è valida in aritmetica esatta, ma non lo è in aritmetica finita. Perché qui
funziona? Perché sia V che U sono matrici ortogonali: moltiplicare o pre-moltiplicare per una
matrice ortogonale non amplifica gli errori, in quanto il numero di condizione è pari a 1. Al
65
contrario, se tentassimo di applicare la stessa logica ai nuovi vettori non perfettamente ortogonali,
ad esempio:
AY2 = A(I − VV⊤ )YT
ci troveremmo a moltiplicare A per una matrice mal condizionata, peggiorando drasticamente la
propagazione degli errori. L’ortogonalizzazione ripetuta rimette invece tutto sotto controllo: si ap-
plica A solo a matrici ortogonali, e poi si moltiplica per matrici ortogonali. Nessuna amplificazione
dell’errore.
Dobbiamo fare ora un’altra considerazione sull’algoritmo usato per l’ortogonalizzazione dei residui.
Sebbene la SVD sia sempre stabile e non esista una perturbazione numerica che non le permetta di
produrre una base realmente ortonormale, essa ha un costo elevato (O(n2 m)) e quindi all’aumentare
della dimensione del problema diventa sempre meno efficiente. Potremmo quindi chiederci perché
non usare una fattorizzazione di Cholesky usando la seguente matrice di Overlap:
O = Y1⊤ Y1
che è idealmente una matrice simmetrica definita positiva, per cui quindi vale prendere la fattoriz-
zazione di Cholesky:
O = LL⊤
da cui:
L−1 Y1⊤ Y1 L−⊤ = I
e definire una trasformazione ortogonalizzante tramite
Y2 = Y1 L−⊤
che produce vettori ortonormali. Questa strategia è pessima per due motivi:
1. Instabilità numerica: se i vettori di Y1 hanno norma 10−8 o 10−9 , allora gli elementi
diagonali di O sono dell’ordine di 10−16 o 10−18 , cioè inferiori o prossimi alla precisione
macchina. Non è garantito che O venga riconosciuta come definita positiva da LAPACK, e
la routine Cholesky potrebbe fallire.
ma poiché L deriva dall’inversione di una matrice quasi singolare, la sua norma sarà enorme.
La norma di L−⊤ prende il nome di growth factor, perché quantifica quanto cresce l’errore di
ortogonalità durante l’ortogonalizzazione. In sintesi esistono due problemi fondamentali:
66
1. Cholesky può fallire. Se l’overlap ha elementi diagonali estremamente piccoli, la matrice
potrebbe non essere riconosciuta come definita positiva in precisione finita, causando il falli-
mento della fattorizzazione.
2. Anche quando funziona, il risultato è pessimo. Il fattore L sarà mal condizionato;
risolvere il sistema triangolare
Y2 = Y1 L−⊤
introdurrà errori rilevanti e Y2⊤ Y2 non sarà affatto vicino all’identità.
Tuttavia vorremmo poter usare Cholesky vista la sua economicità, ma ovviamente non possiamo
pagarne gli svantaggi. L’idea è imbrogliare leggermente:
1. Si consideri l’overlap:
O = Y1⊤ Y1
Oϵ = O + αI
dove α è scelto come un multiplo (ad esempio 100 volte) della precisione macchina moltiplicata
per la norma di O.
3. Se la fattorizzazione di Cholesky fallisce ancora, si aumenta α di un fattore 10 e si riprova.
Nella pratica, un singolo tentativo è quasi sempre sufficiente.
Cosı̀ otteniamo dei fattori di Cholesky L per la matrice leggermente modificata Oϵ , e possiamo
risolvere il sistema triangolare:
Ỹ2 = Y1 L−⊤
I problemi di questa strategia sono:
1. Stiamo ancora invertendo una matrice molto mal condizionata: l’ortogonalizzazione risultante
sarà imprecisa.
2. Inoltre non stiamo neppure ortogonalizzando con la matrice giusta, ma con una versione
”perturbata” dell’overlap.
Tuttavia dobbiamo chiederci, quanto è grave questa imprecisione? Supponiamo di aggiungere
un termine α dell’ordine 10−12 e che il numero di condizionamento diventi 1010 . Ci aspettiamo
allora che l’ortogonalità ottenuta sia accurata a circa 4-5 cifre decimali: ovvero errori dell’ordine di
10−4 − 10−5 . Non è bello, ma non è neanche disastroso. In particolare, dopo questo passo abbiamo
ottenuto vettori Ỹ2 con 4-5 cifre di accuratezza. Un risultato insufficiente per il metodo di David-
son (che amplificherebbe questi errori fino a distruggere l’informazione), ma perfetto come punto
di partenza per un’ortogonalizzazione ulteriore. Perché l’ortogonalizzazione via Cholesky produce
vettori scarsi? Perché il punto di partenza era pessimo: vettori minuscoli, mal condizionati, con-
taminati dall’errore numerico. Cholesky non può fare miracoli.
Una volta ottenuti vettori con norma circa 1 e con errori dell’ordine 10−4 − 10−5 , la situazione
cambia: ora il problema consiste nel correggere un errore piccolo rispetto al vettore stesso, non nel
trasformare vettori di norma 10−9 in vettori di norma 1. A questo punto:
67
Y3 = orto(Ỹ2 )
opera su vettori già ”quasi buoni”, e può correggere l’errore residuo. Dopo il secondo passaggio di
ortogonalizzazione, i vettori risultano ortonormali fino alla precisione macchina:
Y3⊤ Y3 = I + O(ϵmach )
La doppia ortogonalizzazione produce vettori ortonormali con accuratezza comparabile a
quella del SVD, ma a costo molto inferiore. Questa strategia non è immediata né ovvia: richiede
una chiara comprensione delle cause dell’errore, della sua propagazione e degli strumenti numerici
disponibili. Il messaggio finale è che: se conosci bene il tuo problema, se sai stimare gli errori
e prevedere come si propagano, allora puoi progettare soluzioni numeriche efficienti
e robuste. L’idea generale comunque è effettuare l’ortogonalizzazione e l’ortonormalizzazione
in maniera iterativa, fino a quando il numero di condizione della trasformazione ortonormalizzante
risulta minore di 2. Questo garantisce che un errore di ortogonalità su un singolo vettore dell’ordine
della precisione macchina si traduca in un errore massimo di soli due volte tale precisione, ancora
perfettamente accettabile.
Per un generico sistema molecolare con N elettroni e M nuclei possiamo scrivere l’operatore hamil-
toniano in unità atomiche come:
N M N XM N X N M X
M
X 1 X 1 X Za X 1 X ZA zB
H=− ∇2i − ∇2A − + + ([Link])
i=1
2 2MA i=1
riA i=1 j>i rij RAB
A=1 A=1 A=1 B>A
N N XM N X N M X
M
X 1 X Za X 1 X ZA zB
Ĥe = − ∇2i − + + ([Link])
i=1
2 i=1
riA i=1 j>i rij RAB
A=1 A=1 B>A
M
X 1
Ĥn = − ∇2 ([Link])
2MA A
A=1
Dove Ĥe è l’hamiltoniano elettronico e Ĥn è l’operatore di momento nucleare. Praticamente sempre
in chimica si applica la cosidetta approssimazione di Born-Oppenheimer. Questa approssi-
mazione si basa sul fatto che, poiché i nuclei sono molto più massivi degli elettroni, allora possiamo
68
considerare disaccoppiati i moti nucleari con quelli elettronici. Supponiamo di poter scrivere la
funzione totale come:
XK
Ψ(r, R) = χk (r; R)ϕk (R) ([Link])
k
dove K è il basis set e χk sono autofunzioni dell’hamiltoniano elettronico dipendenti solo ”para-
metricamente” dalla geometria molecolare R (ovvero sebbene siano funzioni solo di r, la loro forma
funzionale dipende strettamente da R). Le funzioni ϕk (R) sono invece delle funzioni che tengono
conto solo della geometria molecolare.
Abbiamo sempre che vale:
⟨χk′ (r; R)||χk (r; R)⟩(r) = δk′ k ([Link])
dove il pedice (r) indica che l’integrazione avviene solo sulle coordinate elettroniche. Possiamo
quindi definire la matrice He come la matrice (K × K) con i seguenti elementi:
(He )kk′ = ⟨χk′ (r; R)|Ĥe |χk (r; R)⟩(r) = δk′ k Ek (R) ([Link])
Vediamo quindi che Eq. [Link] forma un insieme di equazioni agli autovalori accoppiate.
L’accoppiamento deriva dai termini non diagonali di Hn (R), che sono detti termini di accoppia-
mento vibronico poiché dipendono dall’accoppiamento tra l’energia cinetica nucleare ed elettronica.
Per completare la trattazione sulla Born-Oppenheimer iniziamo studiando gli elementi:
⟨χk |PAα |χk ⟩(r) ([Link])
⟨χk′ |PAα |χk ⟩(r) ([Link])
ovvero gli elementi rispettivamente diagonali e non diagonali. L’elemento diagonale è sempre nullo
per funzioni d’onda reali; questo deriva dal fatto che esso è un valore di aspettazione, e quindi per
definizione reale, tuttavia per come abbiamo definito PAα esso è immaginario e siccome l’integrale
è anch’esso su variabili reali avremmo:
∂
Z = ⟨χk |PAα |χk ⟩(r) = −i⟨χk | |χk ⟩(r) = −iK ([Link])
∂RAα
69
ma l’unico numero reale e complesso allo stesso tempo è 0 (per una dimostrazione più formale si
può usare l’operatore di inversione temporale e le sue proprietà). Per quanto riguarda gli elementi
non diagonali Eq. [Link] essi sono il centro dell’accoppiamento vibronico e nel caso in cui
questi si annullino allora i moti elettronici e nucleari saranno completamente separabili.
Per capire quando questo è il caso consideriamo il seguente elemento di matrice:
⟨χk′ |[P̂Aα , Ĥe ]|χk ⟩ = ⟨χk′ |P̂Aα Ĥe |χk ⟩ − ⟨χk′ |Ĥe |ˆ⟩PAα χk ([Link])
Sfruttando che l’hamiltoniano è un operatore hermitiano possiamo scrivere:
⟨χk′ |[P̂Aα , Ĥe ]|χk ⟩ = Ek ⟨χk′ |P̂Aα |χk ⟩ − Ek′ ⟨χk′ |P̂Aα |χk ⟩ ([Link])
che riordinando i termini diventa:
⟨χk′ |[P̂Aα , Ĥe ]|χk ⟩
⟨χk′ |P̂Aα |χk ⟩ = ([Link])
Ek − Ek′
Questa relazione lega quindi Eq. [Link] con l’energia di due diversi stati stazionari elettronici
del sistema. All’aumentare della differenza ∆E = Ek − Ek′ il contributo dei termini non diagonali
diventa sempre minore fino ad essere trascurabile. Per quanto riguarda i termini ⟨χk′ |Tn |χk ⟩(r)
possiamo considerare nulli gli elementi non diagonali seguendo lo stesso ragionamento fatto in Eq.
[Link]; per quanto riguarda invece i componenti diagonali essi sono spesso trascurati e possono
essere considerati come correzioni all’energia. Entro queste approssimazione possiamo riscrivere
Eq. [Link] come:
(Hn (R))kk′ = δkk′ Tn ([Link])
che è una matrice diagonale, permettendoci quindi di separare perfettamente le funzioni nucleari
ed elettroniche:
[Tn + Ek (R)]ϕk (R) = Eϕk (R) ([Link])
ovvero possiamo calcolare l’energia elettronica a varie geometrie nucleari e poi calcolare
l’energia totale per ogni geometria. Stati elettronici per cui vale questa approssimazione sono
∼
detti adiabatici. Va da se che nel caso di stati degeneri o quasi degeneri abbiamo Ek (R) = Ek′ (R)
e quindi i termini di Eq. [Link] non sono trascurabili e vanno risolte le equazioni accoppiate. Le
intersezioni tra superfici di energia potenziale (PES) sono dette in generale intersezioni coniche.
+ altri termini
([Link])
70
il primo termine è detto E0 ed è nullo poiché la derivata calcolata nel minimo di potenziale è per
definizione nulla:
3n 2 (el)
(el) 1 X ∂ E
E (R1 , ..., Rn ) = E0 + (Ri − Ri0 )(Rj − Rj0 ) ([Link])
2 i=1,j=1 ∂Ri ∂Rj 0
definiamo quindi gli spostamenti dj = Rj − Rj0 ed otteniamo quindi in forma matriciale Eq. [Link]:
1 †
E (el) ≃ E0 + d Hd ([Link])
2
con H la matrice simmetrica Hessiana di elementi:
2 (el)
∂ E
Hij = ([Link])
∂Ri ∂Rj
dove Mii è un elemento della matrice (3n × 3n) diagonale delle masse nucleari M che ha come
elementi diagonali:
M3n = M3n−1 = M3n−2 = Mn ([Link])
ovvero associamo la massa del nucleo n-esimo ad ognuna delle sue coordinate. Secondo questo
procedimento i nuclei sono considerati distinguibili e quindi non consideriamo effetti quantistici
derivanti dall’indistinguibilità di nuclei identici.
Per semplificare la trattazione ed alleggerire la notazione definiamo delle nuove coordinate pesate
sulle masse:
q = M1/2 d ([Link])
dove abbiamo che M1/2
p
= δij Mij . L’hamiltoniano può essere riscritto come:
ij
3n 3n
1 X ∂2 1 † −1/2 −1/2 1 X ∂2 1 †
ĤBO − E0 ≃ − 2 + q M HM q = − 2 + q Hq ([Link])
2 i=1 ∂qi 2 2 i=1 ∂qi 2
dove la radice si assorbe nella definizione delle nuove coordinate ed H è una matrice che ha come
elementi:
Hij
(H)ij = p ([Link])
Mii Mjj
Per passare ad una somma di operatori indipendenti facciamo una trasformazione unitaria lineare
delle coordinate q. Sappiamo infatti dal teorema spettrale che H sarà diagonalizzabile in quanto
matrice reale e simmetrica:
†
HU = UΛ → H = UΛU ([Link])
71
†
con U = U−1 e Λ una matrice diagonale di autovettori di H. La matrice U ha come colonne
gli autovettori di H. Sostituendo la forma decomposta spettralmente dell’Hessiana in Eq. [Link]
otteniamo:
3n
1 X ∂2 1 † †
ĤBO − E0 ≃ − 2 + q UΛU q ([Link])
2 i=1 ∂qi 2
†
ovvero abbiamo definito un nuovo set di coordinate Q = U q che diagonalizza l’Hessiana e che
sono ottenute per trasformazione unitaria delle vecchie. Queste coordinate sono dette coordinate
normali e, nell’approssimazione armonica, l’operatore di Born-Oppenheimer può essere scritto
come somma di operatori indipendenti:
3n 3n 3n
1 X ∂2 1X X
ĤBO − E0 ≃ − 2 + Λii Q2i = ĥ(Qi ) ([Link])
2 i ∂Qi 2 i i
dove il primo termine rimane formalmente uguale poiché il laplaciano è invariante per trasformazioni
unitarie. La soluzione sarà quindi il prodotto delle funzioni soluzione per i singoli operatori ĥ(Qi ):
Notiamo che dei 3n autovalori 6 (o 5 per molecole lineari) sono nulli. Infatti abbiamo considerato
tutte le coordinate nucleari ma traslazioni e rotazioni infinitesime (diciamo infinitesime poiché in
realtà gli stati rotazionali e quelli vibrazionali sono accoppiati per via della presenza del momento
di inerzia nel primo dei due) non fanno variare l’energia. Per rotazioni e traslazioni sopravvive
quindi solo il contributo cinetico e possiamo esprimere i rispettivi operatori come:
1
T̂cm = − ∇2 ([Link])
2Mtot cm
L̂1 L̂2 L̂3
T̂rot = + + ([Link])
2I1 2I2 2I3
X n
Mtot = mi ([Link])
i
72
Se vogliamo vedere l’effetto di una sola coordinata normale basta porre tutte le coordinate uguali
a 0 e vedere l’effetto sugli spostamenti. Questo procedimento può essere eseguito con Gaussian,il
quale pone la coordinata di interesse come Qk = A cos (ωk t + ϕ) e permette di vedere in modo
oscillante lo spostamento dei nuclei dovuto a quella coordinata. Queste soluzioni che abbiamo
ottenuto sono esattamente analoghe a quelle che si otterrebbero in meccanica lagrangiana per un
sistema di oscillatori armonici accoppiati classici (infatti non abbiamo
√ considerato effetti quantistici
per i nuclei). Secondo questo ragionamento vale che ωk = Λii (che è quello che si ottiene per
un oscillatore armonico classico). Tutte queste radici sono reali e positive. Questo deriva
dal fatto che tendenzialmente le superfici di potenziale hanno concavità verso l’alto, nei casi in
cui questo non vale Gaussian da frequenze negative (ma in realtà sono immaginarie). Casi del
genere possono comunque essere di interesse, per esempio spesso punti di sella rappresentano stati
di transizioni e abbiamo per queste zone della superficie un solo autovalore immaginario (se cerco
zone in cui ho un solo autovalore immaginario posso quindi trovare facilmente delle selle).
Passiamo ora a trattare le soluzioni dei singoli operatori:
1 ∂2 1
ĥ(Qi ) = − + Λii Q2i ([Link])
2 ∂Q2i 2
come abbiamo detto questo è proprio l’operatore dell’oscillatore armonico e le soluzioni possono
essere espresse come:
1p
χki (Qi ) = Hki exp (− Λii Q2i ) ([Link])
2
con Hki polinomi di Hermite di ordine ki . Gli autovalori saranno invece:
p 1
ϵki =
Λii ( + ki ) ([Link])
2
√
e sono tutti equispaziato con differenza ∆E = Λ = ω. Per ogni stato è verificato:
che vale per via del fatto che i polinomi di Hermite di ordine pari sono dispari e quelli di ordine
dispari sono pari.
La funzione nucleare e i suoi autovalori di energia possono quindi essere espressi come:
3n−6(5)
Y
⃗ cm , α, β, γ, Q1 , ..., Q3n−6(5) ) = ϕk (R
Φ(n) (R ⃗ cm )ϕrot
j,m (α, β, γ) χνl (Ql ) ([Link])
l=1
3n−6(5)
k2 X p 1
E = E0 + + ϵrot + Λl ( + νl ) ([Link])
2Mtot 2
l=1
Con queste equazioni possiamo simulare uno spettro vibrazionale (IR o Raman) di una
generica molecola; ovvero possiamo studiare l’interazione di una radiazione elettromagnetica con
la nostra molecola, per poi analizzare la radiazione trasmessa o diffusa. Per via dell’interazione
radiazione-materia nell’hamiltoniano della molecola isolata si aggiunge un termine, detto proprio
di interazione:
Ĥtot = Ĥmol + Ĥint ([Link])
73
e dovremmo usare l’equazione di schroedinger dipendente dal tempo. Noi non comunque non
siamo tanto interessati a risolvere l’equazione di per se ma a vedere se dall’interazione consegue un
assorbimento di fotone e, siccome l’interazione radiazione-materia può in prima approssimazione
essere ”mediata” dal momento di dipolo molecolare, scriviamo il termine di interazione come:
⃗ˆ · E
Ĥint = −µ ⃗ rad ([Link])
ovviamente la radiazione è assunta polarizzata e la direzione è scelta dall’operatore. Il momento di
dipolo invece è solidale ad un s.d.r molecolare (in mezzi dispersi quindi l’orientazione campo-dipolo
è variabile, mentre in mezzi solidi è fissata). L’operatore di momento di dipolo può essere scritto
per n nuclei e N elettroni in unità atomiche come:
n
X N
X
µ̂ = ⃗a −
za R ri ([Link])
a=1 i
che è moltiplicativo poiché siamo nella rappresentazione delle coordinate. La radiazione può essere
espressa come:
⃗ ⃗
Erad = E0 exp i(ωt − k · ⃗r) ≃ E0 exp iωt ([Link])
dove l’approssimazione vale poiché tendenzialmente la lunghezza d’onda è molto maggiore delle
dimensioni molecolari e quindi possiamo trascurare il contributo di ⃗k · ⃗r. Solitamente il grafico
di uno spettro ha sull’asse delle x la frequenza della radiazione e sulle y la probabilità di
transizione Wi→f , la quale è proporzionale al seguente elemento di matrice:
fermandoci solo ai termini lineari. Possiamo quindi ottenere il valore di aspettazione del momento
di dipolo (ovvero il momento che misuriamo sperimentalmente) come:
⟨µ̂⟩mol = ⟨Φ(n) |µ̂|Φ(m) ⟩(R) ([Link])
74
Se ci limitiamo ai moti nucleari l’hamiltoniano di interazione sarà quindi esprimibile come:
ĥint = −⃗ ⃗
µ(Q) · E ([Link])
e la probabilità di transizione sarà:
Wi→f = |⟨Φf |ĥint |Φi ⟩(R) |2 ([Link])
Dobbiamo poi considerare i contributi rotazionali, siccome il campo elettrico è fissato dal laboratorio
mentre la molecola si trova in un mezzo disperso o in soluzione e quindi varia la sua orientazione
rispetto al sistema di riferimento del laboratorio:
µ ⃗ = |⃗
⃗ ·E µ|E0 cos (Ω) ([Link])
con Ω angolo tra il momento di dipolo e la radiazione. Possiamo quindi avere uno stato finale
che differisce da quello iniziale solo per la parte rotazionale (tipicamente nelle microonde) e quindi
avremo:
(rot) (rot)
⟨ϕf ||⃗ µ0 |E0 cos (Ω)|ϕi ⟩(R) ([Link])
dove possiamo non considerare tutti i contributi non rotazionali poiché si integrano ad 1. Abbiamo
usato in Eq. [Link] direttamente µ ⃗ 0 e non lo sviluppo poiché nelle microonde è il contributo
più importante (assumiamo lo spostamento dalla geometria di equilibrio piccolo nella spettroscopia
rotazionale e quindi se µ
⃗ 0 è nullo la molecola non assorbe nelle microonde in prima approssimazione).
L’angolo Ω è legato agli angoli di eulero con cui definiamo la funzione rotazionale. Noi comunque
siamo interessati solo solo alla parte vibrazionale (anche se ricordiamo che non sono del tutto
indipendenti). Nella nostra approssimazione abbiamo che la probabilità di transizione sarà:
3n−6(5) 3n−6(5) 3n−6(5)
Y
⃗
X ∂ µ̂ Y
Wi→f = |⟨ϕel ϕrot χνj (Qj )| − E · µ̂0 + Qi |ϕel ϕrot χνk (Qk )⟩|2
j=1 i
∂Q i 0 k=1
([Link])
ma µ
⃗ 0 è una costante e quindi possiamo portarlo fuori; le funzioni rotazionali e elettroniche sono
anch’esse costanti e si integrano ad uno. L’equazione di riduce quindi ad una somma sulle coordinate
normali, ovvero ad una somma dei seguenti elementi di matrice:
∂ µ̂ ⃗ 0 |χν (Ql )⟩|2
|⟨χν±1 (Ql )| Ql · E ([Link])
∂Ql 0
questo è valido per via delle proprietà dei polinomi di Hermite, infatti affinché l’integrale non sia
nullo lo stato di arrivo deve essere uguale a ν ± 1 quello di partenza. Ovviamente affinché sia
possibile la transizione deve anche essere sempre che:
∂ µ̂
̸= 0 ([Link])
∂Ql 0
Infatti è questa quantità che decide l’intensità del picco di emissione. Se il momento di
dipolo non cambia con la coordinata normale Ql allora la transizione non può avvenire. Tutti i
termini che abbiamo evitato di scrivere sono costanti che semplicemente scalano il risultato.
Un esempio tipico di spettro IR è la molecola di CO2 , la quale possiede 3n − 5 = 9 − 5 = 4 modi
normali di vibrazione. Di questi 2 sono moti di bending degeneri che fanno variare il momento di
dipolo; gli altri due sono moti di stretching di cui uno simmetrico (il dipolo non cambia) e uno
asimmetrico per il quale varia il momento di dipolo. In uno spettro IR quindi l’unico moto che non
vedremo sarà lo stretching asimmetrico.
75
2.3 Spettri Raman e polarizzabilità
In uno spettro Raman analizziamo la radiazione diffusa. La sorgente è un laser con una lunghezza
d’onda nel visibile. In uno spettro Raman abbiamo un picco centrale dovuto alla diffusione
Rayleigh (che a sua volta è una risposta molecolare alla luce incidente legata alla polarizzabilità
della molecola); ed attorno a questo picco abbiamo altre righe disposte in modo simmet-
rico rispetto al picco centrale. A destra abbiamo le bande anti-Stokes e a sinistra le Stokes.
Queste righe dipendono da un processo a due fotoni con uno stato intermedio che è uno stato
virtuale, ovvero uno stato non stazionario e con tempo di vita nullo. L’anti-Stokes dipende da un
fotone entrante ad energia maggiore di quella uscente, viceversa per le bande Stokes. A temperatura
ambiente quindi le righe Stokes saranno molto più intense. Le bande Stokes ed Anti-Stokes sono
simmetriche rispetto allo scattering di Rayleigh, ovvero:
ωAnti−Stokes − ω0 = ω0 − ωStokes ([Link])
Per spiegare questo fenomeno si introduce il concetto di diffusione come risposta della materia a
un campo esterno dinamico, in questo caso il campo elettromagnetico della radiazione incidente.
La proprietà che descrive tale risposta è detta polarizzabilità. Un sistema costituito da particelle
cariche, quando viene sottoposto a un campo esterno, si polarizza. Nel nostro caso, il campo
esterno è oscillante e quindi la polarizzazione stessa risulta variabile nel tempo, generando dipoli
indotti oscillanti. Secondo la teoria di Maxwell, un dipolo oscillante emette radiazione e quindi
la radiazione diffusa osservata sperimentalmente è proprio quella emessa da questi
dipoli indotti; questi dipoli sono proporzionali al campo applicato attraverso la polarizzabilità
dipolare. Assumendo una risposta lineare (ovvero assumendo campi non troppo intensi), il dipolo
indotto è direttamente proporzionale al campo tramite il tensore di polarizzabilità. Formalmente, se
consideriamo una componente del dipolo indotto (lungo x, y o z), essa è legata a una componente
del campo applicato attraverso il tensore di polarizzabilità α, che è di rango 2 (dipende da due
indici, che rappresentano le componenti del campo e del dipolo). Ciò significa che un campo diretto
lungo un asse J può generare una risposta (cioè un dipolo indotto) lungo un asse diverso (x, y o z):
X
(⃗
µind )i = αij ξi ([Link])
i
dove (⃗
µind )i è la componente iesima del momento di dipolo indotto e ξi è la componente iesima del
campo. Per la Raman il termine di interazione nell’hamiltoniano va quindi come:
⃗ ind · ξ⃗rad
ĥint ∝ µ ([Link])
la quale è una forma quadratica rispetto alle componenti del campo. L’interazione sarà esprimibile
come:
1X
αij ξi ξj ([Link])
2 i,j
Questo differenzia la Raman dalle spettroscopie di assorbimento, dove il termine di interazione è
lineare.
Parliamo ora nello specifico della polarizzabilità. Essa è una proprietà specifica della molecola
e dipende dai suoi stati elettronici. Essa può essere anche considerata come proprietà dinamica,
nel qual caso sarà esprimibile come:
X (E el − E (el) )⟨Φ(el) |di |Φ(el) ⟩r ⟨Φ(el) |di |Φ(el) ⟩r
k i 1 k 1 1
αij(ω) = −2 (el) 2
([Link])
k>1 ω 2 − (Ekel − Ei )
76
con di operatore di momento di dipolo elettronico definito come
X
− ri ([Link])
i
con ri posizione dell’elettrone i-esimo. Dunque, risolvendo l’equazione per gli stati adiabatici (cosa
di cui ci occuperemo a breve), saremo in grado di calcolare la polarizzabilità, che ci dice come il
sistema risponde a un campo oscillante con frequenza ω. Abbiamo anche ottenuto che, nel regime
di risposta lineare, quando ωrad si avvicina alla differenza di energia tra il generico stato eccitato e
lo stato fondamentale (Ek − E1 ), la funzione diverge ed otteniamo il fenomeno della risonanza.
Nel caso della spettroscopia Raman, la radiazione incidente del laser è nel visibile, mentre le prime
risonanze elettroniche tipicamente cadono nell’ultravioletto. Pertanto la frequenza ω0 del laser si
trova lontano da queste risonanze. Nonostante ciò, la diffusione avviene ugualmente; per esempio
la diffusione Rayleigh è dovuta proprio alla risposta del sistema alla frequenza ω0 . Un
esempio macroscopico di questo fenomeno è il colore azzurro del cielo: l’azoto e l’ossigeno diffon-
dono selettivamente la radiazione solare bianca, rendendo la luce diffusa prevalente nell’azzurro. A
seconda dell’ora e dell’inclinazione solare la diffusione cambia, modificando il colore [Link]
spettri Raman, oltre alla diffusione Rayleigh, però, osserviamo anche le righe Stokes e anti-Stokes.
Queste dipendono dal fatto che il processo diffusivo coinvolge due fotoni: uno incidente e uno
diffuso, simultanei. Durante questo processo, lo stato nucleare della molecola può cambiare e ciò
comporta che la radiazione diffusa abbia frequenze diverse rispetto alla radiazione incidente. Anche
in questo caso la proprietà fondamentale è il tensore di polarizzabilià, il quale essendo simmetrico,
si può trattare come una matrice e sarà esprimibile come:
3n−6(5)
0
X ∂αij
αij = αij + Qk ([Link])
∂Qk 0
k=1
ed abbiamo che le righe Stokes sono quelle per cui ν → ν + 1 e saranno quelle più intense; le
righe Anti-Stokes sono invece quelle per cui ν → ν − 1 e saranno quelle meno intense visto che
a Tamb solo lo stato vibrazionale fondamentale è popolato.
77
che interagisce con il momento magnetico delle particelle. Quando il fascio attraversa il dispositivo,
si osserva una separazione in due sottofasci distinti. Ciò dimostra che ciascun elettrone possiede un
momento magnetico con due soli possibili valori proiettati lungo la direzione del campo applicato.
In altre parole:
• metà degli elettroni vengono deflessi in una direzione
Attraverso la misura si può determinare il modulo del momento magnetico dell’elettrone e la sua
componente lungo la direzione del campo applicato. Lo spin si comporta quindi come un momento
angolare e lo si può razionalizzare come tale:
⃗ eb = gµb S
µ ([Link])
⃗ eb è
dove g è il fattore di Lande, µb è il magnetone di Bohr, S è l’operatore associato allo spin e µ
il momento magnetico intrinseco dell’elettrone. L’operatore Ŝ è a tutti gli effetti un operatore di
momento angolare, infatti abbiamo che:
1
Ŝz = ± ([Link])
2
in unità atomiche. Siccome noi siamo interessati ad incorporare lo spin nelle nostre funzioni moleco-
lare introduciamo la rappresentazione dello spin. Questa è una rappresentazione nella quale
Ŝz è un operatore moltiplicativo. Definiamo quindi una funzione dello spin η(σ) tale che:
78
siccome gli stati possibili sono solo due il prodotto scalare tra funzioni del genere sarà:
Z σ= 12
X
⟨ηi ||ηj ⟩ = ηi∗ ηj dσ = ηi∗ ηj ([Link])
σ=− 12
β( 12 ) = 0
β(σ) = ([Link])
β(− 12 ) = 1
per cui valgono le seguenti proprietà:
1 1
⟨α(σ)||α(σ)⟩ = α(− )2 + α( )2 = 1 ([Link])
2 2
1 2 1 2
⟨β(σ)||β(σ)⟩ = β(− ) + β( ) = 1 ([Link])
2 2
1 1 1 1
⟨β(σ)||α(σ)⟩ = β(− )α(− ) + β( )α( ) = 0 = ⟨α(σ)||β(σ)⟩ ([Link])
2 2 2 2
le funzioni α, β formano quindi una base ortonormale nella rappresentazione di spin.
Ogni funzione in questa base sarà quindi esprimibile come:
con c1 , c2 coefficienti complessi arbitrari. Possiamo ora vedere come si comportano gli operatori
che ci interessano in questa base. Per costruzione Ŝz è ortogonale:
1
2 0 1 1
Ŝz = (α, β) = ( α, − β) ([Link])
0 − 12 2 2
Siamo quindi interessati ad Ŝx , Ŝy , Ŝ 2 . Sappiamo che le proiezioni non possono essere diagonali
(dalle regole di commutazione del momento angolare) e quindi per scriverli usiamo gli operatori di
salita e discesa:
Ŝ± = Ŝx ± iŜy ([Link])
Dimostriamo che operatori di questa forma agiscono sui ket di Ŝz facendo scendere (o
salire) di un ”gradino” la proiezione del momento di spin:
usando infatti le proprietà di commutazione del momento angolare si vede che Eq. [Link] può
essere riscritta come:
[Ŝz , Ŝ± ] = iŜy ± i(−iŜx )
([Link])
= iŜy ± Ŝx = ±Ŝ±
troviamo ora l’effetto di Ŝ± su un autoket di Ŝz considerando la seguente equazione:
79
questa sarà riscrivibile come:
Ŝz (Ŝ± |s, ms ⟩) = [Ŝz , Ŝ± ] + Ŝ± Ŝz |s, ms ⟩ ([Link])
che è proprio quello che volevamo dimostrare. Dimostrato l’effetto di Ŝ± possiamo definire le
seguenti proprietà di questi operatori:
Ŝ+ α = 0 Ŝ+ β = α
Ŝ− α = β Ŝ− β = 0
le matrici di Ŝx , Ŝy nella rappresentazione dello spin saranno esprimibili come:
0 12
Ŝx = 1 ([Link])
0
2
0 − 2i
Ŝy = i ([Link])
2 0
dove abbiamo sfruttato la relazione biunivoca tra un operatore e la sua matrice di rappresentazione
rispetto ad una certa base. Le tre matrici che rappresentano Ŝx , Ŝy , Ŝz nella rappresen-
tazione dello spin sono dette matrici di pauli:
1
0 12 0 − 2i
0
Ŝz = 2 Ŝx = Ŝy = ([Link])
0 − 12 1
2 0 i
2 0
note queste matrici sfruttando lo stesso principio usato in Eq. [Link] possiamo ottenere la matrice
di Ŝ 2 in questa rappresentazione:
3 1 0
Ŝ 2 (α, β) = (Ŝx2 + Ŝy2 + Ŝz2 )(α, β) = (α, β) ([Link])
4 0 1
80
fondamentale della simmetria permutazionale. In meccanica quantistica non possiamo più
associare ad una particella una traiettoria ben definita e quindi questa diventa fondamentelmente
indistinguibile dalle altre particelle dello stesso tipo nel sistema. Matematicamente deve essere che
per qualsiasi permutazione di due particelle il modulo quadro delle funzioni d’onda che descrive il
sistema deve rimanere invariato:
In teoria, si potrebbe anche pensare anche ad un fattore complesso di modulo unitario davanti alla
funzione, tuttavia usando la teoria dei gruppi di permutazione si riesce ad escludere questa opzione.
Fisicamente esistono quindi due tipi di particelle:
• Fermioni: particelle a spin semi-intero che seguono la statistica di fermi-dirac, queste devono
essere antisimmetriche rispetto allo scambio di due indici.
• Bosoni: particelle a spin intero che seguono la statistica di bose-einstein, queste devono essere
simmetriche rispetto allo scambio di due indici.
introducendo un operatore di scambio possiamo scrivere per una funzione d’onda di fermioni
che:
P̂ Φ(1, 2, 3, ...., n) = (−1)k Φ(1, 2, 3, ...., n) ([Link])
con k numero di scambi eseguiti. La funzione d’onda quindi, per essere fisicamente rilevante, non
solo deve essere autofunzione degli operatori di spin ma deve anche soddisfare Eq. [Link]. For-
tunatamente non sarà difficile far si che la nostra funzione soddisfi questi vincoli. Concentriamoci
prima sulla parte strettamente riguardante lo spin. Per un sistema di n elettroni avremo 2n autofun-
zioni di spin che formano una base ortonormale completa. Queste funzioni possono essere espresse
come combinazione lineare delle funzioni di spin monoelettroniche. Siccome il numero di funzioni
di base aumenta esponenzialmente all’aumentare degli elettroni cerchiamo un modo intelligente di
costruirle. Per farlo definiamo gli operatori associati a stati di spin polielettronici:
⃗ˆtot =
Xˆ
S ⃗i
S ([Link])
i
Xˆ
Ŝz,tot = ⃗z,i
S ([Link])
i
⃗ˆtot · S
⃗ˆtot =
XX
Ŝ 2 = S Ŝi · Ŝj ([Link])
i j
da queste relazioni vediamo che mentre Ŝz,tot è semplicemente una somma di operatori monoelet-
tronici, per quanto riguarda Ŝ 2 abbiamo una somma sia di operatori monoelettronici (i = j) che di
operatori a due elettroni (i ̸= j). Vediamo che Ŝ 2 può essere scritto in funzione degli operatori di
81
scambio:
n
X X
Ŝ 2 = Ŝi2 + 2 Ŝj · Ŝk
i j<k
([Link])
3 X
= nI + 2 Ŝj · Ŝk
4
j<k
ci basta quindi legare il secondo termine all’operatore di scambio, per farlo cerchiamo due valori
tali che:
Ŝj · Ŝk ≡ B P̂jk + CI ([Link])
dove usiamo la forma B P̂jk − C perché sappiamo già che i due operatori non hanno stessi
autovalori e quindi sarebbe inutile cercare solo B. Per trovare le due costanti ci basta
sfruttare che per i due elettroni k, j abbiamo solo due possibili stati: un singoletto e un tripletto.
Calcoliamo quindi gli autostati:
(λ) 1 2 (λ) 1 3 3 3
(Ŝj · Ŝk )Θ0,2 = (Ŝ − Ŝj2 − Ŝk2 )Θ0,2 = (0 − − )=−
2 2 4 4 4 ([Link])
(λ) 1 (λ) 1 3 3 1
(Ŝj · Ŝk )Θ1,2 = (Ŝ 2 − Ŝj2 − Ŝk2 )Θ1,2 = (2 − − )=
2 2 4 4 4
dove abbiamo usato Eq. [Link] nel caso di due elettroni. Siccome sappiamo poi che i due stati
sono rispettivamente antisimmetrici e simmetrici rispetto allo scambio; otteniamo un sistema di
due equazioni in due variabili: (
B + C = 14
([Link])
−B + C = − 34
che tramite semplici sostituzioni ci porta ad:
1 1
Ŝj · Ŝk ≡ − I + Pjk ([Link])
4 2
che, sostituito in Eq. [Link] insieme ad Eq. [Link], ci da:
n
⃗ˆtot · S
⃗ˆtot =
XX 3 X 1 1
Ŝ 2 = S Ŝi · Ŝj = nI + 2 − I + Pjk =
i j
4 4 2
j<k
3 1 X
= n − n(n − 1) I + P̂jk ([Link])
4 4
j<k
n
3 1 X
Ŝ 2 = n − n(n − 1) δil + P̂jk
il 4 4
j<k
permettendoci quindi di riscrivere Ŝ 2 solo in termini di operatori di scambio. Questa forma è molto
utile poiché connette direttamente lo spin con il gruppo delle permutazioni ma in questo
corso non andremo più a fondo di cosı̀. In ogni caso questa relazione semplifica anche i calcoli
a mano:
82
• n = 3, ms = 12 : Abbiamo 3 funzioni possibili:
S12 = 1
+ 12
1
s1 = 2
− 12
S12 = 0
aggiungendo ancora un altro elettrone possiamo ottenere 3 stati diversi di spin; un quartetto
S = 34 e due doppietti S = 12 :
83
3
+ 21 S= 2
S12 = 1
+ 12
− 21 S= 1
2
1
s1 = 2 + 21
− 12
S12 = 0
Gli stati di spin per n = 4 sono ottenibili in modo analogo a quanto visto per n = 2 → n = 3,
partendo da quelli ottenuti per n = 3. In questo modo possiamo ”analiticamente” calcolare
gli stati di spin per n elettroni. Per ottenere invece le funzioni di spin associate ad ogni stato
dobbiamo fare una considerazione, ovvero che preso uno stato di spin S per n elettroni sono due
gli stati che ci interessano:
• (S + 12 ,n + 1)
• (S − 12 ,n + 1)
Possiamo ottenere il primo moltiplicando per α:
che poi andrà normalizzata. Quello che abbiamo fatto è creare, partendo da ΘS,n β, una funzione
ortogonale. Le funzioni cosı̀ ottenute sono quindi già ortonormali e per ogni stato con un
certo D(S, n) abbiamo esattamente λ funzioni con 1 ≤ λ ≤ D(S, n); ovvero:
(λ)
∀(S, n) ∃ ΘS,n , 1 ≤ λ ≤ D(S, n) ([Link])
per via di ciò sono queste funzioni quelle tendenzialmente usate in calcoli di chimica quantistica.
Concludiamo dicendo che vale sempre che:
X
D(S, n)(2S + 1) = 2n ([Link])
S
Concentriamoci ora invece sulla parte riguardante la simmetria permutazionale. Abbiamo visto
in Eq. [Link] come esprimere l’effetto dello scambio di due indici sulla funzione d’onda; siccome lo
scambio di indici ci ridà la funzione d’onda originale (al netto della costante moltiplicativa) allora
deve essere che:
[Ĥ, P̂ ] = 0 ([Link])
84
ovvero l’hamiltoniano e l’operatore di scambio commutano sempre. In ogni caso noi
richiediamo che la nostra funzione sia antisimmetrica rispetto allo scambio di coordinate e se ri-
solvessimo esattamente l’equazione di Schroedinger dovremmo aggiungere anche questo vincolo.
Siccome noi siamo interessati a soluzioni approssimate ottenute come combinazione di una base di
lavoro, semplicemente scriviamo la nostra funzione di prova come:
con Ŝn gruppo simmetrico. Grazie a questa relazione derivabile dalla teoria dei gruppi possiamo
definire anche un operatore di simmetrizzazione come:
1 X
Ŝ = P̂ ([Link])
N!
P̂ ∈Sn
questo operatore può per esempio essere usato per considerare effetti quantistici per nuclei bosonici
identici. Concludiamo il paragrafo dicendo che la comodità d’uso di questi operatori sta nel fatto
che ”proiettano” la combinazione delle funzioni di lavoro nella loro componente antisimmetrica o
simmetrica (ovvero quelle a cui siamo interessati) e, nel caso in cui questa non vi sia, rendono la
funzione nulla, permettendoci quindi di considerare il principio di Pauli direttamente nella
formulazione matematica senza doverne tener conto direttamente.
dove abbiamo usato Eq. [Link] ed abbiamo posto xi = (ri , σi ). Si vede che Eq. [Link] è
equivalente al determinante di una matrice (n × n) (Eq. [Link]) che ha come elementi:
ovvero come indice di riga abbiamo le coordinate dell’i-esimo elettrone, mentre come indice di
colonna lo j-esimo spin-orbitale. Dalle proprietà dei determinanti abbiamo che Ψtrial
k sarà nulla se:
85
• Righe uguali: avremmo in questo caso che xi = xj , che è in teoria permesso ma nullo nel
caso di fermioni. Questa condizione viene definita nodo fermionico.
• Colonne uguali: avremmo in questo caso che due elettroni nello stesso stato con uguali xi .
Questo è strettamente vietato per il principio di esclusione di Pauli.
e poiché tendenzialmente nei calcoli numerici si usano più orbitali di quanti elettroni possieda il
sistema abbiamo dalla combinatoria:
m
numero di possibili determinanti = ([Link])
n
con m ≥ n numero di orbitali della base di lavoro e n numero di elettroni. Che in generale è un
numero grande; vediamo però che mentre vale sempre per costruzione che:
1
Ŝz Dk = (Nα − Nβ )Dk ([Link])
2
in generale non vale che:
Ŝ 2 Dk = S(S + 1)Dk ([Link])
con Dk un generico determinante di Slater e Nα , Nβ numero di spin α, β rispettivamente. Si cercherà
quindi solo funzioni con una certo valore di ms e il numero di determinanti da calcolare sarà minore:
m m
numero di determinanti con ms fissato = ([Link])
Nα Nβ
Tuttavia anche questo è un numero generalmente grande e quindi più spesso si costru-
iscono funzioni usando il secondo metodo introdotto.
Una funzione di configurazione (CFS) è definita come:
Y n
(λ)
Ψtrial,CSF
k = Â ϕ(ri )ki ΘS,n (σ1 , σ2 , ..., σn ) ([Link])
i=1
(λ)
dove ϕ(ri )ki è una funzione spaziale e ΘS,n (σ1 , σ2 , ..., σn ) è una funzione di spin ottenuta per esempio
con Yamanouchi-kotani. Queste funzioni sono comode poiché hanno ”inserito” il principio di
esclusione di pauli nella parte spaziale e sono, per costruzione, autofunzioni sia di Ŝz
che di Ŝ 2 . Da un punto di vista matematico possono essere viste come combinazioni lineari di
determinanti di Slater. Per poter organizzare meglio le CFS ed evitare ripetizioni si possono
usare le tabelle di Weyl:
86
All’interno delle celle si inseriscono gli orbitali secondo un ordine preciso e univoco. Questo ordine
segue la disposizione con cui gli orbitali compaiono nel prodotto complessivo della funzione d’onda:
si procede riga per riga da sinistra verso destra, e successivamente dall’alto verso il basso fino a
completare la tabella. In tal modo, ogni configurazione è definita in maniera univoca. Nel caso di
celle adiacenti, si assegnano degli indici orbitali i, j, k, ecc., rispettando la condizione:
i ≤ j, i < k
i j ([Link])
k
lungo la colonna, cosı̀ da evitare violazioni del principio di Pauli, che proibisce la presenza di due
elettroni con gli stessi numeri quantici. Inoltre, questa regola garantisce anche che non si generino
ripetizioni o combinazioni triviali, assicurando quindi l’indipendenza lineare delle funzioni costruite.
Prendiamo il caso di n = 3, m = 3; saranno possibili un quartetto e un doppietto:
• Quartetto (S = 32 ):
1
2
3
• Doppietto (S = 12 ):
1 1 1 2 1 3 1 1 1 2 1 3
2 2 2 3 3 3
2 2 2 3
3 3
1 3
Ogni tabella è associata ad una CSF. Abbiamo ottenuto 8 CSF per S = 2 e 1 per S = 2,
mentre con un determinante di slater per Sz = 21 avremmo ottenuto:
3 3
Nslater = =9 ([Link])
2 1
Ovvero con le CSF otteniamo direttamente le funzioni con Sz fissato senza bisogno di
separarle dalle altre (con slater invece le otteniamo tutte). Quando il numero di orbitali e quello
di elettroni diventa elevato è quindi sempre preferibile lavorare con le CSF. Il numero di CSF (anche
detto numero di Weyl) è dato da:
2S + 1 M + 1 M +1
NCSF = ([Link])
M + 1 n2 − S M − n2 − S
Che è generalmente minore di Nslater . In conclusione, i determinanti di Slater forniscono i mattoni
di base per costruire funzioni elettroniche, mentre le CSF permettono di ottenere funzioni con spin
ben definito, indispensabili per una modellistica accurata.
Con tutte queste nozioni è possibile scrivere correttamente una funzione di prova che
poi dovrà essere ottimizzata, come vedremo nel possimo paragrafo.
87
2.6 Calcolo degli stati adiabatici
Passiamo quindi al calcolo vero e proprio degli stati adiabatici. Ricordiamo che abbiamo sempre a
disposizione una base completa di funzioni di spin e quindi dovremmo preoccuparci solo della
parte spaziale della funzione d’onda. Scriveremo sempre la funzione polielettronica da calcolare
come:
Xn
Ψ= c k Φk ([Link])
k
dove Φk sono generiche funzioni di lavoro, riscrivibili come una qualche combinazione lineare di
funzioni orbitaliche:
Xm
Φk = tlk χl ([Link])
l
queste funzioni orbitaliche (orbitali atomici o molecolari per esempio) sono a loro volta riscritti come
combinazione di funzioni elementari, ovvero delle funzioni ”semplici” (tipicamente gaussiane) che
permettano di rendere gli integrali da calcolare computazionalmente analitici:
h
X
χl = fe dle ([Link])
e
i parametri da ottimizzare sono quindi molti e dispersi tra i vari livelli di organizzazione. Questi
parametri vanno ottimizzati per ottenere una funzione migliore possibile e per farlo usiamo il teo-
rema variazionale:
Theorem 2.1. Sia Ψ una funzione approssimata normalizzata per l’hamiltoniano, allora vale sem-
pre che Ĥ ⟨Ψ|Ĥ|Ψ⟩ ≥ Egs .
Dimostrazione: prendiamo le autofunzioni esatte dell’hamiltoniano in esame:
˜
Ĥ Φ̃k = Ek Φ̃ ([Link])
che può essere al massimo uguale allo stato fondamentale nel caso k = 1 per via di come abbiamo
ordinato gli stati.
88
Per ottimizzare i parametri basta quindi derivare l’energia rispetto ai parametri delle
funzioni d’onda e trovare i punti stazionari. Tendenzialmente il teorema variazionale è usato
per ottimizzare i parametri delle funzioni elementari che vengono scelte come base per descrivere le
funzioni orbitaliche. Spesso comunque piuttosto che ricercare questi parametri ogni volta, si usano
delle basi già ottimizzate con parametri ad hoc per meglio approssimare il sistema in esame.
I risultati ottenuti per lo stato fondamentale con il teorema variazionale possono essere gener-
alizzati usando il teorema variazionale lineare (o di Ritz). Supponiamo di avere la seguente
funzione di prova: X
Ψ= c n Φn ([Link])
n
dove {Φn } sono generiche funzioni orbitaliche (ottenibili come determinanti di slater o CSF e
calcolate da funzioni elementari). Settiamo il seguente sistema di equazioni che vogliamo risolvere:
min ⟨Ψ|Ĥ|Ψ⟩
([Link])
⟨Ψ||Ψ⟩ = 1
questo è un problema di Lagrange e per risolverlo prendiamo variazioni infinitesime dei coefficienti:
ck → ck + δck ([Link])
la lagrangiana diventerà quindi:
δ ⟨Ψ|Ĥ|Ψ⟩ − λδ ⟨Ψ||Ψ⟩ − 1 = 0
X X ([Link])
∗ ∗ ∗ ∗
(δck cl + ck δcl )Hkl − λ (δck cl + ck δcl )Skl = 0
k,l kl
dove i prodotti tra variazioni infinitesime non sono considerati e abbiamo che Hkl = ⟨Φk |Ĥ|Φl ⟩, Skl =
⟨Φk ||Φl ⟩. Raccogliendo i termini otteniamo:
X X
δc∗k (cl Hkl − λcl Skl ) + δcl (c∗k Hkl − λc∗k Skl ) = 0 ([Link])
kl kl
dove sono state sfruttare le proprietà degli operatori Hermitiani per S, H. I due termini sono
quindi l’uno il complesso coniugato dell’altro e, siccome le variazioni sulla parte reale e su
quella complessa sono indipenti tra loro possiamo scrivere per una generica variazione:
δck (K) = Kδa + Kiδb ([Link])
∗
δck (K ) + δc∗k (K) ∗ ∗ ∗
= K δa + K iδb + Kδa − Kiδb = δa(K + K) + iδb(K − K) ∗
([Link])
Siccome stiamo cercando un minimo le poniamo uguale a zero ed otteniamo un sistema di due
equazioni: (
K ∗ = −K
([Link])
K∗ = K
89
che è valido se e solo se K = K ∗ = 0. Il sistema da risolvere diventa quindi:
(
K=0
([Link])
K∗ = 0
che sono a loro volta le rispettive complesse coniugate. Se una è vera sarà vera anche l’altra
e quindi ci basta quindi risolvere: X
(Hkl − λSkl )cl = 0 ([Link])
kl
o in forma matriciale:
Hc = λSc ([Link])
ma per definizione la matrice S è definita positiva e simmetrica e quindi può essere diagonalizzata
con una trasformazione unitaria:
†
U SU = Λ ([Link])
con (Λ)kk > 0. Eq. [Link] può essere riscritta come:
diagonalizzando H′ trovo il set di autovalori e autovettori {λk , c(k) }. In ogni caso siccome spesso da
un punto di vista computazionale conviene di molto usare funzioni ortogonali d’ora in poi useremo
semplicemente:
Hc = λc ([Link])
Il valor medio dell’energia, per un certo autovettore, sarà esprimibile come:
X (λ)∗ (λ) X (λ)∗ (λ) X (λ)
< E(λ) >= ck cl Hkl = ck cl δkl λ = |ck |λ = λ ([Link])
kl kl k
Partendo da questa relazione possiamo dimostrare un corollario importante del teorema variazionale
lineare, ovvero il teorema di Mcdonald, grazie al quale potremo estendere i risultati prima
ottenuti per stati adiabatici eccitati. Il teorema va dimostrato per induzione; supponiamo quindi
di avere due funzioni approssimate: X (1)
Ψ1 = ck Φk
k
X (2)
([Link])
Ψ2 = ck Φk
k
90
e invece due funzioni esatte non degeneri:
Ĥ Ψ̃1 = E1 Ψ̃1
([Link])
Ĥ Ψ̃2 = E2 Ψ̃2
Abbiamo dimostrato dal teorema variazionale il caso base, ovvero < Egs >= λ1 ≥ Egs . Diciamo
ora di poter scrivere la funzione del primo stato eccitato esatto Ψ̃2 come combinazione lineare di
quelle approssimate e scriviamo i vincoli del sistema:
(
⟨aΨ1 + bΨ2 ||Ψ̃1 ⟩ = 0
([Link])
a2 + b2 = 1
⟨aΨ1 + bΨ2 |Ĥ|aΨ1 + bΨ2 ⟩ = a2 < E1 > +b2 < E2 >= a2 λ1 + b2 λ2 ([Link])
dove i termini misti si annullano per ortogonalità ed abbiamo usato Eq. [Link]. La funzione esatta
sarà quindi approssimabile come:
X
Ψ̃2 = aΨ1 + bΨ2 = Dk Ψ̃k ([Link])
k>1
Ottenendo:
a2 λ1 + b2 λ2 + c2 λ3 ≥ E3
a2 (λ1 − λ3 ) + b2 (λ2 − λ3 ) + λ3 ≥ E3 ([Link])
2 2
λ3 ≥ E3 + a (λ3 − λ1 ) + b (λ3 − λ2 )
Si vede facilmente quindi che nel caso generale otteniamo:
k−1
X
λk ≥ Ek + a2l (λk − λl ) ([Link])
l=1
91
confermando che ciascun autovalore esatto fornisce un limite inferiore al corrispondente autovalore
approssimato. I tre teoremi quà visti: il teorema variazionale ristretto, il teorema vari-
azionale lineare ed il teorema di Mcdonald sono considerabili le basi di ogni metodo di
chimica quantistica. Nel prossimo paragrafo mostreremo quindi il modo più semplice di sfruttare
questi risultati per cercare funzioni d’onda accettabili.
Per analizzare in modo accurato questi processi in realtà bisognerebbe affrontare un problema di
scattering, ma possiamo limitarci a una considerazione più semplice. Se un elettrone ha una certa
energia di legame nella materia, allora possiamo pensare a un modello analogo a quello dell’atomo
di idrogeno; in cui ogni elettrone si comporta come una particella indipendente ”oc-
cupando” un orbitale idrogenoide. Questa approssimazione è infatti il modello strettamente
usato nella chimica per giustificare tutte le considerazioni sulla configurazione fatte per esempio in
chimica organica. Ogni elettrone comunque è soggetto a un potenziale efficace che rappresenta
l’effetto medio degli altri elettroni. Il sistema sarà, dunque, soggetto a un’equazione simile a quella
dell’atomo di idrogeno:
1
[− ∇2 − vSCF ] ϕj (⃗r) = Ej ϕj (⃗r) ([Link])
2
Il potenziale vSCF viene indicato come potenziale autocompatibile. In tale schema, ciascun
elettrone occupa una funzione d’onda analoga a quella dell’orbitale dell’atomo di idrogeno, e la sua
energia Ej risulta pari a meno l’energia di legame corrispondente. Dai profili di fotoemissione si
possono cosı̀ dedurre le energie di questi livelli. La funzione d’onda ϕj di ciascun elettrone rappre-
senta quindi il suo stato indipendente nel campo medio prodotto dagli altri elettroni.
La densità elettronica complessiva del sistema risulta dalla somma dei contributi di tutti gli
elettroni:
N
X
ρ(⃗r) = |ϕi (⃗r)|2 . ([Link])
i=1
92
Il quadrato del modulo dell’orbitale atomico (come nel caso dell’idrogeno) fornisce la distribuzione
spaziale dell’elettrone; la somma dei moduli al quadrato di tutte le funzioni rappresenta la densità
totale del sistema.
Questo modello, tuttavia è nella realtà approssimato poiché gli elettroni interagiscono tra loro.
Questo potenziale vSCF (self-consistent field ) contiene naturalmente il termine di interazione tra
gli elettroni e i nuclei. Possiamo quindi scrivere una somma estesa a tutte le cariche nucleari nel
modo seguente:
n
X Za
vSCF = − ([Link])
⃗ a|
a=1 |⃗
r−R
Tuttavia, in un sistema polielettronico sono presenti anche tutti gli altri elettroni, i quali con-
tribuiscono a un effetto complessivo di interazione elettrostatica classica. Tale effetto può essere
rappresentato come una somma estesa su tutte le posizioni elettroniche ⃗r′ , nella quale compare la
densità elettronica:
n Z
X Za ρ(⃗r)
vSCF = − d⃗r ([Link])
a=1 |⃗
r−R ⃗ a| |⃗r − ⃗r′ |
È importante notare che questa densità elettronica rappresenta una densità numerica, ossia il
numero di elettroni per unità di volume, e non una densità di carica. Per ottenere la densità di
carica, è sufficiente moltiplicare per la carica elementare −e dell’elettrone.
Passiamo ora alla parte di calcolo vera e proprio. Nel metodo di Hartree-Fock, la funzione d’onda
elettronica viene espressa nella forma più semplice possibile adatta a un modello di fermioni in-
dipendenti; ovvero come prodotto di spin-orbitali monoelettronici:
h i
(1) (2) (N )
ΦHF (⃗r1 , σ1 , ⃗r2 , σ2 . . . , ⃗rN , σN ) = N Â ψ1 ψ2 · · · ψN ([Link])
93
DObbiamo quindi ricercare gli spin-orbitali ottimali in cui collocare gli N elettroni, in modo da
ottenere il miglior determinante di Slater possibile, ossia quello che minimizza l’energia totale del
sistema. Questa minimizzazione avviene applicando il teorema variazionale:
dove Ĥe è l’operatore hamiltoniano elettronico. Per semplificare il calcolo, si impone la condizione
che gli spin-orbitali siano ortogonali tra loro:
94
L’uso di spin-orbitali ortogonali risulta a questo punto estremamente vantaggioso: se gli orbitali
non fossero ortogonali, la somma su N ! permutazioni comporterebbe il calcolo di un numero enorme
di integrali, mentre grazie all’ortogonalità molti termini si annullano, riducendo drasticamente il
lavoro computazionale.
Prima di andare avanti ricordiamo ora la struttura di Ĥe :
N n
1 X 2 XX Za X 1 X Za Zb
Ĥe = − ∇i − + + ([Link])
2 i ⃗ a| r Rab
i=1 a=1 |⃗
ri − R i<j ij a<b
dove:
• il primo termine è l’operatore monoelettronico, che comprende il termine cinetico
• il secondo termine rappresneta l’interazione con i nuclei
• 1
rij rappresenta la repulsione coulombiana tra gli elettroni i e j,
• l’ultimo è il termine di repulsione tra i nuclei, che non dipende dalle coordinate elettroniche e
può quindi essere considerato costante e alla fine all’energia, fungendo quindi da riferimento
costante.
È importante notare che la prima parte dell’hamiltoniano (il primo e il secondo termine) contiene
operatori a un corpo, nel senso che è una somma sugli elettroni;
N
X X 1 X Za Zb
Ĥe = ĥi (⃗r) + + ([Link])
i=1
r
i<j ij
Rab
a<b
con:
N n
1 X 2 XX Za
ĥ(⃗r) = − ∇i − ([Link])
2 i ⃗ a|
i=1 a=1 |⃗
ri − R
Si tratta di un operatore a un elettrone nel campo di tutti i nuclei. Al contrario, la seconda parte
contiene operatori a due corpi. In generale, la presenza di operatori a tre corpi o di ordine
superiore renderebbe la trattazione molto più complessa, pertanto in Hartree-Fock ci
si limita a quelli a uno e due corpi.
Possiamo scrivere il valore di aspettazione ottenuto in Eq. [Link] per l’energia come:
X (1) (N )
X X 1 (1) (N )
EHF = (−1)P ⟨ψ1 . . . ψN | ĥ(⃗ri ) + |P̂ ψ1 . . . ψN ⟩ ([Link])
i
r
i<j ij
P̂ ∈SN
Andiamo a studiare per ora soltanto i termini che vengono fuori dall’operatore ĥ1 , avremo:
n
Y n
Y
⟨ψ1 |ĥ(1)|ψ1 ⟩ ⟨Ψi ||Ψi ⟩ − ⟨ψ1 |ĥ(1)|ψ2 ⟩⟨ψ2 |ψ1 ⟩ ⟨Ψi ||Ψi ⟩ + . . . ([Link])
i=2 i=3
Vediamo che per ortogonalità sopravvive solo il primo termine. Questo ragionamento ovviamente
vale per ogni ĥ(⃗ri ) e quindi la parte mono-elettronica nel metodo di Hartree-Fock è:
n
X
EHF = hii ([Link])
i=1
95
dove:
⟨ψ1 |ĥ(1)|ψ1 ⟩ = h11 ([Link])
In altre parole, ciascun elettrone contribuisce indipendentemente con la propria energia.
monoelettronica nel campo dei nuclei.
Passiamo ora alla parte più complessa: il contributo bielettronico associato all’operatore di
interazione coulombiana tra coppie di elettroni. Procediamo in modo simile a quanto fatto per il
caso monoelettronico:
(1) (N ) 1 (1) (N )
⟨ψ1 . . . ψN | |(1 − P̂12 )ψ1 . . . ψN ⟩ ([Link])
r12
Dove abbiamo usato l’ortogonalità delle funzioni d’onda, la quale fa si che tutte le permutazioni
che non scambiano esattamente 1 e 2 si annullano. Otteniamo quindi:
1 1
⟨ψ1 ψ2 | |ψ1 ψ2 ⟩⟨ψ3 |ψ3 ⟩ . . . ⟨ψN |ψN ⟩ − ⟨ψ1 ψ2 | |ψ2 ψ1 ⟩⟨ψ3 |ψ3 ⟩ ([Link])
r12 r12
In questa espressione compaiono due tipi di integrali. Questi integrali, detti integrali bielettron-
ici, coinvolgono le coordinate spaziali e di spin degli elettroni 1 e 2 (dove con 1 e 2 non indichiamo
il primo e il secondo elettrone ma semplicemente i due interessati dalla permutazione). Durante
l’integrazione, si distinguono le variabili integrate: xi = (ri , σi ), con l’integrale completo espresso
come: Z Z
dx1 dx2 . . . dxN = dr1 dσ1 dr2 dσ2 · · · drN dσN . ([Link])
Quando interviene l’operatore di scambio P̂12 , esso modifica l’ordine degli argomenti ma non annulla
l’integrale, poiché la variazione non comporta l’ortogonalità delle funzioni coinvolte (a differenza
dei termini monoelettronici). In generale, si indica un integrale bielettronico nella forma:
ZZ
1 1
(ψa ψb | |ψc ψd ) = (ψa∗ ψc )(1) (ψ ∗ ψd )(2) dx1 dx2 ([Link])
r12 |⃗r1 − ⃗r2 | b
dove indichiamo che le funzioni a sinistra sono integrate sulle coordinate dell’elettrone 1 e quelle
a destra sull’elettrone 2. L’integrale bielettronico è certamente il collo di bottiglia com-
putazionale dei calcoli Hartree–Fock: è infatti il termine che, storicamente, negli anni ’50
limitava severamente i calcoli fattibili. Questi integrali coinvolgono due particelle e, dopo aver
separato la parte di spin (che si gestisce analiticamente), rimangono integrali spaziali su quattro
coordinate: si tratta quindi di integrali in sei dimensioni per ogni coppia, che senza i computer
moderni erano pressoché proibitivi da calcolare in grande numero. Introduciamo ora la notazione
compatta per gli integrali bielettronici (la cosiddetta “doppia barra” ).
1 1
⟨ψi ψj ||ψi ψj ⟩ = ⟨ψi ψj | |ψi ψj ⟩ − ⟨ψi ψj | |ψj ψi ⟩ ([Link])
r12 r12
Con questa notazione compatta l’energia elettronica di Hartree-Fock in forma finale è della forma:
N N
X 1X
EHF = hii + ⟨ψi ψj ||ψi ψj ⟩ ([Link])
i=1
2 i,j
Per comodità di notazione e per evitare doppi conteggi, si introduce spesso un fattore 1/2 nella
somma sulle coppie: si somma su i e j da 1 a N e si moltiplica per 1/2 (oppure si somma su
96
i < j senza il fattore 1/2). Se si includono anche i termini con i = j la loro combinazione si
annulla opportunamente, quindi la scrittura estesa non è problematica purché si tenga presente la
cancellazione formale. In ogni caso il primo termine rappresenta la somma dei contributi
monoelettronici, mentre il secondo incorpora la media delle interazioni coulombiane e
di scambio tra tutte le coppie di elettroni.
Ora concentriamoci sulla procedura variazionale. L’obbiettivo è minimizzare l’energia sotto il
vincolo di ortonormalità degli spin-orbitali. Imponiamo la varizione arbitraria degli spin-orbitali
(che, in generale, possono essere complessi, quindi si varia indipendentemente la parte reale e
immaginaria):
ψi → ψi + δψi ([Link])
Con il metodo dei moltiplicatori di lagrange si ottiene quindi:
" #
X
δ ⟨EHF ⟩ − λij ⟨ψi |ψj ⟩ − 1 = 0, ([Link])
i,j
dove λij sono i moltiplicatori di Lagrange, che saranno tanti quanti sono i vincoli. Nella pratica si
scrive l’espressione della variazione espandendo i termini monoelettronici e bielettronici. La parte
relativa alla repulsione nucleare ENN non dipende dalle coordinate elettroniche, perciò scompare
nella variazione e può essere aggiunta a posteriori come costante di riferimento dell’energia totale.
Esplicitando la variazione otteniamo termini del tipo:
X 1 X
⟨δψi |ĥ|ψi ⟩ + ⟨ψi |ĥ|δψi ⟩ + ⟨δψi ψj ||ψi ψj ⟩ + ⟨ψi δψj ||ψi ψj ⟩ + ⟨ψi ψj ||δψi ψj ⟩
i
2 i,j
X
+ ⟨ψi ψj ||ψi δψj ⟩ − ⟨δψi |ψj ⟩ + ⟨ψi |δψj ⟩ λij = 0
ij
([Link])
A questo punto, con un ragionamento simile a quanto fatto per la dimostrazione del teorema
variazionale lineare, dividiamo i termini in coppie di complessi coniugati:
X X X
⟨δψi |ĥ|ψi ⟩ + ⟨δψi ψj ||ψi ψj ⟩ − ⟨δψi |ψj ⟩λij + c.c. = 0 ([Link])
i j ij
dove abbiamo usato Eq. [Link]; infatti deriva dal fatto che gli elettroni siano interscambiabili che:
ZZ ZZ
∗ 1 ∗ 1
(δψi ψi )(1) (ψj ψj )(2) dx1 dx2 = (δψi∗ ψi )(1) (δψj∗ ψj )(2) dx1 dx2 ([Link])
|⃗r1 − ⃗r2 | |⃗r1 − ⃗r2 |
che si può verificare facilmente anche per gli integrali di scambio. Otteniamo quindi alla fine:
X X X
⟨δψi |ĥ|ψi ⟩ + ⟨δψi ψj ||ψi ψj ⟩ = ⟨δψi |ψj ⟩λij ([Link])
i j ij
Poiché la variazione è arbitraria per ciascuno spin–orbitale, ogni singolo termine deve annullarsi,
non la somma, ottenendo: X
⟨δψi |ĥ + Jˆ − K̂|ψi ⟩ = ⟨δψi |ψj ⟩λij ([Link])
j
97
ˆ K̂. In ogni caso alla fine diventa:
Mostremo tra poco la definizione di J,
X
(ĥ + Jˆ − K̂)ψi = ψj λji ([Link])
j
F̂ = ĥ + Jˆ − K̂ ([Link])
dove:
• Jˆ è l’operatore coulombiano, che rappresenta l’interazione media dell’elettrone con la
densità elettronica totale.
• K̂ è l’operatore di scambio, privo di analogo classico, derivante dall’indistinguibilità degli
elettroni.
possiamo quindi passare a definire in modo preciso i due operatori.
ˆ
Definizione esplicita di J. L’operatore coulombiano è un operatore moltiplicativo tale per
cui:
ˆ i⟩ =
X 1
⟨δψi |J|ψ ⟨δψi ψj | |ψi ψj ⟩ ([Link])
j
r12
Componendo ora le due parti, possiamo definire in modo piuttosto chiaro l’operatore di Coulomb.
A sinistra abbiamo: Z
ˆ i ⟩ = dx δψi∗ (Jψ
⟨δψi |J|ψ ˆ i )(x) ([Link])
A destra abbiamo:
XZ (δψi∗ ψi )(x1 )(ψj∗ ψj )(x2 )
Z
X 1
⟨δψi ψj | |ψi ψj ⟩ = dx1 dx2 ([Link])
j
r12 j
r12
confrontando le due parti si ottiene che l’azione dell’operatore Jˆ può essere espressa come:
X Z (ψj∗ ψj )(x′ ) Z
ρ(⃗r)
J(x) = dx′ = d⃗r ([Link])
⃗′
|⃗r − r | |⃗r − r⃗′ |
j
98
Definizione esplicita di K̂. Passiamo quindi alla definizione di K̂, detto operatore di scam-
bio. Si chiama cosı̀ perché la sua origine risiede nel processo di antisimmetrizzazione della funzione
d’onda. Se non avessimo antisimmetrizzato il prodotto degli spin–orbitali, il termine di scambio non
comparirebbe, e ci saremmo fermati al solo termine coulombiano. Dunque, il contributo di K̂ nasce
direttamente dal principio di antisimmetria e, quindi, dall’indistinguibilità degli elettroni.
Analogamente a quanto fatto per J, ˆ possiamo scrivere il contributo di K̂ come:
X 1
⟨δψi |K̂|ψi ⟩ = ⟨δψi ψj | |ψj ψi ⟩ ([Link])
j
r12
dove però le funzioni d’onda risultano invertite rispetto al caso coulombiano: l’indice i e j scambiano
le loro posizioni tra le coordinate elettroniche 1 e 2. Questo comporta che K̂ non sia un operatore
moltiplicativo, bensı̀ un operatore integrale, ovvero un operatore definito mediante un kernel
K(x, x′ ) tale che: Z
(K̂f )(x) = K(x, x′ ) f (x′ ) dx′ ([Link])
cioè il valore dell’operatore applicato a una funzione f nel punto x è dato da un integrale sullo
spazio x′ pesato dal kernel. Pertanto, il termine di variazione ⟨δψi |K̂|ψi ⟩ si scrive come:
ZZ
⟨δψi |K̂|ψi ⟩ = δψi∗ (x) K(x, x′ ) ψi (x′ ) dx dx′ ([Link])
99
Questo effetto è stabilizzante, perché obbliga gli elettroni con lo stesso spin a disporsi più lontano
gli uni dagli altri.
Al contrario, gli elettroni con spin opposto non sono soggetti a tale effetto e possono quindi sovrap-
porsi spazialmente. Nel formalismo di Hartree–Fock, al di là del termine coulombiano mediato,
l’interazione tra elettroni di spin opposto rimane approssimata. Gli elettroni con spin opposto pos-
sono sovrapporsi nella stessa regione di spazio: la loro interazione rimane puramente coulombiana
e mediata. È proprio da questa differenza che nascerà la correlazione elettronica, che vedremo
successivamente e che non può essere descritta dal solo metodo di Hartree-Fock e richiede approcci
più sofisticati per essere trattata in modo accurato.
F̂ = ĥ + Jˆ − K̂ ([Link])
e dobbiamo ora risolvere la seguente equazione per ognuno degli orbitali che compone la base:
X
F̂ ψi = ψi λji ([Link])
j
I moltiplicatori di Lagrange, indicati con λij , costituiscono una matrice quadrata n × n, poiché
coinvolgono solo gli n spin-orbitali che compaiono nel determinante di Slater. La Matrice di
Lagrange Λ è la rappresentazione matriciale dell’operatore di Hartee-Fock. Possiamo
scrivere:
F̂ (ψ1 , ψ2 , . . . , ψN ) = (ψ1 , ψ2 , . . . , ψN ) Λ ([Link])
L’operatore di Fock F̂ è Hermitianiano, ma non lineare. Infatti, poiché contiene Jˆ e K̂, che
dipendono dalle soluzioni stesse, non soddisfa la proprietà di linearità. Per questo motivo,come
abbiamo già detto, F̂ non corrisponde a un’osservabile fisica. Comunque, anche se non lineare,
vediamo cosa otteniamo diagonalizzando la matrice Λ tramite una trasformazione unitaria:
U † ΛU = ε ([Link])
Otteniamo la matrice ε, la quale è una matrice diagonale i cui elementi εi rappresentano gli autoval-
ori orbitali, ossia le energie degli spin-orbitali canonici. Dobbiamo ora domandarci quale effetto
100
tale trasformazione abbia sulla funzione d’onda totale del sistema. In particolare, possi-
amo chiederci se la funzione d’onda determinante, costruita a partire dagli spin-orbitali originali,
cambi o meno forma dopo la trasformazione.
Sia Ψ il determinante di Slater costruito dagli spin-orbitali {ψi }, e Ψ′ quello costruito dai nuovi
spin-orbitali {ψi′ } ottenuti tramite una trasformazione unitaria U :
X
ψi′ = Uji ψj ([Link])
j
P P P
U ψ (1) U ψ (1) . . . U ψ (1)
Pj j1 j Pj j2 j Pj jN j
1 j Uj1 ψj (2) j Uj2 ψj (2) . . . j UjN ψj (2)
Ψ′ = √ .. .. .. .. ([Link])
N! . . . .
P P P
j Uj1 ψj (N ) j Uj2 ψj (N ) . . . j UjN ψj (N )
Utilizzando le proprietà del determinante, possiamo scrivere il determinante sopra come il prodotto
del determinante della matrice originaria (costruita dagli spin-orbitali ψj ) per il determinante della
matrice unitaria U :
Ψ′ = (det U ) Ψ ([Link])
Poiché U è unitaria, il suo determinante ha modulo unitario:
| det U | = 1 ([Link])
Ne consegue che la densità elettronica e, più in generale, ogni grandezza fisica dipendente
dal modulo quadrato della funzione d’onda, rimane invariata sotto una trasformazione
unitaria. La trasformazione unitaria che diagonalizza Λ (e quindi F̂ ) non altera la funzione d’onda
se non per un fattore di fase globale, privo di significato fisico. In altre parole, la probabilità di
trovare gli elettroni in una determinata configurazione spaziale resta la stessa. Possiamo quindi
concludere che la diagonalizzazione dell’operatore di Fock è una procedura lecita e conveniente, in
quanto consente di esprimere le equazioni di Hartree–Fock in una forma semplificata, detta forma
canonica, senza modificare la descrizione fisica del sistema. Abbiamo quindi dimostrato che nella
101
rappresentazione canonica, l’operatore di Fock assume la forma: con εi elemento della matrice
diagonale:
U † ΛU = ε ([Link])
Inoltre è possibile dare un significato fisico ai moltiplicatori di lagrange λi attraverso il
teorema di Koopmans.
Consideriamo un sistema costituito da N elettroni, la cui energia totale è EN . Rimuovendo un
elettrone da un determinato spinorbitale ϕi , otteniamo un nuovo sistema con N − 1 elettroni, di
(i)
energia EN −1 . L’energia richiesta per tale rimozione è:
(i)
Ii = EN −1 − EN ([Link])
Nel quadro approssimato di Hartree–Fock, tale differenza è pari, in prima approssimazione, al valore
negativo dell’autovalore orbitale corrispondente:
Ii ≈ −εi ([Link])
Questo risultato implica che gli autovalori εi dell’operatore di Fock, nella sua forma canonica,
rappresentano (a meno di un segno) le energie di ionizzazione degli elettroni appartenenti ai rispet-
tivi spin-orbitali. In termini fisici, ciò significa che la variazione di energia associata alla rimozione di
un elettrone da un orbitale ϕi è pari all’opposto del corrispondente autovalore di Fock. Prendiamo
l’energia per N elettroni:
N N
X 1X
EN = ⟨Ψi |ĥ|Ψi ⟩ + ⟨Ψi Ψj ||Ψi Ψj ⟩ + Enucl ([Link])
i
2 i,j
l’energia per N − 1 elettroni ottenuto rimuovendo un elettrone dallo spin-orbitale a invece è:
N N
(a)
X 1 X
EN −1 = ⟨Ψi |ĥ|Ψi ⟩ + ⟨Ψi Ψj ||Ψi Ψj ⟩ + Enucl ([Link])
2
i̸=a i̸=a,j
che è proprio l’energia dello spin-orbitale associato, dimostrando quindi il teorema di Koop-
mans. Nella realtà però una volta tolto l’elettrone il sistema rilassa e gli spin orbitali
che descrivono il sistema con N − 1 elettroni non potranno essere uguali a quelli che
descrivevano il sistema con N elettroni; tale interpretazione quindi, sebbene estrema-
mente utile, non è rigorosamente esatta. L’operatore di Fock, infatti, non rappresenta un
osservabile fisico: esso non è lineare e dipende dalla stessa soluzione che si cerca. Pertanto, gli
autovalori εi non corrispondono a quantità fisiche direttamente misurabili, ma solo ad
approssimazioni di energie di ionizzazione, valide nell’ambito del modello di campo
medio. Nonostante questa limitazione, il teorema di Koopmans fornisce un potente collegamento
tra teoria e osservazione sperimentale: consente infatti di interpretare gli autovalori di Fock come
differenze di energia tra il sistema a N elettroni e quello a N − 1 elettroni, rendendo la teoria di
102
Hartree–Fock un efficace strumento di analisi qualitativa e quantitativa della struttura elettronica
dei sistemi atomici e molecolari. Se volessimo trattare in modo più preciso la ionizzazione dovremmo
usare tuttavia la teoria dello scattering quantistico.
Entriamo ora più nello specifico della teoria di Hartree-Fock. Vi sono due versioni del metodo:
Restricted (RHF) e Unrestricted (UHF). Quando il set di orbitali fornito è lo stesso sia per gli α
che per i β siamo in RHF. Nel UHF invece non abbiamo questa limitazione e possiamo usare due
set diversi di orbitali, uno per gli spin α e uno per i β. Useremo RHF quando abbiamo un numero
pari di elettroni, in questo modo per N orbitali abbiamo che N/2 sono occupati e gli altri N/2
sono invece detti orbitali virtuali (se ci sono degli elettroni spaiati il metodo è open shell-RHF);
ogni orbitale in RHF ha quindi un occupazione doppia. In UHF invece per ogni orbitale
dobbiamo scegliere uno spin da assegnargli. Per via di tutte queste considerazioni tendenzialmente
i calcoli con RHF danno energie maggiori rispetto a quelli fatti in UHF.
Iniziamo trattando il Restricted Hartee-Fock poiché è quello che useremo più spesso. Prendiamo
come funzioni di lavoro degli orbitali di tipo atomico χk , questo poiché si è visto che in generale si
lavora meglio (computazionalmente) con un set di funzioni distribuite sui nuclei (anche se il set non è
completo), rispetto ad un set completo non distribuito (come l’insieme delle funzioni dell’oscillatore
armonico per esempio). Infatti un tale set consente di ridurre il numero di funzioni a un valore
gestibile, limitando l’onere computazionale. Noi comunque cerchiamo degli orbitali molecolari della
forma:
m
X
Φj (⃗r) = χk (⃗r)Tkj ([Link])
k
dove i Tkj sono i parametri variazionali. Sappiamo che l’operatore coulombiano Jˆ ha la forma:
N Z
X (Ψ∗j Ψj )(x′ )
Jˆ = dx′ ([Link])
j
|⃗r − ⃗r′ |
ma nello schema RHF abbiamo metà degli elettroni con spin α e l’altra metà con spin β la somma
si può quindi separare in:
N/2 Z
X (χ∗j χj )(⃗r′ )
ˆ
J =2 d⃗r′ ([Link])
j
|⃗r − ⃗r′ |
Ovvero per ogni elettrone ho due termini riferiti uno ad α e uno a β; per via di ciò e
siccome dx′ = d⃗r′ dσ (usando le proprietà delle funzioni di spin) abbiamo potuto sostiuire dx′ con
d⃗r. In questo caso possiamo scrivere:
JˆRHF = 2Jˆorb ([Link])
dove Jˆorb è il contributo sugli orbitali. Anche la densità sarà esprimibile come:
N/2
X
ρRHF (⃗r) = 2 |χi |2 ([Link])
j
Per quanto riguarda l’operatore di scambio invece esso può essere scritto come:
103
dove i termini tra parentesi graffa sono detti proiettori. Per dimostrare questa scrittura conside-
riamo il Kernel dell’operatore:
N
X Ψj (x)Ψj (x′ )∗
K(x, x′ ) = ([Link])
j
|⃗r − ⃗r′ |
il termine tra le graffe garantisce che ”sopravviva” solo la parte di spin equivalente a quella della
funzione di lavoro e quindi il numero di interazioni di scambio è dimezzato rispetto a quelle coulom-
biane. Comunque vediamo con un esempio generico:
Z
K̂ϕ(⃗r)α(σ) = (K̂orb ϕ)(⃗r) dσ ′ (αα∗ + ββ ∗ )α(σ ′ ) = (K̂orb ϕ)(⃗r)α(σ) ([Link])
ovvero in RHF possiamo esprimere il problema energetico solo in funzione della parte spaziale:
ed ogni orbitale sarà poi occupato da due elettroni con spin opposto. Se usiamo come funzioni di
lavoro degli orbitali atomici ci troviamo all’interno dell’approssimazione LCAO, e gli orbitali trovati
saranno detti orbitali molecolari, sempre esprimibili nella forma:
m
X
Φj (⃗r) = χk (⃗r)Tkj ([Link])
k
104
che è esprimibile in forma matriciale come:
FT = STε ([Link])
questo è un sistema di equazioni differenziali noto come equazioni di Roothaan. T è una matrice
(m×m) e ogni sua colonna rappresenta i coefficienti di un singolo orbitale molecolare e per ognugno
di questi possiamo scrivere un equazione del tipo:
Questa equazione rappresenta dunque un sistema di m equazioni lineari omogenee, le cui incog-
nite sono i coefficienti Tkj :
(F − εi S) Ti = 0
T1j
T2j ([Link])
(F − εi S)
. =0
Tkj
ed avrò soluzioni non banali se e solo se:
det(F − εi S) = 0 ([Link])
L’equazione di Roothaan tuttavia non è realmente un sistema lineare, poiché la matrice di Fock F
dipende dagli stessi coefficienti T attraverso i termini coulombiani e di scambio (J e K). In altre
parole:
F = F(C). ([Link])
Ne consegue che non è possibile risolvere il sistema in un solo passaggio. È necessario
adottare una procedura iterativa:
1. Scelta di un guess iniziale: Si parte da una stima iniziale T(0) per la matrice dei coefficienti
(detta initial guess), che può essere ottenuta con diversi criteri — ad esempio, combinazioni
di orbitali atomici o orbitali di legame preliminari.
2. Costruzione della matrice di Fock: Dalla matrice T(0) si calcola la corrispondente matrice
di Fock F(0) .
105
5. Iterazione: Si ripetono i passi precedenti fino a quando la soluzione non raggiunge l’autoconsistenza,
ovvero fino a quando le matrici di Fock di due iterazioni successive differiscono di meno di
una certa soglia:
F(n−1) T(n) = S T(n) ε(n) ([Link])
∥T(n−1) − T(n) ∥2 ≤ δ ([Link])
con δ tipicamente dell’ordine di 10−8 –10−10 .
Quando la convergenza è raggiunta, le matrici T(n) e F(n) rappresentano rispettiva-
mente i coefficienti orbitali e l’operatore di Fock auto–consistente, e gli autovalori εi
costituiscono le energie orbitali finali del sistema.
Pertanto, si selezionano i primi N/2 orbitali (nel caso restricted ) di energia più bassa:
bensı̀ nella diagonalizzazione diretta della matrice di Fock F. Le basi di funzioni atom-
iche {χk } impiegata per costruire gli orbitali molecolari infatti non sono, in generale, ortonormali.
Funzioni di tipo s, p, d, poste su centri atomici differenti, non risultano ortogonali, in particolare
quando gli atomi sono vicini. In tal caso la matrice di sovrapposizione S che ha come elementi:
106
non sarà quindi diagonale.
Per lavorare in una base ortonormale, si effettua una trasformazione lineare delle funzioni di base:
tale che
⟨χ̄k |χ̄l ⟩ = δkl ([Link])
Ovvero diagonalizza S. Una delle ortogonalizzazioni più comuni e simmetriche è l’ortogonalizzazione
di Löwdin, basata sulla radice quadrata della matrice di sovrapposizione:
Questa trasformazione permette di passare da una base non ortogonale a una base ortonormale,
semplificando notevolmente la forma matriciale dell’equazione di Hartree-Fock. Consideriamo ora
l’equazione di Hartree-Fock nella sua forma matriciale:
FT = STε ([Link])
107
l’equazione diventa:
F̄T̄ = T̄ε ([Link])
Si tratta ora di un’equazione di autovalori nella forma canonica, cioè di una diagonalizzazione di
una matrice hermitiana, in cui la matrice T è unitaria.
dove T̄† = T̄−1 per via dell’unitarietà. In pratica, si utilizza un algoritmo numerico (ad esempio
Jacobi o Householder) che restituisce la matrice unitaria T̄ e il vettore degli autovalori ε.
T = S−1/2 T̄ ([Link])
In questo modo, gli orbitali molecolari risultano espressi come combinazioni lineari di
orbitali atomici ortonormali.
La matrice T cosı̀ definita conserva il carattere atomico delle funzioni di base, poiché la
trasformazione S−1/2 è simmetrica e tratta tutti gli orbitali in modo bilanciato. Dunque, possiamo
scrivere gli MO come:
m
X m
X m
X m
X
ϕj = χ̄k T̄kj = χk Tkj = χk (S −1/2 )kl T̄lj ([Link])
k=1 k=1 k=1 l=1
Per illustrare la procedura, consideriamo il caso più semplice della molecola di idrogeno H2 .
Utilizziamo come funzioni di base due orbitali atomici 1s, uno centrato su ciascun nucleo:
χS e χS ′ ([Link])
ψ = A χS + B χS ′ ([Link])
dove
A = ⟨χS |F̂ |χS ⟩ = ⟨χS ′ |F̂ |χS ′ ⟩, B = ⟨χS |F̂ |χS ′ ⟩ ([Link])
108
La matrice di sovrapposizione è anch’essa 2 × 2:
1 s
S= , ([Link])
s 1
dove s = ⟨χS |χS ′ ⟩ rappresenta l’integrale di sovrapposizione, positivo e dipendente dalla distanza
internucleare. All’aumentare della distanza, s → 0, ma a distanza di legame tipica s assume valori
dell’ordine di 0.1–0.2.
(F − εS)T = 0 ([Link])
ovvero:
A − ε B − εs T1
=0 ([Link])
B − εs A − ε T2
Affinché esistano soluzioni non banali, il determinante deve annullarsi:
Espandendo e semplificando:
ε− < ε+ ([Link])
Pertanto, la soluzione con il segno “+” (combinazione costruttiva) corrisponde allo stato legante,
caratterizzato da una densità elettronica aumentata nella regione interatomica e da un’energia più
bassa. La soluzione “−” (combinazione distruttiva) descrive invece lo stato antilegante, a energia
maggiore.
Analizziamo quindi nello specifico il segno del termine B. Per definizione:
1
B = ⟨χS | − ∇2 + VSCF |χS ′ ⟩ ([Link])
2
dove VSCF (r) è il potenziale self-consistent, contenente il contributo nucleare e quello elettronico.
109
Nella regione di sovrapposizione, il prodotto χS χS ′ è positivo, mentre il potenziale VSCF è negativo,
dominato dal termine attrattivo nucleare. Pertanto, l’integrale complessivo risulta negativo:
B < 0. ([Link])
ψ± = N± (χS ± χS ′ ) ([Link])
dove N± sono costanti di normalizzazione. L’orbitale ψ− (con segno + nella combinazione) è legante
e simmetrico rispetto all’inversione delle coordinate, mentre ψ+ è antilegante e antisimmetrico. Nel
caso di H2 , gli elettroni occupano l’orbitale di energia più bassa, ε− , lasciando vuoto quello antile-
gante.
In questo semplice caso non è necessario alcun procedimento iterativo: la simmetria del sistema
determina completamente la forma delle soluzioni, e la convergenza è immediata.
2.6.2 Scelta delle funzioni di base: orbitali atomici di tipo Slater e Gaussiani
È il momento di introdurre le basi di funzioni atomiche, ossia le funzioni base localizzate sugli
atomi, che sono quelle utilizzate nel metodo LCAO (Linear Combination of Atomic Orbitals) a tutti
i livelli di teoria: non solo nell’approccio di tipo Hartree–Fock, ma anche nei metodi più avanzati,
come il metodo a campo autoconsistente (SCF) e oltre.
Queste basi, chiamate basi di orbitali atomici o atomic orbitals (AO basis sets), sono distribuite sui
centri atomici. Questa scelta è la più conveniente dal punto di vista computazionale e fisico, poiché,
anche quando vogliamo descrivere un materiale o una molecola estesa, il comportamento elettronico
locale ricorda fortemente quello degli orbitali atomici dei singoli atomi costituenti. Pertanto, una
base formata da orbitali localizzati su ciascun nucleo è più adeguata rispetto a una base completa
puramente monocentrica. Nel metodo di Hartree–Fock,per esempio, la scelta delle funzioni atom-
iche di base è un aspetto cruciale che influenza in modo determinante la qualità e l’efficienza del
calcolo. A seconda del livello di approssimazione e del tipo di sistema studiato, è possibile utilizzare
110
diverse basi, che vengono classificate in base al numero di gusci atomici rappresentati. Una base
minima è costituita dal numero minimo di funzioni necessario per descrivere ciascun
atomo, cioè un guscio per ogni livello energetico occupato. Ad esempio, per l’idrogeno o
l’elio, una base minima comprende solo un orbitale 1s. Una base doppia include due gusci per ogni
atomo (ad esempio 1s e 2s, 2p), mentre una base tripla aggiunge anche il terzo guscio (3s, 3p, 3d),
e cosı̀ via.
Il numero di funzioni atomiche cresce rapidamente con il numero dei gusci inclusi. Infatti, ogni
livello quantico n contiene n2 funzioni orbitali, e il numero totale di funzioni fino al livello N cresce
in modo cubico:
N (N + 1)(2N + 1)
12 + 22 + 32 + · · · + N 2 = ∼ O(N 3 ). ([Link])
6
Di conseguenza, l’aumento del numero di gusci porta a una crescita cubica del numero
di funzioni atomiche e, soprattutto, del numero di integrali bielettronici da calcolare.
In passato questo rappresentava un forte limite computazionale; oggi, grazie alle risorse di calcolo
moderne, il problema è molto più gestibile, ma resta comunque un aspetto critico.
Questo comportamento esponenziale è cruciale, perché riproduce fedelmente sia la coda a grandi
distanze sia la forma della funzione vicino al nucleo. Consideriamo ora l’orbitale fondamentale 1s,
corrispondente ai numeri quantici:
n = 1, l = 0.
In questo caso la funzione radiale si riduce a:
111
Analizziamo ora il comportamento della funzione d’onda per r → 0 e per r → ∞.
• Per r → 0: I termini dominanti sono quelli contenenti fattori in 1/r provenienti dall’operatore
di Laplace e dal potenziale coulombiano, poiché sono i più divergenti. Espandendo la funzione
per piccoli valori di r, otteniamo:
1
e−Zr ≃ 1 − Zr + (Zr)2 + · · · .
2
Sostituendo nell’Hamiltoniano:
1 d2
2 d Z
Ĥψ(r) = − + (1 − Zr + · · · ) − (1 − Zr + · · · ) .
2 dr2 r dr r
Si nota che il termine in 1/r proveniente dal potenziale bilancia esattamente la derivata radiale
al primo ordine, garantendo una cuspide all’origine (r = 0), cioè una discontinuità nella
derivata della funzione.
• Per r → ∞: La funzione decresce esponenzialmente come e−Zr , coerentemente con il com-
portamento atteso per un elettrone legato in un potenziale coulombiano.
Questo andamento di cuspide all’origine e il decadimento esponenziale a grandi distanze, mostra
che le funzioni d’onda idrogenoidi (e in particolare le parti esponenziali di quest’ultime) forniscono
un modello fisico accurato degli orbitali atomici reali. Per questo motivo, nelle basi atomiche
utilizzate nei metodi di struttura elettronica, si preferiscono funzioni di tipo esponen-
ziale, distribuite sui centri atomici, poiché descrivono in modo realistico sia la regione
vicina al nucleo che la coda elettronica.
Nel metodo di Hartree–Fock, l’interazione di scambio garantisce che due elettroni con lo stesso
spin non possano trovarsi nello stesso punto dello spazio, riducendo cosı̀ la probabilità di collisione.
Tuttavia, gli elettroni con spin opposto non sono soggetti a tale vincolo: possono trovarsi nella
stessa regione spaziale e quindi “collidere”. Questo effetto non gestito dal solo termine di scambio
porta a un incremento dell’energia del sistema, che viene successivamente corretto nei metodi più
accurati introducendo la cosiddetta energia di correlazione elettronica. Di essa parleremo più
approfonditamente in seguito.
STO e GTO
Partendo dalle considerazioni precedenti sono state sviluppate due categorie di funzioni atomiche
di base principali:
1. le funzioni di tipo Slater (Slater Type Orbitals, STO);
2. le funzioni di tipo Gaussiano (Gaussian Type Orbitals, GTO).
112
Orbitali di tipo Slater
Le funzioni di tipo Slater sono concepite per imitare il comportamento delle soluzioni esatte
dell’atomo di idrogeno o, più in generale, di un sistema idrogenoide (un nucleo di carica Z con
un solo elettrone). Per un atomo idrogenoide, le funzioni d’onda in coordinate polari, hanno la
forma:
zr
ψnℓm (r, θ, ϕ) = rℓ Rn,l (r) e− n Yℓm (θ, ϕ) ([Link])
Le funzioni di tipo Slater adottano una forma semplificata ma fisicamente realistica:
dove l’esponenziale e−ζr riproduce correttamente il decadimento esponenziale delle soluzioni idrogenoidi.
Queste funzioni descrivono accuratamente sia il comportamento vicino al nucleo sia
la coda esponenziale a grande distanza. Infatti, come abbiamo già visto, un aspetto fonda-
mentale delle funzioni esatte dell’atomo di idrogeno è la cosiddetta condizione di cuspide. Nel
potenziale coulombiano V (r) = −Z/r, quando r → 0 il potenziale tende a −∞, generando una
singolarità nella funzione d’onda. Affinché l’energia resti finita, il comportamento della funzione
deve compensare tale divergenza.
cioè la pendenza della funzione all’origine deve essere proporzionale alla carica nucleare.
Per una dimostrazione di Eq. [Link] cercare la condizione di cuspide di Kato.
Questa è la condizione di cuspide, che gli orbitali di tipo Slater rispettano naturalmente, garan-
tendo una rappresentazione fisicamente corretta nei pressi del nucleo. Per via della loro formulazione
aderente con la fisica le funzioni di Slater rispettano di default anche il comportamento asintotico
delle funzioni; a grande distanza (r → ∞), il decadimento esponenziale e−ζr delle funzioni di tipo
Slater riproduce correttamente la diminuzione della densità elettronica, in accordo con il potenziale
di ionizzazione del sistema. In sintesi, gli orbitali di tipo Slater sono fisicamente accurati
sia vicino al nucleo sia nella regione asintotica.
113
2. il loro decadimento a grande distanza è troppo rapido rispetto al comportamento fisico cor-
retto.
Dal punto di vista fisico, quindi, le GTO sono meno adatte a descrivere accuratamente gli or-
bitali atomici. Tuttavia, esse offrono un enorme vantaggio computazionale: tutti gli integrali
molecolari espressi in termini di funzioni gaussiane possono essere calcolati in forma
analitica. Questo significa che gli integrali monoelettronici e bielettronici possono essere trattati
con precisione elevata e in modo estremamente efficiente al calcolatore. Le funzioni di tipo Slater,
al contrario, richiedono l’integrazione numerica, molto più costosa e complessa. Per questo motivo,
nonostante la minore aderenza fisica, le funzioni gaussiane sono ampiamente utiliz-
zate nella chimica quantistica computazionale moderna. Sono infatti le basi più usate nei
principali codici di chimica quantistica (come Gaussian, ORCA, Q-Chem, ecc.).
Tuttavia, poiché una funzione STO descrive meglio la fisica del problema si cerca di scriverle
come combinazione lineari di GTO:
N
X 2
e−ζr ≈ ci e−αi r ,
i=1
dove i coefficienti ci e gli esponenti αi vengono scelti in modo da riprodurre al meglio la forma dello
STO. In generale, serve un numero N > 1 di gaussiane per ottenere una buona approssimazione
(ad esempio, 3, 6 o 9 gaussiane in una “contracted Gaussian”).
Esiste inoltre una rappresentazione integrale esatta che collega una funzione di tipo Slater a una
combinazione continua di gaussiane. In particolare, si può scrivere:
Z +∞
1 ζ ζ2 2
e−ζr = √ √ e− 4α e−αr d(ln α).
π −∞ 2 α
Questa relazione mostra che una funzione di tipo Slater può essere vista come una su-
perposizione continua di gaussiane con pesi opportuni. In pratica, utilizzando un numero
finito di gaussiane si ottiene un’approssimazione discreta molto accurata, base del principio su cui
si fondano i moderni set di basi gaussiane (come STO-3G, 6-31G, cc-pVTZ, ecc.). Questa relazione
114
mostra chiaramente il collegamento tra le due forme funzionali: da un lato il termine e−ζr , tipico
2
delle funzioni di tipo Slater, e dall’altro e−αr , proprio delle gaussiane. Sono esattamente le due
tipologie di funzioni che stiamo confrontando per definire la forma più adatta agli orbitali di base
atomici.
Se analizziamo quindi la funzione:
ζ ζ2
f (α) = √ e− 4α ,
2 α
In questa espressione, i valori di αk sono scelti in progressione geometrica. Definendo una costante
iniziale α0 e un fattore di scala β > 1, si ha:
αk = α0 β k
∆(ln α) = ln αk+1 − ln αk = ln β
α0 ≈ 0.01, β ≈ 3, M ≈ 10,
In altre parole, una funzione di tipo Slater può essere approssimata in modo eccellente
mediante la combinazione di circa 10 funzioni gaussiane, ovvero:
M
X 2
e−ζr ≈ ck e−αk r ,
k=1
115
temperate, o even-tempered Gaussians. Questa formulazione fu introdotta per la prima volta da
Rudenberg, che mostrò come fosse possibile ottenere buone approssimazioni delle funzioni di Slater
mediante un numero limitato di gaussiane con esponenti distribuiti geometricamente.
Il primo a proporre esplicitamente l’uso di funzioni gaussiane nei calcoli molecolari fu però Roy
McWeeny, nel 1957. A quell’epoca, tuttavia, la potenza di calcolo disponibile non era sufficiente
a rendere pratiche tali applicazioni: l’idea era teoricamente valida, ma prematura. Qualche anno
più tardi, John A. Pople riprese e sviluppò il concetto, implementandolo nel celebre programma
Gaussian. La prima versione del codice Gaussian risale al 1970 e rappresentò una svolta epocale
nella chimica computazionale, grazie all’efficienza del calcolo degli integrali molecolari basati su
gaussiane. Le basi sviluppate da Pople (come le famose 6-31G, 6-311G**, ecc.) si fondano proprio
su questa idea di combinare più gaussiane per approssimare un singolo orbitale di tipo Slater. Il
successo del metodo fu tale che nel 1998 Pople ricevette il Premio Nobel per la Chimica, con-
diviso con Walter Kohn — quest’ultimo premiato per lo sviluppo della teoria del funzionale della
densità (DFT). Nel laboratorio vedremo come queste funzioni gaussiane siano effettivamente im-
plementate nei software di calcolo e come vengano utilizzate per descrivere realisticamente orbitali
atomici e molecolari.
Indipendentemente dal livello di complessità della funzione d’onda, l’energia ha sempre una forma
di questo tipo:
X X 1
E = EN N + ⟨χi |ĥ|χj ⟩Pij [Ψ] + ⟨χi χj | |χk χl ⟩ij Γijkl [Ψ]
ij ij
r12
Dove P e Γ dipendono dal tipo di funzione d’onda. Il metodo che usiamo determina i loro valori.
In formule generali, l’operatore monoelettronico agisce su un singolo orbitale χi come:
NX
nuclei
1 ZA
h(i) = − ∇2i − ,
2 |ri − RA |
A=1
116
Gli integrali monoelettronici sono quindi:
Z NX
nuclei
1 ZA
hij = χ∗i (r) (− ∇2i − ) χj (r) dr
2 |ri − RA |
A=1
Gli integrali monoelettronici sono calcolati una volta sola per tutta la base o, in alternativa, pos-
sono essere letti da disco se già disponibili. Oggi, grazie alla velocità di lettura dei dischi e alla
disponibilità di memoria, spesso si preferisce caricarli direttamente in memoria per velocizzare i
calcoli successivi. Essendo una matrice simmetrica (hij = hji ), il numero di integrali distinti da
2
calcolare è m2 , dove m è il numero di orbitali atomici nella base.
Gli integrali bielettronici invece coinvolgono due elettroni distinti e sono espressi come integrali
a sei dimensioni: ZZ ∗
χi (r1 )χ∗j (r2 )χk (r1 )χl (r2 )
⟨ij|kl⟩ = dr1 dr2
|r1 − r2 |
In questi integrali, r1 e r2 rappresentano le coordinate dei due elettroni, mentre le funzioni χi , χj , χk , χl
sono gli orbitali atomici (o contratti) scelti come base. Questi integrali contengono tutti i numeri
quantici necessari per descrivere correttamente gli elettroni, ma in forma generica possono essere
indicati con indici i, j, k, l, senza specificare la loro origine.
Gli integrali bielettronici rappresentano la repulsione coulombiana tra due distribuzioni di carica.
Questi integrali coinvolgono quattro indici e, per una base di m orbitali, il numero totale di inte-
grali bielettronici è dell’ordine di m4 . Tuttavia, sfruttando le simmetrie permutazionali, possiamo
ridurre il numero effettivo di calcoli:
• scambio di i ↔ k;
• scambio di j ↔ l;
• scambio degli elettroni 1 e 2.
Considerando queste simmetrie, il numero effettivo di integrali distinti diventa circa m4 /8. Come
osservato nei calcoli di tipo Hartree-Fock, il costo computazionale principale non deriva dalla
diagonalizzazione della matrice Fock, che scala come m3 per algoritmi tipo Jacobi, ma
dal calcolo degli integrali bielettronici, che rappresenta il vero collo di bottiglia. Negli
anni ’60 e ’70, questo problema spinse allo sviluppo dei metodi semi-empirici, che permettono di
bypassare il calcolo esatto degli integrali bielettronici utilizzando parametri calibrati su molecole
di riferimento. Questi metodi, pur non essendo completamente ab initio, offrono un compromesso
tra precisione e costi computazionali Nei metodi ab initio, tutti gli integrali elettronici (mono- e
bielettronici) sono calcolati esattamente, con un errore trascurabile. Al contrario, nei metodi semi-
empirici (es. CNDO, INDO), alcuni integrali sono sostituiti da valori parametrizzati su molecole
note, introducendo approssimazioni controllate.
Per esempio, il metodo Kohn-Sham della DFT viene spesso considerato ab initio nella pratica,
anche se il funzionale dell’energia contiene parametri e non è rigorosamente derivato dalle equazioni
di Schrödinger. In senso stretto, quindi, sarebbe semi-empirico, ma viene comunque trattato come
ab initio perché non utilizza calcoli di integrali parametrizzati come CNDO.
117
In ogni caso gli integrali bielettronici sono espressi come integrali a sei dimensioni: richiedono
quindi calcoli complessi, potenzialmente analitici solo per alcune scelte di orbitali (come gaussiane
contratte). Gli integrali monoelettronici, al contrario, sono molto più semplici da calcolare e non
rappresentano un problema in termini computazionali. In particolare con gli Slater Type Orbitals
(STO), gli integrali sono analitici solo per sistemi a centro singolo, come gli atomi isolati. In-
fatti, per calcoli accurati sugli atomi, si utilizzano le STO e non le gaussiane, perché con le STO
servono meno funzioni atomiche per ottenere la stessa accuratezza. Al contrario, le gaussiane libere
richiederebbero un numero maggiore di funzioni per descrivere correttamente gli orbitali atomici.
Ad esempio, nelle cosiddette basi di orbitali naturali, il numero di gaussiane generiche necessarie
può essere molto elevato (fino a decine di gaussiane per singolo orbitale), rendendo i calcoli pesanti.
Tuttavia per molecole con più centri (ad esempio l’acqua), gli integrali bielettronici con STO non
possono più essere calcolati in modo completamente analitico. In questi casi si utilizzano griglie o
quadrature numeriche. Le formule di quadratura si basano su polinomi speciali e scelgono punti
di valutazione opportuni in funzione della distanza o della forma della funzione, per approssimare
l’integrale. Tuttavia, se la dimensione della griglia è grande (ad esempio 10–100 punti per ogni coor-
dinata), in uno spazio a sei dimensioni il numero totale di punti da calcolare diventa enorme (1006 ),
rendendo il metodo troppo costoso. Un’alternativa in questi casi è rappresentata dai metodi di
Monte Carlo: gli integrali vengono calcolati in modo stocastico, campionando punti casuali con
una distribuzione adeguata all’interno dell’integrale. Questo metodo è particolarmente competitivo
per integrali a dimensione superiore a quattro. Con questo approccio, non importa se le fun-
zioni siano STO o gaussiane: il calcolo procede nello stesso modo, utilizzando il campionamento
stocastico.
Gli integrali con le gaussiane invece sono sempre analitici. Per mostrarlo consideriamo
una gaussiana tipo s centrata in un punto A:
2
ϕs (r) = e−a|r−A|
dove a è l’esponente della gaussiana. La simmetria di questa funzione è sferica attorno al centro A.
Dalla gaussiana tipo s è possibile generare gaussiane tipo px , py , pz derivando rispetto alle coor-
dinate del centro:
∂ −a|r−A|2 2
ϕPx (r) = e = −2a(x − xA )e−a|r−A|
∂xA
Allo stesso modo, derivando più volte rispetto a x, y, z si ottengono funzioni tipo d, f e cosı̀ via.
Tutte queste operazioni sono analitiche e, di conseguenza, basta sapere come calcolare gli integrali
per le gaussiane tipo S per ottenere in maniera analitica anche gli integrali su tutte le funzioni
generate (P, D, F, ecc.).
Gli integrali di gaussiane sono analitici anche su più centri, infatti:
2 2 2 2
e−a|r1 −A| e−b|r2 −B| e−c|r1 −C| e−d|r2 −D|
Z
⟨ab|cd⟩ = dr1 dr2 ,
|r1 − r2 |
dove a, b, c, d sono gli esponenti delle gaussiane e A, B, C, D sono i centri dei rispettivi orbitali.
Anche dopo aver integrato su r1 e r2 , l’integrale dipende ancora dai centri A, B, C, D, che non
sono variabili di integrazione ma semplicemente le posizioni delle gaussiane. Se vogliamo calcolare,
ad esempio, un integrale contenente una funzione di tipo px centrata in A, possiamo ottenere il
118
risultato derivando l’integrale s rispetto alla coordinata xA . In questo modo, calcolando l’integrale
sulle gaussiane s, possiamo ricavare tutti gli integrali necessari per funzioni atomiche di qualsiasi
tipo. Un’altra proprietà fondamentale delle gaussiane è la riduzione dei centri: un integrale
bielettronico su quattro centri può essere ridotto a un integrale su due centri, che è calcolabile in
forma analitica. Consideriamo il prodotto di due gaussiane centrate in A e B:
2 2
e−a|r−A| e−b|r−B| .
−a|r − A|2 − b|r − B|2 = −(a + b)|r|2 + 2r · (aA + bB) − (a|A|2 + b|B|2 ).
Riorganizzando i termini, possiamo riscrivere il prodotto come una singola gaussiana centrata in
un punto intermedio P definito come media pesata dei centri originali:
aA + bB
P= ,
a+b
con esponente totale a + b. La parte indipendente da r viene raccolta in un fattore costante
2 2 2
γ = e−(a|A| +b|B| −(a+b)|P| ) .
Il prodotto di due gaussiane quindi genera una nuova gaussiana centrata in un punto interme-
dio, moltiplicata per un fattore numerico proveniente dai termini rimanenti. Questo permette di
definire dei centroidi su cui calcolare gli integrali bielettronici in modo analitico. Un integrale
bielettronico generico sarà esprimibile come:
1
⟨ac| |bd⟩ = ΛSac Sbd F (t)
r12
dove Λab,cd dipende dagli esponenti delle gaussiane a + b e c + d, normalizzato da fattori π −1/2
1
(a + b)(c + d) 2
Λ=2
π(a + b + c + d)
I centroidi SB sono integrali di sovrapposizione delle gaussiane, mentre il parametro t dipende dal
rapporto degli esponenti e dalla distanza tra i centroidi.
119
Una volta calcolati gli integrali sulle gaussiane tipo S, tutti gli altri integrali necessari si otten-
gono tramite derivate rispetto alle coordinate dei centri, rendendo il procedimento completamente
analitico. Questo è il motivo per cui le gaussiane sono cosı̀ utili: gli integrali si calcolano con
formule note, veloci da eseguire con i moderni processori.
Per dare un’idea delle prestazioni, un processore moderno capace di 1010 moltiplicazioni al sec-
ondo può calcolare 100 moltiplicazioni in circa un milionesimo di secondo. Negli anni ’60 e ’70,
invece, i computer riuscivano a fare appena 6 × 106 moltiplicazioni al secondo. All’epoca, molti
calcoli venivano fatti praticamente a mano, e i tempi erano molto più lunghi. Oggi, con la memo-
ria disponibile, possiamo conservare gli integrali bielettronici su disco o in memoria, o calcolarli al
bisogno. La scelta dipende dalla dimensione del calcolo: leggere dal disco o calcolare al momento può
essere equivalente in termini di tempo. Un tempo, scrivere integrali bielettronici su disco richiedeva
strategie particolari, combinando quattro indici in un unico numero per risparmiare spazio. Per
riprodurre gli orbitali di tipo Slater (STO) con gaussiane (GTO), si utilizzano tipicamente molte
gaussiane per ogni STO. Ad esempio, per un 1s STO si possono usare 10 gaussiane, per un 2s
altre 10, coprendo un intervallo molto ampio di esponenti, da 0.01 fino a 100 000, in progressione
geometrica. Gli esponenti maggiori producono gaussiane molto acute, necessarie per riprodurre la
cuspide vicino al nucleo.
Tuttavia, usare molte gaussiane per ogni STO è estremamente costoso. Per calcoli più pratici
si ricorre ai già menzionati schemi di contrazione, dove più gaussiane vengono combinate in
blocchi per definire un singolo orbitale atomico. Esempi di queste basi contratte sono quelle di
Dunning, TZVP, NZ o di Pople, come la 6-3-1-1G. Suddividendo le gaussiane in blocchi, possi-
amo creare orbitali atomici compatti e gestibili, riducendo il costo computazionale mantenendo una
buona qualità del calcolo.
Con questi schemi è possibile ottimizzare il numero di funzioni atomiche e ridurre il costo com-
putazionale; non si utilizzano quindi tutte le 10 gaussiane per generare un singolo STO, ma si
preferisce suddividerle in blocchi generati secondo lo schema di contrazione in uso, ad esempio
6-3-1. In questo schema infatti:
• Le prime 6 gaussiane vengono combinate per generare un orbitale atomico più compatto.
6
(1)
X
χ1 = ck gk
k=1
120
• Le successive 3 gaussiane generano un orbitale meno compatto.
9
(2)
X
χ2 = ck gk
k=7
La combinazione finale degli orbitali atomici derivanti da queste contrazioni è determinata dal
calcolo molecolare stesso, adattandosi al contesto chimico dell’atomo. Questo approccio permette
di ottenere una base flessibile, che non è vincolata a un singolo atomo isolato, ma può essere
utilizzata anche in molecole più complesse, come il metano.
Esempio: Carbonio
Per mostrare il metodo facciamo un esempio con il carbonio; le 6 gaussiane più compatte vengono
utilizzate per descrivere il guscio interno (1s), mentre le 3 gaussiane successive e la gaussiana libera
contribuiscono alla formazione degli orbitali valenza. Cosı̀:
• Gli orbitali 1s derivano da una combinazione lineare delle prime 3 gaussiane contratte.
• Gli orbitali 2s derivano da una combinazione lineare delle stesse 3 gaussiane, ma con coeffici-
enti differenti, appropriati per descrivere l’orbitale 2s.
• Gli orbitali di valenza 2p (px , py , pz ) vengono generati utilizzando solo le 4 gaussiane più
diffuse del blocco finale (ad esempio α7 , α8 , α9 , α10 ).
In questo modo, la contrazione dei blocchi di gaussiane permette di ridurre il numero di fun-
zioni necessarie per descrivere sia gli orbitali interni sia quelli di valenza, garantendo una maggiore
flessibilità per applicazioni molecolari. Gli orbitali atomici ottenuti in questo modo riman-
gono flessibili: possono descrivere sia l’atomo isolato sia lo stesso atomo in un contesto
molecolare, mantenendo una buona rappresentazione della densità elettronica. In gen-
erale nel caso delle valenze p, come nel carbonio, si usano solo le gaussiane più diffuse del secondo
blocco, evitando di impiegare le 6 gaussiane più compatte che vengono invece assegnate agli orbitali
1s. In particolare, per un px atomico, l’orbitale sarà una combinazione lineare delle 4 gaussiane fi-
nali con i coefficienti adeguati, eventualmente moltiplicata per x per ottenere la simmetria corretta.
La stessa logica vale per py e pz .
Due tra i principali contributori alla costruzione di basi gaussiane contratte sono Pople e Dunning,
i quali hanno sviluppato approcci differenti:
121
Nelle basi di Pople o Dunning, se si stampano gli esponenti dal programma Gaussian, non si os-
serva una progressione geometrica perfetta come nelle basi di tipo even-tempered. Tuttavia, gli
esponenti più grandi e più piccoli seguono comunque una progressione quasi geometrica, con valori
che possono variare da 10−2 fino a 105 anche per atomi del primo periodo. Gli esponenti più grandi
sono associati ai gusci interni, mentre quelli più piccoli rappresentano le funzioni diffuse necessarie
a descrivere le regioni esterne della densità elettronica.
Vediamo quindi il metodo di Dunning; esso può essere visualizzato chiaramente attraverso un es-
perimento computazionale. Supponiamo di eseguire un calcolo su un atomo di carbonio utilizzando
10 gaussiane di tipo s e, per l’idrogeno, 4 gaussiane di tipo s. Il calcolo completo fornirà diversi
orbitali, ad esempio l’orbitale 1s e l’orbitale 2s, ciascuno espresso come combinazione lineare delle
gaussiane primitive:
X10 X10
ψ1s = cj,1 gj ψ2s = cj,2 gj
j=1 j=1
dove cj,1 e cj,2 sono i coefficienti risultanti dal calcolo variazionale. Per costruire una contrazione,
si analizza il rapporto tra i coefficienti delle gaussiane nei due orbitali:
cj,1
Rj =
cj,2
Se per un certo gruppo di gaussiane il rapporto Rj rimane costante, significa che quelle gaussiane
contribuiscono in modo proporzionale sia all’orbitale 1s sia al 2s. In tal caso, si può
procedere a contrarle, cioè a combinarle in un’unica funzione, poiché entrano nei due orbitali con
lo stesso peso relativo. Ad esempio, se per due gaussiane consecutive si trova che:
c1,1 c2,1
= = 0.2
c1,2 c2,2
allora queste due gaussiane possono essere contratte insieme: la loro combinazione lineare descrive
correttamente il comportamento degli orbitali 1s e 2s senza perdita di precisione. In questo modo,
Dunning ha individuato la possibilità di combinare le funzioni gaussiane per costruire uno schema
di contrazione che può essere visto in forma compatta come:
In questa notazione:
• le parentesi tonde indicano le funzioni gaussiane primitive (GTO);
• le parentesi quadre indicano le funzioni contratte, cioè gli orbitali atomici effettivi
ottenuti dalla combinazione lineare delle primitive.
Naturalmente, laddove i rapporti tra i coefficienti non siano costanti, non ha senso combinare le
gaussiane.
122
parte studiando l’atomo isolato e costruendo una base che mantenga una flessibilità intrinseca, in
modo che, quando l’atomo si trova inserito in un contesto chimico diverso, la base possa adattarsi.
In questo modo si ottiene una base flessibile, utile e riutilizzabile in diversi contesti molecolari.
È importante ricordare che non si costruiscono mai espansioni gaussiane con l’obiettivo
di riprodurre esattamente le funzioni di tipo STO, piuttosto si cerca di trovare orbitali gaus-
siani contratti che sempre meglio le approssimano con il minimo possibile costo computazionale.
Nel programma Gaussian, tuttavia, esistono anche le cosiddette STO-NG, dove N indica il numero
di gaussiane utilizzate per approssimare una funzione di tipo Slater. Ad esempio, se si specifica
STO-6G, il calcolo utilizzerà sei gaussiane per ogni orbitale di tipo Slater. La più semplice tra queste
è la STO-3G, che viene spesso utilizzata dal programma GaussView quando si genera la geometria
di una molecola molto grande. Infatti, quando si salva l’input di una simile molecola per Gaussian,
si troverà nel file la base STO-3G, la più piccola disponibile. Ciò avviene perché il programma, per
default, evita di assegnare basi troppo grandi e computazionalmente pesanti, in modo da rendere
il calcolo almeno eseguibile. Naturalmente, l’utente è libero di modificare la base successivamente.
Tuttavia, le STO-NG hanno un valore scientifico limitato: sono basi minimali, poco flessibili
e quindi incapaci di descrivere accuratamente un legame chimico, tanto meno di rappresentare in
modo realistico variazioni molecolari o reazioni chimiche. Le loro energie risultano pertanto poco
affidabili. Le si può usare solo per molecole molto grandi e solo quando si desidera ottenere una
stima qualitativa, ma non per calcoli accurati.
A differenza di queste, esistono poi le basi di Dunning e quelle di Pople, come ad esempio
6-31G*, oppure le serie cc-pVnZ (tipiche di Dunning) o le ANO-L. Queste ultime sono basi molto
estese e sofisticate, costituite da schemi di contrazione complessi, con anche dieci gaussiane per sin-
golo orbitale atomico. Sono quindi basi pesanti dal punto di vista computazionale, da utilizzare
solo quando strettamente necessario — ad esempio per riprodurre calcoli già presenti in letteratura
o per ottenere risultati di elevata accuratezza — poiché ogni funzione atomica coinvolge numerose
gaussiane e il calcolo degli integrali risulta estremamente oneroso. Infatti, in tali casi, ogni integrale
può richiedere quattro sommatorie interne.
Tutte le basi note possono comunque essere consultate online: esistono siti dedicati che raccol-
gono una vera e propria tavola periodica interattiva delle basi. Cliccando su un elemento, è
possibile visualizzare tutte le basi disponibili per quell’atomo e persino scaricare il formato compat-
ibile con il programma di calcolo scelto. In questo modo non è necessario generare manualmente la
base: basta scaricarla e inserirla nel file di input del calcolo.
Scegliere una base più accurata cambia non solo il risultato numerico, ma anche le
proprietà geometriche della molecola.
Visti i dettagli dei basis set ed alcune regole di contrazione potremmo chiederci, quali sono dunque
i requisiti fondamentali che si chiedono a una base?
• buona descrizione dei gusci interni. I gusci interni sono occupati dagli elettroni più vicini
al nucleo, e l’energia totale del sistema dipende fortemente dalla corretta descrizione di tali
elettroni. Anche se gli elettroni dei gusci interni non partecipano direttamente ai processi
chimici o fotochimici - che coinvolgono principalmente gli elettroni di valenza - essi esercitano
comunque un’influenza sui gusci di valenza, modificandone parzialmente il comportamento e
le energie associate.
123
• Flessibilità della base. Un altro requisito essenziale è la flessibilità. La base deve essere
in grado di adattarsi a diversi contesti chimici: ad esempio, il carbonio e l’idrogeno presenti
nel metano si trovano in un ambiente elettronico molto diverso rispetto a quello del benzene.
Di conseguenza, una buona base deve essere sufficientemente flessibile da descrivere accu-
ratamente entrambi i casi, riproducendo in modo realistico sia le caratteristiche del legame
chimico che le proprietà elettroniche locali.
Le basi minime, come suggerisce il nome, non sono in grado di fornire una descrizione accurata di
tali stati. Possono al massimo essere impiegate per una stima preliminare, come l’ottimizzazione
della geometria di molecole molto grandi, ma non per calcoli accurati di energie o proprietà elet-
troniche.
124
Ad esempio:
• per il carbonio, che possiede orbitali 1s, 2s e 2p, le funzioni di polarizzazione sono di tipo d,
f , g, h, ecc.;
• per l’idrogeno, che possiede solo l’orbitale 1s, le funzioni di polarizzazione sono di tipo p, d,
ecc.
Queste funzioni sono fondamentali per:
1. descrivere con maggiore accuratezza i legami chimici;
2. rappresentare la risposta del sistema a un campo elettrico esterno;
3. migliorare la descrizione della correlazione elettronica.
Per chiarire, consideriamo la molecola di idrogeno. Le funzioni di tipo p allineate lungo l’asse di
legame possono contribuire a descrivere l’orbitale molecolare occupato e aiutare a rappresentare
correttamente la dissociazione a grandi distanze internucleari. In presenza di un campo elettrico
esterno, invece, la polarizzazione dell’atomo d’idrogeno non può essere descritta solo da una funzione
s (simmetrica rispetto al nucleo), ma richiede anche l’aggiunta di una funzione p che permetta lo
spostamento della densità elettronica nella direzione del campo.
Funzioni diffuse Oltre alle funzioni di polarizzazione, sono fondamentali anche le funzioni dif-
fuse. Esse hanno lo stesso momento angolare delle funzioni occupate, ma esponenti più piccoli, il
che le rende più estese nello spazio. Le funzioni diffuse sono importanti per rappresentare corretta-
mente le code della densità elettronica.
Per il carbonio, le funzioni diffuse sono di tipo s e p; per l’idrogeno, sono solo di tipo s. Queste
funzioni si rivelano particolarmente utili nella descrizione delle interazioni deboli e nella rappresen-
tazione accurata di specie anioniche o stati eccitati.
Come riconoscere funzioni diffuse e di polarizzazione Nelle basi di tipo Pople, le funzioni
diffuse e di polarizzazione sono indicate da simboli specifici:
• un + nella sigla (6-31+G*) indica l’aggiunta di funzioni diffuse sugli atomi pesanti;
• due segni ++ (6-31++G**) indicano che le funzioni diffuse sono aggiunte anche agli idrogeni;
• un singolo asterisco * aggiunge funzioni di polarizzazione solo agli atomi pesanti;
• due asterischi ** aggiungono funzioni di polarizzazione anche agli idrogeni.
Nelle basi di tipo Dunning, invece, le funzioni diffuse sono indicate dal prefisso aug- (da
augmented ), mentre le funzioni di polarizzazione sono comprese nella denominazione polarized
valence. Vediamo più nel dettaglio le basi della famiglia di Dunning indicate con la sigla
cc-pVnZ, dove:
• cc sta per correlation consistent;
• pV indica polarized valence, ossia che le funzioni di polarizzazione sono applicate solo agli
elettroni di valenza;
• nZ (dove n = D, T, Q, 5, . . .) rappresenta il numero di zeta, cioè il grado di duplicazione dei
gusci di valenza: doppia, tripla, quadrupla zeta, ecc.
125
Esempio: il caso dell’azoto Consideriamo, ad esempio, l’atomo di azoto. La configurazione
elettronica evidenzia una distinzione tra i gusci di valenza e il guscio interno (o core). In una
base cc-pVDZ (correlation consistent polarized valence double zeta), la denominazione “doppia zeta”
si riferisce esclusivamente alla duplicazione dei gusci nella regione di valenza, mentre il core rimane
invariato. Ciò significa che la base descrive con maggiore accuratezza gli elettroni più esterni, ma
non introduce funzioni aggiuntive per gli elettroni interni.
Basi con core polarizzato Esistono però anche basi che includono la polarizzazione del core,
indicate con la sigla cc-pCVnZ, dove “C” sta per core. In queste basi, vengono introdotte funzioni
supplementari dedicate agli elettroni interni, rendendo il core polarizzabile. Ad esempio, se uti-
lizziamo una base cc-pCVDZ per il carbonio, avremo due gusci dedicati anche al core: un orbitale
1s e un ulteriore 1s′ (più compatto), oltre a un p′ aggiuntivo.
Questa estensione della base risulta utile in contesti particolari, come:
• studi di effetti relativistici, nei quali il core deve essere descritto con maggiore accuratezza;
• simulazioni di spettri Auger o altri processi che coinvolgono la creazione di lacune nel
guscio interno (core holes), dove la polarizzabilità del core è cruciale per descrivere il rilas-
samento elettronico successivo all’eccitazione.
Scelta della base e costo computazionale La scelta della base dipende fortemente dallo scopo
del calcolo e dalle risorse computazionali disponibili. Non esiste una base “migliore in asso-
luto”: una base molto grande fornisce risultati più accurati, ma diventa rapidamente proibitiva
in termini di tempo e memoria. La dimensione della base influisce in modo diretto sui tempi di
calcolo, che aumentano in maniera esponenziale con la complessità del metodo impiegato.
Nel caso del metodo Hartree–Fock, la complessità computazionale scala approssimativamente come
M 4 , dove M rappresenta il numero di funzioni di base. Questo è già un costo significativo, ma nei
metodi post-Hartree–Fock, come la Configuration Interaction (CI) o la Coupled Cluster (CC), la
complessità può arrivare a scalare come M 7 o peggio, rendendo i calcoli molto onerosi.
In questi casi, la base non incide solo sul numero di funzioni, ma anche sulla complessità delle
operazioni, in particolare nella trasformazione degli integrali bielettronici dalla base atom-
ica alla base molecolare. Tale trasformazione coinvolge quattro indici e quindi richiede operazioni
su tensori a quattro dimensioni, con un costo computazionale estremamente elevato. Inoltre, la
costruzione della matrice hamiltoniana nel metodo post-Hartree–Fock dipende dal numero di deter-
minanti di Slater considerati, il quale cresce in modo esponenziale con il numero di orbitali.
126
Tuttavia, quando la molecola diventa molto grande, ad esempio nel caso del benzene, anche una
base di questo tipo può risultare eccessiva per i metodi ab initio più costosi. In tali situazioni, i
metodi di teoria del funzionale della densità (DFT) rappresentano un’ottima alternativa: i DFT,
infatti, hanno una complessità computazionale simile a quella di Hartree–Fock, ma consentono di
utilizzare basi più ampie senza un incremento drammatico del costo computazionale. Naturalmente,
la dimensione della molecola rimane comunque un fattore determinante.
• le funzioni di tipo S sono contratte secondo lo schema 6s, 3s, 1s, 1s, con un’ulteriore funzione
diffusa (indicato dal segno “+”);
• le funzioni di tipo P seguono lo schema 3p, 1p, 1p, anch’esse con una funzione diffusa (“+”);
• la presenza degli asterischi “**” segnala l’aggiunta di funzioni di polarizzazione, che per
l’azoto corrispondono a funzioni di tipo D.
NN = (5 × 1) + (4 × 3) + (1 × 5) = 5 + 12 + 5 = 22
127
Per ciascun atomo di idrogeno abbiamo dunque:
NH = (4 × 1) + (1 × 3) = 4 + 3 = 7
La molecola di ammoniaca contiene un atomo di azoto e tre atomi di idrogeno, quindi:
Ntotale = NN + 3 × NH = 22 + 3 × 7 = 22 + 21 = 43
Il numero totale di funzioni atomiche di base utilizzate nel calcolo della molecola di ammoniaca con
la base 6-311++G** è quindi pari a:
Ntotale = 43
Lo schema di contrazione per questo sistema può quindi essere rappresentato complessivamente nel
modo seguente:
Azoto (N): 12S → 5S, 6P → 4P, 1D → 1D
Idrogeno (H): 6S → 4S, 1P → 1P
Questo esempio evidenzia come l’utilizzo di una base di gaussiane permetta di rappresentare in
modo efficiente gli orbitali atomici.
Conlcudiamo dicendo che quindi le funzioni gaussiane, pur essendo semplici dal punto di vista
matematico, contratte e combinate in blocchi per formare funzioni più complesse, sono in grado
di riprodurre accuratamente gli orbitali atomici effettivi. Lavorare invece con gaussiane non con-
tratte (cioè “libere”) sarebbe estremamente costoso dal punto di vista computazionale, poiché il
numero di funzioni da gestire crescerebbe enormemente. Il processo di contrazione, quindi, riduce
drasticamente il costo computazionale mantenendo comunque un’elevata accuratezza.
Infine, è il calcolo variazionale quantomeccanico stesso che, nel corso della procedura, combina
e rimescola i blocchi contratti nel modo ottimale per fornire l’energia più bassa possibile e quindi
gli orbitali molecolari più accurati. Per questo motivo, non è opportuno cercare di riprodurre
manualmente funzioni del tipo STO tramite gaussiane: è preferibile lasciare che sia il
metodo variazionale a determinare la combinazione ottimale delle funzioni di base.
Supponiamo di avere una base di funzioni atomiche di dimensione m. In questo caso, poiché
la molecola H2 è costituita da due elettroni, la funzione d’onda migliore possibile può essere scritta
in modo relativamente semplice:
X αβ − βα
Ψbest (r1 , r2 ) = Ckl [χk (r1 )χl (r2 ) + χl (r1 )χk (r2 )] √
k≤l
2
128
dove abbiamo assunto lo stato di spin di singoletto (e quindi le combinazioni simmetriche) come
stato a più bassa energia. Abbiamo inoltre che:
• χk e χl sono le funzioni della base atomica
• Ckl sono coefficienti variazionali determinati minimizzando l’energia
Poiché disponiamo di m funzioni di base, il numero totale di coefficienti variazionali Ckl sarà pari
a:
m(m + 1)
Nparam = −1
2
dove il termine (−1) tiene conto della condizione di normalizzazione della funzione d’onda, che
riduce di un’unità i gradi di libertà effettivi. Per la molecola H2 , anche utilizzando basi di dimen-
sioni relativamente grandi, il calcolo rimane computazionalmente accessibile, poiché il numero di
parametri resta modesto.
Sfruttando il metodo RHF, gli elettroni vengono collocati negli orbitali a più bassa energia disponi-
bili. Siccome inoltre abbiamo solo due elettroni e vogliamo una parte spaziale simmetrica la funzione
d’onda RHF può essere espressa semplicemente come il seguente determinante di Slater:
αβ − βα
ΦRHF (r1 , r2 ) = ϕ1 (r1 )ϕ1 (r2 ) √
2
Supponiamo ora di utilizzare una base atomica molto ampia, comprendente funzioni 1s, 2s, 2p
per ciascun atomo. Indichiamo con χ1s , χ2s , χ2p le funzioni centrali sull’atomo 1, e con χ′1s , χ′2s ,
χ′2p quelle centrate sull’atomo 2. Siccome la molecola di idrogeno possiede un centro di simmetria,
allora la funzione orbitale ϕ1 , rispettando tale simmetria, conterrà un contributo dominante dalla
combinazione simmetrica delle funzioni 1s:
ϕ1 ≈ N (χ1s + χ′1s )
e cosı̀ via per tutte le funzioni compatibili con la simmetria del sistema.
129
Per RHH → ∞ l’energia dovrebbe tendere al valore corrispondente a due atomi di idrogeno
isolati, ovvero E = −1 unità atomica (dato che per un singolo atomo di idrogeno E = − 12 ).
Tuttavia, per distanze molto grandi, il metodo RHF può sbagliare anche di 0.2 a.u. . Ciò mostra
come il metodo Hartree-Fock fallisca completamente nel descrivere correttamente la dissociazione
della molecola a grande distanza, dove gli elettroni dovrebbero localizzarsi su ciascun atomo.
Per vedere il motivo di ciò vediamo che succede quando i due atomi di idrogeno si trovano a distanza
molto grande (RHH → ∞) in RHF; in questo caso l’orbitale ϕ1 tende alla forma:
1
lim ϕ1 = √ (χ1s + χ′1s )
RHH →∞ 2
Il fattore √12 compare perché, a grandi distanze, le due funzioni χ1s e χ′1s sono sostanzialmente
ortogonali: la loro sovrapposizione tende a zero, e quindi la combinazione simmetrica deve essere
normalizzata con tale coefficiente. Quindi, la funzione RHF sarà:
1 αβ − βα
ΦRHFRHH →∞ = (χ1s + χ′1s ) (χ1s + χ′1s ) √
2 2
Quando le due funzioni χ1s e χ′1s diventano disgiunte (nessuna sovrapposizione), la fun-
zione ϕ1 risulta estesa su entrambi gli atomi e non riesce più a descrivere correttamente
la localizzazione elettronica. Di conseguenza, l’energia calcolata con il metodo Hartree–Fock
diventa superiore a quella corretta dei due atomi separati. Vediamolo calcolando l’energia Hartree-
Fock in queste condizioni limite. Ricordiamo che entrambi gli elettroni occupano lo stesso orbitale
ϕ1 , e che l’operatore monoelettronico ĥ, in unità atomiche, è dato da:
1 1 1
ĥ = − ∇2 − −
2 riA riB
per ognuno dei due elettroni del sistema. L’energia Hartree-Fock totale risulta quindi:
1
EHF = 2⟨ϕ1 |ĥ|ϕ1 ⟩ + ⟨ϕ1 (r1 )ϕ1 (r2 )| |ϕ1 (r1 )ϕ1 (r2 )⟩
r12
130
Va inoltre ricordato che esiste anche un termine di repulsione nucleare, ma questo, a distanza
infinita, tende a zero. Analizziamo dunque i vari termini della parte monolettronica e bielelettronica
dell’hamiltoniano nelle condizioni di distanza infinita. La parte monolettronica diventa:
1 1 1 1 1
⟨ √ (χs + χs′ ) | − ∇2 − − | √ (χs + χs′ )⟩
2 2 |r − RH | |r − RH ′ | 2
Poiché le due funzioni sono centrate su nuclei molto distanti, la loro sovrapposizione tende a zero:
R →∞
⟨χS |χS ′ ⟩ −−HH
−−−−→ 0
Emono = 2EH
ossia due volte l’energia dell’atomo di idrogeno isolato. Il contributo in eccesso, che distingue
il comportamento del calcolo HF da quello ideale, è quindi rappresentato dal termine
di repulsione interelettronica residuo. Analizziamo dunque questo termine nel dettaglio:
1 1
Eee = ⟨(χs + χs′ )(χs + χs′ )| |(χs + χs′ )(χs + χs′ )⟩
4 r12
dove abbiamo esplicitato i fattori di normalizzazione:
4
1 1
√ =
2 4
Ricordiamo che i termini incrociati contenenti χS χS ′ vanno a zero per RHH → ∞, per cui soprav-
vivono solo i termini puri:
1 1 1
Eee ≈ ⟨χs χs | |χs χs ⟩ + ⟨χs′ χs′ | |χs′ χs′ ⟩
4 r12 r12
I due termini sono identici, poiché ciascuno rappresenta l’interazione elettrostatica tra due
“mezzi elettroni” localizzati sullo stesso atomo. Il fattore 14 deriva dal prodotto 12 × 12 , che
riflette il fatto che, su ogni centro atomico, si ha una densità elettronica corrispondente a metà
elettrone. Questi due mezzi elettroni interagiscono tra loro, generando un contributo
residuo che non può annullarsi nemmeno all’infinito, e che rappresenta l’errore statico
del metodo Hartree–Fock. L’errore che abbiamo evidenziato in precedenza viene definito statico,
poiché non dipende dalla distanza internucleare; esso infatti anche all’aumentare della distanza,
dalla distanza di equilibrio fino a RHH → ∞, permane dando un contributo energetico extra che non
dovrebbe essere presente. Questo termine, che rimane costante anche a distanza infinita, introduce
un errore significativo sull’energia totale e deve essere assolutamente eliminato. Di
conseguenza, non è possibile utilizzare il metodo RHF per descrivere correttamente la dissociazione
della molecola di idrogeno, H2 , o, più in generale, per studiare una dissociazione omolitica di un
131
legame chimico. Nel caso di una dissociazione omolitica, infatti, il legame viene rotto mantenendo
un elettrone su ciascun frammento risultante. In queste condizioni, il metodo HF produce lo stesso
tipo di errore statico, poiché impone che entrambi gli elettroni occupino il medesimo orbitale spaziale
a qualsiasi distanza. Pertanto, il calcolo RHF non è adatto allo studio di processi di dissociazione
omolitica: esso fallisce completamente nel rappresentare correttamente la separazione di un legame
in due frammenti neutri. Viceversa, può fornire risultati accettabili solo in caso di dissociazioni
eterolitiche, dove gli elettroni rimangono localizzati su un unico centro. Nel caso di dissociazione
omolitica, l’errore è grande e non trascurabile, specialmente alle grandi distanze.
Come correggere il problema. Per eliminare questo errore, è necessario consentire agli elettroni
di riorganizzarsi autonomamente quando la distanza internucleare cresce, in modo che a distanza
infinita ciascun elettrone possa localizzarsi su un diverso centro atomico.
Affrontiamo il problema adottando un approccio non ristretto, l’Unrestricted Hartree–Fock
(UHF), che permette agli elettroni di occupare orbitali spaziali differenti. Per semplificare la dis-
cussione, consideriamo una base minima, sufficiente a mettere in evidenza il fenomeno senza
introdurre eccessiva complessità formale. L’errore del calcolo RHF deriva infatti dal vincolo im-
posto agli elettroni di condividere lo stesso orbitale spaziale a tutte le distanze. Poiché la forma
dell’orbitale non cambia, l’errore associato resta invariato e non può essere eliminato. Per ovviare
a ciò, occorrono almeno due orbitali distinti, in modo da consentire agli elettroni di spin opposto
β
di localizzarsi su centri diversi. Definiamo quindi due orbitali, ϕα 1 e ϕ1 , ciascuno espresso come
combinazione lineare delle funzioni atomiche χS e χS ′ :
ϕα
1 = χS + λ χS ′ , ϕβ1 = χS ′ + λχS
dove λ è un parametro variazionale (reale e compreso tra 0 e 1) che controlla il grado di mescolamento
tra i due centri atomici H e H’. Il parametro λ sarà determinato minimizzando l’energia totale in
funzione della distanza internucleare RHH .
β
Figure 2.6.1: Per simmetria, ϕα1 e ϕ1 presentano una distribuzione spaziale analoga: ciascuno è
fortemente localizzato su un centro e presenta un contributo minore sull’altro.
β
La figura qualitativa che ne risulta mostra ϕα
1 (continua) concentrata su H e ϕ1 (tratteggiata)
su H’. Notiamo che i due orbitali non sono ortogonali, e non è affatto necessario che lo siano. Anzi,
132
non vogliamo imporre ortogonalità, poiché ciò limiterebbe la flessibilità del sistema: i due orbitali
possono parzialmente sovrapporsi, contribuendo entrambi alla descrizione del legame.
Possiamo costruire quindi la funzione d’onda UHF come un determinante di Slater formato dai due
spin-orbitali:
β
1 ϕα 1 (r1 )α(1) ϕ1 (r1 )β(1)
ΨUHF = √ α β
2 ϕ1 (r2 )α(2) ϕ1 (r2 )β(2)
che può essere esplicitamente scritto come:
β β
ΦUHF = ϕα α
1 (1)α(1)ϕ1 (2)β(2) − ϕ1 (1)β(1)ϕ1 (2)α(2)
Osserviamo immediatamente che questa funzione d’onda non si fattorizza in una parte spaziale e
β
in una parte di spin. Infatti, a differenza del caso Hartree–Fock ristretto, gli orbitali ϕα
1 e ϕ1 sono
diversi, e ciò impedisce la separazione del termine di spin del tipo:
(α(1)β(2) − β(1)α(2))
Questo aspetto è cruciale: la funzione d’onda UHF non ha una parte di spin separabile e presenta
una correlazione intrinseca tra coordinate spaziali e di spin.
ϕα
1 = χS + λχS ′ , ϕβ1 = χS ′ + λχS
Tuttavia, questa scelta non è univoca: si può infatti costruire un insieme equivalente di orbitali
scambiando il ruolo di α e β, cioè:
ϕ̃α
1 = χS ′ + λχS , ϕ̃β1 = χS + λχS ′
β β
Queste due configurazioni, {ϕα α
1 , ϕ1 } e {ϕ̃1 , ϕ̃1 }, sono completamente equivalenti dal punto di vista
energetico e fisico. Gli indici α e β non determinano una scelta preferenziale: semplicemente
rappresentano due possibili orientamenti della stessa soluzione, ottenuti per scambio di centri. Di
conseguenza, si possono costruire due determinanti UHF equivalenti:
β
ΦUHF = ϕα α
1 (1)α(1)ϕ1 (2)β(2) − ϕ1β (1)β(1)ϕ1 (2)α(2)
β β
Φ̃UHF = ϕ̃α α
1 (1)α(1)ϕ̃1 (2)β(2) − ϕ̃1 (1)β(1)ϕ̃1 (2)α(2)
Entrambe le funzioni d’onda rappresentano soluzioni equivalenti del problema UHF. Possiamo utiliz-
zare una di esse per analizzare il comportamento energetico in funzione del parametro variazionale λ.
133
• a distanza infinita, vogliamo che ciascun elettrone si localizzi su un centro diverso; ciò cor-
risponde a λ = 0, poiché in tal caso ϕ1α = χS e ϕ1β = χS ′ .
Pertanto, λ(R) diminuisce da 1 a 0 con l’aumentare della distanza. Tuttavia, tale variazione non è
continua: a un certo valore di R (circa R ≈ 2.7 bohr) si verifica una brusca discontinuità. Questo
punto è noto come punto di Coulson–Fischer.
Analogamente, la funzione con la tilde è semplicemente rovesciata nei ruoli dei centri:
Anche in questo limite, la funzione non si fattorizza in una parte spaziale e una parte di spin.
Tuttavia, il vantaggio essenziale è che questa soluzione fornisce l’energia corretta: il termine
bi-elettronico di interazione tra i due elettroni scompare, poiché gli orbitali χS e χS ′ non si sovrap-
pongono più. Pertanto, l’energia tende correttamente a quella di due atomi di idrogeno isolati.
Curva energetica e rottura di simmetria.
134
La curva di energia ottenuta con il metodo UHF coincide con quella RHF fino al punto di
Coulson–Fischer; oltre tale punto, l’energia UHF continua a decrescere correttamente verso il limite
asintotico E = −1 (in unità atomiche), mentre la soluzione RHF rimane troppo alta.
Conseguenze sulla natura dello stato di spin. Quando la simmetria si rompe, lo stato di
spin non è più ben definito: la funzione d’onda UHF non è un puro singoletto né un puro tripletto,
ma una loro combinazione. Nel limite di grande distanza, gli elettroni risultano disen-
tangled : non sono più correlati tra loro, e il comportamento di ciascuno è indipendente dall’altro.
L’entanglement, infatti, presuppone uno stato di spin ben definito (singoletto o tripletto puro), in
cui le due particelle rimangono quantisticamente collegate anche a distanza. Nel caso UHF, invece,
gli elettroni sono separati e non presentano più correlazione di spin a lunga distanza.
Come visto, esistono due soluzioni UHF equivalenti: ΦUHF e Φ̃UHF . Consideriamo la loro somma e
la loro differenza. La somma delle due funzioni dà:
Φ+ = ΦUHF + Φ̃UHF ∝ (χS (1)χS ′ (2) + χS ′ (1)χS (2)) (α(1)β(2) − β(1)α(2))
In questo caso, la funzione d’onda si fattorizza perfettamente in una parte spaziale simmetrica e una
parte di spin antisimmmetrica, corrispondente a un singoletto con S = 0 e mS = 0. La differenza
delle due funzioni, invece, dà:
Φ− = ΦUHF − Φ̃UHF ∝ (χS (1)χS ′ (2) − χS ′ (1)χS (2)) (α(1)β(2) + β(1)α(2))
che rappresenta una funzione tripletto con S = 1 e mS = 0.
Pertanto, le due soluzioni UHF equivalenti possono essere combinate linearmente per ottenere
stati puri di spin ben definito: una combinazione simmetrica (singoletto) e una antisimmmetrica
(tripletto). Da ciò si deduce che ciascuna delle soluzioni UHF di partenza è in realtà una mescolanza
di singoletto e tripletto con lo stesso peso. Dunque, avremo:
1
Ψ1 = √ (ΦU HF + Φ̃U HF )
2
135
1
Ψ3 = √ (ΦU HF − Φ̃U HF )
2
Possiamo quindi riscrivere:
1
ΦU HF = √ (Ψ1 + Ψ3 )
2
1
Φ̃U HF = √ (Ψ1 − Ψ3 )
2
Abbiamo dunque una mescolanza di singoletto e tripletto √ con uguale peso, determinato dal
modulo quadrato del coefficiente di mescolamento, che vale 1/ 2. Di conseguenza, la composizione
è 0.5 singoletto e 0.5 tripletto.
136
In entrambi i casi, si mantiene spin totale pari a zero, ma gli elettroni risultano localizzati su un
solo centro.
Se ora guardiamo al caso RHF (Restricted Hartree-Fock), vediamo che anch’essa rappresenta
una mescolanza di questi contributi, ma con un peso fisso che non varia con la distanza tra i due
atomi. Infatti, anche a distanza infinita — dove i contributi ionici dovrebbero sparire completa-
mente — nella RHF questi pesano quanto la parte covalente, il che è fisicamente errato. La funzione
RHF, infatti, può essere scritta come:
1
ΦRHF = √ (χS + χS ′ )(χS + χS ′ )(αβ − βα)
2
Qui, gli orbitali molecolari sono combinazioni simmetriche con uguale peso dei due orbitali atomici,
per cui il contributo covalente e quello ionico risultano sempre bilanciati, indipendentemente dalla
distanza.
αβ − βα
ΦRHF = [(χS χS ′ + χS ′ χS ) + (χS χS + χS ′ χS ′ )] √
2
Questo spiega l’andamento energetico scorretto della curva RHF: a grande distanza, infatti, l’energia
tende a quella di mezzo elettrone per atomo, un risultato non fisico. La meccanica quantistica
consente di scrivere matematicamente tale combinazione, ma il risultato non corrisponde alla realtà
del sistema.
Possiamo interpretare meglio il tutto considerando una base minima composta da due orbitali:
• un orbitale legante, riempito;
• e un orbitale antilegante, virtuale.
Nel calcolo RHF, l’orbitale legante (il primo orbitale, a energia più bassa) è una combinazione sim-
metrica χS + χS ′ , mentre l’orbitale virtuale (a energia più alta) è la combinazione antisimettrica
χS − χS ′ , necessaria per mantenere l’ortogonalità tra gli orbitali.
Se poniamo entrambi gli elettroni nell’orbitale legante, otteniamo lo stato fondamentale RHF,
mentre se li poniamo nell’orbitale virtuale otteniamo uno stato eccitato, più alto in energia.
αβ − βα
ΦRHF = ϕ1 ϕ1 √
2
αβ − βα
Φecc
RHF = ϕ2 ϕ2 √
2
Siamo per il momento abituati a visualizzare la situazione in questi termini. Tuttavia, occorre fare
attenzione: quando la distanza interatomica tende all’infinito (RH → ∞), le energie orbitali
del più alto orbitale occupato (HOMO) e del più basso orbitale non occupato (LUMO), cioè ε1 e
137
ε2 , diventano degeneri.
In altre parole:
ε1 = ε2 per RH → ∞
In tali condizioni, le due funzioni orbitali associate hanno la stessa energia. Di conseguenza, i due
elettroni del sistema potrebbero occupare indifferentemente uno dei due orbitali, poiché l’energia
non cambia.
In questo caso, la funzione d’onda ottenuta nel formalismo RHF (Restricted Hartree-Fock) as-
sume lo stesso valore dell’energia della funzione eccitata, ovvero quella con entrambi gli elettroni
collocati nel secondo orbitale (LU M O).
In queste condizioni, i due determinanti (quello RHF e quello eccitato) non sono autofunzioni
dell’hamiltoniano elettronico, e dunque non diagonalizzano l’operatore. Pertanto, il valore ot-
tenuto è un semplice valore d’aspettazione e non un autovalore esatto.
Il termine variazionale ci viene in aiuto: ricordiamo che nel principio variazionale lineare,
se abbiamo una base di funzioni poli-elettroniche, la funzione d’onda totale può essere espressa
come una combinazione lineare:
Ψ = C1 ΦRHF + C2 Φeccitata
Tale combinazione è sempre un miglioramento rispetto a ciascuna funzione singola, poiché i coef-
ficienti C1 e C2 fungono da parametri variazionali lineari. Per determinarli, diagonalizza una
matrice 2 × 2:
! !
C1 C1
ERHF H12
=E
H21 H22 C2 C2
Diagonalizzando questa matrice, otteniamo gli autovalori (energie) e gli autovettori (i coefficienti
C1 , C2 ) corrispondenti. In particolare:
• ERHF rappresenta l’energia della funzione RHF (energia del determinante di riferimento);
• H22 è l’energia del determinante eccitato;
• H12 = H21 è l’elemento fuori diagonale, responsabile del mescolamento tra i due determi-
nanti.
138
Quando la distanza tra i nuclei diventa molto grande, i due coefficienti C1 e C2 risultano uguali in
modulo (a parte il segno):
1
|C1 | = |C2 | = √
2
Questo significa che, a distanza infinita, i due determinanti contribuiscono in modo equivalente alla
funzione d’onda del sistema.
In definitiva, anche in questo caso ci ritroviamo in una situazione analoga a quella vista preceden-
temente: a grande distanza, i due determinanti (RHF e eccitato) partecipano in egual misura
alla descrizione del sistema. La singola funzione RHF, dunque, non è più sufficiente: occorre com-
binare almeno due determinanti per ottenere una rappresentazione accurata.
Si tratta quindi della descrizione più semplice possibile — una funzione costruita dalla combi-
nazione di due determinanti, che corrisponde sostanzialmente alla funzione di singoletto già
incontrata in precedenza.
139
Il profilo energetico risultante mostra una curva intermedia tra quella ottenuta con RHF e
quella esatta (indicata come best). Tale curva intermedia separa l’errore della funzione RHF in due
contributi distinti che chiamiamo ∆1 (inferiore) e ∆2 (superiore):
Ebest − ERHF = ∆1 + ∆2
140
Contributo statico
Il contributo statico rappresenta l’errore intrinseco della funzione RHF, in particolare legato alla
restrizione della singola determinante. A grande distanza tra i nuclei, il contributo statico
diventa predominante e assume un valore notevole (∼ 0.2 a.u.), mentre il contributo dinamico è
trascurabile (∼ 0).
Alla distanza di equilibrio, invece, il contributo statico e quello dinamico diventano comparabili,
ciascuno con un valore di circa 0.02 a.u., cioè un decimo di quello statico a grande distanza.
Contributo dinamico
Il contributo dinamico è legato alla capacità della funzione d’onda di descrivere correttamente le
correlazioni tra elettroni a breve distanza. In altre parole, mentre il contributo statico corregge
l’errore della singola determinante, il contributo dinamico emerge solo quando si includono ulteriori
determinanti oltre quelli necessari per correggere la degenerazione a grande distanza.
In sintesi nel caso della molecola H2 , la separazione in contributi statico e dinamico può essere
ottenuta utilizzando due determinanti:
• Il determinante RHF singolo rappresenta il caso peggiore, con una forte energia di correlazione
statica a grandi distanze.
• La funzione BEST, costruita includendo tutti i determinanti possibili basati sulla base atom-
ica scelta, rappresenta la descrizione più accurata e consente di distinguere chiaramente i
contributi statico e dinamico.
• Un’approssimazione intermedia può essere realizzata con due determinanti, costruiti in
modo da formare un singoletto, in cui gli elettroni rimangono correttamente correlati (en-
tangled ) a spin opposto.
È fondamentale sottolineare che, per preservare la corretta descrizione della correlazione elettronica,
la funzione deve essere un singoletto, e non un approccio unrestricted con elettroni non entangled.
Questo garantisce che il contributo statico venga corretto e che l’energia dinamica possa essere
successivamente migliorata con l’inclusione di ulteriori determinanti.
Questa funzione ci da la curva di energia migliore possibile; quella invece ottenuta con RHF è
consistente con la best nel minimo ma assume valori troppo elevati nella dissociazione. La differenza
di Energia è:
141