Implementation of non-collinear spin-constrained DFT calculations in SIESTA with a fully relativistic Hamiltonian

An accurate and efficient general method to constrain the magnetization of individual atoms or groups of atoms within a fully relativistic non-collinear spin density functional theory formalism is presented and implemented within the SIESTA code. This approach can be applied to study a variety of complex magnetic configurations and to build effective magnetic Hamiltonians for multiscaling micromagnetic simulations. As an example, the method is applied to obtain constrained magnetic states for a Fe3 structure, and for a S = 1/2 kagome layer (vanadium oxyfluoride V7O6F18). Of paramount importance in spintronics is the control and manipulation of magnetic interactions between constituent species, characterized mainly by the pair-wise magnetic exchange tensor  ij . By constraining the atomic magnetizations of an infinite Fe linear chain, the total selfconsistent energy values are mapped to a generalized Heisenberg model, obtaining not only the diagonal terms of  ij but also the off-diagonal contributions due to the explicit presence of the spin–orbit coupling in the formalism. The diagonal values of  ij promote short ranged ferromagnetic alignment whilst the non-zero off-diagonal values can lead to the formation of the spiral states in the chain, as expected from theory.


Introduction
Magnetic materials are of much interest for a myriad of possible applications. A deep understanding of the magnetic properties is crucial, and requires a combination of experiments, theory, and modeling. In most magnetic materials, the arrangement of the atomic magnetizations adopt a collinear(CL) configuration in which a common direction is shared by all of them(usually taken as the z axis) [1,2]. Specifically, each individual magnetization can be parallel or antiparallel to each other, leading to ferromagnetic(FM) or antiferromagnetic(AFM) phases, respectively. The magnetic exchange tensor  characterizes the strength for the FM or AFM coupling between neighboring atoms.
In contrast to the CL magnetic configurations, there are a number of systems where the magnetization can take a more complex structure, and vary its direction at each point in space. This type of noncollinear(NC) structures have been observed for example in the form of helical spin density waves or spin spirals for the ground state of γ-Fe [3], and in Fe/MgO sandwiches, [4] metastable NC structures in small magnetic clusters, such as in Fe, [5] or in geometrically frustrated systems such as magnetic pyrochlore oxides [6].
The spin-orbit coupling(SOC) plays a critical role in fixing the directionality of the magnetization, and is responsible for some important magnetic properties such as the magnetic anisotropy energy(MAE), that determines the tendency of the magnetization to align along some specific axis in solids, surfaces, nanostructures, and clusters.
Several theoretical approaches have been developed during the past years to access the CL and NC states from first-principles calculations, such as the early work by von Barth and Hedin in which the unpolarized formulation of density functional theory(DFT) [7,8] was extended to include the concept of vector-spin density, treating the magnetization as a vector field [9], or the work by Kübler et al [10], where the exchangecorrelation functional was explicitly obtained for the NC case. A step forward was the generalization to deal with arbitrary magnetic configurations, i.e., configurations where the orientations of the local magnetizations are forced to predefined directions, that are unreachable unless those constraints are explicitly imposed on the selfconsistent (SC) solution. These are part of a wider class of methods in which constraints are imposed on the SC solution [11], termed constrained density functional theory (CDFT). The constrained-spin DFT was first developed by Dederichset al [12]. Over the past two and a half decades a number of interesting problems have been studied successfully using CDFT, such as the effect of Cerium impurities [12], long-range electron transfer [13], spin-dynamics [14], and magnetic exchange couplings [15][16][17] (see [11] for a review). As we will see in the following, these theoretical developments allow calculations of both s and MAEs.
Basically there are two computational strategies to determine  . One is the so-called frozen magnon approximation [18], in which the spin configuration is constrained to a spin wave with periodicity determined by a wave vector q and the energy of this spin wave is calculated through the generalized Bloch theorem for a spin-spiral configuration [19]. In the other approach, the coupling constants are calculated directly from the change of the total energy associated with constrained rotations of the spin-polarization at the sites involved [16,17,20]. The later will be the approach followed in this work. In contrast to most of the previous calculations of the magnetic exchange couplings based on constrained magnetization with a specific orientations, which relied on a scalar-relativistic formalism (i.e., the SOC is not explicitly included in the formalism), we include here a fully relativistic treatment for the Hamiltonian, including SOC and two component (spinor) wavefunctions. By overcoming this limitation, both diagonal and off-diagonal contributions to  become available. This can provide access, for example, to the parameters defining Dzyaloshinkii-Moriya interactions [21,22] (among many others).
The paper is structured as follows. In section 2.1 our DFT formalism is summarized. The mathematical structure of the CDFT as well as its relationship with the SIESTA code are detailed in section 2.2. Section 3 focuses on the optimization algorithm. Specific applications of CDFT are shown in sections 4 and 5. Finally, section 6 summarizes the main results.
2. General formalism of fully relativistic constrained DFT 2.1. Noncollinear-spin DFT Within the usual CL-spin DFT, a common arbitrary quantization magnetization axis is taken for the whole system, and the charge density ρ(r) is defined with two possible spin components r  ( ) r and r  ( ) r . The Kohn-Sham(KS) orbitals [8]  The total charge density is given by ρ(r)=r  (r)+r  (r), and the local magnetization along the quantization axis (chosen to be z) as M(r)=[r  (r)−r  (r)]u z . In the NC case the spin quantization axis varies from point to point in space and the local magnetizations M As we will see in section 2.2 it is possible to obtain the components of the local magnetization M(r) from the elements of the density matrix r mn ¢ ss .
In this work, we use the DFT implemented in the SIESTA method and code [23], which includes the description of NC spin states [24]. The SOC contribution is included selfconsistently in the KS Hamiltonian and the energy functional r [ ] E KS , following the method described in  [25]. Hence, a spinor, fully relativistic treatment of the KS Hamiltonian is used. The matrix formulation of the Hamiltonian(and density matrix) gives off-diagonal terms (i.e., non-zero couplings between the two spin components) instead of only the diagonal parts that one has for the CL case. For NC calculations, these off-diagonal terms are complex, while the diagonal ones are complex only when SOC is included in the calculation.

