Metodi ab initio

Ogni elettrone è rappresentato da una funzione d'onda, il cui quadrato corrisponde ad una densità di probabilità, rappresentata da una funzione di distribuzione delle coordinate. La distribuzione gaussiana è un esempio di distribuzione.

Figure 2.1: Distribuzioni Gaussiane $G(x) = e^{(-x^2/\sigma ^2)}$

In Figura sono rappresentate diverse funzioni di distribuzione Gaussiane in un grafico unidimensionale: Sulle code, la probabilità è piccola. Come mostrato in Figura, la distribuzione Gaussiana può essere piú o meno slargata (cioè delocalizzata) a seconda del parametro $\sigma$. Piú grande é $\sigma$ e maggiore sará la larghezza. Una distribuzione gaussiana tridimensionale da luogo ad un "orbitale" di tipo s. La distribuzione tridimensionale relativa all'orbitale di tipo s viene rappresentata attraverso superfici di equiprobabilitá, che consistono in superfici sferiche. La distribuzione Gaussiana e' un esempio di distribuzione localizzata: la probabilitá di trovare l'elettrone e' minima al di fuori di un certo range (per esempio per la Gaussiana verde di Figura 2.1, ho scarse probabilitá di trovare l'elettrone per $x < -2$ e per $x > 2$).

In una molecola ci sono molti elettroni. L'equazione di Schroedinger è l'equazione cui obbedisce la funzione d'onda globale polielettronica della molecola, indicata con la notazione $\Psi(1,2,3,..n)$. L'equazione di Schroedinger possiamo sinteticamente scriverla cosí

$$H(R,r)\Psi + G(\Psi,R,r) \Psi = E\Psi$$ (2.1)

$H$ e $G$ sono due operatori^2.1(cioè ``cose'' che trasformano la funzione $\psi$ in un altra funzione). L'operatore H dipende da $r$, ovvero dalla totalitá delle coordinate elettroniche, e da R, ovvero dalla configurazione dei nuclei; G dipende sia da R che da $\Psi$ stessa. Come vedremo in seguito, nel caso in cui l'operatore che agisce sulla funzione soluzione dell'equazione, come nella Eq. 2.1, dipende dalla soluzione stessa $\Psi$, la soluzione $\Psi$ puó essere ottenuta solo iterativamente: ad ogni iterazione si ottiene una soluzione $\Psi$ che permette di calcolore un nuovo e piú accurato operatore G e cosí via fino a quando $\Psi$, e dunque neppure G, non cambiano piú. $E$ è l'energia associata a $\Psi$. Nella equazione 2.1 $R$ rappresenta l'insieme delle coordinate dei nuclei. Queste coordinate sono parametri fissi. La soluzione della equazione sopra per $\Psi$ ed $E$ è diversa se cambiamo questi parametri, cioè se cambiamo le coordinate dei nuclei. Per esempio $E$ per una molecola di idrogeno con il legame H-H ad 1 Å , è minore della energia $E$ per la stessa molecola con la distanza H-H a 1.5 Å. La soluzione esatta della equazione 2.1 non è conosciuta esattamente (eccetto che per l'atomo di idrogeno). Per arrivare ad una soluzione approssimata che ci dia E e $\Psi$, si assume che la funzione d'onda polielettronica sia data da un prodotto di funzioni d'onda monoelettroniche, dette orbitali molecolari:

$$\Psi(1,2,3,4,5...n) = \psi_{1}(1)\psi_{2}(2).. \psi_{n}(n)$$ (2.2)

Nella Equazione2.2^2.2 la notazione $\psi_{1}(1)$ sta ad indicare che quella funzione (per il momento sconosciuta) dipende dalle sole coordinate dell'elettrone 1. Gli elettroni, come qualunque altro corpo nello spazio, sono caratterizzati da tre coordinate cartesiane $x,y,z$ e da una quarta coordinata detta di spin che puó assumere solo due valori (su e giú oppure up and down oppure ancora $\alpha$ e $\beta$). In tutti i metodi ab initio si suppone inoltre che l'orbitale molecolare $\psi$ riferito all'elettrone 1 (come ad un qualunque altro elettrone) sia data da una combinazione lineare di funzione di base atomiche:

$$\psi_{i}(x) = c_1 \phi_{1}(x) + c_2 \phi_{2} (x) + ...c_n \phi_{n} (x),$$ (2.3)

dove le $\phi$ sono funzioni centrate sugli atomi chiamate orbitali atomici, ed i coefficienti $c$ sono degli scalari.

la combinazione lineare a destra rappresenta il BASIS SET, cioè le funzioni di distribuzione atomiche con la quale rappresentiamo l'orbitale molecolare monoelettronico generico $\psi_{i}$. Per ottenere altri orbitali molecolari basta variare solo i coefficienti $c$ tenendo fissi gli orbitali atomici $\phi(x)$. Per questo l'insieme degli orbitali atomici $\phi_i$ si chiama BASE. Per esempio nella Figura 2.2, data la base $g_{1}(x), g_{2}(x)$, possiamo ottenere due distribuzioni diverse (cioè due orbitali molecolari diversi) con due set di coefficienti diversi: nel primo set ($c_1 = c_2$) si ottiene una distribuzione delocalizzata, il che significa semplicemente che l'elettrone può essere trovato con eguale probabilità intorno alla posizione del primo nucleo (2Å  in Figura) o del secondo nucleo (4Åin Figura) Nel secondo caso ( $c_ 1 » c_2$) la probabilità di trovare l'elettrone sul nucleo 1 è molto piú alta: l'elettrone è in questo caso localizzato.

Figure 2.2: combinazioni lineare di funzioni Gaussiane

Se andiamo a sostituire le equazione 2.2 e 2.3 nella equazione 2.1, troveremo dopo qualche passaggio complicato che la equazione di Schroedinger per gli elettroni a nuclei fissi 2.1 si riduce ad una equazione per i coefficienti $c_{1} ..... c_n$ [Per maggiori dettagli sulla procedura matematica, consultare il seguente link (in inglese)]. Infatti - come abbiamo detto - le funzioni di BASE non variano. Facciamo variare la $\Psi$ variando solo i coefficienti. Di questi coefficienti ce ne sono $n$ per ogni orbitale molecolare, le funzioni d'onda monoelettroniche $\psi$. Nella risoluzione iterativa della equazione 2.1 si ottengono tanti orbitali molecolari monoelettronici (funzioni $\psi$ (vedi Eq 2.2), quanti sono le funzioni di base atomiche $\phi$ (vedi Eq. 2.3). Siccome gli orbitali molecolari devono essere almeno tanti quanti sono gli elettroni, la base delle funzione atomiche $\phi$ deve avere un numero sufficiente di elementi sí da poter allocare nei risultanti orbitali $\psi$ gli elettroni che compongono il sistema in oggetto.

Come si risolve per i coefficienti l'equazione 2.1? Prima di tutto occorre mettersi in testa che vi sono molte $\Psi$ (ed $E$ associate) che soddisfano alla equazione 2.1. A noi interessa la $\Psi$ che da luogo all'energia piú bassa. Quella $\Psi$ rappresenta la funzione di distribuzione elettronica del cosiddetto ``stato fondamentale''. Le altre $\Psi$, che ci interessano meno, rappresentano le funzioni di distribuzione elettroniche associate ai cosiddetti ``stati eccitati''. Perché ci interessa lo stato fondamentale? Perché nelle condizioni di temperatura e di pressione alle quali ci troviamo, il 99.9999% delle molecole si trova in questo stato. Tornando all'equazione 2.1, esiste un teorema che ci permette di arrivare alla soluzione: il teorema variazionale. Questo teorema stabilisce che ogni $\Psi$ - comunque sia scelta - dá luogo ad un valore di $E$ maggiore o uguale al valore vero di $E$ per lo stato fondamentale. Dunque basta determinare quei particolari coefficienti $c$ che minimizzano $E$. Sulla base del teorema variazionale, infatti, quello sarà il valore piú vicino possibile al valore vero dello stato fondamentale. In termini matematici per trovare il minimo di E occorre calcolare la derivata di E rispetto ai coefficienti ed eguagliarla a zero. Dunque saltano fuori tante equazione quanti sono i coefficienti (tante equazioni quante sono le funzioni del BASIS SET). Il sistema di equazioni è risolto da $n$ set degli $n$ coefficienti c del tipo

$$\psi_1 = c_{11} \phi_{1} + c_{12} \phi_{3} .. c_{1n} \phi_{n}$$

$$\psi_2 = c_{21} \phi_{1} + c_{22} \phi_{3} .. c_{2n} \phi_{n}$$

$$.....$$

$$\psi_3 = c_{n1} \phi_{1} + c_{n2} \phi_{3} .. c_{nn} \phi_{n}$$ (2.4)

