Molecular mechanics and quantum mechanical modeling of hexane soot structure and interactions with pyrene

Abstract

Molecular simulations (energy minimizations and molecular dynamics) of an n-hexane soot model developed by Smith and co-workers (M. S. Akhter, A. R. Chughtai and D. M. Smith, Appl. Spectrosc., 1985, 39, 143; ref. 1) were performed. The MM+ (N. L. Allinger, J. Am. Chem. Soc., 1977, 395, 157; ref. 2) and COMPASS (H. Sun, J. Phys. Chem., 1998, 102, 7338; ref. 3) force fields were tested for their ability to produce realistic soot nanoparticle structure. The interaction of pyrene with the model soot was simulated. Quantum mechanical calculations on smaller soot fragments were carried out. Starting from an initial 2D structure, energy minimizations are not able to produce the observed layering within soot with either force field. Results of molecular dynamics simulations indicate that the COMPASS force field does a reasonably accurate job of reproducing observations of soot structure. Increasing the system size from a 683 to a 2732 atom soot model does not have a significant effect on predicted structures. Neither does the addition of water molecules surrounding the soot model. Pyrene fits within the soot structure without disrupting the interlayer spacing. Polycyclic aromatic hydrocarbons (PAH), such as pyrene, may strongly partition into soot and have slow desorption kinetics because the PAH-soot bonding is similar to soot–soot interactions. Diffusion of PAH into soot micropores may allow the PAH to be irreversibly adsorbed and sequestered so that they partition slowly back into an aqueous phase causing dis-equilibrium between soil organic matter and porewater.

Introduction

The importance of black carbon in the environment has been recognized for some time (see ref. 4 and references therein). Soot is one type of black carbon generated from combustion of fossil fuels. When soot particles form an aerosol, soot surfaces can act as a heterogeneous catalyst in a number of environmentally important atmospheric reactions. [5–7] When deposited in soils and sediments, soot has a significant influence on the transport and bioavailability of some organic contaminants due to their strong partitioning into soot.[8] In addition, organic contaminants, such as polycyclic aromatic hydrocarbons (PAH), can form simultaneously with soot in the combustion process.[9] Therefore, trace amounts of soot can play a disproportionate role in the long-term sequestration of otherwise biodegradable compounds.[10] In addition, soot may represent the more hydrophobic end-member in the distributed reactivity model for soil organic matter of Weber and co-workers[11] or the glassy/rubbery model of Pignatello and co-workers.[12] Both of these models suggest a component of soil organic matter that has a high affinity for hydrophobic organic molecules with a rigid structure that can support voids large enough for molecules to diffuse into. Soot has these qualities although it certainly is not the only component of soil organic matter that may play a role in the long-term adsorption of organic contaminants.

Although the role of soot in the environment and the analytical chemistry of soot have been studied extensively, molecular simulations of soot structure and chemistry are rare. Some studies[13] have modeled activated carbon, but the composition and structure of activated carbon are distinct from those of soot (see ref. 7 and references therein). The goal of this study is to discover if molecular modeling techniques with an appropriate force field will produce structural results consistent with observations of soot structure. If so, molecular modeling may provide a method for understanding the strong partitioning and irreversible adsorption of PAH into soot. (Note: 'irreversible' is defined as a process that does not undergo the back reaction via the same pathway as the forward reaction. This does not mean that PAH can not be desorbed once they are sorbed. For instance, if the sorbate changes the structure of the sorbent, then the desorption reaction will likely occur through a different mechanism and the sorption is 'irreversible'.)

