Voltage prediction of vanadium redox flow batteries from first principles

Global energy demand has been increasing for decades, which has created a necessity for large scale energy storage solutions for renewable energy sources. We studied the voltage of vanadium redox flow batteries (VRFBs) with density functional theory (DFT) and a newly developed technique using ab initio molecular dynamics (AIMD). DFT was used to create cluster models to calculate the voltage of VRFBs. However, DFT is not suited for capturing the dynamics and interactions in a liquid electrolyte, leading to the need for AIMD, which is capable of accurately modeling such things. The molarities and densities of all systems were carefully considered to match experimental conditions. With the use of AIMD, we calculated a voltage of 1.23 V, which compares well with the experimental value of 1.26 V. The techniques developed using AIMD for voltage calculations will be useful for the investigation of potential future battery technologies or as a screening process for additives to make improvements to currently available batteries.


Introduction
As lithium-ion batteries are approaching their theoretical capacity limit [2], and increasing demand is putting a strain on the supply chain [3], alternative energy storage technologies must be investigated.One alternative option for large scale battery storage is vanadium redox flow batteries [4].
Vanadium redox flow batteries (VRFBs) stand out among other flow batteries due to their long life cycle [5], thanks in part to their ability to utilize the four oxidation states of Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.vanadium (V 2+ , V 3+ , V 4+ and V 5+ ).A schematic of a VRFB can be seen in figure 1.Like other flow batteries, liquid solutions are held in large tanks and pumped through cell stacks.VRFBs conventionally use components of vanadium salts dissolved in sulphuric acid.Pumping these solutions through the cell stacks allows for hydrogen transfer between the anolyte (V 2+ /V 3+ ) and catholyte (V 4+ /V 5+ ), generating a current as charge is transferred.The use of all-vanadium components allows for excellent stability, even when cross-contamination occurs.Solutions can be re-balanced even when the battery is active, correcting the anolyte or catholyte to its desired charge state.This stability leads to a high number of maximum cycles, lasting 12 000-14 000 cycles [5].For comparison, Liion batteries have a maximum cycle life of 1000-10 000 cycles.
VRFBs are not without their drawbacks though.They are not suitable for mobile applications due to their large size (and corresponding large weight) [5].They are also expensive, with Past studies of VRFBs involving density functional theory (DFT) have focused on the solvation structure of vanadium, the stability of V 5+ and its precipitation as V 2 O 5 .Gupta et al used force field molecular dynamics and well-tempered metadynamics to investigate the interaction of vanadium ions with the molecules in the solvent [11].Vijayakumar et al used DFT and NMR to study the solubility of vanadium ions in electrolyte [12].The contact ion pair between vanadium and sulphuric acid [13] and the solvation structure of vanadium were investigated using DFT by Sepehr and Paddison [14].Ab initio molecular dynamics (AIMD) has also been used to investigate the solvation structure of vanadium in solutions with sulphuric acid and hydrochloric acid [15].Tussupbayev and Kudaibergenova performed benchmark testing on the standard reduction potential of vanadium compounds, where they were able to obtain results within 0.020 V with their most accurate basis set [10], but we were unable to replicate these results.Caro et al have discussed the need for using ab initio molecular dynamics for transition metals in aqueous solution [16].Other studies have also shown accurate redox potentials and acidity constants can be found using DFT based molecular dynamics [17,18].Yu et al have also discussed the use of AIMD to find redox potentials of quinones and their derivatives [19].Those that investigate vanadium focus on the redox potential of V 2+/3+ and have been shown to be difficult to reliably model, often relying on significant explicit solvation as well as implicit solvation.The inclusion of large amounts of explicit solvation can be difficult when using DFT as many configurations must be considered.
Previous attempts to model the voltage of VRFBs relied on mathematical models, and did not utilize first principles atomic modeling [20].To study their low voltage, we developed a novel computational approach, employing ab initio molecular dynamics to model the VRFB system.Existing techniques for modeling the voltage of conventional 'rocking-chair' batteries involve DFT [21][22][23][24][25], and utilize the equation: For example, Eglitis and Borstel were able to demonstrate the use of DFT to predict a 5 V voltage plateau for a Li-ion cell [26,27].While the mechanisms in 'rocking-chair' batteries are different from those of redox flow batteries, this equation is still relevant as it compares the energies of the charged and discharged states.This equation tells us that the voltage of a battery system can be determined by calculating the energy of the discharged state (E discharged ) and the energy of the charged state (E charged ) and dividing by the number of electrons transferred.DFT has been used on liquid systems in some cases, requiring explicit and implicit solvation along with hybrid models [28][29][30][31][32][33].Extensive testing must also be done on the configuration of any explicit molecules outside the first solvation shell as the dynamics of the system cannot be modeled using DFT.To address this issue, AIMD can be used to model the liquid systems, allowing for the movement of molecules and calculation of the system energy at each time step over a simulation.The average energy of the system in its charged and discharged states, taken from the energy values found at each time step, can be used in equation ( 1) to find the voltage of liquid battery systems.

