Structure of the hexagonal NaYF₄ phase from first-principles molecular dynamics

Abstract

The hexagonal phase of NaYF₄ is one of the most popular hosts for the synthesis of upconverting phosphors. The local structure of the lattice is known to have an impact on the optical properties of the doped NaYF₄, however despite extensive research being conducted in this field, the structure of the crystalline β-NaYF₄ phase is little understood. In this paper, the bulk β-NaYF₄ is investigated by means of Car–Parrinello Molecular Dynamics simulations. Three different space groups proposed in the literature are compared. The results show that models based on P6̄ and P6̄2m space groups converge to the same structure with nine-fold coordination sites only, while the model based on the P6₃/m space group is distinct and contains sodium atoms in six-fold coordination sites. Besides coordination numbers, analysis of the distances and average structure, the paper also presents a study of the dynamics of the lattice. Depending on the model, anisotropic thermal vibrations are observed and quantified, as well as oscillations of the sodium atoms between adjacent sites. A different scale of oscillation is observed, depending on the element (sodium or yttrium) and coordination site. In some cases, a reduced symmetry of the coordination shell is observed.

---

1 Introduction

Sodium yttrium fluoride (NaYF₄) doped with lanthanides has become the most promising and popular material exhibiting upconversion properties.¹ The term upconversion, in non-linear optics, refers to a process where a subsequent absorption of two photons (here, near-infrared or infrared) leads first, to a double excitation and afterwards to the emission of a single photon with higher energy (typically in the visible range). The popularity of NaYF₄ as a host lattice for upconverting materials is due to the high efficiency, chemical stability, sharp emission bandwidth of the resulting material, as well as large shifts that separate the absorption and emission frequency.² The latter one can be easily tuned by the choice of dopant, e.g. NaYF₄:Yb³⁺,Er³⁺ emits red and green light, while NaYF₄:Yb³⁺,Tm³⁺ emits blue.¹

From the practical point of view, these materials can be applied, in the form of nanocrystals (NC), mostly for labeling biomolecules with the purpose of both in vitro and in vivo detection.¹ Some of the recent applications include for example FRET biosensors based on doped NaYF₄ nanocrystals³ and cancer cell imaging.⁴ Many of these applications require prior functionalization of the surface of NC. In the case of the FRET biosensor this is necessary in order to bind the NC with the acceptor molecule, but in other situations it might be required to passivate the surface with a monolayer of molecules in order to increase the biocompatibility or solubility of the NaYF₄ nanocrystals.^5,6

Up to this point two important issues arise, both of them strongly related to the structure of the material. An accurate theoretical description of the upconversion mechanism is possible, applying multireference quantum chemical methods, possibly including electron correlation as in the case of Cs₂ZrCl₆:U⁴⁺.⁷ However, methods such as the complete active space perturbation theory (CASPT2) are computationally extremely demanding and the structure of the crystal phase should be well known beforehand. On the other hand, Density Functional Theory (DFT) has proved to be a valuable tool, when the structure of inorganic materials and their surface has to be determined; examples found in the literature include yttrium-containing materials.^8–10 The second issue is related to the surface of the material: in order to understand the binding mechanism of passivating agents such as PEG-phosphate⁶ the structure of the surface of NaYF₄ should be known, including various aspects, such as the surface reconstruction. The symmetry of sites occupied by the dopant in the NaYF₄ lattice influences the UC properties,¹¹ however the symmetry may be affected by the disorder and defects in the lattice, therefore it is interesting to study the dynamics of the lattice. The present paper aims at investigation of atomic motions in the NaYF₄ crystal by means of Car–Parrinello Molecular Dynamics (CPMD) method.

The phase transitions of NaYF₄ were studied by Thoma et al.¹² They have studied phase diagrams of NaF–MF₃ systems and found a pure hexagonal β-NaYF₄ phase with a melting point of 691 °C. The corresponding lattice constants of this phase were a = 5.95 and c = 3.52 Å. Burns has studied¹³ the structure of the hexagonal phase of NaNdF₄ and proposed that β-NaYF₄ conforms to the same space group, namely P6̄. In the P6̄ group, there is a nine-fold coordinated position occupied by Y³⁺, another nine-fold occupied position occupied by Y³⁺/Na⁺ (in the ratio 3 : 1) and a six-fold coordinated position occupied by Na⁺ and vacancies (ratio 1 : 1). On the other hand, Sobolev has suggested¹⁴ that β-NaYF₄ crystals share the same group as the mineral, gagarinite, i.e. the P6₃/m group. Recently, the same structure was considered by Krämer et al.¹⁵ The P6₃/m group differs from P6̄ by the fact, that the first two positions are symmetry-related and the atoms can be intermixed. Finally, Roy and Roy¹⁶ have proposed that the β-NaYF₄ structure should belong to the P6̄2m group. This structure was much later considered by Grzechnik et al.¹⁷ in the context of high-pressure phase transitions.

