background

PortLab

CMS per la generazione di portali web dedicati a progetti di ricerca

Equazione di Schrödinger

Questa pagina risolve l'equazione di Schrödinger di alcuni sistemi quantistici con il metodo Monte Carlo variazionale (VMC): si campiona la densità |ψ|² di una funzione d'onda di prova con l'algoritmo di Metropolis e si stima l'energia come media dell'energia locale. Ogni sistema è implementato in quattro linguaggi (Javascript, PHP, Python, C++) con una regola di ferro: a parità di seed i quattro backend producono risultati numericamente identici — stesso generatore di numeri casuali (xorshift64*), stesso ordine di estrazione, stessa traiettoria stocastica. Il confronto dei tempi diventa così un esperimento controllato: cambia solo il linguaggio (~2 ordini di grandezza tra C++ e Python). Il backend Javascript gira direttamente nel browser (Web Worker); gli altri girano sul server. Oltre alla suite multi-linguaggio ci sono i sistemi a molti elettroni con determinante di Slater generale (Be, Ne, N₂, CH₄) e il motore generale a funzione d'onda esterna: la chimica sta in un file JSON generato una volta con Hartree-Fock (PySCF), il motore C++ è unico e permette sia il VMC sia il DMC fixed-node, il gradino più alto della gerarchia HF → VMC → DMC. Questa scheda ne raccoglie il complemento teorico, con le derivazioni complete — dall'equazione di Schrödinger al principio variazionale, dal campionamento di Metropolis al DMC fixed-node — e chiude con la bibliografia.

1. L'equazione di Schrödinger

Il punto di partenza è l'equazione di Schrödinger stazionaria [1] $$ \hat H \,\Psi(\mathbf R) = E \,\Psi(\mathbf R), $$ dove \( \mathbf R = (\mathbf r_1, \dots, \mathbf r_N) \) raccoglie le coordinate degli \( N \) elettroni. Tutta la suite lavora in unità atomiche (\( \hbar = m_e = e = 4\pi\varepsilon_0 = 1 \)): le lunghezze sono in bohr (\( a_0 \simeq 0{,}529\,\text{Å} \)) e le energie in hartree (\( 1\,\text{Ha} \simeq 27{,}211\,\text{eV} \)). Nell'approssimazione di Born–Oppenheimer i nuclei sono fissi e l'Hamiltoniana elettronica è $$ \hat H = -\frac{1}{2}\sum_{i=1}^{N} \nabla_i^2 \;-\; \sum_{i,A} \frac{Z_A}{|\mathbf r_i - \mathbf R_A|} \;+\; \sum_{i<j} \frac{1}{r_{ij}} \;+\; \sum_{A<B} \frac{Z_A Z_B}{R_{AB}}, $$ con \( r_{ij} = |\mathbf r_i - \mathbf r_j| \). Il termine di repulsione elettrone–elettrone rende il problema non separabile: la funzione d'onda vive in \( 3N \) dimensioni e una griglia con \( M \) punti per dimensione costerebbe \( M^{3N} \) valutazioni. È la maledizione della dimensionalità, ed è il motivo per cui si campiona con il metodo Monte Carlo [10, 13]: l'errore statistico scala come \( 1/\sqrt{n} \) indipendentemente dalla dimensione.

1.1 Stati stazionari ed evoluzione unitaria (il visualizzatore)

Il resto della pagina risolve l'equazione stazionaria; la scheda “Visualizzatore dell'idrogeno” mostra invece l'equazione dipendente dal tempo per l'unico atomo con soluzione analitica chiusa. Gli autostati dell'idrogeno sono $$ \psi_{nlm}(\mathbf r) = R_{nl}(r)\, Y_l^m(\theta,\varphi), \qquad E_n = -\frac{1}{2 n^2}\ \text{Ha} = -\frac{13{,}6}{n^2}\ \text{eV}, $$ con \( R_{nl} \) parte radiale (polinomi di Laguerre) e \( Y_l^m \) armoniche sferiche. Un singolo autostato è stazionario: \( |\psi|^2 \) non dipende dal tempo. Una loro sovrapposizione evolve invece in modo unitario, $$ \Psi(\mathbf r,t) = c_A\,\psi_A\,e^{-iE_A t} + c_B\,\psi_B\,e^{-iE_B t}, $$ e la densità \( |\Psi(t)|^2 \) oscilla («batte») con pulsazione \( \Delta E = |E_A - E_B| \) e periodo $$ T = \frac{2\pi}{\Delta E} $$ (in unità atomiche; il visualizzatore lo riporta in fs). Il colore di ciascun punto codifica la fase locale \( \arg\Psi \), e la nuvola è campionata da \( |\Psi|^2 \) con lo stesso algoritmo di Metropolis (§2) che, reso quantitativo, diventa il VMC: è il ponte visivo tra il campionamento stocastico e il metodo di calcolo delle altre schede.

2. Principio variazionale e Monte Carlo variazionale (VMC)

