Carbons with Regular Pore Geometry Yield Fundamental Insights into Supercapacitor Charge Storage.

We conduct molecular dynamics simulations of electrical double-layer capacitors (EDLCs) using a library of ordered, porous carbon electrode materials called zeolite templated carbons (ZTCs). The well-defined pore shapes of the ZTCs enable us to determine the influence of pore geometry on both charging dynamics and charge storage mechanisms in EDLCs, also referred to as supercapacitors. We show that charging dynamics are negatively correlated with the pore-limiting diameter of the electrode material and display signatures of both progressive charging and ion trapping. However, the equilibrium capacitance, unlike charging dynamics, is not strongly correlated to commonly used, purely geometric descriptors such as pore size. Instead, we find a strong correlation of capacitance to the charge compensation per carbon (CCpC), a descriptor we define in this work as the average charge of the electrode atoms within the coordination shell of a counterion. A high CCpC indicates efficient charge storage, as the strong partial charges of the electrode are able to screen counterion charge, enabling higher ion loading and thus more charge storage within the electrode at a fixed applied voltage. We determine that adsorption sites with a high CCpC tend to be found within pockets with a smaller radius of curvature, where the counterions are able to minimize their distance with multiple points on the electrode surface, and therefore induce stronger local partial charges.

Introduction

Energy storage is becoming ever more important as society moves away from fossil fuels toward a cleaner energy paradigm. Storage technologies are often compared using the high-level metrics of energy density and power density. Batteries have relatively high energy density and low power density, while traditional capacitors have high power density and low energy density.[1,2] One advantage of capacitors is that, since they store energy through electric charge, they do not degrade as easily as batteries, which store energy through a chemical reaction. Capacitors are thus able to provide millions of charge–discharge cycles, versus several thousand cycles for the best-performing batteries.[3] Motivation to increase the energy density of capacitors has spurred the development in the past few decades of a new type of capacitor, called the electrical double layer capacitor (EDLC), or supercapacitor. EDLCs are differentiated from traditional capacitors by having a liquid electrolyte in place of a solid dielectric. Under an applied voltage, the accumulation of positive and negative ions at the anode and cathode, respectively, leads to the creation of a double layer at each electrode.

The energy density of a capacitor is proportional to the capacitance, defined as the amount of charge, Q, stored by a capacitor at a given voltage, V(3)For an idealized parallel plate capacitor consisting of metallic electrodes enclosing a dielectric medium of uniform dielectric constant, ε, the capacitance can be derived exactly aswhere A is the cross-sectional area of the capacitor and d is the separation distance between the electrodes. If we transpose this expression to an EDLC, A would be the surface area of the electrode which is accessible to electrolyte ions, and d the characteristic thickness of the electrode–electrolyte interface. Equation provides some intuition for why EDLCs have increased capacitance over their traditional counterparts: First, the charge separation distance is smaller in an EDLC than in a traditional capacitor, since ions can approach within less than 0.5 nm of the electrode surface, while in a traditional capacitor d is a few nanometers or higher.[3,4] Second, the accessible surface area of an EDLC can be increased by several orders of magnitude using rough or porous electrode.

Multiple theories have arisen to describe the electrochemical double layer, beginning with the classical theory of Helmoltz,[5,6] which was subsequently improved by Gouy,[7] Chapman,[8] and Stern[9] to consider discrete ions and complex double layer structures. These theories can accurately predict capacitances of EDLCs whose pores are macroscopic. However, when the pores become comparable in size to the electrolyte ions, so-called “anomalously” high capacitances have been observed that break with both existing theories and empirical trends.[10,11] Reports of materials with such impressive capacitances have led to considerable growth in the field of microporous materials for EDLCs, to better understand the mechanisms behind capacitance in small pores.[12−18]

A popular choice of material for EDLC electrodes is porous carbon, due to its stability, ease of synthesis, and low cost. Porous carbons used in EDLCs include activated carbons,[11] carbide-derived carbons,[10,19] carbon onions and nanotubes,[20] carbonized precursors such as metal–organic frameworks,[21] and graphene-based composites.[22,23] Experiments and simulations have shed some light on the charge-storage mechanisms in such materials;[14−16,24,25] however, a major challenge is that most microporous carbon materials have neither a narrow pore size distribution nor a regular structure with long-range order, making it difficult to draw conclusions between structure and performance.

