First-Principles Study of Structure and Magnetism in Copper(II)-Containing Hybrid Perovskites

Abstract

We report a first-principles study of hybrid organic–inorganic perovskites with formula [A]Cu(H₂POO)₃ (A = triazolium (Trz) and guanidinium (Gua), and H₂POO⁻ = hypophosphite), and [HIm]Cu(HCO₂)₃ (HIm = imidazolium cation, ${\mathrm{HCO}}_2^{-}$ = formate). The triazolium hypophosphite and the formate have been suggested as possible ferroelectrics. We study the fully relaxed structures with different magnetic orderings and possible phonon instabilities. For the [Trz]Cu hypophosphite, the Trz cation is shown to induce large octahedral distortions due to the Jahn-Teller effect, with Cu-O long-bond ordering along two perpendicular directions, which is correlated with antiferromagnetic ordering and strongly one-dimensional. We find that the structure is dynamically stable with respect to zone-center distortions, but instabilities appear along high symmetry lines in the Brillouin zone. On the other hand, for the [HIm]Cu formate, large octahedral distortions are found, with large Cu-O bonds present in half of the octahedra, in this case along a single direction, and correspondingly, the magnetism is almost two-dimensional.

1. Introduction

Metal-organic frameworks (MOFs), where metal centers or clusters are linked by organic ligands, have received intense research interest recently. The denser MOFs, in particular with ABX₃ perovskite topology, where the X⁻ and/or A-site of simple perovskites are substituted by polyatomic ions, have shown fascinating ferroelectric, magnetic, and multiferroic properties [1,2,3], similar to their inorganic counterparts [4,5]. Theoretical studies of this family of compounds using density functional theory calculations [6] have matched well with experimental results, which highlights the predictive power of computation [7,8].

Recently, hybrid organic-inorganic perovskites based on the hypophosphite ligand have been synthesized [9]. Both formate and hypophosphite ligands produce intermediate-size cubic interstices, allowing for a choice of several different amine cations on the A site.

A recent group-theoretical analysis [10] showed that the number of possible distortions is greatly enhanced in molecular perovskites due to the additional kinds of symmetry-breaking degrees of freedom allowed, including unconventional octahedral rotations [11,12,13], columnar shifts [11,14], and multipolar order of the organic cations [15]. This allows more flexibility for different functionalities. The symmetry analysis [10] of the distortions present in synthesized Mn compounds, triazolium manganese hypophosphite [9] or imidazolium manganese formate [16], has suggested that the respective Cu analogues could be ferroelectric, due to the propensity of the Cu octahedra to show Jahn–Teller distortions in combination with the distortion modes available for hybrid perovskites already present in the Mn compounds. This motivated us to study Cu-based organic-inorganic perovskites using first-principles density functional theory calculations. For the hypophosphites, besides the suggested triazolium hypophosphite, we have also considered guanidinium as an organic cation. The structures considered in this work are shown in Figure 1.

We find that the different Cu compounds have both 1D and 2D magnetism, with small couplings in the spin exchange paths which involve the long Cu-O bonds, but we do not find significant structural instabilities with respect to the centrosymmetric group of the Mn analogues.

2. Methods

We have used the Vienna Ab-initio Simulation Package, VASP (VASP Software GmbH, Vienna, Austria) [17] for density functional calculations, which uses a plane wave basis set and the projector augmented-wave (PAW) method [18]. We considered the GGA-PBE (generalized gradient approximation, Perdew-Burke-Erzernhof) exchange-correlation functional [19] and added an effective Coulomb repulsion U to Cu 3d orbitals, according to the simplified approach of Dudarev et al. [20]. For Cu oxides, values of U = 4 eV may be most appropriate, but it is difficult to choose an appropriate U value that correctly describes different properties [21]. Here, we studied the variation of U from 0 to 5 eV. For the PBE (U = 0) cases, we also applied the DFT-D3 (BJ) method [22,23] to take van der Waals corrections into consideration. The following settings were used, except when stated otherwise: a 3 × 3 × 2 k-points mesh for the hypophosphites and a 2 × 2 × 1 mesh for the formate; the energy cutoff for the plane waves was 600 eV for the [Trz]Cu(H₂POO)₃ and 520 eV for the other two compounds; the full relaxations (atomic positions, cell shape, and volume) did not use symmetry, and stopped when the maximum residual force components were less than 0.005 eV/A. The FINDSYM program (Brigham Young University, Provo, UT, USA) [24,25] was used to find the symmetry of the optimized structure, the VESTA software for analysis and structural representation [26], and PHONOPY to post-process VASP results and calculate phonons.

