Análisis de fundiciones y fluidos de ab initio Molecular Dynamics Simulations con el paquete UMD

Resumen

Los fundidos y fluidos son vectores ubicuos de transporte de masa en sistemas naturales. Hemos desarrollado un paquete de código abierto para analizar simulaciones de dinámica molecular ab initio de dichos sistemas. Calculamos las propiedades estructurales (unión, clusterización, especiación química), transporte (difusión, viscosidad) y termodinámicas (espectro vibratorio).

Resumen

Hemos desarrollado un paquete de código abierto basado en Python para analizar los resultados derivados de simulaciones de dinámica molecular ab initio de fluidos. El paquete es el más adecuado para aplicaciones en sistemas naturales, como fundidos de silicato y óxido, fluidos a base de agua y varios fluidos supercríticos. El paquete es una colección de scripts de Python que incluyen dos bibliotecas principales que se ocupan de los formatos de archivo y la cristalografía. Todos los scripts se ejecutan en la línea de comandos. Proponemos un formato simplificado para almacenar las trayectorias atómicas y la información termodinámica relevante de las simulaciones, que se guarda en archivos UMD, que significa Dinámica Molecular Universal. El paquete UMD permite el cálculo de una serie de propiedades estructurales, de transporte y termodinámicas. Comenzando con la función de distribución de pares, define las longitudes de enlace, construye una matriz de conectividad interatómica y, finalmente, determina la especiación química. Determinar la vida útil de las especies químicas permite realizar un análisis estadístico completo. Luego, los scripts dedicados calculan los desplazamientos cuadráticos medios para los átomos, así como para las especies químicas. El análisis de autocorrelación implementado de las velocidades atómicas arroja los coeficientes de difusión y el espectro vibratorio. El mismo análisis aplicado sobre las tensiones produce la viscosidad. El paquete está disponible a través del sitio web de GitHub y a través de su propia página dedicada del proyecto ERC IMPACT como paquete de acceso abierto.

Introducción

Los fluidos y los fundidos son vectores activos de transporte químico y físico en entornos naturales. Las elevadas tasas de difusión atómica favorecen los intercambios y reacciones químicas, la baja viscosidad junto con la flotabilidad variable favorecen la gran transferencia de masa y las relaciones de densidad de fusión cristalina favorecen la estratificación dentro de los cuerpos planetarios. La ausencia de una red periódica, las altas temperaturas típicas requeridas para alcanzar el estado fundido y la dificultad para el enfriamiento hacen que la determinación experimental de una serie de propiedades obvias, como la densidad, la difusión y la viscosidad, sea extremadamente desafiante. Estas dificultades hacen que los métodos computacionales alternativos sean herramientas fuertes y útiles para investigar esta clase de materiales.

Con el advenimiento de la potencia de cálculo y la disponibilidad de superordenadores, actualmente se emplean dos grandes técnicas de simulación atomística numérica para estudiar el estado dinámico de un sistema atomístico no cristalino, Monte ^Carlo1 y la dinámica molecular (MD)^1,2. En las simulaciones de Monte Carlo, el espacio de configuración se muestrea aleatoriamente; Los métodos de Monte Carlo muestran una escala lineal en paralelización si todas las observaciones de muestreo son independientes entre sí. La calidad de los resultados depende de la calidad del generador de números aleatorios y de la representatividad del muestreo. Los métodos de Monte Carlo muestran escala lineal en paralelización si el muestreo es independiente entre sí. En dinámica molecular (DM) el espacio configuracional es muestreado por trayectorias atómicas dependientes del tiempo. A partir de una configuración dada, las trayectorias atómicas se calculan integrando las ecuaciones newtonianas de movimiento. Las fuerzas interatómicas se pueden calcular utilizando potenciales interatómicos modelo (en la DM clásica) o utilizando métodos de primeros principios (in ab initio, o primeros principios, MD). La calidad de los resultados depende de la longitud de la trayectoria y su capacidad para no sentirse atraído por los mínimos locales.

Las simulaciones de dinámica molecular contienen una gran cantidad de información, toda relacionada con el comportamiento dinámico del sistema. Las propiedades medias termodinámicas, como la energía interna, la temperatura y la presión, son bastante estándar para calcular. Pueden extraerse de los archivos de salida de las simulaciones y promediarse, mientras que las cantidades relacionadas directamente con el movimiento de los átomos, así como con su relación mutua, deben calcularse después de la extracción de las posiciones y velocidades atómicas.

En consecuencia, se ha dedicado mucho esfuerzo a visualizar los resultados, y varios paquetes están disponibles hoy en día, en diferentes plataformas, de código abierto o no [^Ovito3, ^VMD4, ^Vesta5, ^Travis6, etc.]. Todas estas herramientas de visualización tratan eficientemente las distancias interatómicas, y como tales, permiten el cálculo eficiente de las funciones de distribución de pares y los coeficientes de difusión. Varios grupos que realizan simulaciones de dinámica molecular a gran escala tienen software patentado para analizar varias otras propiedades resultantes de las simulaciones, a veces en shareware u otras formas de acceso limitado a la comunidad, y a veces limitado en alcance y uso a algunos paquetes específicos. Algoritmos sofisticados para extraer información sobre enlaces interatómicos, patrones geométricos y termodinámica se desarrollan e implementan en algunos de estos paquetes3,4,5,6,7, etc.