Per ogni funzione di prova \( \Psi_T \) normalizzabile vale il principio variazionale: $$ E_V \;=\; \frac{\langle \Psi_T | \hat H | \Psi_T \rangle} {\langle \Psi_T | \Psi_T \rangle} \;\ge\; E_0 , $$ con uguaglianza se e solo se \( \Psi_T \) è lo stato fondamentale. Il passo chiave del VMC [3, 10] è riscrivere il funzionale come media pesata sulla densità \( |\Psi_T|^2 \): $$ E_V = \int \pi(\mathbf R)\, E_L(\mathbf R)\, d\mathbf R, \qquad \pi(\mathbf R) = \frac{|\Psi_T(\mathbf R)|^2}{\int |\Psi_T|^2\, d\mathbf R}, \qquad E_L(\mathbf R) = \frac{\hat H \Psi_T(\mathbf R)}{\Psi_T(\mathbf R)} , $$ dove \( E_L \) è l'energia locale. La media si stima campionando \( \pi(\mathbf R) \) con l'algoritmo di Metropolis [2]: da \( \mathbf R \) si propone \( \mathbf R' = \mathbf R + \delta \) (spostamento uniforme di semi-ampiezza step) e si accetta con probabilità $$ A(\mathbf R \to \mathbf R') = \min\!\left(1,\; \left|\frac{\Psi_T(\mathbf R')}{\Psi_T(\mathbf R)}\right|^2\right). $$ La catena così costruita soddisfa il bilancio dettagliato \( \pi(\mathbf R)\,P(\mathbf R\to\mathbf R') = \pi(\mathbf R')\,P(\mathbf R'\to\mathbf R) \) e converge a \( \pi \) senza che serva la costante di normalizzazione. I primi passi (burn_in) si scartano; i campioni successivi sono correlati, quindi l'errore statistico si stima con le medie a blocchi: la corsa è divisa in \( N_b \) blocchi (default 100), e $$ \sigma_{E}^2 \simeq \frac{1}{N_b(N_b-1)} \sum_{b=1}^{N_b} \left( \bar E_b - \bar E \right)^2 . $$

Vale inoltre il principio di varianza zero: se \( \Psi_T \) è un'autofunzione esatta, \( E_L(\mathbf R) = E_0 \) è costante e la varianza \( \sigma^2(E_L) \) si annulla. La varianza, oltre all'energia, misura quindi la qualità della funzione di prova — ed è il test di correttezza assoluto usato dall'atomo di Hooke (§4.4).

3. Fattore di Jastrow e condizioni di cuspide di Kato

Le funzioni di prova della suite hanno la forma prodotto \( \Psi_T = \Phi \cdot J \), dove \( \Phi \) fissa la struttura orbitale (prodotti o determinanti) e \( J = e^{\sum_{i<j} u(r_{ij})} \) è il fattore di Jastrow [4], qui in forma di Padé: $$ u(r) = \frac{k\, r}{1 + \beta r} . $$ Il suo compito è correlare esplicitamente le coppie di elettroni e, in particolare, curare le divergenze del potenziale coulombiano. Quando due particelle cariche coincidono, \( 1/r \to \infty \): perché \( E_L \) resti finita, il termine cinetico deve produrre una divergenza uguale e opposta. Ne seguono le condizioni di cuspide di Kato [5]: $$ \left.\frac{1}{\Psi}\frac{\partial \Psi}{\partial r_{iA}}\right|_{r_{iA}=0} = -Z_A \qquad\text{(elettrone–nucleo)}, $$ $$ \left.\frac{1}{\Psi}\frac{\partial \Psi}{\partial r_{ij}}\right|_{r_{ij}=0} = \frac{1}{2} \;\;\text{(spin antiparalleli)}, \qquad \frac{1}{4} \;\;\text{(spin paralleli)} , $$ dove il valore ridotto per gli spin paralleli discende dall'antisimmetria (la funzione d'onda si annulla già linearmente alla coalescenza). Nel Padé la pendenza all'origine è \( u'(0) = k \): scegliendo \( k = 1/2 \) o \( 1/4 \) per coppia si impone la cuspide elettrone–elettrone esatta, mentre \( \beta \) resta il parametro variazionale che ne regola la portata. Le cuspidi elettrone–nucleo sono invece imposte dagli esponenti degli orbitali (o, nel motore generale, dalla correzione di cuspide, §8).

4. I sei sistemi della suite multi-linguaggio

4.1 Elio

Funzione di prova alla Padé–Jastrow con orbitali \( 1s \) a cuspide esatta (\( Z = 2 \)): $$ \Psi_T = e^{-2 r_1}\, e^{-2 r_2}\, \exp\!\left[ \frac{r_{12}}{2\,(1+\alpha\, r_{12})} \right] . $$ L'energia locale ha forma chiusa (derivazione in Thijssen [12]): $$ E_L = -4 + \frac{(\hat{\mathbf r}_1 - \hat{\mathbf r}_2)\cdot(\mathbf r_1 - \mathbf r_2)} {r_{12}\,(1+\alpha r_{12})^2} - \frac{1}{r_{12}\,(1+\alpha r_{12})^3} - \frac{1}{4\,(1+\alpha r_{12})^4} + \frac{1}{r_{12}} . $$ Con \( \alpha = 0{,}15 \) si ottiene \( E = -2{,}87756 \pm 0{,}00061 \) Ha, contro il valore esatto \( -2{,}90372 \) Ha [15].

