Theoretical characterization of several models of nanoporous carbon

Elastic, electronic and vibrational properties of seven models of nanoporous carbon are reported. The studied structures are periodic graphitic arrangements with heptagonal and octagonal rings of carbon, known as Schwarzites. The calculations were performed within a non-orthogonal tight binding framework which has been shown to be reliable for diamond, graphene layers, fullerenes and carbon nanotubes. In contrast with previous studies, each structure was properly relaxed, so that differences between each model must be assigned to intrinsic properties rather than to differences in their construction. Thermodynamic properties were calculated from the vibrational density of states.


Introduction
The production and characterization of porous forms of carbon have been, for quite a long time now, active fields of research, mainly in connection with their applications as active filters and gas storage devices. While traditional porous carbons exhibit rather disordered pore structures, the recent implementation of template techniques for the fabrication of nanostructured materials [1]- [3] brought a new perspective to the field. In fact, several nanoporous carbons with periodic pore structure have been synthesized during the past few years [2]- [10], and have revealed themselves as very good candidates for electrochemical double-layer capacitors (EDLCs), hydrogen storage and electro-catalytic devices, being, therefore, appealing materials for digital communication devices and fuel-cell technologies.
Zeolites (crystalline aluminosilicate materials with pore sizes ranging from 3 to 13 Å) and mesoporous silicates (amorphous solids with ordered pore structures) have been used as molecular sieves, and it has been shown that their pore sizes can be suitably controlled within the fabrication procedure. Impregnation of these porous structures with carbon seems to be a promising method for synthesizing a great variety of new periodic graphitic arrangements with mainly sp 2 bonding. The template (zeolite or silicate) can be removed using hydrofluoric acid (HF), leaving the carbon network intact [7]. After removal of the template material, two different behaviours have been observed: some materials retain the same structure as before sieve elimination, and others develop structural transformations. In every case, the final result is a thermally and mechanically stable periodic carbon material with ordered pores. Possible modifications induced on carbon networks by zeolite frameworks have been discussed in the literature [19]. In particular, in the work of Terasaki et al on zeolite [7], the very small size of the zeolite pores, the lack of evidence for graphitic planes after thermal treatment up to 900 K, and the strong evidence of sp 2 bond hybridization given by the solid state NMR spectra support the idea of ordered aromatic rings conforming curved surfaces as structural blocks of the resulting nanoporous carbon.
Although a few decades ago carbon chemistry was thought to be completely understood in terms of the sp 2 graphene hybridization and sp 3 diamond hybridization, the discovery of nanostructured carbon forms (fullerenes [11], nanotubes [12], nano-onions [13], etc) brought attention to the rich variety of chemical and physical properties related more strongly to the 123.3 geometry, and particularly to the curvature of carbon structures, which imply intermediate degrees of hybridization. In order to build arrangements of positively curved surfaces, such as carbon fullerene crystals [11] and silicon clathrates [14], from sp 2 and sp 3 correlated networks (i.e. with group IV atoms), the inclusion of pentagonal rings is needed. The existence of pentagonal rings poses its own signature over the electronic and vibrational properties and the vast amount of void space in fullerite crystals, along with their cagelike structure, makes them suitable for endohedral doping procedures, allowing for further tuning of their properties [15].
It is also possible to build up crystals entirely from sp 2 hybridized carbon, showing a negatively curved surface [16]- [18]; this type of material can be obtained by including sevenmembered and eight-membered rings. Following the Mackay-Terrones proposal, several models of structures similar in topology to triply periodic minimal surfaces (i.e. mean curvature equal to zero) have been proposed [22]- [24]. In these periodic arrangements, called Schwarzites, the negative Gaussian curvature is obtained by the inclusion of heptagonal rings of carbon. Schwarzites divide space into two regions possessing interconnected channels with pores similar to those found in zeolites. Relationships between zeolites and carbon networks have also been explored in the literature in the context of fullerite crystals [20].
Previous work on Schwarzites has shown that some of them may be stable under normal experimental conditions, making them good candidates for crystalline carbon materials with porosity at the mesoscopic or nanoscopic level [10]. Carbon Schwarzite phases have already been identified and are thought to play a major role in the conformation of non-graphitizing carbons [25]. Schwarzite-like structures may therefore be present in ordered nanoporous materials such as those obtained by Ma and co-workers [7], and also in non-graphitizing meso and nanoporous carbons, such as those reported in [26].
In an attempt to elucidate the actual structure and physical properties of nano and mesoporous carbons, theoretical investigations are needed. In this work, we report a comprehensive characterization of some of the negative curvature structures, with cubic unit cells, both from the mechanical and electronic points of view. Elastic constants, which provide valuable information about mechanical stability and electronic and vibrational densities of states (DOSs) have been calculated, providing a signature for each model structure. We also report on the thermodynamical properties as calculated from the vibrational DOS. In addition, we have found that, besides the curvature of the system, proper relaxation of the structures is essential for an adequate description of the electronic properties, even at a qualitative level.
Due to the large number of atoms per unit cubic cell in the studied structures, most of our results rely on a tight binding method. But whenever we found it necessary to reinforce our conclusions-and it was computationally feasible-plane wave density functional theory (DFT) calculations in both the local density approximation (LDA) and generalized gradient approximation (GGA) to the exchange-correlation energy were also performed using the PWSCF package [27] and the CASTEP package [28].

