Optimizing configurations for determining the magnetic model based on ab-initio calculations

TL;DR

This paper introduces a new strategy to optimize magnetic coupling determination using ab-initio calculations, Monte Carlo sampling, and model refinement, demonstrated on a specific material with an open-source software implementation.

Contribution

It presents a novel method combining ab-initio calculations, Monte Carlo sampling, and model refinement for efficient magnetic coupling estimation.

Findings

Effective magnetic configurations identified for Bi₃Mn₄O₁₂(NO₃)

Coupling constants and error bounds estimated accurately

Open source software automates the entire process

Abstract

In this paper, it is presented a novel strategy to optimize the determination of magnetic couplings by using ab-initio calculations of the energy. This approach allows determining efficiently, in terms of a proposed effective magnetic spin model, an optimal set of magnetic configurations to be simulated by DFT methods. Moreover, a procedure to estimate the values of the coupling constants and their error bounds from the estimated energies is proposed. This method, based on Monte Carlo sampling, takes into account the accuracy of the ab - initio simulations. A strategy to refine models reusing previously computed configuration energies is also presented. We apply the method to determine a magnetic model for the recently synthesized material Bi$_3$Mn$_4$O$_{12}$(NO$_3$). Finally, an open source software that implements and automatizes the whole process is presented.

Tables

Table 1. Table 1: Coupling constants obtained by ab-initio calculations with fifty-four configurations in ref Alaei  et al. ( 2017 )

$J_i/\left|J_1\right|$	$\frac{\text{11 configurations}}{\left(\left|J_1\right|=0.346\mathrm{meV}\right)}$	$\frac{\text{54 configurations}}{\left(\left|J_1\right|=0.349\mathrm{meV}\right)}$
$J_1/\left|J_1\right|$	$-1.0\pm 0.1$	$-1.0\pm 0.06$
$J_2/\left|J_1\right|$	$-0.12\pm 0.06$	$-0.11\pm 0.04$
$J_3/\left|J_1\right|$	$-0.09\pm 0.07$	$-0.09\pm 0.05$
$J_0/\left|J_1\right|$	$-0.3\pm 0.21$	$-0.3\pm 0.12$
$J_{1c}/\left|J_1\right|$	$-0.1\pm 0.09$	$-0.11\pm 0.06$
$J_{2c}/\left|J_1\right|$	$-0.05\pm 0.07$	$-0.03\pm 0.04$
$J_{3c}/\left|J_1\right|$	$-0.07\pm 0.08$	$-0.06\pm 0.05$

Equations

$$|\beta ⟩=|\{{\phi }_i\},\{s_i^{\beta }\}{⟩}_{SL}=\frac{1}{\sqrt{N}}\left||{\phi }_1,s_1^{\beta }⟩\dots |{\phi }_N,s_N^{\beta }⟩\right|$$

$${\mathbf{H}}_0=\sum _{\beta {\beta }^{\prime }}|\beta ⟩⟨\beta |{\mathbf{H}}_{\mathrm{full}}|{\beta }^{\prime }⟩⟨{\beta }^{\prime }|$$

$${\mathbf{H}}_{eff}[\{J_{\alpha }\}]=\sum _{\alpha =1}^M\frac{J_{\alpha }}{2}\sum _{(i,j)\in B_{\alpha }}{\vec{\mathbf{S}}}_i\cdot {\vec{\mathbf{S}}}_j+J_{0}1$$

$${\mathcal{C}}_{\Delta \epsilon }:=\{\{J_{\alpha }\}/|⟨\beta |{\mathbf{H}}_{eff}[\{J_{\alpha }\}]|\beta ⟩-E_{\beta }|<\Delta \epsilon ,\forall |\{s_i\}⟩\}.$$

$${\mathcal{C}}_{\Delta \epsilon }=\{\vec{J}/\Vert A\cdot \vec{J}-{\vec{E}}_0{\Vert }_{\infty }<\Delta \epsilon \},$$

$$[A{]}_{\beta \alpha }=\frac{1}{2}\sum _{(i,j)\in B_{\alpha }}⟨\{s_i{\}}_{\beta }|{\vec{S}}_i\cdot {\vec{S}}_j|\{s_i{\}}_{\beta }⟩$$

$$\Vert \vec{v}{\Vert }_{\infty }=\underset{i}{\max }|{\vec{v}}_i|$$

