Estrarre la complessità dalle simulazioni
microscopiche
Fabrizio Cleri
ENEA, Unità Materiali e Nuove Tecnologie
Centro Ricerche Casaccia, Roma
I nostri collaboratori:
Manuela Volpe
Sara Letardi
Vittorio Sansalone
Mariella Ippolito
Andrea Moriani (ENEA Frascati)
Alessandro Mattoni (Cagliari)
Alessandra Satta (*)
Enrico Pisanu (*)
Manlio Messina (*)
Sonia Costantini (*)
Luigi Brambilla (*)
Staff Team:
Fabrizio Cleri (ENEA Casaccia)
Gregorio D’Agostino (ENEA Casaccia)
Pier Giuseppe Gabrielli (ENEA Casaccia)
Michele Gusso (ENEA Brindisi)
Principali collaborazioni in corso:
Luciano Colombo (Cagliari)
Giovanni Ciccotti (Roma “La Sapienza”)
Antonio Di Carlo (Roma Tre)
Patrizia Trovalusci (Roma “La Sapienza”)
Pawel Keblinski (Rennselaer, NY)
Susan B. Sinnott (U. Florida, Gainesville)
C. Delerue (ISEN Lille)
UTS Materiali e Nuove Tecnologie, Centro Ricerche Casaccia , Roma - Italy
Reattività di superfici,
catalisi, interazione
gas-matrice
Materiali nanostrutturati per
applicazioni funzionali
Microstrutture, interfacce, proprietà meccaniche
Multiscale modelling
Sistemi disordinati
F. Cleri
F. Cleri
G. D’Agostino,
F. Cleri
P.G. Gabrielli
M. Volpe
S. Letardi
V. Sansalone
M. Ippolito
Collaborazioni:
Cagliari (L. Colombo)
Progetti:
Nanocarbonio (ASI,
Conan)
Membrane H2/O2
FISR Idruri metallici
Nanoclusters
G. D’Agostino,
M. Gusso
Progetti:
FISR Nanomat
Collaborazioni:
Roma1 (G. Ciccotti)
ISEN Lille (C. Delerue)
RPI Troy (P. Keblinski)
Florida Un. (S. Sinnott)
Progetti:
TANO
NIRT/NSF
Danneggiamento
materiali irraggiati
F. Cleri,
A. Moriani (FUS)
Progetti:
ENEA/Euratom
Collaborazioni:
Roma1 (P. Trovalusci)
Roma3 (A. Di Carlo)
Cagliari (L. Colombo)
Collaborazioni:
Roma1 (G. Parisi)
Cagliari (L. Colombo)
Progetti:
5% PROMOMAT
COFIN 2003
ITEM FP5
Progetti:
5% PROMOMAT
COFIN 2003
ITEM FP5
Nel triennio 2001-2003:
~ 30 pubblicazioni e conference proceedings
16 invited talks
2 congressi internazionali organizzati
17 seminari di ospiti esterni
RISORSE HARDWARE/SOFTWARE
(Il “laboratorio di fisica computazionale”…?)
- Cluster di PC Linux (5 Athlon 1-1.8 Ghz + 2 twin-AMD 2.8 Ghz)
+ alcuni Mac G4/G5 OsX
- Sistema operativo Linux Debian Woody 3.1, GNU Fortran/C/C++, OpenOffice
- Grafica di base: Xmgr, GNUplot, PGPlot + pacchetti sviluppati in casa
- Grafica avanzata: JMOL, gOpenMol
Principali software di sviluppo:
- CAST (“Casaccia Atomistic Simulation Tools”)
comprende tutto il necessario per dinamica molecolare classica e tight binding + Monte
Carlo (HELLF), in diversi ensemble termodinamici fino al Gran-Canonico e con diversi
potenziali di interazione (EAM+secondo mom., L-J, Stillinger-Weber, Tersoff …)
- ABINIT, codice di simulazione di struttura elettronica e dinamica molecolare da principi
primi, sviluppato da Univ. Libre Bruxelles
- fhi98md, fhi98PP, codice di simulazione di struttura elettronica da principi primi (Fritz-Haber
- Institut, Berlino) e modulo per la generazione di pseudopotenziali a partire da funzioni
-d’onda all electron
DINAMICA MOLECOLARE - Modelli su scala atomica dei materiali
Gli atomi sono descritti esplicitamente come masse puntiformi che interagiscono mediante
forze su scala atomica.
L’evoluzione del sistema si ottiene integrando le equazioni newtoniane del moto per ciascun
atomo, nel campo di potenziale determinato dall’interazione complessiva tra gli atomi.
Le forze interatomiche si possono ottenere da un
potenziale empirico derivabile analiticamente…
Gli elettroni sono cioè “congelati” nello stato
fondamentale, rendendo il calcolo molto efficiente.
Si possono così studiare sistemi di dimensioni fino
a qualche frazione di micrometro.
Oppure, il legame esplicito tra atomi può essere
descritto mediante gli elettroni di legame. In questo
caso le forze hanno una base quanto-meccanica.
Questo approccio è computazionalmente molto più
pesante ed è quindi limitato a sistemi composti al
più da qualche decina-centinaia di atomi.
(J. Schiotz et al., Nature 391, 561 (1998))
METODI “MONTE CARLO” : Coarse-graining dei gradi di libertà
I metodi Monte Carlo sono basati sulla
stima di integrali multi-dimensionali
mediante tecniche di campionamento
stocastico.
L’integrale da calcolare può essere una
funzione termodinamica, un funzionale di
energia totale, il kernel di una equazione
integrale (ad es., equazione di
Boltzmann), o altro ancora.
Le variabili indipendenti sono solo un
sottoinsieme limitato di tutte le variabili
canoniche del sistema.
Campionando le variabili in base a una
funzione di distribuzione data, si calcola
il valore corrispondente dell’integrando e
se ne ottiene così una media statistica.
Simulazione Grand-Canonical Monte Carlo del
fisisorbimento di H2 in una matrice di nanotubi di
carbonio idrogenati.
(M. Volpe, F. Cleri, Chem. Phys. Lett. 231, 133 (2003))
RICERCA DELLA COMPLESSITA’
Le simulazioni su scala atomica coinvolgono un numero estremamente
grande di gradi di libertà, che si spinge facilmente ai limiti delle capacità di
memoria e di elaborazione dei più avanzati calcolatori disponibili.
- modelli atomici o molecolari con interazioni empiriche, in cui il numero di
atomi o molecole è estremamente elevato e le traiettorie di integrazione o
di campionamento sono molto lunghe,
- modelli ab-initio, in cui le funzioni d'onda degli elettroni sono espanse su
basi ortonormali di dimensione (spazio di Hilbert) molto grande.
L'analisi e l'estrazione di informazioni rilevanti da una tale mole di dati, che
supera facilmente le centinaia di GBytes per ciascuna simulazione, è un
compito, complesso con un onere computazionale confrontabile, o in alcuni
casi addirittura maggiore della stessa simulazione.
RICERCA DI COMPONENTI DI
MOMENTO ANGOLARE DEFINITO
IN UNA FUNZIONE D’ONDA
COMPLICATA
LA TEORIA DEL FUZIONALE DENSITA’ ELETTRONICA
A partire dal problema di Schroedinger per un sistema di ioni ed elettroni:
Introduciamo il concetto di “densità elettronica” come variabile univoca in 3 coordinate
{x,y,z}, a sostituire la funzione d’onda Y in 3N coordinate spaziali:
Mediante una formulazione variazionale, possiamo passare al set di N problemi
equivalenti di singola particella in un potenziale efficace:
Gli orbitali di Kohn e Sham, che entrano nella hamiltoniana efficace della formulazione
DFT del problema di Schroedinger, possono essere rappresentati in una qualsiasi
base ortonormale…come orbitali atomici, gaussiane, polinomi di Legendre… oppure
onde piane:
La scelta delle onde piane è conveniente per molti motivi:
- base molto semplice
- nessuna pre-concezione sulla soluzione
- assenza di errori di sovrapposizione (sovra-completezza
del set di base)
- facilità del calcolo delle forze (derivate prime)
- uso di un set discreto (ma infinito) di elementi grazie al
teorema di Bloch
BENZO-TIOLO/TIOLATO SU CLUSTER DI Au
Nelle simulazioni di struttura elettronica mediante la teoria del Funzionale Densità (DFT)
viene accuratamente descritta la chimica dei legami molecola-metallo, la struttura
energetica dei livelli molecolari, il rilassamento ionico, lo shift dei livelli molecolari…
PROBLEMA: Caratterizzare il legame chimico, alla ricerca delle cause della risposta
microscopica (in questo caso, la rotazione di 90°).
EVOLUZIONE DEGLI ORBITALI MOLECOLARI NEL CHEMISORBIMENTO
RSH-Au(3)
RS*-Au(3)
unrelaxed
RS-Au(3)
LE ARMONICHE SFERICHE Ylm(q,f),
AUTOFUNZIONI DEL MOMENTO ANGOLARE
Una qualsiasi funzione f(q,f) può essere sviluppata in armoniche sferiche:
E in particolare, un’onda piana:
ACCELERAZIONE DELLE
COORDINATE DI REAZIONE PER
EVENTI RARI
Un bacino di energia potenziale intrappola un sistema dinamico in uno stato metastabile.
(Le linee indicano contorni ad energia costante)
“Eventi rari”, sono quei processi attivati che corrispondono ad un percorso che collega
uno stato metastabile del sistema ad uno stato adiacente.
IL PROBLEMA DELLA SCALA DEI TEMPI
Le simulazioni microscopiche mediante algoritmi espliciti richiedono l’integrazione delle
equazioni del moto con dei dt molto più piccoli del tempo tipico di vibrazione t …
E un “evento raro” ha una frequenza di (milioni di t)
-1
…
LA RISPOSTA “TIPICA”: KINETIC MONTE CARLO
Il metodo Monte Carlo “cinetico” si basa su una catalogazione dei possibili “eventi rari”
che collegano due stati metastabili.
Mediante un algoritmo stocastico si campiona la distribuzione degli eventi rari, con
delle costanti di tempo inversamente proporzionali all’altezza delle barriere nei vari
punti di sella della ipersuperficie di energia potenziale.
Modello di diffusione di CH4 nella zeolite Na-Y: (a) identificazione dei siti stabili; (b) spanning di uno dei cammini
di transizione (“evento raro”); (c) barriere energetiche lungo il cammino di transizione.
(S.M. Auerbach et al., J. Phys. Chem. 100, 5923 (1996))
PROBLEMA: chi ci dice “a priori” quali sono gli eventi rari da tenere in conto?
REAZIONI CHIMICHE SU UNA SUPERFICIE CATALITICA
PROBLEMA: trovare il cammino di reazione più
probabile in funzione delle condizioni termodinamiche
per una molecola di acqua che viene adsorbita e
dissociata su una superficie, in questo caso TiO2-rutile
(100).
Il metodo “standard” è di inventarsi (con il buon
senso) i possibili cammini di reazione e far “muovere”
la molecola lungo la coordinata appropriata con la
“constrained minimization”.
Altrimenti si può provare con la dinamica molecolare
a temperatura tale da consentire la diffusione
spontanea…
Habs
H2O(g)
H2O(a)
OH(a)+H(a)
IL METODO DELLE REPLICHE PARALLELE
(A)
(B)
(C)
(D)
Il sistema viene replicato N volte (ad esempio, su N processori)
Le repliche vengono randomizzate (o “de-fasate”)
Ogni replica fa un certo percorso dinamico, finché
Un evento raro si verifica. A questo punto si torna ad (A)
IL METODO DELL’ “IPERDINAMICA”
Consiste nell’inventarsi un potenziale
fittizio (linea tratteggiata) che “allaga”
i bacini (stati metastabili) del potenziale
vero.
La condizione primaria è che il potenziale
fittizio V sommato a quello vero, sia zero
ad ogni superficie di separazione tra due
bacini (cioè ad ogni punto di sella).
Il “tempo” della dinamica accelerata risultante è stimato come:
n
t hyper  dt MD exp
i1

