Exploring Metastable Phases During Lithiation of Organic Battery Electrode Materials.

In this work, the Li-ion insertion mechanism in organic electrode materials is investigated through the lens of atomic-scale models based on first-principles theory. Starting with a structural analysis, the interplay of density functional theory with evolutionary and potential-mapping algorithms is used to resolve the crystal structure of the different (de)lithiated phases. These methods elucidate different lithiation reaction pathways and help to explore the formation of metastable phases and predict one- or multi-electron reactions, which are still poorly understood for organic intercalation electrodes. The cathode material dilithium 2,5-oxyterephthalate (operating at 2.6 V vs. Li/Li+) is investigated in depth as a case study, owing to its rich redox chemistry. When compared with recent experimental results, it is demonstrated that metastable phases with peculiar ring-ring molecular interactions are more likely to be controlling the redox reactions thermodynamics and consequently the battery voltage.

Introduction

Organic‐based electrodes have gained attention in the last few years as a promising pathway to achieve a truly sustainable battery technology. Connected to this, the market price volatility of common materials used in inorganic‐based batteries together with the environmental harms they may impose have raised several concerns about the future of electrical energy storage (EES) devices. Likewise, the advent of novel ultra‐portable gadgets brought the demand for higher‐performing EESs to a new paradigm. Recently, a large set of organic electroactive materials (OEM) have been explored as viable options for novel Li‐ion batteries (LIBs)[ , ] to truly achieve sustainable battery technologies. Among these, there exist several examples of compounds with high energy densities,[ , , , , ] fast reaction kinetics,[ , ] high capacities, good cyclability[ , ] or sustainable synthesis routes,[ , , , ] but to date no candidate has been reported that possesses all of these desired properties. A deeper understanding of the underlying physicochemical properties of these types of materials would be a key component to successfully realize their final technological application. Knowledge about structural changes occurring in the electrode during battery operation – during the Li‐ion insertion/deinsertion reactions – would be of foremost importance to capture the relevant structure‐property relationships that control the battery's performance.

In particular, the molecular and solid‐state structures composing the OEM are directly associated with several key properties that affect the lithiation reaction thermodynamics, such as the Li‐ion coordination environment, electronic structure, or reaction kinetics. Furthermore, the Li‐ion insertion process is known to occur by different mechanisms leading to distinct phenomena of the electrode materials such as phase separation and/or the appearance of metastable phases. In the former case, no stable intermediate phases appear during lithiation, thus resulting in a non‐stepwise process. This phenomenon is known to occur for several inorganic LIB cathodes, such as Li FePO4[ , ] and Li Ti2(PO4)3[ , , ] and has recently been pointed out for a few OEMs.[ , , ] Similarly, the lithiation reaction occurring through metastable phases may also cause undesired effects, such as the appearance of complex structural ordering inside the electrode and corresponding changes on the potential energy landscape for the ion diffusion process, which in turn affects the reaction kinetics. Metastable phases controlling the lithiation mechanism have been investigated for some inorganic electrodes,[ , , ] with special attention given to the Li FePO4 compound.[ , , , ] For organic counterparts, recent studies[ , , ] have been reporting similar phenomena, but these are still poorly understood and require further analysis.

In this context, an evolutionary algorithm interplayed with density functional theory (DFT) has previously been successfully employed to predict the crystal structure of several OEMs during their lithiation processes.[ , , , ] However, this algorithm searches for the most energetically favorable phases – the global energy minima in the configurational space – and is therefore not suitable to analyze metastable structures. To address this issue, we have here developed a new computational methodology to predict metastable phases appearing during the lithium insertion process. This novel approach relies on a successive mapping of the potential energy surface (PES) of the solid‐state structure to perform ‘educated guesses’ on the position of the Li‐ions during lithiation. We present an insightful comparison between these two different pathways (two‐phase and metastable phase) to evaluate the structural changes that dicarboxylate‐based molecular crystals undergo upon battery operation, using dilithium 2,5‐oxyterephthalate (Li2DHT; Figure 1a) as a case study. The new method successfully identifies metastable lithiated phases for this specific OEM, which seems to offer a more reasonable comparison with experimental data. The battery open‐circuit voltage, for instance, shows a discrepancy of only 2.7 % regarding experimental results, which is an improvement over the evolutionary algorithm method which leads to a discrepancy of 4.2 %. In a nutshell, the novel method differs to other traditional approaches to identify metastable phases in two central aspects: (i) due to its conceptual simplicity, the method is less computationally demanding, hence usually offering faster assessments; (ii) it is mainly focused on ion‐insertion processes.

