Sommation d'Ewald

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

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

où ${\displaystyle \phi _{sr}(\mathbf {r} )}$ représente le terme de courte portée dont la sommation est rapide dans l'espace réel et ${\displaystyle \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 ${\displaystyle 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.

${\displaystyle 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 ${\displaystyle \rho _{uc}(\mathbf {r} )}$ est une intégration sur les positions ${\displaystyle \mathbf {r} _{k}}$ des charges ${\displaystyle q_{k}}$ dans la maille unitaire centrale :

${\displaystyle \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 ${\displaystyle \rho _{\text{tot}}(\mathbf {r} )}$ est la même intégration sur les charges de la maille unitaire ${\displaystyle q_{k}}$ et leurs images périodiques :

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

${\displaystyle 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 ${\displaystyle \rho _{\text{tot}}(\mathbf {r} )}$ est un produit :

${\displaystyle {\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 δ :

${\displaystyle {\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 : ${\displaystyle \mathbf {b} _{1}\ {\stackrel {\mathrm {def} }{=}}\ \mathbf {a} _{2}\times \mathbf {a} _{3}/\Omega }$ (aux permutations cycliques près) où ${\displaystyle \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). ${\displaystyle L(\mathbf {r} )}$ et ${\displaystyle {\tilde {L}}(\mathbf {k} )}$ sont des fonctions réelles paires.

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

${\displaystyle 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 :

${\displaystyle {\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 :

${\displaystyle {\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 :

${\displaystyle 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 :

${\displaystyle 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ù ${\displaystyle \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 ${\displaystyle {\tilde {\rho }}_{uc}(\mathbf {k} )}$ est calculé, la sommation/intégration sur ${\displaystyle \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.

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 : ${\displaystyle \phi (\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \phi _{sr}(\mathbf {r} )+\phi _{lr}(\mathbf {r} )}$ - un terme d'interaction à courte portée ${\displaystyle \phi _{sr}(\mathbf {r} )}$ à intégration rapide dans l'espace réel est un terme d'interaction à longue portée ${\displaystyle \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 :

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

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

${\displaystyle 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 :

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

où ${\displaystyle {\tilde {\Phi }}_{lr}}$ et ${\displaystyle {\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 ${\displaystyle {\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.

L'énergie électrostatique d'un cristal polaire (c'est-à-dire un cristal avec un dipôle identifié ${\displaystyle \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² ; deuxièmement, la force d'une interaction dipôle-dipôle simple décroît en ${\displaystyle {\frac {1}{R^{3}}}}$ ; et troisièmement, la sommation mathématique ${\displaystyle \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 ${\displaystyle \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 :

${\displaystyle 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ù ${\displaystyle \mathbf {p} _{uc}}$ et ${\displaystyle V_{uc}}$ sont les moments dipolaires nets et le volume de la maille unitaire, ${\displaystyle 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 ${\displaystyle dU=-\mathbf {p} _{uc}\cdot d\mathbf {E} }$ où ${\displaystyle d\mathbf {E} }$ représente le champ électrique infinitésimale généré par une charge surfacique infinitésimale ${\displaystyle dq\ {\stackrel {\mathrm {def} }{=}}\ \sigma \,dS}$ (loi de Coulomb) :

${\displaystyle 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.

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²), où N est le nombre d'atomes du système. La PME procure elle une scalabilité en O(N logN)[1].

• (en) Cet article est partiellement ou en totalité issu de l’article de Wikipédia en anglais intitulé « Ewald summation » (voir la liste des auteurs).