Towards novel calcium battery electrolytes by efficient computational screening

Abstract The development of Ca conducting electrolytes is key to enable functional rechargeable Ca batteries. The here presented screening strategy is initially based on a combined density functional theory (DFT) and conductor-like screening model for real solvents (COSMO-RS) approach, which allows for a rational selection of electrolyte solvent based on a set of physico-chemical and electrochemical properties: solvation power, electrochemical stability window, viscosity, and flash and boiling points. Starting from 81 solvents, N,N-dimethylformamide (DMF) was chosen as solvent for further studies of cation-solvent interactions and subsequent comparisons vs. cation-anion interactions possibly present in electrolytes, based on a limited set of Ca-salts. A Ca2+ first solvation shell of [Ca(DMF)8]2+ was found to be energetically preferred, even as compared to ion-pairs and aggregates, especially for PF6− and TFSI as the anions. Overall, this points to Ca(TFSI)2 and Ca(PF6)2 dissolved in DMF to be a promising base electrolyte for Ca batteries from a physico-chemical point-of-view. While electrochemical assessments certainly are needed to verify this promise, the screening strategy presented is efficient and a useful stepping-stone to reduce the overall R&D effort.


Introduction
Batteries based on multivalent chemistry have emerged as promising alternatives and complements to the currently dominant lithium-ion battery (LIB) technology, mainly due to foreseen high volumetric capacities and improved safety [1][2][3][4] . Amongst these, calcium (Ca) based batteries show some unique advantages [5] : i) Ca is the 5 th most abundant element in Earth's crust [6] -which makes the chemistry long-term sustainable, ii) Ca tends to grow uniformly during battery cycling -i.e. no dendrites form under "normal " operating conditions, as opposed to both Li and Mg metal batteries [2] , and iii) the negative standard potential of Ca is -2.87 V [ 2 , 3 , 6 ], and thus enables higher voltage cells than other multivalent chemistries, especially as compared to Mg and Al [7] . Foremost, however, very high volumetric energy densities can be achieved by employing Ca metal anodes with a volumetric capacity of 2,073 mAh cm − 3 [2] . Yet, the lack of Ca conducting electrolytes, enabling efficient Ca plating/stripping vs . Ca metal anodes and high ion conductivities at room temperature, has stalled Ca batteries from being demonstrated at large scale and furthermore also efficient Ca cathodes from being developed.
For any modern multivalent rechargeable battery technology the optimization of the electrolyte composition and studies of the resulting performance, must be based on considering several basic physico-chemical and electrochemical properties such as: electrochemical stability window (ESW), ionic conductivity, salt solubility, viscosity, flash and boiling points, etc [ 1 , 8 , 9 ]. These properties are largely determined by the choices made for salt and solvent(s) and the salt concentration, and ultimately also affect the viability of efficient Ca metal plating/stripping. Ponrouch et al. [7] showed that it is feasible to plate/strip Ca at temperatures > 75°C by employing organic electrolytes conceptually similar to those of the LIB technology and Wang et al. [10] achieved reversible electrodeposition even at room temperature. Very recently, a few electrolyte concepts have been demonstrated which indeed showed salt solubility as the main issue towards an efficient Ca based electrolyte. Standard solvents such as ethylene carbonate (EC) and N,Ndimethylformamide (DMF) both displayed reasonable salt solubility and dissociation [11] . There are, however, still challenges to overcome; temperature dependence [7] , cycling efficiency, oxidation stability [10] , etc. , which calls for effective strategies to develop new Ca conducting electrolytes.
Addressing both the wide range of electrolyte properties needed and the large search space of salts and solvents, novel electrolyte development is here suggested to be rationalized by applying a computational screening strategy [12][13][14][15][16][17][18][19]   screening of both modern battery and double layer capacitor electrolytes [ 12 , 14 , 17 ]. Properties such as melting/flash/boiling points, solubility, viscosity and ESW were used as parameters to find new electrolytes for electrochemical double layer capacitors in the work of Schütter et al., where finally cyano esters were identified as promising solvents [14] , and Husch et al. reported a volunteer computing approach targeting the identification of electrolytes for Li-ion batteries [12] . Comparisons to experimental data confirmed the efficiency of this strategy to rank electrolytes. The DFT + COSMO-RS methodology has previously been used to investigate Li-S battery electrolytes, both to predict elemental sulfur solubility [20] and to explain the advantages and the change in transport mechanism resulting from applying fluorinated solvents [21] . While the COSMO-RS approach in general has difficulties in properly describing systems with highly localized charge densities, such as alkali and alkali earth metal cations, e.g. Ca 2 + , this was successfully resolved for Li + in those studies. We are therefore confident to employ the DFT + COSMO-RS approach in the search for novel Ca battery electrolytes.

