Stability of FeVO₄-II under Pressure: A First-Principles Study

Abstract

In this work, we report first-principles calculations to study FeVO${}_4$ in the CrVO${}_4$-type (phase II) structure under pressure. Total-energy calculations were performed in order to analyze the structural parameters, the electronic, elastic, mechanical, and vibrational properties of FeVO${}_4$-II up to 9.6 GPa for the first time. We found a good agreement in the structural parameters with the experimental results available in the literature. The electronic structure analysis was complemented with results obtained from the Laplacian of the charge density at the bond critical points within the Quantum Theory of Atoms in Molecules methodology. Our findings from the elastic, mechanic, and vibrational properties were correlated to determine the elastic and dynamic stability of FeVO${}_4$-II under pressure. Calculations suggest that beyond the maximum pressure covered by our study, this phase could undergo a phase transition to a wolframite-type structure, such as in CrVO${}_4$ and InVO${}_4$.

1. Introduction

High-pressure studies on $AB$O${}_4$ compounds have increased in the last two decades due to the evolution of the diamond-anvil cell (DAC) and its impact on other relevant areas, such as materials sciences [1,2]. Among $AB$O${}_4$ compounds, those with a CrVO${}_4$-type structure (Space Group, SG, $Cmcm$, No. 63, Z = 4, with [6 − 4] coordination for the [A-B cations]) have received less attention. CrVO${}_4$-type compounds are a family formed mainly by phosphates, vanadates, chromates, and sulfates [3]. Given their position in the Bastide diagram [4], it has been reported that they could have a structural phase transition driven by pressure to the wolframite [5,6] or zircon [7,8] structure. This implies a change in the coordination to [6 − 4 + 2] and [8 − 4] in the high-pressure phases, respectively.

According to the literature, AVO${}_4$ vanadates with CrVO${}_4$-type structures have been used in applications as photocatalyst materials [1,9,10], cathodoluminescent materials [11,12], and scintillators [13,14], among others. The properties of these compounds are reported to be highly dependent on the occupation of the d valence electrons of the transition metal cations [9,10,11,12,13,14]. Until now, only four vanadates with a CrVO${}_4$-type structure have been reported: CrVO${}_4$ [6], FeVO${}_4$ [15], InVO${}_4$ [16], and the less studied TlVO${}_4$ [3]. Unlike the other vanadates, FeVO${}_4$ presents a great polymorphism, crystallizing at ambient conditions in the triclinic AlVO${}_4$-type [17] structure (SG: $P\overline{1}$, No. 2, Z = 6) [15], while the CrVO${}_4$-type phase was obtained from high-pressure-high-temperature studies [18] and synthesized at ambient conditions by hydrothermal method [19]. Further, the antiferromagnetic order in FeVO${}_4$-II was observed after the characterization of its magnetic properties [19]. However, until now, no more experimental or theoretical studies have been reported to characterize this phase. It is worth mentioning that most of the studies carried out at high pressure on this compound have focused on characterizing the crystal structure of the different polymorphs [15,18,20,21,22].

In this work, we performed a first-principles characterization of FeVO${}_4$-II up to 9.6 GPa. The main goal is to determine the elastic and dynamic stability of the studied polymorph at ambient conditions and under pressure. In this sense, a complete report of the elastic constants, mechanical properties, elastic anisotropy, and vibrational properties is presented. Furthermore, given the relevance of the electronic structure of this compound in materials science, we computed the density of states and the band structure together with an analysis of the charge density with the QTAIM methodology [23,24].

2. Computational Details

Total energy calculations were performed within the framework of the density functional theory (DFT) [25] and the projector-augmented wave (PAW) [26,27] method as implemented in the Vienna Ab initio Simulation Package (VASP) [28,29,30,31]. We used a plane-wave energy cutoff of 520 eV to ensure high precision in our calculations. The exchange-correlation energy was described within the generalized gradient approximation (GGA) in the Perdew–Burke–Ernzerhof for solids (PBEsol) formulation [32]. The GGA+U was used to account for the strong correlation between the electrons in the d shell, on the basis of Dudarev’s method [33]. In this method, the Coulomb interaction U and the onsite exchange interaction $J_{\mathrm{H}}$ are treated together as $U_{\mathrm{eff}}$ = U - $J_{\mathrm{H}}$. For our GGA+U calculations, we choose U = 6 eV and $J_{\mathrm{H}}$ = 0.95 eV. Similar values were previously used with relative success in the study of other iron and vanadate compounds [15,20,34,35,36,37].

