Grille FFT, marges spatiales et vérification des intégrales \(J\) et \(K\)

Cette page rassemble les précautions à prendre dans les programmes molecule_h2_01_03_integrales_2e.py et molecule_h2_01_04_verification_K_fft.py concernant :

le choix de N_GRILLE pour la FFT ;
le centrage de la grille ;
les paramètres MARGE_BOITE et MARGE_NOYAU ;
les petites et grandes distances internucléaires \(R\) ;
les limites asymptotiques des intégrales coulombienne \(J\) et d’échange \(K\).

1. Taille de grille et efficacité de la FFT

Il n’est pas nécessaire que N_GRILLE soit une puissance de 2. Les FFT sont généralement très efficaces lorsque leur dimension se décompose en petits facteurs premiers :

\[ N = 2^a\,3^b\,5^c\,7^d. \]

Dimension	Décomposition	Appréciation
64	\(2^6\)	Excellent
72	\(2^3\times 3^2\)	Excellent
80	\(2^4\times 5\)	Excellent
81	\(3^4\)	Très bon
96	\(2^5\times 3\)	Excellent
121	\(11^2\)	Moins favorable

Un grand nombre premier est généralement moins favorable. Cependant, il faut regarder la forme exacte du tableau transmis à la FFT.

N_GRILLE = 65

Si la grille contient 65 points bornes comprises, elle possède 64 intervalles :

\[ 65-1=64=2^6. \]

C’est commode pour la géométrie de la grille, mais la FFT n’utilise une taille 64 que si le tableau réellement transmis à fftn contient 64 points par direction.

Si le programme applique directement :

np.fft.fftn(tableau)

à un tableau de forme :

(65, 65, 65)

alors la FFT travaille bien avec :

\[ 65=5\times 13. \]

Cette taille reste acceptable, même si elle est moins favorable que 64, 72 ou 80.

Pourquoi utiliser des grilles impaires ?

La suite :

49, 65, 81, 121

correspond à :

\[ 49=2\times24+1,\qquad 65=2\times32+1,\qquad 81=2\times40+1,\qquad 121=2\times60+1. \]

Une grille impaire possède un point central exactement situé à l’origine. Cela facilite le placement symétrique du milieu de la molécule et du plan médiateur.

Conclusion pratique. Pour une étude de convergence physique, conserver des tailles impaires comme 49, 65, 81 ou 121 est raisonnable. Le gain de vitesse obtenu en passant de 65 à 64 peut être moins important que la perte du point central exact.

2. Définition de la boîte spatiale

Le programme utilise une demi-largeur de boîte de la forme :

L = max(MARGE_BOITE, R / 2 + MARGE_NOYAU)

La boîte s’étend alors approximativement de \(-L\) à \(+L\) dans chaque direction. Le pas de grille vaut :

\[ h=\frac{2L}{N_{\mathrm{GRILLE}}-1}. \]

Avec :

N_GRILLE = 65
MARGE_BOITE = 8.0
MARGE_NOYAU = 6.5

pour :

\[ \frac{R}{2}+6.5<8, \]

c’est-à-dire :

\[ R<3\,a_0, \]

on a :

\[ L=8\,a_0, \qquad h=\frac{16}{64}=0.25\,a_0. \]

Ce pas est raisonnable pour une première étude, mais il reste assez grossier pour décrire précisément le voisinage immédiat des noyaux.

3. Précautions lorsque \(R\) est petit

3.1. Taille de la boîte

Lorsque \(R\to0\), les deux noyaux sont proches du centre. Une boîte allant de \(-8a_0\) à \(+8a_0\) est largement suffisante pour des orbitales hydrogénoïdes \(1s\).

À une distance de \(8a_0\) :

\[ \psi(r)\propto e^{-r}, \qquad e^{-8}\simeq3.35\times10^{-4}. \]

La densité est alors de l’ordre de :

\[ |\psi|^2\propto e^{-16}\simeq1.1\times10^{-7}. \]

Les effets de bord sont donc faibles.

3.2. Résolution de la séparation des noyaux

Le problème principal n’est pas la taille de la boîte, mais le nombre de points entre les deux noyaux.

Avec \(h=0.25a_0\) :

\(R\)	Nombre d’intervalles entre les noyaux	Appréciation
\(1.0a_0\)	4	Acceptable
\(0.5a_0\)	2	Assez grossier
\(0.25a_0\)	1	Insuffisant
\(<0.25a_0\)	moins de 1	Non résolu

Une règle pratique consiste à demander :

\[ \frac{R}{h}\gtrsim4\text{ à }6. \]

Avec \(h=0.25a_0\), cela donne :