Screening protocol
The screening is divided into two main steps: i) Selection of a solvent, based on evaluation of properties such as solubility, ESW, viscosity, flash point and boiling point, and ii) Selection of a Ca-salt, by elucidating the Ca 2 + coordination number (CN) in the selected solvent and evaluating the cation-anion interactions ( Fig. 1 ). For the former step, 81 solvents (Table S1) were chosen based on know-how from the fields of Li-ion, Naion and Mg rechargeable batteries, while the latter applied four Ca-salts readily available commercially and with traditionally battery electrolyte used weakly coordinating anions (WCAs): BF 4 − , BH 4 − , PF 6 − , and TFSI. The first, solvent selection, step started by computing the relative solubility of the four salts in each of the 81 solvents and from this the solvents with reasonably large relative solubilities were selected for further studies. These were assessed for their ESWs by computed reduction and oxidation potentials, to pass the criteria of the reduction potential being lower than the cathodic limit, 0 V vs . Ca 2 + /Ca °, and the oxidation potential higher than the anodic limit, here a bit arbitrarily set to 4.5 V vs . Ca 2 + /Ca °. These potential limits, including computational errors and uncertainties, should allow for medium to high voltage Ca battery cells to be created using the selected solvents. Finally, viscosities, flash and boiling points were evaluated for the solvents having passed the two previous stages. A low viscosity is essential to ensure high power rate capability in the cell as the (vehicular) ion transport is viscosity dependent [22] and for proper wetting of porous electrodes [23] , while the flash and boiling points are safety related properties -both at the cell manufacturing stage and in the battery usage phase. Finally, for the salt selection, various complexes using the different anions (An) of the type [Ca(An) n ] 2-n (n = 1-4) were investigated and the cation-anion interaction strengths used to assess the probability of ion-pair and aggregate formation for the best solvent from the above screening, which provides information on the speciation and thereby indirectly also on the (cat)ion transport.

Computational details 2.2.1. Solubilities
The mole fractional absolute salt solubility (x_S) is in COSMOtherm treated as: where X and S are the chemical potentials of the pure and solvated salt, respectively, ΔG fus is the free energy of fusion of the salt [24] . With ΔG fus being unknown for all the four Ca-salts used, no prediction of absolute solubilities can be performed, but setting ΔG fus = 0 [25] creates a constant shift of the absolute solubilities and the relative solubilities (S_R) can simply be evaluated as ln(S_R) = ( X -S )/RT. The criterion for a reasonable relative solubility is defined as: 0 < Log 10 (S_R) < -5. The chemical potentials in Eq. (1) are, hence, obtained by integrating the chemical potential of a solvent´s surface segment with charge over the screened surface of the solute as: ( ) is the segment chemical potential with screened charge and x is the compound mole fraction in the mixture. Size and shape differences are accounted by adding an extra combinatorial term, , that depends on volume and area of the mixture compounds with some more adjustable parameters [26] .
The DFT calculations were performed as implemented in the TUR-BOMOLE package [27] . The BP86 functional was used together with a TZVP basis set for all geometry optimizations [28][29][30] . Single-point calculations on the converged geometries employed a TZVPD basis set, to generate the -surface for the subsequent COSMO-RS step, with a fine grid cavity. All subsequent COSMO-RS calculations, as implemented in the COSMOtherm package, used a BP_TZVPD_FINE_C30_1701 parametrization at 298.15 K [ 26 , 31-33 ].

Electrochemical stability windows
The thermodynamic ESW for a species is determined by the redox potentials: where n is the number of electrons participating in the reduction/oxidation reaction and F is the Faraday constant. Δ ( ) and Δ ( ) are the Gibbs free energy changes associated with the reduction and oxidation reactions, respectively, and here computed as: The changes in the Gibbs free energies are computed as: where ΔE is the electronic energy change, ΔE ZPE is the zero-point energy change, ΔH and ΔS are the changes in the different enthalpy and entropy contributions, respectively, and ΔG solv accounts for solvation effects. ΔE was here computed by DFT employing the M06-2X [34] functional and a TZVP basis set as implemented in TURBOMOLE. The redox potentials were calculated using a vertical, non-adiabatic, transition between ground and excited state, employing the Franck-Condon approximation. Thermal contributions were evaluated by vibrational frequency calculations. Finally, ΔG solv was evaluated within the COSMO-RS framework (see 2.2.1) and the resulting redox potentials were shifted by 1.53 V to obtain a scale vs . Ca 2 + /Ca °.