A new class of materials called zeolite-templated carbons (ZTCs), which are synthesized using a sacrificial zeolite scaffold,[26−28] has been demonstrated as a promising EDLC material.[29−32] Thus far, ZTCs have been synthesized from just three zeolites (FAU, EMT, and beta) of the 245 frameworks recognized by the IZA Structure Commission.[33−35] Recently, Braun and co-workers reported a method to computationally synthesize ZTCs from a given zeolite structure, and then predict whether the ZTCs are experimentally synthesizable based on stability during finite-temperature molecular dynamics simulations.[36] Their predicted ZTCs are composed of sp2-hybridized carbons that tile a surface that is dual to the zeolite. Templating on a crystalline framework confers well-defined pore geometries that could yield insights into the structure–property relationships of electrode materials, motivating further study of ZTCs for energy storage applications.

In this work, we use molecular dynamics (MD) simulations to screen the ZTC materials of Braun et al. as electrode materials in EDLC cells. We show that the charging time scale of the ZTCs is negatively correlated with pore limiting diameter, and that there is evidence of both progressive charge penetration and kinetic trapping within the ZTCs during charging. We then study the equilibrium capacitance of the ZTCs to investigate the correlation between geometric descriptors, local electrolyte configurations, and charge storage mechanisms within the electrode. Introducing the concept of charge compensation per carbon (CCpC), we find that charge storage is more efficient at ion adsorption sites with high CCpC, which are more likely within pores with a lower radius of curvature. Conversely, charge storage is diminished at high-radius-of-curvature sites and within sites with a mismatch of local pore diameter and ion size.

Results and Discussion

As demonstrated in the Supporting Information, we determined suitable protocols for building an EDLC simulation cell and tested different MD equilibration schemes. Here we present the results of our simulations using the constant-potential method for charge equilibration, discussing first the dynamics of charging and then equilibrium capacitances.

Charging Dynamics

An EDLC cell can be represented macroscopically with an equivalent circuit model,[3] which can range in complexity from an RC circuit, consisting of a capacitor in series with a resistor (representing the electrode capacitance and solution resistance, respectively), to more detailed representations such as a transmission line model, which assumes that the properties of the circuit are distributed continuously throughout the material and thus models the EDLC as a series of individual circuit elements, or electrode “slices.”[37] These more complex models, such as the transmission line model used by multiple groups,[38,39] are especially useful when extrapolating macroscopic charging dynamics from microscopic simulations. Since we only want estimates of the relative charging dynamics among different ZTC materials, we use here the simplest RC equivalent circuit (effectively a one-slice model) to model the ZTC EDLCs.

Under an applied constant potential, the accumulated charge in an RC circuit exponentially approaches its equilibrium value.[40] We therefore fit the transient charge per atom q(t) of the EDLCs to the exponential functionwhere q∞ is the infinitely equilibrated charge per atom of the capacitor and τ is the time constant of the exponential, which for an RC circuit is R × C (the product of the resistance and the capacitance of the circuit). The circuit is considered to be equilibrated after 5τ, at which point the charge is within 1% of the infinitely equilibrated charge.

Examples of the exponential fits for the charging of FAU_1, BEA, and EMT are shown in Figure . For all of the structures, the electrode atom charges are underestimated by the exponential fit during the first 1 to 2 ns of constant-potential charging (SI Figure S10). Defining an “average residual” ε1τ to be the relative root-mean-square error (RMSE) between the exponential fit and the charging profile over the first characteristic time, we plot the characteristic time τ against the pore limiting diameter (PLD) of the material, with the points colored to show ε1τ, in Figure . We observe that τ is negatively correlated with PLD, consistent with work from Vasilyev et al.,[41] on CDCs, indicating that it takes longer for charge to equilibrate in materials with smaller PLD, likely due to larger diffusive barriers. Similarly, ε1τ also increases with smaller PLD. Since ε1τ reflects the quality of the exponential fit of the accumulated charge (higher ε1τ indicates worse fit), the data show that the exponential fit at early times is more severely underestimated when PLD decreases. This underestimation of the initial charging profile by the exponential fit suggests the existence of multiple time scales in EDLC charging, due to both progressive charging of the electrodes from the bulk electrolyte inward (consistent with a transmission line model), as well as “overfilling” and ion trapping analogous to that observed by Kondrat et al.[17,42] Both phenomena, as well as further details regarding equilibrium convergence, are discussed further in the Supporting Information.

Figure 1

Evolution of average absolute charge per atom for FAU_1, BEA, and EMT, with a constant potential drop of 1 V applied to the cell starting at t = 0 ps. Exponential fit is shown in the red dashed line and fit parameters are provided in the label.

