Deep Learning Chimie quantique Graph Neural Networks QSAR · QSPR

PhotoGNN

Pipeline d'encodeurs moléculaires multi-tâche pour la photocatalyse organique. Deux backbones hybrides MACE–GINEConv–SchNet entraînés sur 443k molécules TD-DFT (QCDGE) et 15k potentiels d'oxydation expérimentaux (OxPot), fusionnés dans une architecture en goulot qui prédit simultanément 6 propriétés. Précisions finales (erreurs de prédiction moyennes) : gap HOMO-LUMO 47.9 meV · E_S₁ 43.1 meV · E_T₁ 32.8 meV · E_HOMO 81.7 meV · Eox 54.0 mV.
Ce projet est pensé comme la première étape d'un jumeau numérique de réacteur photochimique en flux continu.

PyTorchRDKitOptuna

I. Propriétés photophysiques : généralités et intérêt

La photocatalyse rédox repose sur un principe fondamental : convertir l’énergie lumineuse en énergie chimique, là où la thermochimie classique exige de la chaleur ou des métaux précieux et/ou toxiques.

Dans ce cadre, le photocatalyseur joue un rôle central via un état excité qui est simultanément un meilleur oxydant et un meilleur réducteur que l’état fondamental.

Prédire le gap HOMO-LUMO seul ne suffit pas. Ce qui gouverne l’efficacité d’un photocatalyseur organique, c’est l’ensemble de ses propriétés à l’état excité :

Propriété Symbole Ce qu’elle contrôle
Gap HOMO-LUMO $E_{gap}$ Longueur d’onde d’absorption, réactivité à l’état excité
Énergie du 1er état singulet $E_{S_1}$ Énergie disponible pour le transfert d’électron
Énergie du 1er état triplet $E_{T_1}$ Probabilité et énergie du croisement intersystème (ISC)
Force d’oscillateur $S_0 \to S_1$ $f_{S_1}$ Intensité d’absorption, rendement de photoexcitation
Rendement quantique $\Phi$ Efficacité de la conversion photons $\to$ produits (ou émission)
Durée de vie $T_1$ $\tau_{T_1}$ Fenêtre temporelle pour la diffusion et la réactivité bimoléculaire
Potentiels rédox fondamental $E_{red}/E_{ox}$ Stabilité thermodynamique et point de référence pour l’excitation
Potentiel rédox excité $E^*$ Pouvoir réducteur/oxydant sous lumière (calculé via Rehm-Weller)

Sachant que :

  1. Potentiel rédox à l’état excité (Approximation de Rehm-Weller) :

    \[E_{ox}^* \approx E_{ox} - E_{0,0}\]

    Où $E_{0,0}$ est l’énergie de transition pure, souvent estimée par $E_{S_1}$ ou $E_{T_1}$ selon le mécanisme.

  2. Rendement quantique ($\Phi$) :

    \[\Phi = \frac{k_{obs}}{k_{rad} + k_{non-rad} + k_{quench}}\]

    Le contrôle de $\Phi$ dépend de la compétition entre les constantes de vitesse de désactivation.

Concevoir ou sélectionner un photocatalyseur organique adapté à une réaction cible suppose donc de connaître ces propriétés avec précision, idéalement avant de synthétiser la molécule.

Obtenir ces grandeurs nécessite des calcules par TD-DFT, prennant de plusieurs minutes à plusieurs heures par molécule (avec un coût computationnel croîssant à une puissance élevée du nombre d’atomes), ou des mesures expérimentales couteuses en temps et en réactifs, et évidemment nécessitant les appareils de mesure.

Pour explorer un espace chimique de millions de candidats, ces limitations sont prohibitives.

L’objectif de PhotoGNN est de prédire ces propriétés simultanément en quelques millisecondes, avec une précision proche du calcul TD-DFT.


Le modèle conçu, nommé dans la suite du document “méta-modèle”, se repose sur l’assemblage en parallèle de deux modèles prédictifs séparés en deux branches, respectivements spécialisés dans la prédiction de propriétés photophysiques et des propriétés rédox.

Ce choix de conception a été motivé par plusieurs raisons :

  1. deux branches = plus simple à optimiser sélectivement
  2. l’entrainement du modèle “de tête”, servant d’entonoir, permet d’améliorer le résultat des deux modèles via {raison liée à la loss fonction}
  3. {raison 3}

