Molecular Modeling of Porous Carbons Using the Hybrid Reverse

We calculate the ring statistics of the molecular models using the shortest path criteria of ... term in the acceptance probability, in addition to th...
0 downloads 0 Views 432KB Size
Download PDF
9942

Langmuir 2006, 22, 9942-9948

Molecular Modeling of Porous Carbons Using the Hybrid Reverse Monte Carlo Method Surendra K. Jain,† Roland J.-M. Pellenq,‡ Jorge P. Pikunic,§ and Keith E. Gubbins*,† Center for High Performance Simulation and Department of Chemical and Biomolecular Engineering, North Carolina State UniVersity at Raleigh, Box 7905, Raleigh, North Carolina 27695-7905, Centre de Recherche en Matie` re Condense´ e et Nanosciences, Campus de Luminy, 13288 Marseille, Cedex 09, France, and Department of Biochemistry, UniVersity of Oxford, South Parks Road, Oxford OX1 3QU, United Kingdom ReceiVed December 15, 2005. In Final Form: September 6, 2006 We apply a simulation protocol based on the reverse Monte Carlo (RMC) method, which incorporates an energy constraint, to model porous carbons. This method is called hybrid reverse Monte Carlo (HRMC), since it combines the features of the Monte Carlo and reverse Monte Carlo methods. The use of the energy constraint term helps alleviate the problem of the presence of unrealistic features (such as three- and four-membered carbon rings), reported in previous RMC studies of carbons, and also correctly describes the local environment of carbon atoms. The HRMC protocol is used to develop molecular models of saccharose-based porous carbons in which hydrogen atoms are taken into account explicitly in addition to the carbon atoms. We find that the model reproduces the experimental pair correlation function with good accuracy. The local structure differs from that obtained with a previous model (Pikunic, J.; Clinard, C.; Cohaut, N.; Gubbins, K. E.; Guet, J. M.; Pellenq, R. J.-M.; Rannou, I.; Rouzaud, J. N. Langmuir 2003, 19 (20), 8565). We study the local structure by calculating the nearest neighbor distribution, bond angle distribution, and ring statistics.

Introduction Due to their excellent surface activity,1 porous carbons find widespread application in industry for purification and separation of gas and liquid streams2,3 and as catalyst supports.4,5 However, the microscopic structure of these materials is still a matter of debate. Due to their disordered nature, the detailed microstructure cannot be deduced from experimental techniques such as X-ray diffraction or high-resolution transmission electron microscopy. Thus, suitable atomistic models need to be developed that realistically describe the atomic structure of these materials to understand their different properties and to help in tailoring them for specific purposes. The most commonly used model for porous carbons is the slit pore model,6 in which the carbon is represented as stacks of infinite graphene layers, separated by slit-shaped pores. The heterogeneity of the material is thus represented by a range of pore widths. The pore size distribution for a particular carbon is then determined by matching the calculated adsorption for this model (obtained from molecular simulations or density functional theory) to the experimental adsorption isotherm. However, the slit pore model fails to account for pore connectivity or tortuosity or for the curved and defective graphene segments that are shown by X-ray and high-resolution transmission electron microscopy * To whom correspondence should be addressed. Phone: (919) 5132262. Fax: (919) 513-2470. E-mail: [email protected]. † North Carolina State University at Raleigh. ‡ Centre de Recherche en Matie ` re Condense´e et Nanosciences. § University of Oxford. (1) Gregg, S. J.; Sing, K. S. W. Adsorption, Surface Area and Porosity, 2nd ed.; Academic Press: London, 1982. (2) Sircar, S.; Golden, T. C.; Rao, M. B. Carbon 1996, 34 (1), 1. (3) Marsh, H.; Heintz, E. A.; Rodriguez-Reinoso, F. Introduction to Carbon Technologies; University of Alicante: Alicante, Spain, 1997. (4) Rodriguez-Reinoso, F. Carbon 1998, 36, 159. (5) Arumugam, B. K.; Wankat, P. C. Adsorption 1998, 4, 345. (6) Pikunic, J.; Lastokie, C. M.; Gubbins, K. E. In Handbook of porous solids; Schuth, F., Sing, K., Weitkamp, J., Eds.; Wiley-VCH: Weinheim, Germany, 2003; p 182.

