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

Sintesi Tecnica: Abilitazione della Dinamica Molecolare Ab-Initio sotto Bias tramite l'interfaccia CP2K+SMEAGOL

Definizione del Problema
Sebbene la Teoria del Funzionale della Densità (DFT) combinata con le Funzioni di Green Fuori dall'Equilibrio (NEGF) sia un approccio standard per modellare il trasporto quantistico sotto potenziali di bias esterni, le metodologie convenzionali sono limitate ai sistemi statici. Le implementazioni esistenti operano tipicamente all'interno di un quadro canonico, calcolando le proprietà di trasporto a posizioni atomiche fisse. Ciò crea un significativo divario tra la comprensione teorica e i risultati macroscopici degli esperimenti, in particolare per i sistemi elettrochimici (EC) dove le trasformazioni sono intrinsecamente grand-canoniche. Per descrivere accuratamente l'ambiente naturale delle trasformazioni EC, le simulazioni richiedono la dinamica molecolare (MD) sotto bias e una descrizione a confini espliciti degli elettroni che possano entrare e uscire dalla cella computazionale. Inoltre, le attuali implementazioni DFT-NEGF spesso mancano della capacità di calcolare le forze indotte dalla corrente necessarie per l'ottimizzazione della geometria e la dinamica, oppure sono limitate a specifici codici DFT (ad esempio, SIESTA) che mancano di funzionalità presenti in pacchetti all'avanguardia ottimizzati per le simulazioni in fase condensata.

Metodologia
Gli autori presentano una nuova interfaccia tra il popolare pacchetto di struttura elettronica a base mista Gaussiana/piano d'onda CP2K e il codice NEGF SMEAGOL. Questa interfaccia, denominata CP2K+SMEAGOL, è progettata per essere una libreria standalone collegabile con qualsiasi software DFT, riutilizzando il più possibile il codice originale di SMEAGOL.

Le principali implementazioni tecniche includono:

• Forze Indotte dalla Corrente: L'interfaccia implementa il calcolo delle forze non conservative che agiscono sugli ioni sotto bias. Queste forze sono derivate dalla derivata temporale dei momenti atomici, utilizzando i gradienti dell'Hamiltoniana completa rispetto alle posizioni atomiche. Fondamentalmente, l'implementazione include il termine di forza derivante dalla derivata della matrice di densità ($\Omega_{ij}$), che tiene conto della componente del "vento elettronico" (trasferimento di momento dagli elettroni agli ioni). Questo termine è spesso trascurato in altre implementazioni ma si dimostra essenziale per la correttezza qualitativa.
• Costruzione della Matrice di Densità: Per occupare correttamente gli stati legati all'interno della finestra di bias, l'interfaccia aggiorna il calcolo della matrice di densità. Invece di utilizzare una singola valutazione del contorno o una semplice media dei contributi sinistra/destra, impiega uno schema a doppio contorno pesato (seguendo Brandbyge et al.) che somma gli integrali da entrambi i lead con pesi determinati dai contributi fuori equilibrio.
• Configurazione della Simulazione: Il metodo consente la dinamica molecolare Born-Oppenheimer sotto bias. Il sistema è modellato con una regione centrale di "molecola estesa" collegata a elettrodi semi-infiniti. Le auto-energie dei lead sono trattate come statiche e indipendenti dalla densità di carica della regione centrale, un'ipotesi valida quando la regione centrale è sufficientemente grande da schermare le variazioni di carica.

Risultati Chiave e Validazione
Gli autori hanno validato l'interfaccia CP2K+SMEAGOL rispetto ai calcoli di riferimento eseguiti con SIESTA+SMEAGOL attraverso diversi sistemi:

1. Forze a Zero-Bias (Filo d'Oro Infinito): L'interfaccia riproduce con successo le forze atomiche a zero bias, corrispondendo ai calcoli solo con CP2K. Dimostra inoltre che trascurare il termine della densità di energia ($\Omega$) porta a forze qualitativamente errate.
2. Condensatore a Piastre Parallele: L'interfaccia riproduce la caduta di potenziale attesa attraverso un condensatore, mostrando un potenziale costante nei lead e una caduta lineare nel vuoto. Conferma che lo schema a doppio contorno pesato è necessario per produrre densità di carica simmetriche e potenziali di Hartree corretti ad alto bias, mentre le valutazioni a singolo contorno producono risultati asimmetrici e errati.
3. Ottimizzazione della Geometria Giunzione Au-H$_2$-Au: Sotto bias applicato (fino a 1,5 V), l'interfaccia predice correttamente l'allungamento del legame H-H e gli spostamenti atomici, in linea con i precedenti risultati di SIESTA+SMEAGOL. I cambiamenti sono attribuiti alla corrente elettrica piuttosto che al campo elettrico.
4. Dinamica Molecolare di un Filo d'Oro Solvatato: Gli autori hanno eseguito le prime simulazioni di dinamica molecolare DFT-NEGF su un sistema in fase condensata su larga scala (un filo monoatomico d'oro solvatato da 166 molecole d'acqua).
   • Stabilità: La dinamica è stata stabile con una deriva energetica a lungo termine minima, nonostante la natura a sistema aperto.
   • Migrazione Elettronica: In una traiettoria di 200 fs a 1 V di bias, la simulazione ha osservato effetti di migrazione elettronica, dove il minimo della densità elettronica locale si è spostato dall'atomo centrale all'atomo più a destra, coerentemente con la schermatura della carica positiva dell'elettrodo destro.
   • Trasmissione: Gli spettri di trasmissione per il filo solvatato hanno mostrato una somiglianza qualitativa con un filo ideale infinito, con differenze attribuite a ulteriori canali di tunneling aperti dalle molecole d'acqua.

Prestazioni
Il principale collo di bottiglia computazionale è la valutazione della densità elettronica come integrale della funzione di Green, che richiede 96 calcoli di punti energetici (64 sull'asse reale, 32 sull'asse complesso) per ogni passo SCF, sostituendo una singola diagonalizzazione dell'Hamiltoniana.

• Parallelizzazione: Il codice utilizza MPI per i calcoli nello spazio reale e per i punti k, e OpenMP per i cicli DO e le routine LAPACK. Le prestazioni scalano bene fino a 32 thread OpenMP.
• Costo: Per un sistema di filo d'oro solvatato, il calcolo CP2K+SMEAGOL è stato circa 3,2 volte più lento per passo SCF rispetto a un calcolo standard CP2K, nonostante l'uso di 8 volte più core (2048 contro 256). Gli autori notano che questo rappresenta un aumento del costo di un ordine di grandezza, rendendo le traiettorie MD a lungo termine impegnative senza ulteriore accelerazione (ad esempio, tramite ScaLAPACK o campi di forza basati sul machine learning).

Significato e Rivendicazioni
L'articolo sostiene che CP2K+SMEAGOL sia l'implementazione DFT-NEGF più ricca di funzionalità e completa disponibile per CP2K ad oggi, superando l'implementazione nativa di CP2K che è limitata al trasporto al punto $\Gamma$ senza forze.

Il significato primario di questo lavoro è la dimostrazione che la DFT-NEGF può essere utilizzata per eseguire simulazioni di dinamica molecolare sotto bias per sistemi in fase condensata su larga scala in condizioni operative realistiche. Per quanto ne sappiano gli autori, si tratta dei primi calcoli di questo tipo eseguiti. L'interfaccia consente lo studio di reazioni elettro-catalitiche e sistemi elettrochimici dove l'interazione tra potenziale applicato, flusso di corrente e dinamica atomica è critica, colmando il divario tra i modelli teorici microscopici e le condizioni sperimentali macroscopiche.