# Calculs ab initio : estimation d'erreur en modèle a posteriori

Denis Aubry, Amélie Fau, Ann-Lenaïg Hamon

► To cite this version:

Denis Aubry, Amélie Fau, Ann-Lenaïg Hamon. Calculs ab initio : estimation d'erreur en modèle a posteriori. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01408466

HAL Id: hal-01408466

https://hal.science/hal-01408466v1

Submitted on 5 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Calculs ab initio :
estimation d'erreur en modèle a posteriori

D. Aubry, A. Fau, A.-L. Hamon

Laboratoire MSS-Mat (UMR CNRS 8579)

École Centrale Paris, Grande Voie des Vignes, 92290 Châtenay-Malabry

{denis.aubry, amelie.fau, ann-lenaig.hamon}@ecp.fr

**Résumé** – Les calculs *ab initio* consistent à résoudre l'équation de Schrödinger via différentes approximations. On développe ici quelques modèles hiérarchisés. Le système le plus simple présenté ici consiste à n'utiliser l'approximation de Hartree-Fock (HF) que pour un nombre restreint d'électrons, dits électrons de valence, qui pilotent les propriétés essentielles du système. Un modèle plus sophistiqué, dit configuration d'interactions (CI), consiste lui en une série infinie de déterminants qu'il faut tronquer pertinemment.

**Mots clés** – Approximation de Hartree-Fock, configuration d'interactions, charge effective de Born, méthode des éléments finis.

# 1 Introduction

Les calculs *ab initio* consistent à calculer les propriétés d'un matériau donné (molécule ou cristal) via son état électronique. En effet, ces propriétés (mécaniques, conductives, optiques,...) découlent de l'état des liaisons dans le matériau donc de l'état des électrons du système et particulièrement des électrons de la couche externe, dits électrons de valence. La modélisation *ab initio* passe par la résolution de l'équation de Schrödinger sous une forme variationnelle pour la fonction d'onde électronique du matériau, équation qui n'est malheureusement pas soluble telle quelle. Il existe donc toute une hiérarchie de modèles plus ou moins sophistiqués qui approchent plus ou moins la solution exacte du problème. Il est donc essentiel de pouvoir évaluer l'erreur commise en fonction de la stratégie de résolution adoptée, par rapport à une quantité d'intérêt choisie, la plus naturelle étant l'énergie.

En unités atomiques, l'équation de Schrödinger s'écrit, pour un système de N particules (les noyaux et les électrons du système) :

$$-1/2 \Delta\Psi + V\Psi = E\Psi$$

où $\Psi: \mathbb{R}^{3N} \rightarrow \mathbb{C}$ est la fonction d'onde du système, $V$ le potentiel d'interaction et $E$ l'énergie totale du système.
Dans le cas classique où le seul type d'interaction entre particules est l'interaction coulombienne, on peut séparer $V$ en trois contributions :

$$V = V_{nn} + V_{en} + V_{ee}$$

où $V_{nn}$, $V_{en}$ et $V_{ee}$ désignent respectivement l'interaction noyaux-noyaux, noyaux-électrons et électrons-électrons.

La première approximation couramment utilisée consiste à figer la position des noyaux atomiques : c'est l'approximation adiabatique ou de Born-Oppenheimer. Ainsi, le nombre de particules contenues dans la fonction d'onde à déterminer devient $N_e$, nombre d'électrons du système, et l'énergie $V_{nn}$ est une constante :

$$-\frac{1}{2} \Delta \psi + (V_{en} + V_{ee}) \psi = E_e \psi$$

où $\psi: \mathbb{R}^{3Ne} \to \mathbb{C}$ est la fonction d'onde *électronique* du système, antisymétrique pour respecter le principe d'exclusion de Pauli, et $E_e = E - V_{nn}$ est l'énergie électronique du système étudié.

# 1.1 Approximation de Hartree-Fock et électrons de valence

## 1.1.1 Système de Fock

Mis à part le cas très simple où le système n'est constitué que d'un seul électron, la résolution de cette équation nécessite de rechercher la solution sous une forme particulière donnée en utilisant une formulation variationnelle du problème. Ainsi, la méthode de Hartree-Fock (HF), consiste à rechercher la solution sous une forme qui vérifie par construction le principe d'exclusion de Pauli :

$$\psi(x_j)_{1 \le j \le Ne} = \text{Det}[\phi_i(x_j)]$$

On résout alors un système de $N_e$ équations de Fock, dans notre cas par la méthode des éléments finis. Cette écriture est assez bien adaptée au cas des molécules, mais prend mal en compte la corrélation électronique :