The Monkhorst-Pack scheme was employed for the Brillouin-zone (BZ) integrations [38] with a mesh 4 × 3 × 3, which corresponds to a set of 8 special k-points in the irreducible BZ. In the relaxed equilibrium configuration, the forces are less than 0.3 meV/Å per atom in each Cartesian direction. The highly converged results on forces are required for the calculations of the dynamical matrix using the direct force constant approach [39], which allows us to identify the irreducible representation and the character of the phonon modes at the zone center ($\mathsf{\Gamma }$ point). The phonon density of states (PDOS) were obtained from the calculation of the phonons in the whole BZ with a supercell 2 × 2 × 2 times the conventional unit cell by using the PHONON software. The calculations of elastic constants were made following the methodology validated in Ref. [40] with a k-points mesh of 8 × 6 × 6, a plane-wave energy cutoff of 570 eV, and a POTIM parameter of 0.016. The elastic tensor is determined by performing six finite distortions of the lattice and deriving the elastic constants from the strain–stress relationship [41]. The electron density was computed with a refined grid to be analyzed within the fundaments of Quantum Theory of Atoms in Molecules (QTAIM) [42] with CRITIC [43,44], and AIM [45] codes.

3. Results and Discussion

3.1. Crystal Structure

According to the literature, FeVO${}_4$ has been successfully synthesized as a metastable polymorph at ambient conditions in the orthorhombic CrVO${}_4$-type structure (SG: $Cmcm$, No. 63, Z = 4) [19], also called the FeVO${}_4$-II phase, see Figure 1a. In this structure, there are four non-equivalent Wyckoff positions (WP): Fe occupies the 4a (0, 0, 0), V the 4c (0, y, 1/4), and O the 8f (0, y, z), and 8g (x, y, 1/4). Internally, the structure consists of edge-sharing FeO${}_6$ distorted octahedra along the c direction. The –FeO${}_6$–FeO${}_6$– chains are linked together by means of the VO${}_4$ distorted tetrahedra, while the VO${}_4$ tetrahedra are not linked to each other, as can be seen in Figure 1b. As shown in Figure 1c, there are two non-equivalent interatomic distances in the FeO${}_6$ octahedra, the apical Fe-O distance ($d_{\mathrm{Fe}-\mathrm{O}1}$) and the equatorial Fe-O distance ($d_{\mathrm{Fe}-\mathrm{O}2}$). In the VO${}_4$ tetrahedra, there are two V–O interatomic distances ($d_{\mathrm{V}-\mathrm{O}1}$ and $d_{\mathrm{V}-\mathrm{O}2}$). In both polyhedra, the respective interatomic distance $d_{X-\mathrm{O}2}$ is larger than $d_{X-\mathrm{O}1}$ (X = Fe, V). Inside the polyhedrons, there are important angles to mention that are relevant for the structural changes in the polyhedra under pressure. The $\alpha$ angle is formed between the equatorial plane of the FeO${}_6$ octahedra and the polar bond, while the $\beta$ ($\gamma$) angle is between the two $d_1$ ($d_2$) interatomic distances of the VO${}_4$ tetrahedra.

The equilibrium volume and unit-cell parameters at $P\approx 0$ GPa were computed by minimizing the total crystal energy for different volumes ranging from −4 to 10 GPa while allowing the internal atomic positions and the lattice parameters to relax. We used the third-order Birch–Murnaghan [47] equation-of-state (EOS) to fit the volume–energy curve. Our results showed that the lowest energy configuration is obtained when an antiferromagnetic coupling between Fe ions is obtained, in accordance with the experimental results from Mössbauer effect [19]. The optimized lattice parameters (a, b, c), volume (V), interatomic distances ($d_{\mathrm{Fe}-\mathrm{O}1}$, $d_{\mathrm{Fe}-\mathrm{O}2}$, $d_{\mathrm{V}-\mathrm{O}1}$, $d_{\mathrm{V}-\mathrm{O}2}$), internal angles ($\alpha$, $\beta$, $\gamma$), bulk modulus (B) and pressure derivative of B (B’) appear in

Table 1. Calculated structural parameters for FeVO${}_4$-II phase. a, b, and c are the lattice parameters, V is the volume, $d_{X-\mathrm{O}n}$ (X = Fe, V; n = 1, 2) and $\alpha$, $\beta$, and $\gamma$, are the interatomic distances and angles inside the crystal structure, respectively (see Figure 1). B is the bulk modulus and $B^{\prime }$ is the pressure derivative of B. Results from this work are denoted with 🟉. We also include experimental results from the literature of AVO${}_4$ compounds for comparison.

	a (Å)	b (Å)	c (Å)	V (Å${}^3$)	${\mathit{d}}_{\mathbf{Fe}-\mathbf{O}1}$ (Å)	${\mathit{d}}_{\mathbf{Fe}-\mathbf{O}2}$ (Å)	${\mathit{d}}_{\mathbf{V}-\mathbf{O}1}$ (Å)	${\mathit{d}}_{\mathbf{V}-\mathbf{O}2}$ (Å)	$\mathit{\alpha }$ (${}^{\circ }$)	$\mathit{\beta }$ (${}^{\circ }$)	$\mathit{\gamma }$ (${}^{\circ }$)	B (GPa)	${\mathit{B}}^{\prime }$	Ref.
