feature(Biaxial Anisotropy): introduce more complex anisotropy term

This implementation of the anisotropy supports complex realisations of
an anisotropy energy landscape. The two anisotropy axes can be specified
independently for each site in the unit cell as well as a list of
coefficients with a set of three exponents for the angular dependencies.
For details see `core/docs/Hamiltonian.md`.

core/CMakeLists.txt

@@ -379,6 +379,7 @@ if( SPIRIT_USE_CUDA )
         ${CMAKE_CURRENT_LIST_DIR}/src/engine/interaction/Zeeman.cpp
         ${CMAKE_CURRENT_LIST_DIR}/src/engine/interaction/Anisotropy.cpp
         ${CMAKE_CURRENT_LIST_DIR}/src/engine/interaction/Cubic_Anisotropy.cpp
+        ${CMAKE_CURRENT_LIST_DIR}/src/engine/interaction/Biaxial_Anisotropy.cpp
         ${CMAKE_CURRENT_LIST_DIR}/src/engine/interaction/Exchange.cpp
         ${CMAKE_CURRENT_LIST_DIR}/src/engine/interaction/DMI.cpp
         ${CMAKE_CURRENT_LIST_DIR}/src/engine/interaction/DDI.cpp

core/docs/Hamiltonian.md

@@ -0,0 +1,79 @@
+Heisenberg Hamiltonian
+====================================================
+
+### Hamiltonian
+The Spin-Hamiltonian is defined as
+$$
+\mathcal{H}[\vec{n}] =
+    \mathcal{H}_{Z}[\vec{n}]
+    + \mathcal{H}_{A}[\vec{n}]
+    + \mathcal{H}_{XC}[\vec{n}]
+    + \mathcal{H}_{DMI}[\vec{n}]
+    + \mathcal{H}_{DDI}[\vec{n}]
+    + \mathcal{H}_{Quad}[\vec{n}]
+$$
+
+---
+
+### Zeeman
+$$ \mathcal{H}_Z[\vec{n}] = - \sum_i \mu_i \vec{B}\cdot\vec{n}_i $$
+
+---
+
+### Anisotropy
+The anisotropy term is implemented in three variants:
+- uniaxial anisotropy
+$$ \mathcal{H}_{A_1}[\vec{n}] = - \sum_i \sum_j K_j (\hat{K}_j\cdot\vec{n}_i)^2 = - \sum_i \sum_j K_j [\cos(\theta_j)]^2 $$
+- cubic anisotropy
+$$ \mathcal{H}_{A_c}[\vec{n}] = - \frac{1}{2} \sum_i \sum_j K_j ([\vec{n}_i]_x^4 + [\vec{n}_i]_y^4 + [\vec{n}_i]_z^4) $$
+- biaxial anisotropy
+$$
+\begin{alignedat}{1}
+\mathcal{H}_{AC} =& \sum_{i,j} K_j^{(n_1, n_2, n_3)} (1 - (\hat{K}_j^{(1)}\cdot\vec{n}_i)^2)^{n_1} \cdot (\hat{K}_j^{(2)}\cdot\vec{n}_i)^{n_2} \cdot ((\hat{K}_j^{(1)} \times \hat{K}_j^{(2)} ) \cdot\vec{n}_i)^{n_3} \\
+=& \sum_{i,j} K_j^{(n_1, n_2, n_3)} \cdot [\sin(\theta_j)]^{2n_1} \cdot [\cos(\varphi_j)\sin(\theta_j)]^{n_2} \cdot [\sin(\varphi_j)\sin(\theta_j)]^{n_3} \\
+\end{alignedat}
+$$
+
+The uniaxial anisotropy is equivalent to the biaxial anisotropy up to an offset to the total energy when setting
+$$
+\begin{alignedat}{1}
+\hat{K}_j^{(1)} =& \hat{K}_j \\
+K_j^{(n_1, n_2, n_3)} =& -K_j \delta_{n_1, 1}\delta_{n_2, 0}\delta_{n_3, 0}.
+\end{alignedat}
+$$
+
+---
+
+### Exchange
+- symmetric exchange interaction
+$$ \mathcal{H}_{XC}[\vec{n}] = - \sum\limits_{\braket{ij}}\, J_{ij} \vec{n}_i\cdot\vec{n}_j $$
+
+- Dzyaloshinskii–Moriya interaction (antisymmetric exchange)
+$$ \mathcal{H}_{DMI}[\vec{n}] = - \sum\limits_{\braket{ij}}\, \vec{D}_{ij} \cdot (\vec{n}_i\times\vec{n}_j) $$
+
+where it is important to note that `<ij>` denotes the unique pairs of interacting spins `i` and `j`.
+
+---
+
+### Dipole-Diplole Interaction
+$$
+\mathcal{H}_{DDI}[\vec{n}]
+    = \frac{1}{2}\frac{\mu_0}{4\pi} \sum_{\substack{i,j \\ i \neq j}} \mu_i \mu_j \frac{(\vec{n}_i \cdot \hat{r}_{ij}) (\vec{n}_j\cdot\hat{r}_{ij}) - \vec{n}_i \vec{n}_j}{{r_{ij}}^3}
+$$
+
+---
+
+### Quadruplet Interaction
+$$ \mathcal{H}_{Quad}[\vec{n}] = - \sum_{i,j,k,l} K_{ijkl} (\vec{n}_i \cdot \vec{n}_j)(\vec{n}_k \cdot \vec{n}_l) $$
+
+---
+
+<!-- ![](https://math.now.sh?from=E_%5Cmathrm%7BQuad%7D%20%3D%20-%20%5Csum%5Climits_%7Bijkl%7D%5C%2C%20K_%7Bijkl%7D%20%5Cleft%28%5Cvec%7Bn%7D_i%5Ccdot%5Cvec%7Bn%7D_j%5Cright%29%5Cleft(%5Cvec%7Bn%7D_k%5Ccdot%5Cvec%7Bn%7D_l%5Cright)) -->
+
+
+### Gaussian (test-) Hamiltonian
+
+The Hamiltonian is defined as
+$$ \mathcal{H}[\vec{n}] = \sum\limits_i \mathcal{H}_i[\vec{n}]  = \sum\limits_i a_i \exp\left( -\frac{(1 - \vec{n}\cdot\vec{c}_i)^2}{2\sigma_i^2} \right) $$
+
+<!-- ![](https://math.now.sh?from=%5Cmathcal%7BH%7D%20%3D%20%5Csum%5Climits_i%20%5Cmathcal%7BH%7D_i%20%20%3D%20%5Csum%5Climits_i%20a_i%20%5Cexp%5Cleft%28%20-%5Cfrac%7B(1%20-%20%5Cvec%7Bn%7D%5Ccdot%5Cvec%7Bc%7D_i%29%5E2%7D%7B2%5Csigma_i%5E2%7D%20%5Cright)) -->

core/include/Spirit/Hamiltonian.h

@@ -88,6 +88,11 @@ PREFIX void Hamiltonian_Set_Anisotropy(
 PREFIX void
 Hamiltonian_Set_Cubic_Anisotropy( State * state, scalar magnitude, int idx_image = -1, int idx_chain = -1 ) SUFFIX;
 
+// Set a global biaxial anisotropy [meV]
+PREFIX void Hamiltonian_Set_Biaxial_Anisotropy(
+    State * state, const scalar * magnitude, const unsigned int exponents[][3], const scalar * primary,
+    const scalar * secondary, int n_terms, int idx_image = -1, int idx_chain = -1 ) SUFFIX;
+
 // Set the exchange interaction in terms of neighbour shells [meV]
 PREFIX void Hamiltonian_Set_Exchange(
     State * state, int n_shells, const scalar * jij, int idx_image = -1, int idx_chain = -1 ) SUFFIX;
@@ -136,6 +141,20 @@ PREFIX void Hamiltonian_Get_Anisotropy(
 PREFIX void
 Hamiltonian_Get_Cubic_Anisotropy( State * state, scalar * magnitude, int idx_image = -1, int idx_chain = -1 ) SUFFIX;
 
+// Retrieves the number of atoms in the magnetic unit cell for which a biaxial anisotropy has been set.
+PREFIX int Hamiltonian_Get_Biaxial_Anisotropy_N_Atoms( State * state, int idx_image = -1, int idx_chain = -1 ) SUFFIX;
+
+// Retrieves the number of terms that contribute to the biaxial anisotropy.
+// The returned value is an ascending array of indices that mark the boundaries between the terms contributing to an
+// atom. It has length `n_atoms + 1`.
+PREFIX void Hamiltonian_Get_Biaxial_Anisotropy_N_Terms(
+    State * state, int n_atoms, int n_terms[], int idx_image = -1, int idx_chain = -1 ) SUFFIX;
+
+// Retrieves the biaxial anisotropy.
+PREFIX void Hamiltonian_Get_Biaxial_Anisotropy(
+    State * state, const int * n_terms, int * indices, scalar primary[][3], scalar secondary[][3], scalar * magnitude,
+    int exponents[][3], int n_indices, int idx_image = -1, int idx_chain = -1 ) SUFFIX;
+
 /*
 Retrieves the exchange interaction in terms of neighbour shells.
 

core/include/engine/Hamiltonian.hpp

@@ -5,6 +5,7 @@
 #include <engine/interaction/ABC.hpp>
 #include <engine/interaction/Anisotropy.hpp>
 #include <engine/interaction/Cubic_Anisotropy.hpp>
+#include <engine/interaction/Biaxial_Anisotropy.hpp>
 #include <engine/interaction/DDI.hpp>
 #include <engine/interaction/DMI.hpp>
 #include <engine/interaction/Exchange.hpp>

core/include/engine/Vectormath_Defines.hpp

@@ -89,6 +89,18 @@ struct Quadruplet
 #pragma omp declare reduction( + : Vector3 : omp_out = omp_out + omp_in ) initializer( omp_priv = Vector3::Zero() )
 #endif
 
+struct PolynomialTerm
+{
+    scalar coefficient;
+    unsigned int n1, n2, n3;
+};
+
+struct AnisotropyPolynomial
+{
+    Vector3 k1, k2, k3;
+    field<PolynomialTerm> terms;
+};
+
 struct Neighbour : Pair
 {
     // Shell index

core/include/engine/interaction/Biaxial_Anisotropy.hpp

@@ -0,0 +1,75 @@
+#pragma once
+#ifndef SPIRIT_CORE_ENGINE_INTERACTION_BIAXIAL_ANISOTROPY_HPP
+#define SPIRIT_CORE_ENGINE_INTERACTION_BIAXIAL_ANISOTROPY_HPP
+
+#include <engine/interaction/ABC.hpp>
+
+namespace Engine
+{
+
+namespace Interaction
+{
+
+class Biaxial_Anisotropy : public Interaction::Base<Biaxial_Anisotropy>
+{
+public:
+    Biaxial_Anisotropy(
+        Hamiltonian * hamiltonian, intfield indices, field<AnisotropyPolynomial> polynomials ) noexcept;
+
+    void setParameters( const intfield & pIndices, const field<AnisotropyPolynomial> & pPolynomials )
+    {
+        assert( pIndices.size() == pPolynomials.size() );
+        this->anisotropy_indices     = pIndices;
+        this->anisotropy_polynomials = pPolynomials;
+        hamiltonian->onInteractionChanged();
+    };
+    void getParameters( intfield & pIndices, field<AnisotropyPolynomial> & pPolynomials ) const
+    {
+        pIndices     = this->anisotropy_indices;
+        pPolynomials = this->anisotropy_polynomials;
+    };
+    [[nodiscard]] std::size_t getN_Atoms() const
+    {
+        return anisotropy_polynomials.size();
+    }
+
+    [[nodiscard]] std::size_t getN_Terms( const std::size_t i ) const
+    {
+        if( i >= anisotropy_polynomials.size() )
+            return 0;
+        return anisotropy_polynomials[i].terms.size();
+    }
+
+    bool is_contributing() const override;
+
+    void Energy_per_Spin( const vectorfield & spins, scalarfield & energy ) override;
+    void Hessian( const vectorfield & spins, MatrixX & hessian ) override;
+    void Sparse_Hessian( const vectorfield & spins, std::vector<triplet> & hessian ) override;
+
+    std::size_t Sparse_Hessian_Size_per_Cell() const override
+    {
+        return anisotropy_indices.size() * 9;
+    };
+
+    void Gradient( const vectorfield & spins, vectorfield & gradient ) override;
+
+    // Calculate the total energy for a single spin to be used in Monte Carlo.
+    //      Note: therefore the energy of pairs is weighted x2 and of quadruplets x4.
+    scalar Energy_Single_Spin( int ispin, const vectorfield & spins ) override;
+
+    // Interaction name as string
+    static constexpr std::string_view name          = "Generalized Anisotropy";
+    static constexpr std::optional<int> spin_order_ = std::nullopt;
+
+protected:
+    void updateFromGeometry( const Data::Geometry * geometry ) override;
+
+private:
+    intfield anisotropy_indices;
+    field<AnisotropyPolynomial> anisotropy_polynomials;
+};
+
+} // namespace Interaction
+
+} // namespace Engine
+#endif

core/include/engine/interaction/CMakeLists.txt

@@ -5,6 +5,7 @@ set(HEADER_SPIRIT_ENGINE_INTERACTION
     ${CMAKE_CURRENT_SOURCE_DIR}/Zeeman.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/Anisotropy.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/Cubic_Anisotropy.hpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/Biaxial_Anisotropy.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/Exchange.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/DMI.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/DDI.hpp

core/include/engine/interaction/Hamiltonian.hpp

@@ -37,6 +37,7 @@ class Gaussian;
 class Zeeman;
 class Anisotropy;
 class Cubic_Anisotropy;
+class Biaxial_Anisotropy;
 class Exchange;
 class DMI;
 class DDI;
@@ -82,6 +83,7 @@ class Hamiltonian
     friend class Interaction::Zeeman;
     friend class Interaction::Anisotropy;
     friend class Interaction::Cubic_Anisotropy;
+    friend class Interaction::Biaxial_Anisotropy;
     friend class Interaction::Exchange;
     friend class Interaction::DMI;
     friend class Interaction::DDI;

core/include/io/Dataparser.hpp

@@ -26,6 +26,18 @@ void Anisotropy_from_File(
     intfield & anisotropy_index, scalarfield & anisotropy_magnitude, vectorfield & anisotropy_normal,
     intfield & cubic_anisotropy_index, scalarfield & cubic_anisotropy_magnitude ) noexcept;
 
+void Biaxial_Anisotropy_Axes_from_File(
+    const std::string & anisotropy_axes_file, const std::shared_ptr<Data::Geometry> geometry, int & n_axes,
+    std::map<int, std::pair<Vector3, Vector3>> & anisotropy_axes ) noexcept;
+
+void Biaxial_Anisotropy_Terms_from_File(
+    const std::string & anisotropy_terms_file, const std::shared_ptr<Data::Geometry> geometry, int & n_terms,
+    std::map<int, std::vector<PolynomialTerm>> & anisotropy_terms ) noexcept;
+
+void Biaxial_Anisotropy_from_File(
+    const std::string & anisotropy_axes_file, const std::string & anisotropy_terms_file, const std::shared_ptr<Data::Geometry> geometry, int & n_indices,
+    intfield & anisotropy_indices, field<AnisotropyPolynomial> & anisotropy_polynomials ) noexcept;
+
 void Basis_from_File(
     const std::string & basis_file, Data::Basis_Cell_Composition & cell_composition, std::vector<Vector3> & cell_atoms,
     std::size_t & n_cell_atoms ) noexcept;