$$-\frac{1}{2} \Delta \mathbf{\Phi} + V_{\text{en}} \mathbf{\Phi} + (\text{Tr} \underline{\mathbf{G}}) \mathbf{\Phi} - \underline{\mathbf{G}} \cdot \mathbf{\Phi} = \begin{pmatrix} \epsilon_1 & 0 & \cdots & 0 \\ 0 & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & \epsilon_{N_e} \end{pmatrix} \mathbf{\Phi}$$

où $\Delta \underline{\mathbf{G}} = -4 \pi \mathbf{\Phi} \otimes \mathbf{\Phi}$ ($g_{ij} = 0$ pour les spins antiparallèles)

et l'énergie électronique totale du système est donnée par :

$$E_e = \frac{1}{2} \sum_i h_i + \epsilon_i \quad \text{où } h_i = \frac{1}{2} \int_{\mathbb{R}^3} |\nabla \varphi_i|^2 + \int_{\mathbb{R}^3} V_{\text{en}} |\varphi_i|^2$$

## 1.1.2 Électrons de valence et charge effective

Par ailleurs, d'autres approximations permettent de simplifier la complexité du calcul, liée essentiellement au nombre d'électrons à traiter. Ainsi, par définition, les électrons de cœur d'un atome ne sont pas perturbés par l'environnement de cet atome, que celui-ci soit isolé, ou lié à d'autres atomes d'espèce quelconques. Il est donc possible, au lieu de les recalculer systématiquement, de les représenter par un pseudopotentiel [1], et de ne calculer que les électrons des couches externes du système, les électrons de valence.
Ainsi, $N = N_v < N_e$.

Cette approximation permet également d'éviter de recalculer tous les termes de **G**. Lors du calcul des orbitales $\phi_i$ contenant les électrons de valence ($N_e-N_v < i \le N_e$), on peut remplacer les termes de **G** faisant intervenir les orbitales des électrons de cœur $\phi_j$ par une orbitale hydrogénoïde paramétrée par $Z_{eff}$, c'est-à-dire solution (analytiquement calculable, et même connue) de :

$$-\frac{1}{2}\Delta\varphi - \frac{Z_{eff}}{r}\varphi = -\frac{Z_{eff}^2}{n^2}\varphi$$

À titre d'exemple, on présente ici le cas de l'hélium à l'état fondamental. Les deux orbitales occupées sont totalement identiques $\phi_1 = \phi_2 = \varphi$, et les états électroniques des deux électrons ne diffèrent que par le spin ($G_{12} = 0$). Cette orbitale est une orbitale 1s (niveau le plus bas, à symétrie sphérique) et la valeur de $Z_{eff}$ est obtenue en minimisant la distance entre les deux fonctions, par exemple au sens de la norme $L^2$. Alors $Z_{eff} = 1,62$, valeur assez éloignée de celle fournie par la littérature ($Z_{eff} = 1,69$) [2].

<figure><img src="image_2.png" /><figcaption>Fig. 1 : Représentation axisymétrique d'un calcul par éléments finis de :
a) $\varphi_1$, orbitale 1s réelle ; b) $|\varphi_1 - \varphi_{1eff}|$ par un calcul de $Z_{eff}$ avec une norme uniforme et c) $|\varphi_1 - \varphi_{1eff}|$ par un calcul de $Z_{eff}$ avec pondération.</figcaption></figure>

On peut comparer les erreurs commises dans la détermination des modes 1s par le calcul de façon exacte et en utilisant l'orbitale approchée $\varphi_{2approx}$ comme exposé ci-dessus (figure 1.b). Par ailleurs, on peut choisir pour déterminer $Z_{eff}$ de minimiser la distance entre l'orbitale réelle et l'orbitale hydrogénoïde la plus proche, mais en pondérant la norme par la probabilité de présence des autres électrons à cet endroit : en effet, pour calculer efficacement et le plus précisément possible les termes $G_{ij}$, il faut que l'orbitale hydrogénoïde $\varphi_{3approx}$ soit une approximation plus particulièrement pertinente dans les régions où $|\varphi_i|$ prend des valeurs importantes. Alors on peut choisir de déterminer $Z_{eff}$ en minimisant :
$$ \int_{\mathbb{R}^3} |\varphi_1|^2 |\varphi_2^2 - \varphi_{2\text{approx}}^2| \quad \text{au lieu de :} \quad \int_{\mathbb{R}^3} |\varphi_2^2 - \varphi_{2\text{approx}}^2| $$

