Density functional theory study of single-molecule ferroelectricity in Preyssler-type polyoxometalates

Ferroelectricity (FE) is a characteristic of some electro-active materials that possess a robust spontaneous polarization, which can be switched under an external electric field.¹ Due to the nonlinear nature,^2,3 hysteretic behavior,^4,5 and catalytic properties,^6,7 ferroelectric materials (FEs) have been widely employed as capacitors,^8,9 sensors,^10,11 tunnel junctions,^12–14 catalysts,^15,16 and nonvolatile memory.^17–20 Up to now, many important FEs belong to the class of perovskites,^21–26 such as barium titanate (BaTiO₃)^27,28 and lead zirconate titanate.^29,30 However, FE has also been described in many other materials,^31–34 and the development of computing techniques and algorithms also drives tremendous efforts to look for higher-performance ferroelectric materials.^35–39 

Polyoxometalates (POMs) are nanoscale anionic metal-oxide clusters of d-block transition metals, which play a very active role in many fields such as homogeneous and heterogeneous catalyses,^40–43 electrocatalysis,^44,45 photocatalysis,^46–48 materials science,^49–51 molecular magnetism and spintronics,^52,53 and macromolecular crystallography⁵⁴ due to their structural diversity, electronic properties, modifiability, tunability, and redox properties.^55,56

In the vast family of POMs, the donut-shaped molecule consisting of five PO₄ tetrahedra surrounded by 30 WO₆ octahedra is termed as Preyssler-type POMs ([Mⁿ⁺P₅W₃₀O₁₁₀]^(15−n)−, abbreviated as [Mⁿ⁺@P₅W₃₀]^(15−n)−), and shown in Fig. 1.^57,58 Compared to the classical Keggin- and Dawson-type POMs, Preyssler anions display distinctive properties, among which the most notable are their (1) high negative charges (−12 to −15), (2) large oxidation potential, and (3) reversible transformation. The diameter of the internal cavity of Preyssler clusters approaches 5 Å, large enough to accommodate at least one metal cation (Mⁿ⁺).^59,60

So far, most studies addressed the (electro)catalytic activities and extending frameworks for functionalization of Preyssler clusters,^61–63 but more recently, Nishihara and co-workers reported a Preyssler cluster, [Tb³⁺@P₅W₃₀]¹²⁻, showing switchable electric polarization with high spontaneous polarization (6 µC/cm²) and coercive electric fields, which is the first discovery of FE in the family of POMs.⁶⁴ It was assumed that a Tb³⁺ cation encapsulated in the Preyssler cluster could be interconverted by thermal energy and electric field between two isoenergetic sites (Site-A and Site-A′ in Fig. 1), leading to the inversion of polarization. Although the bi-stability of the two isoenergetic sites provides the basic ingredients to explain FE, there are still some points that need further clarification for a complete description of the origin of FE in Preyssler clusters. In the first place, it is not clear whether there is a water molecule bound to the Tb³⁺ cation inside the cavity of the POM, which would certainly influence the interconversion of the cation between the two isoenergetic sites. X-ray and EXAFS studies of [Mⁿ⁺@P₅W₃₀]^(15−n)− (M = Na, Ca, Y, Ag, and Bi) have typically indicated the presence of a water molecule inside the cavity linked to the Mⁿ⁺ ion.^62,65,66 Moreover, due to the relatively high charge of the Preyssler anion, quite a few counter-cations are present in the crystal, which is closely associated with the redistribution of charges and may play a role in the polarization switching. While experimental studies have been addressed for Preyssler systems, a detailed analysis based on density functional theory (DFT) for the microscopic mechanism at the origin of this “single-molecule” ferroelectricity is still lacking due to the complexity of the system. This motivated us to close this gap.

