Dataset for article "Ab initio investigation of electronic and lattice properties of Fe4(P2O7)3"

Ab initio investigation of electronic and lattice properties of Fe4(P2O7)3 Svitlana Pastukha, Pawe l T. Jochyma, Jan La˙zewskia, Dominik Legutb,c,∗and Przemys law Piekarza aInstitute of Nuclear Physics Polish Academy of Science, PL-31342 Krak´ow, Poland bDepartment of Condensed Matter Physics, Faculty of Mathematics and Physics, Charles University CZ-121 16, Prague, Czech Republic cIT4Innovations, VSB-Technical University of Ostrava CZ-708 00, Ostrava, Czech Republic In this research, we examine the electronic, magnetic, and lattice properties of the Fe4(P2O7)3 compound using the first principles calculations based on the density functional theory. The crystal lattice has a monoclinic structure, belonging to the P21/n space group. The optimized lattice parameters are a=7.406 ˚ A, b=21.425 ˚ A, c=9.529 ˚ A, and agree very well with the experimental data, thanks to the local Coulomb interactions and van der Waals forces included in the calculations. The investigation considers several magnetic orderings. The lowest total energy was found for the antiferromagnetic configuration with the magnetic moment of ∼4.6µBper Fe atom. The electronic structure calculation shows the Mott insulating state with the energy gap Eg= 2.87 eV. For the relaxed crystal structure, the elastic properties were obtained and analyzed. The phonon dispersion relations and density of states were calculated within the temperature-dependent effective potential methodusing atomic multidisplacements obtained by high efficiency configuration space sampling. I. INTRODUCTION The iron pyrophosphate compounds are of great scientific and technological interest due to their diverse chemical and physical properties and a broad range of applications. Iron pyrophosphate shows high chemical stability in a solution and can be applied as an active sorbent for the removal of radionuclides from water [1]. The recent study revealed that iron pyrophosphate doped with carbon can be successfully used to control tetracycline antibiotic pollution in aqueous solution by peroxymonosulfate activation [2]. An enhanced photoelectrochemical water oxidation activity of nanocrystalline iron pyrophosphate in amorphous iron phosphate overlayer was also reported [3]. A new compound of iron pyrophosphate with added cobalt carbonate was shown to be prospective for environmental applications [4]. These compounds are also intensively studied for biomedical applications. Numerous studies indicate the important role of ferric pyrophosphate in iron deficiency treatment and iron transport control in human body [5–7]. In addition, iron pyrophosphate has attracted attention due to its potential applications in electrochemical energy storage devices, such as lithium-ion batteries. Porous iron pyrophosphate nanostructure was shown to possess an enhanced electrocatalytic activity in Li-O2batteries [8], while lithium iron pyrophosphate compound synthesized by conventional solid-state reaction was reported as promising cathode for a lithium-ion battery system for large-scale applications [9]. Although the above studies demonstrate a broad range of iron pyrophosphate applications, there is a limited number of investigations concerning its crystal structure, electronic, and magnetic properties. There are four types of such compounds reported in the literature: Fe2P2O7, ∗Contact author: [email protected] Fe3(P2O7)2, Fe4(P2O7)3and Fe7(P2O7)4[10–13], varying by the crystal structure and the way the octahedral iron ions link to each other. In particular, to the best of our knowledge, there are no ab initio, theoretical studies of electronic and lattice dynamical properties of these materials. In the current study, we focus on the Fe4(P2O7)3compound with the monoclinic P21/n symmetry. The structural description and chemical characterization of this material was done by Pahdi et al. [12], who reported presence of two distinct Fe2O9groups of face shared octahedra and three kinds of the P2O7groups. More detailed analysis of the crystal structure and magnetic properties of Fe4(P2O7)3was done by Elbouaanani et al. [13], using neutron diffraction and M¨ossbauer spectrometry techniques. The authors reported the existence of four antiferromagnetic iron sublattices, which correspond to the four crystallographically distinct positions of iron atoms. Also, the existence of a paramagnetic to antiferromagnetic transition was observed at about T= 50 K. This phase transition, as a result of Fe3+ magnetic ordering, was also confirmed in the low temperature heat capacity study [14]. In general, iron pyrophosphate research covers a wide range of experimental topics such as synthesis methods, analysis of structure and magnetic properties, and exploration of catalytic applications. However, there is currently a lack of ab initio studies that delve into the electronic band structure and the impact of local Coulomb interactions. Furthermore, the lattice dynamical properties of these compounds remain unexplored. Here, we investigate the structural, electronic, and dynamical properties of Fe4(P2O7)3using density functional theory (DFT), with the local Coulomb and van der Waals (vdW) interactions included within the DFT+U method. The obtained lattice parameters show good agreement with the experimental data. Also, the magnetic configuration for the ground state corresponds well to the results of neutron measurements. The calculated arXiv:2503.15201v1 [cond-mat.mtrl-sci] 19 Mar 2025 2 phonon dispersion relations provide further insight into the dynamical properties of the material. II. CALCULATION METHOD The Vienna Ab initio Simulation Package (VASP) [15, 16] with the projector augmented-wave method [17] and the generalized gradient approximation [18] was used for relaxation of the crystal structure and to calculate the force constants. The geometry optimization was carried out within a unit cell, which is equivalent to the primitive cell, containing a total of 124 atoms. In order to determine the magnetic ground state, three different antiferromagnetic (AFM) arrangements were considered. To effectively sample the Brillouin zone, integration in the reciprocal space was performed over 4×2×4 k-point mesh generated with the Monkhorst-Pack scheme [19]. To ensure high precision of the calculations, a cut-off energy of 500 eV was applied. Computations were carried out using the DFT+U technique, incorporating distinct parameters for iron and oxygen atoms. The primary objective of these calculations was to gain insights into the electronic properties of iron pyrophosphate. The results presented here were obtained using Ud=8 eV and Jd=0.9 eV for the Fe(3d) states and Up=4 eV and Jp=0.5 eV for the O(2p) states. This parameter set resulted in an insulating state with an electronic band gap Eg=2.87 eV. For this specific combination of Uand Jparameters, we obtained a good agreement with the experimental lattice parameters, as shown in Tab. I. It should be noted that calculations performed with smaller Coulomb parameters (Ud= 4 eV, Up= 2 eV) as well as without the local interactions (Ud=Up=Jd=Jp= 0) did not converge to a welldefined ground state. This indicates the significance of local Coulomb interactions in iron pyrophosphate. To overcome the limitations of conventional approximations in accurately representing the long-range dispersive interactions during crystal structure optimization, we included approximate van der Waals (vdW) forces into the process. This involved incorporating vdW corrections proposed by Grimme, namely the D2 [20] and D3 [21] methods, as well as those suggested by Tkatchenko and Scheffler (TS) [22]. Our objective in incorporating these corrections was to refine the treatment of interatomic interactions and elevate the precision of crystal structure optimization. The mechanical stability of Fe4(P2O7)3was investigated by computing its elastic constants at T= 0 K using the AELAS code [23].The mechanical properties, such as Young’s modulus, shear modulus, elastic anisotropy, as well as Poisson’s and Pugh’s ratios, were also derived and analyzed. Phonon dispersion relations and phonon density of states (PDOS) were computed utilizing the temperature dependent effective potential (TDEP) approach [24, 25] and the Alamode software [26]. To obtain the forceconstant matrices, the atomic multidisplacements were generated by the high efficiency configuration space sampling (HECSS) [27], and followed by the calculation of the resulting Hellmann-Feynman forces acting on all atoms. Phonon energies and polarization vectors were subsequently determined through the exact diagonalization of the dynamical matrix. The calculations were performed for the 2 ×1×1 supercell containing 248 atoms. Furthermore, the LO-TO splitting was calculated using the static dielectric tensor and Born effective charges obtained within density functional perturbation theory [28]. III. CRYSTAL STRUCTURE AND MAGNETIC ORDER At room temperature, Fe4(P2O7)3crystallizes in the monoclinic structure with the P21/n (14) space group and its primitive unit cell contains four chemical formula units [13]. There are four different iron crystallographic sites surrounded by deformed octahedra. Six distinct diphosphate groups (P2O7) provide six terminal oxygen atoms, which act as their coordination ligands. As depicted in Fig. 1, two types of octahedra joined by three common oxygen atoms (sharing one face) create double octahedra. This arrangement causes significant shortening of distances between iron atoms in two different Fe2O9dimers, as reported in Ref. [13]: l1=2.967 ˚ A, and l2=3.048 ˚ A. These experimental values were very well reproduced with the DFT calculations, which yielded l1=3.005 ˚ A, and l2=3.048 ˚ A. The atomic configuration also includes three non-equivalent types of diphosphate groups. In each case, the coordination of iron atoms involves all six oxygen atoms, except for the bridging oxygen atom, which remains unengaged. It results from the specific geometric arrangement where a mirror plane, perpendicular to the lattice’s a-axis, includes both phosphorus atoms, the bridging oxygen atom, and one of the terminal oxygen atoms. It’s important to note that in this structure, two variations of P2O7anions exist with two different bridging bond angles φ1=179.81◦and φ2=148.20◦. In our investigation, we performed calculations for the ferromagnetic (FM) and three antiferromagnetic (AFM) states. The results indicate that the energy of the FM state is higher than that of the AFM states, with an energy difference of ∆E= 0.135 eV per formula unit. Consequently, we focused our analysis on the AFM configuration. For the three different AFM configurations (AFM1, AFM2, and AFM3) shown in Fig. 2, the cell parameters and magnetic moments on Fe atoms were computed and compared with the experimental data in Tab. I. The AFM1 pattern reflects the magnetic configuration observed via neutron diffraction measurements below the N´eel temperature TN= 50 K [13]. This arrangement shows reverse magnetic moments on two Fe3+ ions within the dimers (along the adirection), while nearest-neighbor 3 φ1 φ2 l1 l2 FIG. 1. Crystal structure of Fe4(P2O7)3. Iron, phosphorus, and oxygen ions are represented as pink, green, and blue balls, respectively. The lengths l1,l2and angles φ1,φ2denote the intra-dimer Fe-Fe distances and bridging bond P-O-P angles, respectively. The image was rendered using VESTA software [29]. iron atoms along the cdirection exhibit the ferromagnetic (FM) order. The antiparallel orientation of moments within the Fe dimers stems from a direct AFM interaction due to relatively short Fe-Fe distance and indirect superexchange interactions along the Fe-O-Fe path [13]. Interactions between dimers reveal a frustrated nature due to both FM and AFM superexchange pathways via phosphorus ions. AFM2 displays moments in opposite directions along the aand caxes, while AFM3 features moments arranged in antiparallel alignment along the c axis and parallel magnetization within the dimers. The computed magnetic moment (m= 4.6µB/Fe atom) and lattice parameters exhibit negligible sensitivity to the specific AFM configuration. Only the monoclinic angle βfor AFM3 slightly deviates from the experimental value (refer to Tab. I). This configuration also possesses the highest total energy due to the parallel alignment of magnetic moments within the Fe dimers. For each AFM configuration, the magnetic moments at distinct Fe crystallographic sites are slightly different, and their values are presented in Supplementary Material (SM). TABLE I. Lattice parameters and magnetic moments in iron pyrophosphate compared with the experimental data. AFM1 AFM2 AFM3 Experimental [13] a(˚ A) 7.441 7.441 7.438 7.389 b(˚ A) 21.571 21.569 21.571 21.337 c(˚ A) 9.583 9.583 9.585 9.517 β(deg) 90.07 90.08 91.27 90.03 ∆E(eV) 0 0.067 0.07 – m(µB) 4.6 4.6 4.6 4.55(5) AFM1 AFM2 AFM3 FIG. 2. Representation of magnetic moments on the iron atoms of Fe4(P2O7)3for three AFM configurations. The identified ground state with the AFM1 configuration closely corresponds to experimental outcomes, displaying a very good agreement with the measured mag- 4 netic moment of m= 4.55 µBat T= 1.7 K [13]. Additionally, inclusion of the vdW attractive forces reduces the lattice constants, significantly improving their values (refer to Tab. II). The three different types of vdW corrections give similar lattice parameters, and the best agreement with the experimental data was achieved for the D2 method. In the following investigations, we focus on exploring the electronic and phonon properties exclusively for the AFM1 ground state configuration. TABLE II. Lattice parameters of iron pyrophosphate computed for three types of vdW corrections and compared with the experimental data. D2 D3 TS Experimental [13] a(˚ A) 7.406 7.432 7.403 7.389 b(˚ A) 21.425 21.521 21.492 21.337 c(˚ A) 9.529 9.556 9.522 9.517 β(deg) 90.07 90.08 91.27 90.03 IV. ELECTRONIC PROPERTIES The electronic band structure of Fe4(P2O7)3shows the Mott insulating state with the band gap Eg= 2.87 eV (see Fig. 3), which separates a narrow, flat conducting band from a wide valence band spanning approximately 11 eV, and very weakly depends on the wave vector. The AFM1 ordering results in the same total electronic density of states (EDOS) for both spin components (refer to Fig. 4). An examination of the partial EDOS, which showcases the electronic states projected onto various iron atom types, reveals the presence of exchange Y C Z 10.0 7.5 5.0 2.5 0.0 2.5 5.0 7.5 E (eV) FIG. 3. Electronic band structure of iron pyrophosphate. The Fermi level (dotted line) was set to zero. 0 25 50 75 100 (a) Total 0 20 40 60 (b)Fe E (eV) 0.0 2.5 5.0 7.5 10.0 (c)P 20 15 10 5 0 5 10 E (eV) 0 10 20 30 (d)O EDOS (1/eV) FIG. 4. The total and partial electronic density of states of iron pyrophosphate, with the Fermi level (dotted line) set at zero. splitting within the 3dstates (see the EDOS for the four distinct iron atoms shown in Fig. [S2] in SM). Within the energy range spanning from -7 eV to the Fermi energy (EF), weak contributions from both spin components are evident. A main band, centered roughly at -8 eV, primarily consists of majority-spin states, constituting the lower Hubbard band. The unoccupied minority-spin 3d states displays a notable peak at around 4 eV, forming the upper Hubbard band. Within the energy range spanning from -12 eV to EF, the Fe(3d) states demonstrate hybridization with the O(2p) states, with the latter prevailing in the vicinity of EF. The O(2s) states are localized within the energy range between -22 and -17 eV. In the empty bands situated above the forbidden gap, a minor contribution stemming from oxygen sand pstates becomes discernible. Additionally, the 3pstates inherent to the phosphorus atoms remain occupied within the energy span from -11 eV to EF, accompanied by an extra peak discernible around -11 eV. The slender bands below -17 eV primarily accommodate the P(3s) states. 5 V. MECHANICAL STABILITY In order to confirm the intrinsic phase stability of the AFM1 phase, we determine the elastic constants (Cij) at T= 0 K using the AELAS code [23]. For the monoclinic and lower symmetry crystals, one usually inspects the eigenvalues of the elastic constants tensor or the corresponding compliance coefficients tensor. The crystal phase with the monoclinic symmetry is described by thirteen independent elastic constants, see Ref. [30] and references therein. The Cij values for the AFM1 phase are listed in Table III. TABLE III. The elastic constants (Cij ) in GPa for the AFM1 phase of Fe4(P2O7)3. C11 C12 C13 C15 C22 C23 C25 168.51 43.80 58.00 -6.84 162.87 53.47 -8.41 C33 C35 C44 C46 C55 C66 179.01 -5.93 52.18 -0.84 52.74 18.90 Some of the Cij are negative, however, this is allowed for the monoclinic symmetry. The important fact is that all eigenvalues of the elastic constants tensor are positive, and therefore this compound is intrinsically stable in monoclinic structure. Using the single crystal Cij values, we further elaborate on other mechanical characteristics of Fe4(P2O7)3. The polycrystaline elastic behavior was derived using the Hill averaging scheme [31] which combines the Reuss and Voigt[32] limits of bulk (BR,BV) and shear (GR,GV) moduli. This allows us to determine the averaged bulk, shear, and Young’s moduli: E=9BRV H GRV H 3BRV H +GRV H as well as the Poisson ratio: ν=3BRV H −ERV H 6BRV H ), and elastic anisotropy: ARV H =(GV−GR) (GV+GR) see the table IV. TABLE IV. The Hill averaged values of bulk (BRV H ), shear modulus (GRV H ), elastic anisotropy (ARV H ), the Pugh ratio (BRV H /(GRV H ), Young modulus (E), and the Poisson ratio (ν). The units for (BRV H ), (GRV H ), (ERV H ) is GPA, the ARV H , Pugh and Poisson ratios are dimensionless. BRV H GRV H ARV H BRV H /GRV H ERV H ν 90.5 44.2 0.096 2.05 113.9 0.29 The elastic anisotropy within the interval (0–1) indicates the deviation from the elastically isotropic behavior, i.e.,A= 0, which indeed is the case here, as A= 0.096 is much lower than similar low-dimensional frustrated antiferromagnets like KCuF3 and Cu(en)(H2O)2SO4, where A=0.26 and A=0.20, respectively [30]. To assess the ductility or brittleness of Fe4(P2O7)3we use the values of BRV H , GRV H ,ERV H ,and ν. In contrast to KCuF3and Cu(en)(H2O)2SO4, Fe4(P2O7)3is definitively not soft, as BRV H is rather moderate. The Pugh ratio (BRV H /GRV H ) is over factor 2 indicating ductility of the compound similar to the above-mentioned AFM compounds. Similarly, the ductility is correlated with the lower values of ERV H and GRV H as well as having the Poisson ratio (ν) below 0.5, indicating weaker and less directional bonding, i.e. more isotropic behavior [33]. The directional Young moduli are obtained from the compliance tensor elements as Ex=S−1 11 = 145 GPa, Ey=S−1 22 = 142 GPa, and Ez=S−1 33 = 149 GPa, again indicating very small elastic anisotropy despite large difference in length of the lattice vectors. Let us note, that these directional Young moduli are higher than the polycrystalline averaged one E due to the monoclinic structure. Inserting the density, ρ= 3.219g/cm−3from our calculated volume and the averaged BRV H ,GRV H one determines the longitudinal vl=pBRV H + (4/3)GRV H /ρ and transverse vt=pGRV H /ρ sound velocities in this compound[34]. The calculated values amount to vl= 2154.6 m.s−1and vt= 1171.8 m.s−1, respectively. VI. DYNAMICAL PROPERTIES In this section, we delve into the lattice dynamical properties of iron pyrophosphate. The primitive cell contains 124 atoms, which results in a 372 phonon dispersion relations. The phonon dispersion relations and density of states obtained for the primitive cell are presented in SM. In this case, most frequencies are real, except for four phonon modes near the Y and C points. Notably, the lowest soft mode in the primitive cell is identified at the Y point with the wave vector q= (0.5,0,0). Therefore, we doubled the crystal structure along the adirection to investigate the dynamical properties and stability of this material. Fig. 5(a, b) presents the phonon dispersion curves calculated for the 2 ×1×1 supercell obtained within the harmonic approximation and with the third-order anharmonic correction included in the atomic potential at T= 100 K. Comparing to the result obtained for the primitive cell (see SM), we observe the stabilization of the soft modes at the Y and C points, and only weak soft modes close to the Γ and B points are found in the harmonic calculations. However, the third-order anharmonic corrections remove the imaginary modes, showing the significant influence of anharmonicity on dynamical properties. Therefore, the inclusion of anharmonic corrections significantly impacts the overall stability and al- 6 lows for a more accurate representation of the atomic vibrations and interactions within the structure. This result agrees with the elastic data analysis, confirming the mechanical stability of the system. Fig. 5(c) displays the total phonon density of states (PDOS) along with the partial contributions from Fe, O, and P vibrations calculated with the third-order anharmonic correction. In the low-frequency range (0–400 cm−1), O vibrations predominantly govern phonon states, with comparatively smaller contributions from Fe and P atoms. Between 400-700 cm−1, oxygen atom vibrations dominate, while around 720 cm−1, phosphorous atom vibrations prevail. At higher frequencies (approximately 900-1250 cm−1), an elevated density of O atom vibrations is noticeable, with slightly lesser contributions from P atom vibrations. VII. SUMMARY In summary, the presented study of the structural, electronic, and dynamical properties of Fe4(P2O7)3has yielded new knowledge about this compound. The calculations carried out for three different antiferromagnetic orderings allowed to determine the ground state configuration, which closely match the experimental data. By using the DFT+U approach, the Mott insulating state with the electronic band gap Eg= 2.87 eV was predicted. The examination of the crystal structure demonstrated the existence of coordinated terminal oxygen atoms from diphosphate groups, resulting an interconnected double octahedron of iron ions. The accuracy of the computational approach was validated by the calculated lattice parameters, iron-iron separations, and bond angles, which showed good agreement with the experimental data. Additionally, the analysis revealed that Fe4(P2O7)3contains three different types of crystallographically non-equivalent diphosphate groups, each of which adds to the overall structural stability of the compound. The analysis of the obtained elastic parameters confirmed the mechanical stability of the crystal structure. The lattice dynamics of iron pyrophosphate was investigated by calculating phonon dispersion relations and phonon density of states using the harmonic and thirdorder anharmonic approximations within the TDEP approach. The analysis of the phonon dispersion curves identified soft modes; however, these were stabilized by taking into account anharmonicity of the lattice. VIII. ACKNOWLEDGEMENT D.L. acknowledges financial support by the projects No. 22-35410K from the Czech Science Foundation and No. CZ.02.01.01/00/22 008/0004572 by the Ministry of Education, Youth, and Sports of Czechia. [1] X. Liu, Y. Xu, R. Jin, P. Yin, L. Sun, T. Liang, and S. Gao, Facile synthesis of hierarchical Fe4(P2O7)3for removal of u(vi), Journal of Molecular Liquids 200, 311 (2014). [2] Q. Ma, Y. Sun, C. Zhang, Y. Xue, Y. Chen, W. Teng, and J. Fan, Iron pyrophosphate doped carbon nanocomposite for tetracycline degradation by activation of peroxymonosulfate, New J. Chem. 46, 17985 (2022). [3] C. Xia, Y. Li, M. Je, J. Kim, S. M. Cho, C. H. Choi, H. Choi, T. H. Kim, and J. K. Kim, Nanocrystalline iron pyrophosphate-regulated amorphous phosphate overlayer for enhancing solar water oxidation, Nano-Micro Letters 14, 2150 (2022). [4] B. Boonchom and N. Phuvongpha, Synthesis of new binary cobalt iron pyrophosphate CoFeP2O7, Materials Letters 63, 1709 (2009). [5] A. Gupta, N. B. Amin, A. Besarab, S. E. Vogel, G. W. Divine, J. Yee, and J. V. Anandan, Dialysate iron therapy: Infusion of soluble ferric pyrophosphate via the dialysate during hemodialysis, Kidney International 55, 1891 (1999). [6] M. C. Fidler, T. Walczyk, L. Davidsson, C. Zeder, N. Sakaguchi, L. R. Juneja, and R. F. Hurrell, A micronised, dispersible ferric pyrophosphate with high relative bioavailability in man, British Journal of Nutrition 91, 107–112 (2004). [7] R. Trivedi and K. Barve, Delivery systems for improving iron uptake in anemia, International Journal of Pharmaceutics 601, 120590 (2021). [8] G.-H. Lee, Y. S. Kim, and D.-W. Kim, Redox effect of Fe2+/Fe3+ in iron phosphates for enhanced electrocatalytic activity in li-o2 batteries, Chemical Engineering Journal 388, 124294 (2020). [9] S. Nishimura, M. Nakamura, R. Natsui, and A. Yamada, New lithium iron pyrophosphate as 3.5 v class cathode material for lithium ion battery, Journal of the American Chemical Society 132, 13596 (2010), https://doi.org/10.1021/ja106297a. [10] J. Hoggins, J. Swinnea, and H. Steinfink, Crystal structure of Fe2P2O7, Journal of Solid State Chemistry 47, 278 (1983). [11] Z. W. Xiao, G. R. Hu, Z. D. Peng, K. Du, and X. G. Gao, Solid state synthesis and characterization of iron (II) pyrophosphate fe2p2o7, Chinese Chemical Letters 18, 1525 (2007). [12] A. K. Padhi, K. S. Nanjundaswamy, C. Masquelier, S. Okada, and J. B. Goodenough, Effect of structure on the Fe3+/Fe2+ redox couple in iron phosphates, Journal of the Electrochemical Society 144, 1609 (1997). [13] L. K. Elbouaanani, B. Malaman, R. G´erardin, and M. Ijjaali, Crystal structure refinement and magnetic properties of Fe4(P2O7)3studied by neutron diffraction and m¨ossbauer techniques, Journal of Solid State Chemistry 163, 412 (2002). [14] Q. Shi, L. Zhang, M. E. Schlesinger, J. Boerio-Goates, and B. F. Woodfield, Low temperature heat capacity 7 Y C Z B D Wave vector 0 200 400 600 800 1000 1200 Frequency (cm 1) (a) Y C Z B D Wave vector (b) (c) Total Fe P O PDOS (arb. units) Harmonic Anharmonic FIG. 5. The phonon dispersion relations calculated using a) harmonic approximation and b) third-order anharmonic correction, and c) the total and partial element-projected PDOS for anharmonic calculation in Fe4(P2O7)3. study of Fe3PO7and Fe4(P2O7)3, The Journal of Chemical Thermodynamics 62, 86 (2013). [15] G. Kresse and J. Furthm¨uller, Software vasp, Vienna University of Technology, Vienna (1999). [16] G. Kresse and J. Furthm¨uller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comp. Mat. Sci. 6, 15 (1996). [17] P. E. Bl¨ochl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994). [18] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996). [19] H. J. Monkhorst and J. D. Pack, Special points for brillouin-zone integrations, Phys. Rev. B 13, 5188 (1976). [20] S. Grimme, Semiempirical GGA-type density functional constructed with a long-range dispersion correction, J. Comp. Chem. 27, 1787 (2006). [21] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (dft-d) for the 94 elements H-Pu, The Journal of Chemical Physics 132, 154104 (2010). [22] A. Tkatchenko and M. Scheffler, Accurate molecular van der waals interactions from ground-state electron density and free-atom reference data, Phys. Rev. Lett. 102, 073005 (2009). [23] S. H. Zhang and R. F. Zhang, Aelas: Automatic elastic property derivations via high-throughput first-principles computation, Computer Physics Communications 220, 403 (2017). [24] O. Hellman, I. A. Abrikosov, and S. I. Simak, Lattice dynamics of anharmonic solids from first principles, Phys. Rev. B 84, 180301 (2011). [25] O. Hellman, P. Steneteg, I. A. Abrikosov, and S. I. Simak, Temperature dependent effective potential method for accurate free energy calculations of solids, Phys. Rev. B 87, 104111 (2013). [26] T. Tadano, Y. Gohda, and S. Tsuneyuki, Anharmonic force constants extracted from first-principles molecular dynamics: applications to heat transfer simulations, J. Phys.: Condens. Matter 26, 225402 (2014). [27] P. T. Jochym and J. La˙zewski, High Efficiency Configuration Space Sampling – probing the distribution of available states, SciPost Phys. 10, 129 (2021). [28] M. Gajdoˇs, K. Hummer, G. Kresse, J. Furthm¨uller, and F. Bechstedt, Linear optical properties in the projectoraugmented wave methodology, Phys. Rev. B 73, 045112 (2006). [29] K. Momma and F. Izumi, Vesta 3 for three-dimensional visualization of crystal, volumetric and morphology data, Journal of applied crystallography 44, 1272 (2011). [30] D. Legut, R. S`ykora, U. D. Wdowik, and A. Orend´aˇecov´a, Mechanical properties of the quasi-one-dimensional antiferromagnet Cu (en)(H2O) 2SO4, Journal of Nanoscience and Nanotechnology 19, 3016 (2019). [31] R. Hill, The elastic behaviour of a crystalline aggregate, Proceedings of the Physical Society. Section A 65, 349 (1952). [32] G. G¨oran, Thermophysical properties of materials (Elsevier, 1999). [33] F. Chu, T. E. Mitchell, B. Majumdar, D. Miracle, T. K. Nandy, and D. Banerjee, Elastic properties of the O phase in Ti-Al-Nb alloys, Intermetallics 5, 147 (1997). [34] O. L. Anderson, A simplified method for calculating the debye temperature from elastic constants, Journal of Physics and Chemistry of Solids 24, 909 (1963).