Curved-Surface Atomic Modeling of Nanoporous ... - ACS Publications

Applied Physics, RMIT UniVersity, GPO Box 2476V, Melbourne 3001, Australia, and UnileVer ... scope Unit, The University of Sydney, NSW, 2006, Australi...
0 downloads 0 Views 1MB Size
Download PDF
802

J. Phys. Chem. C 2007, 111, 802-812

Curved-Surface Atomic Modeling of Nanoporous Carbon Timothy C. Petersen,*,† Ian K. Snook,† Irene Yarovsky,† Dougal G. McCulloch,† and Brendan O’Malley‡ Applied Physics, RMIT UniVersity, GPO Box 2476V, Melbourne 3001, Australia, and UnileVer Corporate Research, Sharnbrook, Bedford, MK44 1LQ, U.K. ReceiVed: June 25, 2006; In Final Form: October 11, 2006

We present a methodology for modeling atomically layered nanoporous carbonaceous solids. Our method models carbon atoms interacting on sets of analytically defined surfaces, in conjunction with reverse Monte Carlo techniques (RMC). Our approach allows the construction of large, physically realistic nanoporous models of low-density carbonaceous solids that are both physical and statistically consistent with experimental observations. To demonstrate the flexibility of our approach and to investigate different mesostructures, we have modeled two nanoporous structures of glassy carbon. One model is comprised of multiple surfaces with the same center of symmetry, and the other is based on sets of parallel surfaces. The resulting models contain highly curved graphitic regions that also resemble the glassy carbon microstructure observed in high-resolution electron microscopy.

1. Introduction Diffraction-based experiments cannot directly determine the atomic structure of disordered materials and, therefore, the structures of the majority of disordered carbon solids remain difficult to describe. Disordered carbonaceous solids are often referred to in crystallographic contexts as “amorphous” or lacking long-range order. Indeed, solids such as hydrogenated or tetrahedral amorphous carbon certainly possess high degrees of structural disorder. However, other carbon structures such as glassy carbon, carbon black, and turbostratic graphite1,2 are much more ordered yet decidedly noncrystalline. Modeling ordered and disordered structures in the absence of high degrees of symmetry is a complex problem. Aromatic bonding of carbon results in the formation of sheets, which may vary in extent and even possess significant curvature or even closure, as in the case of nanotubes and fullerenes. The structural interaction between aromatic sheets necessarily leads to layered and/or porous domains within low-density carbon solids, which further establishes a complicated mesostructure over length scales intermediate to those probed by conventional diffraction and small-angle scattering. Hence, the mesoscopic and atomic microstructure are necessarily coupled and any experimentally consistent model of such solids must describe both length scales together. In a detailed neutron diffraction study of nanoporous carbons, Zetterstrom et al.3 recently discussed the complexities involved with nanoporous modeling and, in particular, the required notion of scale-dependent microstructural densities. One approach to combine the disparate length scales is to consider models based upon graphite fragments, using which both mesoscopic constraints and conventional diffraction data can be statistically modeled via reverse Monte Carlo (RMC) techniques4 or similar * Corresponding author. E-mail: [email protected]. Tel: +61 2 9351 7679. Fax: +61 2 9351 7682. Address: Australian Key Centre for Microscopy & Microanalysis, Madsen Building F09, Electron Microscope Unit, The University of Sydney, NSW, 2006, Australia. † RMIT University. ‡ Unilever Corporate Research.

algorithms.5 These techniques can also be employed with initially random porous carbon structures.6 Alternatively, chemical constraints can be imposed upon graphitic fragments to create nanoporous models.7 Another approach, which will be explored here, is to consider layered surfaces as defining the carbonaceous mesostructure, with aromatic bonding on the surfaces dictating the microstructure. For a chosen mesostructure with this approach, the realistic design of porous carbons requires novel theoretical procedures for modeling aromatic structures on surfaces. Previous research has shown that it is possible to model fullerene-like surfaces such as negatively curved graphite8 in a physically realistic manner.9-13 Giant Bucky spheres have been constructed by explicitly considering topological carbon bonding defects on spheres,14 and other curved surface structures have been considered as well.15-17 Supporting the existence of such exotic structures, toroidal18-21 and zeolite-like carbon models22-24 have been shown to be energetically favorable and physically reasonable.25-28 For example, Hyde and O’Keefe29 have explicitly considered the energetics associated with the elastic warping of graphitic sheets. Calculations of the electronic,30-32 magnetic,33 mechanical, and vibrational properties34 of some of these models have also been reported in the literature. Experimental and theoretical studies have suggested that lowdensity disordered carbon solids can in fact be modeled by interconnected fullerene structures.8,35,36 In a recent publication,37 we developed a method for modeling these types of structures using a Monte Carlo approach to optimize the energies of atoms interacting on analytic surfaces. In this work, we extend our methodology to realistically describe large nanoporous carbon solids consisting of many layered and curved aromatic sheets. To demonstrate the feasibility of our novel techniques, we have chosen to model a glassy carbon solid. In perhaps the most well-known model for glassy carbon, Jenkins and Kawamura38 represent the mesoscopic structure as consisting of sets of intertwined curved sheets like those depicted in Figure 1. We have used hybrid RMC (HRMC) techniques previously39,40 to model the circled portion of Figure 1, which represents the

10.1021/jp063973f CCC: $37.00 © 2007 American Chemical Society Published on Web 12/13/2006

Curved-Surface Atomic Modeling

J. Phys. Chem. C, Vol. 111, No. 2, 2007 803 the surface labeled by the Cartesian coordinates x1 ) x, x2 ) y, x3 ) z(x,y):

(

T3 cn A1 cos-1 2π A3 A3 A2 2πx 2πy cos - cos kˆ (2) T1 A3 T2

s ) xiˆ + yjˆ + z(x,y)kˆ ) xiˆ + yjˆ +

( )

Figure 1. Depiction of the glassy carbon model proposed by Jenkins and Kawamura.32 The circled portion was modeled in our previous investigations33,34 of the glassy carbon microstructure.

