introduce trivial struct Spin::quantity for an abstract state repre…

…sentation

This allows writing generic code that accesses only the spin part of an
image without a reference to a specific get function. The get function
is still useful when iteration over all fields is required.

core/include/engine/Manifoldmath.hpp

@@ -118,8 +118,9 @@ void hessian_covariant(
 scalar dist_geodesic( const vectorfield & v1, const vectorfield & v2 );
 
 // Calculate the "tangent" vectorfields pointing between a set of configurations
+template<typename StateType>
 void Tangents(
-    std::vector<std::shared_ptr<vectorfield>> configurations, const std::vector<scalar> & energies,
+    const std::vector<std::shared_ptr<StateType>> & configurations, const std::vector<scalar> & energies,
     std::vector<vectorfield> & tangents );
 
 } // namespace Manifoldmath

core/include/engine/Vectormath.inl

@@ -2,6 +2,7 @@
 
 #include <data/Geometry.hpp>
 #include <engine/Backend.hpp>
+#include <engine/StateType.hpp>
 #include <engine/Vectormath_Defines.hpp>
 #include <utility/Constants.hpp>
 
@@ -317,59 +318,59 @@ inline void cross( const vectorfield & vf1, const vectorfield & vf2, vectorfield
 }
 
 template<typename StateType, typename GradientType, typename EnergyFunction>
-void Gradient( const StateType & spins, GradientType & gradient, EnergyFunction && energy, scalar delta )
+void Gradient( const StateType & state, GradientType & gradient, EnergyFunction && energy, scalar delta )
 {
-    static_assert( std::is_convertible_v<decltype( energy( spins ) ), scalar> );
+    static_assert( std::is_convertible_v<decltype( energy( state ) ), scalar> );
 
-    std::size_t nos = spins.size();
+    std::size_t nos = state.spin.size();
 
     // Calculate finite difference
-    vectorfield spins_plus( nos );
-    vectorfield spins_minus( nos );
+    StateType state_plus  = Engine::make_state<StateType>( nos );
+    StateType state_minus = Engine::make_state<StateType>( nos );
 
-    spins_plus  = spins;
-    spins_minus = spins;
+    state_plus  = state;
+    state_minus = state;
 
     for( std::size_t i = 0; i < nos; ++i )
     {
         for( std::uint8_t dim = 0; dim < 3; ++dim )
         {
             // Displace
-            spins_plus[i][dim] += delta;
-            spins_minus[i][dim] -= delta;
+            state_plus.spin[i][dim] += delta;
+            state_minus.spin[i][dim] -= delta;
 
             // Calculate gradient component
-            scalar E_plus    = energy( spins_plus );
-            scalar E_minus   = energy( spins_minus );
+            scalar E_plus    = energy( state_plus );
+            scalar E_minus   = energy( state_minus );
             gradient[i][dim] = 0.5 * ( E_plus - E_minus ) / delta;
 
             // Un-Displace
-            spins_plus[i][dim] -= delta;
-            spins_minus[i][dim] += delta;
+            state_plus.spin[i][dim] -= delta;
+            state_minus.spin[i][dim] += delta;
         }
     }
 }
 
 template<typename StateType, typename HessianType, typename GradientFunction>