V (r(t i )
kB T

DINAMICA ACCELERATA DALLA TEMPERATURA
Il metodo si basa sulla semplice osservazione che per un processo attivato
la frequenza è data dalla legge di Arrhenius:
Rate   0 exp(EA /kBT)
Aumentando la temperatura T, quindi, si abbassa l’altezza delle barriere EA.
l’esponenziale fa sì che questa operazione abbia effetti non
Complicazione:
lineari per le diverse EA di una ipersuperficie di energia.
Processo di compattazione di una superficie di Cu mediante scambio concertato di tre atomi
(J.A. Sprague, F. Montalenti e A.F. Voter)
IL METODO DEI DIMERI, O “NUDGED ELASTIC BANDS”
Un “dimero” è una replica virtuale del sistema che passeggia sulla ipersuperficie
di energia alla ricerca dei vari punti di sella.
Variando l’orientazione e la lunghezza del passo tra il sistema e la sua replica, sulla
base di una stima locale della matrice delle derivate (Hessiano), il “dimero” segue il
percorso di minimo gradiente tra i due stati metastabili e individua “automaticamente” il
punto di sella.
Nell’ipotesi armonica, la “transition state theory” ci dà la frequenze dell’evento raro:
13N iinit (E saddleE init )/ kB T
Rate  3N 1 saddle e
1  i
che dipendono solo dalle caratteristiche del punto di partenza e del punto di sella.
In pratica il metodo stima il più piccolo autovalore e autovettore associato dell’Hessiano
senza mai dover
 calcolare l’intera matrice.
Essendo una ricerca casuale (si parte da molti punti infinitesimi vicino al punto di
partenza) non è dato che TUTTI i saddle point rilevanti vengano trovati.
PROSPETTIVE E PROBLEMI
Kinetic Monte Carlo - metodo molto efficiente, permette di superare barriere di tempi
estremamente lunghe. Richiede una enorme quantità di informazioni a priori sui
percorsi possibili e le relative barriere di energia.
Repliche parallele - metodo “banale”, non richiede alcuna conoscenza a priori.
Legato alla disponibilità di un numero elevatissimo di processori.
Hyperdynamics - metodo concettualmente semplice. La definizione del potenziale di
bias è però molto complessa. Se non è fatta in maniera molto studiata, il guadagno è
molto piccolo.
Time-accelerated MD - metodo piuttosto approssimativo, anche se capace di
elevatissimi speed-up. Richiede alcune ipotesi arbitrarie sull’andamento relativo delle
probabilità dei diversi eventi.
Dimeri - metodo efficace per casi di moti collettivi molto complessi. Basato
sull’approssimazione armonica, è piuttosto laborioso e funziona per sistemi
relativamente piccoli.
RICERCA DI EVENTI CINETICI
“CATASTROFICI” IN FASI
CONDENSATE
… cioè, problemi in cui si integrano le equazioni del moto di un sistema
contenente un numero enorme di particelle per tempi molto lunghi, e si
vanno a cercare fenomeni critici (o “catastrofici”) senza avere alcuna
informazione a priori su DOVE, QUANDO e COME si verificheranno…
Ad esempio:
- Nucleazione ed avanzamento di una microfrattura,
- Spostamento di un fronte di fase (liquido-solido, …)
- Nucleazione di una fase ordinata entro una fase disordinata
- Creazione di uno sciame di particelle secondarie durante irraggiamento
- ………
UN CASO SEMPLICE: AVANZAMENTO DI UN FRONTE SOLIDO-LIQUIDO
COESISTENZA DI
SILICIO CRISTALLINO
E AMORFO A T=1500 K
(L. Brambilla et al.,
Applied Phys. Lett. 2001)
Esiste una proprietà semplice
della fase ordinata, ad
esempio il fattore di struttura
S(k) lungo la direzione di
crescita, che permette di
evidenziare la regione ordinata
rispetto a quella disordinata.
INDIVIDUARE NUCLEI CRISTALLINI IN UN AMORFO
Quando non si ha alcuna informazione preventiva sulla posizione e sulla dinamica del
sistema bisogna inventarsi delle funzioni apposite che ricerchino i cluster di particelle
che soddisfano particolari condizioni:
i(1)= 1 + | ni/4 -1|
i(2)= j [ 1 + (dij/d0 Π1)2 ]/ ni
i(3)=
j k [ 2 - cos (qjik -q0) ]/ni(ni Π1)
2
(L. Brambilla et al., MRS Proceedings vol. 638 (2000))
CALCOLO DI FUNZIONI DI MERITO PER LA
COORDINAZIONE, BOND-LENGTH E ANGOLI
DI LEGAME:
MEDIA PER OGNI ATOMO AD OGNI TIMESTEP
SCALING… N3 !!!
SIMULAZIONE DI CASCATE DI IRRAGGIAMENTO
Anche ad energie moderate (10-100 eV) i primari e i secondari viaggiano per decine di nm
nel bulk. Pochissimi eventi significativi ad alta energia in un mare di atomi “imperturbati”.
Come in altri casi analoghi (frattura,
dislocazioni) le regioni interessanti vengono
identificate in base alla loro differenza
rispetto al bulk (coordinazione, energia…)
Dinamica Molecolare
danneggiamento primario
Sarebbe peraltro necessario
sviluppare algoritmi di tracking che
seguano in modo dinamico solo la
parte interessante della cascata di
spostamenti, riducendo di parecchi
ordini di grandezza il numero dei
gradi di libertà microscopici.
Kinetic Monte Carlo
propagazione su lunga scala
CASCATA PRODOTTA DA UN PRIMARIO DA 5 keV IN SiO2 AMORFO
0.16 ps
2.67 ps
28.6 nm
0.10 ps
I difetti in rosso sono deficienze di O, quelli in blu sono atomi di O che tornano a posto
(recovery), identificati come Si 4-coordinati che hanno cambiato uno o più vicini.
(M.J. Caturla et al., Laser-IFE Workshop, Livermore 2001)
In conclusione…
Le simulazioni al calcolatore offrono una grossa opportunità
alla scienza dei materiali: per la prima volta è infatti possibile
confrontare modelli teorici e dati sperimentali in scala 1:1
(almeno per le lunghezze…).
Sussistono limitazioni intrinseche di scala per le simulazioni
atomistiche, che non potranno essere superate neanche con
un eccezionale aumento della potenza di calcolo.
Viene in aiuto il paradigma multiscala, il quale però richiede un
deciso cambiamento di prospettiva.
Peraltro, l’improvement dei metodi teorici e degli algoritmi
numerici, unito allo sviluppo di tecniche di intelligenza
artificiale, potrà rendere solubili problemi oggi ancora lontani
dallo stadio di fattibilità tecnologica (e.g., materiali ibridi).