Noncollinear constrained spin-DFT
In general, a given SC solution corresponds to the ground state or a metastable state of the system, with a certain arrangement of the magnetization vectors. The goal of constrained-spin DFT is to force the system out from its equilibrium (ground or metastable) magnetic state in order to obtain SC configurations in which the orientation of the atomic magnetizations M i are arbitrary fixed to specific directions t i . This task is possible by the use of the Lagrange multipliers technique in which an additional term is added to the usual KS functional r [ ( )] E r KS to force the solution to satisfy the constraints that the atomic magnetizations point at the preset directions. The solution is then SC but constrained to satisfy those conditions.
The magnetization of the system at each point r, defined as x y z are the Pauli matrices, can be partitioned into local magnetization vectors by defining weight functions ( ) W r i , which can be rather arbitrary, as long as the sum of the atomic magnetizations coincides with the actual total magnetization of the system (integral of equation (7)). This allows us to define atomic magnetization vectors M i by: More generally, we can partition the magnetization on orbitals, atoms, or group of atoms. From now on, let us assume that we have N magnetic centers and their local magnetizations are M i [i=1, 2, K, N]. Following the same procedure as in the usual DFT formalism to obtain the stationary state by means of the minimization of the KS functional E KS [ρ], we can construct a new functional L that includes E KS [ρ] and additional terms that account for the directional constraints of each local magnetization: i the unitary vector that defines each local target magnetization direction and l i = l l l , the Lagrange multipliers. Apart from being simple, the constraint terms in equation (10) are linear in the density, thus avoiding multiple stationary solutions and/or energy discontinuities [26].
Stationary solutions of equation (10) for a given set of fixed t i orientations imply: and the orthonormality condition for the spinors Y ( ) r l yields a set of modified fully relativistic NC constrainedspin KS equations for the system: The value of L at the stationary solutions for each set of t { } i , obtained by the SC solution of equations (11) and (13) is the KS energy of the system in the presence of the constraints, t ({ }) E ct i . SIESTA expands the KS orbitals(equation (4)) as a a linear combination of numerical strictly localized atomic orbitals as in equation (2). The matrix elements of the fully relativistic constrained-spin DFT Hamiltonian are given by: where we can rename the components of t i by t q j = sin cos i and t q = cos z i to perform explicitly the vectorial operations bracketed in the second terms on the right hand side of equation (14): It is worth to mention that the matrix elements mn ( ) W i of the weight functions ( ) W r i only depend on the atoms or group of atoms i for which their spin magnetization has been constrained and hence they do not depend on the spin components that have to be updated when a new set Lagrange multipliers is obtained. So they can be calculated at the beginning of the calculation once for all.

