Divulging the Hidden Capacity and Sodiation Kinetics of NaxC6Cl4O2

Article pubs.acs.org/JPCC

Divulging the Hidden Capacity and Sodiation Kinetics of NaxC6Cl4O2: A High Voltage Organic Cathode for Sodium Rechargeable Batteries Rafael B. Araujo,*,† Amitava Banerjee,† and Rajeev Ahuja†,‡ †

Condensed Matter Theory Group, Department of Physics and Astronomy, Box 516, Uppsala University, S-75120 Uppsala, Sweden Applied Materials Physics, Department of Materials and Engineering, Royal Institute of Technology (KTH), S-100 44 Stockholm, Sweden

‡

S Supporting Information *

ABSTRACT: In the current emerging sustainable organic battery field, quinones are seen as one of the prime candidates for application in rechargeable battery electrodes. Recently, C6Cl4O2, a modified quinone, has been proposed as a high voltage organic cathode. However, the sodium insertion mechanism behind the cell reaction remained unclear due to the nescience of the right crystal structure. Here, the framework of the density functional theory (DFT) together with an evolutionary algorithm was employed to elucidate the crystal structures of the compounds NaxC6Cl4O2 (x = 0.5, 1.0, 1.5 and 2). Along with the usefulness of PBE functional to reflect the experimental potential, also the importance of the hybrid functional to divulge the hidden theoretical capacity is evaluated. We showed that the experimentally observed lower specific capacity is a result of the great stabilization of the intermediate phase Na1.5C6Cl4O2. The calculated activation barriers for the ionic hops are 0.68, 0.40, and 0.31 eV, respectively, for NaC6Cl4O2, Na1.5C6Cl4O2, and Na2C6Cl4O2. These results indicate that the kinetic process must not be a limiting factor upon Na insertion. Finally, the correct prediction of the specific capacity has confirmed that the theoretical strategy used, employing evolutionary simulations together with the hybrid functional framework, can rightly model the thermodynamic process in organic electrode compounds.

1. INTRODUCTION The quest to identify promising positive electrodes is also extended toward organic molecules along with conventional transition metal oxides.1 This is due to their low production cost, environmental friendliness, and lightweight and relatively simple design strategies that are able to tune the required properties.2−4 Despite the scientific community’s vigorous effort to produce the best organic materials for cathode application in secondary batteries, researchers are still experiencing difficulties to reach the desired properties for these compounds. Those hurdles are low conductivity in the bulk form, low energy densities due to the lower deliverable potentials and, also in many compounds, low stability with limited cycle life.5−7 Quinones derivatives have recently gained attention due to their great deliverable specific capacity,7−11 and high Na insertion potential, comparable to the wellestablished cathode compounds.12 One of the prime reasons behind the excellent performance of quinones is the presence of carbonyl groups, which act as an electroactive center and accoun for the high voltage.10 Recently, Haegyeom et al.10 proposed that chloranil, C6Cl4O2, compound as a cathode of sodium rechargeable batteries. It is a modification of the p-benzoquinone with electron-withdrawing elements; in this case, halogen atoms are fully substituted the H atoms, to achieve high potentials. Due to the high electronegativity of the substituents, the electronic © XXXX American Chemical Society

density must be displaced out of the molecular center, resulting in higher values of the ion insertion potential. They reported an average deintercalation voltage of the order of 2.72 V vs Na/ Na+. This reported potential is comparable to the potential of inorganic compounds used in sodium batteries such as NaFeO2, showing flat voltage profile of 3.3 V vs Na/Na+ related to the Fe3+/Fe4+ oxidation couple13 and reasonably higher than the reported values of organic cathode materials already studied for sodium battery application.5,6,14−17 The excellent performance of the chloranil as a cathode has motivated us to point out the missing links in its electrochemical reaction pathway. Indeed, even with the great achievable potential of the concerned compound, the sodium insertion mechanism behind the cell reaction is still unclear. This is due to the lack of knowledge of the actual crystal structure governed by the material after ionic insertion. Haegyeom et al.10 showed that the host material could accommodate approximately 1.4 Na per formula unit (f.u.) yielding an energy density of 405 Wh/kg. They speculated that the reason behind the lower specific capacity with less than 2 electrons per f.u. is owed to lower electronic conductivity as well as the higher Na−Na interaction upon sodium insertion. Received: April 18, 2017 Revised: June 13, 2017 Published: June 14, 2017 A

DOI: 10.1021/acs.jpcc.7b03621 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

crystal lattice and fractional coordinates were relaxed until forces get smaller than 0.01 eV/Å. Moreover, a plane-wave energy cutoff of 600 eV and a k-mesh with the resolution of 2π × 0.08 Å−1 is used. The selected stable crystal structures coming from the USPEX predictions were subsequently submitted to a second local optimization. In this part, two different strategies were used to model the exchange and correlation term. First, the PBE method including the van der Waals correction within the Grimme’s approach (PBE-D2)26 was employed. It is wellknown that van der Waals interactions play a significant role in molecular systems and pure GGA functionals are not able to properly describe them. Here, we used the same strategy employed in refs 22 and 27 to describe the parameters used in the calculation of the energy dispersion in the Grimme’s approach. Second, the hybrid functional of Heyd−Scuseria− Ernzerhof (HSE06)28,29 was employed to treat the electron over delocalization coming from the self-interaction error of the PBE functional. In this approach, a fraction of the exact Hartree−Fock exchanged is used to compensate the selfinteraction error correctly modeling the electronic localization. In this study, we have varied the exact exchange amount of the HSE06, controlled by the α parameter, as 0.25, 0.35, 0.45, and 0.55. A linear relation of the band gap is obtained with the variation of α, as showed in Figure 1. The same behavior is

