Fast Mg²⁺ diffusion in Mo₃(PO₄)₃O for Mg batteries

Abstract

In this work, we identify a new potential Mg battery cathode structure Mo₃(PO₄)₃O, which is predicted to exhibit ultra-fast Mg²⁺ diffusion and relatively high voltage based on first-principles density functional theory calculations. Nudged elastic band calculations reveal that the migration barrier of the percolation channel is only ∼80 meV, which is remarkably low, and comparable to the best Li-ion conductors. This low barrier is verified by ab initio molecular dynamics and kinetic Monte Carlo simulations. The voltage and specific energy are predicted to be ∼1.98 V and ∼173 W h kg⁻¹, respectively. If confirmed by experiments, this material would have the highest known Mg mobility among inorganic compounds.

---

A promising strategy to increase the energy density of rechargeable batteries is to transition from monovalent to multivalent batteries, such as Mg batteries,^1–3 in which Mg²⁺ is reversibly inserted into/extracted from the cathode. Because dendrites are less likely to grow on Mg metal than Li metal during charging, Mg metal can be directly used as the anode, resulting in a substantial improvement of the theoretical volumetric energy density from 760 mA h cm⁻³ for Li–graphite to 3833 mA h cm⁻³ for Mg metal.⁴

One of the scientific challenges for the development of Mg-ion batteries is the limited mobility of Mg²⁺ in cathode structures. Unlike its alkali cation cousins Li⁺ and Na⁺,^5–7 the higher charge of Mg²⁺ makes it much harder to overcome the migration barrier. In addition, the difficulty of developing Mg anodes and electrolytes that are compatible with candidate cathode materials has posed a challenge for the experimental evaluation of the electrochemical performance of interesting cathode candidates.⁸ As a result, the search for Mg cathode structures has only yielded a few materials that exhibit reasonably fast reversible electrochemical Mg²⁺ intercalation, i.e., Chevrel Mo₆S₈ (∼135 mA h g⁻¹ at ∼1.0–1.3 V),³ orthorhombic V₂O₅ (∼150 mA h g⁻¹ at ∼2.3–2.6 V),^5,9 MoO3 (∼220 mA h g⁻¹ at ∼1.7–2.8)⁵ and spinel TiS₂ (∼200 mA h g⁻¹ at ∼1–1.5 V).^10,11 However, the low voltages limit the power densities of sulfides. Oxides and phosphates generally have higher voltages but at the cost of poorer Mg²⁺ diffusivity.

In this work, we identify a new possible phosphate compound Mo₃(PO₄)₃O, which is shown to exhibit ultra-fast Mg²⁺ diffusion and relatively high voltage based on first-principles density functional theory (DFT) calculations. First-principles calculations have proven to be accurate and effective in studying the voltage and mobility of Li-ion^12–19 and multivalent electrode materials.^6,7,20 Our first-principles nudged elastic band (NEB) calculations^21–23 predict that Mo₃(PO₄)₃O has an unusually low Mg migration barrier of ∼80 meV, which is lower than the values previously reported for spinel TiS₂ (∼550 meV)¹⁰ and Chevrel phases (∼360 meV),¹¹ suggesting that this structure may enable very high Mg²⁺ diffusivity.

The crystal structure of MgMo₃(PO₄)₃O is shown in Fig. 1. MgMo₃(PO₄)₃O was derived from known compounds such as CaFe₃(PO₄)₃O,²⁴ SrFe₃(PO₄)₃O²⁵ and Bi_0.4Fe₃(PO₄)₃O²⁶ by substituting Fe by Mo and the other metal ion (Ca, Sr, Bi) by Mg. The structure is relaxed in both lattice parameters and atomic positions after direct substitutions. Apart from placing Mg in the Ca site in CaFe₃(PO₄)₃O (site B in Fig. 1), we also investigated other possible sites for Mg to reside in the empty host structure Mo₃(PO₄)₃O. One site (site A in Fig. 1) is 9.8 meV lower in energy than site B. As shown in Fig. 1, there are one A site and two B sites in the unit cell, and site A and site B are too close to accommodate Mg cations simultaneously. Hence our calculations predicted that Mg resides in site A at Mg_0.5Mo₃(PO₄)₃O composition and in site B in MgMo₃(PO₄)₃O. Chains composed of edge-sharing MoO₆ octahedra along the [010] direction form the backbone of the structure and are interconnected by predominantly corner-sharing MoO₅ trigonal bipyramids, MoO₄ tetrahedra, and PO₄ tetrahedra. Except for the [Mo₆O₂₈] chains, there are only two edge-sharing links in the unit cell, which are between PO₄ tetrahedra and MoO₆ octahedra. All the other links between polyhedra are corner-sharing links, which enables the polyhedra to rotate slightly, creating an adaptive tunnel to facilitate Mg²⁺ diffusion. Similar to LiFePO₄,²⁷ the MgMo₃(PO₄)₃O structure has a 1D diffusion channel along the b-axis, along which Mg²⁺ diffusion is expected to occur.