nonporous microstructure probed by conventional diffraction experiments. By fusing the HRMC technique with that of our surface modeling algorithms, we present here two large nanoporous structures of glassy carbon consisting of layered surfaces that are physically realistic and consistent with experimental observations over a range of length scales. One model consists of parallel surfaces, and the other is composed of layered surfaces with the same center. These models are referred to throughout as the “parallel” and “concentric” models, respectively. The statistical properties of these models are compared to results from our previous studies, which only modeled the (nonporous) atomic microstructure of glassy carbon. 2. Definitions and Theory 2.1. Hybrid Reverse Monte Carlo. We developed the HRMC procedure previously to describe disordered carbon solids34,41,42 by extending the original RMC algorithm of McGreevy et al.43 Our algorithm seeks to statistically minimize interatomic energies to improve the physical nature of the atomic configurations, while also ensuring consistency of the structure with experimental observations. For the simulations reported here, balance between these two requirements is controlled by weighting the RMC χ2 error term with parameters σ and kT; where kT is the metropolis Monte Carlo (MMC) Boltzmann factor and σ is merely considered constant for all values of the modeled experimental structure factor S(k); where k is the magnitude of the scattering vector. The interatomic interaction energies in our HRMC algorithm are calculated using the environment dependent interaction potential (EDIP), which was first developed by Justo et al.44 and later adapted for carbon by Marks.45 Adapting the HRMC algorithm to model experimentally consistent configurations of carbon atoms interacting on surfaces requires modification to the conventional Monte Carlo step procedure. These details are discussed in Section 2.3 2.2. The Surface Definition. The prominent features in Figure 1 are the interconnected and significantly curved sheets which are, in places, layered. Using ideas from our earlier work on carbon nanostructured surfaces,31 it was hypothesised that layered sets of three-dimensionally periodic surfaces could suitably describe these qualitative aspects of the glassy carbon mesostructure. Fourier approximations46,47 to Schwarz’s minimal p-surface48 were therefore chosen for the defining surface functions φ(s) used in this work 3

φn(s) )

Ai cos(2πxi/Ti) ) cn ∑ i)1

(1)

where Ai, Ti, and cn are scalar parameters and s is a vector on

( ))

2.2.1. Concentric Surfaces. Different parameters cn in eq 1 create highly varied surfaces with a common center of symmetry. Juxtaposing these surfaces together naturally creates a layered structure, which forms a template for the mesostructure of the concentric model presented in this work. The concentric surfaces are defined here by setting Ti ) T and Ai ) 1 in eq 1, with coordinates xi chosen to range between 0 and T, thereby generating simple-cubic periodic boundary conditions. Each cn in eqs 1 and 2 is constant for a given surface and will be referred to here as a “level parameter”. The explicit choice of these surface parameters is motivated by experimental observations, as discussed in Section 3 and calculated in Section 4. 2.2.2. Parallel Surfaces. Given a surface equation φ0(s) ) 0, if the gradient normal to φ0(s) exists over all s of interest, then one can define a family of parallel surfaces, φp(s) ) 0, which are displaced from the base φ0(s) by p multiples of a shortest distance d from s; as represented by eq 3.

(

φp(s) ) φ0 s + pd

∇φ0(s) |∇φ0(s)|

)

(3)

Clearly for some surfaces, such as the Mo¨bius strip, eq 3 may fail because of topological requirements such as orientability. However, we will simply note here that many surfaces exist that are amenable to analysis through eq 3 and are therefore suitable within the context of this work. For consistency with the concentric surface model, we chose the function φ0(s) defined by eq 1 with c0 ) 0 as the base surface for the parallel model. The combination of this base surface with other parallel counterparts given by eq 3 defines the parallel model presented here. The particular choices of the parallel surface parameters are discussed in Sections 3 and 4. 2.3. Displacements of Points on the Surface. In accordance with the work of Holyst et al.,49,50 Monte Carlo moves on the surface are achieved through uniform displacements in the tangent plane combined with projections back down to the surface. The particular algorithm employed here ensures that the step size in the tangent plane is randomly selected from a uniform distribution of desired maximum length, with random direction of displacement. Explicit details are given in our recent publication.31 2.3.1. Displacements on Parallel Surfaces. Many parallel surfaces of, for example, planes, spheres, or toroi, are trivially of the same analytic form as the parent surfaces from which they were derived. However, for other surfaces such as sheets of corrugated iron, for example, their parallel counterparts are not of the same form and may not be analytically known. To calculate Monte Carlo displacements of points on such nontrivial parallel surfaces, we have modified the displacement algorithm of Holyst et al.,43,44 as described in the following paragraph. For an analytically defined base surface φ0(s) furnished with a gradient normal vector n over all s of interest, movements of points on surfaces φp(s) parallel to φ0(s) can be effected through projections from φp(s) down to φ0(s) and vice versa. Equation 4 describes this procedure more accurately, where ni and nf are the gradient normal vectors at the respective initial and final

804 J. Phys. Chem. C, Vol. 111, No. 2, 2007

Petersen et al. 2.5.1. Surface Area of Parallel Surfaces. Without knowing the analytic form of a parallel surface, it is nonetheless possible to calculate the surface area by considering the mean H and Gaussian K curvatures of the (known) base surface alone. Following the work of Hyde et al.51 explicitly, the form of the parallel surface area is given in relation to the base area by eq 6

A)

Figure 2. Movements of points on a parallel surface are effected through normal projection ni down to the base surface φ0(s) (I), uniform displacements t in the tangent plane of φ0(s) (II), projection g from the tangent plane back to φ0(s) (III), followed by final normal projection nf (IV) returning to φp(s).

positions on the base surface φ0(s).

sf ) ((si - ni) + t + g) + nf

(4)

