diff --git a/core/include/Spirit/Parameters_EMA.h b/core/include/Spirit/Parameters_EMA.h index bfa8195b4..e0cd6eaf0 100644 --- a/core/include/Spirit/Parameters_EMA.h +++ b/core/include/Spirit/Parameters_EMA.h @@ -19,6 +19,9 @@ This method, if needed, calculates modes (they can also be read in from a file) and perturbs the spin system periodically in the direction of the eigenmode. */ +// Clears all the previously calculated modes from memory +PREFIX void Parameters_EMA_Clear_Modes( State * state, int idx_image = -1, int idx_chain = -1 ) SUFFIX; + /* Set -------------------------------------------------------------------- @@ -42,6 +45,9 @@ Parameters_EMA_Set_Amplitude( State * state, float amplitude, int idx_image = -1 // Set whether to displace the system statically instead of periodically. PREFIX void Parameters_EMA_Set_Snapshot( State * state, bool snapshot, int idx_image = -1, int idx_chain = -1 ) SUFFIX; +// Set whether to use sparse matrices. +PREFIX void Parameters_EMA_Set_Sparse( State * state, bool sparse, int idx_image = -1, int idx_chain = -1 ) SUFFIX; + /* Get -------------------------------------------------------------------- @@ -62,5 +68,8 @@ PREFIX float Parameters_EMA_Get_Amplitude( State * state, int idx_image = -1, in // Returns whether to displace the system statically instead of periodically. PREFIX bool Parameters_EMA_Get_Snapshot( State * state, int idx_image = -1, int idx_chain = -1 ) SUFFIX; +// Returns whether to use sparse matrices. +PREFIX bool Parameters_EMA_Get_Sparse( State * state, int idx_image = -1, int idx_chain = -1 ) SUFFIX; + #include "DLL_Undefine_Export.h" #endif \ No newline at end of file diff --git a/core/include/Spirit/Parameters_GNEB.h b/core/include/Spirit/Parameters_GNEB.h index 9b930e7cc..ff25fd0a2 100644 --- a/core/include/Spirit/Parameters_GNEB.h +++ b/core/include/Spirit/Parameters_GNEB.h @@ -111,6 +111,15 @@ PREFIX void Parameters_GNEB_Set_Spring_Force_Ratio( State * state, float ratio, PREFIX void Parameters_GNEB_Set_Path_Shortening_Constant( State * state, float path_shortening_constant, int idx_chain = -1 ) SUFFIX; +// Set if moving endpoints should be used +PREFIX void Parameters_GNEB_Set_Moving_Endpoints( State * state, bool moving_endpoints, int idx_chain = -1 ) SUFFIX; + +// Set if attracting endpoints should be used +PREFIX void Parameters_GNEB_Set_Translating_Endpoints( State * state, bool translating_endpoints, int idx_chain = -1 ) SUFFIX; + +// Set equilibrium Rx, used for the moving endpoints method +PREFIX void Parameters_GNEB_Set_Equilibrium_Delta_Rx( State * state, float delta_Rx_left, float delta_Rx_right, int idx_chain = -1 ) SUFFIX; + // Set the GNEB image type (see the integers defined above). PREFIX void Parameters_GNEB_Set_Climbing_Falling( State * state, int image_type, int idx_image = -1, int idx_chain = -1 ) SUFFIX; @@ -169,6 +178,15 @@ PREFIX float Parameters_GNEB_Get_Spring_Force_Ratio( State * state, int idx_chai // Return the path shortening constant. PREFIX float Parameters_GNEB_Get_Path_Shortening_Constant( State * state, int idx_chain = -1 ) SUFFIX; +// Return if moving endpoints are used +PREFIX bool Parameters_GNEB_Get_Moving_Endpoints( State * state, int idx_chain = -1 ) SUFFIX; + +// Set if translating endpoints are used +PREFIX bool Parameters_GNEB_Get_Translating_Endpoints( State * state, int idx_chain = -1 ) SUFFIX; + +// Get the equilibrium Rx, used for the moving endpoints method +PREFIX void Parameters_GNEB_Get_Equilibrium_Delta_Rx( State * state, float * delta_Rx_left, float * delta_Rx_right, int idx_chain = -1 ) SUFFIX; + /* Returns the integer of whether an image is regular, climbing, falling, or stationary. diff --git a/core/include/Spirit/Quantities.h b/core/include/Spirit/Quantities.h index 6826d49f8..21e14688c 100644 --- a/core/include/Spirit/Quantities.h +++ b/core/include/Spirit/Quantities.h @@ -14,12 +14,19 @@ Quantities ``` */ +// Average spin direction +PREFIX void Quantity_Get_Average_Spin( State * state, float s[3], int idx_image = -1, int idx_chain = -1 ); + // Total Magnetization PREFIX void Quantity_Get_Magnetization( State * state, float m[3], int idx_image = -1, int idx_chain = -1 ); // Topological Charge PREFIX float Quantity_Get_Topological_Charge( State * state, int idx_image = -1, int idx_chain = -1 ); +// Topological Charge Density +PREFIX int Quantity_Get_Topological_Charge_Density( + State * state, float * charge_density, int * triangle_indices, int idx_image = -1, int idx_chain = -1 ); + // Minimum mode following information PREFIX void Quantity_Get_Grad_Force_MinimumMode( State * state, float * gradient, float * eval, float * mode, float * forces, int idx_image = -1, diff --git a/core/include/Spirit/Simulation.h b/core/include/Spirit/Simulation.h index a62e6d7de..bcfd84bb8 100644 --- a/core/include/Spirit/Simulation.h +++ b/core/include/Spirit/Simulation.h @@ -61,8 +61,17 @@ struct Simulation_Run_Info int total_walltime = 0; float total_ips = 0; float max_torque = 0; + + int n_history_iteration = 0; + int * history_iteration = nullptr; + int n_history_max_torque = 0; + float * history_max_torque = nullptr; + int n_history_energy = 0; + float * history_energy = nullptr; }; +PREFIX void free_run_info(Simulation_Run_Info info) SUFFIX; + /* Start or stop a simulation -------------------------------------------------------------------- @@ -71,27 +80,27 @@ Start or stop a simulation // Monte Carlo PREFIX void Simulation_MC_Start( State * state, int n_iterations = -1, int n_iterations_log = -1, bool singleshot = false, - Simulation_Run_Info * info = NULL, int idx_image = -1, int idx_chain = -1 ) SUFFIX; + Simulation_Run_Info * info = nullptr, int idx_image = -1, int idx_chain = -1 ) SUFFIX; // Landau-Lifshitz-Gilbert dynamics and energy minimisation PREFIX void Simulation_LLG_Start( State * state, int solver_type, int n_iterations = -1, int n_iterations_log = -1, bool singleshot = false, - Simulation_Run_Info * info = NULL, int idx_image = -1, int idx_chain = -1 ) SUFFIX; + Simulation_Run_Info * info = nullptr, int idx_image = -1, int idx_chain = -1 ) SUFFIX; // Geodesic nudged elastic band method PREFIX void Simulation_GNEB_Start( State * state, int solver_type, int n_iterations = -1, int n_iterations_log = -1, bool singleshot = false, - Simulation_Run_Info * info = NULL, int idx_chain = -1 ) SUFFIX; + Simulation_Run_Info * info = nullptr, int idx_chain = -1 ) SUFFIX; // Minimum mode following method PREFIX void Simulation_MMF_Start( State * state, int solver_type, int n_iterations = -1, int n_iterations_log = -1, bool singleshot = false, - Simulation_Run_Info * info = NULL, int idx_image = -1, int idx_chain = -1 ) SUFFIX; + Simulation_Run_Info * info = nullptr, int idx_image = -1, int idx_chain = -1 ) SUFFIX; // Eigenmode analysis PREFIX void Simulation_EMA_Start( State * state, int n_iterations = -1, int n_iterations_log = -1, bool singleshot = false, - Simulation_Run_Info * info = NULL, int idx_image = -1, int idx_chain = -1 ) SUFFIX; + Simulation_Run_Info * info = nullptr, int idx_image = -1, int idx_chain = -1 ) SUFFIX; /* Single iteration of a Method diff --git a/core/include/data/Geometry.hpp b/core/include/data/Geometry.hpp index 9a748950c..0d1628da4 100644 --- a/core/include/data/Geometry.hpp +++ b/core/include/data/Geometry.hpp @@ -81,8 +81,9 @@ class Geometry // ---------- Constructor // Build a regular lattice from a defined basis cell and translations Geometry( - std::vector bravais_vectors, intfield n_cells, std::vector cell_atoms, - Basis_Cell_Composition cell_composition, scalar lattice_constant, Pinning pinning, Defects defects ); + const std::vector & bravais_vectors, intfield n_cells, const std::vector & cell_atoms, + const Basis_Cell_Composition & cell_composition, scalar lattice_constant, const Pinning & pinning, + const Defects & defects ); // ---------- Convenience functions // Retrieve triangulation, if 2D diff --git a/core/include/data/Parameters_Method_EMA.hpp b/core/include/data/Parameters_Method_EMA.hpp index 77612a914..0612c1589 100644 --- a/core/include/data/Parameters_Method_EMA.hpp +++ b/core/include/data/Parameters_Method_EMA.hpp @@ -15,6 +15,7 @@ struct Parameters_Method_EMA : Parameters_Method scalar frequency = 0.02; scalar amplitude = 1; bool snapshot = false; + bool sparse = false; // ----------------- Output -------------- // Energy output settings diff --git a/core/include/data/Parameters_Method_GNEB.hpp b/core/include/data/Parameters_Method_GNEB.hpp index b80a06883..6b4b0165b 100644 --- a/core/include/data/Parameters_Method_GNEB.hpp +++ b/core/include/data/Parameters_Method_GNEB.hpp @@ -33,6 +33,14 @@ struct Parameters_Method_GNEB : Parameters_Method_Solver // Mersenne twister PRNG std::mt19937 prng = std::mt19937( rng_seed ); + bool moving_endpoints = false; + bool translating_endpoints = false; + + scalar equilibrium_delta_Rx_left = 1.0; + scalar equilibrium_delta_Rx_right = 1.0; + + bool escape_first = false; + // ----------------- Output -------------- bool output_energies_step = false; bool output_energies_divide_by_nspins = true; diff --git a/core/include/engine/Eigenmodes.hpp b/core/include/engine/Eigenmodes.hpp index 9b1d0cc32..3e3a4d878 100644 --- a/core/include/engine/Eigenmodes.hpp +++ b/core/include/engine/Eigenmodes.hpp @@ -33,8 +33,16 @@ bool Hessian_Full_Spectrum( // gradient and hessian should be the 3N-dimensional representations without constraints bool Hessian_Partial_Spectrum( const std::shared_ptr parameters, const vectorfield & spins, const vectorfield & gradient, - const MatrixX & hessian, std::size_t n_modes, MatrixX & tangent_basis, MatrixX & hessian_constrained, - VectorX & eigenvalues, MatrixX & eigenvectors ); + const MatrixX & hessian, std::size_t n_modes, MatrixX & tangent_basis, MatrixX & hessian_constrained, VectorX & eigenvalues, + MatrixX & eigenvcetors ); + +// Calculate a partial eigenspectrum of a sparse Hessian +// gradient and hessian should be the 3N-dimensional representations without constraints +bool Sparse_Hessian_Partial_Spectrum( + const std::shared_ptr parameters, const vectorfield & spins, const vectorfield & gradient, + const SpMatrixX & hessian, int n_modes, SpMatrixX & tangent_basis, SpMatrixX & hessian_constrained, VectorX & eigenvalues, + MatrixX & eigenvectors +); } // end namespace Eigenmodes } // end namespace Engine diff --git a/core/include/engine/Hamiltonian.hpp b/core/include/engine/Hamiltonian.hpp index b2bc5e7e6..f4826e53c 100644 --- a/core/include/engine/Hamiltonian.hpp +++ b/core/include/engine/Hamiltonian.hpp @@ -96,7 +96,7 @@ class Hamiltonian virtual std::size_t Number_of_Interactions(); // Hamiltonian name as string - virtual const std::string & Name(); + virtual const std::string & Name() const = 0; // Boundary conditions intfield boundary_conditions; // [3] (a, b, c) diff --git a/core/include/engine/Hamiltonian_Gaussian.hpp b/core/include/engine/Hamiltonian_Gaussian.hpp index 4ed4cb96a..e07269b65 100644 --- a/core/include/engine/Hamiltonian_Gaussian.hpp +++ b/core/include/engine/Hamiltonian_Gaussian.hpp @@ -36,7 +36,7 @@ class Hamiltonian_Gaussian : public Hamiltonian scalar Energy_Single_Spin( int ispin, const vectorfield & spins ) override; // Hamiltonian name as string - const std::string & Name() override; + const std::string & Name() const override; // Parameters of the energy landscape int n_gaussians; diff --git a/core/include/engine/Hamiltonian_Heisenberg.hpp b/core/include/engine/Hamiltonian_Heisenberg.hpp index 3e185aa3c..5ff650bfe 100644 --- a/core/include/engine/Hamiltonian_Heisenberg.hpp +++ b/core/include/engine/Hamiltonian_Heisenberg.hpp @@ -56,7 +56,7 @@ class Hamiltonian_Heisenberg : public Hamiltonian scalar Energy_Single_Spin( int ispin, const vectorfield & spins ) override; // Hamiltonian name as string - const std::string & Name() override; + const std::string & Name() const override; // ------------ Single Spin Interactions ------------ // External magnetic field across the sample diff --git a/core/include/engine/Hamiltonian_Micromagnetic.hpp b/core/include/engine/Hamiltonian_Micromagnetic.hpp index 392ac042c..eee550e22 100644 --- a/core/include/engine/Hamiltonian_Micromagnetic.hpp +++ b/core/include/engine/Hamiltonian_Micromagnetic.hpp @@ -42,7 +42,7 @@ class Hamiltonian_Micromagnetic : public Hamiltonian scalar Energy_Single_Spin( int ispin, const vectorfield & spins ) override; // Hamiltonian name as string - const std::string & Name() override; + const std::string & Name() const override; std::shared_ptr geometry; diff --git a/core/include/engine/Manifoldmath.hpp b/core/include/engine/Manifoldmath.hpp index e7c84bca2..c7f8d0030 100644 --- a/core/include/engine/Manifoldmath.hpp +++ b/core/include/engine/Manifoldmath.hpp @@ -87,6 +87,9 @@ void spherical_to_cartesian_coordinate_basis( const vectorfield & vf, MatrixX & // void spherical_coordinate_christoffel_symbols( const vectorfield & vf, MatrixX & gamma_theta, MatrixX & gamma_phi ); +void sparse_hessian_bordered_3N( + const vectorfield & image, const vectorfield & gradient, const SpMatrixX & hessian, SpMatrixX & hessian_out ); + // Calculate Hessian for a vectorfield constrained to unit length, at any extremum (i.e. where vectors || gradient) void hessian_bordered( const vectorfield & image, const vectorfield & gradient, const MatrixX & hessian, MatrixX & tangent_basis, diff --git a/core/include/engine/Method.hpp b/core/include/engine/Method.hpp index 5bfe31586..0daff1414 100644 --- a/core/include/engine/Method.hpp +++ b/core/include/engine/Method.hpp @@ -158,6 +158,11 @@ class Method // Number of times to save long n_log; + // History of relevant quantities + std::vector history_iteration; + std::vector history_max_torque; + std::vector history_energy; + protected: //////////// Information for Logging and Save_Current //////////////////////// int idx_image; @@ -182,9 +187,6 @@ class Method // Maximum force component per image std::vector force_max_abs_component_all; - // History of relevant quantities - std::map> history; - //////////// General ///////////////////////////////////////////////////////// // Systems the Solver will access diff --git a/core/include/engine/Method_GNEB.hpp b/core/include/engine/Method_GNEB.hpp index 5509b2c7d..572fe383b 100644 --- a/core/include/engine/Method_GNEB.hpp +++ b/core/include/engine/Method_GNEB.hpp @@ -75,8 +75,15 @@ class Method_GNEB : public Method_Solver std::vector F_gradient; std::vector F_spring; vectorfield f_shrink; + + vectorfield F_translation_left; + vectorfield F_translation_right; + // Last calculated tangents std::vector tangents; + + vectorfield tangent_endpoints_left; + vectorfield tangent_endpoints_right; }; } // namespace Engine diff --git a/core/include/engine/Solver_Depondt.hpp b/core/include/engine/Solver_Depondt.hpp index 4abdcd3f0..25ed98c3e 100644 --- a/core/include/engine/Solver_Depondt.hpp +++ b/core/include/engine/Solver_Depondt.hpp @@ -16,7 +16,7 @@ inline void Method_Solver::Initialize() configurations_predictor[i] = std::shared_ptr( new vectorfield( this->nos, { 0, 0, 0 } ) ); this->temp1 = vectorfield( this->nos, { 0, 0, 0 } ); -}; +} /* Template instantiation of the Simulation class for use with the Depondt Solver. @@ -74,16 +74,16 @@ inline void Method_Solver::Iteration() // Get new spin conf n_new = R( (H+H')/2 ) * n Vectormath::rotate( conf, temp1, angle, conf ); } -}; +} template<> inline std::string Method_Solver::SolverName() { return "Depondt"; -}; +} template<> inline std::string Method_Solver::SolverFullName() { return "Depondt"; -}; \ No newline at end of file +} \ No newline at end of file diff --git a/core/include/engine/Solver_Heun.hpp b/core/include/engine/Solver_Heun.hpp index 3f33d89b1..01aa0b654 100644 --- a/core/include/engine/Solver_Heun.hpp +++ b/core/include/engine/Solver_Heun.hpp @@ -16,7 +16,7 @@ inline void Method_Solver::Initialize() configurations_predictor[i] = std::shared_ptr( new vectorfield( this->nos ) ); this->temp1 = vectorfield( this->nos, { 0, 0, 0 } ); -}; +} /* Template instantiation of the Simulation class for use with the Heun Solver. @@ -77,16 +77,16 @@ inline void Method_Solver::Iteration() // Copy out conf = conf_temp; } -}; +} template<> inline std::string Method_Solver::SolverName() { return "Heun"; -}; +} template<> inline std::string Method_Solver::SolverFullName() { return "Heun"; -}; \ No newline at end of file +} \ No newline at end of file diff --git a/core/include/engine/Solver_Kernels.hpp b/core/include/engine/Solver_Kernels.hpp index de55b6bd2..30603755a 100644 --- a/core/include/engine/Solver_Kernels.hpp +++ b/core/include/engine/Solver_Kernels.hpp @@ -17,6 +17,7 @@ namespace Engine { namespace Solver_Kernels { + // SIB void sib_transform( const vectorfield & spins, const vectorfield & force, vectorfield & out ); diff --git a/core/include/engine/Solver_RK4.hpp b/core/include/engine/Solver_RK4.hpp index 6c5cff621..ff1c48ef6 100644 --- a/core/include/engine/Solver_RK4.hpp +++ b/core/include/engine/Solver_RK4.hpp @@ -32,7 +32,7 @@ inline void Method_Solver::Initialize() this->configurations_k4[i] = std::shared_ptr( new vectorfield( this->nos ) ); this->temp1 = vectorfield( this->nos, { 0, 0, 0 } ); -}; +} /* Template instantiation of the Simulation class for use with the 4th order Runge Kutta Solver. @@ -144,16 +144,16 @@ inline void Method_Solver::Iteration() // Copy out conf = conf_temp; } -}; +} template<> inline std::string Method_Solver::SolverName() { return "RK4"; -}; +} template<> inline std::string Method_Solver::SolverFullName() { return "Runge Kutta (4th order)"; -}; \ No newline at end of file +} \ No newline at end of file diff --git a/core/include/engine/Solver_SIB.hpp b/core/include/engine/Solver_SIB.hpp index 7afe1d32d..f95a894aa 100644 --- a/core/include/engine/Solver_SIB.hpp +++ b/core/include/engine/Solver_SIB.hpp @@ -10,7 +10,7 @@ inline void Method_Solver::Initialize() this->configurations_predictor = std::vector>( this->noi ); for( int i = 0; i < this->noi; i++ ) configurations_predictor[i] = std::shared_ptr( new vectorfield( this->nos ) ); -}; +} /* Template instantiation of the Simulation class for use with the SIB Solver. @@ -47,16 +47,16 @@ inline void Method_Solver::Iteration() Solver_Kernels::sib_transform( image, forces_virtual_predictor[i], image ); } -}; +} template<> inline std::string Method_Solver::SolverName() { return "SIB"; -}; +} template<> inline std::string Method_Solver::SolverFullName() { return "Semi-implicit B"; -}; \ No newline at end of file +} \ No newline at end of file diff --git a/core/include/engine/Solver_VP.hpp b/core/include/engine/Solver_VP.hpp index 072d5861b..e7461487e 100644 --- a/core/include/engine/Solver_VP.hpp +++ b/core/include/engine/Solver_VP.hpp @@ -15,7 +15,7 @@ inline void Method_Solver::Initialize() this->forces_previous = velocities; // [noi][nos] this->projection = std::vector( this->noi, 0 ); // [noi] this->force_norm2 = std::vector( this->noi, 0 ); // [noi] -}; +} /* Template instantiation of the Simulation class for use with the VP Solver. @@ -111,16 +111,16 @@ inline void Method_Solver::Iteration() conf[idx] = conf_temp[idx].normalized(); } ); } -}; +} template<> inline std::string Method_Solver::SolverName() { return "VP"; -}; +} template<> inline std::string Method_Solver::SolverFullName() { return "Velocity Projection"; -}; \ No newline at end of file +} \ No newline at end of file diff --git a/core/include/engine/Solver_VP_OSO.hpp b/core/include/engine/Solver_VP_OSO.hpp index 3b4b2333e..f6132fac5 100644 --- a/core/include/engine/Solver_VP_OSO.hpp +++ b/core/include/engine/Solver_VP_OSO.hpp @@ -16,7 +16,7 @@ inline void Method_Solver::Initialize() this->projection = std::vector( this->noi, 0 ); // [noi] this->force_norm2 = std::vector( this->noi, 0 ); // [noi] this->searchdir = std::vector( this->noi, vectorfield( this->nos, { 0, 0, 0 } ) ); -}; +} /* Template instantiation of the Simulation class for use with the VP Solver. @@ -117,10 +117,10 @@ template<> inline std::string Method_Solver::SolverName() { return "VP_OSO"; -}; +} template<> inline std::string Method_Solver::SolverFullName() { return "Velocity Projection using exponential transforms"; -}; \ No newline at end of file +} \ No newline at end of file diff --git a/core/include/engine/Vectormath.hpp b/core/include/engine/Vectormath.hpp index e58979dc9..0a5b403a8 100644 --- a/core/include/engine/Vectormath.hpp +++ b/core/include/engine/Vectormath.hpp @@ -531,7 +531,13 @@ inline int idx_from_pair( //////// Vectorfield Math - special stuff // Calculate the mean of a vectorfield -std::array Magnetization( const vectorfield & vf ); +std::array Magnetization( const vectorfield & vf, const scalarfield & mu_s ); + +// Calculate the topological charge density inside a vectorfield +void TopologicalChargeDensity( + const vectorfield & vf, const Data::Geometry & geometry, const intfield & boundary_conditions, + scalarfield & charge_density, std::vector & triangle_indices ); + // Calculate the topological charge inside a vectorfield scalar TopologicalCharge( const vectorfield & vf, const Data::Geometry & geom, const intfield & boundary_conditions ); diff --git a/core/include/utility/Exception.hpp b/core/include/utility/Exception.hpp index b59c4acd6..344a2502b 100644 --- a/core/include/utility/Exception.hpp +++ b/core/include/utility/Exception.hpp @@ -6,6 +6,9 @@ #include +#include +#include + namespace Utility { @@ -29,18 +32,17 @@ enum class Exception_Classifier }; /* - * Spirit library exception class: - * Derived from std::runtime_error. - * Adds file, line and function information to the exception message. - * Contains an exception classifier and a level so that the handler - * can decide if execution should be stopped or can continue. + * Spirit library exception class, derived from std::runtime_error. + * + * Adds file, line and function information to the exception message. Contains an exception classifier and a level so + * that the handler can decide if execution should be stopped or can continue. */ class Exception : public std::runtime_error { public: Exception( Exception_Classifier classifier, Log_Level level, const std::string & message, const char * file, - unsigned int line, const char * function ) + unsigned int line, const char * function ) noexcept( false ) : std::runtime_error( fmt::format( "{}:{} in function \'{}\':\n{:>49}{}", file, line, function, " ", message ) ), classifier( classifier ), @@ -51,7 +53,17 @@ class Exception : public std::runtime_error { } - ~Exception() throw() override = default; + // Copy-constructor + Exception( const Exception & ) noexcept = default; + // Copy-assignment operator + Exception & operator=( const Exception & ) = delete; + + // Move-constructor + Exception( Exception && ) noexcept = default; + // Move-assignment constructor + Exception & operator=( Exception && ) = delete; + + ~Exception() noexcept override = default; const Exception_Classifier classifier; const Log_Level level; @@ -60,43 +72,58 @@ class Exception : public std::runtime_error const char * function; }; +// To make sure that the noexcept-specifications are correct +static_assert( + std::is_nothrow_copy_constructible::value, + "Utility::Exception is expected to be nothrow copy-constructible" ); +static_assert( + std::is_nothrow_move_constructible::value, + "Utility::Exception is expected to be nothrow move-constructible" ); +static_assert( + std::is_nothrow_destructible::value, "Utility::Exception is expected to be nothrow destructible" ); + /* * Rethrow (creating a std::nested_exception) an exception using the Exception class * to add file and line info */ -void rethrow( const std::string & message, const char * file, unsigned int line, const char * function ); +void rethrow( const std::string & message, const char * file, const unsigned int line, const char * function ) noexcept( + false ); /* * Handle_Exception_API finalizes what should be done when an exception is encountered at the API layer. - * This function should only be used inside API functions, since that is the top level at which an - * exception is caught. + * This function catches any exceptions and should only be used in the catch-statements of API functions, since that is + * the top level at which an exception is caught. */ void Handle_Exception_API( - const char * file, unsigned int line, const char * function = "", int idx_image = -1, int idx_chain = -1 ); + const char * file, const unsigned int line, const char * function = "", const int idx_image = -1, + const int idx_chain = -1 ) noexcept; /* - * Handle_Exception_Core finalizes what should be done when an exception is encountered inside the core. + * Handle_Exception_Core finalizes what should be done when an exception is encountered inside the core library. * This function should only be used inside the core, below the API layer. */ -void Handle_Exception_Core( const std::string & message, const char * file, unsigned int line, const char * function ); +void Handle_Exception_Core( + const std::string & message, const char * file, unsigned int line, const char * function ) noexcept( false ); // Shorthand for throwing a Spirit library exception with file and line info -// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) #define spirit_throw( classifier, level, message ) \ + /* NOLINTNEXTLINE(cppcoreguidelines-macro-usage,hicpp-no-array-decay) */ \ throw Utility::Exception( classifier, level, message, __FILE__, __LINE__, __func__ ) // Rethrow any exception to create a backtraceable nested exception -// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) -#define spirit_rethrow( message ) Utility::rethrow( message, __FILE__, __LINE__, __func__ ) +#define spirit_rethrow( message ) \ + /* NOLINTNEXTLINE(cppcoreguidelines-macro-usage,hicpp-no-array-decay) */ \ + Utility::rethrow( message, __FILE__, __LINE__, __func__ ) // Handle exception with backtrace and logging information on the calling API function -// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) #define spirit_handle_exception_api( idx_image, idx_chain ) \ + /* NOLINTNEXTLINE(cppcoreguidelines-macro-usage,hicpp-no-array-decay) */ \ Utility::Handle_Exception_API( __FILE__, __LINE__, __func__, idx_image, idx_chain ) // Handle exception with backtrace and logging information on the calling core function -// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) -#define spirit_handle_exception_core( message ) Utility::Handle_Exception_Core( message, __FILE__, __LINE__, __func__ ) +#define spirit_handle_exception_core( message ) \ + /* NOLINTNEXTLINE(cppcoreguidelines-macro-usage,hicpp-no-array-decay) */ \ + Utility::Handle_Exception_Core( message, __FILE__, __LINE__, __func__ ) } // namespace Utility diff --git a/core/python/setup.py b/core/python/setup.py index e9f827848..2c5d8d91c 100644 --- a/core/python/setup.py +++ b/core/python/setup.py @@ -135,7 +135,7 @@ def run(self): classifiers = CLASSIFIERS, install_requires = INSTALL_REQUIRES, package_data = { - 'spirit': ['libSpirit.dylib', 'libSpirit.so', 'libSpirit.dll'], + 'spirit': ['libSpirit.dylib', 'libSpirit.so', 'Spirit.dll'], }, cmdclass = {'bdist_wheel': bdist_wheel, 'clean': CleanCommand}, test_suite = 'setup.my_test_suite', diff --git a/core/python/spirit/parameters/ema.py b/core/python/spirit/parameters/ema.py index 6996d3ea9..1085286db 100644 --- a/core/python/spirit/parameters/ema.py +++ b/core/python/spirit/parameters/ema.py @@ -9,10 +9,19 @@ import spirit.spiritlib as spiritlib import ctypes -### Load Library +# Load Library _spirit = spiritlib.load_spirit_library() -## ---------------------------------- Set ---------------------------------- +_EMA_Clear_Modes = _spirit.Parameters_EMA_Clear_Modes +_EMA_Clear_Modes.argtypes = [ctypes.c_void_p, + ctypes.c_int, ctypes.c_int] +_EMA_Clear_Modes.restype = None +def clear_modes(p_state, idx_image=-1, idx_chain=-1): + """Clears the modes.""" + _EMA_Clear_Modes(ctypes.c_void_p(p_state), + ctypes.c_int(idx_image), ctypes.c_int(idx_chain)) + +# ---------------------------------- Set ---------------------------------- _EMA_Set_N_Modes = _spirit.Parameters_EMA_Set_N_Modes _EMA_Set_N_Modes.argtypes = [ctypes.c_void_p, ctypes.c_int, @@ -32,7 +41,16 @@ def set_n_mode_follow(p_state, n_mode, idx_image=-1, idx_chain=-1): _EMA_Set_N_Mode_Follow(ctypes.c_void_p(p_state), ctypes.c_int(n_mode), ctypes.c_int(idx_image), ctypes.c_int(idx_chain)) -## ---------------------------------- Get ---------------------------------- +_EMA_Set_Sparse = _spirit.Parameters_EMA_Set_Sparse +_EMA_Set_Sparse.argtypes = [ctypes.c_void_p, ctypes.c_bool, + ctypes.c_int, ctypes.c_int] +_EMA_Set_Sparse.restype = None +def set_sparse(p_state, sparse, idx_image=-1, idx_chain=-1): + """Set wether to use sparse matrices.""" + _EMA_Set_Sparse(ctypes.c_void_p(p_state), ctypes.c_bool(sparse), + ctypes.c_int(idx_image), ctypes.c_int(idx_chain)) + +# ---------------------------------- Get ---------------------------------- _EMA_Get_N_Modes = _spirit.Parameters_EMA_Get_N_Modes _EMA_Get_N_Modes.argtypes = [ctypes.c_void_p, ctypes.c_int, ctypes.c_int] @@ -46,4 +64,14 @@ def get_n_modes(p_state, idx_image=-1, idx_chain=-1): _EMA_Get_N_Mode_Follow.restype = ctypes.c_int def get_n_mode_follow(p_state, idx_image=-1, idx_chain=-1): """Returns the index of the mode to use.""" - return int(_EMA_Get_N_Mode_Follow(p_state, ctypes.c_int(idx_image), ctypes.c_int(idx_chain))) \ No newline at end of file + return int(_EMA_Get_N_Mode_Follow(p_state, ctypes.c_int(idx_image), ctypes.c_int(idx_chain))) + + +_EMA_Get_Sparse = _spirit.Parameters_EMA_Get_Sparse +_EMA_Get_Sparse.argtypes = [ctypes.c_void_p, ctypes.c_bool, + ctypes.c_int, ctypes.c_int] +_EMA_Get_Sparse.restype = ctypes.c_bool +def get_sparse(p_state, sparse, idx_image=-1, idx_chain=-1): + """Set wether to use sparse matrices.""" + return bool(_EMA_Get_Sparse(ctypes.c_void_p(p_state), + ctypes.c_int(idx_image), ctypes.c_int(idx_chain))) diff --git a/core/python/spirit/parameters/gneb.py b/core/python/spirit/parameters/gneb.py index bbe549498..1c46dc60d 100644 --- a/core/python/spirit/parameters/gneb.py +++ b/core/python/spirit/parameters/gneb.py @@ -135,6 +135,27 @@ def set_path_shortening_constant(p_state, shortening_constant, idx_image=-1, idx _GNEB_Set_Path_Shortening_Constant(ctypes.c_void_p(p_state), ctypes.c_float(shortening_constant), ctypes.c_int(idx_image), ctypes.c_int(idx_chain)) +_GNEB_Set_Moving_Endpoints = _spirit.Parameters_GNEB_Set_Moving_Endpoints +_GNEB_Set_Moving_Endpoints.argtypes = [ctypes.c_void_p, ctypes.c_bool, ctypes.c_int] +_GNEB_Set_Moving_Endpoints.restype = None +def set_moving_endpoints(p_state, moving_endpoints, idx_chain=-1): + """Set if moving endpoints should be used.""" + _GNEB_Set_Moving_Endpoints(ctypes.c_void_p(p_state), ctypes.c_bool(moving_endpoints), ctypes.c_int(idx_chain)) + +_GNEB_Set_Translating_Endpoints = _spirit.Parameters_GNEB_Set_Translating_Endpoints +_GNEB_Set_Translating_Endpoints.argtypes = [ctypes.c_void_p, ctypes.c_bool, ctypes.c_int] +_GNEB_Set_Translating_Endpoints.restype = None +def set_translating_endpoints(p_state, translating_endpoints, idx_chain=-1): + """Set if attracting endpoints should be used.""" + _GNEB_Set_Translating_Endpoints(ctypes.c_void_p(p_state), ctypes.c_bool(translating_endpoints), ctypes.c_int(idx_chain)) + +_GNEB_Set_Equilibrium_Delta_Rx = _spirit.Parameters_GNEB_Set_Equilibrium_Delta_Rx +_GNEB_Set_Equilibrium_Delta_Rx.argtypes = [ctypes.c_void_p, ctypes.c_float, ctypes.c_float, ctypes.c_int] +_GNEB_Set_Equilibrium_Delta_Rx.restype = None +def set_equilibrium_delta_Rx(p_state, delta_Rx_left, delta_Rx_right, idx_chain=-1): + """Set if moving endpoints should be used.""" + _GNEB_Set_Equilibrium_Delta_Rx(ctypes.c_void_p(p_state), ctypes.c_float(delta_Rx_left), ctypes.c_float(delta_Rx_right), ctypes.c_int(idx_chain)) + _GNEB_Set_Climbing_Falling = _spirit.Parameters_GNEB_Set_Climbing_Falling _GNEB_Set_Climbing_Falling.argtypes = [ctypes.c_void_p, ctypes.c_int, ctypes.c_int, ctypes.c_int] _GNEB_Set_Climbing_Falling.restype = None @@ -198,6 +219,34 @@ def get_path_shortening_constant(p_state, idx_image=-1, idx_chain=-1): return float( _GNEB_Get_Path_Shortening_Constant(ctypes.c_void_p(p_state), ctypes.c_int(idx_image), ctypes.c_int(idx_chain))) +_GNEB_Get_Moving_Endpoints = _spirit.Parameters_GNEB_Get_Moving_Endpoints +_GNEB_Get_Moving_Endpoints.argtypes = [ctypes.c_void_p, ctypes.c_int] +_GNEB_Get_Moving_Endpoints.restype = ctypes.c_bool +def get_moving_endpoints(p_state, idx_chain=-1): + """Return if moving endpoints are used.""" + return bool( _GNEB_Get_Moving_Endpoints(ctypes.c_void_p(p_state), + ctypes.c_int(idx_chain)) ) + +_GNEB_Get_Translating_Endpoints = _spirit.Parameters_GNEB_Get_Translating_Endpoints +_GNEB_Get_Translating_Endpoints.argtypes = [ctypes.c_void_p, ctypes.c_int] +_GNEB_Get_Translating_Endpoints.restype = ctypes.c_bool +def get_translating_endpoints(p_state, idx_chain=-1): + """Return if translating endpoints are used.""" + return bool( _GNEB_Get_Translating_Endpoints(ctypes.c_void_p(p_state), + ctypes.c_int(idx_chain)) ) + +_GNEB_Get_Equilibrium_Delta_Rx = _spirit.Parameters_GNEB_Get_Equilibrium_Delta_Rx +_GNEB_Get_Equilibrium_Delta_Rx.argtypes = [ctypes.c_void_p, ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_float), ctypes.c_int] +_GNEB_Get_Equilibrium_Delta_Rx.restype = None +def get_equilibrium_delta_Rx(p_state, idx_chain=-1): + """Return the equilibrium delta_Rx for the moving endpoints.""" + delta_Rx_left = ctypes.c_float() + delta_Rx_right = ctypes.c_float() + + _GNEB_Get_Equilibrium_Delta_Rx(ctypes.c_void_p(p_state), ctypes.byref(delta_Rx_left), ctypes.byref(delta_Rx_right), ctypes.c_int(idx_chain)) + + return [float(delta_Rx_left.value), float(delta_Rx_right.value)] + _GNEB_Get_Climbing_Falling = _spirit.Parameters_GNEB_Get_Climbing_Falling _GNEB_Get_Climbing_Falling.argtypes = [ctypes.c_void_p, ctypes.c_int, ctypes.c_int] _GNEB_Get_Climbing_Falling.restype = ctypes.c_int diff --git a/core/python/spirit/quantities.py b/core/python/spirit/quantities.py index 7e72d9d5b..458d6aa7e 100644 --- a/core/python/spirit/quantities.py +++ b/core/python/spirit/quantities.py @@ -83,4 +83,30 @@ def get_topological_charge(p_state, idx_image=-1, idx_chain=-1): Returns 0 for systems of other dimensionality. """ return float(_Get_Topological_Charge(ctypes.c_void_p(p_state), - ctypes.c_int(idx_image), ctypes.c_int(idx_chain))) \ No newline at end of file + ctypes.c_int(idx_image), ctypes.c_int(idx_chain))) + + +_Get_Topological_Charge_Density = _spirit.Quantity_Get_Topological_Charge_Density +_Get_Topological_Charge_Density.argtypes = [ctypes.c_void_p, ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_int), ctypes.c_int, ctypes.c_int] +_Get_Topological_Charge_Density.restype = ctypes.c_int +def get_topological_charge_density(p_state, idx_image=-1, idx_chain=-1): + """Calculates and returns the total topological charge of 2D systems. + + Note that the charge can take unphysical non-integer values for open boundaries, + because it is not well-defined in this case. + + Returns 0 for systems of other dimensionality. + """ + + # Figure out the number of triangles + num_triangles = _Get_Topological_Charge_Density(ctypes.c_void_p(p_state), None, None, + ctypes.c_int(idx_image), ctypes.c_int(idx_chain)) + + # Allocate the required memory for indices_ptr and charge_density + charge_density_ptr = (ctypes.c_float * num_triangles)() + triangle_indices_ptr = (ctypes.c_int * (3 * num_triangles))() + + _Get_Topological_Charge_Density(ctypes.c_void_p(p_state), charge_density_ptr, triangle_indices_ptr, + ctypes.c_int(idx_image), ctypes.c_int(idx_chain)) + + return [float(c) for c in charge_density_ptr], [ [int(triangle_indices_ptr[3*i]), int(triangle_indices_ptr[3*i + 1]), int(triangle_indices_ptr[3*i + 2])] for i in range(num_triangles)] \ No newline at end of file diff --git a/core/python/spirit/simulation.py b/core/python/spirit/simulation.py index 601a2df35..87e2be3f0 100644 --- a/core/python/spirit/simulation.py +++ b/core/python/spirit/simulation.py @@ -94,9 +94,27 @@ class simulation_run_info(ctypes.Structure): ("total_iterations", ctypes.c_int), ("total_walltime", ctypes.c_int), ("total_ips", ctypes.c_float), - ("max_torque", ctypes.c_float) + ("max_torque", ctypes.c_float), + ("_n_history_iteration", ctypes.c_int), + ("_history_iteration", ctypes.POINTER(ctypes.c_int)), + ("_n_history_max_torque", ctypes.c_int), + ("_history_max_torque", ctypes.POINTER(ctypes.c_float)), + ("_n_history_energy", ctypes.c_int), + ("_history_energy", ctypes.POINTER(ctypes.c_float)) ] + def history_iteration(self): + return [ self._history_iteration[i] for i in range(self._n_history_iteration) ] + + def history_max_torque(self): + return [ self._history_max_torque[i] for i in range(self._n_history_max_torque) ] + + def history_energy(self): + return [ self._history_energy[i] for i in range(self._n_history_energy) ] + + def __del__(self): + _spirit.free_run_info(self) + ### ----- Start methods ### MC _MC_Start = _spirit.Simulation_MC_Start diff --git a/core/python/test/quantities.py b/core/python/test/quantities.py index 0eff6bfe6..3bc3c6691 100644 --- a/core/python/test/quantities.py +++ b/core/python/test/quantities.py @@ -5,13 +5,13 @@ spirit_py_dir = os.path.abspath(os.path.join(os.path.dirname( __file__ ), "..")) sys.path.insert(0, spirit_py_dir) -from spirit import state, quantities, configuration +from spirit import state, quantities, configuration, geometry import unittest ########## -cfgfile = spirit_py_dir + "/../test/input/fd_neighbours.cfg" # Input File +cfgfile = spirit_py_dir + "/../test/input/api.cfg" # Input File p_state = state.setup(cfgfile) # State setup @@ -25,11 +25,21 @@ class Quantities_Get(TestParameters): def test_magnetization(self): configuration.plus_z(self.p_state) + mu_s = 1.34 + geometry.set_mu_s(self.p_state, mu_s) M = quantities.get_magnetization(self.p_state) self.assertAlmostEqual(M[0], 0) self.assertAlmostEqual(M[1], 0) - self.assertAlmostEqual(M[2], 1) - + self.assertAlmostEqual(M[2], mu_s) + + def test_topological_charge(self): + configuration.plus_z(self.p_state) + configuration.skyrmion(self.p_state, radius=5, pos=[1.5, 0, 0]) + Q = quantities.get_topological_charge(self.p_state) + [Q_density, triangles] = quantities.get_topological_charge_density(self.p_state) + self.assertAlmostEqual(Q, -1.0) + self.assertAlmostEqual(Q, sum(Q_density)) + ######### def suite(): diff --git a/core/src/Spirit/IO.cpp b/core/src/Spirit/IO.cpp index 9ad9bac82..976b3f666 100644 --- a/core/src/Spirit/IO.cpp +++ b/core/src/Spirit/IO.cpp @@ -1052,7 +1052,7 @@ try Log( Utility::Log_Level::Warning, Utility::Log_Sender::IO, fmt::format( "Resizing eigenmode buffer because the number of modes in the OVF file ({}) " - "is greater than the buffer size ({})", + "is different from the buffer size ({})", n_eigenmodes, image->modes.size() ) ); image->modes.resize( n_eigenmodes ); image->eigenvalues.resize( n_eigenmodes ); diff --git a/core/src/Spirit/Parameters_EMA.cpp b/core/src/Spirit/Parameters_EMA.cpp index c0178dbbe..922010f0e 100644 --- a/core/src/Spirit/Parameters_EMA.cpp +++ b/core/src/Spirit/Parameters_EMA.cpp @@ -7,6 +7,25 @@ #include +// Clears all the previously calculated modes from memory +void Parameters_EMA_Clear_Modes( State * state, int idx_image, int idx_chain ) noexcept +{ + std::shared_ptr image; + std::shared_ptr chain; + + // Fetch correct indices and pointers + from_indices( state, idx_image, idx_chain, image, chain ); + + Log( Utility::Log_Level::Info, Utility::Log_Sender::API, "Clearing modes", idx_image, idx_chain ); + + image->Lock(); + for( auto & el : image->modes ) + { + el.reset(); + } + image->Unlock(); +} + /*------------------------------------------------------------------------------------------------------ */ /*---------------------------------- Set EMA ---------------------------------------------------------- */ /*------------------------------------------------------------------------------------------------------ */ @@ -123,6 +142,26 @@ catch( ... ) spirit_handle_exception_api( idx_image, idx_chain ); } +void Parameters_EMA_Set_Sparse( State * state, bool sparse, int idx_image, int idx_chain ) noexcept +try +{ + std::shared_ptr image; + std::shared_ptr chain; + + // Fetch correct indices and pointers + from_indices( state, idx_image, idx_chain, image, chain ); + + Log( Utility::Log_Level::Info, Utility::Log_Sender::API, fmt::format( "Setting parameter 'sparse' to {}", sparse ), + idx_image, idx_chain ); + image->Lock(); + image->ema_parameters->sparse = sparse; + image->Unlock(); +} +catch( ... ) +{ + spirit_handle_exception_api( idx_image, idx_chain ); +} + /*------------------------------------------------------------------------------------------------------ */ /*---------------------------------- Get EMA ----------------------------------------------------------- */ /*------------------------------------------------------------------------------------------------------ */ @@ -208,6 +247,23 @@ try return image->ema_parameters->snapshot; } catch( ... ) +{ + spirit_handle_exception_api( idx_image, idx_chain ); + return false; +} + +bool Parameters_EMA_Get_Sparse( State * state, int idx_image, int idx_chain ) noexcept +try +{ + std::shared_ptr image; + std::shared_ptr chain; + + // Fetch correct indices and pointers + from_indices( state, idx_image, idx_chain, image, chain ); + + return image->ema_parameters->sparse; +} +catch( ... ) { spirit_handle_exception_api( idx_image, idx_chain ); return false; diff --git a/core/src/Spirit/Parameters_GNEB.cpp b/core/src/Spirit/Parameters_GNEB.cpp index d2e22118b..27e9a500b 100644 --- a/core/src/Spirit/Parameters_GNEB.cpp +++ b/core/src/Spirit/Parameters_GNEB.cpp @@ -237,6 +237,80 @@ catch( ... ) spirit_handle_exception_api( -1, idx_chain ); } +// Set if moving endpoints should be used +void Parameters_GNEB_Set_Moving_Endpoints( State * state, bool moving_endpoints, int idx_chain ) noexcept +try +{ + int idx_image = -1; + std::shared_ptr image; + std::shared_ptr chain; + + // Fetch correct indices and pointers + from_indices( state, idx_image, idx_chain, image, chain ); + + chain->Lock(); + auto p = chain->gneb_parameters; + p->moving_endpoints = moving_endpoints; + chain->Unlock(); + + Log( Utility::Log_Level::Parameter, Utility::Log_Sender::API, + fmt::format( "Set GNEB moving endpoints = {}", moving_endpoints ), idx_image, idx_chain ); +} +catch( ... ) +{ + spirit_handle_exception_api( -1, idx_chain ); +} + +// Set if translating endpoints should be used +void Parameters_GNEB_Set_Translating_Endpoints( State * state, bool translating_endpoints, int idx_chain ) noexcept +try +{ + int idx_image = -1; + std::shared_ptr image; + std::shared_ptr chain; + + // Fetch correct indices and pointers + from_indices( state, idx_image, idx_chain, image, chain ); + + chain->Lock(); + auto p = chain->gneb_parameters; + p->translating_endpoints = translating_endpoints; + chain->Unlock(); + + Log( Utility::Log_Level::Parameter, Utility::Log_Sender::API, + fmt::format( "Set GNEB translating endpoints = {}", translating_endpoints ), idx_image, idx_chain ); +} +catch( ... ) +{ + spirit_handle_exception_api( -1, idx_chain ); +} + +void Parameters_GNEB_Set_Equilibrium_Delta_Rx( State * state, float delta_Rx_left, float delta_Rx_right, int idx_chain ) noexcept +try +{ + int idx_image = -1; + std::shared_ptr image; + std::shared_ptr chain; + + // Fetch correct indices and pointers + from_indices( state, idx_image, idx_chain, image, chain ); + + chain->Lock(); + auto p = chain->gneb_parameters; + p->equilibrium_delta_Rx_left = delta_Rx_left; + p->equilibrium_delta_Rx_right = delta_Rx_right; + + chain->Unlock(); + + Log( Utility::Log_Level::Parameter, Utility::Log_Sender::API, + fmt::format( "Set equilibrium delta Rx for GNEB with moving endpoints. delta_Rx_left = {}, delta_Rx_right = {}", delta_Rx_left, delta_Rx_right ), idx_image, idx_chain ); +} +catch( ... ) +{ + spirit_handle_exception_api( -1, idx_chain ); +} + + void Parameters_GNEB_Set_Climbing_Falling( State * state, int image_type, int idx_image, int idx_chain ) noexcept try { @@ -507,6 +581,63 @@ catch( ... ) return 0; } +bool Parameters_GNEB_Get_Moving_Endpoints( State * state, int idx_chain ) noexcept +try +{ + int idx_image = -1; + std::shared_ptr image; + std::shared_ptr chain; + + // Fetch correct indices and pointers + from_indices( state, idx_image, idx_chain, image, chain ); + + auto p = chain->gneb_parameters; + return static_cast( p->moving_endpoints ); +} +catch( ... ) +{ + spirit_handle_exception_api( -1, idx_chain ); + return 0; +} + +bool Parameters_GNEB_Get_Translating_Endpoints( State * state, int idx_chain ) noexcept +try +{ + int idx_image = -1; + std::shared_ptr image; + std::shared_ptr chain; + + // Fetch correct indices and pointers + from_indices( state, idx_image, idx_chain, image, chain ); + + auto p = chain->gneb_parameters; + return static_cast( p->translating_endpoints ); +} +catch( ... ) +{ + spirit_handle_exception_api( -1, idx_chain ); + return 0; +} + +void Parameters_GNEB_Get_Equilibrium_Delta_Rx( State * state, float * delta_Rx_left, float * delta_Rx_right, int idx_chain) noexcept +try +{ + int idx_image = -1; + std::shared_ptr image; + std::shared_ptr chain; + + // Fetch correct indices and pointers + from_indices( state, idx_image, idx_chain, image, chain ); + + auto p = chain->gneb_parameters; + *delta_Rx_left = float(p->equilibrium_delta_Rx_left); + *delta_Rx_right = float(p->equilibrium_delta_Rx_right); +} +catch( ... ) +{ + spirit_handle_exception_api( -1, idx_chain ); +} + int Parameters_GNEB_Get_Climbing_Falling( State * state, int idx_image, int idx_chain ) noexcept try { diff --git a/core/src/Spirit/Quantities.cpp b/core/src/Spirit/Quantities.cpp index 8642e6f53..b621104e7 100644 --- a/core/src/Spirit/Quantities.cpp +++ b/core/src/Spirit/Quantities.cpp @@ -14,6 +14,27 @@ namespace C = Utility::Constants; +void Quantity_Get_Average_Spin( State * state, float s[3], int idx_image, int idx_chain ) +try +{ + std::shared_ptr image; + std::shared_ptr chain; + + // Fetch correct indices and pointers + from_indices( state, idx_image, idx_chain, image, chain ); + + // image->Lock(); // Mutex locks in these functions may cause problems with the performance of UIs + + auto mean = Engine::Vectormath::mean( *image->spins ); + + for( int i = 0; i < 3; ++i ) + s[i] = (float)mean[i]; +} +catch( ... ) +{ + spirit_handle_exception_api( idx_image, idx_chain ); +} + void Quantity_Get_Magnetization( State * state, float m[3], int idx_image, int idx_chain ) try { @@ -25,7 +46,7 @@ try // image->Lock(); // Mutex locks in these functions may cause problems with the performance of UIs - auto mag = Engine::Vectormath::Magnetization( *image->spins ); + auto mag = Engine::Vectormath::Magnetization( *image->spins, image->geometry->mu_s ); image->M = Vector3{ mag[0], mag[1], mag[2] }; // image->Unlock(); @@ -65,6 +86,51 @@ catch( ... ) return 0; } +int Quantity_Get_Topological_Charge_Density( + State * state, float * charge_density_ptr, int * triangle_indices_ptr, int idx_image, int idx_chain ) +try +{ + std::shared_ptr image; + std::shared_ptr chain; + + // Fetch correct indices and pointers + from_indices( state, idx_image, idx_chain, image, chain ); + + // image->Lock(); // Mutex locks in these functions may cause problems with the performance of UIs + + scalar charge = 0; + int dimensionality = Geometry_Get_Dimensionality( state, idx_image, idx_chain ); + scalarfield charge_density( 0 ); + std::vector triangle_indices( 0 ); + + if( dimensionality == 2 ) + { + Engine::Vectormath::TopologicalChargeDensity( + *image->spins, *image->geometry, image->hamiltonian->boundary_conditions, charge_density, + triangle_indices ); + } + + if( charge_density_ptr != nullptr && triangle_indices_ptr != nullptr ) + { + for( int i = 0; i < charge_density.size(); i++ ) + { + charge_density_ptr[i] = charge_density[i]; + triangle_indices_ptr[3 * i] = triangle_indices[3 * i]; + triangle_indices_ptr[3 * i + 1] = triangle_indices[3 * i + 1]; + triangle_indices_ptr[3 * i + 2] = triangle_indices[3 * i + 2]; + } + } + + // image->Unlock(); + + return charge_density.size(); // return the number of triangles +} +catch( ... ) +{ + spirit_handle_exception_api( idx_image, idx_chain ); + return 0; +} + void check_modes( const vectorfield & image, const vectorfield & grad, const MatrixX & tangent_basis, const VectorX & eigenvalues, const MatrixX & eigenvectors_2N, const vectorfield & minimum_mode ) diff --git a/core/src/Spirit/Simulation.cpp b/core/src/Spirit/Simulation.cpp index 23bbdc67f..c6e68d7c6 100644 --- a/core/src/Spirit/Simulation.cpp +++ b/core/src/Spirit/Simulation.cpp @@ -11,6 +11,15 @@ #include #include +#include + +void free_run_info(Simulation_Run_Info info) noexcept +{ + delete[] info.history_energy; + delete[] info.history_iteration; + delete[] info.history_max_torque; +}; + // Helper function to start a simulation once a Method has been created void run_method( std::shared_ptr method, bool singleshot, Simulation_Run_Info * info = nullptr ) { @@ -38,6 +47,27 @@ void run_method( std::shared_ptr method, bool singleshot, Simula info->total_iterations = method->getNIterations(); info->total_walltime = method->getWallTime(); info->total_ips = float( info->total_iterations ) / info->total_walltime * 1000.0; + + if( method->history_iteration.size() > 0 ) + { + info->n_history_iteration = method->history_iteration.size(); + info->history_iteration = new int[method->history_iteration.size()]; + std::copy(method->history_iteration.begin(), method->history_iteration.end(), info->history_iteration); + } + + if( method->history_max_torque.size() > 0 ) + { + info->n_history_max_torque = method->history_max_torque.size(); + info->history_max_torque = new float[method->history_max_torque.size()]; + std::copy(method->history_max_torque.begin(), method->history_max_torque.end(), info->history_max_torque); + } + + if( method->history_energy.size() > 0 ) + { + info->n_history_energy = method->history_energy.size(); + info->history_energy = new float[method->history_energy.size()]; + std::copy(method->history_energy.begin(), method->history_energy.end(), info->history_energy); + } } } } @@ -225,10 +255,10 @@ try std::string( "There are still one or more simulations running on the specified chain!" ) + std::string( " Please stop them before starting a GNEB calculation." ) ); } - else if( Chain_Get_NOI( state, idx_chain ) < 3 ) + else if( Chain_Get_NOI( state, idx_chain ) < 2 ) { Log( Utility::Log_Level::Error, Utility::Log_Sender::API, - std::string( "There are less than 3 images in the specified chain!" ) + std::string( "There are less than 2 images in the specified chain!" ) + std::string( " Please insert more before starting a GNEB calculation." ) ); } else diff --git a/core/src/data/Geometry.cpp b/core/src/data/Geometry.cpp index deb95fe39..4d9f5d715 100644 --- a/core/src/data/Geometry.cpp +++ b/core/src/data/Geometry.cpp @@ -21,8 +21,9 @@ namespace Data { Geometry::Geometry( - std::vector bravais_vectors, intfield n_cells, std::vector cell_atoms, - Basis_Cell_Composition cell_composition, scalar lattice_constant, Pinning pinning, Defects defects ) + const std::vector & bravais_vectors, intfield n_cells, const std::vector & cell_atoms, + const Basis_Cell_Composition & cell_composition, scalar lattice_constant, const Pinning & pinning, + const Defects & defects ) : bravais_vectors( bravais_vectors ), n_cells( n_cells ), n_cell_atoms( cell_atoms.size() ), diff --git a/core/src/engine/Eigenmodes.cpp b/core/src/engine/Eigenmodes.cpp index 98d5459dc..442768225 100644 --- a/core/src/engine/Eigenmodes.cpp +++ b/core/src/engine/Eigenmodes.cpp @@ -3,7 +3,9 @@ #include // #include -#include // Also includes +#include // Also includes +#include + #include #include @@ -66,7 +68,6 @@ void Calculate_Eigenmodes( std::shared_ptr system, int idx_im // Calculate the Eigenmodes vectorfield gradient( nos ); - MatrixX hessian( 3 * nos, 3 * nos ); // The gradient (unprojected) system->hamiltonian->Gradient( spins_initial, gradient ); @@ -77,17 +78,39 @@ void Calculate_Eigenmodes( std::shared_ptr system, int idx_im // }); Vectormath::set_c_a( 1, gradient, gradient, system->geometry->mask_unpinned ); - // The Hessian (unprojected) - system->hamiltonian->Hessian( spins_initial, hessian ); - - // Get the eigenspectrum - MatrixX hessian_constrained = MatrixX::Zero( 2 * nos, 2 * nos ); - MatrixX tangent_basis = MatrixX::Zero( 3 * nos, 2 * nos ); VectorX eigenvalues; MatrixX eigenvectors; - bool successful = Eigenmodes::Hessian_Partial_Spectrum( - system->ema_parameters, spins_initial, gradient, hessian, n_modes, tangent_basis, hessian_constrained, - eigenvalues, eigenvectors ); + SpMatrixX tangent_basis = SpMatrixX( 3 * nos, 2 * nos ); + + bool sparse = system->ema_parameters->sparse; + bool successful; + if( sparse ) + { + // The Hessian (unprojected) + SpMatrixX hessian( 3 * nos, 3 * nos ); + system->hamiltonian->Sparse_Hessian( spins_initial, hessian ); + // Get the eigenspectrum + SpMatrixX hessian_constrained = SpMatrixX( 2 * nos, 2 * nos ); + + successful = Eigenmodes::Sparse_Hessian_Partial_Spectrum( + system->ema_parameters, spins_initial, 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 ); + // Get the eigenspectrum + MatrixX hessian_constrained = MatrixX::Zero( 2 * nos, 2 * nos ); + MatrixX _tangent_basis = MatrixX( tangent_basis ); + + successful = Eigenmodes::Hessian_Partial_Spectrum( + system->ema_parameters, spins_initial, gradient, hessian, n_modes, _tangent_basis, hessian_constrained, + eigenvalues, eigenvectors ); + + tangent_basis = _tangent_basis.sparseView(); + } if( successful ) { @@ -217,5 +240,49 @@ bool Hessian_Partial_Spectrum( return ( hessian_spectrum.info() == Spectra::SUCCESSFUL ) && ( nconv > 0 ); } +bool Sparse_Hessian_Partial_Spectrum( + const std::shared_ptr parameters, const vectorfield & spins, const vectorfield & gradient, + const SpMatrixX & hessian, int n_modes, SpMatrixX & tangent_basis, SpMatrixX & hessian_constrained, + VectorX & eigenvalues, MatrixX & eigenvectors ) +{ + int nos = spins.size(); + + // Restrict number of calculated modes to [1,2N) + n_modes = std::max( 1, std::min( 2 * nos - 2, n_modes ) ); + + // Calculate the final Hessian to use for the minimum mode + Manifoldmath::sparse_tangent_basis_spherical( spins, tangent_basis ); + + SpMatrixX hessian_constrained_3N = SpMatrixX( 3 * nos, 3 * nos ); + Manifoldmath::sparse_hessian_bordered_3N( spins, gradient, hessian, hessian_constrained_3N ); + + hessian_constrained = tangent_basis.transpose() * hessian_constrained_3N * tangent_basis; + + // TODO: Pinning (see non-sparse function for) + + hessian_constrained.makeCompressed(); + int ncv = std::min(2*nos, std::max(2*n_modes + 1, 20)); // This is the default value used by scipy.sparse + int max_iter = 20*nos; + + // Create the Spectra Matrix product operation + Spectra::SparseSymMatProd op( hessian_constrained ); + // Create and initialize a Spectra solver + Spectra::SymEigsSolver> hessian_spectrum( + &op, n_modes, ncv); + hessian_spectrum.init(); + + // Compute the specified spectrum, sorted by smallest real eigenvalue + int nconv = hessian_spectrum.compute( max_iter, 1e-10, int( Spectra::SMALLEST_ALGE ) ); + + // Extract real eigenvalues + eigenvalues = hessian_spectrum.eigenvalues().real(); + + // Retrieve the real eigenvectors + eigenvectors = hessian_spectrum.eigenvectors().real(); + + // Return whether the calculation was successful + return ( hessian_spectrum.info() == Spectra::SUCCESSFUL ) && ( nconv > 0 ); +} + } // namespace Eigenmodes } // namespace Engine diff --git a/core/src/engine/Hamiltonian.cpp b/core/src/engine/Hamiltonian.cpp index a3aa9cb05..a596cd13a 100644 --- a/core/src/engine/Hamiltonian.cpp +++ b/core/src/engine/Hamiltonian.cpp @@ -177,12 +177,4 @@ scalar Hamiltonian::Energy_Single_Spin( int ispin, const vectorfield & spins ) "Tried to use Hamiltonian::Energy_Single_Spin() of the Hamiltonian base class!" ); } -static const std::string name = "--"; -const std::string & Hamiltonian::Name() -{ - spirit_throw( - Exception_Classifier::Not_Implemented, Log_Level::Error, - "Tried to use Hamiltonian::Name() of the Hamiltonian base class!" ); -} - } // namespace Engine diff --git a/core/src/engine/Hamiltonian_Gaussian.cpp b/core/src/engine/Hamiltonian_Gaussian.cpp index 39c6f07a2..2239bca8b 100644 --- a/core/src/engine/Hamiltonian_Gaussian.cpp +++ b/core/src/engine/Hamiltonian_Gaussian.cpp @@ -123,7 +123,7 @@ scalar Hamiltonian_Gaussian::Energy_Single_Spin( int ispin, const vectorfield & // Hamiltonian name as string static const std::string name = "Gaussian"; -const std::string & Hamiltonian_Gaussian::Name() +const std::string & Hamiltonian_Gaussian::Name() const { return name; } diff --git a/core/src/engine/Hamiltonian_Heisenberg.cpp b/core/src/engine/Hamiltonian_Heisenberg.cpp index 9d99f1603..2419163bc 100644 --- a/core/src/engine/Hamiltonian_Heisenberg.cpp +++ b/core/src/engine/Hamiltonian_Heisenberg.cpp @@ -855,7 +855,7 @@ void Hamiltonian_Heisenberg::Gradient_DDI_FFT( const vectorfield & spins, vector const int * c_it_bounds_pointwise_mult = it_bounds_pointwise_mult.data(); // Loop over basis atoms (i.e sublattices) -#pragma omp parallel for collapse( 5 ) +#pragma omp parallel for collapse( 4 ) for( int i_b1 = 0; i_b1 < c_n_cell_atoms; ++i_b1 ) { for( int c = 0; c < c_it_bounds_pointwise_mult[2]; ++c ) @@ -1517,7 +1517,7 @@ void Hamiltonian_Heisenberg::Clean_DDI() // Hamiltonian name as string static const std::string name = "Heisenberg"; -const std::string & Hamiltonian_Heisenberg::Name() +const std::string & Hamiltonian_Heisenberg::Name() const { return name; } diff --git a/core/src/engine/Hamiltonian_Heisenberg.cu b/core/src/engine/Hamiltonian_Heisenberg.cu index ed3153fd0..a4b0c4681 100644 --- a/core/src/engine/Hamiltonian_Heisenberg.cu +++ b/core/src/engine/Hamiltonian_Heisenberg.cu @@ -1413,7 +1413,10 @@ namespace Engine // Hamiltonian name as string static const std::string name = "Heisenberg"; - const std::string& Hamiltonian_Heisenberg::Name() { return name; } + const std::string & Hamiltonian_Heisenberg::Name() const + { + return name; + } } #endif diff --git a/core/src/engine/Hamiltonian_Micromagnetic.cpp b/core/src/engine/Hamiltonian_Micromagnetic.cpp index 49125611e..9ce57091a 100644 --- a/core/src/engine/Hamiltonian_Micromagnetic.cpp +++ b/core/src/engine/Hamiltonian_Micromagnetic.cpp @@ -912,7 +912,7 @@ void Hamiltonian_Micromagnetic::Hessian( const vectorfield & spins, MatrixX & he // Hamiltonian name as string static const std::string name = "Micromagnetic"; -const std::string & Hamiltonian_Micromagnetic::Name() +const std::string & Hamiltonian_Micromagnetic::Name() const { return name; } diff --git a/core/src/engine/Hamiltonian_Micromagnetic.cu b/core/src/engine/Hamiltonian_Micromagnetic.cu index 50f351b10..7bb001ba1 100644 --- a/core/src/engine/Hamiltonian_Micromagnetic.cu +++ b/core/src/engine/Hamiltonian_Micromagnetic.cu @@ -1075,7 +1075,7 @@ void Hamiltonian_Micromagnetic::Hessian( const vectorfield & spins, MatrixX & he // Hamiltonian name as string static const std::string name = "Micromagnetic"; -const std::string & Hamiltonian_Micromagnetic::Name() +const std::string & Hamiltonian_Micromagnetic::Name() const { return name; } diff --git a/core/src/engine/Manifoldmath.cpp b/core/src/engine/Manifoldmath.cpp index e9c55f934..9ecad2461 100644 --- a/core/src/engine/Manifoldmath.cpp +++ b/core/src/engine/Manifoldmath.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include @@ -21,12 +22,12 @@ namespace Manifoldmath { void project_parallel( vectorfield & vf1, const vectorfield & vf2 ) { - vectorfield vf3 = vf1; - project_orthogonal( vf3, vf2 ); -// TODO: replace the loop with Vectormath Kernel -#pragma omp parallel for - for( unsigned int i = 0; i < vf1.size(); ++i ) - vf1[i] -= vf3[i]; + scalar proj = Vectormath::dot(vf1, vf2); + Backend::par::apply( vf1.size(), [vf1 = vf1.data(), vf2 = vf2.data(), proj] SPIRIT_LAMBDA (int idx) + { + vf1[idx] = proj * vf2[idx]; + } + ); } void project_orthogonal( vectorfield & vf1, const vectorfield & vf2 ) @@ -73,6 +74,42 @@ scalar dist_geodesic( const vectorfield & v1, const vectorfield & v2 ) return sqrt( dist ); } +/* + Helper function for a more accurate tangent +*/ +void Geodesic_Tangent( vectorfield & tangent, const vectorfield & image_1, const vectorfield & image_2, const vectorfield & image_mid ) +{ + // clang-format off + Backend::par::apply( + image_1.size(), + [ + image_minus = image_1.data(), + image_plus = image_2.data(), + image_mid = image_mid.data(), + tangent = tangent.data() + ] SPIRIT_LAMBDA (int idx) + { + const Vector3 ex = { 1, 0, 0 }; + const Vector3 ey = { 0, 1, 0 }; + scalar epsilon = 1e-15; + + Vector3 axis = image_plus[idx].cross( image_minus[idx] ); + + // If the spins are anti-parallel, we choose an arbitrary axis + if( std::abs(image_minus[idx].dot(image_plus[idx]) + 1) < epsilon ) // Check if anti-parallel + { + if( std::abs( image_mid[idx].dot( ex ) - 1 ) > epsilon ) // Check if parallel to ex + axis = ex; + else + axis = ey; + } + tangent[idx] = image_mid[idx].cross( axis ); + } + ); + Manifoldmath::normalize(tangent); + // clang-format on +}; + /* Calculates the 'tangent' vectors, i.e.in crudest approximation the difference between an image and the neighbouring */ @@ -91,15 +128,13 @@ void Tangents( if( idx_img == 0 ) { auto & image_plus = *configurations[idx_img + 1]; - Vectormath::set_c_a( 1, image_plus, tangents[idx_img] ); - Vectormath::add_c_a( -1, image, tangents[idx_img] ); + Geodesic_Tangent( tangents[idx_img], image, image_plus, image ); // Use the accurate tangent at the endpoints, useful for the dimer method } // Last Image else if( idx_img == noi - 1 ) { auto & image_minus = *configurations[idx_img - 1]; - Vectormath::set_c_a( 1, image, tangents[idx_img] ); - Vectormath::add_c_a( -1, image_minus, tangents[idx_img] ); + Geodesic_Tangent( tangents[idx_img], image_minus, image, image ); // Use the accurate tangent at the endpoints, useful for the dimer method } // Images Inbetween else @@ -161,14 +196,12 @@ void Tangents( Vectormath::set_c_a( 1, t_plus, tangents[idx_img] ); Vectormath::add_c_a( 1, t_minus, tangents[idx_img] ); } - } - - // Project tangents into tangent planes of spin vectors to make them actual tangents - project_tangential( tangents[idx_img], image ); - - // Normalise in 3N - dimensional space - Manifoldmath::normalize( tangents[idx_img] ); + // Project tangents into tangent planes of spin vectors to make them actual tangents + project_tangential( tangents[idx_img], image ); + // Normalise in 3N - dimensional space + Manifoldmath::normalize( tangents[idx_img] ); + } } // end for idx_img } // end Tangents } // namespace Manifoldmath @@ -469,6 +502,36 @@ void spherical_to_cartesian_christoffel_symbols( const vectorfield & vf, MatrixX } } +void sparse_hessian_bordered_3N( + const vectorfield & image, const vectorfield & gradient, const SpMatrixX & hessian, SpMatrixX & hessian_out ) +{ + // Calculates a 3Nx3N matrix in the bordered Hessian approach and transforms it into the tangent basis, + // making the result a 2Nx2N matrix. The bordered Hessian's Lagrange multipliers assume a local extremum. + + int nos = image.size(); + VectorX lambda( nos ); + for( int i = 0; i < nos; ++i ) + lambda[i] = image[i].normalized().dot( gradient[i] ); + + // Construct hessian_out + typedef Eigen::Triplet T; + std::vector tripletList; + tripletList.reserve( hessian.nonZeros() + 3 * nos ); + + // Iterate over non zero entries of hesiian + for( int k = 0; k < hessian.outerSize(); ++k ) + { + for( SpMatrixX::InnerIterator it( hessian, k ); it; ++it ) + { + tripletList.push_back( T( it.row(), it.col(), it.value() ) ); + } + int j = k % 3; + int i = ( k - j ) / 3; + tripletList.push_back( T( k, k, -lambda[i] ) ); // Correction to the diagonal + } + hessian_out.setFromTriplets( tripletList.begin(), tripletList.end() ); +} + void hessian_bordered( const vectorfield & image, const vectorfield & gradient, const MatrixX & hessian, MatrixX & tangent_basis, MatrixX & hessian_out ) @@ -499,6 +562,7 @@ void hessian_bordered( hessian_out = tangent_basis.transpose() * tmp_3N * tangent_basis; } + void hessian_projected( const vectorfield & image, const vectorfield & gradient, const MatrixX & hessian, MatrixX & tangent_basis, MatrixX & hessian_out ) diff --git a/core/src/engine/Manifoldmath.cu b/core/src/engine/Manifoldmath.cu index 7a20cf7c0..3bdd1299a 100644 --- a/core/src/engine/Manifoldmath.cu +++ b/core/src/engine/Manifoldmath.cu @@ -2,6 +2,7 @@ #include #include +#include #include #include @@ -18,9 +19,12 @@ namespace Engine { void project_parallel(vectorfield & vf1, const vectorfield & vf2) { - vectorfield vf3 = vf1; - project_orthogonal(vf3, vf2); - Vectormath::add_c_a(-1, vf3, vf1); + scalar proj = Vectormath::dot(vf1, vf2); + Backend::par::apply( vf1.size(), [vf1 = vf1.data(), vf2 = vf2.data(), proj] SPIRIT_LAMBDA (int idx) + { + vf1[idx] = proj * vf2[idx]; + } + ); } __global__ void cu_project_orthogonal(Vector3 *vf1, const Vector3 *vf2, scalar proj, size_t N) diff --git a/core/src/engine/Method.cpp b/core/src/engine/Method.cpp index 0b9a5a67e..3f70440f2 100644 --- a/core/src/engine/Method.cpp +++ b/core/src/engine/Method.cpp @@ -24,8 +24,9 @@ Method::Method( std::shared_ptr parameters, int idx_img // Sender name for log messages this->SenderName = Log_Sender::All; - // Default history contains max_torque - this->history = std::map>{ { "max_torque", { this->max_torque } } }; + this->history_iteration = std::vector(); + this->history_max_torque = std::vector(); + this->history_energy = std::vector(); // TODO: is this a good idea? this->n_iterations = std::max( long( 1 ), this->parameters->n_iterations ); diff --git a/core/src/engine/Method_GNEB.cpp b/core/src/engine/Method_GNEB.cpp index 28cb05d0f..5d410bc09 100644 --- a/core/src/engine/Method_GNEB.cpp +++ b/core/src/engine/Method_GNEB.cpp @@ -1,4 +1,6 @@ #include +#include +#include #include #include #include @@ -6,13 +8,11 @@ #include #include #include +#include #include #include #include -#include -#include - #include using namespace Utility; @@ -40,10 +40,15 @@ Method_GNEB::Method_GNEB( std::shared_ptr chain this->F_gradient = std::vector( this->noi, vectorfield( this->nos, { 0, 0, 0 } ) ); // [noi][nos] this->F_spring = std::vector( this->noi, vectorfield( this->nos, { 0, 0, 0 } ) ); // [noi][nos] this->f_shrink = vectorfield( this->nos, { 0, 0, 0 } ); // [nos] - this->xi = vectorfield( this->nos, { 0, 0, 0 } ); // [nos] + this->xi = vectorfield( this->nos, { 0, 0, 0 } ); + + this->F_translation_left = vectorfield( this->nos, { 0, 0, 0 } ); + this->F_translation_right = vectorfield( this->nos, { 0, 0, 0 } ); // Tangents this->tangents = std::vector( this->noi, vectorfield( this->nos, { 0, 0, 0 } ) ); // [noi][nos] + this->tangent_endpoints_left = vectorfield( this->nos, { 0, 0, 0 } ); // [nos] + this->tangent_endpoints_right = vectorfield( this->nos, { 0, 0, 0 } ); // [nos] // We assume that the chain is not converged before the first iteration this->max_torque = this->chain->gneb_parameters->force_convergence + 1.0; @@ -55,7 +60,7 @@ Method_GNEB::Method_GNEB( std::shared_ptr chain this->configurations[i] = this->systems[i]->spins; // History - this->history = std::map>{ { "max_torque", { this->max_torque } } }; + // this->history = std::map>{ { "max_torque", { this->max_torque } } }; //---- Initialise Solver-specific variables this->Initialize(); @@ -100,13 +105,10 @@ void Method_GNEB::Calculate_Force( // while the gradient force is manipulated (e.g. projected) auto eff_field = this->chain->images[img]->effective_field.data(); auto f_grad = F_gradient[img].data(); - Backend::par::apply( - image.size(), - [eff_field, f_grad] SPIRIT_LAMBDA( int idx ) - { - eff_field[idx] *= -1; - f_grad[idx] = eff_field[idx]; - } ); + Backend::par::apply( image.size(), [eff_field, f_grad] SPIRIT_LAMBDA( int idx ) { + eff_field[idx] *= -1; + f_grad[idx] = eff_field[idx]; + } ); if( img > 0 ) { @@ -145,11 +147,11 @@ void Method_GNEB::Calculate_Force( scalar range_Rx = interp[0].at( std::distance( interp[0].begin(), std::max_element( interp[0].begin(), interp[0].end() ) ) ) - interp[0].at( - std::distance( interp[0].begin(), std::min_element( interp[0].begin(), interp[0].end() ) ) ); + std::distance( interp[0].begin(), std::min_element( interp[0].begin(), interp[0].end() ) ) ); scalar range_E = interp[1].at( std::distance( interp[1].begin(), std::max_element( interp[1].begin(), interp[1].end() ) ) ) - interp[1].at( - std::distance( interp[1].begin(), std::min_element( interp[1].begin(), interp[1].end() ) ) ); + std::distance( interp[1].begin(), std::min_element( interp[1].begin(), interp[1].end() ) ) ); for( int idx_image = 1; idx_image < this->chain->noi; ++idx_image ) { @@ -251,6 +253,104 @@ void Method_GNEB::Calculate_Force( // Copy out Vectormath::set_c_a( 1, F_total[img], forces[img] ); } // end for img=1..noi-1 + + // Moving endpoints + 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] ); + + // Overall translational force + if( chain->gneb_parameters->translating_endpoints ) + { + // clang-format off + Backend::par::apply(nos, + [ + F_translation_left = F_translation_left.data(), + F_translation_right = F_translation_right.data(), + F_gradient_left = F_gradient[0].data(), + F_gradient_right = F_gradient[noi-1].data(), + spins_left = this->chain->images[0]->spins->data(), + spins_right = this->chain->images[noi-1]->spins->data() + ] SPIRIT_LAMBDA ( int idx) + { + const Vector3 axis = spins_left[idx].cross(spins_right[idx]); + const scalar angle = acos(spins_left[idx].dot(spins_right[idx])); + + // Rotation matrix that rotates spin_left to spin_right + Matrix3 rotation_matrix = Eigen::AngleAxis(angle, axis.normalized()).toRotationMatrix(); + + if ( abs(spins_left[idx].dot(spins_right[idx])) >= 1.0 ) // Angle can become nan for collinear spins + rotation_matrix = Matrix3::Identity(); + + const Vector3 F_gradient_right_rotated = rotation_matrix * F_gradient_right[idx]; + F_translation_left[idx] = -0.5 * (F_gradient_left[idx] + F_gradient_right_rotated); + + const Vector3 F_gradient_left_rotated = rotation_matrix.transpose() * F_gradient_left[idx]; + F_translation_right[idx] = -0.5 * (F_gradient_left_rotated + F_gradient_right[idx]); + } + ); + // clang-format on + + Manifoldmath::project_parallel( F_translation_left, tangents[0] ); + Manifoldmath::project_parallel( F_translation_right, tangents[chain->noi - 1] ); + } + + scalar rotational_coeff = 1.0; + if(chain->gneb_parameters->escape_first) + { + // Estimate the curvature along the tangent and only activate the rotational force, if it is negative + scalar proj_left = Vectormath::dot( F_gradient[0], tangents[0] ); + scalar proj_right = Vectormath::dot( F_gradient[chain->noi-1], tangents[chain->noi-1] ); + if (proj_left > proj_right) + { + rotational_coeff = 0.0; + } + } + + for( int img : { 0, chain->noi - 1 } ) + { + scalar delta_Rx0 = ( img == 0 ) ? chain->gneb_parameters->equilibrium_delta_Rx_left : + chain->gneb_parameters->equilibrium_delta_Rx_right; + scalar delta_Rx = ( img == 0 ) ? Rx[1] - Rx[0] : Rx[chain->noi - 1] - Rx[chain->noi - 2]; + + auto spring_constant = ( ( img == 0 ) ? 1.0 : -1.0 ) * this->chain->gneb_parameters->spring_constant; + auto projection = Vectormath::dot( F_gradient[img], tangents[img] ); + + auto F_translation = ( img == 0 ) ? F_translation_left.data() : F_translation_right.data(); + + // std::cout << " === " << img << " ===\n"; + // fmt::print( "tangent_coeff = {}\n", spring_constant * (delta_Rx - delta_Rx0) ); + // fmt::print( "delta_Rx {}\n", delta_Rx ); + // fmt::print( "delta_Rx0 {}\n", delta_Rx0 ); + + // clang-format off + Backend::par::apply( + nos, + [ + F_total = F_total[img].data(), + F_gradient = F_gradient[img].data(), + forces = forces[img].data(), + tangents = tangents[img].data(), + tangent_coeff = spring_constant * (delta_Rx - delta_Rx0), + F_translation, + projection, + rotational_coeff + ] SPIRIT_LAMBDA ( int idx ) + { + forces[idx] = rotational_coeff * (F_gradient[idx] - projection * tangents[idx]) + + tangent_coeff * tangents[idx] + + F_translation[idx]; + + // std::cout << F_translation[idx].transpose() << "\n"; + F_total[idx] = forces[idx]; + } + ); + // clang-format on + } + } + } // end Calculate template @@ -261,8 +361,14 @@ void Method_GNEB::Calculate_Force_Virtual( using namespace Utility; // Calculate the cross product with the spin configuration to get direct minimization - for( std::size_t i = 1; i < configurations.size() - 1; ++i ) + for( std::size_t i = 0; i < configurations.size(); ++i ) { + + if( !chain->gneb_parameters->moving_endpoints && ( i == 0 || i == configurations.size() - 1 ) ) + { + continue; + } + auto & image = *configurations[i]; auto & force = forces[i]; auto & force_virtual = forces_virtual[i]; @@ -305,7 +411,7 @@ void Method_GNEB::Hook_Post_Iteration() this->max_torque = 0; std::fill( this->max_torque_all.begin(), this->max_torque_all.end(), 0 ); - for( int img = 1; img < chain->noi - 1; ++img ) + for( int img = 0; img < chain->noi; ++img ) { scalar fmax = this->MaxTorque_on_Image( *( this->systems[img]->spins ), F_total[img] ); // Set maximum per image @@ -338,7 +444,7 @@ void Method_GNEB::Hook_Post_Iteration() // Rx chain->Rx = this->Rx; // E - for( int img = 1; img < chain->noi; ++img ) + for( int img = 0; img < chain->noi; ++img ) chain->images[img]->E = this->energies[img]; // Rx interpolated chain->Rx_interpolated = interp[0]; @@ -491,7 +597,9 @@ template void Method_GNEB::Save_Current( std::string starttime, int iteration, bool initial, bool final ) { // History save - this->history["max_torque"].push_back( this->max_torque ); + this->history_iteration.push_back( this->iteration ); + this->history_max_torque.push_back( this->max_torque ); + this->history_energy.push_back( this->systems[0]->E ); // File save if( this->parameters->output_any ) @@ -514,8 +622,8 @@ void Method_GNEB::Save_Current( std::string starttime, int iteration, bo preEnergiesFile = this->parameters->output_folder + "/" + fileTag + "Chain_Energies"; // Function to write or append image and energy files - auto writeOutputChain = [this, preChainFile, preEnergiesFile, iteration]( const std::string& suffix, bool append ) - { + auto writeOutputChain = [this, preChainFile, preEnergiesFile, + iteration]( const std::string & suffix, bool append ) { try { // File name @@ -559,8 +667,7 @@ void Method_GNEB::Save_Current( std::string starttime, int iteration, bo }; Calculate_Interpolated_Energy_Contributions(); - auto writeOutputEnergies = [this, preChainFile, preEnergiesFile, iteration]( const std::string & suffix ) - { + auto writeOutputEnergies = [this, preChainFile, preEnergiesFile, iteration]( const std::string & suffix ) { bool normalize = this->chain->gneb_parameters->output_energies_divide_by_nspins; bool readability = this->chain->gneb_parameters->output_energies_add_readability_lines; diff --git a/core/src/engine/Method_LLG.cpp b/core/src/engine/Method_LLG.cpp index a3a459ebb..39e260550 100644 --- a/core/src/engine/Method_LLG.cpp +++ b/core/src/engine/Method_LLG.cpp @@ -41,11 +41,6 @@ Method_LLG::Method_LLG( std::shared_ptr system, int i this->force_converged = std::vector( this->noi, false ); this->max_torque = system->llg_parameters->force_convergence + 1.0; - // History - this->history = std::map>{ { "max_torque", { this->max_torque } }, - { "E", { this->max_torque } }, - { "M_z", { this->max_torque } } }; - // Create shared pointers to the method's systems' spin configurations this->configurations = std::vector>( this->noi ); for( int i = 0; i < this->noi; ++i ) @@ -304,11 +299,16 @@ template void Method_LLG::Save_Current( std::string starttime, int iteration, bool initial, bool final ) { // History save - this->history["max_torque"].push_back( this->max_torque ); - this->systems[0]->UpdateEnergy(); - this->history["E"].push_back( this->systems[0]->E ); - auto mag = Engine::Vectormath::Magnetization( *this->systems[0]->spins ); - this->history["M_z"].push_back( mag[2] ); + this->history_iteration.push_back( this->iteration ); + this->history_max_torque.push_back( this->max_torque ); + this->history_energy.push_back( this->systems[0]->E ); + + // this->history["max_torque"].push_back( this->max_torque ); + // this->systems[0]->UpdateEnergy(); + // this->history["E"].push_back( this->systems[0]->E ); + // Removed magnetization, since at the moment it required a temporary allocation to compute + // auto mag = Engine::Vectormath::Magnetization( *this->systems[0]->spins ); + // this->history["M_z"].push_back( mag[2] ); // File save if( this->parameters->output_any ) @@ -334,42 +334,40 @@ void Method_LLG::Save_Current( std::string starttime, int iteration, boo // Function to write or append image and energy files auto writeOutputConfiguration - = [this, preSpinsFile, preEnergyFile, iteration]( const std::string & suffix, bool append ) - { - try - { - // File name and comment - std::string spinsFile = preSpinsFile + suffix + ".ovf"; - std::string output_comment = fmt::format( - "{} simulation ({} solver)\n# Desc: Iteration: {}\n# Desc: Maximum torque: {}", - this->Name(), this->SolverFullName(), iteration, this->max_torque ); - - // File format - IO::VF_FileFormat format = this->systems[0]->llg_parameters->output_vf_filetype; - - // Spin Configuration - auto & spins = *this->systems[0]->spins; - auto segment = IO::OVF_Segment( *this->systems[0] ); - 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 = 3; - segment.valuelabels = strdup( "spin_x spin_y spin_z" ); - segment.valueunits = strdup( "none none none" ); - if( append ) - IO::OVF_File( spinsFile ).append_segment( segment, spins[0].data(), int( format ) ); - else - IO::OVF_File( spinsFile ).write_segment( segment, spins[0].data(), int( format ) ); - } - catch( ... ) - { - spirit_handle_exception_core( "LLG output failed" ); - } - }; - - auto writeOutputEnergy - = [this, preSpinsFile, preEnergyFile, iteration]( const std::string & suffix, bool append ) - { + = [this, preSpinsFile, preEnergyFile, iteration]( const std::string & suffix, bool append ) { + try + { + // File name and comment + std::string spinsFile = preSpinsFile + suffix + ".ovf"; + std::string output_comment = fmt::format( + "{} simulation ({} solver)\n# Desc: Iteration: {}\n# Desc: Maximum torque: {}", + this->Name(), this->SolverFullName(), iteration, this->max_torque ); + + // File format + IO::VF_FileFormat format = this->systems[0]->llg_parameters->output_vf_filetype; + + // Spin Configuration + auto & spins = *this->systems[0]->spins; + auto segment = IO::OVF_Segment( *this->systems[0] ); + 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 = 3; + segment.valuelabels = strdup( "spin_x spin_y spin_z" ); + segment.valueunits = strdup( "none none none" ); + if( append ) + IO::OVF_File( spinsFile ).append_segment( segment, spins[0].data(), int( format ) ); + else + IO::OVF_File( spinsFile ).write_segment( segment, spins[0].data(), int( format ) ); + } + catch( ... ) + { + spirit_handle_exception_core( "LLG output failed" ); + } + }; + + auto writeOutputEnergy = [this, preSpinsFile, preEnergyFile, + iteration]( const std::string & suffix, bool append ) { bool normalize = this->systems[0]->llg_parameters->output_energy_divide_by_nspins; bool readability = this->systems[0]->llg_parameters->output_energy_add_readability_lines; diff --git a/core/src/engine/Method_MC.cpp b/core/src/engine/Method_MC.cpp index 923e7037f..6ab8b836a 100644 --- a/core/src/engine/Method_MC.cpp +++ b/core/src/engine/Method_MC.cpp @@ -35,9 +35,9 @@ Method_MC::Method_MC( std::shared_ptr system, int idx_img, in // this->max_torque = system->mc_parameters->force_convergence + 1.0; // History - this->history = std::map>{ { "max_torque", { this->max_torque } }, - { "E", { this->max_torque } }, - { "M_z", { this->max_torque } } }; + // this->history = std::map>{ { "max_torque", { this->max_torque } }, + // { "E", { this->max_torque } }, + // { "M_z", { this->max_torque } } }; this->parameters_mc = system->mc_parameters; diff --git a/core/src/engine/Method_MMF.cpp b/core/src/engine/Method_MMF.cpp index 2005954fa..36c5e4007 100644 --- a/core/src/engine/Method_MMF.cpp +++ b/core/src/engine/Method_MMF.cpp @@ -36,7 +36,7 @@ Method_MMF::Method_MMF( std::shared_ptr system, int i this->SenderName = Utility::Log_Sender::MMF; // History - this->history = std::map>{ { "max_torque", { this->max_torque } } }; + // this->history = std::map>{ { "max_torque", { this->max_torque } } }; // We assume that the systems are not converged before the first iteration this->max_torque = system->mmf_parameters->force_convergence + 1.0; @@ -427,7 +427,8 @@ template void Method_MMF::Save_Current( std::string starttime, int iteration, bool initial, bool final ) { // History save - this->history["max_torque"].push_back( this->max_torque ); + this->history_iteration.push_back( this->iteration ); + this->history_max_torque.push_back( this->max_torque ); // File save if( this->parameters->output_any ) diff --git a/core/src/engine/Sparse_HTST.cpp b/core/src/engine/Sparse_HTST.cpp index 640fca12e..f0a6cf579 100644 --- a/core/src/engine/Sparse_HTST.cpp +++ b/core/src/engine/Sparse_HTST.cpp @@ -26,37 +26,46 @@ namespace Engine namespace Sparse_HTST { -void Sparse_Get_Lowest_Eigenvector( const SpMatrixX & matrix, int nos, scalar & lowest_evalue, VectorX & lowest_evec ) +// TODO: fix the mess of the partial spectrum solvers, should in general be refactored and merged with ema method somehow +enum Partial_Spectrum_Solver { - Log( Utility::Log_Level::All, Utility::Log_Sender::HTST, " Using Spectra to compute lowest eigenmode..." ); + RAYLEIGH, // The "old" way + LANCZOS // Lanczos with sensible parameters +}; +const Partial_Spectrum_Solver spectrum_solver = Partial_Spectrum_Solver::LANCZOS; - VectorX evalues; - MatrixX evectors; +void Sparse_Get_Lowest_Eigenvectors( const SpMatrixX & matrix, scalar max_evalue, scalarfield & evalues, std::vector & evecs ) +{ + Log( Utility::Log_Level::All, Utility::Log_Sender::HTST, " Using Spectra to compute lowest eigenmodes..." ); + + int nos = matrix.rows() / 2; + int n_modes = 6; // Number of lowest modes to be computed (should always be enough) - int n_steps = std::max( 2, nos ); + int ncv = std::min(2*nos, std::max(2*n_modes + 1, 20)); // This is the default value used by scipy.sparse + int max_iter = 20*nos; - // Create a Spectra solver Spectra::SparseSymMatProd op( matrix ); + // Create and initialize a Spectra solver Spectra::SymEigsSolver> matrix_spectrum( - &op, 1, n_steps ); - + &op, n_modes, ncv); matrix_spectrum.init(); - int nconv = matrix_spectrum.compute(); - if( matrix_spectrum.info() == Spectra::SUCCESSFUL ) - { - evalues = matrix_spectrum.eigenvalues().real(); - evectors = matrix_spectrum.eigenvectors().real(); - } - else + int nconv = matrix_spectrum.compute( max_iter, 1e-10, int( Spectra::SMALLEST_ALGE ) ); + + if( matrix_spectrum.info() != Spectra::SUCCESSFUL ) { Log( Utility::Log_Level::All, Utility::Log_Sender::HTST, - " Failed to calculate lowest eigenmode. Aborting!" ); + " Failed to calculate lowest eigenmode. Aborting!" ); return; } - lowest_evalue = evalues[0]; - lowest_evec = evectors.col( 0 ); + for( int i=0; i evecs_sp = std::vector( 0 ); - Sparse_Get_Lowest_Eigenvectors_VP( sparse_hessian_sp_geodesic_2N, epsilon, evalues_sp, evecs_sp ); + + if ( spectrum_solver == Partial_Spectrum_Solver::RAYLEIGH ) + Sparse_Get_Lowest_Eigenvectors_VP( sparse_hessian_sp_geodesic_2N, epsilon, evalues_sp, evecs_sp ); + else if ( spectrum_solver == Partial_Spectrum_Solver::LANCZOS ) + Sparse_Get_Lowest_Eigenvectors( sparse_hessian_sp_geodesic_2N, epsilon, evalues_sp, evecs_sp ); + scalar lowest_evalue = evalues_sp[0]; VectorX & lowest_evector = evecs_sp[0]; @@ -432,7 +446,11 @@ void Calculate( Data::HTST_Info & htst_info ) // Calculate modes at minimum (needed for zero-mode volume) std::vector evecs_min = std::vector( 0 ); - Sparse_Get_Lowest_Eigenvectors_VP( sparse_hessian_geodesic_min_2N, epsilon, evalues_min, evecs_min ); + + if ( spectrum_solver == Partial_Spectrum_Solver::RAYLEIGH ) + Sparse_Get_Lowest_Eigenvectors_VP( sparse_hessian_geodesic_min_2N, epsilon, evalues_min, evecs_min ); + else if ( spectrum_solver == Partial_Spectrum_Solver::LANCZOS ) + Sparse_Get_Lowest_Eigenvectors( sparse_hessian_geodesic_min_2N, epsilon, evalues_min, evecs_min ); // Checking for zero modes at the minimum.. Log( Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Checking for zero modes at the minimum ..." ); diff --git a/core/src/engine/Vectormath.cpp b/core/src/engine/Vectormath.cpp index f81094ad4..623b2ee2a 100644 --- a/core/src/engine/Vectormath.cpp +++ b/core/src/engine/Vectormath.cpp @@ -492,16 +492,21 @@ Vector3 decompose( const Vector3 & v, const std::vector & basis ) ///////////////////////////////////////////////////////////////// -std::array Magnetization( const vectorfield & vf ) +std::array Magnetization( const vectorfield & vf, const scalarfield & mu_s ) { - Vector3 vfmean = mean( vf ); + vectorfield temp_vf = vf; + scale( temp_vf, mu_s ); + Vector3 vfmean = mean( temp_vf ); std::array M{ vfmean[0], vfmean[1], vfmean[2] }; return M; } -scalar -TopologicalCharge( const vectorfield & vf, const Data::Geometry & geometry, const intfield & boundary_conditions ) +void TopologicalChargeDensity( + const vectorfield & vf, const Data::Geometry & geometry, const intfield & boundary_conditions, + scalarfield & charge_density, std::vector & triangle_indices ) { + charge_density.resize( 0 ); + // This implementations assumes // 1. No basis atom lies outside the cell spanned by the basis vectors of the lattice // 2. The geometry is a plane in x and y and spanned by the first 2 basis_vectors of the lattice @@ -567,6 +572,7 @@ TopologicalCharge( const vectorfield & vf, const Data::Geometry & geometry, cons for( int a = 0; a < geometry.n_cells[0]; ++a ) { std::array tri_spins; + std::array tri_indices; // bools to check wether it is allowed to take the next lattice site in direction a, b or a+b bool a_next_allowed = ( a + 1 < geometry.n_cells[0] || boundary_conditions[0] ); bool b_next_allowed = ( b + 1 < geometry.n_cells[1] || boundary_conditions[1] ); @@ -600,13 +606,28 @@ TopologicalCharge( const vectorfield & vf, const Data::Geometry & geometry, cons break; } tri_spins[i] = vf[idx]; + tri_indices[i] = idx; } if( valid_triangle ) - charge += sign * solid_angle_2( tri_spins[0], tri_spins[1], tri_spins[2] ); + { + triangle_indices.push_back( tri_indices[0] ); + triangle_indices.push_back( tri_indices[1] ); + triangle_indices.push_back( tri_indices[2] ); + charge_density.push_back( + sign / ( 4.0 * Pi ) * solid_angle_2( tri_spins[0], tri_spins[1], tri_spins[2] ) ); + } } } } - return charge / ( 4 * Pi ); +} + +// Calculate the topological charge inside a vectorfield +scalar TopologicalCharge( const vectorfield & vf, const Data::Geometry & geom, const intfield & boundary_conditions ) +{ + scalarfield charge_density( 0 ); + std::vector triangle_indices( 0 ); + TopologicalChargeDensity( vf, geom, boundary_conditions, charge_density, triangle_indices ); + return sum( charge_density ); } void get_gradient_distribution( diff --git a/core/src/io/Configparser.cpp b/core/src/io/Configparser.cpp index 9a840c8fc..5574130de 100644 --- a/core/src/io/Configparser.cpp +++ b/core/src/io/Configparser.cpp @@ -530,7 +530,7 @@ try parameter_log.emplace_back( fmt::format( " pinned site[{}]: {} at ({} {} {}) = ({})", i, pinning.sites[i].i, pinning.sites[i].translations[0], pinning.sites[i].translations[1], - pinning.sites[i].translations[2], pinning.spins[0].transpose() ) ); + pinning.sites[i].translations[2], pinning.spins[i].transpose() ) ); } } } @@ -1056,6 +1056,10 @@ std::unique_ptr Parameters_Method_GNEB_from_Config config_file_handle.Read_Single( parameters->n_iterations_log, "gneb_n_iterations_log" ); config_file_handle.Read_Single( parameters->n_iterations_amortize, "gneb_n_iterations_amortize" ); config_file_handle.Read_Single( parameters->n_E_interpolations, "gneb_n_energy_interpolations" ); + config_file_handle.Read_Single( parameters->moving_endpoints, "gneb_moving_endpoints" ); + config_file_handle.Read_Single( parameters->equilibrium_delta_Rx_left, "gneb_equilibrium_delta_Rx_left" ); + config_file_handle.Read_Single( parameters->equilibrium_delta_Rx_right, "gneb_equilibrium_delta_Rx_right" ); + config_file_handle.Read_Single( parameters->translating_endpoints, "gneb_translating_endpoints" ); } catch( ... ) { @@ -1077,8 +1081,11 @@ std::unique_ptr Parameters_Method_GNEB_from_Config parameter_log.emplace_back( fmt::format( " {:<18} = {}", "maximum walltime", str_max_walltime ) ); parameter_log.emplace_back( fmt::format( " {:<18} = {}", "n_iterations", parameters->n_iterations ) ); parameter_log.emplace_back( fmt::format( " {:<18} = {}", "n_iterations_log", parameters->n_iterations_log ) ); - parameter_log.emplace_back( - fmt::format( " {:<18} = {}", "n_iterations_amortize", parameters->n_iterations_amortize ) ); + parameter_log.emplace_back( fmt::format( " {:<18} = {}", "n_iterations_amortize", parameters->n_iterations_amortize ) ); + parameter_log.emplace_back( fmt::format( " {:<18} = {}", "moving_endpoints", parameters->moving_endpoints ) ); + parameter_log.emplace_back( fmt::format( " {:<18} = {}", "equilibrium_delta_Rx_left", parameters->equilibrium_delta_Rx_left ) ); + parameter_log.emplace_back( fmt::format( " {:<18} = {}", "equilibrium_delta_Rx_right", parameters->equilibrium_delta_Rx_right ) ); + parameter_log.emplace_back( fmt::format( " {:<18} = {}", "translating_endpoints", parameters->translating_endpoints ) ); parameter_log.emplace_back( fmt::format( " {:<18} = \"{}\"", "output_folder", parameters->output_folder ) ); parameter_log.emplace_back( fmt::format( " {:<18} = {}", "output_any", parameters->output_any ) ); parameter_log.emplace_back( fmt::format( " {:<18} = {}", "output_initial", parameters->output_initial ) ); diff --git a/core/src/utility/Exception.cpp b/core/src/utility/Exception.cpp index f29918ac7..576c0c2da 100644 --- a/core/src/utility/Exception.cpp +++ b/core/src/utility/Exception.cpp @@ -6,7 +6,8 @@ namespace Utility { -void rethrow( const std::string & message, const char * file, const unsigned int line, const char * function ) +void rethrow( const std::string & message, const char * file, const unsigned int line, const char * function ) noexcept( + false ) try { std::rethrow_exception( std::current_exception() ); @@ -27,7 +28,7 @@ catch( ... ) std::throw_with_nested( ex ); } -void Backtrace_Exception() +void Backtrace_Exception() noexcept( false ) try { if( std::exception_ptr eptr = std::current_exception() ) @@ -38,7 +39,7 @@ try { Log( Log_Level::Severe, Log_Sender::API, "Unknown Exception. Terminating" ); Log.Append_to_File(); - std::exit( EXIT_FAILURE ); // Exit the application. May lead to data loss + std::exit( EXIT_FAILURE ); // Exit the application. May lead to data loss! } } catch( const Exception & ex ) @@ -67,7 +68,8 @@ catch( const std::exception & ex ) } void Handle_Exception_API( - const char * file, const unsigned int line, const char * api_function, const int idx_image, const int idx_chain ) + const char * file, const unsigned int line, const char * function, const int idx_image, + const int idx_chain ) noexcept try { // Rethrow in order to get an exception reference (instead of pointer) @@ -85,8 +87,8 @@ try str_exception = "SEVERE exception"; Log( ex.level, Log_Sender::API, fmt::format( - "Caught {} in API function \'{}\' (at {}:{})\n{}Exception backtrace:", str_exception, api_function, - file, line, Log.tags_space ), + "Caught {} in API function \'{}\' (at {}:{})\n{}Exception backtrace:", str_exception, function, file, + line, Log.tags_space ), idx_image, idx_chain ); // Create backtrace @@ -101,7 +103,7 @@ try { Log( Log_Level::Severe, Log_Sender::API, "TERMINATING!", idx_image, idx_chain ); Log.Append_to_File(); - std::exit( EXIT_FAILURE ); // exit the application. May lead to data loss + std::exit( EXIT_FAILURE ); // Exit the application. May lead to data loss! } // If it was recoverable we now write to Log @@ -113,8 +115,8 @@ try idx_chain ); Log( Log_Level::Severe, Log_Sender::API, fmt::format( - "Caught std::exception in API function \'{}\' (at {}:{})\n{:>49}Exception backtrace:", api_function, - file, line, " " ), + "Caught std::exception in API function \'{}\' (at {}:{})\n{:>49}Exception backtrace:", function, file, + line, " " ), idx_image, idx_chain ); // Create backtrace Backtrace_Exception(); @@ -124,7 +126,7 @@ try // Terminate Log( Log_Level::Severe, Log_Sender::API, "TERMINATING!", idx_image, idx_chain ); Log.Append_to_File(); - std::exit( EXIT_FAILURE ); // Exit the application. May lead to data loss + std::exit( EXIT_FAILURE ); // Exit the application. May lead to data loss! } catch( ... ) { @@ -134,25 +136,25 @@ try fmt::format( "Caught unknown exception in API function \'{}\' (at {}:{})\n{:>49}Cannot backtrace unknown " "exceptions...", - api_function, file, line, " " ), + function, file, line, " " ), idx_image, idx_chain ); Log( Log_Level::Severe, Log_Sender::API, "-----------------------------------------------------", idx_image, idx_chain ); Log( Log_Level::Severe, Log_Sender::API, "Unable to handle! TERMINATING!", idx_image, idx_chain ); Log.Append_to_File(); - std::exit( EXIT_FAILURE ); // Exit the application. May lead to data loss + std::exit( EXIT_FAILURE ); // Exit the application. May lead to data loss! } } catch( ... ) { - std::cerr << "Another exception occurred while handling an exception from \'" << api_function << "\' (at " << file + std::cerr << "Another exception occurred while handling an exception from \'" << function << "\' (at " << file << ":" << line << ")!" << std::endl; std::cerr << "TERMINATING!" << std::endl; - std::exit( EXIT_FAILURE ); // Exit the application. May lead to data loss + std::exit( EXIT_FAILURE ); // Exit the application. May lead to data loss! } void Handle_Exception_Core( - const std::string & message, const char * file, const unsigned int line, const char * function ) + const std::string & message, const char * file, const unsigned int line, const char * function ) noexcept( false ) try { // Rethrow in order to get an exception reference (instead of pointer) @@ -196,47 +198,4 @@ catch( ... ) spirit_rethrow( message ); } -// void Spirit_Exception( const Exception & ex, int idx_image, int idx_chain ) -// { -// switch ( ex ) -// { -// case Exception::File_not_Found : -// Log( Log_Level::Warning, Log_Sender::API, "File not found. Unable to open.", -// idx_image, idx_chain ); -// break; - -// case Exception::System_not_Initialized : -// Log( Log_Level::Severe, Log_Sender::API, "System is uninitialized. Terminating.", idx_image, idx_chain ); -// Log.Append_to_File(); -// std::exit( EXIT_FAILURE ); // exit the application. may lead to data loss -// break; - -// case Exception::Division_by_zero: -// Log( Log_Level::Severe, Log_Sender::API, "Division by zero", idx_image, idx_chain ); -// Log.Append_to_File(); -// std::exit( EXIT_FAILURE ); // exit the application. may lead to data loss -// break; - -// case Exception::Simulated_domain_too_small: -// Log( Log_Level::Error, Log_Sender::API, std::string( "Simulated domain is " ) + -// std::string( "too small. No action taken." ), idx_image, idx_chain ); -// break; - -// case Exception::Not_Implemented: -// Log( Log_Level::Warning, Log_Sender::API, std::string( "Feature not " ) + -// std::string( " implemented. No action taken." ), idx_image, idx_chain ); -// break; - -// case Exception::Non_existing_Image: -// Log( Log_Level::Warning, Log_Sender::API, "Non existing image. No action taken.", -// idx_image, idx_chain ); -// break; - -// case Exception::Non_existing_Chain: -// Log( Log_Level::Warning, Log_Sender::API, "Non existing chain. No action taken.", -// idx_image, idx_chain ); -// break; -// } -// } - } // namespace Utility \ No newline at end of file diff --git a/core/src/utility/Timing.cpp b/core/src/utility/Timing.cpp index fc873fe8f..3f384bb9f 100644 --- a/core/src/utility/Timing.cpp +++ b/core/src/utility/Timing.cpp @@ -77,29 +77,37 @@ scalar HoursPassed( duration dt ) duration DurationFromString( const std::string & dt ) { - std::uint8_t hours = 0, minutes = 0; - std::uint64_t seconds = 0; + std::int32_t hours = 0, minutes = 0; + std::int64_t seconds = 0; std::istringstream iss( dt ); std::string token = ""; // Hours if( std::getline( iss, token, ':' ) ) + { if( !token.empty() ) - iss >> hours; + hours = std::stoi( token ); + } + // Minutes if( std::getline( iss, token, ':' ) ) + { if( !token.empty() ) - iss >> minutes; + minutes = std::stoi( token ); + } + // Seconds if( std::getline( iss, token, ':' ) ) + { if( !token.empty() ) - iss >> seconds; + seconds = std::stol( token ); + } // Convert to std::chrono::seconds seconds += 60 * minutes + 60 * 60 * hours; - std::chrono::seconds chrono_seconds( seconds ); + std::chrono::seconds chrono_seconds( seconds ); // Return duration return duration( chrono_seconds ); } diff --git a/core/test/input/api.cfg b/core/test/input/api.cfg index 106f21f25..8c03cf0f5 100644 --- a/core/test/input/api.cfg +++ b/core/test/input/api.cfg @@ -39,9 +39,13 @@ log_file_level 3 ################### Geometry ##################### ### The bravais lattice type -bravais_lattice sc +bravais_lattice hex2d +basis +2 +0 0 0 +0.333 0.333 0.0 ### Number of basis cells along principal ### directions (a b c) -n_basis_cells 100 100 1 +n_basis_cells 50 50 1 ################# End Geometry ################### \ No newline at end of file diff --git a/core/test/test_solvers.cpp b/core/test/test_solvers.cpp index 06a1947b6..c2d26af87 100644 --- a/core/test/test_solvers.cpp +++ b/core/test/test_solvers.cpp @@ -37,7 +37,7 @@ TEST_CASE( "Solvers testing", "[solvers]" ) // Expected values float energy_expected = -5849.69140625f; - std::vector magnetization_expected{ 0, 0, 0.79977f }; + std::vector magnetization_expected{ 0, 0, 2.0 * 0.79977f }; // Result values scalar energy; @@ -78,7 +78,7 @@ TEST_CASE( "Solvers testing", "[solvers]" ) // Expected values float energy_sp_expected = -5811.5244140625f; - std::vector magnetization_sp_expected{ 0, 0, 0.96657f }; + std::vector magnetization_sp_expected{ 0, 0, 2.0 * 0.96657f }; // Result values scalar energy_sp; diff --git a/core/test/test_vectormath.cpp b/core/test/test_vectormath.cpp index 2c190d492..b7eca4d59 100644 --- a/core/test/test_vectormath.cpp +++ b/core/test/test_vectormath.cpp @@ -7,16 +7,18 @@ TEST_CASE( "Vectormath operations", "[vectormath]" ) int N = 10000; int N_check = std::min( 100, N ); scalarfield sf( N, 1 ); + scalar mu_s = 1.34; + scalarfield mu_s_field( N, mu_s ); vectorfield vf1( N, Vector3{ 1.0, 1.0, 1.0 } ); vectorfield vf2( N, Vector3{ -1.0, 1.0, 1.0 } ); SECTION( "Magnetization" ) { Engine::Vectormath::fill( vf1, { 0, 0, 1 } ); - auto m = Engine::Vectormath::Magnetization( vf1 ); + auto m = Engine::Vectormath::Magnetization( vf1, mu_s_field ); REQUIRE( m[0] == Approx( 0 ) ); REQUIRE( m[1] == Approx( 0 ) ); - REQUIRE( m[2] == Approx( 1 ) ); + REQUIRE( m[2] == Approx( mu_s ) ); } SECTION( "Rotate" ) diff --git a/ui-cpp/ui-imgui/include/parameters_widget.hpp b/ui-cpp/ui-imgui/include/parameters_widget.hpp index eb9f295a7..ea7bcc578 100644 --- a/ui-cpp/ui-imgui/include/parameters_widget.hpp +++ b/ui-cpp/ui-imgui/include/parameters_widget.hpp @@ -140,6 +140,12 @@ struct Parameters // Seed for RNG int rng_seed = 2006; + bool moving_endpoints = false; + bool translating_endpoints = false; + + float delta_Rx_left = 1.0; + float delta_Rx_right = 1.0; + // ----------------- Output -------------- bool output_energies_step = false; bool output_energies_divide_by_nspins = true; diff --git a/ui-cpp/ui-imgui/src/parameters_widget.cpp b/ui-cpp/ui-imgui/src/parameters_widget.cpp index 9d706513f..33203fe30 100644 --- a/ui-cpp/ui-imgui/src/parameters_widget.cpp +++ b/ui-cpp/ui-imgui/src/parameters_widget.cpp @@ -550,6 +550,41 @@ void ParametersWidget::show_content() { Parameters_GNEB_Set_Path_Shortening_Constant( state.get(), parameters_gneb.path_shortening_constant ); } + + ImGui::BeginTable( "MovingEndpointsTable", 2 ); + ImGui::TableNextColumn(); + ImGui::TextUnformatted( "Moving endpoints" ); + ImGui::TableNextColumn(); + if( ImGui::Checkbox( "##gneb_moving_endpoints", ¶meters_gneb.moving_endpoints ) ) + { + Parameters_GNEB_Set_Moving_Endpoints( state.get(), parameters_gneb.moving_endpoints ); + } + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex( 0 ); + ImGui::TextUnformatted( "Translating endpoints" ); + ImGui::TableNextColumn(); + if( ImGui::Checkbox( "##gneb_translating_endpoints", ¶meters_gneb.translating_endpoints ) ) + { + Parameters_GNEB_Set_Translating_Endpoints( state.get(), parameters_gneb.translating_endpoints ); + } + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex( 0 ); + ImGui::TextUnformatted( "Delta Rx left" ); + ImGui::TableNextColumn(); + if( ImGui::InputFloat( "##gneb_delta_Rx_left", ¶meters_gneb.delta_Rx_left ) ) + { + Parameters_GNEB_Set_Equilibrium_Delta_Rx( + state.get(), parameters_gneb.delta_Rx_left, parameters_gneb.delta_Rx_right ); + } + ImGui::TableNextColumn(); + ImGui::TextUnformatted( "Delta Rx right" ); + ImGui::TableNextColumn(); + if( ImGui::InputFloat( "##gneb_delta_Rx_right", ¶meters_gneb.delta_Rx_right ) ) + { + Parameters_GNEB_Set_Equilibrium_Delta_Rx( + state.get(), parameters_gneb.delta_Rx_left, parameters_gneb.delta_Rx_right ); + } + ImGui::EndTable(); } else if( ui_shared_state.selected_mode == GUI_Mode::MMF ) { @@ -768,6 +803,11 @@ void ParametersWidget::update_data() state.get(), ¶meters_gneb.output_energies_step, ¶meters_gneb.output_energies_interpolated, ¶meters_gneb.output_energies_divide_by_nspins, ¶meters_gneb.output_energies_add_readability_lines ); + parameters_gneb.moving_endpoints = Parameters_GNEB_Get_Moving_Endpoints( state.get() ); + parameters_gneb.translating_endpoints = Parameters_GNEB_Get_Translating_Endpoints( state.get() ); + Parameters_GNEB_Get_Equilibrium_Delta_Rx( + state.get(), ¶meters_gneb.delta_Rx_left, ¶meters_gneb.delta_Rx_right ); + // ----------------- MMF Parameters_MMF_Get_N_Iterations( state.get(), ¶meters_mmf.n_iterations, ¶meters_mmf.n_iterations_log ); // Output diff --git a/ui-cpp/ui-qt/include/PlotWidget.hpp b/ui-cpp/ui-qt/include/PlotWidget.hpp index feb62159b..23e6a161a 100644 --- a/ui-cpp/ui-qt/include/PlotWidget.hpp +++ b/ui-cpp/ui-qt/include/PlotWidget.hpp @@ -20,6 +20,8 @@ class PlotWidget : public QtCharts::QChartView bool plot_image_energies; bool plot_interpolated; int plot_interpolated_n; + bool divide_by_nos; + bool renormalize_Rx; private: std::shared_ptr state; diff --git a/ui-cpp/ui-qt/src/InfoWidget.cpp b/ui-cpp/ui-qt/src/InfoWidget.cpp index 5c4d7ee7a..242f1acaf 100644 --- a/ui-cpp/ui-qt/src/InfoWidget.cpp +++ b/ui-cpp/ui-qt/src/InfoWidget.cpp @@ -50,10 +50,10 @@ void InfoWidget::updateData() // Magnetization float M[3]; - Quantity_Get_Magnetization( state.get(), M ); - this->m_Label_Mx->setText( QString::fromLatin1( "Mx: " ) + QString::number( M[0], 'f', 8 ) ); - this->m_Label_My->setText( QString::fromLatin1( "My: " ) + QString::number( M[1], 'f', 8 ) ); - this->m_Label_Mz->setText( QString::fromLatin1( "Mz: " ) + QString::number( M[2], 'f', 8 ) ); + Quantity_Get_Average_Spin( state.get(), M ); + this->m_Label_Mx->setText( QString::fromLatin1( "Sx: " ) + QString::number( M[0], 'f', 8 ) ); + this->m_Label_My->setText( QString::fromLatin1( "Sy: " ) + QString::number( M[1], 'f', 8 ) ); + this->m_Label_Mz->setText( QString::fromLatin1( "Sz: " ) + QString::number( M[2], 'f', 8 ) ); // Force double f_max = Simulation_Get_MaxTorqueNorm( state.get() ); diff --git a/ui-cpp/ui-qt/src/MainWindow.cpp b/ui-cpp/ui-qt/src/MainWindow.cpp index 79e278813..28390f87b 100644 --- a/ui-cpp/ui-qt/src/MainWindow.cpp +++ b/ui-cpp/ui-qt/src/MainWindow.cpp @@ -870,8 +870,8 @@ void MainWindow::updateStatusBar() QString::fromLatin1( "E: " ) + QString::number( E, 'f', 6 ) + QString::fromLatin1( " " ) ); float M[3]; - Quantity_Get_Magnetization( state.get(), M ); - this->m_Label_Mz->setText( QString::fromLatin1( "M_z: " ) + QString::number( M[2], 'f', 6 ) ); + Quantity_Get_Average_Spin( state.get(), M ); + this->m_Label_Mz->setText( QString::fromLatin1( "S_z: " ) + QString::number( M[2], 'f', 6 ) ); float ips; int precision; diff --git a/ui-cpp/ui-qt/src/ParametersWidget.cpp b/ui-cpp/ui-qt/src/ParametersWidget.cpp index 9c8843698..4bfd2e3e0 100644 --- a/ui-cpp/ui-qt/src/ParametersWidget.cpp +++ b/ui-cpp/ui-qt/src/ParametersWidget.cpp @@ -189,6 +189,16 @@ void ParametersWidget::Load_Parameters_Contents() else if( image_type == 3 ) this->radioButton_Stationary->setChecked( true ); + // Moving endpoints + bool moving_endpoints = Parameters_GNEB_Get_Moving_Endpoints( state.get() ); + this->checkBox_moving_endpoints->setChecked( moving_endpoints ); + bool translating_endpoints = Parameters_GNEB_Get_Translating_Endpoints( state.get() ); + this->checkBox_translating_endpoints->setChecked( translating_endpoints ); + float delta_Rx_left, delta_Rx_right; + Parameters_GNEB_Get_Equilibrium_Delta_Rx( state.get(), &delta_Rx_left, &delta_Rx_right); + this->doubleSpinBox_delta_Rx_left->setValue( delta_Rx_left ); + this->doubleSpinBox_delta_Rx_right->setValue( delta_Rx_right ); + // EMA // modes to calculate and visualize i1 = Parameters_EMA_Get_N_Modes( state.get() ); @@ -445,6 +455,16 @@ void ParametersWidget::set_parameters_gneb() Parameters_GNEB_Set_N_Iterations( state.get(), i1, i2 ); std::string folder = this->lineEdit_gneb_output_folder->text().toStdString(); Parameters_GNEB_Set_Output_Folder( state.get(), folder.c_str() ); + + // Moving endpoints + bool moving_endpoints = this->checkBox_moving_endpoints->isChecked(); + bool translating_endpoints = this->checkBox_translating_endpoints->isChecked(); + float delta_Rx_left = this->doubleSpinBox_delta_Rx_left->value(); + float delta_Rx_right = this->doubleSpinBox_delta_Rx_right->value(); + + Parameters_GNEB_Set_Moving_Endpoints(state.get(), moving_endpoints); + Parameters_GNEB_Set_Translating_Endpoints(state.get(), translating_endpoints); + Parameters_GNEB_Set_Equilibrium_Delta_Rx(state.get(), delta_Rx_left, delta_Rx_right); } void ParametersWidget::set_gneb_auto_image_type() @@ -580,7 +600,7 @@ void ParametersWidget::load_Spin_Configuration_Eigenmodes() .c_str() ); auto file = string_q2std( fileName ); - IO_Eigenmodes_Read( this->state.get(), file.c_str(), type ); + IO_Eigenmodes_Read( this->state.get(), file.c_str() ); // n_modes parameter might be change by IO_Eigenmodes_Read so update that first this->spinBox_ema_n_modes->setValue( Parameters_EMA_Get_N_Modes( state.get() ) ); @@ -691,6 +711,7 @@ void ParametersWidget::Setup_Parameters_Slots() connect( this->checkBox_gneb_output_any, SIGNAL( stateChanged( int ) ), this, SLOT( set_parameters_gneb() ) ); connect( this->checkBox_gneb_output_initial, SIGNAL( stateChanged( int ) ), this, SLOT( set_parameters_gneb() ) ); connect( this->checkBox_gneb_output_final, SIGNAL( stateChanged( int ) ), this, SLOT( set_parameters_gneb() ) ); + connect( this->checkBox_gneb_output_energies_step, SIGNAL( stateChanged( int ) ), this, SLOT( set_parameters_gneb() ) ); connect( @@ -699,6 +720,12 @@ void ParametersWidget::Setup_Parameters_Slots() connect( this->checkBox_gneb_output_chain_step, SIGNAL( stateChanged( int ) ), this, SLOT( set_parameters_gneb() ) ); + // Endpoints + connect( this->checkBox_moving_endpoints, SIGNAL( stateChanged( int ) ), this, SLOT( set_parameters_gneb() ) ); + connect( this->checkBox_translating_endpoints, SIGNAL( stateChanged( int ) ), this, SLOT( set_parameters_gneb() ) ); + connect( this->doubleSpinBox_delta_Rx_left, SIGNAL( editingFinished() ), this, SLOT( set_parameters_gneb() ) ); + connect( this->doubleSpinBox_delta_Rx_right, SIGNAL( editingFinished() ), this, SLOT( set_parameters_gneb() ) ); + // EMA connect( this->spinBox_ema_n_modes, SIGNAL( editingFinished() ), this, SLOT( set_parameters_ema() ) ); connect( this->spinBox_ema_n_mode_follow, SIGNAL( editingFinished() ), this, SLOT( set_parameters_ema() ) ); diff --git a/ui-cpp/ui-qt/src/PlotWidget.cpp b/ui-cpp/ui-qt/src/PlotWidget.cpp index 323a55258..8ce4d45cd 100644 --- a/ui-cpp/ui-qt/src/PlotWidget.cpp +++ b/ui-cpp/ui-qt/src/PlotWidget.cpp @@ -177,14 +177,30 @@ void PlotWidget::plotEnergies() // Add data to series int idx_current = System_Get_Index( state.get() ); - current.push_back( QPointF( Rx[idx_current] / Rx_tot, energies[idx_current] / nos ) ); + + scalar Rx_cur = Rx[idx_current]; + if (renormalize_Rx && Rx_tot > 0) + Rx_cur /= Rx_tot; + + scalar E_cur = energies[idx_current]; + if (divide_by_nos) + E_cur /= nos; + + current.push_back( QPointF( Rx_cur, E_cur) ); + if( this->plot_image_energies ) { for( int i = 0; i < noi; ++i ) { if( i > 0 && Rx_tot > 0 ) - Rx[i] = Rx[i] / Rx_tot; - energies[i] = energies[i] / nos; + Rx[i] = Rx[i]; + energies[i] = energies[i]; + + if (renormalize_Rx && Rx_tot > 0) + Rx[i] /= Rx_tot; + + if (divide_by_nos) + energies[i] /= nos; if( Parameters_GNEB_Get_Climbing_Falling( state.get(), i ) == 0 ) normal.push_back( QPointF( Rx[i], energies[i] ) ); @@ -206,8 +222,14 @@ void PlotWidget::plotEnergies() for( int i = 0; i < size_interp; ++i ) { if( i > 0 && Rx_tot > 0 ) - Rx_interp[i] = Rx_interp[i] / Rx_tot; - energies_interp[i] = energies_interp[i] / nos; + Rx_interp[i] = Rx_interp[i]; + energies_interp[i] = energies_interp[i]; + + if (renormalize_Rx && Rx_tot > 0) + Rx_interp[i] /= Rx_tot; + + if (divide_by_nos) + energies_interp[i] /= nos; interp.push_back( QPointF( Rx_interp[i], energies_interp[i] ) ); @@ -268,6 +290,21 @@ void PlotWidget::plotEnergies() float delta = 0.1 * ( ymax - ymin ); if( delta < 1e-6 ) delta = 0.1; + this->chart->axisY()->setMin( ymin - delta ); this->chart->axisY()->setMax( ymax + delta ); + + // Rescale x axis + if (!renormalize_Rx && Rx_tot > 0) + { + delta = 0.04 * Rx_tot; + this->chart->axisX()->setMin( Rx[0] - delta ); + this->chart->axisX()->setMax( Rx[noi-1] + delta ); + } else if(Rx_tot > 0) { + this->chart->axisX()->setMin( -0.04 ); + this->chart->axisX()->setMax( 1.04 ); + } else { + this->chart->axisX()->setMin( -0.04 ); + this->chart->axisX()->setMax( 0.04 ); + } } \ No newline at end of file diff --git a/ui-cpp/ui-qt/src/PlotsWidget.cpp b/ui-cpp/ui-qt/src/PlotsWidget.cpp index 753efe9a9..0546dcaa3 100644 --- a/ui-cpp/ui-qt/src/PlotsWidget.cpp +++ b/ui-cpp/ui-qt/src/PlotsWidget.cpp @@ -24,6 +24,8 @@ PlotsWidget::PlotsWidget( std::shared_ptr state ) connect( this->checkBox_ImageEnergies, SIGNAL( stateChanged( int ) ), this, SLOT( updatePlotSettings() ) ); connect( this->checkBox_InterpolateEnergies, SIGNAL( stateChanged( int ) ), this, SLOT( updatePlotSettings() ) ); connect( this->spinBox_InterpolateEnergies_N, SIGNAL( editingFinished() ), this, SLOT( updatePlotSettings() ) ); + connect( this->checkBox_divide_by_NOS, SIGNAL( stateChanged( int ) ), this, SLOT( updatePlotSettings() ) ); + connect( this->checkBox_normalize_Rx, SIGNAL( stateChanged( int ) ), this, SLOT( updatePlotSettings() ) ); // Update Timer auto timer = new QTimer( this ); @@ -46,5 +48,9 @@ void PlotsWidget::updatePlotSettings() { this->energyPlot->plot_image_energies = this->checkBox_ImageEnergies->isChecked(); this->energyPlot->plot_interpolated = this->checkBox_InterpolateEnergies->isChecked(); + + this->energyPlot->divide_by_nos = this->checkBox_divide_by_NOS->isChecked(); + this->energyPlot->renormalize_Rx = this->checkBox_normalize_Rx->isChecked(); + this->energyPlot->plot_interpolated_n = this->spinBox_InterpolateEnergies_N->value(); } \ No newline at end of file diff --git a/ui-cpp/ui-qt/ui/ParametersWidget.ui b/ui-cpp/ui-qt/ui/ParametersWidget.ui index 1983806fb..892f8f6ce 100644 --- a/ui-cpp/ui-qt/ui/ParametersWidget.ui +++ b/ui-cpp/ui-qt/ui/ParametersWidget.ui @@ -83,15 +83,15 @@ - 0 + 2 0 0 - 513 - 670 + 533 + 684 @@ -110,22 +110,6 @@ LLG - - - - Qt::Vertical - - - QSizePolicy::Fixed - - - - 20 - 30 - - - - @@ -138,42 +122,26 @@ Parameters - - - - - - - 60 - 16777215 - - - - Inclination - - - - - - - - 60 - 16777215 - - + + + + - 0 + Polarisation (x y z) - - + + Qt::Horizontal + + QSizePolicy::Minimum + - 40 + 15 20 @@ -181,6 +149,13 @@ + + + + Temperature + + + @@ -250,6 +225,107 @@ + + + + Convergence at + + + + + + + Qt::Horizontal + + + + 40 + 20 + + + + + + + + Time Step + + + + + + + + + + 45 + 16777215 + + + + 3 + + + -100.000000000000000 + + + 100.000000000000000 + + + + + + + + 45 + 16777215 + + + + 3 + + + -100.000000000000000 + + + 100.000000000000000 + + + + + + + + 45 + 16777215 + + + + 3 + + + -100.000000000000000 + + + 100.000000000000000 + + + + + + + Qt::Horizontal + + + + 40 + 20 + + + + + + @@ -285,6 +361,22 @@ + + + + + 40 + 0 + + + + + 50 + 16777215 + + + + @@ -325,62 +417,6 @@ - - - - Convergence at - - - - - - - - - - 15 - 16777215 - - - - 1E - - - - - - - - 60 - 0 - - - - -20 - - - 20 - - - -8 - - - - - - - Qt::Horizontal - - - - 40 - 20 - - - - - - @@ -442,74 +478,36 @@ - - - - Time Step - - - - - - - - - - 45 - 16777215 - - - - 3 - - - -100.000000000000000 - - - 100.000000000000000 - - - - - + + + + - 45 + 60 16777215 - - 3 - - - -100.000000000000000 - - - 100.000000000000000 + + Inclination - - + + - 45 + 60 16777215 - - 3 - - - -100.000000000000000 - - - 100.000000000000000 + + 0 - - + + Qt::Horizontal @@ -523,55 +521,48 @@ - - - - - 40 - 0 - - - - - 50 - 16777215 - - - - - - - - Qt::Horizontal - - - - 40 - 20 - - - - - - - - + + + + + + + 15 + 16777215 + + - Polarisation (x y z) + 1E - - + + + + + 60 + 0 + + + + -20 + + + 20 + + + -8 + + + + + Qt::Horizontal - - QSizePolicy::Minimum - - 15 + 40 20 @@ -579,13 +570,6 @@ - - - - Temperature - - - @@ -601,7 +585,7 @@ Output - + @@ -626,7 +610,7 @@ - + @@ -660,7 +644,7 @@ - + @@ -706,7 +690,7 @@ - + Qt::Horizontal @@ -719,7 +703,7 @@ - + @@ -763,6 +747,22 @@ + + + + Qt::Vertical + + + QSizePolicy::Fixed + + + + 20 + 30 + + + + @@ -842,7 +842,7 @@ 0 0 483 - 637 + 652 @@ -1193,8 +1193,8 @@ 0 0 - 483 - 715 + 469 + 848 @@ -1213,7 +1213,7 @@ GNEB - + Qt::Vertical @@ -1226,23 +1226,7 @@ - - - - Qt::Vertical - - - QSizePolicy::Fixed - - - - 20 - 30 - - - - - + @@ -1303,31 +1287,31 @@ - - + + - Energies Step + Chain Step - - + + - Energies divide by NOS + Energies Step - - + + - Chain Step + Energies Interpolated - - + + - Energies Interpolated + Energies divide by NOS @@ -1405,65 +1389,57 @@ - - - - - 0 - 0 - + + + + Qt::Vertical - + - 16777215 - 16777215 + 20 + 40 - - GNEB Parameters - - - - - - - - Qt::Horizontal - - - - 40 - 20 - - - - - - - - false - - - Temperature - - - - - - - 4 - - - 1000.000000000000000 - - - - - - - + + + + + + Qt::Vertical + + + QSizePolicy::Fixed + + + + 20 + 30 + + + + + + + + + 0 + 0 + + + + + 16777215 + 16777215 + + + + GNEB Parameters + + + + - + Qt::Horizontal @@ -1476,28 +1452,41 @@ - - - false - - + + 60 - 16777215 + 0 + + -20 + + + 20 + + + -8 + - + - Time Step + Convergence at 1E - + + + + Moving endpoints + + + + @@ -1563,10 +1552,10 @@ - - + + - + Qt::Horizontal @@ -1579,63 +1568,57 @@ - + - 60 + 40 0 - - -20 - - - 20 + + + 50 + 16777215 + - - -8 + + 0.5 - + - Convergence at 1E + Spring force ratio - - + + - + - Spring Constant + Path shortening constant 1E - - - - 40 - 0 - + + + -1000 - - - 50 - 16777215 - + + 1000 - - 1.0 + + -3 - + Qt::Horizontal @@ -1649,30 +1632,80 @@ - - + + + + + + Qt::Horizontal + + + + 40 + 20 + + + + + + + + false + + + Temperature + + + + + + + 4 + + + 1000.000000000000000 + + + + + + + + + Translating endpoints + + + + + - + - Path shortening constant 1E + Spring Constant - - - -1000 + + + + 40 + 0 + - - 1000 + + + 50 + 16777215 + - - -3 + + 1.0 - + Qt::Horizontal @@ -1686,10 +1719,10 @@ - - + + - + Qt::Horizontal @@ -1702,33 +1735,73 @@ - - - - 40 - 0 - + + + false - 50 + 60 16777215 - - 0.5 - - + - Spring force ratio + Time Step + + + + Equilibrium Delta_Rx (right) + + + + + + + + 100 + 0 + + + + + 100 + 16777215 + + + + + + + + + 100 + 0 + + + + + 100 + 16777215 + + + + + + + + Equilibrium Delta_Rx (left) + + + @@ -1740,7 +1813,7 @@ 0 0 483 - 637 + 652 @@ -2055,7 +2128,7 @@ 0 0 483 - 637 + 652 @@ -2463,7 +2536,7 @@ - + diff --git a/ui-cpp/ui-qt/ui/PlotsWidget.ui b/ui-cpp/ui-qt/ui/PlotsWidget.ui index 775ab7276..1e2d45549 100644 --- a/ui-cpp/ui-qt/ui/PlotsWidget.ui +++ b/ui-cpp/ui-qt/ui/PlotsWidget.ui @@ -110,28 +110,18 @@ - - - - Qt::Horizontal - - - - 40 - 20 - - - - - - - - Refresh - - - + + + + 1 + + + 10000 + + + @@ -142,6 +132,13 @@ + + + + Normalize Rx + + + @@ -149,18 +146,35 @@ - - - - 1 - - - 10000 + + + + Divide by NOS + + + + Qt::Horizontal + + + + 40 + 20 + + + + + + + + Refresh + + +