Spin wave excitations of magnetic metalorganic materials

The Organic Materials Database (OMDB) is an open database hosting about 22,000 electronic band structures, density of states and other properties for stable and previously synthesized 3dimensional organic crystals. The web interface of the OMDB offers various search tools for the identification of novel functional materials such as band structure pattern matching and density of states similarity search. In this work the OMDB is extended to include magnetic excitation properties. For inelastic neutron scattering we focus on the dynamical structure factor S(Q, ω) which contains information on the excitation modes of the material. We introduce a new dataset containing atomic magnetic moments and Heisenberg exchange parameters for which we calculate the spin wave spectra and dynamic structure factor with linear spin wave theory and atomistic spin dynamics. We thus develop the materials informatics tools to identify novel functional organic and metalorganic magnets.

The Organic Materials Database (OMDB) is an open database hosting about 22,000 electronic band structures, density of states and other properties for stable and previously synthesized 3dimensional organic crystals. The web interface of the OMDB offers various search tools for the identification of novel functional materials such as band structure pattern matching and density of states similarity search. In this work the OMDB is extended to include magnetic excitation properties. For inelastic neutron scattering we focus on the dynamical structure factor S(Q, ω) which contains information on the excitation modes of the material. We introduce a new dataset containing atomic magnetic moments and Heisenberg exchange parameters for which we calculate the spin wave spectra and dynamic structure factor with linear spin wave theory and atomistic spin dynamics. We thus develop the materials informatics tools to identify novel functional organic and metalorganic magnets.