The primary goal of this work is to unveil the entity of the sodium insertion mechanism in C6Cl4O2 by predicting the crystal structures of the sodiated phases: Na0.5C6Cl4O2, NaC6Cl4O2, Na1.5C6Cl4O2. and Na2C6Cl4O2. This is done by applying the framework of the density functional theory together with an evolutionary algorithm to find out the most stable phases in each Na content. Then, the crystal structure changes upon each sodium insertion step in the lattice of the C6Cl4O2 are clarified. We have presented a logical explanation for the lower specific capacity obtained in the experiment, which is mostly concerned with the high stabilization of an intermediated phase. The prediction of the Na insertion potential reveals good agreement with the experimental results with a minimum discrepancy of 1% for the first observed voltage plateau and 9.5% for the second plateau when compared with the experimental data. Moreover, the application of the hybrid functional method, properly describing the self-interaction error, changed the most stable intermediated phases in the convex hull. Then, the thermodynamics analysis of the battery reaction showed an unlikely process for the Na insertion after the formation of the compound Na1.5C6Cl4O2. The ionic diffusion was also investigated by employing the nudged elastic band method together with the DFT+D2 formalism. The activation barriers for Na hops are 0.68, 0.40, and 0.31 eV for NaC6Cl4O2, Na1.5C6Cl4O2, and Na2C6Cl4O2, respectively. The analysis of the electronic structure displayed the first reduction reaction occurs mostly on the oxygen sites of the chloranil molecules. Then, for the second step of the reaction, Na varying from 1.0 to 1.5, the extra charge is mostly localized in the rings of the molecule. It elucidates the main reason for the obtained twovoltage plateau in the experiment. Finally, the correct prediction of the specific capacity has confirmed that the used theoretical strategy, employing evolutionary simulations together with the hybrid functional framework, is able to rightly model the thermodynamic process in organic electrode compounds.

2. COMPUTATIONAL DETAILS The framework of the density functional theory (DFT) in conjunction with the evolutionary algorithm USPEX18−20 has been used to search the most stable crystal structures for the compounds NaxC6Cl4O2 (x = 0.5, 1.0, 1.5, and 2). This method proved to be a powerful tool for seeking the complex crystal structures of molecular units.21,22 The main idea of this algorithm is to start the search with a set of randomly created crystal structures and then apply specifically designed operators to modify the initial set of structures and evolving them in the direction to the global minimum. In this case, we have started with an initial population size of 300 structures. The variation operators initially applied are heredity (fracGene), softmutation (fracAomsMut), randomly from the space group (fracRand), molecular orientation (fracRotMut). Moreover, we have considered maximum 30 generations with stopping criteria of 10 times appearance of the same best structure. Each structure built by USPEX underwent an energy minimization process. In this work, the local optimizations were performed on the basis of the density functional theory with the projector-augmented wave (PAW) method as implemented in Vienna ab initio simulation package (VASP).23,24 Spin-polarized generalized gradient approximation was applied in the Perdew, Burke, and Ernzerhof (PBE)25 parametrization to describe the exchange and correlation term of the Kohm−Sham Hamiltonian. The

Figure 1. Computed band gap energy vs the HSE06 mixing parameter α. The dashed line is highlighting the band gap energy predicted by the G0W0−PBE framework. Moreover, the band gap evaluated with PBE and PBE+D2 is showed in the down part of the picture.

shown for different materials in refs 30 and 31. In order to emerge with the proper band gap value of the chloranil structure, we have performed calculation by employing many body perturbation theory in the GW approximation. This approach solves the Hedin’s quasi-particle equations by a first order expansion of the self-energy operator in terms of one body greens function and the screened Coulomb interaction potential. The G0W0 (nonself consistent calculation) was employed on top of the PBE wave functions. This calculation has used 272 electronic bands yielding a band gap energy of 3.98 eV for the C6Cl4O2 structure. Then, a direct comparison of the G0W0 band gap is done to set the optimal value of the mixing parameter, α. Figure 1 shows that the best value of α is 0.55. This mixing parameter produces a band gap with a difference of only 0.08 eV concerning to the one computed in the G0W0 approach. It is also important to highlight that this B


Article


Table 1. Experimental Lattice Parameters from Reference 38 as Well as the Predicted Lattice Parameters Considering the Distinct Schemes Uesd to Treat the Exchange and Correlation Terms of the Kohn−Sham Hamiltoniana a (Å) b (Å) c (Å) α (deg) β (deg) γ (deg) volume (Å3) a

expt

PBE

PBE+D2

HSE06 (25%)

HSE06 (55%)

8.67 5.68 8.54 90.00 105.94 90.00 405.14

9.52 (9.8%) 6.16 (8.4%) 9.15 (7.1%) 90.00 109.00 89.99 507.22 (25%)

8.93 (2.7%) 5.66 (−0.3%) 8.70 (1.7%) 90.00 106.71 89.99 421.79 (4%)

9.32 (7.2%) 5.94 (4.5%) 8.82 (3.2%) 89.99 108.09 90.02 465.02 (15%)

8.79 (1.3%) 5.91 (3.8%) 8.58 (0.4%) 89.99 106.68 90.00 428.19 (6%)

The value in the parentheses is the difference with the experimental value in percentage.

Figure 2. Representative picture of the experimental crystal structure of the chloranil (a) together with the most stable phases found by the evolutionary simulation for the compounds (b) Na0.5C6Cl4O2, (c) NaC6Cl4O2, (d) Na1.5C6Cl4O2, and (e) Na2C6Cl4O2. Here, the yellow, red, green, and brown balls represent, respectively, Na atoms, O atoms, Cl atoms, and C atoms. First two columns of the figure depict different orientational views of all structures. The last column presents the changes in the Na−O bonds upon sodium insertion. The bond distances in this picture were computed in the PBE-D2 approach.

method to fit the mixing parameter using the band gap coming from G0W0 has been used in other investigations.31−34 The average (de)intercalation potential vs Na/Na+ is assessed by computing the Gibbs free energy of the reaction (y−x) (Na+ + e−)+ NaxC6Cl4O2 → NayCl4O2; where the difference of y and x being the charge transfers from the Na