Here, we report DFT calculations carried out to explore the effect of a water molecule inside the cavity and the origin and main contributions to the ferroelectric property experimentally observed in the [M³⁺@P₅W₃₀]¹²⁻ system with M = La, Gd, and Lu. Gd³⁺ has a similar ionic radius as Tb³⁺ but is computationally much simpler because of the half-filled 4f shell. DFT has, in general, important problems in treating ions with partially filled f-shells other than 4f⁷. La³⁺ (4f⁰) and Lu³⁺ (4f¹⁴) are studied for comparison, being the largest and smallest ions of the lanthanide series, respectively. Our results are arranged as follows. We first systematically examine the structures of the Preyssler anions [M³⁺@P₅W₃₀]¹²⁻, combined with the energy comparison between ferroelectric and anti-ferroelectric phases. We find that one water molecule linked to the internal cation is necessary for a ferroelectric ground state. Without an extra water molecule in the cavity, there is no indication that the metal will leave the center of the cavity and break the structural symmetry. Then, we explore the switching path connecting Site-A and Site-A′ of the encapsulated cation–water unit (M³⁺–H₂O). Combining the displacement of M³⁺ between the two isoenergetic sites with a 180° rotation of H₂O results in a barrier of 1.15 eV for Gd³⁺, in good agreement with the experimental barrier of 0.96 eV for Tb³⁺.⁶⁴ Finally, we calculate the polarization value to be 4.7 µC/cm² of the Preyssler crystal using the Berry phase method,⁶⁷ which is in fairly good agreement with the experimental value of about 6 µC/cm².⁶⁴ We also observe that the interconversion of M³⁺–H₂O between two isoenergetic sites induces a polarization of 3.4 µC/cm², and it can be considered as the main origin of FE. The polarization contribution of K⁺ counterions is about 1.2 µC/cm². It is certainly non-negligible while the contribution of the [P₅W₃₀]¹⁵⁻ framework is much smaller and can safely be ignored.

Geometry optimizations for the isolated molecule were performed with the ADF 2016 package⁶⁸ using the OPBE functional.^69,70 All-electron triple-ζ-quality basis sets with single polarization functions (TZP) were used for all atoms.^71,72 Relativistic effects were included via the zeroth-order regular approximation (ZORA).^73,74 The solvent effects were treated by the conductor-like screening model (COSMO) with a dielectric constant (ε) of 78 to simulate the aqueous environment, as implemented in ADF 2016.^75,76 Different point-group symmetries were applied considering the complexity of the Preyssler cluster. Open-shell electronic configurations were considered with the spin-unrestricted formalism in the Gd³⁺ case.⁷⁷ For the crystal structure, we started our calculations from the experimental crystallographic data of the Preyssler system with M³⁺ = Bi³⁺.⁶⁶ The structure is orthorhombic, space group Pnma with a = 28.6997 Å, b = 21.5091 Å, and c = 20.8737 Å, and has 114 symmetry unique centers. It is optimized until Hellmann–Feynman forces were smaller than 0.05 eV/Å. Kohn–Shan equations were solved using the projector-augmented-wave (PAW) method,^78,79 as implemented in VASP.^80,81 The standard corresponding pseudopotentials represent the La [Kr-4d]; Gd, Lu, W [Xe]; K (1s²-3s²); P (1s²-2p⁶); and O (1s²) core electrons. We have used the Perdew–Burke–Ernzerhof (PBE) approximations for the exchange–correlation functional.^82,83 The energy cutoff for the plane wave expansion was set to 400 eV in agreement with cutoffs used in other studies.^84–86 Geometry optimizations at the gamma point and a 2 × 2 × 2 grid were used for the calculation of electronic polarization. By establishing a suitable antiferroelectric phase, the Berry phase approach is used to evaluate the ferroelectric polarization.⁶⁷ 

The first issue that needs to be clarified is the possible presence of water molecules inside the cavity. In the experimental study, the crystals were dried before the polarization measurements.⁶⁴ This could imply that in addition to the interstitial water molecules, also the water molecule(s) on the inside of the Preyssler cage (if present at all) were removed from the sample. Experimental studies showed that in some [M@P₅W₃₀] systems, the encapsulated ion is bound to a water molecule⁶⁵ and located off-center, while in others, no sign of water was found.^84,87 In previous studies in our group of a series of [M@P₅W₃₀]^q− structures, it was established that only the inclusion of a water molecule inside the cavity leads to an optimized structure with the Na⁺ ion at the position measured in experiment.⁵⁹ Due to computational limitations, other systems (among which also the here studied compound with La³⁺) were optimized without water on the inside of the cage. Now that the computational limitations are no longer pressing, we have studied the geometry of [M@P₅W₃₀]¹²⁻ and [M–H₂O@P₅W₃₀]¹²⁻ (M = La³⁺, Gd³⁺, and Lu³⁺) only imposing the symmetry restrictions of the C_s spatial group. This ensures that all structural relaxations (off-site movement of the encapsulated cation, W–O–W bond length alternation, etc.) can take place when energetically favorable. Without the water molecule, all three systems converge to a geometry with the cation at the center of the cavity, while the incorporation of a water molecule pushes the cations away from it. The bond lengths are 2.381 Å for La³⁺–H₂O, 2.283 Å for Gd³⁺–H₂O, and 2.203 Å for Lu³⁺–H₂O, as depicted in Fig. 2 (right). Two isoenergetic positions are also found: Site-A in the upper half of the cavity and Site-A′ in the lower half. The position of Site-A (A′) is not exactly the same for the three ions and depends on the radius of the lanthanide, i.e., 1.798 Å for La³⁺, 1.743 Å for Gd³⁺, and 1.723 Å for Lu³⁺ above (or below) the center of the cavity. We roughly estimate the binding energy (ΔE_{binding_H2O}, kcal/mol) between the internal M³⁺ cation and a H₂O molecule and compare the modulus with the case where the internal cation is Na⁺. The structure of [Na⁺@P₅W₃₀]¹⁴⁻ has been experimentally confirmed that Na⁺ is combined with a H₂O molecule encapsulated in the cavity. The results are summarized in Fig. 3 in which the cases of M^q+ = Ca²⁺ and Bi³⁺ are also listed for comparison. The observation that the binding energy of the three lanthanide ions is larger than in the case of Na⁺ indicates that the water molecule can indeed reside inside the cavity.