I. INTRODUCTION
Magnetism and magnetically ordered materials have played a crucial role in the development of the technology used in our every-day life. The identification of novel materials with desired magnetic target properties as well as the investigation of coupling mechanisms, the resulting order and its excitations are therefore of great importance. While basic concepts are usually explored in materials with feasible complexity, materials with complex unit cells and dominant interaction effects quite often exhibit the more desirable properties with respect to technological applications. In particular metal organic frameworks and organic molecular crystals exhibit promising structures for electron-spin-based devices. Magnetic organics have attracted attention with respect to spintronic devices [1][2][3] and magnon spintronics [4], multiferroic phases [5,6], molecular qubits [7,8], and spin-liquid physics [9][10][11]. Local magnetic moments in organic materials can arise due to transition metal and rare earth ions embedded in the molecules or due to local unsaturized bonds as they occur in stable organic radicals [12][13][14]. Some organics establish stable magnetic order of various kind up to room temperature [15][16][17]. For device engineering and technological applications, an important advantage of organic materials over other functional materials such as transition metal oxides, is that organics can be synthesized with relative low cost and large-scale production methods.
The vast increase of experimental and theoretical data obtained over the past century has opened a novel approach to materials research based on computer science methods and the construction of materials databases [18][19][20][21][22][23][24][25]. Such databases were successfully applied in mining for functional materials [21,[26][27][28][29][30][31][32] and as training * hellsvik@kth.se data for machine learning algorithms predicting complex materials properties [33][34][35][36], bypassing computationally demanding ab initio calculations. In this article we report about the implementation of a novel dataset for organic magnets into the organic materials database OMDB [19,37], available at https://omdb.mathub.io. The data was calculated by means of ab initio calculations and comprises information about local magnetic moments, magnetic exchange coupling, expected magnetic ground state, as well as spin-wave excitation spectra and the dynamical structure factor S(Q, ω) for hundreds of previously synthesized organic molecular crystals and metal organic frameworks. This data is embedded into the existing framework of the OMDB and can be explored using interactive statistics and non-trivial search tools such as pattern matching [38,39].
The implementation of such a database is closely related to the recent enhancement of neutron flux [40][41][42] and detector technology, allowing for inelastic neutron scattering experiments to be performed with higher signal to noise ratio, and higher energy resolution [43,44]. For example, for inelastic neutron scattering experiments, the central entity is the dynamical structure factor S(Q, ω), which contains information on the excitation modes of the material. The vast majority of inelastic neutron scattering (INS) measurements of magnetic materials are analyzed by means of fitting the experimental dynamical structure factor S(Q, ω) to the dispersion and structure factor of a spin wave Hamiltonian as provided in our dataset.
The remaining of the article is organized as follows.
In Section II we present the models and approximations used, as well as the scheme developed for high throughput calculation of magnetic properties and excitation spectra. Results are presented in Sec. III, followed by a conclusion and outlook in Sec. IV. The workflow for high throughput calculation of magnetic ground states and spin wave spectra. The blue boxes indicate the steps we use for the results presented in this article, with the figures and equations in the right hand column illustrating central entities. The black arrows on the left hand side indicate steps which can be performed with machine learning (ML) techniques, to replace computationally expensive DFT and ASD simulations, ultimately enabling for prediction of magnon spectra (green TARGET box) directly from the crystal structure (red INPUT box).

II. MODELS AND METHODS
A state-of-the-art approach for ab initio investigations of spin dynamics in magnetic solids is to solve the full electronic dynamics with time-dependent density functional theory (TD-DFT) and calculating the energies and lifetimes of the spin waves without resorting to the adiabatic approximation [45,46]. However, both the real time and the linear response variants of TD-DFT are computationally demanding and currently not feasible for organic materials with large chemical unit cells containing in average ≈ 80 atoms for the materials stored within the OMDB. The calculation of the coupling parameters for magnetic Hamiltonians is nowadays supported by various DFT softwares, such as RSPt [47], SPR-KKR [48], KKRnano [49], and HUTSEPOT [50,51]. This opens the path for a more tractable approach in context of the highthroughput computations performed by us, allowing for an ab initio modeling of spin wave excitation spectra in form of a two-step approach: first, the coupling parameters of an effective magnetic Hamiltonian are calculated from DFT calculations; second, analytical or numerical spin wave theory calculations are performed to obtain the spin wave spectra. Furthermore, thermodynamic properties of the spin Hamiltonian can be investigated with Monte Carlo and atomistic spin dynamics (ASD) simulations [52].
The description of our methodology is contained in a number of sections. The magnetic model Hamiltonian and details about the density functional theory calculations go into Sec. II A, followed by a brief description of linear spin wave theory for collinear magnets in Sec. II B, and the atomistic spin dynamics method in Sec. II C. The high throughput workflow is described in Section II D and is shown schematically in Fig. 1.

A. Magnetic Hamiltonian
For a magnetic solid with N A atoms in the crystallographic unit cell and a total of N c cells, we use the notation R iµ = R i + r µ to specify atomic positions, where r µ is the position of the basis atom µ in the unit cell, and R i is the position of the crystallographic unit cell i. The magnetic model is a Heisenberg Hamiltonian formulated in terms of unit vectorsê iν for the directions of the local magnetic moment at site iµ. Note the sign convention with a leading minus sign for the first term, and that the double summation NcN A jν count each bond twice. The magnetic moment of an atom is µ iµ = µ B m µêiµ , where m µ is the size of the magnetic moment of atomic type µ in terms of multiples of Bohr magnetons µ B .
In order to formulate spin Hamiltonians for magnetic organic materials we have used the full-potential linear muffin-tin method (FP-LMTO) software RSPt [47] to calculate Heisenberg exchange parameters. The fullpotential basis set allows for accurate electronic calculations irrespective of the geometry of the crystal structure. The latter is important for the sparse unit cells typically present in organics. Furthermore, in RSPt Heisenberg exchange interactions can be calculated from Green functions for a reference spin structure, e.g. a ferromagnetic or collinear antiferromagnetic ordering, by means of the Liechtenstein-Katsnelson-Antropov-Gubanov (LKAG) formalism [53], without the need to set up supercells with different spin configurations. Details on the implementation of the LKAG formula in RSPt can be found in [54].
The calculations were made with the local density approximation (LDA) for the exchange-correlation potential. The charge density and the potential inside the muffin-tin spheres were represented using an angular momentum decomposition up to l max = 8. One energy set was used for the valence electrons. For the description of the states in the interstitial region three kinetic energy tails were used: -0.3, -2.3, and -1.5 Ryd.

B. Adiabatic magnon theory
The excitation spectra of a spin Hamiltonian can be obtained by various analytical and semi-analytical methods. The traditional approach relies on the introduction of Holstein-Primakoff (HP) operators for the bosonic operators, and expanding the spin Hamiltonian to quadratic or higher order in these HP operators. The expansion to quadratic order leads to linear spin wave theory which has established itself as the main formalism to model inelastic neutron scattering (INS) data of magnetic materials. This technique and aspects of how it can be implemented in software is described for the general noncollinear case of multisublattice magnets in [55][56][57]. In the following we briefly review the main steps of the so called adiabatic magnon theory, a formalism that naturally connects to the frozen magnon technique for obtaining spin wave energies from ab initio calculations. More details can be found in [58]. The spatial Fourier transform of the exchange interaction is defined according to where R 0µjν = R 0 +r µ −R j −r ν , and i = 0 is the central unit cell. After linearizing the equation of motion for the Heisenberg Hamiltonian in Eq. (1) for small cone angle excitation relative to the magnetization quantization direction, here chosen to be alongẑ, we can introduce the quantitỹ in which e z µ = {1, −1} specify the collinear (anti)ferromagnetic groundstate of the system. The spinwave dispersion as a function of wave vector q can then be obtained by diagonalizingJ µν (Q).

C. Atomistic spin dynamics
In common with spin wave theory a core element in the atomistic spin dynamics method [52,59] is that it is possible to parametrize the energies and dynamics of the magnetic system to a magnetic Hamiltonian formulated in terms of local spins (or magnetic moments) as in Eqn.
(1). The other core element is the stochastic Landau-Lifshitz-Gilbert equation (SLLG), which is a semi-classical equation used to model the motion of the atomic magnetic moments at zero or finite temperature. The first term of the equation is the precessional motion, while the second term describes the damping motion. m i = m i (t) is the magnetic moment and experiences an effective magnetic field B i = B i (t), calculated from the Hamiltonian as B i = − ∂H ∂mi . γ is the gyromagnetic ratio, α is a scalar (isotropic) Gilbert damping constant. Temperature is included in the form of Langevin dynamics with the power of the stochastic magnetic field B fl i = B fl i (t) related to the damping constant α via a fluctuation-dissipation relation [52]. For the ASD simulations we use the UppASD software [60].
In numerical simulations of the stochastic Landau-Lifshitz-Gilbert equation there is no need for a linearization of the equations of motion wherefore the full dynamics of a magnetic system described by a bilinear-inspin-operators Hamiltonian such as Eqn. (1) is retained [52,61]. Furthermore, also the full dynamics of a magnetic Hamiltonian which contains higher order coupling such as biquadratic exchange can be studied without the need for mean-field treatment of these couplings, see e.g. Ref. [62]. Spatial and temporal fluctuations are sampled using a time-dependent correlation function where α and β are the Cartesian components of the magnetic moments. Fourier transforming C(r, t) over space and time yields the dynamic structure factor which can be related to the scattering intensity measured in inelastic neutron or electron scattering on magnetic materials. Only the magnetic excitations perpendicular to the scattering wave vector contribute to the scattering intensity [57,63] S(Q, ω) = α,β For the measurement of the dynamic structure factor ASD simulations were run at T = 1 K using a time step dt = 5 · 10 −16 s and a small Gilbert damping α = 0.0001. The correlation function in Eq. (5) was measured using a sampling step of t samp = 5 · 10 −15 s, over a moving sampling window of t win = 5 · 10 −11 s. The corresponding frequency range for the dynamic structure factors in Eqs. (6) and (7) is ω/(2π) = [0.02, 0.04, . . . , 200] THz (0.0827 meV to 827 meV).

D. High throughput calculation
In this section is described the methodology that we have developed for high throughput calculation of magnetic ground states and spin wave spectra. The main steps are shown in Figure 1.
Input files for the RSPt ground state calculations are prepared from the crystallographic information files (CIF) using the cif2cell program [64]. The calculations are performed for a ferromagnetic ground state configuration, without any a priori consideration of what the magnetic ground state could be. In order to classify the atoms as magnetic or non-magnetic, a simple criteria is used, namely that the spin-polarized charge density integrated over the muffin-sphere of the atomic sites is larger than 0.1 µ B . For these atomic sites Heisenberg exchange interactions are calculated using the approach described in Sec. II A.
The linear spin wave calculation require as input not only the parameters of the magnetic Hamiltonian, but also ground state spin configuration. In order to search for the ground state we use an atomistic spin-dynamics quenching scheme for simulation cells with edge length L, corresponding to a number N atom = N A L 3 of sites, where N A is the number of magnetic atoms in the primitive chemical cell. Starting from having all magnetic moments initially in a ferromagnetic configuration, the system evolves in Langevin dynamics simulation at finite temperature. This is followed by a sequence of simulations at zero temperature but with finite Gilbert damping α. For zero or static external magnetic field, the deterministic Landau-Lifshitz equation has the property that energy is a non-increasing function of time, and consequently that the energy of the spin Hamiltonian will be minimized [52]. For systems with competing magnetic phases at low temperatures, there is a finite probability that the system will get caught in a local minimum so that the true ground state of the magnetic Hamiltonian is not found, an issue that could be of particular concern for systems with highly degenerate (free) energy landscapes such as frustrated magnets and spin glasses [65]. We do not have quenched chemical disorder wherefore a spin glass phase is not possible, furthermore a pure Heisenberg Hamiltonian as in Eq. 1 cannot stabilize a skyrmion spin texture. For the present data set, we classify the systems as ferromagnetic, collinear antiferromagnetic, noncollinear antiferromagnetic, or param-agnetic. An overview and statistics of the dataset will be presented in Section III. Another important distinction is whether the magnetic ordering is commensurate with the primitive chemical cell or not. Knowledge of the shape and size of the primitive magnetic cell, allows for the linear spin wave calculation to be performed for the size of the spin wave Hamiltonian that will give the correct number of magnon dispersion eigenvalues and eigenvectors. In order to obtain from an size N A L 3 simulation cell the smallest possible magnetic cell, we make use of the capability of the Elk software [66] to reduce from an input crystallographic structure to the smallest primitive cell, considering also the magnetic ordering.

III. RESULTS
The main dataset of the OMDB constitute of spinpolarized electron band structures and density of states [19,37]. As a very first step in the production of our current set of magnetic Hamiltonians and excitation spectra, we have chosen a subset of magnetic materials and did for these deploy the high throughput methodology described in Sec. II D. We have so initially obtained the magnetic excitation spectra for slightly more than 100 materials which all have in common that each of them, in the RSPt LDA calculations, displayed one or more sites with a magnetic moment of at least 0.1 µ B . The majority of the materials in this magnetic excitations dataset have noncollinear ground states. Rather remarkably, none of our currently considered materials displays a ferromagnetic groundstate. The ASD quenching simulations to obtain magnetic ground states used simulation supercells with edge length L = 6, and the ASD simulations for sampling S(Q, ω) at finite temperature used supercells with L = 20.
We will in the following present our magnetic excitations dataset (MagExcite), closely connecting to how the results are presented on the Organic Materials Database. We refer to the materials using their ID numbers (COD-ID) within the crystallographic open data base (COD) [67,68], and the OMDB-ID numbers.
Our first example is the compound C 20 H 10 CoN 12 (OMDB-ID 11913, cod-id 2014058) which has very low symmetry and crystallizes in spacegroup number 2. The primitive nonmagnetic cell has 43 atoms, corresponding to one formula unit. The multiplicity of the Co Wyckoff position is 1. For this site we obtained a magnetic moment µ Co = 0.60µ B . Heisenberg interactions were calculated between the Co sites for distances up to a maximum of three multiples of the lattice constants. The outcome of the ASD quenching calculation was a collinear antiferromagnetic spin configuration that is shown in Fig.  2. With regard to the original lattice vectors for the primitive chemical cell, a 2 × 1 × 2 supercell of the primitive cell can accommodate the antiferromagnetic ordering. Calculating the spin wave spectra for this conventional magnetic cell, we obtain two doubly degenerate magnon bands which are shown in the lower left corner of Fig. 2. The conventional magnetic cell can be reduced to a two-atom magnetic cell for which the magnon bands are shown in the upper left corner of Fig. 2.
A snapshot of the OMDB magnetic properties panel for C 20 H 10 CoN 12 is shown in Fig. 3. At the top of the panel is shown the T = 0 K adiabatic magnon spectra along a high symmetry path in the Brillouin zone, with an accompanying interactive figure for the special points in the Brillouin zone. For the same high symmetry path the corresponding T = 1 K dynamical structure factor intensity, sampled using Eqs. (5)- (7), is displayed in a colorplot with logarithmic scaling for the intensity. The exchange couplings are listed, ordered to show the strongest couplings first, in a table in the lower left part of the magnetic properties panel and are also displayed in the JSmol viewer. A listing of the magnetic moments can be found in the lower right part of the panel.
Also the compound C 10 H 10 N 8 Ni (OMDB-ID 19963, cod-id 2216760) has a structure that belongs to space group number 2. The primitive nonmagnetic cell has 58 atoms, corresponding to two formula units. For the Ni site we obtain a magnetic moment µ Ni = 1.05µ B . The magnetic ground state is collinear, and can be accommodated in a 2 × 1 × 2 supercell of the primitive cell. Given that the primitive nonmagnetic cell has two Ni sites, the supercell contains a total of eight magnetic sites. Similarily as for C 20 H 10 CoN 12 , the supercell can be reduced to a primitive magnetic cell, this time with four atoms. The linear spin wave calculations yield two doubly degenerate magnon bands, one pair being an acoustic mode, and the other pair being an optical mode with finite values of the dispersion at the zone center. Our results for C 10 H 10 N 8 Ni are displayed in Fig. 4. The color plot for the dynamic structure factor gives an indication on the relative occupancy of higher and lower energy magnons at the given temperature of T = 1 K [69].
As a third example, we will discuss the noncollinear antiferromagnet material manganese dicarboxylate C 2 H 3 MnO 3 [70] (OMDB-ID 11283, cod-id 1515486), with results displayed in Fig. 5. The material crystallizes in the non-symmorphic space group number 14. The primitive chemical cell contains 36 atoms, corresponding to four formula units. A magnetic moment is carried by the Mn sites for which our DFT calculations gave a magnetic moment µ Mn = 2.20µ B . The outcome of the ASD quenching simulation was a canted antiferromagnetic phase. Our current implementation of linear spin wave theory does not support fully automated calculations of the spin wave spectra for systems with noncollinear ground states. In contrast, the ASD simulation method allows for semiclassical sampling of the dynamic structure factor for magnetic systems with any spin ordering, including also the paramagnetic regime where long range ordering is absent. Figure 5 shows the dynamical structure factor at T = 1 K for C 2 H 3 MnO 3 , revealing two main bands of dispersion.
Dispersion relations for magnons are commonly calcu- lated along lines in reciprocal space between high symmetry points in the Brillouin zone. We have chosen to use for each material the same paths for magnons dispersions as has been used for the creation of the OMDB electronic band structure data set [19,37]. Similar to the previously described pattern matching for electronic band structures [38], we have implemented a pattern matching algorithm for magnon bands. The algorithm is based on representing pairs of bands as vectors within a chosen moving window. The similarity between defined pattern and magnon bandstructure is obtained by calculating the cosine distance. To speed up the approach for online usage, an approximate nearest neighbor algorithm is used. Figure 6 shows its application for the example of identifying Dirac nodes within magnon spectra. Using as a query Dirac-type crossing of linear bands, the pattern matching algorithm here identified such a crossing along a line in reciprocal space between the high symmetry points Z and D, for the space group 14 collinear antiferromagnet C 6 H 14 NiO 8 (OMDB-ID 21042, cod-id 4325324).

IV. CONCLUSIONS AND OUTLOOK
In summary we have developed a method for highthroughput calculation of the magnetic excitation spectra of organometallic materials. To this end we have deployed a multiscale modeling approach in which first the atomic magnetic moments and interaction coupling constants of magnetic Hamiltonians are calculated, followed by ground state determination, and calculation of the spin wave dispersion and the dynamic structure factor. For the present results we have used Heisenberg Hamiltonians with interactions ranging up to three lattice constants. Frustration among the Heisenberg interactions can lead to noncollinear magnetic ground states, and we have among our materials encountered collinear antiferromagnets as well as noncollinear spin orderings. Currently we have used an implementation of linear spin wave theory for collinear magnets, wherefore dispersions were calculated only for the collinear systems, however, we have sampled the dynamic structure factor at T = 1 K for all materials using atomistic spin dynamics.
The method can be naturally extended to include other interactions such as Dzyaloshinskii-Moriya interaction, Ising interaction, and single-site magnetocrystalline anisotropy energy to the magnetic Hamiltonian, of relevance for the low-energy excitation spectra of quantum organometallic materials. Related to this we expect an even higher fraction of noncollinear ground states, for which the need for calculation of the spin wave dispersions with the more general framework of linear spin wave theory for noncollinear magnets is desired [56,57]. Furthermore, recent developments of machine learning techniques for lattice models and spin Hamiltonians, as for instance a profile method for recognition of three-dimensional magnetic structures [71], determination of phase transition temperatures by means of self-organizing maps [72], and a support vector machines based method for multiclassification of phases [73], will be most useful for identification and classification of competing magnetic phases at finite temperature, and the corresponding phase transition temperatures.
A high-throughput database of magnetic properties of organometallic materials opens the path towards establishing general concepts of materials statistics. We combine our developed search tools for materials properties [19,39] with the new data for magnetic properties, enabling us to search for user-specified patterns in the magnon dispersion relations and magnon density of states. Such tools provide a novel approach towards identifying functional magnetic materials, such as magnon Dirac materials [74][75][76], topological magnon insulators [77][78][79], and magnon Hall materials [80,81]. While there is a huge presence of theoretical concepts of such novel topological magnon phases, only very few materials are known exhibiting these bosonic quantum phases. We see strong indications that the tools presented by us during this article will help us to efficiently explore the realm of organic materials in this context. With the database grown to a decent size we furthermore see great opportunities in training machine learning models. We have recently shown that such machine learning tools can be applied to predict basic materials properties for extremely complex organic materials hosting thousands of atoms in the unit cell and with that are outside the realm of materials which can be calculated straightforwardly using density functional theory [36]. In a similar manner, we  see great opportunities in training models towards automatically formulating effective Heisenberg models for arbitrary organic materials, significantly accelerating the accumulation of magnon and S(Q, ω) spectra, as needed, for instance by large-scale experimental facilities [40][41][42].

V. ACKNOWLEDGEMENTS
We acknowledge economic support from the Swedish Research Council (VR) through a neutron project grant (BIFROST, Dnr. 2016-06955), and by VILLUM FONDEN via the Centre of Excellence for Dirac Materials (Grant No. 11744). The authors acknowledge computational resources from the Swedish National Infrastructure for Computing (SNIC) at the Center for High Performance Computing (PDC), the High Performance Computing Center North (HPC2N), and the Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX).