Une erreur de simplification ? (somme d'Ewald)

Sommes et fonctions trigonométriques

L'auteur de ce sujet a trouvé une solution à son problème.
Auteur du sujet

Bonjour,

Pour le moment, je travaille sur l'implémentation d'une somme d'Ewald, À un moment, je dois implémenter la somme suivante,

$$\sum_{\vec{k} \neq 0}\sum_{ij} q_i\,q_j\,\cos(\vec{k}\cdot\vec{r}_{ij})$$

Ou $q_i$ et $q_j$ sont les charges des atomes $i$ et $j$, $r_{ij}$ est la distance entre l'atome $i$ et $j$ ($\vec{r}_{ij} = \vec{r}_j - \vec{r}_i$) et $\vec{k}$ est le vecteur réciproque, tel que $\vec{k} = (n_x \frac{2\pi}{L_x}; n_y \frac{2\pi}{L_y}; n_z \frac{2\pi}{L_z})$ ou $n_i$ est un nombre entier compris entre -1 et 1 (on pourrait choisir des nombre plus grand, mais c'est plus facile pour débugger). On exclu $n_x = n_y = n_z = 0$ de la somme (parce que dans l'expression complète de la somme d'Ewald, il y a un $|\vec{k}|^{-1}$).

Très justement, la publication que je suis me propose la simplification suivante,

$$\sum_{\vec{k} \neq 0}\sum_{ij} q_i\,q_j\,\cos(\vec{k}\cdot\vec{r}_{ij}) = \sum_{\vec{k}\neq 0}\left[\sum_i q_1\,cos{(\vec{k}\cdot\vec{r}_i)}\right]^2 + \left[\sum_i q_1\,sin{(\vec{k}\cdot\vec{r}_i)}\right]^2$$

… Ce qui est tout à fait correct si on prend le temps de l'écrire sur une feuille de papier et en se rappellant que $\cos{(k\cdot\vec{r}_{ij})} = \cos{(k\cdot\vec{r}_{j}-k\cdot\vec{r}_{i})}$ et que $\cos{(a-b)} = \cos{a}\cos{b} + \sin{a}\sin{b}$.

Implémentée (très) grossièrement en python, on a les codes suivant…

  • Pour $\sum_{\vec{k} \neq 0}\sum_{ij} q_i\,q_j\,\cos(\vec{k}\cdot\vec{r}_{ij})$,
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import math

atoms = ([1, 2, 1], [-1, 0, -1], [1, 1, 1], [0, 0, 0])
charges = (.7, -.3, -.2, -.2)  # il faut que ça soit électriquement neutre, que ∑ q_i = 0
reciprocal = [1/3, 1/3, 1/3]

n = [-1, 0, 1]

U1 = 0

for nx in n:
    for ny in n:
        for nz in n:
            if nx == ny == nz == 0:
                continue

            n_ = [nx, ny, nz]
            k = [n_[a] * reciprocal[a] for a in range(3)]

            for i in range(len(atoms)):
                for j in range(len(atoms)):
                    k_dot_r = k[0] * (atoms[j][0]-atoms[i][0]) \
                              + k[1] * (atoms[j][1]-atoms[i][1]) \
                              + k[2] *(atoms[j][2]-atoms[i][2])
                    U1 += charges[i] * charges[j] * math.cos(k_dot_r)

print(U1)
  • Et pour $\sum_{\vec{k}\neq 0}\left[\sum_i q_1\,cos{(\vec{k}\cdot\vec{r}_i)}\right]^2 + \left[\sum_i q_1\,sin{(\vec{k}\cdot\vec{r}_i)}\right]^2$, on aura
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# ... suite du précédent
U2 = 0

for nx in n:
    for ny in n:
        for nz in n:
            if nx == ny == nz == 0:
                continue

            n_ = [nx, ny, nz]
            k = [n_[a] * reciprocal[a] for a in range(3)]

            cos_part = 0.0
            sin_part = 0.0

            for i in range(len(atoms)):
                k_dot_r = k[0] * atoms[i][0] + k[1] * atoms[i][1] + k[2] * atoms[i][2]

                cos_part += charges[i] * math.cos(k_dot_r)
                sin_part += charges[i] * math.sin(k_dot_r)

            U2 += cos_part**2 + sin_part**2