FeVO${}_4$	5.6038	8.2624	6.1726	285.80	1.9734	2.0377	1.6705	1.7842	89.271	106.990	112.547	88.8	2.8	🟉
	5.5941	8.3216	6.2252	289.80	1.9784	2.0521	1.6544	1.7762				88.3	3.2	DFT [15]
	5.6284	8.2724	6.1118	284.57	2.0060	2.0290	1.6520	1.7920	89.600	107.200	112.500			Exp. [19]
CrVO${}_4$	5.5785	8.2830	6.0576	279.90								93.0	4.0	DFT [6]
	5.5680	8.2080	5.9770	273.16								63.0	4.0	Exp. [6]
InVO${}_4$	5.7547	8.6168	6.6751	331.00	2.1877	2.1392	1.6690	1.7824				71.0	4.0	DFT [5]
	5.7380	8.4920	6.5820	320.72								69.0	4.0	Exp. [16]
TlVO${}_4$	5.8390	8.6870	6.8000	344.92										Exp. [48]

together with experimental results from the literature [19]. Furthermore, we have included the reported parameters for other CrVO${}_4$-type structure vanadates for comparison. According to Table 1, the equilibrium volume at ambient pressure of AVO${}_4$ (A = Cr, Fe, In, Tl) compounds increases with the atomic mass of cation A, while the opposite occurs for bulk modulus. The obtained results for FeVO${}_4$–II are in excellent agreement with experiments, with a difference of less than 1% in the lattice parameters. For this volume, the WPs are 4a (0, 0, 0) for Fe, 4c (0, 0.35788, 1/4) for V, 8f (0, 0.76239, 0.96753), and 8g (0.26481, 0.47778, 1/4) for O. Following Table 1, the equilibrium volume at ambient pressure of AVO${}_4$ (A = Cr, Fe, In, Tl) compounds increases with the atomic mass of cation A. In contrast, the opposite occurs for bulk modulus.

According to the most recent high-pressure experiments carried out on FeVO${}_4$ [15,20], this compound undergoes the following phase transition sequence, $P\overline{1}\to P{\overline{1}}^{\prime }\to$$\alpha$-MnMoO${}_4$-type, driven by pressure. However, previous high-pressure studies reported the next phase transition sequence $P\overline{1}\to$ CrVO${}_4$-type → wolframite in a small range of pressure [49]. Therefore, it is unknown if the CrVO${}_4$-type structure of FeVO${}_4$ would transition to wolframite when high-pressure experiments are carried out starting from the CrVO${}_4$-type structure, as occurs with CrVO${}_4$ [6] and InVO${}_4$ [5,16], and was predicted through first-principles calculations [15]. Therefore, from now on, we will only deal with studying the behavior of FeVO${}_4$-II under pressure without considering whether the CrVO${}_4\to$ wolframite transition occurs. Figure 2 shows the pressure evolution of the structural parameters for the FeVO${}_4$-II phase. According to Figure 2a, a volume reduction of $\mathsf{\Delta }V$ = −9% is observed for a range of pressure of ≈10 GPa. Considering the anisotropy of the orthorhombic crystal structure, $\mathsf{\Delta }V$ is more reflected in the [010] and [001] directions, while the change in [100] is practically negligible, see Figure 2b. Indeed the compressibility of the [001] direction is comparable to the linear compressibility of diamond, making the studied material a near zero linear compressibility material along [001]. The anisotropy can also be understood with the help of the other internal parameters from Figure 2c–f. As seen in Figure 2c, the interatomic distances $d_{\mathrm{Fe}-\mathrm{O}n}$ of FeO${}_6$ undergo more changes than the $d_{\mathrm{V}-\mathrm{O}n}$ distances of VO${}_4$, which is reflected in the larger compressibility of $V_{{\mathrm{FeO}}_4}$ octahedra against the minimal change in $V_{{\mathrm{VO}}_4}$, see Figure 2d. Hence, the large compressibility of the y-axis is due primarily to the compressibility of octahedra due to the reduction in $d_{\mathrm{Fe}-\mathrm{O}2}$. Whereas, despite the reduction in $d_{\mathrm{V}-\mathrm{O}2}$, the $\beta$ and $\gamma$ angles adjust in such a way that $V_{{\mathrm{VO}}_4}$ remains almost constant with pressure. On the other hand, the distortion parameter ${\mathsf{\Delta }}_d$ presents an important change at 6.2 GPa, which will be discussed in detail in Section 3.3.

3.2. Electronic Structure