Viscosity, boiling and flash points
Viscosity, boiling and flash points were all computed using COS-MOtherm. The viscosity ( ) was computed based on the Quantitative-Structure-Property-Relationship (QSPR) model [35] as: where is the surface area of the solvent i from its -surface, is the second -moment of the solvent, is the number of ring atoms in the solvent and the pure solvent entropy times temperature. The five constants (C x ) were derived from a set of 175 room temperature viscosities of organic compounds. Boiling points were computed by raising the temperature of the solvent until the difference between the COS-MOtherm prediction of the system total vapor pressure and a preselected 'pressure' (1 atm) was < 10 − 4 mbar [ 12 , 26 ]. Similarly, the flash points were calculated by raising the system temperature, computing the vapor pressure, and a difference vs. the flash point pressure < 10 − 4 mbar [ 12 , 26 ]. The flash point pressure is calculated using the solvent surface areas as the only descriptor, in general a rather satisfying strategy [36] .

Ca 2 + solvation structure and interactions
To gain insight on the interactions possibly present in the electrolytes, which are a function of the stability of the first solvation shell of Ca 2 + , both generic [ Ca ( solvent ) x ] 2+ cation-solvent complexes (x = 5-9) and [Ca(An) n ] 2-n (n = 1-4) cation-anion complexes were optimized using a minima hopping (MH) global optimization, as implemented in the Atomic Simulation Environment (ASE), using temperature to overcome energy barriers and explore the energy landscape [ 37 , 38 ]. Short molecular dynamics (MD) simulations were used to escape local minima, followed by local optimizations using on-the-fly adjusted parameters. The ab initio MD (AIMD) parts of the MH were performed in the microcanonical ensemble (NVE) with a velocity Verlet algorithm and a time step of 1 fs. Local optimizations were performed in the framework of the BFGS quasi-Newton algorithm. For both, the BP86 functional was employed together with a SVP basis set, as implemented in TURBOMOLE. The MH algorithm started at 300 K and an initial energy threshold of 48 kJ mol − 1 . The energetically most stable structures for each complex obtained by the MH algorithm were further geometry optimized at the M06-2X/TZVP level of theory for the binding energy, ΔE b , calculations.
The binding energies, Δ , of the [ Ca ( solvent ) x ] 2+ complexes (x = 5-9) were calculated as: where [ Ca ( solvent ) x ] 2+ and are the electronic energies corrected by the zero point energies (ZPE) and solvation energies. 2+ is the cation electronic energy and is the number of solvents in the complex. The solvation energies were calculated within the COSMO-RS framework following the procedure outlined in 2.2.2 and the ZPEs were obtained from vibrational frequency calculations. The preferred Ca 2 + cation solvation number (SN) is, hence, the x providing the lowest Δ .
Similarly, the dissociation energies, Δ , of the complexes [Ca(An) n ] 2-n were calculated as: where n is the number of anions and , A n − and C a 2+ are the electronic energies corrected by the ZPE and solvation energies of the complexes and components, respectively. The Δ indicates the cationanion interaction strength and, therefore, the "free " ions vs. ion-pairs and higher aggregates balance for a specific choice of anion.

Results and discussion
First the solvents were screened for Ca-salt solvation power and subsequently the ESWs were computed for the more promising solvents, whereafter their viscosities, boiling and flash points were all analyzed. Subsequently, the details of the Ca 2 + coordination and ion-pair and aggregate formation were assessed for the most promising solvent.