In this paper, we investigate the structure of the bulk β-NaYF₄ (hexagonal) phase. The initial models are build using P6̄, P6₃/m and P6̄2m space groups. These structures are then optimized and subjected to CPMD simulations. The resulting data are analysed in order to find the average structure and fluctuations in positions of atoms, as well as to determine the coordination of Na/Y by fluorine atoms.

2 Methods

The simulations have been performed at the level of the Density Functional Theory methodology (DFT). Due to the lack of empirical potentials, classical simulations of NaYF₄ cannot be performed and quantum-chemical methods have to be used. This has two implications: one can expect higher versatility and accuracy, but the methods are expensive and therefore the size of the model and simulation time are limited. The CPMD method^18,19 has been used, which has the advantage over ab initio Molecular Dynamics, that the wavefunction does not need to be reoptimized at every step; instead, classical equations of motion are applied to the wavefunction coefficients.

The crystal structure of NaYF₄ has a certain disorder and is made up of positions that contain mixed Na and Y atoms as well as vacancies. To allow for a very good sampling one should use larger supercell (containing various configurations of atoms). This however can be cumbersome at the DFT level due to the high cost of the calculations. The models used here contain 12 unit cells to account for the aforementioned disorder, nevertheless some caution is required when analysing the results. In order to test various space groups proposed in the literature, three models have been studied. They have been labeled according to their initial space group: P6̄, P6₃/m (gagarinite structure) and P6̄2m. The lattice constants of the cell have been set according to those given by Wang et al.,² namely a = 5.96 Å and c = 3.53 Å. The hexagonal unit cell has been repeated and then truncated in a way preserving the right periodicity, obtaining an orthorhombic cell with roughly equal edges of 11.920 × 10.323 × 10.590 Å. Such a cell contains 108 atoms (Na₁₈Y₁₈F₇₂), which, as already mentioned, corresponds to twelve unit cells containing Na_1.5Y_1.5F₆. The atoms in the supercell were distributed according to the space group. Whenever the position allows for mixing of Na/Y atoms, they were placed in such a way that the cell contained different arrangements of atoms stacked one on top of another (e.g. Na–Y–Na and Na–Na–Y). Fig. 1 shows the actual layer-by-layer placement of atoms in the considered models. These three models have been optimized in CPMD, with fixed lattice parameters. Additionally, in order to investigate the influence of the lattice constant, the P6̄2m model has been optimized using the CP2K program,²⁰ allowing for the lattice constants to be changed, but preserving the orthorhombic system. The optimizations were performed with two different sets of parameters (basis-set and external pressure), which led to cells slightly smaller and bigger that the one reported in the literature. The resulting structures will be denoted as P6̄2m small and P6̄2m big, respectively. Another simulation, using the P6̄ model was set-up in order to test the influence of the functional selection; as an alternative to PBE, the HCTH²¹ functional was used and the results are referred to as P6̄ HCTH. Next, an NVT simulation in CPMD was performed for each model. The optimized atomic positions were further relaxed during MD, therefore the simulation was continued until at least 10 ns of data was collected from the equilibrated system.

Fig. 1 Arrangement of atoms in the models. Orange color – yttrium, blue color – sodium atoms. Symmetry-related positions A, B, C are shown for the bottom layers only.

Calculations in CPMD (both, optimization and MD) have been performed using the Perdew–Burke–Erzenhof (PBE) functional.²² The core shells have been replaced with Troullier–Martins pseudopotentials²³ and the cut-off radius for the wavefunction was set to 80 Ry. The real-space mesh was set to 160 × 140 × 140. The calculations benefited from the Γ-point approximation, which is justified due to the non-metallic character of the system and the larger box used. The Ewald summation was extended to 3 cells in each direction. The geometry optimization was continued until the gradient reached 5 × 10⁻³ a.u.; which is sufficient to prepare the initial structure for MD simulation. After that, the wavefunction was optimized until the maximum gradient reached 1 × 10⁻⁵ a.u. Finally, the MD simulation was performed in the NVT ensemble, using the Nose–Hoover thermostat.^24,25 The thermostat was applied to the nuclear degrees of freedom only, with the temperature set to 300 K, characteristic frequency of 700 cm⁻¹ and a chain length of 3. The fictitious electron mass was set to 600 a.u. and the time step was set to 5 a.u. The trajectory was sampled every 16 steps, namely 1.935 fs.

The P6̄2m small and P6̄2m big models were optimized using the CP2K program. The lattice parameters a, b, c were optimized, while the straight angles, maintaining the orthorhombic point group, were preserved. The PBE functional was used, together with the Gaussian and Plane Wave (GPW) method.²⁶ In this method, the wavefunction is built using GTO's, while the electron density is described with plane waves. For plane waves, a multigrid approach was used with a cut-off of 280 Ry at the finest grid level. The core electrons were replaced by GTH pseudopotentials.²⁷