studies and so fails to account correctly for the adsorption7 and is particularly inadequate for predicting heats of adsorption or diffusion in these materials. Although there have been several models proposed that incorporate departures from the slit pore model,8 most of these do not provide a realistic description of the highly disordered nature of these carbons. Reconstruction methods, in which a 3D structural model is built that is consistent with a set of experimental structure data, offer the most promising route to realistic models of such carbons at the present time. Reverse Monte Carlo (RMC)9 is one such reconstruction method, in which an atomistic model is built that matches experimental structure factor data from X-ray or neutron diffraction. However, for covalently bonded materials, such as carbons, there is a nonuniqueness problem associated with the direct application of RMC.10 According to the uniqueness theorem of statistical mechanics,11,12 if the intermolecular potential for the material is pairwise additive, the pair correlation function, g(r) (and hence the structure factor), and all higher correlation functions are uniquely determined for a given pair potential; thus, in this case RMC should give the exact (within the experimental accuracy of the data) unique structure. However, most materials do not exhibit pairwise additivity. Chemical bonding is intrinsically a multibody process, as in carbons. In such cases, the structures obtained from unconstrained RMC are nonunique; many molecular structures can match the experimental structure factor. To overcome this nonuniqueness problem, constraints are needed in the fitting procedure. These constraints may be physical or chemical and (7) Coasne, B.; Pikunic, J. P.; Pellenq, R. J.-M.; Gubbins, K. E. Proc. Mater. Res. Soc. 2004, 790, 8.5. (8) For a review see: Bandosz, T. J.; Biggs, M. J.; Gubbins, K. E.; Hattori, Y.; Liyama, T.; Kaneko, K.; Pikunic, J.; Thomson, K. In Chemistry and Physics of Carbons; Radovic, L., Ed.; Marcel Dekker: New York, 2003; Vol. 28, p 41. (9) McGreevy, R. L.; Pustzai, L. Mol. Simul. 1988, 1, 359. (10) Evans, R. Mol. Simul. 1990, 4, 409. (11) Henderson, R. L. Phys. Lett. A 1974, 49, 197. (12) Gray, C. G.; Gubbins, K. E. Theory of Molecular Fluids; Clarendon Press: Oxford, 1984; p 178.

10.1021/la053402z CCC: $33.50 © 2006 American Chemical Society Published on Web 10/28/2006

Molecular Modeling of Porous Carbons

