Simulazioni Monte Carlo e Dinamica Molecolare

Simulazioni Molecolari
Minimizzazione dell’energia
Un solo minimo
Ricerca conformazionale
Molti minimi
Simulazione molecolare
Genera 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 sistemamacroscopico
Proprietà energetiche
Proprietà strutturali
Proprietà termodinamiche
Proprietà dinamiche
Ricerca conformazionale
Studio del moto delle molecole
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:
Aave  lim
 
•
1

 t 0
A(P N (t )r N (t )dt
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
• 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:
 A   A(P N , r N ) (P N , r N )d

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
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:
 A   A(P N , r N ) (P N , r N )d

La funzione peso  (P , r ) 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
N
N
•
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)
Per l’insieme canonico (N V T costanti) la probabilità è data
dall’espressione di Boltzmann function:
 (P N , r N )  exp(  E (P N , r N ) / k BT ) / Q
QNVT
 E (P N , r N ) 
 C  exp 
d

k BT 

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:
Aave  lim
 
1

 t 0
A(P N (t )r N (t )dt
L’ ipotesi ergodica stabilisce che:
 A  Aave
o meglio, visto che di necessita dobbiamo usare un tempo finito:
 A  Aave
Propagazione del sistema nel tempo
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)
Dinamica Molecolare
Deterministica
1
 A 
M
M
N N
A
(
P
r )

i 1
Monte Carlo
Stocastica
1
 A 
M
M
N
A
(
r
 )
i 1
Proprietà Termodinamiche
Energia Interna:
1
u  E 
M
M
E
i 1
i
Capacità termica
Cv  { E 2    E  2 } / k BT 2
Pressione:

1
1 N N
P   Nk BT    rij f ij 
V
3 i1 j i1

N
Temperatura:
2
Pi
k BT
K 

(3 N  N c )
2
i 1 2mi
Impostazioni comuni della simulazione
Scelta del modello dell’energia
Meccanica molecolare
Ab initio
Semiempirico
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 confine
Permette il calcolo di proprietà 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 tempo
Simulazioni 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
 (P N , r N )  exp(  E (P N , r N ) / k BT ) / Q
• 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
1
 A 
M
M
N
A
(
r
 )
i 1
• 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
n. di punti sotto la curva
Area sotto la curva = _______________________________
__
n. di punti
totali
Proprietà termodinamiche
Energia potenziale media:  V (r N )   dr NV (r N )  (r N )
N
exp[

V
(
r
) / k BT ]
 (r N ) 
Z
Z NVT   dr N exp[ V (r N ) / k BT ]
Z
Usando l’integrazione Monte Carlo
Si generano Ntrial configurazioni casuali e si calcola:
 V (r
N

) 
N trial
i 1
Vi (r N ) exp[ Vi (r N ) / k BT ]

N trial
i 1
exp[ Vi (r N ) / k BT ]
• 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
1
 A 
M
M
N
A
(
r
 )
i 1
Metropolis Monte Carlo
Si generano selettivamente configurazioni che danno contributi
maggiori all’ integrale
Movimento casuale
“Prova”
DE
NO
YES
Metropolis
Test
DE < 0
or
exp(-DE/RT) < X[0,1]
?
X[0,1] è un numero casuale fra 0 e 1
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
1)
Configurazione iniziale
1
Energia potenziale E1:
r13
r12
3
r23
V (r1 , r2 , r3 )  4 [(

r12
 4 [(
2
2)
)12  (

r23

r12
)12  (
) 6 ]  4 [(

r23

r13
)12  (
)6 ]
Variazione casuale della posizione di
un atomo qualsiasi, ad esempio
l’atomo 2
1
3
2

r13
)6 ]
3)
Nuova configurazione
1
r13
3
r12
r23
Energia potenziale nuova configurazione
E2:
V (r1 , r2 , r3 )  4 [(
2

r12
 4 [(
) (

r23
12

r12
)12  (
) ]  4 [(
6

r23

r13
) (
12

r13
)6 ]
)6 ]
Confronto energia configurazione
precedente:
a) Se DE < 0 o exp(- DE/RT) < X [0,1] la nuova
configurazione è accettata
b) Se exp(- DE/RT) > X [0,1] la nuova
configurazione non è accettata
4)
5)
1
2
3
Nel caso a) si
riparte con una
variazione
casuale della
precedente
configurazione,
la prima
3
1
2
Nel caso b) si
effettua una
variazione
casuale
dell’ultima
configurazione,
la seconda
6)
7)
Si valuta l’energia di questa terza configurazione e si
prosegue nello stesso modo
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:
1
 V 