Figure 2

Characteristic time scale τ of EDLCs during charging at ΔΨ = 1 V, colored with the normalized ε1τ of the exponential fit during the first characteristic time scale.

Capacitance Screening

The capacitance under a constant applied potential is plotted against various geometric descriptors (top row) and average local interfacial properties (bottom row) in Figure . While it would be interesting to compare the values obtained in this work with previous results obtained for disordered porous structures, we limit ourselves to the systems simulated here since the design of the model supercapacitors and the equilibration process seem to affect the simulated capacitances. Indeed, previous simulations done with constant charge or constant potential equilibrations and with slightly different starting densities have shown a difference of 36% between calculated capacitances.[14,38] Geometric descriptors, in particular, average pore size and pore size distributions, are often used to rationalize structure–capacitance relationships in EDLCs.[43,44] For example, In Käärik et al., the authors trained a regression model to link capacitance to three experimental descriptors: surface area, pore size fraction, and density. While they showed overall statistical agreement of their model with most of the 100 amorphous carbons they studied, their regression model has the highest relative error (roughly 20%) for the highest-capacitance materials, suggesting that these geometric descriptors lack a fundamental correlation to the anomalous increases in capacitance achievable with subnanometer pores.[45] Other theoretical studies on slit nanopores even suggest that the capacitance could follow an oscillatory behavior with the average pore size for carbons with well-defined pore sizes.[46−48] However, in our materials we do not see a definite correlation with any of the geometric descriptors alone. Indeed, for all purely geometric descriptors the Pearson correlation coefficient (Pearson’s r) with gravimetric capacitance is less than 0.5 (Table ).

Figure 3

Gravimetric capacitance as a function of geometric descriptors in the top row: density, void fraction, average pore size, and PLD; and local properties in the bottom row: average distance between counterions and electrode (⟨d̅sep⟩), the surface area to charge separation ratio (A/⟨d̅sep⟩), the average degree of confinement (⟨DoC⟩), and the average charge compensation per carbon (⟨CCpC⟩). The inset plots the gravimetric capacitance against total charge compensation in the coordination shell, where the x-axis range goes from 0.15 e to 0.23 e.

Table 1

Pearson Correlation Coefficients Computed between Gravimetric Capacitance and Descriptors Shown in Figure

descriptor	Pearson’s r
Geometric Properties
density	–0.385
void fraction	0.461
dpore	–0.005
PLD	0.039
Local Interfacial Properties
⟨d̅sep⟩	–0.206
A/⟨d̅sep⟩	0.751
⟨DoC⟩	–0.387
⟨CCpC⟩	0.934
total charge compensation	0.252

This suggests that while particular geometric descriptors might be a useful indicator of capacitance within particular families of materials, a clear relationship between capacitance and, for example, pore size is not the rule, but rather the exception for materials that are otherwise geometrically similar. ZTCs, due to their regular templated structures, exhibit a diversity of topologies, pore geometries, and local curvatures, which are not well captured by traditional geometric descriptors, but are known to influence charge storage.[15] Thus, the insights we can glean from local interfacial properties in ZTCs might be better translated to microporous carbon materials in general.

The bottom row of Figure plots capacitance versus quantities related to the electrolyte–electrode interfacial configuration, which are computed for an ion in relation to the electrode atoms within its coordination shell. The charge separation (⟨d̅sep⟩) is the average distance between the counterion and the carbons within its coordination shell. The degree of confinement (DoC) is defined as the fraction of the maximum solid angle around a counterion that is occupied by carbon atoms within the coordination shell cutoff (set to the first minimum in the ion-carbon radial distribution function).[15] And finally, the charge compensation per carbon (CCpC), a quantity introduced in this work, is defined as the magnitude of the net charge of the coordination shell divided by the number of atoms within the coordination cutoff. Since the coordination shell represents the “sphere of influence” of an ion, a high CCpC indicates that a relatively large amount of the ionic charge is being compensated by a relatively small number of electrode atoms. For all quantities, the angle brackets ⟨⟩ denote averaging over all counterions in an electrode.

