Van Vleck analysis of angularly distorted octahedra using VanVleckCalculator

A method and associated Python script, VanVleckCalculator, are described for parameterizing octahedral shear and first-order Jahn–Teller distortions in crystal structures.


Introduction
The van Vleck distortion modes (Van Vleck, 1939) modes describe all possible displacements of octahedrally-coordinated ligands about a core atom.They are particularly useful in the context of the Jahn-Teller effect (Jahn & Teller, 1937), which in general occurs when a high-symmetry coordination is destabilised with respect to a deviation to lower symmetry as a consequence of electronic degeneracy.The Jahn-Teller effect distorts the crystal structure via the Jahn-Teller distortion.While the Jahn-Teller distortion is not unique to octahedra in bulk crystalline materials, it is in octahedra that it was first observed experimentally (Bleaney & Bowers, 1952), and it is in materials with Jahn-Teller-distorted octahedra that colossal magnetoresistance (Millis et al., 1996) and high-temperature superconductivity (Fil et al., 1992;Keller et al., 2008) were discovered.
A transition metal (TM) cation in an octahedral configuration will have its d orbitals split into three t 2g orbitals 1 at lower energy and two e g orbitals at higher energy.It will have a number, n, of electrons in these d orbitals (hereafter described as d n ).
For certain values of n and, where applicable, certain low-or high-spin characters 2 , there will exist multiple orbitals that could be occupied by an electron or an electron hole with equal energy.This degeneracy is destabilising, resulting in the most stable configuration of atomic sites being one in which the ligands distort from their highsymmetry positions in order to rearrange the orbitals into a non-degenerate system with minimised energy.This is shown for a low-spin d 7 TM cation (such as Ni 3+ or Co 2+ ) in Figure 1, though such distortions may occur for any value of n in d n where there is a degenerate occupancy.The stabilisation energy due to the Jahn-Teller effect 1 In this paper, we use the notation that lower case symmetry descriptors (such as eg or t2g) refer to orbitals with this symmetry, and upper case descriptors (such as Eg or T2g) refer to the symmetry more generally. 2In the low-spin case, t2g orbitals fill fully before eg orbitals gain electrons; in the high-spin case, once the t2g orbitals are singly-occupied, the next two electrons will populate the eg orbitals.
IUCr macros version 2.1.10:2016/01/28 is larger for e g degeneracy than t 2g degeneracy, and so the effect is prominent to higher temperatures, and hence more widely-studied, in JT-active materials with e g degeneracy (Castillo-Martínez et al., 2011).In the literature, various techniques for parameterising the Jahn-Teller distortion are used.An often-used example (Kimber, 2012;Lawler et al., 2021;Nagle-Cocco et al., 2022;Genreith-Schriever et al., 2023) is the bond length distortion index, defined by Baur (1974), as: IUCr macros version 2.1.10:2016/01/28 where l i is the distance between the core ion and the ith coordinated ion, and l av is the average of all the distances between the core ion and coordinated ions.
A similar parameter (Shirako et al., 2012;Sarkar et al., 2018;Nagle-Cocco et al., 2022) is the effective coordination number, which for an octahedron deviates from 6 only when there is bond length distortion, defined by Hoppe (1979) as: where l ′ av is a modified average distance defined as: Finally, a third parameter used to quantify the Jahn-Teller distortion (Schofield et al., 1997;Kyono et al., 2015;Mikheykin et al., 2015) is the quadratic elongation, < λ >, defined by Robinson et al. (1971) as: where l 0 is the centre-to-vertex distance of a regular polyhedron of the same volume.
More recently, an alternative approach to modelling polyhedral distortion has been described (Cumby & Attfield, 2017), involving fitting an ellipsoid to the positions of the ligands around a coordination polyhedron, calculating the three principal axes of the ellipsoid, R 1 , R 2 , and R 3 , where R 1 ≤ R 2 ≤ R 3 , and using the variance of these three radii as a metric for the distortion.This has been applied to the first-order Jahn-Teller distortion in Pughe et al. (2023).
These parameterisations each have merits.However, they are not sensitive to the symmetry of the octahedral distortion.The van Vleck modes are conceptually different to each of these for quantifying the Jahn-Teller distortion because they can be used to IUCr macros version 2.1.10:2016/01/28 quantify distortion with the precise symmetry of the transition metal e g orbitals.This is important because Jahn-Teller distortions typically follow a particular symmetry.
When the distortion is due to degeneracy in the e g orbitals it will be of E g symmetry; when it is due to degeneracy in the t 2g -degenerate orbitals it may be either E g or T 2g symmetry (Child & Roach, 1965;Bacci et al., 1975;Holland et al., 2002;Halcrow, 2009;Teyssier et al., 2016;Schmitt et al., 2020;Streltsov et al., 2022), although there is relatively little unambiguous experimental evidence for a Jahn-Teller-induced shear as compared with more typical E g distortion..
In this paper, we present a Python (Van Rossum & Drake, 2009) package, Van-VleckCalculator, for calculating the van Vleck distortion modes.We show that the approach to calculating the modes which is commonly used in the literature is a reasonable approximation for octahedra with negligible angular distortion, but results in the loss of information in other cases.We propose a new metric, the shear fraction η, for understanding the correlation between octahedral shear and angular distortion.
Finally, we re-analyse some previously-published data in terms of the van Vleck modes to show that these can be an effective way of understanding octahedral behaviour.