-void Hessian( const StateType & spins, HessianType & hessian, GradientFunction && gradient, scalar delta )
+void Hessian( const StateType & state, HessianType & hessian, GradientFunction && gradient, scalar delta )
 {
     static_assert( std::is_invocable_v<GradientFunction, const StateType &, vectorfield &> );
 
     // This is a regular finite difference implementation (probably not very efficient)
     // using the differences between gradient values (not function)
     // see https://v8doc.sas.com/sashtml/ormp/chap5/sect28.htm
 
-    std::size_t nos = spins.size();
+    std::size_t nos = state.spin.size();
 
-    vectorfield spins_pi( nos );
-    vectorfield spins_mi( nos );
-    vectorfield spins_pj( nos );
-    vectorfield spins_mj( nos );
+    StateType state_pi = Engine::make_state<StateType>( nos );
+    StateType state_mi = Engine::make_state<StateType>( nos );
+    StateType state_pj = Engine::make_state<StateType>( nos );
+    StateType state_mj = Engine::make_state<StateType>( nos );
 
-    spins_pi = spins;
-    spins_mi = spins;
-    spins_pj = spins;
-    spins_mj = spins;
+    state_pi = state;
+    state_mi = state;
+    state_pj = state;
+    state_mj = state;
 
     vectorfield grad_pi( nos );
     vectorfield grad_mi( nos );
@@ -385,26 +386,26 @@ void Hessian( const StateType & spins, HessianType & hessian, GradientFunction &
                 for( std::uint8_t beta = 0; beta < 3; ++beta )
                 {
                     // Displace
-                    spins_pi[i][alpha] += delta;
-                    spins_mi[i][alpha] -= delta;
-                    spins_pj[j][beta] += delta;
-                    spins_mj[j][beta] -= delta;
+                    state_pi.spin[i][alpha] += delta;
+                    state_mi.spin[i][alpha] -= delta;
+                    state_pj.spin[j][beta] += delta;
+                    state_mj.spin[j][beta] -= delta;
 
                     // Calculate Hessian component
-                    gradient( spins_pi, grad_pi );
-                    gradient( spins_mi, grad_mi );
-                    gradient( spins_pj, grad_pj );
-                    gradient( spins_mj, grad_mj );
+                    gradient( state_pi, grad_pi );
+                    gradient( state_mi, grad_mi );
+                    gradient( state_pj, grad_pj );
+                    gradient( state_mj, grad_mj );
 
                     hessian( 3 * i + alpha, 3 * j + beta )
                         = 0.25 / delta
                           * ( grad_pj[i][alpha] - grad_mj[i][alpha] + grad_pi[j][beta] - grad_mi[j][beta] );
 
                     // Un-Displace
-                    spins_pi[i][alpha] -= delta;
-                    spins_mi[i][alpha] += delta;
-                    spins_pj[j][beta] -= delta;
-                    spins_mj[j][beta] += delta;
+                    state_pi.spin[i][alpha] -= delta;
+                    state_mi.spin[i][alpha] += delta;
+                    state_pj.spin[j][beta] -= delta;
+                    state_mj.spin[j][beta] += delta;
                 }
             }
         }

core/include/engine/common/interaction/Functor_Prototypes.hpp

@@ -56,7 +56,7 @@ struct Reduce_Functor
     scalar operator()( const typename Interaction::state_t & state ) const
     {
         using std::begin, std::end;
-        scalarfield energy_per_spin( state.size() );
+        scalarfield energy_per_spin( state.spin.size() );
         functor( state, energy_per_spin );
         return Backend::cpu::reduce( SPIRIT_CPU_PAR begin( energy_per_spin ), end( energy_per_spin ) );
     };

core/include/engine/spin/Hamiltonian.hpp

