Método simple para estimar parámetros del potencial entre pares: propiedades presión-volumen-temperatura del benceno mediante simulación molecular Monte Carlo

A. Elías-Domínguez* iD, M. Torres-Suárez, F. Pérez-Villaseñor, A. Castro-Agüero, A. Ortiz-Arroyo, M.V. Hernández-Ruíz, E. Reyes-Pérez, M.R. Castillo-Estrada

Universidad Autónoma de Tlaxcala, Facultad de Ciencias Básicas, Ingeniería y Tecnología, Ingeniería Química. Calzada Apizaquito s/n. Apizaco, Tlaxcala, México, C.P. 90401. *aeliasdom@yahoo.com.mx

http://doi.org/10.5281/zenodo.5091165

Bajar cita (RIS): Elías-Domínguez y cols., 2019 AyTBUAP 4(15): 1-15

Editado por: Jesús Muñoz-Rojas (Instituto de Ciencias, Benemérita Universidad Autónoma de Puebla)

Fecha de Publicación: 3 de octubre de 2019

EOI: https://eoi.citefactor.org/10.11235/BUAP.04.15.01

URI: https://hdl.handle.net/20.500.12371/9333

Referencia: Elías-Domínguez A, Torres-Suárez M, Pérez-Villaseñor F, Castro-Agüero A, Ortiz-Arroyo A, Hernández-Ruíz M V, et al. Método simple para estimar parámetros del potencial entre pares: propiedades presión-volumen-temperatura del benceno mediante simulación molecular Monte Carlo. Alianzas y Tendencias BUAP [Internet]. 2019;4(15):1–15. Available from: https://www.aytbuap.mx/publicaciones#h.l2n12mg4jp24

RESUMEN

En este trabajo se determinó la presión, con la ecuación del virial molecular, de benceno puro mediante la técnica de simulación molecular Monte Cario en el ensamble canónico, para un rango de temperatura de 513.15 a 628.15 kelvin considerando densidad constante de 0.9144 mol/L y densidad variable de 0.8829 mol/L a 2.1336 mol/L, empleando un modelo molecular de átomos unidos con interacciones exponencial-6. Las estimaciones obtenidas por simulación molecular mostraron desviaciones respecto a los datos experimentales, presentando errores promedio de 9.12% a densidad constante y de 16.70% considerando densidad variable. Con el fin de mejorar la predicción, se propone un método simple para estimar los parámetros del potencial exponencial-6, basado en extrapolar el valor del error de predicción a cero, a partir de un ajuste lineal del error. Los parámetros así obtenidos, se utilizaron para calcular nuevamente la presión del benceno puro a las mismas condiciones de temperatura y densidad antes mencionadas, con errores promedios de predicción de 0.68% a densidad constante y de 2.29% considerando densidad variable.

Palabras clave: Estimación de parámetros; simulación Monte Carlo; benceno.

ABSTRACT

Pressure of puré benzene is calculated using a molecular Monte Cario simulation technique in the canonical ensemble with the molecular virial equation, data was calculated in a temperatura span of 513.15 to 628.15 Kelvin, considering two cases. In the first instance, constant density of 0.9144 mol/L is used and in the second case density is spanned from 0.8829 mol/L to 2.1336 mol/L, wlth valúes calculated using the molecular model of unlted atoms wlth exponential-6 interactions. The estimations obtained by molecular simulatlon show average errors of 9.12% from the experimental data for constant density and 16.70% for variable density case. To improve data prediction, a simple method to estímate the parameters of the exponential-6 potential is proposed, the estimation method is based on the extrapolation of the error prediction valué to zero, from a linear adjustment of the error. The parameters obtained were used to recalculate the pressure of puré benzene at the same temperature and density conditions, predicted data showed an average error of 0.68% for constant density and 2.29% for variable density.

Keywords: Parameter estimation; Monte Carlo simulation; benzene.

INTRODUCCIÓN

La termodinámica deriva relaciones matemáticas que conectan diferentes propiedades experimentales de los sistemas macroscópicos en sistemas de equilibrio que contienen muchas moléculas, del orden de 1O²⁰ o más. Estas interconexiones de la termodinámica no dan ninguna información sobre la interpretación o explicación, a nivel molecular, de las propiedades experimentales observadas. Tal explicación cae más bien dentro del campo de la mecánica estadística o la termodinámica estadística [1].