Figure 3a,b displays our calculated density of states (DOS) and band structures along the high-symmetry path Z–$\mathsf{\Gamma }$–S–R–Z–T–Y–$\mathsf{\Gamma }$ of the FeVO${}_4$-II phase at ambient pressure and 9.6 GPa, respectively. Red dashed lines at 0 eV indicate the Fermi level. At ambient pressure, the valence band maximum (VBM) is located at an H ($\frac{1}{2}$, $\frac{1}{2}$, $-u$) point (between T and Y), mainly occupied by 2p ($p_x$) states from oxygen. In contrast, the conduction band minima (CBM) is located at the $\mathsf{\Gamma }$ point with contributions of the 3d states from Fe, mostly from $d_{yz}$ sub-orbital. Therefore, FeVO${}_4$-II phase presents an indirect energy band-gap $E_g$ = 2.11 eV. In comparison, the direct $\mathsf{\Gamma }$–$\mathsf{\Gamma }$ energy band-gap is 2.17 eV. In this case, the maximum at $\mathsf{\Gamma }$ below the Fermi level is occupied by 2p ($p_x$ and $p_y$) states from oxygen, as happens in other $AB$O${}_4$ compounds [40]. It has been reported that CBM of many rare-earth AVO${}_4$ vanadates is dominated by vanadium 3d states, which makes these compounds present energy band-gaps close to 3.8 eV [50,51]. However, in the case of vanadates InVO${}_4$ (phase II), CrVO${}_4$ (phase II), and FeVO${}_4$ (phase I), the energy band-gap lowers to 3.2, 2.6, and 2.1 eV, respectively, due to contribution from cation A to the valence and conduction band [5,6,50,52]. Hence, our $E_g$ value follows the observed results for those vanadates. According to partial electronic DOS of Figure 3a, the states below the Fermi level mainly belong to the O atoms. In contrast, the states above the conduction bands are due to Fe, followed by V and O.

As pressure increases up to 9.6 GPa, the VBM shifts to the high-symmetry point Y while the CBM is still at the $\mathsf{\Gamma }$ point with the same band occupation observed at ambient pressure, see Figure 3b. Hence, at this pressure, the band-gap is still indirect with $E_g$ = 2.05 eV, which gives a pressure coefficient of $d{E}_g/dp$ = −7.1 meV/GPa. In contrast, the pressure coefficients for InVO${}_4$ and CrVO${}_4$ are positive, with values of 13.5 [52] (8.9 [5]) meV for InVO${}_4$ and 1.9 (4.1) meV for CrVO${}_4$ [6], obtained by optical absorption measurements (DFT calculations). On the other hand, the forbidden direct $\mathsf{\Gamma }$–$\mathsf{\Gamma }$ band-gap increases with pressure up to 2.21 eV at 9.6 GPa, which gives a value of $d{E}_g/dP$ = 4.2 meV, similar to the obtained values for the other AVO${}_4$ vanadates.

We performed a topological analysis of the charge density ($\rho$) with the QTAIM methodology to complement the study of the electronic structure of FeVO${}_4$-II. Figure 4a,b show the gradient vector field of the electron density ($\nabla \rho$) of FeO${}_6$ and VO${}_4$ polyhedra, respectively (gold for Fe, purple for V, and red for O). Blue points indicate the atomic nucleus, corresponding to the maximum critical points (MCP). The minimal CP (mCP) are in green, and the bond CPs (BCPs) are in black and are located in the bond path (blue dashed lines) between Fe (V) and O. The distance between Fe (V) and the BCP ($d_{\mathrm{BCP}}$) is plotted as a function of pressure in Figure 4c. As seen, the distances $d_{\mathrm{BCP}}$ follow the same trend observed in the interatomic distances of Figure 2c. Therefore, the crystal structure can be interpreted in terms of charge density. In this sense, the Laplacian of the charge density (${\nabla }^2{\rho }_{\mathrm{BCP}}$) at the BCP gives an idea of the ionicity (${\nabla }^2{\rho }_{\mathrm{BCP}}>0$) or covalency (${\nabla }^2{\rho }_{\mathrm{BCP}}<0$) of a system [53]. According to Figure 4d, the Fe–O and V–O bonds have an ionic nature. The obtained values of ${\nabla }^2{\rho }_{\mathrm{BCP}}$ (0.214–0.295 a.u.) for Fe–O bonds in the distorted octahedra are larger than those reported for A–O bonds in AWO${}_4$ (0.086–0.147 a.u.) [54], AMoO${}_4$ (0.085–0.22 a.u.) [53], and ATcO${}_4$ (0.045–0.153 a.u.) [40] compounds. In the case of the VO${}_4$ distorted tetrahedra, there is a significative difference in the values of ${\nabla }^2{\rho }_{{\mathrm{BCP}}_1}$ and ${\nabla }^2{\rho }_{{\mathrm{BCP}}_2}$, which is related to the interatomic distances $d_{\mathrm{V}-\mathrm{O}n}$ (see Table 1) being significantly larger than those obtained for the FeO${}_6$ distorted octahedra. This is the reason for the low compressibility of the VO${}_4$ polyhedra. As the pressure grows, the Laplacian values increase, strengthening the Fe–O and V–O bonds and the resistance of the polyhedrons to compression.

