Structure, bonding and ionic mobility in Na–V–P–O glasses for energy storage applications

Abstract

Na–V–P–O glasses are promising materials for sodium ion batteries, and yet a thorough understanding of their atomic scale behavior has so far been elusive. In this work we leverage structural and electrochemical experiments with first-principles and large-scale machine learning-accelerated molecular dynamics to elucidate quantitatively the interplay among structure, bonding, and ion mobility on space and time scales of unprecedented extensions. We unravel the existence of a broad V coordination distribution together with heterogeneous Na-ion mobility featuring percolation channels. Our results are instrumental in the search for NVP glass optimization for electrochemical applications.

---

Sodium vanadium phosphate (NVP) oxides have long been a focus of research in energy storage, driven by the electrochemical properties of crystalline NVP phases.^1,2 Nowadays glassy and glass-ceramic counterparts are emerging as valuable alternatives, targeting improvements upon the crystalline forms in terms of synthesis, stability, scalability, and cost.^3–5 In particular, transition metal oxide-containing glasses enable multielectron reactions and impact thermal behavior, enhancing their energy storage potential.^6,7 We present here a novel characterization of NVP glasses that integrates experiments with first-principles and machine learning-accelerated molecular dynamics simulations. Our focus is on the link between structural, bonding, and ion mobility properties in such non-crystalline systems. While detailed structural data are available for binary systems such as V_xO_y–P₂O₅ and Na₂O–V_xO_y (e.g., total structure factors and correlation functions), similar insights for ternary NVP glasses remain scarce.^8,9

On the theoretical side, classical MD (CMD) has provided valuable insights,¹⁰ but current force fields struggle to capture the complex V local environment, particularly VO_n polyhedra, as shown in recent studies of binary VP and ternary NVP glasses.^11,12 In these systems, we find that the account of chemical bonding via first-principles molecular dynamics (FPMD) significantly improves the description of the local V environment. Relevant findings are the presence of distinct single and double/vanadyl bonding fingerprints, the stronger network-forming role of V⁵⁺, and the higher degree of phosphate network polymerization, aligning well with experimental results.

Overall, the present study advances the understanding of NVP glasses by providing a multifaceted characterization of two new representative compositions of Na₂O–V_xO_y–P₂O₅ glasses (NVP10 and NVP28; Table 1 and Table S1 in ESI†). Within this purpose, experimental methods include X-ray diffraction, X-ray photoelectron spectroscopy (XPS), and electrochemical impedance spectroscopy (EIS). On the atomic-scale simulation side, we use a full thermal annealing within Born–Oppenheimer MD (f-BOMD) instead of a short equilibration at 300 K (s-BOMD), as previously reported.^11,12 f-BOMD enables full relaxation of the glass constituents during the melting and cooling processes, resulting in a more accurate medium-range structure compared to s-BOMD, which only allows local relaxation while retaining most of the CMD medium-range order. To significantly extend the reach of space and time scale we take full advantage of machine learning interatomic potential (MLIP) trajectories based on Gaussian approximation potential (GAP), trained on f-BOMD data.^13–15

Table 1 Chemical compositions (in mol (%)) and the glass transition temperature (T_g), along with density (d) of Na₂O–V_xO_y–P₂O₅ system NVP0, NVP10, and NVP28 glasses. See ESI for simulated compositions

System	Na2O (%)	VxOy [V2O4 + V2O5]^a (%)	P2O5 (%)	d (g cm^−3)	T g (C)
a Estimated by XPS (see ESI). b Estimated by He pycnometer (see ESI). c Estimated by DSC (see ESI). d From ref. 8. e From ref. 11.
NVP0^d	—	50.0 [17.6 + 32.4]	50.0	2.800	395^e
NVP10	10.0	40.0 [13.2 + 26.8]	50.0	2.777	439
NVP28	28.5	43.00 [5.59 + 37.41]	28.5	2.939	297

---

Indeed, by comparing the performances of the different models, Fig. 1 underpins the significant improvements of f-BOMD over CMD and s-BOMD for NVP10, in terms of X-rays total structure factors S(k) and reduced total pair correlation functions G(r). These changes are exemplified, for both S(k) and G(r) (Table 2 and Table S5 in ESI†), by the substantial reduction in the R_χ parameter measuring the quantitative agreement with respect to experimental data. In reciprocal space, f-BOMD enhances S(k) by providing a broadened first peak and a minor contribution around 1 Å⁻¹, more accurately capturing intermediate-range order. In direct space and referring to G(r), f-BOMD performs better when compared to experimental data for the first and second peaks, by revealing effectively the presence of double/vanadyl and single bonds, respectively. This holds true also for the fourth peak, corresponding to polyhedral connections involving VV and VP distances. Fig. 2 provides in a comparative fashion the partial pair correlation function g_VO(r) for VO pairs obtained from the different computational schemes. The s-BOMD model features well separated single and double V⋯O bonds, with peaks at ∼1.6 Å (V=O) and ∼1.83 Å (VO), consistent with previous results,^11,12 improving upon the single broadened peak obtained with CMD. This separation is better defined within f-BOMD showing a well-defined minimum between peaks, and a VO bond distance (∼1.88 Å) closer to experimental data (∼1.92 Å).