A positive correlation can be observed between the capacitance and A/⟨d̅sep⟩ (r = 0.751), reminiscent of classical theories of capacitance as expressed by eq . This suggests that we can view capacitance in the ZTCs as arising from an “ideal” contribution from a reference electrode with the same A/⟨d̅sep⟩, and a “non-ideal” contribution responsible for the deviations from classical double layer theory, arising from the microporosity. One measure of how micropores influence charge storage is the DoC. Here, we note that we are plotting in Figure the average degree of confinement, ⟨DoC⟩, which obscures differences in the range and distribution of confinement values within a material. We do not observe a strong correlation with capacitance when ⟨DoC⟩ is below 0.25, and when ⟨DoC⟩ is above 0.25 the capacitance seems to be slightly negatively correlated with confinement. This finding adds nuance to the conclusions from previous studies that more confinement generally has a positive influence on charge storage efficiency, i.e., that ions in highly confined sites generate a higher counter-charge on the electrode.[15,18] We discuss confinement effects further in a later section, where we examine charge storage mechanisms in individual pores.

Finally, the local descriptor that appears to have the best correlation with capacitance is ⟨CCpC⟩, for which we observe a positive and nearly linear relationship (r = 0.934) with even less scatter than for A/⟨d̅sep⟩. Capacitance and ⟨CCpC⟩ both aggregate information about the charge stored by the electrode atoms; however, their strong correlation is not trivial because only about 30–45% of the electrode atoms are within the coordination shell of a counterion at a given time step. These coordination shell carbons have a slightly larger-than-proportional share of charge, carrying between 35% and 50% of the net charge in the electrode (SI Figure S13). Perhaps surprisingly, the capacitance does not correlate with the total charge compensation within the coordination shell (inset of Figure ). The observation that per-carbon (rather than total) charge compensation correlates so strongly with the capacitance indicates that electrodes store charge most efficiently when the local environment of an ion can compensate a large fraction of its charge while using relatively little electrode “real estate”.

One complication with comparing materials using local properties is that they are computed with a definition of the coordination shell that uses a cutoff radius, rcut around the ion. rcut radius was chosen following the literature[15] as the first minimum in the ion-carbon RDF. However, we found in our materials that the first minima were not all at the same location in all materials, and some of them did not have a clear “minimum” at all. Therefore, we opted to use the same rcut of 6.3 Å for all materials, as this was the location of most of the RDF first minima and also was consistent with the literature. Further work is needed to determine how to better define a coordination shell and compute local interfacial properties. However, since we were able to observe quite a strong correlation of capacitance with ⟨CCpC⟩ with the existing coordination shell definition, we leave this complication for a future study.

Having investigated geometric descriptors and local interfacial properties of EDLCs, averaged over the entire electrode, we find that almost all of them other than ⟨CCpC⟩ lack a clear correlation with capacitance or, in the case of A/⟨d̅sep⟩, are correlated but exhibit significant scatter. In the following sections we turn our attention to the relationship between pore geometry, local electrolyte properties, and charge storage within individual pores of selected materials. We then move toward a more general framework for rationalizing differences in capacitance among ZTC materials. Due to the structural diversity of the ZTC frameworks, we believe insights drawn from ZTCs are also relevant general design rules for porous carbon EDLC electrodes.

Charge Storage Mechanism in Selected Materials

We begin our examination of individual materials by considering BEA and BEA_beta, which are templated on different polymorphs of the same zeolite (known as zeolite beta) as shown in Figure a.[49] Naturally occurring zeolite beta consists of a mixture of polymorphs A and B, both of which contain layers of the same tertiary building unit which are rotated by 90° with respect to each other. In polymorph A (corresponding to BEA_beta ZTC), the layers are stacked in a chiral fashion, while in polymorph B (corresponding to BEA ZTC), the rotation of the layers alternates. As a result, the pore size distributions of BEA and BEA_beta differ, with slightly larger pore sizes for BEA_beta as shown in Figure b.

Figure 4

(a) BEA and BEA_beta ZTC unit cells. (b) pore size distributions. (c) scatterplot showing joint distribution of CCpC and DoC for counterion adsorption sites, with probability distributions on the corresponding axes. (d) conditional expectation of CCpC for a given DoC, denoted ⟨CCpC⟩DoC.

The capacitances of these ZTCs differ widely, with 34.0 F g–1 (45.8 F g–1) gravimetric and 2.33 μF cm–2 (2.65 μF cm–2) areal capacitances computed for BEA (BEA_beta). The ions within the pores also have different degrees of confinement, possibly arising from the slight differences in the most probable pore sizes. As seen in Figure c, the anions in the anode of BEA_beta have a single peak in their DoC histogram around 0.33, while the anions in the BEA have on average higher DoCs, with one peak at 0.35 and another at 0.42. We might suppose from this that BEA should have the higher charge storage efficiency, since Merlet et al. showed that highly confined ions are able to store more charge in supercapacitors;[15] however, in this case the opposite is true: ⟨CCpC⟩DoC is higher in BEA_beta than in BEA for all DoC values (Figure d). In the cathode, as well, the average charge compensation is lower for BEA than for BEA_beta (SI Figure S14).