Ogni set corrisponde ad un orbitale molecolare (ad un spin-orbitale piú precisamente). Ogni orbitale molecolare $\psi_i$ è caratterizzato da una sua energia $\epsilon_i$, anch'essa ottenuta dalla risoluzione del sistema. La funzione d'onda globale polielettronica sarà data dal prodotto $\psi_{1}(1) ...\psi_{n}(n)$ e l'energia totale (di stato fondamentale) sarà data dalla somma delle energie degli (spin)orbitali molecolari. Per la particolare forma degli operatori $H$ e $G$, ogni orbitale molecolare accoglie in realtà due elettroni. Se nella molecola ci sono N elettroni ci devono pertanto essere almeno N/2 funzioni spaziali di base che corrispondono ad N funzioni spin-orbitaliche di base (la coordinata spin puó assumere solo due valori). In realtà la base cosiddetta minima contiene un numero di orbitali atomici maggiore di N/2 cosicché alcuni degli orbitali molecolari risultanti dalla soluzione della Eq. 2.4 risulteranno non occupati o virtuali. Il passaggio di uno degli elettroni dagli orbitali molecolari occupato a questi orbitali virtuali dá luogo ad uno stato eccitato della molecole. Discutiamo un esempio pratico. Nella molecola di acqua ci sono 10 elettroni (8 dell'ossigeno) e 2 dagli idrogeni. La base minima si compone degli orbitali atomici 1s, 2s, 2px, 2py, 2pz appartenenti all'ossigeno (per un totale di 5 funzioni orbitali) e degli orbitali 1s appartenenti agli idrogeni, per un totale quindi di sette orbitali, 5 sull'ossigeno ed uno ciascuno sui due idrogeni. In questi sette orbitali, nello stato fondamentale (che corrisponde alla soluzione della Eq. 2.4 utilizzando il metodo variazionale), dovremo allocare 10 elettroni. Datosi che per ogni orbitale molecolare vi sono due elettroni (con spin opposti), avremo pertanto 5 orbitali occupati e due orbitali virtuali vuoti. Dalla funzione d'onda polielettronica risultante dal prodotto dei 10 spin-orbitali occupati (o meglio dal suo quadrato che corrisponde alla densità) si possono calcolare tutte le proprietà della molecola (momento di dipolo, modi di vibrazione, plausibilità, etc. ). Resta da dire una cosa importante: come abbiamo accennato, nella equazione $\ref{eq:sch}$ uno dei due operatori dipende da $\Psi$. Dunque una volta calcolata la $\Psi$ in termini dei coefficienti risolvendo il sistema 2.4, si può ricalcolare l'operatore $G$ e ricominciare tutto da capo. Si arriva in tal modo ad una altra $\Psi$ piú precisa, che permette di ricalcolare l'operatore G ancora meglio, fino a quando né $\Psi$ né $G$ cambiano piú. Questa procedure iterativa di risoluzione si chiama ``self-consistente''. Al primo step di questa procedura non si conosce $\Psi$ e dunque neppure $G$. All'inizio si utilizza dunque un $G$ qualunque approssimato. Piú vicina è l'approssimazione di G iniziale al G vero, meno step di iterazione saranno necessari nella procedura self-consistente. Tutte queste operazioni, soluzione del sistema, calcolo iterativo di G, soluzione del nuovo sistema, ricalcolo del G aggiornato e cosí via, sono oggi tutte eseguite su di un computer, possibilmente piuttosto potente. Esperti programmatori (chimici, fisici, matematici, informatici) scrivono i ``codici'', cioè la sequenza di istruzioni che i calcolatori devono eseguire per pervenire al risultato finale (cioè i coefficienti finali, e l'energia ed eventualmente alcune proprietà molecolari che possono essere calcolate conoscendo la distribuzione elettronica, come il dipolo il quadrupolo etc). Come biotecnologi voi, molto probabilmente, potreste un giorno utilizzare codici ab initio all'interno di suite commerciali, per assistenza nella progettazione di farmaci, per conoscere o predire le proprietà molecolari (come il dipolo, la polarizzabilitá, l'idrofobia, la reattività etc.) di un certo farmaco/molecola, per docking, per analisi statistiche o conformazionali etc. A quel punto si porrà il problema di quale metodo utilizzare nella giungla sterminata delle metodologie ab initio. Questa ``giungla'' di approssimazioni può essere razionalizzata introducendo due principali categorie o livelli di approssimazione:

• livello 1. L'equazione 2.1 è un modo simbolico di scrivere l'equazione di Schroedinger. In realtà al posto di H e G si possono mettere espressioni matematiche diverse a seconda di quale approssimazione abbiamo scelto di adottare. Quindi la forma funzionale della Eq. 2.1 determina un primo livello di approssimazione. Per esempio se utilizziamo il semplice metodo Hartree-Fock (in base al quale l'elettrone ``sente'' gli altri elettroni in maniera ``media'') H e G avranno una certa forma. Se utilizziamo il metodo della interazione di configurazione (CI) o quello del funzionale densità (DFT), gli operatori H e G avranno altre forme analitiche, piú complicate. Tutte queste metodologie sostanzialmente approssimano l'equazione vera di Schroedinger in modo tale che si possa arrivare in tempi ragionevoli - tramite l'esecuzione di un codice su un computer - alla soluzione. Invariabilmente i metodi di approssimazione migliore (per esempio l'interazione di configurazione) sono di gran lunga i piú costosi da un punto di vista computazionale. Il metodo da utilizzare in questa pletora di approcci dipende dalle proprietà che occorre calcolare. Se per esempio vogliamo calcolare il momento di dipolo di una molecola, ci basterà utilizzare dei metodi ab initio ``a livello Hartree-Fock''. Se invece vogliamo valutare la stereo-interazione tra due molecole neutre, allora dovremo utilizzare un metodo ``correlato'', nel quale le correlazioni (cioè le interazioni istantanee) tra elettroni siano incluse esplicitamente, come per esempio nel caso della CI (interazione di configurazione), del CCT (Coupled Cluster triple excitation), del MP2(4) (Moeller Plesset perturbation theory) o DFT (Density functional theory). Tutte queste complicate terminologie si riferiscono essenzialmente al modo di rappresentare gli operatori H e G nella Eq. 2.1 in modo da approssimare a vari livelli la vera equazione di Schroedinger polielettronica.