$${\mathcal{C}}_{\Delta \epsilon }({\mathcal{S}}_n):=\{\vec{J}/\Vert A^{\prime }\cdot \vec{J}-{\vec{E}}_0^{\prime }{\Vert }_{\infty }<\Delta \epsilon \}$$

$$\frac{\Vert \vec{x}{\Vert }_{\infty }}{\sqrt{n}}\le \frac{\Vert \vec{x}{\Vert }_2}{\sqrt{n}}\le \Vert \vec{x}{\Vert }_{\infty }\le \Vert \vec{x}{\Vert }_2$$

$${\mathcal{C}}_{\Delta \epsilon }^{(2)}({\mathcal{S}}_n):=\{\vec{J}/\Vert A^{\prime }\cdot \vec{J}-{\vec{E}}_0^{\prime }{\Vert }_2<\Delta \epsilon \},$$

$${\chi }^2(\vec{J})=\Vert A^{\prime }.\vec{J}-{\vec{E}}_0^{\prime }{\Vert }^2$$

$$A^{\prime }=U\Sigma V^t$$

$$\Vert {\mathcal{C}}_{\Delta \epsilon }^{(2)}({\mathcal{S}}_n)\Vert <\Delta \epsilon cn({\mathcal{S}}_n)$$

$$cn(A^{\prime })=\underset{s\in \sigma (A^{\prime })}{\max }\frac{1}{s}$$

$${\mathcal{C}}_{\Delta \epsilon }^{(2)}({\mathcal{S}}_n)\subset {\mathcal{C}}_{\Delta \epsilon }({\mathcal{S}}_n)\subset {\mathcal{C}}_{\Delta \epsilon \sqrt{n}}^{(2)}({\mathcal{S}}_n)$$

$$C({\mathcal{S}}_n):=\sqrt{n}cn(A^{\prime }({\mathcal{S}}_n))$$

$$J_{\alpha }^{LS}={\mathrm{argmin}}_{J_{\alpha }}{\chi }^2(J_{\alpha })$$

$$|1⟩$$

$$|2⟩$$

$$|10⟩$$

$$|17⟩$$

$$|18⟩$$

$$|24⟩$$

$$|28⟩$$

$$|34⟩$$

$$|41⟩$$

$$|45⟩$$

$$|47⟩$$

$$A$$

$$U^{†}U$$

$${\vec{x}}_0={\mathrm{argmin}}_{\vec{x}\in {\mathbb{R}}^n}\Vert A\cdot \vec{x}-\vec{b}{\Vert }^2$$

$$A^{+}=(A^{t}A{)}^{-1}{A}^t=A^t(A{A}^t{)}^{-1}=V{\Sigma }^{-1}{U}^t$$

$$A{A}^{+}$$

$$A^{+}A$$

$$\Vert A(\delta \vec{x}+{\vec{x}}_0)-\vec{b}{\Vert }^2$$

$$\Vert A\cdot \vec{x}-\vec{b}\Vert <\epsilon$$

Full text

I Introduction

The physical description of magnetic degrees of freedom in a broad class of compounds is usually based on the well known Heisenberg-Dirac-van Vleck Hamiltonian Auerbach (1994); Lhuillier and Misguich (2011); Balents (2010); Lamas et al. (2012); Bishop et al. (2012); Lamas, Capponi, and Pujol (2011); Gong et al. (2013). In this context, the use of ab-initio techniques of electronic structure to determine the effective values of the exchange couplings ($J_{i}$) is a very useful tool for theoretical physicists. This theoretical description, despite its simplicity, allows us to understand the basic ingredients that give rise to a wide variety of magnetic phases Cabra, Lamas, and Rosales (2011a); Messio, Bernu, and Lhuillier (2012); Honecker, Mila, and Troyer (2000); Mosadeq, Shahbazi, and Jafari (2011); Mulder et al. (2010); Ceccatto, Gazza, and Trumper (1993); Lamas and Matera (2015); Brenig and Becker (2001); Chubukov and Jolicoeur (1991); Zhang, Arlego, and Lamas (2014); Lamas et al. (2013). The successful determination of the magnetic parameters lies in the appropriate balance between the selection of the model and the spin configuration used in the ab-initio calculations. In contrast, the super-exchange pathways via bridging ligands may cause that interactions between sites separated by a long distance to be large, forcing us to take very large unit cells. The computational cost is greatly increased with the unit cell size and the use of a large number of spin configurations is restrictive. In order to be able to perform calculations of the Heisenberg exchange coupling constants and determine the minimal model we should be able to detect the set of configurations that gives us the best determination of our model. In this paper we present a strategy to determine this set of optimal configurations to be used in the ab-initio estimation of the energy. In this way we are able to obtain a well conditioned system of linear equations to determine the parameters of the magnetic model.