3. Triazolium Cu Hypophosphite

3.1. Magnetism

We started from the structural data for [Trz]Mn hypophosphite given by Wu et al. [9], wherein the compound adopts the space group P2₁/c and subsequently replaces the Mn atoms by Cu. We performed a series of full structural optimizations with the ferromagnetic and the three antiferromagnetic orders obtained by the different combinations of up/down spins in the unit cell. These can be called A, C, and G-type antiferromagnetism (AFM) [27], considering that the Cu framework is cubic-like (even if the framework is distorted from ideal cubes by columnar shifts).

Figure 2 shows the total energies obtained for the optimized structures as a function of U. There is a clear separation between two pairs of states; the A and G-type AFM are always lowest in energy, while the C-type AFM and FM are highest. This separation decreases with U, as expected from the increased localization of states, but it is always much larger than between each pair. These degeneracies between states with different (F/AF) in-plane interactions suggest that the in-plane exchange interaction is very weak. On the other hand, the states where the out-of-plane magnetic ordering switches from F to AF are separated by large energy differences; thus, the out-of-plane exchange must be much stronger (particularly for low U values and AFM). The A-AFM state is slightly lower in energy than the G-AFM state, with the difference slightly larger for lower U values, suggesting a small in-plane ferromagnetic (FM) coupling, while the C and F states remain almost degenerate. Thus, this compound presents a quasi-1D magnetism, particularly for lower U values, with very weakly coupled 1D antiferromagnetic chains. This unusual property has previously been shown in another Cu metal–organic framework perovskite [28], dimethylammonium copper formate ([DMA]Cu[HCO₂]3), and it could be understood on the basis of Jahn–Teller distortions in the CuO₆ octahedra. We will show in the following structural analysis that strong octahedral distortions are also present in our compound. Indications of low-dimensional magnetism were also seen in guanidinium Cu formate [29], which is consistent with first-principles calculations that have shown a dominant 1D interaction [7]. We also used the DFT-D3(BJ) method to check the influence of dispersion corrections, parametrized for GGA (U = 0), shown in the rightmost set of points of Figure 2. The FM/AFM out-of-plane difference corresponds to intermediate values of U, while the A/G and particularly C/F states are now more separated, indicating higher in-plane exchange, but dominant 1D magnetism is still seen.

The calculated Cu moments are in the range between 0.6 and 0.8 μ_B, and the total moment per cell is 1 μ_B/Cu for the FM state. The values are consistent with the Cu²⁺ S = 1/2 Jahn-Teller state. Considering the cubic-like framework of Cu ions, we can use the same model as in Ref. [28] to map the density functional theory (DFT)total energies to a Heisenberg model and obtain three exchange interactions between spins: intra-chain nearest-neighbor (J₁), interchain in-plane (J₂), and interchain out-of-plane (J₃). The interactions are shown in

Table 1. Spin exchange interactions (in K) for H = 1/2∑_ijJ_ij S_i S_j, with S = 1/2, calculated with GGA, GGA+U, GGA+D3(BJ), and comparison with another quasi-1D Cu MOF.

	U = 0 eV	D3	U = 5 eV	[DMA]Cu(HCO2)3 1
J1	780	533	141	214
J2	−18.6	−36	−6.8	−8.8
J3	7.8	5.9	3.1	2.7

¹ U = 5 eV [27].