metal anode to the cathode. The Gibbs free energy is accounted as G = E + PV − TS; where E is the internal energy of the system, V is the volume, T is the temperature, P is the pressure, and S is entropy. It is shown that for solid systems, the PV and TS contributions of the Gibbs free energy are very small.35 C


Article


It is also interesting to notice that the prediction of a and c lattice parameters with the hybrid functional reveals better results than the case of PBE-D2, as compared with the experimental measures, while the b parameter is better described by the PBE-D2 scheme. This is related to the molecular orientation and packing directions that can induce greater or lower effects in the treatment of a specific lattice parameter when long-range dispersion corrections are included. The structure parameters of the compounds NaxC6Cl4O2 with x = 0.5, 1.0, 1.5, and 2.0 are summarized in Table A1. Moreover, the structural stabilities in terms of formation energy as a function of x were evaluated for the concerning structures using the following equation:

Therefore, the driving component of the Gibbs free energy is the internal energy of the system. Then, the average (de)intercalation potential can be computed without loss of accuracy as ⎧ E(Na yC6Cl4O2 ) − E(Na xC6Cl4O2 ) − (y − x)E(Na) ⎫ ⎬ V = −⎨ y−x ⎭ ⎩ (1) ⎪

⎪

⎪

⎪

where, the total energies, E, were computed by using the framework of the DFT as described above. The investigation of the minimum energy paths (MEP), through the climbing nudged elastic band (cNEB)36,37 method, was performed by creating five replicas of the system. In which case, the diffusing Na ions were moved in equidistant steps between equilibrium positions. Calculated activation energies converged within 0.01 eV/A. Moreover, an AIMD simulations were performed with a time step of 2 fs producing a sample equilibration of 30 ps for analysis of the most probable pathways of the Na+ hop. A temperature of 800 K was considered in the canonical ensemble (fixed particle number, volume, and temperature, NVT) for all simulations. This temperature was chosen to accelerate the dynamic process.

Ef = E NaxC6Cl4O2 −

(xE Na 2C6Cl4O2 + (2 − x)EC6Cl4O2)

(2) 2 Here E stands for the total energy of the subscripted compound. Figure 3a shows the formation energies evaluated

3. RESULTS AND DISCUSSIONS 3.1. Crystal Structure Evolution upon Sodium Insertion. The crystal structure of chloranil was first proposed by Chu et al.38 at room temperature. More recently, Kubozono et al.39 has also studied this compound by varying the temperature conditions. They confirm a small phase change when the temperature goes smaller than 90 K. In this work, we consider the phase of chloranil investigated at 150 K reported by Kubozono. This material stabilizes in a monoclinic phase with space group P21/a and lattice parameters as summarized in Table 1. The crystal unit shows two crystallographic nonequivalent molecules as shown in Figure 2a. The molecules display an interaction between C and O atoms of distinct molecules forming a network over all the crystal lattice. Due to this interaction the chloranil presents 2D like structure lying on the (001) plane as depicted in Figure 2a1. Table 1 also summarizes the crystal lattice parameters computed by applying different schemes to treat the exchange and correlation term of the Kohn−Sham Hamiltonian. The reproducibility of the chloranil experimental structure gives confidence over the methods used for this study. The structural energy minimization with PBE leads to great overestimations of lattice parameters a, b, and c in comparison with the experimental values. This behavior is consistent with the reported one by Todorova et al.40 It indicates that this functional does not properly describe this material. The inclusion of dispersion, as suggested by Grimme (PBE-D2), leads to much better results than the ones found with pure PBE functional. A very small, of only 4% volume discrepancy compare to the experimental value reflects the necessity to correctly treat the van der Waals dispersion forces. To ensure a proper description of the electronic structure, the HSE06 hybrid functional is employed. The inclusion of short-range exchange partially overcomes the self-interaction error of PBE that causes delocalization of electrons. That creates remarkable improvement in the structural description when α is 0.55 in comparison to the GGA case, as shown in Table 1. Even without properly treat the dispersion interactions, this method itself commits excellent agreement with the experimental data.

Figure 3. Formation energies as a function of x in eV/f.u. with (a) PBE and PBE+D2 and (b) HSE06 with 25% of exact exchange and HSE06 with 55% of exact exchange. The formation energies were computed considering the ground state structures found in the evolutionary simulation.

with the PBE and PBE+D2 schemes. The compound with x = 0.5 appears in the tie-line of the convex hull for the PBE+D2 picture and out of the tie-line for the pure PBE, but with a small energetic distance from the line of the convex hull. On the other hand, the compound with x = 1 emerged as a stable phase for both methods. Finally, the linear combination of the initial and final x compositions (x = 0 and x = 2) showed to be more energetically favorable than the case of x = 1.5 for both the functional. The formation energies depicted in Figure 3b emerge with another scenario. The HSE06 approach, with different amount of exact exchange, has shown the bottom of the hull in a distinct position if compared with the PBE and PBE+D2 methods. Moreover, it displays the case of x = 0.5 out of the convex hull line while the other considered compositions show favorable formations energies. Since the Grimme method (PBE+D2) was successfully able to reproduce the experimental lattice parameters of the C6Cl4O2 with a relatively good agreement with the experimental data, we decided to follow this methodology in the forthcoming discussions of the found structures Na0.5C6Cl4O2, NaC6Cl4O2, Na1.5C6Cl4O2, and Na2C6Cl4O2. D


Article