We apply the algorithms to determine the magnetic couplings corresponding to the compund Bi3Mn4O12(NO3)Matsuda et al. (2010). In this compound the Mn4+ ions form a honeycomb lattice. Two layers of such honeycomb planar configurations are separated by Bismuth atoms, forming an almost isolated bilayer structure separated by a long distance to the next bilayer structure. Furthermore, the magnetic susceptibility dataSmirnova et al. (2009) and neutron scatteringMatsuda et al. (2010) suggests two-dimensional magnetism, so it seems reasonable to model the system with a bilayer structure of Mn4+ ions. There is some experimental evidence that the interlayer coupling and the first and second neighbors intralayer couplings are the most relevant interactions and could be a strong competition between themMatsuda et al. (2010).

The outline of the paper is the following. In section II we state criteria used in determining a suitable strategy to obtain a well conditioned set of equations for the exchange couplings. This system defines a family of representative models. In section III we apply the method to define the family of magnetic models corresponding to the material Bi3Mn4O12(NO3). In section IV we present the conclusions and perspectives. Finally in the appendices A,B and C we display the single value decomposition, the detailed algorithm to optimize the set of magnetic configurations and the visual interfase for the scripts that implement the method.

II Coupling constants and Relevant Configurations

In this section, we discuss a method based on ab-initio calculation of total energy differences to estimate the coupling constants in an effective magnetic model.

We start by considering a certain atomic lattice involving magnetic atoms. Magnetism in matter is an intrinsically quantum phenomenon, requiring for its description a full quantum treatment of the electronic degrees of freedom. An exact approach of such problem is computationally unfeasible due to the huge size of its associated Hilbert space.

The usual strategy to tackle this problem is based on the relative weakness of the magnetic contribution to the energy, compared to those coming from the spatial degrees of freedom. This allows to approximate the true ground state of the system as a linear combination of those Slater’s determinant-like states

[TABLE]

that minimizes $E_{\beta}=\langle\beta|{\bf H}_{\rm full}|\beta\rangle$ for each fixed spin configuration $\{s^{\beta}_{i}\}$. By construction, these states define an orthogonal basis of the ground multiplet of the system. The evaluation of these optimizations can be performed in a relatively efficient way by means of Hartree Fock or Density Functional Theory (DFT)-like methods111In DFT, the test wave functions are not exactly Slater’s determinats, but the argument holds: GS associated to different spin configurations are orthogonal among them.. This approximation is justified if, for each $\{s_{i}\}$, the corresponding spectrum is gapped. In such a case, the true Ground State (GS) (and its low-lying excited states) can be obtained diagonalizing

[TABLE]

DFT and Hartree-Fock formalisms provide a method to evaluate (individual) diagonal elements of ${\bf H}_{0}$ in an efficient way. Notice that at this point, ${\bf H}_{0}$ is still in principle a huge matrix, in a way that even to evaluate every diagonal entry is a non-affordable task. To go further, we will suppose that ${\bf H}_{0}$ can be approximated by a simpler model, depending on a relatively small number of parameters. The family of Heisenberg’s models

[TABLE]

supplies a very versatile class of models, with a rich phase space, that do not break $SU(2)$ symmetry. Here, $B_{a}$ are sets of equivalent pairs of sites in the lattice, $J_{\alpha}$ the corresponding coupling constants ( $J_{0}$ is a global energy offset). The coupling constants $J_{\alpha}$ can then be choosen in a way that ${\bf H}_{eff}$ has diagonal entries close to the computed diagonal elements on ${\bf H}_{0}$. Notice that the condition on all the diagonal entries in both Hamiltonians defines an overconditioned set of linear equations for $J_{\alpha}$, in general it will not be possible to find $J_{\alpha}$ for a perfect match. On the other hand, as energies in DFT can only be estimated up to a finite accuracy $\Delta\varepsilon$, it makes sense to ask whether the diagonal elements in both matrices differ in less than $\Delta\varepsilon$. We define then the set of $\Delta\varepsilon-$compatible parameters