. (For comparison with the previous work, we follow the convention where AFM interactions are positive.) As it can be seen by comparison of the two rightmost columns, the exchange interactions in our compound are qualitatively similar to [DMA]Cu[HCO₂]₃ for the same U value, with an AFM intrachain larger by far than the interchain interactions J₂ + J₃ (overall weakly ferromagnetic). Therefore, similar experimental signatures of 1D magnetism should be found. The AFM intrachain spin exchange J₁ is slightly weaker in our case, but so is the interchain exchange (J₂ + J₃), such that the dimensionless ratio J₁/(J₂ + J₃) is −38 for [Trz]Cu(H₂POO)₃ and −35 for [DMA]Cu(HCO₂)₃.

3.2. Structure

Concerning the structure, while the Mn analogue at room temperature has a measured volume of 277.22 ų/Mn, for the Cu compound, our calculations (after full relaxation and without use of symmetry) result in 279.29 to 285.91 ų/Cu, depending on U_Cu, with smaller volumes for larger U. For G/A-AFM orders, the changes in volume are similar and come mostly from the change of U. There is a small expansion in any case, as could be expected on the basis of a slightly larger Cu²⁺ ionic radius relative to Mn²⁺ [30]. Concerning the transition metal–oxygen bonds, in the synthesized Mn compound, the Mn-O bonds vary between 2.11 and 2.24 Å, with an average bond length of 2.18 Å, while for the Cu compound, our predictions for U = 0 (U = 5) show a much more distorted octahedron, consisting of a Cu-O square with “equatorial” bonds 1.95–1.99 (1.96–1.99) Å, while the other two “axial” distances are much larger, 2.62 (2.42) and 2.94 (2.76) Å. These distances are almost the same for A and G-type orders. The strong AFM coupling between ab-planes can be related to the Cu-O bonds along the c direction connecting planes, which are all the smaller equatorial ones, while in the ab-plane, half of the bonds are the long axial bonds. The very large Cu-O 2.94 (2.76) Å distance in particular, alternating direction in half of the in-plane framework edges, must considerably weaken the Cu-Cu super-exchange path through the O atoms, disrupting the formation of long-range magnetic order in the plane. For increasing U values, the two axial distances decrease, so the distortion of the Cu-O octahedra is not so large; nevertheless, the in-plane magnetic coupling is lower due to the enhanced localization.

The comparison of structural parameters between GGA with and without D3 corrections and U is shown in

Table 2. Structural parameters a, b, c (Å), β (°), Volume V (ų/fu), and Cu–O equatorial and axial distances (Å) for different approximations: GGA, GGA+D3(BJ), and GGA+U (U_Cu = 5 eV).

Structural Parameter	GGA	D3	U = 5 eV
a	9.88	9.58	9.86
b	9.18	8.54	9.04
c	12.68	12.66	12.73
β	99.47	99.14	99.30
V	285.91	255.45	279.29
Cu-O eq.	1.95–1.99	1.95–1.99	1.96–1.99
Cu-O ax.	2.56 & 2.95	2.39 & 2.63	2.49 & 2.76

. The D3 corrections cause a large decrease in volume (10.6%) corresponding to the decrease of in-plane lattice parameters a and b by 3.0 and 7.0%. The effect of U = 5 eV is much smaller, and the small decrease in volume is mostly caused by a decrease in b by 1.5%. Looking at the Cu-O octahedral distances, the equatorial bonds remain in the range 1.95–1.99 Å, and the difference caused by the dispersion corrections is smaller axial distances, now 2.39 and 2.63 Å, compared with 2.62 and 2.94 Å with GGA, respectively. In this case, GGA+U also produces a considerable decrease in the axial bonds to 2.49 and 2.76 Å, although it is not as great as D3. Thus, the effect of both corrections is a less distorted octahedron.

We looked at the symmetry of the optimized structures, and, except for small (below 0.01 Å) tolerances in atomic positions and lattice parameters, the structure remains with the initial space group, P2₁/c, which is centrosymmetric. Using the structure found at U = 0 with the A-AFM order, we searched for possible dynamical instabilities using the VASP software by calculating the Hessian matrix (with central differences and a step size of 0.015 Å) and corresponding zone-center (Γ point) vibrational frequencies. The 297 (3N−3, where N = 100 is the number of atoms) non-translational modes are all real, showing that the P2₁/c structure is stable with respect to zone-center distortions, and thus, this compound will likely not be a proper ferroelectric. For example, the five lowest mode frequencies are 34, 46, 47, 48, and 51 cm⁻¹. However, this does not exclude other possibilities such as hybrid improper ferroelectricity [31], where the polarization would be induced by other modes. Therefore, we also calculated phonons in other regions of the reciprocal space, using density functional perturbation theory to calculate the force constants [32] and then post-processing using the PHONOPY code [33] to calculate the phonon dispersion in high symmetry lines of the Brillouin zone. The resulting phonon band structure plot is shown in Figure 3. Instabilities are seen in many lines in the Brillouin zone, which are shown as negative frequencies in the plot, but with small magnitudes (≈33 cm⁻¹). Although a more extensive check is in order here, the calculations are already too computationally expensive using the original cell.