@@ -60,7 +60,7 @@ class Hamiltonian : public Common::Hamiltonian<state_type, StandaloneAdaptorType
 
     void Gradient( const state_t & state, vectorfield & gradient )
     {
-        const auto nos = state.size();
+        const auto nos = state.spin.size();
 
         if( gradient.size() != nos )
             gradient = vectorfield( nos, Vector3::Zero() );
@@ -96,7 +96,7 @@ class Hamiltonian : public Common::Hamiltonian<state_type, StandaloneAdaptorType
 
 struct HamiltonianVariantTypes
 {
-    using state_t     = vectorfield;
+    using state_t     = StateType;
     using AdaptorType = Spin::Interaction::StandaloneAdaptor<state_t>;
 
     using Gaussian   = Hamiltonian<state_t, AdaptorType, Interaction::Gaussian>;

core/include/engine/spin/Interaction_Standalone_Adaptor.hpp

@@ -100,8 +100,8 @@ class StandaloneAdaptor_Local : public Common::Interaction::StandaloneAdaptor_Lo
     void Gradient( const state_t & state, vectorfield & gradient ) final
     {
         using std::begin, std::end;
-        auto functor           = typename InteractionType::Gradient( this->data, this->cache );
-        const auto * state_ptr = state.data();
+        auto functor   = typename InteractionType::Gradient( this->data, this->cache );
+        auto state_ptr = state.data();
         Backend::transform(
             SPIRIT_PAR this->indices.begin(), this->indices.end(), gradient.begin(),
             [state_ptr, functor] SPIRIT_LAMBDA( const IndexTuple & index )

core/include/engine/spin/Method_GNEB.hpp

@@ -32,13 +32,13 @@ class Method_GNEB : public Method_Solver<solver>
     std::string_view Name() override;
 
     void Calculate_Force(
-        const std::vector<std::shared_ptr<vectorfield>> & configurations,
+        const std::vector<std::shared_ptr<StateType>> & configurations,
         std::vector<vectorfield> & forces ) override; // Moved to public, because of cuda device lambda restrictions
 
 private:
     // Calculate Forces onto Systems
     void Calculate_Force_Virtual(
-        const std::vector<std::shared_ptr<vectorfield>> & configurations, const std::vector<vectorfield> & forces,
+        const std::vector<std::shared_ptr<StateType>> & configurations, const std::vector<vectorfield> & forces,
         std::vector<vectorfield> & forces_virtual ) override;
 
     // Check if the Forces are converged

core/include/engine/spin/Method_LLG.hpp

@@ -35,9 +35,9 @@ class Method_LLG : public Method_Solver<solver>
     void Prepare_Thermal_Field() override;
     // Calculate Forces onto Systems
     void Calculate_Force(
-        const std::vector<std::shared_ptr<vectorfield>> & configurations, std::vector<vectorfield> & forces ) override;
+        const std::vector<std::shared_ptr<StateType>> & configurations, std::vector<vectorfield> & forces ) override;
     void Calculate_Force_Virtual(
-        const std::vector<std::shared_ptr<vectorfield>> & configurations, const std::vector<vectorfield> & forces,
+        const std::vector<std::shared_ptr<StateType>> & configurations, const std::vector<vectorfield> & forces,
         std::vector<vectorfield> & forces_virtual ) override;
 
 private:

core/include/engine/spin/Method_MC.hpp

@@ -29,7 +29,7 @@ class Method_MC : public Method
     void Iteration() override;
 
     // Metropolis iteration with adaptive cone radius
-    void Metropolis( const vectorfield & spins_old, vectorfield & spins_new );
+    void Metropolis( const StateType & state_old, StateType & state_new );
 
     // Save the current Step's Data: spins and energy
     void Save_Current( std::string starttime, int iteration, bool initial = false, bool final = false ) override;

core/include/engine/spin/Method_MMF.hpp

@@ -28,7 +28,7 @@ class Method_MMF : public Method_Solver<solver>
 private:
     // Calculate Forces onto Systems
     void Calculate_Force(
-        const std::vector<std::shared_ptr<vectorfield>> & configurations, std::vector<vectorfield> & forces ) override;
+        const std::vector<std::shared_ptr<StateType>> & configurations, std::vector<vectorfield> & forces ) override;
 
     // Check if the Forces are converged
     bool Converged() override;
@@ -66,9 +66,9 @@ class Method_MMF : public Method_Solver<solver>
 
     // Functions for getting the minimum mode of a Hessian
     void Calculate_Force_Spectra_Matrix(
-        const std::vector<std::shared_ptr<vectorfield>> & configurations, std::vector<vectorfield> & forces );
+        const std::vector<std::shared_ptr<StateType>> & configurations, std::vector<vectorfield> & forces );
     void Calculate_Force_Lanczos(
-        const std::vector<std::shared_ptr<vectorfield>> configurations, std::vector<vectorfield> & forces );
+        const std::vector<std::shared_ptr<StateType>> configurations, std::vector<vectorfield> & forces );
 };
 
 } // namespace Spin

