Parámetros de los modelos de energía potencial

Las simulaciones de dinámica molecular fueron realizadas usando el ensamble tipo NPT (número de moléculas, presión y temperatura constante) (^{Frenkel & Smit, 2002}). La molécula de n-octano fue descrita con el campo de fuerza all-atom GROMOS53A6 (^{Oostenbrink et al., 2004}), donde la forma del funcional para calcular la energía total, E , es definida por las ecuaciones 1 y 2:

Donde r y r ₀ corresponden a las longitudes del enlace entre los átomos en el tiempo de simulación y en equilibrio. Los términos θ y θ ₀ son los ángulos de enlace entre tres átomos en el tiempo de simulación y en equilibrio. El término φ corresponde al ángulo torsional formado entre cuatro átomos. Aquí, C12_ij y C6_ij son los parámetros de interacción entre pares de átomos i y j, los cuales están asociados a la energía de interacción tipo Lennard-Jones. Estos términos son calculados usando las siguientes relaciones: c12ij=c12i*c12j y c6ij=c6i*c6j. Además, ε es el coeficiente de permitividad en el vacío, las variables q_i y q_j son las cargas definidas para los átomos y r_ij es la distancia entre los átomos.

Los modelos SPC flexible, SPC-E flexible y TIP3P flexible fueron los campos de fuerzas usados para describir las moléculas de agua en los sistemas. En estos modelos de tres centros, los átomos de hidrógeno tienen cargas positivas, el átomo de oxígeno tiene carga negativa y el potencial de interacción describe únicamente la interacción entre los pares de átomos de oxígeno. La energía intermolecular del agua es definida por la ecuación 2. En la tabla1, se presenta los parámetros de los campos de fuerzas del agua.

Particularmente, estos parámetros pueden ser modificados y ajustados por el usuario en función del campo de fuerza que desea utilizar. Estos parámetros han sido reportados en diferentes trabajos y pueden ser agregados o modificados por el usuario en los archivos con extensión itp y top presentes en el programa GROMACS-4.5.4 (^{Malde et al., 2011}; ^{Abraham et al., 2019}).

Construcción del sistema agua/n-octano

Para las simulaciones con el software GROMACS-4.5.4 se requieren de los archivos que contienen la información de la configuración de las moléculas en una celda (archivo.pdb), los parámetros del campo de fuerza (archivo.top) y las condiciones de la simulación (archivo.mdp) (^{Abraham et al., 2019}).

Inicialmente, el programa Avogadro-1.1.1 (^{Hanwell et al., 2012}) fue utilizado para la construcción de las geometrías moleculares de agua y n-octano (ver figura 1(a) y 1(b)). Estas geometrías fueron almacenadas en el formato de archivo Protein Data Bank (PDB). Luego, el software GROMACS-4.5.4 fue empleado en la construcción, preparación y simulación del sistema agua/n-octano. Para ello, usando la herramienta genbox (^{van der Spoel et al., 2010}) se construyeron dos celdas rectangulares con condiciones límites periódicas en todos los ejes de coordenadas y con dimensiones de 4×4×3nm³. Esta construcción fue realizada utilizando las densidades experimentales de 0.997 g/cm³ y 0.699 g/cm³ a 298.15 K para el agua y n-octano, respectivamente (^{Yaws, 2003}). Con estas densidades, una celda contiene 180 moléculas de n-octano y la otra celda 1600 moléculas de agua. Las siguientes instrucciones se utilizaron en una terminal para la construcción de las celdas:

a) $ genbox -ci agua.pdb -nmol 1600 -box 4 4 3 -o aguabox.pdb

b) $ genbox -ci octano.pdb -nmol 180 -box 4 4 3 -o octanobox.pdb

En estas celdas, las moléculas están distribuidas de manera aleatoria (ver figura 1(c) y 1(d)).

La salida en la terminal muestra lo siguiente:

a) Output configuration contains 4800 atoms in 1600 residues

Density: 997.165 (g/l)

Number of SOL molecules: 1600