Following the assumption that the ferroelectric property of the Preyssler cluster originates from the interconversion of the cation trapped in the cavity between the two isoenergetic sites (Fig. 1), the energies related to the process of M³⁺ moving from Site-A to Site-A′ along the vertical Ĉ₅ axis are calculated. The study is carried out as a series of restricted geometry optimizations, fixing the M³⁺ ion at different positions on the Ĉ₅ axis with steps of 0.5 Å. Every structure is re-optimized, keeping the position of M³⁺ frozen, under the C_5v symmetry without H₂O and under the C_s symmetry with one H₂O molecule bound to the cation. Figure 4 displays the schematic energy profiles associated with the interconversion process of M³⁺, where the solid curve in the top panel corresponds to the presence of H₂O inside the cavity, while the dashed curve in the bottom panel corresponds to the absence of H₂O. While the relative energies are different for the three ions, the shape is strictly the same in all cases. When M³⁺ is not connected to H₂O, the energy decreases by 0.28 eV for La³⁺, 0.34 eV for Gd³⁺, and 1.38 eV for Lu³⁺ as the cation moves from Site-A to the center position and then increases again by the same amount moving on to Site-A′. If the cation moves further outward beyond Site-A′, the energy of the system continues to rise, and it will reach a maximum when the ion passes through the window of the framework. Therefore, the energy profile of the whole process is a simple parabola with a single minimum representing as M³⁺ displaces along the Ĉ₅ axis. However, the situation changes dramatically when a H₂O molecule is linked to M³⁺. The energy increases by 0.57 eV for La³⁺, 0.24 eV for Gd³⁺, and 0.03 eV for Lu³⁺ from A to the center, reaching a maximum at the center of the cavity, and then decreases again when M³⁺ moves to Site-A′. We have also verified that the M³⁺ ion located at Site-A shifts to the center of the molecule when the structure is fully relaxed without symmetry constraints. This is important since POM frameworks present what is known as alternating bond distortions that are hampered with the presence of symmetry constraints.⁸⁸ The two completely different profiles suggest that the presence of water is necessary for having the cation off-centered giving rise to two isoenergetic sites of the cation, a necessary condition for (anti-)ferroelectric properties. So far, we have only studied the displacement of the M³⁺–H₂O unit from the center of the cavity to Site-A or A′, but as can be seen in the broken solid curve (the half green part and half blue part) of Fig. 4, the two potential energy surfaces cannot be connected as the orientation of the water molecule is opposite for the two cases.