core/include/engine/spin/Method_Solver.hpp

@@ -101,19 +101,19 @@ class SolverMethods : public Method
     using Method::Method;
 
     virtual void Prepare_Thermal_Field() = 0;
-    virtual void Calculate_Force(
-        const std::vector<std::shared_ptr<vectorfield>> & configurations, std::vector<vectorfield> & forces )
+    virtual void
+    Calculate_Force( const std::vector<std::shared_ptr<StateType>> & configurations, std::vector<vectorfield> & forces )
         = 0;
     virtual void Calculate_Force_Virtual(
-        const std::vector<std::shared_ptr<vectorfield>> & configurations, const std::vector<vectorfield> & forces,
+        const std::vector<std::shared_ptr<StateType>> & configurations, const std::vector<vectorfield> & forces,
         std::vector<vectorfield> & forces_virtual )
         = 0;
     // Actual Forces on the configurations
     std::vector<vectorfield> forces;
     // Virtual Forces used in the Steps
     std::vector<vectorfield> forces_virtual;
     // Pointers to Configurations (for Solver methods)
-    std::vector<std::shared_ptr<vectorfield>> configurations;
+    std::vector<std::shared_ptr<StateType>> configurations;
 };
 
 // default implementation (to be overwritten by class template specialization)
@@ -175,7 +175,7 @@ class Method_Solver : public SolverData<solver>
      * TODO: maybe rename to separate from deterministic and stochastic force functions
      */
     void Calculate_Force(
-        const std::vector<std::shared_ptr<vectorfield>> & configurations, std::vector<vectorfield> & forces ) override
+        const std::vector<std::shared_ptr<StateType>> & configurations, std::vector<vectorfield> & forces ) override
     {
         Log( Utility::Log_Level::Error, Utility::Log_Sender::All,
              "Tried to use Method_Solver::Calculate_Force() of the Method_Solver class!", this->idx_image,
@@ -190,7 +190,7 @@ class Method_Solver : public SolverData<solver>
      * Default implementation: direct minimization
      */
     void Calculate_Force_Virtual(
-        const std::vector<std::shared_ptr<vectorfield>> & configurations, const std::vector<vectorfield> & forces,
+        const std::vector<std::shared_ptr<StateType>> & configurations, const std::vector<vectorfield> & forces,
         std::vector<vectorfield> & forces_virtual ) override
     {
         // Not Implemented!

core/include/engine/spin/Solver_Depondt.hpp

@@ -26,7 +26,7 @@ class SolverData<Solver::Depondt> : public SolverMethods
     // Virtual Forces used in the Steps
     std::vector<vectorfield> forces_virtual_predictor;
 
-    std::vector<std::shared_ptr<vectorfield>> configurations_predictor;
+    std::vector<std::shared_ptr<StateType>> configurations_predictor;
 
     vectorfield temp1;
 };
@@ -44,9 +44,10 @@ inline void Method_Solver<Solver::Depondt>::Initialize()
     this->angle               = scalarfield( this->nos, 0 );
     this->forces_virtual_norm = std::vector<scalarfield>( this->noi, scalarfield( this->nos, 0 ) );
 
-    this->configurations_predictor = std::vector<std::shared_ptr<vectorfield>>( this->noi );
+    this->configurations_predictor = std::vector<std::shared_ptr<StateType>>( this->noi );
     for( int i = 0; i < this->noi; i++ )
-        configurations_predictor[i] = std::make_shared<vectorfield>( this->nos, Vector3{ 0, 0, 0 } );
+        configurations_predictor[i]
+            = std::make_shared<StateType>( StateType{ vectorfield( this->nos, Vector3{ 0, 0, 0 } ) } );
 
     this->temp1 = vectorfield( this->nos, { 0, 0, 0 } );
 }
@@ -72,7 +73,7 @@ inline void Method_Solver<Solver::Depondt>::Iteration()
     for( int i = 0; i < this->noi; ++i )
     {
         Solver_Kernels::depondt_predictor(
-            forces_virtual[i], rotationaxis[i], angle, *configurations[i], *configurations_predictor[i] );
+            forces_virtual[i], rotationaxis[i], angle, configurations[i]->spin, configurations_predictor[i]->spin );
     }
 
     // Calculate_Force for the Corrector
