Enabling Ab-Initio Molecular Dynamics under Bias: The CP2K+SMEAGOL Interface for Integrating Density Functional Theory and Non-Equilibrium Green Functions

Résumé Technique : Permettre la Dynamique Moléculaire Ab-Initio sous Biais via l'Interface CP2K+SMEAGOL

Énoncé du Problème
Bien que la Théorie de la Fonctionnelle de la Densité (DFT) combinée aux Fonctions de Green hors Équilibre (NEGF) soit une approche standard pour modéliser le transport quantique sous potentiels de biais externes, les méthodologies conventionnelles sont limitées aux systèmes statiques. Les implémentations existantes opèrent généralement dans un cadre canonique, calculant les propriétés de transport à des positions atomiques fixes. Cela crée un fossé significatif entre la compréhension théorique et les résultats macroscopiques des expériences, particulièrement pour les systèmes électrochimiques (EC) où les transformations sont intrinsèquement de nature grand canonique. Pour décrire avec précision l'environnement naturel des transformations EC, les simulations nécessitent une dynamique moléculaire (MD) sous biais et une description explicite des frontières ouvertes pour les électrons pouvant entrer et sortir de la cellule de calcul. De plus, les implémentations DFT-NEGF existantes manquent souvent de la capacité de calculer les forces induites par le courant nécessaires à l'optimisation de la géométrie et à la dynamique, ou sont restreintes à des codes DFT spécifiques (par exemple, SIESTA) qui ne possèdent pas les fonctionnalités présentes dans les logiciels de pointe optimisés pour les simulations en phase condensée.

Méthodologie
Les auteurs présentent une nouvelle interface entre le célèbre package de structure électronique mixte Gaussien/onde plane CP2K et le code NEGF SMEAGOL. Cette interface, nommée CP2K+SMEAGOL, est conçue comme une bibliothèque autonome capable d'être liée à n'importe quel logiciel DFT, en réutilisant autant que possible le code original de SMEAGOL.

Les principales implémentations techniques incluent :

• Forces Induites par le Courant : L'interface implémente le calcul des forces non conservatives agissant sur les ions sous biais. Ces forces sont dérivées de la dérivée temporelle des impulsions atomiques, utilisant les gradients du Hamiltonien complet par rapport aux positions atomiques. Crucialement, l'implémentation inclut le terme de force provenant de la dérivée de la matrice de densité ($\Omega_{ij}$), qui rend compte de la composante du « vent électronique » (transfert de quantité de mouvement des électrons vers les ions). Ce terme est souvent négligé dans d'autres implémentations mais s'avère essentiel pour la justesse qualitative.
• Construction de la Matrice de Densité : Pour occuper correctement les états liés au sein de la fenêtre de biais, l'interface met à jour le calcul de la matrice de densité. Au lieu d'utiliser une évaluation de contour unique ou une simple moyenne des contributions gauche/droite, elle emploie un schéma à double contour pondéré (suivant Brandbyge et al.) qui somme les intégrales des deux électrodes avec des poids déterminés par les contributions hors équilibre.
• Configuration de la Simulation : La méthode permet la dynamique moléculaire de Born-Oppenheimer sous biais. Le système est modélisé avec une région centrale de « molécule étendue » attachée à des électrodes semi-infinies. Les auto-énergies des électrodes sont traitées comme statiques et indépendantes de la densité de charge de la région centrale, une hypothèse valide lorsque la région centrale est suffisamment grande pour écranter les changements de charge.

Résultats Clés et Validation
Les auteurs ont validé l'interface CP2K+SMEAGOL par rapport à des calculs de référence effectués avec SIESTA+SMEAGOL sur plusieurs systèmes :

1. Forces à Biais Nul (Fil d'Or Infini) : L'interface reproduit avec succès les forces atomiques à biais nul, correspondant aux calculs CP2K seuls. Elle démontre également que négliger le terme de densité d'énergie ($\Omega$) conduit à des forces qualitativement incorrectes.
2. Condensateur à Plaques Parallèles : L'interface reproduit la chute de potentiel attendue à travers un condensateur, montrant un potentiel constant dans les électrodes et une chute linéaire dans le vide. Elle confirme que le schéma à double contour pondéré est nécessaire pour produire des densités de charge symétriques et des potentiels de Hartree corrects à haut biais, alors que les évaluations à contour unique produisent des résultats asymétriques et incorrects.
3. Optimisation de la Géométrie de Jonction Au-H$_2$-Au : Sous un biais appliqué (jusqu'à 1,5 V), l'interface prédit correctement l'élongation de la liaison H-H et les déplacements atomiques, ce qui est cohérent avec les résultats précédents de SIESTA+SMEAGOL. Ces changements sont attribués au courant électrique plutôt qu'au champ électrique.
4. Dynamique Moléculaire d'un Fil d'Or Solvaté : Les auteurs ont réalisé les premières simulations de dynamique moléculaire DFT-NEGF sur un système de grande échelle en phase condensée (un fil monoatomique d'or solvaté par 166 molécules d'eau).
   • Stabilité : La dynamique est stable avec une dérive d'énergie minimale à long terme, malgré la nature de système ouvert.
   • Migration des Électrons : Dans une trajectoire de 200 fs à 1 V de biais, la simulation a observé des effets de migration électronique, où le minimum de densité électronique locale s'est déplacé de l'atome central du fil vers l'atome le plus à droite, conformément à l'écrantage de la charge positive de l'électrode de droite.
   • Transmission : Les spectres de transmission du fil solvaté ont montré une similitude qualitative avec un fil idéal, les différences étant attribuées à des canaux de tunnel supplémentaires ouverts par les molécules d'eau.

Performance
Le principal goulot d'étranglement computationnel est l'évaluation de la densité électronique comme une intégrale de la fonction de Green, nécessitant 96 calculs de points d'énergie (64 sur l'axe réel, 32 sur l'axe complexe) par étape SCF, remplaçant une seule diagonalisation de Hamiltonien.

• Parallélisation : Le code utilise MPI pour les calculs dans l'espace réel et les points k, et OpenMP pour les boucles DO et les routines LAPACK. Les performances passent bien jusqu'à 32 threads OpenMP.
• Coût : Pour un système de fil d'or solvaté, le calcul CP2K+SMEAGOL était environ 3,2 fois plus lent par étape SCF qu'un calcul CP2K standard, malgré l'utilisation de 8 fois plus de cœurs (2048 contre 256). Les auteurs notent que cela représente une augmentation du coût d'un ordre de grandeur, rendant les longues trajectoires de MD difficiles sans accélération supplémentaire (par exemple, via ScaLAPACK ou des champs de force d'apprentissage automatique).

Signification et Revendications
L'article affirme que CP2K+SMEAGOL est l'implémentation DFT-NEGF la plus riche en fonctionnalités et la plus complète disponible pour CP2K à ce jour, surpassant l'implémentation native de CP2K qui est limitée au transport au point $\Gamma$ sans forces.

La principale signification de ce travail est la démonstration que la DFT-NEGF peut être utilisée pour effectuer des simulations de dynamique moléculaire sous biais pour des systèmes de grande échelle en phase condensée dans des conditions opérationnelles réalistes. À la connaissance des auteurs, il s'agit des premiers calculs de ce type réalisés. L'interface permet l'étude des réactions électrocatalytiques et des systèmes électrochimiques où l'interaction entre le potentiel appliqué, le flux de courant et la dynamique atomique est critique, comblant ainsi le fossé entre les modèles théoriques microscopiques et les conditions expérimentales macroscopiques.