Therefore, we explore the potential energy profiles for the displacement of the internal M³⁺–H₂O unit in more detail, especially the way the unit passes through the center. We study three possibilities for the unit to switch between the two isoenergetic points. First, the orientation of M³⁺ and H₂O (M³⁺ is up and H₂O is down) is assumed to remain constant during the interconversion, as depicted by the blue curve in Fig. 5. Starting from the minimum A, the energy first rises and then falls until reaching the second minimum when M³⁺ is located near the center of the cavity, which is 0.72 eV (in the Gd³⁺ case) higher than A. As M³⁺–H₂O continues to move down, the energy increases again and reaches a maximum when M³⁺ is about 1 Å below the center and H₂O is exactly located at the window, the narrowest part of Preyssler, suggesting a strong interaction between H₂O and the framework here. If the metal is pushed further downward to Site-A′, the M–H₂O bond breaks and water leaves the cavity. As explained above, the bare metal ion tends to occupy the center of the cavity. Considering the experimental observation that the ferroelectric polarization has the same numerical value with the opposite direction at two isoenergetic points, we expect that the structures at these points should have a symmetry operation relating to each other. Hence, this path is not consistent with the experimental fact, although the relatively low barrier makes it possible for Gd³⁺ and Lu³⁺ to move from the upper part of the cavity to the lower one at room temperature. The second hypothesis assumes that the displacement of M³⁺–H₂O is accompanied by the gradual rotation of H₂O and their relative orientation is reversed from up-down at Site-A to down-up at Site-A′, described by the red curve in Fig. 5. In this way, one can guarantee that the two structures at the isoenergetic sites are symmetrical to each other. The intermediate situation with M³⁺–H₂O lying in the σ_h plane (the center of the cavity, θ = 90°) turns out to be a (meta-)stable intermediate connected to the initial position (Site-A) by a transition state (θ = 39°), 0.6 Å from the σ_h plane, at slightly higher energy (0.19 eV for Gd³⁺). The overall energy barriers to go from Site-A to A′ are 1.61 eV, 1.15 eV, and 0.52 eV for La³⁺, Gd³⁺, and Lu³⁺, respectively. Among them, the barrier of 1.15 eV for Gd³⁺ matches the experimental barrier of 0.96 eV for Tb³⁺ quite well, in line with the similar radii of both ions. The final possibility that we explored is that a second H₂O molecule enters the Preyssler framework from the outside and pushes the M³⁺–H₂O unit down until the original H₂O molecule leaves the Preyssler framework (the far right panel of Fig. 5). However, the calculations show that inserting a water molecule requires more than 2.0 eV of energy. Therefore, we also have to overthrow the third hypothesis and mark the second scenario as the most probable mechanism for switching the dipole moment of the POM.

For a more detailed analysis of the ferroelectric properties, we need to go beyond the molecular approach that we have adopted so far. Therefore, we now turn our attention to the crystal structure and analyze the differences between the ferroelectric phase (FP) and the antiferroelectric phase (AFP), as shown in Fig. 6. The full relaxation of the ferroelectric phase (FP) (Fig. 6, left) maintains the orthorhombic nature with four Preyssler molecules per unit cell and reproduces quite accurately the lattice vectors of the x-ray structure of the Preyssler system with the Bi ion inside the cavity.⁶⁶ It also predicts the four cations (La³⁺ or Gd³⁺) to be off-center and arranged in the same direction (all Site-A or all Site-A′) when linked to a water molecule. The displacement of the metal ion from the center is similar to the one found in the molecular calculations. In order to estimate the polarization and to avoid the inclusion of quantum of polarization,⁶⁷ it is necessary to build a centrosymmetric antiferroelectric phase with inversion symmetry as well as a way to connect the FP with the AFP structure through atomic displacements. The AFP is constructed by moving the M³⁺–H₂O unit from Site-A to Site-A′ along the Ĉ₅ axis in the two Preyssler anions, accompanied with the rotation of 180° of the H₂O molecule so that the four encapsulated M³⁺–H₂O units are in a head-to-tail arrangement (head-to-head arrangements are expected to be higher in energy) with two units in Site-A and the other two in Site-A′, as shown in Fig. 6 (right, AFP). In order to ensure the inversion symmetry of the AFP structure, the two Preyssler frameworks and the surrounding 24 K⁺ counter-cations outside the frameworks have also been properly centrosymmetrized. For both structures, we have performed a full geometry relaxation (cell parameters and atomic positions), and the energy difference between the two is 1.67 eV per unit cell in favor of the FP structure. In line with the experimental conditions, we have removed 96 water molecules from the unit cell outside the POM, only leaving the four that are inside the cavity bonded to the metal ion.