One noteworthy feature in the charge compensation distribution of the BEA anode is a minimum in ⟨CCpC⟩DoC at 0.43 DoC (Figure d), the location of the higher peak in the DoC histogram. This drop is significant as it does not appear in BEA_beta, showing that sites with the same DoC can have drastically different contributions to capacitance. The source of this low-CCpC region becomes clear when examining representative configurations of the most highly confined anions of BEA and BEA_beta in Figure . In BEA, the anion is located in a cylindrical, nanotube-like structure, with a coordination shell of electrode atoms encircling the anion on all sides, while the anion in BEA_beta is only confined on two out of four sides by the electrode. The cylindrical pore of BEA is too small to fit another anion or even solvent molecule, but too large to snugly fit BF4–, causing it to be stuck in the middle of the pore where it is not close enough to induce a strong compensating charge on any of the atoms within its coordination shell. As a result, the coordination shell atoms in BEA have a total charge of −0.0263 e, while the coordination shell in BEA_beta has a total charge of 0.223 e. This effect has also been observed in the literature; for example, Kondrat et al. showed a local minimum in capacitance when the pore-to-ion-diameter ratio L/d was around 1.5, and increased capacitance when L/d was near 1 and 2.[13]

Figure 5

Snapshots of the coordination shells of highly coordinated anions in the positive electrodes of (a) BEA and (b) BEA_beta. Carbon atoms are colored according to their charge, with blue indicating negative and red indicating positive charge, scaled from −0.01 e to 0.01 e. The green color corresponds to the anions, BF4–, and the orange color to the cations, BMI+; the transparent linear molecules are the solvent.

The charge storage efficiency also depends on the ion size in relation to the pore size: In the electrolyte studied here, the cation is larger than the anion and its charge is distributed on three sites. As such, when a cation is at the center of the nanotube-like pore, the partial charges of its coarse-grained sites are able to approach more closely to the electrode surface, making the equivalent pore in the cathode more efficient at storing charge. This explains why there is no drop in ⟨CCpC⟩DoC in the BEA cathode (SI Figure S14). These observations highlight one important role of pore geometry in determining charge storage efficiency, by influencing local ion density and electrolyte coordination environment.

We next examine another pair of ZTCs, BEC, and ISV, to test our hypothesis that differences in microporosity leads to deviations from the classical capacitance relation, eq . We chose these two structures because they have virtually the same A/⟨d̅sep⟩ (3.5 × 1012 mg–1), but different capacitances of 50 F g–1 and 43 F g–1 for BEC and ISV, respectively, putting them at high and low ends of the capacitance range for the given A/⟨d̅sep⟩ value, as seen in the plot of Csim – Cfit in (Figure ). Analyzing these structures in a similar fashion as for the BEA polymorphs, we find that BEC has more highly confined sites than ISV and a higher ⟨CCpC⟩DoC for all values of DoC (SI Figure S15).

Figure 6

Deviation of gravimetric capacitance from the linear least-squares fit of capacitance vs A/⟨d̅sep⟩. Structures mentioned in the text are labeled. The inset shows the simulated gravimetric capacitance with the fitted line and the color bar indicates the pore diameter.

We visualize in Figure the instantaneous CCpC and d̅sep of each counterion at snapshots taken every 0.5 ps during the production run, along with probability distributions for CCpC and d̅sep along the corresponding axes. The joint distributions for the two materials are largely overlapping, but the adsorption sites in BEC have their highest probability density at slightly lower d̅sep than the sites in ISV, corresponding to more efficient charge storage (as seen by the higher CCpC distribution) and higher overall capacitance of BEC.

Figure 7

Scatterplot of CCpC and d̅sep of each counterion adsorption site during production run, with probability distributions, for BEC and ISV.

By examining individual pairs of structures, and comparing properties of the individual pores inside those structures, we start to obtain insights into how charge storage mechanisms are related to pore geometries. We generalize this approach in the following section, looking at many structures to extract structure–property trends.

Effect of Pore Geometry on Charge Storage Efficiency

