Abstract
We combine two approaches to describe an activated carbon which contains micro- and mesopores. We propose a simple approach of packing individual three-dimensional rigid structures to represent microporous region of commercially available BAM-P109 activated carbon. Mesopores are represented by planar walls interacting via the standard Steele 10-4-3 potential assuming experimental pore-size distribution. The models are used to predict n-perfluorohexane adsorption at 273 K temperature. We used uncharged united atom model for n-perfluorohexane molecule. The non-bonded interactions were described by the Gordon potential model, with parameters adjusted to reproduce the Mie force field potential, which was designed to accurately represent vapour pressure of alkanes and perfluoroalkanes. Adsorption isotherm was calculated using Grand Canonical Monte Carlo method implemented in MCCCS Towhee. We show that predicted n-perfluorohexane loading at relative pressures of interest are significantly lower than the experimental ones. The underprediction is attributed to the low aromatic carbon interaction strength and sensitive calculations of chemical potential. Adjusted interaction parameter of aromatic carbon and chemical potential led to n-perfluorohexane adsorption that agreed well with the reported experimental values.
Introduction
Activated carbon is one of the most usable materials in modern industry, spanning its applications from storage and purification to catalysis. Activated carbons are relatively inexpensive, have high adsorption capacity, and their properties can be tailored for specific applications. The demand for these materials is increasing, and the ability to predict and tune its performance using virtual techniques would be essential for cost- and time-effective product development. The aim of this work is to predict the adsorption of n-perfluorohexane (PFH) on activated carbon BAM-P109, using grand canonical Monte Carlo simulations and force fields reported in the literature. This work is an entry for 8th Industrial Fluid Properties Simulation Challenge (IFPSC, 2014).
PFH is a hexane derivative where all hydrogen atoms are substituted with fluorine atoms. PFH is used in many industrial and medical applications due to its unique properties such as chemical and thermal stability, biological inertness and ability to dissolve gases from air, just to name a few examples. Until the previous 7th Industrial Fluid Properties Simulation Challenge (Ross et al., 2014) there were no studies on PFH adsorption in porous media. However, PFH is considered to be moderately challenging molecule to model and recognized as potentially one of the strong global warming compounds (Ross et al., 2014). Therefore, it is important to assess currently available computational techniques and force fields for accurate predictions of PFH adsorption.
In this work, we use a simple approach towards representing activated carbons as a packed system of individual molecules (Gonciaruk and Siperstein, 2015). These molecules consist of graphene-like flat arms connected through rigid core. The core prevents molecules from efficient packing leading to porous material. By controlling graphene arm size it is possible to adjust the structural properties of the material, such as accessible surface area and pore-size distribution (PSD). We have also used idealized slit pores where the Steele potential represents interactions, and the main differences between the two approaches are discussed.
Simulation details
Perfluorohexane
The PFH molecules were represented using a united atom approximation with no explicit charges. In a PFH molecule, fluorine atoms are combined with the neighbouring carbon atom to form a single interaction site. The non-bonded and bonded interactions parameters were based on the force field developed by Potoff and Bernard-Brunel which was reported to accurately describe perfluoroalkanes, alkanes, alkenes and their mixtures (Potoff and Bernard-Brunel, 2009; Potoff and Kamath, 2014). In this force field, site–site interactions are described via the Mie potential (a generalized form of Lennard-Jones 12-6 potential), given by:
In this work, simulations are carried out using the MCCCS Towhee v7.0.6, which is a powerful all-purpose Monte Carlo simulation software (Martin, 2013). Unfortunately, the Mie potential is not available in Towhee v.7.0.6. Therefore, we have approximated the Mie force field functional form by using the so-called Gordon n–6 potential (Gordon, 2006) which is implemented in Towhee, and is given by:
Gordon and Mie potential parameters for perfluorohexane units.