print(U2)

… Et on obtient bien 2 fois la même réponse (logique), à l'erreur numérique près (les résultats divergent après la 12ième décimale ou un truc du genre).

La question que je me pose, c'est pourquoi j'obtient la même réponse avec le code suivant,

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
# ... C'est toujours la suite du précédent
U3 = 0

for nx in n:
    for ny in n:
        for nz in n:
            if nx == ny == nz == 0:
                continue

            n_ = [nx, ny, nz]
            k = [n_[a] * reciprocal[a] for a in range(3)]

            sin_p_part = 0.0

            for i in range(len(atoms)):
                k_dot_r = k[0] * atoms[i][0] + k[1] * atoms[i][1] + k[2] * atoms[i][2]
                sin_p_part += charges[i] * math.sin(k_dot_r + math.pi/4)

            U3 += 2 * sin_p_part **2

print(U3)

… Qui correspondrait à la simplification suivante,

$$\sum_{\vec{k}\neq 0}\left[\sum_i q_1\,cos{(\vec{k}\cdot\vec{r}_i)}\right]^2 + \left[\sum_i q_1\,sin{(\vec{k}\cdot\vec{r}_i)}\right]^2 = \sum_{\vec{k}\neq 0} 2 \times\left[\sum_i q_1\,sin{\left(\vec{k}\cdot\vec{r}_i + \frac{\pi}{4}\right)}\right]^2$$

J'ai "trouvé" cette simplification sur base d'une erreur de raisonnement, pensant que $\sin(x)+\cos(x) = \sqrt{2}\,\sin{\left(x+\frac{\pi}{4}\right)}$, ce qui est correct en soit, mais inapplicable à mon cas, puisque je ne peux pas dire que,

$$\sum_{\vec{k}\neq 0}\left[\sum_i q_1\,cos{(\vec{k}\cdot\vec{r}_i)}\right]^2 + \left[\sum_i q_1\,sin{(\vec{k}\cdot\vec{r}_i)}\right]^2 \overset{?}{=} \sum_{\vec{k}\neq 0}\left[\sum_i q_1\,cos{(\vec{k}\cdot\vec{r}_i)} + q_1\,sin{(\vec{k}\cdot\vec{r}_i)}\right]^2$$

Ce qui est pour moi, à première vue, une grosse erreur. Et pourtant, ça fonctionne ! Et ça m'intéresse, puisque en terme de performances ça ne me fait plus qu'un seul cosinus par atome à évaluer, contre un sinus et un cosinus dans le second cas et une pelletée dans le premier (dans cas réel, le nombre de fois qu'on fait cette évaluation est important, le nombre de $\vec{k}$ est important et le nombre d'atomes est important, donc ça commence à compter).

Alors, si quelqu'un pouvait m'expliquer pourquoi ça fonctionne ou l'endroit ou je me suis trompé en implémentant mon code, je vous avoue que ça me rassurerait. Merci d'avance :)

PS: bien entendu, j'ai testé en changeant le nombre de charges, les positions, le nombre de vecteurs réciproques à considérer, etc, et le code me sort toujours trois réponses identiques.

Édité par pierre_24