II. Branche photophysique : modèle prédictif de propriétés à l’état excité

II. A. Preuve de concept : validation sur QM9

Avant d’attaquer QCDGE, l’architecture a été validée sur QM9 (Ramakrishnan et al., 2014) : 133 885 molécules organiques avec gap HOMO-LUMO calculé par B3LYP/6-31G(2df,p). Ce POC a permis de confirmer la valeur de chaque brique architecturale par ablation progressive, et d’identifier le plateau de cette combinaison dataset/méthode. Le checkpoint résultant (90 meV sur le gap) sert de warmstart pour QCDGE.

Étape Architecture Erreur moyenne (Gap HOMO-LUMO, meV)
1a Baseline SVR / Random Forest / MLP 393 / 359 / 348
2 GNN topologique RDKit (GINE/GAT/SAGE) 212
3 GNN RDKit + SchNet 157
4 MACE gelé + pooling moyen + MLP 126
5 MACE gelé + GINEConv per-atome + RDKit 95
6 Architecture hybride 4 branches 90

Chaque étape a révélé la limite de la précédente : les architectures de machine learning classique et réseaux de neurones simples, entrée 1, montrent une précision non utilisable.

L’utilisation de réseaux de neurones graphiques avec des features moléculaires calculées avec RDKit (entrée 2) permet un gain de performence significatif. L’ajout des informations spatiales avec l’incorporation de SchNet montre aussi une amélioration de la performence (entrée 3).

L’incorporation de MACE, modèle fondation entrainé sur le dataset OMOL,[ref] permet d’apporter des informations physiques atomiques résultant sur le saut de précision le plus important. L’entrée 5, qui a permis de passer sous la barre des 100 meV d’erreur prédictive, est obtenu en utilisant une structure à 4 branches : le pooling global perd l’information différentielle inter-atomes, et le parallelisme de branche apporte une complémentarité que chacune seule ne peut fournir. La dernière architecture se conclut par un plateau à 90 meV après une optimisation doussée du modèle, les limites du dataset QM9 et le besoin de plus de données.

II. B. Dataset QCDGE, 443 000 molécules

QCDGE (Zhu et al., 2024) est un dataset de 443 106 molécules organiques (CNOF, ≤10 atomes lourds), calculées au niveau TD-ωB97X-D/6-31G* avec correction de dispersion D3BJ pour l’état fondamental. C’est la base de données qui pilote l’entraînement principal.

Par rapport à QM9 (le dataset de référence classique, utilisé comme POC — voir section suivante), QCDGE apporte trois avantages décisifs :

La répartition des données suit un scaffold split de Murcko (80/10/10), assurant une réparition proportionnelle de la taille des molécules au sein des différents splits.

Alignement DFT des embeddings MACE. Les graphes PyG sont reconstruits depuis les coordonnées DFT extraites du HDF5 (via rdDetermineBonds), garantissant que l’ordre atomique MACE coïncide exactement avec l’ordre du graphe moléculaire. L’information per-atome de MACE est exploitée localement, et non poolée globalement. 98,1 % des molécules bénéficient de cet alignement ; les 1,9 % restants utilisent un fallback sur la géométrie ETKDG.

Les embeddings MACE ont été pré-calculés (~443 000 molécules, 101 Go au total) puis stockés en fichiers memmap numpy pour optimiser l’utilisation de la mémoire système, et ainsi les partager entre les workers lors de l’optimisation des hyperparamètres sans duplication en mémoire vive.


II. C. Modèle fondation MACE comme encodeur atomistique pré-entraîné

Le saut le plus important, de 157 meV à 126 meV puis 95 meV, s’est produit en introduisant MACE-OMOL comme encodeur atomique gelé.

MACE est un modèle atomistique pré-entraîné sur 100 millions de conformères de petites molécules organiques (dataset OMol25, Meta/FAIR, 2024–2025), avec des calculs DFT de haute précision (ωB97M-V). Son architecture exploite des harmoniques sphériques d’ordre élevé pour encoder la géométrie 3D locale de chaque atome de manière équivariante (c’est-à-dire de façon cohérente quelle que soit l’orientation de la molécule dans l’espace).