The symbols t and g in eq 4 denote the respective tangent and projection vectors for the base surface. To implement eq 4, each point on φp(s) is accompanied in our algorithm by its base image on φ0(s), with known normal vector ni. For the derived parallel surfaces φp(s), the base surface φ0(s) therefore acts as sub-space on which analytic displacements are carried out in four separate stages, I f IV, as depicted in Figure 2. It should be noted that the normal (and hence also tangent) vectors of φ0(s) and φp(s) are, by definition of φp(s), equal in direction. 2.4. Imposing Symmetry on the Surface. Disordered materials, by their definition, lack symmetry. Therefore, it seems inappropriate that symmetry should play any role in the modeling of such materials. However, we are concerned here with very large systems which, when subdivided into a suitable number of similar sized pieces or subsystems, can still have sufficient flexibility to describe local structural disorder. It should be noted that, although advantageous with regards to computational efficiency, the modeling algorithms presented here are not dependent upon the use of symmetry. Symmetry operations and algorithms used to subdivide the concentric surface into identical patches are discussed in our previous work.31 In accordance with that work, here we also employ a total of 48 patches. As shown in the Appendix, symmetry can be also imposed upon parallel surfaces, using symmetry operations on the defining base surface. 2.5. Surface Area Calculation. As a function of the level parameters cn and boundary length T, the surface area A(cn,T) can be numerically approximated by discrete elements of area over the surface and summing them all up N

A(cn,T) =

x(

) (

δz(i∆x,j∆y) δx

N

∑ ∑ ∆x∆y i)1 j)1

2

+

)

δz(i∆x,j∆y) δy

(5)

2

+1

where ∆x ) ∆y ) T/N and the form of z(x,y) was given in eq 2. Using large values of N, eq 5 was used to evaluate the areas for all surfaces in this work.

∫S dA + 2nd∫S HdA + (nd)2∫S KdA

(6)

where S is the base surface to be integrated over, dA is an infinitesimal element of the base surface area, and nd is the normal separation between the base surface and its nth parallel counterpart. Holyst et al.52 recently expressed the curvatures H and K in a convenient form for direct numerical integration, the explicit details of which are listed in the Appendix. 3. Experimental Observations Parts a and b of Figure 3 show high-resolution transmission electron microscopy (TEM) images of a glassy carbon solid, synthesized by heat treatment to 2500 °C, which we refer to as V25 glassy carbon. The tangled graphitic microstructure portrayed in Figure 1 can be inferred from parts a and b of Figure 3, which were collected using a JEOL 2010 TEM and zero-loss energy filtering with a Gatan Imaging Filter (GIF2000). Importantly, the sample used to generate these images is the same as that from which the experimental diffraction data was collected. 4. Surface Design 4.1. The Nanoporous Concentric Model. In previous investigations,33,34 it was noted that the V25 glassy carbon microstructure is more appropriately described by a local density equal to that of graphite (2.27 g/cm3), differing considerably from the reported bulk density (1.55 g/cm3).33 The concentric system modeled here was therefore designed to have a bulk density equal to that of glassy carbon, in conjunction with a local microstructural density equal to that of graphite. Even for slit-like porous systems, these two constraints are difficult to implement3 and are dependent, in any case, upon what one means by “local density”. Images presented in Figure 3a and b suggest that a model based on approximately 10 nonparallel concentric surfaces defined by a certain range of the cn in eq 1 would qualitatively describe some of the mesoscopic features. To produce lowenergy configurations and partly satisfy density requirements, preliminary analysis and our previous work31 indicated that it is appropriate to set the density F of atoms on each surface to that of graphite (F ) Fg = 0.3829 atoms/Å2). The bulk density Fb is fixed by the total number of atoms and boundary conditions imposed by the Ti parameters in eq 1. However, to achieve even qualitative agreement with diffraction-based observations, the microstructural density must be further refined by imposing significant correlations between surfaces relating to the graphitic (002) diffraction peak. These requirements imply that neighboring surfaces must, in some regions, be approximately separated by the graphitic c-axis distance d = 3.35 Å. These c-axis correlations were chosen to be enhanced in regions of φ(s) containing small surface curvature, defined by the faces of the cubic boundaries of φ(s). In these regions, by appropriate derivation of the level parameters cn, adjacent surfaces were chosen to be separated by the graphitic c-axis distance in the

Curved-Surface Atomic Modeling

J. Phys. Chem. C, Vol. 111, No. 2, 2007 805

Figure 3. High-resolution TEM images of an ultrathin V25 glassy carbon film show evidence of highly curved, adjacent graphitic layers. The drawn arrows mark regions of interest, where highly curved and layered structures are evident.

ci that satisfy all of the desired design aspects for the mesostructure template. The level parameters are fixed by the distance d between the intersection points in Figure 4 and were derived using proof by induction (see the Appendix), leading to eq 8:

cn ) 1 + 2 cos

Figure 4. P-surface based concentric surfaces in cross section (x, y, or z ) 0), where the small circle marks the projected center of the model and the arrows are along the main diagonals of the faces of the simple cubic cell in which the surfaces are periodically bounded.

manner shown in Figure 4, which shows a cross section of 10 concentric p-surfaces along the cubic face: The 10 concentric surfaces were designed to orthogonally intersect the arrows in Figure 4 at precise intervals. At these intersections, the distances between each pair of surfaces were set equal to the graphitic (002) pair correlation length (d = 3.35 Å). To calculate the correct level parameters for the design suggested in Figure 4, while also satisfying both surface and bulk-density requirements, additional analysis was required as explained in the next section. 4.1.1. Determining the Concentric Model LeVel Parameters cn and Cubic Cell Length T. Irrespective of the separation between neighboring surfaces, another requirement for ensuring local graphitic density is that the surface density of atoms must also equal that of graphite. Thus, the total number of atoms NTotal in the system is controlled by two factors: the bulk density Fb (of glassy carbon) and the surface density Fg (of graphite). These two conditions can be jointly fulfilled by noting that

Fg

∑n A(cn,T) ) NTotal ) T3Fb.

so that minimizing

χ2({ci},T) )

(∑ Fg

Fb

)

A(cn,T) - T3

n

2

(7)

with respect to T allows the determination of both T and the 10