are based upon the understanding of the material being modeled. RMC has been used before in modeling amorphous nonporous carbons incorporating different constraints, such as the ratio of sp2 to sp3 sites, coordination number constraints, bond angle constraints, etc.13-16 Thomson and Gubbins17 modeled a porous carbon as a collection of graphene microcrystals and used the following constraints: (1) an atom can only have two or three neighbors, (2) all interatomic distances are 1.42 Å, (3) all bond angles are 120°. In their final model the pores were slit shaped but connected and the model allowed for the walls to depart from being parallel. Since the model was made of defect-free graphene segments, it did not allow for the formation of five- and sevenmembered rings, which are responsible for the plate curvature observed in TEM images. Pikunic et al.18,19 developed an atomistic model for porous carbons starting from carbon atoms. In addition to fitting to the radial distribution function, they included constraints on the bond angle and the fraction of carbon atoms having three nearest carbon neighbors. All carbon atoms were assumed to be sp2 hybridized, and the hydrogen atoms were implicitly taken into account by this assumption and the neighbor constraint. This simulation protocol, termed constrained RMC (CRMC), was used to model two saccharose-based carbons heat treated at 400 and 1000 °C. The resulting models revealed the disordered nature of the carbons and also contained five- and seven-membered rings in addition to the six-membered rings. Pikunic et al.20 also calculated the isosteric heat of adsorption of argon and nitrogen in their resultant model and obtained excellent agreement with the experiments, suggesting that their model captured the surface heterogeneity present in the sample well. The CRMC procedure does not include any interatomic interaction, so the stability of the resultant models cannot be guaranteed. This procedure was found to yield some unphysical features such as three- and four-membered rings, which are highenergy structures. Recently,21 we studied the stability of the models obtained from CRMC for saccharose-based carbons by relaxing them using two different approaches that realistically describe the interaction between the carbon atoms. The first approach was based on a bond order potential, while the second considered a tight binding model. We found that the local structures obtained from CRMC19 change upon relaxation and the three- and four-membered rings were eliminated upon relaxation. In this work we present a method in which an energy term is included in the acceptance probability for atomic moves, so the simulation seeks to simultaneously minimize the total energy of the system and also the error in the radial distribution function. The presence of the energy term greatly decreases the presence of unrealistic, high-energy structures in the resulting models while simultaneously matching the experimental data. This approach is called hybrid reverse Monte Carlo (HRMC), since it combines the features of the regular Monte Carlo and reverse (13) Bakaev, V. A. J. Chem. Phys. 1995, 102, 1398. (14) Walters, J. K.; Rigden, J. S.; Newport, R. J. Phys. Scr. 1995, T57, 137. (15) O’Malley B.; Snook, I.; McCulloh, D. Phys. ReV. B 1998, 57, 14148. (16) Rosato, V.; Lascovich, J. C.; Santoni, A.; Colombo, L. Int. J. Mod. Phys. C 1998, 9, 917. (17) Thomson, K. T.; Gubbins, K. E. Langmuir 2000, 16, 5761. (18) Pikunic, J.; Clinard, C.; Cohaut, N.; Gubbins, K. E.; Guet, J. M.; Pellenq, R. J.-M.; Rannou, I.; Rouzaud, J. N. Stud. Surf. Sci. Catal. 2002, 144, 19. (19) Pikunic, J.; Clinard, C.; Cohaut, N.; Gubbins, K. E.; Guet, J. M.; Pellenq, R. J.-M.; Rannou, I.; Rouzaud, J. N. Langmuir 2003, 19 (20), 8565. (20) Pikunic, J.; Llewellyn, P.; Pellenq, R. J.-M.; Gubbins, K. E. Langmuir 2005, 21, 4431. (21) Jain, S. K.; Fuhr, J.; Pellenq, R. J.-M.; Pikunic, J.; Bichara, C.; Gubbins, K. E. Stud. Surf. Sci. Catal., in press.

Langmuir, Vol. 22, No. 24, 2006 9943 Table 1. Composition and Density Data for the Carbon Samples sample

H/C

O/C

Hg F (g/mL)

CS400 CS1000 CS1000a

0.53 0.15 0.091

0.123 0.041 0.0087

1.275 1.584 0.722

Table 2. Carbon Density for the Carbon Samples sample

FC (atoms/Å3)

FC (g/cm3)

CS400 CS1000 CS1000a

0.05290 0.07424 0.03621

1.2 1.5 0.9

Monte Carlo methods. Snook and co-workers22-24 used a similar procedure to successfully model amorphous carbons. We use the HRMC method to develop atomic models for three porous saccharose-derived cokes. In contrast to models derived from the CRMC method,18,19 unrealistic high-energy structural features such as three- and four-membered carbon rings are absent in these models. Carbon Samples We develop atomic models for three porous saccharose cokes of various densities and pore structures, namely, CS400, CS1000, and CS1000a; these were previously studied by Pikunic et al.19 (CS400 and CS1000) and by Jain et al.25 (CS1000a) using the CRMC method, thus making a direct comparison of the CRMC and HRMC models here. CS400 and CS1000 were obtained by pyrolizing pure saccharose at 400 and 1000 °C, respectively, under a nitrogen flow. CS1000a, an activated form of CS1000, was obtained by heating CS1000 in an atmosphere of CO2 for 20 h. The experimental (target) radial distribution function, g(r), and the density and composition data for the samples CS400, CS1000, and CS1000a are taken from refs 19 and 25; the density and composition data are shown in Tables 1 and 2.