The basic model used in this study was taken from the work of Smith and co-workers.[1] This structure was generated from experimental studies on hexane soot.[14] The model was derived from compositional analyses, spectroscopic analyses (e.g., FTIR, Raman, EPR and ¹³C CP/MAS NMR) and adsorption experiments.[5] Consequently, the initial composition and structure used in this study should be an excellent approximation for hexane soot. The composition of this model is C₄₈₉H₁₆₀O₃₄ (683 atoms). The structure used in this study is not meant to be a model for all soots, carbon black or char. Soots with different origins are expected to have different compositions and structures,[15] so the results presented in this paper may not be applicable for a wide variety of substances commonly known as black carbon. The purpose of this paper is to determine which computational methods give reasonable results regarding soot structure and behavior. However, many of the physicochemical principles modeled herein may operate to a greater or lesser degree in various types of black carbon.

Methods

Molecular mechanics simulations were performed with two different force fields in two different molecular modeling programs. The MM+ force field[2] was used within the program HyperChem 5.1;[16] and the COMPASS force field[3] was used with the program Cerius[2] (see ref. 17). A major difference between these force fields is that the parameters within COMPASS[3] were optimized to fit data of condensed phases; whereas MM+ parameters were optimized against gas-phase molecules. Thus, COMPASS should be more likely to give accurate results for a soot particle. Both sets of force fields were used within energy minimizations and molecular dynamics (MD) simulations.[18] The main difference between these two methods is that energy minimizations occur at 0 K and result in a local energy minimum structure; whereas MD simulations include temperature and time and allow for a larger sampling of configuration space. This difference is critical when examining possible structures of large, complex molecules.[19] The soot model was first subjected to energy minimization from the 2D structure,[1] then MD simulations were carried out using the minimum energy structure as the starting configuration. MD simulations were performed with a time step of 1 × 10^-15 s at 300 K for 10 000 to 100 000 time steps.

The 683 atom soot model represents a small fragment of a soot particle which can be an aggregate 0.1 to 1 μm in size made up of smaller particles on the order of 30 nm in diameter. Hence, the effect of model system size on the calculated results was investigated. A larger soot model was constructed from the energy-minimized structure of four 683 atom soot molecules. These models were randomly placed in proximity to one another and subjected to energy re-minimization and MD simulations with the COMPASS force field in Cerius[2] (see ref. 17). The final diameter of this 2732 atom soot model is approximately 5 nm, closer to the actual soot particle diameter. Calculations on larger models with 30 nm diameters may be necessary in the future, however, in order to better characterize system-size effects on the soot model structure.

Polycyclic aromatic hydrocarbons (PAH) strongly partition into soots[8] and may form simultaneously with soots.[9] Hence, the interaction of the PAH pyrene (C₁₆H₁₀) with soot was modeled. Pyrene molecules were placed in various positions near the energy minimized soot model. Initial positioning of the four pyrene molecules was made to reflect π–π (London dispersion), H–π (dipole–dipole), H–O (H-bonding) and H–H (van der Waals) interactions between the pyrene and soot. Thus, the starting configuration was given the opportunity to find stable arrangements for a variety of bonding arrangements rather than pre-selecting only one type of bond. The pyrene–soot system was then subjected to energy minimizations with both the MM+ and COMPASS force fields. Following energy minimization, MD simulations of 10 000 to 100 000 time steps were carried out. The system size was then quadrupled as explained above. Energy minimizations and MD simulations (100 000 time steps) were then performed on the 2732 atom soot model with 16 pyrene molecules using only the COMPASS force field.

The above calculations were performed with the model system isolated in vacuo. To test the possible effects of solvation on the soot structure and partitioning of pyrene, MD simulations were performed with water molecules present. Runs of 100 000 time steps with the COMPASS force field in Cerius[2] (ref. 17) were carried out on the 683 atom soot model including four pyrene molecules surrounded by 2000 water molecules in a nanodroplet (i.e., no periodic boundary conditions). The soot plus pyrene model system was taken from the previous MD simulation of this system in vacuo using the same force field.