@@ -84,7 +85,7 @@ inline void Method_Solver<Solver::Depondt>::Iteration()
     for( int i = 0; i < this->noi; i++ )
     {
         Solver_Kernels::depondt_corrector(
-            forces_virtual[i], forces_virtual_predictor[i], temp1, angle, *configurations[i] );
+            forces_virtual[i], forces_virtual_predictor[i], temp1, angle, configurations[i]->spin );
     }
 }
 

core/include/engine/spin/Solver_Heun.hpp

@@ -20,8 +20,8 @@ class SolverData<Solver::Heun> : public SolverMethods
     // Virtual Forces used in the Steps
     std::vector<vectorfield> forces_virtual_predictor;
 
-    std::vector<std::shared_ptr<vectorfield>> configurations_predictor;
-    std::vector<std::shared_ptr<vectorfield>> delta_configurations;
+    std::vector<std::shared_ptr<StateType>> configurations_predictor;
+    std::vector<std::shared_ptr<StateType>> delta_configurations;
 };
 
 template<>
@@ -33,13 +33,13 @@ inline void Method_Solver<Solver::Heun>::Initialize()
     this->forces_predictor         = std::vector<vectorfield>( this->noi, vectorfield( this->nos, { 0, 0, 0 } ) );
     this->forces_virtual_predictor = std::vector<vectorfield>( this->noi, vectorfield( this->nos, { 0, 0, 0 } ) );
 
-    this->configurations_predictor = std::vector<std::shared_ptr<vectorfield>>( this->noi );
+    this->configurations_predictor = std::vector<std::shared_ptr<StateType>>( this->noi );
     for( int i = 0; i < this->noi; i++ )
-        configurations_predictor[i] = std::make_shared<vectorfield>( this->nos );
+        configurations_predictor[i] = std::make_shared<StateType>( make_state<StateType>( this->nos ) );
 
-    this->delta_configurations = std::vector<std::shared_ptr<vectorfield>>( this->noi );
+    this->delta_configurations = std::vector<std::shared_ptr<StateType>>( this->noi );
     for( int i = 0; i < this->noi; i++ )
-        delta_configurations[i] = std::make_shared<vectorfield>( this->nos );
+        delta_configurations[i] = std::make_shared<StateType>( make_state<StateType>( this->nos ) );
 }
 
 /*
@@ -65,8 +65,8 @@ inline void Method_Solver<Solver::Heun>::Iteration()
     {
         // First step - Predictor
         Solver_Kernels::heun_predictor(
-            *this->configurations[i], this->forces_virtual[i], *this->delta_configurations[i],
-            *this->configurations_predictor[i] );
+            this->configurations[i]->spin, this->forces_virtual[i], this->delta_configurations[i]->spin,
+            this->configurations_predictor[i]->spin );
     }
 
     // Calculate_Force for the Corrector
@@ -79,8 +79,8 @@ inline void Method_Solver<Solver::Heun>::Iteration()
     {
         // Second Step - Corrector
         Solver_Kernels::heun_corrector(
-            this->forces_virtual_predictor[i], *this->delta_configurations[i], *this->configurations_predictor[i],
-            *this->configurations[i] );
+            this->forces_virtual_predictor[i], this->delta_configurations[i]->spin,
+            this->configurations_predictor[i]->spin, this->configurations[i]->spin );
     }
 }