2π π cos( ) + ) (2πnd T 8 3

(8)

The 10 concentric surfaces in the bulk model were chosen to be defined by c-5 through to c4. Equation 8 was then used to iteratively minimize eq 7 and numerically arrive at the initial cell parameters listed in Table 1, where the reduced area refers to the surface area for T ) 1. As designed, the parameters in Table 1 yield the following characteristics of the bulk model: T ) 102.3 Å; bulk density, 1.55 g/cm3; surface density for each surface, 0.3829 atoms/Å2; total number of atoms, 83142. To implement 48 imposed symmetry elements on each surface during the surface-modified HRMC calculations of the concentric model, we rounded the number of atoms on each surface in Table 1 to the nearest factor of 48, giving a new total number of atoms: 83 136 (or 1732 unique atoms). 4.2. The Nanoporous Parallel Model. For a carbonaceous system of parallel surfaces, a suitable choice for the local density is to set both the atomic surface density equal to that of graphite Fg = 0.3818 atoms/Å2 and the shortest distance between neighboring surfaces equal to the c-axis distance in graphite (d = 3.35 Å). Appropriately, this definition of local density becomes equal to the bulk density of graphite for systems comprised of infinitely many parallel and planar (non-curved) surfaces. For a fixed number of N parallel surfaces bounded within a cubic cell, merging these bulk and local density requirements together leads to a unique value of the cell boundary length T, which is arrived at by the fact that TABLE 1: Initial Surface Parameters for the Concentric Model level number (n)

cn

-5 -4 -3 -2 -1 0 1 2 3 4

.405 1.117 0.826 0.539 0.262 0.000 -0.241 -0.455 -0.638 -0.787

reduced area 1.279 1.635 2.059 2.235 2.323 2.354 2.329 2.274 2.184 2.081

number of atoms 5125 6549 8250 8953 9306 9433 9330 9110 8749 8337

806 J. Phys. Chem. C, Vol. 111, No. 2, 2007

Petersen et al.

N

FBT 3 ) NTotal ) FBT 2

∑A1n(nd) n)1

so that

T)

Fg FB

N

A1n(nd) ∑ n)1

(9)

where A1n(nd) is the reduced surface area for a surface with T ) 1, NTotal is the total number of atoms inside the cubic cell, N is the number of surfaces, and the dependence of A1n(nd) on nd is given by eq 6. Combining eqs 6 and 9 together, T is therefore completely specified by eq 10:

T)

Fg Fb

Fg

N

∫S dA + F ∑{2nd∫S H dA + (nd)2∫S K dA} bn)1

(10)

In eq 10, S refers to the surface region containing all points in eq 10 that lie in the range 0 < xi e T, with T ) 1. Because the surfaces are all parallel, it was decided that N ) 7 surfaces would give sufficient (002) graphitic correlations and that a smaller number of surfaces would not be consistent with the qualitative high-resolution TEM observations of the projected V25 glassy carbon microstructure shown in Figure 3a and b. Using eq 6 with the number of discrete grid points set to N ) 80 000, the integrals in eq 10 were approximated to give

∫S dA = 2.35, ∫S H dA = 0.00, and ∫S K dA = -25.18 which yields T = 74.23 Å. These values produce equivalent details to those in Table 1, through use of eq 6. Rounding the number of atoms on each surface to the nearest factor of 48, as required for the imposition of symmetry, this value of T gives Natoms ) 31 776 with the bulk density equal to the desired value FB ) 1.55 g/cm3. 5. Simulation Details 5.1. General Simulation Details. 5.1.1. Multisurface Energy Calculations. For carbon-based multiple-surface models, constituent neighboring surfaces should be separated by distances approximately equal to the graphitic c-axis distance. Combined with the fact that this c-axis distance is close to the cutoff distance for potentials such as the EDIP, beyond which interaction energies are defined to be zero, it is therefore sensible to calculate the energy of a given atom using only those atoms embedded in the same surface. With respect to energy calculations of a single atom, decoupling the surfaces in this manner effectively reduces the total number of atoms by a factor equal to the number of surfaces in the bulk model. When applying symmetries on each surface in conjunction with this additional restriction to intrasurface interactions, the number of unique atoms requiring energy evaluation can be vastly reduced. For part of the nanoporous simulations reported here, energetic interactions were therefore calculated for atoms belonging to identical surfaces only by using surface lists in our algorithms. 5.1.2. Inhomogenous Density Corrections. It is implicitly assumed in the calculation of the radial distribution function g(r) within the HRMC algorithm that the density of the system is the same in all regions. For nanoporous systems, the microstructure is necessarily coupled to local density inhomogeneities, caused by the mesostructure, on the length scale to

