Calcul du facteur de fractionnement isotopique 

Télécharger le fichier pdf d’un mémoire de fin d’études

Equilibre thermodynamique

L’equilibre thermodynamique entre deux systemes A et B peut donner lieu a un fractionnement des isotopes correspondant a une modification de la repartition des isotopes entre les deux systemes.
Constante d’echange Le sch´ema 1.2 repr´esente un equilibre entre deux phases s´epar´ees par une surface d’´echange. On consid`ere l’exemple d’un ´equilibre entre un min´eral compos´e de deux isotopes et un liquide compos´e de mol´ecules diatomiques, form´ees par les mˆemes isotopes que le min´eral. Les isotopes peuvent soit ˆetre echang´es au sein d’un mˆeme syst`eme (pas de fractionnement) soit entre les deux syst`emes (ce qui donne lieu au fractionnement isotopique). Bien que cet ´equilibre ne soit pas une r´eaction chimique, on peut lui associer une ´equation d’´equilibre de la forme : YA+YB∗ =YB+YA∗ (1.8).
o`u YA repr´esente l’isotope Y dans la phase A, Y∗ est un isotope rare de Y. Cet ´equilibre peut alors ˆetre d´ecrit par une constante d’´equilibre qui s’´ecrit sous la forme : Kech = [YA][YB∗][YB][Y ∗] A (1.9).
Si le syst`eme n’est pas a` l’´equilibre, alors le rapport de concentration, que l’on appelle quotient de r´eaction et qui est not´e QR, sera diff´erent de Kech. Le flux d’´echange d’isotopes est traduit par la comparaison entre ces deux rapports. Par exemple, si QR <Kech alors la r´eaction (1.8) va dans le sens direct.

Methodes pour calculer le fractionnement isotopique

Pour calculer le fractionnement isotopique, il faut d´eterminer la variation d’´energie libre associ´ee a` un changement de masse dans deux syst`emes a` l’´equilibre. Comme le fractionnement isotopique est un ph´enom`ene quantique, il faut que l’aspect quantique des noyaux soit pris en compte.
En couplant la dynamique mol´eculaire des int´egrales de chemins (PIMD, voir sec. 2.4) avec la methode de l’Int´egration thermodynamique (voir section 1.4.2.b), il est possible d’avoir acc`es `a la variation d’´energie libre d’un syst`eme, pour un changement de masse. Cette m´ethode permet de prendre en compte les aspects quantiques des noyaux sans n´egliger l’anharmonicit´e. Elle est appel´ee TI-PIMD et sera pr´esent´ee dans le chapitre 2.5.1. Cependant, cette m´ethode requiert beaucoup de ressources en terme de temps de calcul.
Une autre m´ethode consiste `a calculer le spectre vibrationnel lorsque l’approxima-tion harmonique est valable. Cela permet de remonter `a l’information sur l’aspect quantique du noyau, et donc de calculer le fractionnement isotopique. Cepen-dant, en n´egligeant l’anharmonicit´e on perd une partie de cette information car les niveaux d’´energie sont modifi´es lorsque l’on se place dans l’approximation har-monique. L’erreur que l’on commet en se pla¸cant dans l’approximation harmonique d´epends des syst`emes etudi´es.

Approximation harmonique