Bond angle distortions
Fig. 2. The 6 van Vleck modes exhibited for an octahedron, with sites labelled using the notation in the Theory section.For the octahedra exhibiting Q 1 , Q 2 , and Q 3 distortions, there is no angular distortion; for the octahedra exhibiting Q 4 , Q 5 , and Q 6 distortions, there is no bond length distortion.For the octahedral shear (Q 4 , Q 5 , and Q 6 ) modes, axes are drawn to show where the bond directions would be if undistorted.An octahedron can exhibit several, or all, of these distortions simultaneously.

Theory
Within an octahedron, we can split the 6 ligand ions into three pairs, where the two ions within the pair are opposite one another.In the absence of angular distortion (i.e., assuming all ligand-core-ligand angles are an integer number of 90 a basis where each of the three axes exist directly along the x-, y-, and z-axis, and where the origin in space is defined as the centre of the octahedron.
Each pair within an octahedron can therefore be assigned to an axis and labelled as the a, b, or c pair respectively.Within a pair, ions can be labelled as − or + depending on whether they occur at a negative or positive displacement from the origin, along the axis, respectively.This notation is demonstrated in Figure 2, where each pair of ions is represented by a different colour.
For each of the 6 ligands, we define a set of coordinates: x α β , y α β , and z α β , where α is a, b, or c denoting the pair in which the ligand is, and β is − or + denoting which ion within the pair.
Using these, we further define a set of van Vleck coordinates (capitalised to distinguish from true coordinates) which is the displacement of the ion within an axis away from its ideal position.For instance, for the ion with α = a and β = −: Using these coordinates, the first six van Vleck modes (Q j ; j = 1 − 6) are defined as follows (Van Vleck, 1939): IUCr macros version 2.1.10:2016/01/28 We only discuss these first six van Vleck modes, which are shown in Figure 2.
IUCr macros version 2.1.10:2016/01/28 shear distortions.Q 1 is a simple expansion/contraction mode which does not affect symmetry and will not be discussed further.
Q 2 and Q 3 are a planar rhombic distortion and a tetragonal distortion respectively; they are considered degenerate due to the Hamiltonian, which is discussed for instance in Kanamori (1960).These two modes form a basis for distortions describing different octahedral configurations with the symmetry of the transition metal e g orbitals (Goodenough, 1998;Khomskii & Streltsov, 2020).These modes are of most relevance for first-order Jahn-Teller distortions occurring due to degenerate e g orbitals.
A phase space of possible octahedral configurations can be constructed using these two parameters (Kanamori, 1960), as shown in Figure 3.Here the magnitude of the distortion ρ 0 can be calculated as follows: and the angle3 ϕ of this distortion from being of purely Q 3 character can be calculated by: All possible combinations of the Q 2 and Q 3 modes correspond to a particular angle ϕ, and hence a particular configuration as shown in Figure 3.The structural effect of a rotation of ϕ within a range of 120 • can be quite significant, as shown in Figure 3;  octahedral complexes.Note that for angles which are not special angles, there will be mixing of the orbital states of the nearest two special angles.
A characteristic of the Jahn-Teller distortion is that, in the absence of external distortive forces, the symmetry of the structure matches the symmetry of the orbitals involved.Typically, any d-orbital Jahn-Teller distortion will have some planar rhombic (Q 2 ) or tetragonal (Q 3 ) character.However, sometimes when the degeneracy occurs in the t 2g orbital, there may instead be a trigonal component to the symmetry of the distortion, which manifests as an angular distortion instead (Child & Roach, 1965;Bacci et al., 1975;Holland et al., 2002;Halcrow, 2009;Teyssier et al., 2016;Schmitt et al., 2020;Streltsov et al., 2022).For the more commonly-studied case of a degeneracy in the e g orbitals, the effect of a rotation of ϕ similarly changes the symmetry of the d orbitals.Figure 1 shows the splitting of the d orbitals in an octahedrally-coordinated d 7 transition metal due to an elongation-type first-order Jahn-Teller distortion, where the tetragonal elongation occurs along the z-axis.Note that the unpaired e g electron occupies the d z 2 orbital.In the opposite case of a compression-type first-order Jahn-Teller distortion along the z axis, the lower-energy, and hence singly-occupied, orbital would be the d x 2 −y 2 ; this would correspond to a rotation in ϕ of 180 • .More generally, as a function of ϕ, there exist a set of special angles separated by a 60 • rotation corresponding to a particular e g orbital being singly-occupied by a d electron.These are tabulated in Table 1.An octahedron for which ϕ does not correspond to one of these special angles exhibits orbital mixing (Rodriguez-Carvajal et al., 1998;Zhou & Goodenough, 2008b).
IUCr macros version 2.1.10:2016/01/28 The Q 4 to Q 6 modes describe shear of the octahedra, i.e. the effect whereby paired ligands at opposite sides of a central ions are displaced in opposite directions, and have trigonal T 2g character.The shear modes may be used to quantify the Jahn-Teller distortion in octahedra where the degeneracy occurs within t 2g orbitals (Child & Roach, 1965;Teyssier et al., 2016).The magnitude of the calculated shear is typically correlated with angular distortion, which is commonly quantified using the σ 2 ζ metric called the Bond Angle Variance (Robinson et al., 1971) (BAV), defined here as: where m is the number of bond angles (i.e. 12 for octahedra), ζ i is the ith bond angle, and ζ 0 is the ideal bond angle for a regular polyhedron (i.e. 90 • for an octahedron).
However, for direct comparison to the shear modes, it is more appropriate to use the standard deviation σ ζ .η=0 η=1 η=0.5  22.In the case for η = 0, the only distortion is anti-shear within a single plane.In the case for η = 0.5, there are two planes in which there is distortion, a shear and anti-shear distortion equal in magnitude.In the case for η = 1, there is a plane with a purely shear distortion.
For an octahedron with non-zero T 2g (Q 4 , Q 5 , Q 6 ) modes, increasing their magnitude will increase the angular distortion, but an octahedron may have angular distortion without exhibiting octahedral shear.To analyse the extent to which angular distortion in an octahedron is due to shear, we propose a shear fraction parameter η, demonstrated in Figure 4 and defined below.
First, we must define a set of shear and "anti-shear" angular indices, which are modifications of Equations 8 to 10 in terms of angles rather than displacements.The indices are represented with ∆ and a subscript corresponding to the plane in which rotation occurs: the ab-plane corresponds to the Q 4 mode, the ac-plane to the Q 5 mode, and the bc-plane to the Q 6 mode.The absence or presence of a prime symbol, ′, designates whether the index represents shear or anti-shear respectively.Finally, the IUCr macros version 2.1.10:2016/01/28 δ angle is the rotation of the ligand from its ideal van Vleck coordinate in a clockwise direction, within the plane in which the corresponding van Vleck shear (Q 4 to Q 6 ) would occur.These are defined thus (see SI, Figure S7): We then quantify the shear and "anti-shear" distortions using the following equations: From here, we define the shear fraction η as follows: IUCr macros version 2.1.10:2016/01/28 This η parameter will be important in interpreting the relation between the angular distortion, σ ζ , and the van Vleck shear modes Q 4 to Q 6 .

Implementation
In this section, the algorithm used to calculate Van Vleck distortion modes is discussed.
It is written using Python 3 (Van Rossum & Drake, 2009) as a package called VanVleckCalculator, with the full code available on GitHub (Nagle-Cocco, 2023), and also presented with annotations in the Supplementary Information.Data handling and some calculations make use of NumPy (Harris et al., 2020), and crystal structures are handled using PyMatGen (Ong et al., 2013).
A flow chart showing the octahedral rotation algorithm can be found in Supplementary Information, Figure S1.
Besides calculating the van Vleck modes and the angular shear modes described in this paper, VanVleckCalculator can also calculate various other parameters as described in Supplementary Information.

Selecting an origin
Selection of the origin is a key step in calculating van Vleck modes.The most common approach, for an M X 6 octahedron, is to take the M ion as the origin.This is a reasonable approach, given that M ions are typically positioned at, or very close to, the centre of an octahedron.This is particularly appropriate for unit cells derived from Rietveld refinement (Rietveld, 1969) of Bragg diffraction data, where the M ion is likely to occur at a high-symmetry Wyckoff site.A third, similar, option would be to choose the average position of the 6 ligands as the origin in space.An example of when this may be a desirable choice would be for systems exhibiting a pseudo Jahn-Teller effect (also called the second-order Jahn-Teller effect), where the central cation IUCr macros version 2.1.10:2016/01/28 is offset from the centre of the octahedron.
In some instances, a crystal structure may be simulated using a supercell.Examples include so-called "big box" Pair Distribution Function (PDF) analysis (Tucker et al., 2007) and Molecular Dynamics (MD) (Bocharov et al., 2020) simulations.Such a supercell typically retains the periodicity which is an axiom of a typical crystallographic unit cell, but will exhibit local variations.For instance, a unit cell obtained by analysis of Bragg diffraction data is typically regarded as an "average" structure, insensitive to local phenomena such as thermally-driven atomic motion or disordered atomic displacements such as a non-cooperative Jahn-Teller distortion.In a crystallographic unit cell, thermal motion of atoms is typically represented by variable Atomic Displacement Parameters (ADPs) (Peterse & Palm, 1966).In contrast, a supercell should reflect local phenomena, for instance exhibiting local Jahn-Teller distortions in a system with a non-cooperative Jahn-Teller distortion, and representing thermal effects not with ADPs but rather by distributing equivalent atoms in adjacent repeating units in slightly different positions.In this regard, a supercell can be considered a "snapshot" of a crystal system at a point in time.It may not be appropriate to set the core ion as the centre of the octahedron in a supercell, therefore, as the positioning of both core and ligand ions is in part due to thermal effects, and so the "centre" of the octahedron will be displaced due to random motion.The alternative option would be to simply use the crystallographic site of the central ion and fix this as independent of the precise motion of the central ion locally.
In VanVleckCalculator, the user has the option to take as the centre of the octahedron either the central ion, the average position of the 6 ligands, or a specified set of coordinates.

Calculating van Vleck modes along bond directions
The calculation of the van Vleck modes, as described in the Theory section, requires that the basis in space be the octahedral axes (i.e. the three orthogonal axes entering the octahedron via one vertex, passing through the central ion, and exiting via the opposite vertex).For a given crystal structure, this may require that an octahedron be rotated about each of the three axes making up the basis, until the octahedral axes perfectly align with the basis.This becomes more complicated when the octahedron exhibits angular distortion (i.e.exhibits ligand-core-ligand angles are not integer multiples of 90 • ).In this case, it is impossible to define octahedral axes according to the strict criteria previously defined.

Calculating van Vleck modes within Cartesian coordinates
In VanVleckCalculator, we have written an algorithm for rotating an octahedron about three Cartesian axes with a defined origin within the octahedron, such that the ligands are as close as possible to the axes (within the constraint that there is angular distortion).This allows for calculation of van Vleck modes in a way that does not artificially constrain the octahedral shear modes (Q 4 , Q 5 , and Q 6 ) to be zero.
First, three orthogonal axes are taken as the x-, y-, and z-axes6 .By default, these are the [1,0,0], [0,1,0], and [0,0,1] axes respectively, but alternative sets of orthogonal vectors can be given by the user; for instance, for regular octahedra rotated 45 • about the x axis, the user would be recommended to give as axes and . This vector is given as a Python list with shape (3,3).For consistency, the cross product of the first two axes should always be parallel with the third given vector; if anti-parallel, the algorithm will automatically multiply all elements in the third vector by -1.The three pairs of the octahedron (as defined in the Theory section) are each then assigned to one of these three axes on the basis of which pair has the largest projection of its displacement (the vector between two on a particular axis, with the z-axis assigned first, then the y-axis from amongst the two pairs not assigned to the z-axis, then the x axis is automatically assigned to the remaining pair).Within each pair, the ligands are then ordered such that the ligand with the negative distance is along the assigned vector first, then the ligand with positive distance occurs second.
Second, the octahedron is rotated about the x-, y-, and z-directions of the basis repeatedly until the orthogonal axes supplied in the previous step match the basis precisely.This is performed in a while loop structure, with the rotation angles about the three axes summed in quadrature and compared with a defined tolerance (by default, 3×10 −4 radians in VanVleckCalculator), and if the total rotation exceeds the tolerance, the step is repeated7 .This step is usually unnecessary, and can be skipped by leaving the default set of orthogonal axes, which are [1,0,0], [0,1,0], and [0,0,1] (meaning no rotation will occur).
Third, an automatic rotation algorithm will further minimise the effect of angular distortion.For each of the three axes, the four ligands not intended to align with that axis are selected.The angle to rotate these four ligands about the origin such that each is aligned with its intended axis within the plane perpendicular to the axis of rotation is calculated.The octahedron is then rotated about this axis by the average of these four angles.This occurs iteratively until, for a given iteration, the sum (in quadrature) of the three rotation angles is less than the already-mentioned defined tolerance.
At this point, the octahedron is optimally aligned with the basis (given the limitation that there may be angular distortion) and the van Vleck modes can be calculated.

Ignoring or including angular distortion: a comparison
To evaluate the utility of calculating the van Vleck modes without disregarding the angular distortion, we perform a comparison between the two approaches.We have calculated the van Vleck distortion modes and associated parameters for octahedra in NaNiO 2 and LaMnO 3 with both a method that ignores angular distortion and calculates modes along bond directions (consistent with the Q 2 and Q 3 equations defined by Kanamori (1960)), and a method that used Cartesian coordinates in order to take angular distortion into account.Table 2 shows this for these two materials.
Firstly, for the van Vleck modes calculated without ignoring angular distortion, we can see the octahedral shear modes (Q 4 , Q 5 , Q 6 ) are larger for the material with higher angular distortion (as quantified using bond angle variance).While the effect of ignoring angular distortion is significant for the Q 4 , Q 5 , and Q 6 modes, it makes negligible difference for the calculation of Q 2 and Q 3 modes, and the associated ρ 0 and ϕ parameters.It is therefore likely a reasonable approximation to take, particularly for calculation of ϕ as is common in literature, even for octahedra which exhibit higher angular distortion.However, there is a definite loss of information in assuming the shear modes Q 4 to Q 6 are zero.The impact of this is assessed in the case studies.

and η values for
NaNiO 2 and LaMnO 3 at room-temperature, calculated using orthogonal axes as described in this report ("Cartesian" method), and alternatively by ignoring angular distortion and calculating van Vleck modes along bond lengths ("Kanamori" method † ).The centre of the octahedron is taken as the central ion position.Values were calculated using crystal structures reported on the Inorganic Crystal Structure Database (ICSD).To demonstrate the difference in angular distortion, the Bond Angle Variance (defined in Equation 13) is also tabulated.BAV is rounded to the third significant figure; Q modes and related parameters are rounded to the 4th decimal place.

Temperature-dependence of octahedral shear in LaAlO 3
Perovskite and perovskite-like crystal structures are amongst the most important and widely-studied crystalline material classes in materials science today.Perovskite crystal structures have ABX 3 chemical formulae, with A and B being ions at the centres of dodecagons and octahedra, respectively, with the X anion constituting the vertices of these polyhedra.The BX 6 octahedra interact via corner-sharing interactions.There are also perovskite-like crystal structures such as the double perovskites, A 2 BB ′ X 6 (King & Woodward, 2010;Koskelo et al., 2023), for which many of the IUCr macros version 2.1.10:2016/01/28 same principles apply.The ideal perovskite system would be cubic, with space group P m 3m, but many related structures with lower symmetry are known.This typically occurs in three situations (Woodward, 1997): 1. when there is a mismatch between the ionic radii of the octahedrally-coordinated BX 6 cation and the dodecagonally-coordinated AX 12 cation, resulting in tilting of the octahedra; see Figure 5(e).
2. when there is displacement of the central cation from the centre of the octahedron, typically due to the pseudo Jahn-Teller effect.
3. when the ligands of the octahedron are distorted by electronic phenomena such as the first-order Jahn-Teller effect.
In this case study, we focus on the first case, where a size mismatch results in octahedral tilting.Octahedra are often modelled as rigid bodies, but in practice they are not rigid in all systems, and the octahedral tilting will often induce strain resulting in angular distortion.This is typically far smaller than that seen in edge-sharing materials such as NaNiO 2 , but it is large enough that it cannot be disregarded when attempting to fully understand the structure of the material.As was noted by Darlington (Darlington, 1996), this angular distortion commonly manifests as shear.LaAlO 3 is a perovskite-like ABX 3 material which is cubic (space group P m 3m) above around ∼830 K, but which exhibits a rhombohedral distortion below this temperature (with space group R 3c) due to octahedral tilting (Hayward et al., 2005), see Here, we calculate the van Vleck shear modes.Due to the symmetry of the octahedral tilting, there is only one independent shear mode, and We compare this with the bond angle standard deviation given in Equation 13, see Figure 5.We see that despite being distinct parameters, the temperature dependence of both is entirely identical.We attribute this to the shear fraction, η, being precisely 1 for all temperatures where there is angular distortion, meaning that shear is completely correlated with angular distortion.

Big box analysis of Pair Distribution Function data on LaMnO 3
The Jahn-Teller distortion in LaMnO 3 , a perovskite-like ABX 3 material which has the crystal structure shown in Figure 6(a), occurs as a consequence of degeneracy in IUCr macros version 2.1.10:2016/01/28 the e g orbitals on the high-spin d 4 Mn 3+ ion.At ambient temperatures, it is a prime example of a cooperative Jahn-Teller distortion, exhibiting long-range orbital order where the elongation of the Jahn-Teller axis alternates between the a and b directions for neighbouring MnO 6 octahedra, never occurring along the c direction (Khomskii & Streltsov, 2020) [Figure 6(b)].With heating through ∼750 K, the Jahn-Teller distortion can no longer be observed in the average structure obtained from Bragg diffraction (Rodriguez-Carvajal et al., 1998).However, the Jahn-Teller distortion persists locally as has been shown by pair distribution function (Qiu et al., 2005) and EXAFS (García et al., 2005;Souza et al., 2005) measurements.This transition is one of the most widely-studied orbital order-disorder transitions for the first-order Jahn-Teller distortion.The high-temperature orbital regime has been described theoretically in terms of a three-state Potts model (Ahmed & Gehring, 2006;Ahmed & Gehring, 2009), a view supported by big box analysis of combined neutron and x-ray pair distribution function data (Thygesen et al., 2017), as performed using RMCProfile (Tucker et al., 2007).
Orthogonal axes  (2017).In (c), orthogonal axes were used (i.e.angular distortion was included in the calculation, using the method described in this manuscript), whereas in (d) the Mn-O bond directions were taken as the axes regardless of orthogonality.

Bond directions
(e) the Mn 3+ octahedra which exhibit a mixed Q 2 -Q 3 type distortion due to the first-order Jahn-Teller effect, manifesting as three different bond lengths, labelled in ascending order of length as s (orange), m (grey), and l (green).(f) a histogram of the smallest to largest Mn-O bond length within each octahedron in the 10 × 10 × 8 supercell, with the blue vertical lines indicating the bond lengths in the average structure.
In this case study, we take a 10 × 10 × 8 supercell of LaMnO 3 , obtained using RMCProfile against total scattering data obtained at room-temperature, and previously published in the aforementioned work (Thygesen et al., 2017).Results are shown in Figure 6.We repeat the analysis of this supercell from the perspective of IUCr macros version 2.1.10:2016/01/28 the E g (Q 2 , Q 3 ) van Vleck distortion modes, using two different approaches: (1) the algorithm for automatically determining a set of orthogonal axes is applied to each octahedron individually, and (2) following the van Vleck equations 23 and 24 proposed by Kanamori (1960) where angular distortion is disregarded.In each of these cases the crystallographic site of the supercell is taken as the origin, and so thermally-driven variations in the Mn position will not affect the result.
As can be seen in Figures 6(d)-(f), there are two clusters of octahedra within the polar plot, occurring at ϕ ≈ ±107 • .This corresponds to occupation of the d y 2 orbitals (+) and of the d x 2 orbitals (−).In both cases, the superposition of perpendicular Q 3 compression and elongation modes results in an octahedron with mixed Q 2 -Q 3 character.This finding is consistent with previous works which placed MnO 6 octahedra from LaMnO 3 into the framework of an E g (Q 2 , Q 3 ) polar plot (Zhou & Goodenough, 2008a;Zhou et al., 2011).The Q 2 contribution to the distortion, as seen from the three different Mn-O bond lengths in LaMnO 3 , is also present in Jahn-Teller-distorted ACuF 3 (A=Na,K,Rb) (Lufaso & Woodward, 2004;Marshall et al., 2013;Khomskii & Streltsov, 2020) and even in some Jahn-Teller-undistorted perovskites (Zhou & Goodenough, 2008a), indicating it is related to the structure.It is not intrinsic to Jahn-Teller-distorted manganates, as it is absent in high-spin d 4 Mn 3+ with edge-sharing octahedral interactions and colinear orbital ordering such as α-NaMnO 2 and LiMnO 2 (checked using ICSD references 15769 and 82993 (Jansen & Hoppe, 1973;Armstrong & Bruce, 1996) respectively).
The Q 2 component to the octahedral distortion is therefore likely intrinsic to the crys-IUCr macros version 2.1.10:2016/01/28 tal structure (Zhou & Goodenough, 2006;Zhou & Goodenough, 2008a), which occurs as a result of octahedral tilting reducing the symmetry from cubic P m 3m to P nma.
In LaMnO 3 , the combination of the Q 2 component to the distortion and the orbital ordering [Figure 6(b)] are a possible distortion of the P nma space group.In this way, the orbital ordering may be coupled to the octahedral tilting, a link previously made by Lufaso & Woodward (2004).
Finally, we also calculate the Q 4 to Q 6 octahedral shear modes for all octahedra in the supercell, presented as a histogram in Figure S2 in Supplementary Information.We present the average and standard deviation, as calculated assuming orthogonal axes and with the automated octahedral rotation: Q 4 = −0.02±0.13Å, Q 5 = 0.02±0.10Å, and Q 6 = −0.00± 0.11 Å.In each case, the magnitude of the distortion is zero within standard deviation, and also contains the value from the average structure presented in Table 2 within the range of error.This low level of shear generally supports the validity of calculating the E g (Q 2 , Q 3 ) van Vleck modes along bond directions rather than a Cartesian coordinate system for a system like LaMnO 3 .It is interesting to note that the standard deviation is higher for Q 4 , which quantifies the shear within the plane in which there is orbital ordering.

Effect of pressure on the JT distortion in NaNiO 2
In recent years, there have been several studies looking at the effect of applied pressure on the Jahn-Teller distortion in crystalline materials ( Åsbrink et al., 1999;Loa et al., 2001;Choi et al., 2006;Zhou et al., 2008;Zhou et al., 2011;Aguado et al., 2012;Mota et al., 2014;Caslin et al., 2016;Zhao et al., 2016;Collings et al., 2018;Bhadram et al., 2021;Lawler et al., 2021;Scatena et al., 2021;Ovsyannikov et al., 2021;Nagle-Cocco et al., 2022).Most of these have shown that, as a general rule, pressure reduces the magnitude of the Jahn-Teller distortion as a consequence of the elongated bond IUCr macros version 2.1.10:2016/01/28 being more compressible than the shorter bonds.Zhou et al. (2011) use van Vleck modes to quantify the effect of pressure on the Jahn-Teller distortion in the cornersharing perovskite-like compounds LaMnO 3 and KCuF 3 .While application of pressure reduces the magnitude of the distortion, as quantified using ρ 0 (Equation 11), they argue that it does not change the orbital mixing ϕ (Equation 12).KCuF  22.Note that for the NiO 6 octahedra, Q 5 = Q 6 , whereas for FeO 6 octahedra, For Fe 2 O 3 , the average position of the O ligands were taken as the centre of the octahedron.
We previously studied the effect of pressure on the Jahn-Teller distortion in NaNiO 2 (Nagle- et al., 2022), by performing Rietveld refinement (Rietveld, 1969) of neutron diffraction data from the PEARL instrument (Bull et al., 2016) at the ISIS Neutron and Muon Source.However, we did not utilise the van Vleck distortion modes, instead quantifying the Jahn-Teller distortion using the bond length distortion index (Baur, 1974) and the effective coordination number (Hoppe, 1979).In that study, we found no deviation from the ambient-pressure space group C2/m (Dick et al., 1997;Sofin & Jansen, 2005), shown in Figure 7(a), for all pressure points at room-temperature up to ∼4.5 GPa.This space group permits only four short{long} and two long{short} bonds or 6 equal bond lengths, depending on the angle β, and so throughout the measured IUCr macros version 2.1.10:2016/01/28

Cocco
pressure range there exists no Q 2 character to the Jahn-Teller distortion, consistent with the principle that hydrostatic pressure does not change orbital mixing (Zhou et al., 2011).
Here, we perform a fresh analysis of the variable-pressure octahedral behaviour as a function of pressure at room temperature in NaNiO 2 in terms of the E g (Q 2 , Q 3 ) van Vleck distortion modes.For a reference we sought a material which does not exhibit a first-order Jahn-Teller distortion but does exhibit bond length distortion; for this purpose, we selected Fe 2 O 3 , the pressure dependence of which was previously studied in Finger & Hazen (1980), and which exhibits bond length distortion due to its face-and edge-sharing octahedral connectivity.Fe 2 O 3 contains high-spin d 5 Fe 3+ cations within octahedra which interact via both face-and edge-sharing interactions.
It should be noted that Fe 2 O 3 likely exhibits some very subtle pseudo Jahn-Teller distortion (related to, but distinct from the first Jahn-Teller effect discussed here) on account of the Fe 3+ ions (Cumby & Attfield, 2017;Bersuker & Polinger, 2020), but this does not impact the discussion in any meaningful way.
In Figure 7(c) we compare (for NaNiO 2 ) ρ 0 with three other parameters (bond length distortion index, quadratic elongation, and effective coordination number) which are often used to parametrise the magnitude of the Jahn-Teller distortion.The trend for each is near identical, although the magnitudes differ greatly, indicating that each is a reasonable parameter for quantifying the magnitude of the Jahn-Teller distortion.This can be compared to Figure 7(e) which shows the same parameters for the Jahn-Teller-undistorted FeO 6 octahedra in Fe 2 O 3 , where it can be seen that ρ 0 remains approximately at zero throughout the measured pressure range, despite a high level of bond length distortion as represented by the bond length distortion index, effective coordination number, and quadratic elongation (a similar plot for KCuF 3 can be seen in SI, Figure S3).This means that, while these parameters are valid for quan-tifying the magnitude of Jahn-Teller distortion, they are also sensitive to other kinds of distortion.ρ 0 is calculated using Q 2 and Q 3 which have E g symmetry, and so ρ 0 will only be non-zero for a distortion with E g symmetry.Thus, it is arguably the ideal choice for parameterising the magnitude of this type of Jahn-Teller distortion.However, while ρ 0 is more reliable than the other parameters shown in Figures 7(c,d) for demonstrating the presence of a Jahn-Teller distortion, it is not always strictly zero for a Jahn-Teller-inactive octahedron, as it will have a non-zero value if the octahedron is distorted with an e g symmetry.For example, the NaO 6 octahedron in C2/m NaNiO 2 has the same symmetry as the NiO 6 octahedron, and so exhibits a value of ρ 0 between 0.065 and 0.05 within the studied pressure range [Figure S4 in Supplementary Information], and Jahn-Teller-inactive FeO 6 octahedra in RFeO 3 perovskites have non-zero ρ 0 due to the E g symmetry of the distorted octahedra, as shown in Zhou & Goodenough (2008a).
Figure 8 shows a polar plot for the behaviour of NaNiO 2 and KCuF 3 in the range 0 to 5 GPa (the measured range for NaNiO 2 ).It can be seen that within this pressure range, the magnitude of the Jahn-Teller distortion decreases far more for KCuF 3 than NaNiO 2 ; this reflects the fact that KCuF 3 is more compressible, with a bulk modulus 57(1) GPa (Zhou et al., 2011) compared with 121(2) GPa for NaNiO 2 (Nagle-Cocco et al., 2022), as obtained by a fit to the third-order Birch-Murnaghan equation-ofstate (Birch, 1947).Within this pressure range we see that ϕ does not change with pressure for either material, and that this property is true regardless of whether ϕ is or is not a special angle (as in Table 1), consistent with the interpretation of Zhou

Conclusion
We present VanVleckCalculator, a code package written in Python 3 for the calculation of octahedral van Vleck distortion modes.These modes are particularly important for understanding the behaviour of the Jahn-Teller distortion, and we have shown that the parameter ρ 0 (which is based on the van Vleck Q 2 and Q 3 modes) is a more reliable way of quantifying the Jahn-Teller distortion than other oft-used parameters such as the bond length distortion index.
We show the importance of using a Cartesian set of coordinates for this calculation, instead of calculating the modes along bond directions, as is often done in the litera-IUCr macros version 2.1.10:2016/01/28 ture.This is because calculating the van Vleck distortion modes along bond directions relies on the assumption that there is no angular distortion or octahedral shear, which is often a false assumption and artificially constrains the Q 4 , Q 5 , and Q 6 modes to be zero.We show that there is value in calculating these later modes, for instance in understanding the effect of octahedral tiling on octahedra in perovskite-like materials.
These shear modes will also be useful for parameterising the Jahn-Teller effect when the degeneracy occurs in the t 2g orbitals and results in a trigonal distortion, because their symmetry matches the distortion.
We also show that octahedral shear correlates with angular distortion for materials under the influence of tuning parameters such as pressure or temperature where there is a continuously-varying distortion, such as octahedral tilting (as in LaAlO 3 ) or firstorder Jahn-Teller distortion (as in NaNiO 2 ).However, there is no correlation when the distortion is due to competing interactions due to face-or edge-sharing octahedra (as in Fe 2 O 3 ).We propose a new parameter, the shear fraction η (defined in Equation 22), which can be used to predict whether there will be correlation between octahedral shear modes and angular distortion.

Fig. 1 .
Fig. 1.The orbital rearrangement due to a tetragonal elongation for an octahedrallycoordinated low-spin d 7 transition metal ion, which typically occurs due to the first-order Jahn-Teller effect.
and Z a − = z a − .See Figure2for clarification of the ion notation.

Fig. 3 .
Fig. 3.The Q 2 -Q 3 phase space for elongated octahedra, with a representation of the values ρ 0 and ϕ.Based on a figure from an article by Goodwin (2017).

Fig. 4 .
Fig. 4. Three possible octahedral shear/anti-shear distortions, with the associated value of the shear fraction η as defined in Equation22.In the case for η = 0, the only distortion is anti-shear within a single plane.In the case for η = 0.5, there are two planes in which there is distortion, a shear and anti-shear distortion equal in magnitude.In the case for η = 1, there is a plane with a purely shear distortion.
IUCr macros version 2.1.10:2016/01/28 see Figure5(e) and (f).Throughout both temperature regimes, there is the absence of bond length distortion; a calculation of the bond length distortion index would yield a value of zero at all temperatures.In the low-temperature regime, the magnitude of the distortion continuously decreases with increasing temperature, reaching zero at the transition temperature.Most commonly in the literature, the tilting angle between the octahedral axis and the c-axis (0 • in the cubic phase) is used to quantify this distortion; for LaAlO 3 , this is shown in Figure5(a).The strain induced by this distortion results in intra-octahedral angular distortion.Hayward et al. (2005) model this in terms of strain tensors, finding a linear temperature dependence below the transition temperature, which differs from the temperature-dependence of the tilting angle (which resembles an exponential decline), implying the two are related but distinct phenomena.Cumby & Attfield (2017) instead model the octahedral distortion for this same dataset using the radii of a minimum-bounding ellipsoid, and also find approximately linear temperature dependence of the long and short radii as they approach convergence (see Figure5(b)).

Fig. 6 .
Fig. 6.The perovskite-like structure of LaMnO 3 , as obtained from ICSD structure 50334, is shown in (a).(b) the orbital ordering at room-temperature in LaMnO 3 , and is reprinted with permission from Khomskii & Streltsov, Chem.Rev. (2021), 121, 5, 2992-3030, copyright 2021 American Chemical Society.(c) and (d) show polar plots with a point representing the calculated ϕ and ρ 0 values for each MnO 6 octahedron in a 10 × 10 × 8 supercell of LaMnO 3 at room-temperature, as obtained from reverse Monte Carlo analysis of neutron Pair Distribution Function data in Thygesen et al.(2017).In (c), orthogonal axes were used (i.e.angular distortion was included in the calculation, using the method described in this manuscript), whereas in (d) the Mn-O bond directions were taken as the axes regardless of orthogonality.(e) the Mn 3+ octahedra which exhibit a mixed Q 2 -Q 3 type distortion due to the first-order Jahn-Teller effect, manifesting as three different bond lengths, labelled in ascending order of length as s (orange), m (grey), and l (green).(f) a histogram of the smallest to largest Mn-O bond length within each octahedron in the 10 × 10 × 8 supercell, with the blue vertical lines indicating the bond lengths in the average structure.

Figure 6
Figure 6(e) shows the MnO 6 octahedron in the average structure of LaMnO 3 at room temperature, with the three different bond lengths plotted in Figure 6(f) along with a histogram of all the bond lengths in the supercell.This shows how the combination of the Q 2 and Q 3 distortion modes manifests in the octahedral distortion.

Fig. 7 .Fig. 8 .Fig. 9 .
Fig. 7.The crystal structures of Jahn-Teller-active C2/m NaNiO 2 and inactive R 3c Fe 2 O 3 are shown in (a) and (b) respectively.(c) and (d) show a comparison of various metrics for quantifying the degree of Jahn-Teller distortion as a function of pressure, for NiO 6 octahedra in NaNiO 2 and FeO 6 octahedra in Fe 2 O 3 respectively.The parameters subject to comparison are the magnitude ρ 0 , bond length distortion index, effective coordination number, and quadratic elongation.Dashed lines indicate a linear fit to the data, whereas solid lines connect data points.
the previous study(Nagle-Cocco et al., 2022), we showed using specific O-Ni-O bond angles that pressure reduces the angular distortion for NaNiO 2 .Here, we show that pressure also reduces the related shear distortion in NaNiO 2 .This is demonstrated in Figure9where we plot the octahedral shear Q 4 , Q 5 , and Q 6 modes for NaNiO 2 and Fe 2 O 3 against the bond angle standard deviation, σ ζ , defined in Equation13.Unlike the AlO 6 octahedra in LaAlO 3 [Figure5], for NiO 6 octahedra in NaNiO 2 there is no perfect correlation between the shear modes and angular distortion despite η ≈ 1, because there is more than one independent shear mode, but we can see that shear distortion and angular distortion are still highly correlated.However, for Fe 2 O 3 the shear fraction η << 1 and there is no correlation between the shear distortion modes and angular distortion.This difference in behaviour likely arises because the main driver of the change is a continuous decrease in the Jahn-Teller distortion in NaNiO 2 , as compared to Fe 2 O 3 where positions of the oxygen anions are determined by the reduced degrees of freedom arising from trying to satisfy multiple face-and edge-sharing interactions.This result could only be achieved by calculating the van Vleck modes in a Cartesian coordinate system as outlined in this paper, as opposed to calculating the distortion modes along bond directions, indicating the relevance of calculating the van Vleck modes in this way, and of the shear fraction η we propose in this work.
scholarship EP/R513180/1 to pursue doctoral research from the UK Engineering and Physical Sciences Research Council (EPSRC).

Table 1 .
The special angles in the Q 2 -Q 3 phase space [Figure3], as a function of ϕ = arctan (Q 2 /Q 3 ), with the associated singly-occupied e g orbital, for d 4 and low-spin d 7