Simulation Details Reverse Monte Carlo Method. In the reverse Monte Carlo method9 Monte Carlo type atomic moves are made until eventually an atomic configuration of the model system is generated that matches the structural properties of the real system obtained by experiment. Throughout the simulation the differences between the simulation and experimental structural properties are minimized. The most commonly used structural property in the RMC method is the structure factor, S(q), and the quantity to be minimized is nexptl

χ2 )

∑[S

sim(qi)

- Sexptl(qi)]2

(1)

i)1

where Ssim is the structure factor for the model material and Sexptl is the experimental structure factor. After determination of Ssim for a given atomic configuration, atoms are moved randomly in a Monte Carlo procedure to obtain a new configuration. The probability of acceptance of a new atomic configuration is given by

[ (

Pacc ) min 1, exp -

1 (χ 2 - χold2) Tχ new

)]

(2)

where Tχ is a weighting parameter. Radial distribution functions can also be used instead of the structure factor in the RMC studies. The (22) Opletal, G.; Petersen, T.; O’Malley, B.; Snook, I.; McCulloch, D. G.; Marks, N. A.; Yarovsky, I. Mol. Simul. 2002, 28 (10-11), 927. (23) Petersen, T.; Yarovsky, I.; Snook, I.; McCulloch, D. G.; Opletal, G. Carbon 2004, 42, 2457. (24) Opletal, G.; Petersen, T. C.; McCulloch, D. G.; Snook, I. K.; Yarovsky, I. J. Phys.: Condens. Matter 2005, 17, 2605. (25) Jain, S. K.; Pikunic, J.; Pellenq, R. J.-M.; Gubbins, K. E. Adsorption 2005, 11, 355.

9944 Langmuir, Vol. 22, No. 24, 2006

Jain et al.

use of radial distribution functions speeds up the simulation, since the overhead of calculating the Fourier transform to convert gsim(r) to Ssim(q) is absent. Radial distribution functions have been used before in RMC studies.19,25 Hybrid Reverse Monte Carlo Method. In the HRMC method, we introduce an energy penalty term in the acceptance criteria. The energy of the system is calculated using the reactive empirical bond order (REBO) potential of Brenner26 developed for hydrocarbons, which is based on Tersoff’s covalent bonding formalism27 Uij ) VijR(rij) + bijVijA(rij)

(3)

It has a pair repulsive, VijR, and a pair attractive, VijA, potential term and a bond order term, bij, which weights the attractive part of the potential with respect to the repulsive part. The bond order term is a many-body term which depends on the local environment of atoms i and j. A variety of chemical effects that affect the strength of the covalent bonding interaction are all accounted for in this term. Coordination numbers, bond angles, and conjugation effects all contribute to the strength of a particular bonding interaction in the REBO potential. REBO is highly parametrized, and the parameters have been optimized for different properties of diamond, graphite, and small hydrocarbons. The REBO potential can model intramolecular chemical bonding of a wide range of carbons and hydrocarbons and has been used to study diamond, graphite, fullerenes, carbon nanotubes, amorphous carbon, hydrocarbons, etc.28-30 The REBO potential is short ranged and does not contain any dispersion interactions. In our method, we simultaneously minimize the error in the radial distribution function, χ2, as well as the energy of the system. We use the simulated annealing minimization procedure31 to make sure our system is not trapped in a local minimum, as introduced in ref 18. Starting from a random initial configuration, we randomly select an atom and displace it randomly to a new position. We make hydrogen atom moves along with carbon atom moves. The acceptance probabilities for the two kinds of moves are