• livello 2. Il secondo livello di approssimazione corrisponde semplicemente alla grandezza del basis set, ovvero al numero di funzioni atomiche di base nella definizione degli orbitali molecolari $\psi$ (vedi Eq. 2.2). Piú grande è il basis set e migliore sarà l'approssimazione. In linea teorica, fissato il livello di approssimazione 1, il risultato migliore per quel livello si ottiene con un basis set composto da un numero infinito di funzioni $\phi$ indipendenti.^2.3 Dovremmo peró attendere infinitamente a lungo anche disponendo di un computer potentissimo. La nomenclatura con cui vengono designate le basi è molto complicata. Si utilizzano incomprensibili acronimi del tipo 6-31G(d) o 3-21G o STO 3G., sigle che si riferiscono sia al tipo di funzioni atomiche $\phi$ utilizzate, sia al numero delle funzioni di base utilizzate. Cerchiamo di chiarire alcuni concetti principali circa questa complessa nomenclatura. In quasi tutti i calcoli ab initio, le funzioni che rappresentano orbitali atomici $\phi$ sono espresse da combinazioni di funzioni Gaussiane (del tipo $\exp((r-r_0)^{2}/\sigma^2)f(x,y,z))$) con diversa larghezza $\sigma$ opportunamente scelta far riprodurre la forma analitica degli orbitali di Slater idrogenoidi che è del tipo $\exp(-a(r-r0)) f(x,y,x)$.^2.4 Queste combinazioni si chiamano contrazioni. Quando gli orbitali atomici sono rappresentati da due o tre contrazioni diverse si parla di base ``doppia zeta'' o ``tripla zeta''. Gli elementi delle contrazioni - le singole Gaussiane - si chiamano primitive. Il numero di primitive utilizzate per una singola contrazione cambia a seconda che gli orbitali atomici siano interni o esterni (o di valenza). Questo tipo di basi, di gran lunga le piú usate nei calcoli, si chiamano basi di Pople o ``split valence''. La loro notazione generale è del tipo n-ij[k]G. G sta per ``Gaussian'' e si riferisce naturalmente alla forma funzionale della primitiva. La base 4-31G sta ad indicare una singola zeta di quattro primitive per gli orbitali $\phi$ interni ed una doppia zeta con due distinte contrazioni di tre ed una primitiva, rispettivamente per gli orbitali di ``valenza''. Sia la 6-31G che la 4-31G sono basi minime. Esse includono solo gli orbitali atomici piú interni (tipicamente gli s ed i p per gli atomi delle prime tre righe del Sistema Periodico) che sono quelli che danno il contributo maggiore all'orbitale molecolare di stato fondamentale. Queste due basi, la 6-31G e la 4-31G differiscono solo nel numero delle primitive. Piú primitive ci sono piú al lungo durerà il calcolo. Dunque la 6-31G è una base piú costosa della 4-31G, poiché include 6 primitive per gli orbitali piú interni contro le quattro della base 4-31G. Per andare oltre la base minima, occorre includere orbitali atomici che non sono occupati negli atomi isolati (per esempio gli orbitali di tipo d dell'ossigeno e dell'azoto o gli orbitali di tipo p dell'idrogeno). Questi orbitali extra si indicano a sinistra o con delle stellette o esplicitandoli: per esempio 6-31G** = 6-31G(d,p).

  Come per il livello 1, relativo alla forma funzionale degli operatori H e G, anche la scelta delle basi dipende dalle proprietà che vogliamo calcolare. Se vogliamo calcolare uno spettro Raman (cioè delle polarizzabilitá elettroniche) servirá una base grande. Per il dipolo basterà invece una base piccola come per esempio la 4-31G o la 6-31G.

Tipicamente nella soluzione del sistema 2.4 ottenuto ponendo a zero le derivate rispetto ai coefficienti, qualunque sia la forma funzionale scelta per H o G, si devono calcolare un grande numero di ``integrali''. Il numero di questi integrali da risolvere è proporzionale alla quarta potenza del numero degli elettroni, poichè in questi integrali ``multicentrici bielettronici '' sono coinvolti al massimo quattro funzioni orbitaliche $\phi$. Ogni orbitale atomico, come detto sopra, corrisponde ad una somma di primitive e dunque il calcolo è anche proporzionale a $p^4$, dove $p$ è il numero di primitive. Se io raddoppio il numero di primitive (cioè raddoppio la base), il mio calcolo diventa circa sedici volte piú pesante.

Questi che ho discusso nel presente capitolo sono i concetti che stanno dietro ad una operazione in realtà semplicissima: premere un bottone. Questo piccolo sforzo fisico consente infatti, dopo aver compilato un semplicissimo input, di lanciare un programma per poi aspettare fiduciosi che il programma finisca. Nella tabella che segue è mostrato l'input relativo al calcolo dell'eenergia e della funzione d'onda dello stato fondamentale dell'Uracile del piú noto e diffuso dei pacchetti commerciali/accademici ab initio, Gaussian:

Figure 2.3: Input al programma Gaussian per la molecola di Uracile

Questo input si riferisce al calcolo delle proprietà elettroniche dell'uracile, una delle basi del DNA, o meglio di un suo tautomero corrispondente alla forma ``lattame'' [Vedi figura 2.4]].

Figure 2.4: Tautomeri dell'Uracile.

Veniamo all'analisi del file 2.3: la prima riga indica quanta memoria il programma può utilizzare. Piú grande è questo numero piú veloce sarà il calcolo, ma piú ``affaticato'' vi sembrerà il vostro computer. Se sceglierete poca memoria, mentre il programma ``gira'', potrete utilizzare il computer per esempio per navigare in rete, ma dovrete attendere piú a lungo per ottenere i risultati. Nella seconda riga ``#SP'' sta per ``single point'', ovvero fai solo un calcolo dei coefficienti senza muovere i nuclei fino ad ottenere l'energia minima. La seconda parola della seconda riga sintetizza i due livelli di cui sopra: rb3lyp3 è un metodo del tipo DFT (density functional theory) di approssimare gli operatori H e G nella Eq. 2.1; mentre la base scelta è di tipo 6-311G(d,p), singola zeta e tripla zeta per gli orbitali atomici interni ed esterni rispettivamente. La keyword ``Polar'' indica che si vuole calcolare, oltre alla energia ed alla distribuzione elettronica [inclusi dipolo, quadrupolo etc] anche la polarizzabilitá della molecola, cioè la tendenza della molecola a modificare il proprio dipolo per effetto di un campo elettrico statico. La terza , quarta e quinta riga sono riservati a dei commenti: potete scriverci quello che volete. Sulla sesta riga sono specificati due numeri: la carica totale, zero in questo caso, e la molteplicitá di spin, 1, come nel 99% dei casi (numero di elettroni pari e tutti gli elettroni appaiati) con cui avrete a che fare. I numeri corrispondenti alle colonne 2,3,4 delle righe 7-18 corrispondono alle coordinate della molecola (in Angstroem). La prima colonna delle righe 7-19 si riferisce ovviamente al tipo di atomo. I risultati delle energie sono dati in unitá atomiche o Hartree. Nel sistema cgs, l'energia, come è noto, ha le dimensioni di carica al quadrato diviso una distanza. Nelle unita' atomiche, gli Hartree, la distanza si misura in $a_0$ = raggio di Bohr o raggio dell'atomo di idrogeno pari a 0.528 Å, cioè 0.0528 nanometri; mentre la carica si esprime in carica dell'elettrone pari a $1.60 \times 10^{-19}$ Coulomb. L'energia negativa di un Hartree corrisponde a quella che si riscontra tra un protone ed un elettrone posti alla distanza di 0.53 Å. Una mole di queste coppie di protone ed elettrone ha una energia pari a 627.5095 Kcal.