b) Output configuration contains 4680 atoms in 180 residues

Density: 699.00 (g/l)

Number of RES molecules: 180

Figura 1 (a) Molécula de agua, (b) Molécula de n-octano, (c) Celda con 1600 moléculas de agua y (d) Celda con 180 moléculas de n-octano. 

Luego, las celdas fueron modificadas con la herramienta editconf (^{van der Spoel et al., 2010}) y usadas para la construcción del sistema agua/n-octano. Las siguientes instrucciones fueron usadas:

a) $ editconf -f aguabox.pdb -box 4 4 10 -center 2.0 2.0 5.0 -o aguabox-centro.pdb

b) $ editconf -f octanobox.pdb -box 4 4 10 -center 2.0 2.0 1.75 -o octanobox-izquierda.pdb

c) $ editconf -f octanobox.pdb -box 4 4 10 -center 2.0 2.0 8.25 -o octanobox-derecha.pdb

En este caso, la celda del sistema agua/n-octano tiene dimensiones de 4×4×10nm³. El centro de la celda está localizado en 2.0 nm en los ejes xy. La tercera coordenada corresponde a la posición a lo largo del eje z. En la figura 2(a) and 2(b), se muestran las celdas de agua y n-octano modificadas.

Para unir las celdas, se utilizaron las siguientes instrucciones que permiten la obtención del sistema agua/n-octano usado en las simulaciones:

a) $ genbox -cp octanobox-izquierda.pdb -cs aguabox-centro.pdb -o sistema1.pdb

b) $ genbox -cp sistema1.pdb -cs octanobox-derecha.pdb -o sistema-agua-n-octano.pdb

Aquí, el comando -cp identifica al soluto y el comando -cs al solvente. El comando -o indica en que archivo se guarda la información de la configuración inicial del sistema (sistema-agua-n-octano.pdb). En este caso, en el formato PDB. La figura 2(c), muestra la configuración del sistema agua/n-octano, el cual contiene dos regiones interfaciales.

Figura 2 (a) Celda de agua localizada en el centro del sistema, (b) Celdas de n-octano localizadas en los extremos del sistema y (c) Sistema agua/n-octano obtenido en la construcción. 

Preparación de las topologías moleculares

Para la simulación del sistema agua/n-octano son necesarios los parámetros asociados a la ecuación 1 y 2 que permiten determinar la energía total del sistema. Este archivo es denominado topología molecular (^{van der Spoel et al., 2010}). Los parámetros que describen a las moléculas de agua y n-octano corresponden a los modelos all-atom GROMOS53A6, SPC, SPC-E y TIP3P. El siguiente texto muestra el archivo de topología “archivo.top” utilizado:

; Include force field parameters

#include “/home/gromacs-4.5.4/share/top/gromos53a6.ff/forcefield.itp”

#include “/home/gromacs-4.5.4/share/top/gromos53a6.ff/octano.itp”

#include “/home/gromacs-4.5.4/share/top/gromos53a6.ff/spc.itp”

[system]

agua-octano

[molecules]

; Compound #mols

oct 180

SOL 1600

oct 180

El comando #include define sirve para leer la dirección completa en la cual se encuentran localizados archivos externos que son necesarios para realizar la simulación. En este caso, nos permite leer la dirección completa donde están localizados los archivos itp. Estos archivos contienen los parámetros enlazantes y no enlazantes de las moléculas (el archivo octano.itp fue obtenido desde el servidor ATB. (^{Malde et al., 2011}). También, en el archivo de topología “archivo.top” es necesario indicar el número de moléculas presentes en el sistema. Esta cantidad debe coincidir con el número de moléculas presentes en el archivo de configuración inicial “sistema-agua-n-octano.pdb”.

Condiciones de las simulaciones

