Méthodes analytiques et numériques de molecule_h2_01_03_integrales_2e.py
- des expressions analytiques pour les grandeurs à un électron ;
- un calcul numérique tridimensionnel par convolution FFT pour les intégrales à deux électrons.
1. Objectif général du programme
Le programme calcule, pour plusieurs distances internucléaires \(R\), les quatre intégrales :
\[ J_{gg}=(gg|gg), \qquad J_{gu}=(gg|uu), \qquad J_{uu}=(uu|uu), \qquad K_{gu}=(gu|ug). \]Il construit ensuite les quatre énergies :
\[ E_{gg}^{S}, \qquad E_{gu}^{S}, \qquad E_{gu}^{T}, \qquad E_{uu}^{S}. \]Les résultats sont affichés, sauvegardés dans un fichier CSV puis tracés.
2. Constantes physiques et unités
Le programme commence par reconstruire explicitement :
\[ a_0 = \frac{4\pi\varepsilon_0\hbar^2}{m_e e^2}, \] \[ E_h = \frac{e^2}{4\pi\varepsilon_0a_0}. \]Il en déduit :
\[ 1\,E_h = \texttt{EV\_PAR\_HARTREE}\ \text{eV}, \] \[ 1\,a_0 = \texttt{ANGSTROM\_PAR\_A0}\ \text{\AA}. \]En revanche, tous les calculs internes sont faits avec :
\[ a_0=1, \qquad E_h=1. \]3. Distances et paramètres numériques
Les distances \(R\) étudiées sont imposées par le tableau :
R_VALEURS_A0 = np.array([
0.80, 1.00, 1.20, 1.40, 1.60,
2.00, 2.50, 3.00, 4.00,
5.00, 6.00, 8.00
])La grille tridimensionnelle contient :
\[ N_{\mathrm{grille}}^3 \]points, avec ici :
\[ N_{\mathrm{grille}}=65. \]Le programme prévoit explicitement des études de convergence avec :
\[ N_{\mathrm{grille}}=49,\ 65,\ 81. \]4. Partie analytique à un électron
4.1 Recouvrement des orbitales \(1s\)
Les orbitales atomiques \(1s\), centrées sur les noyaux \(A\) et \(B\), sont :
\[ \phi_A(\mathbf r) = \frac{e^{-r_A}}{\sqrt{\pi}}, \qquad \phi_B(\mathbf r) = \frac{e^{-r_B}}{\sqrt{\pi}}. \]Le recouvrement est calculé analytiquement :
\[ S(R) = \langle\phi_A|\phi_B\rangle = e^{-R} \left( 1+R+\frac{R^2}{3} \right). \]Cette expression est utilisée dans la fonction :
def recouvrement_1s(R):
return np.exp(-R) * (1.0 + R + R**2 / 3.0)4.2 Intégrale coulombienne à un électron
Le programme utilise l’expression analytique :
\[ J_{1e}(R) = -\frac{1}{R} + e^{-2R} \left( 1+\frac{1}{R} \right). \]Cette quantité correspond à l’interaction d’un électron localisé sur un noyau avec l’autre noyau.
4.3 Intégrale de résonance à un électron
L’intégrale de résonance est calculée par :
\[ K_{1e}(R) = -e^{-R}(1+R). \]integrale_resonance_1s, ne doit pas être confondue avec l’intégrale d’échange à deux électrons \(K_{gu}\).
4.4 Éléments de matrice \(H_{AA}\) et \(H_{AB}\)
Le programme écrit :
\[ H_{AA} = E_{1s}^{\mathrm H} + J_{1e}(R), \] avec : \[ E_{1s}^{\mathrm H} = -\frac12 E_h. \]Puis :
\[ H_{AB} = E_{1s}^{\mathrm H}S(R) + K_{1e}(R). \]4.5 Énergies orbitalaires \(g\) et \(u\)
Les deux valeurs propres LCAO sont calculées analytiquement :
\[ \varepsilon_g(R) = \frac{H_{AA}+H_{AB}}{1+S}, \] \[ \varepsilon_u(R) = \frac{H_{AA}-H_{AB}}{1-S}. \]5. Construction de la grille cartésienne
5.1 Position des noyaux
Les deux noyaux sont placés sur l’axe \(z\) :
\[ z_A=-\frac{R}{2}, \qquad z_B=+\frac{R}{2}. \]5.2 Taille de la boîte
La demi-largeur de la boîte est définie par :
\[ L = \max \left( L_{\mathrm{fixe}}, \frac{R}{2}+L_{\mathrm{noyau}} \right), \] avec : \[ L_{\mathrm{fixe}}=8\,a_0, \qquad L_{\mathrm{noyau}}=6{,}5\,a_0. \]Cela garantit que chaque noyau reste entouré d’une marge suffisante, même lorsque \(R\) augmente.
5.3 Pas de grille
Sur l’intervalle \([-L,L]\), le pas est :
\[ \Delta x = \frac{2L}{N_{\mathrm{grille}}-1}. \]La grille est isotrope :
\[ \Delta x=\Delta y=\Delta z. \]L’élément de volume discret est :
\[ \Delta V = (\Delta x)^3. \]6. Orbitales atomiques sur la grille
À chaque point \((x,y,z)\), le programme calcule :
\[ r_A = \sqrt{x^2+y^2+\left(z+\frac{R}{2}\right)^2}, \] \[ r_B = \sqrt{x^2+y^2+\left(z-\frac{R}{2}\right)^2}. \]Puis :
\[ \phi_A = \frac{e^{-r_A}}{\sqrt{\pi}}, \qquad \phi_B = \frac{e^{-r_B}}{\sqrt{\pi}}. \]7. Orbitales moléculaires sur la grille
Les orbitales moléculaires sont construites avec la valeur analytique du recouvrement :
\[ \phi_g = \frac{\phi_A+\phi_B} {\sqrt{2(1+S)}}, \] \[ \phi_u = \frac{\phi_A-\phi_B} {\sqrt{2(1-S)}}. \]Le programme contrôle numériquement :
\[ \int |\phi_g|^2\,d^3r, \] \[ \int |\phi_u|^2\,d^3r, \] \[ \int \phi_g\phi_u\,d^3r. \]Ces trois quantités sont évaluées par des sommes discrètes :
\[ \int f(\mathbf r)\,d^3r \approx \sum_{ijk}f_{ijk}\,\Delta V. \]Les résultats sont enregistrés dans :
norme_g
norme_u
orthogonalite_gu8. Densités utilisées pour les intégrales à deux électrons
Le programme construit trois densités :
\[ \rho_g(\mathbf r)=\phi_g^2(\mathbf r), \] \[ \rho_u(\mathbf r)=\phi_u^2(\mathbf r), \] \[ \rho_{gu}(\mathbf r)=\phi_g(\mathbf r)\phi_u(\mathbf r). \]Les deux premières sont des densités électroniques positives. La troisième est une densité de transition qui change de signe.
9. Réduction d’une intégrale double à deux intégrales simples
Une intégrale à deux électrons s’écrit :
\[ I = \iint \frac{ \rho_1(\mathbf r_1) \rho_2(\mathbf r_2) }{ |\mathbf r_1-\mathbf r_2| } \,d^3r_1\,d^3r_2. \]Le programme commence par calculer le potentiel créé par \(\rho_2\) :
\[ V_{\rho_2}(\mathbf r) = \int \frac{ \rho_2(\mathbf r') }{ |\mathbf r-\mathbf r'| } \,d^3r'. \]Puis l’intégrale devient :
\[ I = \int \rho_1(\mathbf r) V_{\rho_2}(\mathbf r) \,d^3r. \]10. Convolution coulombienne par FFT
10.1 Forme de convolution
Le potentiel coulombien est une convolution :
\[ V_\rho = \rho * \frac{1}{r}. \]Sous forme discrète :
\[ V_{ijk} \approx \sum_{i'j'k'} \rho_{i'j'k'} \frac{1}{ |\mathbf r_{ijk}-\mathbf r_{i'j'k'}| } \Delta V. \]10.2 Utilisation de fftconvolve
Le programme appelle :
fftconvolve(
rho,
noyau,
mode="same"
) * pas**3La FFT utilise le théorème de convolution :
\[ \mathcal F(\rho*K) = \mathcal F(\rho)\, \mathcal F(K). \]La convolution directe aurait un coût très élevé. Pour une grille contenant :
\[ N^3 \]points, la méthode FFT réduit fortement le temps de calcul par rapport à une somme directe sur toutes les paires de points.
10.3 Mode same
L’option :
mode="same"renvoie un tableau de même taille que la densité d’entrée.
11. Traitement de la singularité du noyau \(1/r\)
Le noyau coulombien discret est :
\[ K(\mathbf r)=\frac{1}{r}. \]Au centre de la grille :
\[ r=0, \]et la valeur \(1/r\) diverge.
Le programme ne met ni zéro ni une valeur arbitraire. Il remplace la cellule cubique centrale par une sphère de même volume.
Si \(a\) est le rayon de cette sphère :
\[ \frac{4\pi a^3}{3} = (\Delta x)^3. \]D’où :
\[ a = \left( \frac{3(\Delta x)^3}{4\pi} \right)^{1/3}. \]La moyenne de \(1/r\) dans cette sphère vaut :
\[ \left\langle\frac1r\right\rangle = \frac{ \int_0^a \frac1r 4\pi r^2\,dr }{ \frac43\pi a^3 } = \frac{3}{2a}. \]La valeur centrale du noyau est donc remplacée par :
\[ K(0) = \frac{3}{2a}. \]12. Potentiels calculés
Le programme ne calcule que trois potentiels :
\[ V_g = \rho_g * \frac1r, \] \[ V_u = \rho_u * \frac1r, \] \[ V_{gu} = \rho_{gu} * \frac1r. \]Cela suffit pour obtenir les quatre intégrales recherchées.
13. Calcul numérique des intégrales \(J\) et \(K\)
13.1 Intégrale \(J_{gg}\)
\[ J_{gg} = \int \rho_g(\mathbf r) V_g(\mathbf r) \,d^3r. \]Numériquement :
\[ J_{gg} \approx \sum_{ijk} \rho_{g,ijk} V_{g,ijk} \Delta V. \]13.2 Intégrale \(J_{uu}\)
\[ J_{uu} = \int \rho_u(\mathbf r) V_u(\mathbf r) \,d^3r. \]13.3 Intégrale \(J_{gu}\)
Deux formes théoriquement identiques sont calculées :
\[ J_{gu}^{(1)} = \int \rho_g V_u\,d^3r, \] \[ J_{gu}^{(2)} = \int \rho_u V_g\,d^3r. \]En théorie :
\[ J_{gu}^{(1)} = J_{gu}^{(2)}. \]Sur une grille finie, une petite différence peut apparaître. Le programme prend donc :
\[ J_{gu} = \frac12 \left( J_{gu}^{(1)} + J_{gu}^{(2)} \right). \]13.4 Intégrale d’échange \(K_{gu}\)
\[ K_{gu} = \int \rho_{gu}(\mathbf r) V_{gu}(\mathbf r) \,d^3r. \]avec :
\[ \rho_{gu}=\phi_g\phi_u. \]Numériquement :
\[ K_{gu} \approx \sum_{ijk} \rho_{gu,ijk} V_{gu,ijk} \Delta V. \]14. Contrôle analytique sur l’atome d’hydrogène
Avant de calculer les intégrales moléculaires, le programme teste la méthode sur une orbitale \(1s\) isolée :
\[ \phi_{1s}(r) = \frac{e^{-r}}{\sqrt{\pi}}. \]L’intégrale exacte connue est :
\[ (1s\,1s|1s\,1s) = \frac58 E_h = 0{,}625\,E_h. \]Le programme calcule :
\[ J_{\mathrm{num}} = \iint \frac{ |\phi_{1s}(\mathbf r_1)|^2 |\phi_{1s}(\mathbf r_2)|^2 }{ r_{12} } \,d^3r_1\,d^3r_2. \]Puis il affiche :
\[ \Delta J = J_{\mathrm{num}}-\frac58, \] \[ \varepsilon_{\mathrm{rel}} = \frac{ J_{\mathrm{num}}-\frac58 }{ \frac58 }. \]15. Construction analytique des énergies totales
Une fois \(J_{gg}\), \(J_{gu}\), \(J_{uu}\) et \(K_{gu}\) calculés numériquement, les énergies sont obtenues algébriquement.
15.1 Configuration \((1\sigma_g)^2\)
\[ E_{gg}^{S} = 2\varepsilon_g + J_{gg} + \frac1R. \]15.2 Configuration \(1\sigma_g1\sigma_u\), singulet
\[ E_{gu}^{S} = \varepsilon_g+\varepsilon_u + J_{gu}+K_{gu} + \frac1R. \]15.3 Configuration \(1\sigma_g1\sigma_u\), triplet
\[ E_{gu}^{T} = \varepsilon_g+\varepsilon_u + J_{gu}-K_{gu} + \frac1R. \]15.4 Configuration \((1\sigma_u)^2\)
\[ E_{uu}^{S} = 2\varepsilon_u + J_{uu} + \frac1R. \]L’écart singulet–triplet est aussi calculé :
\[ E_{gu}^{S} - E_{gu}^{T} = 2K_{gu}. \]16. Sauvegarde des résultats
Pour chaque distance \(R\), le programme sauvegarde notamment :
| Grandeur | Signification |
|---|---|
R_a0 | Distance en rayons de Bohr |
R_A | Distance en ångströms |
pas_a0 | Pas spatial de la grille |
demi_largeur_a0 | Demi-largeur de la boîte |
norme_g | Norme numérique de \(\phi_g\) |
norme_u | Norme numérique de \(\phi_u\) |
orthogonalite_gu | Valeur numérique de \(\langle g|u\rangle\) |
J_gg, J_gu, J_uu | Intégrales coulombiennes |
K_gu | Intégrale d’échange |
epsilon_g, epsilon_u | Contributions à un électron |
V_pp | Répulsion proton–proton \(1/R\) |
E_... | Énergies totales des configurations |
Le fichier produit est :
molecule_h2_01_03_integrales_2e.csv17. Figures diagnostiques
Le programme produit deux figures :
- les quatre intégrales \(J_{gg}\), \(J_{gu}\), \(J_{uu}\) et \(K_{gu}\) en fonction de \(R\) ;
- les quatre énergies totales en fonction de \(R\).
Les intégrales et les énergies sont converties de hartrees en eV uniquement pour ces tracés.
18. Sources principales d’erreur numérique
| Source d’erreur | Effet possible | Contrôle |
|---|---|---|
| Pas de grille trop grand | Mauvaise description près des noyaux et du noyau \(1/r\) | Comparer \(N=49,65,81\) |
| Boîte trop petite | Troncature des queues exponentielles | Augmenter les marges |
| Régularisation centrale | Erreur locale sur la singularité | Contrôle \(5/8\) |
| Annulations dans \(\rho_{gu}\) | Sensibilité accrue de \(K_{gu}\) | Étude de convergence et méthode indépendante |
| Variation de la boîte avec \(R\) | Variation du pas si \(N\) reste fixe | Enregistrer pas_a0 |
19. Nature exacte du modèle
- deux orbitales atomiques \(1s\) fixes ;
- pas d’optimisation de l’exposant orbital ;
- pas de relaxation Hartree–Fock ;
- pas de corrélation électronique complète ;
- pas d’interaction de configurations dans ce programme.
20. Résumé du déroulement du calcul
- Choisir une distance \(R\).
- Calculer analytiquement \(S\), \(\varepsilon_g\) et \(\varepsilon_u\).
- Construire la grille cartésienne.
- Évaluer \(\phi_A\) et \(\phi_B\).
- Construire \(\phi_g\) et \(\phi_u\).
- Contrôler leurs normes et leur orthogonalité.
- Construire \(\rho_g\), \(\rho_u\) et \(\rho_{gu}\).
- Construire le noyau discret régularisé \(1/r\).
- Calculer les trois potentiels par convolution FFT.
- Calculer \(J_{gg}\), \(J_{gu}\), \(J_{uu}\) et \(K_{gu}\).
- Ajouter les contributions analytiques à un électron et \(1/R\).
- Former les quatre énergies totales.
- Sauvegarder et tracer les résultats.