La mecánica estadística es un formalismo que tiene como objetivo explicar las propiedades físicas de la materia del sistema macroscópico sobre la base del comportamiento de sus componentes microscópicos. Por lo tanto, los campos cubiertos por la mecánica estadística y la termodinámica clásica coinciden. El alcance del formalismo es casi tan ilimitado como el rango mismo de los fenómenos naturales, ya que es aplicable a la materia en cualquier estado. De hecho, se ha aplicado, con considerable éxito, al estudio de la materia en estado sólido, estado líquido o estado gaseoso, materia compuesta de varias fases y/o varios componentes, materia bajo condiciones extremas de densidad y temperatura, etc. Además, el formalismo de la mecánica estadística permite investigar los estados de no equilibrio de la materia así como los estados de equilibrio [2].

Las relaciones presión-volumen-temperatura (PVT) tienen gran utilidad en el diseño y mejoramiento de procesos industriales, como la determinación del tamaño de recipientes y tuberías, la medición de fluidos, etc. Además sirven como parámetros para evaluar propiedades termodinámicas, como la entalpia y la energía interna, las cuales controlan los procesos de conversión de energía. La predicción de estas relaciones mediante simulación molecular resulta ser conveniente para procesos que requieren información que involucran presiones o temperaturas muy altas, que de otra manera seria muy complejo medir experimentalmente [3]. La simulación molecular proporciona resultados exactos para un modelo molecular, sin tener que recurrir a teorías aproximadas y proporciona información de las propiedades a nivel molecular, que no se podrían predecir mediante la termodinámica clásica. El método de simulación Monte Cario es muy apropiado para determinar las propiedades PVT de sustancias puras, el cual requiere del modelo molecular de la sustancia, que incluya su geometría y los parámetros de los potenciales de interacción [3]. Por lo tanto, el modelo molecular es un elemento importante en la simulación molecular, pues de la correcta elección de la función de potencial y de los parámetros del potencial, depende la obtención de buenos resultados.

En la literatura [4, 5, 6, 7, 8, 9, 10, 11, 12] se han presentado y utilizado un sinnúmero de modelos moleculares de diversas sustancias. Estos modelos consideran a las sustancia como moléculas rígidas lineales (dióxido de carbono, etano, etc.) y no lineales (agua, benceno, ciclo hexano, etc.), así como también moléculas flexibles (alcanos y polímeros), etc. La forma de obtener los parámetros que definen la función de potencial de interacción, es mediante prueba y error, es decir, se establecen valores de los parámetros y se mide su predicción con respecto a uno o dos valores experimentales de las propiedades de interés [4, 7, 13]. Contreras- Camacho et al. en 2004 [14] propusieron un método para determinar los parámetros del potencial de interacción, basado en fijar cuatro o cinco valores de las propiedades de interés y minimizar los errores al cuadrado mediante la expansión de cada propiedad considerada por una serie de Taylor de primer orden, esto conlleva a evaluar derivadas parciales de la propiedad con respecto a los parámetros a optimizar y resolver un sistema de ecuaciones lineales; la desventaja de este método, es la evaluación de las derivadas y el resolver el sistema de ecuaciones resultantes, que cada vez que se aumente el número de datos experimentales de las propiedades de interés a considerar para la estimación de los parámetros que definen la función potencial, será también más complejo. Con el fin de disponer de un método de complejidad intermedia, en este trabajo se propone un método simple para estimar los parámetros o (representa el tamaño de los sitios de interacción de las moléculas que forman al sistema) y e (representa la profundidad del pozo de la función potencial del modelo molecular en estudio) del potencial exponencial-6, basado en extrapolar el valor del error de predicción a cero, a partir de un ajuste lineal del error. El método propuesto se aplica al modelo molecular del benceno puro, sustancia utilizada a nivel industrial como solvente, como aditivo y principalmente para la producción de etilbenceno, eumeno y ciclo hexano [15].

Metodología