The P6̄2m small model was optimized with a double-ζ DZVP-MOL-OPT-SR-GTH basis-set²⁸ and external pressure tensor of 30 GPa. The obtained cell had a volume of 1271 ų and lattice parameters of a = 11.839, b = 10.225, c = 10.496 Å. The P6̄2m big was optimized with a smaller, single-ζ SZV-GTH basis-set and a default pressure value (10 MPa) which led to an increase of cell volume from 1303 to 1388 ų and an increase of the parameters to a = 12.206, b = 10.595, c = 10.730 Å. The final cells were optimized again using CPMD, but with the lattice parameters fixed.

The equilibrated part of the MD trajectory (ca. 10 ns in each model) was processed statistically as follows. The small drift of the model was removed by shifting each frame so that the center of mass was at (0, 0, 0). The average configuration of atomic position (x̄, ȳ, z̄) was determined, as well as the standard deviation from this average (for each direction, σ_x, σ_y, σ_z, separately). Four characteristic positions of sodium and yttrium were considered: A:Y, B:Y, B:Na and C:Na (shown in Fig. 1). The difference between the average position (x̄, ȳ, z̄) and ideal position (x_id, y_id, z_id) was calculated. The ideal position is the one that would result from the perfect placement of the atoms according to the space group. For each characteristic position, the Radial Distribution Function (RDF) was calculated, showing the distribution of Y–F or Na–F distances. Based on these RDF's, the coordination numbers were calculated (by integrating the RDF between 0 and the first minimum). Additionally, for selected Na/Y atoms, the distance from the coordinating fluorine atoms were plotted along the trajectory.

The methods used (functional, pseudopotentials etc.) were validated by calculating the equilibrium bond lengths of NaF and YF₃ molecules. A single molecule (NaF or YF₃) has been placed in a cubic box of 15 × 15 × 15 Å and its geometry has been optimized in CPMD using similar methodology as described above, applying two different functionals, PBE and HCTH. The equilibrium bond lengths have been compared with values calculated at various levels of theory (DFT, MP2 and CISD+Q) found in the literature or calculated using MOLCAS package.²⁹ In the case of the NaF molecule, the cc-pVQZ basis set were used. Alternatively, for both molecules, NaF and YF₃, a relativistic effective core potential was used and the valence electrons were described using the following basis sets: (8s7p6d)/[6s5p3d] for Y, (4s4p)/[2s2p] for Na and (4s5p1d)/[2s3p1d] for F (henceforth called RECP-Dolg).³⁰ In the MP2 calculations using all-electron basis sets (in NaF), two variants of the frozen core were used: three orbitals (1s of F and 1s2s of Na) or six orbitals (1s of F and 1s2s2p of Na)—designated FC3 and FC6, respectively.

3 Results

3.1 Validation of the methods

The isolated molecules NaF and YF₃ have been used here as benchmark systems, to validate the applied methodology against higher-level calculations. For the isolated YF₃ molecule, the plane-wave DFT using the PBE functional (as implemented in CPMD) predicts a planar structure with a Y–F bond length of 2.005 Å. This is in accordance with another DFT study by Perrin et al.,³¹ who found the value of 2.002 Å using the B3PW91 functional. MP2 and CISD+Q calculations by Solomonik and Marochko³² predict the Y–F bond length to be 2.014 and 2.012 Å, respectively, whereas our own MP2 calculations predict a value of 2.024 Å. The plane-wave DFT using the HCTH functional gives a slightly lower value of 1.996 Å. All data are collected in Table 1.

Table 1 Equilibrium bond lengths (in Å) calculated at various levels of theory

Method	Basis set	NaF	YF3	Ref.
PBE/CPMD	PW	1.969	2.005
HCTH/CPMD	PW	1.970	1.996
MP2/FC3	cc-pVQZ	1.947
MP2/FC6	cc-pVQZ	1.994
MP2	RECP-Dolg	1.976	2.024
MP2	RECP		2.014	32
CISD+Q	RECP		2.012	32
B3PW91	RECP		2.002	31

---

For the NaF molecule the plane-wave DFT with PBE and HCTH functionals predict the Na–F bond length to be 1.969 and 1.970 Å, respectively. The MP2/RECP-Dolg value is slightly higher, 1.976 Å. The MP2 value using the all-electron cc-pVQZ basis set yields bond lengths of 1.947 and 1.994 Å, depending on how many core orbitals are frozen. Clearly, including the 2p orbitals of sodium in the dynamic correlation has a significant impact on the bond length. The RECP-Dolg basis set leaves out only one valence electron to be treated explicitly, which could explain the overestimated bond length. The effective core potentials used in CPMD simulations on the other hand, replace only the 1s shell, leaving out 9 valence electrons.