Implementation
Following the work of Dederichs et al [12], instead of solving the coupled equations for the stationary points of the constrained functional(equations (11) and (13)), we cast the problem in terms of a variational one, which we can solve using minimization techniques. [11] The constraint process begins with an initial set of N atomic magnetization vectors {M i in } and a set of initial Lagrange multipliers {l i in }, and results in a set of spinors and multipliers {l i } such that equations (13) and (11) are satisfied selfconsistently. These are found by maximizing the constrained functional in equation (10) only as a function of the Lagrange multipliers {λ i }, while the KS spinors are obtained solving selfconsistently the KS equations in equation (13) for each set of values {λ i }. At the end of the process, a set of {l i out } is obtained such that the total constrained KS energy t ({ }) E ct i corresponds to a stationary value of the functional given by the equation (10) in which each one of the atomic desired magnetization orientation t i has been reached.
To perform the numerical optimization of the Lagrangian mn ¢ F ss , we follow a simple steepest descent method in which each optimization step n ct consists of a one-dimensional(1D) maximization along the opposite steepest direction L i . The line maximization is done by first bracketing the line maximum, and then by searching it using the golden section search method [27]. For these line maximizations, a tolerance in the energy of d » -10 D 1 5 eV is used. The whole optimization process finishes when each atomic magnetization lies along the target orientations. In practice, the target orientation of each atomic magnetization i is given by its azimuthal and polar angles, θ i and j i , respectively. As convergence criteria we can use the conditions: The first convergence condition accounts for the difference between the target angles and the ones obtained selfconsistently, and the second for how close the atomic magnetization r [ ] M i is with respect to the target orientation. The motivation to decide what tolerance criterium will be the best to achieve an accurate result depends on how the total SC energy varies between two different angular values and how close the torque is to zero. In principle, just one of these parameters can control the whole constraint process, but we have implemented both criteria in SIESTA, making it possible to use only one or both of them.
The additional complexity of the constraint method compared to an unconstrained calculation is that the selfconsistent solution of the constrained KS Hamiltonian must be repeated for all the values of the Lagrange multipliers during the maximization process. Fortunately once one has the first guess for the Hamiltonian or the density matrix, the number of SC steps for each maximization point is usually small, resulting in a workload which is not prohibitive with respect to unconstrained DFT calculations.

Noncollinear magnetic systems
In the following, we illustrate the use of our constrained noncollinear spin-DFT implementation to obtain nonaccessible magnetic configurations of a linear Fe 3 cluster, extending previous work by Oda et al [5] (section 4.1), and study magnetic couplings in V 7 O 6 F 18 , a spin-1/2 kagome layered structure [28], (section 4.2).