3.3. Elastic Constants and Mechanical Properties

The elastic constants ($c_{ij}$) determine the resistance to deformation in the directions and planes of a crystal, making them highly dependent on the system’s symmetry. For the orthorhombic crystalline system, there are nine independent elastic constants: $c_{11}$, $c_{22}$ and $c_{33}$ for longitudinal compression, $c_{12}$, $c_{13}$ and $c_{23}$ for transverse expansion, $c_{44}$, $c_{55}$, and $c_{66}$ for pure shear [55].

Table 2. Elastic constants $c_{ij}$ (in GPa) at 0 GPa and pressure coefficients.

	${\mathit{c}}_{11}$	${\mathit{c}}_{22}$	${\mathit{c}}_{33}$	${\mathit{c}}_{44}$	${\mathit{c}}_{55}$	${\mathit{c}}_{66}$	${\mathit{c}}_{12}$	${\mathit{c}}_{13}$	${\mathit{c}}_{23}$
$c_{ij}$	234.2	125.5	144.8	62.98	52.49	47.04	78.46	68.21	39.81
$d{c}_{ij}/dP$	2.41	−0.87	3.78	0.74	1.04	0.15	2.39	2.18	2.57

shows the elastic constants $c_{ij}$ of FeVO${}_4$-II at 0 GPa and the respective pressure coefficients $d{c}_{ij}/dP$. The necessary and sufficient Born stability criteria for orthorhombic systems are [56]: $c_{11}>0$, $c_{11}{c}_{22}>c_{12}^2$, $c_{11}{c}_{22}{c}_{33}+2c_{12}{c}_{13}{c}_{33}+2c_{12}{c}_{13}{c}_{23}-c_{11}{c}_{23}^2-c_{22}{c}_{13}^2-c_{33}{c}_{12}^2>0$, $c_{44}>0$, $c_{55}>0$, and $c_{66}>0$, which are fulfilled for FeVO${}_4$-II. Under hydrostatic pressure (P), the stability criteria are [57]: $c_{ii}-P>0$ (i = 1 to 6), $c_{11}+c_{22}+c_{33}+2c_{12}+2c_{13}+2c_{23}+3P>0$, $c_{11}+c_{22}-2c_{12}-4P>0$, $c_{11}+c_{33}-2c_{13}-4P>0$, and $c_{22}+c_{33}-2c_{23}-4P>0$, which are also meet. According to Table 2, the elastic constant $c_{11}$ presents the higher value, followed by $c_{33}$ and $c_{22}$, indicating that longitudinally the crystal is more resistant to deformation in the a-axis; this was also observed in the evolution of the lattice parameters under pressure in Figure 2b. As would be expected, the largest transverse expansion constant is $c_{12}$, followed by $c_{13}$ and $c_{23}$, $c_{33}$ being the smallest of the nine elastic constants. At the same time, the shear elastic constants are generally small against those from longitudinal compression and transverse expansion. We found that elastic constants follow a linear trend with pressure, where only $c_{22}$ presents a negative pressure coefficient related to the trends observed in the structural parameters with pressure.

The elastic constants and compliances were used to compute the bulk (B) and shear (G) modulus with the Hill procedure [58], which is the average of the Voigt (V) [59] and Reuss (R) [60] methods. The bulk modulus in the Voigt and Reuss methods is computed as follows: $9B_{\mathrm{V}}$ = ($c_{11}$ + $c_{22}$ + $c_{33}$ + 2$c_{12}$ + 2$c_{13}$ + 2$c_{23}$), and $B_{\mathrm{R}}^{-1}$ = ($s_{11}$ + $s_{22}$ + $s_{33}$ + 2$s_{12}$ + 2$s_{13}$ + 2$s_{23}$), respectively. While the shear modulus is obtained with the following equations: 15$G_{\mathrm{V}}$ = ($c_{11}$ + $c_{22}$ + $c_{33}$ − $c_{12}$ − $c_{23}$ − $c_{13}$ + 3$c_{44}$ + 3$c_{55}$ + 3$c_{66}$), and $15G_{\mathrm{R}}^{-1}$ = [4($s_{11}$ + $s_{22}$ + $s_{33}$) − 4($s_{12}$ + $s_{13}$ + $s_{23}$) + 3($s_{44}$ + $s_{55}$ + $s_{66}$)], respectively. B and G in the Hill average are obtained by means of $B=(B_{\mathrm{V}}+B_{\mathrm{R}})/2$ and $G=(G_{\mathrm{V}}+G_{\mathrm{R}})/2$, respectively. On the other hand, Young’s modulus (E) and Poisson’s ratio ($\nu$) are given by $E=9BG/(3B+G)$ and $\nu =(3B-2G)/(6B+2G)$, respectively. The Vickers hardness was calculated with the approximation of Tian et al. [61]: $H_v$ = 0.92 $k^{1.137}$$G^{0.708}$, $k=G/B$.