Aquí proponemos el paquete UMD, un paquete de código abierto escrito en Python para analizar la salida de simulaciones de dinámica molecular. El paquete UMD permite el cálculo de una amplia gama de propiedades estructurales, dinámicas y termodinámicas (Figura 1). El paquete está disponible a través del sitio web de GitHub (https://github.com/rcaracas/UMD_package) y a través de una página dedicada (http://moonimpact.eu/umd-package/) del proyecto ERC IMPACT como un paquete de acceso abierto.

Para que sea universal y más fácil de manejar, nuestro enfoque es extraer primero toda la información relacionada con el estado termodinámico y las trayectorias atómicas del archivo de salida de la ejecución de dinámica molecular real. Esta información se almacena en un archivo dedicado, cuyo formato es independiente del paquete MD original donde se ejecutó la simulación. Llamamos a estos archivos archivos "umd", que significa Dinámica Molecular Universal. De esta manera, nuestro paquete UMD puede ser utilizado fácilmente por cualquier grupo ab initio con cualquier software, todo con un mínimo esfuerzo de adaptación. El único requisito para utilizar el paquete actual es escribir el analizador apropiado desde la salida del software MD en particular en el formato de archivo umd, si aún no existe. Por el momento, proporcionamos dichos analizadores para los paquetes ^VASP8 y ^QBox9 .

Figura 1: Diagrama de flujo de la biblioteca UMD.
Las propiedades físicas están en azul, y los principales scripts de Python y sus opciones están en rojo. Haga clic aquí para ver una versión más grande de esta figura.

Los archivos umd son archivos ASCII; la extensión típica es "umd.dat" pero no obligatoria. Todos los componentes de análisis pueden leer archivos ASCII del formato umd, independientemente de la extensión del nombre real. Sin embargo, algunos de los scripts automáticos diseñados para realizar estadísticas rápidas a gran escala en varias simulaciones buscan específicamente archivos con la extensión umd.dat. Cada propiedad física se expresa en una línea. Cada línea comienza con una palabra clave. De esta manera, el formato es altamente adaptable y permite agregar nuevas propiedades al archivo umd, al tiempo que conserva su legibilidad en todas las versiones. Las primeras 30 líneas del archivo umd de la simulación de pirolita a 4.6 GPa y 3000 K, utilizadas a continuación en la discusión, se muestran en la Figura 2.

Figura 2: El comienzo del archivo umd que describe la simulación de pirolita líquida a 4.6 GPa y 3000 K.
El encabezado va seguido de la descripción de cada instantánea. Cada propiedad se escribe en una línea, que contiene el nombre de la propiedad física, los valores y las unidades, todas separadas por espacios. Haga clic aquí para ver una versión más grande de esta figura.

Todos los archivos umd contienen un encabezado que describe el contenido de la celda de simulación: el número de átomos, electrones y tipos atómicos, así como detalles para cada átomo, como su tipo, símbolo químico, número de electrones de valencia y su masa. Una línea vacía marca el final del encabezado y lo separa de la parte principal del archivo umd.

Luego se detalla cada paso de la simulación. Primero, se dan los parámetros termodinámicos instantáneos, cada uno en una línea diferente, especificando (i) el nombre del parámetro, como energía, tensiones, presión hidrostática equivalente, densidad, volumen, parámetros de red, etc., (ii) su(s) valor(es), y (iii) sus unidades. A continuación se presenta una tabla que describe los átomos. Una línea de cabecera da las diferentes medidas, como posiciones cartesianas, velocidades, cargas, etc., y sus unidades. Luego, cada átomo se detalla en una línea. Por grupos de tres, correspondientes a los tres ejes x, y, z , las entradas son: las posiciones reducidas, las posiciones cartesianas plegadas en la celda de simulación, las posiciones cartesianas (que tienen debidamente en cuenta el hecho de que los átomos pueden atravesar varias celdas unitarias durante una simulación), las velocidades atómicas y las fuerzas atómicas. Las dos últimas entradas son escalares: carga y momento magnético.

Dos bibliotecas principales garantizan el correcto funcionamiento de todo el paquete. La biblioteca umd_process.py se ocupa de los archivos umd, como la lectura y la impresión. La biblioteca crystallography.py se ocupa de toda la información relacionada con la estructura atómica real. La filosofía subyacente de la biblioteca crystallography.py es tratar la red como un espacio vectorial. Los parámetros de la celda unitaria junto con su orientación representan los vectores base. El "Espacio" tiene una serie de atributos escalares (volumen específico, densidad, temperatura y número específico de átomos), propiedades termodinámicas (energía interna, presión, capacidad calorífica, etc.) y una serie de propiedades tensoriales (tensión y elasticidad). Los átomos pueblan este espacio. La clase "Lattice" define este conjunto, junto con varios cálculos cortos, como volumen específico, densidad, obtención de la red recíproca del directo, etc. La clase "Átomos" define los átomos. Se caracterizan por una serie de propiedades escalares (nombre, símbolo, masa, número de electrones, etc.) y una serie de propiedades vectoriales (la posición en el espacio, ya sea relativa a la base vectorial descrita en la clase Lattice, o relativa a coordenadas cartesianas universales, velocidades, fuerzas, etc.). Aparte de estas dos clases, la biblioteca de crystallography.py contiene una serie de funciones para realizar una variedad de pruebas y cálculos, como distancias atómicas o multiplicación de celdas. La tabla periódica de los elementos también se incluye como diccionario.

Los diversos componentes del paquete umd escriben varios archivos de salida. Como regla general, todos son archivos ASCII, todas sus entradas están separadas por pestañas y se hacen lo más autoexplicativas posible. Por ejemplo, siempre indican claramente la propiedad física y sus unidades. Los archivos umd.dat cumplen plenamente con esta regla.

Protocolo

1. Análisis de las corridas de dinámica molecular

NOTA: El paquete está disponible a través del sitio web de GitHub (https://github.com/rcaracas/UMD_package) y a través de una página dedicada (http://moonimpact.eu/umd-package/) del proyecto ERC IMPACT como un paquete de acceso abierto.

1. Extraiga cada conjunto específico de propiedades físicas utilizando uno o más scripts de Python dedicados del paquete. Ejecute todos los scripts en la línea de comandos; todos emplean una serie de banderas, que son lo más consistentes posible de un guión a otro. Los indicadores, su significado y los valores predeterminados se resumen en la Tabla 1.

Bandera	Significado	Script usándolo	Valor predeterminado
-h	Ayuda breve	todo
-f	Nombre de archivo UMD	todo
-i	Pasos de termalización a descartar	todo	0
-i	Archivo de entrada que contiene los enlaces interatómicos	especiación	bonds.input
-s	Muestreo de la frecuencia	msd, especiación	1 (se considera cada paso)
-a	Lista de átomos o aniones	especiación
-c	Lista de cationes	especiación
-l	Longitud del enlace	especiación	2
-t	Temperatura	vibraciones, reología
-v	Discretización de la anchura de la ventana de muestreo de la trayectoria para el análisis de desplazamiento cuadrático medio	Msd	20
-z	Discretización del inicio de la ventana de muestreo de la trayectoria para el análisis de desplazamiento cuadrático medio	Msd	20

Tabla 1: Banderas más comunes utilizadas en el paquete UMD y su significado más común.

2. Comience transformando la salida de la simulación MD realizada en un código de primeros principios, como ^VASP8 o ^QBox9, en un archivo UMD.
   1. Si las simulaciones de MD se realizaron en VASP, entonces en la línea de comandos escriba:
      VaspParser.py -f -i
      donde el indicador –f define el nombre del archivo VASP OUTCAR y el –i la longitud de termalización.
      NOTA: El paso inicial, definido por –i permite descartar los primeros pasos de las simulaciones, que representan la termalización. En una ejecución típica de dinámica molecular, la primera parte del cálculo representa la termalización, es decir, el tiempo que tarda el sistema en describir una distribución gaussiana de la temperatura, y en que todo el sistema exhiba fluctuaciones de la temperatura, presión, energía, etc. alrededor de los valores de equilibrio. Esta parte de termalización de la simulación no debe tenerse en cuenta al analizar las propiedades estadísticas del fluido.
3. Transforme el archivo . umd archivos en . xyz para facilitar la visualización en varios otros paquetes, como ^VMD4 o ^Vesta5. En la línea de comandos, escriba:
   umd2xyz.py -f -i -s
   donde –f define el nombre del . umd file, –i define el período de termalización a descartar, y –s la frecuencia del muestreo de la trayectoria almacenada en el . archivo umd . Los valores predeterminados son –i 0 –s 1, es decir, considerando todos los pasos de la simulación, sin descartar ninguno.
4. Invierta el archivo umd en archivos POSCAR de tipo VASP utilizando el script umd2poscar.py; las instantáneas de las simulaciones se pueden seleccionar con una frecuencia predefinida. En la línea de comandos, escriba:
   umd2poscar.py -f -i -l -s
   donde –l representa el último paso que se va a transformar en archivo POSCAR. Los valores predeterminados son -i 0 -l 10000000 -s 1. Este valor de –l es lo suficientemente grande como para cubrir una trayectoria completa típica.

2. Realizar el análisis estructural

1. Ejecute el script gofrs_umd.py para calcular la función de distribución de pares (PDF) _gᴀʙ(r) para todos los pares de tipos atómicos A y B (Figura 3). La salida se escribe en un archivo ASCII, separado por tabulaciones, con la extensión gofrs.dat. En la línea de comandos, escriba:
   gofrs_umd.py -f -s < Sampling_Frequency > -d -i
   NOTA: Los valores predeterminados son Sampling_Frequency (la frecuencia para muestrear la trayectoria) = 1 paso; DiscretizationInterval (para trazar el g(r)) = 0,01 Å; InitialStep (número de pasos al principio de la trayectoria que se descartan) = 0. El PDF radial, _gᴀʙ(r) es el número medio de átomos de tipo B a una distancia d_ᴀʙ dentro de una capa esférica de radio r y espesor dr centrada en los átomos de tipo A (Figura 3):
   
   
   con ρ la densidad atómica, _NA y _NB el número de átomos de tipo A y B, y δ(r−_rᴀʙ) la función delta que es igual a 1 si los átomos A y B se encuentran a una distancia entre r y r+dr. La abscisa del primer máximo de _gᴀʙ(r) da la longitud de enlace de mayor probabilidad entre los átomos de tipo A y B, que es la más cercana a una distancia de enlace promedio que podemos determinar. El primer mínimo delimita el alcance de la primera esfera de coordinación. Por lo tanto, la integral sobre el PDF hasta el primer mínimo da el número de coordinación promedio. La suma de las transformadas de Fourier de _gᴀʙ(r) para todos los pares de tipos atómicos A y B produce el patrón de difracción del fluido, tal como se obtiene experimentalmente con un difractómetro. Sin embargo, en realidad, como a menudo las esferas de coordinación de alto orden faltan en el _gᴀʙ(r), el patrón de difracción no se puede obtener en su totalidad.

Figura 3: Determinación de la función de distribución de pares.
a) Para cada átomo de una especie (por ejemplo, rojo), todos los átomos de la especie coordinadora (por ejemplo, gris y/o rojo) se cuentan en función de la distancia. (b) El gráfico de distribución de distancia resultante para cada instantánea, que en esta etapa es solo una colección de funciones delta, se promedia sobre todos los átomos y todas las instantáneas y se pondera mediante la distribución de gas ideal para generar (c) la función de distribución de pares que es continua. El primer mínimo de g(r) es el radio de la primera esfera de coordinación, utilizado posteriormente en el análisis de especiación. Haga clic aquí para ver una versión más grande de esta figura.

2. Extraer las distancias medias de enlace interatómico como los radios de las primeras esferas de coordinación. Para ello, identifique la posición del primer máximo de las funciones _gᴀʙ(r): trace el archivo gofrs.dat en una aplicación de hoja de cálculo y busque los máximos y mínimos para cada par de átomos.
3. Identifique el radio de la primera esfera de coordinación, como el primer mínimo del PDF, _gᴀʙ(r), utilizando un software de hoja de cálculo. Esta es la base para todo el análisis estructural del fluido; el PDF produce el estado de enlace promedio de los átomos en el fluido.
4. Extraiga las distancias de los primeros mínimos, es decir, la abscisa, y escríbalas en un archivo separado, llamado, por ejemplo, bonds.input. Alternativamente, ejecute uno de los scripts analyze_gofr del paquete UMD para identificar los máximos y los mínimos de las funciones _gᴀʙ(r). En la línea de comandos, escriba:
   analyze_gofr_semi_automatic.py
5. Haga clic en la posición del máximo y el mínimo de la función _gᴀʙ(r) que se muestra en el gráfico que abre el programa. El script escanea automáticamente la carpeta actual, identifica todos los archivos gofrs.dat y realiza el análisis para cada uno de ellos. Haga clic de nuevo en el máximo y el mínimo en la ventana cada vez que el script necesite una suposición inicial informada.
6. Abra y mire el archivo generado automáticamente llamado bonds.input que contiene las distancias de enlace interatómico.

3. Realizar el análisis de especiación

1. Calcular la topología de enlace entre los átomos, utilizando el concepto de conectividad dentro de la teoría de grafos: los átomos son los nodos y los enlaces interatómicos son los caminos. El script speciation_umd.py necesita las distancias de enlace interatómico definidas en el archivo bonds.input .
   NOTA: La matriz de conectividad se construye en cada paso temporal: dos átomos que se encuentran a una distancia menor que el radio de su correspondiente primera esfera de coordinación se consideran unidos, es decir, conectados. Varias redes atómicas se construyen tratando los átomos como nodos en un gráfico cuyas conexiones están definidas por este criterio geométrico. Estas redes son las especies atómicas, y su conjunto define la especiación atómica en ese fluido en particular (Figura 4).

Figura 4: Identificación de los cúmulos atómicos.
Los poliedros de coordinación se definen utilizando distancias interatómicas. Todos los átomos a una distancia menor que un radio especificado se consideran enlazados. Aquí el umbral corresponde a la primera esfera de coordinación (los círculos rojos claros), definida en la Figura 1. La polimerización y, por lo tanto, las especies químicas se obtienen de las redes de los átomos unidos. Nótese el cúmulo central _Red1Grey2, que está aislado de los otros átomos, que forman un polímero infinito. Haga clic aquí para ver una versión más grande de esta figura.

2. Ejecute el script de especiación para obtener la matriz de conectividad y obtener los poliedros de coordinación o la polimerización. En la línea de comandos, escriba:
   speciation_umd.py -f -s -i -l -c -a -m -r
   donde el indicador -i da el archivo con las distancias de enlace interatómico, que se produjo por ejemplo en el paso anterior. Alternativamente, ejecute el script con una sola longitud para todos los enlaces definidos por el indicador -l.
   NOTA: La bandera -c especifica los átomos centrales, y la bandera -a los ligandos. Tanto los átomos centrales como los ligandos pueden ser de diferentes tipos; en este caso, deben estar separados por comas. La bandera -m da el tiempo mínimo que una especie debe vivir para ser considerada en el análisis. De forma predeterminada, este tiempo mínimo es cero, todas las ocurrencias se cuentan en el análisis final.
   1. Ejecute el script speciation_umd.py con el indicador –r 0, que muestra el gráfico de conectividad en el primer nivel para identificar los poliedros de coordinación. Por ejemplo, un átomo central, denotado como un catión , puede estar rodeado por uno o más aniones (Figura 4). La secuencia de comandos de especiación identifica cada uno de los poliedros de coordinación. La media ponderada de todos los poliedros de coordinación da el número de coordinación, idéntico al obtenido de la integración del PDF. En la línea de comandos, escriba:
      speciation_umd.py -f -i -c -a -r 0
      NOTA: Los números de coordinación promedio en fluidos son números fraccionarios. Esta fraccionalidad proviene de la característica media de la coordinación. La definición basada en la especiación produce una representación más intuitiva e informativa de la estructura del fluido, donde se cuantifican las proporciones relativas de las diferentes especies, es decir, las coordinaciones.
   2. Ejecute el script speciation_umd.py con el indicador –r 1, que muestrea el gráfico de conectividad en todos los niveles de profundidad para obtener la polimerización. La red a través del gráfico atómico tiene una cierta profundidad, ya que los átomos están unidos más lejos a otros enlaces (por ejemplo, en secuencias de cationes y aniones alternos) (Figura 4).
3. Abra los dos archivos . popul.dat y . stat.dat consecutivamente; estos constituyen la salida del script de especiación. Cada cúmulo está escrito en una línea, especificando su fórmula química, el momento en que se formó, el momento en que murió, su vida útil, una matriz con la lista de los átomos que forman este cúmulo. Trazar la vida útil de cada cúmulo atómico de todas las especies químicas encontradas en la simulación como se encuentra en el archivo .popul.dat (Figura 5).
4. Trazar el análisis poblacional con la abundancia de cada especie, tal y como se encuentra en el . archivo stat.dat . Este análisis, tanto absoluto como relativo, corresponde a las estadísticas reales de los poliedros de coordinación para el caso -r 0; para el caso de la polimerización, con -r 1 esto debe tratarse cuidadosamente, ya que podría ser necesario aplicar cierta normalización sobre el número relativo de átomos. La abundancia corresponde a la integral a lo largo de las vidas. El. stat.dat archivo también enumera el tamaño de cada clúster, es decir, cuántos átomos lo forman.

4. Calcular los coeficientes de difusión

1. Extraer los desplazamientos cuadráticos medios (MSD) de los átomos en función del tiempo para obtener la autodifusividad. La fórmula estándar del MSD es:
   
   donde los prefactores son renormalizaciones. Con la herramienta MSD, hay diferentes formas de analizar los aspectos dinámicos de los fluidos.
   NOTA: T es el tiempo total de la simulación y N _α es el número de átomos de tipo α. El tiempo inicial _t0 es arbitrario y abarca la primera mitad de la simulación. _Ninit es el número de veces iniciales. τ es la anchura del intervalo de tiempo durante el cual se calculan los TME; su valor máximo es la mitad de la duración de la simulación. En las implementaciones típicas de MSD, cada ventana comienza al final de la anterior. Pero un muestreo más escaso puede acelerar el cálculo del MSD, sin alterar la pendiente del MSD. Para esto, la ventana i-ésima comienza en el tiempo _t0(i), pero la ventana (i+1)-ésima comienza en el tiempo _t0(i) + τ + v, donde el valor de v está definido por el usuario. Del mismo modo, el ancho de la ventana se incrementa en pasos discretos definidos por el usuario, como tales: τ(i) = τ(i-1) + z. Los valores de z ("paso horizontal") y v ("paso vertical") son positivos o cero; el valor predeterminado para ambos es 20.
2. Calcule el MSD utilizando la serie de scripts de msd_umd . Su salida se imprime en un archivo . msd.dat archivo, donde el MSD de cada tipo atómico, átomo o clúster se imprime en una columna en función del tiempo.
   1. Calcule el MSD promedio de cada tipo atómico. Los MSD se calculan para cada átomo y luego se promedian para cada tipo atómico. El archivo de salida contiene una columna para cada tipo atómico. En la línea de comandos, escriba:
      msd_umd.py -f -z -v -b
   2. Calcular el MSD de cada átomo. Los MSD se calculan para cada átomo y luego se promedian para cada tipo atómico. El archivo de salida contiene una columna para cada átomo en la simulación y, a continuación, una columna para cada tipo atómico. Esta característica permite identificar átomos que se difunden en dos entornos diferentes, como líquido y gas, o dos líquidos. En la línea de comandos, escriba:
      msd_all_umd.py -f -z -v -b
   3. Calcular el MSD de la especie química. Utilice la población de clústeres identificados con el script de especiación e impresos en el archivo . popul.dat archivo. Los MSD se calculan para cada clúster individual. El archivo de salida contiene una columna para cada clúster. Para evitar considerar polímeros a gran escala, ponga un límite en el tamaño del clúster; su valor predeterminado es 20 átomos. En la línea de comandos, escriba:
      msd_cluster_umd.py -f -p -s -b -c
      NOTA: Los valores predeterminados son: –b 100 –s 1 –c 20.
3. Traza el MSD usando un software basado en hojas de cálculo (Figura 6). En una representación log-log del MSD versus el tiempo, identifique el cambio de pendiente. Separe la primera parte, generalmente corta, que representa el régimen balístico , es decir, la conservación de la velocidad de los átomos después de las colisiones. La segunda parte más larga representa el régimen difusivo , es decir, la dispersión de la velocidad de los átomos después de las colisiones.
4. Calcule los coeficientes de difusión a partir de la pendiente del MSD como:
   
   donde Z es el número de grados de libertad (Z = 2 para la difusión en el plano, Z = 3 para la difusión en el espacio), y t es el paso del tiempo.

5. Funciones de correlación de tiempo

1. Calcule las funciones de correlación de tiempo como una medida de la inercia del sistema utilizando la fórmula general:
   
   A puede ser una variedad de variables dependientes del tiempo, como las posiciones atómicas, las velocidades atómicas, las tensiones, la polarización, etc., cada una de las cuales produce, a través de las relaciones de Green-Kubo12,13, diferentes propiedades físicas, a veces después de una transformación adicional.
2. Analizar las velocidades atómicas para obtener el espectro vibratorio del líquido y la expresión alternativa de los coeficientes de autodifusión atómica.
   1. Ejecute el script vibr_spectrum_umd.py para calcular la función de autocorrelación velocidad atómica (VAC) para cada tipo atómico y realice su transformada de Fourier rápida. En la línea de comandos, escriba:
      vibr_spectrum_umd.py -f -t
      donde –t es la temperatura que debe definir el usuario. El script imprime dos archivos: el archivo . vels.scf.dat archivo con la función VAC para cada tipo atómico y el archivo . vibr.dat archivo con el espectro vibratorio descompuesto en cada especie atómica y el valor total.
   2. Abra y lea el .dat vels.scf. Trace la función VAC desde el archivo vels.scf.dat utilizando un software similar a una hoja de cálculo.
   3. Mantenga la parte real del Fourier VAC. Esto es lo que produce el espectro vibratorio, en función de la frecuencia:
      
      donde m son las masas atómicas.
   4. Traza el espectro vibratorio desde el archivo vibr.dat utilizando un software similar a una hoja de cálculo (Figura 7). Identificar el valor finito en ω=0 que corresponde al carácter difusivo del fluido y los diversos picos del espectro a frecuencia finita. Identificar la participación de cada tipo atómico en el espectro vibratorio.
      NOTA: La descomposición en tipos atómicos muestra que diferentes átomos tienen diferentes contribuciones ω=0, correspondientes a sus coeficientes de difusión. La forma general del espectro es mucho más suave con menos características que para un sólido correspondiente.
   5. En la cáscara, lea la integral sobre el espectro vibratorio, que produce los coeficientes de difusión para cada especie atómica.
      NOTA: Las propiedades termodinámicas se pueden obtener mediante la integración del espectro vibratorio, pero los resultados deben usarse con precaución debido a dos aproximaciones: la integración es válida dentro de la aproximación cuasi armónica, que no necesariamente se mantiene a altas temperaturas; y es necesario descartar la parte gaseosa del espectro correspondiente a la difusión. La integración debe hacerse solo sobre la parte del espectro similar a una red. Pero esta separación generalmente requiere varios pasos y cálculos posteriores al ^{procesamiento14}, que no están cubiertos por el presente paquete UMD.
3. Ejecute el script viscosity_umd.py para analizar la autocorrelación del tensor de tensión de los componentes para estimar la viscosidad de la masa fundida. En la línea de comandos, escriba:
   viscosity_umd.py -f -i -s -o -l
   NOTA: Esta característica es exploratoria y cualquier resultado debe tomarse con precaución. En primer lugar, comprobar a fondo la convergencia de la viscosidad con respecto a la longitud de la simulación.
   1. Derivar la viscosidad del fluido a partir de la autocorrelación del tensor de ^tensión15 como:
      
      donde V y T son el volumen y la temperatura respectivamente, _κB es la constante de Boltzmann y σ _ij la componente ij fuera de diagonal del tensor de tensión, expresada en coordenadas cartesianas.
   2. Utilizar un ajuste más adecuado para obtener una estimación más robusta de la ^{viscosidad15,16} y evitar el ruido de la función de autocorrelación tensor-tensión que pudiera derivarse del tamaño finito y la duración finita de las simulaciones. Para la función de autocorrelación del tensor de tensión, utilice la siguiente forma ^{funcional15,16} que produzca buenos resultados:
      
      donde A, B, _τ1, τ2 y ω son parámetros de ajuste. Después de la integración, la expresión para la viscosidad se convierte en:
      

6. Parámetros termodinámicos derivados de las simulaciones.

1. Ejecute averages.py para extraer los valores promedio y la propagación (como desviación estándar) de la presión, la temperatura, la densidad y la energía interna de los archivos umd . En la línea de comandos, escriba:
   averages.py -f -s
   con –s 0 como valor predeterminado.
2. Calcular el error estadístico del promedio, utilizando métodos de bloqueo.
   NOTA: Hay varios sabores de este método. Siguiendo el trabajo de Allen y Tildesley2, es común promediar sobre secuencias de bloques temporales, de longitud cada vez mayor, y estimar la desviación estándar con respecto a la media ^{aritmética17}. La convergencia puede alcanzarse en el límite de muchos tamaños de bloque y lo suficientemente largos, cuando el muestreo no está correlacionado. Aunque el valor umbral real para la convergencia generalmente debe elegirse manualmente.
   1. Utilice el método de reducción a la ^mitad18: a partir de la muestra de datos inicial, en cada paso κ, reduzca a la mitad el número de muestras promediando cada dos muestras consecutivas correspondientes del paso anterior κ−1:
      
   2. Ejecute el script fullaverages.py para realizar el análisis estadístico completo, incluido el error de la media. En la línea de comandos, escriba:
      fullaverages.py -s -u
      NOTA: El script se automatiza hasta el punto de buscar todos los archivos .umd.dat en el directorio actual y realizar el análisis para todos ellos. Los valores predeterminados son –s 0 –u 0. Para -u 0 la salida es mínima, y para -u 1 la salida es completa, con varias unidades alternativas impresas. Este script requiere soporte gráfico, ya que crea una imagen gráfica para comprobar la convergencia para estimar el error en la media.

Resultados

La pirolita es un modelo de fusión de silicato multicomponente (0.5Na2O 2CaO 1.5Al2O3 4FeO 30MgO 24SiO2) que mejor se aproxima a la composición del silicato a granel de la Tierra, el promedio geoquímico o de todo nuestro planeta, excepto su núcleo a base de ^hierro19. La Tierra primitiva estuvo dominada por una serie de eventos de fusión a gran ^escala20, el último podría haber envuelto a todo el planeta, después de su condensación para el disco protolunar21. La pirolita representa la mejor aproximación a la composición de tales océanos de magma a escala planetaria. En consecuencia, estudiamos ampliamente las propiedades físicas de la fusión de pirolita en el rango de temperatura de 3.000\u20125.000 K y el rango de presión de GPa de 0\u2012150 a partir de simulaciones de dinámica molecular ab initio en la implementación de VASP. Estas condiciones termodinámicas caracterizan completamente las condiciones oceánicas de magma más extremas de la Tierra. Nuestro estudio es un excelente ejemplo de un uso exitoso del paquete UMD para todo el análisis en profundidad de los ^fundidos22. Calculamos la distribución y las longitudes medias de los enlaces, trazamos los cambios en la coordinación catiónico-oxígeno y comparamos nuestros resultados con estudios experimentales y computacionales previos sobre silicatos amorfos de diversas composiciones. Nuestro análisis en profundidad ayudó a descomponer los números de coordinación estándar en sus constituyentes básicos, describir la presencia de poliedros de coordinación exóticos en la masa fundida y extraer vidas útiles para todos los poliedros de coordinación. También destacó la importancia del muestreo en simulaciones en términos de longitud de la trayectoria y también del número de átomos presentes en el sistema que se modela. En cuanto al post-procesamiento, el análisis UMD es independiente de estos factores, sin embargo, deben tenerse en cuenta al interpretar los resultados proporcionados por el paquete UMD. Aquí, mostramos algunos ejemplos de cómo el paquete UMD se puede utilizar para extraer varios rasgos característicos de las fundiciones, con una aplicación a la pirolita fundida.

La función de distribución de pares Si-O obtenida del script gofrs_umd.py muestra que el radio de la primera esfera de coordinación, que es el primer mínimo de la función g(r), se encuentra alrededor de 2,5 angstroms en T = 3000 K y P = 4,6 GPa. El máximo de g(r) se encuentra en 1.635 Å, esta es la mejor aproximación a la longitud de la curva. La cola larga se debe a la temperatura. Usando este límite como la distancia de enlace si-O, el análisis de especiación muestra que las unidades de _SiO4, que pueden durar hasta unos pocos picosegundos, dominan la fusión (Figura 5). Hay una parte importante de la masa fundida que muestra polimerización parcial, como lo refleja la presencia de dímeros como _Si2O7 y trímeros como las unidades _Si3Ox. Su vida útil correspondiente es del orden del picosegundo. Todos los polímeros de orden superior tienen una vida útil considerablemente más corta.

Figura 5: Vida útil de la especie química Si-O.
La especiación se identifica en una fusión multicomponente a 4,6 GPa y 3000 K. Las etiquetas marcan los monómeros _siO3, _SiO4 y _SiO5 y los diversos polímeros _SixOy. Haga clic aquí para ver una versión más grande de esta figura.

Los diferentes valores de los pasos verticales y horizontales, definidos por las banderas –z y –v anteriores, producen varios muestreos del MSD (Figura 6). Incluso grandes valores de z y v son suficientes para definir las pendientes y, por lo tanto, los coeficientes de difusión de los diferentes átomos. La ganancia en el tiempo para el post-procesamiento es notable cuando se va a valores grandes de z y v. El MSD ofrece un criterio de validación muy fuerte para la calidad de las simulaciones. Si la parte de difusión del MSD no es lo suficientemente larga, eso es una señal de que la simulación es demasiado corta y no alcanza el estado fluido en un sentido estadístico. El requisito mínimo para la parte difusiva del MSD depende en gran medida del sistema. Se puede requerir que todos los átomos cambien su sitio al menos una vez en la estructura de la masa fundida para que sea considerada como un ^fluido10. Un excelente ejemplo con aplicaciones en ciencias planetarias son las complejas fundiciones de silicato a altas presiones cercanas o incluso por debajo de su línea ^liquidus11. Los átomos de Si, los principales cationes formadores de redes, cambian de sitio después de más de dos docenas de picosegundos. Las simulaciones más cortas que este umbral estarían considerablemente inframuestreando el posible espacio de configuración. Sin embargo, como los aniones coordinadores, es decir, los átomos de O, se mueven más rápido que los átomos centrales de Si, pueden compensar parte de la lenta movilidad del Si. Como tal, todo el sistema podría cubrir un mejor muestreo del espacio de configuración que el asumido solo a partir de los desplazamientos de Si.

Figura 6: Desplazamientos cuadráticos medios (TME).
Los MSD se ilustran para algunos tipos atómicos de una fusión de silicato de múltiples componentes. El muestreo con varios pasos horizontales y verticales, z y v, produce resultados consistentes. Círculos sólidos: -z 50 –v 50. Círculos abiertos: -z 250 –v 500. Haga clic aquí para ver una versión más grande de esta figura.

Finalmente, las funciones de VCA atómica producen el espectro vibratorio de la masa fundida. La Figura 7 muestra el espectro en las mismas condiciones de presión y temperatura que las anteriores. Representamos las contribuciones de los átomos de Mg, Si y O, así como el valor total. A la frecuencia cero hay un valor finito del espectro, que corresponde al carácter de difusión de la masa fundida. La extracción de las propiedades termodinámicas del espectro vibratorio debe eliminar este carácter difusivo similar al gas de cero, pero también tener en cuenta adecuadamente su desintegración a frecuencias más altas.

Figura 7: El espectro vibratorio de la fundición de pirolita.
La parte real de la transformada de Fourier de la función de autocorrelación velocidad-velocidad atómica produce el espectro vibratorio. Aquí el espectro se calcula para una fusión de silicato de múltiples componentes. Los fluidos tienen un carácter difusivo similar al gas distinto de cero a una frecuencia cero. Haga clic aquí para ver una versión más grande de esta figura.

Discusión

El paquete UMD ha sido diseñado para funcionar mejor con simulaciones ab initio, donde el número de instantáneas generalmente se limita a decenas a cientos de miles de instantáneas, con unos pocos cientos de átomos por unidad de celda. Las simulaciones más grandes también son manejables siempre que la máquina en la que se ejecuta el posprocesamiento tenga suficientes recursos de memoria activos. El código se distingue por la variedad de propiedades que puede calcular y por su licencia de código abierto.

Los archivos umd.dat son apropiados para los conjuntos que conservan el número de partículas sin cambios a lo largo de la simulación. El paquete UMD puede leer archivos derivados de cálculos donde varía la forma y el volumen del cuadro de simulación. Estos cubren los cálculos más comunes, como NVT y NPT, donde el número de partículas, N, temperatura T, volumen, V y / o presión, P, se mantienen constantes.

Para el momento comienzan la función de distribución de pares, así como todas las escrituras que necesitan estimar las distancias interatómicas, como las escrituras de especiación, funcionan solo para células unitarias ortogonales, es decir, para células cúbicas, tetragonales y ortorrómbicas, donde los ángulos entre los ejes son de 90 °.

Las principales líneas de desarrollo para la versión 2.0 son la eliminación de la restricción de ortogonalidad para distancias y la adición de más características para los scripts de especiación: analizar enlaces químicos individuales, analizar los ángulos interatómicos e implementar la segunda esfera de coordinación. Con la ayuda de la colaboración externa, estamos trabajando en la migración del código a una GPU para un análisis más rápido en sistemas más grandes.

Materiales

Name	Company	Catalog Number	Comments
getopt library			open-source
glob library			open-source
matplotlib library			open-source
numpy library			open-source
os library			open-source
Python software	The Python Software Foundation	Version 2 and 3	open-source
random library			open-source
re library			open-source
scipy library			open-source
subprocess library			open-source
sys library			open-source