Mie and Gordon potential curves for (a) CF2 and (b) CF3 units.
Following Potoff's work, CF x units were separated by a fixed bond length of 1.54 Å. The bond angle was described by standard harmonic oscillation, and OPLS fluorocarbon four-term potential was used to describe torsional conformations (Potoff and Bernard-Brunel, 2009).
Microporous carbon
Model carbons were constructed using Material Studio software (Accelrys, Inc.) using the method inspired by (Gonciaruk and Siperstein, 2015). Graphene arms were created in planar form by connecting six-membered carbon rings. The edge carbon atoms of the graphene were saturated with hydrogen atoms. Three different arms were constructed: small and large arms of disk shape and a small-size ribbon-like arm (Figure 2). The arms were then connected through two different centres inspired by the structures of [2.2.2]-propellane (trip) and octahedral propellane (octa) (Figure 3). The trip and octa structures possess a rigid three-fold symmetry that keeps graphene arms separated in three-dimensional space. The carbon model is described fully atomistically. Although it has been reported that BAM-P109 contains trace amounts of sulphur, sodium and oxygen as well as some other elements (IFPSC, 2014), they were not included into the model, but it would be useful to understand the effect of functional group presence to the overall PFH adsorption in the future. For packing procedure, interactions between atoms are described using the Dreiding force field (Mayo et al., 1990). Partial charges are calculated using charge equilibration method QEq implemented in Material Studio (Rappe and Goddard, 1991). The Lennard-Jones 12-6 potential was used to model the van der Waals interactions, while long-range electrostatic interactions were calculated using the Ewald method.
Representation of the graphene arms used to model the carbon adsorbent materials. (Left) Large and (middle) small disk-like arms, and (right) ribbon-like arm. Cores of model carbons.

Composition of the carbon model system.
Gordon potential parameters for carbon model.
The model carbons were characterised and compared in terms of density, surface area, pore volume and PSD. Geometric surface area is defined by a line that the centre of the probe draws while rolling along the van der Waals surface of the adsorbent. The probe was selected with a 3.76 Å diameter which corresponds to van der Waals diameter of argon (Ar).
Structural properties of carbon models.
The density of the model carbons is simply calculated by dividing the total mass of the system by the total volume of the simulation box. However, it is thought that the skeletal volume should be used instead for more correct comparison with experimental data (Abbott and Colina, 2011). Skeletal densities are calculated from the following equation:
Two carbon structures were created: one that reproduces correct surface area accessible to Ar atom and another that has a free volume and PSD similar to the one measured experimentally (Tan and Gubbins, 1992) (Table 4). The PSD in the micropore region is shown in Figure 4. It is obvious that our models are not capturing the region of very small pores. This may result in slightly lower Ar loadings at low pressures. In the case of larger PFH molecules, adsorption can be expected to occur only in the pores wider than about 0.8 nm from where calculated PSD of first model follows very similar trend to experimental PSD. Nevertheless, we expect that it is possible to obtain model carbons which would represent all experimental structural properties by changing the size of the arms and combination of structures with different arm sizes, shapes or cores. Packed model carbons with similar arm size possess similar structural properties (PSD, surface area and free volume) regardless of arm shape (Gonciaruk and Siperstein, 2015). A slight trend was observed with respect to disc-shaped arms’ size: larger arms tend stack better than small arms, reducing the surface area, while smaller arms can move more freely, penetrating ‘internal molecular free volume’ (Long and Swager, 2001) created between two arms due to three-fold structure of triptycene. Therefore, the large graphene arms were introduced to create the overall network and small arm structures were included to reduce graphene stacking, while ribbons were supposed to penetrate the larger free pores and increase microporosity.
Pore-size distribution for the two different carbon models compared with the experimental data taken from Benchmark data provided by (IFPSC, 2014).
Slit pore model
Adsorption of PFH was also modelled using a simple slit pore model. A collection of slit pores with separations narrower than 2 nm was used to represent microporous region of activated carbon. To account for mesoporosity, adsorption was calculated in a slit pores of 13, 16 and 20 nm wide which are around the local maxima for mesoporous region of PSD as can be seen in the experimental data provided in (IFPSC, 2014).
Solid–fluid interaction parameters for equation (8).
Vapour pressure calculation
The vapour pressure for PFH at high temperatures was obtained using the direct coexistence method (Chapela et al., 1977) using Lammps (Plimpton, 1995). A rectangular box with dimensions L
x
= L
y
= 4.786 nm and L
z
= 23.8 nm with total of 604 PFH molecules was used, and periodic boundary conditions in all directions were considered. The equations of motion were with a time step of 2fs, and a Nose-Hoover thermostat was used to control temperature. The system was equilibrated for ns and averages were obtained over 10 ns. A typical two-phase coexistence simulation is shown in Figure 5. The coexistence densities were obtained using the density profiles along the larger dimension used (z). The vapour pressure was calculated by using an independent simulation at the vapour density calculated during the direct coexistence simulations. Results from the simulations are shown in Figure 6. Simulations were carried out at high temperatures, and the saturation pressure at 273 K was obtained by extrapolation using the Clausius-Clapeyron equation. It was estimated that at 273 K, the saturation pressure for PFH is 0.144 bar, which is slightly higher than the values reported by National Institute of Standards and Technology (NIST) (Crowder et al., 1967; Lemmon et al.).
Representative configuration of vapour–liquid equilibrium at 300 K obtained using the direct-coexistence method. Molecular dynamics results for (a) vapour–liquid coexistence and (b) saturation pressure for PFH. Symbols are values calculated in this work and solid curves are data taken from (Potoff and Bernard-Brunel, 2009).