The calculations performed for P6̄2m using literature (P6̄2m) and optimized (P6̄2m small and P6̄2m big) cell parameters, show very similar fluctuations (Table 3). The overall deviation (σ_x + σ_y + σ_z) for yttrium atoms in P6̄2m and P6̄2m small are virtually the same, while in P6̄2m big they increase by 0.01 Å. However the influence of the cell size is anisotropic: standard deviation of σ_z differs between P6̄ 2m small and P6̄2m big by 0.10 and 0.15 Å, for A:Y and B:Y respectively. The lowest fluctuations of sodium atoms are observed in P6̄2m (the literature parameters). Increasing and decreasing cell size leads to higher deviation, especially in the z direction. The average positions of atoms in P6̄2m and P6̄2m big are very close to the idealized structure and to each other. The only difference is in the x coordinate, where A:Y and B:Y positions deviate by 006–0.008 Å and 0.015–0.017 Å, respectively, from the idealized structure. The P6̄2m small model is characterised by a larger displacement from the idealized structure, with the most profound change (−0.028 Å) in ȳ − y_id of the A:Y site. The RDF's of the P6̄2m, P6̄2m small and P6̄2m big models show no influence of the cell size.

Table 2 shows the relative energy of the models and, as it could be anticipated, the optimized model P6̄2m big is characterised by the lowest energy. Table 3 also shows results of the simulation of P6̄, performed using the HCTH functional. The results are similar to those obtained with the PBE functional.

Table 2 Average potential energy of a unit cell (Na_1.5Y_1.5F₆) in each model

Model	Energy [eV]	σ [eV]
P6̄	0.000	0.030
P6̄2m	0.058	0.029
P6̄2m big	−0.103	0.029
P6̄2m small	0.016	0.029
P63/m	0.043	0.030

---

Table 3 Difference between mean (x̄) and ideal (x_id) position, standard deviation (σ) (in Å) and coordination numbers of the distinct positions

Position	# Atoms	x̄ − xid	ȳ − yid	z̄ − zid	σx	σy	σz	Coord. no.
P6̄
A:Y	12	0.003	0.002	−0.049	0.071	0.071	0.054	9.00
B:Y	6	−0.005	0.002	−0.044	0.060	0.065	0.070	9.00
B:Na	6	0.012	0.008	−0.045	0.104	0.109	0.132	9.00
C:Na	12	−0.001	−0.012	0.523	0.105	0.104	0.329	8.86
P6̄ HCTH
A:Y	12	−0.003	−0.001	0.050	0.071	0.073	0.059	9.00
B:Y	6	0.005	−0.001	0.042	0.062	0.068	0.068	9.00
B:Na	6	−0.014	−0.008	0.041	0.101	0.107	0.133	9.00
C:Na	12	0.000	0.013	−0.518	0.100	0.104	0.280	8.92
P6̄2m
A:Y	12	0.006	0.000	0.000	0.074	0.073	0.061	9.00
B:Y	6	−0.015	0.001	0.000	0.065	0.068	0.078	9.00
B:Na	18	0.000	0.000	0.000	0.100	0.100	0.117	8.99
P6̄2m small
A:Y	12	−0.013	−0.028	−0.012	0.070	0.080	0.059	9.00
B:Y	6	0.001	0.023	0.012	0.064	0.073	0.073	9.00
B:Na	18	−0.015	0.021	0.014	0.100	0.110	0.180	8.96
P6̄2m big
A:Y	12	0.008	0.001	0.000	0.077	0.075	0.069	9.00
B:Y	6	−0.017	−0.001	−0.001	0.071	0.073	0.088	8.98
B:Na	18	0.000	0.000	0.000	0.108	0.111	0.133	8.98
P63/m
B:Y	18	−0.002	−0.003	0.040	0.066	0.065	0.073	9.00
B:Na	6	−0.001	0.001	0.042	0.122	0.113	0.114	8.99
C:Na	12	0.001	0.009	−0.447	0.103	0.100	0.480	6.19

---

3.2 Comparison with the ideal bulk structure