4.2 Ione H2+

Un solo elettrone e due centri a distanza \( R \): LCAO \( \phi = e^{-a r_A} + e^{-a r_B} \). La cuspide esatta su entrambi i nuclei richiede che, alla coalescenza su A, il contributo dell'orbitale centrato su B sia tenuto in conto: la condizione dà l'equazione di punto fisso $$ a = 1 + e^{-aR} , $$ risolta per iterazione. L'energia totale (con \( 1/R \)) ha il minimo attorno a \( R \simeq 2{,}0\text{–}2{,}2 \) bohr con \( E \simeq -0{,}576 \) Ha, qualitativamente corretto per una LCAO minimale (il legame a un elettrone è dimostrato).

4.3 Idrogeno molecolare H2

Due elettroni nello stesso orbitale di legame, per ciascuno \( \phi(i) = e^{-a r_{iA}} + e^{-a r_{iB}} \) con la stessa cuspide di punto fisso, più il Jastrow con \( k = 1/2 \) (spin antiparalleli). Detti \( L_i = \nabla_i^2 \phi_i / \phi_i \) e \( \mathbf F_i = \nabla_i \phi_i / \phi_i \), l'energia cinetica locale è $$ \text{KE} = -\frac{1}{2}\Big[ L_1 + L_2 + 2u'' + \frac{4u'}{r_{12}} + 2u'\,(\mathbf F_1 - \mathbf F_2)\cdot \hat{\mathbf r}_{12} + 2u'^2 \Big] , $$ dove i termini in \( u \) vengono da \( \nabla^2 e^{u} / e^{u} = u'' + 2u'/r + u'^2 \) applicato a entrambi gli elettroni e dal doppio prodotto tra gradiente orbitale e gradiente del Jastrow. Con \( \beta = 0{,}7 \) il minimo cade a \( R = 1{,}40 \) bohr con \( E = -1{,}1513 \pm 0{,}0009 \) Ha: il 98% del valore esatto \( -1{,}1745 \) Ha di Kołos e Wolniewicz [16].

4.4 Atomo di Hooke (caso a varianza zero)

Due elettroni in trappola armonica con repulsione coulombiana: $$ \hat H = -\frac{1}{2}(\nabla_1^2 + \nabla_2^2) + \frac{\omega^2}{2}(r_1^2 + r_2^2) + \frac{1}{r_{12}} . $$ Con \( \Psi_T = e^{-\omega (r_1^2+r_2^2)/2}\, e^{u(r_{12})} \) l'energia locale è $$ E_L = 3\omega - u'' - \frac{2u'}{r_{12}} + \omega\, u' r_{12} - u'^2 + \frac{1}{r_{12}} . $$ Taut [9] ha mostrato che per particolari \( \omega \) la soluzione è analitica: per \( \omega = 1/2 \) lo stato fondamentale ha fattore di correlazione \( J = 1 + r_{12}/2 \) ed energia esattamente \( E = 2 \) Ha. In modalità exact la suite riproduce \( E = 2{,}0000000000 \) con varianza nulla: ogni configurazione dà la stessa energia locale, il test di correttezza più severo che un codice VMC possa superare (§2).

4.5 Ione H3+

Due elettroni e tre protoni ai vertici di un triangolo equilatero di lato \( d \) nel piano \( xy \) (circumraggio \( d/\sqrt3 \)). LCAO a tre centri con cuspide dal punto fisso \( a = 1 + 2e^{-ad} \) (alla coalescenza su un nucleo contribuiscono gli altri due), Jastrow con \( k=1/2 \), repulsione nucleare \( 3/d \). Minimo a \( d \simeq 1{,}65\text{–}1{,}9 \) bohr con \( E \simeq -1{,}248 \) Ha (93% del riferimento \( -1{,}3438 \) Ha [17]), sotto l'energia di H2: la stabilità dello ione — uno dei più abbondanti del mezzo interstellare — è dimostrata.

4.6 Idruro di litio (LiH)

Quattro elettroni (2↑, 2↓): il vincolo di antisimmetria è imposto da determinanti di Slater \( 2\times2 \) separati per spin, \( \Phi = D^{\uparrow} D^{\downarrow} \), con orbitali $$ \phi_1 = e^{-a_1 s}, \qquad \phi_2 = e^{-a_2 t} + c\, e^{-a_3 s} , $$ dove \( s \) e \( t \) sono le distanze da Li e H. Il determinante è invariante per \( \phi_2 \to \phi_2 + \lambda\, \phi_1 \) (combinazioni di colonne), una proprietà usata come test. Il Jastrow è a coppie con cuspide spin-dipendente (\( k = 1/4 \) paralleli, \( 1/2 \) antiparalleli, §3). Ai default \( E = -7{,}988 \pm 0{,}003 \) Ha, praticamente il limite Hartree–Fock (\( -7{,}987 \)); l'esatto non relativistico è \( \simeq -8{,}070 \) Ha: la differenza è l'energia di correlazione che il singolo determinante non può catturare.