The mechanical properties and their evolution with pressure are provided in Figure 5. As seen, B increases almost linearly with pressure, which means that the resistance to deformation increases for this pressure range. Whereas G undergoes a change in the linear trend above 3 GPa, reaching a maximum value of 53.9 at 6.2 GPa and falling down to 53.3 GPa at 9.6 GPa. There is no significant change in G for this pressure range, but the change in trend will substantially affect the other properties that depend on G. Interestingly, in other materials, such change in the slope of the pressure dependence of the shear modulus G has been identified as a precursor effect of a phase transition observed at higher pressures [62]. According to the equations of Voigt and Reuss for G, the behavior of G is due to the negative pressure coefficient for $c_{22}$ and the significant values of $d{c}_{ij}/dP$ for the transverse expansion elastic constants. On the other hand, B behaves linearly because $B_{\mathrm{V}}$ ($B_{\mathrm{G}}$) only sums the longitudinal and transverse elastic constants (compliances), while $c_{22}$ ($s_{22}$) decreases with pressure, so the value of $c_{22}$ affects B by a small constant factor.

The behavior of Young’s modulus is similar to G due to its relation with this mechanical property. According to Figure 5c, E increases by 4.3 % in a pressure range of 7.86 GPa; above this pressure, it decreases from 138.1 to 137.6 GPa. This behavior indicates that FeVO${}_4$ could undergo a phase transition from a CrVO${}_4$-type phase to another one at a relatively small pressure since it has been observed that FeVO${}_4$ undergoes several phase transitions below 10 GPa, see Refs. [15,20] and references therein. The Poisson’s ratio follows an increasing trend, taking values between 0.26 and 0.29, which are in the range for ceramic semiconductors (0.25–0.42) [63] and are similar for other $AB$O${}_4$ compounds [40,54,64]. These values indicate that ionic contribution dominates in FeVO${}_4$-II. The Pugh relation, Figure 5e, represents the ductility or brittleness of the system. A value below 1.75 [11] reveals a tendency to brittleness, while a value greater than this indicates that the material behaves in a ductile manner. FeVO${}_4$-II presents a value of 1.745 at ambient pressure and increases to 2.06 at 9.6 GPa, which means that ductility increases with pressure. As was expected, FeVO${}_4$-II, as other $AB$O${}_4$ compounds, presents a small Vickers hardness due to the relatively small values of B and G. Given the relation between $H_{\mathrm{V}}$ and $B/G$, it is natural that the Vickers hardness decreases with pressure, Figure 5f.

On the other hand, elastic anisotropy has important implications in engineering science and crystal physics because it is commonly responsible for the formation of microcracks in materials [65]. In this work, we analyzed several anisotropic factors discussed below. The anisotropic shear factors provide a measure of the degree of anisotropy in the bonding between atoms in different crystallographic planes. The shear anisotropic factor for the {100} shear planes between the $\langle 011\rangle$ and $\langle 010\rangle$ directions is $A_{\left\{100\right\}}=4c_{44}/(c_{11}+c_{33}-2c_{13})$. For the {010} shear planes between $\langle 101\rangle$ and $\langle 001\rangle$ directions it is $A_{\left\{010\right\}}=4c_{55}/(c_{22}+c_{33}-2c_{23})$. For the {001} shear planes between $\langle 110\rangle$ and $\langle 010\rangle$ directions, it is $A_{\left\{001\right\}}=4c_{66}/(c_{11}+c_{22}-2c_{12})$. The factors $A_{\left\{100\right\}}$, $A_{\left\{010\right\}}$, and $A_{\left\{001\right\}}$ must be one for isotropic crystals, and any value greater or smaller than one measures the degree of elastic anisotropy. According to Figure 6a, $A_{\left\{100\right\}}$ presents the smallest deviation from one, followed by $A_{\left\{010\right\}}$ and $A_{\left\{001\right\}}$. However, once pressure is applied, $A_{\left\{010\right\}}$ is more affected than the others. Factor $A_{\left\{010\right\}}$ is related to the equatorial interatomic distances of the distorted octahedra [$d_{\mathrm{Fe}-\mathrm{O}2}$, see Figure 2c], which makes quite a bit of sense since it is $d_{\mathrm{Fe}-\mathrm{O}2}$ that undergoes the greatest change with pressure. While $A_{\left\{001\right\}}$ is more related to $d_{\mathrm{Fe}-\mathrm{O}1}$, with this interatomic distance being the second one most affected by pressure.

