Simulazioni Molecolari
Minimizzazione dell’energiaUn solo minimo
Ricerca conformazionaleMolti minimi
Simulazione molecolareGenera una successione di strutture che non sono minimi
Permette l’esplorazione dello spazio conformazionale
Si riferiscono a un insieme di molecole e non ad una
Applicazioni
Calcolo delle proprietà di un sistemaProprietà energeticheProprietà strutturaliProprietà termodinamicheProprietà dinamiche
Ricerca conformazionale
Studio del moto delle molecole
macroscopico
Generalità
• I metodi di simulazione molecolare permettono di studiare e prevedere il comportamento di sistemi macroscopici costituiti da un numero enorme di molecole ( 1023)
• Ciò è ottenuto tramite l’uso di modelli ridotti del sistema macroscopico che contengono un numero ridotto di molecole (10-1000)
• Le tecniche di simulazione molecolare sono due: Monte Carlo Dinamica molecolare
Meccanica statistica• In generale il valore di una proprietà di un sistema dipende dalla posizione e dai
momenti pi=mivi delle N particelle che variano nel tempo
x1(t), y1(t), z1(t), x2(t), y2(t) , z2(t) ……., xN(t), yN(t), zN(t)
px1(t), py1(t), pz1(t), px2(t), py2(t) , pz2(t) ……., pxN(t), pyN(t), pzN(t)
• Il valore istantaneo della proprietà A può essere scritto: A (PN (t), rN(t)) e fluttua nel tempo come risultato dell’interazione fra le particelle
• Il valore sperimentale della proprieta è la media di A nel tempo di misura:
• E’ quindi necessario conoscere come variano posizioni e velocità di ogni particella
nel tempo, cosa che è possibile in linea di principio ma impossibile in pratica per
un sistema macroscopico con 1023 particelle
0)()((
1lim
t
NNave dtttAA rP
• Per risolvere questo problema è stata sviluppata la meccanica statistica in cui un singolo sistema che si evolve nel tempo è sostituito da un gran numero di repliche del sistema, considerate però in uno stesso istante (insieme)
• La media temporale è quindi sostituita dalla media sull’insieme: in cui l’integrale è esteso a tutte le coordinate e i momenti del sistema
considerato, cioè d= dPNdrN
• Il sistema più comunemente usato in meccanica statistica é quello canonico in cui il numero di particelle N, il volume V e la temperatura T sono costanti
• In meccanica statistica si introduce poi il concetto di spazio delle fasi che
altro non è un’estensione dello spazio conformazionale da una a più
molecole
dAA NNNN ),(),( rPrP
Un sistema di N-particelle è completamente definito dalle 3N coordinate e dai 3N momenti P=mv
Un insieme di questi 6N valori (PN , rN)=(Px1, Py1,.., Pz3N ; x1, y1.., z3N)definisce uno stato del sistema o un punto nello spazio delle fasi.
Ciascuna proprietà sperimentalmente misurabile A è ottenuta come media pesata di quella proprietà sullo spazio delle fasi:
La funzione peso da la probabilità di trovare questo particolare stato del sistema nello spazio delle fasi, cioè la probabilità che il sistema abbia coordinate rN e momenti PN
dAA NNNN ),(),( rPrP
),( NN rP
Per l’insieme canonico (N V T costanti) la probabilità è data dall’espressione di Boltzmann function:
QTkE BNNNN /)/),(exp(),( rPrP
TkE
CQB
NN
NVT
),(exp
rP
• A (PN , rN) è il valore della proprietà A che il sistema assume nello stato (PN , rN) cioè in quel punto delle spazio delle fasi.
Ad esempio l’energia si ha:
E(PN , rN) = i=1,N Pi2/2mi + V(rN)
d
Ipotesi Ergodica
In meccanica statistica il valore di una proprietà A è quindi ottenuto come media sull’insieme <A>, per calcolare la quale dobbiamo conoscere la probabilita che il sistema sia in un certo stato (PN , rN) Alternativamente possiamo propagare un singolo sistema nello spazio delle fasi e ottenere il valore di A come media temporale:
L’ ipotesi ergodica stabilisce che:
o meglio, visto che di necessita dobbiamo usare un tempo finito:
0)()((
1lim
t
NNave dtttAA rP
aveAA
aveAA
Propagazione del sistema nel tempo
Dinamica MolecolareDeterministica
M
i
NNAM
A1
)(1
rP
M
i
NAM
A1
)(1
r
Monte CarloStocastica
In tutte le tecniche di simulazione molecolare si genera una successione nel
tempo di M configurazioni del sistema e per ciascuna si calcola A (PN , rN)
Proprietà Termodinamiche
M
iiEM
Eu1
1Energia Interna:
222 /}{ TkEEC Bv Capacità termica
N
i
N
ijijijB frTNk
VP
1 1311
Pressione:
)3(221
2
cB
N
i i
i NNTk
m
PK
Temperatura:
Impostazioni comuni della simulazione
Scelta del modello dell’energiaMeccanica molecolareAb initioSemiempirico
Scelta del modello di confine
Scelta della configurazione iniziale
Fase di equilibrazione Dalla configurazione iniziale all’equilibrio
Monitoraggio delle proprietà del sistema
Fase di produzione
Analisi della simulazione
Condizioni al contorno periodiche
Rimuove gli effetti di confinePermette il calcolo di proprie-tà macroscopiche usando un piccolo numero di particelle
Quando una particella esce da una cella è sostituita dalla sua immagine che entra dalla cella sul lato opposto
Il problema del campionamento
La successione di stati del sistema che si ottiene in una simulazione molecolare definisce un campionamento dello spazio delle fasi
Per ottenere buoni risultati il campionamento deve essere omogeneo e permettere di esplorare un’ampia parte dello spazio delle fasi
In particolare il campionamento è convergiuto se: le quantità medie rimangono stabili nel tempoSimulazioni da diverse condizioni iniziali danno gli stessi risultati
Metodo Monte Carlo
• Il metodo Monte Carlo genera casualmente delle configurazioni e usa degli particolari criteri per decidere se accettare o no ciascuna nuova configurazione nel set usato per il calcolo della proprietà A
• Questi criteri assicurano che la probabilità di ottenere una nuova configurazione sia uguale al fattore di Boltzmann
• Configurazioni di bassa energia vengono quindi generati con probabilità più alta che configurazioni di alta energia
• Alla fine dellla simulazione il valor medio della proprietà e ottenuto come semplice media
QTkE BNNNN /)/),(exp(),( rPrP
M
i
NAM
A1
)(1
r
• Nel metodo Monte Carlo l’energia è determinata solo dalla posizione degli atomi e non dai loro momenti
E(PN , rN) = V(rN)
• La limitazione dovuta a questo fatto è che non si calcolano i valori assoluti delle proprietà ma solo le deviazione rispetto al valore per il gas ideale (facilmente calcolabile analiticamente)
• Non c’è alcuna relazione temporale fra configurazioni successive ottenute in una simulazione Monte Carlo
Integrazione Monte Carlo
Area sotto la curva = _________________________________
n. di punti sotto la curva
n. di punti totali
Proprietà termodinamiche
trial
trial
N
i BN
i
N
i BN
iN
iN
TkV
TkVVV
1
1
]/)(exp[
]/)(exp[)()(
r
rrr
)()()( NNNN VdV rrrr
ZTkV B
NN ]/)(exp[
)(r
r
]/)(exp[ TkVdZ BNN
NVT rr
Energia potenziale media:
Usando l’integrazione Monte CarloSi generano Ntrial configurazioni casuali e si calcola:
Z
• Questo approccio non è però fattibile per l’alto numero di configurazioni con fattore di Boltzmann molto piccolo
• Questo problema è ovviato nell’approccio Metropolis in cui le configurazioni sono generate con probabilità uguale proprio al fattore di Boltzmann
• L’integrale per Q o per una qualsiasi proprietà A può quindi essere approsimato da una semplice somma
M
i
NAM
A1
)(1
r
Metropolis Monte Carlo
Metropolis Test
“Prova” Movimento casuale
E
NO
YES
E < 0or exp(-E/RT) < X[0,1]
??X[0,1] è un numero casuale fra 0 e 1
Si generano selettivamente configurazioni che danno contributi maggiori all’ integrale
Applicazioni
Simulazioni atomiche (esempio, cluster atomici)
Simulazioni di liquidi
Determinazione delle conformazioni molecolari (sistema di una sola molecola)
Polimeri
Esempio Sistema di tre atomi uguali in una cella cubica Calcolo dell’energia potenziale totale media del sistema
Configurazione iniziale
Energia potenziale E1:
])()[(4])()[(4),,( 6
13
12
13
6
12
12
12321 rrrrrrrV
])()[(4 6
23
12
23 rr
1
2
3
r13
r12
r23
1)
2) Variazione casuale della posizione di un atomo qualsiasi, ad esempio l’atomo 2
1
2
3
Nuova configurazione
Energia potenziale nuova configurazione E2:
])()[(4])()[(4),,( 6
13
12
13
6
12
12
12321 rrrrrrrV
])()[(4 6
23
12
23 rr
3)
5) Nel caso b) si effettua una variazione casuale dell’ultima configurazione, la seconda
1
2
3r23
r13
r12
4) Confronto energia configurazione precedente: a) Se E < 0 o exp(- E/RT) < X [0,1] la nuova configurazione è accettata b) Se exp(- E/RT) > X [0,1] la nuova configurazione non è accettata
1
2 3
1
2
3
Nel caso a) si riparte con una variazione casuale della precedente configurazione, la prima
6) Si valuta l’energia di questa terza configurazione e si prosegue nello stesso modo
7) Si memorizzano le energie potenziali di tutte le configurazioni accettate E1, E2, E3,…..
Dopo un certo numero di step, ad esempio 100, si conclude la simulazione e si calcola il valor medio dell’energia potenziale come semplice sommatoria:
100
1321 ),,(
100
1
iiVV rrr
Dinamica Molecolare
• A differenza del metodo Monte Carlo, nell’approccio della dinamica
molecolare si calcola la dinamica reale del sistema
• Applicando le equazioni del moto di Newton si calcolano le traiettorie di
tutte le N particelle del sistema
• La traiettoria di una particella specifica la posizione e la velocità di tale
particella in funzione del tempo xi(t), yi(t), zi(t) vxi(t), vyi(t), vzi(t) ; basta
conoscere le coordinate perché le velocità sono le derivate vxi(t) = dxi(t)/dt
• Data che la soluzione delle equazioni di Newton è numerica si ottiene una
successione di posizioni e velocita in M istanti di tempo successivi t1, t2,…,tM
: xi(t1), yi(t1), zi(t1) ; xi(t2), yi(t2), zi(t2) ; ………..; xi(tM), yi(tM), zi(tM)
• In pratica si ottiene una successione di M configurazioni che costituiscono un
campionamento dello spazio delle fasi di tutto il sistema (N particelle)
• Il valor medio di una proprietà A è ottenuta come media temporale
in pratica
in cui
A(PN rN) = A( PN(ti ),rN(ti ) )
è il valore istantaneo della proprietà A che dipende dalle posizioni e dai
momenti (velocità) delle N particelle nell’istante ti
• Si noti che a differenza del metodo Monte Carlo è necessario conoscere le
velocità di tutte le particelle e non solo le posizioni
0)()((
1lim
t
NNave dtttAA rP
M
i
NNAM
A1
)(1
rP
Propagazione temporale
Equazioni del moto di Newton:
)()( tmt iii aF
)(1
2
2
tmt ii
i Fr
)(1
2
2
tV
mt ii
i
r
r
Fi(r) = V(r1,r2,..,rN) / ri
Integrazione NumericaEspansione in serie di Taylor
Errttttttt )(21
)()()( 2avrr
Errttttttt )(21
)()()( 2bavv
f(t+t) = f(t) + (df(t)/dt)t + 1/2 (d2f(t)/dt2)t2 + ..
Per il moto di una particella ci interessa la sua posizione r e si ha:
In cui v(t) = r(t)/dt e a(t)= d2r(t)/dt2 e b(t) è la derivata terza
Algoritmo di Verlet
Algoritmo in praticaSi specificano le posizioni a t-t e t.Dalle positioni si calcolano le forze al tempo tSi calcolano le posizioni al tempo t+ t.Si calcolano le velocità
)(21
)()()( 2 ttttttt avrr
)(21
)()()( 2 ttttttt avrr
)()()(2)( 2 ttttttt arrr
)(1
)()(2)( 2 tm
tttttt Frrr
tttttt 2/)]()([)( rrv
i
j
Fij
-Fij
Calcolo della forza al tempo t
713
224
ijiji
ijij rr
VF
r
612
4ijij
ij rrV
Energia di interazione fra gli atomi i e jLennard-Jones
Forza sull’atomo i dovuta all’atomo j
Dinamica Molecolare
Calcola la nuova posizione
Calcola le forza dalla curvatura del potenziale Calcola le velocità dalla forza
T typically 10-15sec
Salva energia e geometria per la media
Scelta dell’intervallo di tempo
Atomi Traslazione 10-14
Molecole rigide Traslazione, rotazione 5x10-15
Molecole flessibiliLegami rigidi
Traslazione, rotazionetorsione
2x10-15
Molecole flessibiliLegami flessibili
Traslazione, rotazioneTorsione, vibrazione
10-15 or 5x10-16
Sistema Tipo di moto presente Intervallo suggerito
Per un sistema di atomi l’intervallo di tempo deve essere piccolo rispetto al tempo medio fra le collisioni
Per molecole flessibili l’ intervallo ottimale è approssimativamente un decimo del tempo per il movimento più veloce
Per molecole flessibili con legami flessibili il moto più veloce è quello di vibrazione, ad esempio per i legami C-H ~10 fs
Gli algoritmi SHAKE e RATTLE sono usati per bloccare i gradi di libertà interni e permettere intervalli di tempo più lunghi
Preparazione e avvio di una simulazione
Definizione configuraz. inizialeDerivata sperimentalmente Calcolata
Assegnazione delle velocità
Inizio della simulazioneFase di equilibrazione
Controllo della temperatura Fase di produzione
Acquisizione dei dati
VelocitàLe velocità sono assegnate casualmente da una distribuzione di Maxwell-Boltzmann alla temperatura di interesse:
Le componenti x, y, z delle velocità hanno distribuzione
2223
42
exp2
)( vTk
mvTk
mvf
BB
Energia Cinetica e Temperatura
Le velocità sono legate all’energia cinetica da:
2
1 21
ii
N
i
mK v
La temperatura è legata alle velocità da:
N
i
iiB mTkN
1
2
22)63( v
Dinamica e insiemi statistici
Una simulazione dinamica è in genere eseguita a N, V e E costanti (insieme microcanonico)
Può essere spesso utile eseguirla a N,V e T costanti (insieme canonico)
Per eseguire una dinamica a T costante dobbiamo ad ogni step al tempo t calcolare T(t) e scalare le velocità, cioè moltiplicarle per un fattore , tale che la temperatura assuma il valore voluto T
Infatti la temperatura è legata alle velocità da:
N
i
iiB mTkN
1
2
22)63( v
Si dimostra che tale valore vale:
)(tTT
ii λvv
Esempio Sistema di tre atomi uguali in una cella cubica Calcolo dell’energia potenziale totale media del sistema
Configurazione iniziale (r1, r2, r2) al tempo t1
Energia interazione V12 V13 e V23 forze F12, F13, F23
])()(2[24 7
12
13
1212 rrF
1
2
3
r13
r12
r23
1)
2)1
2
3
……….
Si calcola la nuova posizione al tempo t2=t1+t per
ogni particella i=1,2,3
)(1
)()(2)( 12
111 tm
tttttt iii Frrr
Solvatazione
• Le molecole di solvente polare si dispongono attorno allo ione o alla molecola con il lato di polarità opposta rispetto alla carica (parziale) ed esercitano un campo elettrico
• Questo campo elettrico da luogo ad un potenziale di interazione che deve essere aggiunto all’Hamiltoniano nell’equazione di Schroedinger
• Nei metodi di meccanica molecolare si aggiunge semplicemente un contributo elettrostatico nell’espressione empirica dell’energia e si considera poi una costante dielettrica nell’interazione elettrostatica tra i vari atomi
E importante poiché molte molecole sono in soluzione
Top Related