View Full Version : Ricerca scientifica con il calcolo distribuito
Calcolo distribuito
Apro questo topic con la collaborazione del moderatore di sezione Lucrezio, per far far conoscere, a chi non vi partecipa già, la ricerca scientifica tramite il calcolo distribuito e per poter discutere insieme ed approfondire la parte teorica e scientifica che sta alla base dei vari progetti.
Il calcolo distribuito utilizza le risorse di milioni di computer collegati alla rete (internet) per risolvere problemi computazionali su larga scala, in diversi ambiti della ricerca. Negli ultimi anni internet si è diffuso ed è diventato sempre più veloce, e computer sempre più potenti sono connessi alla rete. Per sfruttare al meglio tutta questa potenza spesso inutilizzata, sono nati i progetti di calcolo distribuito, che utilizzano la potenza "libera" dei pc dei volontari per portare avanti importanti ricerche che nemmeno con costosissimi supercomputer potrebbero essere eseguite. In questo modo è possibile ottenere potenze di calcolo enormi, tali da superare quelle di qualsiasi supercomputer.
Questi progetti lavorano sui PC come screensaver oppure sfruttando i "tempi morti" di utilizzo della CPU in cui sono installati (in maniera da non ridurne le prestazioni), elaborando ciascuno piccole parti di un progetto di ricerca molto vasto.
BOINC
Per sfruttare al meglio tutta questa potenza, alcuni ricercatori dell'università di Berkeley (insieme a moltissimi volontari sparsi per il globo) hanno sviluppato il progetto B.O.I.N.C. (http://boinc.berkeley.edu/) (Berkeley Open Infrastructure for Network Computing), che è una nuova e flessibile piattaforma per il calcolo distribuito che permette di sfruttare facilmente le risorse dei pc dei partecipanti ai progetti. Infatti i volontari che intendono partecipare non devono fare altro che installare il client BOINC sul loro PC, selezionare i progetti a cui vogliono partecipare e lasciare eseguire l'elaborazione. A tutto il resto penserà il client BOINC. Ovviamente dopo si possono settare molte opzioni per configurare al meglio l'elaborazione.
Ognuno può partecipare nella misura che vuole, non ci sono limiti minimi o massimi. Non si guadagna nulla, se non la soddisfazione di sapere che anche noi abbiamo contribuito, con le nostre risorse, ad un certo progetto. Il bello del calcolo distribuito è proprio quello di poter contribuire a complesse ricerche senza la necessità di nozioni tecniche, ma solamente configurando e facendo girare il programma. I team, le classifiche e le statistiche poi rendono tutto molto divertente :D
Progetti
Su BOINC, lanciato nell'estate 2004, sono oggi attivi numerosi progetti di calcolo distribuito, che spaziano in molti campi della ricerca, dalla matematica pura all'astrofisica, dalla medicina alle nanotecnologie, ecc. Alcuni progetti sono anche open source (come Seti@home, ecc), e come lo è la piattaforma BOINC stessa che è in continuo siluppo.
Ecco una lista dei principali progetti BOINC:
ClimatePrediction.net (http://climateapps2.oucs.ox.ac.uk/cpdnboinc/) - Ricerca sui cambiamenti climatici
Einstein@home (http://einstein.phys.uwm.edu/) - Ricerca delle onde gravitazionali
LHC@home (http://lhcathome.cern.ch/) - Simulazione di scontri tra particelle tramite l'utilizzo di un acceleratore (l'LHC appunto)
Rosetta@home (http://boinc.bakerlab.org/rosetta/) - Ricerca sulle proteine per la cura di molte malattie degenerative
SETI@home (http://setiathome.ssl.berkeley.edu/) - Ricerca di segnali radio di origine extraterrestre
SIMAP (http://boinc.bio.wzw.tum.de/boincsimap/) - Catalogazione di similitudini sulle proteine
Spinhenge@home (http://spin.fh-bielefeld.de/) - Studio per la realizzazione di Magneti Molecolari (Magnetismo su Nanoscala Controllato)
Quantum Monte Carlo@home (http://qah.uni-muenster.de/) - Studio della struttura e delle reazioni delle molecole tramite l'utilizzo della chimica quantistica
Proteins@home (http://biology.polytechnique.fr/proteinsathome) - Studio della struttura delle proteine
uFluids (http://www.ufluids.net/) - studio del comportamento e della stabilita' dei fluidi in condizioni di microgravità e microfluidità
World Community Grid (http://www.worldcommunitygrid.org/) - Diversi progetti di ricerca medica per la cura dell'AIDS, per lo studio del genoma umano e la cura di molte altre malattie oggi incurabili
e molti altri ancora...
Come fare per partecipare
Dopo aver installato il client BOINC (scaricabile da qui: http://boinc.berkeley.edu/download.php) dovete semplicemente eseguire la procedura di aggancio ad un progetto per ogni progetto al quale volete partecipare. Al termine della procedura di attach BOINC provvede a comunicare con il server del progetto, e a scaricare i dati necessari per avviare l'elaborazione. Il tutto è visualizzabile dalle varie finestre del programma (Project, Work, Transfer, Message, ecc). Si può partecipare sia da soli o condividendo il proprio account con altri utenti (formando quindi un miniteam), e si possono utilizzare più computer in un singolo account per aumentare la potenza elaborativa (basta fare l'attach con gli stessi dati selezionando "utente esistente").
BOINC: divertimento oltre la ricerca
Partecipare ai progetti di calcolo distribuito significa anche divertimento. Infatti oltre a conoscere tanti altri amici con interessi simili, ci si diverte molto guardando le statistiche e le classifiche generate dai vari progetti. Infatti la quantità di dati elaborata viene conteggiata con dei punti (crediti), e per i progetti BOINC sono disponibili decine di siti di statistiche mondiali molto dettagliati (es: http://www.boincstats.com http://www.boincsynergy.com, ecc).
BOINC: un valido test per il sistema
L'elaborazione con i progetti BOINC è anche un ottimo strumento per verificare la stabilità di un PC, ed è sfruttato da molti proprio per questo motivo. Infatti i client dei progetti eseguono complessi calcoli, utilizzando molto cpu e ram, mettendo bene sotto stress il sistema. Se un pc non crea problemi con alcuni progetti BOINC, si può definire veramente ROCK SOLID.
>> BOINC.Italy
http://img209.imageshack.us/img209/9793/boincitalylogo400x1451yt1.jpg (http://www.hwupgrade.it/forum/forumdisplay.php?f=10)
Hardware Upgrade partecipa alla ricerca con il team BOINC.Italy (nato e tutt'ora attivo proprio su questo forum), e invita tutti gli utenti a farne parte. Il team è presente nei principali progetti BOINC NO-PROFIT. Per approfondimenti sul team, sulla piattaforma BOINC e per la risoluzione di problemi vi invitiamo a visitare la sezione del forum dedicata cliccando >> QUI << (http://www.hwupgrade.it/forum/showthread.php?t=1380757).
Vi aspettiamo! Se decidete di partecipare non mancate di fare il join nel team BOINC.Italy (tramite l'account sul sito del progetto), e di presentarvi in QUESTO (http://www.hwupgrade.it/forum/showthread.php?t=1464542) apposito topic:
QUI (http://www.hwupgrade.it/forum/showpost.php?p=11882821&postcount=3) potete trovare una guida rapida per iniziare subito in pochi click ad elaborare sui progetti BOINC.
_________________
Questo topic, come scritto all'inizio, è nato dalla volontà dei membri del team di conoscere meglio i progetti di calcolo ai quali partecipiamo e la teoria che ci sta dietro, e di capire in che modo i vari dati vengono elaborati. Chi volesse può scrivere qualcosa relativo a uno o più dei progetti BOINC.
Lucrezio
05-06-2007, 08:44
I sistemi a molti corpi: Introduzione al problema della Chimica Teorica
Questa breve introduzione non si pone, ovviamente, lo scopo di descrivere a fondo il problema, utilizzando l'opportuno formalismo matematico ed analizzando nei dettagli i diversi metodi utilizzati in chimica teorica, ma cercherà di descrivere in modo del tutto qualitativo che tipo di sistemi vanno descritti da un punto di vista quantistico, quali sono le forze in gioco e quali le strategie, molto in generale, che si possono seguire.
Ci sono un po' di formulacce sparse, messe per completezza, ma che possono essere totalmente ignorate!
Descrizione quantistica di un sistema
La meccanica quantistica descrive un qualsiasi sistema servendosi di una funzione di stato, detta funzione d'onda. Essa dipende dalle coordinate di tutte le particelle, nonché dal loro spin, e - una volta nota - determina tutte le proprietà presenti e future del sistema.
Tale funzione deve soddisfare l'equazione di Schroedinger, ovvero deve valere - per un sistema stazionario:
http://operaez.net/mimetex/\mathcal{H}\Psi = E\Psi
dove
http://operaez.net/mimetex/\mathcal{H} è l'operatore Hamiltoniano del sistema, ovvero l'operatore che contiene tutte le interazioni energetiche del sistema
http://operaez.net/mimetex/\Psi è la funzione d'onda
E è l'energia del sistema
Sistemi a molti corpi
La descrizione di un qualsiasi tipo di sistema molecolare o atomico da un punto di vista chimico (e quindi, in particolare, "elettronico") è il grande problema della chimica teorica moderna. La meccanica quantistica fornisce, con le sue leggi, gli strumenti adatti: purtroppo le equazioni che descrivono i sistemi sono difficilmente risolubili in modo analitico.
Già la descrizione degli atomi pone dei problemi notevoli: se l'atomo di idrogeno, composto da un protone e da un unico elettrone, si può trattare in maniera esatta già l'atomo di elio è un problema che va affrontato con approssimazioni numeriche.
L'idrogeno, d'altra parte, è un sistema a singola particella: mettendosi nel sistema di riferimento del nucleo l'unica particella presente è il singolo elettrone, che risente del campo elettrico generato dal nucleo. L'elio, con due elettroni, è il più semplice sistema "a molti corpi": considerando il nucleo come una carica puntiforme che produce un campo elettrostatico Coulombiano e mettendosi in un sistema di riferimento con esso solidale rimangono comunque da descrivere due elettroni, i quali non solo interagiranno con la sorgente di potenziale elettrico dovuta alla carica nucleare, ma subiranno delle interazioni reciproche.
Sono proprio queste ultime che rendono impossibile la soluzione esatta del problema: se già un atomo con due elettroni è intrattabile in modo esatto, un atomo più grosso o in una molecola composta da tanti atomi (e quindi moltissimi elettroni!) rappresenta davvero una sfida!
I sistemi di interesse chimico: la descrizione degli elettroni
Gli elettroni sono particelle cariche che appartengono alla famiglia dei fermioni. Questo significa che sono dotati di momento angolare di spin (una proprietà squisitamente quantistica, che non trova nessun corrispettivo classico e che quindi va assunta come una proprietà intrinseca) semiintero (in particolare possono aver spin +1/2 o -1/2) e che la funzione d'onda[/b] che descrive il sistema deve cambiare segno se si scambiano fra di loro due elettroni, ovvero
http://operaez.net/mimetex/\Psi(1,2) = - \Psi(2,1)
Questo requisito è di particolare importanza, in quanto una funzione d'onda senza tale proprietà non è fisicamente accettabile!
Esaminiamo ora quali sono le forze in gioco:
Come tutto ciò che si muove gli elettroni possiedono un'[I]Energia Cinetica, che allo stesso modo è posseduta anche dai nuclei
Nuclei ed elettroni sono particelle cariche, dunque ci saranno interazioni elettrostatiche fra tutte le particelle.
Scendendo in maggiore dettaglio, le energie di cui tener conto sono:
Energia cinetica degli eletttroni http://operaez.net/mimetex/T_e
Energia cinetica dei nuclei http://operaez.net/mimetex/T_N
Interazione elettrostatica fra Nuclei ed elettroni http://operaez.net/mimetex/V_{N,e}
Interazione elettrostatica fra Nuclei e Nuclei http://operaez.net/mimetex/V_{N,N}
Interazione elettrostatica fra elettroni ed elettroni http://operaez.net/mimetex/V_{e,e}
In modo un po' più tecnico si può dire che l'Hamiltoniano del sistema è la somma di questi cinque termini.
Sono quindi davvero molte le interazioni di cui si deve tener conto!
NOTA:
Il problema fin qui è ancora abbastanza vago: non si capisce esattamente che cosa bisogna fare né come!
Introducendo un prima semplificazione esso però si delinea molto meglio e diventa - a livello intuitivo - già più comprensibile... abbiate fiducia!
L'approssimazione di Born-Oppenheimer
Una prima iniziale semplificazione al problema che si vuole trattare consiste nel "separare", almeno in prima approssimazione, l'equazione dei nuclei da quella degli elettroni. L'idea è che gli elettroni sono molto più leggeri dei nuclei (anche nel caso dell'ideogeno, il cui nucleo è costituito da un solo protone, il rapporto fra le masse è di 1/1836 !) dunque è possibile pensare che gli elettroni si adattino istantaneamente agli spostamenti dei nuclei. Dunque si può pensare che i nuclei siano "fermi" rispetto agli elettroni e che quindi interagiscano con essi semplicemente creando un campo elettrico come se fossero un insieme di cariche fisse. In questo modo la posizione dei nuclei diventa un parametro fissato del problema che si suddivide in due:
Fissate le posizioni dei nuclei, quale sarà la funzione d'onda che descrive gli elettroni e quale la loro energia?
Come cambia l'energia totale del sistema (somma dell'energia degli elettroni e di quella dei nuclei) cambiando le posizioni dei nuclei?
E' chiaro che l'energia degli elettroni dipenderà dalla posizione dei nuclei: i processi di ottimizzazione della geometria (e quindi quelli utilizzati per la descrizione della confromazione più stabile delle proteine e via dicendo) dunque procedono calcolando l'energia degli elettroni e quindi cambiando le coordinate dei nuclei in maniera da minimizzarla, il tutto in modo iterativo fino a raggiungere la condizione di stazionarietà.
Il problema elettronico
Chiaramente il punto più complicato è il primo: si tratta di risolvere l'equazione di Schroedinger per gli elettroni, calcolando così la funzione d'onda e l'energia del sistema.
Come già detto questo non può venir fatto esattamente: esistono diversi metodi più o meno precisi che partono da un livello molto semplificato fino ad introdurre delle notevoli complicazioni per migliorare le stime fatte. In generale i metodi si dividono in due categorie:
Metodi ab initio, ovvero metodi che partono dalla conoscenza delle costanti fondamentali della fisica il cui scopo è risolvere in modo più o meno approssimato l'equazione di Schroedinger. Solitamente sono molto gravosi in termini computazionali, ma danno degli ottimi risultati su qualsiasi sistema.
Fra i metodi ab initio recentemente si inizia ad usare sempre più spesso il metodo "quantum monte carlo" (è quello utilizzato ad esempio qmc@home e da folding@home), sul quale si possono trovare abbondanti informazioni in rete, ad esempio:
http://it.wikipedia.org/wiki/Quantistica_Monte_Carlo
Metodi semiempirici, che si avvalgono di alcuni dati sperimentali, come distanze di legame, angoli, energie di legame e così via per risolvere l'equazione di Schroedinger partendo "avvantaggiati". Pur molto meno gravosi sono ovviamente meno generali e affidabili dei metodi ab initio, in quanto si ancorano ai dati sperimentali per introdurre delle correzioni che migliorino dei modelli molto semplificati (ad esempio: considerano gli elettroni come cariche semplici e ci attaccano un fattore correttivo che permetta di riprodurre meglio i dati sperimentali).
[Il problema geometrico]
Supponiamo di avere scelto un metodo che, data una certa configurazione (ovvero assegnato il valore di tutte le coordinate di tutti i nuclei) ci restituisca l'energia elettronica del sistema. L'energia totale sarà semplicemente la somma di quella elettronica e della repulsione coulombiana dovuta al fatto che i nuclei sono tutti carichi positivamente: quest'ultimo termine si calcola agevolmente, quindi non costituisce un problema. Abbiamo quindi ottenuto l'energia come una funzione
http://operaez.net/mimetex/E(\vec{R_1},\vec{R_2},\cdots ,\vec{R_n})
di tutte le coordinate nucleari.
L'idea è di cambiare le coordinate in modo che l'energia diminuisca: per ottenere questo ci sono diversi metodi, che differiscono fra loro per l'uso o meno delle derivate dell'energia rispetto alle coordinate dei nuclei.
Per avere un'idea qualitativa di come si può procedere un buono schema è il seguente:
Si calcolano le derivate dell'energia rispetto a tutte le coordinate: in questo modo si ottiene un'informazione su come cambia l'energia cambiando singolarmente le coordinate (in particolare si capisce se l'energia aumenta o diminuisce aumentando o diminuendo il valore della coordinata in questione)
Si modificano le coordinate in modo che l'energia diminuisca
Si ottiene una nuova configurazione per la quale si ricalcola l'energia elettronica e quindi si ricomincia dal primo punto
In generale il processo in assoluto più dispendioso in termini di costo computazionale è il calcolo dell'energia elettronica.
Per fare un esempio molto semplice, riporto i tempi di calcolo con diversi modelli (in ordine crescente di accuratezza) per un sistema estemamente piccolo e semplice se confrontato con una proteina: il benzene
Metodo: Hartree - Fock (ab initio, molto approssimato)
Energia Totale: - 230.74600 u.a.
Distanza Carbonio-Carbonio: 138.6 pm
Distanza Carbonio-Idrogeno: 107.5 pm
Tempo di elaborazione: 6' 6.3''
Metodo: DFT - b3lyp (è un metodo piuttosto preciso, ma non completamente ab initio)
Energia Totale: - 232.30070 u.a.
Distanza Carbonio-Carbonio: 139.5 pm
Distanza Carbonio-Idrogeno: 108.6 pm
Tempo di elaborazione: 9' 42.6''
Metodo: mp2 (è un metodo che corregge l'Hartree Fock utilizzando la teoria delle perturbazioni. Confrontabile con DFT per accuratezza ma del tutto ab initio)
Energia Totale: - 231.54048 u.a.
Distanza Carbonio-Carbonio: 140.0 pm
Distanza Carbonio-Idrogeno: 108.7 pm
Tempo di elaborazione: 18' 46.9''
Metodo: CCSD (metodo "Coupled Cluster" con singole e doppie eccitazioni, estremamente preciso anche se molto gravoso)
Energia Totale: - 231.54604 u.a.
Distanza Carbonio-Carbonio: 140.6 pm
Distanza Carbonio-Idrogeno: 109.6 pm
Tempo di elaborazione: 1h 59' 8.1'' (!!!)
Metodo: mp4 (è un metodo che corregge l'Hartree Fock utilizzando la teoria delle perturbazioni fino al 4 ordine. Confrontabile con CCSD per accuratezza)
Energia Totale: - 231.55901 u.a.
Distanza Carbonio-Carbonio: 139.7 pm
Distanza Carbonio-Idrogeno: 108.7 pm
Tempo di elaborazione: 1h 8' 15.7''
La distanza carbonio - carbonio sperimentale è di 139 pm, quindi il metodo più costoso non sempre è quello che dà le previsioni migliori. L'esempio è riportato, in ogni caso, per mostrare come crescendo in ordine di precisione i tempi si allunghino notevolmente anche con molecole molto piccole... figuriamoci quindi con una proteina!
Da qui la necessità di suddividere il lavoro il più possibile e (per quel che ci riguarda) di dare il proprio contributo alla ricerca!
Lucrezio
05-06-2007, 08:44
Questo secondo post è dedicato a chi mastica bene la meccanica quantistica... diamo un po' di soddisfazione agli addetti ai lavori!
Il metodo di Hartree - Fock
L'idea del metodo di Hartree - Fock è di scrivere la funzione d'onda del sistema a molti corpi come un singolo determinante di slater - ovvero come quella di un sistema per definizione non correlato - e quindi di trovare il determinante che minimizzi l'energia totale.
Sia
http://operaez.net/mimetex/\Psi(x_1,\ldots,x_n) = \frac{1}{\sqrt{N!}}\sum_{P\in S_N}(-1)^P\prod_{j=1}^N\varphi_{P(j)}(x_j) = \frac{1}{N!}det|\varphi_{P(i)}(x_j)|
la funzione d'onda (un singolo determinante, appunto) del sistema a N particelle.
L'Hamiltoniano elettronico del sistema è quello coulombiano, ovvero, nelle unità naturali:
http://operaez.net/mimetex/\mathcal{H} = h_1 + \sum_{i\neq j}\frac{1}{2}\frac{1}{r_{ij}} = T + V + \sum_{i\neq j}\frac{1}{2}\frac{1}{r_{ij}} = -\sum \frac{\nabla^2_i}{2} - \sum_{i,N} + \frac{Z_N}{r_{iN}} + \sum_{i\neq j}\frac{1}{2}\frac{1}{r_{ij}}
Tale funzione è normalizzata, come è possibile verificare facilmente, dunque l'energia del sistema sarà:
http://operaez.net/mimetex/E=\langle \Psi |\mathcal{H}|\Psi \rangle.
E' facile vedere, specie utilizzando il formalismo della seconda quantizzazione, che l'elemento di matrice appena scritto, per un determinante di slater vale
http://operaez.net/mimetex/E=\sum_{j=1}^n\langle\varphi_j|h_1|\varphi_j\rangle + \frac{1}{2} \sum_{i,j}\langle \varphi_i\varphi_j||\varphi_i\varphi_j\rangle
dove
http://operaez.net/mimetex/\langle\varphi_i\varphi_j||\varphi_i\varphi_j\rangle = \langle\varphi_i\varphi_j|\frac{1}{r_{ij}}|\varphi_i\varphi_j\rangle - \langle\varphi_i\varphi_j|\frac{1}{r_{ij}}|\varphi_j\varphi_i\rangle
Per trovare il "miglior determinante" basta impostare un problema di minimo vincolato (va imposta la condizione al bordo che gli spinorbitali rimangano ortonormali) con i moltiplicatori di Lagrange, aggiungendo all'energia il termine identicamente nullo
http://operaez.net/mimetex/\sum_{i,j}\Lambda_{ij}(\langle\varphi_i|\varphi_j\rangle - \delta_{ij})
La procedura è lecita grazie al teorema variazionale, che garantisce che l'energia ottenuta con una procedura - appunto - variazionale non può essere inferiore a quella esatta.
In particolare la variazione prima dell'energia cambiando gli spinorbitali (detta anche "derivata funzionale") sarà:
http://operaez.net/mimetex/\delta E[\delta \varphi_k] =\langle\delta \varphi_k|h_1|\varphi_k\rangle + \frac{1}{2} \sum_{j}\{\langle \delta\varphi_k\varphi_j||\varphi_k\varphi_j\rangle + \langle \varphi_j\delta \varphi_k||\varphi_j\varphi_k\rangle\} - \sum_j \Lambda_{jk}\langle\varphi_k|\varphi_j\rangle +cc
ovvero, raccogliendo i due integrali bielettronici, che sono uguali per simmetria:
http://operaez.net/mimetex/\delta E[\delta \varphi_k] =\langle\delta \varphi_k|h_1|\varphi_k\rangle + \sum_{j}\langle \delta\varphi_k\varphi_j||\varphi_k\varphi_j\rangle - \sum_j \Lambda_{jk}\langle\varphi_k|\varphi_j\rangle +cc
dove con cc si intende "complesso coniugato".
In un punto di minimo tale variazione sarà per definizione nulla; poiché la variazione del signolo spin orbitale è del tutto arbitraria dovrà essere - scrivendo esplicitamente gli integrali:
http://operaez.net/mimetex/h_1\varphi_k(\vec{x})+ \sum_{j}\int d^3\xi\varphi^{\dagger}_j(\vec{\xi})\frac{1}{|\vec{x}-\vec{\xi}|}\varphi_j(\vec{\xi})\varphi_k(\vec{x})- \sum_j \int d^3\xi\varphi^{\dagger}_j(\vec{\xi})\frac{1}{|\vec{x}-\vec{\xi}|}\varphi_k(\vec{\xi})\varphi_j(\vec{x}) = \sum_j \Lambda_{jk}\varphi_j(\vec{x})
Definendo tre nuovi operatori:
http://operaez.net/mimetex/J(\vec{x}) = \sum_{j}\int d^3\xi\varphi^{\dagger}_j(\vec{\xi})\frac{1}{|\vec{x}-\vec{\xi}|}\varphi_j(\vec{\xi})
(detto operatore Coulombiano o di Hartree)
http://operaez.net/mimetex/K(\vec{x},\vec{\xi}) = \sum_{j}\int d^3\xi\varphi^{\dagger}_j(\vec{\xi})\frac{1}{|\vec{x}-\vec{\xi}|}\varphi_j(\vec{x})
(detto operatore di scambio)
Ed infine:
http://operaez.net/mimetex/\mathcal{F} = h_1 + J - K
detto operatore di Fock
Si ottiene qualcosa di molto vicino alle equazioni di Hartree - Fock:
http://operaez.net/mimetex/\mathcal{F}|\varphi_k\rangle = \sum_j \Lambda_{jk}\varphi_j
Per ottenere le equazioni di Hartree - Fock basta notare che la matrice dei moltiplicatori di Lagrange é sicuramente autoaggiunta e dunque - per il teorema spettrale - diagonalizzabile in una base ortonormale di autovettori.
Facendo il cambio di base, ovvero passando a quelli che vengono detti orbitali canonici di Hartree - Fock si ottengono finalmente le equazioni di Hartree - Fock:
http://operaez.net/mimetex/\mathcal{F}|\phi_k\rangle = \epsilon_k|\phi_k\rangle
Ora, l'energia totale del sistema, in questa approssimazione, sarà data dal valore di aspettazione dell'Hamiltoniano fatto sul determinante ottenuto con gli orbitali che risolvono le equazioni di Hartree - Fock e dunque:
http://operaez.net/mimetex/E=\langle \Psi |\mathcal{H}|\Psi \rangle = \langle \Psi |\mathcal{F} + W - (J - K)|\Psi \rangle
L'energia può essere scritta in funzione delle energie orbitaliche e degli elementi di matrice dell'operatore monoelettronico, ad esempio come:
http://operaez.net/mimetex/E=\frac{1}{2}\sum_j(h_{1,jj} + \epsilon_j).
Si noti che negli operatori J e K compaiono gli spinorbitali, ovvero le soluzioni dell'equazione!
Dunque le equazioni non sono lineari e vanno risolte con un metodo iterativo: si parte da un set di orbitali di prova - con i quali si ottiene l'operatore di Fock - e si itera finché il potenziale non è stazionario (nel senso numerico, ovviamente, del termine). Il campo così ottenuto viene detto autoconsistente, o, nel più tecnico inglese, "Self Consistent Field".
Si noti infine che il metodo di Hartree Fock tiene conto in modo esatto dell'energia di scambio, ma trascura del tutto la correlazione, motivo per ricorrere a metodi più sofisticati come quelli utilizzati nel post precedente per ottenere energia e geometria del benzene!
Posso postare ? :D
domanda : voglio iscrivermi a worldgrid , posso scegliere se scaricare il software united device o boinc. secondo voi che devo fare ? perchè mi dice che uno di ste 4 ricerche va solo con united device ...
Coapzpaza
Posso postare ? :D
domanda : voglio iscrivermi a worldgrid , posso scegliere se scaricare il software united device o boinc. secondo voi che devo fare ? perchè mi dice che uno di ste 4 ricerche va solo con united device ...
Coapzpaza
Il progetto per la distrofia muscolare dispone del client per UD e non per BOINC tutti gli altri hanno un client anche per BOINC.
Comunque per partecipare scaricati e installa BOINC poi fai l'attach ai progetti che ti interessano...
Il progetto per la distrofia muscolare dispone del client per UD e non per BOINC tutti gli altri hanno un client anche per BOINC.
Comunque per partecipare scaricati e installa BOINC poi fai l'attach ai progetti che ti interessano...
uhm allora vado con boinc ...
Coapzpaza
Per la ricerca uso spesso il metodo di Nelder-Mead (detto dell'Ameba o del simplesso, credo) per la minimizzazione delle funzioni di più variabili... E' stato testato?
Lucrezio
05-06-2007, 16:52
In arrivo un'ultima tornata di dati (sto finendo la simulazione, sempre che non mi scoppi il pc come all'ultimo tentativo :muro: )!
In seguito qualcosina su Hartree Fock mi sa che ve la puppate!
:asd:
stbarlet
05-06-2007, 23:45
ci voleva un po di QM :D
In arrivo un'ultima tornata di dati (sto finendo la simulazione, sempre che non mi scoppi il pc come all'ultimo tentativo :muro: )!
In seguito qualcosina su Hartree Fock mi sa che ve la puppate!
:asd:
Ti è scoppiato il PC? :D
Lucrezio
07-06-2007, 19:01
Ti è scoppiato il PC? :D
Ho dovuto usare una base piccina perché la mia partizione linux era troppo piccola (a 35GB di scratch il programma si è piantato...)
:stordita:
Comunque ho fatto l'ultima simulazione (mp4) e l'ho riportata ;)
Lucrezio
08-06-2007, 15:15
UP per Lucrezio :D
Datemi tempo :cry:
Sono sotto esami!
Non vorrei sbagliarmi ma il quarto termine che corregge l'energia nel metodo delle perturbazioni viene positivo in alcuni casi (:rolleyes:) e quindi poi l'energia totale ti aumenta. Cioè in soldoni utilizzando anche un metodo + spinto (quarto ordine) ottieni un'energia maggiore rispetto a quella con un metodo che si ferma alla correzione del terzo ordine. Quindi invece di migliorare il valore dell'energia totale... ti allontani dal valor vero. Ma su che molecola lo stai provando?
Prendiamo x esempio i tuoi calcoli sul benzene:
I sistemi a molti corpi: Introduzione al problema della Chimica Teorica
Il benzene:
Metodo: Hartree - Fock (ab initio, molto approssimato)
Energia Totale: - 230.74600 u.a.
Distanza Carbonio-Carbonio: 138.6 pm
Distanza Carbonio-Idrogeno: 107.5 pm
Tempo di elaborazione: 6' 6.3''
Metodo: DFT - b3lyp (è un metodo piuttosto preciso, ma non completamente ab initio)
Energia Totale: - 232.30070 u.a.
Distanza Carbonio-Carbonio: 139.5 pm
Distanza Carbonio-Idrogeno: 108.6 pm
Tempo di elaborazione: 9' 42.6''
Metodo: mp2 (è un metodo che corregge l'Hartree Fock utilizzando la teoria delle perturbazioni. Confrontabile con DFT per accuratezza ma del tutto ab initio)
Energia Totale: - 231.54048 u.a.
Distanza Carbonio-Carbonio: 140.0 pm
Distanza Carbonio-Idrogeno: 108.7 pm
Tempo di elaborazione: 18' 46.9''
Metodo: CCSD (metodo "Coupled Cluster" con singole e doppie eccitazioni, estremamente preciso anche se molto gravoso)
Energia Totale: - 231.54604 u.a.
Distanza Carbonio-Carbonio: 140.6 pm
Distanza Carbonio-Idrogeno: 109.6 pm
Tempo di elaborazione: 1h 59' 8.1'' (!!!)
Metodo: mp4 (è un metodo che corregge l'Hartree Fock utilizzando la teoria delle perturbazioni fino al 4 ordine. Confrontabile con CCSD per accuratezza)
Energia Totale: - 231.34575 u.a.
Distanza Carbonio-Carbonio: 139.7 pm
Distanza Carbonio-Idrogeno: 108.7 pm
Tempo di elaborazione: 1h 8' 15.7''
Se ci concentriamo sui valori di energia totale della molecola di benzene vediamo come il metodo mp2 sia molto più in linea col risultato del CCSD (molto fine) di quanto non sia un metodo mp4. Sarebbe interessante quindi analizzare per bene i risultati che Gaussian (o altro programma) ti ha dato per individuare le varie correzioni per i vari ordini di perturbazione. Ovviamente tutto il discorso vale a parità di set di base preso per descrivere gli orbitali.
Un altra cosa che si potrebbe fare... vedere i risultati per un mp3. Io direi che un mp3 sicuramente migliorerà i valori di energia rispetto al mp2, senza l'inconveniente della correzione (positiva!) dell'energia che può dare un mp4.
Lucrezio
10-06-2007, 12:12
Prendiamo x esempio i tuoi calcoli sul benzene:
Se ci concentriamo sui valori di energia totale della molecola di benzene vediamo come il metodo mp2 sia molto più in linea col risultato del CCSD (molto fine) di quanto non sia un metodo mp4. Sarebbe interessante quindi analizzare per bene i risultati che Gaussian (o altro programma) ti ha dato per individuare le varie correzioni per i vari ordini di perturbazione. Ovviamente tutto il discorso vale a parità di set di base preso per descrivere gli orbitali.
Un altra cosa che si potrebbe fare... vedere i risultati per un mp3. Io direi che un mp3 sicuramente migliorerà i valori di energia rispetto al mp2, senza l'inconveniente della correzione (positiva!) dell'energia che può dare un mp4.
Ciao!
Ho effettivamente usato Gaussian, ma purtroppo non sono riuscito a fare i conti con basi equivalenti. Ho usato 6-311+G(d) tranne che per mp4 (non mi ha fatto usare le diffuse) e per CCSD (non ho abbastanza spazio su disco per gli scratch generati dal conto con una base TZ e ho dovuto usare quella doppia...)
I risultati quindi non sono direttamente confrontabili, per ovvi motivi, ma ci tenevo semplicemente a dare un´idea dei costi a livello di tempo macchina andando ad utilizzare metodi via via piu' accurati (ho recentemente fatto un seminario sul CC e considerando solo le doppie eccitazioni, per i 6 elettroni pi-greco del benzene con base minima (12 SO) vengono 594 equazioni con un numero che non mi ricordo compreso fra 26000 e 27000 termini l´una... e i termini sono integrali bielettronici! :eek: ) senza introdurre un confronto diretto ¨ad armi pari¨! Alla fine il parametro quanto si avvicina alla distanza sperimentale C-C rispetto a quanto tempo ci mette e'qualitativamente indicativo e spiega perche'andando a livelli di precisione molto alti, che sono quelli che ci si prefigge di raggiungere nei progetti di calcolo distribuito, ci sia bisogno della collaborazione di tutti :D
Per mp4 verifico, ho ancora tutti i log!
Ma sei un computazionale?
EDIT! Avevo sbagliato a riportare l´energia, comunque qui se vuoi c´è la fine del log:
1\1\GINC-FILIPPONB\FOpt\RMP4DQ-FC\6-311G(d)\C6H6\FILIPPO\05-Jun-2007\0
\\#p opt mp4(DQ)/6-311g(d)\\Geometria dft del benzene\\0,1\C,-1.651665
6405,-1.571753493,0.0000622823\C,-0.2541961142,-1.5718231686,0.0005047
613\C,0.4445826826,-0.3616424139,-0.0000457729\C,-0.2540915099,0.84863
68538,-0.0010466366\C,-1.6515281778,0.8487065215,-0.0014861017\C,-2.35
03233726,-0.3615026487,-0.0009418524\H,-2.1955392333,-2.5136962982,0.0
004869318\H,0.2895820054,-2.5138209969,0.0012845076\H,1.5322659483,-0.
3616770827,0.0002934028\H,0.2898154757,1.7905603852,-0.0014730229\H,-2
.1953394217,1.7906852299,-0.0022656724\H,-3.438006642,-0.3614298884,-0
.001284827\\Version=AM64L-GDVRevE.01\State=1-A\HF=-230.741645\MP2=-231
.5342744\MP3=-231.5604688\MP4D=-231.5815274\MP4DQ=-231.559011\RMSD=3.6
19e-09\RMSF=2.688e-05\Dipole=0.,0.,-0.0000015\PG=C01 [X(C6H6)]\\@
The archive entry for this job was punched.
Beh non sono un computazionale, anche perchè sto alla quinquennale però ho fatto l'esame di chimica computazionale e qualcosina l'ho vista :rolleyes:. Peccato che non sei riuscito a usare la stessa base, sarebbe stato interessante per confrontare i vari metodi. Il caso del benzene inizia ad essere bello tosto se si inizia ad affinare il metodo e ad incrementare la "grandezza" del set di base. Per i metodi di perturbazione infatti bisogna essere un pò cauti e valutare bene le cose soprattutto quelli al di sopra del terzo ordine. Sarebbe interessante vedere la mole dei calcoli per un piccolo peptide, a dire la verità non ci ho mai provato. Penso che un HF con una base minima si potrebbe anche tentare :sofico:. Se trovo del materiale te lo mando. :D
Lucrezio
12-06-2007, 04:05
Beh non sono un computazionale, anche perchè sto alla quinquennale però ho fatto l'esame di chimica computazionale e qualcosina l'ho vista :rolleyes:. Peccato che non sei riuscito a usare la stessa base, sarebbe stato interessante per confrontare i vari metodi. Il caso del benzene inizia ad essere bello tosto se si inizia ad affinare il metodo e ad incrementare la "grandezza" del set di base. Per i metodi di perturbazione infatti bisogna essere un pò cauti e valutare bene le cose soprattutto quelli al di sopra del terzo ordine. Sarebbe interessante vedere la mole dei calcoli per un piccolo peptide, a dire la verità non ci ho mai provato. Penso che un HF con una base minima si potrebbe anche tentare :sofico:. Se trovo del materiale te lo mando. :D
Oddio... con cinque o sei amminoacidi secondo me anche a livello dft dovrebbe venire un buon risultato!
Se guardi quello che fanno in qmc@home in fondo riescono ad andare a livelli di precisione molto buona con tre amminoacidi (e quantum montecarlo è davvero molto preciso, anche per le interazioni deboli che richiedono almeno l'mp2 e una base con funzioni diffuse!)
Comunque appena torno a casa (per ora i conti li ho fatti con il portatile!) provo il CCSD con TZ!
Vediamo che succede...
Ma sei al vecchio ordinamento? Che indirizzo hai preso?
Io sto finendo la triennale, poi pensavo di fare chimica fisica...
Si QMC riesce a fare davvero calcoli molto buoni (sia per metodo usato sia per set di base)... ma con tutti quei pc a lavorare vanno come missili :sofico:. Un DFT si potrebbe pure provare per un peptide anche se penso che senza una Z-matrix buona di partenza ci vuole un casino di tempo x ottimizzare l'energia. Io sto alla quinquennale era l'ultimo anno del vecchio ordinamento prima dell'avvento della triennale, indirizzo inorganico a Roma. La specialistica chimico-fisica è molto interessante, complimenti x la scelta.
Lucrezio
07-07-2007, 16:54
C'è voluto un po' di tempo ma ho mantenuto la mia promessa :D
...
sarebbe questo il famigerato messaggio specifico?
non ho capito un paio di crocette sopra le variabili, sono morte? :D
Lucrezio
14-07-2007, 12:58
sarebbe questo il famigerato messaggio specifico?
non ho capito un paio di crocette sopra le variabili, sono morte? :D
:Prrr:
Ho usato il dagger al posto dell'asterisco per dire complesso coniugato... il latex su forum ha qualche problema tecnico con queste cose...
C'è voluto un po' di tempo ma ho mantenuto la mia promessa :D
Ma hai finito di far frullare il pc con Gaussian? :read: che molecole hai testato?
Lucrezio
16-07-2007, 01:30
Ma hai finito di far frullare il pc con Gaussian? :read: che molecole hai testato?
Ci ho fatto la tesina con gaussian, ho fatto girare un b3lyp su un cluster di 10 molecole di pentacene simulando un intorno anisotropo con pcm :sofico:
Altro che benchmark :asd:
Ci ha messo 12 ore :eek:
Ci ho fatto la tesina con gaussian, ho fatto girare un b3lyp su un cluster di 10 molecole di pentacene simulando un intorno anisotropo con pcm :sofico:
Altro che benchmark :asd:
Ci ha messo 12 ore :eek:
Ma il modello pcm è implementato su Gaussian, perchè io non ho mai fatto calcoli con il pcm ma con il scrf (Self Consistent Reaction Field o modello di Onsager, meno sofisticato del tuo :rolleyes:) e ho sempre usato Gamess. Cmq non ci hai messo neanche tanto come tempo, diciamo il giusto. Certo qualche risultato potevi pure farcelo vedere ;), che solvente hai utilizzato?
Lucrezio
16-07-2007, 18:11
Ma il modello pcm è implementato su Gaussian, perchè io non ho mai fatto calcoli con il pcm ma con il scrf (Self Consistent Reaction Field o modello di Onsager, meno sofisticato del tuo :rolleyes:) e ho sempre usato Gamess. Cmq non ci hai messo neanche tanto come tempo, diciamo il giusto. Certo qualche risultato potevi pure farcelo vedere ;), che solvente hai utilizzato?
Self made, perché doveva simulare un reticolo cristallino.
Abbiamo impostato l'eptano per via delle caratteristiche chimico - fisiche, poi abbiamo cambiato il tensore di inerzia.
Comunque PCM è implementato in Gaussian dalla versione G03, inoltre l'implementazione (e il modello) sono stati fatti in gran parte dalla mia relatrice, quindi mi aspettavo che me lo facesse usare :D
Ah, i calcoli erano su una macchina biprocessore con 4 giga di ram e dischi scsi in raid... quindi 12 ore sono effettivamente tantine :D
Self made, perché doveva simulare un reticolo cristallino.
Abbiamo impostato l'eptano per via delle caratteristiche chimico - fisiche, poi abbiamo cambiato il tensore di inerzia.
Comunque PCM è implementato in Gaussian dalla versione G03, inoltre l'implementazione (e il modello) sono stati fatti in gran parte dalla mia relatrice, quindi mi aspettavo che me lo facesse usare :D
Ah, i calcoli erano su una macchina biprocessore con 4 giga di ram e dischi scsi in raid... quindi 12 ore sono effettivamente tantine :D
Capito... interessante! Ma ha scopi pratici questo studio o era puramente teorico? Vedo cmq che sei un esperto in materia :D, se ho qualche intoppo "computazionale" ora so chi tartassare :p. Certo 12 ore per quella macchina non sono poi pochine, io pensavo che lo avessi fatto girare su un computer + "casalingo" ;)
Ciao a tutti. Scusate l'intromissione...
Io faccio elaborare il mio pc su QMC ormai da parecchio tempo e sono convinto della sua utilità, perlomeno dal punto di vista della ricerca di base.
Volevo chiedere a voi esperti se i risultati ottenuti da QMC vengono in qualche modo utilizzati da qualcuno in qualche ricerca o cosa simile.
Cioè... mi piacerebbe capire in che modo i calcoli che esegue il mio computer diventano poi utili all'atto pratico.
Lucrezio
20-09-2007, 21:20
Ciao a tutti. Scusate l'intromissione...
Io faccio elaborare il mio pc su QMC ormai da parecchio tempo e sono convinto della sua utilità, perlomeno dal punto di vista della ricerca di base.
Volevo chiedere a voi esperti se i risultati ottenuti da QMC vengono in qualche modo utilizzati da qualcuno in qualche ricerca o cosa simile.
Cioè... mi piacerebbe capire in che modo i calcoli che esegue il mio computer diventano poi utili all'atto pratico.
Per ora non credo, ma il progetto è di fare in modo che il qmc diventi davvero utile!
Ora come ora ha dei problemi piuttosto gravi: la funzione d'onda che si ottiene deve essere sottoposta a dei vincoli esterni per essere fisicamente accettabile, il che impone dei limiti all'accuratezza... inoltre non si fanno ancora le derivate dell'energia analiticamente.
D'altra parte il metodo scala molto bene con le dimensioni del sistema e rappresenta una grande speranza!
Qmc@home fornisce dati per poter investigare meglio sulle caratteristiche del metodo stesso... ovvero: ricerca pura ;)
Interessante..
Questi problemi di cui parli sono risolvibili o sono proprio teoricamente impossibili da inserire nel metodo?
Lucrezio
21-09-2007, 16:37
Interessante..
Questi problemi di cui parli sono risolvibili o sono proprio teoricamente impossibili da inserire nel metodo?
Non me ne intendo abbastanza per risponderti, l'idea è che non si è ancora riusciti a far produrre a qmc una funzione d'onda antisimmetrica se non imponendo dei vincoli (ovvero i nodi della funzione)... e questo provoca dei casini.
Spero di affrontare un po' meglio questo argomento l'anno prossimo :D
Stando ai dati riportati dal noto sito di statistiche BOINCStats, BOINC ha raggiunto e superato la soglia del PetaFLOP di potenza elaborativa.
Nella giornata di ieri, 30 gennaio, sono stati assegnati in totale ben 106 millioni di crediti, equivanenti a 1.06 PetaFLOP/sec. Dalla news sulla homepage del progetto BOINC arrivano i complimenti a tutti i progetti ed a tutti i partecipanti per il traguardo raggiunto.
Il "flop" in informatica è l'unità di misura della capacità di calcolo, ed quivale ad un'operazione in virgola mobile (in inglese Floating Point). Tali opearazioni, abbastanza complesse, sono usate in tutte le simulazioni scientifiche. "Peta" è un prefisso che significa 10^15, cioè 1.000.000.000.000.000. Quindi un computer da un petaflop al secondo può effettuare 1.000.000.000.000.000 operazioni in virgola mobile al secondo. O, per renderla semplice, 1 petaflop = un milione di miliardi di istruzioni/calcoli al secondo.
E' un traguardo notevole ed impressionante, se confrontato con quello del più potente supercomputer attualmente in funzione, l'IBM Blue Gene/L che ha una potenza di 478 TeraFLOP al secondo (circa mezzo PetaFlop), e soprattutto se si pensa che tutta quella potenza è offerta da normali PC di volontari sparsi su tutta la terra. La rete di computer creata da BOINC permette ai ricercatori di usufruire di una grande potenza di calcolo a costo praticamente nullo e poter investire soldi nella ricerca piuttosto che nell'acquisto e nel mantenimento si un supercomputer a prezzi proibitivi (che praticamente nessuno potrebbe permettersi!).
La potenza totale dedicata ai progetti BOINC, come si può vedere sempre da BOINCStats è in costante crescita, e attualmente la media è di circa 900milioni di crediti al giorno (0,9 PetaFLOP).
Giusto per curiosità: la potenza attuale del nostro team, BOINC.Italy, è circa lo 0,25% del totale della potenza BOINC. Praticamente servirebbero (solo) 400 BOINC.Italy per fare tutto quello che viene fatto attualmente nel mondo!
(autore: Ghz fonte: www.boincitaly.org)
Lucrezio
04-02-2008, 22:55
:eek:
Ah, preparatevi perché quest'anno faccio QMC in un corso quindi... post pieno di formule in arrivo :D
Lucrezio
16-06-2008, 18:15
In arrivo un post di spiegazione sul Quantum Montecarlo!
E' una cosa molto più semplice di quanto si possa immaginare :D
Se riesco posto anche un programmino in fortran per calcolare l'energia dell'atomo di elio commentandolo, in modo da chiarire bene il funzionamento della tecnica.
E' uno dei primissimo programmi che scrivo, quindi i programmatori esperti sono pregati di non sbranarmi e - possibilmente - di darmi qualche buon suggerimento ;)
se a qualcuno interessa ci sonbo un pò di calcoletti da fare qua
http://www.hwupgrade.it/forum/showthread.php?t=1838179
thanks ;)
killercode
15-10-2008, 13:43
se a qualcuno interessa ci sonbo un pò di calcoletti da fare qua
http://www.hwupgrade.it/forum/showthread.php?t=1838179
thanks ;)
In università da voi non c'è un elaboratore per queste cose?
Tipo questo da noi http://www.fis.unipr.it/dokuwiki/doku.php?id=lca:servizi:calcolo
si ma è troppo complicato accedere, bisognerebbe fare dei passaggi burocratici, dovrei installare del software, troppe spiegazioni da dare, sarebbe troppo complicato.
mi pareva la via più breve ;)
Lucrezio
19-10-2008, 19:39
Quantum Montecarlo
Introduzione
Il metodo "Quantum Montecarlo" e' uno dei metodi piu' recenti della chimica teorica e prevede un approccio completamente differente da quello "classico".
Uno dei grandi problemi della chimica teorica e' sempre stato associato alla difficolta' di fare velocemente integrali in molte dimensioni. La funzione d'onda elettronica di un sistema come una molecola, dotata di n elettroni, e' un oggetto che dipende dalle coordinate di ognuno di essi: gia' per un atomo piccolo come il Berillio (4 elettroni) ci si trova a dover lavorare con una funzione di 12 variabili!
Calcolare l'energia di un simile sistema richiederebbe, in linea di principio, di calcolare un integrale di questo genere:
http://operaez.net/mimetex/\int d^3r_1d^3r_2d^3r_3d^3r_4 \bar{\Psi}(\vec{r}_1,\vec{r}_2,\vec{r}_3,\vec{r}_4)\mathscr{H}\Psi(\vec{r}_1,\vec{r}_2,\vec{r}_3,\vec{r}_4)
dove con
http://operaez.net/mimetex/\mathscr{H}=-\sum_{j=1}^4 \frac{\nabla^2_j}{2} - \sum_{j=1}^4\frac{4}{r_j}+\frac{1}{2}\sum_{i\neq j}^4\frac{1}{r_{ij}}
si intende l'Hamiltoniano elettronico (in unita' atomiche), in questo caso particolare dell'atomo di berillio.
Supponendo di voler applicare il principio variazionale (http://en.wikipedia.org/wiki/Variational_principle) e di proporre una certa forma della funzione d'onda dipendente da alcuni parametri, per ottenere la miglior stima dell'energia bisognerebbe calcolare l'energia e minimizzarla rispetto ai parametri stessi, ripetendo l'integrale riportato un numero considerevole di volte.
Il modo "classico" per calcolare un integrale del genere (non si fa a mano, purtroppo difficilmente e' mai analitico!) consiste nel fissare una certa griglia di punti e quindi stimare l'integrale con la regola dei trapezi (http://it.wikipedia.org/wiki/Regola_del_trapezio), o di Cavalieri - Simpson (http://it.wikipedia.org/wiki/Regola_di_Cavalieri-Simpson) o altre forme di quadratura.
Per un integrale di questo calibro (12 dimensioni), supponendo di prendere anche solo 10 punti per dimensione si dovrebbero fare 10^12 valutazioni di funzione (mille miliardi!!!), impresa laboriosa anche per un calcolatore piuttosto potente... e parliamo di un solo atomo con solo quattro elettroni!!!
Se si considera poi che i primi calcoli di chimica teorica sono stati effettuati piu' di 50 anni fa... beh, non c'e' bisogno di dire che tale strada e' del tutto impercorribile!
L'approccio tradizionale a questo problema introduce la cosiddetta approssimazione orbitalica: le funzioni d'onda variazionali vengono scritte come prodotto di funzioni ognuna delle quali dipende dalle coordinate di un solo elettrone (questa funzione viene detta "orbitale", per chi vuole avventurarsi in qualcosa di un po' "tecnico" il terzo post di questo thread contiene un esempio di tecnica variazionale basata sugli orbitali). Si puo' dimostrare che questo comporta il calcolo di integrali al piu' in sei variabili; una scelta opportuna delle funzioni di base (gli orbitali, appunto!) permette di ridurre la dimensionalita' degli integrali fino a 2.
Questa strategia - e' tutt'ora in assoluto la piu' utilizzata! - presenta alcuni vantaggi notevoli:
Permette di ricondurre la soluzione di un'equazione differenziale alla soluzione di un sistema di equazioni lineari
Comporta il calcolo di integrali "facili"
E' esatta per elettroni non interagenti (il cosiddetto gas ideale di fermi)
Si presta con grande facilita' ad interpretazioni intuitive
Per chiarire il punto dolente di questa approssimazione e' necessario scendere un po' piu' nel dettaglio della struttura elettronica.
Gli elettroni di un atomo o di una molecola, infatti, si muovono in modo non indipendente fra di loro. Banalizzando un po' il discorso la correlazione elettronica e' conseguenza di due caratteristiche fisiche degli elettroni: sono fermioni (particelle dotate di spin semintero) indistinguibili e sono carichi.
E' possibile dimostrare che l'indistinguibilita' degli elettroni (non e' possibile pensare ad un esperimento che permetta di "etichettare" gli elettroni di una molecola, in modo da sapere qual e' il primo, qual e' il secondo e cosi' via) si riflette sulla funzione d'onda che li descrive tramite una proprieta' di simmetria: scambiando le coordinate di due elettroni essa deve cambiare segno e quindi essere antisimmetrica per permutazioni dispari (per i bosoni, invece, essa rimarrebbe invariata). La dimostrazione di questo teorema (cosiddetto Teorema Spin-Statistica (http://en.wikipedia.org/wiki/Spin-Statistics_theorem)) richiede di entrare nell'ambito della teoria dei campi; l'antisimmetria della funzione d'onda viene quindi solitamente presa come un postulato (principio di esclusione di Pauli) in chimica quantistica.
Una conseguenza di questo teorema e' che elettroni con spin parallelo tendono ad evitarsi fra di loro: in questo modo possono avvicinarsi di piu' al nucleo carico positivamente e dunque si ha un guadagno di energia (detto barbaramente: se uno sta da una parte e uno dall'altra si danno meno "fastidio" fra di loro che se stanno dalla stessa parte, dunque possono avvicinarsi di piu' al nucleo!). Questo contributo - favorevole - all'energia viene detto energia di scambio, ed e' relativamente facile da calcolare anche nell'ambito dell'approssimazione orbitalica in quanto e' una conseguenza dello spin che viene imposta semplicemente scegliendo funzioni d'onda antisimmetriche.
Il vero "busillis" della chimica quantistica e' rappresentato dal contributo all'energia dovuto alla correlazione istantanea fra i oti degli elettroni. Fisicamente quello che avviene e' che gli elettroni, essendo carichi, tendono ad evitarsi a prescindere dal loro stato di spin proprio grazie alla repulsione elettrostatica. Questo fenomeno non puo' essere descritto bene in termine di un modello a particelle indipendenti (come il modello orbitalico: ogni orbitale, infatti, dipende solo dalle coordinate di un elettrone!), in quanto e' proprio la fisica che gli sta dietro a gridare che si tratta di un contributo che dipende dalla posizione reciproca delle particelle e non da quella assoluta!
Di conseguenza, per riuscire a calcolare l'energia di correlazione, e' necessario ricorrere a tecniche piuttosto complicate.
Se, infatti, si scrive la funzione d'onda come un signolo prodotto antisimmetrizzato (detto anche singolo determinante, in quanto il modo piu' semplice di scrivere tale prodotto e' di vederlo come un determinante di Slater (http://it.wikipedia.org/wiki/Determinante_di_Slater)) non si puo' avere energia di correlazione, dato che gli elettroni sono trattati l'uno indipendentemente dall'altro.
La chimica teorica classica utilizza un teorema di completezza che afferma che i determinanti di slater formano una base completa di un certo spazio di funzioni al quale deve appartenere la funzione d'onda esatta: dunque questa sara' esprimibile come una combinazione lineare di determinanti (metodo Full-CI).
Detta in poche parole, due determinanti saranno meglio di uno, tre meglio di due e cosi' via. L'idea e' sicuramente buona, ma e' praticamente impercorribile: infatti il numero di determinanti che si possono costruire a partire da N funzioni di base scala piu' o meno come il fattoriale di N: gia' sistemi piccoli come l'acqua sono alla portata solo di computer estremamente potenti.
E' quindi necessario troncare questa espansione ed utilizzare un numero ridotto di funzioni, oppure adottare altre tecniche (come la teoria delle perturbazioni).
In tutti i casi si ottengono dei metodi che scalano male e che sono di difficile applicazione a sistemi dell'estensione di una molecola di importanza biologica, pur fornendo dei risultati eccezionali sui sistemi dove sono applicabili - mi occupo, per la tesi, proprio di uno di questi metodi (Coupled-Cluster), posso garantire che e' uno spettacolo, ma per fare dei conti su un sistema con un centinaio d'atomi ci vuole un supercomputer!
Una strategia diversa per calcolare gli integrali
L'idea che sta alla base del metodo montecarlo e' piuttosto semplice. Si supponga di voler calcolare il valore di pigreco: se si considera una circonferenza di raggio unitario, pigreco sara' proprio l'aera della circonferenza stessa. Si inscriva tale circonferenza in un quadrato: esso avra' area pari a 4.
Dunque il rapporto fra le aree sara' pigreco/4
Si supponga di lanciare delle freccette dentro al quadrato in modo che ogni punto abbia la stessa probabilita' di essere raggiunto: il numero di frecce all'interno del cerchio diviso il totale delle frecce dovrebbe tendere a pigreco quarti, se le frecce lanciate sono molte!
E' possibile fare un test con il seguente programmino (scrivo in fortran 77, scusate ma e' l'unico linguaggio che conosco :strodita: ):
program pigreco
implicit real*8 (a-h,o-z)
parameter (Nmax=100000000,nseed=12345)
n=0
call srand(nseed)
do j=1,Nmax
x=2*(0.5d0-rand())
y=2*(0.5d0-rand())
r=x**2+y**2
if (r.le.1) then
n=n+1
endif
enddo
pigr=4.d0*n/Nmax
write(6,*) pigr
end
Variando il parametro Nmax (numero delle frecce) si ottengono questi risultati:
Nmax=1.000 ----> pigr=3.148000
Nmax=1.0000 ----> pigr=3.115600
Nmax=1.00000 ----> pigr=3.133360
Nmax=1000.000 ----> pigr=3.141240
Nmax=10.000.000 ----> pigr=3.141419
Nmax=100.000.000 ----> pigr=3.141656
Come si vede, andando ad aumentare il numero di "frecce", il risultato migliora.
A garanzia di quanto ottenuto c'e' un teorema di analisi matematica, il teorema del limite centrale (http://it.wikipedia.org/wiki/Teorema_del_limite_centrale). A partire dallo stesso (se ci fosse qualche matematico di buona volonta' disposto a farlo sarebbe molto apprezzato :D) e' possibile dimostrare che i risultati ottenuti con il metodo montecarlo convergono, inoltre si puo' stimarne l'errore in quanto la varianza va come 1/N.
In un modo un po' piu' formale si puo' dire che, andando a generare dei punti distribuiti in maniera uniforme, si va a campionare l'area del quadrato in modo tanto migliore quanto piu' sono i punti generati: dunque, disponendo di un generatore di numeri casuali come rand() del fortran (che genera numeri distribuiti in modo uniforme fra 0 e 1) e' possibile calcolare pigreco in modo approssimato.
Lo stesso metodo puo' essere utilizzato per integrare una funzione. Si supponga ad esempio di voler calcolare l'integrale
http://operaez.net/mimetex/I=\int_a^bf(x)dx
La "ricetta" del metodo montecarlo prescrive di procurarsi una distribuzione uniforme definita sull'intervallo [a,b] (usero' la lettera greca "mu" per le distribuzioni), che deve essere non negativa e normalizzata ad uno:
http://operaez.net/mimetex/\int d\mu = 1
A questo punto e' possibile scrivere:
http://operaez.net/mimetex/\int_a^bf(x)dx = \int f(x) d\mu \simeq \sum f(x_i)
dove i punti x_i sono distribuiti secondo mu.
D'altra parte mu e' la distribuzione uniforme, ovvero quella che afferma che tutti i punti appartenenti all'intervallo [a,b] hanno la stessa probabilita' di essere campionati: una buona funzione di distribuzione potrebbe essere la funzione caratteristica di [a,b]:
http://operaez.net/mimetex/\mu(x) = \left \{ \begin{array} 1 \quad x\ \in\ [a,b]\\ 0 \quad\mathrm{altrimenti}\end{array}\right .
Tale funzione, pero', non e' normalizzata:
http://operaez.net/mimetex/\int \mu(x)dx = b-a dunque va divisa per (b-a). In sintesi l'integrale sara':
http://operaez.net/mimetex/\int_a^bf(x)dx = \int f(x) d\mu \simeq \frac{1}{b-a} \sum f(x_i)
con x_i distribuiti in modo uniforme.
Si puo' andare a calcolare la varianza come
http://operaez.net/mimetex/\sigma^2 = \frac{\langle f^2 \rangle - \langle f\rangle^2}{N} = \frac{\frac{1}{N}\sum_{i=1}^N f(x_i)^2 - \left ( \frac{1}{N}\sum x_i \right )^2}{N}
Vediamo cosa succede:
subroutine mc(func,n,ss,err)
implicit real*8 (a-h,o-z)
dimension x(n)
external func
parameter (L=1d0)
parameter (NITMAX=1000000)
C***Inizializzazione
fmedia=0.d0
f2media=0.d0
call srand(1832)
do j=1,nitmax
do i=1,n
x(i)=2*L*(rand()-0.5)
enddo
f=func(x,n)
fmedia=fmedia+f
f2media=fmedia+f**2
enddo
ss=fmedia/NITMAX
err=sqrt(dabs(f2media/NITMAX-ss**2)/NITMAX)
ss=ss*(2*L)**n
return
end
Questa routine integra una funzione func (di n variabili: come e' facile vedere questo non porta alcuna complicazione!) fra -L ed L (eventualmente nell'ipercubo di lato 2L centrato nell'origine, se le dimensioni sono piu' di una!).
Prendendo come funzione f(x) = x^2 (l'integrale fa 2/3):
Nmax=1.0000 ----> I=0.6668659, err=5E-3
Nmax=1.00000 ----> I=0.6662768, err=1E-3
Nmax=1000.000 ----> I=0.6669116, err=5E-4
Nmax=10.000.000 ----> I=0.6664661, err=1E-4
Nmax=100.000.000 ----> I=0.6666737, err=5E-5
Come si vede l'integrale converge, anche se lentamente.
Il metodo, tuttavia, e' poco efficiente, in quanto si vanno a campionare allo stesso modo regioni dove la funzione e' grande e contribuisce molto all'integrale e regioni dove la funzione e' piccola e contribuisce poco. E' di Metropolis l'idea di rimediare a questo difetto scegliendo una funzione di distribuzione non piu' uniforme, ma che concentri il campionamento nelle regioni importanti!
L'idea e' sempre la stessa:
http://operaez.net/mimetex/\int_a^bf(x)dx=\int_a^b\frac{f(x)}{\mu(x)} \mu(x) dx = \int_a^b h(x) d\mu
Il problema e', a questo punto, come fare a generare una successione di punti distribuiti secondo la distribuzione mu.
Metropolis risolse il problema proponendo il seguente algoritmo (detto appunto algoritmo di metropolis):
Si calcolano la funzione e la distribuzione in un punto i
Ci si muove in modo casuale da punto i al punto i+1
Si accetta il movimento solo se soddisfa il test di metropolis, altrimenti si rimane nel punto vecchio
Ad ogni step si aggiornano le somme e si procede come nel montecarlo primitivo
Il test di metropolis e' molto semplice. Si calcola la funzione di distribuzione nel punto nuovo: se
http://operaez.net/mimetex/I=\frac{\mu(x_{i+1})}{\mu(x_i)}\geq1
si accetta il movimento, altrimenti si estrae un numero casuale fra zero ed uno (distribuito uniformemente, come quelli che restituisce il generatore rand() di fortran per intendersi!) e si confronta con il rapporto: se il rapporto e' piu' grande del numero estratto si accetta il movimento, altrimenti lo si respinge.
Dunque, in sistesi, la probabilita' di accettazione e'
http://operaez.net/mimetex/acc=\max\left (\frac{\mu(x_{i+1})}{\mu(x_i)},1\right )
Il codice e' autoesplicativo:
subroutine metropolis(func,dis,n,ss,err,acc)
implicit real*8 (a-h,o-z)
parameter (delta=1.d0, L=2.d0, NITMAX=1000000)
parameter (pigr=3.1415d0)
dimension xnew(n), xold(n)
external func, distr
fmedia=0.d0
f2media=0.d0
nacc=0
call srand(183)
do j=1,n
xold(j)=2*L*(0.5-rand())
enddo
fold=func(xnew,n)
disold=dis(xold,n)
do i=1,NITMAX
do j=1,n
xnew(j)=xold(j)+(rand()-0.5d0)*delta
if (xnew(j).lt.-L) xnew(j)=xnew(j)+2*L
if (xnew(j).gt.L) xnew(j)=xnew(j)-2*L
enddo
disnew=dis(xnew,n)
test=disnew/disold
if ((test.ge.1.d0).or.(test.gt.rand())) then
fold=func(xnew,n)
disold=disnew
nacc=nacc+1
do j=1,n
xold(j)=xnew(j)
enddo
endif
fmedia=fmedia+fold
f2media=f2media+fold**2
enddo
ss=fmedia/NITMAX
var=f2media/NITMAX-ss**2
err=sqrt(var/NITMAX)
acc=nacc*100.d0/NITMAX
return
end
La routine integra fra -L e L una funzione func (quella che prima ho chiamato h(x)!) rispetto ad una distribuzione dis(x), entrambe di n variabili, e restituisce il valore dell'integrale, la deviazione standard e il numero di mosse accettate.
Il montecarlo primitivo non e' altro, quindi, che un caso particolare del montecarlo metropolis!
In particolare, il metodo di metropolis si adatta perfettamente al problema della chimica teorica. Si supponga di voler calcolare in modo variazionale l'energia: Se si parametrizza con un insieme di parametri p la funzione d'onda, si dovra' minimizzare il funzionale
http://operaez.net/mimetex/E[\Psi] = \int d^3r_1d^3r_2\ldots d^3r_n \bar{\Psi}\mathscr{H}\Psi
dove ho supposto che la funzione d'onda fosse normalizzata ad uno.
Moltiplicando sopra e sotto per la funzione d'onda:
http://operaez.net/mimetex/E[\Psi] = \int d^3r_1d^3r_2\ldots d^3r_n \bar{\Psi}\Psi\frac{\mathscr{H}\Psi}{\Psi}=\int d^3r_1d^3r_2\ldots d^3r_n |\Psi|^2E_l
dove
http://operaez.net/mimetex/E_l=\frac{\mathscr{H}\Psi}{\Psi} viene definita ''energia locale''.
A questo punto, essendo il modulo quadrato della funzione d'onda una buona funzione di distribuzione, il gioco e' fatto e si puo' procedere all'integrazione, minimizzando il funzionale con le tecniche usuali!
Lucrezio
19-10-2008, 19:41
Questa e' un'idea di come funziona un Variational-Montecarlo.
Ho cercato di evitare i dettagli e di dare solo qualche idea, spero che si capisca!
Appena ho un po' di tempo posto qualche risultato per elio, litio e berillio (che sono i sistemi che ho provato a trattare)!
Suggerimenti e contestazioni sono benvenuti!
Dal forum di Rosetta@Home, giornale di David Baker (capo ricercatore)
While the results are still preliminary, it appears that Rosetta@home has produced an extremely exciting result! As I described a few posts ago, many of you through rosetta@home contributed to the design of proteins predicted to bind very tightly to the influenza flu virus. We have now completed the first round of testing of the designed proteins, and one of them in the experiments conducted thus far clearly binds very tightly to the virus. Our data also indicate that the binding is at a site critical to the virus invasion of our cells, and so the protein neutralize the virus. I will keep you posted over the next couple of months as the picture becomes clearer--but for now--thank you all for making this possible!!
SerPaguroSniffa³
07-10-2012, 19:02
ciao ho appena scaricato il sw, sto usando grid republic.
Guardando i log dei messaggi vedo che ci sono dei messaggi rossi apparentemente di errore:
http://i48.tinypic.com/14obrps.jpg
non vorrei lasciar lavorare il pc per niente.. voi sapete dirmi cosa sono?
vBulletin® v3.6.4, Copyright ©2000-2025, Jelsoft Enterprises Ltd.