Fig. 1 X-ray S(k) (top) and G(r) (bottom) of NVP10 glass, comparing calculated data (CMD, s-BOMD, f-BOMD, MLIP with varying model sizes (ESI†)) to experiments.

Table 2 Comparison of the agreement between CMD, s-BOMD, f-BOMD, and MLIP models (small (∼400 atoms) and large (∼3200 atoms)) and experimental data using goodness-of-fit R_χ parameters for the X-ray total structure factor S(k) of NVP0, NVP10, and NVP28 glasses at 300 K, as well as for the neutron total structure factor for NVP0

Method	Size	NVP0		NVP10	NVP28
CMD	Small	11.3 ± 0.1	9.8 ± 0.1	17.2 ± 0.1	17.1 ± 1.0
s-BOMD	Small	8.4 ± 0.2	8.4 ± 0.2	14.4 ± 0.2	9.8 ± 1.2
f-BOMD	Small	—	—	12.8 ± 0.1	—
MLIP	Small	8.7 ± 0.1	9.9 ± 0.1	11.5 ± 0.5	10.4 ± 0.5
	Large	7.1 ± 0.1	8.1 ± 0.1	11.9 ± 0.1	10.3 ± 0.1

---

Fig. 2 Left: Partial pair correlation function g_VO(r) for NVP10 glass obtained by CMD, s-BOMD, f-BOMD, and MLIP with varying model sizes (ESI). Right: Atomistic models by f-BOMD (Na (yellow), V (blue), P (orange), and O (red)) with (top) local spin density isosurfaces (0.05 a.u., in purple) and (bottom) Wannier centers (green) within a V local environment in a VO_n polyhedron.

Fig. 2 also shows a spin topology analysis of the NVP10 model obtained from f-BOMD, characterized by spin localization on V sites corresponding to paramagnetic V⁴⁺ and allowing for a precise speciation between V⁵⁺ and V⁴⁺ sites. Additionally, we provide an atomistic perspective on the unique bonding characteristics of VO_n polyhedra, based on maximally localized Wannier functions (WFCs).¹¹ Within the spin-unrestricted DFT-BOMD framework, this method concurs to identify different bonding types: single VO bonds, defined by two WFCs shared between connected atoms, and double/vanadyl V=O bonds, characterized by more than two WFCs shared between connected atoms. Accurate assessment of the local V environment enables better partitioning bridging and non-bridging oxygen contents (BO and NBO, respectively) in the glass network, a key factor impacting Na ion dynamics.^16–18 In Table 3, we also report the V coordination number (n_V), calculated by integrating the first peak of the g_VO(r), as well as the distribution of individual n_V(l) structural units for V with l-fold coordination in the NVP10 glass. The f-BOMD data provide a novel perspective, distinct from CMD results and building upon s-BOMD data, showing a higher n_V (4.6) and a broad coordination distribution (4 to 6), with a dominant presence of VO₅ polyhedra (∼50%). These findings strongly underscore the role of V as a network former in the glass structure.^11,12

Table 3 Average V coordination number (n_V), bond lengths (r_ij (in Å), taken as the position of the first maximum of the g_αβ(r)), and distribution of individual n_V(l) structural units of V of l-fold coordination computed for NVP10 glass, comparing CMD, s-BOMD, f-BOMD, and MLIP models

	Exp.	CMD	s-BOMD	f-BOMD	MLIP
n V	4.20–5.40	4.33	4.73	4.60	4.93
r VO	1.59–1.78	—	1.59	1.57	1.56
	1.8–2.3 (1.92)	1.78	1.82	1.86	1.86
Structural units
l = 3	VO3	5.1 ± 2	—	—	—
l = 4	VO4	60.6 ± 2	38.4 ± 2	44.7 ± 0.1	19.4 ± 3.3
l = 5	VO5	30.1 ± 2	50.2 ± 3	50.3 ± 0.2	64.8 ± 1.6
l = 6	VO6	4.2 ± 1	11.4 ± 5	5.0 ± 0.1	15.9 ± 2.5

---