5. Determinante di Slater generale (Be, Ne, N2, CH4)

Per i sistemi a molti elettroni la suite usa una macchina generale a determinante \( N\times N \) per specie di spin: costruita la matrice di Slater $$ A_{ij} = \phi_j(\mathbf r_i) , $$ il modulo del determinante si valuta come \( \ln|\det A| \) tramite decomposizione LU con pivoting (stabile anche vicino ai nodi) e l'inversa con Gauss–Jordan. Derivare \( \ln D \) è algebra dell'inversa: $$ \nabla_i \ln D = \sum_j (A^{-1})_{ji}\, \nabla \phi_j(\mathbf r_i), \qquad \frac{\nabla_i^2 D}{D} = \sum_j (A^{-1})_{ji}\, \nabla^2 \phi_j(\mathbf r_i) . $$ Con il Jastrow a coppie, l'energia cinetica locale per elettrone si assembla come $$ \text{KE}_i = T_i + \sum_{j \ne i}\Big( u'' + \frac{2u'}{r_{ij}} \Big) + 2\, \mathbf d_i \cdot \mathbf j_i + |\mathbf j_i|^2 , $$ dove \( T_i \) è il termine determinantale puro, \( \mathbf d_i = \nabla_i \ln D \) e \( \mathbf j_i = \nabla_i \ln J \). La struttura del codice è identica per i quattro sistemi; cambiano solo gli orbitali \( \phi_j \) (tipo Slater con esponenti variazionali, simmetrie molecolari fissate: D∞h per N2, Td per CH4). I motori riportano anche i riferimenti Hartree–Fock ed esatto per il confronto.

6. Il motore generale a funzione d'onda esterna

Il salto architetturale dello Stadio C è separare il motore dalla chimica: la funzione d'onda diventa un file JSON generato una volta da un calcolo Hartree–Fock ristretto (RHF) con PySCF [19], e il motore C++ è unico: $$ \Psi_T = D^{\uparrow}[\phi_k]\; D^{\downarrow}[\phi_k]\; J, \qquad \phi_k = \sum_{\mu} C_{k\mu}\, \chi_\mu , $$ dove i \( \chi_\mu \) sono orbitali atomici gaussiani contratti di tipo \( s \) e \( p \): $$ \chi_\mu(\mathbf r) = P_\mu(\mathbf u) \sum_i c_i\, e^{-a_i u^2}, \qquad \mathbf u = \mathbf r - \mathbf R_c , $$ con \( P_\mu = 1 \) (tipo \( s \)) oppure \( u_x, u_y, u_z \) (tipo \( p \)). Definendo \( S_n = \sum_i c_i\, a_i^n\, e^{-a_i u^2} \), le derivate necessarie all'energia locale sono in forma chiusa: $$ s:\;\; \chi = S_0, \quad \nabla\chi = -2 S_1 \mathbf u, \quad \nabla^2\chi = 4 S_2 u^2 - 6 S_1 ; $$ $$ p:\;\; \chi = u_m S_0, \quad \nabla\chi = -2 S_1 u_m \mathbf u + S_0\, \hat{\mathbf e}_m, \quad \nabla^2\chi = u_m (4 S_2 u^2 - 10 S_1) . $$ La validazione è a tre livelli, dal più forte al più debole: (i) puntuale: l'energia locale del motore coincide fino a \( \sim 10^{-6} \) con il laplaciano numerico di \( \ln\Psi \) calcolato in modo indipendente; (ii) statistica: con Jastrow spento (\( \beta = 0 \)) e cuspidi disattivate il VMC riproduce l'energia RHF di PySCF entro l'errore; (iii) fisica: col Jastrow acceso l'energia scende sotto l'HF restando sopra l'esatta (bound variazionale, §2).

6.1 Il parco molecolare

La libreria comprende 18 funzioni d'onda: le prime 7 di base (H2O e CH4 in STO-3G e 6-31G, etano, etanolo, benzene) e 11 aggiunte per coprire, ciascuna, un concetto chimico distinto. Tutte a determinante singolo RHF a guscio chiuso; se non indicato diversamente la base è STO-3G.

