Jsymm: A Python package for symmetry analysis of exchange tensors in magnetic Hamiltonians

To find the constraints imposed by symmetry on the exchange tensors of a bond, we first need to understand how these tensors are transformed under the action of symmetry elements. Recall that spins transform as axial vectors:

where $g$ is a symmetry element, $R_{g}$ is its matrix, and $\mathbf{S}=(S_{x},S_{y},S_{z})^{T}$. We define the corresponding action on the Hamiltonian term by

An exchange term in the Hamiltonian can be written as

where $J$ is the exchange matrix. In particular, the DMI matrix has the form

where $D_{\alpha}$ are the components of the vector $\mathbf{D}$. Since the matrix $D$ is antisymmetric, we have

which one can interpret as

In other words, the vector $\mathbf{D}$ must change sign upon the permutation of the spin labels to compensate the sign change of the cross product and leave the energy contribution invariant.

Suppose that the group element $g$ maps each atom of the bond to itself. Then it follows from Eqs. (4), (5) and (6) that the matrix $J$ transforms as

where we used that $R_{g}$ is orthogonal, so $R_{g^{-1}}=R_{g}^{-1}=R_{g}^{T}$. In the case when $g$ interchanges two atoms, the DMI matrix also changes its sign, according to Eq. (9). We conclude that the symmetry element that maps the bond to itself acts on the exchange matrices as follows:

where the minus sign corresponds to the symmetries that flip the bond.

One can also derive the transformation rule for the vector $\mathbf{D}$. Since the cross product transforms as an axial vector,

Thus, $\mathbf{D}$ also transforms as an axial vector, which additionally changes its sign when $g$ flips the bond.

Our next step is to determine which symmetry elements contribute to the constraints of the first kind introduced above and which bonds are affected by the constraints of the second kind. To this end, we define the notions of a stabilizer subgroup and an orbit of a bond, taking into account the lattice periodicity.

Recall that the point group $G$ of a crystal is defined as a quotient $G=S/T$ of the space group $S$ by the subgroup $T$ of lattice translations. The point group consists of cosets

where $R$ is the reflection or rotation matrix and $\mathbf{t}$ is the translation vector. For each coset $[g]$, we choose a representative element $g=\{R|\mathbf{t}\}$, where $\mathbf{t}$ has non-negative fractional components.

Our goal is to define the action of the point group $G$ on the bonds formed by pairs of atoms. As an intermediate step, consider the action of the representative $g$ on an atom $a$ with the position vector $\mathbf{r}$:

We define a bond as an ordered pair of atoms and denote it

There is no need to distinguish bonds that are related by a uniform lattice translation of both atoms. Formally, one can describe this by an equivalence relation $\sim$ on the set of all bonds:

where $\mathbf{T}$ is a translation by a lattice vector. Now we can consider the equivalence classes of bonds $[\mathbf{b}]$ associated with $\sim$. One can think of the class $[\mathbf{b}]$ as an infinite lattice of all bonds related to $\mathbf{b}$ by a lattice translation.

We define the action of the point group on the class $[\mathbf{b}]=[(a_{1},a_{2})]$ as

One checks that the action is well-defined in both arguments: the result does not depend neither on the choice of the representative $g$ of the coset, nor on the choice of the bond $\mathbf{b}$ in the equivalence class. This happens because the lattice translations in $[g]=gT$ are absorbed by the equivalence of bonds. Note, however, that it is crucial to “mod out” the lattice translations only as a last step, not for the symmetry action on individual atoms. Keeping this in mind, from now on we will not distinguish between the coset $[g]$ and its representative $g$, to simplify the notation.

The action (17) allows us to introduce the standard notions of a stabilizer and an orbit. Note that the equivalence (16) distinguishes the bonds that have the opposite orientation: $\mathbf{b}=(a_{1},a_{2})$ and $-\mathbf{b}=(a_{2},a_{1})$ belong to different classes. However, an important constraint on the exchange tensors of the bond comes from the symmetries that flip the bond. Taking this into account, we define the stabilizer of the bond $[\mathbf{b}]$ as

In words, the stabilizer subgroup consists of the symmetry elements that, up to a lattice translation, either leave the bond fixed or interchange its atoms. This set of symmetry operations determines constraints on the form of exchange tensors for a given bond class. The orbit

contains all the bonds whose exchange tensors are symmetry-related to those of the starting bond. Note that, for the orbit-stabilizer theorem to hold, the elements of the orbit should also be considered “modulo orientation” to avoid double-counting.