Linear Fe 3 cluster
The lowest energy structure for Fe 3 is an equilateral triangle with a collinear spin. However, a metastable linear configuration (see figure 1(a)) was shown by Oda et al [5] to have a NC ground state where the atomic magnetizations for edge atoms are slightly tilted, while for the central Fe atom it points in the direction perpendicular to the molecular axis ( figure 1(a)). This NC solution is obtained even without SOC.
Using unconstrained calculations one cannot access magnetic configurations in which the orientations of the local magnetizations deviate from this NC state, since self-consistency will align the magnetic moments according to the NC ground state solution. However, using CDFT we can fix the orientations for each M i allowing studies of magnetic excitations around this point. In the following, we illustrate this approach.
Following Oda et al [5] our calculations are based on the local density approximation under the Ceperley-Alder [29] parametrization for the exchange-correlation functional. Pseudopotentials with nonlinear core correction and a doubled zeta polarized strictly localized basis set are used to describe the valence electrons. The cluster, with an optimized bond value of 1.965 Å is placed in a cubic unit cell of 20 Å which is sufficient to avoid interactions with neighboring cluster images. Our results are in good agreement with Oda's work, and the NC ground state is obtained by means of unconstrained DFT calculation. The components of M i of the atoms at the edge are 3.17 μ B and 0.24 μ B for z and x, respectively (compared to 2.88 μ B and 0.29 μ B ). Regarding the magnetization of the central atom we obtained 1.42 μ B (1.27 μ B in Oda's work). Figure 1(b) shows the results for the total SC energy of the system, obtained by modifying the tilt angle for the edge atoms, θ, from 0°up to 10°. Any deviation from the NC ground state configuration results in higher values of the energy, as expected.

Vanadium oxyfluoride kagome slab
We have used the CDFT formalism to study possible metastable frustrated magnetic states of a V 4+ kagome slab by comparison of the total SC energy between different states. The kagome lattice is a two-dimensional network of corner-sharing triangles that is known to host exotic quantum magnetic states [30].  2(A) shows a scheme of the V atoms arrangement in the kagome structure). This spin-1/2 network is a candidate for quantum spin liquid, and the compound exhibits a high degree of magnetic frustration. Different bi-layers in the bulk are decoupled by the large organic cations, and in our model only one inorganic double-layer is considered, with compensating charges to ensure the correct number of electrons in the system.
Experimentally, it has been observed that this compound presents an alternating kagome network, with two types of triangles (one with short edges, and the other one with slightly longer edges). We have found that the orientation of the magnetization vectors of the + V 3 determines whether the + V 4 points towards the center of the large, or the small triangles, and induces a small canting out-of-plane for the kagome spins. Figure 2(B) plots the evolution of the out-of-plane vanadium + V 4 magnetization orientation θ with the number of SC steps for V 7 O 6 F 18 slab. The initial magnetization was taken in the xy plane(θ=90°) for each of the double kagome planes (left (blue) panel), and as the SC runs it freely evolves to its ground magnetic state that it is achieved after 105 SC steps. The final arrangement for the + V 4 magnetization is tilted 5°out of the xy plane, due to a magnetic interaction with the inter-layered + V 3 , aligned along the z axis. To parametrize the magnetic interactions with a Heisenberg Hamiltonian for the in-plane kagome magnetic moments, a coplanar configuration can be enforced by performing a constrained SC calculation. On the right(red) panel of the figure, the + V 4 magnetization is constrained to lay in-plane, i.e., the final magnetic configuration corresponds to a target θ angle of 90°. Starting from the ground state achieved in the unconstrained SC calculation the magnetization evolves until it reaches the target θ value, at around 140 SC additional steps. It is worth mentioning that although not depicted in the figure 2(B), several Lagrange multipliers optimization steps were performed up to achieve the desired target angle of 90°. The difference between both energies is ΔE≈5meV which means that the configuration with inplane magnetization represents a metastable magnetic state of the system, and cannot be obtained without a constraint formalism.