En minimisant cette erreur pondérée, la valeur $Z_{eff}$ obtenue (1,71) est proche de celles de la littérature. La figure 1.c) montre alors que la différence entre les modes 1s calculés exactement et à l'aide de la fonction effective $\varphi_{2approx}$ diminue sensiblement au sens de la norme infinie, mais est légèrement augmentée au sens de la norme $L^2$ (6,6.10⁻⁵ au lieu de 3,4.10⁻⁵). Comme perspective de travail, on chercherait donc à systématiser l'évaluation de l'erreur commise par cette approximation.

## 2 Configuration d'interactions

Parmi les modèles plus sophistiqués, dits post-HF, la recherche d'une solution sous la forme d'une série de tels déterminants (interaction de configurations, CI) :

$$ \psi(x_j)_{1 \le j \le Ne} = \sum_k \alpha_k \det[\varphi_i^k(x_j)] $$

permet d'obtenir une meilleure solution. En effet, l'espace des solutions s'écrivant sous la forme d'une telle série (éventuellement infinie) est dense dans l'espace des fonctions $\mathbb{R}^{3Ne} \to \mathbb{C}$. Le premier déterminant est le déterminant HF initial, et les premiers termes suivants ne diffèrent du premier terme que par une ou deux orbitales dites virtuelles. L'enjeu est alors de déterminer le nombre d'orbitales virtuelles et de termes dans la série à prendre en compte.

En choisissant l'énergie comme quantité d'intérêt, on utilise la stratégie de Rannacher [2]. Le lagrangien :

$$ L(\psi, E, z, E_z) = E \langle \psi | \psi \rangle - E \langle \psi | z \rangle + \frac{1}{2} \langle \nabla \psi | \nabla z \rangle + \langle V \psi | z \rangle + E_z (1 - \langle \psi | \psi \rangle) $$

tient compte des différentes contraintes du problème, à savoir une énergie $E$ minimale respectant l'équation de Schrödinger et une fonction d'onde $\psi$ normalisée. $E_z$ est ainsi le multiplicateur de Lagrange associé à la condition sur la norme du problème, et est homogène à une énergie.

$\langle \cdot | \cdot \rangle$ désigne la notation habituellement utilisée pour le produit hermitien :

$$ \langle \varphi | \psi \rangle = \int \varphi^* \psi $$

Après optimisation du lagrangien, on peut alors représenter l'erreur de la façon suivante :

$$ 2 (E - E_h) = \langle \delta\psi | (H - E_h)z_h \rangle + \langle \psi_h | (H - E_h)\delta z \rangle + R_h, $$

où $H$ est le hamiltonien du système, c'est-à-dire l'opérateur intervenant dans le premier membre de l'équation de Schrödinger, et les grandeurs indicées $_h$ ont été déterminées par le modèle grossier (typiquement HF, ou CI à n termes) et celles non indicées par le modèle plus sophistiqué (typiquement CI à n + 1 termes).

L'erreur commise est donc regroupée dans les produits hermitiens :