Studied models
Regarding the notation of the structures studied, interested readers are referred to the original references [16]- [18], [22]- [24]; specially, for a full review, see [30] and [31,33,32]. As stated above, several periodic graphitic structures with negative Gaussian curvature have been proposed in the literature [16]- [18], [22]- [24], [29]. These structures are classified as diamondlike (D), gyroid (G), and primitive (P), according to the symmetry of the network formed by the 123.4 interconnected channels [30,16]. The D, G and P, triply periodic minimal surfaces (TPMSs) divide the space in two equivalent subspaces, and for this reason they are also referred as balanced (BAL) structures. These TPMS can be decorated with carbon atoms, thus obtaining the D8BAL, G8BAL, and P8BAL structures. On the other hand, there are periodic graphitic structures in which each carbon atom is shared by two octagons and one hexagon (type 688 vertices), thus generating the arrangements D688, G688, and P688. Finally, it is also possible to build up periodic surfaces which are parallel (PAR) to a balanced one and, therefore, divide the space into two inequivalent subspaces.
In the present contribution, we studied three structures of type 688 (D688, G688, and P688), three balanced structures containing octagonal rings (D8BAL, G8BAL, and P8BAL), and one parallel structure of the primitive network containing seven-membered rings (P7PAR).

Calculation method
In this work we use a tight binding model [34] to calculate the relaxed structures, electronic structure, mechanical properties and vibrational DOS of seven Schwarzite models. Specifically, we use the non-orthogonal tight binding model due to Porezag et al [35], as implemented in the TROCADERO [36] code.
First of all, we determined the equilibrium structures relaxing the lattice constant and atomic coordinates of each model studied. The k-point meshes used in this step were chosen such that the energy changes for the initial unrelaxed models associated with increasing k-point mesh sizes lay below 1 meV/atom, which is enough to follow accurately the changes in energy with volume and geometry. For the smallest cells (D688, P688, G688, and D8BAL) this accuracy threshold was reached by using four k points, corresponding to a 2 × 2 × 2 mesh in the Monkhorst-Pack generation scheme [37]. For the largest cells -point sampling proved to be accurate enough. Fitting the energy versus volume curves to the Birch-Murnaghan formula [38] we found the equilibrium energy, E 0 , lattice parameter, a = (V 0 ) 1/3 , bulk modulus, B 0 , and the derivative of the bulk modulus with pressure, B p , for each structure. For the studied crystals the first three terms in (1) are enough to obtain a good fit. For cubic systems there are three independent elastic constants: C 11 , C 12 , and C 44 , this last one usually referred to as the trigonal shear modulus. The first two elastic constants can be written in terms of the bulk modulus, B 0 = (C 11 + 2C 12 )/3, and the tetragonal shear modulus, (C 11 −C 12 ). As stated above, the bulk modulus is extracted from the Birch-Murnaghan equation. The remaining elastic moduli were calculated by straining the lattice vectors and relaxing the internal coordinates. For the tetragonal shear modulus we used volume conserving orthorhombic deformations [38], which for cubic structures produce energy changes given by where x is the strain. For the trigonal shear modulus, volume conserving monoclinic deformations with energy changes of the form 123.5 were used. The deformation parameter, x, used in our calculations ranges between 0 and 0.06 for both moduli. Electronic DOSs were calculated from the eigenvalues of the tight binding Hamiltonian, using a Lorentzian broadening of 0.1 eV. The number of k points used in the DOS calculations was larger than the one used for the geometrical optimizations, in order to obtain smooth curves. 256 k points, corresponding to the 8 × 8 × 8 mesh in the Monkhorst-Pack scheme, were used for the D688, P688, G688, and D8BAL structures, and 32 points (4 × 4 × 4 Monkhorst-Pack mesh) for the remaining models.
As discussed below, we have also carried out DFT calculations on the energetics and DOSs for P688 and G688 to test the tight binding results. Within a plane wave pseudopotential framework [27] we performed both LDA calculations using the Perdew-Zunger [39] exchangecorrelation functional with a Von Barth-Car [40] pseudopotential, and GGA calculations with the Perdew-Wang [41] exchange-correlation formula and an ultrasoft [42] pseudopotential. Due to the increase in the computational effort in the DFT calculations we have chosen an accuracy goal of 0.05 eV/atom in the total energy, which is sufficient as long as one is not interested in detailed geometrical optimization. This accuracy threshold was reached with a kinetic energy cut-off of 58 Ryd and four special k points for the Brillouin zone integration (the Monkhorst-Pack 2 × 2 × 2 mesh) for the LDA calculation, and with an energy cut-off of 20 Ryd and ten k points for the GGA calculation.
Vibrational DOSs were calculated constructing and diagonalizing the dynamical matrix and considering replicas of the unit cell in the case of small cells and a single unit cell for larger cells (which give only the vibrational frequencies at the point). From these v-DOSs, a full set of thermodynamical values was also computed.