Las condiciones de las simulaciones están contenidas en el archivo con la extensión mdp (ver archivo agua-octano.mdp en anexo 1). El ensamble NPT (número de moléculas, presión y temperatura constante) fue usado en las simulaciones. Los sistemas son periódicos en las coordenadas xyz. La temperatura fue de 300 K y la presión de 1 atm. El método de escalamiento de la velocidad (^{Bussi, Donadio & Parrinello, 2007}) fue usado para controlar la temperatura y el método de ^{Berendsen (Berendsen et al., 1984}) para la presión. La constante del tiempo para el termostato fue de 0.1 ps y para el barostato fue de 1.0 ps. Las velocidades iniciales fueron generadas usando una distribución Maxweliana a 300 K y las ecuaciones de movimiento fueron integradas usando el algoritmo de Verlet con un paso del tiempo de 1 fs. Las interacciones de Lennard-Jones fueron calculadas con un radio de 1.50 nm y las interacciones electrostáticas con un radio de 1.40 nm con el método Particle Mesh Ewald (PME) (^{Essmann et al., 1995}). Además, correcciones de dispersión de largo alcance fueron aplicados para la presión y la energía. Al principio, los sistemas fueron optimizados con el método de gradiente conjugado, el cual consiste en la minimización de la energía variando las coordenadas de posición de los átomos en la celda. Luego, las configuraciones finales obtenidas en el proceso de minimización fueron utilizadas en las simulaciones de dinámica molecular. Para obtener los sistemas en equilibrio, las simulaciones fueron realizadas por un tiempo de 20 ns. Los últimos 10 ns fueron usados para estimar las propiedades como la tensión interfacial y el espesor de la región interfacial agua/n-octano. Las trayectorias fueron almacenadas cada 1000 fs. Finalmente, la visualización del sistema agua/n-octano fue realizada con el programa pymol y las gráficas de los perfiles de densidad requirieron el uso del programa xmgrace.

El programa GROMACS-4.5.4 utiliza la herramienta grompp (^{van der Spoel, 2010}) para hacer el pre-procesamiento de los archivos con las extensiones mdp, top y pdb. La instrucción se muestra a continuación:

a) $ grompp -f agua-octano.mdp -p agua-octano.top -c sistema-agua-n-octano.pdb -o topol.tpr

Esto produce el archivo “topol.tpr” que almacena la información de los archivos “agua-octano.mdp”, “agua-octano.top” y “sistema-agua-n-octano.pdb”, la cual es necesaria para realizar la simulación. Luego, se ejecuta la siguiente instrucción para realizar la optimización ó la dinámica molecular:

b) $ mdrun -v -s topol.tpr -np 1 -nt 8

Donde -np indica el número de procesadores y -nt el número de hilos a utilizar.

Estimación de la tensión interfacial y el espesor de película interfacial

Para el sistema inmiscible agua/n-octano que posee una región interfacial, la tensión interfacial (γ) fue determinada usando los componentes tangencial (P _t = (P _xx + P _yy)/2) y normal (P _N = P _zz) al plano interfacial del tensor presión. En este caso, la tensión interfacial es definida por la ecuación 3:

Donde, el término L _z es la longitud del eje z de la celda. Este modelo fue propuesto por ^{Kirkwood y Buff (Kirkwood & Buff, 1949}) y depende de los de los parámetros de los campos de fuerza utilizados. Los elementos del tensor presión molecular están definidos por la ecuación 4:

Aquí I representa el tensor unitario, T es la temperatura, y ρ = N/V es la densidad en número. Los términos i y j representan las direcciones x, y, z. Además, r_kl representa el vector entre los centros de masa de las moléculas k y l. En la ecuación 4, F_kl es la fuerza intermolecular entre las moléculas k y l.

Por otro parte, la teoría onda-capilar (TCW) (^{Rowlinson & Widom, 2002}) evalúa las características de las interfaces líquido-líquido. La TCW asume que el espesor de la película interfacial (σ) contiene dos componentes asociados entre sí: la amplitud intrínseca (σ₀) de la interfaz y la amplitud capilar (σ_cw) asociada a las fluctuaciones de la interfaz debido a la temperatura (ondas capilares). El espesor de la película interfacial de acuerdo con la TCW viene dado por la ecuación 5:

