diff --git a/README.md b/README.md index 6517ea937..96c2a04e2 100644 --- a/README.md +++ b/README.md @@ -83,7 +83,8 @@ control of parameters. ### *Physics Features* -- Atomistic Spin Lattice Heisenberg Model including also DMI and dipole-dipole +- Atomistic spin lattice Heisenberg model including DMI and dipole-dipole interaction +- Micromagnetic model including DM and dipolar interactions - **Spin Dynamics simulations** obeying the [Landau-Lifschitz-Gilbert equation](https://en.wikipedia.org/wiki/Landau%E2%80%93Lifshitz%E2%80%93Gilbert_equation) - Direct **Energy minimisation** with different solvers diff --git a/core/README.md b/core/README.md index 389d6487b..8f2ab9fb3 100644 --- a/core/README.md +++ b/core/README.md @@ -4,7 +4,7 @@ Spirit This is the core library of the **Spirit** framework. It is meant to provide useful and easy API functions to enable productive work -with Atomistic Dynamics Simulations and Optimizations. +with Atomistic and Micromagnetic Dynamics Simulations and Optimizations. The current implementation is specific to atomistic spin models, but it may easily be generalised. diff --git a/core/docs/Input.md b/core/docs/Input.md index 1618088d2..45de2581f 100644 --- a/core/docs/Input.md +++ b/core/docs/Input.md @@ -6,11 +6,12 @@ The following sections will list and explain the input file keywords. 1. [General Settings and Log](#General) 2. [Geometry](#Geometry) 3. [Heisenberg Hamiltonian](#Heisenberg) -4. [Gaussian Hamiltonian](#Gaussian) -5. [Method Output](#MethodOutput) -6. [Method Parameters](#MethodParameters) -7. [Pinning](#Pinning) -8. [Disorder and Defects](#Defects) +4. [Micromagnetic Hamiltonian](#Micromagnetic) +5. [Gaussian Hamiltonian](#Gaussian) +6. [Method Output](#MethodOutput) +7. [Method Parameters](#MethodParameters) +8. [Pinning](#Pinning) +9. [Disorder and Defects](#Defects) General Settings and Log @@ -169,7 +170,7 @@ For more details, such as the notation used here, see [Phys. Rev. B **99** 22441 **General Parameters:** ```Python -### Hamiltonian Type (heisenberg_neighbours, heisenberg_pairs, gaussian) +### Hamiltonian Type (heisenberg_neighbours, heisenberg_pairs, micromagnetic, gaussian) hamiltonian heisenberg_neighbours ### Boundary conditions (in a b c) = 0(open), 1(periodical) @@ -306,6 +307,80 @@ Pairwise interactions are specified in meV per unique pair \, while quadruplets are specified in meV per unique quadruplet \. +Micromagnetic Hamiltonian +-------------------------------------------------- + +To use a micromagnetic Hamiltonian, use `micromagnetic` as input parameter after the `hamiltonian` keyword. + +Note that Spirit only supports rectilinear geometries when using this Hamiltonian. + +**Units:** + +In the micromagnetic model, Spirit uses SI units: +- `m` (meter) for distances +- `J` (Joule) for energies +- `s` (second) for time + +Therefore, +- `A [J/m]` for exchange stiffness +- `D [J/m^2]` for DMI +- `K [J/m^3]` for anisotropy +- `Ms [A/m]` for saturation magnetisation + +**General Parameters:** + +```Python +### boundary_conditions (in a b c) = 0(open), 1(periodical) +boundary_conditions 0 0 0 + +### The order of the finite difference approximation of the spatial gradient +spatial_gradient_order 2 +``` + +**Static:** + +```Python +# Saturation magnetisation [A/m] +Ms 1.3e6 +``` +while the magnetocrystalline anisotropy can be specified as an axis, +```Python +# Anisotropy [J/m^3] +anisotropy 0.3e6 +``` +or as a tensor +```Python +# Anisotropy [J/m^3] +tensor_anisotropy +0.3e6 0 0 +0 0.3e6 0 +0 0 0.3e6 +``` + +**Interactions:** + +The exchange interaction and DMI can each be set either as a constant, +```Python +# Stiffness [J/m] +exchange 10e-12 +# DMI [J/m^2] +dmi 6e-3 +``` +or as a tensor, +```Python +# Stiffness [J/m] +tensor_exchange +10e-12 0 0 + 0 10e-12 0 + 0 0 10e-12 +# DMI [J/m^2] +tensor_dmi + 0 -6e-3 6e-3 + 6e-3 0 -6e-3 +-6e-3 6e-3 0 +``` + + Gaussian Hamiltonian ---------------------------------------------------- diff --git a/core/include/Spirit/Constants.h b/core/include/Spirit/Constants.h index 5ca7ddb63..58edc1f30 100644 --- a/core/include/Spirit/Constants.h +++ b/core/include/Spirit/Constants.h @@ -28,17 +28,26 @@ PREFIX scalar Constants_k_B() SUFFIX; // Planck constant [meV*ps / rad] PREFIX scalar Constants_hbar() SUFFIX; -// Millirydberg [mRy / meV] -PREFIX scalar Constants_mRy() SUFFIX; - // Gyromagnetic ratio of electron [rad / (s*T)] PREFIX scalar Constants_gamma() SUFFIX; // Electron g-factor [unitless] PREFIX scalar Constants_g_e() SUFFIX; +// Millirydberg [mRy / meV] +PREFIX scalar Constants_mRy() SUFFIX; + +// Erg [erg / meV] +PREFIX scalar Constants_erg() SUFFIX; + +// Joule [J / meV] +PREFIX scalar Constants_Joule() SUFFIX; + // Pi [rad] PREFIX scalar Constants_Pi() SUFFIX; +// Pi/2 [rad] +PREFIX scalar Constants_Pi_2() SUFFIX; + #include "DLL_Undefine_Export.h" #endif \ No newline at end of file diff --git a/core/include/Spirit/Geometry.h b/core/include/Spirit/Geometry.h index 8ea539732..a135f107a 100644 --- a/core/include/Spirit/Geometry.h +++ b/core/include/Spirit/Geometry.h @@ -62,7 +62,7 @@ Positions are in units of the bravais vectors (scaled by the lattice constant). PREFIX void Geometry_Set_Cell_Atoms( State * state, int n_atoms, float ** atoms ) SUFFIX; /* -Set the magnetic moments of basis cell atoms. +Set the magnetic moments of basis cell atoms in units of [mu_B]. */ PREFIX void Geometry_Set_mu_s( State * state, float mu_s, int idx_image = -1, int idx_chain = -1 ) SUFFIX; @@ -130,7 +130,7 @@ Retrieve dimensionality of the system (0, 1, 2, 3). PREFIX int Geometry_Get_Dimensionality( State * state, int idx_image = -1, int idx_chain = -1 ) SUFFIX; /* -Get the magnetic moments of basis cell atoms. +Get the magnetic moments of basis cell atoms in units of [mu_B]. */ PREFIX void Geometry_Get_mu_s( State * state, float * mu_s, int idx_image = -1, int idx_chain = -1 ) SUFFIX; diff --git a/core/include/Spirit/Hamiltonian.h b/core/include/Spirit/Hamiltonian.h index b713fad1b..dcdc932b9 100644 --- a/core/include/Spirit/Hamiltonian.h +++ b/core/include/Spirit/Hamiltonian.h @@ -56,11 +56,33 @@ Dipole-Dipole method // Use a direct summation with a cutoff radius #define SPIRIT_DDI_METHOD_CUTOFF 3 +/* +Definition of Hamiltonian types +-------------------------------------------------------------------- +*/ + +typedef enum +{ + Hamiltonian_Heisenberg = 0, + Hamiltonian_Micromagnetic = 1, + Hamiltonian_Gaussian = 2 +} Hamiltonian_Type; + /* Setters -------------------------------------------------------------------- */ +/* +Set the kind of Hamiltonian to be used by all systems. +Can be (case is ignored): + +- Heisenberg +- Micromagnetic +- Gaussian +*/ +PREFIX void Hamiltonian_Set_Kind( State * state, Hamiltonian_Type type, int idx_chain = -1 ) SUFFIX; + // Set the boundary conditions along the translation directions [a, b, c] PREFIX void Hamiltonian_Set_Boundary_Conditions( State * state, const bool * periodical, int idx_image = -1, int idx_chain = -1 ) SUFFIX; @@ -162,7 +184,7 @@ Writes the 3Nx3N embedding Hessian to a file. If triplet_format is set to true the hessian is written as a list of triplets, recommended for large and sparse Hessians. */ PREFIX void Hamiltonian_Write_Hessian( - State * state, const char * filename, bool triplet_format = true, int idx_image = -1, int idx_chain = -1) SUFFIX; + State * state, const char * filename, bool triplet_format = true, int idx_image = -1, int idx_chain = -1 ) SUFFIX; #include "DLL_Undefine_Export.h" #endif \ No newline at end of file diff --git a/core/include/data/Geometry.hpp b/core/include/data/Geometry.hpp index 8e33f023d..0d1628da4 100644 --- a/core/include/data/Geometry.hpp +++ b/core/include/data/Geometry.hpp @@ -141,10 +141,17 @@ class Geometry // Dimensionality of the points int dimensionality; int dimensionality_basis; - // Center and Bounds + // Center and bounds [Å] Vector3 center, bounds_min, bounds_max; - // Unit Cell Bounds + // Unit cell bounds [Å] Vector3 cell_bounds_min, cell_bounds_max; + // Unit cell size [m] + Vector3 cell_size; + // Unit cell volume [m^3] + scalar cell_volume; + // Computes the saturation magnetisation density [A/m] + // from the atomistic basis cell composition + scalar getMs(); private: // Generate the full set of spin positions diff --git a/core/include/engine/CMakeLists.txt b/core/include/engine/CMakeLists.txt index 606eceb17..b73e76a6a 100644 --- a/core/include/engine/CMakeLists.txt +++ b/core/include/engine/CMakeLists.txt @@ -3,6 +3,7 @@ set(HEADER_SPIRIT_ENGINE ${CMAKE_CURRENT_SOURCE_DIR}/Neighbours.hpp ${CMAKE_CURRENT_SOURCE_DIR}/Hamiltonian.hpp ${CMAKE_CURRENT_SOURCE_DIR}/Hamiltonian_Heisenberg.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/Hamiltonian_Micromagnetic.hpp ${CMAKE_CURRENT_SOURCE_DIR}/Hamiltonian_Gaussian.hpp ${CMAKE_CURRENT_SOURCE_DIR}/Eigenmodes.hpp ${CMAKE_CURRENT_SOURCE_DIR}/HTST.hpp diff --git a/core/include/engine/Demagnetization_Tensor.hpp b/core/include/engine/Demagnetization_Tensor.hpp new file mode 100644 index 000000000..27cfa25b5 --- /dev/null +++ b/core/include/engine/Demagnetization_Tensor.hpp @@ -0,0 +1,299 @@ +#pragma once +#ifndef DEMAGNETIZATION_TENSOR_HPP +#define DEMAGNETIZATION_TENSOR_HPP + +#include + +#include +#include + +namespace Engine +{ +namespace Demagnetization_Tensor +{ + +namespace Exact +{ + +// Helper functions +template +scalar kappa( const scalar & x, const scalar & y, const scalar & z, const scalar & R ) +{ + // If the argument of the log gets large, it is zeroed out by the prefactor in f or g. Therefore we just set it to + // zero here, to avoid division by zero. + auto res = std::log( ( x + R ) / ( std::sqrt( y * y + z * z ) ) ); + if( std::isnan( res ) || std::isinf( res ) ) + res = 0; + return res; +} + +template +scalar delta( const scalar & x, const scalar & y, const scalar & z, const scalar & R ) +{ + const auto & Pi = Utility::Constants::Pi; + + auto arg = x * y / ( z * R ); + auto res = std::atan( arg ); + + // If the arg is infinite atan(arg) will give +- Pi/2 depending on sign + // The std::atan function should know about this, but we do not rely on it here + if( std::isinf( arg ) ) + { + if( arg < 0 ) + return -Pi / 2; + else + return Pi / 2; + } + + // If arg is nan it most likely means a division 0/0 ocurred, + // we just return 0 because the delta function is cancelled by prefactors in that case + if( std::isnan( arg ) ) + return 0; + + return res; +} + +// Helper function for Nxx. Symmetric in z and y. +template +scalar f( const scalar & x, const scalar & y, const scalar & z ) +{ + scalar R = std::sqrt( x * x + y * y + z * z ); + return ( + ( y / 2.0 ) * ( z - x ) * ( z + x ) * kappa( y, x, z, R ) + + ( z / 2.0 ) * ( y - x ) * ( y + x ) * kappa( z, x, y, R ) - x * y * z * delta( y, z, x, R ) + + 1.0 / 6.0 * ( ( x - y ) * ( x + y ) + ( x - z ) * ( x + z ) ) * R ); +} + +// Helper function for Nxy. Symmetric in x and y. +template +scalar g( const scalar & x, const scalar & y, const scalar & z ) +{ + scalar R = std::sqrt( x * x + y * y + z * z ); + return ( + ( x * y * z ) * kappa( z, x, y, R ) + ( y / 6.0 ) * ( 3 * z * z - y * y ) * kappa( x, y, z, R ) + + ( x / 6.0 ) * ( 3 * z * z - x * x ) * kappa( y, x, z, R ) - ( z * z * z / 6.0 ) * delta( x, y, z, R ) + - ( z * y * y / 2.0 ) * delta( x, z, y, R ) - ( z * x * x / 2.0 ) * delta( y, z, x, R ) - x * y * R / 3.0 ); +} + +template +scalar gamma( const scalar & e1, const scalar & e2, const scalar & e3 ) +{ + return 8.0 / std::pow( -2, std::abs( e1 ) + std::abs( e2 ) + std::abs( e3 ) ); +} + +// Exact term for the demagnetization tensor Nxx and Nxy components, see Newell 1993 +// These formulas suffer from loss of significant digits as the distance increases +// Therefore they also return an estimate of their error, which helps to guide the use of asymptotic formulas. +// See Donahue "Accurate computation of thedemagnetization tensor" + +template +scalar +Nxx( const scalar & X, const scalar & Y, const scalar & Z, const scalar & dx, const scalar & dy, const scalar & dz, + scalar & abs_error ) +{ + scalar cur_max = 0; + scalar res = 0; + for( int e1 = -1; e1 <= 1; e1++ ) + { + for( int e2 = -1; e2 <= 1; e2++ ) + { + for( int e3 = -1; e3 <= 1; e3++ ) + { + auto tmp = gamma( e1, e2, e3 ) / ( 4 * Utility::Constants::Pi * dx * dy * dz ) + * f( X + e1 * dx, Y + e2 * dy, Z + e3 * dz ); + res += tmp; + if( std::abs( tmp ) > cur_max ) + cur_max = std::abs( tmp ); + if( std::abs( res ) > cur_max ) + cur_max = std::abs( res ); + } + } + } + // The main sources of error are temporary values with large magnitudes, while the final result of the sum is small + // in magnitude. Therefore, we approximate the absolute error of the sum as the abosolute error of the largest summand. + abs_error = std::abs( cur_max ) / std::pow( 10, std::numeric_limits::digits10 ); + return res; +} + +// See Nxx comment. +template +scalar +Nxy( const scalar & X, const scalar & Y, const scalar & Z, const scalar & dx, const scalar & dy, const scalar & dz, + scalar & abs_error ) +{ + scalar cur_max = 0; + scalar res = 0; + for( int e1 = -1; e1 <= 1; e1++ ) + { + for( int e2 = -1; e2 <= 1; e2++ ) + { + for( int e3 = -1; e3 <= 1; e3++ ) + { + auto tmp = gamma( e1, e2, e3 ) / ( 4 * Utility::Constants::Pi * dx * dy * dz ) + * g( X + e1 * dx, Y + e2 * dy, Z + e3 * dz ); + res += tmp; + if( std::abs( tmp ) > cur_max ) + cur_max = std::abs( tmp ); + if( std::abs( res ) > cur_max ) + cur_max = std::abs( res ); + } + } + } + abs_error = std::abs( cur_max ) / std::pow( 10, std::numeric_limits::digits10 ); + return res; +} + +} // namespace Exact + +namespace Asymptote +{ + +// The exact formula can be rewritten as +// Nxx = 1/(4*pi*dx*dy*dz) * 2(cosh(dx * del_x)-1) * 2(cosh(dy * del_y)-1) * 2(cosh(dz * del_z)-1) * f(x,y,z) +// where del_x is to be understood as the partial derivative wrt x etc. For Nxy replace f with g. +// The cosh terms compute the finite difference terms, in the exact formula, via taylor series +// e.g: h^2 f''(x) ~ f(x+h) - 2*f(x) + f(x-h) = h^2 * 2(cosh(h * d/dx) - 1) +// We compute the asymptotes by expanding 2(cosh(h * d/dx) - 1) up to finitely many terms +// The first term we get is: +// 1/(4 * pi) * dx dy dz (del_x^2 del_y^2 del_z^2 f(x,y,z)), +// which turns out to be just the dipole approximation. +// The following asymptotes also include the next higher terms. +// See Donahue "Accurate computation of the demagnetization tensor". + +// Implements del_x^2 del_y^2 del_z^2 f(x,y,z) +template +scalar f2( const scalar & x, const scalar & y, const scalar & z, const scalar & R ) +{ + scalar R2 = R * R; + scalar R5 = R * R * R * R * R; + return ( 3.0 * x * x - R * R ) / R5; +} + +// Implements del_x^4 del_y^2 del_z^2 f(x,y,z) +template +scalar f2xx( const scalar & x, const scalar & y, const scalar & z, const scalar & R ) +{ + scalar R2 = R * R; + scalar R5 = R * R * R * R * R; + return ( -40.0 * x * x / R2 - 5.0 * ( 7.0 * x * x / R2 - 1 ) * ( -2.0 * x * x + y * y + z * z ) / R2 + 4.0 ) / R5; +} + +// Implements del_x^2 del_y^4 del_z^2 f(x,y,z) +template +scalar f2yy( const scalar & x, const scalar & y, const scalar & z, const scalar & R ) +{ + scalar R2 = R * R; + scalar R5 = R * R * R * R * R; + return ( 20.0 * y * y / R2 - 5.0 * ( 7.0 * y * y / R2 - 1.0 ) * ( -2.0 * x * x + y * y + z * z ) / R2 - 2.0 ) / R5; +} + +// Implements del_x^2 del_y^2 del_z^4 f(x,y,z) +template +scalar f2zz( const scalar & x, const scalar & y, const scalar & z, const scalar & R ) +{ + return f2yy( x, z, y, R ); // Swap y and z, since f(x,y,z) is symmetric in these +} + +// Implements del_x^2 del_y^2 del_z^2 g(x,y,z) +template +scalar g2( const scalar & x, const scalar & y, const scalar & z, const scalar & R ) +{ + scalar R2 = R * R; + scalar R5 = R * R * R * R * R; + return 3.0 * x * y / R5; +} + +// Implements del_x^4 del_y^2 del_z^2 g(x,y,z) +template +scalar g2xx( const scalar & x, const scalar & y, const scalar & z, const scalar & R ) +{ + scalar R2 = R * R; + scalar R7 = R * R * R * R * R * R * R; + return 15.0 * x * y * ( 7.0 * x * x / R2 - 3.0 ) / R7; +} + +// Implements del_x^2 del_y^4 del_z^2 g(x,y,z) +template +scalar g2yy( const scalar & x, const scalar & y, const scalar & z, const scalar & R ) +{ + return g2xx( y, x, z, R ); // Swap x and y, since g(x,y,z) is symmetric in these +} + +// Implements del_x^2 del_y^2 del_z^4 g(x,y,z) +template +scalar g2zz( const scalar & x, const scalar & y, const scalar & z, const scalar & R ) +{ + scalar R2 = R * R; + scalar R7 = R * R * R * R * R * R * R; + return 15.0 * x * y * ( 7.0 * z * z / R2 - 1.0 ) / R7; +} + +template +scalar Nxx_asym( + const scalar & X, const scalar & Y, const scalar & Z, const scalar & dx, const scalar & dy, const scalar & dz ) +{ + scalar R = std::sqrt( X * X + Y * Y + Z * Z ); + return -1.0 / ( 4.0 * Utility::Constants::Pi ) * dx * dy * dz + * ( f2( X, Y, Z, R ) + + 1.0 / 12.0 + * ( dx * dx * f2xx( X, Y, Z, R ) + dy * dy * f2yy( X, Y, Z, R ) + dz * dz * f2zz( X, Y, Z, R ) ) ); +} + +template +scalar Nxy_asym( + const scalar & X, const scalar & Y, const scalar & Z, const scalar & dx, const scalar & dy, const scalar & dz ) +{ + scalar R = std::sqrt( X * X + Y * Y + Z * Z ); + return -1.0 / ( 4.0 * Utility::Constants::Pi ) * dx * dy * dz + * ( g2( X, Y, Z, R ) + + 1.0 / 12.0 + * ( dx * dx * g2xx( X, Y, Z, R ) + dy * dy * g2yy( X, Y, Z, R ) + dz * dz * g2zz( X, Y, Z, R ) ) ); +} + +} // namespace Asymptote + +namespace Automatic +{ + +// These functions implement an automatic switching between the asymptotic expression and the exact one. +// Based on the floating point precision and the reported estimated error +template +scalar +Nxx( const scalar & X, const scalar & Y, const scalar & Z, const scalar & dx, const scalar & dy, const scalar & dz ) +{ + scalar abs_error = 0; + auto nxx_analytical = Exact::Nxx( X, Y, Z, dx, dy, dz, abs_error ); + auto nxx_asym = Asymptote::Nxx_asym( X, Y, Z, dx, dy, dz ); + + // If the asymptote is within the error due to loss of significance we use it instead of the exact formula + if( std::abs( nxx_analytical - nxx_asym ) < 20.0 * abs_error ) + return nxx_asym; + + return nxx_analytical; +} + +template +scalar +Nxy( const scalar & X, const scalar & Y, const scalar & Z, const scalar & dx, const scalar & dy, const scalar & dz ) +{ + scalar abs_error = 0; + auto nxy_analytical = Exact::Nxy( X, Y, Z, dx, dy, dz, abs_error ); + auto nxy_asym = Asymptote::Nxy_asym( X, Y, Z, dx, dy, dz ); + + // std::cout << "nxy_analytical " << nxy_analytical << "\n"; + // std::cout << "nxy_asym " << nxy_asym << "\n"; + // std::cout << "abs_error " << abs_error << "\n"; + + // If the asymptote is within the error due to loss of significance we use it instead of the exact formula + if( std::abs( nxy_analytical - nxy_asym ) < 20.0 * abs_error ) + return nxy_asym; + + return nxy_analytical; +} + +} // namespace Automatic + +} // namespace Demagnetization_Tensor +} // namespace Engine + +#endif \ No newline at end of file diff --git a/core/include/engine/Hamiltonian.hpp b/core/include/engine/Hamiltonian.hpp index 0a5519c4e..f4826e53c 100644 --- a/core/include/engine/Hamiltonian.hpp +++ b/core/include/engine/Hamiltonian.hpp @@ -5,12 +5,22 @@ #include #include +#include + #include "Spirit_Defines.h" #include namespace Engine { +enum class DDI_Method +{ + FFT = SPIRIT_DDI_METHOD_FFT, + FMM = SPIRIT_DDI_METHOD_FMM, + Cutoff = SPIRIT_DDI_METHOD_CUTOFF, + None = SPIRIT_DDI_METHOD_NONE +}; + /* * The Hamiltonian contains the interaction parameters of a System. * It also defines the functions to calculate the Effective Field and Energy. diff --git a/core/include/engine/Hamiltonian_Heisenberg.hpp b/core/include/engine/Hamiltonian_Heisenberg.hpp index 056ee5aea..5ff650bfe 100644 --- a/core/include/engine/Hamiltonian_Heisenberg.hpp +++ b/core/include/engine/Hamiltonian_Heisenberg.hpp @@ -2,9 +2,6 @@ #ifndef SPIRIT_CORE_ENGINE_HAMILTONIAN_HEISENBERG_HPP #define SPIRIT_CORE_ENGINE_HAMILTONIAN_HEISENBERG_HPP -#include -#include - #include "FFT.hpp" #include "Spirit_Defines.h" #include @@ -12,16 +9,11 @@ #include #include -namespace Engine -{ +#include +#include -enum class DDI_Method +namespace Engine { - FFT = SPIRIT_DDI_METHOD_FFT, - FMM = SPIRIT_DDI_METHOD_FMM, - Cutoff = SPIRIT_DDI_METHOD_CUTOFF, - None = SPIRIT_DDI_METHOD_NONE -}; /* The Heisenberg Hamiltonian using Pairs contains all information on the interactions between spins. diff --git a/core/include/engine/Hamiltonian_Micromagnetic.hpp b/core/include/engine/Hamiltonian_Micromagnetic.hpp new file mode 100644 index 000000000..eee550e22 --- /dev/null +++ b/core/include/engine/Hamiltonian_Micromagnetic.hpp @@ -0,0 +1,152 @@ +#pragma once +#ifndef SPIRIT_CORE_HAMILTONIAN_MICROMAGNETIC_HPP +#define SPIRIT_CORE_HAMILTONIAN_MICROMAGNETIC_HPP + +#include "Spirit_Defines.h" +#include +#include +#include +#include + +#include +#include + +namespace Engine +{ + +/* + The Micromagnetic Hamiltonian +*/ +class Hamiltonian_Micromagnetic : public Hamiltonian +{ +public: + Hamiltonian_Micromagnetic( + scalar Ms, scalar external_field_magnitude, Vector3 external_field_normal, Matrix3 anisotropy_tensor, + Matrix3 exchange_tensor, Matrix3 dmi_tensor, DDI_Method ddi_method, intfield ddi_n_periodic_images, + scalar ddi_radius, std::shared_ptr geometry, int spatial_gradient_order, + intfield boundary_conditions ); + + void Update_Interactions(); + + void Update_Energy_Contributions() override; + + void Hessian( const vectorfield & spins, MatrixX & hessian ) override; + void Gradient( const vectorfield & spins, vectorfield & gradient ) override; + void Energy_Contributions_per_Spin( + const vectorfield & spins, std::vector> & contributions ) override; + void Energy_Update( + const vectorfield & spins, std::vector> & contributions, + vectorfield & gradient ); + // 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; + + // Hamiltonian name as string + const std::string & Name() const override; + + std::shared_ptr geometry; + + // ------------ ... ------------ + int spatial_gradient_order; + + scalar Ms; + + // ------------ Single Spin Interactions ------------ + // External magnetic field across the sample + scalar external_field_magnitude; + Vector3 external_field_normal; + Matrix3 anisotropy_tensor; + + // ------------ Pair Interactions ------------ + // Exchange interaction + Matrix3 exchange_tensor; + // DMI + Matrix3 dmi_tensor; + pairfield neigh; + field spatial_gradient; + bool A_is_nondiagonal = true; + + // Dipole-dipole interaction + DDI_Method ddi_method; + intfield ddi_n_periodic_images; + scalar ddi_cutoff_radius; + pairfield ddi_pairs; + scalarfield ddi_magnitudes; + vectorfield ddi_normals; + bool ddi_pb_zero_padding = true; + +private: + // ------------ Effective Field Functions ------------ + // Calculate the Zeeman effective field of a single Spin + void Gradient_Zeeman( vectorfield & gradient ); + // Calculate the Anisotropy effective field of a single Spin + void Gradient_Anisotropy( const vectorfield & spins, vectorfield & gradient ); + // Calculate the exchange interaction effective field of a Spin Pair + void Gradient_Exchange( const vectorfield & spins, vectorfield & gradient ); + // Calculate the DMI effective field of a Spin Pair + void Gradient_DMI( const vectorfield & spins, vectorfield & gradient ); + void Spatial_Gradient( const vectorfield & spins ); + // Calculates the Dipole-Dipole contribution to the effective field of spin ispin within system s + void Gradient_DDI( const vectorfield & spins, vectorfield & gradient ); + void Gradient_DDI_Cutoff( const vectorfield & spins, vectorfield & gradient ); + void Gradient_DDI_Direct( const vectorfield & spins, vectorfield & gradient ); + void Gradient_DDI_FFT( const vectorfield & spins, vectorfield & gradient ); + + // ------------ Energy Functions ------------ + // Indices for Energy vector + int idx_zeeman, idx_anisotropy, idx_exchange, idx_dmi, idx_ddi; + void E_Update( const vectorfield & spins, scalarfield & Energy, vectorfield & gradient ); + // Calculate the Zeeman energy of a Spin System + void E_Zeeman( const vectorfield & spins, scalarfield & Energy ); + // Calculate the Anisotropy energy of a Spin System + void E_Anisotropy( const vectorfield & spins, scalarfield & Energy ); + // Calculate the exchange interaction energy of a Spin System + void E_Exchange( const vectorfield & spins, scalarfield & Energy ); + // Calculate the DMI energy of a Spin System + void E_DMI( const vectorfield & spins, scalarfield & Energy ); + // Dipolar interactions + void E_DDI( const vectorfield & spins, scalarfield & Energy ); + void E_DDI_Direct( const vectorfield & spins, scalarfield & Energy ); + void E_DDI_Cutoff( const vectorfield & spins, scalarfield & Energy ); + void E_DDI_FFT( const vectorfield & spins, scalarfield & Energy ); + + // Preparations for DDI-Convolution Algorithm + void Prepare_DDI(); + void Clean_DDI(); + + // Plans for FT / rFT + FFT::FFT_Plan fft_plan_spins; + FFT::FFT_Plan fft_plan_reverse; + + field transformed_dipole_matrices; + bool save_dipole_matrices = true; + field dipole_matrices; + + // Number of inter-sublattice contributions + int n_inter_sublattice; + // At which index to look up the inter-sublattice D-matrices + field inter_sublattice_lookup; + + // Lengths of padded system + field n_cells_padded; + // Total number of padded spins per sublattice + int sublattice_size; + + FFT::StrideContainer spin_stride; + FFT::StrideContainer dipole_stride; + + // Calculate the FT of the padded D matriess + void FFT_Demag_Tensors( FFT::FFT_Plan & fft_plan_dipole, int img_a, int img_b, int img_c ); + // Calculate the FT of the padded spins + void FFT_Spins( const vectorfield & spins ); + + // Bounds for nested for loops. Only important for the CUDA version + field it_bounds_pointwise_mult; + field it_bounds_write_gradients; + field it_bounds_write_spins; + field it_bounds_write_dipole; +}; + +} // namespace Engine + +#endif \ No newline at end of file diff --git a/core/include/io/Configparser.hpp b/core/include/io/Configparser.hpp index 2dc3f5320..4c1ae6281 100644 --- a/core/include/io/Configparser.hpp +++ b/core/include/io/Configparser.hpp @@ -10,6 +10,7 @@ #include #include #include +#include namespace IO { @@ -45,6 +46,9 @@ std::unique_ptr Hamiltonian_Heisenberg_from_Conf const std::string & config_file_name, const std::shared_ptr geometry, const std::string & hamiltonian_type ); +std::unique_ptr Hamiltonian_Micromagnetic_from_Config( + const std::string & config_file_name, const std::shared_ptr geometry ); + std::unique_ptr Hamiltonian_Gaussian_from_Config( const std::string & config_file_name, const std::shared_ptr geometry ); diff --git a/core/include/io/Dataparser.hpp b/core/include/io/Dataparser.hpp index 700b2c2c1..da009e295 100644 --- a/core/include/io/Dataparser.hpp +++ b/core/include/io/Dataparser.hpp @@ -1,4 +1,3 @@ - #pragma once #ifndef SPIRIT_CORE_IO_DATAPARSER_HPP #define SPIRIT_CORE_IO_DATAPARSER_HPP diff --git a/core/include/utility/Constants.hpp b/core/include/utility/Constants.hpp index f5d0b2d4a..e8d6b4925 100644 --- a/core/include/utility/Constants.hpp +++ b/core/include/utility/Constants.hpp @@ -34,11 +34,14 @@ double const gamma = 0.1760859644; double const g_e = 2.00231930436182; // Millirydberg [mRy/meV] -double const mRy = 1.0 / 13.605693009; +double const mRy = 1.0 / 13.605693122994; // erg [erg/meV] double const erg = 6.2415091 * 1e14; +// Joule [Joule/meV] +double const Joule = 6.2415091 * 1e+21; + // Pi [rad] double const Pi = 3.141592653589793238462643383279502884197169399375105820974; @@ -74,13 +77,50 @@ double const gamma = 0.1760859644; // Electron g-factor [unitless] double const g_e = 2.00231930436182; +} // namespace Constants_mRy + +/* +Constants by micromagnetic convention (SI units): + - Spatial scale: meters + - Energy scale: Joule + - Time scale: seconds + - Magnetic fields scale: Tesla +*/ +namespace Constants_Micromagnetic +{ +// The Bohr Magneton [Joule/T] +double const mu_B = Constants::mu_B / Constants::Joule; + +// The vacuum permeability [T^2 m^3 / Joule] +double const mu_0 = Constants::mu_0 * Constants::Joule; + +// The Boltzmann constant [J/K] +double const k_B = Constants::k_B / Constants::Joule; + +// Planck constant [J*s/rad] +double const hbar = Constants::hbar / Constants::Joule * 1e-12; + +// Gyromagnetic ratio of electron [rad/(s*T)] +double const gamma = Constants::gamma * 1e+12; + +// Electron (Landé) g-factor = gamma * hbar / mu_B [unitless] +double const g_e = Constants::g_e; + +// meV [meV/Joule] +double const meV = 1.0 / Constants::Joule; + +// Millirydberg [mRy/Joule] +double const mRy = Constants::mRy / Constants::Joule; + +// erg [erg/Joule] +double const erg = Constants::erg / Constants::Joule; + // Pi [rad] -double const Pi = 3.141592653589793238462643383279502884197169399375105820974; +double const Pi = Constants::Pi; // Pi/2 [rad] -double const Pi_2 = 1.570796326794896619231321691639751442098584699687552910487; - -} // namespace Constants_mRy +double const Pi_2 = Constants::Pi_2; +} // namespace Constants_Micromagnetic } // namespace Utility diff --git a/core/python/spirit/constants.py b/core/python/spirit/constants.py index 203710c0f..23435c808 100644 --- a/core/python/spirit/constants.py +++ b/core/python/spirit/constants.py @@ -35,12 +35,6 @@ hbar = _hbar() """Planck's constant [meV*ps / rad]""" -_mRy = _spirit.Constants_mRy -_mRy.argtypes = None -_mRy.restype = scalar -mRy = _mRy() -"""Millirydberg [mRy / meV]""" - _gamma = _spirit.Constants_gamma _gamma.argtypes = None _gamma.restype = scalar @@ -57,4 +51,22 @@ _Pi.argtypes = None _Pi.restype = scalar pi = _Pi() -"""Pi [rad]""" \ No newline at end of file +"""Pi [rad]""" + +_mRy = _spirit.Constants_mRy +_mRy.argtypes = None +_mRy.restype = scalar +mRy = _mRy() +"""MilliRydberg [mRy / meV]""" + +_erg = _spirit.Constants_erg +_erg.argtypes = None +_erg.restype = scalar +erg = _erg() +"""Erg [erg / meV]""" + +_Joule = _spirit.Constants_Joule +_Joule.argtypes = None +_Joule.restype = scalar +Joule = _Joule() +"""Joule [J / meV]""" \ No newline at end of file diff --git a/core/python/test/constants.py b/core/python/test/constants.py index 7efee4716..2b1337c66 100644 --- a/core/python/test/constants.py +++ b/core/python/test/constants.py @@ -18,27 +18,35 @@ class TestConstants(unittest.TestCase): Energy: milli-eV Time: pico-sec Magnetic field: Tesla """ - + # TODO: Find a way to easily switch between units system (eg milli-eV to milli-Ry) and test it def test_Bohr_magneton(self): self.assertEqual( scalar(constants.mu_B).value, scalar(0.057883817555).value ) - + + def test_vacuum_permeability(self): + self.assertEqual( scalar(constants.mu_0).value, scalar(2.0133545*1e-28).value ) + def test_Boltzmann_Constant(self): self.assertEqual( scalar(constants.k_B).value, scalar(0.08617330350).value ) - + def test_Planck_constant(self): self.assertEqual( scalar(constants.hbar).value, scalar(0.6582119514).value ) - - def test_millirydberg(self): - self.assertEqual( scalar(constants.mRy).value, scalar(1.0/13.605693009).value ) - + def test_gyromagnetic_ratio_of_electron(self): self.assertEqual( scalar(constants.gamma).value, scalar(0.1760859644).value ) - + def test_electron_g_factor(self): self.assertEqual( scalar(constants.g_e).value, scalar(2.00231930436182).value ) - + + def test_milliRydberg(self): + self.assertEqual( scalar(constants.mRy).value, scalar(1.0/13.605693122994).value ) + + def test_erg(self): + self.assertEqual( scalar(constants.erg).value, scalar(6.2415091*1e+14).value ) + + def test_Joule(self): + self.assertEqual( scalar(constants.Joule).value, scalar(6.2415091*1e+21).value ) ######### diff --git a/core/src/Spirit/Constants.cpp b/core/src/Spirit/Constants.cpp index 6c35d1615..9e1516d70 100644 --- a/core/src/Spirit/Constants.cpp +++ b/core/src/Spirit/Constants.cpp @@ -37,7 +37,22 @@ scalar Constants_g_e() noexcept return static_cast( Utility::Constants::g_e ); } +scalar Constants_erg() noexcept +{ + return static_cast( Utility::Constants::erg ); +} + +scalar Constants_Joule() noexcept +{ + return static_cast( Utility::Constants::Joule ); +} + scalar Constants_Pi() noexcept { return static_cast( Utility::Constants::Pi ); +} + +scalar Constants_Pi_2() noexcept +{ + return static_cast( Utility::Constants::Pi_2 ); } \ No newline at end of file diff --git a/core/src/Spirit/Geometry.cpp b/core/src/Spirit/Geometry.cpp index bdd321533..35685c2e3 100644 --- a/core/src/Spirit/Geometry.cpp +++ b/core/src/Spirit/Geometry.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -31,6 +32,10 @@ void Helper_System_Set_Geometry( Data::Spin_System & system, const Data::Geometr // Update the Heisenberg Hamiltonian if( system.hamiltonian->Name() == "Heisenberg" ) std::static_pointer_cast( system.hamiltonian )->Update_Interactions(); + + // Micromagnetic Hamiltonian + if( system.hamiltonian->Name() == "Micromagnetic" ) + std::static_pointer_cast( system.hamiltonian )->Update_Interactions(); } void Helper_State_Set_Geometry( @@ -166,8 +171,8 @@ try // The new geometry const auto & old_geometry = *state->active_image->geometry; auto new_geometry = Data::Geometry( - bravais_vectors, old_geometry.n_cells, old_geometry.cell_atoms, old_geometry.cell_composition, - old_geometry.lattice_constant, old_geometry.pinning, old_geometry.defects ); + bravais_vectors, old_geometry.n_cells, old_geometry.cell_atoms, old_geometry.cell_composition, + old_geometry.lattice_constant, old_geometry.pinning, old_geometry.defects ); // Update the State Helper_State_Set_Geometry( *state, old_geometry, new_geometry ); @@ -198,8 +203,8 @@ try // The new geometry auto & old_geometry = *state->active_image->geometry; auto new_geometry = Data::Geometry( - old_geometry.bravais_vectors, n_cells, old_geometry.cell_atoms, old_geometry.cell_composition, - old_geometry.lattice_constant, old_geometry.pinning, old_geometry.defects ); + old_geometry.bravais_vectors, n_cells, old_geometry.cell_atoms, old_geometry.cell_composition, + old_geometry.lattice_constant, old_geometry.pinning, old_geometry.defects ); // Update the State Helper_State_Set_Geometry( *state, old_geometry, new_geometry ); @@ -413,8 +418,8 @@ try // The new geometry auto & old_geometry = *state->active_image->geometry; auto new_geometry = Data::Geometry( - bravais_vectors, old_geometry.n_cells, old_geometry.cell_atoms, old_geometry.cell_composition, - old_geometry.lattice_constant, old_geometry.pinning, old_geometry.defects ); + bravais_vectors, old_geometry.n_cells, old_geometry.cell_atoms, old_geometry.cell_composition, + old_geometry.lattice_constant, old_geometry.pinning, old_geometry.defects ); // Update the State Helper_State_Set_Geometry( *state, old_geometry, new_geometry ); @@ -438,8 +443,8 @@ try // The new geometry auto & old_geometry = *state->active_image->geometry; auto new_geometry = Data::Geometry( - old_geometry.bravais_vectors, old_geometry.n_cells, old_geometry.cell_atoms, old_geometry.cell_composition, - lattice_constant, old_geometry.pinning, old_geometry.defects ); + old_geometry.bravais_vectors, old_geometry.n_cells, old_geometry.cell_atoms, old_geometry.cell_composition, + lattice_constant, old_geometry.pinning, old_geometry.defects ); // Update the State Helper_State_Set_Geometry( *state, old_geometry, new_geometry ); diff --git a/core/src/Spirit/Hamiltonian.cpp b/core/src/Spirit/Hamiltonian.cpp index a65056ff7..fcabbb2c6 100644 --- a/core/src/Spirit/Hamiltonian.cpp +++ b/core/src/Spirit/Hamiltonian.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -12,6 +13,7 @@ #include #include +#include using namespace Utility; @@ -19,6 +21,103 @@ using namespace Utility; /*---------------------------------- Set Parameters ---------------------------------------------------- */ /*------------------------------------------------------------------------------------------------------ */ +void Hamiltonian_Set_Kind( State * state, Hamiltonian_Type type, int idx_chain ) noexcept +try +{ + // TODO + if( type != Hamiltonian_Heisenberg && type != Hamiltonian_Micromagnetic && type != Hamiltonian_Gaussian ) + { + Log( Utility::Log_Level::Error, Utility::Log_Sender::API, + fmt::format( "Hamiltonian_Set_Kind: unknown type index {}", int( type ) ), -1, idx_chain ); + return; + } + + std::shared_ptr image; + std::shared_ptr chain; + + // Fetch correct indices and pointers + int idx_image = -1; + from_indices( state, idx_image, idx_chain, image, chain ); + + idx_image = 0; + std::string kind_str = ""; + if( type == Hamiltonian_Heisenberg ) + { + kind_str = "Heisenberg"; + + if( kind_str == image->hamiltonian->Name() ) + { + Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, + fmt::format( "Hamiltonian is already of {} kind. Not doing anything.", kind_str ), -1, idx_chain ); + return; + } + } + else if( type == Hamiltonian_Micromagnetic ) + { + kind_str = "Micromagnetic"; + + if( kind_str == image->hamiltonian->Name() ) + { + Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, + fmt::format( "Hamiltonian is already of {} kind. Not doing anything.", kind_str ), -1, idx_chain ); + return; + } + } + else if( type == Hamiltonian_Gaussian ) + { + kind_str = "Gaussian"; + + if( kind_str == image->hamiltonian->Name() ) + { + Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, + fmt::format( "Hamiltonian is already of {} kind. Not doing anything.", kind_str ), -1, idx_chain ); + return; + } + } + + for( auto & image : chain->images ) + { + image->Lock(); + try + { + if( type == Hamiltonian_Heisenberg ) + { + // TODO: are these the desired defaults? + image->hamiltonian = std::shared_ptr( new Engine::Hamiltonian_Heisenberg( + 0, Vector3{ 0, 0, 1 }, {}, {}, {}, {}, {}, SPIRIT_CHIRALITY_NEEL, Engine::DDI_Method::None, + { 0, 0, 0 }, false, 0, {}, {}, image->geometry, image->hamiltonian->boundary_conditions ) ); + } + else if( type == Hamiltonian_Micromagnetic ) + { + // TODO: are these the desired defaults? + image->hamiltonian = std::shared_ptr( new Engine::Hamiltonian_Micromagnetic( + 0, 0, Vector3{ 0, 0, 1 }, Matrix3::Zero(), Matrix3::Zero(), Matrix3::Zero(), + Engine::DDI_Method::None, { 0, 0, 0 }, 0, image->geometry, 2, + image->hamiltonian->boundary_conditions ) ); + } + else if( type == Hamiltonian_Gaussian ) + { + // TODO: are these the desired defaults? + image->hamiltonian + = std::shared_ptr( new Engine::Hamiltonian_Gaussian( {}, {}, {} ) ); + } + } + catch( ... ) + { + spirit_handle_exception_api( idx_image, idx_chain ); + } + image->Unlock(); + ++idx_image; + } + + Log( Utility::Log_Level::All, Utility::Log_Sender::API, fmt::format( "Set Hamiltonian kind to {}", kind_str ), -1, + idx_chain ); +} +catch( ... ) +{ + spirit_handle_exception_api( -1, idx_chain ); +} + void Hamiltonian_Set_Boundary_Conditions( State * state, const bool * periodical, int idx_image, int idx_chain ) noexcept try @@ -61,6 +160,13 @@ try // Fetch correct indices and pointers from_indices( state, idx_image, idx_chain, image, chain ); + if( image->hamiltonian->Name() != "Heisenberg" && image->hamiltonian->Name() != "Micromagnetic" ) + { + Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, + "External field cannot be set on " + image->hamiltonian->Name(), idx_image, idx_chain ); + return; + } + // Lock mutex because simulations may be running image->Lock(); @@ -81,15 +187,22 @@ try // Update Energies ham->Update_Energy_Contributions(); + } + else if( image->hamiltonian->Name() == "Micromagnetic" ) + { + auto ham = (Engine::Hamiltonian_Micromagnetic *)image->hamiltonian.get(); - Log( Utility::Log_Level::Info, Utility::Log_Sender::API, - fmt::format( - "Set external field to {}, direction ({}, {}, {})", magnitude, normal[0], normal[1], normal[2] ), - idx_image, idx_chain ); + // Normals + Vector3 new_normal{ normal[0], normal[1], normal[2] }; + new_normal.normalize(); + + // Into the Hamiltonian + ham->external_field_magnitude = magnitude; + ham->external_field_normal = new_normal; + + // Update Energies + ham->Update_Energy_Contributions(); } - else - Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, - "External field cannot be set on " + image->hamiltonian->Name(), idx_image, idx_chain ); } catch( ... ) { @@ -98,6 +211,11 @@ try // Unlock mutex image->Unlock(); + + Log( Utility::Log_Level::Info, Utility::Log_Sender::API, + fmt::format( + "Set external field to {} [T], direction ({}, {}, {})", magnitude, normal[0], normal[1], normal[2] ), + idx_image, idx_chain ); } catch( ... ) { @@ -114,12 +232,22 @@ try // Fetch correct indices and pointers from_indices( state, idx_image, idx_chain, image, chain ); + if( image->hamiltonian->Name() != "Heisenberg" && image->hamiltonian->Name() != "Micromagnetic" ) + { + Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, + "Anisotropy cannot be set on " + image->hamiltonian->Name(), idx_image, idx_chain ); + return; + } + image->Lock(); try { + std::string units = "N/A"; + if( image->hamiltonian->Name() == "Heisenberg" ) { + units = "meV"; auto * ham = dynamic_cast( image->hamiltonian.get() ); int nos = image->nos; int n_cell_atoms = image->geometry->n_cell_atoms; @@ -144,15 +272,28 @@ try // Update Energies ham->Update_Energy_Contributions(); + } + else if( image->hamiltonian->Name() == "Micromagnetic" ) + { + units = "J/m^3"; + auto ham = (Engine::Hamiltonian_Micromagnetic *)image->hamiltonian.get(); - Log( Utility::Log_Level::Info, Utility::Log_Sender::API, - fmt::format( - "Set anisotropy to {}, direction ({}, {}, {})", magnitude, normal[0], normal[1], normal[2] ), - idx_image, idx_chain ); + Vector3 Kn{ normal[0], normal[1], normal[2] }; + Kn.normalize(); + + ham->anisotropy_tensor << Kn[0] * Kn[0], Kn[0] * Kn[1], Kn[0] * Kn[2], Kn[1] * Kn[0], Kn[1] * Kn[1], + Kn[1] * Kn[2], Kn[2] * Kn[0], Kn[2] * Kn[1], Kn[2] * Kn[2]; + ham->anisotropy_tensor *= magnitude; + + // Update Energies + ham->Update_Energy_Contributions(); } - else - Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, - "Anisotropy cannot be set on " + image->hamiltonian->Name(), idx_image, idx_chain ); + + Log( Utility::Log_Level::Info, Utility::Log_Sender::API, + fmt::format( + "Set anisotropy to {} [{}], direction ({}, {}, {})", magnitude, units, normal[0], normal[1], + normal[2] ), + idx_image, idx_chain ); } catch( ... ) { @@ -175,6 +316,13 @@ try // Fetch correct indices and pointers from_indices( state, idx_image, idx_chain, image, chain ); + if( image->hamiltonian->Name() != "Heisenberg" && image->hamiltonian->Name() != "Micromagnetic" ) + { + Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, + "Exchange interaction cannot be set on " + image->hamiltonian->Name(), idx_image, idx_chain ); + return; + } + image->Lock(); try @@ -190,12 +338,34 @@ try std::string message = fmt::format( "Set exchange to {} shells", n_shells ); if( n_shells > 0 ) - message += fmt::format( " Jij[0] = {}", jij[0] ); + message += fmt::format( " Jij[0] = {} [meV/bond]", jij[0] ); Log( Utility::Log_Level::Info, Utility::Log_Sender::API, message, idx_image, idx_chain ); } - else - Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, - "Exchange cannot be set on " + image->hamiltonian->Name(), idx_image, idx_chain ); + else if( image->hamiltonian->Name() == "Micromagnetic" ) + { + if( n_shells > 1 ) + { + Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, + "Cannot set more than one shell of Exchange interaction on micromagnetic Hamiltonian.", idx_image, + idx_chain ); + image->Unlock(); + return; + } + + auto ham = (Engine::Hamiltonian_Micromagnetic *)image->hamiltonian.get(); + auto Jij = jij[0]; + ham->exchange_tensor << Jij, 0, 0, 0, Jij, 0, 0, 0, Jij; + ham->Update_Interactions(); + + Log( Utility::Log_Level::Info, Utility::Log_Sender::API, fmt::format( "Set tensor_exchange to:" ), + idx_image, idx_chain ); + Log( Utility::Log_Level::Info, Utility::Log_Sender::API, fmt::format( "{}", ham->exchange_tensor.row( 0 ) ), + idx_image, idx_chain ); + Log( Utility::Log_Level::Info, Utility::Log_Sender::API, fmt::format( "{}", ham->exchange_tensor.row( 1 ) ), + idx_image, idx_chain ); + Log( Utility::Log_Level::Info, Utility::Log_Sender::API, fmt::format( "{}", ham->exchange_tensor.row( 2 ) ), + idx_image, idx_chain ); + } } catch( ... ) { @@ -218,7 +388,13 @@ try // Fetch correct indices and pointers from_indices( state, idx_image, idx_chain, image, chain ); - image->Lock(); + + if( image->hamiltonian->Name() != "Heisenberg" && image->hamiltonian->Name() != "Micromagnetic" ) + { + Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, + "Dzyaloshinskii-Moriya interaction cannot be set on " + image->hamiltonian->Name(), idx_image, idx_chain ); + return; + } if( chirality != SPIRIT_CHIRALITY_BLOCH && chirality != SPIRIT_CHIRALITY_NEEL && chirality != SPIRIT_CHIRALITY_BLOCH_INVERSE && chirality != SPIRIT_CHIRALITY_NEEL_INVERSE ) @@ -228,6 +404,8 @@ try return; } + image->Lock(); + try { if( image->hamiltonian->Name() == "Heisenberg" ) @@ -243,12 +421,51 @@ try std::string message = fmt::format( "Set dmi to {} shells", n_shells ); if( n_shells > 0 ) - message += fmt::format( " Dij[0] = {}", dij[0] ); + message += fmt::format( " Dij[0] = {} [meV/bond]", dij[0] ); Log( Utility::Log_Level::Info, Utility::Log_Sender::API, message, idx_image, idx_chain ); } - else - Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, - "DMI cannot be set on " + image->hamiltonian->Name(), idx_image, idx_chain ); + else if( image->hamiltonian->Name() == "Micromagnetic" ) + { + if( n_shells > 1 ) + { + Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, + "Cannot set more than one shell of DM interaction on micromagnetic Hamiltonian.", idx_image, + idx_chain ); + image->Unlock(); + return; + } + + auto ham = (Engine::Hamiltonian_Micromagnetic *)image->hamiltonian.get(); + auto Dij = dij[0]; + + if( chirality == SPIRIT_CHIRALITY_BLOCH ) + { + ham->dmi_tensor << Dij, 0, 0, 0, Dij, 0, 0, 0, Dij; + } + else if( chirality == SPIRIT_CHIRALITY_BLOCH_INVERSE ) + { + ham->dmi_tensor << -Dij, 0, 0, 0, -Dij, 0, 0, 0, -Dij; + } + else if( chirality == SPIRIT_CHIRALITY_NEEL ) + { + ham->dmi_tensor << 0, Dij, 0, -Dij, 0, 0, 0, 0, 0; + } + else + { + ham->dmi_tensor << 0, -Dij, 0, Dij, 0, 0, 0, 0, 0; + } + + ham->Update_Interactions(); + + Log( Utility::Log_Level::Info, Utility::Log_Sender::API, fmt::format( "Set tensor_dmi to:" ), idx_image, + idx_chain ); + Log( Utility::Log_Level::Info, Utility::Log_Sender::API, fmt::format( "{}", ham->dmi_tensor.row( 0 ) ), + idx_image, idx_chain ); + Log( Utility::Log_Level::Info, Utility::Log_Sender::API, fmt::format( "{}", ham->dmi_tensor.row( 1 ) ), + idx_image, idx_chain ); + Log( Utility::Log_Level::Info, Utility::Log_Sender::API, fmt::format( "{}", ham->dmi_tensor.row( 2 ) ), + idx_image, idx_chain ); + } } catch( ... ) { @@ -272,6 +489,14 @@ try // Fetch correct indices and pointers from_indices( state, idx_image, idx_chain, image, chain ); + + if( image->hamiltonian->Name() != "Heisenberg" && image->hamiltonian->Name() != "Micromagnetic" ) + { + Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, + "Dipolar interactions cannot be set on " + image->hamiltonian->Name(), idx_image, idx_chain ); + return; + } + image->Lock(); try @@ -295,9 +520,24 @@ try pb_zero_padding ), idx_image, idx_chain ); } - else - Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, - "DDI cannot be set on " + image->hamiltonian->Name(), idx_image, idx_chain ); + else if( image->hamiltonian->Name() == "Micromagnetic" ) + { + auto ham = (Engine::Hamiltonian_Micromagnetic *)image->hamiltonian.get(); + ham->ddi_method = Engine::DDI_Method( ddi_method ); + ham->ddi_n_periodic_images[0] = n_periodic_images[0]; + ham->ddi_n_periodic_images[1] = n_periodic_images[1]; + ham->ddi_n_periodic_images[2] = n_periodic_images[2]; + ham->ddi_cutoff_radius = cutoff_radius; + ham->ddi_pb_zero_padding = pb_zero_padding; + ham->Update_Interactions(); + + Log( Utility::Log_Level::Info, Utility::Log_Sender::API, + fmt::format( + "Set ddi to method {}, periodic images {} {} {}, cutoff radius {} and pb_zero_padding {}", + ddi_method, n_periodic_images[0], n_periodic_images[1], n_periodic_images[2], cutoff_radius, + pb_zero_padding ), + idx_image, idx_chain ); + } } catch( ... ) { @@ -363,23 +603,26 @@ try { auto * ham = dynamic_cast( image->hamiltonian.get() ); - if( ham->external_field_magnitude > 0 ) - { - // Magnitude - *magnitude = (float)( ham->external_field_magnitude / Constants::mu_B ); + *magnitude = (float)( ham->external_field_magnitude / Constants::mu_B ); + normal[0] = (float)ham->external_field_normal[0]; + normal[1] = (float)ham->external_field_normal[1]; + normal[2] = (float)ham->external_field_normal[2]; + } + else if( image->hamiltonian->Name() == "Micromagnetic" ) + { + auto ham = (Engine::Hamiltonian_Micromagnetic *)image->hamiltonian.get(); - // Normal - normal[0] = (float)ham->external_field_normal[0]; - normal[1] = (float)ham->external_field_normal[1]; - normal[2] = (float)ham->external_field_normal[2]; - } - else - { - *magnitude = 0; - normal[0] = 0; - normal[1] = 0; - normal[2] = 1; - } + *magnitude = (float)ham->external_field_magnitude; + normal[0] = (float)ham->external_field_normal[0]; + normal[1] = (float)ham->external_field_normal[1]; + normal[2] = (float)ham->external_field_normal[2]; + } + else + { + *magnitude = 0; + normal[0] = 0; + normal[1] = 0; + normal[2] = 1; } } catch( ... ) @@ -403,13 +646,10 @@ try if( !ham->anisotropy_indices.empty() ) { - // Magnitude *magnitude = (float)ham->anisotropy_magnitudes[0]; - - // Normal - normal[0] = (float)ham->anisotropy_normals[0][0]; - normal[1] = (float)ham->anisotropy_normals[0][1]; - normal[2] = (float)ham->anisotropy_normals[0][2]; + normal[0] = (float)ham->anisotropy_normals[0][0]; + normal[1] = (float)ham->anisotropy_normals[0][1]; + normal[2] = (float)ham->anisotropy_normals[0][2]; } else { @@ -419,6 +659,17 @@ try normal[2] = 1; } } + else if( image->hamiltonian->Name() == "Micromagnetic" ) + { + auto ham = (Engine::Hamiltonian_Micromagnetic *)image->hamiltonian.get(); + } + else + { + *magnitude = 0; + normal[0] = 0; + normal[1] = 0; + normal[2] = 1; + } } catch( ... ) { @@ -447,6 +698,17 @@ try jij[i] = (float)ham->exchange_shell_magnitudes[i]; } } + else if( image->hamiltonian->Name() == "Micromagnetic" ) + { + auto ham = (Engine::Hamiltonian_Micromagnetic *)image->hamiltonian.get(); + + *n_shells = 1; + jij[0] = (float)ham->exchange_tensor( 0, 0 ); + } + else + { + *n_shells = 0; + } } catch( ... ) { @@ -529,6 +791,17 @@ try dij[i] = (float)ham->dmi_shell_magnitudes[i]; } } + else if( image->hamiltonian->Name() == "Micromagnetic" ) + { + auto ham = (Engine::Hamiltonian_Micromagnetic *)image->hamiltonian.get(); + + *n_shells = 1; + dij[0] = (float)ham->dmi_tensor( 0, 0 ); + } + else + { + *n_shells = 0; + } } catch( ... ) { @@ -580,6 +853,23 @@ try *cutoff_radius = (float)ham->ddi_cutoff_radius; *pb_zero_padding = ham->ddi_pb_zero_padding; } + else if( image->hamiltonian->Name() == "Micromagnetic" ) + { + auto ham = (Engine::Hamiltonian_Micromagnetic *)image->hamiltonian.get(); + + *ddi_method = (int)ham->ddi_method; + n_periodic_images[0] = (int)ham->ddi_n_periodic_images[0]; + n_periodic_images[1] = (int)ham->ddi_n_periodic_images[1]; + n_periodic_images[2] = (int)ham->ddi_n_periodic_images[2]; + *cutoff_radius = (float)ham->ddi_cutoff_radius; + } + else + { + Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, + image->hamiltonian->Name() + + " Hamiltonian: fetching dipole-dipole interaction parameters is not available...", + idx_image, idx_chain ); + } } catch( ... ) { @@ -607,9 +897,9 @@ void saveTriplets( std::string fname, const SpMatrixX & matrix ) std::ofstream file( fname ); if( file && file.is_open() ) { - for (int k=0; k < matrix.outerSize(); ++k) + for( int k = 0; k < matrix.outerSize(); ++k ) { - for (SpMatrixX::InnerIterator it(matrix,k); it; ++it) + for( SpMatrixX::InnerIterator it( matrix, k ); it; ++it ) { file << it.row() << "\t"; // row index file << it.col() << "\t"; // col index (here it is equal to k) @@ -623,9 +913,8 @@ void saveTriplets( std::string fname, const SpMatrixX & matrix ) } } - void Hamiltonian_Write_Hessian( - State * state, const char * filename, bool triplet_format, int idx_image, int idx_chain) noexcept + State * state, const char * filename, bool triplet_format, int idx_image, int idx_chain ) noexcept { std::shared_ptr image; std::shared_ptr chain; @@ -635,11 +924,11 @@ void Hamiltonian_Write_Hessian( // Compute hessian auto nos = image->geometry->nos; - SpMatrixX hessian(3*nos, 3*nos); - image->hamiltonian->Sparse_Hessian(*image->spins, hessian); + SpMatrixX hessian( 3 * nos, 3 * nos ); + image->hamiltonian->Sparse_Hessian( *image->spins, hessian ); - if (triplet_format) - saveTriplets(std::string(filename), hessian); + if( triplet_format ) + saveTriplets( std::string( filename ), hessian ); else - saveMatrix(std::string(filename), hessian); + saveMatrix( std::string( filename ), hessian ); } \ No newline at end of file diff --git a/core/src/data/Geometry.cpp b/core/src/data/Geometry.cpp index 5b9cf2834..4d9f5d715 100644 --- a/core/src/data/Geometry.cpp +++ b/core/src/data/Geometry.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include @@ -87,6 +88,12 @@ Geometry::Geometry( // For updates of triangulation and tetrahedra this->last_update_n_cell_step = -1; this->last_update_n_cells = intfield( 3, -1 ); + + for( int i = 0; i < 3; i++ ) + this->cell_size[i] = 1e-10 * lattice_constant * bravais_vectors[i].norm(); + + this->cell_volume + = 1e-30 * pow( lattice_constant, 3 ) * bravais_vectors[0].dot( bravais_vectors[1].cross( bravais_vectors[2] ) ); } void Geometry::generatePositions() @@ -161,6 +168,15 @@ void Geometry::generatePositions() } } +scalar Geometry::getMs() +{ + // Saturation Magnetisation Density in A/m + scalar Ms = 0; + for( auto & mu_s : cell_composition.mu_s ) + Ms += mu_s * Utility::Constants_Micromagnetic::mu_B / cell_volume; + return Ms; +} + std::vector compute_delaunay_triangulation_3D( const std::vector & points ) try { diff --git a/core/src/data/Spin_System.cpp b/core/src/data/Spin_System.cpp index 03b728530..31e853f2d 100644 --- a/core/src/data/Spin_System.cpp +++ b/core/src/data/Spin_System.cpp @@ -68,6 +68,11 @@ try this->hamiltonian = std::make_shared( static_cast( *other.hamiltonian ) ); } + else if( other.hamiltonian->Name() == "Micromagnetic" ) + { + this->hamiltonian = std::make_shared( + static_cast( *other.hamiltonian ) ); + } else if( other.hamiltonian->Name() == "Gaussian" ) { this->hamiltonian = std::make_shared( @@ -113,6 +118,11 @@ try this->hamiltonian = std::make_shared( *(Engine::Hamiltonian_Heisenberg *)( other.hamiltonian.get() ) ); } + else if( other.hamiltonian->Name() == "Micromagnetic" ) + { + this->hamiltonian = std::make_shared( + *(Engine::Hamiltonian_Micromagnetic *)( other.hamiltonian.get() ) ); + } else if( other.hamiltonian->Name() == "Gaussian" ) { this->hamiltonian = std::make_shared( diff --git a/core/src/engine/CMakeLists.txt b/core/src/engine/CMakeLists.txt index fc23e0533..5ec135283 100644 --- a/core/src/engine/CMakeLists.txt +++ b/core/src/engine/CMakeLists.txt @@ -4,6 +4,8 @@ set(SOURCE_SPIRIT_ENGINE ${CMAKE_CURRENT_SOURCE_DIR}/Hamiltonian.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Hamiltonian_Heisenberg.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Hamiltonian_Heisenberg.cu + ${CMAKE_CURRENT_SOURCE_DIR}/Hamiltonian_Micromagnetic.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/Hamiltonian_Micromagnetic.cu ${CMAKE_CURRENT_SOURCE_DIR}/Hamiltonian_Gaussian.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Eigenmodes.cpp ${CMAKE_CURRENT_SOURCE_DIR}/HTST.cpp diff --git a/core/src/engine/Hamiltonian_Micromagnetic.cpp b/core/src/engine/Hamiltonian_Micromagnetic.cpp new file mode 100644 index 000000000..9ce57091a --- /dev/null +++ b/core/src/engine/Hamiltonian_Micromagnetic.cpp @@ -0,0 +1,922 @@ +#ifndef SPIRIT_USE_CUDA + +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include + +#include + +using namespace Data; +using namespace Utility; +namespace C = Utility::Constants; +using Engine::Vectormath::check_atom_type; +using Engine::Vectormath::idx_from_pair; +using Engine::Vectormath::idx_from_tupel; + +namespace Engine +{ + +Hamiltonian_Micromagnetic::Hamiltonian_Micromagnetic( + scalar Ms, scalar external_field_magnitude, Vector3 external_field_normal, Matrix3 anisotropy_tensor, + Matrix3 exchange_tensor, Matrix3 dmi_tensor, DDI_Method ddi_method, intfield ddi_n_periodic_images, + scalar ddi_radius, std::shared_ptr geometry, int spatial_gradient_order, + intfield boundary_conditions ) + : Ms( Ms ), + Hamiltonian( boundary_conditions ), + spatial_gradient_order( spatial_gradient_order ), + geometry( geometry ), + external_field_magnitude( external_field_magnitude ), + external_field_normal( external_field_normal ), + anisotropy_tensor( anisotropy_tensor ), + exchange_tensor( exchange_tensor ), + dmi_tensor( dmi_tensor ), + ddi_method( ddi_method ), + ddi_n_periodic_images( ddi_n_periodic_images ), + ddi_cutoff_radius( ddi_radius ), + fft_plan_reverse( FFT::FFT_Plan() ), + fft_plan_spins( FFT::FFT_Plan() ) +{ + // Generate interaction pairs, constants etc. + this->Update_Interactions(); +} + +void Hamiltonian_Micromagnetic::Update_Interactions() +{ +#if defined( SPIRIT_USE_OPENMP ) + // When parallelising (cuda or openmp), we need all neighbours per spin + const bool use_redundant_neighbours = true; +#else + // When running on a single thread, we can ignore redundant neighbours + const bool use_redundant_neighbours = false; +#endif + + // TODO: make sure that the geometry can be treated with this model: + // - rectilinear, only one "atom" per cell + // if( geometry->n_cell_atoms != 1 ) + // Log(...) + + // TODO: generate neighbour information for pairwise interactions + + // TODO: prepare dipolar interactions + + // Update, which terms still contribute + + neigh = pairfield( 0 ); + neigh.push_back( { 0, 0, { 1, 0, 0 } } ); + neigh.push_back( { 0, 0, { -1, 0, 0 } } ); + neigh.push_back( { 0, 0, { 0, 1, 0 } } ); + neigh.push_back( { 0, 0, { 0, -1, 0 } } ); + neigh.push_back( { 0, 0, { 0, 0, 1 } } ); + neigh.push_back( { 0, 0, { 0, 0, -1 } } ); + + this->spatial_gradient = field( geometry->nos, Matrix3::Zero() ); + this->Prepare_DDI(); + this->Update_Energy_Contributions(); +} + +void Hamiltonian_Micromagnetic::Update_Energy_Contributions() +{ + this->energy_contributions_per_spin = std::vector>( 0 ); + + // External field + if( std::abs( this->external_field_magnitude ) > 0 ) + { + this->energy_contributions_per_spin.push_back( { "Zeeman", scalarfield( 0 ) } ); + this->idx_zeeman = this->energy_contributions_per_spin.size() - 1; + } + else + this->idx_zeeman = -1; + + if( anisotropy_tensor.norm() > 0.0 ) + { + this->energy_contributions_per_spin.push_back( { "Anisotropy", scalarfield( 0 ) } ); + this->idx_anisotropy = this->energy_contributions_per_spin.size() - 1; + } + else + this->idx_anisotropy = -1; + + if( exchange_tensor.norm() > 0.0 ) + { + this->energy_contributions_per_spin.push_back( { "Exchange", scalarfield( 0 ) } ); + this->idx_exchange = this->energy_contributions_per_spin.size() - 1; + } + else + this->idx_exchange = -1; + + if( dmi_tensor.norm() > 0.0 ) + { + this->energy_contributions_per_spin.push_back( { "DMI", scalarfield( 0 ) } ); + this->idx_dmi = this->energy_contributions_per_spin.size() - 1; + } + else + this->idx_dmi = -1; + + if( this->ddi_method != DDI_Method::None ) + { + this->energy_contributions_per_spin.push_back( { "DDI", scalarfield( 0 ) } ); + this->idx_ddi = this->energy_contributions_per_spin.size() - 1; + } + else + this->idx_ddi = -1; + + // printf("idx_zeeman %i\n", idx_zeeman); + // printf("idx_exchange %i\n", idx_exchange); + // printf("idx_dmi %i\n", idx_dmi); + // printf("idx_anisotropy %i\n", idx_anisotropy); + // printf("idx_ddi %i\n", idx_ddi); + // std::cout << exchange_tensor << "\n"; + // std::cout << dmi_tensor << "\n ===== \n"; +} + +void Hamiltonian_Micromagnetic::Energy_Contributions_per_Spin( + const vectorfield & spins, std::vector> & contributions ) +{ + if( contributions.size() != this->energy_contributions_per_spin.size() ) + { + contributions = this->energy_contributions_per_spin; + } + + int nos = spins.size(); + for( auto & contrib : contributions ) + { + // Allocate if not already allocated + if( contrib.second.size() != nos ) + contrib.second = scalarfield( nos, 0 ); + // Otherwise set to zero + else + Vectormath::fill( contrib.second, 0 ); + } + + // External field + if( this->idx_zeeman >= 0 ) + E_Zeeman( spins, contributions[idx_zeeman].second ); + + // Anisotropy + if( this->idx_anisotropy >= 0 ) + E_Anisotropy( spins, contributions[idx_anisotropy].second ); + + // Exchange + if( this->idx_exchange >= 0 ) + E_Exchange( spins, contributions[idx_exchange].second ); + + // DMI + if( this->idx_dmi >= 0 ) + E_DMI( spins, contributions[idx_dmi].second ); + + // DDI + if( this->idx_ddi >= 0 ) + E_DDI( spins, contributions[idx_ddi].second ); +} + +void Hamiltonian_Micromagnetic::E_Zeeman( const vectorfield & spins, scalarfield & Energy ) +{ +#pragma omp parallel for + for( int icell = 0; icell < geometry->n_cells_total; ++icell ) + { + if( check_atom_type( this->geometry->atom_types[icell] ) ) + Energy[icell] -= C::Joule * geometry->cell_volume * Ms * this->external_field_magnitude + * this->external_field_normal.dot( spins[icell] ); + } +} + +void Hamiltonian_Micromagnetic::E_Update( const vectorfield & spins, scalarfield & Energy, vectorfield & gradient ) +{ +#pragma omp parallel for + for( int icell = 0; icell < geometry->n_cells_total; ++icell ) + { + Energy[icell] -= 0.5 * Ms * gradient[icell].dot( spins[icell] ); + } +} + +void Hamiltonian_Micromagnetic::E_Anisotropy( const vectorfield & spins, scalarfield & Energy ) +{ +#pragma omp parallel for + for( int icell = 0; icell < geometry->n_cells_total; ++icell ) + { + Energy[icell] -= geometry->cell_volume * C::Joule * spins[icell].dot( anisotropy_tensor * spins[icell] ); + } +} + +void Hamiltonian_Micromagnetic::E_Exchange( const vectorfield & spins, scalarfield & Energy ) +{ + auto delta = geometry->cell_size; + for( unsigned int icell = 0; icell < geometry->n_cells_total; ++icell ) + { + Vector3 grad_n; + for( unsigned int alpha = 0; alpha < 3; ++alpha ) + { + int icell_plus = idx_from_pair( + icell, boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, + neigh[2 * alpha] ); + + int icell_minus = idx_from_pair( + icell, boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, + neigh[2 * alpha + 1] ); + + if( icell_plus >= 0 || icell_minus >= 0 ) + { + if( icell_plus == -1 ) + icell_plus = icell; + if( icell_minus == -1 ) + icell_minus = icell; + + grad_n = ( spins[icell_plus] - spins[icell_minus] ) / ( 2 * delta[alpha] ); + // meV/J * J/m * 1/m * 1/m * m^3 + Energy[icell] += C::Joule * geometry->cell_volume * ( grad_n.dot( exchange_tensor * grad_n ) ); + } + } + } +} + +void Hamiltonian_Micromagnetic::E_DMI( const vectorfield & spins, scalarfield & Energy ) +{ + // TODO: This implementation is very likely far from optimal (performance wise) + const auto & delta = geometry->cell_size; + + auto epsilon = []( int i, int j, int k ) { return -0.5 * ( j - i ) * ( k - j ) * ( i - k ); }; + + scalar mult = C::Joule * geometry->cell_volume; + +// Implements: epsilon_{mu,alpha,beta} * D_{mu, nu} * [ n_{alpha} dn_{beta}/dr_{nu} - n_{beta} dn_{alpha}/dr_{nu} ] +#pragma omp parallel for + for( int icell = 0; icell < geometry->n_cells_total; ++icell ) + { + for( int nu = 0; nu < 3; ++nu ) + { + int icell_plus = idx_from_pair( + icell, boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, + neigh[2 * nu] ); + + int icell_minus = idx_from_pair( + icell, boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, + neigh[2 * nu + 1] ); + + if( icell_plus >= 0 || icell_minus >= 0 ) + { + if( icell_plus == -1 ) + icell_plus = icell; + if( icell_minus == -1 ) + icell_minus = icell; + + // Todo: Why is there a factor of 2 difference to OOMMF? + Vector3 grad_n = ( spins[icell_plus] - spins[icell_minus] ) / ( 2 * delta[nu] ); + + for( int alpha = 0; alpha < 3; alpha++ ) + { + for( int beta = 0; beta < 3; beta++ ) + { + for( int mu = 0; mu < 3; mu++ ) + { + // meV/J * J/m^2 * 1/m + Energy[icell] + += mult * epsilon( mu, alpha, beta ) * dmi_tensor( mu, nu ) + * ( spins[icell][alpha] * grad_n[beta] - spins[icell][beta] * grad_n[alpha] ); + } + } + } + } + } + } +} + +void Hamiltonian_Micromagnetic::E_DDI( const vectorfield & spins, scalarfield & Energy ) +{ + if( this->ddi_method == DDI_Method::FFT ) + this->E_DDI_FFT( spins, Energy ); + else if( this->ddi_method == DDI_Method::Cutoff ) + { + if( ddi_cutoff_radius < 0 ) + this->E_DDI_Direct( spins, Energy ); + } +} + +void Hamiltonian_Micromagnetic::E_DDI_Direct( const vectorfield & spins, scalarfield & Energy ) +{ + vectorfield gradients_temp; + gradients_temp.resize( geometry->nos ); + Vectormath::fill( gradients_temp, { 0, 0, 0 } ); + this->Gradient_DDI_Direct( spins, gradients_temp ); + +#pragma omp parallel for + for( int ispin = 0; ispin < geometry->nos; ispin++ ) + { + Energy[ispin] += 0.5 * spins[ispin].dot( gradients_temp[ispin] ); + } +} + +scalar Hamiltonian_Micromagnetic::Energy_Single_Spin( int ispin, const vectorfield & spins ) +{ + scalar Energy = 0; + return Energy; +} + +void Hamiltonian_Micromagnetic::Gradient( const vectorfield & spins, vectorfield & gradient ) +{ + + // Set to zero + Vectormath::fill( gradient, { 0, 0, 0 } ); + this->Spatial_Gradient( spins ); + + // External field + this->Gradient_Zeeman( gradient ); + + // Anisotropy + this->Gradient_Anisotropy( spins, gradient ); + + // Exchange + this->Gradient_Exchange( spins, gradient ); + + // DMI + this->Gradient_DMI( spins, gradient ); + + // DDI + this->Gradient_DDI( spins, gradient ); + + // double energy=0; + // #pragma omp parallel for reduction(-:energy) + // for( int icell = 0; icell < geometry->n_cells_total; ++icell ) + // { + // energy -= 0.5 * Ms * gradient[icell].dot(spins[icell]); + // } + // printf("Energy total: %f\n", energy/ geometry->n_cells_total); +} + +void Hamiltonian_Micromagnetic::Gradient_Zeeman( vectorfield & gradient ) +{ + +#pragma omp parallel for + for( int icell = 0; icell < geometry->n_cells_total; ++icell ) + { + if( check_atom_type( this->geometry->atom_types[icell] ) ) + gradient[icell] + -= C::Joule * Ms * geometry->cell_volume * this->external_field_magnitude * this->external_field_normal; + } +} + +void Hamiltonian_Micromagnetic::Gradient_Anisotropy( const vectorfield & spins, vectorfield & gradient ) +{ + Vector3 temp1{ 1, 0, 0 }; + Vector3 temp2{ 0, 1, 0 }; + Vector3 temp3{ 0, 0, 1 }; +#pragma omp parallel for + for( int icell = 0; icell < geometry->n_cells_total; ++icell ) + { + gradient[icell] -= 2.0 * geometry->cell_volume * C::Joule * anisotropy_tensor * spins[icell]; + } +} + +void Hamiltonian_Micromagnetic::Gradient_Exchange( const vectorfield & spins, vectorfield & gradient ) +{ + auto & delta = geometry->cell_size; + +#pragma omp parallel for + for( unsigned int icell = 0; icell < geometry->n_cells_total; ++icell ) + { + for( int alpha = 0; alpha < 3; ++alpha ) + { + int icell_plus = idx_from_pair( + icell, boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, + neigh[2 * alpha] ); + + int icell_minus = idx_from_pair( + icell, boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, + neigh[2 * alpha + 1] ); + + if( icell_plus >= 0 ) + { + Vector3 grad_n = spatial_gradient[icell_plus].col( alpha ); + gradient[icell] + -= C::Joule * geometry->cell_volume * 2 * ( exchange_tensor * grad_n ) / ( 2 * delta[alpha] ); + } + else + { + Vector3 grad_n = spatial_gradient[icell].col( alpha ); + gradient[icell] + += C::Joule * geometry->cell_volume * 2 * ( exchange_tensor * grad_n ) / ( 2 * delta[alpha] ); + } + + if( icell_minus >= 0 ) + { + Vector3 grad_n = spatial_gradient[icell_minus].col( alpha ); + gradient[icell] + += C::Joule * geometry->cell_volume * 2 * ( exchange_tensor * grad_n ) / ( 2 * delta[alpha] ); + } + else + { + Vector3 grad_n = spatial_gradient[icell].col( alpha ); + gradient[icell] + -= C::Joule * geometry->cell_volume * 2 * ( exchange_tensor * grad_n ) / ( 2 * delta[alpha] ); + } + } + } +} + +void Hamiltonian_Micromagnetic::Spatial_Gradient( const vectorfield & spins ) +{ + auto & delta = geometry->cell_size; + + /* + dn1/dr1 dn1/dr2 dn1/dr3 + dn2/dr1 dn2/dr2 dn2/dr3 + dn3/dr1 dn3/dr2 dn3/dr3 + */ + +#pragma omp parallel for + for( unsigned int icell = 0; icell < geometry->n_cells_total; ++icell ) + { + spatial_gradient[icell] = Matrix3::Zero(); + for( unsigned int i = 0; i < 3; ++i ) + { + int icell_plus = idx_from_pair( + icell, boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, + neigh[2 * i] ); + int icell_minus = idx_from_pair( + icell, boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, + neigh[2 * i + 1] ); + + if( icell_plus >= 0 || icell_minus >= 0 ) + { + if( icell_plus == -1 ) + icell_plus = icell; + if( icell_minus == -1 ) + icell_minus = icell; + + spatial_gradient[icell].col( i ) += ( spins[icell_plus] - spins[icell_minus] ) / ( 2 * delta[i] ); + } + } + } +} + +void Hamiltonian_Micromagnetic::Gradient_DMI( const vectorfield & spins, vectorfield & gradient ) +{ + const auto & delta = geometry->cell_size; + + auto epsilon = []( int i, int j, int k ) { return -0.5 * ( j - i ) * ( k - j ) * ( i - k ); }; + + scalar mult = C::Joule * geometry->cell_volume; + +#pragma omp parallel for + for( unsigned int icell = 0; icell < geometry->n_cells_total; ++icell ) + { + for( unsigned int nu = 0; nu < 3; ++nu ) + { + int icell_plus = idx_from_pair( + icell, boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, + neigh[2 * nu] ); + + int icell_minus = idx_from_pair( + icell, boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, + neigh[2 * nu + 1] ); + + // Energy[icell] += mult * epsilon(mu, alpha, beta) * dmi_tensor(mu, nu) * ( spins[icell][alpha] * + // grad_n[beta] - spins[icell][beta] * grad_n[alpha] ); + + Vector3 grad_n = spatial_gradient[icell].col( nu ); + + for( int alpha = 0; alpha < 3; alpha++ ) + { + for( int beta = 0; beta < 3; beta++ ) + { + for( int mu = 0; mu < 3; mu++ ) + { + if( icell_plus >= 0 ) + { + gradient[icell][beta] -= mult * epsilon( mu, alpha, beta ) * dmi_tensor( mu, nu ) + * spins[icell_plus][alpha] / ( 2 * delta[nu] ); + gradient[icell][alpha] += mult * epsilon( mu, alpha, beta ) * dmi_tensor( mu, nu ) + * spins[icell_plus][beta] / ( 2 * delta[nu] ); + } + + if( icell_minus >= 0 ) + { + gradient[icell][beta] += mult * epsilon( mu, alpha, beta ) * dmi_tensor( mu, nu ) + * spins[icell_minus][alpha] / ( 2 * delta[nu] ); + gradient[icell][alpha] -= mult * epsilon( mu, alpha, beta ) * dmi_tensor( mu, nu ) + * spins[icell_minus][beta] / ( 2 * delta[nu] ); + } + + gradient[icell][alpha] + += mult * epsilon( mu, alpha, beta ) * dmi_tensor( mu, nu ) * ( grad_n[beta] ); + if( icell_plus < 0 ) + { + gradient[icell][beta] += mult * epsilon( mu, alpha, beta ) * dmi_tensor( mu, nu ) + * spins[icell][alpha] / ( 2 * delta[nu] ); + gradient[icell][alpha] -= mult * epsilon( mu, alpha, beta ) * dmi_tensor( mu, nu ) + * spins[icell][beta] / ( 2 * delta[nu] ); + } + + gradient[icell][beta] + -= mult * epsilon( mu, alpha, beta ) * dmi_tensor( mu, nu ) * ( grad_n[alpha] ); + if( icell_minus < 0 ) + { + gradient[icell][beta] -= mult * epsilon( mu, alpha, beta ) * dmi_tensor( mu, nu ) + * spins[icell][alpha] / ( 2 * delta[nu] ); + gradient[icell][alpha] += mult * epsilon( mu, alpha, beta ) * dmi_tensor( mu, nu ) + * spins[icell][beta] / ( 2 * delta[nu] ); + } + } + } + } + } + } +} + +void Hamiltonian_Micromagnetic::Gradient_DDI( const vectorfield & spins, vectorfield & gradient ) +{ + if( this->ddi_method == DDI_Method::FFT ) + this->Gradient_DDI_FFT( spins, gradient ); + + else if( this->ddi_method == DDI_Method::Cutoff ) + { + if( ddi_cutoff_radius < 0 ) + this->Gradient_DDI_Direct( spins, gradient ); + } +} + +void Hamiltonian_Micromagnetic::Gradient_DDI_Direct( const vectorfield & spins, vectorfield & gradient ) +{ + Vector3 delta = geometry->cell_size; + scalar mult = Constants_Micromagnetic::mu_0 * geometry->cell_volume * ( Ms ) * (Ms)*C::Joule; + + int img_a = boundary_conditions[0] == 0 ? 0 : ddi_n_periodic_images[0]; + int img_b = boundary_conditions[1] == 0 ? 0 : ddi_n_periodic_images[1]; + int img_c = boundary_conditions[2] == 0 ? 0 : ddi_n_periodic_images[2]; + + for( int idx1 = 0; idx1 < geometry->nos; idx1++ ) + { + for( int idx2 = 0; idx2 < geometry->nos; idx2++ ) + { + scalar Dxx = 0, Dxy = 0, Dxz = 0, Dyy = 0, Dyz = 0, Dzz = 0; + auto diff = this->geometry->positions[idx2] - this->geometry->positions[idx1]; + + for( int a_pb = -img_a; a_pb <= img_a; a_pb++ ) + { + for( int b_pb = -img_b; b_pb <= img_b; b_pb++ ) + { + for( int c_pb = -img_c; c_pb <= img_c; c_pb++ ) + { + scalar X = 1e-10 * diff[0] + geometry->n_cells[0] * a_pb * delta[0]; + scalar Y = 1e-10 * diff[1] + geometry->n_cells[1] * b_pb * delta[1]; + scalar Z = 1e-10 * diff[2] + geometry->n_cells[2] * c_pb * delta[2]; + + scalar dx = delta[0]; + scalar dy = delta[1]; + scalar dz = delta[2]; + + Dxx += mult * Demagnetization_Tensor::Automatic::Nxx( X, Y, Z, dx, dy, dz ); + Dxy += mult * Demagnetization_Tensor::Automatic::Nxy( X, Y, Z, dx, dy, dz ); + Dxz += mult * Demagnetization_Tensor::Automatic::Nxy( X, Z, Y, dx, dz, dy ); + Dyy += mult * Demagnetization_Tensor::Automatic::Nxx( Y, X, Z, dy, dx, dz ); + Dyz += mult * Demagnetization_Tensor::Automatic::Nxy( Z, Y, X, dz, dy, dx ); + Dzz += mult * Demagnetization_Tensor::Automatic::Nxx( Z, Y, X, dz, dy, dx ); + } + } + } + + gradient[idx1][0] -= ( Dxx * spins[idx2][0] + Dxy * spins[idx2][1] + Dxz * spins[idx2][2] ); + gradient[idx1][1] -= ( Dxy * spins[idx2][0] + Dyy * spins[idx2][1] + Dyz * spins[idx2][2] ); + gradient[idx1][2] -= ( Dxz * spins[idx2][0] + Dyz * spins[idx2][1] + Dzz * spins[idx2][2] ); + } + } +} + +void Hamiltonian_Micromagnetic::Gradient_DDI_FFT( const vectorfield & spins, vectorfield & gradient ) +{ + // Size of original geometry + int Na = geometry->n_cells[0]; + int Nb = geometry->n_cells[1]; + int Nc = geometry->n_cells[2]; + + auto cell_volume = geometry->cell_size[0] * geometry->cell_size[1] * geometry->cell_size[2]; + + FFT_Spins( spins ); + + auto & ft_D_matrices = transformed_dipole_matrices; + auto & ft_spins = fft_plan_spins.cpx_ptr; + + auto & res_iFFT = fft_plan_reverse.real_ptr; + auto & res_mult = fft_plan_reverse.cpx_ptr; + + int idx_s, idx_d; + + // Workaround for compability with intel compiler + const int c_n_cell_atoms = geometry->n_cell_atoms; + const int * c_it_bounds_pointwise_mult = it_bounds_pointwise_mult.data(); + +// Loop over basis atoms (i.e sublattices) +#pragma omp parallel for collapse( 3 ) + + for( int c = 0; c < c_it_bounds_pointwise_mult[2]; ++c ) + { + for( int b = 0; b < c_it_bounds_pointwise_mult[1]; ++b ) + { + for( int a = 0; a < c_it_bounds_pointwise_mult[0]; ++a ) + { + idx_s = a * spin_stride.a + b * spin_stride.b + c * spin_stride.c; + idx_d = a * dipole_stride.a + b * dipole_stride.b + c * dipole_stride.c; + + auto & fs_x = ft_spins[idx_s]; + auto & fs_y = ft_spins[idx_s + 1 * spin_stride.comp]; + auto & fs_z = ft_spins[idx_s + 2 * spin_stride.comp]; + + auto & fD_xx = ft_D_matrices[idx_d]; + auto & fD_xy = ft_D_matrices[idx_d + 1 * dipole_stride.comp]; + auto & fD_xz = ft_D_matrices[idx_d + 2 * dipole_stride.comp]; + auto & fD_yy = ft_D_matrices[idx_d + 3 * dipole_stride.comp]; + auto & fD_yz = ft_D_matrices[idx_d + 4 * dipole_stride.comp]; + auto & fD_zz = ft_D_matrices[idx_d + 5 * dipole_stride.comp]; + + FFT::addTo( + res_mult[idx_s + 0 * spin_stride.comp], FFT::mult3D( fD_xx, fD_xy, fD_xz, fs_x, fs_y, fs_z ), + true ); + FFT::addTo( + res_mult[idx_s + 1 * spin_stride.comp], FFT::mult3D( fD_xy, fD_yy, fD_yz, fs_x, fs_y, fs_z ), + true ); + FFT::addTo( + res_mult[idx_s + 2 * spin_stride.comp], FFT::mult3D( fD_xz, fD_yz, fD_zz, fs_x, fs_y, fs_z ), + true ); + } + } + } // end iteration over padded lattice cells + + // Inverse Fourier Transform + FFT::batch_iFour_3D( fft_plan_reverse ); + + // Workaround for compability with intel compiler + const int * c_n_cells = geometry->n_cells.data(); + + // Place the gradients at the correct positions and mult with correct mu + for( int c = 0; c < c_n_cells[2]; ++c ) + { + for( int b = 0; b < c_n_cells[1]; ++b ) + { + for( int a = 0; a < c_n_cells[0]; ++a ) + { + int idx_orig = a + Na * ( b + Nb * c ); + int idx = a * spin_stride.a + b * spin_stride.b + c * spin_stride.c; + gradient[idx_orig][0] -= res_iFFT[idx] / sublattice_size; + gradient[idx_orig][1] -= res_iFFT[idx + 1 * spin_stride.comp] / sublattice_size; + gradient[idx_orig][2] -= res_iFFT[idx + 2 * spin_stride.comp] / sublattice_size; + } + } + } // end iteration sublattice 1 +} + +void Hamiltonian_Micromagnetic::E_DDI_FFT( const vectorfield & spins, scalarfield & Energy ) +{ + scalar Energy_DDI = 0; + vectorfield gradients_temp; + gradients_temp.resize( geometry->nos ); + Vectormath::fill( gradients_temp, { 0, 0, 0 } ); + this->Gradient_DDI_FFT( spins, gradients_temp ); + + // === DEBUG: begin gradient comparison === + // vectorfield gradients_temp_dir; + // gradients_temp_dir.resize(this->geometry->nos); + // Vectormath::fill(gradients_temp_dir, {0,0,0}); + // Gradient_DDI_Direct(spins, gradients_temp_dir); + + // //get deviation + // Vector3 deviation = {0,0,0}; + // scalar max_deviation = 0; + + // std::array avg = {0,0,0}; + // for(int i = 0; i < this->geometry->nos; i++) + // { + // for(int d = 0; d < 3; d++) + // { + // deviation[d] += std::pow(gradients_temp[i][d] - gradients_temp_dir[i][d], 2); + // avg[d] += gradients_temp_dir[i][d]; + // } + // max_deviation = std::max( (gradients_temp_dir[i] - gradients_temp[i]).norm(), max_deviation ); + // } + // std::cerr << "Avg. Gradient = " << avg[0]/this->geometry->nos << " " << avg[1]/this->geometry->nos << " " << + // avg[2]/this->geometry->nos << std::endl; std::cerr << "Avg. Deviation = " << deviation[0]/this->geometry->nos + // << " " << deviation[1]/this->geometry->nos << " " << deviation[2]/this->geometry->nos << std::endl; + // std::cerr << "Max. Deviation = " << max_deviation << "\n"; + //==== DEBUG: end gradient comparison ==== + + // TODO: add dot_scaled to Vectormath and use that + for( int ispin = 0; ispin < geometry->nos; ispin++ ) + { + Energy[ispin] += 0.5 * spins[ispin].dot( gradients_temp[ispin] ); + // Energy_DDI += 0.5 * spins[ispin].dot(gradients_temp[ispin]); + } +} + +void Hamiltonian_Micromagnetic::FFT_Demag_Tensors( FFT::FFT_Plan & fft_plan_dipole, int img_a, int img_b, int img_c ) +{ + auto delta = geometry->cell_size; + + // Prefactor of DDI + // The energy is proportional to spin_direction * Demag_tensor * spin_direction + // The 'mult' factor is chosen such that the cell resolved energy has + // the dimension of total energy per cell in meV + + // mult has the units of [N / A^2] [m^3] [(A/m)^2] [mev/J] = [J] [meV/J] = [meV] + scalar mult = Constants_Micromagnetic::mu_0 * geometry->cell_volume * ( Ms ) * (Ms)*C::Joule; + + std::cout << "cell_size " << geometry->cell_size.transpose() << "\n"; + std::cout << "cell_volume " << geometry->cell_volume << "\n"; + std::cout << "mult " << mult << "\n"; + std::cout << "Ms " << Ms << "\n"; + + // Size of original geometry + int Na = geometry->n_cells[0]; + int Nb = geometry->n_cells[1]; + int Nc = geometry->n_cells[2]; + + auto & fft_dipole_inputs = fft_plan_dipole.real_ptr; + + // Iterate over the padded system + const int * c_n_cells_padded = n_cells_padded.data(); + +#pragma omp parallel for collapse( 3 ) + for( int c = 0; c < c_n_cells_padded[2]; ++c ) + { + for( int b = 0; b < c_n_cells_padded[1]; ++b ) + { + for( int a = 0; a < c_n_cells_padded[0]; ++a ) + { + int a_idx = a < Na ? a : a - n_cells_padded[0]; + int b_idx = b < Nb ? b : b - n_cells_padded[1]; + int c_idx = c < Nc ? c : c - n_cells_padded[2]; + + scalar Dxx = 0, Dxy = 0, Dxz = 0, Dyy = 0, Dyz = 0, Dzz = 0; + + // Iterate over periodic images + for( int a_pb = -img_a; a_pb <= img_a; a_pb++ ) + { + for( int b_pb = -img_b; b_pb <= img_b; b_pb++ ) + { + for( int c_pb = -img_c; c_pb <= img_c; c_pb++ ) + { + scalar X = ( a_idx + a_pb * Na ) * delta[0]; + scalar Y = ( b_idx + b_pb * Nb ) * delta[1]; + scalar Z = ( c_idx + c_pb * Nc ) * delta[2]; + + scalar dx = delta[0]; + scalar dy = delta[1]; + scalar dz = delta[2]; + + Dxx += mult * Demagnetization_Tensor::Automatic::Nxx( X, Y, Z, dx, dy, dz ); + Dxy += mult * Demagnetization_Tensor::Automatic::Nxy( X, Y, Z, dx, dy, dz ); + Dxz += mult * Demagnetization_Tensor::Automatic::Nxy( X, Z, Y, dx, dz, dy ); + Dyy += mult * Demagnetization_Tensor::Automatic::Nxx( Y, X, Z, dy, dx, dz ); + Dyz += mult * Demagnetization_Tensor::Automatic::Nxy( Z, Y, X, dz, dy, dx ); + Dzz += mult * Demagnetization_Tensor::Automatic::Nxx( Z, Y, X, dz, dy, dx ); + } + } + } + + int idx = a * dipole_stride.a + b * dipole_stride.b + c * dipole_stride.c; + + fft_dipole_inputs[idx] = Dxx; + fft_dipole_inputs[idx + 1 * dipole_stride.comp] = Dxy; + fft_dipole_inputs[idx + 2 * dipole_stride.comp] = Dxz; + fft_dipole_inputs[idx + 3 * dipole_stride.comp] = Dyy; + fft_dipole_inputs[idx + 4 * dipole_stride.comp] = Dyz; + fft_dipole_inputs[idx + 5 * dipole_stride.comp] = Dzz; + } + } + } + FFT::batch_Four_3D( fft_plan_dipole ); +} + +void Hamiltonian_Micromagnetic::FFT_Spins( const vectorfield & spins ) +{ + // size of original geometry + int Na = geometry->n_cells[0]; + int Nb = geometry->n_cells[1]; + int Nc = geometry->n_cells[2]; + + auto cell_volume = geometry->cell_size[0] * geometry->cell_size[1] * geometry->cell_size[2]; + + auto & fft_spin_inputs = fft_plan_spins.real_ptr; + +// iterate over the **original** system +#pragma omp parallel for collapse( 3 ) + for( int c = 0; c < Nc; ++c ) + { + for( int b = 0; b < Nb; ++b ) + { + for( int a = 0; a < Na; ++a ) + { + int idx_orig = a + Na * ( b + Nb * c ); + int idx = a * spin_stride.a + b * spin_stride.b + c * spin_stride.c; + + fft_spin_inputs[idx] = spins[idx_orig][0]; + fft_spin_inputs[idx + 1 * spin_stride.comp] = spins[idx_orig][1]; + fft_spin_inputs[idx + 2 * spin_stride.comp] = spins[idx_orig][2]; + } + } + } + + FFT::batch_Four_3D( fft_plan_spins ); +} + +void Hamiltonian_Micromagnetic::Prepare_DDI() +{ + Clean_DDI(); + + if( ddi_method != DDI_Method::FFT ) + return; + + // We perform zero-padding in a lattice direction if the dimension of the system is greater than 1 *and* + // - the boundary conditions are open, or + // - the boundary conditions are periodic and zero-padding is explicitly requested + n_cells_padded.resize( 3 ); + for( int i = 0; i < 3; i++ ) + { + n_cells_padded[i] = geometry->n_cells[i]; + bool perform_zero_padding = geometry->n_cells[i] > 1 && ( boundary_conditions[i] == 0 || ddi_pb_zero_padding ); + if( perform_zero_padding ) + n_cells_padded[i] *= 2; + } + sublattice_size = n_cells_padded[0] * n_cells_padded[1] * n_cells_padded[2]; + + FFT::FFT_Init(); + +// Workaround for bug in kissfft +// kissfft_ndr does not perform one-dimensional FFTs properly +#ifndef SPIRIT_USE_FFTW + int number_of_one_dims = 0; + for( int i = 0; i < 3; i++ ) + if( n_cells_padded[i] == 1 && ++number_of_one_dims > 1 ) + n_cells_padded[i] = 2; +#endif + + sublattice_size = n_cells_padded[0] * n_cells_padded[1] * n_cells_padded[2]; + + // We dont need to transform over length 1 dims + std::vector fft_dims; + for( int i = 2; i >= 0; i-- ) // notice that reverse order is important! + { + if( n_cells_padded[i] > 1 ) + fft_dims.push_back( n_cells_padded[i] ); + } + + // Create FFT plans + FFT::FFT_Plan fft_plan_dipole = FFT::FFT_Plan( fft_dims, false, 6, sublattice_size ); + fft_plan_spins = FFT::FFT_Plan( fft_dims, false, 3, sublattice_size ); + fft_plan_reverse = FFT::FFT_Plan( fft_dims, true, 3, sublattice_size ); + +#ifdef SPIRIT_USE_FFTW + field temp_s = { &spin_stride.comp, &spin_stride.basis, &spin_stride.a, &spin_stride.b, &spin_stride.c }; + field temp_d + = { &dipole_stride.comp, &dipole_stride.basis, &dipole_stride.a, &dipole_stride.b, &dipole_stride.c }; + + FFT::get_strides( + temp_s, { 3, this->geometry->n_cell_atoms, n_cells_padded[0], n_cells_padded[1], n_cells_padded[2] } ); + FFT::get_strides( temp_d, { 6, 1, n_cells_padded[0], n_cells_padded[1], n_cells_padded[2] } ); + it_bounds_pointwise_mult = { ( n_cells_padded[0] / 2 + 1 ), // due to redundancy in real fft + n_cells_padded[1], n_cells_padded[2] }; +#else + field temp_s = { &spin_stride.a, &spin_stride.b, &spin_stride.c, &spin_stride.comp, &spin_stride.basis }; + field temp_d + = { &dipole_stride.a, &dipole_stride.b, &dipole_stride.c, &dipole_stride.comp, &dipole_stride.basis }; + + FFT::get_strides( + temp_s, { n_cells_padded[0], n_cells_padded[1], n_cells_padded[2], 3, this->geometry->n_cell_atoms } ); + FFT::get_strides( temp_d, { n_cells_padded[0], n_cells_padded[1], n_cells_padded[2], 6, 1 } ); + it_bounds_pointwise_mult = { n_cells_padded[0], n_cells_padded[1], n_cells_padded[2] }; + ( it_bounds_pointwise_mult[fft_dims.size() - 1] /= 2 )++; +#endif + + // Perform FFT of dipole matrices + int img_a = boundary_conditions[0] == 0 ? 0 : ddi_n_periodic_images[0]; + int img_b = boundary_conditions[1] == 0 ? 0 : ddi_n_periodic_images[1]; + int img_c = boundary_conditions[2] == 0 ? 0 : ddi_n_periodic_images[2]; + + FFT_Demag_Tensors( fft_plan_dipole, img_a, img_b, img_c ); + transformed_dipole_matrices = std::move( fft_plan_dipole.cpx_ptr ); +} + +void Hamiltonian_Micromagnetic::Clean_DDI() +{ + fft_plan_spins = FFT::FFT_Plan(); + fft_plan_reverse = FFT::FFT_Plan(); +} + +void Hamiltonian_Micromagnetic::Hessian( const vectorfield & spins, MatrixX & hessian ) {} + +// Hamiltonian name as string +static const std::string name = "Micromagnetic"; +const std::string & Hamiltonian_Micromagnetic::Name() const +{ + return name; +} + +} // namespace Engine + +#endif \ No newline at end of file diff --git a/core/src/engine/Hamiltonian_Micromagnetic.cu b/core/src/engine/Hamiltonian_Micromagnetic.cu new file mode 100644 index 000000000..7bb001ba1 --- /dev/null +++ b/core/src/engine/Hamiltonian_Micromagnetic.cu @@ -0,0 +1,1085 @@ +#ifdef SPIRIT_USE_CUDA + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +using namespace Data; +using namespace Utility; +namespace C = Utility::Constants_Micromagnetic; +using Engine::Vectormath::check_atom_type; +using Engine::Vectormath::cu_check_atom_type; +using Engine::Vectormath::cu_idx_from_pair; +using Engine::Vectormath::cu_tupel_from_idx; +using Engine::Vectormath::idx_from_pair; + +namespace Engine +{ + +Hamiltonian_Micromagnetic::Hamiltonian_Micromagnetic( + scalar Ms, scalar external_field_magnitude, Vector3 external_field_normal, Matrix3 anisotropy_tensor, + Matrix3 exchange_tensor, Matrix3 dmi_tensor, DDI_Method ddi_method, intfield ddi_n_periodic_images, + scalar ddi_radius, std::shared_ptr geometry, int spatial_gradient_order, + intfield boundary_conditions ) + : Hamiltonian( boundary_conditions ), + spatial_gradient_order( spatial_gradient_order ), + geometry( geometry ), + external_field_magnitude( external_field_magnitude ), + external_field_normal( external_field_normal ), + anisotropy_tensor( anisotropy_tensor ), + exchange_tensor( exchange_tensor ), + dmi_tensor( dmi_tensor ) +{ + // Generate interaction pairs, constants etc. + this->Update_Interactions(); +} + +void Hamiltonian_Micromagnetic::Update_Interactions() +{ +#if defined( SPIRIT_USE_OPENMP ) + // When parallelising (cuda or openmp), we need all neighbours per spin + const bool use_redundant_neighbours = true; +#else + // When running on a single thread, we can ignore redundant neighbours + const bool use_redundant_neighbours = false; +#endif + + // TODO: make sure that the geometry can be treated with this model: + // - rectilinear, only one "atom" per cell + // if( geometry->n_cell_atoms != 1 ) + // Log(...) + + // TODO: generate neighbour information for pairwise interactions + + // TODO: prepare dipolar interactions + neigh = pairfield( 0 ); + Neighbour neigh_tmp; + neigh_tmp.i = 0; + neigh_tmp.j = 0; + neigh_tmp.idx_shell = 0; + // order x -x y -y z -z xy (-x)(-y) x(-y) (-x)y xz (-x)(-z) x(-z) (-x)z yz (-y)(-z) y(-z) (-y)z results in 9 parts of Hessian + neigh_tmp.translations[0] = 1; + neigh_tmp.translations[1] = 0; + neigh_tmp.translations[2] = 0; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = -1; + neigh_tmp.translations[1] = 0; + neigh_tmp.translations[2] = 0; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = 0; + neigh_tmp.translations[1] = 1; + neigh_tmp.translations[2] = 0; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = 0; + neigh_tmp.translations[1] = -1; + neigh_tmp.translations[2] = 0; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = 0; + neigh_tmp.translations[1] = 0; + neigh_tmp.translations[2] = 1; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = 0; + neigh_tmp.translations[1] = 0; + neigh_tmp.translations[2] = -1; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = 1; + neigh_tmp.translations[1] = 1; + neigh_tmp.translations[2] = 0; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = -1; + neigh_tmp.translations[1] = -1; + neigh_tmp.translations[2] = 0; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = 1; + neigh_tmp.translations[1] = -1; + neigh_tmp.translations[2] = 0; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = -1; + neigh_tmp.translations[1] = +1; + neigh_tmp.translations[2] = 0; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = 1; + neigh_tmp.translations[1] = 0; + neigh_tmp.translations[2] = 1; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = -1; + neigh_tmp.translations[1] = 0; + neigh_tmp.translations[2] = -1; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = 1; + neigh_tmp.translations[1] = 0; + neigh_tmp.translations[2] = -1; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = -1; + neigh_tmp.translations[1] = 0; + neigh_tmp.translations[2] = 1; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = 0; + neigh_tmp.translations[1] = 1; + neigh_tmp.translations[2] = 1; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = 0; + neigh_tmp.translations[1] = -1; + neigh_tmp.translations[2] = -1; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = 0; + neigh_tmp.translations[1] = 1; + neigh_tmp.translations[2] = -1; + neigh.push_back( neigh_tmp ); + + neigh_tmp.translations[0] = 0; + neigh_tmp.translations[1] = -1; + neigh_tmp.translations[2] = 1; + neigh.push_back( neigh_tmp ); + this->spatial_gradient = field( geometry->nos, Matrix3::Zero() ); + + // Dipole-dipole + this->Prepare_DDI(); + + // Update, which terms still contribute + this->Update_Energy_Contributions(); +} + +void Hamiltonian_Micromagnetic::Update_Energy_Contributions() +{ + this->energy_contributions_per_spin = std::vector>( 0 ); + CU_CHECK_AND_SYNC(); + // External field + if( this->external_field_magnitude > 0 ) + { + this->energy_contributions_per_spin.push_back( { "Zeeman", scalarfield( 0 ) } ); + this->idx_zeeman = this->energy_contributions_per_spin.size() - 1; + } + else + this->idx_zeeman = -1; + // TODO: Anisotropy + // if( ... ) + // { + // this->energy_contributions_per_spin.push_back({"Anisotropy", scalarfield(0) }); + // this->idx_anisotropy = this->energy_contributions_per_spin.size()-1; + // } + // else + this->idx_anisotropy = -1; + // TODO: Exchange + // if( ... ) + // { + // this->energy_contributions_per_spin.push_back({"Exchange", scalarfield(0) }); + // this->idx_exchange = this->energy_contributions_per_spin.size()-1; + // } + // else + this->idx_exchange = -1; + // TODO: DMI + // if( ... ) + // { + // this->energy_contributions_per_spin.push_back({"DMI", scalarfield(0) }); + // this->idx_dmi = this->energy_contributions_per_spin.size()-1; + // } + // else + this->idx_dmi = -1; + // TODO: DDI + // if( ... ) + // { + // this->energy_contributions_per_spin.push_back({"DDI", scalarfield(0) }); + // this->idx_ddi = this->energy_contributions_per_spin.size()-1; + // } + // else + this->idx_ddi = -1; +} + +void Hamiltonian_Micromagnetic::Energy_Contributions_per_Spin( + const vectorfield & spins, std::vector> & contributions ) +{ + if( contributions.size() != this->energy_contributions_per_spin.size() ) + { + contributions = this->energy_contributions_per_spin; + } + + int nos = spins.size(); + for( auto & contrib : contributions ) + { + // Allocate if not already allocated + if( contrib.second.size() != nos ) + contrib.second = scalarfield( nos, 0 ); + // Otherwise set to zero + else + Vectormath::fill( contrib.second, 0 ); + } + + // External field + if( this->idx_zeeman >= 0 ) + E_Zeeman( spins, contributions[idx_zeeman].second ); + + // Anisotropy + if( this->idx_anisotropy >= 0 ) + E_Anisotropy( spins, contributions[idx_anisotropy].second ); + + // Exchange + if( this->idx_exchange >= 0 ) + E_Exchange( spins, contributions[idx_exchange].second ); + // DMI + if( this->idx_dmi >= 0 ) + E_DMI( spins, contributions[idx_dmi].second ); +} + +__global__ void CU_E_Zeeman1( + const Vector3 * spins, const int * atom_types, const int n_cell_atoms, const scalar * mu_s, + const scalar external_field_magnitude, const Vector3 external_field_normal, scalar * Energy, size_t n_cells_total ) +{ + for( auto icell = blockIdx.x * blockDim.x + threadIdx.x; icell < n_cells_total; icell += blockDim.x * gridDim.x ) + { + for( int ibasis = 0; ibasis < n_cell_atoms; ++ibasis ) + { + int ispin = icell + ibasis; + if( cu_check_atom_type( atom_types[ispin] ) ) + Energy[ispin] -= mu_s[ispin] * external_field_magnitude * external_field_normal.dot( spins[ispin] ); + } + } +} +void Hamiltonian_Micromagnetic::E_Zeeman( const vectorfield & spins, scalarfield & Energy ) +{ + int size = geometry->n_cells_total; + CU_E_Zeeman1<<<( size + 1023 ) / 1024, 1024>>>( + spins.data(), this->geometry->atom_types.data(), geometry->n_cell_atoms, geometry->mu_s.data(), + this->external_field_magnitude, this->external_field_normal, Energy.data(), size ); + CU_CHECK_AND_SYNC(); +} + +void Hamiltonian_Micromagnetic::E_Anisotropy( const vectorfield & spins, scalarfield & Energy ) {} + +void Hamiltonian_Micromagnetic::E_Exchange( const vectorfield & spins, scalarfield & Energy ) {} + +void Hamiltonian_Micromagnetic::E_DMI( const vectorfield & spins, scalarfield & Energy ) {} + +void Hamiltonian_Micromagnetic::E_DDI( const vectorfield & spins, scalarfield & Energy ) {} + +scalar Hamiltonian_Micromagnetic::Energy_Single_Spin( int ispin, const vectorfield & spins ) +{ + scalar Energy = 0; + return Energy; +} + +void Hamiltonian_Micromagnetic::Gradient( const vectorfield & spins, vectorfield & gradient ) +{ + // Set to zero + Vectormath::fill( gradient, { 0, 0, 0 } ); + this->Spatial_Gradient( spins ); + // External field + this->Gradient_Zeeman( gradient ); + + // Anisotropy + this->Gradient_Anisotropy( spins, gradient ); + + // Exchange + this->Gradient_Exchange( spins, gradient ); + + // DMI + this->Gradient_DMI( spins, gradient ); + scalar Ms = 1.4e6; + double energy = 0; +#pragma omp parallel for reduction( - : energy ) + for( int icell = 0; icell < geometry->n_cells_total; ++icell ) + { + // energy -= 0.5 *Ms* gradient[icell].dot(spins[icell]); + } + // printf("Energy total: %f\n", energy/ geometry->n_cells_total); +} + +__global__ void CU_Gradient_Zeeman1( + const int * atom_types, const int n_cell_atoms, const scalar * mu_s, const scalar external_field_magnitude, + const Vector3 external_field_normal, Vector3 * gradient, size_t n_cells_total ) +{ + for( auto icell = blockIdx.x * blockDim.x + threadIdx.x; icell < n_cells_total; icell += blockDim.x * gridDim.x ) + { + for( int ibasis = 0; ibasis < n_cell_atoms; ++ibasis ) + { + int ispin = icell + ibasis; + if( cu_check_atom_type( atom_types[ispin] ) ) + gradient[ispin] -= mu_s[ispin] * C::mu_B * external_field_magnitude * external_field_normal; + } + } +} +void Hamiltonian_Micromagnetic::Gradient_Zeeman( vectorfield & gradient ) +{ + int size = geometry->n_cells_total; + CU_Gradient_Zeeman1<<<( size + 1023 ) / 1024, 1024>>>( + this->geometry->atom_types.data(), geometry->n_cell_atoms, geometry->mu_s.data(), + this->external_field_magnitude, this->external_field_normal, gradient.data(), size ); + CU_CHECK_AND_SYNC(); +} + +__global__ void CU_Gradient_Anisotropy1( + const Vector3 * spins, const int * atom_types, const int n_cell_atoms, Vector3 * gradient, size_t n_cells_total, + Matrix3 anisotropy_tensor ) +{ + scalar Ms = 1.4e6; + Vector3 temp1{ 1, 0, 0 }; + Vector3 temp2{ 0, 1, 0 }; + Vector3 temp3{ 0, 0, 1 }; + for( auto icell = blockIdx.x * blockDim.x + threadIdx.x; icell < n_cells_total; icell += blockDim.x * gridDim.x ) + { + int ispin = icell; + gradient[ispin] -= 2.0 * C::mu_B * anisotropy_tensor * spins[ispin] / Ms; + // gradient[ispin] -= 2.0 * this->anisotropy_magnitudes[iani] / Ms * ((pow(temp2.dot(spins[ispin]),2)+ + // pow(temp3.dot(spins[ispin]), 2))*(temp1.dot(spins[ispin])*temp1)+ (pow(temp1.dot(spins[ispin]), 2) + + // pow(temp3.dot(spins[ispin]), 2))*(temp2.dot(spins[ispin])*temp2)+(pow(temp1.dot(spins[ispin]),2)+ + // pow(temp2.dot(spins[ispin]), 2))*(temp3.dot(spins[ispin])*temp3)); gradient[ispin] += 2.0 * 50000 / Ms * + // ((pow(temp2.dot(spins[ispin]), 2) + pow(temp3.dot(spins[ispin]), 2))*(temp1.dot(spins[ispin])*temp1) + + // (pow(temp1.dot(spins[ispin]), 2) + pow(temp3.dot(spins[ispin]), 2))*(temp2.dot(spins[ispin])*temp2)); + } +} + +void Hamiltonian_Micromagnetic::Gradient_Anisotropy( const vectorfield & spins, vectorfield & gradient ) +{ + int size = geometry->n_cells_total; + CU_Gradient_Anisotropy1<<<( size + 1023 ) / 1024, 1024>>>( + spins.data(), this->geometry->atom_types.data(), this->geometry->n_cell_atoms, gradient.data(), size, + this->anisotropy_tensor ); + CU_CHECK_AND_SYNC(); +} + +__global__ void CU_Gradient_Exchange1( + const Vector3 * spins, const int * atom_types, const int * boundary_conditions, const int * n_cells, + int n_cell_atoms, int n_pairs, const Pair * neigh, Vector3 * gradient, size_t size, bool A_is_nondiagonal, + Matrix3 exchange_tensor, const scalar * delta, const scalar Ms ) +{ + int bc[3] = { boundary_conditions[0], boundary_conditions[1], boundary_conditions[2] }; + + int nc[3] = { n_cells[0], n_cells[1], n_cells[2] }; + + for( auto icell = blockIdx.x * blockDim.x + threadIdx.x; icell < size; icell += blockDim.x * gridDim.x ) + { + // int ispin = icell;//basically id of a cell + for( unsigned int i = 0; i < 3; ++i ) + { + + int icell_plus = cu_idx_from_pair( icell, bc, nc, n_cell_atoms, atom_types, neigh[2 * i] ); + int icell_minus = cu_idx_from_pair( icell, bc, nc, n_cell_atoms, atom_types, neigh[2 * i + 1] ); + + if( icell_plus >= 0 || icell_minus >= 0 ) + { + if( icell_plus == -1 ) + icell_plus = icell; + if( icell_minus == -1 ) + icell_minus = icell; + + gradient[icell] -= 2 * C::mu_B * exchange_tensor + * ( spins[icell_plus] - 2 * spins[icell] + spins[icell_minus] ) + / ( Ms * delta[i] * delta[i] ); + } + } + /*if (A_is_nondiagonal == true) { + //xy + int ispin_right = idx_from_pair(ispin, boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, + geometry->atom_types, neigh[0]); int ispin_left = idx_from_pair(ispin, boundary_conditions, geometry->n_cells, + geometry->n_cell_atoms, geometry->atom_types, neigh[1]); int ispin_top = idx_from_pair(ispin, + boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, neigh[2]); int + ispin_bottom = idx_from_pair(ispin, boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, + geometry->atom_types, neigh[3]); + + if (ispin_right == -1) { + ispin_right = ispin; + } + if (ispin_left == -1) { + ispin_left = ispin; + } + if (ispin_top == -1) { + ispin_top = ispin; + } + if (ispin_bottom == -1) { + ispin_bottom = ispin; + } + gradient[ispin][0] -= 2 * exchange_tensor(0, 1) / Ms * ((spatial_gradient[ispin_top](0, 0) - + spatial_gradient[ispin_bottom](0, 0)) / 4 / delta[1] + (spatial_gradient[ispin_right](0, 1) - + spatial_gradient[ispin_left](0, 1)) / 4 / delta[0]); gradient[ispin][0] -= 2 * exchange_tensor(1, 0) / Ms * + ((spatial_gradient[ispin_top](0, 0) - spatial_gradient[ispin_bottom](0, 0)) / 4 / delta[1] + + (spatial_gradient[ispin_right](0, 1) - spatial_gradient[ispin_left](0, 1)) / 4 / delta[0]); gradient[ispin][1] + -= 2 * exchange_tensor(0, 1) / Ms * ((spatial_gradient[ispin_top](1, 0) - spatial_gradient[ispin_bottom](1, 0)) + / 4 / delta[1] + (spatial_gradient[ispin_right](1, 1) - spatial_gradient[ispin_left](1, 1)) / 4 / delta[0]); + gradient[ispin][1] -= 2 * exchange_tensor(1, 0) / Ms * ((spatial_gradient[ispin_top](1, 0) - + spatial_gradient[ispin_bottom](1, 0)) / 4 / delta[1] + (spatial_gradient[ispin_right](1, 1) - + spatial_gradient[ispin_left](1, 1)) / 4 / delta[0]); gradient[ispin][2] -= 2 * exchange_tensor(0, 1) / Ms * + ((spatial_gradient[ispin_top](2, 0) - spatial_gradient[ispin_bottom](2, 0)) / 4 / delta[1] + + (spatial_gradient[ispin_right](2, 1) - spatial_gradient[ispin_left](2, 1)) / 4 / delta[0]); gradient[ispin][2] + -= 2 * exchange_tensor(1, 0) / Ms * ((spatial_gradient[ispin_top](2, 0) - spatial_gradient[ispin_bottom](2, 0)) + / 4 / delta[1] + (spatial_gradient[ispin_right](2, 1) - spatial_gradient[ispin_left](2, 1)) / 4 / delta[0]); + + //xz + ispin_right = idx_from_pair(ispin, boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, + geometry->atom_types, neigh[0]); ispin_left = idx_from_pair(ispin, boundary_conditions, geometry->n_cells, + geometry->n_cell_atoms, geometry->atom_types, neigh[1]); ispin_top = idx_from_pair(ispin, boundary_conditions, + geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, neigh[4]); ispin_bottom = idx_from_pair(ispin, + boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, neigh[5]); + + if (ispin_right == -1) { + ispin_right = ispin; + } + if (ispin_left == -1) { + ispin_left = ispin; + } + if (ispin_top == -1) { + ispin_top = ispin; + } + if (ispin_bottom == -1) { + ispin_bottom = ispin; + } + gradient[ispin][0] -= 2 * exchange_tensor(0, 2) / Ms * ((spatial_gradient[ispin_top](0, 0) - + spatial_gradient[ispin_bottom](0, 0)) / 4 / delta[2] + (spatial_gradient[ispin_right](0, 2) - + spatial_gradient[ispin_left](0, 2)) / 4 / delta[0]); gradient[ispin][0] -= 2 * exchange_tensor(2, 0) / Ms * + ((spatial_gradient[ispin_top](0, 0) - spatial_gradient[ispin_bottom](0, 0)) / 4 / delta[2] + + (spatial_gradient[ispin_right](0, 2) - spatial_gradient[ispin_left](0, 2)) / 4 / delta[0]); gradient[ispin][1] + -= 2 * exchange_tensor(0, 2) / Ms * ((spatial_gradient[ispin_top](1, 0) - spatial_gradient[ispin_bottom](1, 0)) + / 4 / delta[2] + (spatial_gradient[ispin_right](1, 2) - spatial_gradient[ispin_left](1, 2)) / 4 / delta[0]); + gradient[ispin][1] -= 2 * exchange_tensor(2, 0) / Ms * ((spatial_gradient[ispin_top](1, 0) - + spatial_gradient[ispin_bottom](1, 0)) / 4 / delta[2] + (spatial_gradient[ispin_right](1, 2) - + spatial_gradient[ispin_left](1, 2)) / 4 / delta[0]); gradient[ispin][2] -= 2 * exchange_tensor(0, 2) / Ms * + ((spatial_gradient[ispin_top](2, 0) - spatial_gradient[ispin_bottom](2, 0)) / 4 / delta[2] + + (spatial_gradient[ispin_right](2, 2) - spatial_gradient[ispin_left](2, 2)) / 4 / delta[0]); gradient[ispin][2] + -= 2 * exchange_tensor(2, 0) / Ms * ((spatial_gradient[ispin_top](2, 0) - spatial_gradient[ispin_bottom](2, 0)) + / 4 / delta[2] + (spatial_gradient[ispin_right](2, 2) - spatial_gradient[ispin_left](2, 2)) / 4 / delta[0]); + + //yz + ispin_right = idx_from_pair(ispin, boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, + geometry->atom_types, neigh[2]); ispin_left = idx_from_pair(ispin, boundary_conditions, geometry->n_cells, + geometry->n_cell_atoms, geometry->atom_types, neigh[3]); ispin_top = idx_from_pair(ispin, boundary_conditions, + geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, neigh[4]); ispin_bottom = idx_from_pair(ispin, + boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, neigh[5]); + + if (ispin_right == -1) { + ispin_right = ispin; + } + if (ispin_left == -1) { + ispin_left = ispin; + } + if (ispin_top == -1) { + ispin_top = ispin; + } + if (ispin_bottom == -1) { + ispin_bottom = ispin; + } + gradient[ispin][0] -= 2 * exchange_tensor(1, 2) / Ms * ((spatial_gradient[ispin_top](0, 1) - + spatial_gradient[ispin_bottom](0, 1)) / 4 / delta[2] + (spatial_gradient[ispin_right](0, 2) - + spatial_gradient[ispin_left](0, 2)) / 4 / delta[0]); gradient[ispin][0] -= 2 * exchange_tensor(2, 1) / Ms * + ((spatial_gradient[ispin_top](0, 1) - spatial_gradient[ispin_bottom](0, 1)) / 4 / delta[2] + + (spatial_gradient[ispin_right](0, 2) - spatial_gradient[ispin_left](0, 2)) / 4 / delta[0]); gradient[ispin][1] + -= 2 * exchange_tensor(1, 2) / Ms * ((spatial_gradient[ispin_top](1, 1) - spatial_gradient[ispin_bottom](1, 1)) + / 4 / delta[2] + (spatial_gradient[ispin_right](1, 2) - spatial_gradient[ispin_left](1, 2)) / 4 / delta[0]); + gradient[ispin][1] -= 2 * exchange_tensor(2, 1) / Ms * ((spatial_gradient[ispin_top](1, 1) - + spatial_gradient[ispin_bottom](1, 1)) / 4 / delta[2] + (spatial_gradient[ispin_right](1, 2) - + spatial_gradient[ispin_left](1, 2)) / 4 / delta[0]); gradient[ispin][2] -= 2 * exchange_tensor(1, 2) / Ms * + ((spatial_gradient[ispin_top](2, 1) - spatial_gradient[ispin_bottom](2, 1)) / 4 / delta[2] + + (spatial_gradient[ispin_right](2, 2) - spatial_gradient[ispin_left](2, 2)) / 4 / delta[0]); gradient[ispin][2] + -= 2 * exchange_tensor(2, 1) / Ms * ((spatial_gradient[ispin_top](2, 1) - spatial_gradient[ispin_bottom](2, 1)) + / 4 / delta[2] + (spatial_gradient[ispin_right](2, 2) - spatial_gradient[ispin_left](2, 2)) / 4 / delta[0]); + + }*/ + } +} +void Hamiltonian_Micromagnetic::Gradient_Exchange( const vectorfield & spins, vectorfield & gradient ) +{ + int size = geometry->n_cells_total; + scalar * delta = geometry->cell_size.data(); + CU_Gradient_Exchange1<<<( size + 1023 ) / 1024, 1024>>>( + spins.data(), this->geometry->atom_types.data(), boundary_conditions.data(), geometry->n_cells.data(), + geometry->n_cell_atoms, this->neigh.size(), this->neigh.data(), gradient.data(), size, A_is_nondiagonal, + exchange_tensor, delta, this->Ms ); + CU_CHECK_AND_SYNC(); +} + +__global__ void CU_Spatial_Gradient( + const Vector3 * spins, const int * atom_types, const int * boundary_conditions, const int * n_cells, + int n_cell_atoms, int n_pairs, const Pair * neigh, Matrix3 * spatial_gradient, size_t size, scalar * delta, + scalar Ms ) +{ + /* + dn1/dr1 dn1/dr2 dn1/dr3 + dn2/dr1 dn2/dr2 dn2/dr3 + dn3/dr1 dn3/dr2 dn3/dr3 + */ + int bc[3] = { boundary_conditions[0], boundary_conditions[1], boundary_conditions[2] }; + int nc[3] = { n_cells[0], n_cells[1], n_cells[2] }; + + for( auto icell = blockIdx.x * blockDim.x + threadIdx.x; icell < size; icell += blockDim.x * gridDim.x ) + { + for( unsigned int i = 0; i < 3; ++i ) + { + int icell_plus = cu_idx_from_pair( icell, bc, nc, n_cell_atoms, atom_types, neigh[2 * i] ); + int icell_minus = cu_idx_from_pair( icell, bc, nc, n_cell_atoms, atom_types, neigh[2 * i + 1] ); + + if( icell_plus >= 0 || icell_minus >= 0 ) + { + if( icell_plus == -1 ) + icell_plus = icell; + if( icell_minus == -1 ) + icell_minus = icell; + + spatial_gradient[icell].col( i ) += ( spins[icell_plus] - spins[icell_minus] ) / ( 2 * delta[i] ); + } + } + } +} + +void Hamiltonian_Micromagnetic::Spatial_Gradient( const vectorfield & spins ) +{ + int size = geometry->n_cells_total; + CU_Spatial_Gradient<<<( size + 1023 ) / 1024, 1024>>>( + spins.data(), this->geometry->atom_types.data(), boundary_conditions.data(), geometry->n_cells.data(), + geometry->n_cell_atoms, this->neigh.size(), this->neigh.data(), spatial_gradient.data(), size, + geometry->cell_size.data(), this->Ms ); + CU_CHECK_AND_SYNC(); +} + +__global__ void CU_Gradient_DMI1( + const Vector3 * spins, Vector3 * gradient, Matrix3 * spatial_gradient, size_t size, Matrix3 dmi_tensor, scalar Ms ) +{ + for( auto icell = blockIdx.x * blockDim.x + threadIdx.x; icell < size; icell += blockDim.x * gridDim.x ) + { + for( unsigned int i = 0; i < 3; ++i ) + { + gradient[icell][0] -= 4 * C::mu_B + * ( dmi_tensor( 1, i ) * spatial_gradient[icell]( 2, i ) + - 2 * dmi_tensor( 2, i ) * spatial_gradient[icell]( 1, i ) ) + / Ms; + gradient[icell][1] -= 4 * C::mu_B + * ( dmi_tensor( 2, i ) * spatial_gradient[icell]( 0, i ) + - 2 * dmi_tensor( 0, i ) * spatial_gradient[icell]( 2, i ) ) + / Ms; + gradient[icell][2] -= 4 * C::mu_B + * ( dmi_tensor( 0, i ) * spatial_gradient[icell]( 1, i ) + - 2 * dmi_tensor( 1, i ) * spatial_gradient[icell]( 0, i ) ) + / Ms; + } + } +} +void Hamiltonian_Micromagnetic::Gradient_DMI( const vectorfield & spins, vectorfield & gradient ) +{ + int size = geometry->n_cells_total; + CU_Gradient_DMI1<<<( size + 1023 ) / 1024, 1024>>>( + spins.data(), gradient.data(), spatial_gradient.data(), size, dmi_tensor, this->Ms ); + CU_CHECK_AND_SYNC(); +} + +__global__ void CU_FFT_Pointwise_Mult1( + FFT::FFT_cpx_type * ft_D_matrices, FFT::FFT_cpx_type * ft_spins, FFT::FFT_cpx_type * res_mult, + int * iteration_bounds, int i_b1, int * inter_sublattice_lookup, FFT::StrideContainer dipole_stride, + FFT::StrideContainer spin_stride, const scalar Ms ) +{ + int n = iteration_bounds[0] * iteration_bounds[1] * iteration_bounds[2] * iteration_bounds[3]; + int tupel[4]; + int idx_b1, idx_b2, idx_d; + + for( int ispin = blockIdx.x * blockDim.x + threadIdx.x; ispin < n; ispin += blockDim.x * gridDim.x ) + { + cu_tupel_from_idx( ispin, tupel, iteration_bounds, 4 ); // tupel now is {i_b2, a, b, c} + + int & b_inter = inter_sublattice_lookup[i_b1 + tupel[0] * iteration_bounds[0]]; + + idx_b1 + = i_b1 * spin_stride.basis + tupel[1] * spin_stride.a + tupel[2] * spin_stride.b + tupel[3] * spin_stride.c; + idx_b2 = tupel[0] * spin_stride.basis + tupel[1] * spin_stride.a + tupel[2] * spin_stride.b + + tupel[3] * spin_stride.c; + idx_d = b_inter * dipole_stride.basis + tupel[1] * dipole_stride.a + tupel[2] * dipole_stride.b + + tupel[3] * dipole_stride.c; + + auto & fs_x = ft_spins[idx_b2]; + auto & fs_y = ft_spins[idx_b2 + 1 * spin_stride.comp]; + auto & fs_z = ft_spins[idx_b2 + 2 * spin_stride.comp]; + + auto & fD_xx = ft_D_matrices[idx_d]; + auto & fD_xy = ft_D_matrices[idx_d + 1 * dipole_stride.comp]; + auto & fD_xz = ft_D_matrices[idx_d + 2 * dipole_stride.comp]; + auto & fD_yy = ft_D_matrices[idx_d + 3 * dipole_stride.comp]; + auto & fD_yz = ft_D_matrices[idx_d + 4 * dipole_stride.comp]; + auto & fD_zz = ft_D_matrices[idx_d + 5 * dipole_stride.comp]; + + if( tupel[0] == 0 ) + { + res_mult[idx_b1].x = FFT::mult3D( fD_xx, fD_xy, fD_xz, fs_x, fs_y, fs_z ).x; + res_mult[idx_b1].y = FFT::mult3D( fD_xx, fD_xy, fD_xz, fs_x, fs_y, fs_z ).y; + res_mult[idx_b1 + 1 * spin_stride.comp].x = FFT::mult3D( fD_xy, fD_yy, fD_yz, fs_x, fs_y, fs_z ).x; + res_mult[idx_b1 + 1 * spin_stride.comp].y = FFT::mult3D( fD_xy, fD_yy, fD_yz, fs_x, fs_y, fs_z ).y; + res_mult[idx_b1 + 2 * spin_stride.comp].x = FFT::mult3D( fD_xz, fD_yz, fD_zz, fs_x, fs_y, fs_z ).x; + res_mult[idx_b1 + 2 * spin_stride.comp].y = FFT::mult3D( fD_xz, fD_yz, fD_zz, fs_x, fs_y, fs_z ).y; + } + else + { + atomicAdd( &res_mult[idx_b1].x, FFT::mult3D( fD_xx, fD_xy, fD_xz, fs_x, fs_y, fs_z ).x ); + atomicAdd( &res_mult[idx_b1].y, FFT::mult3D( fD_xx, fD_xy, fD_xz, fs_x, fs_y, fs_z ).y ); + atomicAdd( + &res_mult[idx_b1 + 1 * spin_stride.comp].x, FFT::mult3D( fD_xy, fD_yy, fD_yz, fs_x, fs_y, fs_z ).x ); + atomicAdd( + &res_mult[idx_b1 + 1 * spin_stride.comp].y, FFT::mult3D( fD_xy, fD_yy, fD_yz, fs_x, fs_y, fs_z ).y ); + atomicAdd( + &res_mult[idx_b1 + 2 * spin_stride.comp].x, FFT::mult3D( fD_xz, fD_yz, fD_zz, fs_x, fs_y, fs_z ).x ); + atomicAdd( + &res_mult[idx_b1 + 2 * spin_stride.comp].y, FFT::mult3D( fD_xz, fD_yz, fD_zz, fs_x, fs_y, fs_z ).y ); + } + } +} + +__global__ void CU_Write_FFT_Gradients1( + const FFT::FFT_real_type * resiFFT, Vector3 * gradient, FFT::StrideContainer spin_stride, int * iteration_bounds, + int n_cell_atoms, scalar * mu_s, int sublattice_size, const scalar Ms ) +{ + int nos = iteration_bounds[0] * iteration_bounds[1] * iteration_bounds[2] * iteration_bounds[3]; + int tupel[4]; + int idx_pad; + for( int idx_orig = blockIdx.x * blockDim.x + threadIdx.x; idx_orig < nos; idx_orig += blockDim.x * gridDim.x ) + { + + cu_tupel_from_idx( idx_orig, tupel, iteration_bounds, 4 ); // tupel now is {ib, a, b, c} + idx_pad = tupel[0] * spin_stride.basis + tupel[1] * spin_stride.a + tupel[2] * spin_stride.b + + tupel[3] * spin_stride.c; + // printf("%d %f %f\n", idx_orig, resiFFT[idx_pad],gradient[idx_orig][0]); + gradient[idx_orig][0] -= C::mu_B * resiFFT[idx_pad] * Ms * 1e-7 / ( sublattice_size ); + gradient[idx_orig][1] -= C::mu_B * resiFFT[idx_pad + 1 * spin_stride.comp] * Ms * 1e-7 / ( sublattice_size ); + gradient[idx_orig][2] -= C::mu_B * resiFFT[idx_pad + 2 * spin_stride.comp] * Ms * 1e-7 / ( sublattice_size ); + } +} + +void Hamiltonian_Micromagnetic::Gradient_DDI( const vectorfield & spins, vectorfield & gradient ) +{ + // this->Gradient_DDI_Direct(spins, gradient); + this->Gradient_DDI_FFT( spins, gradient ); + /* + if (this->ddi_method == DDI_Method::FFT) + { + printf("sasas"); + this->Gradient_DDI_FFT(spins, gradient); + } + else if (this->ddi_method == DDI_Method::Cutoff) + { + // TODO: Merge these implementations in the future + if (this->ddi_cutoff_radius >= 0) + this->Gradient_DDI_Cutoff(spins, gradient); + else + this->Gradient_DDI_Direct(spins, gradient); + } +*/ +} +void Hamiltonian_Micromagnetic::Gradient_DDI_Cutoff( const vectorfield & spins, vectorfield & gradient ) +{ + // TODO +} +void Hamiltonian_Micromagnetic::Gradient_DDI_FFT( const vectorfield & spins, vectorfield & gradient ) +{ + auto & ft_D_matrices = transformed_dipole_matrices; + + auto & ft_spins = fft_plan_spins.cpx_ptr; + + auto & res_iFFT = fft_plan_reverse.real_ptr; + auto & res_mult = fft_plan_reverse.cpx_ptr; + + int number_of_mults = it_bounds_pointwise_mult[0] * it_bounds_pointwise_mult[1] * it_bounds_pointwise_mult[2] + * it_bounds_pointwise_mult[3]; + + FFT_Spins( spins ); + + // TODO: also parallelize over i_b1 + // Loop over basis atoms (i.e sublattices) and add contribution of each sublattice + for( int i_b1 = 0; i_b1 < geometry->n_cell_atoms; ++i_b1 ) + CU_FFT_Pointwise_Mult1<<<( number_of_mults + 1023 ) / 1024, 1024>>>( + ft_D_matrices.data(), ft_spins.data(), res_mult.data(), it_bounds_pointwise_mult.data(), i_b1, + inter_sublattice_lookup.data(), dipole_stride, spin_stride, Ms ); + CU_CHECK_AND_SYNC(); + FFT::batch_iFour_3D( fft_plan_reverse ); + // scalar * delta = geometry->cell_size.data(); + int sublattice_size = it_bounds_write_dipole[0] * it_bounds_write_dipole[1] * it_bounds_write_dipole[2]; + CU_Write_FFT_Gradients1<<<( geometry->nos + 1023 ) / 1024, 1024>>>( + res_iFFT.data(), gradient.data(), spin_stride, it_bounds_write_gradients.data(), geometry->n_cell_atoms, + geometry->mu_s.data(), sublattice_size, Ms ); + CU_CHECK_AND_SYNC(); +} // end Field_DipoleDipole + +void Hamiltonian_Micromagnetic::Gradient_DDI_Direct( const vectorfield & spins, vectorfield & gradient ) +{ + int tupel1[3]; + int tupel2[3]; + int sublattice_size = it_bounds_write_dipole[0] * it_bounds_write_dipole[1] * it_bounds_write_dipole[2]; + // prefactor of ddi interaction + // scalar mult = 2.0133545*1e-28 * 0.057883817555 * 0.057883817555 / (4 + // * 3.141592653589793238462643383279502884197169399375105820974 * 1e-30); + scalar mult = 1 / ( 4 * 3.141592653589793238462643383279502884197169399375105820974 ); + scalar m0 = ( 4 * 3.141592653589793238462643383279502884197169399375105820974 ) * 1e-7; + int img_a = boundary_conditions[0] == 0 ? 0 : ddi_n_periodic_images[0]; + int img_b = boundary_conditions[1] == 0 ? 0 : ddi_n_periodic_images[1]; + int img_c = boundary_conditions[2] == 0 ? 0 : ddi_n_periodic_images[2]; + scalar * delta = geometry->cell_size.data(); + for( int idx1 = 0; idx1 < geometry->nos; idx1++ ) + { + double kk = 0; + for( int idx2 = 0; idx2 < geometry->nos; idx2++ ) + { + int a1 = idx1 % ( it_bounds_write_spins[1] ); + int b1 = ( (int)( idx1 / it_bounds_write_spins[1] ) ) % ( it_bounds_write_spins[2] ); + int c1 = (int)idx1 / ( it_bounds_write_spins[1] * it_bounds_write_spins[2] ); + int a2 = idx2 % ( it_bounds_write_spins[1] ); + int b2 = ( (int)( idx2 / it_bounds_write_spins[1] ) ) % ( it_bounds_write_spins[2] ); + int c2 = (int)idx2 / ( it_bounds_write_spins[1] * it_bounds_write_spins[2] ); + /*int a_idx = a < n_cells[0] ? a : a - iteration_bounds[0]; + int b_idx = b < n_cells[1] ? b : b - iteration_bounds[1]; + int c_idx = c < n_cells[2] ? c : c - iteration_bounds[2];*/ + int a_idx = a1 - a2; + int b_idx = b1 - b2; + int c_idx = c1 - c2; + if( ( a_idx == 0 ) && ( b_idx == 0 ) && ( c_idx == 0 ) ) + { + // printf("test\n"); + // continue; + } + // printf("%d %d %d\n", a_idx,b_idx,c_idx); + /*if ((a_idx==20) || (b_idx==20) || (c_idx==1)){ + //printf("test1\n"); + //if (c_idx!=1) + //printf("%d %d %d %d\n", a_idx, b_idx, c_idx, dipole_stride.comp); + continue; + }*/ + // scalar delta[3] = { 3,3,0.3 }; + // int idx = b_inter * dipole_stride.basis + a * dipole_stride.a + b * dipole_stride.b + c * dipole_stride.c; + scalar Dxx = 0, Dxy = 0, Dxz = 0, Dyy = 0, Dyz = 0, Dzz = 0; + + Vector3 cell_sizes = { geometry->lattice_constant * geometry->bravais_vectors[0].norm(), + geometry->lattice_constant * geometry->bravais_vectors[1].norm(), + geometry->lattice_constant * geometry->bravais_vectors[2].norm() }; + // asa + for( int i = 0; i < 2; i++ ) + { + for( int j = 0; j < 2; j++ ) + { + for( int k = 0; k < 2; k++ ) + { + double r = sqrt( + ( a_idx + i - 0.5f ) * ( a_idx + i - 0.5f ) * cell_sizes[0] * cell_sizes[0] + + ( b_idx + j - 0.5f ) * ( b_idx + j - 0.5f ) * cell_sizes[1] * cell_sizes[1] + + ( c_idx + k - 0.5f ) * ( c_idx + k - 0.5f ) * cell_sizes[2] * cell_sizes[2] ); + Dxx += mult * pow( -1.0f, i + j + k ) + * atan( + ( ( c_idx + k - 0.5f ) * ( b_idx + j - 0.5f ) * cell_sizes[1] * cell_sizes[2] + / cell_sizes[0] / r / ( a_idx + i - 0.5f ) ) ); + // fft_dipole_inputs[idx + 1 * dipole_stride.comp] += -mult * pow(-1.0f, i + j + k) * + // log(abs(((c_idx + k - 0.5f)* cell_sizes[2] + r)/((c_idx + k - 0.5f)* cell_sizes[2] - r))); + // fft_dipole_inputs[idx + 2 * dipole_stride.comp] += -mult * pow(-1.0f, i + j + k) * + // log(abs(((b_idx + j - 0.5f)* cell_sizes[1] + r)/((b_idx + j - 0.5f)* cell_sizes[1] - r))); + Dxy -= mult * pow( -1.0f, i + j + k ) * log( ( ( ( c_idx + k - 0.5f ) * cell_sizes[2] + r ) ) ); + Dxz -= mult * pow( -1.0f, i + j + k ) * log( ( ( ( b_idx + j - 0.5f ) * cell_sizes[1] + r ) ) ); + + Dyy += mult * pow( -1.0f, i + j + k ) + * atan( + ( ( a_idx + i - 0.5f ) * ( c_idx + k - 0.5f ) * cell_sizes[2] * cell_sizes[0] + / cell_sizes[1] / r / ( b_idx + j - 0.5f ) ) ); + // fft_dipole_inputs[idx + 4 * dipole_stride.comp] += -mult * pow(-1.0f, i + j + k) * + // log(abs(((a_idx + i - 0.5f)* cell_sizes[0] + r)/((a_idx + i - 0.5f)* cell_sizes[0] - r))); + Dyz -= mult * pow( -1.0f, i + j + k ) * log( ( ( ( a_idx + i - 0.5f ) * cell_sizes[0] + r ) ) ); + Dzz += mult * pow( -1.0f, i + j + k ) + * atan( + ( ( b_idx + j - 0.5f ) * ( a_idx + i - 0.5f ) * cell_sizes[0] * cell_sizes[1] + / cell_sizes[2] / r / ( c_idx + k - 0.5f ) ) ); + } + } + } /* + Dxx=Nii(a_idx*delta[0],b_idx*delta[1],c_idx*delta[2],delta[0],delta[1],delta[2]); + Dxy=Nij(a_idx*delta[0],b_idx*delta[1],c_idx*delta[2],delta[0],delta[1],delta[2]); + Dxz=Nij(a_idx*delta[0],c_idx*delta[2], b_idx*delta[1],delta[0],delta[2],delta[1]); + Dyy=Nii(b_idx*delta[1],a_idx*delta[0],c_idx*delta[2],delta[1],delta[0],delta[2]); + Dyz=Nij(b_idx*delta[1],c_idx*delta[2], b_idx*delta[1],delta[1],delta[2],delta[0]); + Dzz=Nii(c_idx*delta[2],a_idx*delta[0],b_idx*delta[1],delta[2],delta[0],delta[1]);*/ + if( idx1 == 42 ) + { + if( ( a_idx == 0 ) && ( b_idx == 0 ) && ( c_idx == 0 ) ) + { + printf( "000 Dxx=%f Dxy=%f Dxz=%f Dyy=%f Dyz=%f Dzz=%f\n", Dxx, Dxy, Dxz, Dyy, Dyz, Dzz ); + } + if( ( a_idx == 1 ) && ( b_idx == 0 ) && ( c_idx == 0 ) ) + { + printf( "100 Dxx=%f Dxy=%f Dxz=%f Dyy=%f Dyz=%f Dzz=%f\n", Dxx, Dxy, Dxz, Dyy, Dyz, Dzz ); + } + if( ( a_idx == 0 ) && ( b_idx == 1 ) && ( c_idx == 0 ) ) + { + printf( "010 Dxx=%f Dxy=%f Dxz=%f Dyy=%f Dyz=%f Dzz=%f\n", Dxx, Dxy, Dxz, Dyy, Dyz, Dzz ); + } + if( ( a_idx == -1 ) && ( b_idx == 1 ) && ( c_idx == 0 ) ) + { + printf( "-110 Dxx=%f Dxy=%f Dxz=%f Dyy=%f Dyz=%f Dzz=%f\n", Dxx, Dxy, Dxz, Dyy, Dyz, Dzz ); + } + if( ( a_idx == 1 ) && ( b_idx == 1 ) && ( c_idx == 0 ) ) + { + printf( "110 Dxx=%f Dxy=%f Dxz=%f Dyy=%f Dyz=%f Dzz=%f\n", Dxx, Dxy, Dxz, Dyy, Dyz, Dzz ); + } + if( ( a_idx == 2 ) && ( b_idx == 0 ) && ( c_idx == 0 ) ) + { + printf( "200 Dxx=%f Dxy=%f Dxz=%f Dyy=%f Dyz=%f Dzz=%f\n", Dxx, Dxy, Dxz, Dyy, Dyz, Dzz ); + } + if( ( a_idx == 0 ) && ( b_idx == 2 ) && ( c_idx == 0 ) ) + { + printf( "020 Dxx=%f Dxy=%f Dxz=%f Dyy=%f Dyz=%f Dzz=%f\n", Dxx, Dxy, Dxz, Dyy, Dyz, Dzz ); + } + if( ( a_idx == 2 ) && ( b_idx == 2 ) && ( c_idx == 0 ) ) + { + printf( "220 Dxx=%f Dxy=%f Dxz=%f Dyy=%f Dyz=%f Dzz=%f\n", Dxx, Dxy, Dxz, Dyy, Dyz, Dzz ); + } + if( ( a_idx == 2 ) && ( b_idx == -2 ) && ( c_idx == 0 ) ) + { + printf( "2-20 Dxx=%f Dxy=%f Dxz=%f Dyy=%f Dyz=%f Dzz=%f\n", Dxx, Dxy, Dxz, Dyy, Dyz, Dzz ); + } + // printf("x=%f y=%f z=%f\n",spins[idx2][0],spins[idx2][1],spins[idx2][2]); + } + kk += Dxx; + gradient[idx1][0] -= ( Dxx * spins[idx2][0] + Dxy * spins[idx2][1] + Dxz * spins[idx2][2] ) * Ms * m0; + gradient[idx1][1] -= ( Dxy * spins[idx2][0] + Dyy * spins[idx2][1] + Dyz * spins[idx2][2] ) * Ms * m0; + gradient[idx1][2] -= ( Dxz * spins[idx2][0] + Dyz * spins[idx2][1] + Dzz * spins[idx2][2] ) * Ms * m0; + } + if( idx1 == 30 ) + { + // printf("x=%f y=%f z=%f\n",spins[idx1][0],spins[idx1][1],spins[idx1][2]); + // printf("kk=%f gx=%f gy=%f gz=%f\n",kk, gradient[idx1][0]/8e5/m0,gradient[idx1][1],gradient[idx1][2]); + } + } +} +__global__ void CU_Write_FFT_Spin_Input1( + FFT::FFT_real_type * fft_spin_inputs, const Vector3 * spins, int * iteration_bounds, + FFT::StrideContainer spin_stride, scalar * mu_s ) +{ + int nos = iteration_bounds[0] * iteration_bounds[1] * iteration_bounds[2] * iteration_bounds[3]; + int tupel[4]; + int idx_pad; + for( int idx_orig = blockIdx.x * blockDim.x + threadIdx.x; idx_orig < nos; idx_orig += blockDim.x * gridDim.x ) + { + cu_tupel_from_idx( idx_orig, tupel, iteration_bounds, 4 ); // tupel now is {ib, a, b, c} + idx_pad = tupel[0] * spin_stride.basis + tupel[1] * spin_stride.a + tupel[2] * spin_stride.b + + tupel[3] * spin_stride.c; + fft_spin_inputs[idx_pad] = spins[idx_orig][0]; + fft_spin_inputs[idx_pad + 1 * spin_stride.comp] = spins[idx_orig][1]; + fft_spin_inputs[idx_pad + 2 * spin_stride.comp] = spins[idx_orig][2]; + // printf("%f %f\n",fft_spin_inputs[idx_pad], fft_spin_inputs[idx_pad+30]); + } +} + +void Hamiltonian_Micromagnetic::FFT_Spins( const vectorfield & spins ) +{ + CU_Write_FFT_Spin_Input1<<<( geometry->nos + 1023 ) / 1024, 1024>>>( + fft_plan_spins.real_ptr.data(), spins.data(), it_bounds_write_spins.data(), spin_stride, + geometry->mu_s.data() ); + CU_CHECK_AND_SYNC(); + FFT::batch_Four_3D( fft_plan_spins ); +} +__global__ void CU_Write_FFT_Dipole_Input1( + FFT::FFT_real_type * fft_dipole_inputs, int * iteration_bounds, const Vector3 * translation_vectors, + int n_cell_atoms, Vector3 * cell_atom_translations, int * n_cells, int * inter_sublattice_lookup, int * img, + FFT::StrideContainer dipole_stride, const Vector3 cell_lengths ) +{ + int tupel[3]; + int sublattice_size = iteration_bounds[0] * iteration_bounds[1] * iteration_bounds[2]; + // prefactor of ddi interaction + // scalar mult = 2.0133545*1e-28 * 0.057883817555 * 0.057883817555 / (4 + // * 3.141592653589793238462643383279502884197169399375105820974 * 1e-30); scalar mult = 1 / (4 + // * 3.141592653589793238462643383279502884197169399375105820974); + scalar mult = 1; + for( int i = blockIdx.x * blockDim.x + threadIdx.x; i < sublattice_size; i += blockDim.x * gridDim.x ) + { + cu_tupel_from_idx( i, tupel, iteration_bounds, 3 ); // tupel now is {a, b, c} + auto & a = tupel[0]; + auto & b = tupel[1]; + auto & c = tupel[2]; + /*if ((a>198)||(b>198)||(c>198)){ + printf("%d %d %d\n", a,b,c); + }*/ + /*int a_idx = a < n_cells[0] ? a : a - iteration_bounds[0]; + int b_idx = b < n_cells[1] ? b : b - iteration_bounds[1]; + int c_idx = c < n_cells[2] ? c : c - iteration_bounds[2];*/ + /*int a_idx = a +1 - (int)iteration_bounds[0]/2; + int b_idx = b +1- (int)iteration_bounds[1]/2; + int c_idx = c +1- (int)iteration_bounds[2]/2;*/ + int a_idx = a < n_cells[0] ? a : a - iteration_bounds[0]; + int b_idx = b < n_cells[1] ? b : b - iteration_bounds[1]; + int c_idx = c < n_cells[2] ? c : c - iteration_bounds[2]; + + int idx = a * dipole_stride.a + b * dipole_stride.b + c * dipole_stride.c; + + for( int i = 0; i < 2; i++ ) + { + for( int j = 0; j < 2; j++ ) + { + for( int k = 0; k < 2; k++ ) + { + double r = sqrt( + ( a_idx + i - 0.5f ) * ( a_idx + i - 0.5f ) * cell_lengths[0] * cell_lengths[0] + + ( b_idx + j - 0.5f ) * ( b_idx + j - 0.5f ) * cell_lengths[1] * cell_lengths[1] + + ( c_idx + k - 0.5f ) * ( c_idx + k - 0.5f ) * cell_lengths[2] * cell_lengths[2] ); + fft_dipole_inputs[idx] += mult * pow( -1.0f, i + j + k ) + * atan( + ( ( c_idx + k - 0.5f ) * ( b_idx + j - 0.5f ) * cell_lengths[1] + * cell_lengths[2] / cell_lengths[0] / r / ( a_idx + i - 0.5f ) ) ); + // fft_dipole_inputs[idx + 1 * dipole_stride.comp] += -mult * pow(-1.0f, i + j + k) * + // log(abs(((c_idx + k - 0.5f)* cell_lengths[2] + r)/((c_idx + k - 0.5f)* cell_lengths[2] - r))); + // fft_dipole_inputs[idx + 2 * dipole_stride.comp] += -mult * pow(-1.0f, i + j + k) * + // log(abs(((b_idx + j - 0.5f)* cell_lengths[1] + r)/((b_idx + j - 0.5f)* cell_lengths[1] - r))); + fft_dipole_inputs[idx + 1 * dipole_stride.comp] + -= mult * pow( -1.0f, i + j + k ) * log( ( ( ( c_idx + k - 0.5f ) * cell_lengths[2] + r ) ) ); + fft_dipole_inputs[idx + 2 * dipole_stride.comp] + -= mult * pow( -1.0f, i + j + k ) * log( ( ( ( b_idx + j - 0.5f ) * cell_lengths[1] + r ) ) ); + + fft_dipole_inputs[idx + 3 * dipole_stride.comp] + += mult * pow( -1.0f, i + j + k ) + * atan( + ( ( a_idx + i - 0.5f ) * ( c_idx + k - 0.5f ) * cell_lengths[2] * cell_lengths[0] + / cell_lengths[1] / r / ( b_idx + j - 0.5f ) ) ); + // fft_dipole_inputs[idx + 4 * dipole_stride.comp] += -mult * pow(-1.0f, i + j + k) * + // log(abs(((a_idx + i - 0.5f)* cell_lengths[0] + r)/((a_idx + i - 0.5f)* cell_lengths[0] - r))); + fft_dipole_inputs[idx + 4 * dipole_stride.comp] + -= mult * pow( -1.0f, i + j + k ) * log( ( ( ( a_idx + i - 0.5f ) * cell_lengths[0] + r ) ) ); + fft_dipole_inputs[idx + 5 * dipole_stride.comp] + += mult * pow( -1.0f, i + j + k ) + * atan( + ( ( b_idx + j - 0.5f ) * ( a_idx + i - 0.5f ) * cell_lengths[0] * cell_lengths[1] + / cell_lengths[2] / r / ( c_idx + k - 0.5f ) ) ); + } + } + } + + // if (fft_dipole_inputs[idx]<-0.03) + } +} + +/* +void Hamiltonian_Micromagnetic::FFT_Dipole_Matrices( FFT::FFT_Plan & fft_plan_dipole, int img_a, int img_b, int img_c ) +{ + auto & fft_dipole_inputs = fft_plan_dipole.real_ptr; + + field img = { img_a, img_b, img_c }; + + // Work around to make bravais vectors and cell_atoms available to GPU as they are currently saves as std::vectors + // and not fields ... + auto translation_vectors = field(); + auto cell_atom_translations = field(); + + for( int i = 0; i < 3; i++ ) + translation_vectors.push_back( geometry->lattice_constant * geometry->bravais_vectors[i] ); + + for( int i = 0; i < geometry->n_cell_atoms; i++ ) + cell_atom_translations.push_back( geometry->positions[i] ); + + Vector3 cell_sizes = { geometry->lattice_constant * geometry->bravais_vectors[0].norm(), + geometry->lattice_constant * geometry->bravais_vectors[1].norm(), + geometry->lattice_constant * geometry->bravais_vectors[2].norm() }; + + CU_Write_FFT_Dipole_Input1<<<( sublattice_size + 1023 ) / 1024, 1024>>>( + fft_dipole_inputs.data(), it_bounds_write_dipole.data(), translation_vectors.data(), geometry->n_cell_atoms, + cell_atom_translations.data(), geometry->n_cells.data(), inter_sublattice_lookup.data(), img.data(), + dipole_stride, cell_sizes ); + CU_CHECK_AND_SYNC(); + FFT::batch_Four_3D( fft_plan_dipole ); +} +*/ + +void Hamiltonian_Micromagnetic::Prepare_DDI() +{ + Clean_DDI(); + + n_cells_padded.resize( 3 ); + n_cells_padded[0] = ( geometry->n_cells[0] > 1 ) ? 2 * geometry->n_cells[0] : 1; + n_cells_padded[1] = ( geometry->n_cells[1] > 1 ) ? 2 * geometry->n_cells[1] : 1; + n_cells_padded[2] = ( geometry->n_cells[2] > 1 ) ? 2 * geometry->n_cells[2] : 1; + sublattice_size = n_cells_padded[0] * n_cells_padded[1] * n_cells_padded[2]; + // printf("111 %d %d %d\n", n_cells_padded[0],n_cells_padded[1],n_cells_padded[2]); + + inter_sublattice_lookup.resize( geometry->n_cell_atoms * geometry->n_cell_atoms ); + + // we dont need to transform over length 1 dims + std::vector fft_dims; + for( int i = 2; i >= 0; i-- ) // notice that reverse order is important! + { + if( n_cells_padded[i] > 1 ) + fft_dims.push_back( n_cells_padded[i] ); + } + + // Count how many distinct inter-lattice contributions we need to store + n_inter_sublattice = 0; + for( int i = 0; i < geometry->n_cell_atoms; i++ ) + { + for( int j = 0; j < geometry->n_cell_atoms; j++ ) + { + if( i != 0 && i == j ) + continue; + n_inter_sublattice++; + } + } + printf( "lex%d %d %d\n", n_inter_sublattice, fft_dims[0], fft_dims[1] ); + // Set the iteration bounds for the nested for loops that are flattened in the kernels + it_bounds_write_spins + = { geometry->n_cell_atoms, geometry->n_cells[0], geometry->n_cells[1], geometry->n_cells[2] }; + + it_bounds_write_dipole = { n_cells_padded[0], n_cells_padded[1], n_cells_padded[2] }; + + it_bounds_pointwise_mult = { geometry->n_cell_atoms, + ( n_cells_padded[0] / 2 + 1 ), // due to redundancy in real fft + n_cells_padded[1], n_cells_padded[2] }; + + it_bounds_write_gradients + = { geometry->n_cell_atoms, geometry->n_cells[0], geometry->n_cells[1], geometry->n_cells[2] }; + + FFT::FFT_Plan fft_plan_dipole = FFT::FFT_Plan( fft_dims, false, 6 * n_inter_sublattice, sublattice_size ); + fft_plan_spins = FFT::FFT_Plan( fft_dims, false, 3 * geometry->n_cell_atoms, sublattice_size ); + fft_plan_reverse = FFT::FFT_Plan( fft_dims, true, 3 * geometry->n_cell_atoms, sublattice_size ); + + field temp_s = { &spin_stride.comp, &spin_stride.basis, &spin_stride.a, &spin_stride.b, &spin_stride.c }; + field temp_d + = { &dipole_stride.comp, &dipole_stride.basis, &dipole_stride.a, &dipole_stride.b, &dipole_stride.c }; + ; + FFT::get_strides( + temp_s, { 3, this->geometry->n_cell_atoms, n_cells_padded[0], n_cells_padded[1], n_cells_padded[2] } ); + FFT::get_strides( temp_d, { 6, n_inter_sublattice, n_cells_padded[0], n_cells_padded[1], n_cells_padded[2] } ); + /* + //perform FFT of dipole matrices + int img_a = boundary_conditions[0] == 0 ? 0 : ddi_n_periodic_images[0]; + int img_b = boundary_conditions[1] == 0 ? 0 : ddi_n_periodic_images[1]; + int img_c = boundary_conditions[2] == 0 ? 0 : ddi_n_periodic_images[2]; + + FFT_Dipole_Matrices(fft_plan_dipole, img_a, img_b, img_c); */ + /* + FFT_Dipole_Matrices( fft_plan_dipole, 0, 0, 0 ); + */ + + transformed_dipole_matrices = std::move( fft_plan_dipole.cpx_ptr ); +} // end prepare + +void Hamiltonian_Micromagnetic::Clean_DDI() +{ + fft_plan_spins = FFT::FFT_Plan(); + fft_plan_reverse = FFT::FFT_Plan(); +} + +void Hamiltonian_Micromagnetic::Hessian( const vectorfield & spins, MatrixX & hessian ) {} + +// Hamiltonian name as string +static const std::string name = "Micromagnetic"; +const std::string & Hamiltonian_Micromagnetic::Name() const +{ + return name; +} + +} // namespace Engine + +#endif \ No newline at end of file diff --git a/core/src/io/Configparser.cpp b/core/src/io/Configparser.cpp index 588749b33..f90a9f1ef 100644 --- a/core/src/io/Configparser.cpp +++ b/core/src/io/Configparser.cpp @@ -2,12 +2,14 @@ #include #include #include +#include #include #include #include #include #include +#include #include #include @@ -800,7 +802,8 @@ std::unique_ptr Parameters_Method_LLG_from_Config( parameter_log.emplace_back( fmt::format( " {:<17} = {}", "maximum walltime", str_max_walltime ) ); parameter_log.emplace_back( fmt::format( " {:<17} = {}", "n_iterations", parameters->n_iterations ) ); parameter_log.emplace_back( fmt::format( " {:<17} = {}", "n_iterations_log", parameters->n_iterations_log ) ); - parameter_log.emplace_back( fmt::format( " {:<17} = {}", "n_iterations_amortize", parameters->n_iterations_amortize ) ); + parameter_log.emplace_back( + fmt::format( " {:<17} = {}", "n_iterations_amortize", parameters->n_iterations_amortize ) ); parameter_log.emplace_back( fmt::format( " {:<17} = \"{}\"", "output_folder", parameters->output_folder ) ); parameter_log.emplace_back( fmt::format( " {:<17} = {}", "output_any", parameters->output_any ) ); parameter_log.emplace_back( fmt::format( " {:<17} = {}", "output_initial", parameters->output_initial ) ); @@ -981,7 +984,8 @@ std::unique_ptr Parameters_Method_MC_from_Config( co parameter_log.emplace_back( fmt::format( " {:<17} = {}", "maximum walltime", str_max_walltime ) ); parameter_log.emplace_back( fmt::format( " {:<17} = {}", "n_iterations", parameters->n_iterations ) ); parameter_log.emplace_back( fmt::format( " {:<17} = {}", "n_iterations_log", parameters->n_iterations_log ) ); - parameter_log.emplace_back( fmt::format( " {:<17} = {}", "n_iterations_amortize", parameters->n_iterations_amortize ) ); + parameter_log.emplace_back( + fmt::format( " {:<17} = {}", "n_iterations_amortize", parameters->n_iterations_amortize ) ); parameter_log.emplace_back( fmt::format( " {:<17} = \"{}\"", "output_folder", parameters->output_folder ) ); parameter_log.emplace_back( fmt::format( " {:<17} = {}", "output_any", parameters->output_any ) ); parameter_log.emplace_back( fmt::format( " {:<17} = {}", "output_initial", parameters->output_initial ) ); @@ -1163,7 +1167,8 @@ std::unique_ptr Parameters_Method_MMF_from_Config( parameter_log.emplace_back( fmt::format( " {:<17} = {}", "maximum walltime", str_max_walltime ) ); parameter_log.emplace_back( fmt::format( " {:<17} = {}", "n_iterations", parameters->n_iterations ) ); parameter_log.emplace_back( fmt::format( " {:<17} = {}", "n_iterations_log", parameters->n_iterations_log ) ); - parameter_log.emplace_back( fmt::format( " {:<17} = {}", "n_iterations_amortize", parameters->n_iterations_amortize ) ); + parameter_log.emplace_back( + fmt::format( " {:<17} = {}", "n_iterations_amortize", parameters->n_iterations_amortize ) ); parameter_log.emplace_back( fmt::format( " {:<17} = \"{}\"", "output_folder", parameters->output_folder ) ); parameter_log.emplace_back( fmt::format( " {:<17} = {}", "output_any", parameters->output_any ) ); parameter_log.emplace_back( fmt::format( " {:<17} = {}", "output_initial", parameters->output_initial ) ); @@ -1227,9 +1232,13 @@ Hamiltonian_from_Config( const std::string & config_file_name, std::shared_ptr Hamiltonian_Heisenberg_from_Conf return hamiltonian; } // end Hamiltonian_Heisenberg_From_Config +std::unique_ptr Hamiltonian_Micromagnetic_from_Config( + const std::string & config_file_name, const std::shared_ptr geometry ) +{ + if( geometry->classifier != Data::BravaisLatticeType::Rectilinear + && geometry->classifier != Data::BravaisLatticeType::SC ) + { + spirit_throw( + Utility::Exception_Classifier::System_not_Initialized, Log_Level::Severe, + fmt::format( + "Hamiltonian: Cannot use micromagnetic Hamiltonian on non-rectilinear geometry (type {})", + int( geometry->classifier ) ) ); + } + //-------------- Insert default values here ----------------------------- + // Boundary conditions (a, b, c) + std::vector boundary_conditions_i = { 0, 0, 0 }; + intfield boundary_conditions = { false, false, false }; + + scalar Ms; + + // The order of the finite difference approximation of the spatial gradient + int spatial_gradient_order = 1; + + // External Magnetic Field + scalar field = 0; + Vector3 field_normal = { 0.0, 0.0, 1.0 }; + + scalar anisotropy_magnitude; + Vector3 anisotropy_normal; + Matrix3 anisotropy_tensor; + scalar exchange_magnitude = 0; + Matrix3 exchange_tensor; + scalar dmi_magnitude = 0; + Matrix3 dmi_tensor; + + // Dipolar + std::string ddi_method_str = "none"; + auto ddi_method = Engine::DDI_Method::None; + intfield ddi_n_periodic_images = { 4, 4, 4 }; + scalar ddi_radius = 0.0; + + //------------------------------- Parser -------------------------------- + Log( Log_Level::Debug, Log_Sender::IO, "Hamiltonian_Micromagnetic: building" ); + try + { + IO::Filter_File_Handle config_file_handle( config_file_name ); + try + { + // Boundary conditions + config_file_handle.Read_3Vector( boundary_conditions_i, "boundary_conditions" ); + boundary_conditions[0] = static_cast( boundary_conditions_i[0] != 0 ); + boundary_conditions[1] = static_cast( boundary_conditions_i[1] != 0 ); + boundary_conditions[2] = static_cast( boundary_conditions_i[2] != 0 ); + } + catch( ... ) + { + spirit_handle_exception_core( + fmt::format( "Unable to read boundary conditions from config file \"{}\"", config_file_name ) ); + } + + // Precision of the spatial gradient calculation + config_file_handle.Read_Single( spatial_gradient_order, "spatial_gradient_order" ); + + if( config_file_handle.Find( "Ms" ) ) + { + config_file_handle.Read_Single( Ms, "Ms" ); + } + else + { + Log( Log_Level::Warning, Log_Sender::IO, + "Input for 'Ms' has not been found. Inferring from atomistic cell instead." ); + Ms = geometry->getMs(); + } + + // Field + config_file_handle.Read_Single( field, "external_field_magnitude" ); + config_file_handle.Read_Vector3( field_normal, "external_field_normal" ); + field_normal.normalize(); + if( field_normal.norm() < 1e-8 ) + { + field_normal = { 0, 0, 1 }; + Log( Log_Level::Warning, Log_Sender::IO, + "Input for 'external_field_normal' had norm zero and has been set to (0,0,1)" ); + } + + // TODO: anisotropy + if( config_file_handle.Find( "tensor_anisotropy" ) ) + { + for( int dim = 0; dim < 3; ++dim ) + { + config_file_handle.GetLine(); + config_file_handle >> anisotropy_tensor( dim, 0 ) >> anisotropy_tensor( dim, 1 ) + >> anisotropy_tensor( dim, 2 ); + } + } + else + { + // Read parameters from config + config_file_handle.Read_Single( anisotropy_magnitude, "anisotropy_magnitude" ); + config_file_handle.Read_Vector3( anisotropy_normal, "anisotropy_normal" ); + anisotropy_normal.normalize(); + auto & Kn = anisotropy_normal; + anisotropy_tensor << Kn[0] * Kn[0], Kn[0] * Kn[1], Kn[0] * Kn[2], Kn[1] * Kn[0], Kn[1] * Kn[1], + Kn[1] * Kn[2], Kn[2] * Kn[0], Kn[2] * Kn[1], Kn[2] * Kn[2]; + anisotropy_tensor *= anisotropy_magnitude; + } + + // TODO: exchange + if( config_file_handle.Find( "tensor_exchange" ) ) + { + for( int dim = 0; dim < 3; ++dim ) + { + config_file_handle.GetLine(); + config_file_handle >> exchange_tensor( dim, 0 ) >> exchange_tensor( dim, 1 ) + >> exchange_tensor( dim, 2 ); + } + } + else + { + config_file_handle.Read_Single( exchange_magnitude, "exchange" ); + exchange_tensor << exchange_magnitude, 0, 0, 0, exchange_magnitude, 0, 0, 0, exchange_magnitude; + } + + // TODO: dmi + if( config_file_handle.Find( "tensor_dmi" ) ) + { + for( int dim = 0; dim < 3; ++dim ) + { + config_file_handle.GetLine(); + config_file_handle >> dmi_tensor( dim, 0 ) >> dmi_tensor( dim, 1 ) >> dmi_tensor( dim, 2 ); + } + } + else + { + config_file_handle.Read_Single( dmi_magnitude, "dmi" ); + // dmi_tensor << 0, dmi_magnitude/std::sqrt(3), 0, + // -dmi_magnitude/std::sqrt(3), 0,0, + // 0, 0, 0; + dmi_tensor << dmi_magnitude / std::sqrt( 3 ), 0, 0, 0, dmi_magnitude / std::sqrt( 3 ), 0, 0, 0, + dmi_magnitude / std::sqrt( 3 ); + // dmi_tensor << 0, dmi_magnitude, -dmi_magnitude, + // -dmi_magnitude, 0, dmi_magnitude, + // dmi_magnitude, -dmi_magnitude, 0; + } + + try + { + IO::Filter_File_Handle myfile( config_file_name ); + + // DDI method + myfile.Read_String( ddi_method_str, "ddi_method" ); + if( ddi_method_str == "none" ) + ddi_method = Engine::DDI_Method::None; + else if( ddi_method_str == "fft" ) + ddi_method = Engine::DDI_Method::FFT; + else if( ddi_method_str == "fmm" ) + ddi_method = Engine::DDI_Method::FMM; + else if( ddi_method_str == "cutoff" ) + ddi_method = Engine::DDI_Method::Cutoff; + else + { + Log( Log_Level::Warning, Log_Sender::IO, + fmt::format( + "Hamiltonian_Heisenberg: Keyword 'ddi_method' got passed invalid method \"{}\". Setting to " + "\"none\".", + ddi_method_str ) ); + ddi_method_str = "none"; + } + + // Number of periodical images + myfile.Read_3Vector( ddi_n_periodic_images, "ddi_n_periodic_images" ); + // myfile.Read_Single(ddi_n_periodic_images, "ddi_n_periodic_images"); + + // Dipole-dipole cutoff radius + myfile.Read_Single( ddi_radius, "ddi_radius" ); + } // end try + catch( ... ) + { + spirit_handle_exception_core( + fmt::format( "Unable to read DDI radius from config file \"{}\"", config_file_name ) ); + } + + } // end try + catch( ... ) + { + spirit_handle_exception_core( fmt::format( + "Unable to parse all parameters of the Micromagnetic Hamiltonian from \"{}\"", config_file_name ) ); + } + // Return + std::vector parameter_log; + parameter_log.emplace_back( "Hamiltonian_Micromagnetic:" ); + parameter_log.emplace_back( fmt::format( " {:<24} = {}", "discretisation order", spatial_gradient_order ) ); + parameter_log.emplace_back( fmt::format( " {:<24} = {}", "Ms [A/m]", Ms ) ); + parameter_log.emplace_back( fmt::format( + " {:<24} = {} {} {}", "boundary conditions", boundary_conditions[0], boundary_conditions[1], + boundary_conditions[2] ) ); + parameter_log.emplace_back( fmt::format( " {:<24} = {}", "external field", field ) ); + parameter_log.emplace_back( fmt::format( " {:<24} = {}", "field normal", field_normal.transpose() ) ); + parameter_log.emplace_back( fmt::format( " {:<24} = {}", "anisotropy tensor", anisotropy_tensor.row( 0 ) ) ); + parameter_log.emplace_back( fmt::format( " {:<24} {}", " ", anisotropy_tensor.row( 1 ) ) ); + parameter_log.emplace_back( fmt::format( " {:<24} {}", " ", anisotropy_tensor.row( 2 ) ) ); + parameter_log.emplace_back( fmt::format( " {:<24} = {}", "exchange tensor", exchange_tensor.row( 0 ) ) ); + parameter_log.emplace_back( fmt::format( " {:<24} {}", " ", exchange_tensor.row( 1 ) ) ); + parameter_log.emplace_back( fmt::format( " {:<24} {}", " ", exchange_tensor.row( 2 ) ) ); + parameter_log.emplace_back( fmt::format( " {:<24} = {}", "dmi tensor", dmi_tensor.row( 0 ) ) ); + parameter_log.emplace_back( fmt::format( " {:<24} {}", " ", dmi_tensor.row( 1 ) ) ); + parameter_log.emplace_back( fmt::format( " {:<24} {}", " ", dmi_tensor.row( 2 ) ) ); + parameter_log.emplace_back( fmt::format( " {:<21} = {}", "ddi_method", ddi_method_str ) ); + parameter_log.emplace_back( fmt::format( + " {:<21} = ({} {} {})", "ddi_n_periodic_images", ddi_n_periodic_images[0], ddi_n_periodic_images[1], + ddi_n_periodic_images[2] ) ); + parameter_log.emplace_back( fmt::format( " {:<21} = {}", "ddi_radius", ddi_radius ) ); + Log.SendBlock( Log_Level::Parameter, Log_Sender::IO, parameter_log ); + + auto hamiltonian = std::make_unique( + Ms, field, field_normal, anisotropy_tensor, exchange_tensor, dmi_tensor, ddi_method, ddi_n_periodic_images, + ddi_radius, geometry, spatial_gradient_order, boundary_conditions ); + + Log( Log_Level::Debug, Log_Sender::IO, "Hamiltonian_Micromagnetic: built" ); + return hamiltonian; + +} // end Hamiltonian_Micromagnetic_from_Config + std::unique_ptr Hamiltonian_Gaussian_from_Config( const std::string & config_file_name, std::shared_ptr geometry ) { diff --git a/core/test/input/fd_micromagnetic.cfg b/core/test/input/fd_micromagnetic.cfg new file mode 100644 index 000000000..596f5b958 --- /dev/null +++ b/core/test/input/fd_micromagnetic.cfg @@ -0,0 +1,90 @@ +############## Spirit Configuration ############## + + +### Output Folders +output_file_tag test_pairs_hamiltonian +log_output_folder . +llg_output_folder output +mc_output_folder output +gneb_output_folder output +mmf_output_folder output + + +################## Hamiltonian ################### + +### Hamiltonian Type (heisenberg_neighbours, heisenberg_pairs, gaussian) +hamiltonian micromagnetic + +### boundary_conditions (in a b c) = 0(open), 1(periodical) +boundary_conditions 0 1 1 + +## external magnetic field vector[T] +external_field_magnitude 1.0 +external_field_normal 0.0 0.0 1.0 + +### µSpin +Ms 1.4e6 + +## Uniaxial anisotropy constant [meV] +anisotropy_magnitude 1.0 +anisotropy_normal 0.0 1.0 1.0 + +bravais_vectors +1 0 0 +0 2 0 +0 0 3 + +tensor_exchange +1e-12 0 0 +0 1e-12 0 +0 0 1e-12 + +tensor_dmi +1e-2 0 0 +0 1e-2 0 +0 0 1e-2 + +################ End Hamiltonian ################# + + +############### Logging Parameters ############### +### Save input parameters on creation of State +log_input_save_initial 0 +### Save input parameters on deletion of State +log_input_save_final 0 +### Levels of information +# 0 = ALL - Anything +# 1 = SEVERE - Severe error +# 2 = ERROR - Error which can be handled +# 3 = WARNING - Possible unintended behaviour etc +# 4 = PARAMETER - Input parameter logging +# 5 = INFO - Status information etc +# 6 = DEBUG - Deeper status, eg numerical + +### Print log messages to the console +log_to_console 1 +### Print messages up to (including) log_console_level +log_console_level 5 + +### Save the log as a file +log_to_file 1 +### Save messages up to (including) log_file_level +log_file_level 3 +############# End Logging Parameters ############# + + + +################### Geometry ##################### +### The bravais lattice type + +### Number of basis cells along principal +### directions (a b c) +n_basis_cells 4 2 2 +################# End Geometry ################### + +################ EMA Parameters ################## +### Number of n_modes (for this system max is 6) +ema_n_modes 7 +### Mode to visualize (for this system max is 5) +ema_n_mode_follow 12 +############## End EMA Parameters ################ diff --git a/core/test/input/physics_ddi_micromagnetic.cfg b/core/test/input/physics_ddi_micromagnetic.cfg new file mode 100644 index 000000000..4367884ce --- /dev/null +++ b/core/test/input/physics_ddi_micromagnetic.cfg @@ -0,0 +1,260 @@ +############## Spirit Configuration ############## + + + +### Output Folders +output_file_tag