O. Coulaud Projet Numath IECN/ Inria-Lorraine

1. Couplage de modèles pour la simulation moléculaire complexe SIMBIO O. Coulaud Projet Numath IECN/ Inria-Lorraine

2. Objectifs de ARC SIMBIO • Simuler des phénomènes moléculaires de grande taille pour • suivi de réaction chimique • étude de propriétés sur de gros systèmes • transfert d’un proton • mécanisme d’ouverture de fermeture d’une membrane • Valider les algorithmes en terme : • convergence des méthodes utilisées • erreurs d’approximation : choix des paramètres • Développer de nouveaux algorithmes • Partenaires : Apache (Rhône-Alpes), LCTN, CEA (Dpt. Sciences du vivant, Grenoble) P.E. Bernard , S. Crousy, T. Gautier, R. Krenenou, B. Maigret, B. Pinçon, B. Plateau

3. Motivation • Interactions intra et inter molécules sont décrites par l’équation de schrodinger équation non linéaire dans R3N , N est le nombre d’électrons, K le nombre de noyaux, j est une fonction anti-symétrique normée, E énergie. Limitation à une 100 d ’atomes (H, C, N, O, S, …) • mécanique classique : champ de force empirique différents champs de forces (AMBER, CHARMm, Xplor, …) on traite jusqu’à 100 000 atomes, impossibilité pour suivre les réactions chimiques, ... • On introduit un champ de force hybride

4. Zone 3 Continuum Mécanique quantique Zone 1 QM ~100 atomes MM Mécanique classique Zone 2 Molécule+solvant ~10 à 50 000 atomes Système moléculaire hybride Couplage de modèles : quantique, mécanique moléculaire, continuum

5. Le champ de force hybride est F = FQ + FM + FC + FMC + FMQ où • la contribution du continuum est FC (xat) = -qat Esolvant(xat), • FQ est issu d’un problème de valeurs propres (Hartree-Fock). • FMQ terme de couplage mécanique classique - quantique. • FMC terme de couplage mécanique classique - continuum. On utilise ce champ de force en • statique : optimisation de géométrie (minimisation), • dynamique : suivi de réaction, … (dynamique newtonienne). On se concentre sur le couplage mécanique classique-continuum pour la dynamique moléculaire.

6. Energie hybride • Le potentiel mécanique classique avec et

7. Energie de réaction du solvant où est calculé par une méthode du continuum. • Energie de couplage • un terme de répulsion-attraction de Leonard-Jones • contrainte de la simulation (volume constant, …) • …

8. solvant S e1 x x x (ak,qk) x x e0 x x molécules Potentiel issu du solvant Surface exclue au solvant Modèle du Continuum • représentation explicite de la molécule + quelques couches du solvant • Le solvant est un milieu continu, et a une constante diélectrique e1. • Le solvant n’est pas ionique. Equations

9. Formulation intégrale On décompose u en avec ou sest la solution de On utilise une méthode de collocation • approximation P2 de la surface et une approximation P1 de la solution, • quadrature de degré 4-3-2 pour évaluer l'intégrale. Matrice non symétrique - méthode directe pour résoudre le système linéaire (LU), Modèle à mémoire partagée pour le parallélisme (OpenMP).

10. Dynamique moléculaire Equations de Newton Schéma de Leapfrog le schéma conserve l’énergie (simulation NVE)

11. Takakaw Code développé dans le projet Apache • Fonctionnalités • simulation NVE, • rayon de coupure sur les forces non liées (VDW, électrostatique), • contraintes : harmoniques (position), shake sur les atomes d’hydrogène (distance), • système borné ou conditions périodiques, • dynamique de Newton ou de Langevin. • Parallélisme • décomposition géométrique du domaine - découpage en boites de taille ~Rc, • Athapascan : environnement de programmation (MPI + threads Posix), • modèle maître esclaves ; mais multi-threads (recouvrement calculs-communications), • équilibrage de charge statique et dynamique.

12. Dt = N Dt Algorithme de couplage Différentes échelles de temps • temps caractéristique des vibrations << temps lié à la force de Coulomb On utilise un schéma à pas multiples la force est définie Dt

13. Les informations provenant de la CM • Énergie provenant du solvant • Force provenant du solvant Les informations provenant de la DM

14. Échanges d’informations • Méthode du continuum • Position et charges des atomes pour le second membre et la surface S • Dynamique moléculaire • valeur du gradient du potentiel sur les atomes, • énergie de réaction du solvant. coût : échange d’un vecteur de taille 3N soit pour 100 000 atomes 2x2,4 Mo à échanger par macro pas.

15. Une alternative • Ne pas reconstruire la surface S à chaque macro pas ensemble moléculaire est contenu dans un ellipsoide. • Utiliser les structures internes des applications • pavage de l’espace en cube pour la dynamique • maillage de la surface pour le continuum

16. 4 x le nombre de sommets de la grille • On duplique ces deux structures maillage et grille pour diminuer les échanges • évaluation du second membre sur le maillage par la DM, • évaluation du potentiel et de son gradient sur la grille par la CM. • Le travail supplémentaire • interpolation sur les atomes à partir des valeurs sur la grille. • Coût DM CM Indépendant du nombre d’atomes mais va dépendre de la géométrie et du Cut-off. 3 x le nombre de nœuds du maillage

17. Conclusion & perspectives • Prototype qui marche mais il reste beaucoup à faire • valider le modèle sur un cas réaliste, • liens entre le pas de la grille (Rcoupure) et le pas du maillage (h), • méthode à pas multiples choix entre Dt et Dt, • ... • Mettre en place un algorithme qui conserve • l ’énergie totale du système (en cours), • le centre de masse. • Introduire le couplage quantique