High-Pressure Elastic , Vibrational and Structural Study of Monazite-Type GdPO 4 from Ab Initio Simulations

The GdPO4 monazite-type has been studied under high pressure by first principles calculations in the framework of density functional theory. This study focuses on the structural, dynamical, and elastic properties of this material. Information about the structure and its evolution under pressure, the equation of state, and its compressibility are reported. The evolution of the Raman and Infrared frequencies, as well as their pressure coefficients are also presented. Finally, the study of the elastic constants provides information related with the elastic and mechanical properties of this compound. From our results, we conclude that monazite-type GdPO4 becomes mechanically unstable at 54 GPa; no evidence of soft phonons has been found up to this pressure at the zone center.


Introduction
The APO 4 orthophosphates compounds (A = trivalent cation) are materials analogous to orthoarsenates, orthovanadates, and orthosilicates.The size of the ionic radio of the A cation defines the two different structures where this family of compounds crystallizes.APO 4 compounds with an ionic radio lower than Gd crystallize in the tetragonal zircon structure (I4 1 /amd, space group 141 with Z = 4).All other orthophosphates (A = La to Gd) crystallize into a monoclinic lower symmetry phase, that is, a monazite structure (P2 1 /n, space group 14 with Z = 4) [1].The monazite structure is isostructural to cerium phosphate mineral (monazite).This structure can be viewed as being composed by alternating edge-sharing AO 9 polyhedra, with a structural distortion derived from a rotation of the PO 4 tetrahedra and a lateral shift of the (100) planes that reduce the symmetry from the tetragonal of the zircon structure to the monoclinic symmetry.
Due to their numerous potential applications, the orthophosphate family has attracted a lot of interest in the last decades.The study of the mechanical properties of these materials can be relevant, since some orthophosphates have shown to be promising materials for oxidation-resistant ceramic toughening [2] and to have interesting luminescent and optical properties with a wide range of potential applications [3][4][5][6].
This family of minerals appears in nature as accessory mineral in rhyolites and granitoid rocks [1].Therefore, the study of orthophosphates may be of importance in petrology [7], chemistry, and mineral physics.These materials can be employed as geochronometers for geological dating [8], because of their doping with radioactive elements that possess a very large half-life.Additionally, since orthophosphates incorporate rare earth elements, they can also be used to determine the distribution of rare earth elements in our Planet [9].Another interesting application of these compounds is related to their use as potential matrix for nuclear waste disposal [5,[10][11][12].GdPO 4 and other orthophosphates have been suggested as inclusion compounds to made nuclear fuel canisters due to their high neutron absorption [13].Therefore, the study of these orthophosphates under high pressure is very promising for future practical applications.Although there are several studies focused on pressure-induced transition and on the structural and vibrational properties of monazite phosphates under pressure [7,[14][15][16][17][18], these studies were carried out only up to 30 GPa.A systematic understanding of their properties has not yet been achieved; the study of their elastic and mechanical properties is scarce [5,19,20].Clearly, additional researches are needed to deepen our knowledge of this group of materials.A systematic investigation of the properties of monazite phosphates (particularly at high pressure) is necessary, especially if we want to use them in efficient and safe applications.
The ab initio methods have been successfully applied to study the thermodynamic properties of monazite rare earth orthophosphates [21][22][23][24].Moreover, the use of ab initio simulations has become a strong and very well supported method to investigate the properties of materials under high pressures [25][26][27][28].In particular, in the case of monazite-type orthophosphates under pressure, the ab initio simulations have been employed with good results [29,30].
The aim of the present paper is to contribute to the knowledge of the properties of orthophosphates at high pressure.In this work, we perform a detailed study of the structural, dynamic, and elastic properties of monazite-type GdPO 4 under pressures up to 60 GPa, going from first principles simulations.This method has proven to be quite efficient, in order to study a variety of properties on related compounds.

Simulations Details
The influence of pressure on the crystal structure, mechanical and vibrational properties of monazite-type GdPO 4 has been analyzed using ab initio simulations.The study was based on the Density Functional Theory (DFT) [31] employing the Vienna Ab initio Simulation Package (VASP) [32][33][34] and pseudopotentials generated with the Projector-Augmented Wave scheme (PAW) [35,36], with the f electrons included into the core.This has proven to give good results in the study of monazite-type structures [20,30,37].To reach accurate results, the set of plane waves was extended up to a 520-eV cutoff energy and the exchange-correlation energy was expressed by means of the Generalized-Gradient Approximation (GGA) with the Perdew-Burke-Ernzerhof for solids (PBEsol) functional [38].A dense grid of Monkhorst-Pack [39] k-special points was utilized to perform the integrations on the Brillouin Zone (BZ).The achieved convergence was 1 meV per formula unit in the total energy.All the structural parameters were optimized by minimizing the forces on the atoms and the stress tensor, at selected volumes.This method has been successfully applied to study the phase stability, and structural properties of semiconductors under high pressure [40].
To carry out the study of the mechanical properties of GdPO 4 , we have evaluated the elastic constants, which describe the mechanical properties of a material in the region of small deformations.The elastic constants can be obtained by computing the macroscopic stress of a small strain with the use of the stress theorem [41,42].In the present study, we employ the method implemented in the VASP code: the ground state and fully optimized structures were strained in different directions taking into account their symmetry [43].
The phonons study was performed using the supercell method [44].The diagonalization of the dynamical matrix provides the frequencies of the normal modes.Furthermore, these calculations also provide the symmetry and eigenvectors of the phonon modes at the Γ point.The calculation of the phonon dispersion curves along the main directions within the BZ was carried out with a supercell of size 2 × 2 × 2. Following the conclusions of our previous studies and those of Blanca-Romero et al. [45], in all the calculations, we neglect the spin-orbit interaction, because the structural properties are barely affected [15,17,26,29,30,37].

Structural Properties
The monoclinic structure of monazite, space group P2 1 /n [1] has been described in the introduction.A schematic view of the crystal structure is presented in Figure 1.It can be seen as an alternating chain of phosphorus-oxygen PO 4 tetrahedra and trivalent cation-oxygen GdO 9 polyhedra.

Structural Properties
The monoclinic structure of monazite, space group P21/n [1] has been described in the introduction.A schematic view of the crystal structure is presented in Figure 1.It can be seen as an alternating chain of phosphorus-oxygen PO4 tetrahedra and trivalent cation-oxygen GdO9 polyhedra.The calculated structural parameters and atomic coordinates at ambient pressure are given in Tables 1 and 2, and are in a very close agreement with the available theoretical and experimental results [11,14,18,20,21].Previous theoretical results [20] differ slightly due to the exchange-correlation functional used in the calculations.[18,20] were obtained in a different setting than the one employed in the present work.The calculated structural parameters and atomic coordinates at ambient pressure are given in Tables 1 and 2, and are in a very close agreement with the available theoretical and experimental results [11,14,18,20,21].Previous theoretical results [20] differ slightly due to the exchange-correlation functional used in the calculations.[18,20] were obtained in a different setting than the one employed in the present work.The evolution of the structural parameters with pressure is displayed in Figure 2a.The calculations have been already compared with experiments [14,18], showing an excellent agreement.The obtained linear axial compressibilities of each axis are: Therefore, from the present results and previous studies [14,18], it can be concluded that the compression of GdPO 4 , as in other monazite-type phosphates, is not isotropic.The a-axis is the most compressible one and the c-axis the least compressible one.As shown in Figure 2a, there is a tendency for the unit-cell parameter a to approach the value of c.Moreover, in Figure 2b, it can be seen that the monoclinic β angle decreases under compression.As a result of these two facts, there is a gradual symmetrization of the monazite structure under compression.The evolution of the structural parameters with pressure is displayed in Figure 2a.The calculations have been already compared with experiments [14,18], showing an excellent agreement.The obtained linear axial compressibilities of each axis are:   = .Therefore, from the present results and previous studies [14,18], it can be concluded that the compression of GdPO4, as in other monazitetype phosphates, is not isotropic.The a-axis is the most compressible one and the c-axis the least compressible one.As shown in Figure 2a, there is a tendency for the unit-cell parameter a to approach the value of c.Moreover, in Figure 2b, it can be seen that the monoclinic β angle decreases under compression.As a result of these two facts, there is a gradual symmetrization of the monazite structure under compression.The bulk modulus, B 0 , and its first pressure derivative were obtained by fitting the set of theoretical energy-volume data with a third-order Birch-Murnaghan (BM) Equation of State (EOS) [46].The results (B 0 = 138.3GPa, and B 0 = 4.02) are summarized in Table 3 and are compared with previous results of B 0 reported for GdPO 4 : the agreement between experiments is good.From our calculations, we have determined the pressure dependence of the polyhedral volumes and distortions for monazite GdPO 4 .Figure 3a presents the relative volume of the PO 4 and GdO 9 polyhedra and that of the unit-cell.As it can be observed, the PO 4 tetrahedron is highly incompressible, while the GdO 9 polyhedron is much more compressible.Actually, in the monazite-type GdPO 4 , the volume change of the GO 9 polyhedra seems to be responsible for most of the volume reduction induced by pressure within the structure.In fact, the phosphate tetrahedra can be treated as a rigid unit in the monazite-type phosphate [17,18,48].The P-O distances do not vary greatly within the tetrahedra.At zero pressure, two distances are 1.55 Å, and the two other 1.56 Å, and 1.54 Å, while at 53 GPa these distances decrease to 1.51 Å × 2, 1.52 Å, and 1.49 Å, respectively.Indeed, the distortion index, calculated with VESTA [49], for the PO 4 polyhedron (see Figure 3b) is very small, and it increases slowly with compression.Regarding the GdO 9 polyhedra at zero pressure, the central gadolinium atom is connected to nine oxygen atoms with eight bonds whose lengths vary from 2.34 to 2.55 Å, with one length at 2.78 Å.At 53 GPa the lengths range from 2.16 to 2.46 Å and the largest one changes to 2.53 Å.The GdO 9 distortion index do not vary linearly, as seen in Figure 3b.There is a decrease around 4 GPa, followed by a considerable increase in pressure, in a similar manner than in the isostructural LaPO 4 and CePO 4 .The GdO 9 polyhedron is more distorted than the PO 4 tetrahedron: its distortion index changes from approximately 0.0381 at ambient pressure to 0.0440 at 54 GPa.The bulk modulus, B0, and its first pressure derivative were obtained by fitting the set of theoretical energy-volume data with a third-order Birch-Murnaghan (BM) Equation of State (EOS) [46].The results (B0 = 138.3GPa, and B′0 = 4.02) are summarized in Table 3 and are compared with previous results of B0 reported for GdPO4: the agreement between experiments is good.g Reference [48], h Reference [24].
From our calculations, we have determined the pressure dependence of the polyhedral volumes and distortions for monazite GdPO4.Figure 3a presents the relative volume of the PO4 and GdO9 polyhedra and that of the unit-cell.As it can be observed, the PO4 tetrahedron is highly incompressible, while the GdO9 polyhedron is much more compressible.Actually, in the monazitetype GdPO4, the volume change of the GO9 polyhedra seems to be responsible for most of the volume reduction induced by pressure within the structure.In fact, the phosphate tetrahedra can be treated as a rigid unit in the monazite-type phosphate [17,18,48].The P-O distances do not vary greatly within the tetrahedra.At zero pressure, two distances are 1.55 Å, and the two other 1.56 Å, and 1.54 Å, while at 53 GPa these distances decrease to 1.51 Å × 2, 1.52 Å, and 1.49 Å, respectively.Indeed, the distortion index, calculated with VESTA [49], for the PO4 polyhedron (see Figure 3b) is very small, and it increases slowly with compression.Regarding the GdO9 polyhedra at zero pressure, the central gadolinium atom is connected to nine oxygen atoms with eight bonds whose lengths vary from 2.34 to 2.55 Å, with one length at 2.78 Å.At 53 GPa the lengths range from 2.16 to 2.46 Å and the largest one changes to 2.53 Å.The GdO9 distortion index do not vary linearly, as seen in Figure 3b.There is a decrease around 4 GPa, followed by a considerable increase in pressure, in a similar manner than in the isostructural LaPO4 and CePO4.The GdO9 polyhedron is more distorted than the PO4 tetrahedron: its distortion index changes from approximately 0.0381 at ambient pressure to 0.0440 at 54 GPa.

Elastic Properties
GdPO4 belongs to the monoclinic space group n° 14; it therefore has 13 independent elastic constants [50]: C11, C12, C13, C15, C22, C23, C25, C33, C35, C44, C46, C55, and C66 in the Voigt notation.Table 4 summarizes our results for the set of Cij elastic constants at 0 GPa.The values here reported are in overall good agreement with those of Feng et al. [20].The diagonal elastic constants C11, C22, C33, are the highest ones, while C25, C35, C46 are the lowest ones.It must be noted that in [20], the cell parameters are underestimated, as commented above.When a non-zero uniform stress is applied to a crystal, the relevant magnitudes that describe their elastic properties are the elastic stiffness coefficients Bij [51].In the special case of a hydrostatic pressure, P, applied to a monoclinic crystal, the elastic stiffness coefficients are related to the elastic constants through the following relationships [51,52]:  Figure 4 shows the pressure dependence of the elastic stiffness coefficients, Bij.It can be seen that, B11, B22, B33, B12, B13, and B23 increase with pressure in the whole pressure range studied, especially in the diagonal components, while B44, B55 and B66 decrease under compression.On the other hand, B15 is almost unaffected by pressure.It should be noted that B25, B35, B46 have similar negatives values and that they decrease under compression.4 summarizes our results for the set of C ij elastic constants at 0 GPa.The values here reported are in overall good agreement with those of Feng et al. [20].The diagonal elastic constants C 11 , C 22 , C 33 , are the highest ones, while C 25 , C 35 , C 46 are the lowest ones.It must be noted that in [20], the cell parameters are underestimated, as commented above.When a non-zero uniform stress is applied to a crystal, the relevant magnitudes that describe their elastic properties are the elastic stiffness coefficients B ij [51].In the special case of a hydrostatic pressure, P, applied to a monoclinic crystal, the elastic stiffness coefficients are related to the elastic constants through the following relationships [51,52]:   At zero pressure, a crystal is mechanically stable when the Born stability criteria are fulfilled [53].However, when a hydrostatic pressure is applied to the crystal, the generalized Born stability criteria [52,54] must be employed to study the mechanical stability.This means that the matrix Bij must be positive definite.Consequently, the stability conditions for a monoclinic crystal are given by the following conditions:

Elastic Properties
The simulations of the elastic constants have been extended up to 60 GPa, well above the maximum experimental pressure reported [14,15].The generalized stability criteria for GdPO4 as function of the pressure are presented in Figure 5.At 54 GPa, one of the M6 criterion is violated, making the system mechanically unstable at this pressure.No evidence of soft phonons was found at the Γ point in this pressure range.These results are consistent with the experimental and theoretical reports that confirm that there is no phase transition induced by pressure in the monazite structure up to the highest pressure reached in those studies (30 GPa) [14,17].At zero pressure, a crystal is mechanically stable when the Born stability criteria are fulfilled [53].However, when a hydrostatic pressure is applied to the crystal, the generalized Born stability criteria [52,54] must be employed to study the mechanical stability.This means that the matrix B ij must be positive definite.Consequently, the stability conditions for a monoclinic crystal are given by the following conditions: The simulations of the elastic constants have been extended up to 60 GPa, well above the maximum experimental pressure reported [14,15].The generalized stability criteria for GdPO 4 as function of the pressure are presented in Figure 5.At 54 GPa, one of the M 6 criterion is violated, making the system mechanically unstable at this pressure.No evidence of soft phonons was found at the Γ point in this pressure range.These results are consistent with the experimental and theoretical reports that confirm that there is no phase transition induced by pressure in the monazite structure up to the highest pressure reached in those studies (30 GPa) [14,17].The bulk modulus (B), the shear modulus (G), the Young modulus (E), and the Poisson's ratio (ν) describe the major elastic properties of a material.Analytical expressions can be obtained for these moduli in the Voigt, Reuss and Hill approximations [55][56][57] from the elastic stiffness coefficients, Bij.Hill has proved that the Voigt and Reuss approximations are limits and pointed out that the arithmetic mean of the two bounds can be considered as the actual B and G elastic moduli.In the case of a monoclinic crystal, the bulk and shear moduli can be expressed in the three approximations as [ where Sij refers to components of the elastic compliances tensor, and subscripts V, R, and H stand for Voigt, Reuss and Hill respectively.The Young modulus, E, and the Poisson's ratio, ν, are obtained with the expressions [59]: The bulk modulus (B), the shear modulus (G), the Young modulus (E), and the Poisson's ratio (ν) describe the major elastic properties of a material.Analytical expressions can be obtained for these moduli in the Voigt, Reuss and Hill approximations [55][56][57] from the elastic stiffness coefficients, B ij .Hill has proved that the Voigt and Reuss approximations are limits and pointed out that the arithmetic mean of the two bounds can be considered as the actual B and G elastic moduli.In the case of a monoclinic crystal, the bulk and shear moduli can be expressed in the three approximations as [58]: where S ij refers to components of the elastic compliances tensor, and subscripts V, R, and H stand for Voigt, Reuss and Hill respectively.The Young modulus, E, and the Poisson's ratio, ν, are obtained with the expressions [59]: where the subscript X refers to the symbols V, R, and H.
In order to discuss the elastic properties of GdPO 4 at ambient pressure, the elastic moduli at 0 GPa are summarized in Table 5.The bulk modulus, B H = 134.47GPa is in rather good agreement with our theoretical value of B 0 = 138.3GPa, obtained from a fit with a third order Birch-Murnaghan equation of state [46], and with the experimental values of B 0 = 137 GPa previously reported [47].As for E, G, and B/G elastic moduli, as well as the Poisson's ratio, ν, at zero applied pressure, it should be noted once more that there is good agreement between our calculations and the experimental results.These facts demonstrate the quality of our simulations.The pressure dependence of the elastic moduli up to 56.1 GPa is presented in Figure 6a-d.It can be observed that B H , and ν H increase with pressure.However, E H and G H increase under compression reaching maximum values at 11 GPa and 9 GPa, respectively; above that pressure they decrease.All the elastic moduli change dramatically around 54 GPa, since the monazite structure becomes mechanically unstable.The Poisson's ratio gives information about the characteristics of the bonding forces and the chemical bonding.The value of the Poisson's ratio in the Hill approximation is ν H = 0.29 at 0 GPa.A value of ν > 0.25 indicates that the interatomic bonding forces are predominantly central and that ionic bonding is predominant against covalent bonding [60,61].The increase of ν upon pressure can be interpreted as an increase of the metallization.
where the subscript X refers to the symbols V, R, and H.
In order to discuss the elastic properties of GdPO4 at ambient pressure, the elastic moduli at 0 GPa are summarized in Table 5.The bulk modulus, BH = 134.47GPa is in rather good agreement with our theoretical value of B0 = 138.3GPa, obtained from a fit with a third order Birch-Murnaghan equation of state [46], and with the experimental values of B0 = 137 GPa previously reported [47].As for E, G, and B/G elastic moduli, as well as the Poisson's ratio, ν, at zero applied pressure, it should be noted once more that there is good agreement between our calculations and the experimental results.These facts demonstrate the quality of our simulations.The pressure dependence of the elastic moduli up to 56.1 GPa is presented in Figure 6a-d.It can be observed that BH, and νH increase with pressure.However, EH and GH increase under compression reaching maximum values at 11 GPa and 9 GPa, respectively; above that pressure they decrease.All the elastic moduli change dramatically around 54 GPa, since the monazite structure becomes mechanically unstable.The Poisson's ratio gives information about the characteristics of the bonding forces and the chemical bonding.The value of the Poisson's ratio in the Hill approximation is νH = 0.29 at 0 GPa.A value of ν > 0.25 indicates that the interatomic bonding forces are predominantly central and that ionic bonding is predominant against covalent bonding [60,61].The increase of ν upon pressure can be interpreted as an increase of the metallization.The B/G ratio is 2.05 and it increases with pressure (see Figure 7).Therefore, GdPO4 monazite is more resistant to volume compression than to shear deformation (B > G).Moreover, according to Pugh criterion, monazite-type GdPO4 behaves like a ductile material, since B/G is smaller than 1.75 [62], in the whole pressure range studied (materials with B/G < 1.75 behave as brittle materials).Our B/G value is consistent with the experimental result of Du et al. [47], although it is somewhat larger than the theoretical one of reference [24].This is probably due to the small value of B0 reported in that reference.However, it is clear that monazite GdPO4 is a ductile material.The B/G ratio is 2.05 and it increases with pressure (see Figure 7).Therefore, GdPO 4 monazite is more resistant to volume compression than to shear deformation (B > G).Moreover, according to Pugh criterion, monazite-type GdPO 4 behaves like a ductile material, since B/G is smaller than 1.75 [62], in the whole pressure range studied (materials with B/G < 1.75 behave as brittle materials).Our B/G value is consistent with the experimental result of Du et al. [47], although it is somewhat larger than the theoretical one of reference [24].This is probably due to the small value of B 0 reported in that reference.However, it is clear that monazite GdPO 4 is a ductile material.The B/G ratio is 2.05 and it increases with pressure (see Figure 7).Therefore, GdPO4 monazite is more resistant to volume compression than to shear deformation (B > G).Moreover, according to Pugh criterion, monazite-type GdPO4 behaves like a ductile material, since B/G is smaller than 1.75 [62], in the whole pressure range studied (materials with B/G < 1.75 behave as brittle materials).Our B/G value is consistent with the experimental result of Du et al. [47], although it is somewhat larger than the theoretical one of reference [24].This is probably due to the small value of B0 reported in that reference.However, it is clear that monazite GdPO4 is a ductile material.

Vibrational Properties
The phonon dispersion of GdPO 4 along the main paths on the Brillouin zone is presented in Figure 8. Three frequency regions are separated by two gaps: a little one between the low-and the medium-frequency regions, and a bigger one between the medium-and the high-frequencies regions.The Partial Density Of States (PDOS) for phonons (see Figure 9) allows us to analyze these vibrations.The PDOS presents three zones separated by two gaps.The first zone, corresponding to the low-frequencies, down to 311 cm −1 , is composed of vibrations of Gd atoms and GdO 9 polyhedra, as well as a small contribution of PbO 4 tetrahedra.The intermediate-frequency zone, from 381 cm −1 to 600 cm −1 , is mainly due to vibrations of the PbO 4 tetrahedra with a very small contribution from the GdO 9 polyhedra.A large gap of ≈340 cm −1 separates this zone with the high-frequency region, which corresponds to frequencies between 942 cm −1 and 1073 cm −1 that are only due to vibrations from PbO 4 units.Therefore, the vibrational spectra of the GdPO 4 monazite could be interpreted in term of the modes of the PbO 4 tetrahedra, which can be considered as independent units in the monazite structure, and GdO 9 polyhedra.Indeed, the vibrational spectrum of monazites ABO 4 has been interpreted in terms of internal modes (bending and stretching modes) associated with the tetrahedron BO 4 and external modes (translational and rotational) which involves movement of both A and BO 4 ions [17,63].

Vibrational Properties
The phonon dispersion of GdPO4 along the main paths on the Brillouin zone is presented in Figure 8. Three frequency regions are separated by two gaps: a little one between the low-and the medium-frequency regions, and a bigger one between the medium-and the high-frequencies regions.The Partial Density Of States (PDOS) for phonons (see Figure 9) allows us to analyze these vibrations.The PDOS presents three zones separated by two gaps.The first zone, corresponding to the lowfrequencies, down to 311 cm −1 , is composed of vibrations of Gd atoms and GdO9 polyhedra, as well as a small contribution of PbO4 tetrahedra.The intermediate-frequency zone, from 381 cm −1 to 600 cm −1 , is mainly due to vibrations of the PbO4 tetrahedra with a very small contribution from the GdO9 polyhedra.A large gap of ≈340 cm −1 separates this zone with the high-frequency region, which corresponds to frequencies between 942 cm −1 and 1073 cm −1 that are only due to vibrations from PbO4 units.Therefore, the vibrational spectra of the GdPO4 monazite could be interpreted in term of the modes of the PbO4 tetrahedra, which can be considered as independent units in the monazite structure, and GdO9 polyhedra.Indeed, the vibrational spectrum of monazites ABO4 has been interpreted in terms of internal modes (bending and stretching modes) associated with the tetrahedron BO4 and external modes (translational and rotational) which involves movement of both A and BO4 ions [17,63].Through group theory analysis, it can be established that the monazite structure has 72 vibrational modes at the zone center: 36 optical Raman-active modes (18Ag + 18Bg), 33 optical IRactive modes (17Au + 16Bu), and 3 acoustic modes (1Au + 2Bu).These vibrational modes can be interpreted as 36 internal (ν1, ν2, ν3 and ν4) and 36 external (translational (T) and rotational (R)) modes [17,63].
The Raman spectrum of many monazite-type phosphates (e.g., BiPO4 LaPO4, CePO4, PrPO4 and GdPO4) at ambient conditions has been previously studied by different authors [18,29,64,65].BiPO4 LaPO4, CePO4, and PrPO4 have also been studied under compression [17,[66][67][68].The frequency distribution of Raman modes is very similar for all of them.As in the case of GdPO4, phonon gaps are observed between the three frequency regions in monazite phosphates; they are also observed in the monazite structure of chromates and selenates [30].
GdPO4 has a distribution of Raman modes similar to other monazite phosphates [17].The Raman mode frequencies and their pressure coefficients at room pressure are summarized in Table 6.The mode Grüneisen parameters [69], that provide a dimensionless representation of the response to compression are also reported.
It can be observed that the Raman spectrum of GdPO4 monazite could be divided into three regions: the low-frequency region, up to 310.91 cm −1 , corresponds to the eighteen external or lattice T and R modes (9Ag + 9Bg); the medium-frequency region, between 381.30 and 603.46 cm −1 , corresponds to the ten internal bending modes (5Ag + 5Bg); and the high-frequency region, above 942.72 cm −1 , corresponds to the eight internal stretching modes (4Ag + 4Bg).
As seen in Figure 10, which shows the theoretical pressure dependence of the Raman-active mode frequencies, most modes harden under compression, as is usual in most materials.However, there are a few modes between 84.78 and 153.94 cm −1 whose frequency changes very slightly under compression.Through group theory analysis, it can be established that the monazite structure has 72 vibrational modes at the zone center: 36 optical Raman-active modes (18A g + 18B g) , 33 optical IR-active modes (17A u + 16B u ), and 3 acoustic modes (1A u + 2B u ).These vibrational modes can be interpreted as 36 internal (ν 1 , ν 2 , ν 3 and ν 4 ) and 36 external (translational (T) and rotational (R)) modes [17,63].
GdPO 4 has a distribution of Raman modes similar to other monazite phosphates [17].The Raman mode frequencies and their pressure coefficients at room pressure are summarized in Table 6.The mode Grüneisen parameters [69], that provide a dimensionless representation of the response to compression are also reported.
It can be observed that the Raman spectrum of GdPO 4 monazite could be divided into three regions: the low-frequency region, up to 310.91 cm −1 , corresponds to the eighteen external or lattice T and R modes (9A g + 9B g ); the medium-frequency region, between 381.30 and 603.46 cm −1 , corresponds to the ten internal bending modes (5A g + 5B g ); and the high-frequency region, above 942.72 cm −1 , corresponds to the eight internal stretching modes (4A g + 4B g ).
As seen in Figure 10, which shows the theoretical pressure dependence of the Raman-active mode frequencies, most modes harden under compression, as is usual in most materials.However, there are a few modes between 84.78 and 153.94 cm −1 whose frequency changes very slightly under compression.External modes (low frequency region) involve movements of the trivalent Gd atoms, and in phosphate monazites ABO 4 , their frequencies greatly depend on the mass of the A atom [17,70].These external modes show very different pressure coefficients, since they involve different Gd-O bonds, with different compressibility, as seen in Section 2 [14].Due to the marked differences between the pressure dependence of the different modes, see Table 6, crossing and anti-crossing phenomena are observed in this region.In particular, one of the modes most affected by compression is the A g mode at 264.93 cm −1 ; The lowest-frequency B g mode is almost not affected by pressure.In previous studies of monazite chromates [37] and phosphates [17], softening of the lower Raman modes was related to a pressure-driven instability of the monazite structure, which could result in a phase transition above 30 GPa [17].In the case of the monazite-type GdPO 4 , the lower modes are little affected by compression.As commented in a previous section, the extension of the simulation up to 60 GPa does not show evidence of soft phonons mode at the Γ point.The pressure-driven instability seems to be related to mechanical instability, rather than dynamical instability.External modes (low frequency region) involve movements of the trivalent Gd atoms, and in phosphate monazites ABO4, their frequencies greatly depend on the mass of the A atom [17,70].These external modes show very different pressure coefficients, since they involve different Gd-O bonds, with different compressibility, as seen in Section 2 [14].Due to the marked differences between the pressure dependence of the different modes, see Table 6, crossing and anti-crossing phenomena are observed in this region.In particular, one of the modes most affected by compression is the Ag mode at 264.93 cm −1 ; The lowest-frequency Bg mode is almost not affected by pressure.In previous studies of monazite chromates [37] and phosphates [17], softening of the lower Raman modes was related to a pressure-driven instability of the monazite structure, which could result in a phase transition above 30 GPa [17].In the case of the monazite-type GdPO4, the lower modes are little affected by compression.As commented in a previous section, the extension of the simulation up to 60 GPa does not show evidence of soft phonons mode at the Γ point.The pressure-driven instability seems to be related to mechanical instability, rather than dynamical instability.Internal bending motions of the PO 4 tetrahedron (medium frequencies region) have smaller pressure coefficients.Two modes with frequencies A g = 506.02and B g = 530 cm −1 at room pressure are the less affected by pressure, while the mode most sensitive to pressure in this region is the B g mode with a frequency 381.30cm −1 at room pressure.The crossover of B g (499.32 cm −1 ) and A g (506.02cm −1 ) modes is observed at ≈4 GPa in Figure 10, due to the different pressure dependences.This behavior could be related to the non-isotropic compression of monazite.This can cause the lower-frequency B g mode to move faster towards high frequencies than the higher-frequency mode.
The internal stretching modes (high frequency region) of the PO 4 tetrahedron have similar pressure coefficients, and they are among the modes whose frequency increases faster under compression.
Regarding the infrared modes of monazite GdPO 4 , low-, medium-, and high frequencies modes have the same general behavior with pressure than the Raman modes have.This also happens for monazite chromates [34].In particular, the modes less affected by pressure is the first A u mode (91.057 cm −1 ); see Table 7 and Figure 11.
Table 7. Infrared frequencies at ambient pressure, their pressure coefficients, and Grüneisen parameters (γ).The pressure dependence of frequency has been fitted with a second order polynomial:  The Grüneisen parameters in the high frequency region are very similar, while there are large differences between the Grüneisen parameters for the lower frequencies (Tables 6 and 7).Therefore, there are large differences between the restoring force on the atoms associated with the lowest modes.

Conclusions
The effects of high pressure on the structure of monazite-type GdPO4 were studied from first principles, ab initio calculations.The simulations were extended up to 60 GPa, well above the maximum experimental pressure previously reported (30 GPa).The evolution of the lattice parameters under compression, the equations of state, and the axial and polyhedral compressibilities are reported.Moreover, the under-pressure distortion of the polyhedral units PO4 and GdO9 is shown.The compression of monazite-type GdPO4 is not isotropic, and the GdO9 polyhedra account for almost all the volume reduction under pressure.
The elastic constants and elastic stiffness coefficients were accurately determined, through which the elastic behavior of monazite-type GdPO4 at high pressure was analyzed.The evolution with pressure of the B, G, and E elastic moduli, ν Poisson's ratio, and the B/G ratio was reported.At all pressures, this compound has shown to be ductile and more resistive to volume compression than to shear deformation (B > G).Furthermore, the high-pressure mechanical stability of this monoclinic compound was studied.It was found that monazite-type GdPO4 becomes mechanically unstable above 54 GPa.
Finally, through the total and partial density of states, the contribution of each polynomial units to the vibrational modes was discussed.This contribution is similar to that found in other monazite phosphates.The evolution of the Raman-and infrared active modes of the monazite-type GdPO4 as The Grüneisen parameters in the high frequency region are very similar, while there are large differences between the Grüneisen parameters for the lower frequencies (Tables 6 and 7).Therefore, there are large differences between the restoring force on the atoms associated with the lowest modes.

Conclusions
The effects of high pressure on the structure of monazite-type GdPO 4 were studied from first principles, ab initio calculations.The simulations were extended up to 60 GPa, well above the maximum experimental pressure previously reported (30 GPa).The evolution of the lattice parameters under compression, the equations of state, and the axial and polyhedral compressibilities are reported.Moreover, the under-pressure distortion of the polyhedral units PO 4 and GdO 9 is shown.The compression of monazite-type GdPO 4 is not isotropic, and the GdO 9 polyhedra account for almost all the volume reduction under pressure.
The elastic constants and elastic stiffness coefficients were accurately determined, through which the elastic behavior of monazite-type GdPO 4 at high pressure was analyzed.The evolution with pressure of the B, G, and E elastic moduli, ν Poisson's ratio, and the B/G ratio was reported.At all pressures, this compound has shown to be ductile and more resistive to volume compression than to shear deformation (B > G).Furthermore, the high-pressure mechanical stability of this monoclinic compound was studied.It was found that monazite-type GdPO 4 becomes mechanically unstable above 54 GPa.
Finally, through the total and partial density of states, the contribution of each polynomial units to the vibrational modes was discussed.This contribution is similar to that found in other monazite phosphates.The evolution of the Raman-and infrared active modes of the monazite-type GdPO 4 as function of pressure was analyzed.No soft-modes were found up to 60 GPa at the center of the Brillouin zone.
A possible pressure-driven phase transition or amorphization might occur due to mechanical instability at the high pressure of 54 GPa.No dynamical instability was found up to this pressure.
The results of the presents study confirm the high stability under pressure of the monazite-type GdPO 4 , which is very promising for practical applications.
Author Contributions: Both authors contributed equally to this work: simulations, analysis and interpretation of results and writing of the manuscript.

Figure 1 .
Figure 1.(Color online) Crystal structure of monazite-type GdPO4.Oxygen, phosphorus, and gadolinium atoms are shown in red, orange, and blue, respectively.The PO4 tetrahedral units and the GdO9 polyhedral units are also shown.

Figure 1 .
Figure 1.(Color online) Crystal structure of monazite-type GdPO 4 .Oxygen, phosphorus, and gadolinium atoms are shown in red, orange, and blue, respectively.The PO 4 tetrahedral units and the GdO 9 polyhedral units are also shown.

Figure 2 .
Figure 2. (a) Pressure dependence of the lattice parameters.(b) Pressure dependence of the β angle.Figure 2. (a) Pressure dependence of the lattice parameters.(b) Pressure dependence of the β angle.

Figure 2 .
Figure 2. (a) Pressure dependence of the lattice parameters.(b) Pressure dependence of the β angle.Figure 2. (a) Pressure dependence of the lattice parameters.(b) Pressure dependence of the β angle.

Figure 3 .
Figure 3. (Color online) (a) Relative variation with pressure of the polyhedral and unit cell volumes of monazite-type GdPO 4 .(b) Distortion index of the PO 4 tetrahedron and the GdO 9 polyhedron.

Figure 4 .
Figure 4. (Color online) Pressure evolution of the elastic stiffness coefficients, Bij.

Figure 4 .
Figure 4. (Color online) Pressure evolution of the elastic stiffness coefficients, B ij .

Figure 7 .
Figure 7. (Color online) Pressure dependence of the B/G ratio.

Figure 6 .
Figure 6.(Color online) (a) Pressure dependence of the bulk modulus, B. (b) Pressure dependence of the Young modulus, E. (c) Pressure dependence of the shear modulus, G.(d) Pressure dependence of the Poisson's, ν, in the approximation of Voigt, Reuss and Hill.

Figure 6 .
Figure 6.(Color online) (a) Pressure dependence of the bulk modulus, B. (b) Pressure dependence of the Young modulus, E. (c) Pressure dependence of the shear modulus, G.(d) Pressure dependence of the Poisson's, ν, in the approximation of Voigt, Reuss and Hill.

Figure 7 .
Figure 7. (Color online) Pressure dependence of the B/G ratio.

Figure 7 .
Figure 7. (Color online) Pressure dependence of the B/G ratio.

Figure 8 .
Figure 8. Phonon dispersion along the Brillouin zone of the of monazite-type GdPO4.

Figure 8 .
Figure 8. Phonon dispersion along the Brillouin zone of the of monazite-type GdPO 4.

Figure 9 .
Figure 9. (Color online) Partial and total phonon density of monazite-type GdPO4 states.

Figure 10 .
Figure 10.(Color online) Pressure dependence of the Raman modes of monazite-type GdPO4.

Figure 11 .
Figure 11.(Color online) Pressure dependence of the infrared modes of monazite-type GdPO4.

Figure 11 .
Figure 11.(Color online) Pressure dependence of the infrared modes of monazite-type GdPO 4 .

Table 3 .
Volume at 0 GPa, bulk modulus, B 0 , and its pressure derivative, B 0 , from a fit with a 3th order Birch Murnaghan EOS.

Table 3 .
Volume at 0 GPa, bulk modulus, B0, and its pressure derivative, B′0, from a fit with a 3th order Birch Murnaghan EOS.

11 C 22 C 33 C 44 C 55 C 66 C 12 C 13 C 15 C 23 C 25 C 25 C 46
Figure 4 shows the pressure dependence of the elastic stiffness coefficients, B ij .It can be seen that, B 11 , B 22 , B 33 , B 12 , B 13 , and B 23 increase with pressure in the whole pressure range studied, especially in the diagonal components, while B 44 , B 55 and B 66 decrease under compression.On the other hand, B 15 is almost unaffected by pressure.It should be noted that B 25 , B 35 , B 46 have similar negatives values and that they decrease under compression.

Table 5 .
Elastic moduli B, E and G, the B/G ratio, and the Poisson's ratio, ν, in the Hill approximation, at 0 GPa.

Table 5 .
Elastic moduli B, E and G, the B/G ratio, and the Poisson's ratio, ν, in the Hill approximation, at 0 GPa.

Table 6 .
Raman frequencies at ambient pressure, their pressure coefficients, and Grüneisen parameters (γ).The pressure dependence of frequency has been fitted with a second order polynomial: ω = ω 0 + �