To enhance computational efficiency without compromising the accuracy achieved through BOMD, we leveraged accelerated-MD by an MLIP of GAP type, fitted over NVP10 f-BOMD datasets (Fig. S2 in ESI†). Our new MLIP shows a remarkable overall performance in terms of energy, forces, and virials amounting to 5.8 meV per atom, 0.4 eV Å⁻¹ and 16.9 meV per atom, respectively, when calculated with respect to DFT accuracy for testing sets. Close behaviors are found on training and testing sets. Fig. 1 compares the performance of the newly developed MLIP for NVP10 glass, showing that the MLIP model closely matches the f-BOMD results for both S(k) and G(r), accurately describing the NVP10 structure and bonding distances. This is further substantiated by the comparable R_χ values for the MLIP and f-BOMD in 388-atom models, with additional improvements observed when considering a significantly larger system (3104 atoms, clearly a dimension beyond the capability of FPMD/BOMD when aiming at the production of a significant time trajectory). MLIP also reproduces accurately VO and V=O bond distances, n_V, and n_V(l) when compared to f-BOMD data. The MLIP approach allows reducing drastically the computational cost from about 220 days for f-BOMD to 4 days with MLIP (Table S6 in ESI†).

Fig. 3 allows a comparison between the total correlation functions T(r) for NVP0 and G(r) for NVP28 glasses, highlighting the robustness and transferability of our MLIP to compositions beyond the training set (see also S(k) comparison in Fig. S4 and S5 in ESI†). For NVP0, key improvements include the precise description of the first peak near 1.8 Å⁻¹ in S(k) (Fig. S4, ESI†) and all four characteristic peaks in X-rays T(r). Notably, MLIP captures the second peak in T(r) at ∼1.9 Å (V–O bond) and the low-k region of S(k), reflecting its improved ability to model short-range V–O bonds and intermediate-range order. For NVP28, a high-Na₂O composition relevant to energy applications, MLIP outperforms both s-BOMD and CMD, particularly for the first peak in S(k) (Fig. S5, ESI†), the low-k region, a better resolved minimum between the 1st and 2nd peaks of G(r), and the 4th peak of G(r) near 3.2 Å, thereby highlighting realistic VO_n/PO_n polyhedral interconnections.

Fig. 3 Comparison of experimental and calculated X-rays total correlation functions data (CMD, s-BOMD, and MLIP with varying model sizes (ESI)) for NVP0 (T(r), top) and NVP28 (G(r), bottom) glasses.

Fig. S6 in the ESI† is indicative of the performance of our MLIP when describing NVP thermal relaxation, leading to an estimated glass transition temperature (T_g) in the range of 564–665 K when comparing 405- and 3240-atom models, a large improvement with respect to the typical overestimation of CMD (∼1500 K), and fairly in line with respect to the experimental value of ∼570 K.^18,19 The accurate structural and thermal description of NVP28 by our MLIP allows the reliable assessment of the ionic conductivity (σ_ion) via Na diffusion coefficients obtained from mean square displacements (MSD).^20–22 MSD analysis (Fig. 4) reveals a highly heterogeneous dynamical pattern for Na ions, with a temperature- and composition-dependent fraction of them exhibiting negligible mobility and vibrating within local cages (trapped ions), while some atoms do feature significantly higher mobility (highly mobile ions). This heterogeneity correlates with partial Na ordering, forming Na-rich regions and percolation channels. Such ordering, especially at intermediate range distances, is reflected by a peak in the Faber–Ziman partial structure factor at ∼1 Å (Fig. S7, ESI†).²³ At T = 473 K, the NVP28 glass exhibits a total experimental conductivity of 3.3 × 10⁻⁵ S cm⁻¹ (Fig. 4), determined by EIS. The nearly ideal semicircle of the Nyquist plot (Fig. S8, ESI†), though slightly depressed,²⁴ suggests a coexistence of minor ionic (Na⁺) and dominant electronic conductivity, driven by small polaron hopping of V⁴⁺ and V⁵⁺ sites. This is consistent with the relatively low ionic conductivity contribution (∼5%) from Na⁺ ions estimated via MLIP (Fig. 4 and Fig. S6 in ESI†). Notably, Fig. 4 also includes data for two additional NVP glasses with similar Na₂O content but differing by V/P ratios and V speciation, where mixed ionic and electronic conductivity has been reported, with ionic contributions reaching up to ∼50%.^25,26 While the comparable total conductivity across these results is consistent, the structural differences underscore the need for deeper investigations. The new MLIP has the ability to provide an accurate description of the structure and a detailed partitioning of BO/NBO roles in Na-ion dynamics. This study establishes a comprehensive understanding of the structural and dynamic properties of Na–V–P–O glasses, leveraging a novel combination of experiments, first-principles simulations, and ML-accelerated MD. By providing unprecedented insights into the medium-range order, vanadium speciation, and Na-ion transport pathways, this work lays the foundation for the design of high-performance glass-based materials for energy storage applications. The demonstrated accuracy and scalability of the MLIP approach open avenues for extending these methods to other complex amorphous systems, accelerating innovation in glass science and energy technology.