[ ( (

Pacc ) min 1, exp -

[ ( (

Pacc ) min 1, exp -

))]

1 1 (χ 2 - χold2) + (Unew - Uold) Tχ new w (carbon atom move) (4)

))] (hydrogen atom move)

1 1 (U - Uold) Tχ w new

(5)

where Unew and Uold are the energies of the new and old configurations, respectively, and w is a weighting parameter used to weight the energy term with respect to the structure one. Tχ is a fictitious temperature. We note that the first term on the right-hand side of this equation is the usual reverse Monte Carlo term, while the second term on the right is the usual (canonical) Monte Carlo term. The procedure based on eq 4 is referred to as HRMC. In our simulation method, we start at a high value of Tχ. The system runs at this temperature for a number of moves and is allowed to come to a steady state. The temperature is then lowered by multiplying Tχ with a constant between 0 and 1. The simulation stops when no significant change in χ2 or energy is observed. Setting w to a large number reduces this HRMC to conventional RMC, with just the error in the radial distribution function being minimized. Similarly, setting w to a very small value reduces HRMC to the conventional metropolis Monte Carlo method, with the temperature being reduced slowly. The structural model thus obtained is influenced by the value of the weighting parameter w. We choose the weighting factor to give a (26) Brenner, D. W. Phys. ReV. B 1990, 42 (15), 9458. (27) Tersoff, J. Phys. ReV. Lett. 1988, 61, 2879. (28) Xia, Y.; Ma, Y.; Xing, Y.; Tan, C.; Mei, L. Phys. ReV. B 1999, 61 (16), 11088. (29) Yamaguchi, Y.; Maruyama, S. Chem. Phys. Lett. 1998, 286, 336. (30) Sinnott, S. B.; Colton, R. J.; White, C. T.; Shendeova, O. A.; Brenner, D. W.; Harrison, J. A. J. Vac. Sci. Technol., A 1997, 15, 936. (31) Kirkpatrick, S.; Gelatt, C. D.; Vecchi M. P. Science 1983, 220, 671.

Figure 1. (a) χ2 vs the weighting parameter, (b) energy/atom vs the weighting parameter for CS1000, and (c) number of three-membered rings vs the weighting factor for CS1000a. good fit to the experimental data while yielding structures that are physically realistic, e.g., that have no or very few high-energy configurations, such as three- or four-membered rings. We performed a series of runs with different weighting factors. We show in Figure 1a plot of the final energy, the final error in the radial distribution function χ2, and the total number of three-membered rings in the model as a function of the weighting parameter w. We see that the energy increases and χ2 decreases with increasing weighting parameter. This is expected as a large value of the weighting parameter gives more weight to χ2 and vice versa. Moreover, the number of three-membered rings also increases with increasing weighting factor. Increasing the weighting factor gives less importance to the energy term, and thus, the three-membered rings (high energy structures) are not penalized. It has been pointed out in previous RMC studies of amorphous carbon16,22 that the presence of these small rings (threeand four-membered rings) is an artifact of the RMC method when used without constraints. Thus, we have to choose a weighting parameter that gives a good fit to the radial distribution function and at the same time contains few, if any, three- and four-membered rings. We used weighting factor w ) 25 for CS400, w ) 25 for CS1000, and w ) 15 for CS1000a. Note that the acceptance probability for the carbon atom moves contains both the deviation in the radial distribution function and the energy penalty term, whereas the acceptance probability of hydrogen atom moves contains only the energy penalty term. The experimental radial distribution function is for C-C, so we are only concerned with matching the carboncarbon radial distribution function. In our simulation method, we move a hydrogen atom after every five carbon atom moves. We developed molecular models in a cubic box of 25 Å. Periodic boundary conditions and the minimum image convention were used.

Results Structural Models. Using our HRMC approach, we built molecular models by taking into account hydrogen explicitly. We neglect other heteroatoms, such as oxygen, that might be present in the sample. Carbons prepared from saccharose contain very little of such heteroatoms. It was this feature that led us to choose such carbons for study. We start with a random

Molecular Modeling of Porous Carbons

