DOI: http://dx.doi.org/10.22201/fq.18708404e.2021.3.76027
José Gregorio Parra Figueredo[a], Peter Iza[b] y Elizabeth Perozo[c]
Resumen |
La dinámica molecular es una herramienta muy útil para estudiar sistemas complejos y determinar propiedades de líquidos inmiscibles entre sí. En este trabajo se presenta un procedimiento computacional usando el programa GROMACS-4.5.4, el cual fue empleado en un laboratorio de fisicoquímica, para estimar la tensión interfacial y el espesor de película interfacial del sistema agua/n-octano con el modelo de Kirkwood-Buff y la teoría onda-capilar. En estas simulaciones, los componentes de la presión fueron utilizados para estimar la tensión interfacial del sistema agua/n-octano. En este estudio, los mejores resultados para la tensión interfacial fueron obtenidos con los campos de fuerza TIP3P y all-atom GROMOS53A6 usados en la descripción del agua y n-octano, respectivamente. |
Palabras clave: |
Dinámica molecular, película interfacial, tensión interfacial. |
A guide for the estimation of the interfacial properties of a water/hydrocarbon system using the software gromacs-4.5.4 |
|
Abstract |
Molecular dynamics is a very useful tool for studying complex systems and determining properties of liquids that are immiscible with each other. In this work, a computational procedure is presented using the GROMACS-4.5.4 program, which was used in a physicochemical laboratory, to estimate the interfacial tension and the interfacial film thickness of the water/n-octane system with the Kirkwood-Buff model and wave-capillary theory. In the simulations, the pressure components were used to estimate the interfacial tension of the water/n-octane system. In this study, the best results for interfacial tension were obtained with the TIP3P and all-atom GROMOS53A6 force fields used in the description of water and n-octane, respectively. |
Keywords: |
Molecular dynamics, interfacial film, interfacial tension. |
Introducción
La determinación de la tensión interfacial y el espesor de la región interfacial de un sistema formado por dos líquidos inmiscibles es fundamental para el entendimiento de fenómenos involucrados en los procesos de ingeniería química; como la extracción con disolventes, separación física de compuestos mediante membranas y preparación de emulsiones (Holmberg, 2002). Un sistema formado por dos líquidos inmiscibles posee una región interfacial que se caracteriza por tener propiedades termodinámicas como la tensión interfacial (Lyklema, 2000). Experimentalmente, la determinación de la tensión interfacial de un sistema agua/hidrocarburo puede ser medido usando un tensiómetro de Du Noüy, el cual consiste en un anillo de platino colgado en el brazo de una balanza de torsión (Castellan, 1987). Para ello, el anillo se sumerge hasta la zona donde se forma la región interfacial y los líquidos inmiscibles se encuentran en contacto con la parte interior y exterior del anillo (Castellan, 1987; Moore, 1972). La fuerza necesaria para remover el anillo desde la región interfacial está determinada por la relación f=4πrγ, donde f es la fuerza aplicada; γ la tensión interfacial y r el radio del anillo (Butt, Graf, & Kappl, 2003).
En el laboratorio de fisicoquímica, los estudiantes realizan la medición de la tensión interfacial de un sistema agua/aceite con un tensiómetro de Du Noüy. Sin embargo, ellos tienen problemas para colocar el anillo de platino en la región interfacial y determinar la ruptura del contacto que existe entre los dos líquidos inmiscibles y el anillo de platino. Este hecho ocasiona errores en las mediciones y deben ser repetidas varias veces para obtener resultados aceptables de los sistemas evaluados.
A nivel atomístico, la dinámica molecular permite estimar las propiedades macroscópicas de sistemas moleculares (Allen & Tildesley, 1989). En esta técnica, las coordenadas y velocidades de los átomos son evaluadas en el tiempo con las ecuaciones de movimiento de Newton (Frenkel & Smit, 2002). A su vez, las moléculas son descritas por un campo de fuerza (Force Field en inglés), los cuales incluyen parámetros asociados a los términos enlazantes y no enlazantes de la energía total del sistema (Lewars, 2016). En este sentido, las propiedades estimadas mediante simulaciones moleculares pueden ser comparadas con las mediciones obtenidas experimentalmente con la finalidad de mostrar a los estudiantes nuevos recursos para evaluar los fenómenos interfaciales. Particularmente, la dinámica molecular ha sido utilizada para explorar y evaluar las propiedades interfaciales como la tensión interfacial, el espesor de la región interfacial y la orientación molecular del agua y del hidrocarburo en la interface con fines de investigación especializada en diferentes trabajos (Zhang et al., 1995; Nicolas & de Souza, 2004; Riedleder et al., 2007; Wick et al., 2011; Neyt et al., 2014). Sin embargo, como herramienta para fortalecer los conocimientos básicos de estudiantes de pre-grado en fisicoquímica no ha sido utilizada.
Por tal motivo, en este trabajo se presenta un procedimiento de dinámica molecular donde el estudiante de química de educación superior sea capaz de determinar la tensión interfacial y el espesor de la región interfacial de un sistema agua/hidrocarburo con la finalidad de complementar sus conocimientos a nivel molecular de los sistemas agua/hidrocarburo. En este caso, las propiedades fueron estimadas para el sistema agua/n-octano.
El software de código abierto GROMACS-4.5.4(http://www.gromacs.org ) (Hess et al., 2008) fue utilizado en la construcción, preparación y simulación de los sistemas. Los potenciales de carga puntual simple (SPC) flexible (Berendsen et al., 1981), carga puntual simple extendida (SPC-E) flexible (Berendsen, Grigera, & Straatsma, 1987), modelo TIP3P flexible (Jorgensen et al., 1983) y all-atom GROMOS53A6 (Oostenbrink et al., 2004) fueron utilizados para estudiar las interacciones moleculares entre el agua y n-octano.
Metodología empleada
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:
|
(1) |
|
(2) |
Donde r y r0 corresponden a las longitudes del enlace entre los átomos en el tiempo de simulación y en equilibrio. Los términos θ y θ0 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í, C12ij y C6ij 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: y . Además, ε es el coeficiente de permitividad en el vacío, las variables qi y qj son las cargas definidas para los átomos y rij 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.
Tabla 1. Parámetros no enlazantes, geométricos y electrostáticos de los modelos flexibles de agua usados en las simulaciones.
Parámetros y unidades |
Modelo SPC |
Modelo SPC-E |
Modelo TIP3P |
C600 (kJ mol−1 nm6) |
2.6173×10-3 |
2.6173×10-3 |
2.4897×10-3 |
C1200 (kJ mol−1nm12) |
2.6341×10-6 |
2.6341×10-6 |
2.4351×10-6 |
Enlace O-H (nm) |
0.1000 |
0.1000 |
0.0957 |
KenlaceO-H(kJ nm-4) |
345000 |
345000 |
502416 |
Ángulo H-O-H |
109.47° |
109.47° |
104.52° |
KanguloH-O-H(kJ) |
383.0 |
383.0 |
628.02 |
qH |
-0.8200 |
-0.8476 |
-0.8340 |
q0 |
0.4100 |
0.4238 |
0.4170 |
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).
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×3nm3. Esta construcción fue realizada utilizando las densidades experimentales de 0.997 g/cm3 y 0.699 g/cm3 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×10nm3. 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.
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”.
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.
Para el sistema inmiscible agua/n-octano que posee una región interfacial, la tensión interfacial (γ) fue determinada usando los componentes tangencial (Pt = (Pxx + Pyy)/2) y normal (PN = Pzz) al plano interfacial del tensor presión. En este caso, la tensión interfacial es definida por la ecuación 3:
|
(3) |
Donde, el término Lz 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:
|
(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, rkl representa el vector entre los centros de masa de las moléculas k y l. En la ecuación 4, Fkl 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 (σ0) 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:
|
(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):
|
(6) |
Donde kB es la constante de Boltzmann, Lii 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):
|
(7) |
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 (σ0) 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:
|
(8) |
Donde Z0agua y Z0hidrocarburo 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:
|
(9) |
σ = 1.82 × √2 × σCW |
(10) |
Donde Rg es el radio de giro de la molécula de n-octano obtenido de la simulación.
Resultados y discusión
El software GROMACS-4.5.4 genera los archivos traj.trr y ener.edr, los cuales son necesarios para realizar el análisis de los sistemas estudiados. El archivo ener.edr contiene los componentes de presión usados para determinar la tensión interfacial usando el modelo de Kirwood-Buff y el archivo traj.trr contiene las posiciones y velocidades de los átomos en el tiempo de simulación. Con las posiciones atómicas se obtienen los perfiles de densidad a lo largo del eje z, los cuales son usados para la estimación de los espesores de película interfacial. A continuación, se muestra el análisis de los resultados obtenidos.
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/m3 (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/m3 y 1040 kg/m3 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/m3), 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 -CH2 y -CH3 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).
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/m3) y n-octano (699 kg/m3). 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.
Tabla 1. Densidad de los sistemas simulados y posición del plano divisorio de Gibbs (GDS) a lo largo del eje z.
Capa de agua |
Capa de n-octano |
||||
Sistemas |
Densidad (kg/m3) |
GDS(1) (nm) |
GDS(2) (nm) |
GDS(1) (nm) |
GDS(2) (nm) |
agua(SPC)/n-octano |
823.82±0.17 |
2.81 |
6.00 |
2.75 |
6.07 |
agua(SPC-E)/n-octano |
831.09±0.35 |
2.82 |
5.95 |
2.75 |
6.02 |
agua(TIP3P)/n-octano |
814.48±0.39 |
2.73 |
6.12 |
2.80 |
6.05 |
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.
Tabla 2. Valores de σCW y σo obtenidos para los sistemas simulados, los cuales fueron determinados con la teoría onda-capilar (TCW).
Sistemas |
Capa de agua |
Capa de n-octano |
agua/n-octano |
|||
σCW(1) (nm) |
σCW(2) (nm) |
σCW(1) (nm) |
σCW(2) (nm) |
σo(1) (nm) |
σo(2) (nm) |
|
agua(SPC)/n-octano |
0.133 |
0.133 |
0.131 |
0.131 |
0.06 |
0.07 |
agua(SPC-E)/n-octano |
0.119 |
0.118 |
0.118 |
0.120 |
0.07 |
0.07 |
agua(TIP3P)/n-octano |
0.146 |
0.147 |
0.145 |
0.146 |
0.08 |
0.07 |
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.
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.
Tabla 3. Tensiones interfaciales determinadas con el modelo de Kirwood-Buff (γKB) y la TCW (γcw). La ecuación 6 fue utilizada para estimar la longitud de correlación molecular del hidrocarburo.
Sistemas |
γKB (mN/m) |
ξ (γKB) (nm) |
γcw (mN/m) |
ξ(γcw) (nm) |
agua(SPC)/n-octano agua(SPC-E)/n-octano agua(TIP3P)/n-octano |
55.510.50 64.851.17 48.610.66 |
0.90 0.91 0.78 |
50.60 50.30 50.34 |
0.98 1.27 0.74 |
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.
Conclusiones
La dinámica molecular ha sido usada en el laboratorio de fisicoquímica como herramienta complementaria en la estimación de la tensión interfacial y ancho de la región interfacial del sistema agua/n-octano usando diferentes combinaciones de modelos del agua con el campo de fuerza all-atom GROMOS53A6 usado para el n-octano. Los mejores resultados fueron obtenidos con el modelo TIP3P para el agua. En este caso, la tensión interfacial estimada con la ecuación de Kirwood-Buff fue de 48.610.66 mN/m. A su vez, la teoría Onda-Capilar nos permitió predecir la tensión interfacial de 50.34 mN/m con un valor de una longitud de correlación del hidrocarburo de 0.74 nm. Este valor es muy cercano a la longitud de la molécula de n-octano de forma extendida (0.72 nm). Finalmente, los resultados estimados fueron aceptables al comparar con datos experimentales.
Anexo 1
agua-octano.mdp
define=-DFLEXIBLE
integrator=steep ó md
dt=0.001
nsteps=20000000
nstcomm=1000
nstxout=1000
nstvout=1000
nstfout=1000
nstlog=1000
nstenergy=1000
nstlist=50
energygrps=octanoagua
ns_type=grid
pbc=xyz
vdwtype=cut-off
rlist=1.40
rvdw=1.50
coulombtype=PME
rcoulomb=1.40
pme_order = 6
ewald_rtol=1e-06
optimize_fft = yes
DispCorr=EnerPres
Tcoupl=v-rescale
tau_t=0.1 0.1
tc-grps=octanoagua
ref_t=300 300
Pcoupl=berendsen
pcoupltype=isotropic
tau_p=2
compressibility=4.5e-5
ref_p=0.98
gen_vel=yes
gen_temp=300
Referencias
Allen, M., Tildesley, D., Computer simulation of liquids, Oxford, England: Oxford university press, 1989.
Abraham, M., Van der Spoel, V., Lindahl, E., Hess, B., and the GROMACS development team, GROMACS User Manual version 2019. https://ftp.gromacs.org/pub/manual/manual-5.0.4.pdf.
Berendsen, H., Grigera, J., y Straatsma, T., (1987). The missing term in effective pair potentials, The Journal of Physical Chemistry, 91, 6269-6271. https://doi.org/10.1021/j100308a038.
Berendsen, H., Postma, J., van Gunsteren, W., y Hermans, J., (1981). Intermolecular forces: Interaction models for water in relation to protein hydration, Netherlands: Springer.
Berendsen, H., Postma, J., van Gunsteren, W., DiNola, A., y Haak, J., (1984). Molecular dynamics with coupling to an external bath, Journal of Chemical Physics, 81(8), 3684-3690. https://doi.org/10.1063/1.448118.
Butt, H., Graf, K., y Kappl, M., (2003). Physics and Chemistry of Interfaces, Weinheim, Germany: Wiley-VCH Verlag & Co. KGaA.
Bussi, G., Donadio, D., y Parrinello, M. (2007). Canonical sampling through velocity rescaling. Journal of Chemical Physics, 126(1), 01410. https://doi.org/10.1063/1.2408420.
Castellan, G., Fisicoquímica, México, D.F., México: Editorial Pearson, 1987.
Cordeiro, N., (2003). Interfacial tension behaviour of water/hydrocarbon liquid-liquid interfaces: A molecular dynamics simulation, Molecular Simulation, 29(12), 817-827. https://doi.org/10.1080/0892702031000121905.
Essmann, U., Perera, L., Berkowitz, M., Darden, T., Lee, H., y Pedersen, L., (1995). A smooth particle mesh ewald method, Journal of Chemical Physics, 103(19), 8577-8593. https://doi.org/10.1063/1.470117.
Freitas, A., Quina, F., y Carroll, F., (1997). Estimation of water-organic interfacial tensions: A Linear Free Energy Relationship Analysis of Interfacial Adhesion, Journal of Physical Chemistry B, 101, 7488-7493. https://doi.org/10.1021/jp970927u.
Frenkel, D., y Smit, B. (2002). Understanding molecular simulation: From algorithms to applications, New York: U.S.A.: Academic Press. https://bit.ly/3vFJybY.
Hanwell, M., Curtis, D., Lonie, D., Vandermeersch, T., Zurek, E., & Hutchison, G. (2012). Avogadro: an advanced semantic chemical editor, visualization, and analysis platform, Journal of Cheminformatics, 4(1), 17. https://jcheminf.biomedcentral.com/articles/10.1186/1758-2946-4-17.
Hess, B., Kutzner, C., Van der Spoel, D., y Lindahl, E. (2008). GROMACS 4: algorithms for highly efficient, load-balanced, and scalable molecular simulation, Journal of Chemical Theory and Computation, 4(3), 435-447. https://doi.org/10.1021/ct700301q.
Holmberg, K. (2002). Handbook of applied surface and colloid chemistry (Vol. I), New York, U.S.A.: Van Wiley. https://doi.org/10.1021/ja025281k.
Jorgensen, W., Chandrasekhar, J., Madura, J., Impey, R., y Klein, M. (1983). Comparison of simple potential functions for simulating liquid water, Journal of Chemical Physics, 79, 926-935. https://doi.org/10.1063/1.445869.
Kirkwood, J., y Buff, F. (1949). The statistical mechanical theory of surface tension, Journal of Chemical Physics, 17(3), 338-343. https://doi.org/10.1063/1.1747248.
Lewars, E. (2016). Computational chemistry: Introduction to the theory and applications of molecular and quantum mechanics, Switzerland: Springer.
Lyklema, J. (2000). Liquid-fluid interfaces, volume III, San Diego, U.S.A.: Academic Press.
Malde, A., Zuo, L., Breeze, M., Stroet, M., Poger, D., Nair, P., y Mark, A. (2011). An automated force field topology builder (ATB) and repository: version 1.0, Journal of Chemical Theory and Computation,7(12), 4026-4037. https://doi.org/10.1021/ct200196m.
Mitrinovic, D., Tikhonov, A., Li, M., Huang, Z., y Schlossman, M. (2000). Noncapillary-wave structure at the water-alkane interface, Physical Review Letters, 85(3),582-585.
Mitrinovic, D., Zhang, Z., Williams, S., Huang, Z., y Schlossman, M. (1999). X-ray Reflectivity Study of the Water-Hexane Interface, The Journal of Physical Chemistry B, 103(11), 1779-1783. https://doi.org/10.1021/jp984640o.
Moore, W., Physical chemistry, London, England: Longman Group Limited, 1972.
Neyt, J., Wender, A., Lachet, V., Ghoufi, A., y Malfreyt, P. (2014). Quantitative Predictions of the Interfacial Tensions of Liquid-Liquid Interfaces through Atomistic and Coarse-Grained Models, Journal of Chemical Theory and Computation, 10, 1887-1899. https://doi.org/10.1021/ct500053c.
Nicolas, J., de Souza, N. (2004). Molecular dynamics study of the n-hexane/water interface: Towards a better understanding of the liquid-liquid interfacial broadening, Journal of Chemical Physics, 120, 2464-2469. http://dx.doi.org/10.1063/1.1629278.
Oostenbrink, C., Villa, A., Mark, A., & Van Gunsteren, W. (2004). A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force field parameter sets 53A5 and 53A6, Journal of Computational Chemistry, 25(13), 1656-1676. https://doi.org/10.1002/jcc.20090.
Riedleder, A., Kentish, S., Perera, J., y Stevens, G. (2007). Structural investigation of a water/n-heptane Interface: A Molecular Dynamics Study, Solvent Extraction and Ion Exchange, 25(1), 41-52. https://doi.org/10.1080/07366290601067424.
Rivera, L., McCabe, C., Cummings, P. (2003). Molecular Simulations of Liquid-Liquid Interfacial properties: Water-n-alkane and water-methanol-n-alkane systems, Physical Review E, 67, 011603. https://doi.org/10.1103/PhysRevE.67.011603.
Rowlinson, J., y Widom, B. (2002). Molecular theory of capillarity, NewYork, U.S.A.: Dover publications.
Senapati, S., y Berkowitz, M. (2001). Computer simulation study of the interface width of the liquid-liquid interface, Physical Review Letters, 87(17), 176101-1-4. https://doi.org/10.1103/PhysRevLett.87.176101.
Underwood, T., y Greenwell, C. (2018). The water-alkane interfaces at various NaCl salt concentrations: A molecular dynamics study of the readily available force fields. Scientific Reports, 8, 352. https://www.nature.com/articles/s41598-017-18633-y.
van der Spoel, D., Lindahl, E., Hess, B., van Buuren, E., Apol, E.,Meulenhoff, P.,Tieleman, D.,Sijbers, A.,Feenstra, K., van Drunen, R.,yBerendsen, H. (2010). Gromacs User Manual version 4.5.4, www.gromacs.org.
Wick, C., Chang, T., Slocum, J. y Cummings, O. (2012). Computational Investigation of the n-alkane/Water Interface with Many-Body Potentials: The Effect of Chain Length and Ion Distributions, Journal of Physical Chemistry C, 116, 783-790. https://doi.org/10.1021/jp208459g.
Yaws, C. (2003). Handbook of Thermodynamic and Physical Properties of Chemical Compounds, Knovel.
Zhang, Y., Feller, S., Brooks, B., y Pastor, R. (1995). Computer simulation of liquid/liquid interfaces. I. Theory and application to octane/water, Journal of Chemical Physics, 103, 10252-10266. https://doi.org/10.1063/1.469927.
Recepción: 16-06-2020
Aceptación: 22-03-2021
Cómo citar
Parra Figueredo, J. G., Iza, P. y Perozo, E. (2021, julio-septiembre). Una guía para la estimación de la tensión interfacial y el espesor de la interface de un sistema agua/hidrocarburo usando el programa gromacs-4.5.4. Educación Química, 32(3). http://dx.doi.org/10.22201/fq.18708404e.2021.3.76027.
[a] Universidad de Carabobo, Facultad de Ciencias y Tecnología, Dpto. de Química, Lab. de Fisicoquímica, Lab. de Química Computacional (QUIMICOMP), Bárbula-Venezuela. Contacto: josegregorioparra@gmail.com.
[b] Escuela Superior Politécnica del Litoral, ESPOL, Departamento de Física, Campus Gustavo Galindo, Guayaquil, Ecuador.
[c] Universidad de Carabobo, Facultad de Ciencias y Tecnología, Dpto. de Química, Dirección del Dpto. de Química, Bárbula-Venezuela.