La methode la plus couramment utilis´ee pour calculer le fractionnement isotopique consiste a` se placer dans l’approximation harmonique. On consid`ere que cette approximation est valable si les atomes sont dans leur position d’´equilibre, c’est-a`-dire que le syst`eme est a` un minimum d’energie (voir figure 1.3). Dans cette approximation, il est possible de calculer les fr´equences des modes de phonon d’un syst`eme. Notons que si tous les atomes du syst`eme ne sont pas dans leur minimum d’´energie, on peut obtenir des fr´equences n´egatives, correspondant a` des modes imaginaires, signalant l’instabilit´e du syst`eme. Dans cette section, nous allons montrer comment il est possible de prendre en compte l’aspect quantique des noyaux a` partir du spectre vibrationnel en calculant la fonction de partition du syst`eme quantique qui est reli´ee au β-facteur par l’´equation (1.23).
Dans l’approximation harmonique, dans ce travail, nous avons uniquement etudi´ des syst`emes assimil´es `a des solides cristallins (soit des min´eraux, soit des con-figurations de liquides relax´ees `a 0 K). Notons qu’il est possible d’utiliser cette m´ethode pour d’autres syst`emes tels que des mol´ecules isol´ees. Pour parvenir a` calculer le fractionnement, et donc la variation d’´energie libre, on consid`ere un syst`eme cristallin qui est a` son minimum d’´energie. Dans un syst`eme cristallin, seuls les ´etats vibrationnels contribuent au fractionnement isotopique (Bottinga, 1968). On peut exprimer le rapport des fonctions de partitions intervenant dans l’´equation du β-facteur (eq. (1.23)) : Z∗ = Z∗ = Nq 3N e− µi,q /2 (1 − e−µi,q∗ ) 1/Nq (1.31)
Z Z vib µi,q∗ /2 (1 µi,q )q=1 i=1 e− e− Y Y − o`u µi = hνi/kBT , avec νi les fr´equences des modes de vibration du cristal (phonons), Nq est le nombre de points q utilis´e pour ´echantillonner la zone de Brillouin et N est le nombre d’atomes dans la maille el´ementaire. Dans le produit sur i, les trois premiers modes correspondent aux branches acoustiques. Dans la suite, les fr´equences `a des modes de translation du syst`eme sont implicitement ignor´ees puisqu’elles sont th´eoriquement nulles.
En rempla¸cant (1.31) dans l’´equation (1.23), le β-facteur s’´ecrit sous la forme : βYA∗ = 3N e− µi,q /2 (1 − e−µi,q∗ ) 1/Nq (1.32).
Dans l’approximation harmonique, il est donc possible de calculer la fonction de partition quantique du syst`eme a` partir du spectre vibrationnel calcul´e sur le syst`eme de type solide cristallin. Pour cela, on peut utiliser les m´ethodes pr´esent´ees dans la section 1.4.2.a.

R´esum´e de la PIMD

La PIMD permet, `a partir de trajectoires classiques, de reproduire les propri´et´es d’un syst`eme dont les noyaux sont quantiques. Puisqu’il s’agit d’une dynamique mol´eculaire, on est `a temp´erature finie et on prend en compte l’anharmonicit´e du syst`eme (via la forme des potentiels empiriques ou par le calcul de la densit´e ´electronique par des m´ethodes ab initio). Dans cette m´ethode, il est possible de parall´eliser, entre chaque pas, le calcul des forces pour chacune des P r´epliques du syst`eme ´evoluant ind´ependamment. Dans un second temps, on ajoute les oscillateurs harmoniques pour faire ´evoluer la position des noyaux dans chacune des r´epliques. Il est possible de coupler la PIMD soit avec des potentiels empiriques soit avec des m´ethodes ab initio telles que la CPMD. Pour plus de d´etails, on peut se r´ef´erer au travail de Marx et Parrinello (1996) qui a fond´e l’id´ee du couplage entre ces deux m´ethodes.

Calcul du fractionnement

Calcul du facteur de fractionnement avec la m´ethode TI-PIMD

Int´egration thermodynamique L’int´egration thermodynamique (Thermodynamics Integration, TI) permet de calculer la variation d’´energie libre associ´ee `a un changement de masse (voir sec. 1.4.2.b). Dans cette section, nous allons montrer comment coupler cette m´ethode avec les PIMD dans le but de calculer le fractionnement isotopique.

Choix des masses