Adsorption of argon
All adsorption simulations were performed using MCCCS Towhee v7.0.6 (Martin, 2013) in the grand canonical ensemble. Adsorption of Ar was simulated using universal force field (Rappe et al., 1992). Interaction parameters between unlike atoms were calculated using Lorentz-Berthelot rule. In the simulation, a combination of adsorbate insertion, deletion and translation steps were used. Each simulation was typically carried out for minimum of 10 × 106 Monte Carlo moves, including equilibration and production steps. Framework atoms were frozen in space during adsorption simulation. Ar adsorption was simulated in both carbon models (Model 1 and Model 2). The thermodynamic inputs were temperature of 87.45 K, simulation unit cell volume and chemical potentials corresponding to target relative pressures between 0.2 and 0.7. The relative pressure is defined as the relative to the bulk saturation pressure for 87.45 K temperature and was calculated using Antoine equation with parameters taken from NIST (Drii and Rabinovich, 1966; Lemmon et al.). Units of adsorption data are provided in number of molecules per unit cell which is then converted to volume per mass of adsorbent in order to compare calculated and experimental results. For the conversion, molar volume (22.4 L mol−1) of an ideal gas at standard temperature and pressure (STP) was used.
Monte Carlo simulations in the NPT ensemble for a system containing only Ar molecules were used to establish correlation between chemical potential and pressure of interest.
Adsorption of perfluorohexane
Adsorption PFH was simulated using Gordon force field (Gordon, 2006), with the adjusted parameters as described previously. A combination of adsorbate insertion, deletion and motion (translation, rotation and confirmation change using configurational-bias) steps were allowed. Typically a minimum total of 6 × 106 Monte Carlo moves were used for equilibration and production steps in the model disordered carbon. Shorter simulations were run for the slit pore systems. Adsorption was carried out at the fixed simulation cell volume at temperature of 273 K temperature and constant chemical potentials corresponding to target relative pressures of 0.1, 0.3 and 0.6, taking reported saturation pressure value of 0.0855 bar (IFPSC, 2014). As described earlier, units of PFH loading in carbon models are also converted using 22.4 L mol−1 ideal gas molar volume as a requirement for the challenge.
Total amount of PFH adsorbed n(P) in the slit pores of sizes between 0.5 and 2.2 nm was calculated assuming that the slit pore have the same PSD as carbon model:
Results
Argon adsorption
Calculated Ar adsorption at 87.45 K is presented in Figure 7 along with experimental benchmark isotherm (Tan and Gubbins, 1992). Ar adsorption in Model 1 is overestimated; however, the disagreement can be attributed to higher surface area of the system comparing to experimental one. The ratio of surface areas is the same as the ratio between adsorbed Ar volumes in simulated system and experimental sample. Therefore, we created the system with the satisfactory surface area (Model 2), which is very close to the experimental surface area. In this case, we were able to obtain the same Ar loading in the system as in the experimental sample.
Argon adsorption in experimental sample (solid line), in carbon model with correct surface area (open circles) and carbon model with similar free argon volume and pore-size distribution (open diamonds).
Increased slope of experimental isotherm at relative pressure above 0.9 is attributed to multilayer formation in the large pores. However, mesoporosity of the sample is not considered for Ar adsorption calculations. Therefore, carbon model does not exhibit increased Ar adsorption at pressures close to Ar saturation pressure. Slightly lower PFH loading at 0.9 relative pressure is due to the short simulation runs and is expected to level off after longer simulation.
Perfluorohexane adsorption
Predicted adsorption of PFH in carbon Model 1 at 273 K is given in Figure 8 and Table 6. Loadings in mmol g−1 units are additionally given in Table 6. Standard deviation was computed directly in Towhee over the five blocks that divide the simulation. From the isotherm figure it is seen that PFH loading increase is insignificant between 0.06 and 0.1 relative pressures. Between 0.1 and 0.6 relative pressures, the model steadily adsorbs PFH and saturation is not expected. A snapshot of PFH molecules adsorbed in carbon model can be seen in Figure 9. It shows that PFH molecules tend to align parallel on the surface of graphene arms.
Total amount of perfluorohexane adsorbed in carbon model at 273 K. Adsorbed amount of PFH in carbon model at 273 K. Representative in-cell image of perfluorohexane in carbon model. Colour code: carbon atoms of carbon model are dark grey and units of perfluorohexane are pink (light grey).