[TABLE]

in a way that any element in ${\cal C}_{\Delta\varepsilon}$ generates a ${\bf H}_{eff}$ with diagonal elements compatible with ${\bf H}_{0}$ upto the tolerance $\varepsilon$. If ${\cal C}_{\Delta\varepsilon}$ is small enough, we can expect that ${\bf H}_{eff}$ leads to the same physical predictions for any choice of $\{J_{\alpha}\}$ in ${\cal C}_{\Delta\varepsilon}$. Once we have a representative choice for $\{J_{\alpha}\}$, we can deal with ${\bf H}_{eff}$ by means of different analytical and numerical methodsAuerbach (1994), ranging from boson mapsSachdev and Bhatt (1990); Trumper et al. (1997); Flint and Coleman (2009); Lamas et al. (2013); Cabra, Lamas, and Rosales (2011a); Ceccatto, Gazza, and Trumper (1993); Mulder et al. (2010); Zhang and Lamas (2013); Zhang et al. (2016), path integralsLamas, Capponi, and Pujol (2011); Lamas et al. (2015), exact diagonalizationHonecker, Mila, and Troyer (2000), DMRGHonecker, Mila, and Troyer (2000); Elias, Arlego, and Lamas (2017), etc.

If the system involves just a very small number of magnetic atoms, the set ${\cal C}_{\Delta\varepsilon}$ can be characterized by the evaluation of the energies of all the possible spin configurations. Since ${\bf H}_{eff}$ is linear in the coupling constants $J_{\alpha}$, we can rewrite Eq. 1 as

[TABLE]

where $\vec{J}$ is a vector with components $J_{\alpha}$, $\vec{E}_{0}$ is a vector with components $E_{\beta}$, and $A$ is a matrix with coefficients

[TABLE]

and

[TABLE]

is the maximum norm or infinity normBhatia (1997). As a result, if ${\cal C}_{\Delta\varepsilon}$ is not empty, it is a convex polytope, i.e., a convex set coming from the intersection of many hemi-spaces (see Fig 1).

As the number of magnetic atoms becomes larger (let’s say, $>6$), the number of configurations grows exponentially, and therefore the evaluation of all the constraints defining the compatibility polytope ${\cal C}_{\Delta\varepsilon}$ is not feasible anymore. However, we can bound ${\cal C}_{\Delta\varepsilon}$ just looking at a smaller set of configurations: defining ${\cal S}_{n}$ as a subset of $n$ elements from the set of all possible configurations, we can define

[TABLE]

where $A^{\prime}=A^{\prime}({\cal S}_{n})$ and $E_{0}^{\prime}=E_{0}^{\prime}({\cal S}_{n})$ are built in such a way those rows correspoding to the spin configurations in ${\cal S}_{n}$ are preserved.

These sets satisfies ${\cal C}_{\Delta\varepsilon}({\cal S}_{n})\subset{\cal C}_{\Delta\varepsilon}({\cal S}_{n^{\prime}})$ if ${\cal S}_{n^{\prime}}\subset{\cal S}_{n}$ and ${\cal C}_{\Delta\varepsilon}\subset{\cal C}_{\Delta\varepsilon}({\cal S}_{n})$, so that as we increase the number of evaluated configurations, the compatibility set does not increase its size. Typically, assuming that the spectrum of ${\bf H}_{eff}$ can be accurately represented by a spin Hamiltonian with a moderately small number of couplings, most of the configurations provide no information or just redundant information. In this way, a very tight bound for ${\cal C}_{\Delta\varepsilon}$ can be obtained by just considering a small set of configurations, with a size of the order of the number of free parameters in the model. For this, it is crucial to pick the set of configurations in an optimal way: otherwise, ${\cal C}_{\Delta\varepsilon}({\cal S}_{n})$ can stay covering a much larger region than ${\cal C}_{\Delta\varepsilon}$, even for quite large $n$ (see Fig 2).

II.1 Choosing an optimal set of configurations