Geometry and energetics
Sketches of the structures obtained after relaxation can be viewed in figure 1. Considering that the studied models have been taken from several reported structures, it is important to relax them using the tight binding energy and forces, and to compare the resulting structures with the original ones. The relaxation procedure, implemented using the conjugate gradients algorithm, preserved the structure in all cases, though bond lengths and angles were modified, as can be seen in table 1, where changes both in values and in distribution of bond lengths/angles can be easily seen.
There is no systematic trend from the unrelaxed to the relaxed structures, as should be expected, given that the unrelaxed structures are constructed by several different methods. The noticeable changes seen in table 1 reflect the differences between those methods, specially the relaxation procedure employed-or lack of it, making each model sensitive at a different level to the tight binding relaxation.
The tight binding model we have used produces rather good geometries both for diamond and graphene, the calculated lattice parameters being 3.57 and 2.456 Å respectively, comparing very well with the experimental data. The bonds and angles in the carbon Schwarzite structures are rather close to that of a graphene layer, as can be seen from the bond lengths and angles  table 2 are measured with respect to the graphene structure. As can be seen in the same table, diamond is only slightly above graphene in energy. The Schwarzite structures lie between 0.23 and 0.53 eV above graphene or diamond. The energies calculated for the P688 and polybenzene (D688) structures compare well with those reported in the work of O'Keeffe et al [22].
Of the studied structures, P8BAL shows the lowest energy, 0.236 eV/atom, although the difference with polybenzene (D688) is only about 0.026 eV/atom. The C 60 buckyball lies about 0.44 eV/atom above the graphene layer (see for example [22]), so that the studied models are at least energetically viable. The P688 structure exhibits the larger energy difference with respect to graphene (0.52 eV/atom) and this and the G688 structures are the only ones with energies above that of the C 60 structure. It is also remarkable that the two studied D structures, despite their different curvatures, atomic densities and topology, lie quite near in energy. This behaviour is not shared by the G and P structures, which exhibit rather different energetic behaviour when the curvature and density are changed. Nevertheless there does not seem to be a direct relationship between the angular dispersion and the energy of the structures. It is worth mentioning that the D8BAL, the G8BAL, and the P8BAL possess relatively close energies, so possible phase transformations involving these arrangements must not be ruled out (activation barriers will be discussed in subsequent work). In fact, there is a cyclic mathematical transformation, called the Bonnet transformation-the transformation that maps a catenoid into a helicoid-which transforms the P minimal surface into the G minimal surface and the G minimal surface into the D minimal surface without changing the metric [17,18,21] and angles, thus preserving the Gaussian curvature and the topology (all these surfaces possess genus three per primitive unit cell). Therefore, it can be proposed that atoms on these surfaces (D, G, and P) can follow this transformation, changing the D8BAL into the G8BAL and the G8BAL into the P8BAL.

