Une erreur de simplification ? (somme d'Ewald)

Bonjour,

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

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,

… 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 12^iè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,

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,

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.

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

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 :

Pour moi, ça, c'est pas logique

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

Si c'est le cas, avec n = [-1, 0., 1., 2.] ça ne marche plus.

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

… 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}$).

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.

À 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$¹ à la position $x_2$ et que $k \in [-1; 1]$ (je rapelle que $\vec{k}\neq 0$). On trouve assez facilement que,

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

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

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

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

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

… 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.

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' ¹.

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

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.

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