AIMD method
To determine the full contents of our simulation cell, we began with the initial components used in an experimental cell.VOSO 4 is mixed with a solution of sulphuric acid, which then goes through a full charge and discharge cycle to form the solution for the other charge states [34].The first reaction equation we consider is the conversion from the V 4+ to V 5+ , written in equation ( 2), where the left side of the equation shows the system in the V 4+ charge state and the right side shows the system in the V 5+ charge state.This reaction shows that water gives an O atom to VO 2+ .Half of the H + are transferred across the proton exchange membrane, while the other half will react with water to become H 3 O + in our simulation cell.The V 4+ to V 3+ reaction is written in equation ( 3), where the aH + term on the left side of the reaction comes from the electrolysis needed for the reaction to occur.On the right side of the equations we can see this introduction of aH + causes VO 2+ to lose an O atom, and H 2 SO 4 to lose H protons to form water.The final reaction equation we must consider is V 3+ to V 2+ , written in equation ( 4), where the introduction of aH + from the charging process converts SO 2− 4 to H 2 SO 4 .Now that we have determined the ratio of the contents in each charge state, we must determine the density and the value of coefficients a, b and c, which are determined by the molarity of the components.The molarities chosen were 2 M for VO + 2 , 2 M for SO 2− 4 and 3.5 M for H 2 SO 4 for our starting V 4+ system. (3) These molarities were chosen based on the experimental values given by Gupta et al [11], and were within the operational concentrations listed in other works [7,8,14,34,35].We chose to use a 15 × 15 × 15 Å 3 simulation cell.This cell size was selected as it provides sufficient volume to include four vanadium ions and is not so large as to be prohibitively expensive to perform AIMD.With the molarities listed and the volume of our simulation cell, a = 4 and b = 7.The value for c was then determined by considering the density of the system.The density of our system was based on equation (5), where c V 3+ is the molarity of V 3+ (2 M) and c SO 2− 4 is the molarity of SO 2− 4 (5.5 M) [8].With the values of a and b already determined, the value of c was now selected to reach the correct density of the system, leading us to c = 68.We prepared AIMD simulations were run with the VASP program.We used the Nosé-Hoover thermostat [37][38][39][40] to control the velocity of atoms, a plane wave cutoff energy of 400 eV, and DFT-D3 dispersion correction [41][42][43].All AIMD trajectories were run with a 0.5 fs time step for half of a picosecond to help with initial equilibration.After this setup phase, the time step was increased to 1 fs, allowing for the trajectories to be run for longer timescales without using excessive computational resources.An example system can be seen in figure 2.
Since V 5+ is not stable in molarities above 2 M, and is only stable within the temperature range of 10 • C-40 • C [7,8,11,14,34,35], a temperature of 27 • C (300 K) was used.Outside these conditions V 5+ precipitates out of the solution, resulting in capacity loss for the system.

Cluster model
Before employing this computationally expensive method using AIMD, a cluster model was used to attempt to model the voltage of the VRFB anolyte with DFT.All cluster models are performed in a vacuum.
The experimental voltage of the anolyte has been shown [6][7][8][9][10] as: By only considering the charged vanadium atoms in each of their respective states, we adjusted equation (1) to the form: where E V 3+ is the discharged state of the anolyte (hydrated vanadium ion with a charge of +3), E V 2+ is the charged state of the anolyte (hydrated vanadium ion with a charge of +2), and it is divided by one since the charge is changed by one electron.This method allows us to use only small molecular clusters with a finite amount of explicit solvent.