Fig. 4 Left: MSD of Na ions for the NVP28 3240-atom model at 1200 K (MLIP), showing highly mobile (purple) and trapped ions (orange). Center: Snapshot with Na ions color-coded by MSD; PO_n/VO_n units as transparent polyhedra. Insets: Trajectories of a trapped and a highly mobile Na ion, color-coded by MSD (MSD_final ∼ 2 Ų and ∼85 Ų, respectively). Right: Arrhenius plot of log(σ) vs. 1000/T, comparing this work (green, sim. and exp.) with other NVP glasses (grey²⁵ and orange²⁶).

Work supported by Agence Nationale de la Recherche (AMSES ANR-20-CE08-0021, EUR-QMAT) and by Région Nouvelle Aquitaine via the CANaMIAS project (AAPR2021-2020-1177911). Computational resources were provided by HPC GENCI (A0160807670, A0160912441, A0160910832, AD010914978), HPC Mésocentres of the University of Strasbourg (Équipe@Meso, Alsacalcul/Big Data), and MCIA.

Data availability

The data supporting this article are included in the ESI.† Raw and curated simulation data, along with the developed MLIP potential, are available from the NOMAD²⁷ and CNRS Research Data²⁸ repositories.

Conflicts of interest

There are no conflicts to declare.

References

1. F. Kong, et al. , Energy, 2021, 219, 119513.
2. E. Boivin, et al. , Molecules, 2021, 26, 1428.
3. S. Gandi, et al. , New J. Chem., 2020, 44, 2897–2906.
4. Z. Wei, et al. , Adv. Mater. Interfaces, 2018, 5, 1800639.
5. J. Zhang, et al. , Energy Environ. Mater., 2023, 6, e12573.
6. M. S. Whittingham, Chem. Rev., 2014, 114, 11414–11443.
7. R. Giovanardi, et al. , J. Mater. Chem., 2010, 20, 308–313.
8. U. Hoppe, et al. , J. Non-Cryst. Solids, 2012, 358, 328–336.
9. U. Hoppe, et al. , J. Non-Cryst. Solids, 2021, 572, 121120.
10. G. Ori, et al. , J. Non-Cryst. Solids, 2011, 357, 2571–2579.
11. S. Wansi Wendji, et al. , J. Non-Cryst. Solids, 2024, 634, 122967.
12. S. Wansi Wendji, et al. , J. Non-Cryst. Solids, 2025, 655, 123420.
13. A. P. Bartók, et al. , Phys. Rev. Lett., 2010, 104, 136403.
14. F. Shuaib, et al. , J. Am. Ceram. Soc., 2024, 108, e20128.
15. T.-L. Pham, et al., arXiv, 2024, preprint, arXiv:2404.11442 DOI:10.48550/arXiv.2404.11442.
16. Y. Onodera, et al. , NPG Asia Mater., 2019, 11, 75.
17. J. Habasaki, C. Leon and K. L. Ngai, Dynamics of glassy, crystalline and liquid ionic conductors, Springer, Cham, 2016.
18. A. Pedone, et al. , J. Non-Cryst. Solids: X, 2022, 100115.
19. F. Lodesani, et al. , Sci. Rep., 2020, 10, 2906.
20. T.-L. Pham, et al. , J. Mater. Chem. A, 2023, 11, 22922.
21. Theory and simulation in physics for materials applications, ed. E. Levchenko, Y. J. Dappe and G. Ori, Springer, Cham, 2020.
22. C. Massobrio, The structure of amorphous materials using molecular dynamics, Inst. of Phys. Pub., London, 2022.
23. S. Sørensen, et al. , J. Phys. Chem. B, 2023, 127, 10179–10188.
24. M. Saad, et al. , Mater. Res. Bull., 2017, 89, 224–231.
25. A. Sharma, et al. , J. Phys. Chem. C, 2024, 128, 9793–9801.
26. M. Wasiucionek, et al. , Solid State Ionics, 1994, 70–71, 346–349.
27. Data available at the NOMAD data repository, 2025.
28. Data and MLIP potential available at the CNRS Research Data repository (https://recherche.data.gouv), 2025.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5cc00443h

‡ These authors contributed equally to this work.

---