which g(r) is sensitive to. These local density variations modulate g(r) to an extent which creates large, unphysical ripples in the resulting structure factor S(k). To overcome these issues and to more appropriately model the microstructural parts of the bulk system, g(r) was truncated in the surface-modified HRMC algorithm to a suitably small value of r about which g(r) approached unity. The density in the Fourier relation between S(k) and g(r) was also set to the microstructural density of graphite, rather than the bulk glassy carbon value with which g(r) was calculated. The use of the graphitic density for the microstructure is in accordance with previous RMC and HRMC investigations of glassy carbon.33,34 5.2. Concentric Model Details. Atoms were randomly assigned positions with Cartesian x, y, and z coordinates ranging from 0 to T ) 102.3 Å in order to build the initial configuration for the concentric nanoporous model. Equation 2 was then tested for each randomly chosen position so that atoms were uniformly distributed on the surface within a small tolerance level of |φn(s) - cn| < 0.005. Hard core constraints were also imposed upon candidate surface positions so that no two atoms were separated by less than 0.5 Å. Beginning from the resulting structure and using the 48 symmetry elements defined by τ ) 0 in eq 7, the surface-modified HRMC code was run with only EDIP energy constraints switched on. These initial calculations were performed to minimize the energy of atoms on the surfaces down to an average energy of -6.80 eV per atom. The energy minimization was performed in three successive stages, using temperatures T ) 4300 K, 2300 K, and 300 K, while the decoupling of the interactions between surfaces allowed each surface to be simulated separately and consequently executed in parallel, thereby further increasing computational efficiency. After the energy optimization stage, the surface modeling code was executed serially, using both S(k) and EDIP constraints (but not g(r)) with the weightings σ ) 0.02 and T ) 2500 K. In accordance with previous investigations, the HRMC inputs used the experimental diffraction data of McCulloch et al.53 According to Section 5.1.2., the maximum r in g(r) was set to 29.5 Å and the local density was set to 2.27 g/cm3. Experimental consistency was deemed satisfactory after these simulations, producing configurations with reasonable average energies of approximately -6.35 eV per atom. 5.3. Parallel Model Details. To construct the initial coordinates of the parallel nanoporous system, we tested eq 2 for randomly chosen atom positions so that all 31 776 atoms were uniformly distributed on the base surface φ0(s), within a small tolerance level of |φ0(s)| < 0.005. Groups of these atoms were then promoted to the six parallel surfaces on either side of the base surface (i.e., φ-3(s) through to φ+3(s)), and the initial S(k) was calculated by setting the maximum r in g(r) to 23.0 Å according to Section 5.1.2., the results of which are shown in Figure 5. Imposing 48 symmetry operations throughout, the parallel surface algorithm was serially executed using intrasurface EDIP energy minimization only in three separate stages defined by the temperatures T ) 4000 K, 2000 K, and then 300 K. Annealing the surface atoms in this manner reduced the average energy per atom down to -6.63 eV per atom, which was deemed an acceptable value for the purpose of preparing an initial configuration for HRMC simulation. From inspection of the initial S(k) shown in Figure 5, it was clear that the surface constraints had to be removed in order to reduce the large degree of intersurface structural correlations arising from the perfectly parallel surfaces. To also facilitate buckling of the surfaces and allow for the occurrence of atoms

Curved-Surface Atomic Modeling

J. Phys. Chem. C, Vol. 111, No. 2, 2007 807

Figure 5. S(k) calculated from the initial bulk configuration (before energy optimization), which was generated by placing atoms at virtually random positions on each of the seven surfaces.

between surfaces, as reported in previous work,34 all surface constraints and symmetry operations were therefore switched off. Optimizing the EDIP energy and S(k) consistency only, several HRMC runs of the 31 776 now unique atoms were then conducted on the bulk system with the aforementioned maximum in g(r), energy weighting T ) 2500 K and various values for the S(k) weighting parameter σ. It was found that convergence of the theoretical S(k) toward that from experiment was impractical for values of σ greater than σ = 0.01. The weightings T ) 2500 K and σ ) 0.01 were thus used for extended calculations until the average energy converged to around -6.30 eV per atom and the consistency with the experimental S(k) became effectively unchanged by further simulation.

Figure 6. (a) Theoretical structure factor S(k) for the concentric model compared to that from experiment. (b) Theoretical structure factor S(k) for the parallel model compared to that from experiment.

6. Results and Discussion From the S(k) for the (initially random) parallel model shown in Figure 5, with sharp peaks separated by intervals of ∆k ) 2π/3.35 Å-1, it is clear that the use of just seven parallel surfaces creates significant graphitic (002) correlations; irrespective of the large structural disorder of the initial surface atoms. As Figure 5 serves to demonstrate, it would not therefore be possible to optimize these features using surface-based HRMC modeling, without modifying the analytic form, or at least parameters, of each of the constituent surfaces in the bulk model. Use of fewer surfaces would also reduce these correlations but would then produce inconsistencies with the experimental images in Figure 3a and b. Theoretical S(k) and g(r) functions from the surface-modified HRMC calculations of the concentric and parallel nanoporous models are compared against experiment in Figure 6a and b, respectively. For both models, Figure 6a and b show some discrepancies between k ) 5 Å-1 and 7 Å-1; however, the overall consistency between experiment and theory is good. Figure 7a and b shows two different theoretical radial distribution functions. The labels “HRMC” in Figure 7a and b refer to the g(r) values calculated by assuming that the local density is constant in all regions of the bulk model. The “HRMC Inverted” g(r) values were corrected for local density variations by first calculating S(k) from the Fourier transform of g(r), then setting to zero all S(k) fluctuations below the corresponding minimum k value in the experimental S(k), and finally Fourier inverting to retrieve g(r). This correction procedure is effective here because the modulation of the “HRMC” g(r) can be considered as the superposition of a small number of slowly varying trigonometric functions. When Fourier transformed, it follows that these low-frequency backgrounds are effectively

Figure 7. Theoretical g(r) values for the concentric model (a) and parallel model (b) after surface-modified HRMC simulation. The “HRMC inverted g(r)” has been produced by removing the effects of inhomogenous local density variations from the HRMC g(r).

removed by high-pass filtering. After the correction procedure, it is seen that overall agreement with the experimental g(r) has been achieved where, like previous investigations,34 there is a small overshoot in the first nearest-neighbor peak. There also remains a small discrepancy for the second peak of the concentric model g(r). Figure 8 compares the atomic coordination number distributions for the parallel and concentric models with the HRMC

808 J. Phys. Chem. C, Vol. 111, No. 2, 2007

Petersen et al.

Figure 8. Coordination distributions for the nanoporous V25 glassy carbon models compared to that of a nonporous model from previous work.34

Figure 10. Bond angle distribution resolved into weighted contributions from atoms at the vertices of each bond with differing coordination types.

Figure 9. Ring size distributions for the parallel and concentric nanoporous models compared to that of the nonporous graphitic34 model of V25 glassy carbon.