Appliqué à chaque atome d’une molécule, MACE produit un vecteur de 3 072 dimensions capturant son environnement chimique local jusqu’à plusieurs angströms de rayon, en intégrant explicitement les distances, angles et torsions.

L’idée centrale est de ne pas ré-apprendre ce que MACE sait déjà. Le modèle pré-entraîné est gelé ; seul un réseau aval (plus petit, plus rapide, ciblé sur les propriétés photophysiques) est entraîné. Cette stratégie de frozen encoder permet de transférer une connaissance atomistique apprise sur des centaines de millions de points vers une tâche spécifique avec beaucoup moins de données.


II. D. Architecture parallèle à quatre branches

Le cœur du modèle repose sur quatre branches parallèles, chacune exploitant une source d’information distincte sur la même molécule :

MACE per-atome [N, 3072]  ──►  Branche GINEConv   →  interactions covalentes      ─┐
MACE per-atome [N, 3072]  ──►  Branche SchNet      →  géométrie 3D continue       ─┤─► fusion → [gap, E_S₁, E_T₁, f_S₁]
Graphe topologique [N,15D]──►  Branche TopoGNN     →  connectivité pure           ─┤
Descripteurs RDKit [29D]  ──►  Branche RDKit MLP   →  propriétés globales         ─┘

Sous-branche GINEConv : interactions covalentes avec les embeddings MACE

Les embeddings MACE per-atome sont projetés dans un espace plus compact, puis traités par des couches GINEConv (Graph Isomorphism Network with Edge features) qui propagent l’information le long des liaisons covalentes de la molécule. Les features des liaisons (type SINGLE/DOUBLE/TRIPLE/AROMATIQUE, conjugaison) sont explicitement intégrées à chaque étape de message passing. Un mécanisme résiduel préserve l’information à travers les couches.

Cette branche capture les effets de délocalisation électronique, les conjugaisons, les systèmes aromatiques — les grandes déterminantes du gap et des niveaux d’énergie excités.

Sous-branche SchNet : géométrie 3D continue

Parallèlement, les mêmes embeddings MACE per-atome sont traités par un réseau de type SchNet (Schütt et al., NeurIPS 2017), qui opère sur un graphe spatial construit à partir des distances interatomiques. Toutes les paires d’atomes séparées par moins de 4 Å sont connectées ; la distance est encodée dans une base de fonctions gaussiennes radiales.

SchNet utilise des filtres continus : pour chaque interaction, un petit réseau calcule un filtre dépendant de la distance, qui module le message passé entre atomes. Cette architecture encode les interactions à moyenne portée non capturées par les liaisons covalentes seules.

Le graphe spatial (radius_graph 3D) est pré-calculé une seule fois dans le cache memmap numpy — aucun recalcul de KD-tree pendant l’entraînement.

Sous-branche TopoGNN : topologie pure

Une troisième branche opère directement sur les features atomiques brutes (type d’atome, degré, charge formelle, aromaticité, appartenance à un cycle : total de 15 dimensions), sans passer par MACE. Elle utilise plusieurs sous-réseaux GNN parallèles dont les sorties sont moyennées, et joue le rôle de régularisateur topologique complémentaire aux branches MACE.

Sous-branche RDKit : descripteurs moléculaires globaux

29 descripteurs calculés par RDKit encodent des propriétés globales de la molécule : charges partielles de Gasteiger, masse moléculaire, aromaticité, logP, TPSA, nombre de cycles, stéréocentres, indices de complexité. Ces descripteurs capturent ce que le message passing ne peut pas facilement exprimer — des propriétés qui dépendent de la molécule entière plutôt que d’interactions locales.

Fusion et tête de prédiction

Les quatre sorties sont concaténées dans un vecteur de fusion, puis traitées par une tête à trois couches linéaires pour produire les quatre propriétés prédites simultanément. La dénormalisation finale (z-score inversé appris sur le train set) donne les résultats en eV (ou sans unité pour $f_{S_1}$).

[gine_out | schnet_out | topo_out | rdkit_out]
            ↓
    Linear → LayerNorm → ELU → Dropout
    Linear → ELU
    Linear → [E_gap, E_S1, E_T1, f_S1]

Choix techniques structurants

