Appunti per le esercitazioni di Chimica Teorica

Giorgio F. Signorini

Date: 2007-2008

Indice

• Metodi Monte Carlo
  • Importance Sampling
    • dimostrazione del valore della varianza della media del campione
  • Catene di Markov
  • Schema del metodo Metropolis Monte Carlo
  • Punti salienti di un programma MC
  • Calcolo di $g(r)$
    • in pratica
  
  
• Bibliografia

Metodi Monte Carlo

Importance Sampling

(basato su ())

Cerchiamo il valore di aspettazione di una funzione $g(x)$ distribuita secondo una probabilità $p(x)$ sul dominio delle $x$ :

$$E\left[g(x)\vert p\right]=\int g(x)p(x)dx$$

(assumiamo $\int p(x)dx=1$ ).

Una stima empirica di $p(x)$ , generando $x_{1},x_{2},\ldots$ secondo $p$ , è

$$P_{n}(x)=\frac{1}{n}\sum_{P}\delta(x-x_{i})$$

La corrispondente stima empirica di $g$ è

$$\hat{E}\left[g(x)\vert p\right] = \int g(x)P_{n}(x)dx$$

$$= \frac{1}{n}\sum_{P}g(x_{i})$$

e

$$\lim_{n\rightarrow\infty}\hat{E}\left[g(x)\vert p\right]\equiv\left\langle g\right\rangle _{P}$$

Se però non siamo in grado di generare punti secondo $p$ , ma secondo un'altra distribuzione $q$ (nel caso banale, una distribuzione uniforme), si può scrivere il valore di aspettazione come

$$E\left[g(x)\vert p\right] = \int g(x)\frac{p(x)}{q(x)}q(x)dx$$

$$= \int g(x)r(x)q(x)dx$$

$$= E\left[g(x)\cdot r(x)\vert q\right]$$

con

$$r(x)=\frac{p(x)}{q(x)}$$

La stima empirica su $q$ del valore di aspettazione di $g(x)$ è

$$\hat{E}\left[g(x)r(x)\vert q\right] = \frac{1}{n}\sum_{Q}g(x_{i})r(x_{i})$$

e

$$\lim_{n\rightarrow\infty}\hat{E}\left[g(x)r(x)\vert q\right]\equiv\left\langle g\cdot r\right\rangle _{Q}$$

Quanto è accurata una certa stima di $E$ fatta con $n$ punti? Cioè, quanto si discosta $\hat{E}$ dal valore vero $E$ ? Dipende da $n$ ma anche dalle proprietà di $g(x)$ e di $p$ (o $q$ ). Ad esempio, se $g(x)=G=cost$ , il valore di aspettazione stimato con qualunque $n$ è $G$ , e l'errore è nullo per qualunque $n$ . Negli altri casi, a parità di $n$ , possiamo aspettarci che sarà tanto minore quanto minore è la variazione di $g(x)$ in un campionamento secondo $p$ .

L'accuratezza della stima è misurata dalla varianza della media $s^{2}$ (da non confondere con la varianza dei dati, $\sigma^{2}$ ).

Si dimostra (vedi ()) che

$$s_{P,n}^{2} \approx \frac{1}{n}\left(\left\langle g^{2}\right\rangle _{P}-\left\langle g\right\rangle _{P}^{2}\right)$$ (1.1)

$$= \frac{\sigma_{P}^{2}}{n}$$

L'equazione () conferma che a parità di $n$ l'errore sulla stima è dato dalla variazione quadratica della $g$ , e che se $g=cost$ allora $s^{2}=0$ .

Se invece che su $P$ , campioniamo su $Q$ , l'errore sulla stima è

$$s_{Q,n}^{2} = \frac{1}{n}\left(\left\langle \left(g\cdot r\right)^{2}\right\rangle _{Q}-\left\langle g\cdot r\right\rangle _{Q}^{2}\right)$$ (1.2)

cioè dipende dalla variazione quadratica di $g(x)\cdot r(x)$ .

L'equazione () corrisponde esattamente alla (3.1.9) in (), con $f=g\cdot p$ e $w=q$ .

Consideriamo due casi limite:

1. Nel caso banale di una distribuzione $U$ uniforme, $q=\frac{1}{x_{max}-x_{min}}=\frac{1}{L}$ , e

   $$s_{U,n}^{2}=\frac{1}{n}\left(\left\langle \left(g\cdot p\right)^{2}\right\rangle _{U}-\left\langle g\cdot p\right\rangle _{U}^{2}\right)$$ (1.3)

   
   

   Il confronto tra quest'ultima equazione e la () ci dice che la differenza tra fare una stima del valore di aspettazione di $g(x)$ campionando $x$ in modo uniforme e campionandola secondo la sua distribuzione $p(x)$ è data dalla differenza nella variazione quadratica di $g\cdot p$ e di $g$ rispettivamente. Se $p(x)$ è molto ``piccata'' (come è la distribuzione nell'insieme canonico) la differenza può essere enorme. Se al contrario il prodotto $g\cdot p$ varia meno di $g$ , cioè $g$ varia molto proprio dove $p$ è bassa, l'errore nella stima può essere più basso se la stima viene calcolata con campionatura uniforme.