MolecolaConcetto illustrato
F2Il caso emblematico del DMC: a livello Hartree–Fock F2 risulta quasi non legato; VMC e DMC recuperano l'energia di legame corretta — mostra dal vivo perché serve la correlazione.
C2H4 (etene)Doppio legame C=C (ibridazione sp2).
C2H2 (acetilene)Triplo legame C≡C (sp); con l'etano (C–C singolo, sp3) completa la serie sp3→sp2→sp.
(H2O)2 (dimero)Legame a idrogeno: un'interazione debole (~5 kcal/mol) che l'HF sottostima e che il DMC recupera.
NH3 (ammoniaca)Doppietto solitario e inversione a ombrello.
H2CO (formaldeide)Gruppo carbonile C=O.
CH3OH (metanolo)Compagno più leggero dell'etanolo, per confronto di serie.
CO2Molecola lineare con doppi legami cumulati (simmetria D∞h).
HF (acido fluoridrico)Diatomica fortemente polare, con riferimento "esatto" ben noto.
C5H5N (piridina)Eteroaromatico: come il benzene ma con un azoto nell'anello — stesso costo di calcolo, un eteroatomo in più.
C10H8 (naftalene)Estende la serie aromatica oltre il benzene (68 elettroni): banco di prova della scalabilità del metodo.
6.2 Cosa questo motore non può generare

La funzione di prova è un singolo determinante RHF a guscio chiuso per il Jastrow, all-electron con correzione di cuspide. Questa scelta — semplice e robusta — lascia fuori quattro famiglie di sistemi:

  • Multireferenziali / biradicalici (C2, O3, molte specie a legame in rottura): un solo determinante dà nodi qualitativamente sbagliati e quindi un grande errore di nodo fisso (§10). Servirebbe una funzione d'onda multi-determinante.
  • A guscio aperto (radicali, O2 nel suo tripletto fondamentale): il generatore impone spin totale nullo (RHF); servirebbero orbitali ROHF/UHF e determinanti di dimensione diversa per i due spin.
  • Basi con funzioni d o f (6-31G*, cc-pVDZ, …): il motore valuta solo orbitali di tipo \( s \) e \( p \) (§6).
  • Elementi pesanti (terza riga in poi, metalli di transizione): in all-electron i core profondi fanno esplodere costo ed energia locale; servirebbero pseudopotenziali, non implementati.

In tutti questi casi la strada non è "aggiungere una molecola" ma estendere il formato del wf e il motore (multi-determinante, ROHF/UHF, pseudopotenziali): un lavoro di ben altra portata rispetto alla semplice esportazione di un JSON.

7. Mosse a singolo elettrone e formula di Sherman–Morrison

Muovere tutti gli elettroni insieme e rifare la LU costerebbe \( O(N^3) \) per mossa. Il motore muove un elettrone alla volta: spostando \( \mathbf r_i \to \mathbf r_i' \) cambia una sola riga di \( A \), e il rapporto dei determinanti si legge dall'inversa mantenuta senza ricalcolare nulla: $$ \rho = \frac{\det A'}{\det A} = \sum_j \phi_j(\mathbf r_i')\,(A^{-1})_{ji} . $$ Se la mossa è accettata, l'inversa si aggiorna in \( O(N^2) \) con la formula di Sherman–Morrison [18]: $$ S_k = \sum_l \phi_l(\mathbf r')\,(A^{-1})_{lk} \qquad (S_i = \rho), $$ $$ (A^{-1})'_{jk} = (A^{-1})_{jk} - \frac{(A^{-1})_{ji}\, S_k}{\rho} \;\; (k \ne i), \qquad (A^{-1})'_{ji} = \frac{(A^{-1})_{ji}}{\rho} . $$ Un "campione" diventa uno sweep (una mossa proposta per elettrone) più una valutazione dell'energia dallo stato mantenuto, con cache per elettrone di valori, gradienti e laplaciani degli orbitali molecolari: nel loop caldo non c'è nessuna decomposizione. Per igiene numerica l'inversa viene ricostruita da zero ogni refresh sweep (default 250). Guadagno misurato di efficienza (metrica \( \sigma^2 \times t \)) su H2O: 19× rispetto alle mosse all-electron.

8. Correzione di cuspide per basi gaussiane

Le gaussiane hanno derivata nulla al nucleo: non possono soddisfare la cuspide di Kato (§3), quindi \( E_L \) diverge come \( -Z/r \) vicino ai nuclei e la varianza esplode. Il rimedio è lo schema di Ma, Towler, Drummond e Needs [11] (qui in versione semplificata): dentro un raggio \( r_c \) attorno a ogni nucleo, la parte \( s \) di ciascun orbitale molecolare viene sostituita da $$ R(r) = \operatorname{sgn} \cdot e^{\,p(r)}, \qquad p(r) = \alpha_0 + \alpha_1 r + \alpha_2 r^2 + \alpha_3 r^3 + \alpha_4 r^4 , $$ con cinque vincoli: la cuspide esatta \( p'(0) = -Z_{\text{eff}} \), il raccordo \( C^2 \) a \( r_c \) (continuità di \( p, p', p'' \)) e il valore centrale \( \alpha_0 \). Il punto delicato — scoperto a spese della varianza di H2O, prima peggiorata — è che \( \alpha_0 \) non va fissato al valore HF: va lasciato libero e ottimizzato minimizzando la deviazione dell'energia locale orbitale $$ E_L^{\text{orb}}(r) = -\frac{1}{2}\left( p'' + p'^2 + \frac{2p'}{r} \right) - \frac{Z_{\text{eff}}}{r} $$ dentro la sfera, così che il polinomio resti "domato" (sul core dell'ossigeno \( \alpha_4 \) passa da \( +2\cdot10^4 \) a \( -78 \)). Con la ricetta completa la varianza di H2O crolla da 966 a 39,5 Ha² (24×). Il fitting avviene nell'esportatore; il motore legge i coefficienti pronti dal JSON. Correzione di cuspide e Jastrow elettrone–nucleo imporrebbero la cuspide due volte: il motore rifiuta la combinazione.