Figure 1

a) Lewis representation of the Li2DHT molecule. b) Average voltage of the two‐step reaction (vs. Li/Li+) for all the investigated prediction schemes and theory levels. The dashed line represents the experimentally reported value.[ , ] c) Electronic energy of each lithiated phase relative to pristine Li2DHT. d) Formation energy of the intermediate Li3DHT with respect to Li2DHT and Li4DHT phases.

Results and Discussion

The crystal structures of the Li DHT phases were initially obtained following two different methodologies: (i) an evolutionary algorithm (EA) approach, as employed in some of our previous works,[ , , ] and (ii) a novel PES mapping framework that is performing educated guesses on the Li‐ion positions for lithiated phases (hereafter named MAP). Although the GGA‐PBE functional was used to obtain the initial geometries for all phases, the hybrid functional HSE06 was employed further using two different pathways: (i) a single‐point calculation without any further geometry optimizations (HSE06‐O), and (ii) a new geometry optimization considering this second functional (HSE06‐O). These methodologies are further described in the Computational Methods section and all the structures are provided in the Supporting Information. The use of hybrid functionals aims to improve the overestimated electron delocalization present in the GGA‐PBE scheme. Thus, through these two paths we can analyze if the hybrid functional should be used to further optimize the crystal structure or, alternatively, if merely correcting its electronic structure would suffice.

The Li2DHT molecule is represented in Figure 1a while Figure 1b shows the average voltage for the two‐step lithiation reaction of Li2DHT. While the initial assessment provided by the EA‐PBE method of 2.91 V exhibits a tolerable discrepancy to the experimental value of roughly 2.6 V,[ , ] the MAP‐PBE presents a much closer result of 2.59 V. For both EA and MAP cases, employing the hybrid functional without structural reoptimization (HSE06‐O) results in a voltage increase, thereby worsening the comparison with the experimental outcome. However, reoptimizing the geometries with the hybrid functional (HSE06‐O) improves the agreement with the experimental outcome considerably, with MAP showing a lithiation voltage of 2.53 V and EA one of 2.71 V. Therefore, the use of hybrid functionals is also relevant to obtain molecular geometries of the respective molecular building blocks. This may be due to the fact that the exchange interaction plays an important role in the molecular bonding[ , , ] and the addition of the exact exchange energy present in the hybrid functional could offer an improvement when studying molecules. To further analyze this outcome, Figure 1c presents the relative energy as a function of lithiation step, that is, the electronic energy of each phase ( ) in comparison to the respective delithiated ( ) phase . The HSE06‐O offers structures with higher energies for both EA and MAP cases, with an upward shift in energy when compared to the HSE06‐O counterparts. Furthermore, the good agreement with the experimental voltage provided by this method suggests that the lithiation reaction may occur through metastable phases. It is also important to note that the EA and MAP approaches differ fundamentally on how the structures are obtained, with the former tending to achieve more stable phases. This issue is shown by comparing energy levels of both schemes in Figure 1c and could be the reason for the better agreement of MAP voltages with the experimental value of around 2.6 V. In fact, the higher energies offered by the MAP approach further reinforce the conclusions about metastable phases as the lithiation mechanism. When analyzing the formation energy of Li3DHT, Figure 1d shows that the EA‐PBE/HSE06‐O supports the appearance of this intermediate phase, whereas for HSE06‐O, this phase would not be observed. Moreover, MAP points to this latter result regardless of the formalism. Nevertheless, this thermodynamical analysis suggests the occurrence of a two‐electron process.