Langmuir, Vol. 22, No. 24, 2006 9945

Table 3. Energy (eV/atom) Obtained Using HRMC (This Work) and CRMC19,25 sample

energy (eV/atom) for the HRMC model

energy (eV/atom) for the CRMC model

CS400 CS1000 CS1000a

-5.443 -6.386 -6.323

-5.214 -4.752 -4.068

configuration of an appropriate number of carbon and hydrogen atoms. In Table 3, we show the energy (eV/atom) of the final resulting structures. The energy of the structural models obtained from CRMC is also shown for comparison. As expected, the energy of the HRMC models is less than that of the CRMC models. This difference is particularly dramatic for the CS1000 and CS1000a carbons, reflecting the prevalence of high-energy local structures in the CRMC models. We also built molecular models by considering only carbon atoms and ignoring the hydrogen atoms (not shown here). As expected we found that the energy of the different models obtained by taking into account both carbon and hydrogen atoms is more positive as compared to that obtained if hydrogen atoms are not taken into account. This can be explained by the fact that the presence of hydrogen reduces the reactivity of carbon atoms; i.e., the number of carboncarbon double and triple bonds that are expected to form in the absence of hydrogen is reduced. Instead some less energetic C-H bonds are formed. Also, in the absence of hydrogen, carbon atoms are brought closer together (more C-C bonds are made or the C-C bonds become shorter), thus decreasing the energy of the system. Figure 2 presents snapshots of the three carbon materials. The resulting models are very disordered. Many of the carbon atoms are arranged to form rings, and the rings are connected together in a complex way. The models reveal structures that are far from crystalline graphite. The radial distribution functions as computed from the simulation are compared with the experimental results in Figure 3. We note that the experimental and simulated pair correlation functions for all three samples are in good agreement. The HRMC procedure leads to an overestimate of the first peak for CS400 and underestimates the first peak for CS1000a. The agreement is very good for CS1000, where the first and second peaks are reproduced well. The pair correlation functions obtained from these models match the experimental data more closely than those for the models obtained using only carbon atoms. This improvement is less marked for CS1000a, since the amount of hydrogen present is quite small. Furthermore, it can be seen that CS1000a has more structure than the other two samples; the peaks are more pronounced, and longer range correlations are more pronounced. All the carbon samples have a first peak at 1.4 Å with a somewhat broad distribution, suggesting there are some sp3 and sp sites along with the sp2 sites. The second peak is at around 2.4 Å for all the carbon samples. The second peak for graphite is at 2.46 Å. There is a shoulder in CS400 after the second peak at around 2.8 Å, which becomes more pronounced in CS1000 and develops into a small peak for CS1000a. The third neighbor peak in graphite is at 2.86 Å. However, neither the experimental nor simulated g(r) shows a peak at 3.4 Å, which is the layer separation in graphite. This suggests that the samples and our models consist of highly defective graphene segments which are more randomly distributed. Structural Characterization. We characterize the carbon samples by probing the average local environment in the form of nearest neighbor distributions and bond angle distributions. We also calculate the ring statistics, which gives information about the medium-range order present in the system. Neighbor Distribution. We assume that two carbon atoms can

Figure 2. Snapshots of models of (a) CS400, (b) CS1000, and (c) CS1000a obtained from the HRMC method. The gray rods represent C-C bonds, and the red rods represent C-H or H-H bonds. The simulation box size used is 25 Å for all three samples.

have a minimum distance of approach of 1.2 Å, and two carbon atoms are considered bonded if the distance between them is less than 1.6 Å. This prescription for bonding19 was obtained from the end points of the first peak in the radial distribution function (C-C).

9946 Langmuir, Vol. 22, No. 24, 2006

Jain et al.

Figure 4. Bond angle distributions (2-fold-coordinated, carbon atoms with two carbon neighbors) for CS400, CS1000, and CS1000a.