maximum Na insertion of x = 1.5 against the case of PBE and PBE+D2 with maximum Na concentration of x = 2. The utmost consequence of Na addition into molecular stacking sequence of each crystal structure is also thoroughly investigated. In the case of C6Cl4O2, two different kinds molecular rotation with respect to the a-axis is observed. That suggests two types of stacking: A-type (the blue plane in Figure A1(a)), where molecules are tilted on the left side with respect to the a-axis and in the B-type (the pink plane in Figure A1(a)), molecules are tilted in the right side. Therefore, emanating “AABBAABB...” type of zigzag molecular packing in the c-axis. These two types of plane intersect with a ∼ 89.4° angle. With the addition of Na0.5 per f.u., it is observed a transformation from zigzag to layer types packing. So that, the molecular packing can be either “AA...” or “BB...” type. However, there is an in-plane rotation of the molecules of ∼31° with respect to the c-axis, which is observed in the alternative stacking plane (Figure A1). This type of plane is titled as A′ or B′. Therefore, the molecular packing of Na0.5C6Cl4O2 assumes as “AA′A...” or “BB′B...” packing toward c-axis. The stacking plane distances of either AA′ or BB′ type are on the order of 3.30 Å. The NaC6Cl4O2 compound gets back to “ABAB...” zigzag molecular packing along the b-axis. But, the intrapacking plane distance of A-type and B-type is same as 3.47 Å. The intersection angle of A- and B-type packing plane is much less than 89.4°, so in this case, planes are packed in more compact way. It is also noteworthy to mention that, here, there is no in-plane molecular rotation like in Na0.5C6Cl4O2. In the case of the Na1.5C6Cl4O2 crystal structure, it is obtained a change back to “AA′A...” or “BB′B...” packing sequence along the a-axis. But the stacking plane distance of A−A′ or B−B′ is changed to 3.41 Å. Also the in-plane molecular rotation is lower than the previous one, on the order of 22° with respect to a-axis. The most interesting part is the Na2C6Cl4O2 crystal packing. It does not change back to either “ABAB...” or “AA′A...”. Although it is showing similar layer type of molecular packing along c-axis, as in Na0.5C6Cl4O2, there is no in-plane molecular rotation as previous. Therefore, the packing can be titled as “AAA...” or “BBB...”. So, in general, it is observed that, with the addition of Na, the molecular packing in the crystal structure is governed either by “ABAB...” or “AA′A...” in a flip-flop way except for Na2C6Cl4O2. 3.2. Na Intercalation Voltage. The evaluation of the potential profile is represented in Figure 4. The calculations took into account only the compositions lying on the tie-line of the convex-hull. The results appear in quite good agreement with the experimental measures of the Na insertion potential, as reported in ref.10 First of all, the trend of two voltage plateaus upon the cell discharge is well reproduced by theoretical predictions in the different flavors. For the reaction, where x in Nax is varying from x = 0 to x = 1 the reported experimental potential is 2.9 V vs Na/Na+. The computed values are showing a range from 2.88 to 3.11 V vs Na/Na+ depending on the used exchange and correlation functional. This displays a maximum discrepancy of 7% for the HSE06 with 55% of exact exchange and less than 1% regarding the PBE functional. The second voltage plateau is obtained during the experiments on a potential of 2.6 V vs Na/Na+. In this region of ionic concentration, the predicted Na insertion potential presents a discrepancy of about 9.6% when PBE is employed as the exchange and correlation functional in comparison to the measured data. The addition of the long-range interactions with the PBE+D2 should emerge with better results since it can

The simulations have predicted important changes in the molecular arrangement over sodium insertion. The initial step reaction, where x in Nax vary from 0 to 0.5, has yielded a breaking of C−O linkage between different chloranil molecules of the desodiated phase. Then, the molecules of the desodiated phase have experienced small molecular rotations upon 0.5 Na insertion, producing a distinct layered structure for Na0.5C6Cl4O2, as depicted in Figure 2, parts b and b1. The observed rotation is a product of the Na coordination features. Due to the higher electronegativity of O atoms than the Cl atoms, there exist a tendency of sodium ions to be coordinated by oxygens. Therefore, the observed structural change from x = 0 to x = 0.5 is driven by the quest of Na fractional positions with the highest coordination number possible. Figure 2b2 displays that, in Na0.5C6Cl4O2, all Na atoms are coordinated by four oxygens. It is also observed, from Table A1, a volume change of the order of 1.4% with the addition of Na ions in the crystal lattice when considered the PBE+D2 scheme. The inclusion of more Na ions, varying x from 0.5 to 1, is accompanied by a volume change of about 7%, as observed in Table A1 in the PBE+D2 method. The main observable change in NaC6Cl4O2 compared to the stable phase of Na0.5C6Cl4O2 is the Na−Na distances. As expected, with the insertion of more Na ions, the Na−Na bonds will be lowered, so in this case is showing 3.47 Å compared to 8.57 Å in Na0.5C6Cl4O2. Now, this value is already smaller than the Na−Na bond distance (3.72 Å) in sodium metal.41 Moreover, there are not many discrepancies in the Na−O bond length in both of these compounds, evoking comparable value of 2.29 and 2.33 Å, respectively for Na0.5C6Cl4O2 and NaC6Cl4O2. It is important to highlight that both of these two stable phases reveal Na coordinated by four oxygens. The next step of the reaction, with x going from 1.0 to 1.5, reveals a break in the trend of Na coordination. As depicted in Figure 2d2, two of the nonequivalent Na atoms, at fractional coordinate positions Na1 and Na2, are surrounded by only three O atoms while the Na ions at position Na3 are coordinated by four O atoms. Na1 and Na2 display a bonding length distance of 3.04 Å what leads to a great electrostatic repulsion between them due to this very small bond distance. Even though, it is still energetically favorable to maintain this structure with these small Na distances. The final step of the reaction takes place with x varying from 1.5 to 2.0. Figure 2e depicts the global minimum structure found in the evolutionary simulation for Na2C6Cl4O2. It is expressing a volume change of about 7.4% as compared to the compound with 1 Na per f.u. in the PBE+D2 scheme. This volume variation comes due to two main reasons. First, the distance between the chloranil molecules centeres increases from 3.47 to 3.82 Å when x varies from 1.0 up to 2.0. The higher distance is a consequence of the extra electron being localized in the chloranil molecule, increasing the electrostatic repulsion among them. The second comes due to the Na coordination number. To make possible all the Na atoms being surrounded by three oxygen atoms, as shown by the most stable phase of Na2C6Cl4O2 in Figure 2e2, it becomes necessary to have some volume expansion. Otherwise, the Na−Na bonds distance would be very small with the insertion of the extra Na per f.u. and the electrostatic repulsion among them will be considerably high to stabilize the compound. Even with this volume change, the bonds among the Na−Na atoms go to 3.20 and 3.29 Å. We have discussed more on top of this issued in the next section where we show that the HSE06 function predicts a E


