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

Resumen Técnico: Habilitación de la Dinámica Molecular Ab-Initio bajo Sesgo mediante la Interfaz CP2K+SMEAGOL

Planteamiento del Problema
Si bien la Teoría de la Funcional de la Densidad (DFT) combinada con las Funciones de Green de No Equilibrio (NEGF) es un enfoque estándar para modelar el transporte cuántico bajo potenciales de sesgo externos, las metodologías convencionales están limitadas a sistemas estáticos. Las implementaciones existentes operan típicamente dentro de un marco canónico, calculando propiedades de transporte en posiciones atómicas fijas. Esto crea una brecha significativa entre la comprensión teórica y los resultados macroscópicos de los experimentos, particularmente para sistemas electroquímicos (EC) donde las transformaciones son intrínsecamente gran canónicas. Para describir con precisión el entorno natural de las transformaciones EC, las simulaciones requieren dinámica molecular (MD) bajo sesgo y una descripción de frontera abierta explícita de los electrones que puedan entrar y salir de la celda computacional. Además, las implementaciones de DFT-NEGF existentes a menudo carecen de la capacidad de calcular las fuerzas inducidas por la corriente necesarias para la optimización de la geometría o la dinámica, o están restringidas a códigos DFT específicos (por ejemplo, SIESTA) que carecen de características presentes en paquetes de última generación optimizados para simulaciones de fase condensada.

Metodología
Los autores presentan una nueva interfaz entre el popular paquete de estructura electrónica de mezcla Gaussiano/onda plana CP2K y el código de NEGF SMEAGOL. Esta interfaz, denominada CP2K+SMEAGOL, está diseñada para ser una biblioteca independiente vinculable con cualquier software de DFT, reutilizando el código original de SMEAGOL tanto como sea posible.

Las implementaciones técnicas clave incluyen:

• Fuerzas Inducidas por la Corriente: La interfaz implementa el cálculo de fuerzas no conservativas que actúan sobre los iones bajo sesgo. Estas fuerzas se derivan de la derivada temporal de los momentos atómicos, utilizando los gradientes del Hamiltoniano completo con respecto a las posiciones atómicas. Crucialmente, la implementación incluye el término de fuerza que surge de la derivada de la matriz de densidad ($\Omega_{ij}$), que da cuenta del componente de "viento electrónico" (transferencia de momento de electrones a iones). Este término suele omitirse en otras implementaciones, pero se demuestra que es esencial para la corrección cualitativa.
• Construcción de la Matriz de Densidad: Para ocupar correctamente los estados ligados dentro de la ventana de sesgo, la interfaz actualiza el cálculo de la matriz de densidad. En lugar de utilizar una evaluación de contorno único o un promedio simple de las contribuciones izquierda/derecha, emplea un esquema de doble contorno ponderado (siguiendo a Brandbyge et al.) que suma integrales de ambos electrodos con pesos determinados por las contribuciones de no equilibrio.
• Configuración de la Simulación: El método permite la dinámica molecular de Born-Oppenheimer bajo sesgo. El sistema se modela con una región de "molécula extendida" central conectada a electrodos semi-infinitos. Las auto-energías de los electrodos se tratan como estáticas e independientes de la densidad de carga de la región central, una suposición válida cuando la región central es lo suficientemente grande como para blindar los cambios de carga.

Resultados Clave y Validación
Los autores validaron la interfaz CP2K+SMEAGOL contra cálculos de referencia realizados con SIESTA+SMEAGOL a través de varios sistemas:

1. Fuerzas de Sesgo Cero (Alambre de Au Infinito): La interfaz reproduce con éxito las fuerzas atómicas a sesgo cero, coincidiendo con los cálculos de solo CP2K. También demuestra que omitir el término de la matriz de densidad de energía ($\Omega$) conduce a fuerzas cualitativamente incorrectas.
2. Capacitor de Placas Paralelas: La interfaz reproduce la caída de potencial esperada a través de un capacitor, mostrando un potencial constante en los electrodos y una caída lineal en el vacío. Confirma que el esquema de doble contorno ponderado es necesario para producir densidades de carga simétricas y potenciales de Hartree correctos a alto sesgo, mientras que las evaluaciones de contorno único producen resultados asimétricos e incorrectos.
3. Optimización de Geometría de la Unión Au-H$_2$-Au: Bajo un sesgo aplicado (hasta 1.5 V), la interfaz predice correctamente la elongación del enlace H-H y los desplazamientos atómicos, siendo consistentes con los resultados previos de SIESTA+SMEAGOL. Los cambios se atribuyen a la corriente eléctrica en lugar que al campo eléctrico.
4. Dinámica Molecular de un Alambre de Au Solvatado: Los autores realizaron las primeras simulaciones de dinámica molecular DFT-NEGF en un sistema de fase condensada a gran escala (un alambre monoatómico de Au solvatado por 166 moléculas de agua).
   • Estabilidad: La dinámica fue estable con una deriva de energía mínima a largo plazo, a pesar de la naturaleza de sistema abierto.
   • Migración de Electrones: En una trayectoria de 200 fs a 1 V de sesgo, la simulación observó efectos de migración de electrones, donde el mínimo de densidad electrónica local se desplazó del átomo central del alambre al átomo más a la derecha, consistente con el blindaje de la carga positiva del electrodo derecho.
   • Transmisión: Los espectros de transmisión para el alambre solvatado mostraron una similitud cualitativa con un alambre infinito ideal, con diferencias atribuidas a canales de tunelamiento adicionales abiertos por las moléculas de agua.

Rendimiento
El principal cuello de botella computacional es la evaluación de la densidad electrónica como una integral de la función de Green, lo que requiere 96 computaciones de puntos de energía (64 en el eje real, 32 en el eje complejo) por cada paso de SCF, reemplazando una única diagonalización del Hamiltoniano.

• Paralelización: El código utiliza MPI para computaciones en el espacio real y puntos k, y OpenMP para bucles DO y rutinas LAPACK. El rendimiento escala bien hasta 32 hilos de OpenMP.
• Costo: Para un sistema de alambre de Au solvatado, el cálculo de CP2K+SMEAGOL fue aproximadamente 3.2 veces más lento por paso de SCF que un cálculo estándar de CP2K, a pesar de usar 8 veces más núcleos (2048 vs. 256). Los autores señalan que esto representa un aumento de un orden de magnitud en el costo, lo que hace que las trayectorias de MD largas sean desafiantes sin mayor aceleración (por ejemplo, mediante ScaLAPACK o campos de fuerza de aprendizaje automático).

Significancia y Reivindicaciones
El artículo afirma que CP2K+SMEAGOL es la implementación de DFT-NEGF más completa y rica en funciones disponible para CP2K hasta la fecha, superando la implementación nativa de CP2K que está limitada al transporte en el punto $\Gamma$ sin fuerzas.

La principal significancia de este trabajo es la demostración de que la DFT-NEGF puede utilizarse para realizar simulaciones de dinámica molecular bajo sesgo para sistemas de fase condensada a gran escala bajo condiciones operativas realistas. Según conocimiento de los autores, estas representan las primeras tales cálculos realizados. La interfaz permite el estudio de reacciones electrocatalíticas y sistemas electroquímicos donde la interacción entre el potencial aplicado, el flujo de corriente y la dinámica atómica es crítica, cerrando la brecha entre los modelos teóricos microscópicos y las condiciones experimentales macroscópicas.