One problem about optimizing the size of ${\cal C}_{\Delta\varepsilon}({\cal S}_{n})$ is that, in order to be evaluated, it is necesary to know the energy of each configuration, which in general requires a lot of computing time. One strategy to avoid this issue is to make use of the inequality

[TABLE]

being $\|\vec{x}\|_{2}=\sqrt{\sum_{i=1}^{n}\vec{x}_{i}^{2}}$ the euclidean norm and $n$ the number of components of the vectorBhatia (1997). With this in mind, we define

[TABLE]

where again, $A^{\prime}=A^{\prime}({\cal S}_{n})$ and $E_{0}^{\prime}=E_{0}^{\prime}({\cal S}_{n})$ are the restrictions of $A$ and $\vec{E}_{0}$ to the rows associated to ${\cal S}_{n}$. These sets define ellipsoids centered at the minimum of the quadratic form

[TABLE]

and with main axes defined by the Singular Value Decomposition (SVD) of $A^{\prime}$:

[TABLE]

where $U^{t}U=V^{t}V={\bf 1}_{M+1}$ and $\Sigma={\rm diag}(\sigma(A^{\prime}))\in\mathbb{R}^{(M+1)\times(M+1)}$ a diagonal rectangular matrix with the singular values of the matrix $A^{\prime}$. The size of ${\cal C}^{(2)}_{\Delta\varepsilon}({\cal S}_{n})$ is then bound by

[TABLE]

where $cn$ is the condition number of $A^{\prime}$, i.e,

[TABLE]

in a way that the size of the set depends not on $\vec{E}_{0}^{\prime}$ but just on $\vec{A}^{\prime}$. This allows us to evaluate it before any DFT/ab-initio expensive simulation.

The usefulness of Def. (3) to bound the size of (1) comes from Eq. 2, that leads up to

[TABLE]

in a way that ${\cal C}_{\Delta\varepsilon}\neq\O$ if ${\cal C}^{(2)}_{\Delta\varepsilon}({\cal S}_{n})\neq\O$, and ${\cal C}^{(2)}_{\Delta\varepsilon\sqrt{n}}({\cal S}_{n})$ bounds ${\cal C}_{\Delta\varepsilon\sqrt{n}}({\cal S}_{n})$ (see Fig 3).

From the previous discussion, the strategy to obtain a good set of configurations to bound the compatibility zone is to look for ${\cal S}_{n}$ that minimizes the cost function

[TABLE]

Formally, the problem of finding the absolute minimum of $C({\cal S}_{n})$ is hard, since it typically presents many relative minima with approximately the same cost. However, what we actually need is just one these relative minima, which can be efficiently achieved by the algorithm presented in the appendix B.

II.2 Estimation of $J_{\alpha}$ and its uncertancies

Once an optimal set of configurations is determined, the corresponding magnetic energies can be estimated by means of DFT simulations. The next step is then to find the representative value $J^{(0)}_{\alpha}$ according to them. In the standard approach, $J^{(0)}_{\alpha}$ is estimated by

[TABLE]

the least square condition. This approach is valid if ${\cal C}_{\Delta\varepsilon}^{2}({\cal S}_{n})\neq\O$ for the considered tolerance, since in that case $J^{0}_{\alpha}$ belongs to ${\cal C}_{\Delta\varepsilon}({\cal S}_{n})$. This can always be achieved by choosing a large enough value for tolerance. However, in that case, the uncertancy in the estimated coupling constants could result overestimated regarding the true accuracy of the ab-initio simulation. This is a problem because the accuracy of the simulations are usually close to the scale of energy of the coupling constants which are being calculated. As a result, the estimated values for $J_{\alpha}$ are smaller in magnitud than the uncertancies, in a way that at the end of the day we are not able to state even the sign of the couplings.

To get a more realistic estimation, a Monte Carlo sampling of the region ${\cal C}_{\Delta\varepsilon}({\cal S}_{n})$ can be performed in order to get the limit values of the compatible $J_{\alpha}$. An efficient stragegy consists on explore a set of random points with a gaussian distribution around $J^{LS}_{\alpha}$, with a correlation matrix $\lambda(A^{t}A)^{-1}$ (see Fig. 4).

II.3 Improving a model