Figure 3. Pair correlation function (C-C) obtained from experiment and from the model: (a) CS400, (b) CS1000, (c) CS1000a. Table 4. Neighbor Distribution for CS400, CS1000, and CS1000a sample CS400 CS1000 CS1000a

zero one two three four neighbors neighbor neighbors neighbors neighbors 0.00242 0 0.00177

0.05925 0.01983 0.04240

0.51149 0.33707 0.39399

0.42443 0.62845 0.55830

0.00242 0.01465 0.00353

The neighbor distribution is shown in Table 4. We calculate the neighbor distribution only for the carbon atoms and consider only those neighbors which are carbon atoms. It is seen that 2-fold- and 3-fold-coordinated atoms dominate the structures of all three samples. As expected, all the hydrogen atoms in the models are singly coordinated; i.e., each hydrogen atom has a hydrogen neighbor or a carbon neighbor. Also, the carbonhydrogen bond length distribution is centered around 1.09 Å with a narrow distribution, which is the optimum bond length for the C-H bond.

Bond Angle Distributions. We report the 2-fold-coordinated bond angle distribution for the carbon atoms, i.e., the C-C-C bond angle distribution, in Figure 4. The bond angle distribution of the 2-fold-coordinated atoms is centered around 120°, with a tail toward longer angles. This is in contrast to the results of simulations obtained by considering only carbon atoms, where the 2-fold bond angle distribution was concentrated around larger angles. The results in Figure 4 show that most of the carbon atoms with two carbon neighbors do have a hydrogen as a third neighbor. However, there are still some bond angles having larger values; this is due to the fact that there are some carbon atoms that have only two carbon neighbors and no hydrogen neighbors. The bond angles for 3-fold-coordinated atoms (not shown here) are centered around 120°, as expected for a carbon atom with three neighbors. The distribution is rather broad about 120°, due to the different local environments of different carbon atoms. In graphite, the distribution is almost a δ function at 120°, since all atoms have the same local environment. We note that in the CRMC method all carbon atoms are assumed to be sp2 hybridized, so the bond angle distribution (for both 2-fold- and 3-foldcoordinated atoms) is centered around 120°. Ring Statistics. We calculate the ring statistics of the molecular models using the shortest path criteria of Franzblau.32 The ring statistics for the HRMC models are shown in Table 5. We also include the ring statistics for the CRMC models for comparison. The rings contain only carbon atoms. Note that hydrogen cannot be part of a ring since hydrogen is singly coordinated, and to be a part of a ring, an atom needs to have at least two neighbors. The six-membered rings dominate for all of the structures. Upon comparing these results with those from simulations ignoring hydrogen atoms, we found that the number of rings (five, six, and seven) decreases upon including hydrogen for all the carbon samples. This is to be expected, since the presence of hydrogen reduces the number of carbon-carbon bonds, and also the number of carbon atoms having three neighbors decreases. The rings reflect the amount of connectivity present in the models. The decrease in the number of rings also results in a decrease in the connectivity. We also calculated the bond angle distribution of the atoms present in a ring. In Figure 5 the bond angle distribution of the five-, six- and seven-membered rings is shown. The fivemembered rings have a distribution centered around 100°. The six-membered rings have a distribution around 120°, as expected. The seven-membered rings have a comparatively broad distribu(32) Franzblau, D. S. Phys. ReV. B 1991, 44 (10), 4925.

Molecular Modeling of Porous Carbons

Langmuir, Vol. 22, No. 24, 2006 9947

Table 5. Ring Statistics (Rings Contain Only Carbon Atoms) for CS400, CS1000, and CS1000a no. of three-membered rings no. of four-membered rings no. of five-membered rings no. of six-membered rings no. of seven-membered rings sample

HRMC

CRMC

HRMC

CRMC

HRMC

CRMC

HRMC

CRMC

HRMC

CRMC

CS400 CS1000 CS1000a

0 1 1

0 22 48

0 0 1

3 58 33

8 9 7

2 56 28

19 65 58

18 55 46

14 44 23

3 45 26