Article


Haegyom et al.10 have speculated that the lower specific capacity is attributed to the high Na−Na electrostatic repulsion in the last stage of the discharge together with the lower conductivity of the chloranil compound. However, it seems that the responsibility for the lower capacity it is not only related to the Na−Na repulsion, but there also exist an influence of the reduced chloranil molecules repelling each other and the great stabilization of the Na1.5C6Cl4O2 compound. 3.3. Electronic Structure and Charge Transfer upon Sodium Insertion. Figure 5 displays the projected density of states (pDOS) of NaxC6Cl4O2 (x = 0, 1.0, and 1.5) with the HSE06 approach considering 55% of exact exchange. The pDOS of the C6Cl4O2 compound confirms the insulator to semiconductor behavior of the material exhibiting a band gap of 4.0 eV. In the Fermi level, negligible contributions of the Figure 4. Calculated voltage profile of NaxC6Cl4O2 vs Na/Na+ for x varying from 0 to 2. The dashed line represents the experimental values of 2.9 and 2.6 V.9 The red line is the calculations with PBE, the blue line represents the calculations with PBE+D2, and the green line is the calculations with HSE06 with 25% of exact exchange and orange is the HSE06 with 55%.

properly describe the related structural properties. However, no real improvement is obtained for the computation of the ion insertion potential with this method in comparison with the pure PBE. A similar result is also reported by Tomoki Yamashita et al.21 The evaluation of the intercalation potential employing the HSE06 method revealed features that the pure PBE and the PBE+D2 were not able to show. In fact, the discrepancies between the experimental potentials and the ones computed with the HSE method yield greater values than the case with the pure PBE. On the other hand, the inclusion of a percentage of exact exchange in the calculation induces the Na insertion to a maximum concentration of x = 1.5 (since the insertion potential goes to a negative value for x going to 1.5 to 2). This result is in really good agreement with the experimental measurement, which prescribes a maximum reversible Na concentration of x = 1.42 during the cell reaction.10 In order to have insights of the reason why the HSE06 method ends up with a maximum concentration of x = 1.5 while the PBE predicts a reaction occurring up to x = 2, the total DOS of the Na1.5C6Cl4O2 compound in the PBE+D2 and HSE06 for α = 0.55, have been analyzed and plotted in Figure A2. The addition of the exact exchange partially cancels the selfinteraction error inherent from the PBE method. This is reflected, for instance, in the metal behavior obtained for the PBE+D2 method against the semiconductor properties showed by the HSE06 scheme increasing the stability of the referent intermediated phase, as showed in the convex hull picture. Finally, this effect leads to a negative voltage value, for x varying from 1.5 to 2 in Na2C6Cl4O2, with HSE06 while the PBE calculations predict a positive Na intercalation potential. It is also interesting to notice that the chloranil molecules distances, in the compound Na2C6Cl4O2, present a value of 3.82 Å with the PBE+D2 approach while the HSE06, at 55% of exact exchange, shows a value of 3.93 Å. It is speculated that the greater electron localization in the molecules, reflected by the HSE06 method, induces to higher electrostatic repulsion among them, finally lowering its stability in comparison to the case with PBE+D2.

Figure 5. Projected density of states (pDOS) for C6Cl4O2, NaC6Cl4O2 and Na1.5C6Cl4O2, respectively. Here, gray background stands by the total DOS, red represents the summed states of O atoms, green is the added Cl contributions and brawn represents the summed C states. F


Article