Pour r´ealiser l’int´egration thermodynamique, nous avons vu dans la section 1.4.2.b qu’il ´etait n´ecessaire d’´echantillonner dF/d, pour plusieurs valeurs de , o`u est un param`etre permettant de faire ´evoluer le syst`eme. L’Integration Thermodynamique peut ˆetre r´ealis´ee sur un chemin quelconque pris par entre 0 et 1, voir eq. (1.37). Par exemple pour un changement de masse entre deux isotopes, il est possible de choisir une variation lin´eaire telle que : M() = Ma (1 − ) +Mb (2.121).

Calcul du facteur de fractionnement dans l’approximation harmonique

Comme on l’a vu en section (voir sec. 1.4.2.a), en calculant le spectre de vibrations du syst`eme, dans l’approximation harmonique, il est possible de calculer la fonction de partition quantique du syst`eme et donc la variation d’´energie libre pour un changement de masse. Pour cela, il est necessaire que le syst`eme soit dans un minimum d’´energie.
Dans un cristal form´e d’un ensemble de N cellules unitaires, on peut ´ecrire les forces agissant sur un atome comme : F(Il) = − @V @R(Il) (2.141) o`u est un indice cart´esien indiquant la direction de la force, I est l’indice de l’atome dans une cellule et l correspondant `a l’indice de la cellule o`u est plac´e l’atome. `A partir d’un ensemble de vecteurs forces, on peut calculer la matrice des constantes de forces : (Il, Jl0) = − @F(Jl0)  @R(Il) (2.142) o`u est un autre axe cart´esien, J est la position d’un autre atome dans la cellule l’.
Les coefficients de la matrice correspondent donc `a la variation de la force agissant sur un atome J ayant ´et´e d´eplac´e selon l’axe par rapport `a une variation d’un autre atome I se d´epla¸cant dans une autre direction . Elle peut ˆetre calcul´ee en utilisant la m´ethode des d´eplacements finis (voir sec. 2.5.2) ou la m´ethode DFPT (voir sec. 2.5.2). `A partir de la matrice des constantes de forces, on peut calculer la matrice dynamique : D(IJ, q) = 1 p MIMJ X l0 (I0, Jl0) exp(iq · [R(Jl0) − R(I0)]) (2.143).

Parametres des simulations

Details des parametres pour les liquides