Quantum mechanical calculations on fragments of the soot structure were performed with the programs GAUSSIAN 98[20] and DMol[3] (ref. 21). Molecular mechanics force fields describe atom–atom interactions in terms of parameters derived for each atom type. Hence, bonding arrangements not included in the original parameterization may not be modeled adequately. Intermolecular bonding, such as H-bonding or π–π bonding, may be particularly problematic for analytical force fields because the parameters are typically optimized to reproduce properties of isolated molecules. For example, neither of the force fields used in this study contain specific H-bonding parameters. Instead, the force fields rely on interactions between partial charges to describe H-bonding even though H-bonding involves electron transfer.[22] Molecular orbital (MO) or density functional theory (DFT) quantum mechanical calculations are generally more reliable than molecular mechanics results because the quantum methods solve for the electron density of a given molecular configuration. Thus, the quantum mechanical calculations are more flexible and accurate, especially when it comes to modeling unusual bonding or intermolecular forces. On the other hand, the quantum methods are much more computationally intensive than molecular mechanics. Consequently, it is impractical to model systems of hundreds or thousands of atoms with quantum mechanics. Similarly, long MD simulations of molecules are not practical with quantum methods. (Exceptionally demanding calculations are possible, however.[23]) The solution to this problem in this study was to isolate fragments of the soot for quantum mechanical modeling. In this manner, long-range configurations are studied with molecular mechanics, and details of short-range bonding are studied with quantum mechanics.

Two tests were carried out. First, intermolecular H-bonding between phenolic and carboxylic anhydride functional groups thought to exist in hexane soot[1] was computed using MP2/6-31G* calculations in GAUSSIAN 98.[20] The 'MP2' indicates that electron correlation is handled with Møller-Plesset second-order perturbation theory.[24] The "6-31G" indicates that 6 Gaussian functions were used to describe the electron density in each atomic orbital, and the valence orbitals were split into two sets – one of 3 Gaussians and one of 1 Gaussian. The "*" indicates that d-orbitals were included on all atoms except H to account for polarization.[25] This energy minimization was performed to test the accuracy of intermolecular H-bonding in soot as predicted by MD simulations using the COMPASS force field. Second, DMol[3] calculations using double-numeric basis sets with p-polarization functions and gradient-corrected exchange and correlation functionals[26, 27] were performed on a soot fragment and a pyrene molecule. The accuracy of this method should be similar to the MP2/6-31G* calculations above but computationally more efficient. These static energy calculations were used to test the reliability of interactions simulated between soot and pyrene with the COMPASS force field.

Results

Energy minimization

The initial hexane soot model structure[1] consists of aromatic regions from 2 to 64 fused rings that are terminated with H atoms and carbonyl groups. The aromatic fragments are linked via ether linkages or intermolecular H-bonding. Sergides et al.[28] constructed the 3D structure of hexane soot by layering the fragments of the 2D model.[1] Layering is observed in high-resolution transmission electron microscopy (HRTEM) studies.[29] The aromatic fragments of soot are similar to graphite, so this 3D structure is a reasonable approximation. However, the planar configuration of these soot models[1, 28] requires that the ether linkages have C–O–C angles of approximately 180°, whereas normal C–O–C bonds in ethers are nearly tetrahedral (e.g., 112°).

Given this information, it is not surprising that the energy-minimized structures of the soot model do not give an overall planar configuration (Fig. 1a). Both the MM+ and COMPASS energy minimizations predict that the ether linkages bend toward angles significantly smaller than 180° which causes the model soot to contort into a more random configuration. However, the aromatic regions remain relatively planar except for some of the larger fragments that tended to curve significantly. Such curvature is not unexpected because large PAH, such as coronene with 7 fused rings, can form bowl-like structures.

Figure 1

(a) Energy minimized structure of the 683 atom soot model of ref. 1 with the COMPASS force field in Cerius[2]. (ref. 17). C atoms are blue, H are red and O are green, (b) Resulting structure of MD simulation on 683 atom soot model with four pyrene molecules at 300 K for 100 ps with the COMPASS force field in Cerius[2] (ref. 17) starting with energy minimized structure in (a). Graphics generated with CrystalMaker.[34]