results from our previous investigation,34 which simulated the graphitic microstructure starting from an initial graphite lattice. The coordination numbers were calculated from the HRMC atomic networks using a first nearest-neighbor radius of 1.85 Å. Apart from the small number of 5-fold coordinated atoms, a close match between the two distributions is apparent. Given that the atoms in the concentric model simulation were constrained to move on surfaces only, the production of 4-fold coordinated atoms was not expected and is therefore surprising. Because the overall percentages of 4-fold atoms are similar in all distributions, the data implies the presence of 4-fold atoms is perhaps dictated by the diffraction data. Upon these considerations, it is important to note that g(r) was not fitted in either nanoporous simulation. Figure 9 displays the ring size distribution for both the parallel and concentric nanoporous models compared to the microstructural model from our previous work.34 The ring distributions were calculated using the shortest path criterion of Franzblau.54 Although there is a dominant number of six member rings in the concentric model, it is clear that the pure restriction of atoms to surfaces has resulted in the existence of other rings that were barely represented in the previous nonporous V25 simulations. Although the presence of five, seven, and even eight member rings is to be expected from the curvature of the surfaces in the concentric model, on chemical grounds, the small number of three and four member rings is a nonphysical attribute, likely related to the existence of 4-fold coordinated surface atoms. Supporting these conclusions, the parallel model results in Figure 9 show no significant proportions of small rings.

Figure 10 shows the bond angle subdistributions for atoms at the vertices of each bond, labeled by particular coordination types. These subdistributions are weighted by the respective coordination percentages and normalized such that the sum of all subdistributions gives the total bond angle distribution. Expectedly, the graphitic 3-fold vertices in both nanoporous models possess bond angles peaked at θ = 120°; however, the 4-fold coordinated vertices for the concentric model are only very broadly centered about the diamond angle θ = 110°. It appears that the EDIP energy constraint has attempted to optimize the geometry of these atoms without much success, the failure of which was anticipated from the fact that the movement of concentric model atoms was constrained to surfaces, whereas diamond-like bonding requires geometry available in strictly three-dimensions. The surface-relaxed parallel model distributions support this notion, where the peak at θ = 110° is much more pronounced. Given the appearance of the significant percentage of 4-fold atoms in Figure 9, Figure 10 further implies that inclusion of bonded atoms between graphitic surfaces, as observed in previous investigations,34 is perhaps a necessary feature of glassy carbon. Atomic configurations generated after the initial EDIP energy minimization stage of the concentric model simulations are shown in Figures 11, 12, and 13. Figures 11 and 12 represent small sections of the total bulk structure and display features qualitatively similar to the experimental high-resolution TEM images of Figure 3a and b. The curvaceous and porous aspects of these configurations are a direct result of the surface design, which remained effectively unchanged throughout the simulations. Figure 11 also exhibits both graphitic bonding and symmetry on the surfaces, whereas Figure 12 shows outer graphitic layers that appear to break away from the other concentric surfaces, similar to the marked features in Figure 3b. Re-imaging the periodic boundary conditions revealed that these two outermost surfaces are in fact completely closed pores and are almost spherical in shape. Figure 13 more clearly demonstrates the graphitic bonding on the surfaces induced by

Curved-Surface Atomic Modeling

J. Phys. Chem. C, Vol. 111, No. 2, 2007 809

Figure 11. Interior section of the concentric model after energy minimization of atoms on the surfaces.

Figure 12. Different view of the concentric model after EDIP energy optimization. The two outer layers are in fact completely closed pores, almost spherical in shape.

the EDIP potential and the effect of the imposed symmetries upon the bonding. Different sections of the parallel model structure are shown in Figures 14 and 15, where lighter shaded features indicate 3-fold coordinated atomic bonds and thick dark lines have been used to represent 4-fold atoms bonded to other 4-fold atoms. Figure 14 shows the curvature of the roughly parallel surfaces, which have been buckled in the HRMC optimization process (throughout which the surface constraints were removed). Agreeing with observations from our previous work,34 both Figures 14 and 15 show that diamond-like atoms are located between surfaces, apparently forming cross-links that pin the surfaces together. Created by joining 3 × 3 × 3 periodic blocks of the parallel model together, Figure 16 displays an inner section of the

nanoporous architecture, where all bonds have been assigned light shades. Some graphitic bonds are evident near the edges; however, Figure 16 more prominently shows the manner in which the buckled carbonaceous surfaces are interconnected about the nanoporous domains. 7. Conclusions By confining carbon atoms to curved surfaces, we have shown that it is possible to realistically model low-density carbonaceous solids such as glassy carbon. It has been demonstrated that, by careful design of these surfaces, large nanoporous models that simultaneously describe both micro- and mesostructure can be constructed to be consistent with experimental observations over a wide range of length scales. Practical application of these ideas

810 J. Phys. Chem. C, Vol. 111, No. 2, 2007

Petersen et al. Inverting the above expression gives yn

1 T cos-1 (cn - 1) 2π 2

(

yn )

)

If the separation between aromatic surfaces along the boundary diagonals is defined as d, then the following relationship is sought:

yn+1 - yn ) d cos

(π8)

This gives

T 1 1 π T cos-1 (cn+1 - 1) cos-1 (cn - 1) ) d cos 2π 2 2π 2 8

(

)

(

)

()

or, upon rearranging

1 2π π 1 (c - 1) ) cos d cos + cos-1 (cn - 1) 2 n+1 T 8 2

(

Figure 13. Thin section of the concentric model after energy minimization, more clearly showing both the imposed symmetries and desired graphitic bonding on the curved surfaces.

in conjunction with HRMC was demonstrated here via the construction of large nanoporous models of glassy carbon. As designed, the mesostructure of these models also qualitatively supports high-resolution TEM observations. To extend this work, it would be beneficial for one to consider simulating low angle scattering in addition to conventional diffraction-based data, so as to better cope with density inhomogeneities as well as to achieve further consistency with experiment. By adjusting the parameters in eq 1, or any parametrized surface equation, it is deemed possible to computationally model surface structure within our surface-modified HRMC algorithm for curved graphitic structures, rather than rely on qualitative interpretations of the mesostructure from high-resolution TEM images. The inclusion of some unconstrained atoms in fixed-surface structures (like the concentric model) would also allow for the possibility of diamond-like atoms to bond between surfaces without relaxing the surface constraints of the majority. Other improvements could be made by performing these simulations with different local densities and different energetic annealing schemes.55 Ultimately, rather than using qualitative interpretations from high-resolution TEM image projections or statistical small-angle diffraction data, electron tomographic techniques could be used to extract the mesostructure within certain regions, using which certain surface functions φ could be analytically fitted and subsequently modeled. Acknowledgment. We thank Dr. Nigel Marks (University of Sydney) for use of his EDIP energy code. We also thank the Australian Partnership for Advanced Computing (APAC) and the Victorian Partnership for Advanced Computing (VPAC) for providing computer time for the calculations presented here. Appendix A.1. Derivation of the Level Parameters. Choose the projection x ) 0 (say) so that along the line defined by y ) z