We show in Figure the structures of the materials with highest and lowest Csim – Cfit, indicating the average charge of each crystallographically unique electrode carbon atom during the equilibrated constant potential run, along with the probability density isosurfaces of counterion locations within the electrodes. Isosurfaces for adsorption sites with more than 0.1 e total charge compensation are shown in purple, while sites with less than 0.1 e charge compensation are shown in green. This allows us to visually associate geometry with average contribution to capacitance for an individual pore and see that the adsorption site isosurfaces that have more than 0.1 e coordination-shell charge compensation are close to the surface of the frameworks, while the isosurfaces associated with less than 0.1 e charge compensation tend to be in the middle of the pores.

Figure 8

Average electrode atom charges and isosurfaces of ion probability density, computed for anions within the anode (a–e), probability density of average charges in the anode of materials discussed in the text, colored by deviation of gravimetric capacitance from the linear fit (f). Electrode atoms are colored on a scale from 0 (white) to 0.01 e (red). Purple isosurfaces indicate total charge compensation within the coordination shell cutoff radius greater than 0.1 e, while green isosurfaces indicate total charge compensation less than 0.1 e.

Inspecting the average atomic charges of BEC and h91 visually (Figure a and b) as well as quantitatively (Figure f), we find that these two materials have more individual carbon atoms with high charge, corresponding to high CCpC. As seen in Figure , BEC and h91 are also the two materials with the highest enhancement in capacitance compared to materials of similar A/⟨d̅sep⟩. Comparing these two frameworks also illustrates the diversity of factors contributing to overall capacitance: h91 appears to have more atoms with a low average charge than BEC, but from Table one can also see that h91 has a larger fraction of atoms with an average charge of greater than 0.008 e, and from the charge distribution in Figure f, it is apparent that h91 even has a peak in the charge distribution around 0.017 e, whereas BEC has no atoms with such a high average charge.

Table 2

Fraction of Atoms with Large Positive or Negative Average Charge in Selected ZTCs

	percentage (%) of atoms with average charge
material	>0.008 in anode	<−0.008 e in cathode
h91	12	18
BEC	12	6.5
BEA_beta	15	14
ISV	3.5	5.0
h49	3.5	4.0
h18	3.5	4.0
BEA	0.5	5.5

In contrast, h18, h49, and ISV, which have A/⟨d̅sep⟩ similar or greater than that of BEC, but lower capacitance, have fewer highly charged atoms (Figure c–f). Figure f and Table provide further quantitative evidence of the correlation between the probability of highly charged atoms and Csim – Cfit: structures with a higher probability density of average atomic charge magnitude around 0.1 e or higher tend to have a higher Csim – Cfit.

In order to rationalize the differences in charge storage between these materials, we focus on the local radius of curvature of the materials, as this roughly determines the distribution of ion–electrode distances at a particular adsorption site. In BEC, which has square-shaped windows with right-angle “corners”, we see electrode atoms with large average partial charges at two locations for each adsorption isosurface, corresponding to the positions at which an adsorbed ion can be in close proximity with two “walls” of the framework simultaneously. In h91, cylindrical pores adjoining with rounded beams create small radius of curvature sites where ions can again approach the electrode surface closely at multiple sites, leading to more electrode atoms with large partial charges.

In contrast, adsorption sites that are near large radius of curvature sites, such as in h18, h49, and ISV (Figure c, d, and e), tend to be associated with materials with lower capacitances relative to their respective A/⟨d̅sep⟩. In adsorption sites with a large radius of curvature, an ion is not able to induce as many favorable Coulombic interactions with the electrode surface, leading to lower charge compensation for ions within those materials. ISV merits particular mention, as it does contain some adsorption sites with low radius of curvature and high charge compensation (dashed circles in Figure e), but also has high radius of curvature/low charge compensation sites (solid circles); ISV still has a relatively low capacitance considering its high A/⟨d̅sep⟩ (Figure ).

Overall, our results demonstrate that pore geometries that are capacitance-enhancing tend to facilitate the close approach of counterions to multiple carbons within the electrode via low radius of curvature adsorption sites, so that the compensating charge from the electrodes can be localized and large in magnitude to most efficiently screen countercharges and allow for higher counterion loading (and thus a large magnitude of charge storage) in the pores. Conversely, capacitance-diminishing properties include pores with high radius of curvature and cylindrical and ill-fitting pores (whose diameters are not commensurate with the ion diameters), as these types of sites have inefficient charge storage and therefore decrease the overall capacitance of the material.

Conclusions