La amplitud capilar (σ_cw) está relacionada con la tensión interfacial y las dimensiones de la celda del sistema usada en la simulación por la ecuación 6 (^{Nicolas & de Souza, 2004}):

Donde k_B es la constante de Boltzmann, L_ii es la longitud del eje x ó y de la celda de simulación; y ξ es la longitud de correlación molecular del hidrocarburo.

Adicionalmente, el perfil de densidad de los sistemas fue determinado de las simulaciones. En este caso, el perfil de densidad de cualquier componente en el sistema puede ser descrito por la ecuación 7 (^{Senapati & Berkowitz, 2001}):

En la ecuación 7, el parámetro es la posición del plano divisorio de Gibbs (GDS, por sus siglas en inglés) en el eje z usando el componente A y ρ^A(z) y es la densidad del componente A. La amplitud intrínseca de la interfaz (σ₀) se genera cuando las dos fases se entremezclan y constituyen la distancia entre los dos planos GDS. Esto es definido por la ecuación 8:

Donde Z₀^agua y Z₀^hidrocarburo son los planos divisorios de Gibbs del agua y del n-octano, respectivamente.

También, se han propuesto otros modelos para calcular este espesor usando ciertas propiedades estructurales (^{Mitrinovic et al., 1999}). Esto se muestra en las ecuaciones 9 y 10:

Donde R_g es el radio de giro de la molécula de n-octano obtenido de la simulación.

Perfiles de densidad

Los perfiles de densidad son obtenidos utilizando la herramienta g_density (^{van der Spoel et al., 2010}) con la siguiente instrucción:

$g_density -sl 600 -ng 3 -d z -n index.ndx -b 10000 -f traj.trr -o densidad.xvg

Esta herramienta g_density lee la información de los últimos 10 ns contenidos en el archivo de traj.trr, divide la celda en 600 porciones a lo largo del eje z y determina el perfil de densidad para tres componentes en el sistema. En este caso, el perfil de densidad es obtenido para el agua, n-octano y el sistema completo. Esto es almacenado en el archivo densidad.xvg, el cual es graficado usando el programa xmgrace. En la figura 3(a) se muestra el perfil de densidad del sistema agua/n-octano donde se identifican las dos regiones interfaciales. En esta figura 3(a), se puede observar la variación de las densidades de las capas de agua y n-octano en la región interfacial.

Las figuras 3(b), 3(c) y 3(d) presentan los perfiles de densidad del sistema agua/n-octano usando las diferentes modelos que describen las moléculas de agua. Cada figura muestra una capa de agua y dos capas de n-octano. Estas capas están en contacto delimitando dos regiones interfaciales bien definidas. Esto indica que los campos de fuerza utilizados producen la separación de fases entre líquidos inmiscibles y, a su vez, la presencia de la región interfacial entre el agua y el n-octano. Para los modelos SPC y SPC-E, los valores de densidad del agua a lo largo del eje z están por encima de la densidad reportada experimentalmente de 997 kg/m³ (ver figura 3b y 3c). Los valores de densidad de la capa de agua pura en el sistema agua/n-octano están alrededor de 1020 kg/m³ y 1040 kg/m³ con los modelos SPC y SPC-E, respectivamente. Esto indica una contracción de la capa de agua ocasionada por las dos capas de n-octano a lo largo del eje z. En cambio, el modelo TIP3P muestra una tendencia muy cercana a la densidad experimental del agua (alrededor de 998 kg/m³), en donde no se observa el efecto de contracción de capa (ver figura 3d).

Figura 3 Perfiles de densidad del sistema agua/n-octano usando diferentes combinaciones de modelos de energía potencial. (a) Perfil de densidad que muestra las regiones interfaciales, (b) agua descrita con el modelo SPC, (c) agua descrita con el modelo SPC-E y (d) agua descrita con el modelo TIP3P. 