φn(x,y,z) ) 1 + 2 cos

( )

2πyn ) cn T

()

(

))

Now, the separation between neighboring surfaces must also obey

T T 1 1 π cos-1 (cn+2 - 1) cos-1 (cn+1 - 1) ) d cos 2π 2 2π 2 8

(

)

(

)

()

and inserting the previous expression into the above yields

T 1 1 π T cos-1 (cn+2 - 1) cos-1 (cn - 1) ) 2d cos 2π 2 2π 2 8

(

)

(

)

()

The same procedure can be applied using yn+3 - yn+2 to give

T T 1 1 π cos-1 (cn+3 - 1) cos-1 (cn - 1) ) 3d cos 2π 2 2π 2 8

(

)

(

)

()

so that in general

T 1 1 π T cos-1 (cn+m - 1) cos-1 (cn - 1) ) md cos 2π 2 2π 2 8

(

)

(

)

()

By definition, c ) 0 when n ) 0 so that

T 1 π T 1 cos-1 (cm - 1) cos-1 - ) md cos 2π 2 2π 2 8

(

)

( )

()

Solving for cm finally gives the desired equation A1 (eq 8 in the text):

cm ) 1 + 2 cos

2π π cos( ) + ) (2πmd T 8 3

(A1)

A.2. Imposing Symmetry on Parallel Surfaces. For parallel surfaces defined by eq 4, the symmetries of the base surface φ0(s) are preserved because if φ0(Ss) ) φ0(s′) ) φ0(s), where the vectors s and s′ are positions on the surface φ0(s) and S is a symmetric operation that takes s to s′ then we also have

(

φn(Sˆ s) ) φ0 Sˆ s + nd

∇φ0(Sˆ s) |∇φ0(Sˆ s)|

) (

) φ0 s′ + nd

)

∇φ0(s′) ) |∇φ0(s′)| φn(s′) (A2)

Although it is tempting to interpret the far-left-hand side of equation A2 as involving symmetry operations on points lying on the parallel surface φn(s), the next term in the equation demonstrates that S must act on points lying on the base surface φ0(s) only. Using these concepts in our parallel algorithm,

Curved-Surface Atomic Modeling

J. Phys. Chem. C, Vol. 111, No. 2, 2007 811

Figure 14. Interior section of the parallel surface after HRMC modeling without surface constraints; darkest features highlight the 4-fold coordinated bonds.

Figure 15. Interior section of the parallel model after HRMC modeling without surface constraints; darkest features highlight the 4-fold coordinated bonds.

Figure 16. Periodically expanded section of the parallel model more clearly showing the amalgamation of nanoporous domains with the microstructure.

symmetries are therefore imposed upon the parallel surfaces φn(s) by applying operations such as eq 7 to the base surface φ0(s) only, which are implemented after executing step III in Figure 2. A.3. Gaussian and Mean Curvature. Denoting the partial Cartesian derivatives of a surface function φ as φx, φy, and φz, the Gaussian K and mean H curvatures are given by equations A3 and A4

C ) φ2x (φ2yz - φ2y φ2z ) + φ2y (φ2xz - φ2x φ2z ) + φ2z (φ2xy - φ2x φ2y ) + 2φxφz(φxzφyy - φxyφyz) + 2φxφy(φxyφzz - φxzφyz) + 2φyφz (φyzφxx - φxyφxz)

H)-

K) where

1

B 2A

2xφ2x + φ2y + φz φ2x

C 1 2 2 A + φy + φz

B ) φ2x (φ2y + φ2z ) + φ2y (φ2x + φ2z ) + φ2z (φ2x + φ2y ) 2φxφyφxy - 2φxφzφxz - 2φyφzφyz and

(A3) A ) -φ2x - φ2y - φ2z (A4)

References and Notes (1) Millward, G. R.; Jefferson, D. A. Lattic Resolution of Carbons by Electron Microscopy. In Chemistry and Physics of Carbon; Walker, P. L.; Ed.; Marcel Dekker: New York, 1970; pp 125-190.