Fig. 1 Crystal structure of MgMo₃(PO₄)₃O in (a) b–c and (b) a–c planes. The structure is built on [Mo₆O₂₈] chains interconnected by predominantly corner-sharing polyhedra. Similar to LiFePO₄, a 1D diffusion channel exists along the b-axis for Mg²⁺.

Fig. 2 shows the minimum energy paths for Mg²⁺ migration in Mo₃(PO₄)₃O, as calculated using the NEB method. Mg²⁺ can follow two paths for migration from one stable site to the nearest equivalent site, an inner- or inter-channel path. The inner-channel path involves migration along the b-axis direction and across the unit cell boundary (the unit cell is shown in Fig. 1), with a very low activation barrier of ∼80 meV (Fig. 2(a1 and a2)). The inter-channel path involves migration along the c-axis direction with a much higher activation barrier of ∼1200 meV (Fig. 2(b1 and b2)). Because diffusivity scales as the inverse exponential of the activation barrier, migration along the inter-channel path is unlikely as at room temperature hopping along the c-direction should be ∼10¹⁸ times less frequent than along the inner-channel path (an increase of 60 meV in the migration barrier corresponds to a decrease of one order of magnitude in the diffusivity at room temperature).⁷ The percolation channel for Mg²⁺ intercalation is therefore along the inner-channel path. The inner-channel path is divided into two segments, as shown in Fig. 2(a1) and (a2), with migration lengths and activation energies of 5.32 and 7.00 Å and ∼80 and ∼70 meV, respectively. Multiple minima exist along both paths because the PO₄ groups can easily rotate to accommodate the Mg at different positions. The flexibility of the PO₄ groups is a result of the corner-sharing connection between polyhedrons in the structure.

Fig. 2 The minimum energy paths for Mg²⁺ migration in Mo₃(PO₄)₃O. (a1) Migration path and (a2) minimum energy path for inner-channel diffusion. (b1) Migration path and (b2) minimum energy path for inter-channel diffusion. A–D in (a2) are markers for the migration processes. A and B corresponds to site A and B in Fig. 1. The inner-channel diffusion in (a1) and (a2) is the percolation channel for Mg²⁺ intercalation.

Tables 1 and 2 display the calculated thermodynamic electrochemical properties of MgMo₃(PO₄)₃O. The redox couple during charging and discharging is Mo³⁺/Mo⁴⁺. The average voltage calculated with the GGA(PBE)+U functional is 1.98 V, resulting in a specific energy of 173 W h kg⁻¹. The voltage data and specific energy data were verified using SCAN³¹ and Heyd–Scuseria–Ernzerhof (HSE)³² functionals. An average voltage of 1.69 V was obtained using the HSE functional, which is generally recognized as the most reliable voltage assessment method.³³ The voltage from the SCAN functional was close to that of HSE.

Table 1 Average voltages of MgMo₃(PO₄)₃O calculated using different levels of theory

Functionals	Average voltage (V)	Specific energy (W h kg^−1)
GGA+U (Mo U = 4.38 eV)^28–30	1.98	173
SCAN^31	1.52	133
HSE06^32	1.69	148

---

Table 2 Thermodynamic electrochemical properties of MgMo₃(PO₄)₃O

