# Quantum Monte Carlo

Il quantum Monte Carlo (QMC) consiste in una grande famiglia di algoritmi sfruttati per simulazioni di sistemi quantistici nei campi di studio della fisica della materia condensata e della chimica computazionale. Questi algoritmi, pur differenziandosi tra loro per il diverso approccio che possono sfruttare, si basano tutti sul metodo Monte Carlo per la risoluzione degli integrali multidimensionali implicati.
I quantum Monte Carlo permette una rappresentazione diretta degli effetti delle repulsioni tra elettroni nella funzione d'onda, con una incertezza statistica che può essere ridotta aumentando la durata della simulazione. Per i bosoni esistono algoritmi numericamenti esatti e che variano in modo polinomiale con la dimensione del sistema oggetto di studio. Per i fermioni esistono invece ottime approssimazioni e algoritmi Monte Carlo numericamente esatti che variano in modo esponenziale, costituendo due differenti approcci risolutivi.

## Fondamenti

Un sistema fisico quantistico, con applicazioni che spazino in modo interdisciplinare, può essere descritto con l'equazione di Schrödinger a molti corpi. Ma questa, nella sua risoluzione matematica, risulta abbastanza complessa e richiede una certa quantità di tempo per essere risolta (in particolare per i fermioni, quali sono gli elettroni).
Tradizionalmente i metodi teorici approssimano la funzione d'onda a molti corpi in una funzione antisimmetrica di orbitali a singola particella (approssimazione degli orbitali):
    $\Psi (x_{1},x_{2},\dots ,x_{n})=f(\Phi _{1}(x_{1}),\Phi _{2}(x_{1}),\dots ,\Phi _{n}(x_{1});\Phi _{1}(x_{2})\Phi _{2}(x_{2}),\dots )$
come nel caso del metodo di Hartree-Fock. Questa formulazione semplifica molto i calcoli, ma introduce limitazioni, soprattutto legate al modo di valutare l'effetto delle repulsioni elettroniche, che per molte applicazioni sono inaccettabili.
Il quantum Monte Carlo rappresenta un insieme di metodi post-Hartree-Fock che permette di superare tali limitazioni sfruttando una funzione d'onda a molti corpi scelta da chi opera la simulazione. In particolare, si può utilizzare la funzione d'onda di Hartree-Fock come punto di partenza, procedendo poi a moltiplicarla per una specifica funzione simmetrica - tipiche sono le funzioni di Jastrow - che apporti un contributo correttivo. La maggior parte dei metodi mira al calcolo della funzione d'onda relativa allo stato fondamentale del sistema, con eccezione del path integral Monte Carlo e del auxiliary field Monte Carlo che invece calcolano la matrice densità.
Sebbene il metodo Monte Carlo dia risultati più accurati rispetto agli altri metodi ab initio di chimica quantistica, il suo costo computazionale non lo rende un metodo di routine, quali la teoria del funzionale di densità e il metodo Coupled Cluster. Di seguito sono descritti i metodi variational Monte Carlo, diffusion Monte Carlo e path integral Monte Carlo.

## Variational Monte Carlo

Il variational Monte Carlo (VMC) è un metodo quantum Monte Carlo che applica l'approccio variazionale per l'approssimazione dello stato fondamentale del sistema.
Il valore di aspettazione può essere scritto nella forma rappresentativa in x come
    ${\frac {\langle \Psi (a)|H|\Psi (a)\rangle }{\langle \Psi (a)|\Psi (a)\rangle }}={\frac {\int \Psi (X,a)^{2}{\frac {H\Psi (X,a)}{\Psi (X,a)}}dX}{\int \Psi (X,a)^{2}dX}}$.
Seguendo l'applicazione del metodo Monte Carlo per la valutazione degli integrali, il rapporto
    ${\frac {\Psi (X,a)^{2}}{\int \Psi (X,a)dX}}$