To connect the FP and AFP, we have followed the interconversion path that is studied in the molecular calculations in which the metal displacement is accompanied by a gradual rotation of the H₂O molecule bonded to it. The parameter λ collects all the geometrical changes along this path and interpolates between the AFP (λ = 0) and the FP (λ = 1), increasing the angle of the M–H₂O bond and the Ĉ₅ rotation axis of the Preyssler anion in steps of 10°. The total polarization of the system can be split into two parts, one that arises from the Preyssler molecule with the M–H₂O unit on the inside and the other part that is caused by changes in the positions of the K⁺ counter-cations. The polarization value of the AFP structure is 0 due to the presence of an inversion point at the center of the unit cell by construction. Figure 7 shows that the polarization increases along the b-axis, which coincides with the direction of the Ĉ₅ axis of the isolated Preyssler cluster, resulting from the increase of asymmetry relative to the non-polar reference structure. The polarization along the a and c directions remains very close to zero for all λ values. The calculated polarization value at the FP is 4.7 µC/cm², only 1.3 µC/cm² different from the experimental value. To trace the origin of the polarization and clarify the role of K⁺ counterions, we consider the individual contribution of the Gd³⁺–H₂O unit, [P₅W₃₀]¹⁵⁻, and K⁺ cations to the total polarization by displacing only appropriate functional units while keeping the other units at the centrosymmetric atomic position. For instance, to calculate the contribution of the Gd³⁺–H₂O unit, only Gd³⁺–H₂O units are displaced from the AFP to FP structure while [P₅W₃₀]¹⁵⁻ and K⁺ cations are kept in the non-polar structure. The contributions of [P₅W₃₀]¹⁵⁻ and K⁺ cations are obtained in the same way. The contributions of the Gd³⁺–H₂O unit, [P₅W₃₀]¹⁵⁻, and K⁺ cations are 3.4 µC/cm², 0.5 µC/cm², and 1.2 µC/cm², respectively. The sum of the three individual polarizations is 5.1 µC/cm², approximately the same as the calculated total polarization of 4.7 µC/cm², proving the consistency of decomposition of the polarization into three separate contributions. According to the order of the contributions of 3.4 µC/cm² > 1.2 µC/cm² > 0.5 µC/cm², the interconversion of the Gd³⁺–H₂O unit along the Ĉ₅ axis between two isoenergetic sites dominates the generation of overall polarization and hence can be considered to be the main origin of the FE in the Preyssler cluster. The redistribution of the overall charge caused by the movement of Gd³⁺ relocates the K⁺ counterions, resulting in an additional contribution that cannot be ignored while the contribution of [P₅W₃₀]¹⁵⁻ is quite small.

In summary, the single-molecule ferroelectric property of the Preyssler cluster, [M³⁺@P₅W₃₀]¹²⁻ (M = La, Gd and Lu), is studied by DFT calculations. By examining the structure of the single anion, we find that the cation enclosed in the Preyssler framework needs to bind one water molecule to go off-center and create the possibility of a polar structure. When the metal cation occupies either Site-A or Site-A′, the inversion of the cluster is broken and a permanent dipole moment is created, which is the basis of ferroelectric properties. Combining the interconversion of M³⁺ between Site-A and Site-A′ with a 180° rotation of H₂O results in a barrier of 1.15 eV for Gd³⁺, in good agreement with the experimental barrier for Tb³⁺, leading to the switching of polarization. Moreover, we simulated the polarization of the Preyssler crystal using the Berry phase method to be 4.7 µC/cm², close to the experimentally measured value. The interconversion M³⁺–H₂O contributes more than 66% and the relocation of K⁺ counter cations contributes about 24%, while the framework [P₅W₃₀]^15- contributes less than 10% to the total polarization. Clearly the interconversion of M³⁺–H₂O is the main reason for the generation of the ferroelectric polarization, which agrees with the hypothesis proposed by Nishihara and co-workers⁴¹ that the ferroelectricity of the Preyssler POMs origins from the bi-stability of the two isoenergetic sites of M³⁺.

See the supplementary material for computed alternating long and short W–O bond lengths and binding energy of different internal cations (M^q+) with a H₂O molecule of the [M³⁺@P₅W₃₀]¹²⁻ clusters.

This project received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie (Grant Agreement No. 713679) and from the Universitat Rovira i Virgili (URV). Financial support has been provided by the Spanish Administration (Grant Nos. CTQ2017-87269-P and CTQ2017-83566-P) and the Generalitat de Catalunya (Grant No. 2017-SGR629). F.W. acknowledges the kind hospitality by the Institute of Functional Material Chemistry, Northeast Normal University (China) during the visiting period from November 5, 2018 to January 13, 2019. F.W. also thanks the kind hospitality by CNR-SPIN c/o Department of Chemical and Physical Science of University of L'Aquila (Italy) during the vising period from November 2, 2019 to December 14, 2019.

A dataset collection of computational results is available in the ioChem-BD repository and can be accessed via https://doi.org/10.19061/iochem-bd-2-49.⁸⁹