Mo3(PO4)3O energy above hull	MgMo3(PO4)3O energy above hull	Volume change in charging/discharging
42 meV per atom	36 meV per atom	2%
Volumetric capacity	Energy density	Gravimetric capacity
330 A h L^−1	651 W h L^−1	87 mA h g^−1

---

Since MgMo₃(PO₄)₃O is a hypothetical compound derived by substituting ions in other compounds, we also evaluate its relative phase stability, by constructing the Mg–Mo–P–O phase diagram using available compounds in the Materials Project Database.³⁴ Both the charged and discharged structures are metastable, with a moderate energy above the energy hull of 42 and 36 meV per atom, respectively. These energies measure the driving force to decompose into other phases. A broad analysis of the known compounds in ICSD (Inorganic Crystal Structure Database) has indicated that this energy range for metastability is quite common, indicating the possibility that this compound may be synthesizable.³⁵ The charged material is unstable against decomposition to MoO₂ and Mo₂P₃O₁₁, and the discharged material is metastable with respect to MoO₂, Mo₂P₃O₁₁, and Mg₃(PO₄)₂. Based on the energy above the hull for experimentally synthesized Chevrel phases, i.e., Mo₆S₈ (67 meV per atom) and Mg₂Mo₆S₈ (49 meV per atom), the stability of MgMo₃(PO₄)₃O is considered reasonable for synthesis. Since MgMo₃(PO₄)₃O is a compound that is derived by a double substitution from CaFe₃(PO₄)₃O, and is metastable with the composition given, a carefully constructed route to synthesize it will be necessary. We suggest that one starts with a more stable [X]Mo₃(PO₄)₃O compound, where X is another working ion such as Ca²⁺, Zn²⁺, Li⁺, etc. and then exchange X with Mg²⁺ at elevated temperature. Finally, the volume change during charging and discharging is very small (∼2%), which is highly favorable for the reliability of electrodes.

To understand the unusually high mobility of Mg²⁺ in this compound, and large difference in activation barriers between the inner- and inter-channel paths, we analyzed the coordination number for Mg²⁺ in every image of the two paths. Mg²⁺ almost always maintains 4-fold coordination along the inner-channel path, which differs greatly from its behavior along the inter-channel path, where Mg²⁺ experiences a coordination number change of 4 → 2 → 4. A larger coordination number change has been previously shown to lead to larger site energy differences along the migration path and ultimately to a larger activation barrier.⁷

Moreover, the corner-sharing connections of most of the polyhedra in the structure facilitate Mg²⁺ migration by enabling rotation to accommodate the presence of a local Mg ion. From this perspective, the inner-channel along the b-axis is more advantageous than the inter-channel along the c-axis because rotation of the MoO₆ octahedra along the edge-sharing [010] [Mo₆O₂₈] chains is much easier around the b-axis than the c-axis. In addition, the void in the middle of the inter-channel path is too open, and moderate polyhedral rotations cannot mediate the coordination number decrease.

To verify the migration barriers obtained from zero-K NEB calculations, we performed ab initio molecular dynamics (AIMD) simulations.^36,37 The mean square displacement (MSD) at 650 K is plotted in Fig. 3. The main contribution of the displacement is along the b-axis, which further confirms the 1D diffusion channel, as indicated by the inner-channel path. The diffusivity calculated using AIMD is 2.82 × 10⁻⁵ cm² s⁻¹.

Fig. 3 MSD of Mg²⁺ at 650 K from AIMD simulations. The main contribution of the displacement is along b-axis, which shows a one-dimensional diffusion channel, as indicated by the inner-channel path in Fig. 2(a1) and (a2).

For comparison, a kinetic Monte Carlo (kMC)^38,39 simulation was conducted based on the inner-channel path in Fig. 2(a2). Hopping rates between local minima were calculated using harmonic transition state theory:⁴⁰

where N is the number of atoms; υ_i and υ_i* are the positive normal mode frequencies at the local minimum and transition state, respectively; and ΔE is the energy barrier, as determined from NEB calculations. The pre-factors and barriers between minima A, B, C, and D in Fig. 2(a2) are listed in Table 3. The other half of the path is symmetrically equivalent to ABCD.

Table 3 Summary of hopping processes