Un ensamble es una colección de un gran número de sistemas, cada uno construido para ser una réplica del sistema macroscópico estudiado.

Aunque todos los sistemas en el ensamble son idénticos desde un punto de vista macroscópico, no son todos idénticos en el nivel molecular. De hecho, en general, hay un número extremadamente grande de estados cuánticos consistentes con un estado termodinámico dado [1]. Los métodos de simulación Monte Cario se basan en generar una serie de configuraciones en forma aleatoria, mediante ciertas reglas de aceptación, en aquellas regiones del espado fase que tengan una contribución importante en el promedio del ensamble. Las reglas de aceptación se seleccionan tal que estas configuraciones ocurran con una frecuencia descrita mediante una distribución de probabilidad dada por el ensamble. Existen diferentes tipos de ensambles: el ensamble canónico, ensamble isotérmico-isobárico, ensamble grand canónico, ensamble microcanónico, etc. En la literatura [3, 16, 17] existe una descripción completa sobre los métodos de simulación molecular.

En el presente trabajo, la técnica de simulación Monte Cario en el ensamble canónico se utilizó para el cálculo de la presión del benceno, mediante la ecuación del virial molecular, a diferentes temperaturas y densidades.

En el ensamble canónico, el número de partículas (N), la temperatura (T) y el volumen (V) se fijan. Este también se conoce como ensamble (NVT). La técnica que se aplica al ensamble canónico se realiza en una caja o región de simulación, tal como se muestra en la Figura 1, con la finalidad de hacer un muestreo de la distribución de probabilidad de este ensamble [3, 16]. En la caja o región de simulación, se establecen los parámetros del ensamble de interés, el modelo molecular y las condiciones iniciales de las moléculas, como posiciones y orientaciones (inicialización). Las simulaciones incluyen un número de jugadas iniciales, que llevan el sistema a un estado estacionario en el valor de las propiedades (equilibración), luego de las cuales se realiza un conjunto de jugadas de producción, con el fin de promediar las propiedades de interés. El algoritmo [3, 16] incluye los siguientes pasos:

- Seleccionar aleatoriamente una partícula y calcular la energía potencial de esta configuración U(o).

- Una vez seleccionada la partícula, se efectúa un desplazamiento aleatorio, tal como se muestra en la Figura 3.1, en el intervalo (-Δ, Δ). La nueva configuración se denota por n y su energía potencial por U(n).

- El movimiento se acepta de acuerdo con el criterio de la Ecuación (1):

si el movimiento se rechaza, se mantiene la configuración anterior.

La Ecuación (1) representa la distribución de probabilidad del ensamble canónico [3,16], donde acc(o—>«) denota aceptar o rechazar la configuración nueva (n) y β = (k_BT)^-1 con k_B la constante de Boltzmann equivalente a 3.29946- 10^-27 Kcal/K molec.

Modelos de potenciales moleculares

La predicción del comportamiento macroscópico de los fluidos a partir del comportamiento de sus componentes moleculares es proporcionada por la mecánica estadística y la simulación por computadora. La base de tal predicción es la disponibilidad de modelos moleculares. Estos modelos deben contener información sobre la posición de los átomos y los parámetros de los potenciales de interacción entre ellos. Los modelos moleculares se formulan en términos de la energía de un sistema de moléculas. Esta energía molecular depende de las estructuras geométricas de las moléculas y del campo de fuerza en el que se mueven. La relación entre la estructura geométrica de un cuerpo y su energía se define en términos mecánicos. El campo de fuerza es el resultado de las propiedades eléctricas de las moléculas. En este nivel, los modelos se basan en mecánica clásica y electrostática [18].

En este trabajo se utiliza el modelo molecular del benceno propuesto por Errington y Panagiotopoulos [7], el cual tiene seis sitios de interacción y se considera como una molécula rígida y planar, con la distancia de enlace de 1.40 Å, el ángulo de enlace θ de 120° y los ángulos de diedro 0 iguales a cero debido a que el modelo es hexagonal planar. Con esta geometría las coordenadas cartesianas (x,y,z) de los átomos de la molécula de benceno con respecto al cuerpo se obtienen definiendo uno de los sitios de interacción como origen, asignándole las coordenadas r_t = [0,0,0] y usando trigonometría plana se establecen los 5 sitios de interacción restantes (Tabla 1) [7].