The charge density difference between the compounds Na 1.5 C 6 Cl 4 O 2 and C 6 Cl 4 O 2 ([ρ(Na 1.5 C 6 Cl 4 O 2 ) − ρ(C6Cl4O2)]) has followed the same procedure explained above for the case of one Na insertion per f.u., and it is shown in Figure 6b. The insertion of the extra 0.5 Na per f.u. evokes the distinct characteristics when compared with the case where x varies from 0 to 1.0. In this structure, the extra charge coming from the 0.5 Na ends up in only one of the two adjacent C6Cl4O2 molecules of the unit cell structure. The fingerprint of this localization appears in the plotted charge density difference as well as in the C−O bonds of the carbonyl groups. It is clearly noticeable in Figure 6b that one of the chloranil molecules of the Na1.5C6Cl4O2 unit cell presents greater electronic density in the 6C ring in comparison to the other. Moreover, the C−O bond in extra-reduced molecule goes to 1.30 Å while, for the other case, the C−O bonds remains with 1.25 Å as in the NaC6Cl4O2 compounds. This feature is also reflected in the following fingerprints: the sharp peak inside of the band gap of the DOS revealed by NaC6Cl4O2 gets more delocalized in Na1.5C6Cl4O2 since the extra 0.5 electron mostly reduces the ring of half of the chloranil molecules in the crystal. The second is the obtained two voltage plateaus in the experiment. It means that in the first stage of the cell reaction, O states are being reduced producing a potential of the order of 2.9 V vs Na/Na+. The second part of the reaction, the reduction starts to occur mostly on the ring of the molecule, up to the point where almost half of the chloranil are double reduced and the other half is single reduced. Thus, this effect gives rises to a slightly small voltage of 2.6 V vs Na/Na+. 3.4. Sodium Insertions Kinetic. To further understand the intercalation of Na ions onto the C6Cl4O2 crystal framework, the cNEB method was employed together with the DFT+D2 formalism to get insights of the ionic migration mechanisms. Three compounds were considered for the investigation of the Na migration, NaC6Cl4O2, Na1.5C6Cl4O2 and Na2C6Cl4O2 where the predicted crystal structures of each case were employed. The most likely pathways for Na insertion are determined based on the IAMD simulations. For instance, Figure A3 depicts the Na trajectories during the simulation time scale for NaC6Cl4O2. The ionic diffusion mechanism revealed by this picture indicates that Na ions will hop between equilibrium positions moving through the formed layers of Na−O atoms producing a 1 D like diffusion process. This diffusion mechanism is also observed for the compounds Na1.5C6Cl4O2 and Na2C6Cl4O2. It is worthy to mention that there is no case where Na ions migrate from one Na−O layer to other. Hence, the migration pathways chosen for the computation of the activation barriers will consider only the ionic diffusion through the Na−O layers. Figure 7 displays the most likely paths for the Na+ diffusion in the materials (a) NaC6Cl4O2, (b) Na1.5C6Cl4O2, and (c) Na2C6Cl4O2 together with each of the energy profiles coming from the cNEB calculations. In the first case, NaC6Cl4O2, single diffusion pathway is investigated (path 01 shown in brown) resulting in an activation barrier of 0.68 eV. Three different paths are studied in the compound Na 1.5C 6Cl 4O2 , as represented by red, green and brown colors in Figure 7b. The barrier for the Na hop through path 03, red color, goes to 0 eV. However, the ionic migration mechanism in this compound cannot occur only by hops though the Path 03. Then, the diffusion mechanism must come with a composition of pathways. Path 01 and path 02, represented by brown and green, respectively, emerge with activation barriers amounting

oxygen states are obtained with a greater contribution of C and Cl states. On the other hand, the bottom of the conduction band is mostly formed by oxygen and carbon states. This analysis indicates that upon reduction reaction, the electron donated by the extra Na, must partially occupy oxygen states as well as carbon states with minor contributions of Cl atoms. The pDOS of NaC6Cl4O2 is depicted in Figure 5b. The greater contributions to the DOS, in the Fermi level, come from O and C states, as expected. It also shows the appearance of a new localized state inside the band gap, at around 1.7 eV. This sharp peak present features of a deep impurity. Thus, this localized state would not induce to high conductivity in NaC6Cl4O2 since the sharp state at 1.7 eV must produce a high-localized band with great electron effective mass. Therefore, even after the insertion of one Na per f.u. of the material would still suffer from low electronic conductivity. Finally, Figure 5c presents the pDOS for the Na1.5C6Cl4O2 compound. The extra inserted 0.5 Na per f.u. does not lead to great changes of the pDOS, in comparison to the NaC6Cl4O2. It is still observed the extra state inside the band gap, but, for this material, the sharp state becomes more dispersive. It means that the electrons donated by the extra addition of 0.5 Na would be less localized than the case of NaC6Cl4O2. Insights of the electronic transfer process during the cell reaction can be obtained by the subtraction of the charge densities of the sodiated and desodiated phases. Figure 6a

Figure 6. Isosurface of the charge density difference with sodium insertion computed by the HSE06 method with a mixing parameter of 0.55: (a) the difference [ρ(NaC6Cl4O2) − ρ(C6Cl4O2)]; (b) the difference [ρ(Na1.5C6Cl4O2) − ρ(C6Cl4O2)]. The yellow and blue colors represent the positive and negative isosurface at +0.008 and −0.008 eV/Å3.

shows the charge density difference between NaC6Cl4O2 and the respective desodiated compound, C6Cl4O4, computed in the HSE06 approximation with 55% of exact exchange. In this case, we have performed a new electronic energy minimization of the crystal framework of the sodiated material, NaC6Cl4O4, by removing the Na ions. It permitted us to calculate the appropriate charge density difference. This procedure leads to good acumen of the localization of the electron donated by the Na atom upon the insertion.42 It is very clear from the electronic density difference that the first sodium insertion is accompanied by a reduction over all O atoms of the C6Cl4O2 molecules considered in the unit cell confirming the insights derived from the pDOS analysis. Moreover, the C−O bonds in the carbonyl groups undergo an expansion from 1.20 to1.25 Å, due to the antibonding character of the newly occupied orbital. G


Article


Figure 7. Activation barriers and corresponding diffusion pathway of Na+ of (a) NaC6Cl4O2, (b) Na1.5C6Cl4O2, and (c) Na2C6Cl4O2.

knowledge of these crystal structures allowed us to predict the sodium intercalation voltages for x going from 0 to 2 with different functionals. The voltage computed with PBE shows excellent agreement with the first voltage plateau, only 1% of discrepancy with the experimental data, whereas showing 9.5% discrepancy in the second plateau. The addition of part of the exact exchange, through the hybrid functional approach, showed the importance of correctly treat the electronic localization to properly predict the maximum allowable Na insertion in the chloranil cathode. The analysis of electronic structure showed that the first part of the reduction reaction, with one sodium insertion per f.u., is mostly followed by the population of O states while in the second part of the reduction, where x varies from 1.0 to 1.5, the extra electron mostly reduces the 6C ring of the chloranil molecules. This explains the reason without any speculation of the two plateaus obtained in the experiment. The calculations showed that there exist two main effects as the responsible to the lower specific capacity in comparison to the theoretical one. The first is the great stability presented by the Na1.5C6Cl4O2 compound. The second is the existence of a higher electrostatic repulsion among the reduced chloranil molecules and among the Na ions in the Na2 C6 Cl 4O2 compound. It finally produces an unlikely reaction from the Na1.5C6Cl4O2 to the Na2C6Cl4O2. The computation of the Na+ energy barriers for the Na+ hop has revealed large to small activations energies with the increase of the Na content. It indicates that the ionic diffusion is not a limiting process in the electrode reaction. Finally, the correct prediction of the specific capacity has confirmed that the used theoretical strategy, employing evolutionary simulations together with the hybrid

