Diagnostic des intégrales \(J\) et \(K\) dans la base moléculaire \(g/u\)

Le calcul considéré utilise :

R = 40.0              # en a0
N_GRILLE = 243        # 3^5
MARGE_BOITE = 10.0
MARGE_NOYAU = 10.0

Les quatre intégrales \(J\) et \(K\) obtenues valent environ \(8{,}49259\ \mathrm{eV}\), tandis que les énergies associées aux orbitales \(g\) et \(u\) deviennent identiques à grande distance.

La valeur \(8{,}49259\ \mathrm{eV}\) est très proche de \[\frac{5}{16}E_h=0{,}3125E_h\simeq8{,}5036\ \mathrm{eV}.\] Cela suggère que le calcul conserve la moitié de l’intégrale coulombienne atomique sur un même centre, mais perd une partie de l’interaction entre les deux centres.

1. Intégrale coulombienne sur un même atome

Pour une orbitale hydrogénoïde \(1s\) normalisée :

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

Numériquement :

\[U=0{,}625E_h\simeq17{,}007\ \mathrm{eV},\qquad \frac U2=\frac{5}{16}E_h\simeq8{,}504\ \mathrm{eV}.\]

2. Ce qui doit réellement tendre vers zéro

Pour deux orbitales atomiques centrées sur les deux noyaux :

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

À grande distance :

\[J_{AB}(R)\sim\frac1R,\qquad K_{AB}(R)\longrightarrow0\]

exponentiellement. Pour \(R=40a_0\) :

\[J_{AB}\simeq\frac1{40}E_h=0{,}025E_h\simeq0{,}6803\ \mathrm{eV}.\]

Ce ne sont donc pas toutes les intégrales exprimées dans la base \(g/u\) qui doivent tendre vers zéro. Ce sont surtout \(K_{AB}\), et plus lentement \(J_{AB}\sim1/R\).

3. Orbitales moléculaires \(g\) et \(u\) à grande distance

\[\phi_g\simeq\frac{\phi_A+\phi_B}{\sqrt2},\qquad \phi_u\simeq\frac{\phi_A-\phi_B}{\sqrt2}.\] \[|\phi_g|^2\simeq|\phi_u|^2\simeq\frac12\left(|\phi_A|^2+|\phi_B|^2\right).\]

Posons :

\[U=(AA|AA)=\frac58E_h,\qquad V=(AA|BB)\simeq\frac1R.\]

À grande distance :

\[(gg|gg)\simeq(uu|uu)\simeq(gg|uu)\simeq\frac{U+V}{2}.\]

Donc :

\[J_{\mathrm{MO}}\simeq\frac{5}{16}+\frac{1}{2R}.\]

Pour \(R=40a_0\) :

\[J_{\mathrm{MO}}\simeq\frac{5}{16}+\frac1{80}=0{,}325E_h\simeq8{,}844\ \mathrm{eV}.\]

4. Intégrale d’échange entre \(g\) et \(u\)

\[\phi_g\phi_u\simeq\frac12\left(|\phi_A|^2-|\phi_B|^2\right).\]

Il en résulte :

\[K_{gu}=(gu|ug)\simeq\frac{U-V}{2}=\frac{5}{16}-\frac{1}{2R}.\]

Pour \(R=40a_0\) :

\[K_{gu}\simeq0{,}300E_h\simeq8{,}163\ \mathrm{eV}.\]

5. Pourquoi les quatre intégrales ne devraient pas être identiques

À \(R=40a_0\), on attend environ :

\[\frac{U+V}{2}\simeq8{,}844\ \mathrm{eV},\qquad \frac{U-V}{2}\simeq8{,}163\ \mathrm{eV}.\]

Leur différence vaut :

\[V\simeq\frac1R\simeq0{,}680\ \mathrm{eV}.\]

Si les quatre intégrales donnent toutes \(8{,}49259\ \mathrm{eV}\), le programme semble obtenir \(U/2\) au lieu de \((U\pm V)/2\).

Cela suggère que la contribution locale sur chaque centre est représentée, mais que la contribution coulombienne entre les deux centres, de l’ordre de \(\pm V/2\), est perdue ou annulée dans le traitement FFT.

6. La résolution spatiale semble suffisante

\[L=\frac R2+\mathrm{MARGE\_NOYAU}=20+10=30a_0,\] \[h=\frac{2L}{N_{\mathrm{GRILLE}}-1}=\frac{60}{242}\simeq0{,}248a_0.\]

Cette résolution est raisonnable pour des orbitales \(1s\). Le problème semble donc davantage lié à la convolution FFT qu’au pas de grille.

7. Causes numériques probables

centrage du noyau coulombien \(1/r\) ;
utilisation cohérente de fftshift et ifftshift ;
zéro-padding ;
découpage du résultat après convolution ;
distance maximale représentée par le noyau ;
facteurs de volume \(h^3\).

7.1. Centrage du noyau coulombien