Elastic properties
Bulk and shear moduli calculated in the present work are shown in table 3. We note that the agreement between the calculated values for diamond and the experimental data (B 0 = 543 GPa, C 11 − C 12 = 951 GPa, C 44 = 576 GPa) lies well within the accuracy of typical LDA results. In polycrystalline samples neither the tetragonal nor the trigonal shear modulus can be directly measured, but the isometric shear modulus G can. Although there is no simple relationship between the elastic constants and G, bounds can be calculated from them [38] and we have also included those bounding values in table 3. It is also customary to characterize the mechanical properties of a given material in terms of the Young modulus E y and Poisson ratio ν, which can be related to the bulk and shear moduli [38]: and Taking the value of G as the mean between the obtained lower and upper bounds, we have estimated those parameters, as well as the first Lame constant I , and included them in table 3. Compared to diamond all the Schwarzites are, of course, softer, but still strong: the bulk moduli are larger than, for example, the value for pure silicon (98 GPa and, typically, between 80 and 100 GPa in LDA calculations). All the shear moduli are larger than zero, meaning the structures are mechanically stable.

123.9
As expected, the bulk modulus increases with the atomic density. As a matter of fact, we have seen that B 0 follows an approximate power law B 0 = C * ρ λ , where ρ is the atomic density of each model, and the exponent λ equals 1.071. Nevertheless, D688 and D8BAL exhibit B 0 which are slightly larger than calculated from the above-mentioned power law.
One of the main features is the rather small shear moduli, compared to bulk moduli, meaning that these materials would be much softer under shear than under compression. This can be thought of as a natural feature of porous and layered materials, as some parts of the structure can slide easily thanks to the large amount of void space within them.
Of all the studied models, polybenzene exhibits the largest bulk and shear moduli and, therefore, it would be the most stable under mechanical stress. For this reason, and its relatively low energy, discussed above, polybenzene could be thought of as the major candidate to be synthesized, in agreement with the conclusions of others [22].