Orca results/discussion
We started by considering small cluster models of the solvated vanadium ions with the Orca program [44].Various solvation structures were optimized with the B3PW91/6-31G method [45], beginning with an empty first solvation shell, followed by the addition of water going from 1-10 molecules.
The results from these hydrated vanadium ions were then used in equation ( 7) to calculate their redox potentials.Example systems can be seen in figure 3. Single point calculations were then run with the optimized coordinates with the method B3PW91/6-311++G(d,p) with a dispersion correction (D3BJ) [46][47][48].Single point calculations were also run with a continuum solvent model [49], where the the dielectric constant for water was used to implicitly surround each structure with water.The vibrational frequencies were also calculated with the continuum solvent model with V 3+ as a triplet (as there are two electrons in the 3d-orbital), and with V 2+ as either a doublet (if we consider the 3d-orbital to have one electron and the 4s-orbital to be full) or a quartet (if there are three electrons on the d-orbital, and none in the 4s-orbital) with six water molecules explicitly in the solvation structure.These frequency calculations were performed to account for the entropy contributions to the Gibbs free energy of the systems.
We began our investigation testing the effects of having water in the solvation structure and the basis set used to perform our calculations.Figure 4 shows us that the difference in the voltage between a smaller basis (6-31G, red) and a larger basis (6-311G++(d,p), blue) was negligible, which may be due to some cancellation of error, since we are using the difference in energy for the voltage calculation.The inclusion of the continuum solvent model (yellow) showed significant improvement, but still did not provide an accurate result when compared to the experimental value of −0.255 V (horizontal dashed green line).
The presence of water in the solvation shell also had a profound effect on the calculated voltage, demonstrating the importance of the solvation structure in liquid systems.The effects of the number of water molecules becomes less important once the system has reached six water molecules, which is consistent with the full first solvation shell of V 3+ [11].We ran our frequency calculations only with six waters in the solvation structure, as we no longer observed significant improvement after six waters.These frequency calculations included the entropy term of the Gibbs free energy, but we were still unable to obtain significant improvements to the calculated voltage.
To better account for the full contents of VRFB anolyte, multiple terms were added to the voltage equation to consider the effects of each of the following: sulfate, bisulfate, and sulphuric acid.As the charge transfer is the result of hydrogen transfer across a proton exchange membrane, sulphuric acid (H 2 SO 4 ) becomes bisulfate (HSO − 4 ) or bisulfate becomes sulfate (SO 2− 4 ) during the discharge process.Our new voltage equation is now, where the first three terms in the numerator represent the discharged state and the last three terms represent the charged state.The value of six water molecules in each structure was chosen based on the results shown in figure 4.
The variable x can be adjusted while maintaining the balanced charge of the system and stoichiometry, but changes the amount the solvent molecules (sulphuric acid, bisulfate and water) are contributing to the voltage equation.Table 1 shows how the voltage changes based on the variable 'x' in equation (8).A value of x = 4 was able to bring our voltage within about 0.1 V of the experimental value.While this result is encouraging, this value for x is arbitrary and justification for its selection as our final result would not be fair.To get a more complete understanding of the voltage, experimental values for the molarity of the system should be used.This consideration leads to our AIMD analysis of the voltage.Table 1.The voltage of the anolyte system with various sulphuric acid considerations.The variable 'x' from equation ( 8) is altered to show the effect of the concentration of sulphuric acid. x Voltage (V) 1.14