Application to magnetic exchange coupling
In this section we will show how to use our fully relativistic implementation to map the total SC ab initio energy to a generalized Heisenberg model. For N interacting spins, it has the following form: where S i are unitary vectors in the direction of the magnetization for each site i, and  ij are 3×3 real and symmetric matrices, usually known as the magnetic exchange tensor between spins at sites i and j. Including the full spinor formalism allows the calculation of the off-diagonal terms of  ij , which vanish for non-relativistic or scalar-relativistic calculations, since in the absence of spin-orbit interaction there is no preferred orientation of the magnetization(see the work of Udvardi et al for further [15]). In general, for an arbitrary pair of magnetic moments the six elements of the exchange tensor  ij will be different. However, symmetry considerations can frequently be employed to reduce the number of independent elements to calculate.
Prior to expanding the Heisenberg equation (16), one has to determine the range of the magnetic interactions, or in other words how many neighbours, j, of site i will be relevant. For most magnetic materials, micromagnetic calculations have shown that nearest-neighbour (NN) magnetic interactions are one order of magnitude higher than for next nearest neighbours [4]. Consequently, effective models in the literature consider no more than second neighbour interactions, and in few exceptional cases up to the third neighbours. In DFT this decision can be of paramount practical importance because the number of atoms involved in the calculation can dramatically increase the computational cost. Although SIESTA is particularly well suited to work with relatively large systems, due to the use of strictly localized numerical atomic orbitals, here we will only take into account the interaction of each atom with its NN to show with a simple example (an infinite linear Fe chain) how the fully relativistic CDFT formalism can be applied to study magnetic interactions.

Infinite linear Fe chain
As an illustrative example of how to calculate the complete set of  ij elements for a specific atomic arrangement here we obtain the magnetic exchange tensor for a free standing infinite Fe linear chain, placed along z axis.
As previously mentioned, the presence of the SOC term within the DFT Hamiltonian will bring us the possibility to calculate not only the J J J , , where the (i, j) subscripts represent pair of Fe atoms of the atomic chain. To deal with magnetic interaction between pairs of atoms under periodic boundary conditions, we need to include a sufficiently large number of atoms in the simulation cell. Here we will only consider NN interactions, with a supercell that includes four atoms. In the following we will simplify the notation by removing the subscript (i, j). Inspecting the atomic configuration in the figure 3(a), we observe that = ¹ J J J xx yy zz and = ¹ J J J xz yz xy due to the symmetry of the linear chain, reducing the total number of elements to four.
Prior to the calculation of the  elements we optimized the geometry (interatomic distances) with a conjugate gradient method under the GGA exchange-correlation functional using the PBE parametrization [31]. We used Troullier-Martins [32] norm-conserving pseudopotentials in the separable Kleinman-Bylander [33] form, and for a better description of the magnetic behavior, nonlinear core corrections were included in the exchange-correlation terms [34]. We used a double-zeta polarized basis set to describe the pseudo-valence electronic levels. The mesh cutoff to perform accurate real space integrals was fixed to 800 Ry, the electronic temperature was set to 10 meV and SCF convergence tolerance for the density matrix was set to 10 −6 . The simulation box included 20 Å of vacuum around the 1D chain to avoid any interaction with neighboring periodic replicas. The structural relaxation was performed at the scalar-relativistic level until the forces on the atoms were below 0.01 eV Å -1 , without any geometric constraint on the atomic positions nor the unit cell vectors. In agreement with previous theoretical works [35] out relaxed structure shows a dimerization with bond distances of 2.53 and 2.12 Å. In this case, there are two different magnetic exchange tensors at NN,  intra and  inter , corresponding to the intra-dimer and inter-dimer magnetic interactions, respectively.
A set of convenient magnetic configurations is chosen to map the SC total energy values q j ({ }) E , n i n i n to a Heisenberg model, equation (16). Subscript i represents a specific atomic label and n each magnetic configuration considered. The usual way to proceed is to start from a metastable magnetic state and calculate several out-ofequilibrium states that correspond to different magnetization orientations q j { } ,  With this particular arrangement we have access to ( intra + inter )/2. In the following we will refer to this tensor as  .
Considering an alternative magnetic configuration where atoms 1 and 2 are fixed and atoms 3 and 4 have a different magnetic orientation we could obtain each of the NN exchange tensors (intra or inter), separately.
For periodic systems, and to obtain the correct Heisenberg description of the total energy, one has to take into account the neighbours of the atoms on the edges of the unit cell, i.e., the Fe atoms on the edges will interact with those replicated in the adjacent cells. For instance, atom 1 has two neighbours, one in the unit cell(atom 2) and other belonging to the upper cell replica located in the positive sense of x that would correspond to atom 4. The same rule is followed for all the atoms in the unit cell. Equation (16) for the FM alignments along z, S i z =(0, 0, 1) for any i gives: On the other hand, when the magnetization of atoms 1 and 3 are kept aligned to the z axis, but those of 2 and 4 have j=0°and an arbitrary θ angle, i.e.
sin , 0, cos 2,4 ), our generalized Heisenberg model gives: J zz and J xz can be obtained by mapping this function to the set of (q n , q ( ) E n n ) pairs obtained by means of CDFT calculations. Figure 4 plots the results. A least squares fit gives » J zz 196meV and » J xz -26.5meV(solid line in the figure).  ), the angular dependence becomes: Subtracting E x,FM the following expression is obtained: which give access to J xx and J xy by fitting (j n , j ( ) E n n ). Figure 5 shows the energy dependency with j for θ=90°. The corresponding fitting gives J xx ≈192meV and J xy ≈−100meV.
Alternatively, the diagonal elements for the magnetic exchange tensor can be obtained through the difference between FM and antiferromagnetic (AFM) configurations, along x or z. The AFM configuration corresponds to antiparallel orientations among neighbouring sites, - , are in good agreement (within numerical errors), with those coming from the fitting procedure. The number of points in both graphs is the minimal necessary to obtain an accurate value of magnetic exchange tensor terms in which the fitting error is below 8%. Notice that the diagonal terms are around one order of magnitude larger than the off-diagonal ones and favor the FM order with NN. The off-diagonal terms contribute to the emergence of spiral states along the chain.