812 J. Phys. Chem. C, Vol. 111, No. 2, 2007 (2) Harris, P. J. F.; Burian, A.; Duber S. Philos. Mag. Lett. 2000, 80, 381. (3) Zetterstrom, P.; Urbonaite, S.; Lindberg, F.; Delaplane, R. G.; Leis, J.; Svensson, G. J. Phys.: Condens. Matter 2005, 17, 3509. (4) Thomson, K. T.; Gubbins, K. E. Langmuir 2000, 16, 5761. (5) Petkov, V.; DiFrancesco, R. G.; Billinge, S. J. L.; Acharya, M.; Foley, H. C. Philos. Mag. B 1999, 79, 1519. (6) Pikunic, J.; Clinard, C.; Cohaut, N.; Gubbins, K. E.; Guet, J. M.; Pellenq, R. J. M.; Rannou, I.; Rouzaud, J. N. Langmuir 2003, 19, 8565. (7) Acharya, M.; Strano, M. S.; Mathews, J. P.; Billinge, S. J. L.; Petkov, V.; Subramoney, S.; Foley, H. C. Philos. Mag. B 1999, 79, 1499. (8) Mackay, A.L.; Terrones, H. Nature 1991, 352, 762. (9) Lenosky, T.; Gonze, X.; Teter, M.; Elser, V. Nature 1992, 355, 333. (10) Terrones, H.; Mackay, A. L. Prog. Cryst. Growth 1997, 34, 25. (11) Mackay, A. L.; Terrones, H. Philos. Trans. R. Soc. London, Ser. A 1993, 343, 113. (12) O’Keeffe, M.; Adams, G. B.; Sankey, O. F. Phys. ReV. Lett. 1992, 68, 2325. (13) Terrones, H.; Fayos, J.; Aragon, J. L. Acta Metall. Mater. 1994, 42, 2687. (14) Perez-Garrido, A. Phys. ReV. B 2000, 62, 6979. (15) Terrones, M.; Hsu, W. K.; Hare, J. P.; Walton, D. R. M.; Kroto, H. W.; Terrones, H. Philos. Trans. R. Soc. London, Ser. A 1996, 354, 2025. (16) Terrones, H.; Terrones, M. Phys. ReV. B 1997, 55, 9969. (17) Terrones, H.; Terrones, M. Fullerene Sci. Technol. 1998, 6, 751. (18) Johnson, J. K.; Davidson, B. N.; Pederson, M. R.; Broughton, J. Q. Phys. ReV. B 1994, 50, 17575. (19) Yang, L.; Jiang, J.; Dong, J. Phys. Status Solidi B 2003, 238, 115. (20) Oh, D. H.; Park, J. M.; Kim, K. S. Phys. ReV. B 2000, 62, 1600. (21) Itoh, S.; Ordejon, P.; Drabold, D. A.; Martin, R. M. Phys. ReV. B 1996, 53, 2132. (22) Nesper, R.; Vogel, K.; Blochl, P. E. Angew. Chem., Int. Ed. Engl. 1993, 32, 701. (23) Nesper, R.; Leoni, S. Chemphyschem 2001, 2, 413. (24) Hyde, S. T.; Ramsden, S.; Di Matteo, T.; Longdell, J. J. Solid State Sci. 2003, 5, 35. (25) Maiti, A.; Brabec, C. J.; Bernholc, J. Phys. ReV. Lett. 1993, 70, 3023. (26) Townsend, S. J.; Lenosky, T. J.; Muller, D. A.; Nichols, C. S.; Elser, V. Phys. ReV. Lett. 1992, 69, 921. (27) Rosato, V.; Celino, M.; Gaito, S.; Benedek, G. Comput. Mater. Sci. 2001, 20, 387. (28) Rosato, V.; Celino, M.; Benedek, G.; Gaito, S. Phys. ReV. B 1999, 60, 16928.

Petersen et al. (29) Hyde, S. T.; O’Keefe M. Philos. Trans. R. Soc. London, Ser. A 1996, 354, 1999. (30) Phillips, R.; Drabold, D. A.; Lenosky T.; Adams G. B.; Sankey O. F. Phys. ReV. B 1992, 46, R1941. (31) Huang, M. Z.; Ching W. Y.; Lenosky, T. Phys. ReV. B 1993, 47, 1593. (32) Aoki, H.; Koshino, M.; Takeda, D.; Morise, H.; Kuroki, K. Phys. ReV. B 2001, 65, 035102. (33) Park, N.; Yoon, M.; Berber, S.; Ihm, J.; Osawa, E.; Tomanek, D. Phys. ReV. Lett. 2003, 91, 237204. (34) Valencia, F.; Romero A. H.; Herna´ndez, E.; Terrones, M.; Terrones, H. New J. Phys. 2003, 5, 123.1. (35) Bourgeois, L. N.; Bursill, L. A. Philos. Mag. A 1999, 79, 1155. (36) Harris, P. J. F. Philos. Mag. 2004, 84, 3159. (37) Petersen, T. C.; Snook, I. K.; Yarovsky, I.; McCulloch, D. G. Phys. ReV. B 2005, 72, 125417. (38) Jenkins, G. M.; Kawamura, K. Nature 1971, 231, 175. (39) O’Malley, B.; Snook, I.; McCulloch, D. Phys. ReV. B 1998, 57, 14148. (40) Petersen, T. C.; Yarovsky, I.; Snook, I. K.; McCulloch D. G.; Opletal, G. Carbon 2003, 41, 2403. (41) Opletal, G.; Petersen, T. C.; O’Malley, B.; Snook, I.; McCulloch, D. G.; Marks, N.; Yarovsky, I. Mol. Simul. 2002, 28, 927. (42) Petersen, T.; Yarovsky, I.; Snook, I.; McCulloch, D. G.; Opletal, G. Carbon 2004, 42, 2004. (43) McGreevy, R. L.; Pusztai, L. Mol. Simul. 1988, 1, 359. (44) Justo, J. F.; Bazant, M. Z.; Kaxiras, E.; Bulatov, V. V.; Yip, S. Phys. ReV. B 1998, 58, 2539. (45) Marks, N. A. Phys. ReV. B. 2000, 63, 035401. (46) Wohlgemuth, M.; Yufa, N.; Hoffman, J.; Thomas, E. L. Macromolecules, 2001, 34, 6083. (47) Harper, P. E.; Gruner, S. M. Eur. Phys. J. E 2000, 2, 217. (48) Schwartz, H. A. Gesamemelte Mathematische Abhandlungen; Springer: Berlin, 1890. (49) Holyst, R.; Plewczynski, D.; Aksimentiev, A.; Burdzy, K. Phys. ReV. E 1999, 60, 302. (50) Plewczyn´ski, D.; Holyst, R. Physica A 2001, 295, 371. (51) Hyde, S. T. J. Phys. Chem. 1989, 93, 1458. (52) Gozdz, W. T.; Holyst, R. Phys. ReV. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1996, 54, 5012. (53) McCulloch, D. G.; McKenzie, D. R.; Goringe, C. M.; Cockayne, D. J. H.; McBride, W.; Green, D. C. Acta Crystallogr. 1999, A55, 178. (54) Franzblau, D. S. Phys. ReV. B 1991, 44, 4925. (55) Opletal, G.; Petersen, T. C.; McCulloch, D. G.; Snook, I.; Yarovsky, I. J. Phys.: Condens. Matter 2005, 17, 2605.