Three models (P6̄ P6₃/m and P6̄2m) of the hexagonal phase of NaYF₄ have been built according to their space groups. These initial (ideal) structures were optimized and subjected to NVT molecular dynamics. Table 3 illustrates how much the average structure from the MD (x̄, ȳ, z̄) run differs from the initial (ideal) structure (x_id, y_id, z_id). The data are averaged overall atoms occupying a particular position. As can be seen, there is almost no difference in the position of P6̄2m – atoms fluctuate, but do not diverge from the initial position. In the case of P6̄, there is almost no change in A:Y, B:Y and B:Na positions. In the C:Na position, the x and y coordinates do not change significantly, but z changes by 0.52 Å, on average. The initial z coordinate of C:Na position was 1.18 Å (or 0.33 in fractional coordinates); the average position amounts to 1.77 Å (0.5 in fractional coordinates). In other words, this change indicates either the displacement of the Na atom from a six-fold coordination site to a nine-fold coordination site (typical for P6̄2m) or oscillation between the two six-fold coordination sites (the other one being a vacancy). Indeed, the coordination numbers for the P6̄ model indicate that atoms A:Y, B:Y and B:Na occupy nine-fold coordination sites. The atoms in the C:Na position have a coordination number slightly lower than 9, which suggests some oscillations in the direction of six-fold coordination sites. After the displacement of C:Na atoms to the nine-fold coordination site, the average P6̄ model becomes similar to P6̄2m. The superposition of these two averaged models is shown in Fig. 2. The gagarinite (P6₃/m) model has been built by swapping two yttrium atoms in the A:Y position with two sodium atoms in the B:Na position of the P6̄ model (hence, positions A and B become a single position designated as B). Such a small change has quite a significant impact on the behaviour of the structure. Comparing the average P6₃/m structure with an ideal one, there is no significant difference in B:Y and B:Na positions, however atoms in the C:Na position are shifted along the z axis by −0.45 Å. This again could mean a displacement to a nine-fold coordination site or the oscillation between two adjacent six-fold coordination sites, similarly to P6̄. As the coordination number 6.2 indicates (Table 3), it is the latter case. To summarize, P6̄ and P6̄2m models appear to represent the same structure, with all atoms in the nine-fold coordination site and are better described by the P6̄2m space group. The P6₃/m model is substantially different in the respect that the C:Na atoms occupy six-fold coordination sites.

Fig. 2 Superposition of the averaged structures P6̄ and P6̄2m. Colors: pink and yellow – Y, blue and magenta – Na, cyan and green – F (P6̄ and P6̄2m respectively).

3.3 Thermal vibrations

Further insight in the structure of the investigated models can be gained from the root-mean-square fluctuations (σ) of atoms calculated along the trajectories. To facilitate the analysis, fluctuations have been averaged over groups of atoms occupying the same crystallographic site in the cell (Table 3). These σ parameters are readily accessible from the trajectory and are related to the thermal factors obtained from X-ray scattering. For the yttrium atoms, for all models, σ varies in the range 0.05 to 0.08 Å, indicating small oscillations around the average position. For sodium atoms, in all the models, σ_x and σ_y amount to 0.10–0.11 Å, which indicates stronger vibrations. Most interestingly, the σ_z value is much larger, especially in the C:Na position of P6̄ and P6₃/m (0.33 and 0.48 Å, respectively). When this information is combined with the coordination numbers, it becomes clear that sodium atoms in the P6̄ model occupy nine-fold coordination numbers, but oscillate toward six-fold coordination sites. The C:Na atoms in P6₃/m occupy six-fold coordination sites, jumping back and forth between adjacent sites. The σ values can be conveniently viewed in Fig. 3, which shows them in the form of semi-transparent ellipsoids, which have been stretched to the size of 6σ_x × 6σ_y × 6σ_z. Assuming a normal distribution, they correspond to 99.7% probability of finding a particular atom in a trajectory frame. It can be seen that except for C:Na atoms, the vibrations are almost isotropic. The C:Na atoms form “channels” and ellipsoids are elongated in the z direction. The magnitude of the elongation is not uniform, however. In P6₃/m, two of the four C:Na “channels” have larger ellipsoids, indicating more freedom to displace atoms. The picture for the P6̄ model indicates moderate oscillation of the C:Na atoms along the z direction. Finally, in the P6̄2m model, vibrations along the z-axis are diminished. This may be attributed to insufficient sampling (in terms of the cell size and trajectory length) and to differences in the local symmetry of the coordination site.

Fig. 3 Standard deviation of atomic positions (σ). The ellipsoids span the dimensions 6σ_x, 6σ_y, 6σ_z (99.7% probability for a normal distribution).

3.4 Analysis of distances