En todos los sistemas simulados, el perfil de densidad para la molécula de n-octano presenta considerables fluctuaciones cerca de la región interfacial. Estas fluctuaciones son producidas por las interacciones repulsivas (según parámetros de Lennard-Jones) entre grupos metilos y metilenos presentes en las moléculas de n-octano con el oxígeno del agua en la región interfacial. Esto produce un reacomodo de las moléculas de n-octano en la región interfacial con el fin de minimizar la repulsión entre los grupos presentes en las moléculas. Estas fluctuaciones son más evidentes y extendidas en la región interfacial de los sistemas donde el agua es descrita con los modelos SPC y SPC-E. En estos modelos SPC y SPC-E, las interacciones repulsivas entre átomos de oxígeno con los grupos -CH₂ y -CH₃ son considerables y controlan el arreglo espacial del n-octano con las moléculas de agua en la interface (ver tabla 1). Finalmente, en el seno del agua y n-octano, los perfiles de densidad son uniformes. Los resultados obtenidos en esta parte deben ser consistentes con los reportados en otros trabajos (^{Zhang et al., 1995}; Wick et al., 2011; ^{Neyt et al., 2014}).

Plano divisorio de Gibbs y espesores de película interfacial

Las densidades de los sistemas fueron obtenidas usando la herramienta g_energy (^{van der Spoel et al., 2010}). Para ello, se ejecuta la siguiente instrucción en un terminal; seleccionando la propiedad que será determinada a partir de los 10000 ps:

$g_energy -f ener.edr -b 10000 -o densidad.xvg

Los valores de densidad para cada sistema son mostrados en la tabla 1. Las densidades de los sistemas están comprendidas entre los valores del agua (997 kg/m³) y n-octano (699 kg/m³). En los sistemas estudiados, las desviaciones de la densidad de los sistemas se reducen a lo largo de la simulación hasta hacerse constantes (ver tabla 1). Por lo tanto, estos sistemas se encuentran en equilibrio dinámico. El sistema donde el modelo SPC es usado para el agua presenta la menor desviación estándar de la densidad.

A su vez, la posición más probable de los planos divisorios de Gibbs (GDS, por sus siglas en inglés) fueron determinados ajustando los datos de los perfiles de densidad (contenidos en el archivo density.xvg) a la ecuación 6 mediante la herramienta Solver contenida en el programa libreoffice-6.0. En este caso, por la presencia de dos regiones interfaciales 1 y 2 (ver figura 3(a)) tenemos dos GDS en el sistema estudiado. Los valores para los GDS de la capa de agua y n-octano son mostrados en la tabla 1.

Las posiciones de los planos GDS fueron iguales cuando las moléculas de agua son descritas con los modelos SPC y SPC-E (ver tabla 1). Esto indica que no existen diferencias entre estos modelos para describir la región interfacial de los sistemas estudiados en una celda de simulación con condiciones límites periódicas en todos los ejes de coordenadas. En cambio, la posiciones de las GDS cuando el agua es descrita con el modelo TIP3P están ligeramente desplazados a lo largo del eje z.

La diferencia entre las posiciones de los GDS del agua y n-octano permiten determinar la amplitud intrínseca. Este término tiene un valor promedio de 0.07 nm para todas las regiones interfaciales. A su vez, se observa que esta propiedad no varía con el modelo de energía potencial usado para la descripción de las moléculas de agua. Los resultados encontrados son consistentes con los obtenidos por Rivera et al (^{Rivera, McCabe & Cummings, 2003}).