MD simulations

The soot structures predicted in MD simulations are dramatically different, depending on which force field was used. MD simulation of the 683 atom soot model for 10 000 time steps at 300 K using the MM+ force field in HyperChem results in little change from the energy minimized structure. (Note: MD simulations were terminated after l0 ps because configurational changes were not observed.) The resulting structure does not reproduce the layered structures observed in HRTEM analyses.[29] Thus, the structure resulting from MM+ is unreliable. MD simulations with the COMPASS force field cause the structure to collapse rapidly. Fig. 1b shows the model soot structure generated with COMPASS-based MD simulations after 100 000 time steps. A layered structure is clearly evident between the sections of the soot model.

Differences in the two force fields were also evident regarding soot–pyrene interactions. COMPASS MD simulations resulted in configurations with the pyrene molecules aligned parallel to the aromatic fragments of the soot model. The distance between the pyrene and soot layers is similar to that of the soot interlayers. The parallel configuration is an indication that π–π interactions between aromatic groups are dominating the long-range energetics. This result is consistent with recent results on solvent-soluble components of coal.[30] Pyrene orientations are nearly perpendicular with respect to the soot layers in the MM+ force field MD simulations (Fig. 2). Bonding between the H atoms of the pyrene and the electrons in the aromatic regions or oxygen atoms on the soot edge are the controlling factors in this case.

Figure 2

MM+ energy-minimized structure of pyrene molecules interacting with the 683 atom soot model from ref. 1. Orientation of pyrene molecules is predominantly end-on indicating that the H atoms of the pyrene are bonding with the oxygen atoms and aromatic groups of the soot model. Graphics generated with CrystalMaker.[34]

COMPASS appears to produce a more realistic soot structure in these MD simulations, so the remainder of the simulations were performed with this force field only. The energy-minimized structure of soot with four pyrene molecules was duplicated to produce the 2732 atom soot model mentioned in the Methods section. An energy minimization was then performed on this larger model. The energy-minimized structure of the 2732 atom model was similar to that of the energy minimized 683 atom system. The molecular fragments were distributed in a random, open configuration with the pyrene molecules loosely associated with the soot layers. A 100 000 time step MD simulation caused the structure to collapse with a regular layering pattern (Fig. 3) as was the case for the 683 atom MD simulation with COMPASS. Again, the pyrene molecules were incorporated into the soot structure without significant change to me interlayer spacing. Significant curvature of some larger aromatic regions is obvious (Fig. 3) suggesting that the nanoparticle is beginning to take a spherical shape. However, larger simulated systems may be necessary to achieve the size and density of real soot particles.

Figure 3

Configuration of the 2732 atom soot model after 100 ps of simulation with the COMPASS force field at 300 K. Sixteen pyrene molecules included in the simulation with the soot oriented themselves in parallel to soot layering with indistinguishable interlayer spacings suggesting that pyrene could fit into voids of the soot without disrupting the structure. Note: the rendering has been changed to a stick configuration to better highlight the layering. Graphics generated with CrystalMaker.[34]

Density analyses[31] of the energy-minimized soot model gives a density of approximately 1.36 g cm^-3. The estimated density of the more compact, layered structure produced using COMPASS (Fig. 1b) is 1.37 g cm ^-3. Soot densities can be as low as 1.57 g cm^-3.[32] The lower density of the model structure may be partially an effect of the small model size because the 2732 atom model has a somewhat higher calculated density of 1.44 g cm ^-3. Larger model systems may be necessary to achieve the density of real soot particles.

Solvation

