From 604cfbb2f2071de7668c9b005015aa9e282e6100 Mon Sep 17 00:00:00 2001 From: puerling Date: Thu, 7 Nov 2024 13:10:32 +0100 Subject: [PATCH] introduce trivial `struct Spin::quantity` for an abstract state representation 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 | 3 +- core/include/engine/Vectormath.inl | 71 ++++++++++--------- .../common/interaction/Functor_Prototypes.hpp | 2 +- core/include/engine/spin/Hamiltonian.hpp | 4 +- .../spin/Interaction_Standalone_Adaptor.hpp | 4 +- core/include/engine/spin/Method_GNEB.hpp | 4 +- core/include/engine/spin/Method_LLG.hpp | 4 +- core/include/engine/spin/Method_MC.hpp | 2 +- core/include/engine/spin/Method_MMF.hpp | 6 +- core/include/engine/spin/Method_Solver.hpp | 12 ++-- core/include/engine/spin/Solver_Depondt.hpp | 11 +-- core/include/engine/spin/Solver_Heun.hpp | 20 +++--- .../engine/spin/Solver_LBFGS_Atlas.hpp | 14 ++-- core/include/engine/spin/Solver_LBFGS_OSO.hpp | 4 +- core/include/engine/spin/Solver_RK4.hpp | 40 +++++------ core/include/engine/spin/Solver_SIB.hpp | 12 ++-- core/include/engine/spin/Solver_VP.hpp | 3 +- core/include/engine/spin/Solver_VP_OSO.hpp | 10 +-- core/include/engine/spin/StateType.hpp | 56 +++++++++++---- .../engine/spin/interaction/Anisotropy.hpp | 17 ++--- .../spin/interaction/Biaxial_Anisotropy.hpp | 27 +++---- .../spin/interaction/Cubic_Anisotropy.hpp | 14 ++-- core/include/engine/spin/interaction/DDI.hpp | 7 +- core/include/engine/spin/interaction/DMI.hpp | 17 ++--- .../engine/spin/interaction/Exchange.hpp | 17 ++--- .../spin/interaction/Functor_Prototypes.hpp | 15 ++-- .../engine/spin/interaction/Gaussian.hpp | 15 ++-- .../engine/spin/interaction/Quadruplet.hpp | 29 ++++---- .../engine/spin/interaction/Zeeman.hpp | 11 +-- core/include/io/Dataparser.inl | 7 +- core/src/Spirit/Chain.cpp | 6 +- core/src/Spirit/Configurations.cpp | 57 +++++++-------- core/src/Spirit/Geometry.cpp | 2 +- core/src/Spirit/IO.cpp | 14 ++-- core/src/Spirit/Quantities.cpp | 20 +++--- core/src/Spirit/State.cpp | 2 +- core/src/Spirit/System.cpp | 2 +- core/src/Spirit/Transitions.cpp | 6 +- core/src/engine/Manifoldmath.cpp | 24 ++++--- core/src/engine/spin/Eigenmodes.cpp | 12 ++-- core/src/engine/spin/HTST.cpp | 20 +++--- core/src/engine/spin/Method_EMA.cpp | 8 +-- core/src/engine/spin/Method_GNEB.cpp | 65 +++++++++-------- core/src/engine/spin/Method_LLG.cpp | 16 ++--- core/src/engine/spin/Method_MC.cpp | 26 +++---- core/src/engine/spin/Method_MMF.cpp | 36 +++++----- core/src/engine/spin/Sparse_HTST.cpp | 16 ++--- core/src/engine/spin/interaction/DDI.cpp | 18 ++--- core/src/io/Dataparser.cpp | 15 ++-- core/src/utility/Configuration_Chain.cpp | 10 +-- core/test/test_aggregation.cpp | 6 +- core/test/test_anisotropy.cpp | 60 +++++++++------- core/test/test_physics.cpp | 30 ++++---- 53 files changed, 493 insertions(+), 436 deletions(-) diff --git a/core/include/engine/Manifoldmath.hpp b/core/include/engine/Manifoldmath.hpp index ca610ad76..159280a63 100644 --- a/core/include/engine/Manifoldmath.hpp +++ b/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 void Tangents( - std::vector> configurations, const std::vector & energies, + const std::vector> & configurations, const std::vector & energies, std::vector & tangents ); } // namespace Manifoldmath diff --git a/core/include/engine/Vectormath.inl b/core/include/engine/Vectormath.inl index 5b7551211..d9e46d7f4 100644 --- a/core/include/engine/Vectormath.inl +++ b/core/include/engine/Vectormath.inl @@ -2,6 +2,7 @@ #include #include +#include #include #include @@ -317,41 +318,41 @@ inline void cross( const vectorfield & vf1, const vectorfield & vf2, vectorfield } template -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 ); + static_assert( std::is_convertible_v ); - 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( nos ); + StateType state_minus = Engine::make_state( 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 -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 ); @@ -359,17 +360,17 @@ void Hessian( const StateType & spins, HessianType & hessian, GradientFunction & // 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( nos ); + StateType state_mi = Engine::make_state( nos ); + StateType state_pj = Engine::make_state( nos ); + StateType state_mj = Engine::make_state( 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; } } } diff --git a/core/include/engine/common/interaction/Functor_Prototypes.hpp b/core/include/engine/common/interaction/Functor_Prototypes.hpp index c6ae63cdb..3344c82c7 100644 --- a/core/include/engine/common/interaction/Functor_Prototypes.hpp +++ b/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 ) ); }; diff --git a/core/include/engine/spin/Hamiltonian.hpp b/core/include/engine/spin/Hamiltonian.hpp index e23afa82c..8c92fc365 100644 --- a/core/include/engine/spin/Hamiltonian.hpp +++ b/core/include/engine/spin/Hamiltonian.hpp @@ -60,7 +60,7 @@ class Hamiltonian : public Common::Hamiltonian; using Gaussian = Hamiltonian; diff --git a/core/include/engine/spin/Interaction_Standalone_Adaptor.hpp b/core/include/engine/spin/Interaction_Standalone_Adaptor.hpp index 5ab108856..1515d12bf 100644 --- a/core/include/engine/spin/Interaction_Standalone_Adaptor.hpp +++ b/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 ) diff --git a/core/include/engine/spin/Method_GNEB.hpp b/core/include/engine/spin/Method_GNEB.hpp index 3c715c130..3a972b0c2 100644 --- a/core/include/engine/spin/Method_GNEB.hpp +++ b/core/include/engine/spin/Method_GNEB.hpp @@ -32,13 +32,13 @@ class Method_GNEB : public Method_Solver std::string_view Name() override; void Calculate_Force( - const std::vector> & configurations, + const std::vector> & configurations, std::vector & forces ) override; // Moved to public, because of cuda device lambda restrictions private: // Calculate Forces onto Systems void Calculate_Force_Virtual( - const std::vector> & configurations, const std::vector & forces, + const std::vector> & configurations, const std::vector & forces, std::vector & forces_virtual ) override; // Check if the Forces are converged diff --git a/core/include/engine/spin/Method_LLG.hpp b/core/include/engine/spin/Method_LLG.hpp index 575f62b1b..8881e6be5 100644 --- a/core/include/engine/spin/Method_LLG.hpp +++ b/core/include/engine/spin/Method_LLG.hpp @@ -35,9 +35,9 @@ class Method_LLG : public Method_Solver void Prepare_Thermal_Field() override; // Calculate Forces onto Systems void Calculate_Force( - const std::vector> & configurations, std::vector & forces ) override; + const std::vector> & configurations, std::vector & forces ) override; void Calculate_Force_Virtual( - const std::vector> & configurations, const std::vector & forces, + const std::vector> & configurations, const std::vector & forces, std::vector & forces_virtual ) override; private: diff --git a/core/include/engine/spin/Method_MC.hpp b/core/include/engine/spin/Method_MC.hpp index 94e5db31c..4098458dc 100644 --- a/core/include/engine/spin/Method_MC.hpp +++ b/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; diff --git a/core/include/engine/spin/Method_MMF.hpp b/core/include/engine/spin/Method_MMF.hpp index b371b3f04..5311af340 100644 --- a/core/include/engine/spin/Method_MMF.hpp +++ b/core/include/engine/spin/Method_MMF.hpp @@ -28,7 +28,7 @@ class Method_MMF : public Method_Solver private: // Calculate Forces onto Systems void Calculate_Force( - const std::vector> & configurations, std::vector & forces ) override; + const std::vector> & configurations, std::vector & forces ) override; // Check if the Forces are converged bool Converged() override; @@ -66,9 +66,9 @@ class Method_MMF : public Method_Solver // Functions for getting the minimum mode of a Hessian void Calculate_Force_Spectra_Matrix( - const std::vector> & configurations, std::vector & forces ); + const std::vector> & configurations, std::vector & forces ); void Calculate_Force_Lanczos( - const std::vector> configurations, std::vector & forces ); + const std::vector> configurations, std::vector & forces ); }; } // namespace Spin diff --git a/core/include/engine/spin/Method_Solver.hpp b/core/include/engine/spin/Method_Solver.hpp index 2f2bd03a2..cff40b9f7 100644 --- a/core/include/engine/spin/Method_Solver.hpp +++ b/core/include/engine/spin/Method_Solver.hpp @@ -101,11 +101,11 @@ class SolverMethods : public Method using Method::Method; virtual void Prepare_Thermal_Field() = 0; - virtual void Calculate_Force( - const std::vector> & configurations, std::vector & forces ) + virtual void + Calculate_Force( const std::vector> & configurations, std::vector & forces ) = 0; virtual void Calculate_Force_Virtual( - const std::vector> & configurations, const std::vector & forces, + const std::vector> & configurations, const std::vector & forces, std::vector & forces_virtual ) = 0; // Actual Forces on the configurations @@ -113,7 +113,7 @@ class SolverMethods : public Method // Virtual Forces used in the Steps std::vector forces_virtual; // Pointers to Configurations (for Solver methods) - std::vector> configurations; + std::vector> configurations; }; // default implementation (to be overwritten by class template specialization) @@ -175,7 +175,7 @@ class Method_Solver : public SolverData * TODO: maybe rename to separate from deterministic and stochastic force functions */ void Calculate_Force( - const std::vector> & configurations, std::vector & forces ) override + const std::vector> & configurations, std::vector & 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 * Default implementation: direct minimization */ void Calculate_Force_Virtual( - const std::vector> & configurations, const std::vector & forces, + const std::vector> & configurations, const std::vector & forces, std::vector & forces_virtual ) override { // Not Implemented! diff --git a/core/include/engine/spin/Solver_Depondt.hpp b/core/include/engine/spin/Solver_Depondt.hpp index f835e678c..770a35f06 100644 --- a/core/include/engine/spin/Solver_Depondt.hpp +++ b/core/include/engine/spin/Solver_Depondt.hpp @@ -26,7 +26,7 @@ class SolverData : public SolverMethods // Virtual Forces used in the Steps std::vector forces_virtual_predictor; - std::vector> configurations_predictor; + std::vector> configurations_predictor; vectorfield temp1; }; @@ -44,9 +44,10 @@ inline void Method_Solver::Initialize() this->angle = scalarfield( this->nos, 0 ); this->forces_virtual_norm = std::vector( this->noi, scalarfield( this->nos, 0 ) ); - this->configurations_predictor = std::vector>( this->noi ); + this->configurations_predictor = std::vector>( this->noi ); for( int i = 0; i < this->noi; i++ ) - configurations_predictor[i] = std::make_shared( this->nos, Vector3{ 0, 0, 0 } ); + configurations_predictor[i] + = std::make_shared( StateType{ vectorfield( this->nos, Vector3{ 0, 0, 0 } ) } ); this->temp1 = vectorfield( this->nos, { 0, 0, 0 } ); } @@ -72,7 +73,7 @@ inline void Method_Solver::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::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 ); } } diff --git a/core/include/engine/spin/Solver_Heun.hpp b/core/include/engine/spin/Solver_Heun.hpp index 87940d2b8..ed79e18d6 100644 --- a/core/include/engine/spin/Solver_Heun.hpp +++ b/core/include/engine/spin/Solver_Heun.hpp @@ -20,8 +20,8 @@ class SolverData : public SolverMethods // Virtual Forces used in the Steps std::vector forces_virtual_predictor; - std::vector> configurations_predictor; - std::vector> delta_configurations; + std::vector> configurations_predictor; + std::vector> delta_configurations; }; template<> @@ -33,13 +33,13 @@ inline void Method_Solver::Initialize() this->forces_predictor = std::vector( this->noi, vectorfield( this->nos, { 0, 0, 0 } ) ); this->forces_virtual_predictor = std::vector( this->noi, vectorfield( this->nos, { 0, 0, 0 } ) ); - this->configurations_predictor = std::vector>( this->noi ); + this->configurations_predictor = std::vector>( this->noi ); for( int i = 0; i < this->noi; i++ ) - configurations_predictor[i] = std::make_shared( this->nos ); + configurations_predictor[i] = std::make_shared( make_state( this->nos ) ); - this->delta_configurations = std::vector>( this->noi ); + this->delta_configurations = std::vector>( this->noi ); for( int i = 0; i < this->noi; i++ ) - delta_configurations[i] = std::make_shared( this->nos ); + delta_configurations[i] = std::make_shared( make_state( this->nos ) ); } /* @@ -65,8 +65,8 @@ inline void Method_Solver::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::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 ); } } diff --git a/core/include/engine/spin/Solver_LBFGS_Atlas.hpp b/core/include/engine/spin/Solver_LBFGS_Atlas.hpp index 216c53c91..16cc9e0d0 100644 --- a/core/include/engine/spin/Solver_LBFGS_Atlas.hpp +++ b/core/include/engine/spin/Solver_LBFGS_Atlas.hpp @@ -62,7 +62,7 @@ inline void Method_Solver::Initialize() // Choose atlas3 coordinates for( int i = 0; i < this->nos; i++ ) { - this->atlas_coords3[img][i] = ( *this->configurations[img] )[i][2] > 0 ? 1.0 : -1.0; + this->atlas_coords3[img][i] = this->configurations[img]->spin[i][2] > 0 ? 1.0 : -1.0; // Solver_Kernels::ncg_spins_to_atlas( *this->configurations[i], this->atlas_coords[i], this->atlas_coords3[i] ); } } @@ -77,14 +77,14 @@ template<> inline void Method_Solver::Iteration() { int noi = configurations.size(); - int nos = ( *configurations[0] ).size(); + int nos = configurations[0]->spin.size(); // Current force this->Calculate_Force( this->configurations, this->forces ); for( int img = 0; img < this->noi; img++ ) { - auto & image = *this->configurations[img]; + auto & image = this->configurations[img]->spin; auto & grad_ref = this->atlas_residuals[img]; Backend::transform( @@ -121,22 +121,22 @@ inline void Method_Solver::Iteration() // Rotate spins Solver_Kernels::atlas_rotate( - *this->configurations[img], this->atlas_coords3[img], this->atlas_directions[img] ); + this->configurations[img]->spin, this->atlas_coords3[img], this->atlas_directions[img] ); } if( std::any_of( Backend::cpu::make_zip_iterator( this->configurations.begin(), this->atlas_coords3.begin() ), Backend::cpu::make_zip_iterator( this->configurations.end(), this->atlas_coords3.end() ), Backend::cpu::make_zip_function( - []( const std::shared_ptr & spins, const scalarfield & a3 ) -> bool - { return Solver_Kernels::ncg_atlas_check_coordinates( *spins, a3, -0.6 ); } ) ) ) + []( const std::shared_ptr & state, const scalarfield & a3 ) -> bool + { return Solver_Kernels::ncg_atlas_check_coordinates( state->spin, a3, -0.6 ); } ) ) ) { const auto m_inverse = [] SPIRIT_LAMBDA( const scalar s ) { return scalar( 1.0 ) / s; }; Backend::transform( SPIRIT_PAR rho.begin(), rho.end(), rho.begin(), m_inverse ); for( int img = 0; img < noi; ++img ) Solver_Kernels::lbfgs_atlas_transform_direction( - *this->configurations[img], this->atlas_coords3[img], this->atlas_updates[img], + this->configurations[img]->spin, this->atlas_coords3[img], this->atlas_updates[img], this->grad_atlas_updates[img], this->atlas_directions[img], this->atlas_residuals_last[img], this->rho ); diff --git a/core/include/engine/spin/Solver_LBFGS_OSO.hpp b/core/include/engine/spin/Solver_LBFGS_OSO.hpp index d5511b790..39c3eb1fb 100644 --- a/core/include/engine/spin/Solver_LBFGS_OSO.hpp +++ b/core/include/engine/spin/Solver_LBFGS_OSO.hpp @@ -69,7 +69,7 @@ inline void Method_Solver::Iteration() // calculate gradients for OSO for( int img = 0; img < this->noi; img++ ) { - auto & image = *this->configurations[img]; + auto & image = this->configurations[img]->spin; auto & grad_ref = this->grad[img]; const auto * f = this->forces[img].data(); @@ -97,7 +97,7 @@ inline void Method_Solver::Iteration() { Vectormath::scale( searchdir[img], scaling ); // rotate spins - Solver_Kernels::oso_rotate( *this->configurations[img], this->searchdir[img] ); + Solver_Kernels::oso_rotate( this->configurations[img]->spin, this->searchdir[img] ); } } diff --git a/core/include/engine/spin/Solver_RK4.hpp b/core/include/engine/spin/Solver_RK4.hpp index d5bfbb74a..6c7635814 100644 --- a/core/include/engine/spin/Solver_RK4.hpp +++ b/core/include/engine/spin/Solver_RK4.hpp @@ -21,11 +21,11 @@ class SolverData : public SolverMethods // Virtual Forces used in the Steps std::vector forces_virtual_predictor; - std::vector> configurations_predictor; + std::vector> configurations_predictor; - std::vector> configurations_k1; - std::vector> configurations_k2; - std::vector> configurations_k3; + std::vector> configurations_k1; + std::vector> configurations_k2; + std::vector> configurations_k3; }; template<> @@ -37,21 +37,21 @@ inline void Method_Solver::Initialize() this->forces_predictor = std::vector( this->noi, vectorfield( this->nos, { 0, 0, 0 } ) ); this->forces_virtual_predictor = std::vector( this->noi, vectorfield( this->nos, { 0, 0, 0 } ) ); - this->configurations_predictor = std::vector>( this->noi ); + this->configurations_predictor = std::vector>( this->noi ); for( int i = 0; i < this->noi; i++ ) - this->configurations_predictor[i] = std::make_shared( this->nos ); + this->configurations_predictor[i] = std::make_shared( make_state( this->nos ) ); - this->configurations_k1 = std::vector>( this->noi ); + this->configurations_k1 = std::vector>( this->noi ); for( int i = 0; i < this->noi; i++ ) - this->configurations_k1[i] = std::make_shared( this->nos ); + this->configurations_k1[i] = std::make_shared( make_state( this->nos ) ); - this->configurations_k2 = std::vector>( this->noi ); + this->configurations_k2 = std::vector>( this->noi ); for( int i = 0; i < this->noi; i++ ) - this->configurations_k2[i] = std::make_shared( this->nos ); + this->configurations_k2[i] = std::make_shared( make_state( this->nos ) ); - this->configurations_k3 = std::vector>( this->noi ); + this->configurations_k3 = std::vector>( this->noi ); for( int i = 0; i < this->noi; i++ ) - this->configurations_k3[i] = std::make_shared( this->nos ); + this->configurations_k3[i] = std::make_shared( make_state( this->nos ) ); } /* @@ -71,8 +71,8 @@ inline void Method_Solver::Iteration() for( int i = 0; i < this->noi; ++i ) { Solver_Kernels::rk4_predictor_1( - *this->configurations[i], this->forces_virtual[i], *this->configurations_k1[i], - *this->configurations_predictor[i] ); + this->configurations[i]->spin, this->forces_virtual[i], this->configurations_k1[i]->spin, + this->configurations_predictor[i]->spin ); } // Calculate_Force for the predictor @@ -84,8 +84,8 @@ inline void Method_Solver::Iteration() for( int i = 0; i < this->noi; ++i ) { Solver_Kernels::rk4_predictor_2( - *this->configurations[i], this->forces_virtual_predictor[i], *this->configurations_k2[i], - *this->configurations_predictor[i] ); + this->configurations[i]->spin, this->forces_virtual_predictor[i], this->configurations_k2[i]->spin, + this->configurations_predictor[i]->spin ); } // Calculate_Force for the predictor (k3) @@ -97,8 +97,8 @@ inline void Method_Solver::Iteration() for( int i = 0; i < this->noi; ++i ) { Solver_Kernels::rk4_predictor_3( - *this->configurations[i], this->forces_virtual_predictor[i], *this->configurations_k3[i], - *this->configurations_predictor[i] ); + this->configurations[i]->spin, this->forces_virtual_predictor[i], this->configurations_k3[i]->spin, + this->configurations_predictor[i]->spin ); } // Calculate_Force for the predictor (k4) @@ -110,8 +110,8 @@ inline void Method_Solver::Iteration() for( int i = 0; i < this->noi; i++ ) { Solver_Kernels::rk4_corrector( - forces_virtual_predictor[i], *configurations_k1[i], *configurations_k2[i], *configurations_k3[i], - *configurations_predictor[i], *configurations[i] ); + forces_virtual_predictor[i], configurations_k1[i]->spin, configurations_k2[i]->spin, + configurations_k3[i]->spin, configurations_predictor[i]->spin, configurations[i]->spin ); } } diff --git a/core/include/engine/spin/Solver_SIB.hpp b/core/include/engine/spin/Solver_SIB.hpp index 96e6751db..ab937e680 100644 --- a/core/include/engine/spin/Solver_SIB.hpp +++ b/core/include/engine/spin/Solver_SIB.hpp @@ -20,7 +20,7 @@ class SolverData : public SolverMethods // Virtual Forces used in the Steps std::vector forces_virtual_predictor; - std::vector> configurations_predictor; + std::vector> configurations_predictor; }; template<> @@ -32,9 +32,9 @@ inline void Method_Solver::Initialize() this->forces_predictor = std::vector( this->noi, vectorfield( this->nos, { 0, 0, 0 } ) ); this->forces_virtual_predictor = std::vector( this->noi, vectorfield( this->nos, { 0, 0, 0 } ) ); - this->configurations_predictor = std::vector>( this->noi ); + this->configurations_predictor = std::vector>( this->noi ); for( int i = 0; i < this->noi; i++ ) - configurations_predictor[i] = std::make_shared( this->nos ); + configurations_predictor[i] = std::make_shared( make_state( this->nos ) ); } /* @@ -54,8 +54,8 @@ inline void Method_Solver::Iteration() this->Calculate_Force_Virtual( this->configurations, this->forces, this->forces_virtual ); for( int i = 0; i < this->noi; ++i ) { - auto & image = *this->configurations[i]; - auto & predictor = *this->configurations_predictor[i]; + auto & image = this->configurations[i]->spin; + auto & predictor = this->configurations_predictor[i]->spin; Solver_Kernels::sib_transform( image, forces_virtual[i], predictor ); Backend::transform( @@ -70,7 +70,7 @@ inline void Method_Solver::Iteration() this->configurations_predictor, this->forces_predictor, this->forces_virtual_predictor ); for( int i = 0; i < this->noi; ++i ) { - auto & image = *this->configurations[i]; + auto & image = this->configurations[i]->spin; Solver_Kernels::sib_transform( image, forces_virtual_predictor[i], image ); } diff --git a/core/include/engine/spin/Solver_VP.hpp b/core/include/engine/spin/Solver_VP.hpp index 3829a8e1e..7ed76a484 100644 --- a/core/include/engine/spin/Solver_VP.hpp +++ b/core/include/engine/spin/Solver_VP.hpp @@ -77,7 +77,8 @@ inline void Method_Solver::Iteration() Solver_Kernels::VP::projected_velocity( projections, forces[i], velocities[i] ); // Apply the projected velocity - Solver_Kernels::VP::apply_velocity( velocities[i], forces[i], this->llg_parameters[i]->dt, *configurations[i] ); + Solver_Kernels::VP::apply_velocity( + velocities[i], forces[i], this->llg_parameters[i]->dt, configurations[i]->spin ); } } diff --git a/core/include/engine/spin/Solver_VP_OSO.hpp b/core/include/engine/spin/Solver_VP_OSO.hpp index 5bfa82245..e1ea3127b 100644 --- a/core/include/engine/spin/Solver_VP_OSO.hpp +++ b/core/include/engine/spin/Solver_VP_OSO.hpp @@ -25,7 +25,7 @@ class SolverData : public SolverMethods std::vector grad_pr; std::vector searchdir; - std::vector> configurations_temp; + std::vector> configurations_temp; std::vector> llg_parameters; }; @@ -36,9 +36,9 @@ inline void Method_Solver::Initialize() this->forces = std::vector( this->noi, vectorfield( this->nos, { 0, 0, 0 } ) ); this->forces_virtual = std::vector( this->noi, vectorfield( this->nos, { 0, 0, 0 } ) ); - this->configurations_temp = std::vector>( this->noi ); + this->configurations_temp = std::vector>( this->noi ); for( int i = 0; i < this->noi; i++ ) - configurations_temp[i] = std::make_shared( this->nos ); + configurations_temp[i] = std::make_shared( make_state( this->nos ) ); this->velocities = std::vector( this->noi, vectorfield( this->nos, Vector3::Zero() ) ); // [noi][nos] this->forces_previous = forces; // [noi][nos] @@ -78,7 +78,7 @@ inline void Method_Solver::Iteration() for( int img = 0; img < this->noi; img++ ) { - auto & image = *this->configurations[img]; + auto & image = this->configurations[img]->spin; auto & grad = this->grad[img]; Solver_Kernels::oso_calc_gradients( grad, image, this->forces[img] ); Vectormath::scale( grad, -1.0 ); @@ -101,7 +101,7 @@ inline void Method_Solver::Iteration() Solver_Kernels::VP::set_step( velocities[img], grad[img], this->llg_parameters[img]->dt, searchdir[img] ); // rotate spins - Solver_Kernels::oso_rotate( *this->configurations[img], this->searchdir[img] ); + Solver_Kernels::oso_rotate( this->configurations[img]->spin, this->searchdir[img] ); } } diff --git a/core/include/engine/spin/StateType.hpp b/core/include/engine/spin/StateType.hpp index 970117d54..5837705f7 100644 --- a/core/include/engine/spin/StateType.hpp +++ b/core/include/engine/spin/StateType.hpp @@ -7,6 +7,9 @@ namespace Engine { +template +struct state_traits; + namespace Spin { @@ -16,35 +19,62 @@ enum struct Field }; template -using quantity = T; +struct quantity +{ + T spin; +}; + +template +struct quantity> +{ + field spin; + + auto data() -> typename state_traits>>::pointer + { + return { spin.data() }; + } + + auto data() const -> typename state_traits>>::const_pointer + { + return { spin.data() }; + } +}; + +template +constexpr auto make_quantity( T && value ) -> quantity> +{ + return { std::forward( value ) }; +} template -T & get( T & q ) +T & get( quantity & q ) { - return q; + if constexpr( field == Field::Spin ) + return q.spin; } template const T & get( const T & q ) { - return q; + if constexpr( field == Field::Spin ) + return q.spin; } -using StateType = vectorfield; -using StatePtr = Vector3 *; -using StateCPtr = const Vector3 *; +using StateType = quantity; +using StatePtr = quantity; +using StateCPtr = quantity; } // namespace Spin template struct state_traits; -template<> -struct state_traits +template +struct state_traits> { - using type = Spin::StateType; - using pointer = Spin::StateType::pointer; - using const_pointer = Spin::StateType::const_pointer; + using type = Spin::quantity; + using pointer = Spin::quantity; + using const_pointer = Spin::quantity; }; template @@ -53,7 +83,7 @@ state_t make_state( int nos ); template<> inline Spin::StateType make_state( int nos ) { - return vectorfield( nos ); + return { vectorfield( nos ) }; }; } // namespace Engine diff --git a/core/include/engine/spin/interaction/Anisotropy.hpp b/core/include/engine/spin/interaction/Anisotropy.hpp index 2225fd719..1845f3d71 100644 --- a/core/include/engine/spin/interaction/Anisotropy.hpp +++ b/core/include/engine/spin/interaction/Anisotropy.hpp @@ -3,6 +3,7 @@ #define SPIRIT_CORE_ENGINE_INTERACTION_ANISOTROPY_HPP #include +#include #include #include @@ -18,7 +19,7 @@ namespace Interaction struct Anisotropy { - using state_t = vectorfield; + using state_t = StateType; struct Data { @@ -116,31 +117,31 @@ struct Functor::Local::DataRef }; template<> -inline scalar Anisotropy::Energy::operator()( const Index & index, const Vector3 * spins ) const +inline scalar Anisotropy::Energy::operator()( const Index & index, quantity state ) const { return -1.0 * Backend::transform_reduce( index.begin(), index.end(), scalar( 0 ), Backend::plus{}, - [this, spins] SPIRIT_LAMBDA( const Interaction::IndexType & idx ) -> scalar + [this, state] SPIRIT_LAMBDA( const Interaction::IndexType & idx ) -> scalar { - const auto d = normals[idx.iani].dot( spins[idx.ispin] ); + const auto d = normals[idx.iani].dot( state.spin[idx.ispin] ); return magnitudes[idx.iani] * d * d; } ); } template<> -inline Vector3 Anisotropy::Gradient::operator()( const Index & index, const Vector3 * spins ) const +inline Vector3 Anisotropy::Gradient::operator()( const Index & index, quantity state ) const { return -2.0 * Backend::transform_reduce( index.begin(), index.end(), Vector3( Vector3::Zero() ), Backend::plus{}, - [this, spins] SPIRIT_LAMBDA( const Interaction::IndexType & idx ) -> Vector3 - { return magnitudes[idx.iani] * normals[idx.iani].dot( spins[idx.ispin] ) * normals[idx.iani]; } ); + [this, state] SPIRIT_LAMBDA( const Interaction::IndexType & idx ) -> Vector3 + { return magnitudes[idx.iani] * normals[idx.iani].dot( state.spin[idx.ispin] ) * normals[idx.iani]; } ); } template<> template -void Anisotropy::Hessian::operator()( const Index & index, const vectorfield &, Callable & hessian ) const +void Anisotropy::Hessian::operator()( const Index & index, const StateType &, Callable & hessian ) const { Backend::cpu::for_each( index.begin(), index.end(), diff --git a/core/include/engine/spin/interaction/Biaxial_Anisotropy.hpp b/core/include/engine/spin/interaction/Biaxial_Anisotropy.hpp index e4eb5e124..cafc4883d 100644 --- a/core/include/engine/spin/interaction/Biaxial_Anisotropy.hpp +++ b/core/include/engine/spin/interaction/Biaxial_Anisotropy.hpp @@ -3,6 +3,7 @@ #define SPIRIT_CORE_ENGINE_INTERACTION_BIAXIAL_ANISOTROPY_HPP #include +#include #include #include @@ -24,7 +25,7 @@ namespace Interaction */ struct Biaxial_Anisotropy { - using state_t = vectorfield; + using state_t = StateType; struct Data { @@ -127,7 +128,7 @@ struct Functor::Local::DataRef }; template<> -inline scalar Biaxial_Anisotropy::Energy::operator()( const Index & index, const Vector3 * spins ) const +inline scalar Biaxial_Anisotropy::Energy::operator()( const Index & index, quantity state ) const { using Utility::fastpow; scalar result = 0; @@ -135,9 +136,9 @@ inline scalar Biaxial_Anisotropy::Energy::operator()( const Index & index, const return result; const auto & [ispin, iani] = *index; - const scalar s1 = bases[iani].k1.dot( spins[ispin] ); - const scalar s2 = bases[iani].k2.dot( spins[ispin] ); - const scalar s3 = bases[iani].k3.dot( spins[ispin] ); + const scalar s1 = bases[iani].k1.dot( state.spin[ispin] ); + const scalar s2 = bases[iani].k2.dot( state.spin[ispin] ); + const scalar s3 = bases[iani].k3.dot( state.spin[ispin] ); const scalar sin_theta_2 = 1 - s1 * s1; @@ -151,7 +152,7 @@ inline scalar Biaxial_Anisotropy::Energy::operator()( const Index & index, const } template<> -inline Vector3 Biaxial_Anisotropy::Gradient::operator()( const Index & index, const Vector3 * spins ) const +inline Vector3 Biaxial_Anisotropy::Gradient::operator()( const Index & index, quantity state ) const { using Utility::fastpow; Vector3 result = Vector3::Zero(); @@ -161,9 +162,9 @@ inline Vector3 Biaxial_Anisotropy::Gradient::operator()( const Index & index, co const auto & [ispin, iani] = *index; const auto & [k1, k2, k3] = bases[iani]; - const scalar s1 = k1.dot( spins[ispin] ); - const scalar s2 = k2.dot( spins[ispin] ); - const scalar s3 = k3.dot( spins[ispin] ); + const scalar s1 = k1.dot( state.spin[ispin] ); + const scalar s2 = k2.dot( state.spin[ispin] ); + const scalar s3 = k3.dot( state.spin[ispin] ); const scalar sin_theta_2 = 1 - s1 * s1; @@ -188,7 +189,7 @@ inline Vector3 Biaxial_Anisotropy::Gradient::operator()( const Index & index, co template<> template -void Biaxial_Anisotropy::Hessian::operator()( const Index & index, const vectorfield & spins, Callable & hessian ) const +void Biaxial_Anisotropy::Hessian::operator()( const Index & index, const StateType & state, Callable & hessian ) const { using Utility::fastpow; if( !is_contributing || index == nullptr ) @@ -197,9 +198,9 @@ void Biaxial_Anisotropy::Hessian::operator()( const Index & index, const vectorf const auto & [ispin, iani] = *index; const auto & [k1, k2, k3] = bases[iani]; - const scalar s1 = k1.dot( spins[ispin] ); - const scalar s2 = k2.dot( spins[ispin] ); - const scalar s3 = k3.dot( spins[ispin] ); + const scalar s1 = k1.dot( state.spin[ispin] ); + const scalar s2 = k2.dot( state.spin[ispin] ); + const scalar s3 = k3.dot( state.spin[ispin] ); const scalar st2 = 1 - s1 * s1; diff --git a/core/include/engine/spin/interaction/Cubic_Anisotropy.hpp b/core/include/engine/spin/interaction/Cubic_Anisotropy.hpp index 692622dee..fcb43f9ca 100644 --- a/core/include/engine/spin/interaction/Cubic_Anisotropy.hpp +++ b/core/include/engine/spin/interaction/Cubic_Anisotropy.hpp @@ -3,6 +3,7 @@ #define SPIRIT_CORE_ENGINE_INTERACTION_CUBIC_ANISOTROPY_HPP #include +#include #include #include @@ -19,7 +20,7 @@ namespace Interaction struct Cubic_Anisotropy { - using state_t = vectorfield; + using state_t = StateType; struct Data { @@ -106,7 +107,7 @@ struct Functor::Local::DataRef }; template<> -inline scalar Cubic_Anisotropy::Energy::operator()( const Index & index, const Vector3 * spins ) const +inline scalar Cubic_Anisotropy::Energy::operator()( const Index & index, quantity state ) const { using Utility::fastpow; scalar result = 0; @@ -115,11 +116,12 @@ inline scalar Cubic_Anisotropy::Energy::operator()( const Index & index, const V const auto & [ispin, iani] = *index; return -0.5 * magnitudes[iani] - * ( fastpow( spins[ispin][0], 4u ) + fastpow( spins[ispin][1], 4u ) + fastpow( spins[ispin][2], 4u ) ); + * ( fastpow( state.spin[ispin][0], 4u ) + fastpow( state.spin[ispin][1], 4u ) + + fastpow( state.spin[ispin][2], 4u ) ); } template<> -inline Vector3 Cubic_Anisotropy::Gradient::operator()( const Index & index, const Vector3 * spins ) const +inline Vector3 Cubic_Anisotropy::Gradient::operator()( const Index & index, quantity state ) const { using Utility::fastpow; Vector3 result = Vector3::Zero(); @@ -130,14 +132,14 @@ inline Vector3 Cubic_Anisotropy::Gradient::operator()( const Index & index, cons for( int icomp = 0; icomp < 3; ++icomp ) { - result[icomp] = -2.0 * magnitudes[iani] * fastpow( spins[ispin][icomp], 3u ); + result[icomp] = -2.0 * magnitudes[iani] * fastpow( state.spin[ispin][icomp], 3u ); } return result; } template<> template -void Cubic_Anisotropy::Hessian::operator()( const Index & index, const vectorfield & spins, Callable & hessian ) const +void Cubic_Anisotropy::Hessian::operator()( const Index & index, const StateType & spins, Callable & hessian ) const { // TODO: Not yet implemented } diff --git a/core/include/engine/spin/interaction/DDI.hpp b/core/include/engine/spin/interaction/DDI.hpp index 70cfd51f1..4bae5b26e 100644 --- a/core/include/engine/spin/interaction/DDI.hpp +++ b/core/include/engine/spin/interaction/DDI.hpp @@ -3,6 +3,7 @@ #define SPIRIT_CORE_ENGINE_INTERACTION_DDI_HPP #include +#include #include #include @@ -25,7 +26,7 @@ namespace Interaction struct DDI { - using state_t = vectorfield; + using state_t = StateType; struct Data { @@ -110,7 +111,7 @@ struct DDI template<> template -void DDI::Hessian::operator()( const vectorfield & spins, Callable & hessian ) const +void DDI::Hessian::operator()( const StateType & state, Callable & hessian ) const { namespace C = Utility::Constants; if( !is_contributing ) @@ -121,7 +122,7 @@ void DDI::Hessian::operator()( const vectorfield & spins, Callable & hessian ) c return; const auto & geometry = *cache.geometry; - const auto nos = spins.size(); + const auto nos = state.spin.size(); // Tentative Dipole-Dipole (only works for open boundary conditions) if( data.method != DDI_Method::None ) diff --git a/core/include/engine/spin/interaction/DMI.hpp b/core/include/engine/spin/interaction/DMI.hpp index e4b568594..724fbcdee 100644 --- a/core/include/engine/spin/interaction/DMI.hpp +++ b/core/include/engine/spin/interaction/DMI.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include @@ -20,7 +21,7 @@ namespace Interaction struct DMI { - using state_t = vectorfield; + using state_t = StateType; struct Data { @@ -167,35 +168,35 @@ struct Functor::Local::DataRef }; template<> -inline scalar DMI::Energy::operator()( const Index & index, const Vector3 * spins ) const +inline scalar DMI::Energy::operator()( const Index & index, quantity state ) const { // don't need to check for `is_contributing` here, because `transform_reduce` will short circuit properly return Backend::transform_reduce( index.begin(), index.end(), scalar( 0.0 ), Backend::plus{}, - [this, spins] SPIRIT_LAMBDA( const DMI::IndexType & idx ) -> scalar + [this, state] SPIRIT_LAMBDA( const DMI::IndexType & idx ) -> scalar { const auto & [ispin, jspin, i_pair, inverse] = idx; return ( inverse ? 0.5 : -0.5 ) * magnitudes[i_pair] - * normals[i_pair].dot( spins[ispin].cross( spins[jspin] ) ); + * normals[i_pair].dot( state.spin[ispin].cross( state.spin[jspin] ) ); } ); } template<> -inline Vector3 DMI::Gradient::operator()( const Index & index, const Vector3 * spins ) const +inline Vector3 DMI::Gradient::operator()( const Index & index, quantity state ) const { // don't need to check for `is_contributing` here, because `transform_reduce` will short circuit properly return Backend::transform_reduce( index.begin(), index.end(), Vector3{ 0.0, 0.0, 0.0 }, Backend::plus{}, - [this, spins] SPIRIT_LAMBDA( const DMI::IndexType & idx ) -> Vector3 + [this, state] SPIRIT_LAMBDA( const DMI::IndexType & idx ) -> Vector3 { const auto & [ispin, jspin, i_pair, inverse] = idx; - return ( inverse ? 1.0 : -1.0 ) * magnitudes[i_pair] * spins[jspin].cross( normals[i_pair] ); + return ( inverse ? 1.0 : -1.0 ) * magnitudes[i_pair] * state.spin[jspin].cross( normals[i_pair] ); } ); } template<> template -void DMI::Hessian::operator()( const Index & index, const vectorfield &, Callable & hessian ) const +void DMI::Hessian::operator()( const Index & index, const StateType &, Callable & hessian ) const { // don't need to check for `is_contributing` here, because `for_each` will short circuit properly Backend::cpu::for_each( diff --git a/core/include/engine/spin/interaction/Exchange.hpp b/core/include/engine/spin/interaction/Exchange.hpp index 99abaaf63..d04e5bad9 100644 --- a/core/include/engine/spin/interaction/Exchange.hpp +++ b/core/include/engine/spin/interaction/Exchange.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include @@ -20,7 +21,7 @@ namespace Interaction struct Exchange { - using state_t = vectorfield; + using state_t = StateType; struct Data { @@ -147,34 +148,34 @@ struct Functor::Local::DataRef }; template<> -inline scalar Exchange::Energy::operator()( const Index & index, const Vector3 * spins ) const +inline scalar Exchange::Energy::operator()( const Index & index, quantity state ) const { // don't need to check for `is_contributing` here, because the `transform_reduce` will short circuit correctly return Backend::transform_reduce( index.begin(), index.end(), scalar( 0.0 ), Backend::plus{}, - [this, spins] SPIRIT_LAMBDA( const Exchange::IndexType & idx ) -> scalar + [this, state] SPIRIT_LAMBDA( const Exchange::IndexType & idx ) -> scalar { const auto & [ispin, jspin, i_pair] = idx; - return -0.5 * magnitudes[i_pair] * spins[ispin].dot( spins[jspin] ); + return -0.5 * magnitudes[i_pair] * state.spin[ispin].dot( state.spin[jspin] ); } ); } template<> -inline Vector3 Exchange::Gradient::operator()( const Index & index, const Vector3 * spins ) const +inline Vector3 Exchange::Gradient::operator()( const Index & index, quantity state ) const { // don't need to check for `is_contributing` here, because the `transform_reduce` will short circuit correctly return Backend::transform_reduce( index.begin(), index.end(), Vector3{ 0.0, 0.0, 0.0 }, Backend::plus{}, - [this, spins] SPIRIT_LAMBDA( const Exchange::IndexType & idx ) -> Vector3 + [this, state] SPIRIT_LAMBDA( const Exchange::IndexType & idx ) -> Vector3 { const auto & [ispin, jspin, i_pair] = idx; - return -magnitudes[i_pair] * spins[jspin]; + return -magnitudes[i_pair] * state.spin[jspin]; } ); } template<> template -void Exchange::Hessian::operator()( const Index & index, const vectorfield &, Callable & hessian ) const +void Exchange::Hessian::operator()( const Index & index, const StateType &, Callable & hessian ) const { // don't need to check for `is_contributing` here, because `for_each` will short circuit properly Backend::cpu::for_each( diff --git a/core/include/engine/spin/interaction/Functor_Prototypes.hpp b/core/include/engine/spin/interaction/Functor_Prototypes.hpp index 1249eaab3..65b1f5c43 100644 --- a/core/include/engine/spin/interaction/Functor_Prototypes.hpp +++ b/core/include/engine/spin/interaction/Functor_Prototypes.hpp @@ -7,6 +7,7 @@ #include #include #include +#include #include @@ -78,7 +79,7 @@ struct Energy_Functor : public DataRef using Data = typename Interaction::Data; using Cache = typename Interaction::Cache; - void operator()( const vectorfield & spins, scalarfield & energy ) const; + void operator()( const StateType & state, scalarfield & energy ) const; using DataRef::DataRef; }; @@ -90,7 +91,7 @@ struct Energy_Single_Spin_Functor : public DataRef using Data = typename Interaction::Data; using Cache = typename Interaction::Cache; - scalar operator()( int ispin, const vectorfield & spins ) const; + scalar operator()( int ispin, const StateType & state ) const; using DataRef::DataRef; }; @@ -102,7 +103,7 @@ struct Gradient_Functor : public DataRef using Data = typename Interaction::Data; using Cache = typename Interaction::Cache; - void operator()( const vectorfield & spins, vectorfield & gradient ) const; + void operator()( const StateType & state, vectorfield & gradient ) const; using DataRef::DataRef; }; @@ -115,7 +116,7 @@ struct Hessian_Functor : public DataRef using Cache = typename Interaction::Cache; template - void operator()( const vectorfield & spins, Callable & hessian ) const; + void operator()( const StateType & state, Callable & hessian ) const; using DataRef::DataRef; }; @@ -150,7 +151,7 @@ struct Energy_Functor : public DataRef using Cache = typename Interaction::Cache; using Index = typename Interaction::Index; - SPIRIT_HOSTDEVICE scalar operator()( const Index & index, const Vector3 * spins ) const; + SPIRIT_HOSTDEVICE scalar operator()( const Index & index, quantity state ) const; using DataRef::DataRef; }; @@ -163,7 +164,7 @@ struct Gradient_Functor : public DataRef using Cache = typename Interaction::Cache; using Index = typename Interaction::Index; - SPIRIT_HOSTDEVICE Vector3 operator()( const Index & index, const Vector3 * spins ) const; + SPIRIT_HOSTDEVICE Vector3 operator()( const Index & index, quantity state ) const; using DataRef::DataRef; }; @@ -177,7 +178,7 @@ struct Hessian_Functor : public DataRef using Index = typename Interaction::Index; template - void operator()( const Index & index, const vectorfield & spins, Callable & hessian ) const; + void operator()( const Index & index, const StateType & state, Callable & hessian ) const; using DataRef::DataRef; }; diff --git a/core/include/engine/spin/interaction/Gaussian.hpp b/core/include/engine/spin/interaction/Gaussian.hpp index f963b1b49..9046491e8 100644 --- a/core/include/engine/spin/interaction/Gaussian.hpp +++ b/core/include/engine/spin/interaction/Gaussian.hpp @@ -2,6 +2,7 @@ #ifndef SPIRIT_CORE_ENGINE_INTERACTION_GAUSSIAN_HPP #define SPIRIT_CORE_ENGINE_INTERACTION_GAUSSIAN_HPP +#include #include #include @@ -23,7 +24,7 @@ E = sum_i^N a_i exp( -l_i^2(m)/(2sigma_i^2) ) where l_i(m) is the distance of m */ struct Gaussian { - using state_t = vectorfield; + using state_t = StateType; struct Data { @@ -113,7 +114,7 @@ struct Functor::Local::DataRef }; template<> -inline scalar Gaussian::Energy::operator()( const Index & index, const Vector3 * spins ) const +inline scalar Gaussian::Energy::operator()( const Index & index, quantity state ) const { scalar result = 0; @@ -125,14 +126,14 @@ inline scalar Gaussian::Energy::operator()( const Index & index, const Vector3 * for( unsigned int igauss = 0; igauss < n_gaussians; ++igauss ) { // Distance between spin and gaussian center - scalar l = 1 - center[igauss].dot( spins[ispin] ); + scalar l = 1 - center[igauss].dot( state.spin[ispin] ); result += amplitude[igauss] * std::exp( -l * l / ( 2.0 * width[igauss] * width[igauss] ) ); }; return result; } template<> -inline Vector3 Gaussian::Gradient::operator()( const Index & index, const Vector3 * spins ) const +inline Vector3 Gaussian::Gradient::operator()( const Index & index, quantity state ) const { Vector3 result = Vector3::Zero(); @@ -144,7 +145,7 @@ inline Vector3 Gaussian::Gradient::operator()( const Index & index, const Vector for( unsigned int i = 0; i < n_gaussians; ++i ) { // Scalar product of spin and gaussian center - scalar l = 1 - center[i].dot( spins[ispin] ); + scalar l = 1 - center[i].dot( state.spin[ispin] ); // Prefactor scalar prefactor = amplitude[i] * std::exp( -l * l / ( 2.0 * width[i] * width[i] ) ) * l / ( width[i] * width[i] ); @@ -156,7 +157,7 @@ inline Vector3 Gaussian::Gradient::operator()( const Index & index, const Vector template<> template -void Gaussian::Hessian::operator()( const Index & index, const vectorfield & spins, Callable & hessian ) const +void Gaussian::Hessian::operator()( const Index & index, const StateType & state, Callable & hessian ) const { if( !is_contributing || index == nullptr ) return; @@ -166,7 +167,7 @@ void Gaussian::Hessian::operator()( const Index & index, const vectorfield & spi for( unsigned int igauss = 0; igauss < n_gaussians; ++igauss ) { // Distance between spin and gaussian center - scalar l = 1 - center[igauss].dot( spins[ispin] ); + scalar l = 1 - center[igauss].dot( state.spin[ispin] ); // Prefactor for all alpha, beta scalar prefactor = amplitude[igauss] * std::exp( -std::pow( l, 2 ) / ( 2.0 * std::pow( width[igauss], 2 ) ) ) / std::pow( width[igauss], 2 ) * ( std::pow( l, 2 ) / std::pow( width[igauss], 2 ) - 1 ); diff --git a/core/include/engine/spin/interaction/Quadruplet.hpp b/core/include/engine/spin/interaction/Quadruplet.hpp index 1639dcfe5..cfa87efea 100644 --- a/core/include/engine/spin/interaction/Quadruplet.hpp +++ b/core/include/engine/spin/interaction/Quadruplet.hpp @@ -4,6 +4,7 @@ #include #include +#include #include namespace Engine @@ -17,7 +18,7 @@ namespace Interaction struct Quadruplet { - using state_t = vectorfield; + using state_t = StateType; struct Data { @@ -142,53 +143,55 @@ struct Functor::Local::DataRef }; template<> -inline scalar Quadruplet::Energy::operator()( const Index & index, const Vector3 * spins ) const +inline scalar Quadruplet::Energy::operator()( const Index & index, quantity state ) const { // don't need to check for `is_contributing` here, because the `transform_reduce` will short circuit correctly return Backend::transform_reduce( index.begin(), index.end(), scalar( 0.0 ), Backend::plus{}, - [this, spins] SPIRIT_LAMBDA( const Quadruplet::IndexType & idx ) -> scalar + [this, state] SPIRIT_LAMBDA( const Quadruplet::IndexType & idx ) -> scalar { const auto & [ispin, jspin, kspin, lspin, iquad] = idx; - return -0.25 * magnitudes[iquad] * ( spins[ispin].dot( spins[jspin] ) ) - * ( spins[kspin].dot( spins[lspin] ) ); + return -0.25 * magnitudes[iquad] * ( state.spin[ispin].dot( state.spin[jspin] ) ) + * ( state.spin[kspin].dot( state.spin[lspin] ) ); } ); } template<> -inline Vector3 Quadruplet::Gradient::operator()( const Index & index, const Vector3 * spins ) const +inline Vector3 Quadruplet::Gradient::operator()( const Index & index, quantity state ) const { // don't need to check for `is_contributing` here, because the `transform_reduce` will short circuit correctly return Backend::transform_reduce( index.begin(), index.end(), Vector3{ 0.0, 0.0, 0.0 }, Backend::plus{}, - [this, spins] SPIRIT_LAMBDA( const Quadruplet::IndexType & idx ) -> Vector3 + [this, state] SPIRIT_LAMBDA( const Quadruplet::IndexType & idx ) -> Vector3 { const auto & [ispin, jspin, kspin, lspin, iquad] = idx; - return spins[jspin] * ( -magnitudes[iquad] * ( spins[kspin].dot( spins[lspin] ) ) ); + return state.spin[jspin] * ( -magnitudes[iquad] * ( state.spin[kspin].dot( state.spin[lspin] ) ) ); } ); } template<> template -void Quadruplet::Hessian::operator()( const Index & index, const vectorfield & spins, Callable & hessian ) const +void Quadruplet::Hessian::operator()( const Index & index, const StateType & state, Callable & hessian ) const { Backend::cpu::for_each( index.begin(), index.end(), - [this, &index, &spins, &hessian]( const Quadruplet::IndexType & idx ) + [this, &index, &state, &hessian]( const Quadruplet::IndexType & idx ) { const auto & [ispin, jspin, kspin, lspin, iquad] = idx; for( int alpha = 0; alpha < 3; ++alpha ) { - hessian( 3 * ispin + alpha, 3 * jspin + alpha, -magnitudes[iquad] * spins[kspin].dot( spins[lspin] ) ); + hessian( + 3 * ispin + alpha, 3 * jspin + alpha, + -magnitudes[iquad] * state.spin[kspin].dot( state.spin[lspin] ) ); for( int beta = 0; beta < 3; ++beta ) { hessian( 3 * ispin + alpha, 3 * kspin + beta, - -magnitudes[iquad] * spins[jspin][alpha] * spins[lspin][beta] ); + -magnitudes[iquad] * state.spin[jspin][alpha] * state.spin[lspin][beta] ); hessian( 3 * ispin + alpha, 3 * lspin + beta, - -magnitudes[iquad] * spins[jspin][alpha] * spins[kspin][beta] ); + -magnitudes[iquad] * state.spin[jspin][alpha] * state.spin[kspin][beta] ); } } } ); diff --git a/core/include/engine/spin/interaction/Zeeman.hpp b/core/include/engine/spin/interaction/Zeeman.hpp index a110adde8..ac5a8a627 100644 --- a/core/include/engine/spin/interaction/Zeeman.hpp +++ b/core/include/engine/spin/interaction/Zeeman.hpp @@ -3,6 +3,7 @@ #define SPIRIT_CORE_ENGINE_INTERACTION_ZEEMANN_HPP #include +#include #include namespace Engine @@ -16,7 +17,7 @@ namespace Interaction struct Zeeman { - using state_t = vectorfield; + using state_t = StateType; struct Data { @@ -109,19 +110,19 @@ struct Functor::Local::DataRef }; template<> -inline scalar Zeeman::Energy::operator()( const Index & index, const Vector3 * spins ) const +inline scalar Zeeman::Energy::operator()( const Index & index, quantity state ) const { if( is_contributing && index != nullptr && *index >= 0 ) { const auto & ispin = *index; - return -mu_s[ispin] * external_field_magnitude * external_field_normal.dot( spins[ispin] ); + return -mu_s[ispin] * external_field_magnitude * external_field_normal.dot( state.spin[ispin] ); } else return 0; } template<> -inline Vector3 Zeeman::Gradient::operator()( const Index & index, const Vector3 * ) const +inline Vector3 Zeeman::Gradient::operator()( const Index & index, quantity ) const { if( is_contributing && index != nullptr && *index >= 0 ) { @@ -134,7 +135,7 @@ inline Vector3 Zeeman::Gradient::operator()( const Index & index, const Vector3 template<> template -void Zeeman::Hessian::operator()( const Index &, const vectorfield &, Callable & ) const {}; +void Zeeman::Hessian::operator()( const Index &, const StateType &, Callable & ) const {}; } // namespace Interaction diff --git a/core/include/io/Dataparser.inl b/core/include/io/Dataparser.inl index ac94bcb2f..3fbdab299 100644 --- a/core/include/io/Dataparser.inl +++ b/core/include/io/Dataparser.inl @@ -1,5 +1,6 @@ #pragma once +#include #include namespace IO @@ -18,7 +19,8 @@ class Buffer using pointer = std::conditional_t, const T *, T *>; public: - explicit Buffer( T & state ) : m_state( &state ) {}; + explicit Buffer( StateType & state ) : m_state( &state.spin ) {}; + explicit Buffer( const StateType & state ) : m_state( &state.spin ) {}; scalar * data() { @@ -57,6 +59,9 @@ private: pointer m_state; }; +Buffer( StateType & ) -> Buffer; +Buffer( const StateType & ) -> Buffer; + } // namespace State } // namespace Spin diff --git a/core/src/Spirit/Chain.cpp b/core/src/Spirit/Chain.cpp index 6e26b8e2d..65ebad5c2 100644 --- a/core/src/Spirit/Chain.cpp +++ b/core/src/Spirit/Chain.cpp @@ -682,9 +682,9 @@ try { chain->images[i]->UpdateEnergy(); if( i > 0 ) - chain->Rx[i] - = chain->Rx[i - 1] - + Engine::Manifoldmath::dist_geodesic( *chain->images[i - 1]->state, *chain->images[i]->state ); + chain->Rx[i] = chain->Rx[i - 1] + + Engine::Manifoldmath::dist_geodesic( + chain->images[i - 1]->state->spin, chain->images[i]->state->spin ); } catch( ... ) { diff --git a/core/src/Spirit/Configurations.cpp b/core/src/Spirit/Configurations.cpp index 6a0973dd4..48d4ed301 100644 --- a/core/src/Spirit/Configurations.cpp +++ b/core/src/Spirit/Configurations.cpp @@ -120,7 +120,7 @@ try // Fetch correct indices and pointers auto [image, chain] = from_indices( state, idx_image, idx_chain ); - state->clipboard_spins = std::make_shared( get( *image->state ) ); + state->clipboard_spins = std::make_shared( image->state->spin ); Log( Utility::Log_Level::Info, Utility::Log_Sender::API, "Copied spin configuration to clipboard.", idx_image, idx_chain ); } @@ -149,8 +149,8 @@ try // Apply configuration image->lock(); - Utility::Configurations::Insert( get( *image->state ), geometry, *state->clipboard_spins, 0, filter ); - geometry.Apply_Pinning( get( *image->state ) ); + Utility::Configurations::Insert( image->state->spin, geometry, *state->clipboard_spins, 0, filter ); + geometry.Apply_Pinning( image->state->spin ); image->unlock(); auto filterstring = filter_to_string( position, r_cut_rectangular, r_cut_cylindrical, r_cut_spherical, inverted ); @@ -197,9 +197,8 @@ try auto filter = get_filter( vpos, r_cut_rectangular, r_cut_cylindrical, r_cut_spherical, inverted ); image->lock(); - Utility::Configurations::Insert( - get( *image->state ), geometry, *state->clipboard_spins, delta, filter ); - geometry.Apply_Pinning( get( *image->state ) ); + Utility::Configurations::Insert( image->state->spin, geometry, *state->clipboard_spins, delta, filter ); + geometry.Apply_Pinning( image->state->spin ); image->unlock(); auto filterstring @@ -243,8 +242,8 @@ try // Apply configuration Vector3 vdir{ direction[0], direction[1], direction[2] }; image->lock(); - Utility::Configurations::Domain( get( *image->state ), geometry, vdir, filter ); - geometry.Apply_Pinning( get( *image->state ) ); + Utility::Configurations::Domain( image->state->spin, geometry, vdir, filter ); + geometry.Apply_Pinning( image->state->spin ); image->unlock(); auto filterstring = filter_to_string( position, r_cut_rectangular, r_cut_cylindrical, r_cut_spherical, inverted ); @@ -299,8 +298,8 @@ try // Apply configuration Vector3 vdir{ 0, 0, 1 }; image->lock(); - Utility::Configurations::Domain( get( *image->state ), geometry, vdir, filter ); - geometry.Apply_Pinning( get( *image->state ) ); + Utility::Configurations::Domain( image->state->spin, geometry, vdir, filter ); + geometry.Apply_Pinning( image->state->spin ); image->unlock(); auto filterstring = filter_to_string( position, r_cut_rectangular, r_cut_cylindrical, r_cut_spherical, inverted ); @@ -333,8 +332,8 @@ try // Apply configuration Vector3 vdir{ 0, 0, -1 }; image->lock(); - Utility::Configurations::Domain( get( *image->state ), geometry, vdir, filter ); - geometry.Apply_Pinning( get( *image->state ) ); + Utility::Configurations::Domain( image->state->spin, geometry, vdir, filter ); + geometry.Apply_Pinning( image->state->spin ); image->unlock(); auto filterstring = filter_to_string( position, r_cut_rectangular, r_cut_cylindrical, r_cut_spherical, inverted ); @@ -378,9 +377,9 @@ try { prng = &image->llg_parameters->prng; } - Utility::Configurations::Random_Sphere( get( *image->state ), geometry, *prng, filter ); + Utility::Configurations::Random_Sphere( image->state->spin, geometry, *prng, filter ); } - geometry.Apply_Pinning( get( *image->state ) ); + geometry.Apply_Pinning( image->state->spin ); image->unlock(); auto filterstring = filter_to_string( position, r_cut_rectangular, r_cut_cylindrical, r_cut_spherical, inverted ); @@ -413,9 +412,8 @@ try // Apply configuration image->lock(); std::mt19937 prng{ 123456789 }; - Utility::Configurations::Add_Noise_Temperature_Sphere( - get( *image->state ), geometry, temperature, prng, filter ); - geometry.Apply_Pinning( get( *image->state ) ); + Utility::Configurations::Add_Noise_Temperature_Sphere( image->state->spin, geometry, temperature, prng, filter ); + geometry.Apply_Pinning( image->state->spin ); image->unlock(); auto filterstring = filter_to_string( position, r_cut_rectangular, r_cut_cylindrical, r_cut_spherical, inverted ); @@ -453,7 +451,7 @@ try { image->lock(); - auto & spins = get( *image->state ); + auto & spins = image->state->spin; auto & mode = *image->modes[idx_mode]; int nos = spins.size(); @@ -508,8 +506,8 @@ try // Apply configuration image->lock(); Utility::Configurations::Hopfion( - get( *image->state ), geometry, vpos, r, order, { normal[0], normal[1], normal[2] }, filter ); - geometry.Apply_Pinning( get( *image->state ) ); + image->state->spin, geometry, vpos, r, order, { normal[0], normal[1], normal[2] }, filter ); + geometry.Apply_Pinning( image->state->spin ); image->unlock(); auto filterstring = filter_to_string( position, r_cut_rectangular, r_cut_cylindrical, r_cut_spherical, inverted ); @@ -550,8 +548,8 @@ try // Apply configuration image->lock(); Utility::Configurations::Skyrmion( - get( *image->state ), geometry, vpos, r, order, phase, upDown, achiral, rl, false, filter ); - geometry.Apply_Pinning( get( *image->state ) ); + image->state->spin, geometry, vpos, r, order, phase, upDown, achiral, rl, false, filter ); + geometry.Apply_Pinning( image->state->spin ); image->unlock(); auto filterstring = filter_to_string( position, r_cut_rectangular, r_cut_cylindrical, r_cut_spherical, inverted ); @@ -600,9 +598,8 @@ try // Apply configuration image->lock(); Utility::Configurations::DW_Skyrmion( - get( *image->state ), geometry, vpos, dw_radius, dw_width, order, phase, upDown, achiral, rl, - filter ); - geometry.Apply_Pinning( get( *image->state ) ); + image->state->spin, geometry, vpos, dw_radius, dw_width, order, phase, upDown, achiral, rl, filter ); + geometry.Apply_Pinning( image->state->spin ); image->unlock(); auto filterstring = filter_to_string( position, r_cut_rectangular, r_cut_cylindrical, r_cut_spherical, inverted ); @@ -653,9 +650,8 @@ try Vector3 vq{ q[0], q[1], q[2] }; Vector3 vaxis{ axis[0], axis[1], axis[2] }; image->lock(); - Utility::Configurations::SpinSpiral( - get( *image->state ), geometry, dir_type, vq, vaxis, theta, filter ); - geometry.Apply_Pinning( get( *image->state ) ); + Utility::Configurations::SpinSpiral( image->state->spin, geometry, dir_type, vq, vaxis, theta, filter ); + geometry.Apply_Pinning( image->state->spin ); image->unlock(); auto filterstring = filter_to_string( position, r_cut_rectangular, r_cut_cylindrical, r_cut_spherical, inverted ); @@ -699,8 +695,7 @@ try Vector3 vq2{ q2[0], q2[1], q2[2] }; Vector3 vaxis{ axis[0], axis[1], axis[2] }; image->lock(); - Utility::Configurations::SpinSpiral( - get( *image->state ), geometry, dir_type, vq1, vq2, vaxis, theta, filter ); + Utility::Configurations::SpinSpiral( image->state->spin, geometry, dir_type, vq1, vq2, vaxis, theta, filter ); image->unlock(); auto filterstring = filter_to_string( position, r_cut_rectangular, r_cut_cylindrical, r_cut_spherical, inverted ); @@ -738,7 +733,7 @@ try // Apply configuration image->lock(); auto geometry = image->hamiltonian->get_geometry(); - Utility::Configurations::Set_Pinned( geometry, get( *image->state ), pinned, filter ); + Utility::Configurations::Set_Pinned( geometry, image->state->spin, pinned, filter ); image->hamiltonian->set_geometry( geometry ); image->unlock(); diff --git a/core/src/Spirit/Geometry.cpp b/core/src/Spirit/Geometry.cpp index 9cbd04b12..dcaeb86ed 100644 --- a/core/src/Spirit/Geometry.cpp +++ b/core/src/Spirit/Geometry.cpp @@ -42,7 +42,7 @@ void Helper_System_Set_Geometry( State::system_t & system, const Data::Geometry if( !same_size( old_geometry, new_geometry ) ) { - Helper_Change_Dimensions( get( *system.state ), old_geometry, new_geometry, { 0, 0, 1 } ); + Helper_Change_Dimensions( system.state->spin, old_geometry, new_geometry, { 0, 0, 1 } ); Helper_Change_Dimensions( system.M.effective_field, old_geometry, new_geometry, { 0, 0, 0 } ); } // Update the system geometry diff --git a/core/src/Spirit/IO.cpp b/core/src/Spirit/IO.cpp index bc996bd13..62d338226 100644 --- a/core/src/Spirit/IO.cpp +++ b/core/src/Spirit/IO.cpp @@ -256,11 +256,9 @@ try } // Read data - file.read_segment_data( idx_image_infile, segment, system_state[0].data() ); - using Engine::Field; - using Engine::get; + file.read_segment_data( idx_image_infile, segment, system_state.spin[0].data() ); - vectorfield & spins = get( system_state ); + auto & spins = system_state.spin; for( std::size_t ispin = 0; ispin < spins.size(); ++ispin ) { if( spins[ispin].norm() < 1e-5 ) @@ -564,8 +562,8 @@ try } // Read data - file.read_segment_data( start_image_infile, segment, system_state[0].data() ); - auto & spins = Engine::get( system_state ); + auto & spins = system_state.spin; + file.read_segment_data( start_image_infile, segment, spins[0].data() ); for( unsigned int ispin = 0; ispin < spins.size(); ++ispin ) { if( spins[ispin].norm() < 1e-5 ) @@ -1012,7 +1010,7 @@ try { const std::string extension = Get_Extension( filename ); - auto & spins = *image->state; + auto & system_state = *image->state; // Open auto file = IO::OVF_File( filename, true ); @@ -1065,7 +1063,7 @@ try { // If the mode buffer is created by resizing then it needs to be allocated if( !image->modes[idx].has_value() ) - image->modes[idx].emplace( vectorfield( spins.size(), Vector3{ 1, 0, 0 } ) ); + image->modes[idx].emplace( vectorfield( system_state.spin.size(), Vector3{ 1, 0, 0 } ) ); // Read header file.read_segment_header( idx + 1, segment ); diff --git a/core/src/Spirit/Quantities.cpp b/core/src/Spirit/Quantities.cpp index ea1577081..e107a7c35 100644 --- a/core/src/Spirit/Quantities.cpp +++ b/core/src/Spirit/Quantities.cpp @@ -29,7 +29,7 @@ try // image->lock(); // Mutex locks in these functions may cause problems with the performance of UIs - auto mean = Engine::Vectormath::mean( get( *image->state ) ); + auto mean = Engine::Vectormath::mean( image->state->spin ); for( int i = 0; i < 3; ++i ) s[i] = mean[i]; @@ -49,9 +49,8 @@ try // image->lock(); // Mutex locks in these functions may cause problems with the performance of UIs - const auto mag = Engine::Vectormath::Magnetization( - get( *image->state ), image->hamiltonian->get_geometry().mu_s ); - image->M.mean = mag; + const auto mag = Engine::Vectormath::Magnetization( image->state->spin, image->hamiltonian->get_geometry().mu_s ); + image->M.mean = mag; // image->unlock(); @@ -76,8 +75,7 @@ try int dimensionality = Geometry_Get_Dimensionality( state, idx_image, idx_chain ); if( dimensionality == 2 ) charge = Engine::Vectormath::TopologicalCharge( - get( *image->state ), image->hamiltonian->get_geometry(), - image->hamiltonian->get_boundary_conditions() ); + image->state->spin, image->hamiltonian->get_geometry(), image->hamiltonian->get_boundary_conditions() ); // image->unlock(); @@ -106,8 +104,8 @@ try if( dimensionality == 2 ) { Engine::Vectormath::TopologicalChargeDensity( - get( *image->state ), image->hamiltonian->get_geometry(), - image->hamiltonian->get_boundary_conditions(), charge_density, triangle_indices ); + image->state->spin, image->hamiltonian->get_geometry(), image->hamiltonian->get_boundary_conditions(), + charge_density, triangle_indices ); } if( charge_density_ptr != nullptr && triangle_indices_ptr != nullptr ) @@ -273,7 +271,7 @@ try int mode_positive = 0; mode_positive = std::max( 0, std::min( n_modes - 1, mode_positive ) ); - Eigen::Ref image_3N = Eigen::Map( image[0].data(), 3 * nos ); + Eigen::Ref image_3N = Eigen::Map( image.spin[0].data(), 3 * nos ); Eigen::Ref grad_3N = Eigen::Map( grad[0].data(), 3 * nos ); // The gradient (unprojected) @@ -305,7 +303,7 @@ try VectorX eigenvalues; MatrixX eigenvectors; bool successful = Engine::Spin::Eigenmodes::Hessian_Partial_Spectrum( - geometry, image, grad, hess, n_modes, basis_3Nx2N, hessian_final, eigenvalues, eigenvectors ); + geometry, image.spin, grad, hess, n_modes, basis_3Nx2N, hessian_final, eigenvalues, eigenvectors ); if( successful ) { @@ -333,7 +331,7 @@ try scalar mode_grad_angle = std::abs( mode_grad / ( mode_3N.norm() * grad_3N.norm() ) ); // Make sure there is nothing wrong - check_modes( image, grad, basis_3Nx2N, eigenvalues, eigenvectors, minimum_mode ); + check_modes( image.spin, grad, basis_3Nx2N, eigenvalues, eigenvectors, minimum_mode ); // If the lowest eigenvalue is negative, we follow the minimum mode if( eigenvalues[0] < -1e-6 && mode_grad_angle > 1e-8 ) // -1e-6)// || switched2) diff --git a/core/src/Spirit/State.cpp b/core/src/Spirit/State.cpp index 4a6ea996b..b8fb14df6 100644 --- a/core/src/Spirit/State.cpp +++ b/core/src/Spirit/State.cpp @@ -136,7 +136,7 @@ try state->active_image = IO::Spin_System_from_Config( state->config_file ); auto & image = state->active_image; Configurations::Random_Sphere( - get( *image->state ), image->hamiltonian->get_geometry(), image->llg_parameters->prng ); + image->state->spin, image->hamiltonian->get_geometry(), image->llg_parameters->prng ); //------------------------------------------------------------------------------------------ //----------------------- Initialize spin system chain ------------------------------------- diff --git a/core/src/Spirit/System.cpp b/core/src/Spirit/System.cpp index 0fc686e95..dfdd41d59 100644 --- a/core/src/Spirit/System.cpp +++ b/core/src/Spirit/System.cpp @@ -43,7 +43,7 @@ try // Fetch correct indices and pointers auto [image, chain] = from_indices( state, idx_image, idx_chain ); - return (scalar *)get( *image->state )[0].data(); + return (scalar *)image->state->spin[0].data(); } catch( ... ) { diff --git a/core/src/Spirit/Transitions.cpp b/core/src/Spirit/Transitions.cpp index c1c574207..7e848b732 100644 --- a/core/src/Spirit/Transitions.cpp +++ b/core/src/Spirit/Transitions.cpp @@ -41,8 +41,7 @@ try { Utility::Configuration_Chain::Homogeneous_Rotation( *chain, idx_1, idx_2 ); for( int img = 0; img < chain->noi; ++img ) - chain->images[img]->hamiltonian->get_geometry().Apply_Pinning( - get( *chain->images[img]->state ) ); + chain->images[img]->hamiltonian->get_geometry().Apply_Pinning( chain->images[img]->state->spin ); Log( Utility::Log_Level::Info, Utility::Log_Sender::API, fmt::format( "Set homogeneous transition between images {} and {}", idx_1 + 1, idx_2 + 1 ), -1, @@ -105,8 +104,7 @@ try { Utility::Configuration_Chain::Add_Noise_Temperature( *chain, idx_1, idx_2, temperature ); for( int img = 0; img < chain->noi; ++img ) - chain->images[img]->hamiltonian->get_geometry().Apply_Pinning( - get( *chain->images[img]->state ) ); + chain->images[img]->hamiltonian->get_geometry().Apply_Pinning( chain->images[img]->state->spin ); Log( Utility::Log_Level::Info, Utility::Log_Sender::API, fmt::format( "Added noise with temperature T={} to images {} - {}", temperature, idx_1 + 1, idx_2 + 1 ), diff --git a/core/src/engine/Manifoldmath.cpp b/core/src/engine/Manifoldmath.cpp index fe4aaf2dd..6aee70978 100644 --- a/core/src/engine/Manifoldmath.cpp +++ b/core/src/engine/Manifoldmath.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include #include @@ -111,20 +112,21 @@ void Geodesic_Tangent( /* Calculates the 'tangent' vectors, i.e.in crudest approximation the difference between an image and the neighbouring */ +template void Tangents( - std::vector> configurations, const std::vector & energies, + const std::vector> & configurations, const std::vector & energies, std::vector & tangents ) { const auto noi = configurations.size(); - const auto nos = ( *configurations[0] ).size(); + const auto nos = configurations[0]->spin.size(); if( noi < 2 ) return; // first image { - const auto & image = *configurations[0]; - const auto & image_plus = *configurations[1]; + const auto & image = configurations[0]->spin; + const auto & image_plus = configurations[1]->spin; Geodesic_Tangent( tangents[0], image, image_plus, image ); // Use the accurate tangent at the endpoints, useful for the dimer method @@ -133,9 +135,9 @@ void Tangents( // Images Inbetween for( unsigned int idx_img = 1; idx_img < noi - 1; ++idx_img ) { - const auto & image = *configurations[idx_img]; - const auto & image_plus = *configurations[idx_img + 1]; - const auto & image_minus = *configurations[idx_img - 1]; + const auto & image = configurations[idx_img]->spin; + const auto & image_plus = configurations[idx_img + 1]->spin; + const auto & image_minus = configurations[idx_img - 1]->spin; // Energies scalar E_mid = 0, E_plus = 0, E_minus = 0; @@ -200,14 +202,18 @@ void Tangents( // Last Image { - const auto & image = *configurations[noi - 1]; - const auto & image_minus = *configurations[noi - 2]; + const auto & image = configurations[noi - 1]->spin; + const auto & image_minus = configurations[noi - 2]->spin; Geodesic_Tangent( tangents[noi - 1], image_minus, image, image ); // Use the accurate tangent at the endpoints, useful for the dimer method } } // end Tangents +template void Tangents( + const std::vector> & configurations, const std::vector & energies, + std::vector & tangents ); + scalar norm( const vectorfield & vf ) { scalar x = Vectormath::dot( vf, vf ); diff --git a/core/src/engine/spin/Eigenmodes.cpp b/core/src/engine/spin/Eigenmodes.cpp index d40f9b578..d0976271f 100644 --- a/core/src/engine/spin/Eigenmodes.cpp +++ b/core/src/engine/spin/Eigenmodes.cpp @@ -68,7 +68,7 @@ void Calculate_Eigenmodes( system_t & system, int idx_img, int idx_chain ) auto & n_modes = system.ema_parameters->n_modes; // vectorfield mode(nos, Vector3{1, 0, 0}); - auto spins_initial = *system.state; + auto state_initial = *system.state; Log( Log_Level::Info, Log_Sender::EMA, fmt::format( "Started calculation of {} Eigenmodes ", n_modes ), idx_img, idx_chain ); @@ -77,7 +77,7 @@ void Calculate_Eigenmodes( system_t & system, int idx_img, int idx_chain ) vectorfield gradient( nos ); // The gradient (unprojected) - system.hamiltonian->Gradient( spins_initial, gradient ); + system.hamiltonian->Gradient( state_initial, gradient ); // auto mask = system.geometry->mask_unpinned.data(); // auto g = gradient.data(); // Backend::for_each_n( SPIRIT_PAR Backend::make_counting_iterator( 0 ), gradient.size(), [g, mask] @@ -96,25 +96,25 @@ void Calculate_Eigenmodes( system_t & system, int idx_img, int idx_chain ) { // The Hessian (unprojected) SpMatrixX hessian( 3 * nos, 3 * nos ); - system.hamiltonian->Sparse_Hessian( spins_initial, hessian ); + system.hamiltonian->Sparse_Hessian( state_initial, hessian ); // Get the eigenspectrum SpMatrixX hessian_constrained = SpMatrixX( 2 * nos, 2 * nos ); successful = Eigenmodes::Sparse_Hessian_Partial_Spectrum( - system.hamiltonian->get_geometry(), spins_initial, gradient, hessian, n_modes, tangent_basis, + system.hamiltonian->get_geometry(), state_initial.spin, gradient, hessian, n_modes, tangent_basis, hessian_constrained, eigenvalues, eigenvectors ); } else { // The Hessian (unprojected) MatrixX hessian( 3 * nos, 3 * nos ); - system.hamiltonian->Hessian( spins_initial, hessian ); + system.hamiltonian->Hessian( state_initial, hessian ); // Get the eigenspectrum MatrixX hessian_constrained = MatrixX::Zero( 2 * nos, 2 * nos ); MatrixX _tangent_basis = MatrixX( tangent_basis ); successful = Eigenmodes::Hessian_Partial_Spectrum( - system.hamiltonian->get_geometry(), spins_initial, gradient, hessian, n_modes, _tangent_basis, + system.hamiltonian->get_geometry(), state_initial.spin, gradient, hessian, n_modes, _tangent_basis, hessian_constrained, eigenvalues, eigenvectors ); tangent_basis = _tangent_basis.sparseView(); diff --git a/core/src/engine/spin/HTST.cpp b/core/src/engine/spin/HTST.cpp index 4500262d5..350953e8f 100644 --- a/core/src/engine/spin/HTST.cpp +++ b/core/src/engine/spin/HTST.cpp @@ -43,7 +43,7 @@ void Calculate( Data::HTST_Info & htst_info, int n_eigenmodes_keep ) auto & image_minimum = *htst_info.minimum->state; auto & image_sp = *htst_info.saddle_point->state; - int nos = image_minimum.size(); + int nos = image_minimum.spin.size(); if( n_eigenmodes_keep < 0 ) n_eigenmodes_keep = 2 * nos; @@ -68,7 +68,7 @@ void Calculate( Data::HTST_Info & htst_info, int n_eigenmodes_keep ) Log( Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Checking if initial configuration is an extremum..." ); Vectormath::set_c_a( 1, gradient_minimum, force_tmp ); - Manifoldmath::project_tangential( force_tmp, image_minimum ); + Manifoldmath::project_tangential( force_tmp, image_minimum.spin ); scalar fmax_minimum = Vectormath::max_norm( force_tmp ); if( fmax_minimum > epsilon_force ) { @@ -90,7 +90,7 @@ void Calculate( Data::HTST_Info & htst_info, int n_eigenmodes_keep ) Log( Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Checking if transition configuration is an extremum..." ); Vectormath::set_c_a( 1, gradient_sp, force_tmp ); - Manifoldmath::project_tangential( force_tmp, image_sp ); + Manifoldmath::project_tangential( force_tmp, image_sp.spin ); scalar fmax_sp = Vectormath::max_norm( force_tmp ); if( fmax_sp > epsilon_force ) { @@ -119,8 +119,8 @@ void Calculate( Data::HTST_Info & htst_info, int n_eigenmodes_keep ) htst_info.eigenvalues_sp = VectorX::Zero( 2 * nos ); htst_info.eigenvectors_sp = MatrixX::Zero( 2 * nos, 2 * nos ); Geodesic_Eigen_Decomposition( - image_sp, gradient_sp, hessian_sp, hessian_geodesic_sp_3N, hessian_geodesic_sp_2N, htst_info.eigenvalues_sp, - htst_info.eigenvectors_sp ); + image_sp.spin, gradient_sp, hessian_sp, hessian_geodesic_sp_3N, hessian_geodesic_sp_2N, + htst_info.eigenvalues_sp, htst_info.eigenvectors_sp ); // Print some eigenvalues block = std::vector{ "10 lowest eigenvalues at saddle point:" }; @@ -163,12 +163,12 @@ void Calculate( Data::HTST_Info & htst_info, int n_eigenmodes_keep ) // Calculation of the 'a' parameters... htst_info.perpendicular_velocity = VectorX::Zero( 2 * nos ); MatrixX basis_sp = MatrixX::Zero( 3 * nos, 2 * nos ); - Manifoldmath::tangent_basis_spherical( image_sp, basis_sp ); + Manifoldmath::tangent_basis_spherical( image_sp.spin, basis_sp ); // Manifoldmath::tangent_basis(image_sp, basis_sp); // Calculate_Perpendicular_Velocity_2N(image_sp, hessian_geodesic_sp_2N, basis_sp, htst_info.eigenvectors_sp, // perpendicular_velocity_sp); Calculate_Perpendicular_Velocity( - image_sp, htst_info.saddle_point->hamiltonian->get_geometry().mu_s, hessian_geodesic_sp_3N, basis_sp, + image_sp.spin, htst_info.saddle_point->hamiltonian->get_geometry().mu_s, hessian_geodesic_sp_3N, basis_sp, htst_info.eigenvectors_sp, htst_info.perpendicular_velocity ); // Reduce the number of saved eigenmodes @@ -223,8 +223,8 @@ void Calculate( Data::HTST_Info & htst_info, int n_eigenmodes_keep ) htst_info.eigenvalues_min = VectorX::Zero( 2 * nos ); htst_info.eigenvectors_min = MatrixX::Zero( 2 * nos, 2 * nos ); Geodesic_Eigen_Decomposition( - image_minimum, gradient_minimum, hessian_minimum, hessian_geodesic_minimum_3N, hessian_geodesic_minimum_2N, - htst_info.eigenvalues_min, htst_info.eigenvectors_min ); + image_minimum.spin, gradient_minimum, hessian_minimum, hessian_geodesic_minimum_3N, + hessian_geodesic_minimum_2N, htst_info.eigenvalues_min, htst_info.eigenvectors_min ); // Print some eigenvalues block = std::vector{ "10 lowest eigenvalues at minimum:" }; @@ -324,7 +324,7 @@ void Calculate( Data::HTST_Info & htst_info, int n_eigenmodes_keep ) scalar Calculate_Zero_Volume( const system_t & system ) { - const auto & spins = *system.state; + const auto & spins = system.state->spin; const auto & geometry = system.hamiltonian->get_geometry(); const int nos = geometry.nos; const auto & bravais_vectors = geometry.bravais_vectors; diff --git a/core/src/engine/spin/Method_EMA.cpp b/core/src/engine/spin/Method_EMA.cpp index 45929d4a2..cb0ffc5ec 100644 --- a/core/src/engine/spin/Method_EMA.cpp +++ b/core/src/engine/spin/Method_EMA.cpp @@ -38,7 +38,7 @@ Method_EMA::Method_EMA( std::shared_ptr system, int idx_img, int idx_c this->angle = scalarfield( this->nos ); this->angle_initial = scalarfield( this->nos ); this->axis = vectorfield( this->nos ); - this->spins_initial = *this->system->state; + this->spins_initial = this->system->state->spin; Eigenmodes::Check_Eigenmode_Parameters( *system ); @@ -64,7 +64,7 @@ void Method_EMA::Iteration() // Reset members this->following_mode = this->parameters_ema->n_mode_follow; // Restore the initial spin configuration - ( *this->system->state ) = this->spins_initial; + this->system->state->spin = this->spins_initial; // Set the new mode this->mode = *this->system->modes[following_mode]; @@ -91,7 +91,7 @@ void Method_EMA::Iteration() Vectormath::scale( this->angle, t_angle ); // Rotate the spins - Vectormath::rotate( this->spins_initial, this->axis, this->angle, image ); + Vectormath::rotate( this->spins_initial, this->axis, this->angle, image.spin ); ++this->counter; std::this_thread::sleep_for( std::chrono::milliseconds( 50 ) ); @@ -109,7 +109,7 @@ void Method_EMA::Finalize() { this->lock(); // The initial spin configuration must be restored - ( *this->system->state ) = this->spins_initial; + this->system->state->spin = this->spins_initial; this->unlock(); } diff --git a/core/src/engine/spin/Method_GNEB.cpp b/core/src/engine/spin/Method_GNEB.cpp index cd23dcd91..e2f62ec7d 100644 --- a/core/src/engine/spin/Method_GNEB.cpp +++ b/core/src/engine/spin/Method_GNEB.cpp @@ -57,7 +57,7 @@ Method_GNEB::Method_GNEB( std::shared_ptr chain, int idx_chain this->max_torque_all = std::vector( this->noi, 0 ); // Create shared pointers to the method's systems' spin configurations - this->configurations = std::vector>( this->noi ); + this->configurations = std::vector>( this->noi ); for( int i = 0; i < this->noi; ++i ) this->configurations[i] = this->systems[i]->state; @@ -82,7 +82,7 @@ std::vector Method_GNEB::getTorqueMaxNorm_All() template void Method_GNEB::Calculate_Force( - const std::vector> & configurations, std::vector & forces ) + const std::vector> & configurations, std::vector & forces ) { // We assume here that we receive a vector of configurations that corresponds to the vector of systems we gave the // Solver. @@ -102,7 +102,7 @@ void Method_GNEB::Calculate_Force( auto * eff_field = this->chain->images[img]->M.effective_field.data(); auto * f_grad = F_gradient[img].data(); Backend::for_each_n( - SPIRIT_PAR Backend::make_counting_iterator( 0 ), image.size(), + SPIRIT_PAR Backend::make_counting_iterator( 0 ), image.spin.size(), [eff_field, f_grad] SPIRIT_LAMBDA( const int idx ) { eff_field[idx] *= -1; @@ -111,7 +111,7 @@ void Method_GNEB::Calculate_Force( if( img > 0 ) { - Rx[img] = Rx[img - 1] + Manifoldmath::dist_geodesic( image, *configurations[img - 1] ); + Rx[img] = Rx[img - 1] + Manifoldmath::dist_geodesic( image.spin, configurations[img - 1]->spin ); if( Rx[img] - Rx[img - 1] < 1e-10 ) { Log( Log_Level::Error, Log_Sender::GNEB, @@ -167,10 +167,10 @@ void Method_GNEB::Calculate_Force( // Get the total force on the image chain // Loop over images to calculate the total force on each Image - std::size_t nos = configurations[0]->size(); + std::size_t nos = configurations[0]->spin.size(); for( int img = 1; img < chain->noi - 1; ++img ) { - auto & image = *configurations[img]; + auto & image = configurations[img]->spin; // The gradient force (unprojected) is simply the effective field // this->chain->images[img]->hamiltonian->Gradient(image, F_gradient[img]); @@ -203,10 +203,10 @@ void Method_GNEB::Calculate_Force( // Calculate finite difference secants vectorfield t_plus( nos ); vectorfield t_minus( nos ); - Vectormath::set_c_a( 1, *this->chain->images[img + 1]->state, t_plus ); - Vectormath::add_c_a( -1, *this->chain->images[img]->state, t_plus ); - Vectormath::set_c_a( 1, *this->chain->images[img]->state, t_minus ); - Vectormath::add_c_a( -1, *this->chain->images[img - 1]->state, t_minus ); + Vectormath::set_c_a( 1, this->chain->images[img + 1]->state->spin, t_plus ); + Vectormath::add_c_a( -1, this->chain->images[img]->state->spin, t_plus ); + Vectormath::set_c_a( 1, this->chain->images[img]->state->spin, t_minus ); + Vectormath::add_c_a( -1, this->chain->images[img - 1]->state->spin, t_minus ); Manifoldmath::normalize( t_plus ); Manifoldmath::normalize( t_minus ); // Get the finite difference (path shrinking) direction @@ -258,16 +258,16 @@ void Method_GNEB::Calculate_Force( if( chain->gneb_parameters->moving_endpoints ) { int noi = chain->noi; - Manifoldmath::project_tangential( F_gradient[0], *configurations[0] ); - Manifoldmath::project_tangential( F_gradient[noi - 1], *configurations[noi - 1] ); + Manifoldmath::project_tangential( F_gradient[0], configurations[0]->spin ); + Manifoldmath::project_tangential( F_gradient[noi - 1], configurations[noi - 1]->spin ); // Overall translational force if( chain->gneb_parameters->translating_endpoints ) { const auto * F_gradient_left = F_gradient[0].data(); const auto * F_gradient_right = F_gradient[noi - 1].data(); - const auto * spins_left = this->chain->images[0]->state->data(); - const auto * spins_right = this->chain->images[noi - 1]->state->data(); + const auto * spins_left = this->chain->images[0]->state->spin.data(); + const auto * spins_right = this->chain->images[noi - 1]->state->spin.data(); auto * F_translation_left_ptr = F_translation_left.data(); auto * F_translation_right_ptr = F_translation_right.data(); // clang-format off @@ -348,7 +348,7 @@ void Method_GNEB::Calculate_Force( template void Method_GNEB::Calculate_Force_Virtual( - const std::vector> & configurations, const std::vector & forces, + const std::vector> & configurations, const std::vector & forces, std::vector & forces_virtual ) { using namespace Utility; @@ -362,7 +362,7 @@ void Method_GNEB::Calculate_Force_Virtual( continue; } - auto & image = *configurations[i]; + auto & image = configurations[i]->spin; const auto & force = forces[i]; auto & force_virtual = forces_virtual[i]; auto & parameters = *this->systems[i]->llg_parameters; @@ -407,7 +407,7 @@ void Method_GNEB::Hook_Post_Iteration() for( int img = 0; img < chain->noi; ++img ) { - Manifoldmath::project_tangential( F_total[img], *( this->systems[img]->state ) ); + Manifoldmath::project_tangential( F_total[img], this->systems[img]->state->spin ); const scalar fmax = Vectormath::max_norm( F_total[img] ); // Set maximum per image this->max_torque_all[img] = fmax; @@ -416,7 +416,7 @@ void Method_GNEB::Hook_Post_Iteration() this->max_torque = fmax; // Set the effective fields - Manifoldmath::project_tangential( this->forces[img], *this->systems[img]->state ); + Manifoldmath::project_tangential( this->forces[img], this->systems[img]->state->spin ); // Vectormath::set_c_a(1, this->forces[img], this->systems[img]->effective_field); } @@ -457,7 +457,7 @@ void Method_GNEB::Calculate_Interpolated_Energy_Contributions() Log( Utility::Log_Level::Debug, Utility::Log_Sender::GNEB, std::string( "Calculating interpolated energy contributions" ), -1, -1 ); - std::size_t nos = this->configurations[0]->size(); + std::size_t nos = this->configurations[0]->spin.size(); int noi = this->chain->noi; if( chain->images[0]->hamiltonian->Name() != "Heisenberg" ) @@ -515,21 +515,24 @@ void Method_GNEB::Finalize() template void Method_GNEB::Message_Block_Start( std::vector & block ) { - scalar length = Manifoldmath::dist_geodesic( *this->configurations[0], *this->configurations[this->noi - 1] ); + scalar length + = Manifoldmath::dist_geodesic( this->configurations[0]->spin, this->configurations[this->noi - 1]->spin ); block.emplace_back( fmt::format( " Total path length: {}", length ) ); } template void Method_GNEB::Message_Block_Step( std::vector & block ) { - scalar length = Manifoldmath::dist_geodesic( *this->configurations[0], *this->configurations[this->noi - 1] ); + scalar length + = Manifoldmath::dist_geodesic( this->configurations[0]->spin, this->configurations[this->noi - 1]->spin ); block.emplace_back( fmt::format( " Total path length: {}", length ) ); } template void Method_GNEB::Message_Block_End( std::vector & block ) { - scalar length = Manifoldmath::dist_geodesic( *this->configurations[0], *this->configurations[this->noi - 1] ); + scalar length + = Manifoldmath::dist_geodesic( this->configurations[0]->spin, this->configurations[this->noi - 1]->spin ); block.emplace_back( fmt::format( " Total path length: {}", length ) ); } @@ -587,21 +590,23 @@ void Method_GNEB::Save_Current( std::string starttime, int iteration, bo segment.title = strdup( title.c_str() ); std::string output_comment = fmt::format( "{}\n# Desc: Image {} of {}", output_comment_base, 0, chain->noi ); - segment.comment = strdup( output_comment.c_str() ); - segment.valuedim = IO::Spin::State::valuedim; - segment.valuelabels = strdup( IO::Spin::State::valuelabels.data() ); - segment.valueunits = strdup( IO::Spin::State::valueunits.data() ); - auto & spins = *this->chain->images[0]->state; - IO::OVF_File( chainFile ).write_segment( segment, spins[0].data(), static_cast( format ) ); + segment.comment = strdup( output_comment.c_str() ); + segment.valuedim = IO::Spin::State::valuedim; + segment.valuelabels = strdup( IO::Spin::State::valuelabels.data() ); + segment.valueunits = strdup( IO::Spin::State::valueunits.data() ); + const auto & system_state = *this->chain->images[0]->state; + IO::OVF_File( chainFile ) + .write_segment( segment, system_state.spin[0].data(), static_cast( format ) ); } // Append all the others for( int i = 1; i < this->chain->noi; i++ ) { - auto & spins = *this->chain->images[i]->state; + const auto & system_state = *this->chain->images[i]->state; std::string output_comment = fmt::format( "{}\n# Desc: Image {} of {}", output_comment_base, i, chain->noi ); segment.comment = strdup( output_comment.c_str() ); - IO::OVF_File( chainFile ).append_segment( segment, spins[0].data(), static_cast( format ) ); + IO::OVF_File( chainFile ) + .append_segment( segment, system_state.spin[0].data(), static_cast( format ) ); } } catch( ... ) diff --git a/core/src/engine/spin/Method_LLG.cpp b/core/src/engine/spin/Method_LLG.cpp index ff0298e50..10ec8b85f 100644 --- a/core/src/engine/spin/Method_LLG.cpp +++ b/core/src/engine/spin/Method_LLG.cpp @@ -42,7 +42,7 @@ Method_LLG::Method_LLG( std::shared_ptr system, int idx_img, i this->max_torque = system->llg_parameters->force_convergence + 1.0; // Create shared pointers to the method's systems' spin configurations - this->configurations = std::vector>( this->noi ); + this->configurations = std::vector>( this->noi ); for( int i = 0; i < this->noi; ++i ) this->configurations[i] = this->systems[i]->state; @@ -72,7 +72,7 @@ void Method_LLG::Prepare_Thermal_Field() template void Method_LLG::Calculate_Force( - const std::vector> & configurations, std::vector & forces ) + const std::vector> & configurations, std::vector & forces ) { // Loop over images to calculate the total force on each Image for( std::size_t img = 0; img < this->systems.size(); ++img ) @@ -92,7 +92,7 @@ void Method_LLG::Calculate_Force( template void Method_LLG::Calculate_Force_Virtual( - const std::vector> & configurations, const std::vector & forces, + const std::vector> & configurations, const std::vector & forces, std::vector & forces_virtual ) { for( int i = 0; i < this->noi; ++i ) @@ -100,7 +100,7 @@ void Method_LLG::Calculate_Force_Virtual( const auto & sys = *this->systems[i]; common_methods[i].Virtual_Force_Spin( *sys.llg_parameters, sys.hamiltonian->get_geometry(), sys.hamiltonian->get_boundary_conditions(), - *configurations[i], forces[i], forces_virtual[i] ); + configurations[i]->spin, forces[i], forces_virtual[i] ); } } @@ -133,7 +133,7 @@ void Method_LLG::Hook_Post_Iteration() for( std::size_t img = 0; img < this->systems.size(); ++img ) { this->force_converged[img] = false; - Manifoldmath::project_tangential( this->forces_virtual[img], *( this->systems[img]->state ) ); + Manifoldmath::project_tangential( this->forces_virtual[img], this->systems[img]->state->spin ); const scalar fmax = Vectormath::max_norm( this->forces_virtual[img] ); if( fmax > 0 ) @@ -153,7 +153,7 @@ void Method_LLG::Hook_Post_Iteration() // ToDo: How to update eff_field without numerical overhead? // systems[0]->effective_field = Gradient[0]; // Vectormath::scale(systems[0]->effective_field, -1); - Manifoldmath::project_tangential( this->forces[0], *this->systems[0]->state ); + Manifoldmath::project_tangential( this->forces[0], this->systems[0]->state->spin ); Vectormath::set_c_a( 1, this->forces[0], this->systems[0]->M.effective_field ); // systems[0]->UpdateEffectiveField(); @@ -263,7 +263,7 @@ void Method_LLG::Save_Current( std::string starttime, int iteration, boo IO::VF_FileFormat format = sys.llg_parameters->output_vf_filetype; // Spin Configuration - auto & spins = *sys.state; + auto & system_state = *sys.state; auto segment = IO::OVF_Segment( sys.hamiltonian->get_geometry() ); std::string title = fmt::format( "SPIRIT Version {}", Utility::version_full ); segment.title = strdup( title.c_str() ); @@ -272,7 +272,7 @@ void Method_LLG::Save_Current( std::string starttime, int iteration, boo segment.valuelabels = strdup( IO::Spin::State::valuelabels.data() ); segment.valueunits = strdup( IO::Spin::State::valueunits.data() ); - const IO::Spin::State::Buffer buffer( spins ); + const IO::Spin::State::Buffer buffer( system_state ); if( append ) IO::OVF_File( spinsFile ).append_segment( segment, buffer.data(), static_cast( format ) ); else diff --git a/core/src/engine/spin/Method_MC.cpp b/core/src/engine/spin/Method_MC.cpp index 0c753407a..983139626 100644 --- a/core/src/engine/spin/Method_MC.cpp +++ b/core/src/engine/spin/Method_MC.cpp @@ -54,20 +54,20 @@ Method_MC::Method_MC( std::shared_ptr system, int idx_img, int idx_cha void Method_MC::Iteration() { // Temporaries - auto & spins_old = *this->system->state; - auto spins_new = spins_old; + auto & state_old = *this->system->state; + auto state_new = state_old; // Generate randomly displaced spin configuration according to cone radius // Vectormath::get_random_vectorfield_unitsphere(this->parameters_mc->prng, random_unit_vectors); // TODO: add switch between Metropolis and heat bath // One Metropolis step - Metropolis( spins_old, spins_new ); - Vectormath::set_c_a( 1, spins_new, spins_old ); + Metropolis( state_old, state_new ); + Vectormath::set_c_a( 1, state_new.spin, state_old.spin ); } // Simple metropolis step -void Method_MC::Metropolis( const vectorfield & spins_old, vectorfield & spins_new ) +void Method_MC::Metropolis( const StateType & state_old, StateType & state_new ) { auto distribution = std::uniform_real_distribution( 0, 1 ); auto distribution_idx = std::uniform_int_distribution<>( 0, this->nos - 1 ); @@ -115,9 +115,9 @@ void Method_MC::Metropolis( const vectorfield & spins_old, vectorfield & spins_n if( this->parameters_mc->metropolis_step_cone ) { // Calculate local basis for the spin - if( spins_old[ispin].z() < 1 - 1e-10 ) + if( state_old.spin[ispin].z() < 1 - 1e-10 ) { - local_basis.col( 2 ) = spins_old[ispin]; + local_basis.col( 2 ) = state_old.spin[ispin]; local_basis.col( 0 ) = ( local_basis.col( 2 ).cross( e_z ) ).normalized(); local_basis.col( 1 ) = local_basis.col( 2 ).cross( local_basis.col( 0 ) ); } @@ -138,7 +138,7 @@ void Method_MC::Metropolis( const vectorfield & spins_old, vectorfield & spins_n Vector3 local_spin_new{ sintheta * std::cos( phi ), sintheta * std::sin( phi ), costheta }; // New spin orientation in regular basis - spins_new[ispin] = local_basis * local_spin_new; + state_new.spin[ispin] = local_basis * local_spin_new; } // Sample the entire unit sphere else @@ -152,12 +152,12 @@ void Method_MC::Metropolis( const vectorfield & spins_old, vectorfield & spins_n phi = 2 * Constants::Pi * distribution( this->parameters_mc->prng ); // New spin orientation in local basis - spins_new[ispin] = Vector3{ sintheta * std::cos( phi ), sintheta * std::sin( phi ), costheta }; + state_new.spin[ispin] = Vector3{ sintheta * std::cos( phi ), sintheta * std::sin( phi ), costheta }; } // Energy difference of configurations with and without displacement - scalar Eold = this->system->hamiltonian->Energy_Single_Spin( ispin, spins_old ); - scalar Enew = this->system->hamiltonian->Energy_Single_Spin( ispin, spins_new ); + scalar Eold = this->system->hamiltonian->Energy_Single_Spin( ispin, state_old ); + scalar Enew = this->system->hamiltonian->Energy_Single_Spin( ispin, state_new ); scalar Ediff = Enew - Eold; // Metropolis criterion: reject the step if energy rose @@ -166,7 +166,7 @@ void Method_MC::Metropolis( const vectorfield & spins_old, vectorfield & spins_n if( this->parameters_mc->temperature < 1e-12 ) { // Restore the spin - spins_new[ispin] = spins_old[ispin]; + state_new.spin[ispin] = state_old.spin[ispin]; // Counter for the number of rejections ++this->n_rejected; } @@ -181,7 +181,7 @@ void Method_MC::Metropolis( const vectorfield & spins_old, vectorfield & spins_n if( exp_ediff < x_metropolis ) { // Restore the spin - spins_new[ispin] = spins_old[ispin]; + state_new.spin[ispin] = state_old.spin[ispin]; // Counter for the number of rejections ++this->n_rejected; } diff --git a/core/src/engine/spin/Method_MMF.cpp b/core/src/engine/spin/Method_MMF.cpp index 7cc0d8f6e..15fefa39c 100644 --- a/core/src/engine/spin/Method_MMF.cpp +++ b/core/src/engine/spin/Method_MMF.cpp @@ -59,7 +59,7 @@ Method_MMF::Method_MMF( std::shared_ptr system, int idx_chain this->mode_follow_previous = 0; // Create shared pointers to the method's systems' spin configurations - this->configurations = std::vector>( 1 ); + this->configurations = std::vector>( 1 ); this->configurations[0] = this->system->state; //---- Initialise Solver-specific variables @@ -68,7 +68,7 @@ Method_MMF::Method_MMF( std::shared_ptr system, int idx_chain template void Method_MMF::Calculate_Force( - const std::vector> & configurations, std::vector & forces ) + const std::vector> & configurations, std::vector & forces ) { // if (this->mm_function == "Spectra Matrix") // { @@ -177,7 +177,7 @@ void check_modes( template void Method_MMF::Calculate_Force_Spectra_Matrix( - const std::vector> & configurations, std::vector & forces ) + const std::vector> & configurations, std::vector & forces ) { auto & image = *configurations[0]; auto & force = forces[0]; @@ -206,7 +206,7 @@ void Method_MMF::Calculate_Force_Spectra_Matrix( // The Hessian (unprojected) hamiltonian.Hessian( image, hessian ); - Eigen::Ref image_3N = Eigen::Map( image[0].data(), 3 * nos ); + Eigen::Ref image_3N = Eigen::Map( image.spin[0].data(), 3 * nos ); Eigen::Ref gradient_3N = Eigen::Map( gradient[0].data(), 3 * nos ); //////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -218,7 +218,7 @@ void Method_MMF::Calculate_Force_Spectra_Matrix( VectorX eigenvalues; MatrixX eigenvectors; bool successful = Eigenmodes::Hessian_Partial_Spectrum( - geometry, image, gradient, hessian, n_modes, basis_3Nx2N, hessian_final, eigenvalues, eigenvectors ); + geometry, image.spin, gradient, hessian, n_modes, basis_3Nx2N, hessian_final, eigenvalues, eigenvectors ); if( successful ) { @@ -289,9 +289,9 @@ void Method_MMF::Calculate_Force_Spectra_Matrix( scalar mode_grad_angle = std::abs( mode_grad / ( mode_3N.norm() * gradient_3N.norm() ) ); // Make sure there is nothing wrong - check_modes( image, gradient, basis_3Nx2N, eigenvalues, eigenvectors, minimum_mode ); + check_modes( image.spin, gradient, basis_3Nx2N, eigenvalues, eigenvectors, minimum_mode ); - Manifoldmath::project_tangential( gradient, image ); + Manifoldmath::project_tangential( gradient, image.spin ); // Some debugging prints if( mode_evalue < -1e-6 && mode_grad_angle > 1e-8 ) // -1e-6)// || switched2) @@ -419,7 +419,7 @@ void Method_MMF::Hook_Post_Iteration() // Loop over images to calculate the maximum torques for( unsigned int img = 0; img < this->systems.size(); ++img ) { - Manifoldmath::project_tangential( this->forces_virtual[img], *( this->systems[img]->state ) ); + Manifoldmath::project_tangential( this->forces_virtual[img], this->systems[img]->state->spin ); const scalar fmax = Vectormath::max_norm( this->forces_virtual[img] ); if( fmax > 0 ) this->max_torque = fmax; @@ -476,16 +476,16 @@ void Method_MMF::Save_Current( std::string starttime, int iteration, boo IO::VF_FileFormat format = sys.mmf_parameters->output_vf_filetype; // Spin Configuration - const auto & spins = *sys.state; - auto segment = IO::OVF_Segment( sys.hamiltonian->get_geometry() ); - std::string title = fmt::format( "SPIRIT Version {}", Utility::version_full ); - segment.title = strdup( title.c_str() ); - segment.comment = strdup( output_comment.c_str() ); - segment.valuedim = IO::Spin::State::valuedim; - segment.valuelabels = strdup( IO::Spin::State::valuelabels.data() ); - segment.valueunits = strdup( IO::Spin::State::valueunits.data() ); - - const IO::Spin::State::Buffer buffer( spins ); + const auto & system_state = *sys.state; + auto segment = IO::OVF_Segment( sys.hamiltonian->get_geometry() ); + std::string title = fmt::format( "SPIRIT Version {}", Utility::version_full ); + segment.title = strdup( title.c_str() ); + segment.comment = strdup( output_comment.c_str() ); + segment.valuedim = IO::Spin::State::valuedim; + segment.valuelabels = strdup( IO::Spin::State::valuelabels.data() ); + segment.valueunits = strdup( IO::Spin::State::valueunits.data() ); + + const IO::Spin::State::Buffer buffer( system_state ); if( append ) IO::OVF_File( spinsFile ).append_segment( segment, buffer.data(), int( format ) ); else diff --git a/core/src/engine/spin/Sparse_HTST.cpp b/core/src/engine/spin/Sparse_HTST.cpp index 8ad37460f..7da50272c 100644 --- a/core/src/engine/spin/Sparse_HTST.cpp +++ b/core/src/engine/spin/Sparse_HTST.cpp @@ -264,7 +264,7 @@ void Calculate( Data::HTST_Info & htst_info ) auto & image_minimum = *htst_info.minimum->state; auto & image_sp = *htst_info.saddle_point->state; - std::size_t nos = image_minimum.size(); + std::size_t nos = image_minimum.spin.size(); Log( Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Saving NO eigenvectors." ); @@ -284,7 +284,7 @@ void Calculate( Data::HTST_Info & htst_info ) Log( Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Checking if initial configuration is an extremum..." ); Vectormath::set_c_a( 1, gradient_minimum, force_tmp ); - Manifoldmath::project_tangential( force_tmp, image_minimum ); + Manifoldmath::project_tangential( force_tmp, image_minimum.spin ); scalar fmax_minimum = Vectormath::max_norm( force_tmp ); if( fmax_minimum > epsilon_force ) { @@ -306,7 +306,7 @@ void Calculate( Data::HTST_Info & htst_info ) Log( Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Checking if transition configuration is an extremum..." ); Vectormath::set_c_a( 1, gradient_sp, force_tmp ); - Manifoldmath::project_tangential( force_tmp, image_sp ); + Manifoldmath::project_tangential( force_tmp, image_sp.spin ); scalar fmax_sp = Vectormath::max_norm( force_tmp ); if( fmax_sp > epsilon_force ) { @@ -327,7 +327,7 @@ void Calculate( Data::HTST_Info & htst_info ) Log( Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Evaluate tangent basis ..." ); SpMatrixX tangent_basis = SpMatrixX( 3 * nos, 2 * nos ); - Manifoldmath::sparse_tangent_basis_spherical( image_sp, tangent_basis ); + Manifoldmath::sparse_tangent_basis_spherical( image_sp.spin, tangent_basis ); // Evaluation of the Hessian... Log( Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Evaluate the Hessian..." ); @@ -337,7 +337,7 @@ void Calculate( Data::HTST_Info & htst_info ) // Transform into geodesic Hessian Log( Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Transform Hessian into geodesic Hessian..." ); SpMatrixX sparse_hessian_sp_geodesic_3N( 3 * nos, 3 * nos ); - sparse_hessian_bordered_3N( image_sp, gradient_sp, sparse_hessian_sp, sparse_hessian_sp_geodesic_3N ); + sparse_hessian_bordered_3N( image_sp.spin, gradient_sp, sparse_hessian_sp, sparse_hessian_sp_geodesic_3N ); SpMatrixX sparse_hessian_sp_geodesic_2N = tangent_basis.transpose() * sparse_hessian_sp_geodesic_3N * tangent_basis; @@ -394,7 +394,7 @@ void Calculate( Data::HTST_Info & htst_info ) Log( Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Evaluate the dynamical matrix" ); SpMatrixX velocity( 3 * nos, 3 * nos ); Sparse_Calculate_Dynamical_Matrix( - image_sp, htst_info.saddle_point->hamiltonian->get_geometry().mu_s, sparse_hessian_sp_geodesic_3N, + image_sp.spin, htst_info.saddle_point->hamiltonian->get_geometry().mu_s, sparse_hessian_sp_geodesic_3N, velocity ); SpMatrixX projected_velocity = tangent_basis.transpose() * velocity * tangent_basis; @@ -432,7 +432,7 @@ void Calculate( Data::HTST_Info & htst_info ) Log( Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Evaluate tangent basis ..." ); SpMatrixX tangent_basis = SpMatrixX( 3 * nos, 2 * nos ); - Manifoldmath::sparse_tangent_basis_spherical( image_minimum, tangent_basis ); + Manifoldmath::sparse_tangent_basis_spherical( image_minimum.spin, tangent_basis ); // Evaluation of the Hessian... Log( Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Evaluate the Hessian..." ); @@ -443,7 +443,7 @@ void Calculate( Data::HTST_Info & htst_info ) Log( Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Transforming Hessian into geodesic Hessian..." ); SpMatrixX sparse_hessian_geodesic_min_3N = SpMatrixX( 3 * nos, 3 * nos ); sparse_hessian_bordered_3N( - image_minimum, gradient_minimum, sparse_hessian_minimum, sparse_hessian_geodesic_min_3N ); + image_minimum.spin, gradient_minimum, sparse_hessian_minimum, sparse_hessian_geodesic_min_3N ); SpMatrixX sparse_hessian_geodesic_min_2N = tangent_basis.transpose() * sparse_hessian_geodesic_min_3N * tangent_basis; diff --git a/core/src/engine/spin/interaction/DDI.cpp b/core/src/engine/spin/interaction/DDI.cpp index 9b1372c22..12ecbac41 100644 --- a/core/src/engine/spin/interaction/DDI.cpp +++ b/core/src/engine/spin/interaction/DDI.cpp @@ -108,7 +108,7 @@ void DDI::applyGeometry( } template<> -void DDI::Energy::operator()( const vectorfield & spins, scalarfield & energy ) const +void DDI::Energy::operator()( const StateType & state, scalarfield & energy ) const { if( !is_contributing ) return; @@ -118,19 +118,19 @@ void DDI::Energy::operator()( const vectorfield & spins, scalarfield & energy ) return; if( data.method == DDI_Method::FFT ) - Energy_per_Spin_FFT( *cache.geometry, *cache.boundary_conditions, cache, spins, energy ); + Energy_per_Spin_FFT( *cache.geometry, *cache.boundary_conditions, cache, state.spin, energy ); else if( data.method == DDI_Method::Cutoff ) { // TODO: Merge these implementations in the future if( data.cutoff_radius >= 0 ) - Energy_per_Spin_Cutoff( *cache.geometry, *cache.boundary_conditions, cache, spins, energy ); + Energy_per_Spin_Cutoff( *cache.geometry, *cache.boundary_conditions, cache, state.spin, energy ); else - Energy_per_Spin_Direct( *cache.geometry, *cache.boundary_conditions, data, spins, energy ); + Energy_per_Spin_Direct( *cache.geometry, *cache.boundary_conditions, data, state.spin, energy ); } }; template<> -void DDI::Gradient::operator()( const vectorfield & spins, vectorfield & gradient ) const +void DDI::Gradient::operator()( const StateType & state, vectorfield & gradient ) const { if( !is_contributing ) return; @@ -140,21 +140,21 @@ void DDI::Gradient::operator()( const vectorfield & spins, vectorfield & gradien return; if( data.method == DDI_Method::FFT ) - Gradient_FFT( *cache.geometry, *cache.boundary_conditions, cache, spins, gradient ); + Gradient_FFT( *cache.geometry, *cache.boundary_conditions, cache, state.spin, gradient ); else if( data.method == DDI_Method::Cutoff ) { // TODO: Merge these implementations in the future if( data.cutoff_radius >= 0 ) - Gradient_Cutoff( *cache.geometry, *cache.boundary_conditions, cache, spins, gradient ); + Gradient_Cutoff( *cache.geometry, *cache.boundary_conditions, cache, state.spin, gradient ); else - Gradient_Direct( *cache.geometry, *cache.boundary_conditions, data, spins, gradient ); + Gradient_Direct( *cache.geometry, *cache.boundary_conditions, data, state.spin, 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. template<> -scalar DDI::Energy_Single_Spin::operator()( const int ispin, const vectorfield & spins ) const +scalar DDI::Energy_Single_Spin::operator()( const int ispin, const StateType & state ) const { if( !is_contributing ) return 0; diff --git a/core/src/io/Dataparser.cpp b/core/src/io/Dataparser.cpp index ba8bcd86c..fb9661b1c 100644 --- a/core/src/io/Dataparser.cpp +++ b/core/src/io/Dataparser.cpp @@ -29,8 +29,7 @@ namespace Spin // Reads a non-OVF spins file with plain text and discarding any headers starting with '#' void Read_NonOVF_System_Configuration( - vectorfield & spins, Data::Geometry & geometry, const int nos, const int idx_image_infile, - const std::string & file ) + StateType & state, Data::Geometry & geometry, const int nos, const int idx_image_infile, const std::string & file ) { IO::Filter_File_Handle file_handle( file, "#" ); @@ -40,13 +39,13 @@ void Read_NonOVF_System_Configuration( for( int i = 0; i < nos && file_handle.GetLine( "," ); i++ ) { - file_handle >> spins[i][0]; - file_handle >> spins[i][1]; - file_handle >> spins[i][2]; + file_handle >> state.spin[i][0]; + file_handle >> state.spin[i][1]; + file_handle >> state.spin[i][2]; - if( spins[i].norm() < 1e-5 ) + if( state.spin[i].norm() < 1e-5 ) { - spins[i] = { 0, 0, 1 }; + state.spin[i] = { 0, 0, 1 }; // In case of spin vector close to zero we have a vacancy #ifdef SPIRIT_ENABLE_DEFECTS geometry.atom_types[i] = -1; @@ -55,7 +54,7 @@ void Read_NonOVF_System_Configuration( } // normalize read in spins - Engine::Vectormath::normalize_vectors( spins ); + Engine::Vectormath::normalize_vectors( state.spin ); } } // namespace Spin diff --git a/core/src/utility/Configuration_Chain.cpp b/core/src/utility/Configuration_Chain.cpp index a51d698e5..a69cc10e7 100644 --- a/core/src/utility/Configuration_Chain.cpp +++ b/core/src/utility/Configuration_Chain.cpp @@ -27,14 +27,14 @@ void Add_Noise_Temperature( State::chain_t & c, int idx_1, int idx_2, scalar tem for( int img = idx_1 + 1; img <= idx_2 - 1; ++img ) { Configurations::Add_Noise_Temperature_Sphere( - get( *c.images[img]->state ), c.images[img]->hamiltonian->get_geometry(), temperature, prng ); + c.images[img]->state->spin, c.images[img]->hamiltonian->get_geometry(), temperature, prng ); } } void Homogeneous_Rotation( State::chain_t & c, int idx_1, int idx_2 ) { - auto & spins_1 = get( *c.images[idx_1]->state ); - auto & spins_2 = get( *c.images[idx_2]->state ); + auto & spins_1 = c.images[idx_1]->state->spin; + auto & spins_2 = c.images[idx_2]->state->spin; scalar angle, rot_angle; Vector3 rot_axis; @@ -63,7 +63,7 @@ void Homogeneous_Rotation( State::chain_t & c, int idx_1, int idx_2 ) for( int img = idx_1 + 1; img < idx_2; ++img ) { angle = rot_angle * scalar( img - idx_1 ) / scalar( idx_2 - idx_1 ); - Engine::Vectormath::rotate( spins_1[i], rot_axis, angle, get( *c.images[img]->state )[i] ); + Engine::Vectormath::rotate( spins_1[i], rot_axis, angle, c.images[img]->state->spin[i] ); } } // Otherwise we simply leave the spin untouched @@ -71,7 +71,7 @@ void Homogeneous_Rotation( State::chain_t & c, int idx_1, int idx_2 ) { for( int img = idx_1 + 1; img < idx_2; ++img ) { - get( *c.images[img]->state )[i] = spins_1[i]; + c.images[img]->state->spin[i] = spins_1[i]; } } } diff --git a/core/test/test_aggregation.cpp b/core/test/test_aggregation.cpp index 943b292ea..5d782cd3e 100644 --- a/core/test/test_aggregation.cpp +++ b/core/test/test_aggregation.cpp @@ -58,7 +58,7 @@ TEST_CASE( "Ensure that Hamiltonian is really just an aggregator", "[aggregation Configuration_Random( state.get() ); const auto & spins = *state->active_image->state; auto & hamiltonian = state->active_image->hamiltonian; - auto nos = get( spins ).size(); + auto nos = spins.spin.size(); if( hamiltonian->active_count() == 0 ) { @@ -86,7 +86,7 @@ TEST_CASE( "Ensure that Hamiltonian is really just an aggregator", "[aggregation scalarfield( nos, 0 ), [&spins]( const scalarfield & v, const auto & interaction ) -> scalarfield { - const auto nos = get( spins ).size(); + const auto nos = spins.spin.size(); auto energy_per_spin = scalarfield( nos, 0 ); interaction->Energy_per_Spin( spins, energy_per_spin ); #pragma omp parallel for @@ -111,7 +111,7 @@ TEST_CASE( "Ensure that Hamiltonian is really just an aggregator", "[aggregation vectorfield( nos, Vector3::Zero() ), [&spins]( const vectorfield & v, const auto & interaction ) -> vectorfield { - auto gradient = vectorfield( get( spins ).size(), Vector3::Zero() ); + auto gradient = vectorfield( spins.spin.size(), Vector3::Zero() ); interaction->Gradient( spins, gradient ); Engine::Vectormath::add_c_a( 1.0, v, gradient ); return gradient; diff --git a/core/test/test_anisotropy.cpp b/core/test/test_anisotropy.cpp index e8d07587e..d33619604 100644 --- a/core/test/test_anisotropy.cpp +++ b/core/test/test_anisotropy.cpp @@ -76,31 +76,32 @@ TEST_CASE( "Uniaxial nisotropy", "[anisotropy]" ) REQUIRE_THAT( normal[2], within_digits( init_normal[2], 12 ) ); } + using Engine::Spin::StateType; SECTION( "Total energies for different orientations should match expected values" ) { - vectorfield spins( state->nos ); + StateType system_state{ vectorfield( state->nos ) }; - for( auto & spin : spins ) + for( auto & spin : system_state.spin ) spin = { 1.0, 0.0, 0.0 }; - scalar energy_x = state->active_image->hamiltonian->Energy( spins ); + scalar energy_x = state->active_image->hamiltonian->Energy( system_state ); // X and Z orientations energies should differ by NOS*init_magnitude - for( auto & spin : spins ) + for( auto & spin : system_state.spin ) spin = { 0.0, 0.0, 1.0 }; - scalar energy_z = state->active_image->hamiltonian->Energy( spins ); + scalar energy_z = state->active_image->hamiltonian->Energy( system_state ); REQUIRE_THAT( energy_x - energy_z, within_digits( init_magnitude * state->nos, digits_a ) ); // X and XY orientations energies should have equal energies scalar sqrt2_2 = std::sqrt( 2 ) / 2; - for( auto & spin : spins ) + for( auto & spin : system_state.spin ) spin = { sqrt2_2, sqrt2_2, 0.0 }; - scalar energy_xy = state->active_image->hamiltonian->Energy( spins ); + scalar energy_xy = state->active_image->hamiltonian->Energy( system_state ); REQUIRE_THAT( energy_x - energy_xy, within_digits( 0, 12 ) ); } SECTION( "All individual energy gradients should match the expected value" ) { - vectorfield spins( state->nos, { 0.0, 0.0, 1.0 } ); + StateType spins{ vectorfield( state->nos, { 0.0, 0.0, 1.0 } ) }; auto gradients = vectorfield( state->nos ); state->active_image->hamiltonian->Gradient( spins, gradients ); @@ -118,8 +119,10 @@ TEST_CASE( "Uniaxial nisotropy", "[anisotropy]" ) TEST_CASE( "Cubic anisotropy", "[anisotropy]" ) { + using Engine::Spin::StateType; + auto state = std::shared_ptr( State_Setup(), State_Delete ); - vectorfield spins( state->nos ); + StateType system_state{ vectorfield( state->nos ) }; // Set uniaxial anisotropy to zero scalar init_normal_uniaxial[3] = { 0.0, 0.0, 1.0 }; @@ -138,35 +141,35 @@ TEST_CASE( "Cubic anisotropy", "[anisotropy]" ) SECTION( "Total energies for different orientations should match expected values" ) { scalar sqrt2_2 = std::sqrt( 2 ) / 2; - for( auto & spin : spins ) + for( auto & spin : system_state.spin ) spin = { sqrt2_2, sqrt2_2, 0.0 }; - scalar energy_xy = state->active_image->hamiltonian->Energy( spins ); + scalar energy_xy = state->active_image->hamiltonian->Energy( system_state ); // X and XY orientations energies should differ by NOS*init_magnitude/4 - for( auto & spin : spins ) + for( auto & spin : system_state.spin ) spin = { 1.0, 0.0, 0.0 }; - scalar energy_x = state->active_image->hamiltonian->Energy( spins ); + scalar energy_x = state->active_image->hamiltonian->Energy( system_state ); REQUIRE_THAT( energy_x - energy_xy, within_digits( -init_magnitude / 4 * state->nos, digits_a ) ); // Y and XY orientations energies should differ by NOS*init_magnitude/4 - for( auto & spin : spins ) + for( auto & spin : system_state.spin ) spin = { 0.0, 1.0, 0.0 }; - scalar energy_y = state->active_image->hamiltonian->Energy( spins ); + scalar energy_y = state->active_image->hamiltonian->Energy( system_state ); REQUIRE_THAT( energy_y - energy_xy, within_digits( -init_magnitude / 4 * state->nos, digits_a ) ); // Y and Z orientations should have equal energies - for( auto & spin : spins ) + for( auto & spin : system_state.spin ) spin = { 0.0, 0.0, 1.0 }; - scalar energy_z = state->active_image->hamiltonian->Energy( spins ); + scalar energy_z = state->active_image->hamiltonian->Energy( system_state ); REQUIRE_THAT( energy_y - energy_z, within_digits( 0, 12 ) ); } SECTION( "All individual energy gradients should match the expected value" ) { - for( auto & spin : spins ) + for( auto & spin : system_state.spin ) spin = { 0.0, 0.0, 1.0 }; auto gradients = vectorfield( state->nos ); - state->active_image->hamiltonian->Gradient( spins, gradients ); + state->active_image->hamiltonian->Gradient( system_state, gradients ); Vector3 gradient_expected{ 0.0, 0.0, scalar( -2.0 * init_magnitude ) }; @@ -183,19 +186,19 @@ TEST_CASE( "Cubic anisotropy", "[anisotropy]" ) { Hamiltonian_Set_Anisotropy( state.get(), 0.1, init_normal_uniaxial ); - for( auto & spin : spins ) + for( auto & spin : system_state.spin ) spin = { 0.0, 0.0, 1.0 }; // Direct energy calculation and energy calculated from gradient should be equal auto gradients_a = vectorfield( state->nos ); scalar energy_from_gradient{}; - state->active_image->hamiltonian->Gradient_and_Energy( spins, gradients_a, energy_from_gradient ); - scalar energy_direct = state->active_image->hamiltonian->Energy( spins ); + state->active_image->hamiltonian->Gradient_and_Energy( system_state, gradients_a, energy_from_gradient ); + scalar energy_direct = state->active_image->hamiltonian->Energy( system_state ); REQUIRE_THAT( energy_from_gradient, within_digits( energy_direct, 12 ) ); // Direct gradient calculation and gradient out of gradient-and-energy calculation should be equal auto gradients_b = vectorfield( state->nos ); - state->active_image->hamiltonian->Gradient( spins, gradients_b ); + state->active_image->hamiltonian->Gradient( system_state, gradients_b ); for( int idx = 0; idx < state->nos; ++idx ) { INFO( @@ -209,6 +212,8 @@ TEST_CASE( "Cubic anisotropy", "[anisotropy]" ) TEST_CASE( "Biaxial anisotropy", "[anisotropy]" ) { + using Engine::Spin::StateType; + auto state = std::shared_ptr( State_Setup(), State_Delete ); int idx_image = -1, idx_chain = -1; @@ -329,9 +334,9 @@ TEST_CASE( "Biaxial anisotropy", "[anisotropy]" ) Hamiltonian_Set_Biaxial_Anisotropy( state.get(), magnitude.data(), array_cast( exponents ), init_primary.data(), init_secondary.data(), N ); - vectorfield spins( state->nos, make_spherical( theta, phi ) ); + StateType system_state{ vectorfield( state->nos, make_spherical( theta, phi ) ) }; scalarfield energy( state->nos, 0.0 ); - interaction->Energy_per_Spin( spins, energy ); + interaction->Energy_per_Spin( system_state, energy ); // reference energy const scalar ref_energy @@ -351,10 +356,11 @@ TEST_CASE( "Biaxial anisotropy", "[anisotropy]" ) INFO( "trial: " << idx << ", theta=" << theta << ", phi=" << phi ); INFO( term_info( terms... ) ); REQUIRE_THAT( - interaction->Energy( spins ) / static_cast( state->nos ), WithinAbs( ref_energy, epsilon_3 ) ); + interaction->Energy( system_state ) / static_cast( state->nos ), + WithinAbs( ref_energy, epsilon_3 ) ); vectorfield gradient( state->nos, Vector3::Zero() ); - interaction->Gradient( spins, gradient ); + interaction->Gradient( system_state, gradient ); const auto k1 = Vector3{ init_primary[0], init_primary[1], init_primary[2] }; const auto k2 = Vector3{ init_secondary[0], init_secondary[1], init_secondary[2] }; diff --git a/core/test/test_physics.cpp b/core/test/test_physics.cpp index c23995669..e491a7ffa 100644 --- a/core/test/test_physics.cpp +++ b/core/test/test_physics.cpp @@ -219,8 +219,8 @@ TEST_CASE( "Finite difference and regular Hamiltonian should match", "[physics]" auto state = std::shared_ptr( State_Setup( input_file ), State_Delete ); Configuration_Random( state.get() ); - const auto & spins = *state->active_image->state; - auto & hamiltonian = state->active_image->hamiltonian; + const auto & system_state = *state->active_image->state; + auto & hamiltonian = state->active_image->hamiltonian; // Compare gradients auto grad = vectorfield( state->nos, Vector3::Zero() ); @@ -229,10 +229,10 @@ TEST_CASE( "Finite difference and regular Hamiltonian should match", "[physics]" { Engine::Vectormath::fill( grad, Vector3::Zero() ); Engine::Vectormath::fill( grad_fd, Vector3::Zero() ); - interaction->Gradient( spins, grad ); + interaction->Gradient( system_state, grad ); Engine::Vectormath::Gradient( - spins, grad_fd, - [&interaction]( const auto & spins ) -> scalar { return interaction->Energy( spins ); } ); + system_state, grad_fd, + [&interaction]( const auto & state ) -> scalar { return interaction->Energy( state ); } ); INFO( "Interaction: " << interaction->Name() << "\n" ); for( int i = 0; i < state->nos; i++ ) { @@ -252,9 +252,9 @@ TEST_CASE( "Finite difference and regular Hamiltonian should match", "[physics]" hessian_fd.setZero(); Engine::Vectormath::Hessian( - spins, hessian_fd, - [&interaction]( const auto & spins, auto & gradient ) { interaction->Gradient( spins, gradient ); } ); - interaction->Hessian( spins, hessian ); + system_state, hessian_fd, + [&interaction]( const auto & state, auto & gradient ) { interaction->Gradient( state, gradient ); } ); + interaction->Hessian( system_state, hessian ); INFO( "Interaction: " << interaction->Name() << "\n" ); INFO( "epsilon = " << epsilon_3 << "\n" ); INFO( "Hessian (FD) =\n" << hessian_fd << "\n" ); @@ -271,19 +271,19 @@ TEST_CASE( "Dipole-Dipole Interaction", "[physics]" ) auto state = std::shared_ptr( State_Setup( input_file ), State_Delete ); Configuration_Random( state.get() ); - auto & spins = *state->active_image->state; + auto & system_state = *state->active_image->state; auto ddi_interaction = state->active_image->hamiltonian->getInteraction(); // FFT gradient and energy auto grad_fft = vectorfield( state->nos, Vector3::Zero() ); - ddi_interaction->Gradient( spins, grad_fft ); - auto energy_fft = ddi_interaction->Energy( spins ); + ddi_interaction->Gradient( system_state, grad_fft ); + auto energy_fft = ddi_interaction->Energy( system_state ); { auto grad_fft_fd = vectorfield( state->nos, Vector3::Zero() ); Engine::Vectormath::Gradient( - spins, grad_fft_fd, - [&ddi_interaction]( const auto & spins ) -> scalar { return ddi_interaction->Energy( spins ); } ); + system_state, grad_fft_fd, + [&ddi_interaction]( const auto & state ) -> scalar { return ddi_interaction->Energy( state ); } ); INFO( "Interaction: " << ddi_interaction->Name() << "\n" ); for( int i = 0; i < state->nos; i++ ) @@ -298,8 +298,8 @@ TEST_CASE( "Dipole-Dipole Interaction", "[physics]" ) auto n_periodic_images = std::vector{ 4, 4, 4 }; Hamiltonian_Set_DDI( state.get(), SPIRIT_DDI_METHOD_CUTOFF, n_periodic_images.data(), -1 ); auto grad_direct = vectorfield( state->nos, Vector3::Zero() ); - ddi_interaction->Gradient( spins, grad_direct ); - auto energy_direct = ddi_interaction->Energy( spins ); + ddi_interaction->Gradient( system_state, grad_direct ); + auto energy_direct = ddi_interaction->Energy( system_state ); // Compare gradients for( int i = 0; i < state->nos; i++ )