Other important factors are the universal anisotropy $A_U=5G_{\mathrm{V}}/G_{\mathrm{R}}+B_{\mathrm{V}}/B_{\mathrm{R}}-6$, the percentage of elastic anisotropy in bulk [$A_B=(B_{\mathrm{V}}-B_{\mathrm{R}})/(B_{\mathrm{V}}+B_{\mathrm{R}})$] and shear [$A_G=(G_{\mathrm{V}}-G_{\mathrm{R}})/(G_{\mathrm{V}}+G_{\mathrm{R}})$] modulus, where a value of zero is associated with isotropic elastic constants, which means that $B_{\mathrm{V}}=B_{\mathrm{R}}$ and $G_{\mathrm{V}}=G_{\mathrm{R}}$. For the entire range of pressure, $A_B$ is larger than $A_G$ because the difference between $B_{\mathrm{V}}$ and $B_{\mathrm{R}}$ is bigger than between $G_{\mathrm{V}}$ and $G_{\mathrm{R}}$. However, the pressure affects more $A_G$ than $A_B$ as happens with G. This difference increases above 6.2 GPa as $G_{\mathrm{R}}$ presents the same behavior of G, while $G_{\mathrm{V}}$ is always increasing. $A_U$ is larger than $A_G$ and $A_B$ since the anisotropies of both modules are added.

3.4. Vibrational Properties

According to the group theory, the space group $Cmcm$ has the following Raman, infrared (IR), and silent phonon modes at the zone center $\mathsf{\Gamma }$ = 5$A_g$ + 4$B_{1g}$ + 2$B_{2g}$ + 4$B_{3g}$, $\mathsf{\Gamma }$ = 6$B_{1u}$ + 7$B_{2u}$ + 5$B_{3u}$, and $\mathsf{\Gamma }$ = 3$A_u$, respectively. These modes can be further classified into internal (symmetric stretching, ${\nu }_1$, symmetric bending, ${\nu }_2$, asymmetric stretching, ${\nu }_3$, and asymmetric bending, ${\nu }_4$) and external (translational, T, and rotational, R) modes of VO${}_4$ units. The calculated Raman, IR, and silent phonon frequencies, the pressure coefficients ($d\omega /dP$), and the Grüneisen parameters ($\gamma$ = -∂ln($\omega$)/∂lnV) at the $\mathsf{\Gamma }$ point appear in

Table 3. Calculated Raman, infrared, and silent phonon frequencies ($\omega$ in cm${}^{-1}$), pressure coefficients ($d\omega /dp$), and Grüneisen parameters ($\gamma$) of FeVO${}_4$-II.

Mode	Raman			Mode	Infrared			Mode	Silent
	$\mathit{\omega }$	$\mathit{d}\mathit{\omega }/\mathit{dp}$	$\mathit{\gamma }$		$\mathit{\omega }$	$\mathit{d}\mathit{\omega }/\mathit{dp}$	$\mathit{\gamma }$		$\mathit{\omega }$	$\mathit{d}\mathit{\omega }/\mathit{dp}$	$\mathit{\gamma }$
$T\left(B_{3g}\right)$	153.7	1.37	0.57	$T\left(B_{1u}\right)$	153.7	0.24	0.15	$T\left(A_u\right)$	164.4	1.29	0.77
$T\left(B_{1g}\right)$	177.4	−1.05	−0.82	$T\left(B_{1u}\right)$	180.8	1.68	0.86	${\nu }_2\left(A_u\right)$	220.3	0.55	0.23
$R\left(B_{1g}\right)$	219.6	0.52	0.22	$R\left(B_{3u}\right)$	201.2	−0.01	−0.03	${\nu }_2\left(A_u\right)$	404.9	7.45	1.68
$T\left(A_g\right)$	221.8	3.23	1.29	${\nu }_2\left(B_{2u}\right)$	274.9	1.12	0.37
$R\left(B_{2g}\right)$	261.8	4.25	1.47	${\nu }_4\left(B_{3u}\right)$	298.8	1.51	0.46
${\nu }_2\left(A_g\right)$	355.1	0.30	0.07	$T\left(B_{3u}\right)$	315.1	3.82	0.25
${\nu }_4\left(B_{3g}\right)$	361.7	1.36	0.33	${\nu }_4\left(B_{2u}\right)$	327.9	−0.88	−1.06
${\nu }_2\left(B_{2g}\right)$	366.4	3.78	0.94	$R\left(B_{1u}\right)$	363.1	2.41	0.48
${\nu }_4\left(B_{1g}\right)$	374.0	3.58	0.88	${\nu }_4\left(B_{2u}\right)$	381.5	6.61	1.61
${\nu }_4\left(A_g\right)$	394.2	1.55	0.36	${\nu }_4\left(B_{1u}\right)$	398.3	4.93	1.13
${\nu }_4\left(B_{3g}\right)$	431.3	5.90	1.25	${\nu }_2\left(B_{2u}\right)$	420.6	4.38	0.92
${\nu }_3\left(B_{1g}\right)$	652.8	5.71	0.78	${\nu }_3\left(B_{3u}\right)$	654.2	5.64	0.77
${\nu }_3\left(A_g\right)$	751.0	5.74	0.72	${\nu }_3\left(B_{2u}\right)$	735.8	5.47	0.67
${\nu }_1\left(A_g\right)$	918.9	1.33	0.12	${\nu }_3\left(B_{1u}\right)$	866.8	2.58	0.26
${\nu }_3\left(B_{3g}\right)$	922.0	2.44	0.24	${\nu }_3\left(B_{2u}\right)$	924.3	1.38	0.13