Para describir las interacciones entre sitios de benceno, y calcular la energía potencial del sistema Uij (r) se usa la función de potencial exponencial-6 de Buckingham dado por la Ecuación (2):

con r_m= 1.09623 σ_ij^ab, α =20, Ɛ_ij^ab/K_B =74.06 Kelvin(K), σ_ij^ab =3.71Å y r_ij^ab es la distancia; los superíndices ab denotan la interacción entre los sitios a y b localizados en las molécula i y j,. Al potencial se le aplica un truncamiento a la distancia r_cut=12 Å, por lo que es necesario añadir a la u_ij(r) un término de corrección de largo alcance U_corr de acuerdo a Frenkel y Smit [3].

La presión (P) se calcula mediante la Ecuación (3):

Donde p es la densidad del sistema, V es el volumen de la región de simulación, w(r) es el virial molecular y se obtiene mediante la Ecuación (4):

Utilizando el potencial de pares Exp-6 de Buckingham, w(r_ij) queda expresada mediante la Ecuación (5):

Para generar una configuración inicial de los sistemas moleculares y de los movimientos siguientes, es de suma importancia describir los procedimientos para asignar y cambiar la orientación de las moléculas, que para el caso de estudio corresponden a las moléculas del benceno. La orientación de un cuerpo rígido especifica la relación de un sistema de ejes fijo en el espacio y otro sistema fijo con respecto al cuerpo. Cualquier vector unitario e se puede expresar en términos de los componentes con respecto al cuerpo o con respecto al espacio: se utilizara la notación e^b y e^s respectivamente. Estos vectores están relacionados mediante la matriz de rotación A, Ecuación (6).

donde la operación indicada es un producto punto o escalar.

La matriz de rotación describe completamente la orientación relativa de los dos sistemas, es decir, define la orientación molecular. Como e es un vector fijo en el marco molecular, entonces e^b queda definido por la geometría del modelo molecular (cuerpo rígido) y no cambia con el tiempo. Sin embargo, en las coordenadas del espacio fijo, los componentes de e^s variarán. La matriz de rotación A contiene tres cantidades (ángulos) independientes. Por conveniencia se utilizan los ángulos de Euler ϕθψ para describir la orientación de una molécula no lineal [16] y se generan aleatoriamente. La matriz de rotación A en términos de los ángulos de Euler está dada por la Ecuación (7)

La Ecuación (7) se utiliza para describir la rotación de las moléculas del benceno.

Detalles de la simulación

Se realizaron simulaciones moleculares mediante el método Monte Cario en el ensamble canónico NVT para calcular la presión de benceno puro en un rango de temperatura de 513.15 K a 628.15 K, considerando densidad constante de 0.9144 mol/L, y considerando densidad variable de 0.8829 mol/L a 2.1336 mol/L. Se utilizaron propiedades reducidas. Las simulaciones se llevaron a cabo utilizando el modelo molecular propuesto por Errington y Panagiotopoulos [7] siguiendo el procedimiento descrito en la metodología; se inició definiendo una caja de simulación con 125 moléculas de benceno, con un truncamiento esférico al potencial a la distancia de 12 Å y se establecieron condiciones de frontera periódicas y de mínima imagen [3,16]. La configuración de arranque se estableció asignando a las moléculas en un enrejado cúbico simple y definiendo su orientación de manera aleatoria mediante ángulos de Euler [19]. Las simulaciones se realizaron en ciclos, cada uno con 100 movimientos (desplazamiento y rotación simultánea de las moléculas). Para equilibrar el sistema se requirieron de 10 millones de configuraciones, seguidas de 10 millones de configuraciones para promediar las propiedades.

Todos los parámetros utilizados en las simulaciones se deben expresar en unidades reducidas, para facilitar los cálculos durante las simulaciones. Para ello se eligen una unidad de energía (Ɛ), longitud (σ) y masa (m) y se expresan todas las cantidades en términos de estas unidades. La densidad (ρ), temperatura (T), energía (U) y presión (P) se expresan en unidades reducidas (denotadas por el superíndice *) como:

Método de estimación de parámetros del potencial de interacción

Para la estimación de los parámetros, se requieren valores iniciales de σ y Ɛ/k_B del potencial de interacción del benceno, con estos valores iniciales se realizan cálculos de la presión, a cierta temperatura y volumen fijas, mediante el método de simulación Monte Cario y evaluar el error de predicción, mediante la Ecuación (12).

Después variar el valor de σ (ó Ɛ/k_B) hacia arriba y hacia abajo, y evaluar la presión mediante simulación manteniendo constante el valor de Ɛ/k_B (ó σ) y su %error; en una gráfica de %error en función de o (ó Ɛ/k_B) representar los tres resultados realizados, y observar el comportamiento, si es lineal, obtener el ajuste por mínimos cuadrados. Para obtener la ecuación lineal, por ejemplo, si se varía Ɛ/k_B con o constante, entonces se obtendrá por linearización la Ecuación (13).

Para estimar un mejor valor del parámetro Ɛ/k_B que minimice el %error, sustituir el %error = 0 en la Ecuación (13) y despejar Ɛ/k_B que será Ɛ/k_B^opt. Ahora variar σ y manteniendo el valor de Ɛ/k_B^opt, evaluar la presión mediante simulación molecular y su %error, y en una gráfica de %error en función de σ representar los tres resultados realizados, y observar el comportamiento, si es lineal, obtener el ajuste por mínimos cuadrados, para obtener la Ecuación (14).

Con esta ecuación se puede calcular el valor σ^0pt haciendo %error = 0.

En el caso de no observar un comportamiento lineal de los %error, se deberá eliminar un punto para lograrlo; esto para extrapolar al valor cero en el %error, pues puede darse el caso que con los tres puntos representados en el gráfico, nunca se extrapole al cero en el %error, debido a un comportamiento creciente-decreciente o viceversa.

Resultados y Discusión

Se realizaron simulaciones en el ensamble NVT a 513.15 K y 0.9144 mol/L del benceno puro, con valores iniciales de los parámetros σ = 3.710Å y Ɛ/k_B = 74.06 K correspondientes a los reportados por Errington y Panagiotopoulos [7]. Para aplicar el método propuesto y estimar los parámetros σ y s/k_B que mejoren la predicción de las propiedades, en este caso la presión obtenida mediante simulación molecular (P^sim), se propone modificar el valor del parámetro Ɛ/k_B en 2 unidades hada arriba y hada abajo del valor ¡nidal, en la Tabla 2 se presentan los valores de los parámetros y los %error en la estimación de la presión mediante la Ecuación (12), donde el valor de la presión experimental (P^exp) se tomó de la literatura [20].

En la Figura 2, se muestra al %error como una función del parámetro Ɛ/k_B, observando que el valor de Ɛ/k_B determina el valor de la P^sim; los datos graficados se ajustan a un modelo lineal, Ecuación (15).

con coeficiente de correlación de 0.9993; mediante extrapolación a un valor de %error = 0 se obtiene el valor de Ɛ/k_B^opt = 82.16 K. De manera similar, se obtiene el valor para σ^opt. En la Tabla 3 se presentan los valores modificados para el parámetro σ, P^sim y el %error con Ɛ/k_B°^pt = 82.16 K.

los resultados (Figura 3) se ajustaron a una línea recta con coeficiente de correlación de 1.0, obteniendo la Ecuación (16).

y haciendo %error = 0, se obtiene σ^opt= 3.76 Å. Notar que en la Figura 3, se representan dos puntos, los correspondientes a a de 3.51 y 3.71 Å, eliminando el tercer punto correspondiente a σ = 3.91 A una vez que hace perder la linealidad.

Para probar la predicción del modelo molecular del benceno con Ɛ/k_B°^pt = 82.16 K y σ^opt= 3.76 Å, se realizaron simulaciones moleculares mediante el método Monte Cario en el ensamble canónico NVT para calcular la presión de benceno puro, en un rango de temperatura de 513.15 K a 628.15 K, considerando densidad constante de 0.9144 mol/L, y considerando densidad variable de 0.8829 mol/L a 2.1336 mol/L. Para efectos de comparación, también se realizaron simulaciones con los valores iniciales de Ɛ/k_B y σ.