There are imaginary frequencies (negative in the plot) in all the studied lines. In addition to the three lowest (acoustic) modes, two more modes (optic) get imaginary frequencies near the E point (A–E and E–Z directions), and near C₂ (Z–C₂ and C₂–Y₂ directions) there is one more, with six total imaginary modes.

To further explore these instabilities, we have looked at the eigenvectors of the phonons in some high-symmetry points (Z, D, E, C₂), and constructed modulated supercell structures with the displacements given by weighted eigenvectors. Then, these structures were optimized, as before, and the symmetry and energies of the optimized structures were analyzed. For the Z point, the modulated structure, as well as the optimized final structure, have P1 symmetry, although it remains P2₁/c within tolerances of 10⁻³ Å in lattice and 10⁻² Å in atomic positions or larger, which shows that the distortions are small. The energy lowering with respect to the initial P2₁/c structure is only a 0.30 meV/formula unit. For the D and E points, the initial modulated and optimized structures have $\mathrm{P}\overline{1}$ symmetry (for D, the structure remains P2₁/c symmetry for tolerances 10⁻³ Å in lattice and 10⁻² Å in atomic positions or larger; for E, the structure remains P2₁/c symmetry for tolerances 10⁻⁴ Å in lattice and 10⁻² Å in atomic positions or larger), and the energies are lower by 1.00 and 0.82 meV/Cu, respectively. These energy and structural differences are indeed small but due to the computational cost of the calculations, we did not explore this point further.

4. Guanidinium Cu Hypophosphite

We studied another compound of the Cu hypophosphite family with the Guanidinium cation (C(NH₂)₃^±) at the A site. For this system, we also started with the structure from the Mn analogue synthesized by Wu et al. [9]. In this case, we used the low temperature $P\overline{1}$ structure [9]. Due to the low symmetry, the only possible symmetry-lowering distortion is the break of inversion symmetry. We considered the FM and the three AFM states supported by the four Cu atoms in the triclinic unit cell. In this case, the + + − − state (where + (−) refers to the sign of projection of the spin moment on a given Cu atom along the quantization axis, considered as a z-axis, and the order corresponds to Cu atoms with z coordinates approximately 0.25, 0.75, 0.5, and 0), corresponds to the G-AFM alternating state for the pseudo-cubic framework, while the + − + − and − + + − cases correspond to double-layered AFM states for ferromagnetic (001) triclinic planes (corresponding to (111) pseudo-cubic planes), as shown in Figure 4. Furthermore, the initial calculations in some cases converged to a ferrimagnetic state of the type + + − +, suggesting that this is also a competing state; therefore, we also included it in our analysis.

Figure 5 shows the total energies obtained for the optimized structures of these five states as a function of U, where FM is taken as the zero-energy state. Indeed, for low values of U (<4 eV), the ferrimagnetic state (+ + − +) turns out to be the lowest energy state among those considered. In the range U > 3 eV, which may be most appropriate, the G-AFM state becomes the ground state.

5. Imidazolium Cu Formate

5.1. Magnetism

We started from the Mn analogue structure, which was also previously synthesized and characterized [16]. The low-temperature ordered structure has the same space group (P2₁/c) as [Trz]Mn(H₂PO₂)₃, but the HIm cations order such that there are now two crystallographically inequivalent cations and metal centers in 4e Wyckoff positions, and thus, the cell is approximately doubled with respect to the previous hypophosphites.