9. Ottimizzazione variazionale di β con campionamento correlato

Il parametro \( \beta \) del Jastrow va ottimizzato: il principio variazionale garantisce che il minimo di \( E_V(\beta) \) è la scelta migliore. Rifare un VMC per ogni \( \beta \) sarebbe costoso e rumoroso (minimi diversi per il rumore statistico); si usa invece il campionamento correlato: si campiona una sola volta a \( \beta_{\text{rif}} \) e si riappesa $$ E(\beta) \simeq \frac{\sum_n w_n E_L^{(\beta)}(\mathbf R_n)}{\sum_n w_n}, \qquad w_n = \left| \frac{\Psi_\beta(\mathbf R_n)}{\Psi_{\beta_{\text{rif}}}(\mathbf R_n)} \right|^2 = e^{\,2[U_\beta(\mathbf R_n) - U_{\beta_{\text{rif}}}(\mathbf R_n)]} , $$ dove le parti determinantali si semplificano nel rapporto: si riusano i termini \( T \) e \( \mathbf d \) già calcolati e cambiano solo i termini di Jastrow — l'intera curva \( E(\beta) \) esce da una passata. La frazione di campioni effettivi \( w_{\text{eff}} = (\sum w_n)^2 / (n \sum w_n^2) \) segnala quando l'overlap tra le due distribuzioni degrada. Due accorgimenti pratici: il riferimento si sposta iterativamente sull'argmin (con rifinitura parabolica), e l'esponente dei pesi va ricentrato sul primo campione — nei sistemi grandi \( \Delta U \) ha una parte costante enorme che satura il clamp numerico e renderebbe i pesi tutti uguali (un \( w_{\text{eff}} = 1 \) fasullo), ma i pesi sono definiti a meno di una costante che nel rapporto si cancella. Un risultato non ovvio: più la base è buona, più \( \beta \) piccoli fanno danno (il Jastrow a lungo raggio distorce la densità monoelettronica), quindi il \( \beta \) ottimo cresce con la qualità della base (H2O: 7,51 in STO-3G, 4,66 in 6-31G).

10. Diffusion Monte Carlo (DMC) fixed-node

Il VMC è limitato dalla forma di \( \Psi_T \). Il DMC [6, 7, 20] proietta via le eccitazioni risolvendo l'equazione di Schrödinger in tempo immaginario \( \tau = it \): $$ -\frac{\partial \Phi}{\partial \tau} = (\hat H - E_T)\,\Phi \qquad\Longrightarrow\qquad \Phi(\tau) = e^{-(\hat H - E_T)\tau}\, \Phi(0) \;\xrightarrow[\tau\to\infty]{}\; c_0\, \Psi_0 , $$ perché ogni autocomponente decade come \( e^{-(E_n - E_T)\tau} \) e la più lenta a morire è lo stato fondamentale. Campionare \( \Phi \) direttamente è inefficiente (e impossibile dove cambia segno); si campiona invece la distribuzione mista \( f = \Psi_T\, \Phi \), che soddisfa $$ -\frac{\partial f}{\partial \tau} = -\frac{1}{2}\nabla^2 f + \nabla\cdot\left[ \mathbf v(\mathbf R)\, f \right] + \left[ E_L(\mathbf R) - E_T \right] f , \qquad \mathbf v = \nabla \ln |\Psi_T| . $$ I tre termini si interpretano come processi stocastici su una popolazione di walker: diffusione gaussiana di ampiezza \( \sqrt{\tau} \), drift lungo \( \mathbf v \) (che spinge i walker dove \( \Psi_T \) è grande), e branching con peso $$ w = \exp\!\left[ -\tau \left( \frac{E_L(\mathbf R) + E_L(\mathbf R')}{2} - E_T \right) \right] , $$ realizzato con \( \lfloor w + u \rfloor \) copie (\( u \) uniforme). Il riferimento \( E_T \) viene aggiornato per stabilizzare la popolazione (default 400 walker) e l'energia DMC è la media pesata di \( E_L \). Un'accettazione di Metropolis sulla funzione di Green approssimata corregge l'errore di passo temporale finito.