Fig. 4 shows propagation of Na/Y–F distances of selected pairs. The panel of P6̄ A:Y shows a nine-fold coordination site of yttrium, typical not only for the P6̄ model, but also for P6₃/m and P6̄2m. The coordination site is highly symmetric, with all nine contacts of equal length. In the P6₃/m model, A and B coordination sites are unique, what introduces small asymmetry to some coordination sites (for example the one shown in panel P6₃/m B:Y). Panel P6̄ B:Na shows one of the nine-fold coordination sites occupied partially by Na and Y (here, occupied by Na). This site possesses lower symmetry, with three atoms in plane with Na having slightly closer contact than the remaining six fluorines. The P6̄ C:Na panel shows a coordination site occupied by Na or a vacancy. It is clear, that this is also a nine-fold coordination site, however the Na atom oscillates out-of-plane, which is manifested by the opposite changes in the blue and green line (atoms coordinating from top and bottom), while the red line (one of the three fluorines in plane with Na) is more fixed. In the P6₃/m model, sodium occupies the B coordination site together with yttrium, however the position is less rigid (compare panels P6₃/m B:Y and P6₃/m B:Na). The panel P6₃/m C:Na shows that the C coordination site in gagarinite is a six-fold coordination site, with Na atom jumping between adjacent sites. For about 5 ps the red and blue lines show equal distances, but after that the atom jumps to the vacant site and the green line becomes similar to the red one. Two bottom panels in Fig. 4 show that in the P6̄2m model sodium atoms in the B position occupy nine-fold coordination sites, but the magnitude of the out-of-plane motion depends on the site (larger in the bottom-left panel and smaller in the bottom-right panel). This is because in our model, the supercell spans three unit cells along the z axis and, as it happens, two “channels” occupied exclusively by Na are formed. For those “channels” oscillations are larger.

Fig. 4 Distance between Na/Y and F from the first coordination sphere, measured along the trajectory for selected pairs. In each case, three distances are shown: between Na/Y and F from the middle (red), bottom (blue) and top (green) triangle of the nine-fold coordination site.

Fig. 4 presents details of the coordination of selected atoms; an alternative view is presented in the form of RDF's, averaged over the whole trajectory and over all atoms in a particular position (Fig. 5 and 6). The RDF's for yttrium feature a sharp and narrow peak indicating the well-defined Y–F distance in the first coordination sphere. In contrast, RDF's for sodium show broader and shorter peaks, which indicate more variability of distances. The RDF's of P6̄ C:Na and P6̄2m B:Na share a characteristic broadening of the first peak, which may result from oscillations towards the adjacent six-fold coordination sites. In general, RDF's of P6̄ and P6̄2m look very similar. The first peak of the C:Na RDF in the P6₃/m model has no broadening at the base, since the atoms occupy the six-fold coordination sites.

Fig. 5 Radial distribution functions of Y–F and Na–F contacts in P6̄ and P6₃/m.

Fig. 6 Radial distribution functions of Y–F and Na–F contacts in P6̄2m, P6̄2m small and P6̄2m big.

4 Discussion

The NaYF₄ hexagonal crystals, being a useful host crystal for non-linear optics applications, pose a difficult system for computational studies. The crystal contains positions occupied by mixed Na and Y atoms and possibly also vacancies. A number of X-ray studies on this system was performed to date,^13–17,33 however they were inconclusive regarding the space group and structure of the material. The early papers by Roy and Roy¹⁶ and Sobolev et al.¹⁴ proposed the material to share the structure of a mineral, gagarinite (P6₃/m space group). Burns proposed that the hexagonal NaYF₄ phase has the P6̄ space group.¹³ The powder X-ray diffraction study by Krämer et al.¹⁵ yielded a P6̄, however, a later single-crystal study by Aebischer³³ produced a structure conforming to the P6₃/m space group. The discrepancy was attributed to the missing information in the powder diagram. Finally, a third structure, according to the P6̄2m space group, was proposed for the high-pressure structure of NaYF₄ by Grzechnik et al.¹⁷ The situation can be further complicated by doping with lanthanides (used in the upconversion application), which breaks the local symmetry of the lattice.¹¹

The three space groups have the following characteristics: in the P6̄ structure there is a nine-fold coordination site at (0, 0, 0) occupied solely by yttrium atoms (denoted here as A:Y). In P6̄ there is another nine-fold coordination site at (2/3, 1/3, 1/2) occupied by Y and Na in the 1 : 1 ratio (here denoted B:Y and B:Na). The third site in P6̄ (C:Na), located more or less at (1/3, 2/3, 1/3) and (1/3, 2/3, 2/3) is occupied by Na and vacancies in the ratio 1 : 1. In the P6₃/m structure, the sites A and B become unique and occupied by Y and Na atoms (denoted B:Y and B:Na here). In the P6̄2m structure, there is a nine-fold coordination site at (0, 0, 0) occupied solely by yttrium atoms, just like in P6̄ (denoted A:Y). Second a nine-fold coordination site, located at (1/3, 2/3, 1/2) and (2/3, 1/3, 1/2) is occupied by Y and Na with a 1 : 3 ratio (B:Y and B:Na).