AIMD results/discussion
With the use of AIMD, we are able to model a more complete and accurate electrolyte system accounting for the interactions and reactions that occur in a liquid system.The energies of these systems can then be used to calculate the voltage of the total system.To calculate the total voltage of the VRFB system, both the anolyte and catholyte need to be considered.The voltage determined in experiment [6][7][8][9] is as shown in equations ( 9)- (11), where the left side of the equation shows the charged state and the right side shows the discharged state.
For simplicity, VO + 2 will be referred to as V 5+ and VO 2+ will be referred to as V 4+ .Our overall voltage equation simply uses equation ( 1) in a more system-specific form, shown in equation (12), where the terms in the numerator describe the energies of the electrolytes with vanadium in each of its four charge states.For our systems, N electrons = 4.The average energy of the AIMD trajectories over two picoseconds were used to calculate the voltage, shown in figure 5.It should also be noted that the V 2+ system required some restriction on bond distances to prevent an unphysical interaction between a vanadium ion and a hydrogen atom.This was not observed in any of the other simulations, and no literature on V-H bonding was found.The bond distance between the hydrogen atom and the water it was bonded with was constrained to prevent the H atom from dissociating and forming a H-V bond.The energy of the V 2+ system without these restrictions can be seen in the supplementary material.
As we can see in figure 5, the relative energies undergo large changes for the first 5 ps, while the system is equlibrating.After this initial phase, the relative energies change much less, which has a noticeable impact on the variability of the calculated voltage.While the voltage begins at a value of 0.7 V, this value increases to above 1 V after the first 5 ps of the trajectory.The voltage reaches a value of 1.2 V by the time the trajectory has been run for 8 ps, and remains within 0.1 V of the experimental voltage value for the rest of the trajectory.At 10 ps, we were able find a value of 1.23 V for the VRFB system.Our calculated voltage is only 0.03 V away from the experimental voltage of 1.26 V. Despite the excellent agreement, we acknowledge that this result may be fortuitous and that there may be some cancellation of errors.One possible source of error is the limited AIMD trajectory accessible to our calculations.We also tested other values for the rolling average which can be seen in the supplementary material.The use of the 2 ps rolling average was based on our observations of the variation in the energy and voltage values.A 1 ps rolling average has high variation, as slight changes are noticeable over the smaller time-scale.A 5 ps rolling average does not react as quickly to changes, requiring longer trajectories before variations have a noticeable effect and also require longer trajectories before the equilibrium phase is contained within an average.We believe that the 2 ps rolling average offers a nice balance to capture the average value, while not requiring the system to be run excessively past the point of the systems reaching their equilibrium phase.

Summary
Vanadium redox flow batteries are a promising large scale energy storage technology, but low voltage and narrow operating conditions are obstacles that must be overcome before they can become a viable large scale energy solution.We are able to model these systems using ab initio atomistic modeling to better understand the interactions occurring and how each component affects the voltage.From our results using DFT on cluster models, we can see the importance of considering all components of a liquid system.Only once all components and their concentrations were considered, were we able to obtain a voltage within 0.1 V from the experimental value for the anolyte using DFT, but could not justify this value as the final result due to the arbitrary inclusion of the variable 'x' in equation (8).To consider all components in the system while maintaining accurate concentrations, AIMD can be used to model the interactions and find an average energy value for the liquid system.We calculated a voltage of 1.23 V, only 0.03 V from the experimental voltage value.With the help of this new technique for determining voltages, we hope this process can be applied to other systems, or as a screening process for any additives that could help improve existing redox flow batteries, including the testing of higher ion concentrations, as well as the voltage at various states of charge.Other flow batteries could be modeled using this same technique, or other battery types that require liquid electrolyte contributions to system voltage, such as dual-ion intercalation batteries [22,[50][51][52].

Figure 1 .
Figure 1.Schematic showing the components of a VRFB.Blue arrows indicate the direction of flow for the catholyte/anolyte as they are pumped from external tanks into the cell stacks.Note that the catholyte (anolyte) is a single solution containing both types of ions, V 4+ and V 5+ (V 2+ and V 3+ ).They are shown here separately to illustrate the redox conversions that occur, made possible by the proton transfer across the proton exchange membrane and the electron current from the external circuit.

Figure 3 .
Figure 3. Example structures of water being added to the solvation structure of a vanadium ion.Vanadium ion with (a) one, (b) six, (c) ten water molecules.

Figure 4 .
Figure 4.The voltage of the anolyte system as a results of DFT structure optimizations for V 3+ and V 2+ surrounded by water.The dashed green line at −0.255 V represents the experimental value [6-9].

Figure 5 .
Figure 5. Average energy and voltage calculated over the preceding 2 ps, as a function of time.The dashed green line at 1.26 V represents the experimental voltage.