LayerNorm plutôt que BatchNorm. La normalisation par batch est instable avec des graphes moléculaires de tailles très variables : les statistiques running calculées sur un batch de graphes hétérogènes ne convergent pas, et le mode évaluation devient pathologique. LayerNorm normalise chaque sample indépendamment — comportement identique en entraînement et en évaluation.

Huber Loss masquée par cible. La loss Huber est robuste aux outliers (grands systèmes conjugués, hétéroatomes multiples). Sur QCDGE, certaines molécules ne disposent pas de calcul TD-DFT valide pour toutes les cibles : la loss est masquée par cible, excluant les NaN sans supprimer la molécule. Chaque cible contribue proportionnellement à la disponibilité des données.

Scaffold split de Murcko avec correction grands groupes. Les splits aléatoires contaminent le jeu de test avec des analogues proches des molécules d’entraînement. Sur QCDGE, 92 % des molécules sont acycliques (scaffold ""), ce qui crée un déséquilibre massif avec un split naïf. La solution : distribuer proportionnellement les grands groupes de scaffold entre les trois splits, garantissant une répartition 80/10/10 sur l’ensemble du dataset.

Optimisation bayésienne des hyperparamètres (Optuna). L’espace de recherche couvre +40 hyperparamètres (dimensions des branches, nombre de couches, dropout, learning rate, scheduler, grad_clip, batch size). L’optimisation tourne en parallèle sur deux GPU via deux workers indépendants, chacun responsable d’un trial complet sur un GPU dédié.


II. E. Résultats

Le modèle final, entraîné sur 300 époques avec les hyperparamètres optimaux (~4,5 M de paramètres), atteint les résultats suivants sur le jeu de test QCDGE :

Métrique Standalone final
Score combiné (mean MAE/std) 0.0339
$E_{gap}$ MAE 47.9 meV
$E_{S_1}$ MAE 41.8 meV
$E_{T_1}$ MAE 32.1 meV
$f_{S_1}$ MAE 8.0

Sur la qualité des labels. QM9 est calculé au niveau B3LYP, connu pour sous-estimer le gap HOMO-LUMO de plusieurs centaines de meV. QCDGE utilise TD-ωB97X-D, une fonctionnelle à séparation de portée corrigeant ces biais sur les états excités, et le seul niveau de calcul adapté à la prédiction de E_S₁ et E_T₁. La précision d’un modèle est toujours plafonnée par la méthode qui a généré ses labels : le passage de B3LYP à TD-ωB97X-D est donc une progression fondamentale, indépendamment de l’architecture.

Ces résultats sont obtenus dans le domaine de QCDGE : des molécules organiques de 1 à 10 atomes lourds (C, N, O, F). C’est là que s’arrête la validité directe du modèle … et que commence le vrai défi !


III. Branche redox : potentiels d’oxydation

Prédire les propriétés à l’état excité ne suffit pas. Dans un transfert d’électron photoinduit, c’est le potentiel d’oxydation à l’état excité ($E_{ox}^*$) qui détermine si la réaction est thermodynamiquement favorable — via la relation de Rehm-Weller :

\[E^*_{ox} = E_{ox} - E_{S_1}\]

Un second backbone indépendant est donc entraîné sur deux nouvelles cibles :

Propriété Symbole Ce qu’elle contrôle
Énergie de l’orbitale HOMO $E_{HOMO}$ Niveau d’énergie du donneur d’électrons à l’état fondamental
Potentiel d’oxydation $E_{ox}$ Mesure expérimentale directe de la capacité réductrice

III. A. Dataset OxPot

OxPot (Goldammer et al., JCIM 2025) regroupe 15 238 molécules organiques avec leur énergie HOMO calculée par DFT (PBE0/cc-pVDZ) et leur potentiel d’oxydation mesuré expérimentalement (V vs. SCE). Il couvre un espace chimique plus large que QCDGE (hétérocycles, fonctions carbonylées, amines aromatiques), groupements présents dans les photocatalyseurs organiques réels.

L’absence de coordonnées DFT pour ces molécules implique que toutes les géométries 3D sont générées par ETKDG (RDKit). L’information per-atome MACE est donc perdue : le modèle travaille avec l’embedding MACE moyenné sur tous les atomes, répété sur les nœuds du graphe RDKit. Malgré cette limitation, les résultats montrent que la branche SchNet et les descripteurs RDKit compensent efficacement.