Screening
First, the relative solubilities for the four Ca-salts in each of the 81 solvents (Table S1) were calculated using COSMO-RS and to further assist the interpretation of the results they were visually divided into five families: boron, fluorine, nitrogen, sulfur containing, and carbonates/ethers (colour coded in Fig. 2 a).
The relative salt solubilites were normalized vs. DMSO and the lowest were obtained for the fluorine containing solvents, which is quite natural as fluorine is an electron-withdrawing group that decreases the solvent donor number (DN) and hence lowers the cation solvation ability [39] . Donor interaction solvents act as Lewis bases towards cations and the DN is a measure of the Lewis basicity (the ability to donate electrons). By using, the limited set of, DN data available in the literature, the higher solubilities obtained for the N and S containing solvents can vice versa be explained by their higher DNs ( Fig. 2 b). Hence, to increase the solubility of Ca-salts electron-donating groups to further increase the DNs is instrumental. All the solvents show the same trends and no significant changes with respect to the different Ca-salts used, which is not surprising given the much higher solvation energy of the Ca 2 + cation as compared to that of the anions (Table S2).
Second, the ESWs of these solvents were calculated. All four solvents have reduction potentials below the cathodic limit, i.e. 0 V vs. Ca 2 + /Ca °. They also have oxidation potentials above 4.5 V (the by us set limit). Thus they are all, in principle, applicable to high voltage Ca battery cells ( Fig. 3 ).
Overall, the nitrogen-based solvents (DMA and DMF) have wider ESWs than the two sulfur-based solvents (DMSO and TMSO), and DMF has the highest both reductive and oxidative stabilities. Kumar et al. pointed out that interfacial interactions between solvents and electrodes may alter the solvent HOMO and LUMO energy levels, hence, changing the ESW [40] . However, DMF has a stability of 2.2 V below the cathodic limit and 1.4 V above the anodic limit, and these large margins should minimize the risk of solvent decomposition on the electrode surface, even if no interfacial interactions were explicitly taken into account.
Next, the viscosities, boiling and flash points of the selected solvents were assessed and these calculated physico-chemical properties all qualitatively agree with the experimental data ( Table 1 ).
Notably, a property like viscosity is notoriuosly difficult to predict as it is a dynamic property which depends on multiple phenomena, which largely explains why the quantitative deviations vs. experimental data are the largest. Nevertheless, the two nitrogen-based solvents display the lower viscosities, in agreement with the experimental data, due to weaker intermolecular interactions [ 49 , 50 ], and DMF the very lowest viscosity due to its lower molecular weight. On the other hand, the safety assessments, i.e. the boiling and flash points, render the sulfur-based solvents better choices, whereas the flash points of DMF and DMA might be an issue vs. the operation/storage temperature of the battery and/or in the cell production processes.
Overall, the screening results in DMF to be the preferred solvent based on the combination of a high solvation power, a wide ESW, and a relatively low viscosity, but again attention must be given to its low flash point.