xk = (np.arange(Npad) - Npad // 2) * h
Xk, Yk, Zk = np.meshgrid(xk, xk, xk, indexing="ij")
r = np.sqrt(Xk**2 + Yk**2 + Zk**2)
noyau = np.zeros_like(r)
masque = r > 0.0
noyau[masque] = 1.0 / r[masque]
noyau_fft = np.fft.fftn(np.fft.ifftshift(noyau))

Si ifftshift manque, l’origine physique du noyau \(1/r\) ne correspond pas à l’origine attendue par la FFT.

7.2. Zéro-padding

Pour une convolution linéaire entre deux tableaux de taille \(N\), il faut une taille au moins égale à :

\[2N-1.\]

from scipy.fft import next_fast_len
Npad = next_fast_len(2 * N_GRILLE - 1)

Pour \(N_{\mathrm{GRILLE}}=243\) :

\[2N_{\mathrm{GRILLE}}-1=485.\]

Si la FFT est faite directement sur un tableau \(243^3\), elle réalise une convolution circulaire.

8. Test décisif avec les densités atomiques

rho_A = phi_A**2
rho_B = phi_B**2

U_A = interaction_coulombienne(rho_A, rho_A)
U_B = interaction_coulombienne(rho_B, rho_B)
V_AB = interaction_coulombienne(rho_A, rho_B)

À \(R=40a_0\), on attend :

\[U_A\simeq U_B\simeq\frac58E_h=0{,}625E_h\simeq17{,}007\ \mathrm{eV},\] \[V_{AB}\simeq\frac1{40}E_h=0{,}025E_h\simeq0{,}6803\ \mathrm{eV}.\]

print("U_A =", U_A)
print("U_B =", U_B)
print("V_AB =", V_AB)
print("1 / R =", 1.0 / R)
print("R * V_AB =", R * V_AB)

\[R\,V_{AB}\longrightarrow1.\]

Si le programme donne \(V_{AB}\simeq0\), l’erreur est localisée : la convolution FFT ne restitue pas correctement l’interaction entre deux densités éloignées.

9. Singlet et triplet \(gu\)

\[E_{\mathrm{sing}}=h_g+h_u+J_{gu}+K_{gu},\] \[E_{\mathrm{trip}}=h_g+h_u+J_{gu}-K_{gu}.\]

Leur différence vaut :

\[E_{\mathrm{sing}}-E_{\mathrm{trip}}=2K_{gu}.\]

Il ne faut cependant pas attendre que \(K_{gu}\) tende vers zéro, car les orbitales \(g\) et \(u\) restent délocalisées :

\[K_{gu}\longrightarrow\frac{5}{16}E_h\]

à une correction en \(1/R\) près.

Cela ne signifie pas que le vrai singulet et le vrai triplet physiques de deux atomes d’hydrogène restent séparés de près de \(17\ \mathrm{eV}\). Les configurations simples dans la base \(g/u\) ne décrivent pas correctement toutes les limites de dissociation.

Pour la bonne dissociation du singulet fondamental de \(\mathrm H_2\), il faut mélanger au moins :

\[|\sigma_g^2\rangle\qquad\text{et}\qquad|\sigma_u^2\rangle.\]

10. Pourquoi les courbes \(g\) et \(u\) se confondent

\[S\to0,\qquad H_{AB}\to0,\qquad E_g-E_u\to0.\]

Une coïncidence visuelle à partir de \(R\simeq5{,}86\ \text{Å}\), soit environ \(11{,}1a_0\), est plausible selon l’échelle verticale.

La valeur \(-28{,}57196\ \mathrm{eV}\) est proche de :

\[2\left(-\frac12-\frac1R\right)E_h\]

pour \(R=40a_0\) :

\[2(-0{,}5-0{,}025)E_h=-1{,}05E_h\simeq-28{,}572\ \mathrm{eV}.\]

La courbe appelée « énergie d’orbitale » représente donc peut-être en réalité \(2h_g\) ou \(2h_u\), et non l’énergie d’une seule orbitale.

11. Absence de minimum pour \(\sigma_g^2\)

Il faut vérifier que la répulsion internucléaire est bien ajoutée :

\[E_{\mathrm{noyaux}}=\frac1R E_h.\]

En électronvolts, si \(R\) est exprimé en \(a_0\) :

E_repulsion_noyaux_eV = EV_PAR_HARTREE / R

Il ne faut pas utiliser directement :

EV_PAR_HARTREE / R_angstrom

sans conversion cohérente des unités.

12. Bloc de diagnostic minimal

rho_A = phi_A**2
rho_B = phi_B**2

U_A = coulomb_fft(rho_A, rho_A)
U_B = coulomb_fft(rho_B, rho_B)
V_AB = coulomb_fft(rho_A, rho_B)

print()
print("Diagnostic des intégrales atomiques")
print("-----------------------------------")
print("U_A en Eh =", U_A)
print("U_B en Eh =", U_B)
print("Valeur attendue de U =", 5.0 / 8.0)
print("V_AB en Eh =", V_AB)
print("Valeur attendue 1/R =", 1.0 / R)
print("R * V_AB =", R * V_AB)

Grandeur	Valeur attendue en \(E_h\)	Valeur attendue en eV
\(U_A\)	0,625	environ 17,007
\(U_B\)	0,625	environ 17,007
\(V_{AB}\)	0,025	environ 0,6803
\(R V_{AB}\)	environ 1	sans objet

13. Conclusion

La valeur \(8{,}49259\ \mathrm{eV}\) est proche de \(5/16E_h\).
Les intégrales dans la base \(g/u\) ne doivent pas toutes tendre vers zéro.
À \(R=40a_0\), on attend deux familles proches de \(8{,}844\) et \(8{,}163\ \mathrm{eV}\), et non quatre valeurs strictement égales.
Le pas \(h\simeq0{,}248a_0\) est raisonnable ; le problème semble plutôt lié à la convolution FFT.
Le test prioritaire est \(V_{AB}\simeq1/R\), donc \(R V_{AB}\simeq1\).
Il faut vérifier le zéro-padding, le centrage du noyau \(1/r\), les décalages FFT et le découpage du résultat.
La répulsion internucléaire \(1/R\) doit être ajoutée avec des unités cohérentes.

Les distances sont exprimées en rayons de Bohr \(a_0\), les énergies en hartrees \(E_h\) ou en électronvolts selon l’indication.