III. B. Architecture et résultats

Le squelette de cette branche est identique à la branche précédente (même backbone MACE–GINEConv–SchNet–TopoGNN–RDKit), avec une tête à deux sorties [$E_{HOMO}$, $E_{ox}$]. Le backbone est initialisé depuis les poids de de la branche précédente, avec seulement une réinitialisation de la tête de préditcion.

Métrique Standalone final (201 époques)
Score combiné 0.1078
$E_{HOMO}$ MAE 80.8 meV
$E_{ox}$ MAE 53.4 mV
Objectif $E_{ox}$ < 70 mV

53,4 mV représente une précision remarquable : la relation DFT→expérimental dans OxPot présente elle-même 64 mV de RMSE résiduel. Le modèle compense partiellement les erreurs systématiques du solvant implicite (COSMO), ce qui place son erreur en dessous de la limite théorique de la méthode de référence.


IV. Fusion des branches : structure en Y

Avec deux backbones indépendants, l’un spécialisé sur les propriétés à l’état excité (QCDGE) et l’autre sur les potentiels redox (OxPot), la question centrale est la méthode combinatoire pour prédire simultanément les six propriétés qui gouvernent un photocatalyseur.

La solution retenue a été architecture en Y à fusion tardive : les deux backbones sont gelés, et seule une tête MLP légère est entraînée sur la concaténation de leurs embeddings.

                  QCDGE (443k mol., TD-DFT)
                         ↓
                   GAP_BASE (gelé)
                   ──────────────  →  emb_A [B, 1344]    ──┐
                                                           ├──► concat [B, 2592]  ──► MLP  ──► [6 cibles]
                   OxPot (15k mol., exp.)                  │
                         ↓                                 │
                   REDOX_BASE (gelé)                       │
                   ────────────────  →  emb_B [B, 1248]  ──┘

Les six cibles sont ainsi prédites simultanément : $E_{gap}$, $E_{S_1}$, $E_{T_1}$, $f_{S_1}$, $E_{HOMO}$, $E_{ox}$.

IV. A. Défi technique : décalage d’échelle des espaces d’embedding

Les deux espaces d’embedding présentent des distributions très différentes : les écarts-types de emb_A varient dans [0.009, 7.830], ceux de emb_B dans [0.029, 10.017], soit un facteur 10 entre les deux espaces. Sans correction, le réseau apprendrait essentiellement à travers emb_B.

La solution est une normalisation z-score séparée de chaque embedding (statistiques calculées sur le train uniquement), avant concaténation :

emb_A_norm = (emb_A − μ_A) / σ_A       emb_B_norm = (emb_B − μ_B) / σ_B
fused = concat([emb_A_norm, emb_B_norm])   →   [B, 2592]

IV. B. Loss masquée sur les cibles mixtes

QCDGE ne contient pas de potentiel d’oxydation, et OxPot ne contient pas de propriétés à l’état excité. La loss Huber est masquée par cible : les NaN sont exclus du gradient molécule par molécule. Le déséquilibre de dataset (4 cibles QCDGE vs. 2 cibles OxPot) est compensé par un sampler pondéré (WeightedRandomSampler).

IV. C. Architecture de la tête MLP (résultat HPO sur 44 + 33 trials)

Linear(2592 → 1024) → LayerNorm → ELU → Dropout(0.028)
Linear(1024 → 256)  → ELU
Linear(256  → 128)  → ELU
Linear(128  → 6)    → [gap, E_S₁, E_T₁, f_S₁, E_HOMO, Eox]

IV. D. Résultats

Cible Backbone seul Structure en Y Dégradation
$E_{gap}$ 47.9 meV 49.2 meV +2.6 % ✓
$E_{S_1}$ 41.8 meV 43.1 meV +3.1 % ✓
$E_{T_1}$ 32.1 meV 32.8 meV +2.3 % ✓
$E_{HOMO}$ 80.8 meV 81.7 meV +1.1 % ✓
$E_{ox}$ 53.4 mV 54.0 mV +1.1 % ✓

