Forme cyclique du butadiyne

Réactivité en milieu atmosphérique et analyse Monte Carlo quantique de la localisation électronique

L’optimisation de la fonction d’onde par la minimisation de la variance

 La variance de l’énergie locale, définie par σ 2 [ΨT ] = Z h HbΨT (R) i2 dR − Z ΨT (R)HbΨT (R) dR 2 (1.72) est un estimateur de la qualité de la fonction d’onde. Plus la fonction d’onde est loin d’un état propre, plus la variance de l’énergie locale est grande. Donc, pour la fonction d’onde exacte la variance est nulle. Dans les calculs utilisant les méthodes traditionnelles (Hartree-Fock, DFT,. . .), optimiser la fonction d’onde consiste à faire varier des paramètres (coefficients des orbitales moléculaires, poids des configurations,. . .) dans le but de minimiser l’énergie totale du système. Avec les méthodes QMC, on peut considérer la variance de l’énergie locale comme un estimateur de la qualité de la fonction d’onde.

Donc l’optimisation des paramètres variationnels peut se faire en réduisant le plus possible la variance de l’énergie locale vers sa limite inférieure, zéro. Nature de la fonction d’onde d’essai Dans les calculs qui suivent, nous avons utilisé des fonctions d’onde HartreeFock sur des bases de fonctions gaussiennes, calculées à l’aide du programme GAMESS. à ces fonctions d’essai, on affecte un facteur de Jastrow qui prend en compte la corrélation : ΨT = DαDβe J (1.73) Dα et Dβ sont les déterminants de Slater de la fonction Hartree-Fock, et e J est le facteur de Jastrow. Une forme minimale typique pour le facteur de Jastrow est la suivante : J = N Xelec i N Xelec j>i C(rij) − N Xelec i N Xnucl M PM(riM) (1.74) C(rij ) est le terme de cusp électron-électron : C(rij ) = aσ rij 1 + bσrij (1.75) Le couple (aσ, bσ) prend la valeur (a↑↑, b↑↑) si les électrons i et j sont de même spin (αα ou ββ), sinon il prend la valeur (a↑↓, b↑↓) si les électrons i et j sont de spins opposés (αβ ou βα). Pim est le terme de pénétration électron-noyau qui permet essentiellement de corriger la densité ρ(r) : PM (rim) = pMriM (1.76) Nous avons également considéré des formes plus sophistiquées incluant des développements systématiques en puissances de rij et de riM, ainsi que des termes à trois corps couplant rij et riM.

Nous n’entrerons pas ici dans les détails[16]. Ainsi, les paramètres variationnels sont les coefficients des orbitales moléculaires, et les paramètres du facteur de Jastrow. En pratique, les coefficients des orbitales moléculaires ne sont pas optimisés, à cause de leur nombre important (contrôle du calcul plus délicat). On dispose d’un nombre constant de marcheurs obtenus d’après l’échantillonnage de la densité résultant de la fonction d’essai, c’est-a-dire de la fonction d’onde Hartree-Fock : les paramètres du facteur de Jastrow (cusp et pénétration) sont tous nuls. Au cours de l’optimisation, on ne déplace pas les marcheurs.

Le calcul de l’énergie et de la variance est toujours effectué sur le même jeu de marcheurs. On utilise un optimiseur de type gradient conjugué ou simplexe pour minimiser la variance en faisant varier les paramètres du facteur de Jastrow. Lorsque le minimiseur a convergé, les paramètres obtenus sont ceux qui minimisent la variance pour un jeu de configurations qui échantillonne bien la densité de la fonction d’essai, mais moins bien celle de la fonction optimisée. Il est donc nécessaire de réaliser une simulation avec la nouvelle fonction d’onde, puis de l’optimiser à nouveau, jusqu’a ce que la variance et l’énergie ne varient plus. Cette méthode, dite de l’échantillonnage corrélé[17] ou Monte Carlo différentiel, nous permet d’une part de réduire énormément le temps de calcul en évitant de refaire une simulation pour chaque calcul de l’énergie, et d’autre part, grˆace à la corrélation entre les marcheurs, la différence hE (Ψi) L i − hE (Ψj ) L i a une variance plus faible que si les énergies étaient calculées à partir de deux jeux de marcheurs indépendants, ce qui guide mieux les paramètres vers le minimum de la variance. Au préalable, nous avons ajouté des fonctions de Slater à l’ensemble des fonctions base, initialement constitué de fonctions gaussiennes (base mixte STO et CHAPITRE 1. METHODES ´ 25 GTO). Les orbitales moléculaires (OM) 1s de cœur autour de chaque atome sont décomposées sur la base STO de Koga et al, et les autres OM sont exprimées sur la base GTO.