The Li‐ion uptake occurring via a two‐step reaction and metastable phases may lead to intrinsic changes in the material's structure. Figure 2a shows the crystal structures of the pristine Li2DHT (EA), the first reduced Li3DHT (EA/MAP) and Li4DHT (EA/MAP) as obtained from the HSE06‐O approach. These are the structures assumed henceforth in this analysis and their lattice parameters are given in Table 1, whereas all obtained structures are provided in the Supporting Information. Additionally, an ab initio molecular dynamics (AIMD) simulation was employed to investigate the stability of the predicted Li2DHT phase. The AIMD indicates that the compound is stable at room temperature during the simulation time length, maintaining the initial crystal structure characteristics (see the Supporting Information for details). At first glance, the EA seems to provide a well‐defined layer of Li‐ions coordinated by oxygen atoms and intercalated by the organic counterpart, for both Li2DHT and Li4DHT phases, whereas MAP shows that the inserted Li‐ions (Li4DHT) are found both in the salt layer and intercalated between the organic parts of the structure. The contrast regarding the positioning of the inserted Li‐ions appears to originate from the lack of space within the delithiated material's salt layer to accommodate the two extra lithium atoms. Therefore, MAP places some ions amidst the organic layers while the EA solves this issue by finding a structure with a thicker salt layer containing all the Li‐ions, which may justify the difference in the lattice parameter between the two methodologies (Table 1).

Figure 2

a) Two‐dimensional (planar) representations of the complex 3D crystal structures (EA and MAP) from (de)lithiated phases of Li2DHT. b) Li−O bond length distribution considering a cutoff radius of 2.8 Å around the respective Li atom for different prediction schemes.

Table 1

Lattice parameters for pristine Li2DHT and reduced Li4DHT for the EA and MAP prediction schemes.

	a [Å]	b [Å]	c [Å]	α [°]	β[ °]	γ [°]
Li2DHT EA	9.15	9.12	4.33	92.59	96.76	65.69
Li4DHT EA	4.25	6.05	12.97	94.76	75.03	102.61
Li4DHT MAP	8.83	9.46	4.03	87.68	86.59	70.05

[°]

°]

[°]

Figure 2b presents the Li−O bond length (BL) distribution for a cutoff radius of 2.8 Å around the Li atom. The BL analysis of these atom pairs further improves our understanding of these materials as oxygen plays a key role in coordinating Li‐ions, hence affecting the thermodynamics of the reaction. Furthermore, Figure 2b also offers a direct comparison between the EA and MAP structure prediction schemes for the Li4DHT compound. The Li‐ion coordination (inferred qualitatively in conjunction with the structure analysis) in this phase changes based on the chosen scheme: the EA shows coordination numbers of five and four while the MAP display six, four and two. The peak near 1.8 Å that appears for both methods is due to different reasons. For EA, this BL is related to the quinone's oxygen and Li‐ions with a coordination number of four, whereas for MAP this is due to the lithium found between molecules and being coordinated by two different COO oxygens (i. e., with a coordination number of two). The EA peaks around 2.0 and 2.4 Å are due to COO oxygens while quinone oxygens are present at 1.9 and 2.6 Å. MAP shows a broader peak around 2.1 Å related to COO and quinone oxygens, with a small peak near 1.9 Å due to carboxylates.