To find the most general form of the exchange matrices for the bonds in the class $[\mathbf{b}]$ and its orbit, we use basic tools from the representation theory of finite groups. Since the DMI matrix $D$, Eq. (7), has three independent matrix elements, one can consider it as an element of a three-dimensional vector space $V_{D}$:

Here, $e_{\alpha}$ are the basis elements in the space of antisymmetric rank 3 matrices. For example,

In a similar way, the anisotropic symmetric matrix $\Gamma$ belongs to the six-dimensional vector space $V_{\Gamma}$ of rank 3 symmetric matrices. Since the action (11) preserves (anti)symmetry of the exchange matrix $J$, it can be interpreted as a linear transformation of the corresponding vector space $V_{J}$:

where $J_{j}$ is a component of a $d$-dimensional vector of matrix elements ($d=3$ for $J=D$ and $d=6$ for $J=\Gamma$). The set of such linear transformations

forms a representation of the stabilizer $\operatorname{St}$ of the bond class $[\mathbf{b}]$, which is a subgroup of the point group $G$.

An exchange matrix $J$ for the bond in the class $[\mathbf{b}]$ is symmetry-compatible if for all elements in the stabilizer we have

This means that $J$ belongs to the invariant subspace of the trivial representation $W_{1}\subset V_{J}$. Thus, our task of finding all symmetry-compatible matrices $J$ amounts to specifying the subspace $W_{1}$. First, we determine its dimension:

If $\dim W_{1}=0$, the anisotropic exchange interaction $J$ is not allowed by symmetry. If $\dim W_{1}=d$, that is, the whole space transforms according to the trivial representation, then the stabilizer of the bond does not put any symmetry constraints on the matrix $J$ (but there may be symmetry relations between the matrices for different bonds in the orbit).

In the intermediate case $0<\dim W_{1}<d$ we need to find $n=\dim W_{1}$ linearly independent vectors in the subspace $W_{1}$. This is done by the projection operator

which projects any vector onto the invariant subspace $W_{1}$ of the trivial representation. Applying $\operatorname{Proj}_{1}$ to the basis elements $\{e_{i}\}$, one can find the desired basis for the space $W_{1}$.

Once the subspace $W_{1}$ is found for the bond class $[\mathbf{b}]$, one can compute symmetry-compatible exchange matrices for all bonds in the orbit of $[\mathbf{b}]$ using the requirement that $H$ be equal to $g\cdot H$ for all pairs of sites. Note that the bond $[\mathbf{b}]$ may be mapped to $g\cdot[\mathbf{b}]$ by more than one group element. One checks that the resulting Hamiltonian $g\cdot H$ does not depend on the choice of this element, by virtue of using $g^{-1}$ in Eq. (5).

For a specific example, let $g$ map the bond $(a_{1},a_{2})$ to the bond $(a_{3},a_{4})$. Then the invariance $g\cdot H=H$ implies that

for each exchange matrix.

Information on chemical composition, cell vectors, and coordinates of atoms is obtained from a cif file, which is processed using the ase library [17]. We also check that the site occupancy is an integer for all atoms with the parsnip-cif package [18].

Then the list of atomic coordinates is processed with the spglib library [19], which finds the space group of the crystal and provides the list of all symmetry elements. Each symmetry element consists of the matrix part and the translation part, both expressed in the basis of the cell vectors. The only control parameter for the symmetry analysis is the numerical precision sym_tolerance, which is used to determine whether coordinates of atoms coincide after applying symmetry operations. Its default value is $10^{-5}$ Å. Note that this parameter might influence the resulting symmetry group in some cases.

For each symmetry element, we compute the Cartesian version of the matrix part. To this end, we construct the basis transformation matrix from the cell vectors. After the transformation, we convert the resulting matrix to a SymPy matrix object, whose elements are exact expressions, such as $\tfrac{\sqrt{3}}{2}$. The list of possible exact values is determined by the fact that a crystal can only have rotation axes of orders $1,2,3,4$ and $6$. We assume that the crystal structure is defined using a conventional unit cell in the standard crystallographic setting. Then the matrix elements can only take values $0$, $\pm 1$, $\pm\tfrac{1}{2}$ and $\pm\tfrac{\sqrt{3}}{2}$.

All geometric data, including coordinates of atoms and symmetry elements, is stored in an instance of class Compound.

Atoms are represented by instances of AtomicSite class. An atom is specified by an index $n$, which enumerates atoms in the unit cell, and by the integer cell coordinates $(m_{x},m_{y},m_{z})$. The value of $n$ encodes both the fractional part of atom’s scaled coordinates $\mathbf{r}_{f}$, which take value in the interval $[0,1)$, and its chemical element. We will refer to the cell with coordinates $(0,0,0)$ as the home unit cell.