Doctorant et assistant en chimie à l'Université de NamurEx-dev' pour ZdS (a aidé à réaliser la ZEP-12 !) • Carniste cis (y parait que c'est une injure)

+0 -0

Quand on a le nez dessus, c'est souvent plus difficile mais:

$$ 2 * [sin(x + \frac{\pi}{4})]^{2} = 2 * [(sin(x) + cos(x)) * \frac{\sqrt{2}}{2}]^{2} = [cos(x) + sin(x)]^{2} $$
Je pense qu'utiliser l'identité:
$$ 1 + sin(2*x) $$
serait encore meilleure.

Édité par Gawaboumga

+0 -0
Auteur du sujet

Tout à fait d'accord (et effectivement, je n'avais pas pensé à l'identité $1+\sin{2x}$).

… Mais c'est à l'étape avant que je pige pas :

$$\sum_{\vec{k}\neq 0}\left[\sum_i q_1\,cos{(\vec{k}\cdot\vec{r}_i)}\right]^2 + \left[\sum_i q_1\,sin{(\vec{k}\cdot\vec{r}_i)}\right]^2 \overset{?}{=} \sum_{\vec{k}\neq 0}\left[\sum_i q_1\,cos{(\vec{k}\cdot\vec{r}_i)} + q_1\,sin{(\vec{k}\cdot\vec{r}_i)}\right]^2$$

Pour moi, ça, c'est pas logique ^^

Doctorant et assistant en chimie à l'Université de NamurEx-dev' pour ZdS (a aidé à réaliser la ZEP-12 !) • Carniste cis (y parait que c'est une injure)

+0 -0

Cette réponse a aidé l'auteur du sujet

Tu dois pouvoir faire ton calcul avec n'importe quel ensemble de $k$ ?

Normalement c'est quelque chose que l'on fait pour n'importe quel ensemble de $k$ de la forme $[-k, k]$.

Sinon je n'ai pas une réponse à ta question, mais j'ai utilisé une autre astuce pour Ewald. Si tu passe le cosinus en complexe, tu as une définition de ta somme plus haut qui utilise une définition récursive de la phase du cosinus, et qui ne fait donc intervenir qu'une seule évaluation du cos. Mais quand j'essaye de transformer ta version de l'équation pour la faire correspondre à cette définition récursive, je me plante au milieu. Je vais ré-essayer quand j'aurai un moment devants moi, avec un papier et un crayon.

Mon Github — Tuto Homebrew — Article Julia

+0 -0
Auteur du sujet

Luthaf a tout à fait raison (j'aurai du le préciser), l'ensemble doit être "symétrique". Je me doute d'ailleurs bien que la raison de cette équivalence doit provenir de là (des histoires du genre $\cos{-x} = \cos{x}$ et $\sin{-x} = -\sin{x}$), mais là, ça ne me saute pas aux yeux (et le fait que ça soit en 3D n'aide pas vraiment, puisque même avec $n \in [-1; 0; 1]$, ça me fait déjà 26 termes, donc c'est pas marrant à écrire sur une feuille de papier).

Sinon, en notation exponentielle, j'ai

$$\sum_{\vec{k} \neq 0}\sum_{ij} q_i\,q_j\,\cos(\vec{k}\cdot\vec{r}_{ij}) = \sum_{\vec{k} \neq 0} \left[q_i\,e^{i\vec{k}\cdot\vec{r}_i}\right]\,\left[q_i\,e^{-i\vec{k}\cdot\vec{r}_i}\right]$$

… Mais je préfère éviter, y'a pas d'exponentielle complexe native en C (même si ça revient à $\cos{x}+i\sin{x}$).

Édité par pierre_24

Doctorant et assistant en chimie à l'Université de NamurEx-dev' pour ZdS (a aidé à réaliser la ZEP-12 !) • Carniste cis (y parait que c'est une injure)

+0 -0

J'arrive jusque :

$\sum_k\left[\left\{\sum_i q_i\sin\left(k.r_i+\frac{\pi}{4}\right)\right\}^2+\left\{\sum_i q_i\sin\left(k.r_i-\frac{\pi}{4}\right)\right\}^2\right]$

En forçant l'introduction des $\frac{\pi}{4}$ depuis le début.

J'aurais envie de conclure en invoquant un réordonnement des $k$ qui permettrait d'avoir ce qu'il faut (grâce aux symétries de ton espace des $k$), mais ca me semble pas évident pour le moment.

Édité par Freedom

Auteur du sujet

À noter que, même si je m'en doutais un peu, ça fonctionne en 1D.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
# version 1D
U1 = 0
for nx in n:
    if nx == 0:
        continue

    k = nx * reciprocal[0]

    for i in range(len(atoms)):
        for j in range(len(atoms)):
            k_dot_r = k * (atoms[j][0]-atoms[i][0])
            U1 += charges[i] * charges[j] * math.cos(k_dot_r)

print(U1)

U2 = 0

for nx in n:
    if nx == 0:
        continue

    k = nx * reciprocal[0]

    cos_part = 0.0
    sin_part = 0.0

    for i in range(len(atoms)):
        k_dot_r = k * atoms[i][0]

        cos_part += charges[i] * math.cos(k_dot_r)
        sin_part += charges[i] * math.sin(k_dot_r)

    U2 += cos_part**2 + sin_part**2

print(U2)

U3 = 0

for nx in n:
    if nx == 0:
        continue

    k = nx * reciprocal[0]

    sin_p_part = 0.0

    for i in range(len(atoms)):
        k_dot_r = k * atoms[i][0]
        sin_p_part += charges[i] * math.sin(k_dot_r + math.pi/4)

    U3 += 2 * sin_p_part **2

print(U3)

Du coup, je crois que j'ai trouvé, parce que ça fait tout de suite moins de terme à traiter, c'est déjà plus facile à écrire sur le papier.

Mettons qu'on est en une dimension, qu'on a un atome 1 avec une charge $q_1$ à la position $x_1$ et un atome 2 avec une charge $q_2$1 à la position $x_2$ et que $k \in [-1; 1]$ (je rapelle que $\vec{k}\neq 0$). On trouve assez facilement que,

$$\sum_{-1, 1}^{k}\sum_{ij} q_i\,q_j\,\cos(k\,r_{ij}) = q_1^2 + q_2^2 + 2\,q_1\,q_2\,[\cos{(x_1-x_2)+\cos{(x_2-x_1)}}]$$

Bon, maintenant, si on prend l'expression avec $\sin{\left(\vec{k}\cdot\vec{r}_i + \frac{\pi}{4}\right)}$ et qu'on joue avec

$$\begin{align} &\sum_{-1,1}^{k} 2 \times\left[\sum_i q_i\,\sin{\left(k\,x_i + \frac{\pi}{4}\right)}\right]^2 \\ % &= \sum_{-1,1}^{k} 2 \times\left[q_1\,\sin{\left(k\,x_1 + \frac{\pi}{4}\right)} + q_2\,\sin{\left(k\,x_2 + \frac{\pi}{4}\right)}\right]^2 \\ % &= \sum_{-1,1}^{k} 2\,q_1^2\,\sin^2{\left(k\,x_1 + \frac{\pi}{4}\right)} + 2\,q_2^2\,\sin^2{\left(k\,x_2 + \frac{\pi}{4}\right)} + 4\,q_1\,q_2\,\sin{\left(k\,x_1 + \frac{\pi}{4}\right)}\,\sin{\left(k\,x_2 + \frac{\pi}{4}\right)} \end{align}$$

Or, $\sin^2{\left(x + \frac{\pi}{4}\right)} = \frac{1}{2}\,[\sin{(2x)}+1]$, donc on obtient,

$$q_1^2+ q_2^2+\sum_{-1,1}^{k} q_1^2\,\sin{(2k\,x_1)} +q_2^2\,\sin{(2k\,x_2)} + 4\,q_1\,q_2\,\sin{\left(k\,x_1 + \frac{\pi}{4}\right)}\,\sin{\left(k\,x_2 + \frac{\pi}{4}\right)}$$

On voit déjà presque ou je veux en venir, mais je vais quand même pousser une étape ou deux plus loin. On sait que $\sin{\left(x + \frac{\pi}{4}\right)} = \frac{1}{\sqrt{2}}\,[\cos{x} + \sin{x}]$. Du coup, on trouve que

$$\begin{align} &q_1^2+ q_2^2+\sum_{-1,1}^{k} q_1^2\,\sin{(2k\,x_1)} +q_2^2\,\sin{(2k\,x_2)} + 4\,q_1\,q_2\,\sin{\left(k\,x_1 + \frac{\pi}{4}\right)}\,\sin{\left(k\,x_2 + \frac{\pi}{4}\right)} \\ &= q_1^2+ q_2^2+\sum_{-1,1}^{k} q_1^2\,\sin{(2k\,x_1)} +q_2^2\,\sin{(2k\,x_2)} + 2\,q_1\,q_2\,\left\{[\sin{(k\,x_1)}+\cos{(k\,x_1)}]\,[\sin{(k\,x_2)}+\cos{(k\,x_2)}]\right\} \\ &= q_1^2+ q_2^2+\sum_{-1,1}^{k} q_1^2\,\sin{(2k\,x_1)} +q_2^2\,\sin{(2k\,x_2)} + 2\,q_1\,q_2\,\cos{(k\,x_1-k\,x_2)} + 2\,q_1\,q_2\,\sin{(k\,x_1+k\,x_2)} \end{align}$$

Et là ou ça fonctionne, c'est évidement parce que $\sin{(-x)} = -\sin{x}$ et que l'espace des $k$ est symétrique, donc les termes $q_1^2\,\sin{(2k\,x_1)}$, $q_2^2\,\sin{(2k\,x_2)}$ et $2\,q_1\,q_2\,\sin{(k\,x_1+k\,x_2)}$ se simplifient ! Seul le terme en $\cos$ reste, et on retombe bien sur le résultat du début.

Bon, c'est pas très très rigoureux et surtout, c'est une "démonstration" basée sur un exemple, mais j'ai le sentiment, comme Freedom, que la clé, c'est $\sin{(-x)} = -\sin{x}$, et que comme on a toujours $\vec{k}$ et $-\vec{k}$ dans la somme, c'est ce qui fait marcher le machin. À voir comment utiliser ça proprement.

EDIT: se prouve à peu près de la même manière en écrivant que

$$\left[\sum_i^N a_i\right]^2 = \sum_i^N a_i^2 + \sum_{j>i}^N 2\,a_i\,a_j$$

Cette formule porte probablement un nom, mais comme c'est moi qui vient d'y penser, je le connais pas.


  1. qui, pour le coup, est égal à $-q_1$, mais c'est pas la question, et ça marche de toute façon. 

Édité par pierre_24

Doctorant et assistant en chimie à l'Université de NamurEx-dev' pour ZdS (a aidé à réaliser la ZEP-12 !) • Carniste cis (y parait que c'est une injure)

+1 -0

Cette réponse a aidé l'auteur du sujet

Voila le détail d'où j'étais arrivé :

$$\begin{aligned}\sum_{k,i,j}q_i q_j \cos(k.r_{ij}) &= \sum_{k,i,j}\frac{q_i q_j}{2} \left(e^{\imath k.r_{ij}}+e^{-\imath k.r_{ij}}\right) \\ &= \sum_{k,i,j}\frac{q_i q_j}{2} \left(e^{\imath k.r_j}e^{-\imath k.r_i}+e^{-\imath k.r_j}e^{\imath k.r_i}\right) \\ &= \sum_k \left[\sum_i q_i e^{\imath k.r_i} \right] \left[\sum_i q_i e^{-\imath k.r_i} \right] \\ &= \sum_k \left|\sum_i q_i e^{\imath k.r_i} \right|^2 \\ &= \sum_k \left|\sum_i q_i e^{\imath \left(k.r_i+\alpha\right)} \right|^2 \forall\alpha\in\mathbb{R} \\ \sum_{k,i,j}q_i q_j \cos(k.r_{ij}) &= \sum_k \left\{ \left[\sum_i q_i \cos\left(k.r_i+\alpha\right)\right]^2 + \left[\sum_i q_i \sin\left(k.r_i+\alpha\right)\right]^2\right\} \forall\alpha\in\mathbb{R} \end{aligned}$$

Avec $\alpha=\frac{\pi}{4}$ tu fais apparaître la moitié de ta somme, il reste à montrer que :

$$ \sum_k \left[\sum_i q_i \cos\left(k.r_i+\frac{\pi}{4}\right)\right]^2 = \sum_k \left[\sum_i q_i \sin\left(k.r_i+\frac{\pi}{4}\right)\right]^2 $$

Ce qui est vrai grâce à tes symétries :

$$ \begin{aligned} \sum_k \left[\sum_i q_i \cos\left(k.r_i+\frac{\pi}{4}\right)\right]^2 &= \sum_k \left[\sum_i q_i \cos\left(-k.r_i+\frac{\pi}{4}\right)\right]^2 \\ &= \sum_k \left[\sum_i q_i \cos\left(k.r_i-\frac{\pi}{4}\right)\right]^2 \\ \sum_k \left[\sum_i q_i \cos\left(k.r_i+\frac{\pi}{4}\right)\right]^2 &= \sum_k \left[\sum_i q_i \sin\left(k.r_i+\frac{\pi}{4}\right)\right]^2 \end{aligned}$$

Édité par Freedom

… Mais je préfère éviter, y'a pas d'exponentielle complexe native en C (même si ça revient à cosx+isinx

Normalement avec cette méthode, tu n'as besoin que d'une seule évaluation de l'exponentielle, que tu peut justement faire à la main. Et si tu peut faire du C99, tu as le droit aux exponentielles complexes ! MSVC ne supporte pas cette partie du C99 ça par contre.

Mon Github — Tuto Homebrew — Article Julia

+0 -0
Auteur du sujet

Ah, j'avais pas vu que Freedom avait posté la réponse entre temps, merci ;)

Normalement avec cette méthode, tu n'as besoin que d'une seule évaluation de l'exponentielle, que tu peut justement faire à la main. Et si tu peut faire du C99, tu as le droit aux exponentielles complexes !

Ah, je savais pas (on en apprend tout les jours). En fait, je partais du principe que la méthode ici ne me donnait qu'un seul cosinus à évaluer, contre un sinus et un cosinus si je l'évaluais "à la main" ($\cos{x}+i\sin{x}$). Après, si il y a une fonction native, j'imagine que les mecs qui l'ont codé derrière ont fait quelque chose de plus intelligent que ça et que ça doit ce valoir question perf' 1.

MSVC ne supporte pas cette partie du C99 ça par contre.

Luthaf

… Loin de moi l'idée d'utiliser ce truc :-°


  1. après, on est plus dans les années 80 ou la différence entre l'évaluation d'une ou deux fonction trigonométrique était dramatique, mais quand même, si je peux gagner du temps, ça m'intéresse. Les temps de calculs sont limités sur le supercalculateur ;) 

Édité par pierre_24

Doctorant et assistant en chimie à l'Université de NamurEx-dev' pour ZdS (a aidé à réaliser la ZEP-12 !) • Carniste cis (y parait que c'est une injure)

+0 -0

Ah, je savais pas (on en apprend tout les jours). En fait, je partais du principe que la méthode ici ne me donnait qu'un seul cosinus à évaluer, contre un sinus et un cosinus si je l'évaluais "à la main" (cosx+isinx)

Normalement, un appel à cos et un appel à sin avec le même argument proche l'un de l'autre sont optimisé en un unique appel à sincos qui calcule les deux d'un coup. Et il y a des versions optimisées si ton processeur supporte des instructions SIMD. Donc on a bien une seule fonction trigonométrique à évaluer.

Mon Github — Tuto Homebrew — Article Julia

+0 -0
Auteur du sujet

On a déjà vu que mon processeur était pas le plus récent :-° (et pourtant, le PC à peut être genre 3 ans).

Par contre, qu'on soit bien d'accord, sincos(), c'est pas ANSI, c'est un truc ajouté, non ? (je l'ai vu dans GNU C, par exemple). Après, c'est ma faute, j'ai choisi de coder en ANSI C un peu par défaut, je peux toujours switcher en GNU C. Ou aller emprunter le code à quelqu'un (si la licence le permet).

… J'ai aussi vu un truc bizarre chez Google à ce sujet ^^

Doctorant et assistant en chimie à l'Université de NamurEx-dev' pour ZdS (a aidé à réaliser la ZEP-12 !) • Carniste cis (y parait que c'est une injure)

+0 -0
Vous devez être connecté pour pouvoir poster un message.
Connexion

Pas encore inscrit ?

Créez un compte en une minute pour profiter pleinement de toutes les fonctionnalités de Zeste de Savoir. Ici, tout est gratuit et sans publicité.
Créer un compte