Following experimental findings from Chen and colleagues, the lithiation of Li2DHT to Li4DHT appears to occur through stepwise one‐electron processes, in contrast with what our thermodynamical analysis suggests. Nonetheless, the cyclic voltammogram (CV) Chen and colleagues have presented shows a single oxidation peak at 2.7 V for the first cycle (one single two‐electron step) and the appearance of two peaks near 2.5 V and 2.6 V for the next cycle (two one‐electron steps). The fact that the Li4DHT was initially obtained through a chemical synthesis suggests that its crystal structure is more likely to be represented by the EA method, shown in Figure 2a. As discussed above, the evolutionary algorithm renders more stable phases and this could be on par with the synthesis process. During the first cycle, the oxidation of Li4DHT to Li2DHT follows a two‐electron process, and the calculated EA voltage of 2.71 V shown in Figure 1a agrees well with the reported CV single peak of 2.7 V. In the next cycle, however, the CV peaks of 2.5 V/2.6 V have a better agreement with the MAP voltage of 2.53 V. Thereafter, it is highly likely that the pristine (as‐synthesized) Li4DHT phase is not fully recovered after the first cycle and the lithiated phases may therefore follow the metastable structures presented by the MAP scheme. Furthermore, the appearance of two CV peaks could be related to the two different sites for inserted Li‐ions, that is, within the lithium salt layer and amidst the organic molecules (Figure 2a).

The charge density distribution of the two added electrons after reduction from Li2DHT to Li4DHT is illustrated in Figure 3 for both EA and MAP schemes when following the HSE06‐O and HSE06‐O approaches. This is helpful when analyzing the redox active center and where in the structure the charge tends to localize. The charge resulting from MAP appears to be more equally distributed over molecular rings and quinone oxygens (Figure 3a, b), in contrast with EA where the charge is unevenly distributed over these units. However, when optimizing geometries with the hybrid functional (HSE06‐O), the evolutionary algorithm presents a charge density uniformly localized over rings and quinoid oxygens (Figure 3c). Incidentally, Li‐ions are also more evenly placed within the salt layer. MAP, in contrast, shows a slightly delocalized charge density over the organic part (Figure 3d), which could be the reason for the decreased lithiation voltage as compared to EA. Surprisingly, a strong interaction between the molecular rings seems to appear (Figure 3e) with a shared charge density between two different rings. The fact that the two displayed carbons are 1.73 Å apart from each other seems to further characterizes this interaction, although not supporting that a bond is actually forming between the rings. This could be a direct consequence of the interstitial Li+ present near the molecules, which may induce an unbalanced electronic charge distribution over the molecular ring.

Figure 3

a–d) Charge densities of the two inserted electrons after Li2DHT reduction to Li4DHT as obtained by following a) EA‐HSE06‐NO, b) MAP‐HSE06‐NO, c) EA‐HSE06‐O and d) MAP‐HSE06 O. e) A closer look at the ring‐ring interaction resulting from MAP‐HSE06‐O. The isosurface values are 0.008 in all cases.

Figure 4 shows the fragment projected density of states (pDOS) normalized by the number of atoms in each fragment. The pDOS for structures obtained using the EA methodology are shown in Figure 4a, b, and d for Li2DHT, Li3DHT and Li4DHT, respectively. For the MAP framework, Figure 4c shows the pDOS for Li3DHT and that for Li4DHT is shown in Figure 4e. The first unoccupied band (Figure 4a) is mainly dominated by the quinone oxygen and the molecular ring (C6H2), with a minor contribution from the carboxylate unit (COO). This result, in conjunction with Figure 3, strengthens the proposed reaction mechanism for this compound, in which the double bonding quinone oxygen is expected to be reduced upon lithiation. As a result, charge is transferred to the molecular rings, which gains aromaticity in the process. For the Li3DHT phase (Figure 4b, c), both EA and MAP methods show that the first band is being partially filled by the electron inserted upon lithiation. Nevertheless, our thermodynamics analysis suggests that this first lithiated phase would not be stabilized. The Li4DHT phase (Figure 4d, e) presents a sensible contrast in the electronic structure between both methods. The two inserted electrons after the reduction of Li2DHT into Li4DHT are found from −2.5 to 0 eV for EA and from −1.75 to 0 eV for the MAP, which corroborates the fact that the former provides a more stable structure and better stabilizes the injected electron upon lithiation. Additionally, MAP shows more perceptive changes on the quinone oxygen states compared to the pristine Li2DHT and a higher contribution from the carboxylate states in accommodating the upcoming electrons. This COO contribution is mainly due to the oxygens coordinating the interstitial Li ion (Figure 3d).