We define the Bond object as an ordered pair of AtomicSite objects. Since we are interested in bonds only up to a lattice translation, we define the standard set as a subset of the set of all possible bonds. The bond belongs to the standard set if the following condition is satisfied: in each pair of cell coordinates, $(m_{x}^{1},m_{x}^{2}),(m_{y}^{1},m_{y}^{2}),(m_{z}^{1},m_{z}^{2})$, the smallest coordinate is $0$. In other words, the bonds in the standard set belong to the sector with non-negative cell coordinates and lie as close to the origin as possible. In the language of Sec. 2.2, the standard set is a set of representatives for equivalence classes that correspond to the equivalence relation (16). Moving the bond to the standard set after any transformation allows us to work with the equivalence classes $[\mathbf{b}]$.

Symmetry elements are represented by instances of SymElement class, which contain rotation/reflection matrices $R$ and translation vectors $\mathbf{t}$ in scaled and Cartesian coordinates.

An exchange tensor is represented by an instance of ExchangeTensor class. Its properties include $3\times 3$ SymPy matrix and a list of symbols used in the matrix. There are two types of exchange tensors, Symm and DMI, whose matrices are symmetric and antisymmetric, respectively. Each tensor can be represented as a “vector” of independent matrix elements. For DMI tensor, this vector is $(D_{x},D_{y},D_{z})$, while for the symmetric exchange it reads $(\Gamma_{xx},\Gamma_{yy},\Gamma_{zz},\Gamma_{xy},\Gamma_{xz},\Gamma_{yz})$.

Our code provides functionality for a common task: for a given crystal, find all symmetry constraints on the exchange tensors for bonds between magnetic elements $A$ and $B$, whose length does not exceed a given value $l$.

We start by listing all atoms and their positions for a $3\times 3\times 3$ cluster of cells centered at the home unit cell. Then we use the cKDTree method of the SciPy spatial module to build a list of pairs of atoms $a_{A},a_{B}$ such that

The first atom always lies in the home unit cell, while the second atom can be anywhere in the cluster. Next, we create Bond objects from all pairs, forming a pool of bonds. Our goal is to organize these bonds by their symmetry and length. The result will be the list of starting bonds, each being the first in its orbit. The orbit of each starting bond will be stored as an attribute of the Bond object. This is done as follows:

1. 1.

   Pick any bond $\mathbf{b}$ from the pool and compute its symmetry orbit. Compute the exchange tensors for all bonds in the orbit. Add $\mathbf{b}$ to the list of bonds.

2. 2.

   Exclude the orbit elements from the pool of bonds.

3. 3.

   Repeat until the pool is empty.

Finally, we sort the bonds in the list by their length.

Given a cif file of a crystal, one can quickly analyze the exchange tensors using the package in a Python console. Below we show a typical workflow, using hematite $\alpha$-Fe₂O₃ as an example.

We start by initializing the Compound object from a cif file.

The output indicates that the Compound was initialized successfully.

Now one can find all bonds between Fe atoms, whose length does not exceed 5Å:

Here, an atom is represented by a label, which indicates the atom’s position in the unit cell, and a triple of integer unit cell coordinates. The list contains starting bonds of the orbits. Also shown are the length of each bond and the number of bonds in its orbit. The columns DMI and Symm show the numbers of independent matrix elements in the exchange matrices $D$ and $\Gamma$. One can readily conclude that the DMI exchange in the second and in the last orbit is prohibited by symmetry.

Exchange tensors for the bonds from the first orbit are

The first table indicates that the DMI vector of the bonds belongs to the vertical axis. To interpret the signs correctly, one should take into account the orientation of the bonds:

To orient the bonds uniformly, we need to flip the first and fourth bonds, which will change the signs of $\mathbf{D}$.

Now consider the fourth orbit. It contains 36 bonds, which equals the order of the group. This means that each of these bonds is stabilized only by the identity transformation. Still, symmetry puts constraints on the exchange matrices:

Thus, while each DMI vector is described by three independent variables, the vectors are not at all arbitrary, since they are related by the symmetries that interchange the bonds in the orbit.

The package provides a web interface that can be run locally [20] or accessed on the web [21]. The user can upload a cif file and analyze it with a chosen numerical precision. Alternatively, one can load an example file from the Help page. Once the file is loaded, the interface shows crystal geometric data, such as position of atoms and representation of symmetry elements. The user can then choose a pair of chemical elements and a maximal bond length $L$, and obtain a list of all bonds shorter than $L$ classified into symmetry orbits.