Deux liquides — l’un fait de 60 mol´ecules d’eau et une mol´ecule de H4SiO4 et l’autre fait de 60 mol´ecules d’eau et une mol´ecule de H3SiO−4 — ont ´et´e simul´es par MD ab initio en utilisant la m´ethode de Car-Parrinello (CPMD, Car et Parrinello (1985a)). Les d´etails des param`etres utilis´es pour les calculs DFT utilis´es en CPMD (cutoff en energie, fonctionnelle d’´echange corr´elation, pseudopotentiels et ´echantillonnage de la zone de Brillouin) sont donn´es plus loin. Nous sommes partis d’une configuration de 64 mol´ecules d’eau, ´equilibr´ee en utilisant une dynamique mo´eculaire classique `a 300 K dans une boˆıte cubique d’arrˆete 12.41°A(Lee et Tuckerman, 2006), de densit´e 1. Dans cette boˆıte, 4 mol´ecules d’eau ont ´et´e remplac ´ees par l’esp`ece H4 (resp. H3), la densit´e de ces deux syst`emes est 1.020 (resp. 1.019). Avant de d´emarrer la dynamique mol´eculaire ab initio, nous avons r´ealis´e une optimisation g´eom´etrique pour minimiser les forces agissant sur les atomes, dues au remplacement de 4 mol´ecules d’eau par une mol´ecule de taille moyenne.
En ce qui concerne la dynamique, le pas de temps choisi est de 0.12 fs pour toutes les dynamiques CPMD. Pour le syst`eme H3, qui est charg´e, nous avons utilis´e une charge de fond uniforme pour assurer la neutralit´e du syst`eme. Pour chaque syst`eme, la proc´edure suivante a ´et´e utilis´ee : premi`erement, une br`eve dynamique mol´eculaire est r´ealis´ee dans l’ensemble (N, V, E) durant 1 ps jusqu’`a ce qu’une temp´erature stable soit atteinte. Ensuite, le syst`eme est doucement amen´e `a la temp´erature de notre ´etude (300 K) durant 2 ps en incr´ementant la temp´erature `a chaque pas. Comme cela a ´et´e pr´esent´e dans le chapitre 2, les thermostats utilis´es pour contrˆoler la temp´erature sont des thermostats de Nos´e-Hoover massifs (Martyna et al., 1992). La fr´equence des thermostats a ´et´e choisie `a 3000 cm−1 et leur masse fictive est de 800 a.u. Les deux liquides ont ´et´e ´equilibr´es durant 15 ps dans l’ensemble (N, V, T) et les donn´ees ont ´et´e collect´ees durant les 40 ps suivantes de la trajectoire. Toutes les 10 ps, les ´electrons ont ´et´e refroidis `a 0 K pour atteindre la surface d’´energie de Born-Oppenheimer et assurer que les orbitales de Kohn et Sham restent proches de celles de l’´etat fondamental. Dans le but d’acc`eder `a la pression dynamique des liquides, d’autres trajectoires ind´ependantes ont ´et´e r´ealis´ees `a diff´erents cutoffs en ´energie (80, 140 et 160 Ry) et le tenseur des contraintes a ´et´e calcul´e. La pression dans les liquides a ´et´e calcul´ee comme la valeur moyenne de la trace de ce tenseur. En particulier, ces r´esultats seront utilis´es pour ´etudier l’erreur due `a la diff´erence de pression dans les calculs entre le liquide et le solide (voir section 3.3.3.c).
D´etails des param`etres des calculs des modes de vibration Les -facteurs ont ´et´e obtenus `a partir des ´equations (1.32) et (1.35) en utilisant les fr´equences phonon calcul´ees par la m´ethode DFPT (Baroni et al., 2001). Pour les solides, ces fr´equences de phonons sont calcul´ees sur la structure cristalline relax´ee `a 0 K et `a pression nulle. Pour les liquides, des configurations ont ´et´e extraites et les positions des atomes ont ´et´e relax´ees `a 0 K puis les phonons ont ´et´e calcul´es sur ces structures dite “structures inh´erentes” (Inherent Structures, IS).
Nous avons utilis´e l’approximation GGA (generalized-gradient approximation) et la fonctionnelle d’´echange-corr´elation de Becke, Lee et Par (BLYP, Becke, 1988; Lee et al., 1988). La fonctionnelle BLYP est connue pour donner une bonne description de l’eau liquide (structure et dynamique, Sprik et al. (1996); Lin et al. (2012)).
Comme cela est montr´e par ces auteurs, la fonctionnelle PBE donne un moins bon accord avec les mesures exp´erimentales de la structure de l’eau liquide. Pour les fr´equences de vibration, Silvestrelli et al. (1997) et Lee et Tuckerman (2007) ont montr´e la capacit´e de la fonctionnelle BLYP `a reproduire de fa¸con correcte le spectre infrarouge de l’eau. Cependant, pour l’´etude des silicates, par exemple les silicates hydrat´es (M´eheut et al., 2007), les liquides silicat´es (P¨ohlmann et al., 2004), ou mˆeme pour les esp`eces du silicium dissoutes (Spiekermann et al., 2012a,b), la fonctionnelle PBE est souvent pr´ef´er´ee. Comme cela est discut´e dans les sections 3.3.2 (Table 3.4) et 3.3.3.b, la BLYP sur´estime les distances Si-O davantage que la PBE par rapport aux valeurs exp´erimentales. N´eanmoins, le spectre de vibrations harmoniques (Figure 3.15) et les propri´et´es de fractionnement calcul´ees pour des solides avec les fonctionnelles BLYP et PBE sont en tr`es bon accord (Figure 3.16). Afin d’´etudier le lien entre la structure des liquides et les propri´et´es de fractionnement, la fonctionnelle BLYP a ´et´e pr´ef´er´ee `a la PBE pour cette ´etude.
Les ´electrons de coeur sont d´ecris par des pseudopotentiels de type norm conserving (Troullier et Martins, 1991) selon l’´equation de Kleinman-Bylander (Kleinman et Bylander, 1982). Les d´etails de simulation pour le quartz et la kaolinite sont ceux utilis´es par M´eheut et al. (2007). Pour les liquides, les fonctions d’ondes´electroniques sont projet´ees sur une base d’ondes planes en utilisant une ´energie de cutoff de 80 Ry, et le cutoff pour la densit´e de charge est de 4Ecut. Notons que seul le point ? est utilis´=e pour l’´echantillonnage de la zone de Brillouin, la dispersion est donc negligee.

Table des matières