100
100
V (r , r , r )
i
i 1
1
2
3
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
Aave  lim
 
1

 t 0
A(P N (t )r N (t )dt
in pratica
1
 A 
M
M
N N
A
(
P
r )

i 1
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
Propagazione temporale
Equazioni del moto di Newton:
Fi (t )  mi ai (t )
 2ri 1
 Fi (t )
2
t
mi
Fi(r) =  V(r1,r2,..,rN) /  ri
 2ri 1 V

2
t
mi ri (t )
Integrazione Numerica
Espansione in serie di Taylor
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:
1 2
r (t  t )  r (t )  tv(t )  t a(t )  Err
2
1 2
v(t  t )  v(t )  ta(t )  t b(t )  Err
2
In cui v(t) = r(t)/dt e a(t)= d2r(t)/dt2 e b(t) è la derivata terza
Algoritmo di Verlet
1 2
r (t  t )  r (t )  tv(t )  t a(t )
2
1 2
r (t  t )  r (t )  tv(t )  t a(t )
2
r(t  t )  2r(t )  r(t  t )  t 2a(t )
1
r (t  t )  2r (t )  r (t  t )  t
F(t )
m
v(t )  [r (t  t )  r (t  t )] / 2t
2
Algoritmo in pratica
Si specificano le posizioni a t-t e t.
Dalle positioni si calcolano le forze al tempo t
Si calcolano le posizioni al tempo t+ t.
Si calcolano le velocità
Calcolo della forza al tempo t
j
-Fij
Energia di interazione fra gli atomi i e j
Lennard-Jones
Fij


Vij  4 
 rij

i
Forza sull’atomo i dovuta all’atomo j
Vij 24   
2
Fij 

ri
   rij

13
7

 
    

 rij  

  
12
6

  
   

 rij  

  
Dinamica Molecolare
Calcola le forza dalla curvatura del potenziale
Calcola le velocità dalla forza
DT
Salva energia e geometria per la media
typically 10-15sec
Calcola la nuova posizione
Scelta dell’intervallo di tempo
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
Sistema
Tipo di moto presente
Intervallo suggerito
Atomi
Traslazione
10-14
Molecole rigide
Traslazione, rotazione
5x10-15
Molecole flessibili
Legami rigidi
Traslazione, rotazione
torsione
2x10-15
Molecole flessibili
Legami flessibili
Traslazione, rotazione
Torsione, vibrazione
10-15 or 5x10-16
Preparazione e avvio di una
simulazione
Definizione configuraz. iniziale
Derivata sperimentalmente
Calcolata
Assegnazione delle velocità
Inizio della simulazione
Fase 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:
 m 
f (v )  

 2k BT 
Le componenti x, y, z delle
velocità hanno
distribuzione
32
 mv 2 
2
exp 
4

v

2
k
T

B 
Energia Cinetica e Temperatura
Le velocità sono legate all’energia cinetica da:
N
1
 K   mi v i2
i 1 2
La temperatura è legata alle velocità da:
(3N  6)k BT N mi v i2

2
2
i 1
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:
(3N  6)k BT N mi v i2

2
2
i 1
Si dimostra che tale valore vale:
  T T (t )
vi  λvi
Esempio
Sistema di tre atomi uguali in una cella cubica
Calcolo dell’energia potenziale totale media del sistema
1)
1
Configurazione iniziale (r1, r2, r2) al tempo t1
r13
r12
Energia interazione V12 V13 e V23 forze F12, F13, F23
3
r23
F12 
24

[2(

r12
)13  (

r12
)7 ]
……….
2
2)
Si calcola la nuova posizione al tempo t2=t1+t
per ogni particella i=1,2,3
1
3
2
1
ri (t1  t )  2ri (t1 )  ri (t1  t )  t
F(t1 )
m
2
Solvatazione
E importante poiché molte molecole sono in soluzione
•
•
•
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