Figure 6. Pore size distribution of the three carbon models. The curves for CS1000 and CS1000a have been shifted upward by 0.1 and 0.2, respectively, for clarity.

We also calculate the geometric pore size distribution (PSD) of the three models using the method of Gelb and Gubbins.34 As we can see from Figure 6 both CS400 and CS1000 have a narrow PSD and are comprised of narrow micropores (widths up to 7.5 Å). CS1000a, however, has a wide PSD with pore sizes reaching 12 Å or more. Upon comparing the PSD of these models with those of the models obtained using only carbon atoms, we found that, for every carbon sample, the PSD for the models obtained here shift toward smaller pores as compared to the PSD obtained for the models simulated with no hydrogen. Adding hydrogen reduces the porosity, causing the PSD to shift to smaller pores. This is especially evident for CS1000. For CS1000a the change in PSD is small, reflecting the small amount of hydrogen present in this carbon.

Conclusion

Figure 5. Bond angle distribution for the bond angles between neighboring atoms in a ring: (a) five-membered rings, (b) sixmembered rings, (c) seven-membered rings.

tion centered around 115°. It has been reported that sevenmembered rings can have a bond angle distribution around 120° and be stable, whereas the five-membered rings do not have a bond angle distribution around 120°.33 This is also observed in the models of our disordered carbons. This again suggests that the local environment of the carbon atoms is modeled correctly by our simulation protocol. Pore Size Distribution. We calculated the porosity of the three models using argon as a test particle. The porosities calculated for the models of porous carbons are 0.232 for CS400, 0.123 for CS1000, and 0.583 for CS1000a. The interaction parameters used in this work are the same as those used in a previous work.25 (33) Lenosky, T. J.; Gonze, X.; Teter, M.; Elser, V. Nature (London) 1992, 355, 333.

In this work we propose and apply a hybrid reverse Monte Carlo method for modeling of porous carbons. The method introduces an energy minimization term in the acceptance probability, in addition to the usual reverse Monte Carlo term, thus combining features of RMC with the conventional Monte Carlo method. The method used here differs from the CRMC method19 in that we replace the bond angle and neighbor constraints by a realistic interatomic potential and account for the presence of hydrogen atoms. The effect of this is to eliminate unrealistic high-energy local structures, such as three- and fourmembered carbon rings, which are a common artifact of structures obtained using the CRMC and other RMC methods. The new HRMC method has been used to build atomic models of three microporous carbons obtained by pyrolizing saccharose. The energies of the HRMC structures are considerably less than those for the CRMC models. Moreover, the many three- and fourmembered carbon rings present in the CRMC models are absent in the HRMC models, in which six-membered rings dominate, with a significant number of five- and seven-membered rings. Thus, these carbons are largely made up of highly defective and (34) Gelb, L. D.; Gubbins, K. E. Langmuir 1999, 15, 305.

9948 Langmuir, Vol. 22, No. 24, 2006

curved graphene segments. These results, together with bond angle and neighbor distributions and the good agreement with experimental radial distribution functions, suggest that the HRMC method is able to give a realistic description of the local atomic structure in microporous carbons. To further characterize the different carbon models, we are currently studying the adsorption of argon in these materials. Calculating the adsorption isotherms and isosteric heats of adsorption in our improved models will provide further insight into the adsorption behavior in these materials.

Jain et al.

Acknowledgment. We thank Isabelle Rannou (CRMD, Universite´ d’Orle´ans, France) for providing us the X-ray results, Nathalie Cohaut and Jean Michael Guet (CRMD, Universite´ d’Orleans, France) for providing us with SAXS and density results, and Dr. S. J. Stuart (Department of Chemistry, Clemson University, South Carolina) for providing us with the REBO code. This work was funded by the National Science Foundation through Grant CTS-0211792. This research used supercomputing resources from HPC-NCSU, SDSC (NSF/MRAC Grant CHE050047S), and NERSC (DOE Grant DE-FGO2-98ER14847). LA053402Z