Les isotopes stables en geologie
Historique et recents developpements des methodes de calcul du fractionnement isotopique
Objectifs et notions de calculs realistes
Organisation du manuscrit
Chapitre 1. Calcul du facteur de fractionnement isotopique 
1.1. Introduction
1.2. Cause du fractionnement isotopique
1.3. Equilibre thermodynamique
1.4. Calcul du facteur de fractionnement
1.4.1. Notations
1.4.2. M´ethodes pour calculer le fractionnement isotopique
1.4.2.a. Approximation harmonique
1.4.2.b. Integration Thermodynamique (TI)
1.5. Temps de calcul
1.6. Conclusion
Chapitre 2. M´ethodologie 
2.1. Introduction
2.2. Description des interactions
2.2.1. Potentiels empiriques
2.2.1.a. Potentiels d’interactions pour Li2O
2.2.1.b. Potentiels d’interactions pour le lithium en solution
2.2.2. Calcul de la structure ´electronique
2.2.2.a. Hamiltonien du syst`eme ´electronique
2.2.2.b. Matrice densit´e
2.2.2.c. Th´eorie de la fonctionnelle de la densit´e
2.2.2.d. Calcul des forces
2.2.2.e. ´Energie d’´echange et corr´elation
2.2.2.f. Impl´ementation
2.2.2.g. Pseudo-potentiels utilis´es dans cette ´etude
2.2.2.h. Pseudo-potentiel du lithium
2.3. Dynamique mol´eculaire
2.3.1. Thermostats
2.3.2. Dynamique mol´eculaire empirique
2.3.3. Dynamique mol´eculaire ab initio
2.3.3.a. CPMD
2.3.3.b. R´esum´e
2.4. Dynamique mol´eculaire des Int´egrales de Chemin (Path Integrals Molecular Dynamics, PIMD)
2.4.1. Int´egrales de Chemin
2.4.2. PIMD
2.4.3. Difficult´es d’´echantillonnage en PIMD et changement de variables
2.4.4. Formation des anneaux
2.4.5. R´esum´e de la PIMD
2.5. Calcul du fractionnement
2.5.1. Calcul du facteur de fractionnement avec la m´ethode TI-PIMD
2.5.1.a. Choix des masses
2.5.1.b. Estimateurs pour la variation d’´energie libre
2.5.2. Calcul du facteur de fractionnement dans l’approximation harmonique
2.6. Codes de simulation
2.7. Conclusion
Chapitre 3. Calcul du fractionnement du silicium entre deux especes en solution 
3.1. Introduction
3.2. Param`etres des simulations
3.2.1. D´etails des param`etres pour les liquides
3.3. R´esultats
3.3.1. Propri´et´es structurales et dynamiques des liquides
3.3.1.a. Fonctions de distribution radiale
3.3.1.b. D´eplacement carr´e moyen
3.3.1.c. Couches de solvatation
3.3.2. Propri´et´es structurales et vibrationnelles du quartz et de la kaolinite
3.3.3. Incertitudes et -facteurs
3.3.3.a. Tests statistiques et valeur p
3.3.3.b. -facteurs
3.3.3.c. Incertitudes
3.3.3.d. Modes de phonon partiels
3.4. Discussion
3.4.1. Dispersion statistique et corr´elation entre la structure et les propri´et´es de fractionnement
3.4.1.a. La probl´ematique du d´esordre configurationnel
3.4.1.b. Utilisation de la corr´elation entre ¯dSiO et ln 30Si
3.4.1.c. Relation entre les distances Si-O et la premi`ere couche de solvatation
3.4.1.d. Vibrational analysis
3.4.2. Implications en g´eologie
3.5. Conclusion
Chapitre 4. Calcul du fractionnement du lithium entre une espece en solution et un mineral 
4.1. Introduction
4.2. Param`etres des simulations
4.3. R´esultats
4.3.1. Structure Li2O et polylithionite
4.3.2. Spectre vibrationnel de Li2O
4.3.3. Fonctions de distribution radiales du liquide
4.3.4. Notations
4.3.5. Calcul des estimateurs en fonction du temps
4.3.6. Convergence des estimateurs
4.3.6.a. Convergence temporelle
4.3.6.b. Convergence en fonction du nombre de r´epliques
4.3.6.c. ´Echantillonnage de la masse
4.3.6.d. Effets de taille et de dispersion
4.3.7. Calcul des propri´et´es de fractionnement isotopique
4.3.7.a. Liquides dans l’approximation harmonique
4.3.7.b. Solides dans l’approximation harmonique
4.3.7.c. Propri´et´es de partage isotopique
4.4. Discussion
4.4.1. Incertitude de la m´ethode TI-PIMD
4.4.2. Anharmonicite de la liaison Li-O
4.4.3. Facteur de fractionnement min´eral-Li+ calcule et experimental
4.4.4. Perspectives
4.5. Conclusion
Conclusion generale 

Télécharger le rapport complet