Conclusions
We have developed and implemented the fully relativistic constrained noncollinear spin-DFT formalism in the SIESTA code using the Lagrange multipliers technique. An additional term is added to the total Hamiltonian that accounts for a penalty energy contribution when an out-of-equilibrium magnetic configuration is computed self-consistently. The method has been validated first by obtaining magnetic configurations that are not accessible in conventional, unconstrained DFT calculations: (1) The linear Fe 3 cluster presents a noncollinear ground state which can be obtained by conventional selfconsistent calculations. In this configuration, the magnetization of the edges atoms are tilted by ∼4°from the molecular axis and the magnetization of the central atom points perpendicular to the linear structure [5]. Spin excitations can be explored by applying the constrained DFT to access different magnetic states where the tilt angle of the edge atoms is changed. (2) The vanadium oxyfluoride kagome slab has also been studied and we found that one of its possible magnetic state corresponds to the + V 4 atomic magnetizations slightly tilted out-of-plane. Forcing the V magnetizations to lie inplane through the CDFT formalism we obtained an increased total energy of 5meV. Figure 5. Total SC energy variation as a function of j angle for the infinite linear Fe chain. The zero of energy is set to the minimum value E 0 . Empty blue squares represent E n values when the constraint formalism is used. The blue solid line is the least squares fit of (j n , E n ) to equation (23).
The fully relativistic constrained spin-DFT formalism can further be used to parametrize effective spin Hamiltonians. By including the SOC in the total Kohn-Sham Hamiltonian, the off-diagonal terms of the magnetic exchange tensor  ij can be directly accessed through a mapping of the ab initio total energy to a generalized Heisenberg model. As an illustrative example, the full NN exchange tensor  ij of an infinite linear chain of Fe atoms has been obtained. For the NN atoms the magnetic order corresponds to a FM state and the non-zero values of the J xy and J zx can trigger spin-spiral phases.