Los resultados de la presión de benceno, tanto experimentales [20] como los obtenidos mediante simulación molecular con los valores iniciales de Ɛ/k_B = 74.06 K y σ = 3.710Å, se muestran en la Figura 4 y en la Tabla 4. Para un rango de temperatura de 513.15 K a 628.15 K considerando densidad constante de 0.9144 mol/L (Figura 4) el error promedio respecto a los datos experimentales fue de 9.12% y considerando densidad variable de 0.8829 mol/L a 2.1336 mol/L (Tabla 4) el error promedio fue de 16.70%. Los errores de predicción tienden a aumentar cuando las condiciones de temperatura y densidad se incrementan.

Con Ɛ/k_B^opt=82.16 K y σ^opt=3.76 Å se realizaron nuevamente las simulaciones moleculares para el benceno. En la Figura 5 se comparan los resultados de la presión obtenida por simulación con los valores experimentales, para un rango de temperatura de 513.15 K a 628.15 K con densidad constante de 0.9144 mol/L, presentando un error promedio de 0.68% y en la Tabla 5, considerando densidad variable de 0.8829 mol/L a 2.3024 mol/L, presentando un error promedio de 2.29%. Comparado la predicción de la presión del benceno mediante simulación molecular, con valores iniciales y optimizados de Ɛ/k_B y σ, se logró reducir el error promedio en un 92.47% considerando densidad constante y en un 86.28% considerando densidad variable.

Los valores de Ɛ/k_B^op=82.16 K y σ^0pt=3.76 Å para el modelo molecular del benceno obtenidos en este trabajo, se comparan con los reportados por otros autores [7, 14, 21, 22, 23, 24, 25, 26] donde los valores de Ɛ/k_B se encuentran entre 48 K a 89.415 K y los de ct entre 3.24 Å y 3.75 Å. El parámetro Ɛ/k_B representa la profundidad del pozo referida a k_B de la función potencial dada por la ecuación (2) y el Ɛ/k_B^opt cae dentro del rango de valores observados, mientras que, el parámetro a representa el tamaño de los sitios de interacción CH de la molécula del benceno, el valor de σ^opt=3.76 Å sugiere que la molécula del benceno sea apenas ligeramente mayor al tamaño máximo reportado en la literatura [7, 14, 21, 22, 23, 24, 25, 26]. Finalmente, con los valores de Ɛ/k_B^opt=82.16 K y σ^opt=3.76 Å se predice de mejor manera el comportamiento presión-volumen-temperatura del benceno en el rango de temperaturas estudiado.

CONCLUSIONES

El método propuesto para la estimación de los parámetros σ y Ɛ/k_B del potencial exponencial-6, para el modelo molecular del benceno puro, fue satisfactorio, pues comparando la predicción de la presión del benceno mediante simulación molecular, con valores iniciales y optimizados de Ɛ/k_B y σ, se logró reducir el error promedio hasta en un 92.47%. El método es simple, no requiere de evaluaciones de derivadas y de la solución de sistema de ecuaciones, pues se basa en extrapolar el valor del error de predicción a cero, a partir de un ajuste lineal del error.

CONFLICTO DE INTERESES

Los autores declaran no tener conflictos de intereses.

REFERENCIAS

[1]. HUI TL. An introduction to statistical thermodynamics. Dover Publications: New York, 1986.

[2]. Pathria RK., Paul DB. Statistical mechanics. Third edition. Academic Press: Oxford, 2011.

[3]. Frenkel D., Smit B. Understanding Molecular Simulation. Academic Press: San Diego, California 2002.

[4]. Jorgensen WL., Chandrasekhar J., Madura JD., Impey RW., Klein ML. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983; 79(2): 926-935.

[5]. Harris JG., Yung KH. Carbón dioxide's liquid- vapor coexistence curve and critical properties as predicted by a simple molecular model. J. Phys. Chem. 1995; 99(31): 12021-12024.