In summary, we investigated both the charging dynamics and equilibrium capacitance in zeolite-templated carbon electrodes. We found that the equilibration time at constant potential varied depending on the material, from as little as 3 ns to over 10 ns, increasing as the pore size decreased. This indicated that the limiting process in EDLC charge convergence is solvent diffusion and reorganization, which becomes slower as the PLD decreases. We also found evidence of ion trapping and progressive charging of the electrode, which leads to an underestimation of the initial charging profile by a single-exponential fit.

We observed that the equilibrium capacitance of the ZTCs was not correlated with geometric pore descriptors such as PLD, and only weakly correlated with most “globally” averaged local properties such as ⟨DoC⟩. Instead, we found that capacitance was correlated strongly with ⟨CCpC⟩, a measure of how much of the ionic charge is compensated per electrode atom in the coordination shell. Turning our focus next to microscopic properties, we identified subtle geometric differences between otherwise similar ZTC structures that give rise to large differences in capacitance. By comparing probability distributions for charge compensation, degree of confinement, and charge separation distances, we determined specific pore geometries that corresponded to enhanced or diminished charge storage efficiency. For example, we illustrate how a pore with a diameter that is not commensurate with the ion size barely compensates the ionic charge at all, while a slightly smaller pore in a similar material responds to a local ion with a strong compensating charge. Finally, by leveraging local CCpC to characterize ion adsorption sites within the electrode, and analyzing the average charge distribution around low- and high-charge-compensation sites, we observed that the capacitance is higher in materials that have pockets with sharp angles or a low radius of curvature, compared to pockets with rounded corners or a larger radius of curvature. Low radius of curvature pockets allow ions to approach close to the pore at multiple points, inducing strong opposing charges in the electrode which efficiently screen the ion charge.

Our work contributes clear molecular insights into the effect of local pore geometry on the charge storage mechanisms within EDLCs. The efficient charge storage of “pockets” versus “hollow” or “surface” adsorption sites has been previously reported, but our work shows that not all pockets are equivalent; those with diameters that match the ion size, or with sharper corners, are able to compensate substantially more of an adsorbed ion’s charge. The relevance of our results rests on the well-defined and yet diverse pore shapes within the ZTC materials that we study, which allow us to reach conclusions that are not dependent on synthesis conditions or other material-specific factors in amorphous carbons. One promising experimental application of our conclusions is in salt-soft templating,[50] where commensurate ion pore diameters can be engineered by templating a porous carbon with the same salt used in the electrolyte. Our work also signals the potential for improvement in the computational study of EDLC electrodes, in order to develop more sophisticated methods for geometric characterization of pore shapes so that potential electrode materials may be screened with less computational expense, enabling high-throughput screening. Some options for geometric analysis include using Voronoi tessellation or mesh-based approaches to characterize pore curvature, or applying a method such as the Smooth Overlap of Atomic Positions (SOAP),[51] a descriptor of local atomic geometry. Further research into more efficient molecular simulation protocols would also be interesting, either to shorten the equilibration time needed or to determine whether a descriptor such as the CCpC calculated at a constant zero potential is relevant to capacitance at a nonzero potential. These experimental and computational research directions can build upon the insights we present in this work, to advance our fundamental understanding of and design capabilities for improved supercapacitors.

Methods

Model and Force Field

In our simulations we used an organic electrolyte composed of a mixture of 1-butyl-3-methylimidazolium tetrafluoroborate ([BMI+][BF4–]) and acetonitrile (ACN) with the concentration of ions equal to 1 M. We modeled the organic electrolyte using a coarse-grained description consisting of a three-site model for BMI+ and ACN, and a single-site model for BF4–, as shown in SI Figure S1. Nonbonded interactions were described by a pairwise Lennard-Jones potential with Lorentz–Berthelot mixing rules, and electrostatic interactions by a Coulombic potential. Long-range electrostatic interactions were computed using the Ewald summation for systems with slab geometry.[52]

For the nonbonded parameters of BMI+ we used those developed by Roy and Maroncelli.[53] The nonbonded parameters for BF4– and ACN were taken from Merlet et al.[14] and from Edwards et al.,[54] respectively. Bonds and angles of BMI+, and bonds of ACN, were kept rigid using the SHAKE algorithm.[55,56] For the angles of ACN we used a harmonic potential with a stiff spring constant of 400 kcal2 rad–1 mol–1 to keep the molecule close to linear. The carbon atoms of the electrodes were modeled as rigid. During the constant applied potential simulations, the charges of the electrode atoms were computed at each time step according to the constant-potential method.[57,58] All force field parameters and further details regarding the constant-potential method are provided in the Supporting Information.