Fixed-node. Per i fermioni \( \Phi \) deve cambiare segno, ma \( f = \Psi_T \Phi \) deve restare positiva: si impone che \( \Phi \) abbia gli stessi nodi di \( \Psi_T \) [20, 7], rigettando le mosse che attraversano un nodo (il segno del determinante è tracciato tramite i pivot della LU). L'energia FN-DMC è esatta per i nodi dati e resta un limite superiore: l'errore residuo è solo l'errore nodale, come si vede confrontando H2O con nodi STO-3G (\( -76{,}311 \)) e 6-31G (\( -76{,}426 \pm 0{,}008 \), a 12 mHa dall'esatto \( -76{,}438 \)).

Stabilità: il drift limitato. Vicino ai nodi \( \mathbf v \sim 1/d \) diverge: un walker propone salti enormi, viene sempre rigettato e, restando fermo con \( E_L \) enorme, si clona fino a far collassare la popolazione (la patologia osservata ai \( \tau \) piccoli). La cura di Umrigar, Nightingale e Runge [8] sostituisce il drift con la versione limitata $$ \tilde{\mathbf v} = \mathbf v\; \frac{-1 + \sqrt{1 + 2\tau |\mathbf v|^2}}{\tau |\mathbf v|^2} , $$ che per \( \tau |\mathbf v|^2 \to 0 \) tende a \( \mathbf v \) e per \( |\mathbf v| \to \infty \) resta \( O(\sqrt{2/\tau}) \), usata coerentemente nella proposta e in entrambe le funzioni di Green; nel branching \( \tau \) è sostituito da \( \tau_{\text{eff}} = \tau \cdot \langle \text{accettazione} \rangle \) e le code di \( E_L \) sono tagliate a \( \pm 2/\sqrt{\tau} \). Dopo la cura, \( \tau = 0{,}005 \) e \( 0{,}0025 \) concordano (l'errore di time-step è sotto controllo) e l'accettazione cresce monotonamente al calare di \( \tau \).

11. Identità numerica cross-linguaggio e generatore di numeri casuali

La regola di ferro della suite — a parità di seed, risultati identici in C++, Python, Javascript e PHP — poggia sul generatore xorshift64* di Marsaglia e Vigna [14]: $$ s \mathrel{\oplus{=}} s \gg 12, \qquad s \mathrel{\oplus{=}} s \ll 25, \qquad s \mathrel{\oplus{=}} s \gg 27, \qquad \text{out} = s \times \text{0x2545F4914F6CDD1D} , $$ con l'uniforme in \( [0,1) \) costruita dai 53 bit alti divisi per \( 2^{53} \) (esattamente la mantissa di un double). Le insidie sono nei linguaggi senza interi a 64 bit nativi: Python maschera con & 0xFFFFFFFFFFFFFFFF a ogni passo; Javascript tiene lo stato come coppia di parole a 32 bit e fa la moltiplicazione a "limb" di 16 bit; PHP emula lo shift logico senza segno. Stesso ordine di estrazione ovunque ⇒ stessa traiettoria stocastica, fino all'MD5 della lista di walker. Nel motore parallelo i seed di catene e walker sono derivati con splitmix64 da (seed, passo, indice): il risultato è identico bit-per-bit con qualunque numero di thread — la stessa proprietà che rende riproducibile il futuro porting su GPU (schema counter-based alla Philox, come in QMCPACK [21]).

12. Bibliografia

  1. E. Schrödinger, Quantisierung als Eigenwertproblem (Erste Mitteilung), Annalen der Physik 384 (79), 361–376 (1926). doi:10.1002/andp.19263840404
  2. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller, Equation of State Calculations by Fast Computing Machines, J. Chem. Phys. 21, 1087–1092 (1953). doi:10.1063/1.1699114
  3. W. L. McMillan, Ground State of Liquid He⁴, Phys. Rev. 138, A442–A451 (1965). doi:10.1103/PhysRev.138.A442
  4. R. Jastrow, Many-Body Problem with Strong Forces, Phys. Rev. 98, 1479–1484 (1955). doi:10.1103/PhysRev.98.1479
  5. T. Kato, On the eigenfunctions of many-particle systems in quantum mechanics, Commun. Pure Appl. Math. 10, 151–177 (1957). doi:10.1002/cpa.3160100201
  6. D. M. Ceperley, B. J. Alder, Ground State of the Electron Gas by a Stochastic Method, Phys. Rev. Lett. 45, 566–569 (1980). doi:10.1103/PhysRevLett.45.566
  7. P. J. Reynolds, D. M. Ceperley, B. J. Alder, W. A. Lester Jr., Fixed-node quantum Monte Carlo for molecules, J. Chem. Phys. 77, 5593–5603 (1982). doi:10.1063/1.443766
  8. C. J. Umrigar, M. P. Nightingale, K. J. Runge, A diffusion Monte Carlo algorithm with very small time-step errors, J. Chem. Phys. 99, 2865–2890 (1993). doi:10.1063/1.465195
  9. M. Taut, Two electrons in an external oscillator potential: Particular analytic solutions of a Coulomb correlation problem, Phys. Rev. A 48, 3561–3566 (1993). doi:10.1103/PhysRevA.48.3561
  10. W. M. C. Foulkes, L. Mitas, R. J. Needs, G. Rajagopal, Quantum Monte Carlo simulations of solids, Rev. Mod. Phys. 73, 33–83 (2001). doi:10.1103/RevModPhys.73.33
  11. A. Ma, M. D. Towler, N. D. Drummond, R. J. Needs, Scheme for adding electron–nucleus cusps to Gaussian orbitals, J. Chem. Phys. 122, 224322 (2005). doi:10.1063/1.1940588
  12. J. M. Thijssen, Computational Physics, 2ª ed., Cambridge University Press (2007). ISBN 978-0-521-83346-2.
  13. B. L. Hammond, W. A. Lester Jr., P. J. Reynolds, Monte Carlo Methods in Ab Initio Quantum Chemistry, World Scientific, Singapore (1994). doi:10.1142/1170
  14. G. Marsaglia, Xorshift RNGs, J. Stat. Softw. 8(14), 1–6 (2003). doi:10.18637/jss.v008.i14; S. Vigna, An experimental exploration of Marsaglia's xorshift generators, scrambled, ACM Trans. Math. Softw. 42(4), 30 (2016). doi:10.1145/2845077
  15. C. L. Pekeris, 1¹S and 2³S States of Helium, Phys. Rev. 115, 1216–1221 (1959). doi:10.1103/PhysRev.115.1216
  16. W. Kołos, L. Wolniewicz, Potential-Energy Curves for the X¹Σg⁺, b³Σu⁺, and C¹Πu States of the Hydrogen Molecule, J. Chem. Phys. 43, 2429–2441 (1965). doi:10.1063/1.1697142
  17. W. Cencek, J. Rychlewski, R. Jaquet, W. Kutzelnigg, Sub-microhartree accuracy potential energy surface for H₃⁺ … I. Calculation of the potential points, J. Chem. Phys. 108, 2831–2836 (1998). doi:10.1063/1.475702
  18. J. Sherman, W. J. Morrison, Adjustment of an Inverse Matrix Corresponding to a Change in One Element of a Given Matrix, Ann. Math. Statist. 21, 124–127 (1950). doi:10.1214/aoms/1177729893
  19. Q. Sun et al., PySCF: the Python-based simulations of chemistry framework, WIREs Comput. Mol. Sci. 8, e1340 (2018). doi:10.1002/wcms.1340
  20. J. B. Anderson, A random-walk simulation of the Schrödinger equation: H₃⁺, J. Chem. Phys. 63, 1499–1503 (1975). doi:10.1063/1.431514
  21. J. Kim et al., QMCPACK: an open source ab initio quantum Monte Carlo package…, J. Phys.: Condens. Matter 30, 195901 (2018). doi:10.1088/1361-648X/aab9c3

Nota sulla numerazione: i richiami nel testo usano [1]=Schrödinger, [2]=Metropolis, [3]=McMillan, [4]=Jastrow, [5]=Kato, [6]=Ceperley–Alder, [7]=Reynolds et al., [8]=Umrigar et al., [9]=Taut, [10]=Foulkes et al., [11]=Ma et al., [12]=Thijssen, [13]=Hammond et al., [14]=Marsaglia/Vigna, [15]=Pekeris, [16]=Kołos–Wolniewicz, [17]=Cencek et al., [18]=Sherman–Morrison, [19]=PySCF, [20]=Anderson, [21]=QMCPACK.

Sistema e parametri
bohr
bohr
bohr
Risultati

Nessun calcolo eseguito.

Nuvola di walker

Ogni punto è la posizione di un elettrone in una configurazione campionata da |ψ|² (colore per elettrone; sfere chiare = nuclei). Trascina per ruotare, rotella per lo zoom.

Confronto di prestazioni

Tempi di calcolo normalizzati a 1 milione di campioni (secondi), aggiornati a ogni esecuzione: stessa traiettoria stocastica, cambia solo il linguaggio.

Esegui lo stesso sistema con backend diversi per popolare la tabella.

Visualizzatore dell'atomo di idrogeno: sovrapposizioni di autostati |n,l,m〉 (fino a n=5), evoluzione temporale tramite le fasi e−iEt, campionamento Metropolis della densità |ψ|² con il colore che codifica la fase. È lo stesso campionamento che, reso quantitativo, diventa il VMC della prima scheda.

Idrogeno / evoluzione unitaria

t = 0.000 fs
Stato A · |n ℓ m⟩
Stato B · |n ℓ m⟩
Miscela |c_B|²50%
Velocità1.0×
Punti
Dimensione punto1.0
E_A = · E_B =
Periodo battimento T =
densità stazionaria (ΔE = 0)
±π 0
arg ψ

Keywords: Schrödinger, Monte Carlo, VMC, quantistica, benchmark