En effet, les fonctions gaussiennes n’ont pas de cusp en leur centre, donc toute la région proche des noyaux est très mal décrite par celles-ci. Ce choix de fonctions nous a permis de réduire la variance d’un facteur 5 avant de commencer l’optimisation du facteur de Jastrow. 1.4.5 Diffusion Monte Carlo (DMC) Dans cette partie, nous présentons la méthode DMC avec nœuds fixés. Cette méthode est la méthode de référence pour effectuer des calculs au-dela de la méthode variationnelle (VMC). Elle permet d’éliminer presque entièrement la dépendance des calculs d’énergie par rapport à la fonction d’essai, l’erreur résiduelle étant uniquement due à sa structure nodale approchée. 

Méthode de projection 

On commence par écrire l’équation de Schrodinger en temps imaginaire − ∂Φ(R,t) ∂t =  − 1 2 ∇ 2 + V (R) − ET  Φ(R,t) (1.77) ou V (R) est l’énergie potentielle de la molécule et ET est une estimation de l’énergie de l’état fondamental. En exprimant la fonction d’onde Φ(R,t) sur une base complète de fonctions propres de l’opérateur hamiltonien, Φ(R,t) = X i ci(t)Φi(R) (1.78) avec ci(t) = hΦi(R)|Φ(R, 0)ie −t(Ei−ET ) (1.79) on a : Φ(R,t) = X i cie −(Ei−ET )tΦi(R) (1.80) ou Ei sont les valeurs propres correspondant aux états propres Φi(R). Pour un temps imaginaire t suffisamment long, parmi les termes avec un coefficient ci non-nul, seul celui de plus basse énergie subsiste, tous les autres termes  tendent vers zéro. La solution asymptotique de l’équation (1.77) est Φ(R,t) = c0e −(E0−ET )tΦ0(R) (1.81) Si ET est ajusté comme la valeur de l’énergie exacte de l’état fondamental, la solution asymptotique est l’état fondamental exact. En revanche, si la fonction d’onde d’essai est orthogonale à Φ0(R), on a c0 = 0, et la solution asymptotique devient l’état excité de plus basse énergie avec un coefficient ci non-nul.

 Interprétation en terme de diffusion/branchement 

On peut réécrire l’équation (1.77) comme ∂Φ(R,t) ∂t = D∇ 2Φ(R,t) − [V (R) − ET ] Φ(R,t) (1.82) Le premier terme de cette équation peut être assimilé à une équation de diffusion dans l’espace des 3N coordonnées ou Φ(R,t) joue le rôle de la densité de particules qui diffusent, et le facteur D = 1/2 devant l’opérateur laplacien correspond au coefficient de diffusion. Une telle équation peut être facilement simulée avec une marche aléatoire des particules dans l’espace des configurations. Le second terme ressemble, quant à lui, à une équation de vitesse décrivant un processus de mort et de naissance d’individus dans une population. Ainsi, l’équation (1.82) peut être simulée par la combinaison d’un processus de diffusion et d’un processus de branchement, dans lequel le nombre de particules diffusant augmente ou diminue de façon à réduire la densité de probabilité dans les régions ou V (R) est grand et à l’augmenter dans les zones d’énergie potentielle favorable. 