.

According to Table 3, one Raman, $R\left(B_{1g}\right)$, and two IR modes, $R\left(B_{3u}\right)$ and ${\nu }_4\left(B_{2u}\right)$, have negative pressure coefficients and Grüneisen parameters (soften with pressure), which are related to the instability of the CrVO${}_4$-type structure under pressure [2]. Such behavior was already observed in InVO${}_4$ [5] and CrVO${}_4$ [20]. The soft translational Raman $T\left(B_{1g}\right)$ mode involves a translation of the VO${}_4$ polyhedra accompanied by a rotation of the FeO${}_6$ units around the z-axis as rigid units, see Figure 7, whose freedom of movement is affected by the decrease in the lattice parameter b, which is the most compressible under pressure. In contrast, the remaining Raman modes harden with pressure. The IR soft $R\left(B_{3u}\right)$ mode is quite similar to the $R\left(B_{1g}\right)$, but in this case, the FeO${}_6$ vibrates asymmetrically. While the other soft $B_{2u}$ mode implies the asymmetric bending in the VO${}_4$ polyhedra with a translation of the Fe cation in the c direction, which produces a stretching of the apical $d_{\mathrm{Fe}-\mathrm{O}1}$ interatomic distances. Therefore, the softening of R and IR modes is related to the compression of the y-axis. Interestingly, phosphates with CrVO${}_4$-type structure, such as InPO${}_4$ and TiPO${}_4$, only present some IR soft modes. The main difference being that Ti(In)PO${}_4$ is more compressible in the c(b)-axis [7].

The phonon dispersion relation and the phonon density of states (DOS) are displayed in Figure 8 at (a) ambient pressure and (b) 9.6 GPa. As can be seen in Figure 8a, the CrVO${}_4$-type structure of FeVO${}_4$ complains with the dynamic stability criteria at P = 0, which implies that ${\omega }^2$ (q, s)> 0 for all wave vectors q and polarization s (longitudinal and transverse modes) [66]. Therefore, FeVO${}_4$-II complains with the elastic and dynamic stability criteria to be experimentally synthesized [19], such as InVO${}_4$ [5] and CrVO${}_4$ [20]. We can see no significant differences in the phonon spectrum at 9.6 GPa, which ensures that the FeVO${}_4$-II phase would be stable in a pressure range where the isostructural vanadates InVO${}_4$ and CrVO${}_4$ undergo a phase transition to the wolframite-type structure [5,20].

The phonon DOS of FeVO${}_4$-II is more similar to CrVO${}_4$ than InVO${}_4$ because In is more than double the mass of Cr and Fe. The partial phonon DOS of the first lower frequency zone (up to ≈ 505 cm${}^{-1}$) is due to external (R and T) and internal (${\nu }_2$ and ${\nu }_4$) phonon modes, with Fe and V similar contributions. In contrast, the two high-frequency zones are due to pure internal vibrations ${\nu }_1$ and ${\nu }_3$, which indicates that phonon DOS is only due to O and V.

4. Conclusions

First-principles calculations were performed to study FeVO${}_4$ in the CrVO${}_4$-type structure up to 9.6 GPa. Our results indicate that this phase is elastically and dynamically stable for this pressure range. FeVO${}_4$ presents an indirect band-gap $E_g$ = 2.11 eV, similar to the value reported for the triclinic phase of FeVO${}_4$ at ambient pressure. The Laplacian charge density at BCPs in both polyhedra is larger than those reported for AMoO${}_4$ and AWO${}_4$ compounds and increases with pressure. It was observed that the y-axis is the most compressible, related to the reduction in the interatomic distances of the FeO${}_6$ polyhedra. These changes in the structural parameters have strong effects on the phonon modes with negative pressure coefficients. We found that the shear modulus G softens with pressure, which suggests that the FeVO${}_4$-II phase could undergo a phase transition below 10 GPa as was predicted in a previous study and happen in CrVO${}_4$ and InVO${}_4$, where these compounds have a phase transition from CrVO${}_4$-type to the wolframite structure. We hope this work encourages other research groups to perform high-pressure characterization studies with the diamond-anvil cell device in order to increase the knowledge of less studied orthovanadates with CrVO${}_4$-type structure, such as FeVO${}_4$ and TlVO${}_4$.