If the accuracy in the DFT energy estimation is high enough, it could happen that the proposed model becomes incompatible. In that case, a more sofisticate model can be required, for instance, by considering different two sets of coupligs that initially were considered with the same value, or by adding interactions with more distant neighbors. For this new model, the optimal set of configurations to determine the new coupling constants can be different. Although, it can be expected that there is a set of configurations, including those in the optimal set for the simpler model, that is also a relative minimum. This allows us to reuse the energies obtained in the simpler model, reducing the computational cost of estimating couplings in the new model. A way to achieve this consists on, first, optimize $C({\cal S}_{n}\cup{\cal S}^{\prime}_{n^{\prime}})$ regarding an ${\cal S}_{n^{\prime}}$ with a fixed $n^{\prime}$, being ${\cal S}_{n}$ the optimal set for the simpler model. Due to the factor $\sqrt{n}$ in $C({\cal S}_{n})$, this result can be improved by reducing the size of ${\cal S}_{n}\cup{\cal S}^{\prime}_{n^{\prime}}$, which can be accomplish by dropping one by one elements from the set.

II.4 Visualbond Spectrojotometer

The method proposed here is suitable for being automatized. In this aim, an open source project has been developedMatera and Lamas (2018) to provide to the community a software tool that performs each step in the analysis. The tool consists on a set of python libraries that helps to build magnetic models from structural data of the target magnetic compound (in cif format), and once the model is defined, to determine the optimal magnetic configurations to be evaluated, and finally, when the ab-initio simulations are performed, to estimate the corresponding coupling constants including their error bounds.

III Aplication: Magnetic model for Bi3Mn4O12(NO3)

The synthesis of the material bismuth oxynitrate, Bi3Mn4O12(NO3), obtained by Smirnova et al.Smirnova et al. (2009), has given a great impulse to the study of two dimensional effective antiferromagnetic models on the bilayer honeycomb lattice. Here, the Mn4+ ions form a honeycomb lattice, grouped in pairs and separated by a large distance by bismuth atoms. For this material the bilayer honeycomb lattice brings an appropriate geometry to build an effective Hamiltonian capable to describe its magnetic properties. The Heisenberg model on the bilayer Honeycomb lattice has been attracted a plenty of theoretical studies in the last yearsGong et al. (2013); Arlego, Lamas, and Zhang (2014); Cabra, Lamas, and Rosales (2011a, b); Fouet, Sindzingre, and Lhuillier (1975); Kandpal and van den Brink (2011); Krokhmalskii et al. (2017); Mosadeq, Shahbazi, and Jafari (2011); Oitmaa and Singh (2012). The richness of this model makes it very interesting from the theoretical point of view.

This two dimensional nature of the effective model is reinforced by magnetic susceptibility data. Also, no evidence of long-range ordering has been observed down to 0.4 KSmirnova et al. (2009); Matsuda et al. (2010); Okubo et al. (2010); Azuma et al. (2011). On the other hand, recent experimental studiesMatsuda et al. (2010) have suggested that the inter-layer coupling, as well as on-layer nearest and, to a lesser extent, next-nearest couplings, are dominant. This last characteristic of the model makes the ab-initio calculation a good tool to understand the mechanism involved by determining an appropriated effective model.

RecentlyAlaei et al. (2017) the magnetic model of Bi3Mn4O12(NO3) has been estimated by using fifty-four different spin configurations and determined the energy by ab-initio calculations with an error $\Delta E=0.5$ meV. As the DMFT calculation of the energy is hard, it is convenient to reduce the number of spin configurations needed to obtain the magnetic couplings. Using the strategy presented in the previous sections we determine a set of optimal configurations and calculate the magnetic couplings. In table 1 we present the couplings obtained inAlaei et al. (2017) using fifty-four calculations of the energy and the result using our optimal eleven configurations. All the couplings agree up to the error. Using our strategy to select the optimal configurations before the ab-initio determination of energy can save many hours of machine work and allows us to work with larger unit cells.

The optimal subset of eleven configurations that we have found with the method developed in Sec II is

[TABLE]

where the configurations are labeled according to those presented in ref. Alaei et al., 2017 and sites are labeled as in Fig. 5.

IV Conclusions