Figure 4

The fragment projected density of states (pDOS) normalized by the fragment's number of atoms: a) Li2DHT; b) Li3DHT EA; c) Li3DHT MAP; d) Li4DHT EA; e) Li4DHT MAP. The vertical dotted line in these plots represents the Fermi level.

Conclusions

In this work, two distinct computational methods to analyze the lithiation mechanism on organic‐based electrodes for Li‐ion batteries are compared. The first of these methods (EA), based on an evolutionary algorithm, stands out in finding the most stable phase of the compound; the structure with the global energy minimum. However, this method is less suitable when the lithiation reaction happens through metastable phases. In this sense, a second method (MAP) is developed, which is based on a sequential mapping of the structure's potential energy surface to provide educated guesses of the most suitable place to insert a Li+.

When both methods are employed to study the Li2DHT compound, the analysis suggests that the lithiation in this system indeed occurs through metastable phases, that is, the thermodynamics presented by the MAP method correspond better to the reported experimental values. In addition, the relevance of using the correct functional within the theoretical framework to achieve an improved description of the system is also demonstrated. Moreover, when highlighting the main differences between the metastable phase from MAP and the more stable phase from EA, it can be concluded that the MAP method has proven to be a suitable choice to identify metastable phases during the Li+ insertion in organic‐based battery electrodes.

Computational Methods

The electrochemical lithiation reaction can be written as Equation 1:

where is the host material and lithium metal is chosen as the reference electrode. To obtain the voltage profile for this reaction, the Nernst equation is employed [Eq. (2)]. Henceforth, all presented voltage values are referenced to the lithium metal electrode (V vs. Li+/Li). Furthermore, the total electronic energy is expected to be the dominant term in the reaction free energy, hence it is used as an approximation for the Gibbs free energy.

The total electronic energy depends directly on the crystal structure of the considered (de)lithiated phases. We here chose the first two lithiation steps for the host material (the first two Li‐ion insertion reactions) and therefore the and phases are needed. Additionally, the formation energy of the intermediate phase after the first lithiation reaction was evaluated through Equation 3:

where represents the same total electronic energy as in Equation (2) and the formation energy is given with respect to and . This method for evaluating the formation energy is a convenient tool when analyzing alloys and is useful to estimate the stability of certain such compositions. For the lithiated phases in organic electrodes, this type of analysis provides valuable insights regarding if the electrochemical process is happening through a one‐step or a two‐step reaction (i. e., if the intermediate phase stabilizes or not).

Evolutionary algorithm (EA)

The evolutionary algorithm has shown to be a powerful tool to predict the crystal structure of a given material by following a process where an initial population of proposed candidates goes through several generations of solid‐state structures, and then finally “evolving” into the most energetically suited candidate. For OEM crystal structures, the first generation candidates can be randomly generated from a list of possible space group symmetries, having the molecular units and the desired number of Li‐ions as building blocks. For instance, the Li2DHT molecule here served as the building block for the delithiated phase, whereas the lithiated phases were obtained by combining this building block with the proper amount of Li atoms (Li3DHT and Li4DHT). This strategy of using the entire molecule as a building block significantly reduces the effort in searching the combinatorial space necessary to find the desired crystal structure. As for the evolution process, a number of genetic operations called heredity, mutation and permutation are employed to systematically create new generations of candidates, forwarding the best characteristics from one generation to the next and thereby ensuring the evolution scheme. For the compounds studied here, the evolution started with a group of 200 randomly generated candidates, followed by a set of 20 candidates per generation for a maximum of 40 generations.