può essere assunto come una funzione di distribuzione della probabilità, valutando il valore dell'energia di aspettazione $E(a)$ come media della funzione locale ${\tfrac {H\Psi (X,a)}{\Psi (X,a)}}$, e quindi si procederà alla minimizzazione di $E(a)$.
Il metodo VMC è sostanzialmente simile a ogni classico metodo variazionale, ad eccezione del fatto che mentre gli integrali multidimensionali vengono valutati numericamente, nel VMC occorre invece calcolare solamente il valore della funzione d'onda, fattore che conferisce grande flessibilità al metodo in questione. Uno dei più grandi vantaggi in termini di accuratezza, derivanti dallo scrivere la funzione d'onda separatamente, deriva dall'introduzione del cosiddetto fattore Jastrow che esprime la funzione d'onda come
    $\exp(\sum {u(r_{ij})})$,
dove $r_{ij}$ è la distanza tra una coppia di particelle quantistiche.
Per i sistemi chimici, fattori leggermente più sofisticati permettono di ottenere un 80-90% dell'energia di correlazione utilizzando meno di 30 parametri. A confronto, a titolo di esempio, un calcolo di interazione di configurazione può richiedere circa 50.000 parametri per raggiungere lo stesso livello di accuratezza, sebbene questa dipenda fortemente dalla tipologia di sistema studiato.
Esistono essenzialmente tre metodi di ottimizzazione delle funzioni d'onda nel variational Monte Carlo:
* utilizzo di metodi deterministici (Stochastic Gradient Approximation, SGA);
* valutazione del costo computazionale della funzione e valutazione delle derivate;
* tecniche iterative.

## Diffusion Monte Carlo

Il metodo diffusion Monte Carlo (DMC) utilizza una funzione di Green per risolvere l'equazione di Schrödinger. Il DMC è potenzialmente un metodo numericamente esatto, il che significa che è in grado di trovare l'esatto valore dell'energia dello stato fondamentale all'interno di un dato intervallo di tolleranza definito per un qualsiasi sistema quantistico. Effettuando i calcoli, si trova che per i bosoni l'algoritmo fornisce dei valori che variano in modo polinomiale con la dimensione del sistema, mentre la variazione nel caso dei fermioni è di tipo esponenziale. Questo fatto rende impossibili le simulazioni DMC per i sistemi fermionici su larga scala; comunque è possibile ottenere risultati molto accurati applicando opportune approssimazioni (fixed-node approximation). Di seguito verrà presentato l'algoritmo basilare.

### Metodo del proiettore

Si prenda in considerazione l'equazione di Schrödinger per una particella nella sola dimensione della coordinata x:
    $i{\frac {d\Psi (x,t)}{dt}}={\frac {d^{2}\Psi (x,t)}{dx^{2}}}+V(x)\Psi (x,t).$
È possibile ottenere una forma contratta di questa equazione introducendo l'operatore hamiltoniano H, con $H={\frac {d^{2}}{dx^{2}}}+V(x)$. Quindi l'equazione assume la forma
    $i{\frac {d\Psi (x,t)}{dt}}=H\Psi (x,t)$.
Esistono particolari funzioni, chiamate autofunzioni, per le quali $H\Psi (x)=E\Psi (x)$, dove E rappresenta un valore numerico. Queste funzioni sono speciali perché nel valutare l'effetto dell'operatore H sulla funzione d'onda, si ottiene sempre lo stesso valore di E (conservazione dell'energia totale della particella quantistica). Queste funzioni sono dette stati stazionari, perché la loro derivata in funzione del tempo, in ogni punto x, è sempre la stessa e quindi l'ampiezza della funzione d'onda non varia al variare del tempo.
Solitamente si è interessati alla funzione d'onda caratterizzata dal più basso autovalore di energia, ovvero alla determinazione dello stato fondamentale. Scrivendo l'equazione di Schrödinger in modo leggermente diverso è possibile ottenere sempre lo stesso autovalore di energia ma questo, piuttosto che avere carattere oscillatorio, sarà convergente. Tale forma, nota anche come equazione di Schrödinger tempo immaginario, è
    ${\frac {-d\Psi (x,t)}{dt}}=(H-E_{0})\Psi (x,t)$,