ZTC materials were synthesized in silico as described in Braun et al.[36] Carbide-derived carbon (CDC) materials, which are studied in depth computationally by Merlet and co-workers,[14,59] are used here as a reference material. CDC structures were taken from Palmer et al.,[60] who generated them using Quench Molecular Dynamics. In this work, we consider 27 ZTC and 2 CDC materials for the constant-charge simulations and a subset of 19 of the ZTCs for constant-potential simulations (CDCs were not simulated at constant potential due to computational expense arising from their large unit cell sizes). CDCs are named as in the original article by Palmer et al. ZTCs are referred to using the name of the templating zeolite. We indicate hypothetical zeolites using the prefix “h” and the last 2 digits of their 7-digit identifier (e.g., h37 for h8326837). Complete names for all the zeolites referenced in the text can be found in the Supporting Information, along with information on framework properties (pore size, void fraction, accessible surface area).

We used a semiautomated protocol to build two-electrode EDLC cells using the Zeo++ software suite[61] and the VMD script interface with the TopoTools package.[62,63] Further details are provided in the Supporting Information. This protocol was designed to fill the EDLC cell with an amount of electrolyte such that when the capacitor is equilibrated at either constant charge or constant voltage, the electrolyte density and composition in the bulk region matches the experimental values. An example of the simulation setup for FAU_1 ZTC is provided in Figure .

Figure 9

Example of simulation setup for FAU_1, with box lengths a = 34.4 Å, b = 30.7 Å, and c = 194.2 Å. BMI+ ions are in red, BF4– ions in blue, and ACN molecules in tan and translucent. The ZTC electrode and graphene caps are in gray. The length of the electrodes in the z direction is 43.4 Å and they are separated by 106.6 Å. Graphene sheets cap the electrodes on both sides of the capacitor cell. The electrolyte consists of 144 [BMI+][BF4–] ion pairs and 1344 ACN molecules. In total, the system contains 9072 atoms.

Molecular Dynamics Simulations

MD simulations were done using the LAMMPS simulation package[64] with Velocity-Verlet time integration using a time step of 1 fs, and a Nosé–Hoover thermostat[65] to maintain a temperature of 300 K. The initial EDLC cell was equilibrated with a constant charge of 0 e on each carbon atom for 4 ns. Final configurations from the zero-charge simulations were used to initialize all further simulation steps.

Molecular simulation is a powerful tool to study EDLCs, as it allows for precise determination of the microscopic properties, such as the structure of the electrolyte within the pores, which can be difficult to access experimentally but which play an important role in determining the capacitance of the material.[13,14] At the same time, simulation of EDLCs presents technical challenges due to long equilibration times[38] and the need to compute the fluctuating charges in the electrode in response to a constant applied potential. The constant potential approach is more accurate but also much more computationally expensive than simulating an EDLC with constant charges on the electrode atoms.[66,67]

We tested multiple protocols for equilibrating the simulation cells, one with a constant-charge equilibration followed by a short constant-potential equilibration step,[14] and the other with a long constant-potential equilibration. In the constant-charge equilibration method, partial charges of ±0.01 e were applied to all the electrode atoms, positive charges for anode atoms, and negative charges for cathode atoms. The EDLC cell was equilibrated with these fixed charges for 8 ns. Then, the effective potential across the cell was calculated using either the 1-D Poisson equation or the averaged local potentials at each electrode atom, and this potential was applied for the constant-potential equilibration and production runs.

In the constant-potential equilibration, a constant potential difference of 1 V was applied across the EDLC cell (±0.5 V for each electrode), and the constant potential simulation was run for at least 10 ns. During the constant-potential run, the average absolute charge on the electrodes was monitored and fit to an exponential. The equilibration step was considered completed when the simulation was at least as long as 5τ, where τ is the time constant of the exponential. This equilibration process was found to be the best following a number of tests which are described in the Supporting Information section, “Development of Computational Protocol”.

Production runs for simulation of capacitance were carried out after equilibration at constant potential. The length of production runs was at least 2 ns. Capacitance was computed using eq , where V = ΔΨ is the voltage drop applied across the cell, and Q is the average absolute value of the charge stored on a single electrode. From the production run we also computed local properties of interest using an in-house software package developed for this study,[68] such as the degree of confinement (DoC)[15] and the charge compensation per carbon (CCpC), in order to understand the mechanisms of charge storage and gain physical insights into differences in capacitances between the materials. The definitions of these local properties, as well as further details regarding the probability density isosurfaces for charge-compensated ion adsorption sites, are provided in the Supporting Information.