In this paper we have presented a method to find an optimal set of configurations in order to determine the couplings in a magnetic model by means of ab-initio calculations. This strategy allows to enhance the use of ab-initio calculations to establish the parameters of an effective magnetic model to be consistent with the calculation of the energies. We apply the method to calculate the family of coupling constants consistent with the ab-initio energies. We show that, while taking an optimal set of magnetic configurations, it is possible to reduce the number of ab-initio calculations to determine the model with a given error and we obtain results for the coupling constants in agreement to previous calculations. Finally we make available a free software to implement the algorithms described in the present work.

Acknowledgments

This research was partially supported by CONICET (Grant No. PIP 0691, PIP0747, PIP0720 and P-UE 22920170100066CO), UNLP (Grant No. 11/X678, 11/X680, 11/X708, 11/X788, 11/X792, 11/X806), ANPCyT (Grant No. PICT 2012-1724, 2013-2616, 2016-4083), UNNOBA (Grant No. SIB0176/2017), and Proyecto Acelerado de Cálculo 2017 de la Red Nacional de Computación de Alto Desempeño (SNCAD-MINCyT) - HPC Cluster, Rosario - Argentina.

Appendix A Singular Value Decomposition (SVD)

Let $A\in{\mathbb{C}}^{n\times m}$ be a rectangular matrix with complex entries. Then, the theorem of the Singular Value Decomposition states that there exist two rectangular matrices $U\in{\mathbb{C}}^{n\times k}$, $V\in{\mathbb{C}}^{m\times k}$ with $k\leq m,n$, and a diagonal positive matrix ${\Sigma}\in{\mathbb{C}}^{k\times k}$ such that

[TABLE]

in a way that the columns in $U$ ($V$) are orthonormals (see, for instance,Bhatia (1997)). The diagonal elements of $\Sigma$ are called the singular values of $A$, denoted by $\sigma(A)={\rm diag}(\Sigma)={\sigma_{1},\ldots\sigma_{k}}$, and the corresponding columns in $U$ ($V$) the left (right) singular vectors.

Since $AA^{\dagger}=U\Sigma^{2}U^{\dagger}$ and $A^{\dagger}A=V\Sigma^{2}V^{\dagger}$ are two semi-definite positive matrices, the left (right) singular vectors of a matrix correspond to the eigenvectors of $AA^{\dagger}$ ($A^{\dagger}A$) with non vanishing eigenvalues. If $A\in\mathbb{R}^{n\times m}$, $U$ and $V$ can be choosen as real matrices.

Among many other applications, the singular value decomposition allows to solve the linear least square problem

[TABLE]

given $A\in\mathbb{C}^{m\times n}$ ($A\in\mathbb{R}^{m\times n}$) and $\vec{b}\in\mathbb{C}^{m}$ ($\vec{b}\in\mathbb{R}^{m}$ ), and where $\|\ldots\|\equiv\|\ldots\|_{2}$ stands here for the Euclidean norm. To solve this, the Moore-Penrose pseudo inverse is introduced

[TABLE]

satisfying

[TABLE]

being $\Pi_{L}$ ($\Pi_{R}$) a projector over the subspace generated by the left (right) singular vectors of $A$. Now, choosing $\vec{x}_{0}=A^{+}\vec{b}$ we obtain

[TABLE]

In the same way, it follows that the region

[TABLE]

(if it is not empty) corresponds to an ellipsoid center at $\vec{x}_{0}=A^{+}\vec{b}$ with main $i$-esim axis parallel to the directions of the $i$-esim left singular vectors, with a size $d_{i}=\frac{\sqrt{\varepsilon^{2}-|{A}\vec{x}_{0}-\vec{b}|^{2}}}{\sigma_{i}}$.

Appendix B The algorithm

A first problem to find an optimal set of equations is related to how to produce a set of inequivalent spin configurations. This is important since for a system large enough, symmetries present on it imply that two different spin configurations result in the same condition over the couplings. To check the equivalence under symmetries is non trivial for general lattices. A more direct way to reject equivalent configurations is just by comparing the coefficients of the linear equations associated to them. The algorithm is then

Pick a random integer number between [math] and $2^{N-1}$, and assign a spin configuration according to its binary digits. 2. 2.

Build the coefficients $A_{ik}$ from this configuration. 3. 3.

If the size of the set of configurations is of the right size, return the set. 4. 4.

Pick a new random configuration in the same way than 1, and add it to the set of configurations. 5. 5.

Build $A_{ik}$. If the last row has the same coefficients that another row, remove the added configuration from the set and go to 4. Otherwise, go to 3.