Electronic properties
In table 4 the calculated bandgaps are shown. Those gaps are calculated as the difference between the highest occupied energy level and the lowest unoccupied energy level in the k-point mesh used to sample the Brillouin zone. The band gap for diamond is remarkably overestimated in our calculations (about 42% over the actual value of 5.47 eV), a feature that is related to the minimal basis set used in the tight binding model; the same trend is to be expected for the remaining values. It should be noted here that even first principles DFT calculations can lead to significant errors in bandgap estimation (up to 30% underestimation or more), and therefore similar errors in more simplistic methodologies such as tight binding (even if of the opposite sign) are to be considered in this context. Proper determination of the electronic characteristics is better achieved using the gap information in conjunction with the electronic DOS.
Electronic DOSs calculated in this work are presented in figure 3. We can identify four structures as semiconductors: D688, G688, P8BAL, and P7PAR. The remaining three structures, P688, G8BAL, and D8BAL, can be identified as gapless systems with low electronic DOS at the Fermi level. It is also interesting to mark the similarities between the electronic structure of the D8BAL and G8BAL. Besides being gapless systems with low e-DOS at the Fermi level, further correspondence can be traced between some of the major characteristics. There is a strong and narrow peak around 8-9 eV, a broader peak near 7 eV and a third one around 3 eV, although in D8BAL this third peak appears to be segmented. These three peaks also seem to be present in the P8BAL structure, but the correspondence is less clear, except for the 9 eV peak which appears clearly. Interestingly, from the electronic properties we also see similarities among the D8BAL, G8BAL, and P8BAL structures, which might be explained by the Bonnet transformation [17,18].
A comparison of the tight binding e-DOS with the GGA e-DOS for the G688 structure is shown in figure 4. Despite its deficiencies, mainly the clear overestimation of the bandgap, the tight binding model describes reasonably well the shape of the e-DOS both in the valence and in the conduction region and, specially, the major changes in the density coming from the geometry change: the deepening of the valence band valley around 3 eV, the reduction in the valence state density near the Fermi level and the rightward shift of the first hump in the conduction state density. It is also to be noted that our GGA results underestimate the diamond bandgap by about 25%, so that DFT results in figure 4 may also underestimate the semiconductor character of the relaxed G688 structure. Similar results were obtained for the P688 structure.    As stated by Huang and Ching [29], the electronic properties are more strongly dependent upon global geometrical properties than on the exact atomic positions. In fact, most of the changes led to by the relaxation procedure are small: besides the reinforcement of the semiconductor character (or weakening of the metallic character) the relaxation does not usually bring about dramatic changes. Still, proper relaxation plays a relevant role in the theoretical characterization of these structures because even small atomic rearrangement can induce changes in the global geometries. In the P688 structure after relaxation the hexagonal rings, originally occurring in quite flat planes, get distorted. Subtler changes take place in G688: there is a change in the alternance of the bond orders (see table 1

Vibrational properties and thermodynamics
The vibrational DOSs computed in this work are shown in figure 5. The most remarkable feature in these v-DOS spectra is the peaks-and-gaps structure beyond 1500 cm −1 in the smaller cells, although the 10 cm −1 Lorentzian broadening used tends to somewhat obscure this behaviour. For D688 there is a clear gap centred at 1670 and of 85 cm −1 width. For P688 more gaps appear: the first one centred at 1535.5 and 50 cm −1 wide, another one at 1642.5 and 88 cm −1 wide, a third one centred at 1767 and 141 cm −1 wide, the fourth is centred at 1817 and 41 cm −1 of width, and finally a fifth, very narrow, gap appearing centred at 1863 cm −1 and only 3 cm −1 wide. For G688 there are no longer gaps in the investigated region. For the larger structures it is difficult to establish the existence of gaps as long as we are limited to the vibrations at the gamma point. These vibrations, nevertheless, seem to fill rather densely the range of allowed wavevectors and these structures are unlikely to develop vibrational gaps. As is to be expected for sp 2 hybridized networks, all these structures exhibit rather high v-DOS near the graphene peak (around 1540 cm −1 ); nevertheless, this feature is no longer the most relevant in the v-DOS spectra. Now we turn our attention to the thermodynamic potentials for the Schwarzite structures. In table 5 we have reported the calculated values of the ionic (vibrational) contribution to the total energy, E, free energy, F, entropy, S, and specific heat, C, at 300 K, showing very similar values. In particular, the G688 and D688 have the same values, while the P7PAR and the P8BAL have almost the same, exhibiting the highest entropies and specific heats and the lowest free energies.     Figure 6. Temperature dependence of the thermodynamic variables for the polybenzene (D688) structure. E is the vibrational contribution to the total energy, F is the free energy, S is the entropy and C is the specific heat.
The dependence of the thermodynamic potentials on temperature for the polybenzene structure is sketched in figure 6. For the remaining structures the behaviour of the potentials as a function of temperature is very similar. All thermodynamic quantities of the studied structures behave as functions of temperature very much like other carbon structures such as diamond and graphite. The heat capacity indicates that the D688 structure behaves as a Debye crystal at about 1000 K.

Summary and conclusions
We have made a thorough study of seven Schwarzite structures. In particular, we have relaxed the studied models using a tight binding method, the results of which are confirmed by plane wave LDA and GGA calculations for the smallest cells. Five of the seven models are closer in energy to graphene than the C 60 fullerene and the remaining two exhibit the same order of energy, meaning that these models are at least energetically viable. The mechanical stability of the structures is demonstrated by the fact that the bulk and shear moduli are all positive and of the same order as those of silicon.
Regarding electronic properties, we have found four of the studied models to behave as semiconductors, although the calculated gaps are likely to be overestimated, so they could be considered small gap semiconductors. The other three are gapless systems with relatively low e-DOS at the Fermi level. In the studied cases, the relaxation procedure results in the reinforcement of the semiconductor character or weakening of the metallic behaviour. The vibrational spectra of the polybenzene and P688 structures shows rather clear gaps, a feature that does not seem to be shared by the remaining models. Despite this difference, the thermodynamical variables studied (total energy, free energy, entropy and specific heat) are very similar for all structures. The differences in the peak structure beyond 1500 cm −1 between D688, P688 and G688 structures may be used to discriminate between them via Raman measurements. For the larger structures the absence of clearly defined peaks in the spectra would make it difficult to discern between those Schwarzite structures and amorphous phases.