- Email: [email protected]
water interfaces
Fuel 267 (2020) 117252
Contents lists available at ScienceDirect
Fuel journal homepage: www.elsevier.com/locate/fuel
Full Length Article
Molecular dynamics investigation on n-alkane-air/water interfaces a
b
a,⁎
George Rucker , Xiong Yu , Liqun Zhang a b
T
Department of Chemical Engineering, Tennessee Technological University, Cookeville, TN, 38505, United States Department of Civil and Environmental Engineering, Case Western Reserve University, Cleveland, OH, 44106, United States
ARTICLE INFO
ABSTRACT
Keywords: Surface tension Alkane Interface Diffusion coefficient Order parameters Molecular dynamics simulation
Alkanes are the important molecules in crude oil that have applications in different fields. In this project, molecular dynamics simulations have been applied to investigate the n-alkane-air and n-alkane-water interfaces. The surface tension, diffusion coefficients, order parameters, radial distribution function, density profiles and radius of gyration of alkanes in the bulk and on the interfaces with air/water were calculated and compared. It was found that with the n-alkane chain length increasing, the surface tension of the interface system increases. However, with the temperature increasing, the surface tension of the alkanes decreases. The alkanes on the interface can diffuse faster than in the bulk in both kinds of interface systems. While n-alkane diffusion rate in the bulk of interface systems is faster than in the pure alkane systems, the radius of gyration of alkanes in the bulk, on the interface and in the pure alkane systems is consistent. Decomposing the diffusion coefficients into the lateral direction and z-direction, the lateral diffusion rate of alkanes on the interface is about 2–10 times faster than in the z-direction in both kinds of interface systems. Although molecules prefer to pack randomly in the bulk, they are more orderly packed on the interface. The result can help to explain the interface effect in water/air-alkane interface systems. The project also demonstrates that using NAMD program and CHARMM forcfield can predict reasonable structure and properties of n-alkane-air/water systems.
1. Introduction Alkanes are an important class of molecules that have various applications in different fields. They are natural nonpolar solvents, and are the major components in crude oil [1]. The structure and dynamics of the interface of alkane with air and with water are of great interest to researchers [2–4] and of great importance to its application in chemistry, physics, engineering and biology, such as the application of water flooding and surfactant flooding in enhanced crude oil recovery [5]. Besides experimental methods, molecular dynamics simulation methods have been applied to address the liquid/air and liquid/liquid interface properties, and satisfying agreements between simulation prediction and experimental data were reached [6–11]. Carpenter and Hehre [12] studied the n-hexane/water interface system using molecular dynamics simulation method, and found differences between the interfacial and bulk liquids. Zhang et al [4] did CHARMM all-atom molecular dynamics simulations on neat n-octane/water and n-octane/ air interfaces. They predicted the surface tension of both kinds of interfaces using different simulation ensembles. Consistent agreement was reached from different ensembles and the predictions were also consistent with experimental work. Van Buuren et al. [13] worked on the decane/water interface, and found that decane molecules on the ⁎
water interface preferred to orient laterally. Rivera et al. [14] worked on the n-alkane/water interfaces and found that with the alkane chain length increasing, the alkane-water surface tension increases. Nicolas and Smit computed the surfaced tension of linear alkanes including nhexane, n-decane and n-hexadecane using united atom OPLS forcefield and SKS forcefield [15]. They found that both models can predict the surface tension of three alkanes consistently agreeing with experimental work at high temperatures. In order to predict the surface tension correctly, it is critical that the model can predict the density correctly. However, up to now, no research on the alkane-air and alkanewater interfaces has been done in a wide carbon number range and in a wide temperature range using CHARMM forcefield, and systematically comparing results from both kinds of systems. Thus, this project studies alkanes from n-C8 to n-C24 in both water and air interfaces, and in the temperature range of 300–400 K. Generally, there are four different kinds of liquid/air, liquid/liquid and liquid–liquid–air interfaces as the sketches shown in Fig. 1. In order to calculate the surface tension of liquid–air interfaces, it is straightforward to set up the interface system as Fig. 1(a), in which the air is represented by vacuum in simulations. In order to calculate the surface tension of liquid–liquid interfaces, it is straightforward to set up the interface system as Fig. 1(b) in which the same amount of liquid I is
Corresponding author at: 1020 Stadium Dr., Cookeville, TN 38505, United States. E-mail address: [email protected] (L. Zhang).
https://doi.org/10.1016/j.fuel.2020.117252 Received 16 November 2019; Received in revised form 27 January 2020; Accepted 28 January 2020 Available online 08 February 2020 0016-2361/ © 2020 Elsevier Ltd. All rights reserved.
Fuel 267 (2020) 117252
G. Rucker, et al.
Fig. 1. Sketches of the typical liquid–liquid/vapor interface systems: Liquid-vacuum interface (a); liquid–liquid interface (b); and two kinds of liquid–liquid-vapor interfaces (c) and (d). The vacuum has been applied to represent the air in simulations.
above and below liquid II in the simulation box. The amount of liquid I should also be large enough to avoid any PBC effect on liquid II. In order to account for the liquid–liquid–air system, such as the froth oilsalt water–air interface [10], or the surfactant-water–air interfaces [9], setting up the interface systems as shown in Fig. 1(c) or (d) is appropriate. Those systems in Fig. 1(c) and (d) include two kinds of liquids, liquid I and liquid II, plus air which is represented as vacuum in the surface tension calculations. The difference in model(c) and model(d) is that in model(d) only slight amount of liquid I is mixed with the major liquid II on the interface, thus the thickness of liquid I in model (d) usually is around 3–5 Å. In the past decade, updated CHARMM forcefield parameters have become available to describe n-alkanes with different chain lengths using the renown online program CGENFF [16,17]. More complicated organic molecules widely applied in industry can be described by CHARMM forcefield [16,18,19]. In order to find out if CHARMM forcefield can predict the surface tension, dynamics and structure of those molecules correctly, CHARMM forcefield and all-atom Nanoscale Molecular Dynamics (NAMD) simulations [20] were applied to work on nalkanes with the chain length from 8 to 24. Two kinds of interface systems were investigated, one is the alkane-air interface, and the other is the alkane-water interface. In order to study the interface effect on the dynamic properties of n-alkanes, simulations on pure alkane systems were also conducted. A wide temperature range from 300 K to 400 K was investigated. The surface tension of 9 n-alkanes with air and with water interface systems was calculated. The density profiles, structure properties including radial distribution function, radius of gyration, order parameters, surface tension, dynamics properties including diffusion coefficient of alkane in the bulk and on the interface were calculated and compared. The result can help to understand the chain length and temperature effect on the dynamics and properties of alkane-air/water interfaces.
simulation was continued using NAMD program and an NVT ensemble for another 11 ns. In total 9 alkane molecules were focused on, including C8H18, C9H20, C12H26, C13H28, C16H34, C17H36, C20H42, C21H44, C24H50. The mol2 files for those alkane molecules were built using ChemDraw Professional 17.0 edition. Using those mol2 files as input, the CGENFF online program was applied to build the CHARMM forcefield parameters. After that, CHARMM simulations were conducted to set up 9 pure alkane systems each having 148 n-alkane molecules packed randomly and with long enough distances between each other to avoid any initial confliction. After a brief energy minimization, the equilibration runs in an NPT ensemble, which means the number of atoms, pressure and temperature were kept constant during simulations, were performed for 1 ns at 1 atm and 300 K. Then the z-direction of the pure alkane system was extended to be 3 times of the original height to build the alkane-air interface system as shown in Fig. 3 model(a), which has the volume of empty space above and below the liquid being at least the alkane box size so that the two alkane-vacuum interfaces do not interact with each other even though periodic boundary conditions (PBCs) are employed in all the interface systems in x, y and z directions. Then, the NAMD simulation using an NVT (constant atom numbers, volume and temperature during simulations) ensemble was applied to continue the simulation for another 10 ns. The box sizes of the alkane-air interface system at 300 K are shown in Table 1. Besides that, in order to compare the dynamic properties of n-alkanes in pure alkane and in the bulk of interface systems, the pure n-alkane simulations on an NPT ensemble were also continued for another 10 ns for different alkanes starting from the initially equilibrated systems shown above at different temperatures. The size of the pure n-alkane simulations is shown in Table 1. In order to build the alkane-water interface system, TIP3P water models were applied to represent water molecules, with the bond length and angle restrained using SHAKE algorithm [21]. Two TIP3P water boxes each having the same size of the alkane system were appended to the alkane box, with one above and one below. The sketch of the alkane-water interface system is shown in Fig. 1 model(b). Thus the height of the box was tripled after adding water layers. For each water mixed system, another 1–2 ns simulations were performed in NPAT ensemble (constant atom number, pressure, surface area and temperature during simulation) at P = 1 atm and with the x and y dimensions of the surface area being equal to build the stable water-organic mixture interfaces, then simulations in NVT ensemble were continued for 8–10 ns. In all the NVT simulations, the long-range electrostatic interactions were calculated using the particle mesh Ewald method. The long-range cutoff for electrostatic interaction is 12 Å, and the van der Waals
2. Simulation method In order to test the CHARMM forcefield and the simulation method, one box of water molecule system which has in total 1981 TIP3P water molecules was set up first using CHARMM program. After the initial energy minimization and equilibration at 300 K and 1 atm for 1 ns using NAMD all-atom molecular dynamics simulations, the box size of the water system was 40.0×40.0×36.9 Å3. After that, empty spaces above and below the water box were appended, each of which has the height close to the height of the original water box. Such kind of watervacuum system in simulation is to represent the water–air system. The size of the simulation box became 40.0×40.0×110.7 Å3, and the 2
Fuel 267 (2020) 117252
G. Rucker, et al.
Table 1 Simulations performed and system details at 300 K. Box size of interface systems (Å3)
Alkanes
n-C8 n-C9 n-C12 n-C13 n-C16 n-C17 n-C20 n-C21 n-C24 water
Box size of pure alkane systems (Å3)
Alkane-air interface
Alkane-water interface
Pure n-alkane system
34.0 35.3 38.3 39.1 41.6 42.4 44.5 45.4 47.1 40.0
34.8 35.1 38.7 35.9 42.5 41.1 46.5 46.5 47.0
34.0 35.3 38.2 39.1 41.6 42.4 44.5 45.4 47.1
× × × × × × × × × ×
34.0 35.3 38.3 39.1 41.6 42.4 44.5 45.4 47.1 40.0
× × × × × × × × × ×
104.6 105.8 114.0 117.9 124.6 127.0 132.1 132.8 140.1 110.7
potential energy is switched to zero over the range of 10–12 Å. A time step of 1 fs during the energy minimization, equilibration, and sampling runs was applied. The frame output frequency is 1 ps. The box sizes for different systems at 300 K are shown in Table 1. In order to consider the influence of temperature on the alkane properties, simulations were performed at 5 different temperatures: 300 K, 340 K, 360 K, 380 K, and 400 K. Based on the sampling run output, the surface tension was calculated from the pressure tensor output in x, y and z directions (denoted using PXX, PYY, and PZZ) [22,23], following Eq. (1):
=
1 LZ Pzz 2
× × × × × × × × ×
34.8 35.1 38.7 35.9 42.5 41.1 46.5 46.5 47.0
(Pxx + Pyy )
106.5 107.3 118.3 109.8 110.1 109.7 114.2 114.3 125.4
× × × × × × × × ×
34.0 35.3 38.2 39.1 41.6 42.4 44.5 45.4 47.1
× × × × × × × × ×
34.8 35.3 38.0 39.3 41.5 42.3 44.1 44.3 46.7
chosen to be consistent with Cummings et al [14] since we work on alkanes with the carbon number from 8 to 24. Both the upper layer and lower layer are labelled as shown in Fig. 2 (Left) and (Right) for both kinds of interface systems. Besides that, the diffusion coefficient of n-alkanes in the bulk and in the interface region (including the upper layer and lower layer regions) in both alkane-water and alkane-air interface systems were calculated based on the mean squared displacement of molecule center of mass (COM) [24], with Eq. (2):
D= 0.5
× × × × × × × × ×
(1)
d |r (t ) r (0)|2 1 lim 2d t dt
(2)
Here, d is the dimensionality, which is 3 for the bulk diffusion coefficient calculation; t is the simulation time, while r(t) is the position of molecule COM at time t, and r(0) the position of molecule at time zero. < … > means ensemble average over the simulation time. When applying Eq. (2) to calculate the diffusion coefficient of n-alkanes in the interface region, only alkanes inside the region are included. Once an alkane moves out of the interface region by 1 ps, it will not be counted. The diffusion coefficient of alkanes in the interface region can be further divided into the lateral direction (x and y directions) Dxy, and z direction Dz, and can be calculated separately using Eqs. (3) and (4), based on the COM displacement of alkanes inside the interface slabs (upper layer slab and lower layer slab) as a function of time t in the lateral and z directions individually. The sum is over all the molecule i that are in slab during the time interval [0, t+ ]. N(t) is the number of
Here, ϒ is the surface tension, LZ is the length of the simulation box in z-direction, as shown in the snapshots of the simulation box for the alkane-air interface system and the alkane-water interface in Fig. 2(Left) and Fig. 2(Right). The angular brackets denote the time average, and ½ accounts for the slab having two interfaces. In Fig. 2, the LX and LY are the length of the simulation box in x and y directions. In all the systems at different temperatures, the Lz is around 3 times of the original alkane box height in order to avoid any interaction between interfaces. In order to consider the interface effect on the dynamics and structure of alkane on the interfaces, alkanes in the distances below the top surface of alkane up to 6 Å (called the upper layer), and above the bottom surface of alkane up to 6 Å (called the lower layer), are defined as alkanes in the interface region. A interface thickness of 6 Å was
Fig. 2. Simulation box of 148 n-C8 molecules with air interface (Left), and 148 n-C9 molecules slab with water interface (Right) at 300 K. The n-alkane molecules are shown in cyan sticks, while TIP3P water molecules in red and white sticks. The Lz, Lx, Ly are the dimensions of the simulation box in z, x and y directions. The upper layer and lower layer are the alkanes within 6 Å to the top and bottom of the alkane slab respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
3
Fuel 267 (2020) 117252
G. Rucker, et al.
data points that contribute to the sum in the interface slab at time t. Once the alkane's COM is outside of the interface slab for more than 1 ps, the molecule’s t and are set to zero before it moves into the interface slab again.
Dxy = lim t
Dz = lim t
1 4N (t ) t
1 2N ( t ) t
[x i (t + )
x i ( )]2 + [yi (t + )
yi ( )]2
i,
[z i (t + )
z i ( )]2
i,
(3) (4)
The radial distribution function (g(r)) was calculated based on molecule COM distances between molecule pairs during the sampling runs. The density profiles of n-alkanes at different temperatures were calculated based on computing the 1-D projections of atomic densities by averaging result over all the frames generated during the sampling runs using the VMD plugin [25]. The order parameter [26] which represents the orientational preference of the alkane chain packing with each other was calculated based on the vector built from the alkane head carbon atom to the tail carbon atom (vector built from C1 to CN) following the same method as our previous work [27]. The order parameters of alkane on the alkaneair/water interface regions and in the bulk which includes all the alkane molecules were calculated separately and compared. Besides that, the radius of gyration (Rg) of n-alkanes was calculated for alkanes in the bulk and in the interfacial region at different temperatures as well.
Fig. 4. Density profiles of n-alkane-water interfaces at 300 K for 9 systems. The mass density profile of water is from the n-C8-water system at 300 K.
result in both kinds of interface systems. Symmetric but having more fluctuations in density profiles were observed in the water-interface systems especially for long-chain alkanes. However, more smooth by not necessarily symmetric density profiles were observed in the airinterface systems. That should relate to the interaction between the interface media and alkane. Because of the interaction between water and alkanes but not in vacuum-alkane system, the mass distribution of n-alkanes on the water-interface is more symmetric.
3. Results and discussion 3.1. Simulation system and structure result
3.1.2. Radial distribution function result at different temperatures Calculating the radial distribution function of alkane molecule pairs inside the system, the result is shown in Fig. 5 for the even-carbonnumbered alkanes and in Figure S1 for the odd-carbon-numbered alkanes. The packing of n-C8 and n-C12 can reach a random packing at 5 different temperatures. However, n-C16 crystalized at 300 K in the airinterface system by showing a very sharp peak, although not in the water-interface system. Because of the long-chain alkane crystallization at 300 K, n-C16 has a preferred packing order and packing distance instead of a random packing. Similarly, n-C20 and n-C24, and n-C17 and n-C21 crystalized at 300 K in both the air-interface and water-interface systems. The closest distance between alkane-alkane pair is around 5 Å with the carbon number increasing from C16 to C24. The crystallization of long-chain alkanes predicted in the radial distribution function result can agree with the order parameter result shown in the Section 3.1.4. More comparison of g(r) for all alkanes is shown in Figure S2 in the Supplemental Material.
3.1.1. Density profiles of alkane-air/water interface system comparison The mass density profiles of the different n-alkane systems along the z-direction at 300 K were calculated, with the result for the alkane-air shown in Fig. 3 and alkane-water systems shown in Fig. 4. With the carbon number increasing, the mass density of n-alkanes increases. At 20 °C, the density of octane is 0.703 g/ml, 0.7176 g/ml for n-C9H20, 0.75 g/ml for n-C12H26, 0.7564 g/ml for n-C13H28 and 0.7701 g/ml for n-C16H34 [28]. Because of the extension of the alkane-air interface system, the densities of n-alkanes are slightly lower than the experimental density data at 300 K. Because of the partial crystallization of long-chain alkanes at 300 K, the density profiles of n-C16, n-C17, n-C20 to n-C24 are not smooth in the height range from −20 Å to 20 Å, but show a valley. In the n-alkane-water interface systems, the densities of n-alkanes at 300 K are very close to the experimental data as listed above. At 300 K, the density of water is 1.0 g/ml, which is consistent with the experimental data. Differences in density profiles were observed comparing
3.1.3. Radius of gyration comparison In order to quantitatively measure the structure change of n-alkanes during the simulation time, Root Mean Squared deviation (RMSD) of an example n-alkane molecule was calculated based on the position of heavy carbon atoms after aligning the trajectory on the initial structure. The RMSD result for the first n-alkane in the simulation box in the airinterface systems and the water-interface systems are shown in Fig. S3 (Left) and (Right), while the result in the pure alkane systems in Fig. S4. As can be seen, different alkane molecules reached equilibrated state in 1 ns and the structure of n-alkane fluctuated evenly during the simulation time. Calculating the size of alkanes staying in the bulk, on the interface and the standard deviations, we have the example Rg result at 300 K shown in Fig. 6(Left). The standard deviations of n-alkanes in the bulk and on the interface are up to 10% of the size of Rg. To avoid overlapping and show data clearly, only the standard deviations of the Rg in the bulk of water-interface at 300 K are shown as an example. Goodsaid-Zalduondo and Engelman measured the radius of gyration of several n-alkanes using neutral scattering method, and got Rg of n-C8 to be
Fig. 3. Density profiles of n-alkane-air interfaces at 300 K. 4
Fuel 267 (2020) 117252
G. Rucker, et al.
Fig. 5. Radial distribution function of n-C8, n-C12, n-C16, n-C20, and n-C24 in the air-interface (labelled as air in the legend) and water-interface (labelled as water in the legend) systems at five different temperatures. The magenta solid line is y = 1.0. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
3.7 ± 0.2 Å, n-C12 to be 4.5 ± 0.2 Å, n-C16 to be 4.8 ± 0.2 Å, and 5.1 ± 0.2 Å for n-C20 [29] in dilute solution of perdeuterated chains in hydrogenated chain solvents of the same length at 45˚C. Our simulation prediction on n-C8, n-C12 are smaller than the experimental data by up to 1 Å, the Rg of n-C16 predicted is very close, while Rg of n-C20 predicted is higher than experimental data by no more than 1 Å. However, very consistent result was reached when comparing our simulation prediction with other simulation prediction such as using Monte Carlo simulation method [30]: Rg of n-C8 from the MC prediction is in the range of 2.74 to 2.78 Å while 4.7 to 5.1 Å for n-C16 in the temperature range of 303 K to 373 K. Those prove that our simulation prediction is reasonable. Fig. 6 (Left) also shows that with the alkane carbon number increasing from 8 to 24, the Rg of alkanes increases almost linearly at 300 K. No matter for alkanes in the upper/lower interface or in the
bulk, the alkane sizes are very similar. The alkane size also has negligible dependence on temperatures as shown in Fig. 6(Right), which also agrees with the MC simulation prediction [30]. The size of alkane in the alkane-air interface systems is very close to the alkane-water systems, and to pure alkane systems (shown in green open symbols in Fig. 6(Right)). 3.1.4. Order parameter result Calculating the order parameters of alkanes in the bulk, and on the interface regions of the alkane-air, alkane-water systems, the result in the temperature range of 340 K to 400 K is shown in Fig. 7(Left) and (Right). As can be seen, the bulk of alkane molecules packed almost randomly in both kinds of systems by showing very small order parameter
Fig. 6. Rg comparison for 9 alkanes in the alkane-water interface systems at 300 K (Left) in the bulk (black circles), on the upper interface (red squares), and on the lower interface (blue diamonds) with experimental data [29] and MC simulation data[30]; and for alkanes in the water-interface (in black) and air-interface (in red), and pure alkane (in green) systems at 340 K, 360 K, 380 K and 400 K (Right). Only the standard deviations of Rg in the bulk of water-interface at 300 K and pure alkane systems at different temperatures are shown to be clear. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 5
Fuel 267 (2020) 117252
G. Rucker, et al.
Fig. 7. Order parameter result comparison for alkanes in the alkane-air interface system (Left) and alkane-water system (Right) in the bulk (in black open symbols), on the interface layer (in solid symbols) (in blue, and orange) at four different temperatures: 340 K (shown by red circles), 360 K (by green squares), 380 K (by blue diamonds), and 400 K (by orange upside-triangles). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
result for n-alkanes at different chain lengths and at different temperatures. Similarly, the alkanes packed randomly in the pure alkane systems at the same temperature range with result shown in Figure S5. However, much larger order parameter result for n-C17, n-C20, n-C21, nC24 at 300 K was observed because of crystallization in both the pure alkane system (shown in Figure S5) and air/water-interface system. That agrees with the result from our previous work, in which n-C22 packed randomly in liquid phase [27] at high temperatures, but crystallized at room temperature. Similarly, we calculated the order parameters of alkane in the interface regions in both kinds of systems. The order parameters of alkanes on the interface are larger than the ones in the bulk; thus molecules in the interface region are more orderly packed. Comparing the result in Fig. 7(Left) and (Right), alkanes are orderly packed on the alkane-air interface and on the alkane-water interface similarly. 3.2. Surface tension result 3.2.1. Surface tension result on pure water box Applying Eq. (1) to calculate the surface tension of the water–air system, it is 52.48 dyn/cm, which is lower than the experimental data which is around 71.7 dyn/cm at 300 K [31]. However, it is within the surface tension range from the seven flexible water models tested by Yuet and Blankschtein [32].
Fig. 8. Surface tension comparison of n-alkanes-air interface systems at five different temperatures.
3.2.2. Surface tension result for n-alkane-air interfaces For a total of 10 ns sampling runs, the surface tension from each 2 ns was calculated. Based on that, the average surface tension and the standard deviations were calculated accordingly, with the surface tension and the standard deviation comparison shown in Fig. 8. With the temperature increasing, the surface tension of alkane-air interface decreases almost linearly. Also, at the same temperature with the alkane chain length increasing, the surface tension increases. Only modest standard deviations of surface tension at different temperatures and for different n-alkanes are observed. With the alkane chain length increasing and temperature increasing, the standard deviation increases slightly. Comparing the surface tension of n-C8 to n-C16 predicted with the experimental data from Jasper et al. who measured the surface tension of n-alkanes using the capillary rise method [2], and from Somayajulu [33] who predicted the surface tension of n-alkanes using three-parameter equation, results are shown in Fig. 9. A satisfying agreement was reached in terms of surface tension dependence on temperature and on
Fig. 9. Comparison of alkane-air interface surface tension predicted from simulations with the available experimental data at different temperatures. Only the surface tension of n-C8H18 to n-C17H36 at five different temperatures is shown, since they have experimental data available from ref[2] in solid symbol and from ref[34] in half-solid symbol. 6
Fuel 267 (2020) 117252
G. Rucker, et al.
chain length in the alkane-air and alkane-water systems at different temperatures was investigated, as the result shown in Figure S6 in the Supplemental Material. With the carbon number increasing, the size of n-alkanes increases as shown in Fig. 6, thus the alkane diffusion coefficient decreases. Generally, the higher the temperatures, the higher the diffusion coefficients. 3.3.2. Diffusion coefficient results for alkanes in the interface region Calculating the diffusion coefficients of n-alkanes on the interface and comparing them with those in the bulk, Fig. 12(Left) shows that alkanes diffuse faster in the alkane-air interface region than in the bulk, and Fig. 12(Right) shows that the alkanes in the alkane-water interface systems also can diffuse faster than in the bulk. Since the diffusional movement can be decomposed into the lateral direction (x and y directions) and z-direction, the diffusion coefficients of molecules in the upper and lower interface regions in the x-y plane and in the z-direction were calculated separately, and the result at different temperatures are shown in Fig. 13(Left) for the alkane-air system and Fig. 13(Right) for the alkane-water system. The original data are shown in Table S1 and Table S2 in the Supplemental Material. As can be seen, the molecules can diffuse much faster on the x-y plane (shown as black circles) than in the z-direction (shown in red-squares). In general, the diffusion of alkanes in the lateral direction is 2 to 10 times faster than in the z-direction in both the air-interface and waterinterface systems.
Fig. 10. Surface tension and standard deviation of 9 different alkane-water interfaces from 300 K to 400 K. The experimental data from ref [35] for n-C8, nC9, and n-C12 are shown in solid symbols and connected using solid lines; while simulation data are shown in open symbols and connected using dashed lines.
the carbon number. The absolute surface tension data from the simulation predictions at the temperature range of 300–400 K are slightly lower than the experimental data by around 2 dyn/cm for each alkane. 3.2.3. Surface tension result for alkane-water interfaces Similarly, the surface tension and standard deviation of alkane with water at different temperatures were calculated using Eq. (1). As the result shown in Fig. 10, the surface tension of alkane-water decreases almost linearly with the temperature increasing. With the alkane chain length increasing, the surface tension of alkane-water interface increases. Comparing the simulation prediction with experimental data from ref [34] (shown as solid symbols) in Fig. 10, the predicted surface tension of n-C8, n-C9, and n-C12 is very close to the experimental data, especially for n-C12. Comparing the surface tension of the alkane on the interface with air and water shown in Figs. 8 and 10, the surface tension of each system almost doubled after the air was replaced by water. However, the molecule-air interface has a stronger temperature dependence than the corresponding alkane-water interfaces. Comparing to the surface tension result of alkane-air interface, the alkane-water interface surface tension showed a larger standard deviation and a larger discrepancy between the simulation prediction and experimental data. The reason should come from the TIP3P water model, which showed a large deviation from the experimental data in this work as shown in Section 3.2.1.
3.4. Discussion The surface tension of 9 n-alkanes was calculated on the interface with air and with water. Based on Nicolas and Smit [15], predicting the correct density profiles is necessary to predict the correct surface tension. Based on the density profile results shown in Figs. 3 and 4 which look reasonable in comparison with available experimental data, it is expected that the surface tension prediction should be reasonable. The surface tensions of n-alkanes on the air-interface are close to experimental data as shown in Fig. 9, but consistently smaller by around 2 dyn/cm at different temperatures and for all six n-alkanes from n-C8 to n-C17. Similarly, in the alkane-water systems, the surface tension predicted is smaller than the available experimental data by 4–5 dyn/cm. The small deviation in surface tension for alkane-air and alkane-water proves that the CHARMM forcefield parameters are reasonable for the surface tension calculation, although the surface tension of water predicted at 300 K using the same input parameters and procedure is smaller than the experimental data by around 20 dyn/cm. Based on our work, the surface tension of water-alkane and air-alkane interfaces is sensitive to the alkane chain length and temperature. With the alkane chain length increasing, the surface tension of the interface increases; with the temperature increasing, the surface tension of the interface decreases. In this paper, the order parameter between n-alkane molecules was calculated by connecting the head carbon to the tail carbon (C1 to CN) to build the unit vector. Other researchers built the unit vector using the carbon atom pair of CK to CK+2 in the order parameter calculation [36]. In order to check the consistency of the result, different carbon atom pairs such as C2-C4, C3-C5, and C6-C8 in the water-n-C8 system were tried at different temperatures, which gave the order parameter result close to the result as shown in Fig. 7 using the carbon atom pair of C1C8. On the interface, which is defined as the alkanes within 6 Å of thickness to the air/water phase in this project, alkanes prefer to pack more orderly than in the bulk. Michael and Benhamin [37] found that alkanes on the alkane-water interface prefer to pack laterally with the interface based on the order parameter result in lateral and z-directions. Because of the shape of n-alkanes, such kind of packing can explain why alkanes move faster in the lateral direction than in the z-direction in both kinds of interface systems found out in this work.
3.3. Diffusion coefficient result 3.3.1. Bulk diffusion coefficient result at different temperatures The bulk diffusion coefficients of alkanes in both the air-alkane systems and water-alkane systems were calculated which include all the alkanes molecules in the system, and result was compared with pure nalkane systems at different temperatures, and with available self-diffusion coefficient data from experimental work [35]. In order to show the result clearly, the even-carbon-numbered alkane results are shown in Fig. 11(Left), while the odd-carbon-numbered alkanes in Fig. 11(Right). In general, n-alkanes diffuse slower in the water-alkane systems than in the air-alkane systems, but alkanes diffuse much faster in both kinds of interface systems than in the pure n-alkane systems. The logD of alkanes has an almost linear relationship with 1/T in the interface systems and in the pure alkane systems. The simulation data on n-C8, n-C9, n-C12, n-C13 and n-C16 at 300 K are very close to experimental data at 303 K using the MRI method[36] (shown as patterned symbols). That proves that our diffusion coefficients predicted using the simulation method are reasonable. Similarly, the dependence of the diffusion coefficient on the alkane 7
Fuel 267 (2020) 117252
G. Rucker, et al.
Fig. 11. Bulk diffusion coefficient of n-C8, n-C12, n-C16, n-C20 and n-C24 (Left) and n-C9, n-C13, n-C17, and n-C21 (Right) in the air-interface systems (in empty symbols connected with solid line), in the water-interface systems (in solid symbols connected with dashed line) and in the pure alkane systems (in half-solid symbols connected with dot-dashed line) in comparison with experimental data [36] (in patterned symbols). Data at 300 K for n-C17, n-C20, n-C21, and n-C24 are not shown because of crystallization in pure-alkane, air-interface and water-interface systems.
Fig. 12. Diffusion coefficient of alkanes in the interface region (by squares) compared with in the bulk (by circles) for the alkane-air (Left) and alkane-water (Right) interface systems.
Fig. 13. Diffusion coefficients of 9 n-alkanes in the lateral direction and in the z-direction at five different temperatures in the air-alkane (Left) and water-alkane (Right) interface systems.
8
Fuel 267 (2020) 117252
G. Rucker, et al.
Up to now, some researchers have worked on the liquid–vapor interface and found the enhanced surface reactions in droplets and the air–water interface [38,39]. Liu et al. [40] calculated the diffusion tensor of water in confined fluids and interface and found that the lateral diffusion tensor of water was around 4 times of that in the bulk, and the normal diffusion tensor was around 2 times of that in the bulk. In this work, the diffusion coefficients of alkane at the interface region with air and with water were calculated. The diffusion coefficient was enhanced in both normal and lateral directions in both kinds of interface systems. More specifically, the lateral direction diffusion rate was around 2 ~ 5 times of that in the bulk, while the normal direction increased moderately. When comparing the diffusion coefficients of n-alkanes in the bulk of air/water-interfaces and in the pure alkane systems, faster diffusion coefficients in the system having interfaces were observed as shown in Fig. 11. That emphasized the effect of the interface interaction to the dynamic properties of alkanes. Comparing our diffusion coefficients predicted with experimental result using the MRI method [36], very consistent agreement was found. Comparing the radius of gyration result and radial distribution function result for both kinds of interface systems with available experimental data, consistent agreements were reached for both individually. With the n-alkane chain length increasing, the size of n-alkane increases, which has negligible dependence on temperature and on the position of alkanes in the interface systems.
resources were provided through the high-performance computers at TTU and resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.fuel.2020.117252. References [1] Roberts FL, Kandhal PS, Brown ER, Lee D-Y, Kennedy TW. Hot Mix Asphalt Materials, Mixture Design, and Construction, 2nd; National Asphalt Pavement Association: Lanham, MD; 1996, Chapter2. [2] Jasper JJ, Kerr ER, Gregorich F. The orthobaric surface tensions and thermodynamic properties of the liquid surfaces of the n-aIkanes, C5 to C28. JACS 1953;75:5252–4. [3] Ismail A, Tsige M, Veld in PJ, Grest G. Surface tension of normal and branched alkanes. Mol Phys 2007;105(23):3155–63. [4] Zhang Y, Feller SE, Brooks BR, Pastor RW. Computer simulation of liquid/liquid interfaces. I. Theory and application to octane/Water. J Chem Phys 1995;103:10252. https://doi.org/10.1063/1.469927. [5] Morrow NR. Interfacial Phenomena in Petroleum Recovery. New York: Dekker; 1991. p. 1. [6] Thomas JL, Roeselová M, Dang LX, Tobias DJ. Molecular dynamics simulations of the solution−air interface of aqueous sodium nitrate. J Phys Chem A 2007;111(16):3091–8. [7] Nguyen CV, Phan CM, Ang HM, Nakahara H, Shibata O, Moroi Y. Molecular dynamics investigation on adsorption layer of alcohols at the air/brine interface. Langmuir 2015;31(1):50–6. [8] Yuet PK, Blankschtein D. Molecular dynamics simulation study of water surfaces: comparison of flexible water models. J Phys Chem B 2010;114(43):13786–95. [9] Zhang L, Liu Z, Ren T, Wu P, Shen J-W, Zhang W, et al. Understanding the structure of hydrophobic surfactants at the air/ water interface from molecular level. Langmuir 2014;30:13815–22. [10] Chong L, Lai Y, Gray M, Soong Y, Shi F, Duan Y. Effects of frothers and oil at saltwater–air interfaces for oil separation: molecular dynamics simulations and experimental measurements. J Phys Chem B 2017;121(27):6699–707. [11] Alejandrea J, Tildesley DJ. Molecular dynamics simulation of the orthobaric densities and surface tension of water. J Chem Phys 1995;102(11):4574–83. [12] Carpenter IL, Hehre WJ. A molecular dynamics study of the hexane/water interface. J Phys Chem 1990;94:531–6. [13] van Buuren AR, Marrink SJ, Berendsen HJC. A molecular dynamics study of the decane/water interface. J Phys Chem 1993;97:9206–12. [14] Rivera JL, McCabe C, Cummings PT. Molecular simulations of liquid-liquid interfacial properties: water–n-alkane and water-methanol–n-alkane systems. Phys Rev E 2003;67. 011603. [15] Nicolas JP, Smit B. Molecular dynamics simulations of the surface tension of nhexane, n-decane and n-hexadecane. Mol Phys 2002;100(15):2471–5. [16] Vanommeslaeghe K, Hatcher E, Acharya C, Kundu S, Zhong S, Shim J, et al. CHARMM general force field: a force field for drug-like molecules compatible with the CHARMM all-atom additive biological force field. J Comput Chem 2010;31:671. [17] Yu W, He X, Vanommeslaeghe K, MacKerell Jr. AD. Extension of the CHARMM general force field to sulfonyl-containing compounds and its utility in biomolecular simulations. J Comput Chem 2012;33:2451. [18] Vanommeslaeghe K, MacKerell Jr. AD. Automation of the CHARMM general force field (CGenFF) I: bond perception and atom typing. J Chem Inf Model 2012;52:3144. [19] Vanommeslaeghe K, Raman EP, MacKerell Jr. AD. Automation of the CHARMM general force field (CGenFF) II: assignment of bonded parameters and partial atomic charges. J Chem Inf Model 2012;52:3155. [20] Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, et al. Scalable molecular dynamics with NAMD. J Comput Chem 2005;26:1781. [21] Ryckaert J-P, Ciccotti G, Berendsen HJC. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J Comput Phys 1977;23:327. [22] Hill TL. Introduction to Statistical Mechanics. New York: Dover; 1986. p. 314. [23] Ono S, Kondo S. Molecular Theory of Surface Tension in Liquids. Berlin: Springer; 1960. p. 134. [24] Frenkel D, Smit B. Understanding Molecular Simulation. 2nd San Diego: Academic; 2002. [25] Giorgino T. Computing 1-D atomic densities in macromolecular simulations: the density profile tool for VMD. Comput Phys Commun 2014;185(1):317–22. [26] Cecchini M, Rao F, Seeber M, Caflisch A. Replica exchange molecular dynamics simulations of amyloid peptide aggregation. J Chem Phys 2004;121:10748. [27] Sonibare K, Rathnayaka L, Zhang L. Comparison of CHARMM and OPLS-aa forcefield predictions for components in one model asphalt mixture. Constr Build Mater 2020;236. 117577. [28] Griesbaum K, Behr A, Biedenkapp D, Voges H-W, Garbe D, Paetz C, et al. “Hydrocarbons” in Ullmann's Encyclopedia of Industrial Chemistry. Weinheim: Wiley-VCH; 2005.
4. Conclusion This paper works on 9 n-alkanes with the air and with the water interface systems at five different temperatures, ranging from 300 K to 400 K. The density profiles, surface tension, diffusion coefficients, radius of gyration, radial distribution function, and order parameters of nalkanes in the bulk and at the interface region in both kinds of interface systems were calculated at different temperatures. It was found that the density profiles of n-alkanes can agree with available experimental data at the room temperature. The higher the chain length of n-alkanes, the higher the surface tension in both kinds of interface systems. The wateralkane interface usually has a higher surface tension than the air-alkane interface. At the same temperature, the alkanes can diffuse faster in the air-alkane interface system than in the water-alkane interface system, and faster in the interface system than in the pure-alkane system. Focusing on the alkanes on the interface in both kinds of systems, molecules are more orderly packed on the interface than in the bulk, and they also diffuse much faster than in the bulk. Decomposing the diffusion coefficients of alkanes into the lateral direction and the z-direction, alkanes on the interface can diffuse much faster in the lateral direction than in the z-direction. CRediT authorship contribution statement George Rucker: Data curation. Xiong Yu: Conceptualization. Liqun Zhang: Conceptualization, Methodology, Software, Investigation. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements We acknowledge TTU student Yilun Lee for helping with the initial setting up and analyzing some simulation trajectories. This work was supported by the Center for Energy Systems Research at TTU and by the Faculty Research Fund TTU provided to Dr. Zhang. The computer 9
Fuel 267 (2020) 117252
G. Rucker, et al. [29] Goodsaid-Zalduondo F, Engelman DM. Conformation of liquid N-alkanes. Biophys J 1981;35(3):587–94. [30] Bessières D, Piñeiro MM, De Ferron G, Plantier F. Analysis of the orientational order effect on n-alkanes: evidences on experimental response functions and description using Monte Carlo molecular simulation. 074507 J Chem Phys 2010;133. https:// doi.org/10.1063/1.3472283. [31] Vargaftik NB, Volkov BN, Voljak LD. International tables of the surface tension of water. J Phys Chem Ref Data 1983;12:817–20. [32] Yuet PK, Blankschtein D. Molecular dynamics simulation study of water surfaces: comparison of flexible water models. J Phys Chem B 2010;114:13786–95. [33] Somayajulu GR. A generalized equation for surface tension from the triple point to the critical point. Int J Thermophys 1988;9(4):559–66. [34] Zeppieri S, Rodrı́guez J, López de Ramos AL. Interfacial tension of alkane + water systems. J Chem Eng Data 2001;46:1086–8. [35] Tofts PS, Lloyd D, Clark CA, Barker GJ, Parker GJM, McConville P, et al. Test liquids
[36] [37] [38] [39] [40]
10
for quantitative MRI measurements of self-diffusion coefficient in vivo. Magn Reson Med 2000;43:368–74. Harris JG. Liquid-vapor interfaces of alkane oligomers: structure and thermodynamics from molecular dynamics simulations of chemically realistic models. J Phys Chem 1992;96:5077–86. Michael D, Benjamin I. Solute orientational dynamics and surface roughness of water/hydrocarbon interfaces. J Phys Chem 1995;99:1530–6. Liu P, Harder E, Berne BJ. On the calculation of diffusion coefficients in confined fluids and interfaces with an application to the liquid−vapor interface of water. J Phys Chem B 2004;108(21):6595–602. Valsaraj KT. Trace gas adsorption thermodynamics at the air−water interface: Implications in atmospheric chemistry. Pure Appl Chem 2009;81(10):1889–901. Nissenson P, Knox CJH, Finlayson-Pitts BJ, Phillips LF, Dabdub D. Enhanced photolysis in aerosols: evidence for important surface effects. Phys Chem Chem Phys 2006;8(40):4700–10.
Contents lists available at ScienceDirect
Fuel journal homepage: www.elsevier.com/locate/fuel
Full Length Article
Molecular dynamics investigation on n-alkane-air/water interfaces a
b
a,⁎
George Rucker , Xiong Yu , Liqun Zhang a b
T
Department of Chemical Engineering, Tennessee Technological University, Cookeville, TN, 38505, United States Department of Civil and Environmental Engineering, Case Western Reserve University, Cleveland, OH, 44106, United States
ARTICLE INFO
ABSTRACT
Keywords: Surface tension Alkane Interface Diffusion coefficient Order parameters Molecular dynamics simulation
Alkanes are the important molecules in crude oil that have applications in different fields. In this project, molecular dynamics simulations have been applied to investigate the n-alkane-air and n-alkane-water interfaces. The surface tension, diffusion coefficients, order parameters, radial distribution function, density profiles and radius of gyration of alkanes in the bulk and on the interfaces with air/water were calculated and compared. It was found that with the n-alkane chain length increasing, the surface tension of the interface system increases. However, with the temperature increasing, the surface tension of the alkanes decreases. The alkanes on the interface can diffuse faster than in the bulk in both kinds of interface systems. While n-alkane diffusion rate in the bulk of interface systems is faster than in the pure alkane systems, the radius of gyration of alkanes in the bulk, on the interface and in the pure alkane systems is consistent. Decomposing the diffusion coefficients into the lateral direction and z-direction, the lateral diffusion rate of alkanes on the interface is about 2–10 times faster than in the z-direction in both kinds of interface systems. Although molecules prefer to pack randomly in the bulk, they are more orderly packed on the interface. The result can help to explain the interface effect in water/air-alkane interface systems. The project also demonstrates that using NAMD program and CHARMM forcfield can predict reasonable structure and properties of n-alkane-air/water systems.
1. Introduction Alkanes are an important class of molecules that have various applications in different fields. They are natural nonpolar solvents, and are the major components in crude oil [1]. The structure and dynamics of the interface of alkane with air and with water are of great interest to researchers [2–4] and of great importance to its application in chemistry, physics, engineering and biology, such as the application of water flooding and surfactant flooding in enhanced crude oil recovery [5]. Besides experimental methods, molecular dynamics simulation methods have been applied to address the liquid/air and liquid/liquid interface properties, and satisfying agreements between simulation prediction and experimental data were reached [6–11]. Carpenter and Hehre [12] studied the n-hexane/water interface system using molecular dynamics simulation method, and found differences between the interfacial and bulk liquids. Zhang et al [4] did CHARMM all-atom molecular dynamics simulations on neat n-octane/water and n-octane/ air interfaces. They predicted the surface tension of both kinds of interfaces using different simulation ensembles. Consistent agreement was reached from different ensembles and the predictions were also consistent with experimental work. Van Buuren et al. [13] worked on the decane/water interface, and found that decane molecules on the ⁎
water interface preferred to orient laterally. Rivera et al. [14] worked on the n-alkane/water interfaces and found that with the alkane chain length increasing, the alkane-water surface tension increases. Nicolas and Smit computed the surfaced tension of linear alkanes including nhexane, n-decane and n-hexadecane using united atom OPLS forcefield and SKS forcefield [15]. They found that both models can predict the surface tension of three alkanes consistently agreeing with experimental work at high temperatures. In order to predict the surface tension correctly, it is critical that the model can predict the density correctly. However, up to now, no research on the alkane-air and alkanewater interfaces has been done in a wide carbon number range and in a wide temperature range using CHARMM forcefield, and systematically comparing results from both kinds of systems. Thus, this project studies alkanes from n-C8 to n-C24 in both water and air interfaces, and in the temperature range of 300–400 K. Generally, there are four different kinds of liquid/air, liquid/liquid and liquid–liquid–air interfaces as the sketches shown in Fig. 1. In order to calculate the surface tension of liquid–air interfaces, it is straightforward to set up the interface system as Fig. 1(a), in which the air is represented by vacuum in simulations. In order to calculate the surface tension of liquid–liquid interfaces, it is straightforward to set up the interface system as Fig. 1(b) in which the same amount of liquid I is
Corresponding author at: 1020 Stadium Dr., Cookeville, TN 38505, United States. E-mail address: [email protected] (L. Zhang).
https://doi.org/10.1016/j.fuel.2020.117252 Received 16 November 2019; Received in revised form 27 January 2020; Accepted 28 January 2020 Available online 08 February 2020 0016-2361/ © 2020 Elsevier Ltd. All rights reserved.
Fuel 267 (2020) 117252
G. Rucker, et al.
Fig. 1. Sketches of the typical liquid–liquid/vapor interface systems: Liquid-vacuum interface (a); liquid–liquid interface (b); and two kinds of liquid–liquid-vapor interfaces (c) and (d). The vacuum has been applied to represent the air in simulations.
above and below liquid II in the simulation box. The amount of liquid I should also be large enough to avoid any PBC effect on liquid II. In order to account for the liquid–liquid–air system, such as the froth oilsalt water–air interface [10], or the surfactant-water–air interfaces [9], setting up the interface systems as shown in Fig. 1(c) or (d) is appropriate. Those systems in Fig. 1(c) and (d) include two kinds of liquids, liquid I and liquid II, plus air which is represented as vacuum in the surface tension calculations. The difference in model(c) and model(d) is that in model(d) only slight amount of liquid I is mixed with the major liquid II on the interface, thus the thickness of liquid I in model (d) usually is around 3–5 Å. In the past decade, updated CHARMM forcefield parameters have become available to describe n-alkanes with different chain lengths using the renown online program CGENFF [16,17]. More complicated organic molecules widely applied in industry can be described by CHARMM forcefield [16,18,19]. In order to find out if CHARMM forcefield can predict the surface tension, dynamics and structure of those molecules correctly, CHARMM forcefield and all-atom Nanoscale Molecular Dynamics (NAMD) simulations [20] were applied to work on nalkanes with the chain length from 8 to 24. Two kinds of interface systems were investigated, one is the alkane-air interface, and the other is the alkane-water interface. In order to study the interface effect on the dynamic properties of n-alkanes, simulations on pure alkane systems were also conducted. A wide temperature range from 300 K to 400 K was investigated. The surface tension of 9 n-alkanes with air and with water interface systems was calculated. The density profiles, structure properties including radial distribution function, radius of gyration, order parameters, surface tension, dynamics properties including diffusion coefficient of alkane in the bulk and on the interface were calculated and compared. The result can help to understand the chain length and temperature effect on the dynamics and properties of alkane-air/water interfaces.
simulation was continued using NAMD program and an NVT ensemble for another 11 ns. In total 9 alkane molecules were focused on, including C8H18, C9H20, C12H26, C13H28, C16H34, C17H36, C20H42, C21H44, C24H50. The mol2 files for those alkane molecules were built using ChemDraw Professional 17.0 edition. Using those mol2 files as input, the CGENFF online program was applied to build the CHARMM forcefield parameters. After that, CHARMM simulations were conducted to set up 9 pure alkane systems each having 148 n-alkane molecules packed randomly and with long enough distances between each other to avoid any initial confliction. After a brief energy minimization, the equilibration runs in an NPT ensemble, which means the number of atoms, pressure and temperature were kept constant during simulations, were performed for 1 ns at 1 atm and 300 K. Then the z-direction of the pure alkane system was extended to be 3 times of the original height to build the alkane-air interface system as shown in Fig. 3 model(a), which has the volume of empty space above and below the liquid being at least the alkane box size so that the two alkane-vacuum interfaces do not interact with each other even though periodic boundary conditions (PBCs) are employed in all the interface systems in x, y and z directions. Then, the NAMD simulation using an NVT (constant atom numbers, volume and temperature during simulations) ensemble was applied to continue the simulation for another 10 ns. The box sizes of the alkane-air interface system at 300 K are shown in Table 1. Besides that, in order to compare the dynamic properties of n-alkanes in pure alkane and in the bulk of interface systems, the pure n-alkane simulations on an NPT ensemble were also continued for another 10 ns for different alkanes starting from the initially equilibrated systems shown above at different temperatures. The size of the pure n-alkane simulations is shown in Table 1. In order to build the alkane-water interface system, TIP3P water models were applied to represent water molecules, with the bond length and angle restrained using SHAKE algorithm [21]. Two TIP3P water boxes each having the same size of the alkane system were appended to the alkane box, with one above and one below. The sketch of the alkane-water interface system is shown in Fig. 1 model(b). Thus the height of the box was tripled after adding water layers. For each water mixed system, another 1–2 ns simulations were performed in NPAT ensemble (constant atom number, pressure, surface area and temperature during simulation) at P = 1 atm and with the x and y dimensions of the surface area being equal to build the stable water-organic mixture interfaces, then simulations in NVT ensemble were continued for 8–10 ns. In all the NVT simulations, the long-range electrostatic interactions were calculated using the particle mesh Ewald method. The long-range cutoff for electrostatic interaction is 12 Å, and the van der Waals
2. Simulation method In order to test the CHARMM forcefield and the simulation method, one box of water molecule system which has in total 1981 TIP3P water molecules was set up first using CHARMM program. After the initial energy minimization and equilibration at 300 K and 1 atm for 1 ns using NAMD all-atom molecular dynamics simulations, the box size of the water system was 40.0×40.0×36.9 Å3. After that, empty spaces above and below the water box were appended, each of which has the height close to the height of the original water box. Such kind of watervacuum system in simulation is to represent the water–air system. The size of the simulation box became 40.0×40.0×110.7 Å3, and the 2
Fuel 267 (2020) 117252
G. Rucker, et al.
Table 1 Simulations performed and system details at 300 K. Box size of interface systems (Å3)
Alkanes
n-C8 n-C9 n-C12 n-C13 n-C16 n-C17 n-C20 n-C21 n-C24 water
Box size of pure alkane systems (Å3)
Alkane-air interface
Alkane-water interface
Pure n-alkane system
34.0 35.3 38.3 39.1 41.6 42.4 44.5 45.4 47.1 40.0
34.8 35.1 38.7 35.9 42.5 41.1 46.5 46.5 47.0
34.0 35.3 38.2 39.1 41.6 42.4 44.5 45.4 47.1
× × × × × × × × × ×
34.0 35.3 38.3 39.1 41.6 42.4 44.5 45.4 47.1 40.0
× × × × × × × × × ×
104.6 105.8 114.0 117.9 124.6 127.0 132.1 132.8 140.1 110.7
potential energy is switched to zero over the range of 10–12 Å. A time step of 1 fs during the energy minimization, equilibration, and sampling runs was applied. The frame output frequency is 1 ps. The box sizes for different systems at 300 K are shown in Table 1. In order to consider the influence of temperature on the alkane properties, simulations were performed at 5 different temperatures: 300 K, 340 K, 360 K, 380 K, and 400 K. Based on the sampling run output, the surface tension was calculated from the pressure tensor output in x, y and z directions (denoted using PXX, PYY, and PZZ) [22,23], following Eq. (1):
=
1 LZ Pzz 2
× × × × × × × × ×
34.8 35.1 38.7 35.9 42.5 41.1 46.5 46.5 47.0
(Pxx + Pyy )
106.5 107.3 118.3 109.8 110.1 109.7 114.2 114.3 125.4
× × × × × × × × ×
34.0 35.3 38.2 39.1 41.6 42.4 44.5 45.4 47.1
× × × × × × × × ×
34.8 35.3 38.0 39.3 41.5 42.3 44.1 44.3 46.7
chosen to be consistent with Cummings et al [14] since we work on alkanes with the carbon number from 8 to 24. Both the upper layer and lower layer are labelled as shown in Fig. 2 (Left) and (Right) for both kinds of interface systems. Besides that, the diffusion coefficient of n-alkanes in the bulk and in the interface region (including the upper layer and lower layer regions) in both alkane-water and alkane-air interface systems were calculated based on the mean squared displacement of molecule center of mass (COM) [24], with Eq. (2):
D= 0.5
× × × × × × × × ×
(1)
d |r (t ) r (0)|2 1 lim 2d t dt
(2)
Here, d is the dimensionality, which is 3 for the bulk diffusion coefficient calculation; t is the simulation time, while r(t) is the position of molecule COM at time t, and r(0) the position of molecule at time zero. < … > means ensemble average over the simulation time. When applying Eq. (2) to calculate the diffusion coefficient of n-alkanes in the interface region, only alkanes inside the region are included. Once an alkane moves out of the interface region by 1 ps, it will not be counted. The diffusion coefficient of alkanes in the interface region can be further divided into the lateral direction (x and y directions) Dxy, and z direction Dz, and can be calculated separately using Eqs. (3) and (4), based on the COM displacement of alkanes inside the interface slabs (upper layer slab and lower layer slab) as a function of time t in the lateral and z directions individually. The sum is over all the molecule i that are in slab during the time interval [0, t+ ]. N(t) is the number of
Here, ϒ is the surface tension, LZ is the length of the simulation box in z-direction, as shown in the snapshots of the simulation box for the alkane-air interface system and the alkane-water interface in Fig. 2(Left) and Fig. 2(Right). The angular brackets denote the time average, and ½ accounts for the slab having two interfaces. In Fig. 2, the LX and LY are the length of the simulation box in x and y directions. In all the systems at different temperatures, the Lz is around 3 times of the original alkane box height in order to avoid any interaction between interfaces. In order to consider the interface effect on the dynamics and structure of alkane on the interfaces, alkanes in the distances below the top surface of alkane up to 6 Å (called the upper layer), and above the bottom surface of alkane up to 6 Å (called the lower layer), are defined as alkanes in the interface region. A interface thickness of 6 Å was
Fig. 2. Simulation box of 148 n-C8 molecules with air interface (Left), and 148 n-C9 molecules slab with water interface (Right) at 300 K. The n-alkane molecules are shown in cyan sticks, while TIP3P water molecules in red and white sticks. The Lz, Lx, Ly are the dimensions of the simulation box in z, x and y directions. The upper layer and lower layer are the alkanes within 6 Å to the top and bottom of the alkane slab respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
3
Fuel 267 (2020) 117252
G. Rucker, et al.
data points that contribute to the sum in the interface slab at time t. Once the alkane's COM is outside of the interface slab for more than 1 ps, the molecule’s t and are set to zero before it moves into the interface slab again.
Dxy = lim t
Dz = lim t
1 4N (t ) t
1 2N ( t ) t
[x i (t + )
x i ( )]2 + [yi (t + )
yi ( )]2
i,
[z i (t + )
z i ( )]2
i,
(3) (4)
The radial distribution function (g(r)) was calculated based on molecule COM distances between molecule pairs during the sampling runs. The density profiles of n-alkanes at different temperatures were calculated based on computing the 1-D projections of atomic densities by averaging result over all the frames generated during the sampling runs using the VMD plugin [25]. The order parameter [26] which represents the orientational preference of the alkane chain packing with each other was calculated based on the vector built from the alkane head carbon atom to the tail carbon atom (vector built from C1 to CN) following the same method as our previous work [27]. The order parameters of alkane on the alkaneair/water interface regions and in the bulk which includes all the alkane molecules were calculated separately and compared. Besides that, the radius of gyration (Rg) of n-alkanes was calculated for alkanes in the bulk and in the interfacial region at different temperatures as well.
Fig. 4. Density profiles of n-alkane-water interfaces at 300 K for 9 systems. The mass density profile of water is from the n-C8-water system at 300 K.
result in both kinds of interface systems. Symmetric but having more fluctuations in density profiles were observed in the water-interface systems especially for long-chain alkanes. However, more smooth by not necessarily symmetric density profiles were observed in the airinterface systems. That should relate to the interaction between the interface media and alkane. Because of the interaction between water and alkanes but not in vacuum-alkane system, the mass distribution of n-alkanes on the water-interface is more symmetric.
3. Results and discussion 3.1. Simulation system and structure result
3.1.2. Radial distribution function result at different temperatures Calculating the radial distribution function of alkane molecule pairs inside the system, the result is shown in Fig. 5 for the even-carbonnumbered alkanes and in Figure S1 for the odd-carbon-numbered alkanes. The packing of n-C8 and n-C12 can reach a random packing at 5 different temperatures. However, n-C16 crystalized at 300 K in the airinterface system by showing a very sharp peak, although not in the water-interface system. Because of the long-chain alkane crystallization at 300 K, n-C16 has a preferred packing order and packing distance instead of a random packing. Similarly, n-C20 and n-C24, and n-C17 and n-C21 crystalized at 300 K in both the air-interface and water-interface systems. The closest distance between alkane-alkane pair is around 5 Å with the carbon number increasing from C16 to C24. The crystallization of long-chain alkanes predicted in the radial distribution function result can agree with the order parameter result shown in the Section 3.1.4. More comparison of g(r) for all alkanes is shown in Figure S2 in the Supplemental Material.
3.1.1. Density profiles of alkane-air/water interface system comparison The mass density profiles of the different n-alkane systems along the z-direction at 300 K were calculated, with the result for the alkane-air shown in Fig. 3 and alkane-water systems shown in Fig. 4. With the carbon number increasing, the mass density of n-alkanes increases. At 20 °C, the density of octane is 0.703 g/ml, 0.7176 g/ml for n-C9H20, 0.75 g/ml for n-C12H26, 0.7564 g/ml for n-C13H28 and 0.7701 g/ml for n-C16H34 [28]. Because of the extension of the alkane-air interface system, the densities of n-alkanes are slightly lower than the experimental density data at 300 K. Because of the partial crystallization of long-chain alkanes at 300 K, the density profiles of n-C16, n-C17, n-C20 to n-C24 are not smooth in the height range from −20 Å to 20 Å, but show a valley. In the n-alkane-water interface systems, the densities of n-alkanes at 300 K are very close to the experimental data as listed above. At 300 K, the density of water is 1.0 g/ml, which is consistent with the experimental data. Differences in density profiles were observed comparing
3.1.3. Radius of gyration comparison In order to quantitatively measure the structure change of n-alkanes during the simulation time, Root Mean Squared deviation (RMSD) of an example n-alkane molecule was calculated based on the position of heavy carbon atoms after aligning the trajectory on the initial structure. The RMSD result for the first n-alkane in the simulation box in the airinterface systems and the water-interface systems are shown in Fig. S3 (Left) and (Right), while the result in the pure alkane systems in Fig. S4. As can be seen, different alkane molecules reached equilibrated state in 1 ns and the structure of n-alkane fluctuated evenly during the simulation time. Calculating the size of alkanes staying in the bulk, on the interface and the standard deviations, we have the example Rg result at 300 K shown in Fig. 6(Left). The standard deviations of n-alkanes in the bulk and on the interface are up to 10% of the size of Rg. To avoid overlapping and show data clearly, only the standard deviations of the Rg in the bulk of water-interface at 300 K are shown as an example. Goodsaid-Zalduondo and Engelman measured the radius of gyration of several n-alkanes using neutral scattering method, and got Rg of n-C8 to be
Fig. 3. Density profiles of n-alkane-air interfaces at 300 K. 4
Fuel 267 (2020) 117252
G. Rucker, et al.
Fig. 5. Radial distribution function of n-C8, n-C12, n-C16, n-C20, and n-C24 in the air-interface (labelled as air in the legend) and water-interface (labelled as water in the legend) systems at five different temperatures. The magenta solid line is y = 1.0. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
3.7 ± 0.2 Å, n-C12 to be 4.5 ± 0.2 Å, n-C16 to be 4.8 ± 0.2 Å, and 5.1 ± 0.2 Å for n-C20 [29] in dilute solution of perdeuterated chains in hydrogenated chain solvents of the same length at 45˚C. Our simulation prediction on n-C8, n-C12 are smaller than the experimental data by up to 1 Å, the Rg of n-C16 predicted is very close, while Rg of n-C20 predicted is higher than experimental data by no more than 1 Å. However, very consistent result was reached when comparing our simulation prediction with other simulation prediction such as using Monte Carlo simulation method [30]: Rg of n-C8 from the MC prediction is in the range of 2.74 to 2.78 Å while 4.7 to 5.1 Å for n-C16 in the temperature range of 303 K to 373 K. Those prove that our simulation prediction is reasonable. Fig. 6 (Left) also shows that with the alkane carbon number increasing from 8 to 24, the Rg of alkanes increases almost linearly at 300 K. No matter for alkanes in the upper/lower interface or in the
bulk, the alkane sizes are very similar. The alkane size also has negligible dependence on temperatures as shown in Fig. 6(Right), which also agrees with the MC simulation prediction [30]. The size of alkane in the alkane-air interface systems is very close to the alkane-water systems, and to pure alkane systems (shown in green open symbols in Fig. 6(Right)). 3.1.4. Order parameter result Calculating the order parameters of alkanes in the bulk, and on the interface regions of the alkane-air, alkane-water systems, the result in the temperature range of 340 K to 400 K is shown in Fig. 7(Left) and (Right). As can be seen, the bulk of alkane molecules packed almost randomly in both kinds of systems by showing very small order parameter
Fig. 6. Rg comparison for 9 alkanes in the alkane-water interface systems at 300 K (Left) in the bulk (black circles), on the upper interface (red squares), and on the lower interface (blue diamonds) with experimental data [29] and MC simulation data[30]; and for alkanes in the water-interface (in black) and air-interface (in red), and pure alkane (in green) systems at 340 K, 360 K, 380 K and 400 K (Right). Only the standard deviations of Rg in the bulk of water-interface at 300 K and pure alkane systems at different temperatures are shown to be clear. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 5
Fuel 267 (2020) 117252
G. Rucker, et al.
Fig. 7. Order parameter result comparison for alkanes in the alkane-air interface system (Left) and alkane-water system (Right) in the bulk (in black open symbols), on the interface layer (in solid symbols) (in blue, and orange) at four different temperatures: 340 K (shown by red circles), 360 K (by green squares), 380 K (by blue diamonds), and 400 K (by orange upside-triangles). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
result for n-alkanes at different chain lengths and at different temperatures. Similarly, the alkanes packed randomly in the pure alkane systems at the same temperature range with result shown in Figure S5. However, much larger order parameter result for n-C17, n-C20, n-C21, nC24 at 300 K was observed because of crystallization in both the pure alkane system (shown in Figure S5) and air/water-interface system. That agrees with the result from our previous work, in which n-C22 packed randomly in liquid phase [27] at high temperatures, but crystallized at room temperature. Similarly, we calculated the order parameters of alkane in the interface regions in both kinds of systems. The order parameters of alkanes on the interface are larger than the ones in the bulk; thus molecules in the interface region are more orderly packed. Comparing the result in Fig. 7(Left) and (Right), alkanes are orderly packed on the alkane-air interface and on the alkane-water interface similarly. 3.2. Surface tension result 3.2.1. Surface tension result on pure water box Applying Eq. (1) to calculate the surface tension of the water–air system, it is 52.48 dyn/cm, which is lower than the experimental data which is around 71.7 dyn/cm at 300 K [31]. However, it is within the surface tension range from the seven flexible water models tested by Yuet and Blankschtein [32].
Fig. 8. Surface tension comparison of n-alkanes-air interface systems at five different temperatures.
3.2.2. Surface tension result for n-alkane-air interfaces For a total of 10 ns sampling runs, the surface tension from each 2 ns was calculated. Based on that, the average surface tension and the standard deviations were calculated accordingly, with the surface tension and the standard deviation comparison shown in Fig. 8. With the temperature increasing, the surface tension of alkane-air interface decreases almost linearly. Also, at the same temperature with the alkane chain length increasing, the surface tension increases. Only modest standard deviations of surface tension at different temperatures and for different n-alkanes are observed. With the alkane chain length increasing and temperature increasing, the standard deviation increases slightly. Comparing the surface tension of n-C8 to n-C16 predicted with the experimental data from Jasper et al. who measured the surface tension of n-alkanes using the capillary rise method [2], and from Somayajulu [33] who predicted the surface tension of n-alkanes using three-parameter equation, results are shown in Fig. 9. A satisfying agreement was reached in terms of surface tension dependence on temperature and on
Fig. 9. Comparison of alkane-air interface surface tension predicted from simulations with the available experimental data at different temperatures. Only the surface tension of n-C8H18 to n-C17H36 at five different temperatures is shown, since they have experimental data available from ref[2] in solid symbol and from ref[34] in half-solid symbol. 6
Fuel 267 (2020) 117252
G. Rucker, et al.
chain length in the alkane-air and alkane-water systems at different temperatures was investigated, as the result shown in Figure S6 in the Supplemental Material. With the carbon number increasing, the size of n-alkanes increases as shown in Fig. 6, thus the alkane diffusion coefficient decreases. Generally, the higher the temperatures, the higher the diffusion coefficients. 3.3.2. Diffusion coefficient results for alkanes in the interface region Calculating the diffusion coefficients of n-alkanes on the interface and comparing them with those in the bulk, Fig. 12(Left) shows that alkanes diffuse faster in the alkane-air interface region than in the bulk, and Fig. 12(Right) shows that the alkanes in the alkane-water interface systems also can diffuse faster than in the bulk. Since the diffusional movement can be decomposed into the lateral direction (x and y directions) and z-direction, the diffusion coefficients of molecules in the upper and lower interface regions in the x-y plane and in the z-direction were calculated separately, and the result at different temperatures are shown in Fig. 13(Left) for the alkane-air system and Fig. 13(Right) for the alkane-water system. The original data are shown in Table S1 and Table S2 in the Supplemental Material. As can be seen, the molecules can diffuse much faster on the x-y plane (shown as black circles) than in the z-direction (shown in red-squares). In general, the diffusion of alkanes in the lateral direction is 2 to 10 times faster than in the z-direction in both the air-interface and waterinterface systems.
Fig. 10. Surface tension and standard deviation of 9 different alkane-water interfaces from 300 K to 400 K. The experimental data from ref [35] for n-C8, nC9, and n-C12 are shown in solid symbols and connected using solid lines; while simulation data are shown in open symbols and connected using dashed lines.
the carbon number. The absolute surface tension data from the simulation predictions at the temperature range of 300–400 K are slightly lower than the experimental data by around 2 dyn/cm for each alkane. 3.2.3. Surface tension result for alkane-water interfaces Similarly, the surface tension and standard deviation of alkane with water at different temperatures were calculated using Eq. (1). As the result shown in Fig. 10, the surface tension of alkane-water decreases almost linearly with the temperature increasing. With the alkane chain length increasing, the surface tension of alkane-water interface increases. Comparing the simulation prediction with experimental data from ref [34] (shown as solid symbols) in Fig. 10, the predicted surface tension of n-C8, n-C9, and n-C12 is very close to the experimental data, especially for n-C12. Comparing the surface tension of the alkane on the interface with air and water shown in Figs. 8 and 10, the surface tension of each system almost doubled after the air was replaced by water. However, the molecule-air interface has a stronger temperature dependence than the corresponding alkane-water interfaces. Comparing to the surface tension result of alkane-air interface, the alkane-water interface surface tension showed a larger standard deviation and a larger discrepancy between the simulation prediction and experimental data. The reason should come from the TIP3P water model, which showed a large deviation from the experimental data in this work as shown in Section 3.2.1.
3.4. Discussion The surface tension of 9 n-alkanes was calculated on the interface with air and with water. Based on Nicolas and Smit [15], predicting the correct density profiles is necessary to predict the correct surface tension. Based on the density profile results shown in Figs. 3 and 4 which look reasonable in comparison with available experimental data, it is expected that the surface tension prediction should be reasonable. The surface tensions of n-alkanes on the air-interface are close to experimental data as shown in Fig. 9, but consistently smaller by around 2 dyn/cm at different temperatures and for all six n-alkanes from n-C8 to n-C17. Similarly, in the alkane-water systems, the surface tension predicted is smaller than the available experimental data by 4–5 dyn/cm. The small deviation in surface tension for alkane-air and alkane-water proves that the CHARMM forcefield parameters are reasonable for the surface tension calculation, although the surface tension of water predicted at 300 K using the same input parameters and procedure is smaller than the experimental data by around 20 dyn/cm. Based on our work, the surface tension of water-alkane and air-alkane interfaces is sensitive to the alkane chain length and temperature. With the alkane chain length increasing, the surface tension of the interface increases; with the temperature increasing, the surface tension of the interface decreases. In this paper, the order parameter between n-alkane molecules was calculated by connecting the head carbon to the tail carbon (C1 to CN) to build the unit vector. Other researchers built the unit vector using the carbon atom pair of CK to CK+2 in the order parameter calculation [36]. In order to check the consistency of the result, different carbon atom pairs such as C2-C4, C3-C5, and C6-C8 in the water-n-C8 system were tried at different temperatures, which gave the order parameter result close to the result as shown in Fig. 7 using the carbon atom pair of C1C8. On the interface, which is defined as the alkanes within 6 Å of thickness to the air/water phase in this project, alkanes prefer to pack more orderly than in the bulk. Michael and Benhamin [37] found that alkanes on the alkane-water interface prefer to pack laterally with the interface based on the order parameter result in lateral and z-directions. Because of the shape of n-alkanes, such kind of packing can explain why alkanes move faster in the lateral direction than in the z-direction in both kinds of interface systems found out in this work.
3.3. Diffusion coefficient result 3.3.1. Bulk diffusion coefficient result at different temperatures The bulk diffusion coefficients of alkanes in both the air-alkane systems and water-alkane systems were calculated which include all the alkanes molecules in the system, and result was compared with pure nalkane systems at different temperatures, and with available self-diffusion coefficient data from experimental work [35]. In order to show the result clearly, the even-carbon-numbered alkane results are shown in Fig. 11(Left), while the odd-carbon-numbered alkanes in Fig. 11(Right). In general, n-alkanes diffuse slower in the water-alkane systems than in the air-alkane systems, but alkanes diffuse much faster in both kinds of interface systems than in the pure n-alkane systems. The logD of alkanes has an almost linear relationship with 1/T in the interface systems and in the pure alkane systems. The simulation data on n-C8, n-C9, n-C12, n-C13 and n-C16 at 300 K are very close to experimental data at 303 K using the MRI method[36] (shown as patterned symbols). That proves that our diffusion coefficients predicted using the simulation method are reasonable. Similarly, the dependence of the diffusion coefficient on the alkane 7
Fuel 267 (2020) 117252
G. Rucker, et al.
Fig. 11. Bulk diffusion coefficient of n-C8, n-C12, n-C16, n-C20 and n-C24 (Left) and n-C9, n-C13, n-C17, and n-C21 (Right) in the air-interface systems (in empty symbols connected with solid line), in the water-interface systems (in solid symbols connected with dashed line) and in the pure alkane systems (in half-solid symbols connected with dot-dashed line) in comparison with experimental data [36] (in patterned symbols). Data at 300 K for n-C17, n-C20, n-C21, and n-C24 are not shown because of crystallization in pure-alkane, air-interface and water-interface systems.
Fig. 12. Diffusion coefficient of alkanes in the interface region (by squares) compared with in the bulk (by circles) for the alkane-air (Left) and alkane-water (Right) interface systems.
Fig. 13. Diffusion coefficients of 9 n-alkanes in the lateral direction and in the z-direction at five different temperatures in the air-alkane (Left) and water-alkane (Right) interface systems.
8
Fuel 267 (2020) 117252
G. Rucker, et al.
Up to now, some researchers have worked on the liquid–vapor interface and found the enhanced surface reactions in droplets and the air–water interface [38,39]. Liu et al. [40] calculated the diffusion tensor of water in confined fluids and interface and found that the lateral diffusion tensor of water was around 4 times of that in the bulk, and the normal diffusion tensor was around 2 times of that in the bulk. In this work, the diffusion coefficients of alkane at the interface region with air and with water were calculated. The diffusion coefficient was enhanced in both normal and lateral directions in both kinds of interface systems. More specifically, the lateral direction diffusion rate was around 2 ~ 5 times of that in the bulk, while the normal direction increased moderately. When comparing the diffusion coefficients of n-alkanes in the bulk of air/water-interfaces and in the pure alkane systems, faster diffusion coefficients in the system having interfaces were observed as shown in Fig. 11. That emphasized the effect of the interface interaction to the dynamic properties of alkanes. Comparing our diffusion coefficients predicted with experimental result using the MRI method [36], very consistent agreement was found. Comparing the radius of gyration result and radial distribution function result for both kinds of interface systems with available experimental data, consistent agreements were reached for both individually. With the n-alkane chain length increasing, the size of n-alkane increases, which has negligible dependence on temperature and on the position of alkanes in the interface systems.
resources were provided through the high-performance computers at TTU and resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.fuel.2020.117252. References [1] Roberts FL, Kandhal PS, Brown ER, Lee D-Y, Kennedy TW. Hot Mix Asphalt Materials, Mixture Design, and Construction, 2nd; National Asphalt Pavement Association: Lanham, MD; 1996, Chapter2. [2] Jasper JJ, Kerr ER, Gregorich F. The orthobaric surface tensions and thermodynamic properties of the liquid surfaces of the n-aIkanes, C5 to C28. JACS 1953;75:5252–4. [3] Ismail A, Tsige M, Veld in PJ, Grest G. Surface tension of normal and branched alkanes. Mol Phys 2007;105(23):3155–63. [4] Zhang Y, Feller SE, Brooks BR, Pastor RW. Computer simulation of liquid/liquid interfaces. I. Theory and application to octane/Water. J Chem Phys 1995;103:10252. https://doi.org/10.1063/1.469927. [5] Morrow NR. Interfacial Phenomena in Petroleum Recovery. New York: Dekker; 1991. p. 1. [6] Thomas JL, Roeselová M, Dang LX, Tobias DJ. Molecular dynamics simulations of the solution−air interface of aqueous sodium nitrate. J Phys Chem A 2007;111(16):3091–8. [7] Nguyen CV, Phan CM, Ang HM, Nakahara H, Shibata O, Moroi Y. Molecular dynamics investigation on adsorption layer of alcohols at the air/brine interface. Langmuir 2015;31(1):50–6. [8] Yuet PK, Blankschtein D. Molecular dynamics simulation study of water surfaces: comparison of flexible water models. J Phys Chem B 2010;114(43):13786–95. [9] Zhang L, Liu Z, Ren T, Wu P, Shen J-W, Zhang W, et al. Understanding the structure of hydrophobic surfactants at the air/ water interface from molecular level. Langmuir 2014;30:13815–22. [10] Chong L, Lai Y, Gray M, Soong Y, Shi F, Duan Y. Effects of frothers and oil at saltwater–air interfaces for oil separation: molecular dynamics simulations and experimental measurements. J Phys Chem B 2017;121(27):6699–707. [11] Alejandrea J, Tildesley DJ. Molecular dynamics simulation of the orthobaric densities and surface tension of water. J Chem Phys 1995;102(11):4574–83. [12] Carpenter IL, Hehre WJ. A molecular dynamics study of the hexane/water interface. J Phys Chem 1990;94:531–6. [13] van Buuren AR, Marrink SJ, Berendsen HJC. A molecular dynamics study of the decane/water interface. J Phys Chem 1993;97:9206–12. [14] Rivera JL, McCabe C, Cummings PT. Molecular simulations of liquid-liquid interfacial properties: water–n-alkane and water-methanol–n-alkane systems. Phys Rev E 2003;67. 011603. [15] Nicolas JP, Smit B. Molecular dynamics simulations of the surface tension of nhexane, n-decane and n-hexadecane. Mol Phys 2002;100(15):2471–5. [16] Vanommeslaeghe K, Hatcher E, Acharya C, Kundu S, Zhong S, Shim J, et al. CHARMM general force field: a force field for drug-like molecules compatible with the CHARMM all-atom additive biological force field. J Comput Chem 2010;31:671. [17] Yu W, He X, Vanommeslaeghe K, MacKerell Jr. AD. Extension of the CHARMM general force field to sulfonyl-containing compounds and its utility in biomolecular simulations. J Comput Chem 2012;33:2451. [18] Vanommeslaeghe K, MacKerell Jr. AD. Automation of the CHARMM general force field (CGenFF) I: bond perception and atom typing. J Chem Inf Model 2012;52:3144. [19] Vanommeslaeghe K, Raman EP, MacKerell Jr. AD. Automation of the CHARMM general force field (CGenFF) II: assignment of bonded parameters and partial atomic charges. J Chem Inf Model 2012;52:3155. [20] Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, et al. Scalable molecular dynamics with NAMD. J Comput Chem 2005;26:1781. [21] Ryckaert J-P, Ciccotti G, Berendsen HJC. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J Comput Phys 1977;23:327. [22] Hill TL. Introduction to Statistical Mechanics. New York: Dover; 1986. p. 314. [23] Ono S, Kondo S. Molecular Theory of Surface Tension in Liquids. Berlin: Springer; 1960. p. 134. [24] Frenkel D, Smit B. Understanding Molecular Simulation. 2nd San Diego: Academic; 2002. [25] Giorgino T. Computing 1-D atomic densities in macromolecular simulations: the density profile tool for VMD. Comput Phys Commun 2014;185(1):317–22. [26] Cecchini M, Rao F, Seeber M, Caflisch A. Replica exchange molecular dynamics simulations of amyloid peptide aggregation. J Chem Phys 2004;121:10748. [27] Sonibare K, Rathnayaka L, Zhang L. Comparison of CHARMM and OPLS-aa forcefield predictions for components in one model asphalt mixture. Constr Build Mater 2020;236. 117577. [28] Griesbaum K, Behr A, Biedenkapp D, Voges H-W, Garbe D, Paetz C, et al. “Hydrocarbons” in Ullmann's Encyclopedia of Industrial Chemistry. Weinheim: Wiley-VCH; 2005.
4. Conclusion This paper works on 9 n-alkanes with the air and with the water interface systems at five different temperatures, ranging from 300 K to 400 K. The density profiles, surface tension, diffusion coefficients, radius of gyration, radial distribution function, and order parameters of nalkanes in the bulk and at the interface region in both kinds of interface systems were calculated at different temperatures. It was found that the density profiles of n-alkanes can agree with available experimental data at the room temperature. The higher the chain length of n-alkanes, the higher the surface tension in both kinds of interface systems. The wateralkane interface usually has a higher surface tension than the air-alkane interface. At the same temperature, the alkanes can diffuse faster in the air-alkane interface system than in the water-alkane interface system, and faster in the interface system than in the pure-alkane system. Focusing on the alkanes on the interface in both kinds of systems, molecules are more orderly packed on the interface than in the bulk, and they also diffuse much faster than in the bulk. Decomposing the diffusion coefficients of alkanes into the lateral direction and the z-direction, alkanes on the interface can diffuse much faster in the lateral direction than in the z-direction. CRediT authorship contribution statement George Rucker: Data curation. Xiong Yu: Conceptualization. Liqun Zhang: Conceptualization, Methodology, Software, Investigation. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements We acknowledge TTU student Yilun Lee for helping with the initial setting up and analyzing some simulation trajectories. This work was supported by the Center for Energy Systems Research at TTU and by the Faculty Research Fund TTU provided to Dr. Zhang. The computer 9
Fuel 267 (2020) 117252
G. Rucker, et al. [29] Goodsaid-Zalduondo F, Engelman DM. Conformation of liquid N-alkanes. Biophys J 1981;35(3):587–94. [30] Bessières D, Piñeiro MM, De Ferron G, Plantier F. Analysis of the orientational order effect on n-alkanes: evidences on experimental response functions and description using Monte Carlo molecular simulation. 074507 J Chem Phys 2010;133. https:// doi.org/10.1063/1.3472283. [31] Vargaftik NB, Volkov BN, Voljak LD. International tables of the surface tension of water. J Phys Chem Ref Data 1983;12:817–20. [32] Yuet PK, Blankschtein D. Molecular dynamics simulation study of water surfaces: comparison of flexible water models. J Phys Chem B 2010;114:13786–95. [33] Somayajulu GR. A generalized equation for surface tension from the triple point to the critical point. Int J Thermophys 1988;9(4):559–66. [34] Zeppieri S, Rodrı́guez J, López de Ramos AL. Interfacial tension of alkane + water systems. J Chem Eng Data 2001;46:1086–8. [35] Tofts PS, Lloyd D, Clark CA, Barker GJ, Parker GJM, McConville P, et al. Test liquids
[36] [37] [38] [39] [40]
10
for quantitative MRI measurements of self-diffusion coefficient in vivo. Magn Reson Med 2000;43:368–74. Harris JG. Liquid-vapor interfaces of alkane oligomers: structure and thermodynamics from molecular dynamics simulations of chemically realistic models. J Phys Chem 1992;96:5077–86. Michael D, Benjamin I. Solute orientational dynamics and surface roughness of water/hydrocarbon interfaces. J Phys Chem 1995;99:1530–6. Liu P, Harder E, Berne BJ. On the calculation of diffusion coefficients in confined fluids and interfaces with an application to the liquid−vapor interface of water. J Phys Chem B 2004;108(21):6595–602. Valsaraj KT. Trace gas adsorption thermodynamics at the air−water interface: Implications in atmospheric chemistry. Pure Appl Chem 2009;81(10):1889–901. Nissenson P, Knox CJH, Finlayson-Pitts BJ, Phillips LF, Dabdub D. Enhanced photolysis in aerosols: evidence for important surface effects. Phys Chem Chem Phys 2006;8(40):4700–10.