The same routine can be used if we want to enlarge a previous set of independent configurations.

The second subrutine that we need for the optimization consists on a method for given a set of configurations, look for a smaller set, with a fixed size or not, in a way that the cost function is maximized. To perform this task, we implement the following algorithm:

Build the matrix $A_{ik}$. 2. 2.

Compute its SVD $A_{ik}=U\Sigma V^{t}$. 3. 3.

Assign weights to the configurations according to the (Euclidean) norm of the corresponding row in $U$. 4. 4.

Take the first $L\geq M$ configurations with larger weight. 5. 5.

If the set leads to a singular set of equations, new configurations are added one by one until the full set is saturated.

If we want to obtain a set of optimal size, we start with this method picking $L=M$ configurations in the step 4. Then, in the step 5, we continue adding configurations until the cost function starts to grow. In both cases, it is convenient to decouple the determination of the coupling constants from the determination of the non magnetic energy. This can be achieved by adding a row to ${A}_{ik}$ built by averaging the rows associated to each configuration, and afterwards, substracting it to the other rows.

To find out the set of $N>M$ optimal magnetic configurations for a given model, the algorithm starts from a random set of configuration, of size larger than the number $N^{\prime}>N$ of configurations to be determined. To generate it, we pick integer random numbers between 0 and $2^{l-1}$, and assign the spin orientations according to its binary digits. From this set, the matrix $A[{\cal S}]$ is built. If two configurations $i$ and $j$ are equivalent, in a way that $A[{\cal S}]_{ik}=A[{\cal S}]_{jk}\forall_{k}$, the configuration $j$ is dropped and a new configuration is picked. Since the last column in $A[{\cal S}]$, corresponding to the coefficient of the energy offset has all its elements equal to 1, the matrix is normalized by averaging all the rows, substracting it to each row, and adding the average row at the end. This procedure leads to an equivalent set of equations, but where the equation for the coupling constants are decoupled from the equation for the non-magnetic contribution $E_{0}$. From this initial set, the corresponding SVD

[TABLE]

and the cost function 5 are evaluated. If the cost function is $\infty$, then a new set of random inequivalent configurations is added to the set, and the procedure is repeated until a finite cost function is achieved, or a maximal size is attained. In the last case, the algorithm fails. Once an initial set of configurations with finite cost function is found, a weight to each configuration is assigned according the (euclidean) norm of the i-esim column of $U$. Using this weights, the subset of $N$ configurations of larger weights is kept. From these configurations, the normalized matrix $A[{\cal S}]$ is evaluated, and from it, the new cost function.

Appendix C Visualbond: graphical environment

In this section we are going to discuss the basics for the use of VisualBond, the graphic interface for the spectrojotometer libraryMatera and Lamas (2018). A more detailed description of functionalities and a set of examples can be found in the application documentation.

The graphic interface is composed by three tabs, corresponding to each step of the process, a menu and a status region.

The process starts by providing structural data for the compound to be simulated to the application. The current version accepts CIF files for this input, as well as partial support for struct files (Wien2k). Once the model is loaded, a reduced .cif file is shown in the edition box of the first tab. The reduced CIF includes just the lattice parameters, the set of magnetic atoms with their positions, and the set of magnetic bonds. In the first tab, at the left of the edition box, there is a dialog for automatically add new bonds in terms of the interatomic distances. By modifying the CIF file we can define different bonds among magnetic atoms.

Once the model is loaded, the other two tabs become enabled. We can then use the second tab for looking for a set of optimal configurations, loading a previously found set, or adding by hand a particular set. Once the set of configurations is loaded, the corresponding set of equations can be calculated and shown in different formats. Notice that by default, the energy associated to each configuration is established to “nan”.

Once the simulations have been performed, these “nan” values must be replaced by the energies obtained from the simulations. After doing that, in the third tab, the configurations edit box is found again. At its left, we can set the parameters for the fitting - the tolerance in the computed energies, the method for bounding the compatibility zone and its parameters - and a button to build the estimation of the coupling constants. At the right side, the corresponding equations and the estimation of the parameters, and the estimation of the relative discrepances between the values of the energy from the simulation and from the magnetic model, can be set. This tolerances can be also presented as a plot, helping to identify which configurations are best represented by the model.