Fig. 4 presents a close-up view of the soot and pyrene model surrounded by 2000 water molecules. This MD simulation was run at 300 K for 100 000 time steps (100 ps). Although the water molecules surround the soot model and H-bond to the oxygen atoms terminating the soot layers, water does not penetrate between the soot layers. Within the time frame of this simulation, no observable change in soot structure was detected. The possibility existed in this simulation for the H-bonding between soot layers to be disrupted as the polar groups on the edge of the soot model bond to H₂O rather than to neighbors within the soot structure. However, the magnitude of this effect was not significant enough to change the interlayer spacing compared to the simulation in vacuo. No desorption of pyrene from the model soot into the aqueous phase was observed during the simulation. The time scale of this simulation is long enough to allow desorption because 100 ps MD simulation in vacuo allowed a pyrene molecule on the surface to detach from the soot and then re-attach to the soot. Pyrene is a relatively insoluble PAH, however, so simulations with more water soluble PAH (e.g., naphthalene) should be run in the future to examine partitioning behavior.

Figure 4

Configuration of the 683 atom soot model with 2000 water molecules after 100 ps of simulation with the COMPASS force field at 300 K. The soot model structure is indistinguishable from the MD simulation performed without waters of solvation. No movement of the four pyrene molecules from the soot substrate was observed suggesting that partitioning will strongly favor soot in agreement with analytical results.[8] Graphics generated with CrystalMaker.[34]

As was just mentioned above, H-bonding between polar functional groups terminating the aromatic fragments is possible. It has been suggested[1] that such bonding arrangements would be common; and, indeed, the MD simulations in this study predict that H-bonding should occur between soot layers. For instance, the H atom of a phenolic moiety rotates toward an oxygen atom in a carboxylic anhydride edge site in the MD simulation of the 683 atom soot model with the COMPASS force field (Fig. 5a). The H-bond distance in this case is 1.8 A, a moderately strong H-bond.

Figure 5

(a) A close-up view of H-bonding between a phenol group and a carboxylic anhydride group as predicted from the MD simulations of soot with the COMPASS force field. The interlayer spacing is controlled by π–π interactions, but this result indicates that H-bonding between layers may contribute to the bonding energy as well, (b) MP2/6-31G* calculation[20] of H-bonding between phenol and carboxylic anhydride used to test the result of the COMPASS-based MD simulation shows similar H-bonding behavior. High-level molecular orbital theory calculations such as these are not possible on large molecules such as the soot model, but small fragments of the larger molecules may be modeled using this type of quantum mechanical calculations. Graphics generated with CrystalMaker.[34]

Molecular orbital calculations on H-bonding within soot

To test whether this type of H-bond is realistic or just an artifact of the simulation, quantum mechanical calculations on a phenol–carboxylic anhydride dimer were carried out. An energy minimization using the MP2/6-31G* method in GAUSSIAN 98[20] was used to determine the H-bonding distance between these two molecules. The MP2/6-31G* method should give a reasonably accurate prediction in this case because the 6-31G* basis set does a reliable job of reproducing structures for these types of organic molecules, and MP2[24] electron correlation increases the accuracy of the H-bond description. The optimized structure in Fig. 5b shows that the predicted H-bond distance is 2.0 A, somewhat longer than the distance in the soot model. The 0.2 A difference is not dramatic, but somewhat surprising because the edge-on orientation of Fig. 5b should promote stronger H-bonding compared to the parallel orientation of the aromatic rings in Fig. 5a. Thus, one would expect the H-bond distance from the MP2/6-31G* calculations to be smaller rather than larger than that found in the MD simulations. One reason for this increased H-bonding in the MD simulation could be that π–π bonding determines the interlayer spacing such that the H-bonding is optimized within a constrained structure. This interpretation is consistent with the observation that the interlayer spacing does not change significantly in the soot simulation with water present, as discussed above.

Density functional theory calculations on soot fragment and pyrene