$$ \langle \text{Det}[\varphi_i^k(x_j)] | \text{Det}[\varphi_{i'}^{k'}(x_j')] \rangle, $$

où $k \neq k'$. Dans le cas où les déterminants ne diffèrent que d'une ou deux orbitales, ce produit hermitien se limite à des intégrales sur les seules fonctions $\varphi_i$ qui diffèrent entre les deux déterminants considérés.
À titre d'exemple, on traite ici le cas de l'hélium, atome composé d'un noyau de charge Z = 2 et de deux électrons. On se restreint au calcul de l'orbitale 1s réellement occupée par deux électrons de spins antiparallèles à l'état fondamental, et d'une orbitale virtuelle 2s (figure 2). Le développement en série de déterminantes est limité à 4 termes :

- le premier déterminant comprenant les deux orbitales 1s dites occupées ;

- les deuxième et troisième déterminants correspondant aux cas d'une excitation simple;

- le quatrième déterminant correspondant au cas d'une double excitation.

$$
\begin{align*}
\Psi(x_1, x_2) &= \alpha_0 \Psi_0 + \alpha_{1\uparrow} \Psi_{1\uparrow} + \alpha_{1\downarrow} \Psi_{1\downarrow} + \alpha_2 \Psi_2 \\
&= \frac{\alpha_0}{\sqrt{2}} \begin{vmatrix} \varphi_{1s\uparrow}(x_1) & \varphi_{1s\uparrow}(x_2) \\ \varphi_{1s\downarrow}(x_1) & \varphi_{1s\downarrow}(x_2) \end{vmatrix} + \frac{\alpha_{1\uparrow}}{\sqrt{2}} \begin{vmatrix} \varphi_{2s\uparrow}(x_1) & \varphi_{2s\uparrow}(x_2) \\ \varphi_{1s\downarrow}(x_1) & \varphi_{1s\downarrow}(x_2) \end{vmatrix} \\
&\quad + \frac{\alpha_{1\downarrow}}{\sqrt{2}} \begin{vmatrix} \varphi_{1s\uparrow}(x_1) & \varphi_{1s\uparrow}(x_2) \\ \varphi_{2s\downarrow}(x_1) & \varphi_{2s\downarrow}(x_2) \end{vmatrix} + \frac{\alpha_2}{\sqrt{2}} \begin{vmatrix} \varphi_{2s\uparrow}(x_1) & \varphi_{2s\uparrow}(x_2) \\ \varphi_{2s\downarrow}(x_1) & \varphi_{2s\downarrow}(x_2) \end{vmatrix}
\end{align*}
$$

<figure><img src="image_3.png" /><figcaption>Fig. 2 : Représentations axisymétriques d'un calcul par éléments finis des orbitales virtuelles a) 2s et b) 3s, utilisées pour la configuration d'interactions.</figcaption></figure>

La recherche des coefficients ($\alpha_i$) qui minimisent l'énergie de $\Psi$ tout en respectant la condition de normalisation impose de rechercher les vecteurs propres de la matrice hermitienne :

$$
\underline{\mathbf{E}} = \begin{pmatrix} E_{00} & E_{01\uparrow} & E_{01\downarrow} & E_{02} \\ E_{1\uparrow1\uparrow} & E_{1\uparrow1\downarrow} & E_{1\downarrow2} & E_{1\downarrow1\downarrow} \\ E_{1\downarrow1\downarrow} & E_{1\downarrow1\uparrow} & E_{1\downarrow2} & E_{22} \end{pmatrix} = \begin{pmatrix} E_{00} & 0 & 0 & E_{02} \\ E_{11} & -E_{02} & E_{12} & 0 \\ 0 & E_{11} & E_{12} & 0 \\ E_{22} & 0 & 0 & E_{22} \end{pmatrix}
$$

Soit $E = -2.890320095184268 < E_{00}$, valeur de l'énergie obtenue par le modèle HF.

L'erreur commise e avec notre choix de lagrangien correspond directement aux termes
extradiagonaux de la matrice **E** ci-dessus : pour la norme $L^2$,

$$e = \sqrt{2E_{02}^2 + 2E_{12}^2} \approx 0,0114.$$

Le nombre chiffres significatifs mis pour **E** semble déraisonnable, mais la minimisation de **E** par ce biais ne modifie que le $10^e$ chiffre significatif. Par ailleurs, on remarque alors que les deux termes $\alpha_1$ sont très petits devant les autres, et correspondent plutôt à un 0 numérique, situation en accord avec la théorie. Un développement en série sur deux termes uniquement ($\Psi_0$ et $\Psi_2$) est tout aussi efficace en étant moins lourd, puisqu'on ne considère plus alors qu'une matrice 2x2.

De façon plus surprenante il apparaît que le mode virtuel 3s contribue infiniment plus à l'approche de l'énergie fondamentale du système dans le cas présent :

$$\text{Alors} \qquad E = -2.890320143115093.$$

Cela s'explique par le fait que $\epsilon_{2s} < \epsilon_{3s}$ (conformément à nos attentes), mais que $E_e(3s) < E_e(2s)$ car $h_{3s} < h_{2s}$. Ce mode 3s étant plus proche en énergie du mode fondamental au sens HF, les interactions entre ces deux modes sont plus importantes, et ce terme contribue plus à approcher la solution exacte que le mode 2s.

## Références

[1] M.C. Payne, M.P. Teter, D.C. Allan, T.A. Arias, J.D. Joannopoulos. Iterative minimization techniques for ab initio total-energy calculations: molecular dynamics and conjugate gradient, *Review of Modern Physics* **64** (4), 1045-1097, 1992 (article de revue).

[2] E. Clementi and D.L. Raimondi, *J. Chem. Phys.* **38**,2686, 1963 ; E. Clementi, D.L. Raimondi, and W.P. Reinhardt, *J. Chem. Phys.* **47**,1300, 1967.

[3] W. Bangerth, R. Rannacher, *Adaptive finite element methods for differential equations*, éd. Birkhäuser Verlag, 2003.