Figure 6 shows the energies for the three simple AFM orders, relative to FM, as a function of U. In this case, due to the larger cell, other orders could have been investigated, but we will focus on the G/C/A/F orders considering c as the out-of-plane axis. For example, the A-AFM order consists of ferromagnetic ab-planes with antiferromagnetic stacking along the c axis. In this case, the C/G states are much lower in energy than A/F; both pairs near degenerate, meaning that the out-of-plane coupling must be much weaker than the in-plane coupling, and there is a quasi-two-dimensional antiferromagnetism, particularly for low values of U.

We can use the same model as before to calculate the exchange interactions: J₂ is the in-plane interaction, while J₁ is the vertical out-of-plane interaction, and J₃ is the diagonal one.

Table 3. Spin exchange interactions (in K) for H = 1/2∑_ijJ_ij S_i S_j, with S = 1/2, calculated with GGA and GGA+U.

	U = 0 eV	U = 5 eV
J1 (out-of-plane)	−11.7	−1.5
J2 (in-plane)	+307.3	+77.9
J3 (out-of-plane)	2.5	−1.9

shows the calculated exchange interactions. As seen for [Trz]Cu(H₂PO₂)₃, the main exchange interaction decreases with U. In the dimethylammonium Cu formate case, the DFT values were much closer to experiment for the higher U values calculations [28], which indicates that for this Cu compound, the U = 5 eV results, with overall ≈75 K AF interaction, are also more reasonable than U = 0 eV. In this case, the degree of 2D magnetism can be estimated by J₂/(J₁ + J₃) = −22 for U = 5 eV, i.e., the AFM in-plane interaction is much larger than the out-of-plane ones (overall weakly ferromagnetic). Note that there are two times more octahedral in-plane than out-of-plane bonds, so the overall exchange interactions, here dominated by J₂ (in-plane), are of similar magnitude to the ones in the hypophosphite (where the dominant out-of-plane J₁ = 780 K for U = 0 eV and 141 K for U = 5 eV, approximately double the in-plane interactions here), and to the ones in the previously studied DMA Cu formate (J₁ = 214 K for U = 5 eV) [28].

5.2. Structure

In this case, there are two inequivalent octahedra. Half of the octahedra have four bonds with 2 Å, one with 2.5 Å, and one with 2.9 Å, while the other half have four bonds with 2 Å, one with 2.4 Å, and one with 2.6 Å. As in the [Trz] Cu hypophosphite, long 2.9 Å bonds should reduce the magnetic interactions. In the [Him]Cu formate, these bonds are only along the c out-of-plane direction, decreasing the coupling in this direction.

Using the FINDSYM program to check the symmetry group of the optimized structures, we found the P2₁ polar space group for the ferromagnetic calculation, U = 4 eV for small tolerances (10⁻³ Å for lattice and position parameters), while increasing the positions tolerance (to 10⁻² Å), the centrosymmetric P2₁/c space group is found instead. Nevertheless, this indicated the possible stability of the polar phase. In order to quantify a possible energy lowering, we performed relaxations with the symmetry constraints of both P2₁ and P2₁/c space groups (using the initial structures found with different tolerances) to well converged values for the more stable magnetic orders (C and G cases), with U = 4 eV. However, after converging the calculations to more stringent parameters (energy cutoff of 600 eV and forces threshold of 2 meV/fu), the energy difference is found to be 0.03 meV/fu for the C magnetic order, meaning that P2₁ and P2₁/c are almost degenerate in energy, and it is also almost vanishing for the G order, −0.006 meV/fu. In this case, we could not try to calculate phonons, due to the larger size of the system.

6. Conclusions

Starting from the Mn analogues, we have explored the possible compounds triazolium/guanidinium Cu hypophosphites and imidazolium Cu formate. We could not find significant energy-lowering distortions with respect to the P2₁/c space group of the Mn analogues, which suggests these compounds are not ferroelectric. We have shown that [Trz]Cu hypophosphite and [HIm]Cu formate exhibit strong, almost ideal 1D and 2D magnetism, respectively, highlighting the role of strong octahedral distortions in combination with different periodic orderings of the organic cations. We hope that our study motivates the synthesis and characterization of these compounds.