Process	Calculated prefactor (THz)	Barrier (meV)
A → B	1.59	70.4
B → A	0.58	60.6
B → C	0.46	66.5
C → B	2.75	16.5
C → D	8.18	29.5
D → C	2.24	46.3

---

The diffusivity at 650 K determined using kMC is 13.70 × 10⁻⁵ cm² s⁻¹. The fact that the diffusivities from AIMD and kMC are within one order of magnitude confirms that the effective barrier of diffusion is low. The room temperature (300 K) Mg diffusivity estimated using kMC is 4.68 × 10⁻⁵ cm² s⁻¹.

In summary, the key property for developing Mg battery cathode materials is the Mg²⁺ cation mobility in the host structure. In this work, we show that Mo₃(PO₄)₃O exhibits extraordinarily fast Mg²⁺ cation mobility based on NEB (activation barrier ∼80 meV), AIMD (diffusivity ∼2.82 × 10⁻⁵ cm² s⁻¹ at 650 K), and kMC (diffusivity ∼13.70 × 10⁻⁵ cm² s⁻¹ at 650 K, ∼4.68 × 10⁻⁵ cm² s⁻¹ at 300 K) simulations. This is to our knowledge the lowest migration barrier ever predicted for Mg²⁺ in an oxide. Its voltage is slightly higher than previously reported sulfides based on GGA+U (1.98 V), SCAN (1.52 V), and HSE06 (1.69 V) calculations, but the capacity of 91 mA h g⁻¹ is relatively low. Our systematic first-principles studies indicate that Mo₃(PO₄)₃O may be a promising 1D cathode material for Mg batteries and is worthy of possible experimental investigation. In addition, the unusually high predicted mobility indicates that while Mg²⁺ diffusion generally is slow in inorganic compounds, there may be notable exceptions.

This work was supported by the Joint Center for Energy Storage Research (JCESR), an Energy Innovation Hub funded by the U.S. Department of Energy, Office of Science and Basic Energy Sciences (Subcontract 3F-31144). This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357.

Notes and references