Adsorption in slit pores shows clear layering formation as shown in Figure 10. All adsorbed molecules preferably adsorb parallel to the pore walls and do not extend perpendicular between two surfaces. No adsorption was observed for pores smaller than 0.7 nm (taking the pore size as the distance between the smooth solid walls).
Representative configurations showing the adsorption of PFH in microporous slit pores at P/P0 = 0.1 and 273 K.
Although no significant difference was observed for the amount adsorbed in the 16 nm slit pore at the pressures of interest, at a pressure of P/P0 = 0.8 condensation in the whole pore was observed. Therefore, the contribution to the amount adsorbed of the mesopore region for the range of pressures of interest is negligible, but an increase in about 1.8 mmol g−1 (or 40 STP cm3 g−1) can be expected at high pressures.
Amount of PFH adsorbed in the microporous slit pore model is 0.74 mmol g−1 (16.5 STP cm3 g−1) at 0.1 relative pressure, which overestimates the adsorption compared to PFH loading in carbon model. It is obvious that at even very low relative pressure PFH molecules start to align parallel to each other forming an ordered arrangement, as shown in Figure 11(a). In the carbon models with irregular complex pore structure, it is more difficult for PFH molecules to form ordered structures; therefore, less PFH is adsorbed. However, it is possible that both models can reach the same PFH loading at higher pressures. It is interesting to note that in the slit pores of the sizes that commensurate with the length of PFH molecule, some of the molecules form bridges between the walls (Figure 11(b)). This suggests that packing of the molecules can be significantly higher than what would be expected from the formation of parallel layers to the surface.
Alignment of perfluorohexane molecules in the slit pore of (a) 0.9 nm and (b) 1.5 nm width at P/P0 = 0.6 and 273 K.
Comparing the results with the experimental data (Ross et al., 2016) it is obvious that our calculations significantly underpredict PFH adsorption, although the slope of the isotherms is very similar (Figure 12). We revised our work and believe that the proposed carbon model is still valid and would be beneficial for further adsorption studies as it is simple to construct and maintain the geometrical complexity of the pores. It is evident from visualisation of the simulation cell with adsorbed PFH that there is still space for additional molecules (Figure 9). Also, at higher chemical potential (pressure) the system was able to accommodate greater amount of PFH, suggesting that the selected interaction parameters or incorrect chemical potential causes such a weak adsorption. By careful examination of the implemented approach, we found several points that could affect the underprediction:
Solid–fluid interaction parameters. The value of epsilon (ε) used for aromatic carbon was different when simulating argon and PFH adsorption. The PFH adsorption simulations used the ε for aromatic carbons in the Steele wall. However, for Ar adsorption, the UFF force field was used where ε is 52.84 K for aromatic carbon. The ε from the Steele potential is 1.75 times smaller than UFF value (Table 3). Since the calculated Ar adsorption match the experimental isotherm well, it is reasonable to assume that larger value assigned by UFF force field describes aromatic carbon-adsorbate interactions better. Charges are not used in the calculations. Knowing that fluorine atoms are highly electronegative, electrostatic interaction might be significant, but whether this would lead to higher adsorption is questionable since the negative ‘cloud’ created by fluorine atoms on each molecule may result in stronger repulsion between the adsorbate molecules. However, long-range interactions as well as their magnitude are already incorporated in the parameters reported by Potoff and Bernard-Brunel (Potoff and Bernard-Brunel, 2009). On the other hand, interactions between electron-rich atom and electron-deficient aromatic ring have been recognized in many biomolecules (Mooibroek et al., 2008). Consideration of lone pair–π interactions may help to better understand and represent adsorption of fluorinated compounds on activated carbons. The default Martin and Frischknecht 2006 algorithm in Towhee (Martin and Frischknecht, 2006) was used for the chemical potential calculation. The Widom insertion method (Widom, 1963, 1982) is recommended to ensure correct chemical potential calculation as it performs only one trial for any non-bonded interaction selection step and not 10 trials. Unfortunately, due to limited time, we did not establish an accurate pressure-chemical potential correlation. Experimental and previous entry results for 8th Fluid Challenge of perfluorohexane adsorption isotherms at 273 K.