in cui è stato eliminato il numero immaginario dalla derivata temporale e se ne tiene conto nel secondo membro dove compare la grandezza $E_{0}$, che è l'energia dello stato fondamentale. Questa equazione modificata possiede alcune utili proprietà. La prima cosa da notare è che quando si vuole determinare la funzione d'onda dello stato fondamentale si ottiene $H\Phi _{0}(x)=E_{0}\Phi _{0}(x)$ e la derivata temporale assume valore uguale a zero. Supponendo adesso di utilizzare un'altra funzione d'onda di partenza ($\Psi$), che non rappresenta lo stato fondamentale ma è non ortogonale a esso. Quindi questa funzione può essere scritta come somma lineare di autofunzioni:
    $\Psi =c_{0}\Phi _{0}+\sum _{i=1}^{\infty }c_{i}\Phi _{i}$.
Trattandosi di una equazione differenziale lineare, si può considerare il contributo di ogni parte separatamente. Si è già detto che $\Phi _{0}$ è stazionaria. Supponiamo di considerare $\Phi _{1}$: essendo $\Phi _{0}$ l'autofunzione a più bassa energia, l'autovalore associato a $\Phi _{1}$ soddisfa la relazione $E_{1}>E_{0}$. Così la derivata temporale di $c_{1}$ è negativa ed eventualmente tenderà a zero, fornendo solamente lo stato fondamentale. Da questa osservazione trae origine un modo di determinare $E_{0}$: osservando l'ampiezza della funzione d'onda nella sua propagazione temporale, un suo incremento o decremento permette di stimare il valore dell'energia.

### Implementazione stocastica

Precedentemente si è arrivati a ottenere una equazione che, valutandone l'evoluzione temporale e utilizzando un valore appropriato di $E_{0}$, permette di ottenere lo stato fondamentale di un qualsivoglia hamiltoniano. Questo rappresenta ancora un grosso problema nella meccanica classica, perché piuttosto che considerare la propagazione di una singola posizione delle particelle occorre invece considerare la propagazione di intere funzioni. Nella meccanica classica il moto delle particelle si potrebbe simulare sulla base dell'equazione $x(t+\tau )=x(t)+\tau v(t)+0.5F(t)\tau ^{2}$, assumendo che la forza sia costante.
Per l'equazione di Schrödinger tempo immaginario, la propagazione temporale viene invece effettuata utilizzando un integrale di convoluzione con una funzione speciale definita funzione di Green. In tal modo si ottiene
    $\Psi (x,t+\tau )=\int G(x,x',\tau )\Psi (x',t)dx'$.
Similmente al caso della meccanica classica si possono effettuare propagazioni solamente per piccole porzioni di tempo, altrimenti la funzione di Green diviene non accurata. Aumentando il numero di particelle di cui è composto il sistema quantistico, aumenta anche la dimensione dell'integrale, dovendo integrare in tutte le coordinate di tutte le particelle. L'integrazione può essere effettuata sfruttando il metodo Monte Carlo.

## Integrale Monte Carlo

L'integrale sui cammini Monte Carlo sfrutta la formulazione dell'integrale sui cammini della meccanica quantistica, sviluppata da Richard Feynman generalizzando il concetto dell'azione della meccanica classica. In pratica una particella quantistica può passare da un punto A nel tempo $t=0$ a un punto B nel tempo $t_{1}$ seguendo differenti possibili percorsi.
L'integrale sui cammini, generalizzato per i problemi quantomeccanici, è
    $Z=\int e^{i{\mathcal {S}}/\hbar }\operatorname {d} \!x$
in cui l'azione viene definita a partire dal tempo t=0 a t=T, x la posizione.