Ca 2 + solvation structure and ion-pair and aggregate formation
With the aim to in more detail understand the nature of the charge carriers in DMF, the binding energies, ΔE b , between Ca 2 + and DMF in the complexes [ Ca ( DMF ) x ] 2+ (x = 5-9) were calculated. The more stable species are those with x = 7 and x = 8, i.e. [ Ca ( DMF ) 7 ] 2+ and [ Ca ( DMF ) 8 ] 2+ ( Fig. 4 ), with a slight preference for the latter. However, the ΔE b difference of 5.9 kJ mol − 1 between these species is insignifi-    [51] . Assuming no ion-pairing, which as a first order approximation is reasonable as the perchlorate anion is a WCA, and mono-dentate solvent coordination, the CN is equal to the SN for Ca 2 + . Furthermore, Asada also showed Mg 2 + to have a CN/SN = 6 in DMF, which supports our general argument that Ca 2 + has a higher CN/SN, much due to its larger radius providing the space to fit more DMF ligands, without causing too much steric hindrance. The latter is indeed the reason to as to why SN = 9 becomes less energetically beneficial.
The possible formation of ion-pairs and higher aggregates directly influences the details of the Ca 2 + charge carriers, which is addressed by modeling a wide range of [Ca(An) n ] 2-n complexes for n = 1-4 ( Fig. 5 , Table 2 ).
From these simple models, all complexes with n = 1 and n = 2 have the anions coordinating Ca 2 + tri-dentately, excepted TFSI for n = 2. This agrees well with previous DFT studies that showed TFSI to prefer to coordinate bi-dentately in various [M(TFSI) 2 ] complexes (M = Li, Mg, Ca, Ba, Zn and Cu) [ 52 , 53 ]. For n = 3 there is a wide distribution of bi- dentate and tri-dentate coordination and correspondingly also a wide variety in CNs; 6, 9, 7 and 6, for the anions BF 4 − , BH 4 − , PF 6 − and TFSI, respectively. Moving to n = 4 the complex [ Ca ( B F 4 ) 4 ] 2− has three of the anions in bi-dentate coordination and one in mono-dentate coordination, hence CN = 7. In the literature, MD simulations have indeed shown the BF 4 − coordination to Ca 2 + to strongly depend on the SN [54] . Yet, the "isocompositional " [ Ca ( B H 4 ) 4 ] 2− has all its anions coordinating tri-dentately (CN = 12), and [ Ca ( P F 6 ) 4 ] 2− in contrast have all its anions coordinating bi-dentately (CN = 8). Altogether, this suggests the BH 4 − anions to coordinate Ca 2 + stronger, which could be detrimental to the (cat)ion tranport and transfer, due to the high dissociation energies of the [Ca(BH 4 ) n ] 2-n complexes ( Table 2 ). Indeed, Hahn et al. attributed the very low ionic conductivity of Ca(BH 4 ) 2 in THF ( ca. 4 S.cm − 1 ) to the very strong cation-anion association, even at 0.1 M salt concentration [55] . Finally, the [ Ca ( TFSI ) 4 ] 2− resulted in a first Ca 2 + solvation shell with SN = 3 and CN = 6, i.e. similar to n = 3, most likely due to steric hindrance effects, and is not treated further here.
In general, the n = 3 complexes have larger total cation-anion interaction strengths, as per their dissociation energies, but some of them only marginally so as compared to n = 2 ( Table 2 ). The exception is the BH 4 − anion where n = 4 is energetically preferred, in line with the strikingly different CN of this complex. Monti et al. showed n = 3 to be the energetically preferable for Na + as the central atom and TFSI as anion [56] , and this is corroborated here, even considering the multivalent nature of the Ca 2 + cation, which makes the complex [ Ca ( TFSI ) 3 ] − "only " single-negatively charged. Overall the TFSI and PF 6 − anions display the lower dissociation energies, which also could be expected from previous studies using Li + and Na + as cations [ 57 , 58 ].
The ion-pair and aggregate formation is a function of the Ca 2 + first solvation shell stability with solvents and/or anions. Hence there exists a competition which we here simply assess as |ΔE b | − |ΔE d |. For the TFSI anion, the highest Δ is 1285 kJ mol − 1 and the binding energy of the most stable complex [ Ca ( DMF ) 8 ] 2+ is -1592 kJ mol − 1 producing a |ΔE b | − |ΔE d | = 307 kJ mol − 1 , which thus indicates that the formation of ion-pairs and aggregates in DMF is an unlikely process due to its large solvation power. This is in complete agreement with the experimental Raman spectroscopy observations by Forero-Saboya et al. wherein the SN of DMF did not change up to 1.2 M Ca(TFSI) 2 concentration, suggesting that ion-pair and aggregate formation both are negligible [11] . Additionally, the Ca(TFSI) 2 based EC electrolytes had the highest ionic conductivities. We note that the difference |ΔE b | − |ΔE d | between [Ca(DMF) 8 ] 2 + and [Ca(PF 6 ) 3 ] − (458 kJ mol − 1 ) is even larger, but for this Ca-salt there are no corresponding experimental data.

Concluding remarks
By employing a joint DFT + COSMO-RS screening approach we have enabled a rational selection of electrolyte solvent(s) and Ca-salt(s) based on the solvent solvation power, ESWs, viscosity, flash and boiling points. The computed relative solubilities point out sulfur and nitrogen containing solvents as the more promising alternatives and this is successfully correlated with their high DNs that arguably result in a high cation solvation capability. While no silver bullet has been found, the purpose of demonstrating the strategy of computational screening as an efficient step-wise process has been successful -and can be recommended as a time and effort saving tool prior to undertaking any experimental studies to fast evaluate Ca conducting electrolyte candidates. Here, as a compromise with a beneficially wide ESW and a low viscosity, but low boiling and flash points, DMF was selected as the most promising solvent and demonstrated to readily solvate Ca 2 + cations. Furthermore, both Ca(PF 6 ) 2 and Ca(TFSI) 2 are identified as favourable salts to be dissolved in DMF due to their low dissociation energies. While the overall stability of the cation-solvent interactions leads to large desolvation energies, which unfavorably correlates with large experimental activation barriers [59][60][61][62] affecting the rate capabilities of batteries negatively, a step-wise desolvation might be easier as each solvent is bound some-what more weakly, but again the action of the Ca metal anode surface on a partially desolvated Ca 2 + to render a complete plating is very difficult to speculate upon. However, the stable cation-solvent interactions may lower the flash point of DMF based electrolytes, especially at high salt concentrations [22] . Again, however, this alters both structure, dynamics, and transport mechanisms [63] .

Declaration of Competing 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.