The hexagonal phase of NaYF₄ was studied by means of DFT by Yao et al.^34,35 They performed geometry optimizations of periodic models with varying number of atoms in the unit cell and subsequent electronic structure analysis, which indicated that charge transfer from F to Y centers should occur during excitation.³⁴ This study was followed by DFT investigation of the hexagonal NaYF₄ phase doped with lanthanide ions, showing that the Hund's rule is obeyed and the lanthanide contraction with increasing atomic number is observed.³⁵ The CPMD method exploiting DFT formalism has been successfully used to study the bulk structure of yttria-stabilized zirconia,³⁶ adsorption on the surface of YVO₄ (ref. 9) or hydration of Y³⁺ ion.³⁷ A similar host lattice, namely LiYF₄ was studied using first-principles methods, however only static calculations (geometry optimization) were performed,^38,39 whereas here we focus on the dynamics of the lattice.

The comparison of the bond lengths in YF₃ and NaF using the methodology applied for crystal cell simulations shows that the structure of YF₃ is predicted correctly, while the Na–F bond might be slightly overestimated (by ca. 0.02 Å).

Based on the analysis of the trajectory, one can say that the P6̄ and P6̄2m structures are indistinguishable, at least at the studied temperature of 300 K and the assumed lattice constant, however comparison with the optimized models (P6̄2m small and P6̄2m big) shows that the results are insensitive to the cell size. Both models, P6̄ and P6̄2m converge to the same structure. They do differ in the magnitude of the vibrations of sodium atoms, but this can be attributed to the discrepancies in the initial positions of atoms: in the P6̄2m structure, sodium atoms do not have to form mixed Na-vacancy “channels”. In fact, two such “channels” are present in our P6̄2m due to the fact that the supercell spans only three unit cells along the z dimension, but in a real crystal they would have a limited length. In both models, the yttrium and sodium atoms occupy nine-fold coordination sites only. Surprisingly, the gagarinite structure (model P6₃/m) behaves quite differently from P6̄ despite only a small structural difference. In this structure, as predicted by Krämer et al., the six-fold coordination sites are occupied by sodium atoms and vacancies and thermal vibrations at 300 K are sufficient for the atom to jump from one site to the other.¹⁵

The analysis of the atomic fluctuations reveals a few interesting facts. First of all, the magnitude of the oscillation is larger for sodium atoms than for yttrium, regardless of the coordination site (see Fig. 3 and 4). For the sodium in the C position (mixed Na-vacancy sites), there is a clear anisotropy with larger oscillations along the z coordinate. These oscillations result probably from the presence of neighbouring six-fold coordination sites, which are empty.

5 Conclusions

Based on the presented results, it is difficult to decide on the right space group for the hexagonal NaYF₄ crystal. In fact, due to the high variability of the structure, one could suppose that an unequivocal assignment is not possible at all. Our MD study therefore, is more likely showing three possible “scenarios” – placement and behaviour of the atoms, depending on the local constitution of the lattice. Two of the models, P6̄ and P6̄2m seem to represent the same structure, with all Na/Y atoms in the nine-fold coordination sites. However, the existence of six-fold coordination sites is not completely excluded, as shown in the P6₃/m model. It seems that it is the second coordination sphere that decides between the six- or nine-fold coordination site of sodium. In the P6̄ model, the C:Na position possesses in its second coordination sphere three A:Y positions (therefore occupied by yttrium atoms only), while in the P6₃/m model, the C:Na position has in its second coordination sphere three B:Na/Y positions (mixed Na/Y atoms).

The sites occupied by yttrium are characterised by more fixed Y–F distances, while the same sites occupied by sodium are characterised by larger fluctuations of the Na–F distance. The A:Y site of P6̄ and P6̄2m are more symmetrical (equal Y–F distances), while the coordination in other sites, especially those occupied by sodium, are less symmetrical and characterised by stronger vibrations in the z-direction.

Acknowledgements

The work has been financed from the state grant no. 0543/IP3/2011/71 awarded by the Ministry of Science and Higher Education (Poland) within the project Iuventus Plus. This research was carried out with the support of the HPC Infrastructure for Grand Challenges of Science and Engineering Project, co-financed by the European Regional Development Fund under the Innovative Economy Operational Programme. The calculations have been performed at the Interdisciplinary Centre for Mathematical and Computational Modelling (ICM, Warsaw, Poland) and Wroclaw Centre for Networking and Supercomputing (WCSS, Wroclaw, Poland).

References