1. P. Canepa, S. Gautam, D. C. Hannah, R. Malik, M. Liu, K. G. Gallagher, K. A. Persson and G. Ceder, Chem. Rev., 2017, 117, 4287–4341.
2. J. Muldoon, C. B. Bucur and T. Gregory, Chem. Rev., 2014, 114, 11683–11720.
3. D. Aurbach, Z. Lu, A. Schechter, Y. Gofer, H. Gizbar, R. Turgeman, Y. Cohen, M. Moshkovich and E. Levi, Nature, 2000, 407, 724–727.
4. H. D. Yoo, I. Shterenberg, Y. Gofer, G. Gershinsky, N. Pour and D. Aurbach, Energy Environ. Sci., 2013, 6, 2265–2279.
5. G. Gershinsky, H. D. Yoo, Y. Gofer and D. Aurbach, Langmuir, 2013, 29, 10964–10972.
6. M. Liu, Z. Rong, R. Malik, P. Canepa, A. Jain, G. Ceder and K. A. Persson, Energy Environ. Sci., 2015, 8, 964–974.
7. Z. Rong, R. Malik, P. Canepa, G. Sai Gautam, M. Liu, A. Jain, K. Persson and G. Ceder, Chem. Mater., 2015, 27, 6016–6021.
8. P. Canepa, G. S. Gautam, R. Malik, S. Jayaraman, Z. Rong, K. R. Zavadil, K. Persson and G. Ceder, Chem. Mater., 2015, 27, 3317–3325.
9. G. Amatucci, F. Badway, A. Singhal, B. Beaudoin, G. Skandan, T. Bowmer, I. Plitz, N. Pereira, T. Chapman and R. Jaworski, J. Electrochem. Soc., 2001, 148, A940–A950.
10. X. Sun, P. Bonnick, V. Duffort, M. Liu, Z. Rong, K. A. Persson, G. Ceder and L. F. Nazar, Energy Environ. Sci., 2016, 9, 2273–2277.
11. M. Liu, A. Jain, Z. Rong, X. Qu, P. Canepa, R. Malik, G. Ceder and K. A. Persson, Energy Environ. Sci., 2016, 9, 3201–3209.
12. S. P. Ong, V. L. Chevrier, G. Hautier, A. Jain, C. Moore, S. Kim, X. Ma and G. Ceder, Energy Environ. Sci., 2011, 4, 3680–3688.
13. Y. Mo, S. P. Ong and G. Ceder, Chem. Mater., 2014, 26, 5208–5214.
14. G. K. P. Dathar, D. Sheppard, K. J. Stevenson and G. Henkelman, Chem. Mater., 2011, 23, 4032–4037.
15. T. Song, H. Cheng, H. Choi, J.-H. Lee, H. Han, D. H. Lee, D. S. Yoo, M.-S. Kwon, J.-M. Choi and S. G. Doo, ACS Nano, 2011, 6, 303–309.
16. K.-S. Park, P. Xiao, S.-Y. Kim, A. Dylla, Y.-M. Choi, G. Henkelman, K. J. Stevenson and J. B. Goodenough, Chem. Mater., 2012, 24, 3212–3218.
17. D. A. Tompsett and M. S. Islam, Chem. Mater., 2013, 25, 2515–2526.
18. A. Urban, D.-H. Seo and G. Ceder, npj Comput. Mater., 2016, 2, 16002.
19. A. Jain, Y. Shin and K. A. Persson, Nat. Rev. Mater., 2016, 1, 15004.
20. G. Sai Gautam, P. Canepa, A. Abdellahi, A. Urban, R. Malik and G. Ceder, Chem. Mater., 2015, 27, 3733–3742.
21. G. Mills and H. Jónsson, Phys. Rev. Lett., 1994, 72, 1124.
22. G. Mills, H. Jónsson and G. K. Schenter, Surf. Sci., 1995, 324, 305–337.
23. Z. Rong, D. Kitchaev, P. Canepa, W. Huang and G. Ceder, J. Chem. Phys., 2016, 145, 074112.
24. M. Hidouri and M. Ben Amara, Acta Crystallogr., Sect. E: Struct. Rep. Online, 2009, 65, i66.
25. V. Morozov, K. Pokholok, B. Lazoryak, A. Malakho, A. Lachgar, O. Lebedev and G. Van Tendeloo, J. Solid State Chem., 2003, 170, 411–417.
26. A. Benabad, K. Bakhous, F. Cherkaoui and E. M. Holt, Acta Crystallogr., Sect. C: Cryst. Struct. Commun., 2000, 56, 1292–1293.
27. R. Malik, D. Burch, M. Bazant and G. Ceder, Nano Lett., 2010, 10, 4123–4127.
28. Z. Wu and R. E. Cohen, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 235116.
29. A. Jain, G. Hautier, S. P. Ong, C. J. Moore, C. C. Fischer, K. A. Persson and G. Ceder, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 045115.
30. L. Wang, T. Maxisch and G. Ceder, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 195107.
31. D. A. Kitchaev, H. Peng, Y. Liu, J. Sun, J. P. Perdew and G. Ceder, Phys. Rev. B: Condens. Matter Mater. Phys., 2016, 93, 045132.
32. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215.
33. G. Ceder, MRS Bull., 2010, 35, 693–701.
34. S. Ping Ong, L. Wang, B. Kang and G. Ceder, Chem. Mater., 2008, 20, 1798–1807.
35. W. Sun, S. T. Dacek, S. P. Ong, G. Hautier, A. Jain, W. D. Richards, A. C. Gamst, K. A. Persson and G. Ceder, Sci. Adv., 2016, 2, e1600225.
36. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558.
37. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 48, 13115.
38. A. Chatterjee and D. G. Vlachos, J. Comput.-Aided Mater. Des., 2007, 14, 253–308.
39. G. Henkelman and H. Jónsson, J. Chem. Phys., 2001, 115, 9657–9666.
40. A. F. Voter and J. D. Doll, J. Chem. Phys., 1984, 80, 5832–5838.

---

Footnotes

† Electronic supplementary information (ESI) available: Computational methods details for the results presented. See DOI: 10.1039/c7cc02903a

‡ Ziqin Rong and Penghao Xiao contributed equally to this work.

---