to 0.72 and 0.40 eV. It is most likely, then, that the Na ions migrate between the Path 02 and Path 03 resulting in the global activation energy of 0.40 eV. Finally, we also have investigated the diffusion process in the compound Na2C6Cl4O2, even considering that it will not be observed in the experiment. This compound comes with two likely path-accumulating barriers of 0.18 and 0.31 eV. Such as for Na1.5C6Cl4O2, the diffusion mechanism has to come as a composition of both paths and the hindering process appears with the barrier of 0.31 eV. One remarkable observation is that the higher Na concentration leads to lower activation energies. The general barriers amounts to 0.68 eV, 0.40 and 0.31 eV for NaC6Cl4O2, Na1.5C6Cl4O2, and Na2C6Cl4O2, respectively. This means that the increasing Na content in the host compound, indeed, leads to higher diffusion coefficients. Therefore, the observed lower specific capacity, as compared with the theoretical specific capacity, obtained by the experiments must have a small influence on the ionic diffusion process.

4. CONCLUSIONS First-principles calculations have been successfully employed together with an evolutionary algorithm to study key properties of chloranil as the cathode for the next generation of organic rechargeable batteries. We have predicted the ground state phases for the compounds NaxC6Cl4O2 with x = 0, 0.5, 1.0, 1.5, and 2.0. The evolutionary simulations have predicted great structural changes under the sodium insertion reaction with a volume change going to about 14% when x varies from 0 to 2 in the PBE+D2 approach. The analysis of the most stable structures reveals a tendency of Na ions being coordinated by 4 to 3 oxygen atoms. Moreover, this feature seems to drive the main changes occurring during the Na insertion process. The H


Article


(12) Meng, Y. S.; Arroyo-de Dompablo, M. E. First principles computational materials design for energy storage materials in lithium ion batteries. Energy Environ. Sci. 2009, 2, 589−609. (13) Okada, S.; Takahashi, Y.; Kiyabu, T.; Doi, T.; Yamaki, J.-i.; Nishida, T. Layered transition metal oxides as cathodes for sodium secondary battery. Meet. Abstr. 2006, MA2006-02, 201−201. (14) Wang, S.; Wang, L.; Zhang, K.; Zhu, Z.; Tao, Z.; Chen, J. Organic Li4C8H2O6 Nanosheets for Lithium-Ion Batteries. Nano Lett. 2013, 13, 4404−4409. (15) Wang, S.; Wang, L.; Zhu, Z.; Hu, Z.; Zhao, Q.; Chen, J. All Organic Sodium-Ion Batteries with Na4C8H2O6. Angew. Chem. 2014, 126, 6002−6006. (16) Luo, W.; Allen, M.; Raju, V.; Ji, X. An Organic Pigment as a High-Performance Cathode for Sodium-Ion Batteries. Adv. Energy Mater. 2014, 4, 1400554. (17) Lee, M.; Hong, J.; Seo, D.-H.; Nam, D. H.; Nam, K. T.; Kang, K.; Park, C. B. Redox Cofactor from Biological Energy Transduction as Molecularly Tunable Energy-Storage Compound. Angew. Chem., Int. Ed. 2013, 52, 8322−8328. (18) Oganov, A. R.; Glass, C. W. Crystal structure prediction using ab initio evolutionary techniques: Principles and applications. J. Chem. Phys. 2006, 124, 244704. (19) Oganov, A. R.; Lyakhov, A. O.; Valle, M. How evolutionary crystal structure prediction works and why. Acc. Chem. Res. 2011, 44, 227−237. (20) Lyakhov, A. O.; Oganov, A. R.; Stokes, A. R.; Zhu, Q. New developments in evolutionary structure prediction algorithm USPEX. Comput. Phys. Commun. 2013, 184, 1172−1182. (21) Yamashita, T.; Momida, H.; Oguchi, T. Crystal structure predictions of NaxC6O6 for sodium-ion batteries: First-principles calculations with an evolutionary algorithm. Electrochim. Acta 2016, 195, 1−8. (22) Zhu, Q.; Oganov, A. R.; Glass, C. W.; Stokes, H. T. Constrained evolutionary algorithm for structure prediction of molecular crystals: methodology and applications. Acta Crystallogr., Sect. B: Struct. Sci. 2012, 68, 215−226. (23) Kresse, G.; Hafner, J. Ab initio molecular dynamics for liquid metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 47, 558−561. (24) Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169−11186. (25) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (26) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 1787−1799. (27) Frayret, C.; Izgorodina, E. I.; MacFarlane, D. R.; Villesuzanne, A.; Barrès, A.-L.; Politano, O.; Rebeix, D.; Poizot, P. Eletrochemical properties of crystallized dilithium squarate: insight from dispersioncorrected density functional theory. Phys. Chem. Chem. Phys. 2012, 14, 11398−11412. (28) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. J. Chem. Phys. 2003, 118, 8207− 8215. (29) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Erratum: “Hybrid functionals based on a screened Coulomb potential” [J. Chem. Phys. 2003, 118, 8207 ]. J. Chem. Phys. 2006, 124, 219906. (30) Seo, D.-H.; Urban, A.; Ceder, G. Calibrating transition-metal energy levels and oxygen bands in first-principles calculations: Accurate prediction of redox potentials and charge transfer in lithium transition-metal oxides. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 92, 115118. (31) Araujo, R. B.; Chakraborty, S.; Ahuja, R. Unveiling the charge migration mechanism in Na2O2: implications for sodium−air batteries. Phys. Chem. Chem. Phys. 2015, 17, 8203−8209. (32) Choi, M.; Janotti, A.; Van de Walle, C. G. Native point defects and dangling bonds in α-Al2O3. J. Appl. Phys. 2013, 113, 044501.