La dégradation est inférieure à 5 % sur toutes les cibles physiques critiques. La tête MLP (~1,5 M paramètres) capture les corrélations inter-cibles, $E_{S_1}$ > $E_{T_1}$ est respecté pour 100 % des prédictions, sans réentraîner les encodeurs.


V. Généralisation aux photocatalyseurs : limitations et domain shift

Les performances sur QCDGE et OxPot masquent un problème fondamental : ces datasets ne contiennent pas de photocatalyseurs réels. QCDGE couvre des molécules à 10 atomes lourds au maximum ; OxPot est dominé par des structures aliphatiques simples. Un photocatalyseur organique typique (acridinium, éosine Y, thioxanthone, rhodamine …) compte 20 à 50 atomes lourds avec des systèmes π étendus fusionnés.

27 photocatalyseurs issus des travaux de Nicewicz et Tlili (acridiniums, pyryliums, xanthones, rhodamines, phénothiazines…) ont été testés avec leurs valeurs expérimentales de référence.

Propriété MAE obtenu Critère Molécules dans le critère
$E_{S_1}$ 868 meV < 300 meV 3 / 20
$E_{T_1}$ 708 meV < 300 meV 6 / 22
$E_{ox}$ 547 mV < 300 mV 4 / 10

Ce résultat bloque la suite du pipeline : les calculs DFT sur les photocatalyseurs cibles ne seront pas lancés avant que le problème de domain shift soit résolu. Deux axes sont en cours d’exploration : le remplacement de l’encodeur de base par MACE_polar-1 (pré-entraîné sur des molécules polaires de taille variable), et un curriculum de pré-entraînement sur 17 000 chromophores organiques expérimentaux (λ_abs, λ_em, ϕ_f) pour adapter le backbone aux grands systèmes aromatiques avant de retourner aux données TD-DFT.


VI. Limites structurelles et perspectives

Ces résultats, autant pour PhotoGNN que pour les modèles d’état de l’art sur QCDGE, sont à mettre en perspective avec la taille réelle des molécules photoactives étudiées en photochimie :

Molécule Formule Atomes lourds
QCDGE (maximum) ≤ 10
Azobenzène (photocommutateur) C₁₂H₁₀N₂ 14
Benzophénone (photosensibilisateur) C₁₃H₁₀O 14
Éosine Y (photocatalyseur) C₂₀H₆Br₄O₅ 30
Ir(ppy)₃ (OLED) C₃₃H₂₄IrN₃ 37

Une imbrication d’imprécisions. L’erreur du modèle (~40–50 meV par rapport aux labels DFT) s’ajoute à l’erreur intrinsèque de TD-ωB97X-D elle-même (200–400 meV par rapport aux valeurs expérimentales selon la molécule). La précision réelle est plafonnée par la méthode qui a généré les données : améliorer l’architecture au-delà d’un certain seuil n’a de sens que si les labels sont eux-mêmes plus précis (CASPT2, données expérimentales directes), ou si la distribution de test est couverte par la distribution d’entraînement.

L’architecture en Y, le curriculum de pré-entraînement, et le fine-tuning ciblé sur les photocatalyseurs réels constituent les trois leviers actifs pour repousser cette limite. La suite du pipeline, un LSTM pour la dynamique temporelle et un PINN pour les contraintes physiques (Beer-Lambert, Rehm-Weller, Stern-Volmer), ne pourra démarrer qu’une fois le domain shift résolu avec pour objectif final la conception d’un jumeau numérique de réacteur photochimique en flux continu.


VII. Stack technique

Outil Rôle
PyTorch Entraînement, différentiation automatique
PyTorch Geometric Couches GNN (GINEConv, GATConv, SAGEConv, GCNConv), pooling, radius_graph
MACE Encodeur atomistique pré-entraîné — embeddings per-atome gelés (3072D)
RDKit Parsing SMILES, construction des graphes moléculaires, 29 descripteurs globaux
Optuna Optimisation bayésienne des hyperparamètres (TPESampler, MedianPruner)
NumPy memmap Accès aux embeddings per-atome QCDGE (~101 Go) sans saturer la RAM
Matériel 2× RTX 5070 Ti (32 Go VRAM totale)

L’intégralité de l’optimisation et de l’entraînement a été réalisée sur infrastructure personnelle.


VIII. Références


← Retour au portfolio