In practice, this process has been performed through an interplay of the evolutionary algorithm as implemented in the USPEX code[ , ] with DFT calculations performed in the Vienna Ab‐initio Simulation Package (VASP).[ , , ] The criterium for the evolution process was the total electronic energy obtained from DFT, in which calculations have been carried out by employing the projector augmented wave (PAW) method and the GGA‐PBE as exchange‐correlation functional. For each generated structure candidate, a two‐step geometry optimization was performed without any constraints. A k‐mesh resolution of 0.12  and plane‐wave cut‐off energy of 400 eV were used for the first optimization step while a higher k‐mesh resolution of 0.08  and a plane‐wave cut‐off of 500 eV was used for the second step. After the evolution scheme, a new geometry optimization was performed for the best candidate with a 600 eV plane‐wave cut‐off and a ‐mesh.

Potential mapping (MAP)

This alternative method, which this study explores, offers the possibility to guess the inserted ion's position inside the electrode during the electrochemical reaction. In order to do so, the first aspect that must be considered is the charge state of the ion (a positively charged state in the case of Li‐ion). Thereafter, the potential energy surface (PES) inside the crystal structure, generated by electrons and atomic nuclei forming the electrode, must be considered. The PES can be evaluated by employing ab‐initio methodology and solving the corresponding quantum mechanics equations, for example, through DFT. In the first approximation considered, the Li+ feels this potential energy surface when being inserted, and a logical guess for the new cation's position would be at the PES maximum. However, this approach can only provide an initial educated guess and the interactions between the inserted ion and other components of the system must be considered further. In fact, the presence of a new Li+ is going to systematically affect the PES itself. This can be remedied by carrying out a complete DFT geometry optimization after the insertion, relaxing the structure to a local energy minimum. The updated PES may then be used to find the position of a new Li+ being inserted, and then repeating the process as many times as necessary to represent the desired lithiation steps. The entire process can be summarized by the workflow presented in Figure 5. One drawback of this method is that it relies on the crystal structure of the delithiated phase to sequentially anticipate the structure of lithiated phases. Therefore, it is necessary to rely on either experimental results or use an initial evolutionary algorithm prediction for the delithiated crystal structure. The novel method derived is here also introduced as a new software package, named Mapion, coded in the python language and is freely distributed on an open‐source basis at GitLab. Although the software is optimized to work with VASP, it can be easily tuned to work with other platforms. In addition, it is also possible to investigate different ion‐insertion elements (e. g., Na+, K+).

Figure 5

Workflow illustrating how the potential mapping methodology identifies the lithiated phases. PES=potential energy surface.

In short, the process starts by calculating the PES of the delithiated phase, carried out in VASP by employing the PAW method and the GGA‐PBE as exchange‐correlation functional with a ‐mesh of and a plane‐wave cut‐off energy of 600 eV. A new Li‐ion is then placed at the PES maximum, generating a new structure. This new lithiated phase is further fully relaxed to a new energy minimum, following a two‐step calculation: first using a ‐mesh and 500 eV energy cut‐off, and then a ‐mesh with 600 eV as energy cut‐off. From this final geometry, a new PES can be calculated, repeating the process to obtain the next lithiation step.

Post‐processing step

After finding the desired lithiated phases for both EA and MAP methods, the hybrid functional HSE06 was employed for two different pathways: (i) a single‐point calculation without any further geometry optimizations with a ‐mesh (HSE06‐O) and (ii) a new structure optimization considering only the gamma‐point in the reciprocal space, followed by a second calculation with the new relaxed geometry and a k‐mesh of (HSE06‐O). Finally, the lithiation voltages and formation energies were evaluated following Equations (2) and (3).

For all performed calculations, using both EA and MAP methods, the Grimme‐D2 semi‐empirical corrections for dispersion interactions were considered, except for Li atoms. The increasing amount of lithium in the lithiated phases may result in error accumulation due to limitations of the DFT‐D theory. In addition, the convergence criterion for all geometry relaxations was the norms of all forces being smaller than 0.01 eV  .

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supporting Information

Click here for additional data file.