functional framework, is able to rightly model the thermodynamic process in organic electrode compounds.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.7b03621. Lattice parameters with different functionals and also stacking of each structure (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*(R.B.A.) E-mail: [email protected]; rafael.barros@physics. uu.se. ORCID

Rafael B. Araujo: 0000-0002-3964-2807 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS A.B. and R.B.A. are thankful to Erasmus Mundus for the doctoral fellowship. We would like to acknowledge the Swedish Research Council (VR), the Swedish Energy Agency and StandUP for financial support. SNIC, HPC2N, and UPPMAX are acknowledged for providing computing time.

■

REFERENCES

(1) Vadehra, G. S.; Maloney, R. P.; Garcia-Garibay, M. A.; Dunn, B. Naphthalene Diimide Based Materials with Adjustable Redox Potentials: Evaluation for Organic Lithium-Ion Batteries. Chem. Mater. 2014, 26, 7151−7157. (2) Liang, Y.; Tao, Z.; Chen, J. Organic electrode materials for rechargeable lithium batteries. Adv. Energy Mater. 2012, 2, 742−769. (3) Burkhardt, S. E.; Lowe, M. A.; Conte, S.; Zhou, W.; Qian, H.; Rodríguez-Calero, G. G.; Gao, J.; Hennig, R. G.; Abruña, H. D. Tailored redox functionality of small organics for pseudocapacitive electrodes. Energy Environ. Sci. 2012, 5, 7176−7187. (4) Häupler, B.; Wild, A.; Schubert, U. S. Carbonyls: powerful organic materials for secondary batteries. Adv. Energy Mater. 2015, 5, 1402034. (5) Armand, A.; Grugeon, S.; Vezin, H.; Laruelle, S.; Ribière, P.; Poizot, P.; Tarascon, J.-M. Conjugated Dicarboxylate Anodes for Li Ion Batteries. Nat. Mater. 2009, 8, 120−125. (6) Chen, H.; Armand, M.; Courty, M.; Jiang, M.; Grey, C. P.; Dolhem, F.; Tarascon, J.-M.; Poizot, P. Lithium Salt of Tetrahydroxybenzoquinone: Toward the Development of a Sustainable Li-Ion Battery. J. Am. Chem. Soc. 2009, 131, 8984−8988. (7) Yao, M.; Yamazaki, S.; Senoh, H.; Sakai, T.; Kiyobayashi, T. Crystalline Polycyclic Quinone Derivatives as Organic PositiveElectrode Materials for Use in Rechargeable Lithium Batteries. Mater. Sci. Eng., B 2012, 177, 483−487. (8) Chihara, K.; Chujo, N.; Kitajou, A.; Okada, S. Cathode properties of Na2C6O6 for sodium-ion batteries. Electrochim. Electrochim. Acta 2013, 110, 240−246. (9) Zhu, Z.; Li, H.; Liang, J.; Tao, Z.; Chen, J. The disodium salt of 2,5-dihydroxy-1,4-benzoquinone as anode material for rechargeable sodium ion batteries. Chem. Commun. 2015, 51, 1446−1448. (10) Kim, H.; Kwon, J. E.; Lee, B.; Hong, J.; Lee, M.; Park, S. Y.; Kang, K. High energy organic cathode for sodium rechargeable batteries. Chem. Mater. 2015, 27, 7258−7264. (11) Li, A.; Feng, Z.; Sun, Y.; Shang, L.; Xu, L. Porous organic polymer/RGO composite as high performance cathode for half and full sodium ion batteries. J. Power Sources 2017, 343, 424−430. I


Article

The Journal of Physical Chemistry C (33) Alkauskas, A.; Broqvist, P.; Pasquarello, A. Defect levels through hybrid density functionals: Insights and applications. Phys. Status Solidi B 2011, 248, 775−789. (34) Radin, M. D.; Siegel, D. J. Charge transport in lithium peroxide: relevance for rechargeable metal-air batteries. Energy Environ. Sci. 2013, 6, 2370−2379. (35) Aydinol, M. K.; Kohan, A. F.; Ceder, G.; Cho, K.; Joannopoulos, J. Ab Initio Study of Lithium Intercalation in Metal Oxides and Metal Dichalcogenides. Phys. Rev. B: Condens. Matter Mater. Phys. 1997, 56, 1354−1365. (36) Henkelman, G.; Jónsson, H. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. J. Chem. Phys. 2000, 113, 9978. (37) Jónsson, H.; Mills, G.; Jacobsen, K. W.; Nudged elastic band method for finding minimum energy paths of transitions, in Classical and Quantum Dynamics in Condensed Phase Simulations; World Scientific: 1998; pp 385−404. (38) Chu, S. S. C.; Jeffrey, G. A.; Sakurai, T. The crystal structure of tetrachloro-p-benzoquinone (chloranil). Acta Crystallogr. 1962, 15, 661−671. (39) Kubozono, Y.; Yoshida, T.; Maeda, H.; Kashino, S.; Terauchi, H.; Ishii, T. X-ray structure analyses of chloranil above and below phase transition temperature. J. Phys. Chem. Solids 1997, 58, 1375− 1381. (40) Todorova, T.; Delley, B. Molecular Crystals: A Test System for Weak Bonding. J. Phys. Chem. C 2010, 114, 20523−20530. (41) Aruja, E.; Perlitz, H. Neubestimmung der Gitterkonstante von Natrium. Z. Kristallogr. - Cryst. Mater. 1939, 100, 195. (42) Chevrier, V. L.; Ong, S. P.; Armiento, R.; Chan, M. K. Y.; Ceder, G. Hybrid density functional calculations of redox potentials and formation energies of transition metal compounds. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 075122.

J


Divulging the Hidden Capacity and Sodiation Kinetics of NaxC6Cl4O2

Recommend Documents