[6]. Smit B., Karaborni S., Siepmann JI. Computer simulation of vapor-liquid phase equilibria of n- alkanes. J. Chem. Phys. 1995; 102(5): 2126- 2140.

[7]. Errington JR., Panagiotopoulos AZ. New intermolecular potential models for benzene and cyclohexane. J. Chem. Phys., 1999; 111(21): 9731-9738.

[8]. Martin MG., Siepmann JI. Novel configurational- bias Monte Cario method for branched molecules. Transferable potentials for phase equilibria. 2. United-atom description of branched alkanes. J. Phys. Chem. B. 1999; 103(21): 4508-4517.

[9]. Jindal KS., Maginn El Monte Cario Simulations of gas solubility ¡n the ionic liquid l-n-butyl-3- methylimidazolium hexafluorophosphate. J. Phys. Chem. B. 2005; 109(2): 10395-10404.

[10]. Singh JK., Errington JR. Calculation of phase coexistence properties and surface tensión of n- alkanes with grand-canonical transition-matrix Monte Cario simulation and finite-size scaling. J. Phys. Chem. B. 2006; 110(3): 1369-1376.

[11]. Chong EL., Myshakin EM. Molecular simulation of competitive adsorption of carbón dioxide- methane mixture on illitic day surfaces. Fluid Phase Equilibria. 2018; 472: 185-195.

[12]. Spyriouni T., Economou IG., Theodorou DN. Molecular simulation of the puré n-hexadecane vapor-liquid equilibria at elevated temperatura. Macromolecules 1998; 31(4): 1430-1431.

[13]. Jorgensen WL., Madura JD., Swenson CJ. Optimizad intermolecular potential function for liquid hydrocarbons. J. Am. Chem. Soc. 1984; 106(22): 6638-6646.

[14]. Contreras-Camacho O., Ungerer P., Boutin A., Mackie AD. Optimizad intermolecular potential for aromatic hydrocarbons based on anisotropic united atoms. 1. Benzene. 1 Phys. Chem. B. 2004; 108(37): 14109-14114.

[15]. Occupational Safety 8i Health Administration (OSHA), Occupational Health Guideline for Benzene, 1988. Available at: www.cdc.gov/niosh/pdfs/0049.pdf. (Último acceso: 15 Agosto 2018).

[16]. Alien MP., Tildesley DJ. Computer simulation of liquids. ed. Oxford University Press Inc.: New York 1987.

[17]. Rapaport DC. The art of molecular dynamics simulation. Second edition. Cambridge University press: New York 1995.

[18]. Lucas K. Molecular models for fluid. First edition. Cambridge University Press: 2007.

[19]. Essén H., Svensson M. Calculation of coordinates from molecular geometric parameters and the concept of a geometric calculator. Computers Chem. 1996; 20(4): 389- 395.

[20]. Goodwin R. Benzene termophysical properties from 279 to 900 K at pressures up to 1000 bar, Journal of Physical and Chemical Reference Data 1988; 17(4): 1541-1636.

[21]. Evans DJ., Watts RO. Interactions between benzene molecules, Molecular Physics, 1976; 31(1): 83-96.

[22]. Friedrich A., Lustig RJ. Thermodynamics of fluid benzene from molecular dynamics simulations, Mol. Liq. 2002; 98-99, 243-261.

[23]. Jorgensen WL., Severance DL. Aromatic- aromatic interactions: free energy profiles for the benzene dimer in water, chloroform, and liquid benzene, Journal of the American Chemical Society 1990; 112(12): 4768-4774.

[24]. Linse, P. Thermodynamic and structural aspects of liquid and solid benzene. Monte Cario study. Journal of the American Chemical Society 1984; 106(19): 5425-5430.

[25]. Claessens M., Ferrarlo M., Ryckaert JP. The structure of liquid benzene. Molecular Physics 1983; 50(1): 217-227.

[26]. Wick CD., Martin MG., Siepmann JI. Transferable Potentials for Phase Equilibria. United-Atom Description of Linear and Branched Alkenes and Alkylbenzenes. The Journal of Physical Chemistry B 2000; 104(33): 8008-8016.