For a chosen orbit, the interface shows a table of symmetry-compatible exchange tensors. For any bond, the user can further address its geometric details (coordinates of atoms), its stabilizer subgroup, and interact with a 3D model showing the placement of the bond inside the unit cell.

We start with the famous La₂CuO₄, which is the parent compound for cuprate high-temperature superconductors. In the overdoped regime, it exhibits magnetic order with spins directed along orthorhombic $a$ axis, but canted out of the CuO₂ plane (see Fig. 1 in  [22]). A detailed structural study has shown that the crystal is orthorhombic in this phase, with $Cmcm$ space group [23]. Experimentally spins are mostly lying in the $ab$ plane and point nearly two the second nearest neighbor Cu ions, but have a small canting out of the plane.

Analysis of the exchange tensor using Jsymm immediately reveals a non-zero Dzyaloshinskii–Moriya interaction (DMI), with the $\mathbf{D}$ vector having two independent components for the nearest-neighbor Cu–Cu bonds lying in the $ab$ plane (distance is 3.8Å), as shown in Table 1. This result is in agreement with microscopic calculations of the DMI, such as those presented in [24, 25].

The presence of a $C_{2}$ rotational axis directed along the $c$ axis (i.e. perpendicular to each of these bonds) and passing through their centers forces the $z$ component of the DM vectors to vanish. Consequently, the DM vectors for nearest-neighbor Cu ions possess only $x$ and $y$ components, which naturally explains the canting of spins out of the $ab$ plane observed experimentally [22].

Jsymm also demonstrates relations between exchange tensors corresponding to different bonds that belong to the same orbit. For example, consider the DM vectors $\mathbf{D}^{1}$ and $\mathbf{D}^{3}$ of the bonds 3-$2_{1}$ and $2_{3}$-3 formed by Cu ions, as shown in Fig. 1. Their components are related as

This agrees with the results of Ref. [25] (if one flips the bond $2_{3}$-3, so that $\mathbf{D}^{3}$ changes its sign) and Ref. [24] (coordinate system must be rotated by $\tfrac{\pi}{4}$ about the $z$ axis for comparison).

Another classical example of a material with anisotropic exchange coupling is the mineral hematite, $\alpha$-Fe₂O₃, which has a corundum structure described by the $R\bar{3}c$ space group [27]. Below $T_{N}\sim 950$K, but above the so-called Morin temperature $T_{M}\sim 250$K — the temperature of the spin reorientation transition — it orders antiferromagnetically with moments perpendicular to the [111] direction, but with a small canting angle of $\sim 10^{-3}$ within the basal plane [28, 29]. Using phenomenological theory and symmetry arguments, Dzyaloshinskii was the first to explain this canting [30]. This theory was later elaborated by Moriya using microscopic considerations [16]. In this work, he also introduced the symmetry principles known as the “Moriya rules”, which correspond in our terms to the constraints put on an individual bond by its stabilizer subgroup.

Our analysis of the crystal structure taken from [31] demonstrates that for nearest-neighbor Fe ions forming a short Fe–Fe bond (distance is 2.89Å), the DM vector is $(0,0,D_{z})$, while $\Gamma$ is allowed to have only two terms: $\Gamma_{xx}=\Gamma_{yy}$ and $\Gamma_{zz}$. DMI with a nonzero $z$ component of the DM vector causes spins lying perpendicular to the $c$ axis to cant within the (111) plane in rhombohedral notation. Exchange matrices for the next neighbors are given in Eq. (33), which was generated by the code listed in Sec. 4.3.5:

The degree and character of the symmetry constraints vary among the bonds. While the anisotropic symmetric exchange $\Gamma$ can be present to certain extent in all bonds, only the first, third, and fourth bonds are allowed to have DM exchange interaction. The fourth bond has the lowest possible symmetry, so its matrix has a generic form. Still, there are the symmetry restrictions for the exchange tensors of the bonds in its orbit, as we saw in Sec. 4.1.

These matrices are in complete agreement with the recent results presented in Ref. [32] (see Eq. (25) there). The apparent discrepancy for the exchange matrices of the third bond can be removed by the rotation $U_{\alpha}$ of the Cartesian coordinate system by $\alpha=\tfrac{\pi}{12}$ about the $z$ axis and by selecting an appropriate bond. In our coordinates, its bond vector in Angstroms reads $[2.5173,1.4533,1.6957]$. Conjugation of its exchange matrix by $U_{\alpha}$ produces a matrix whose elements are related to each other precisely as in Ref. [32]. This illustrates that the form of the symmetry-adapted exchange tensors may vary depending on the relative orientation of the bond vector and the axes of the Cartesian coordinate system.