\[ R\gtrsim1\text{ à }1.5\,a_0. \]

Pour étudier des distances plus petites, on peut augmenter la résolution :

N_GRILLE = 129

donne, pour \(L=8a_0\) :

\[ h=\frac{16}{128}=0.125\,a_0. \]

Et :

N_GRILLE = 257

donne :

\[ h=\frac{16}{256}=0.0625\,a_0. \]

3.3. Singularités coulombiennes sur un point de grille

Si le programme calcule explicitement des termes comme :

1.0 / r_A
1.0 / r_B

il faut éviter qu’un point de grille conduise directement à \(r_A=0\) ou \(r_B=0\), car cela produit une divergence numérique.

Une protection élémentaire est :

epsilon = 1.0e-12
r_A = np.maximum(r_A, epsilon)
r_B = np.maximum(r_B, epsilon)

La singularité coulombienne est intégrable dans une intégrale tridimensionnelle, mais sa représentation ponctuelle peut produire inf ou nan. Pour une grande précision, une moyenne analytique du potentiel sur la maille centrale est préférable à un simple découpage par epsilon.

4. Précautions lorsque \(R\) est grand

4.1. Marge autour de chaque noyau

Lorsque \(R>3a_0\), on a :

\[ L=\frac{R}{2}+6.5. \]

Le noyau droit est placé en \(+R/2\), tandis que la frontière droite est en \(+L\). La distance entre le noyau et la frontière est donc :

\[ L-\frac{R}{2}=6.5\,a_0. \]

Chaque noyau conserve donc une marge extérieure constante égale à MARGE_NOYAU.

À \(6.5a_0\) :

\[ e^{-6.5}\simeq1.50\times10^{-3}, \qquad e^{-13}\simeq2.26\times10^{-6}. \]

Pour une vérification grossière, MARGE_NOYAU = 6.5 peut convenir. Pour vérifier des intégrales très petites, il est préférable de comparer avec :

MARGE_NOYAU = 8.0

puis éventuellement :

MARGE_NOYAU = 10.0

4.2. Le pas augmente si \(N_{\mathrm{GRILLE}}\) reste fixe

Même si la boîte reste assez large, la résolution se dégrade lorsque \(R\) augmente.

\(R\)	\(L=R/2+6.5\)	\(h=2L/64\)
\(10a_0\)	\(11.5a_0\)	\(0.359a_0\)
\(20a_0\)	\(16.5a_0\)	\(0.516a_0\)
\(40a_0\)	\(26.5a_0\)	\(0.828a_0\)

À \(R=40a_0\), une orbitale \(1s\), dont l’échelle caractéristique est \(a_0\), n’est plus décrite que par quelques points.

Il faut donc augmenter N_GRILLE avec \(R\). Pour conserver un pas cible \(h_{\mathrm{cible}}\) :

\[ N_{\mathrm{GRILLE}} \simeq 1+\frac{2L}{h_{\mathrm{cible}}}. \]

Avec \(h_{\mathrm{cible}}=0.25a_0\) :

\(R\)	\(L\)	\(N_{\mathrm{GRILLE}}\) minimal approximatif
\(10a_0\)	11.5	93
\(20a_0\)	16.5	133
\(30a_0\)	21.5	173
\(40a_0\)	26.5	213

5. Limite à grande distance de l’intégrale coulombienne \(J\)

On considère :

\[ J=(AA|BB) = \iint \frac{ |\phi_A(\mathbf r_1)|^2 |\phi_B(\mathbf r_2)|^2 }{ r_{12} } \,d^3r_1\,d^3r_2. \]

À grande distance, chaque densité normalisée se comporte comme une charge ponctuelle vue depuis l’autre centre. On obtient :

\[ J(R)\sim\frac{1}{R}. \]

\(R\) en \(a_0\)	\(J\simeq1/R\) en hartree
5	0.20
10	0.10
20	0.05
50	0.02
100	0.01
1000	0.001

L’intégrale \(J\) tend donc bien vers zéro, mais très lentement. À \(R=20a_0\), elle vaut encore environ \(0.05E_h\).

Le meilleur test asymptotique n’est pas seulement de vérifier que \(J\) diminue, mais de vérifier que : \[ R\,J(R)\longrightarrow1. \] Ce test est accessible pour des valeurs comme \(R=10,\ 15,\ 20,\ 30a_0\).

6. Limite à grande distance de l’intégrale d’échange \(K\)

On considère :

\[ K=(AB|BA) = \iint \frac{ \phi_A(\mathbf r_1)\phi_B(\mathbf r_1) \phi_B(\mathbf r_2)\phi_A(\mathbf r_2) }{ r_{12} } \,d^3r_1\,d^3r_2. \]