2. Se $q(x)$ è una buona approssimazione di $p(x)$ , allora $r(x)\approx1$ , e
   $$s^{2}\approx\frac{1}{n}\left(\left\langle g^{2}\right\rangle _{Q}-\left\langle g\right\rangle _{Q}^{2}\right)$$

dimostrazione del valore della varianza della media del campione

Si immagina^1.1 di fare $N\rightarrow\infty$ stime ciascuna basata su $n$ punti: $\hat{E}_{1}$ , $\hat{E}_{2}$ , $\ldots$ , e di fare la varianza di queste stime:

$$s^{2} = \frac{1}{N}\sum_{k}\left(\hat{E}_{k}-\hat{E}\right)^{2}$$

con

$$\hat{E}_{k} = \frac{1}{n}\sum_{i\in P}g(x_{i}^{k})$$

$$\hat{E} = \frac{1}{N}\sum_{k}\hat{E}_{k}=\frac{1}{Nn}\sum_{k}\sum_{i}g(x_{i}^{k})$$

Con notazione sintetica ( $g_{i}^{k}=g(x_{i}^{k})$ e $i=\left\{ i\in P\right\}$ ) si ha:

$$s^{2} = \frac{1}{N}\sum_{k}\left(\hat{E}-\frac{1}{n}\sum_{i}g_{i}^{k}\right)^{2}$$

$$= \frac{1}{N}\sum_{k}\left(\frac{1}{n}\sum_{i}\left[g_{i}^{k}-\hat{E}\right]\right)^{2}$$

$$= \frac{1}{N}\sum_{k}\frac{1}{n^{2}}\left(\sum_{i}\left[g_{i}^{k}-\hat{E}\right]\right)\left(\sum_{j}\left[g_{j}^{k}-\hat{E}\right]\right)$$

se i punti $i,j$ sono estratti in maniera casuale sul dominio $P$ , i termini misti si annullano:

$$s^{2} \approx \frac{1}{N}\sum_{k}\frac{1}{n^{2}}\sum_{i}\left[g_{i}^{k}-\hat{E}\right]^{2}$$

$$= \frac{1}{n}\frac{1}{N}\sum_{k}\frac{1}{n}\sum_{i}\left[\left(g_{i}^{k}\right)^{2}-\hat{E}^{2}\right]$$

$$\approx \frac{1}{n}\left[\left\langle g^{2}\right\rangle _{P}-\left\langle g\right\rangle _{P}^{2}\right]$$

$$= \frac{\sigma^{2}}{n}$$

dove il passaggio al valore di aspettazione vero $\left\langle \right\rangle$ si ha nel limite $N\rightarrow\infty$ ; l'ultima uguaglianza può essere presa come una definizione di $\sigma^{2}$

$$\sigma^{2}=\left\langle g^{2}-\left\langle g\right\rangle _{P}^{2}\right\rangle _{P}$$

In sostanza, nel limite $N\rightarrow\infty$

$$s^{2} = \left\langle \left(g-\left\langle g\right\rangle _{P}\right)^{2}\right\rangle _{P}$$

$$= \frac{1}{n}\left\langle g^{2}-\left\langle g\right\rangle _{P}^{2}\right\rangle _{P}$$

$$= \frac{1}{n}\sigma^{2}$$

Catene di Markov

(basato su ())

Il problema è generare punti con una distribuzione di probabilità che approssimi la distribuzione canonica.

In una catena di Markov

• c'è un numero finito di stati
• ogni stato dipende solo dal precedente.

In qs condizioni si può costruire una matrice di transizione, $\pi$ . Si ha

$$\rho^{(2)} = \rho^{(1)}\pi$$

$$\rho^{(3)} = \rho^{(2)}\pi=\left[\rho^{(1)}\pi\right]\pi$$

$$\ldots$$

$$\rho^{(n)} = \rho^{(1)}\pi^{n}$$

La catena ammette un limite

$$\rho=\lim_{\tau\rightarrow\infty}\rho^{(1)}\pi^{\tau}$$

tale che

$$\rho\pi=\rho$$

Dall'ultima equazione si capisce che lo stato limite $\rho$ esiste se $\pi$ ammette autovalore $1$ . Se il numero di stati è finito, lo è anche il numero di vettori $\rho^{(i)}$ . (approfondire: come si può essere sicuri che si arriva comunque al limite $\rho$ ?)

Nell'esempio del sistema a due stati (il calolatore), si suppone che $\pi$ sia nota.

Nelle simulazioni MC, si conosce $\rho$ (ogni elemento $\rho_{m}=\rho_{NVT}(\Gamma_{m})$ ), e si vuole trovare $\pi$ .

Vi sono varie possibili soluzioni.