L’approximation des nœuds fixés

 Pour que l’interprétation en terme de diffusion reste valide, il faut que Φ(R,t), qui représente une densité de particules dans l’équation (1.82), reste d’un signe constant (choisi positif) dans tout l’espace. Cette méthode paraît donc restreinte aux systèmes bosoniques dont la fonction d’onde fondamentale n’admet pas de zéros (nœuds). Compte tenu du caractère antisymétrique de la fonction d’onde des fermions, la fonction d’onde est nécessairement négative dans certaines régions de l’espace à 3N dimensions. La méthode des nœuds fixés consiste à séparer la simulation des régions positives et négatives de la fonction d’essai et à interdire la diffusion à travers une hypersurface nodale. On fait une approximation, appelée approximation des nœuds fixés (FixedNode DMC), qui consiste à effectuer la simulation DMC en utilisant les nœuds de la fonction d’essai. Cette approximation nous donne la meilleure énergie possible pour l’hypersurface nodale de la fonction d’essai. On peut donc voir la méthode FN-DMC comme une méthode variationnelle qui donne la valeur exacte de l’énergie si l’hypersurface nodale est exacte. Dans la majorité des cas, l’hypersurface nodale de la fonction d’essai (∼ Hartree-Fock) est peu différente celle de la fonction exacte, donc les énergies obtenues sont très proches des valeurs exactes. 

Reformulation avec la fonction de Green

 L’équation de Schrodinger en temps imaginaire est la suivante : − ∂ ∂t Φ(R,t) = (Hb − ET )Φ(R,t) (1.83) ou t est une variable réelle qui mesure la propagation en temps imaginaire. On cherche une solution de l’équation de Schr¨odinger sous forme intégrale : Φ(R,t + τ ) = Z G(R0 → R, τ )Φ(R0 ,t) dR0 (1.84) ou G(R0 → R, τ ) est une fonction de Green. L’équation à résoudre est donc − ∂ ∂t G(R0 → R,t) = (Hb − ET )G(R0 → R,t) (1.85) et l’on trouve G(R0 → R, τ ) = e −τHb = X i |Φiie −τEi hΦi | (1.86) {Φi} et {Ei} dénotent les ensembles complets des fonctions propres et des valeurs propres de Hb. Lorsque τ → ∞, l’opérateur e −τ (Hb−ET ) projette l’état fondamental |Φ0i qui a un recouvrement non-nul avec la fonction d’onde d’essai choisie. limτ→∞ hR|e −τ (Hb−ET ) |Φ(R, 0)i = limτ→∞ Z G(R0 → R, τ )Φ(R, 0)(R0 ) dR0 = limτ→∞ X i |Φi(R)ie −τ (Ei−ET ) hΦi |Φ(R, 0)i = limτ→∞ Φ0(R)e −τ (E0−ET ) hΦ0|Φ(R, 0)i (1.87) En ajustant ET égal à E0, on rend constant le facteur exponentiel et on fait disparaître les termes supérieurs, car leur énergie est supérieure à E0.

Considérons l’hamiltonien complet qui inclut les termes d’énergie cinétique et potentielle. Lorsque les particules interagissent, dans la majorité des cas on ne connaît pas l’expression exacte de la fonction de Green, et il faut avoir recours à des approximations, comme la formule de Trotter-Suzuki. agit comme une renormalisation de la fonction de Green en fonction du temps. On utilise l’algorithme de branchement mort/naissance ou P détermine le nombre de marcheurs qui survivent jusqu’au pas suivant[22] selon la règle suivante : le nombre de marcheurs au pas suivant est M = E[P + η] ou E est la partie entière et η un nombre aléatoire entre 0 et 1. D’après l’équation (1.91), on remarque que les marcheurs vont être multipliés dans les régions d’énergie potentielle favorable, et qu’ils vont disparaître dans les régions d’énergie potentielle défavorable. L’énergie de référence ET détermine la renormalisation asymptotique et contrôle donc la population totale des marcheurs.