De la misma manera, usando la ecuación 7 fueron obtenidos los valores de σ_cw (ver tabla 2). El valor de σ_cw varía en función del modelo de energía potencial usado para describir las moléculas de agua. De hecho, no existe diferencia en los valores de σ_cw con respecto a la capa de líquido usada en la estimación. Los valores más bajos de fueron obtenidos cuando las moléculas de agua fueron descritas con el modelo SPC-E. En este trabajo, los valores más precisos fueron obtenidos cuando las moléculas de agua están descritas por el modelo TIP3P. Por ejemplo, Rivera et al (obtuvieron un valor de 0.148 nm para el sistema agua/n-octano.

El espesor de película interfacial total fue determinado con las ecuaciones 9 y 10. Para ello, el radio de giro de la molécula de n-octano (valor promedio de 0.282 nm) obtenido de las simulaciones fue empleado. Con la ecuación 9, los valores obtenidos fueron de 0.312 nm, 0.306 nm y 0.318 nm para los sistemas agua(SPC)/n-octano, agua(SPC-E)/n-octano y agua(TIP3P)/n-octano, respectivamente. En cambio, los valores de 0.339 nm, 0.307 nm y 0.307 nm fueron obtenidos para los sistemas usando la ecuación 10. Esta diferencia entre los valores obtenidos es debido al valor de radio de giro usada en la ecuación 10. Los resultados son consistentes con los valores de 0.344 nm y 0.322 nm reportados por Mitrinovic et al (^{Mitrinovic et al, 2000}) y Zhang et al (^{Zhang et al, 1995}), respectivamente.

Tensión interfacial

La tensión interfacial es una de las propiedades más importantes medidas en el laboratorio de fisicoquímica. En este trabajo, la tensión interfacial del sistema agua/n-octano fue estimada usando la ecuación 3 y los componentes diagonales del tensor presión obtenidos de la simulación. Aquí, la siguiente instrucción fue utilizada (^{van der Spoel et al., 2010}):

$ g_energy -b 10000 -f ener.edr -o tension-interfacial.xvg

En el terminal se seleccionan las opciones que generan los componentes de presión y la tensión interfacial. Además, con los valores σcw de fueron determinadas las tensiones interfaciales usando la ecuación 5. Los valores de tensión interfacial son mostrados en la tabla 3.

La tensión interfacial experimental del sistema agua/n-octano tiene un valor de 50.80mN/m a 298.15 K (^{Freitas, Quina & Carroll, 1997}). En este trabajo, con el modelo de Kirwood-Buff, la combinación de SPC y SPC-E con el modelo all-atom GROMOS53A6 produce una sobreestimación de la tensión interfacial, lo cual puede ser debido a la combinación de los parámetros de Lennard-Jones de las moléculas involucradas en la simulación (^{Neyt et al., 2014}; ^{Underwood & Greenwell, 2018}). Aquí, el error es más considerable cuando se utiliza el modelo SPC-E. En cambio, con el modelo TIP3P se obtiene el mejor valor de tensión interfacial.

Además, la tensión interfacial fue calculada usando la TCW. En este caso, la tensión interfacial depende de la longitud de correlación del hidrocarburo (ξ), la cual fue estimada usando la ecuación 6. Para obtener una buena tensión interfacial se requieren de altos valores de la longitud de correlación del hidrocarburo (ver tabla 3) que exceden la longitud molecular del octano (~0.72 nm). En este caso, estos valores se desvían considerablemente de los reportados en otros trabajos. Por ejemplo, Cordeiro (^{Cordeiro, 2002}) y Riedleder et al (^{Riedleder et al., 2007}) reportaron valores de 0.650 nm y 0.53 nm para los sistemas agua/isoctano y agua/n-heptano. Aquí, la mejor estimación de la longitud de correlación del hidrocarburo fue obtenida con el modelo TIP3P para el agua. Además, si modificamos su valor a 0.65 nm se obtiene un valor promedio de tensión interfacial igual a 54.30 mN/m para el sistema agua/n-octano donde el agua fue descrita con el modelo TIP3P. Finalmente, los resultados obtenidos indican que las simulaciones con parámetros a nivel molecular son adecuadas para estimar la tensión interfacial y el espesor de la región interfacial de un sistema formado por dos líquidos inmiscibles. Sin embargo, para obtener valores más exactos de las propiedades interfaciales se requiere de la evaluación de la regla de combinación y un reajuste de los parámetros de Lennard-Jones asociados a la interacción intermolecular entre las moléculas de agua y n-octano en el sistema.