Si può porre la condizione (non necessaria) della reversibilità microscopica:

$$\rho_{m}\pi_{mn}=\rho_{n}\pi_{nm}$$

La soluzione asimmetrica di Metropolis prevede ($\alpha$ è una matrice simmetrica qualsiasi)

$$\begin{array}{lcc} \pi_{mn}=\alpha_{mn} & \rho_{n}>\rho_{m} ... ...o_{m} & m\neq n\\ \pi_{mm}=1-\sum_{n\neq m}\pi_{mn}\end{array}$$

Schema del metodo Metropolis Monte Carlo

1. selezionare a caso una particella e calcolare la sua energia, $U(\mathbf{r^{N}})$
2. spostarla a caso, e calcolare la sua nuova energia, $U(\mathbf{\mathbf{r'^{N}}})$
3. calcolare $p=e^{-\beta\left(U(r')-U(r)\right)}$ :
   • se $p\geq1$ , accettare la mossa
   • se $p<1$ , estrarre un numero $\zeta$ a caso nell'intervallo $[0,1]$ :
     • se $\zeta<p$ , accettare la mossa
     • altrimenti, lasciare la particella nella vecchia posizione
4. accumulare le medie e ripartire dal punto 1.

Punti salienti di un programma MC

1. algoritmo base (v. punto precedente)()
2. a ogni mossa tentata non occorre calcolare tutta E, ma solo il $\Delta E$ dovuto allo spostamento della molecola spostata, $i$ .

   N.B. che questa variazione è $\Delta E(i)=\sum_{j}\Delta E(i,j)$ , e se la mossa viene accettata, va aggiornata sia la $E(i)$ che tutte le altre $E(j)=\sum_{k}E(j,k)$ , di una quantità $\Delta E(j)=\Delta E(j,i).$ In sostanza si deve conservare in memoria ogni singola interazione $E(i,j)$ .

3. Bilancio dettagliato: qualunque schema si adotti, la probabilità a priori di tornare indietro deve essere uguale a quella di andare avanti: $\alpha_{o\rightarrow n}=\alpha_{n\rightarrow o}$ . Ad esempio: evitare di eseguire spostamenti casuali con un ordine fisso, tipo $T_{x},T_{y},T_{z},R_{x},R_{y},R_{z},T_{x},\ldots$ , perché la mossa inversa è impossibile.
4. campionamento efficiente: v Frenkel, una buona misura di efficienza è la quantità di spazio delle fasi esplorato nell'unità di tempo, e di questo una buona misura è $\frac{\Delta r^{2}}{T_{CPU}}$ . In particolare:
   1. conviene muovere una particella per volta
   2. non è detto che un rapporto di accettazione $=50\%$ sia ottimale.
5. Muovere $T_{x},T_{y},T_{z}$ insieme o solo uno per volta? Traslazioni e rotazioni insieme o separate?
6. Quanto deve essere grande $\Delta$ ?
7. PBC-MinimumImage: occhio a mantenere molecole intere

Calcolo di $g(r)$

in pratica

In pratica il $g(r)$ si calcola come:

$$g_{AB}\left(r+\frac{\Delta}{2}\right)=\frac{\left\langle n^{AB}\left(r,\Delta\right)\right\rangle }{n_{ideale}^{AB}\left(r,\Delta\right)}$$

cioè il rapporto tra il numero medio di particelle B a distanza da una particella A compresa tra $r$ e $\left(r+\Delta\right)$

$$\left\langle n^{AB}\left(r,\Delta\right)\right\rangle =\frac{1}{TN_{A}}\sum_{t}^{T}\sum_{i}^{N_{A}}n_{it}^{AB}\left(r,\Delta\right)$$

e lo stesso in un sistema ideale

$$n_{ideale}^{AB}\left(r,\Delta\right) = \rho_{B}\Delta V\left(r,\Delta\right)$$

$$n_{ideale}^{AA}\left(r,\Delta\right) = \rho_{A}\frac{N_{A}-1}{N_{A}}\Delta V\left(r,\Delta\right)$$

Quindi:

$$g_{AB}\left(r+\frac{\Delta}{2}\right)=\frac{\sum_{t}^{T}\sum_{i}^... ...a\right)}{T\left(N_{A}-\delta_{AB}\right)\rho_{B}\Delta V\left(r,\Delta\right)}$$

Note pratiche:

1. Il calcolo di $n_{it}^{AB}$ comprende comunque un loop sugli atomi di tipo B (per contarli)
2. Per le molecole $i,j$ si contano solo le coppie $j\geq i$ , ma all'interno di ciascuna delle due, tutti gli atomi $A,B$ : $\left\{ i_{A}j_{A},i_{A}j_{B},i_{B}j_{A},i_{B}j_{B}\right\}$
3. Le interazioni di tipo $A-A$ e $B-B$ sono la metà di quelle $A-B$ (in cui vanno anche le $B-A$ ); ogni volta, vanno contate due volte tanto

No References!

... immagina^1.1

cfr. () e ()