Contrairement à \(J\), l’intégrale \(K\) dépend du recouvrement des orbitales centrées sur les deux noyaux. Pour deux orbitales \(1s\), le recouvrement vaut :

\[ S(R)=e^{-R} \left( 1+R+\frac{R^2}{3} \right). \]

\(R\)	\(S(R)\)	\(S^2(R)\)
6	\(4.71\times10^{-2}\)	\(2.22\times10^{-3}\)
8	\(1.02\times10^{-2}\)	\(1.04\times10^{-4}\)
10	\(2.01\times10^{-3}\)	\(4.05\times10^{-6}\)
12	\(3.75\times10^{-4}\)	\(1.40\times10^{-7}\)
15	\(2.78\times10^{-5}\)	\(7.75\times10^{-10}\)

L’intégrale \(K\) n’est pas exactement égale à \(S^2/R\), mais elle possède une décroissance exponentielle comparable.

vers \(R=8a_0\), \(K\) doit déjà être très petit ;
vers \(R=10a_0\), il peut être de l’ordre de quelques \(10^{-6}\) à \(10^{-5}E_h\), selon la convention exacte ;
vers \(R=12\) à \(15a_0\), il peut devenir plus petit que l’erreur numérique de la grille.

Pour tester la décroissance de \(K\), les valeurs \[ R=6,\ 8,\ 10,\ 12,\ 15\,a_0 \] sont particulièrement utiles.

7. Précaution particulière pour une convolution par FFT

Une FFT traite naturellement le tableau comme périodique. Sans remplissage par des zéros, une densité proche d’un bord peut interagir avec son image périodique située de l’autre côté de la boîte.

Le programme doit donc effectuer un zéro-padding suffisant, par exemple :

forme_doublee = (2 * Nx, 2 * Ny, 2 * Nz)

ou utiliser une convolution linéaire correctement rembourrée :

scipy.signal.fftconvolve(densite, noyau, mode="same")

Le paramètre mode="same" ne suffit pas à lui seul à garantir l’absence d’interactions périodiques : il faut vérifier la manière dont les tableaux sont prolongés et centrés avant la FFT.

Sans zéro-padding adapté, \(J\) ou \(K\) peuvent conserver une valeur parasite lorsque \(R\) devient grand.

8. Série de tests recommandée

R_TESTS = [2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 15.0, 20.0]

Pour conserver approximativement un pas constant :

PAS_CIBLE = 0.25

L = max(MARGE_BOITE, R / 2 + MARGE_NOYAU)

N_GRILLE = int(np.ceil(2 * L / PAS_CIBLE)) + 1

# On impose une grille impaire pour conserver un point central.
if N_GRILLE % 2 == 0:
    N_GRILLE += 1

Afficher systématiquement :

h = 2 * L / (N_GRILLE - 1)

print("R =", R)
print("L =", L)
print("N_GRILLE =", N_GRILLE)
print("h =", h)

print("J =", J)
print("R * J =", R * J)

print("K_fft =", K_fft)
print("K_direct =", K_direct)

Puis calculer :

erreur_absolue = abs(K_fft - K_direct)

seuil = 1.0e-14
erreur_relative = (
    abs(K_fft - K_direct)
    / max(abs(K_direct), seuil)
)

Lorsque \(K\) est presque nul, l’erreur relative peut devenir énorme ou perdre sa signification. Il faut alors examiner surtout l’erreur absolue.

9. Résultats attendus

Les comportements physiques attendus sont :

\[ R\,J(R)\longrightarrow1, \]

et :

\[ K(R)\longrightarrow0 \]

beaucoup plus rapidement, approximativement de façon exponentielle.

Résumé.

N_GRILLE = 65 est un choix acceptable et bien centré.
Pour \(R\) petit, il faut surtout vérifier le nombre de points entre les noyaux.
Pour \(R\) grand, il faut augmenter N_GRILLE afin de conserver un pas presque constant.
MARGE_NOYAU = 6.5 convient pour une première étude, mais 8 ou 10 permet une meilleure vérification des petites intégrales.
Le test asymptotique pertinent pour \(J\) est \(R\,J\to1\).
Le test pertinent pour \(K\) est sa décroissance exponentielle et la comparaison entre calcul FFT et calcul direct.
Un zéro-padding correct est indispensable pour éviter les interactions périodiques artificielles.

Toutes les distances sont exprimées en rayons de Bohr \(a_0\), et les intégrales énergétiques en hartrees \(E_h\), sauf indication contraire.