Sommation d'Ewald
# Sommation d'Ewald

La sommation d'Ewald (ou parfois somme d'Ewald) est une méthode de calcul des énergies d'interaction de systèmes périodiques (et particulier des cristaux), et tout particulièrement les énergies électrostatiques. La sommation d'Ewald est un cas particulier de la formule sommatoire de Poisson, avec le remplacement de la sommation des énergies d'interaction dans l'espace réel par une sommation équivalente dans un espace de Fourier. L'avantage de cette approche est la convergence rapide de la sommation dans l'espace de Fourier comparée à son équivalent dans l'espace réel lorsque les interactions se font à longue distance. Les énergies électrostatiques comprenant à la fois des termes d'interactions de courtes et de longues portées, il est très intéressant de décomposer le potentiel d'interaction en termes de courte portée - dont la sommation se fait dans l'espace réel - et de longue portée - dont la sommation se fait dans l'espace de Fourier. La méthode fut développée par Paul Peter Ewald en 1921 (voir références) afin de déterminer l'énergie électrostatique (et, par là, la constante de Madelung) des cristaux ioniques.

## Principe

La sommation d'Ewald réécrit le potentiel d'interaction comme la somme de deux termes :

    $\phi (\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \phi _{sr}(\mathbf {r} )+\phi _{lr}(\mathbf {r} )$

où $\phi _{sr}(\mathbf {r} )$ représente le terme de courte portée dont la sommation est rapide dans l'espace réel et $\phi _{lr}(\mathbf {r} )$ le terme de longue portée dont la sommation est rapide dans l'espace de Fourier. La partie à longue portée doit être finie pour tous les arguments (en particulier pour $r=0$) mais peuvent avoir une forme mathématique adéquate quelconque, le plus généralement une distribution gaussienne. La méthode postule que la partie à courte portée peut être sommée facilement, ainsi, le problème est de traiter la sommation du terme à longue portée. En raison de l'utilisation d'une intégrale de Fourier, la méthode postule implicitement que le système étudié est infiniment périodique (postulat pertinent pour un cristal parfait - i.e. pour l'intérieur d'un cristal réel). L'unité de base reproduite par périodicité est appelée maille unitaire. Une telle maille est choisie comme « maille centrale » de référence et les cellules reproduites par périodicité sont appelées images.

L'énergie d'interaction à longue portée est la somme des énergies d'interaction entre les charges d'une maille unitaire centrale et toutes les charges du réseau. Ainsi, elle peut être décrite comme une double intégrale deux champs de densités de charges correspondant aux champs de la maille unitaire et de celle du réseau cristallin.

    $E_{lr}=\int \int d\mathbf {r} d\mathbf {r} ^{\prime }\rho _{\text{tot}}(\mathbf {r} )\rho _{uc}(\mathbf {r} ^{\prime })\ \phi _{lr}(\mathbf {r} -\mathbf {r} ^{\prime })$

où le champ de densité de charge de la maille unitaire $\rho _{uc}(\mathbf {r} )$ est une intégration sur les positions $\mathbf {r} _{k}$ des charges $q_{k}$ dans la maille unitaire centrale :

    $\rho _{uc}(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{\mathrm {charges} \ k}q_{k}\delta (\mathbf {r} -\mathbf {r} _{k})$

et le champ de densité de charge totale $\rho _{\text{tot}}(\mathbf {r} )$ est la même intégration sur les charges de la maille unitaire $q_{k}$ et leurs images périodiques :

    $\rho _{\text{tot}}(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{n_{1},n_{2},n_{3}}\sum _{\mathrm {charges} \ k}q_{k}\delta (\mathbf {r} -\mathbf {r} _{k}-n_{1}\mathbf {a} _{1}-n_{2}\mathbf {a} _{2}-n_{3}\mathbf {a} _{3})$

Ici, $\delta (\mathbf {x} )$ est la fonction δ de Dirac, $\mathbf {a} _{1}$, $\mathbf {a} _{2}$ et $\mathbf {a} _{3}$ sont les vecteurs de maille et $n_{1}$, $n_{2}$ et $n_{3}$ sont des entiers. Le champ total $\rho _{\text{tot}}(\mathbf {r} )$ peut être décrit comme une convolution de $\rho _{uc}(\mathbf {r} )$ avec une fonction de réseau $L(\mathbf {r} )$ :

    $L(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{n_{1},n_{2},n_{3}}\delta (\mathbf {r} -n_{1}\mathbf {a} _{1}-n_{2}\mathbf {a} _{2}-n_{3}\mathbf {a} _{3})$

Puisque l'on a une convolution, la transformation de Fourier de $\rho _{\text{tot}}(\mathbf {r} )$ est un produit :

    ${\tilde {\rho }}_{\text{tot}}(\mathbf {k} )={\tilde {L}}(\mathbf {k} ){\tilde {\rho }}_{uc}(\mathbf {k} )$

dans laquelle la transformée de Fourier de la fonction de réseau est une autre intégration sur des fonctions δ :

    ${\tilde {L}}(\mathbf {k} )={\frac {\left(2\pi \right)^{3}}{\Omega }}\sum _{m_{1},m_{2},m_{3}}\delta (\mathbf {k} -m_{1}\mathbf {b} _{1}-m_{2}\mathbf {b} _{2}-m_{3}\mathbf {b} _{3})$

où les vecteurs de l'espace réciproque sont définis comme il suit : $\mathbf {b} _{1}\ {\stackrel {\mathrm {def} }{=}}\ \mathbf {a} _{2}\times \mathbf {a} _{3}/\Omega$ (aux permutations cycliques près) où $\Omega \ {\stackrel {\mathrm {def} }{=}}\ \mathbf {a} _{1}\cdot \left(\mathbf {a} _{2}\times \mathbf {a} _{3}\right)$ est le volume de la maille unitaire centrale (si c'est géométriquement un parallélépipède, ce qui est souvent mais pas nécessairement le cas). $L(\mathbf {r} )$ et ${\tilde {L}}(\mathbf {k} )$ sont des fonctions réelles paires.

Pour la brièveté, on définit un potentiel effectif à une seule particule :

    $v(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \int d\mathbf {r} ^{\prime }\rho _{uc}(\mathbf {r} ^{\prime })\ \phi _{lr}(\mathbf {r} -\mathbf {r} ^{\prime })$

Cette équation étant aussi une convolution, la transformation de Fourier de cette équation est aussi un produit :

    ${\tilde {V}}(\mathbf {k} )\ {\stackrel {\mathrm {def} }{=}}\ {\tilde {\rho }}_{uc}(\mathbf {k} ){\tilde {\Phi }}(\mathbf {k} )$

dans laquelle la transformée de Fourier est définié de la manière suivante :

    ${\tilde {V}}(\mathbf {k} )=\int d\mathbf {r} \ v(\mathbf {r} )\ e^{-i\mathbf {k} \cdot \mathbf {r} }$

L'énergie peut maintenant être écrite comme une intégrale d'un unique champ :

    $E_{lr}=\int d\mathbf {r} \ \rho _{\text{tot}}(\mathbf {r} )\ v(\mathbf {r} )$

En utilisant le théorème de Parseval, l'énergie peut aussi être sommée dans l'espace de Fourier :

    $E_{lk}=\int {\frac {d\mathbf {k} }{\left(2\pi \right)^{3}}}\ {\tilde {\rho }}_{\text{tot}}^{*}(\mathbf {k} ){\tilde {V}}(\mathbf {k} )=\int {\frac {d\mathbf {k} }{\left(2\pi \right)^{3}}}{\tilde {L}}^{*}(\mathbf {k} )\left|{\tilde {\rho }}_{uc}(\mathbf {k} )\right|^{2}{\tilde {\Phi }}(\mathbf {k} )={\frac {1}{\Omega }}\sum _{m_{1},m_{2},m_{3}}\left|{\tilde {\rho }}_{uc}(\mathbf {k} )\right|^{2}{\tilde {\Phi }}(\mathbf {k} )$

où $\mathbf {k} =m_{1}\mathbf {b} _{1}+m_{2}\mathbf {b} _{2}+m_{3}\mathbf {b} _{3}$ dans la dernière intégrale.

C'est le résultat principal. Une fois que ${\tilde {\rho }}_{uc}(\mathbf {k} )$ est calculé, la sommation/intégration sur $\mathbf {k}$ est simple et devrait rapidement converger. La raison la plus courante d'absence de convergence est une mauvaise définition de la maille unitaire, qui doit être électriquement neutre afin d'éviter des sommes infinies.

## Méthode du maillage particulier d'Ewald

La sommation d'Ewald fut développée comme méthode de physique théorique, bien avant la venue des ordinateurs et de l'informatique. Cependant, elle s'est très largement répandue depuis les années 1970 dans les simulations numériques de systèmes de particules, et tout particulièrement celles interagissant par des lois de forces en carré inverse, comme la gravité ou l'électrostatique. Les applications de ces simulations comprennent des objets aussi divers que les plasmas, les galaxies ou les biomolécules.

Comme dans une sommation d'Ewald « normale », un potentiel d'interaction générique est séparé en deux termes : $\phi (\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \phi _{sr}(\mathbf {r} )+\phi _{lr}(\mathbf {r} )$ - un terme d'interaction à courte portée $\phi _{sr}(\mathbf {r} )$ à intégration rapide dans l'espace réel est un terme d'interaction à longue portée $\phi _{lr}(\mathbf {r} )$ qui s'intègre rapidement dans l'espace de Fourier. L'idée de base de la sommation sur un maillage particulier d'Ewald est de remplacer la sommation directe des énergies d'interaction entre les particules ponctuelles :

    $E_{TOT}=\sum _{i,j}\phi (\mathbf {r} _{j}-\mathbf {r} _{i})=E_{sr}+E_{lr}$

avec deux sommations, une intégration directe $E_{sr}$ du potentiel à courte portée dans l'espace réel :

    $E_{sr}=\sum _{i,j}\phi _{sr}(\mathbf {r} _{j}-\mathbf {r} _{i})$

(partie particulière du maillage particulier d'Ewald) et une intégration dans l'espace de Fourier de la partie à longue portée :

    $E_{lr}=\sum _{\mathbf {k} }{\tilde {\Phi }}_{lr}(\mathbf {k} )\left|{\tilde {\rho }}(\mathbf {k} )\right|^{2}$

où ${\tilde {\Phi }}_{lr}$ et ${\tilde {\rho }}(\mathbf {k} )$ représentent les transformées de Fourier du potentiel électrique et la densité de charge (partie Ewald). En raison de leurs convergences rapides dans leurs espaces respectifs (réel et de Fourier), elles peuvent être tronquées sans perte réelle de précision et un gain important du temps de calcul requis. Afin d'évaluer efficacement la transformée de Fourier ${\tilde {\rho }}(\mathbf {k} )$ du champ de densité de charge, on peut utiliser la transformée de Fourier rapide, qui nécessite que le champ de densité soit évaluée sur une maille discrète de l'espace (partie maillage). En raison du postulat de périodicité implicite de la sommation d'Ewald, les applications de la méthode du maillage particulier d'Ewald (particular mesh Ewald en anglais, abrégé par l'acronyme PME) aux systèmes physiques nécessitent d'imposer une symétrie par périodicité. Pour cette raison, la méthode est bien plus efficace pour les systèmes pouvant être simulé comme étant infinis dans l'espace. Dans les simulations de dynamique moléculaire, cette configuration est normalement atteinte en construisant de manière délibérée une maille unitaire de charge neutre qui peut être infiniment reproduite par pavage afin de former des images; cependant, afin de prendre correctement en compte les effets de cette approximation, ces images sont réincorporées dans la boîte de simulation originale. L'effet global est similaire à une version tridimensionnelle du jeu Asteroids, dans laquelle chaque dimension se « replie » sur elle-même. Cette propriété de la maille est appelée condition périodique aux limites. Afin de mieux visualiser cela, on peut imaginer un cube unitaire ; la face du dessus est en contact avec la face du dessous, la face de droite avec celle de gauche, et celle de devant avec celle de derrière. Il en résulte que la taille de la boîte doit être choisie avec précaution pour éviter de générer des corrélations de mouvement impropres entre deux faces « en contact », mais également assez petite pour pouvoir être traitée numériquement. La définition du rayon de coupure, distance à laquelle on passe de la partie d'interaction à courte distance de la partie interaction à longue distance, peut aussi introduire des artéfacts. La restriction du champ de densité à un maillage rend la méthode PME encore plus efficace pour des systèmes avec des variations de densité « douces », ou des fonctions de potentiels continues. Les systèmes localisés (par exemple des systèmes radicalaires) ou avec de fortes variations de densité pourront être traités plus efficacement par la méthode des multipôles rapides de Greengard et Rokhlin.

## Terme dipolaire

L'énergie électrostatique d'un cristal polaire (c'est-à-dire un cristal avec un dipôle identifié $\mathbf {p} _{uc}$ dans la maille unitaire) est semi-convergente, i.e. elle dépend de l'ordre de la sommation. Par exemple, si les interactions dipôle-dipôle d'une maille unitaire centrale avec des mailles unitaires sont localisées sur un cube toujours croissant, l'énergie converge vers une valeur différente de celle obtenue pour une intégration sphérique. En d'autres termes, cette semi-convergence se produit car : premièrement, le nombre de dipôles interagissants dans une couche de rayon R augmente en $R^2$ ; deuxièmement, la force d'une interaction dipôle-dipôle simple décroît en ${\frac {1}{R^{3}}}$ ; et troisièmement, la sommation mathématique $\sum _{n=1}^{\infty }{\frac {1}{n}}$ diverge. Ce résultat relativement surprenant peut être concilié avec le fait que l'énergie des cristaux réels est finie car ceux-ci ne sont pas infinis, et donc des limites particulières. Plus spécifiquement, la frontière d'un cristal polaire montre une densité de charge effective surfacique sur sa surface $\sigma =\mathbf {P} \cdot \mathbf {n}$ où n est le vecteur normal à la surface et P représente le moment dipolaire net volumique. L'énergie d'interaction du dipôle dans une maille unitaire centrale avec cette densité de charge surfacique peut être écrite :

    $U={\frac {1}{2V_{uc}}}\int {\frac {\left(\mathbf {p} _{uc}\cdot \mathbf {r} \right)\left(\mathbf {p} _{uc}\cdot \mathbf {n} \right)dS}{r^{3}}}$

où $\mathbf {p} _{uc}$ et $V_{uc}$ sont les moments dipolaires nets et le volume de la maille unitaire, $dS$ une aire infinitésimale sur la surface du cristal et r est le vecteur de la maille unitaire à cette aire infinitésimale. Cette formule provient de l'intégration de l'énergie $dU=-\mathbf {p} _{uc}\cdot d\mathbf {E}$ où $d\mathbf {E}$ représente le champ électrique infinitésimale généré par une charge surfacique infinitésimale $dq\ {\stackrel {\mathrm {def} }{=}}\ \sigma \,dS$ (loi de Coulomb) :

    $d\mathbf {E} \ {\stackrel {\mathrm {def} }{=}}\ \left({\frac {-1}{4\pi \epsilon }}\right){\frac {dq\ \mathbf {r} }{r^{3}}}=\left({\frac {-1}{4\pi \epsilon }}\right){\frac {\sigma dS\ \mathbf {r} }{r^{3}}}$

Le signe négatif provient de la définition de r, qui pointe vers la charge, et non de la charge.

## Scalabilité de la méthode

De manière générale, différentes méthodes de sommation d'Ewald donnent des scalabilités différentes. Ainsi, les calculs directs donnent une scalabilité en $O(N^2),$ où N est le nombre d'atomes du système. La PME procure elle une scalabilité en O(N logN).