We performed additional PFH adsorption simulations at larger chemical potentials and are using the same ε as the one used for the Ar adsorption isotherms. The results are summarised in Figure 13. The increase in ε resulted in a larger amount adsorbed of PFH in proposed model carbons. Nevertheless, adsorbed PFH adsorbed volume is too low compared to experimental values. Simulations at higher chemical potential showed larger PFH adsorption that are closer to the experimental values suggesting that initially calculated chemical potential does not correspond to pressures of interest. Adsorption in mesoporous region represented by the slit pore model at the same increased chemical potential can additionally contribute about 10 to 11 STP cm3 g−1 of adsorbed PFH.
Perfluorohexane adsorption in model carbon as a function of chemical potential.
Conclusions
Perfluorohexane adsorption was predicted in commercially available BAM-P109 activated carbon at 273 K and relative pressures of 0.1, 0.3 and 0.6. To represent microporous region, activated carbon model was constructed as a collection of randomly packed individual molecules that possess three-dimensional rigid structures. Prior to PFH adsorption, the model system was validated by calculating Ar isotherm and comparing the results with the experimental benchmark data. The model was able to accommodate the same volume of Ar as in experiment up to 0.9 relative pressures. PFH adsorption in carbon model was modelled using Gordon force field, and non-bonded interactions were described using modified Gordon n-6 potential, with the parameters adapted from the work of Potoff and Bernard-Brunel (2009). We found that the virtual model of activated carbon can adsorb 4.1 ± 0.6 to 11.8 ± 1.3 cm3 g−1 of PFH in the range of relative pressures between 0.1 and 0.6, and at higher pressures when the mesoporous region is filled, the amount adsorbed can reach over 40 cm3 g−1. The calculated amounts of PFH loadings in model carbons are significantly underpredicted. Errors in the estimation of the chemical potential are considered to be the main problem in the estimation of the experimental adsorption isotherm.
Footnotes
Acknowledgments
The authors would like to acknowledge the assistance given by IT Services and the use of the Computational Shared Facility at The University of Manchester. AG is grateful for the postgraduate scholarship from the School of Chemical Engineering and Analytical Sciences from the University of Manchester.
Declaration of Conflicting Interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) received no financial support for the research, authorship, and/or publication of this article.