1. M. Haase and H. Schäfer, Angew. Chem., Int. Ed., 2011, 50, 5808–5829.
2. F. Wang, Y. Han, C. S. Lim, Y. Lu, J. Wang, J. Xu, H. Chen, C. Zhang, M. Hong and X. Liu, Nature, 2010, 463, 1061–1065.
3. D. Tu, L. Liu, Q. Ju, Y. Liu, H. Zhu, R. Li and X. Chen, Angew. Chem., Int. Ed., 2011, 50, 6306–6310.
4. M. Deng, Y. Ma, S. Huang, G. Hu and L. Wang, Nano Res., 2011, 4, 685–694.
5. F. Wang, D. K. Chatterjee, Z. Li, Y. Zhang, X. Fan and M. Wang, Nanotechnology, 2006, 17, 5786–5791.
6. J.-C. Boyer, M.-P. Manseau, J. I. Murray and F. C. J. M. van Veggel, Langmuir, 2010, 26, 1157–1164.
7. L. Seijo and Z. Barandiarán, J. Chem. Phys., 2003, 118, 5335–5346.
8. D. Schiffbauer, C. Wickleder, G. Meyer, M. Kirm, M. Stephan and P. C. Schmidt, Z. Anorg. Allg. Chem., 2005, 631, 3046–3052.
9. M. Oshikiri, M. Boero, A. Matsushita and J. Ye, J. Chem. Phys., 2009, 131, 034701.
10. J. Wen, L. Ning, C.-K. Duan, Y. Chen, Y. Zhang and M. Yin, J. Phys. Chem. C, 2012, 116, 20513–20521.
11. D. Tu, Y. Liu, H. Zhu, R. Li, L. Liu and X. Chen, Angew. Chem., Int. Ed., 2013, 52, 1128–1133.
12. R. E. Thoma, G. M. Herbert, H. Insley and C. F. Weaver, Inorg. Chem., 1963, 2, 1005–1012.
13. J. H. Burns, Inorg. Chem., 1965, 4, 881–886.
14. V. P. P. B. P. Sobolev and D. A. Mineev, Dokl. Akad. Nauk SSSR, 1963, 150, 791.
15. K. W. Krämer, D. Biner, G. Frei, H. U. Güdel, M. P. Hehlen and S. R. Lüthi, Chem. Mater., 2004, 16, 1244–1251.
16. D. M. Roy and R. Roy, J. Electrochem. Soc., 1964, 111, 421–429.
17. A. Grzechnik, P. Bouvier, M. Mezouar, M. D. Mathews, A. K. Tyagi and J. Köhler, J. Solid State Chem., 2002, 165, 159–164.
18. http://www.cpmd.org/, Copyright IBM Corp 1990-2008, Copyright MPI für Festkörperforschung Stuttgart 1997–2001.
19. R. Car and M. Parrinello, Phys. Rev. Lett., 1985, 55, 2471–2474.
20. J. Hutter, M. Iannuzzi, F. Schiffmann and J. VandeVondele, WIREs Comput. Mol. Sci., 2014, 4, 15–25.
21. F. A. Hamprecht, A. J. Cohen, D. J. Tozer and N. C. Handy, J. Chem. Phys., 1998, 109, 6264–6271.
22. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
23. N. Troullier and J. L. Martins, Phys. Rev. B: Condens. Matter, 1991, 43, 1993–2006.
24. S. Nosé, Mol. Phys., 1984, 52, 255–268.
25. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695–1697.
26. G. Lippert, J. Hutter and M. Parrinello, Mol. Phys., 1997, 92, 477–488.
27. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter, 1996, 54, 1703–1710.
28. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105.
29. F. Aquilante, L. De Vico, N. Ferré, G. Ghigo, P.-Å. Malmqvist, P. Neogrády, T. B. Pedersen, M. Pitonak, M. Reiher, B. O. Roos, L. Serrano-Andrés, M. Urban, V. Veryazov and R. Lindh, J. Comput. Chem., 2010, 31, 224–247.
30. D. Andrae, U. Häußermann, M. Dolg, H. Stoll and H. Preuß, Theor. Chim. Acta, 1990, 77, 123–141.
31. L. Perrin, L. Maron and O. Eisenstein, Faraday Discuss., 2003, 124, 25–39.
32. V. G. Solomonik and O. Y. Marochko, J. Struct. Chem., 2000, 41, 725–732.
33. A. Aebischer, M. Hostettler, J. Hauser, K. Krämer, T. Weber, H. U. Güdel and H.-B. Bürgi, Angew. Chem., Int. Ed., 2006, 45, 2802–2806.
34. G. Yao, M. T. Berry, P. S. May and D. S. Kilin, Int. J. Quantum Chem., 2012, 112, 3889–3895.
35. G. Yao, M. T. Berry, P. S. May and D. S. Kilin, J. Phys. Chem. C, 2013, 117, 17177–17185.
36. G. Stapper, M. Bernasconi, N. Nicoloso and M. Parrinello, Phys. Rev. B: Condens. Matter, 1999, 59, 797–810.
37. T. Ikeda, M. Hirata and T. Kimura, J. Chem. Phys., 2005, 122, 024510.
38. W. Y. Ching, Y.-N. Xu and B. K. Brickeen, Phys. Rev. B: Condens. Matter, 2001, 63, 115101.
39. J. Yin, Q. Zhang, T. Liu, X. Guo, M. Song, X. Wang and H. Zhang, Phys. B, 2009, 404, 1053–1057.

---