Calculations on a fragment of the soot model[1] and one pyrene molecule (Fig. 6) were carried out with the program DMol[3] (ref. 21). The structure of the soot fragment and pyrene were optimized with respect to energy within a constrained minimum 3.5 A distance between C atoms of the soot fragment and pyrene. This spacing was chosen based on the results of the MD simulations of the larger soot system. Two static energy calculations were then carried out. One with the spacing constraint decreased to 3.25 A and the other with the spacing set to 4.0 A. The binding energy of the system decreased from 3.25 to 3.50 to 4.00 A with values of-2338.443, -2338.479, and -2338.490 hartrees, respectively. The monotonic decrease in energy suggests that there is not a minimum potential energy between these two molecules between 3.25 and 4.0 A. Rather, the two molecules appear to prefer separation rather than intermolecular bonding. These results highlight the weakness of density functional calculations in modeling intermolecular dispersion forces.

Figure 6

BLYP/DNP calculation[21] of pyrene and a fragment of the soot model.[1] The plotted spheres define the atomic radii of each element. Note that the two do not overlap, but there is not room for another molecule between the sheets. The plan view illustrates that the aromatic rings are offset to place carbon atoms of the pyrene molecule above the holes in the aromatic rings of the soot fragment. The density functional theory calculations do not show an energy minimum commensurate with the observed spacing in soot or with the MD simulations. Graphics generated with CrystalMaker.[34]

Discussion

The observation that the COMPASS force field performs better for modeling soot than the MM+ force field is not surprising. A force field developed for modeling condensed phase systems should include a better representation of the long-range forces (e.g., van der Waals bonding, π–π bonding, H-bonding) operative in a solid phase like soot compared to a force field developed for gas-phase molecules. The point of this comparison is not to criticize the MM+ force field but to make the point that selection of a force field is of paramount importance in molecular simulations. The availability of user-friendly molecular modeling programs has made MD simulations possible for researchers whose main interest is not molecular simulation. However, the results of simulations will not be accurate unless the force field used to produce the results has been tested on similar systems. Another source of error is that simplifying a system to create a molecular model may eliminate several important physical or chemical factors involved with the process of interest. Thus, in this study, the simulations were geared to increase the level of complexity by increasing the system size and adding waters of solvation. Furthermore, the results were compared against experimental observations and tested by higher-level theory calculations. Perhaps most importantly, the initial model used for the simulations was generated by exhaustive experimental and analytical research.

The soot structure produced and its interactions with pyrene are consistent with experimental observations. This is not to say that more accurate agreement between simulation and experiment is not possible. On the contrary, this study suggests that larger model systems simulated for longer periods of time may produce more accurate soot densities and a more realistic structure. Higher temperature simulations may be necessary to anneal the soot model and overcome long relaxation times in the soot–water simulations. Similarly, longer simulations at elevated temperatures may be useful for studying PAH partitioning and diffusion within soots.

With the above limitations, the current simulation results help explain the observed PAH–soot partitioning and the sluggish desorption of PAH from soot. There are two sources of PAH in soots found in sediments: PAH co-formed with the soot[9] and PAH partitioned into the soot from an aqueous phase.[10] The former group is more strongly sequestered than the latter.[10] Based on the position of pyrene within the model soot, the conclusion of this research is that the pyrene is part of the soot structure and does not expand the interlayer spacing. These structures were generated in MD simulations of soot and pyrene simultaneously, so they illustrate how a PAH could be incorporated with the normal soot structure. Removing the pyrene molecule not only requires that it diffuse through the tortuous micropores of the soot,[33] but the energy lost due to formation of a large "vacancy" in the soot structure must also be regained. If there is no species diffusing into the soot to replace the PAH, then the soot structure must undergo considerable relaxation to anneal out this "vacancy". Under ambient temperature conditions found in soils and surface sediments, large-scale relaxation of the rigid soot structure does not seem likely. Hence, the bioavailability of PAH within soot, especially those PAH that co-formed with the soot, would be limited. Consequently, the environmental risk associated with these sequestered PAH would be decreased.