Merge branch 'develop' into fixed-micromagnetic

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

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.
 

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,

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

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<Vector3> bravais_vectors, intfield n_cells, std::vector<Vector3> cell_atoms,
-        Basis_Cell_Composition cell_composition, scalar lattice_constant, Pinning pinning, Defects defects );
+        const std::vector<Vector3> & bravais_vectors, intfield n_cells, const std::vector<Vector3> & cell_atoms,
+        const Basis_Cell_Composition & cell_composition, scalar lattice_constant, const Pinning & pinning,
+        const Defects & defects );
 
     // ---------- Convenience functions
     // Retrieve triangulation, if 2D

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

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;

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<Data::Parameters_Method> 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<Data::Parameters_Method> 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

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)

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;

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

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<Data::Geometry> geometry;
 

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,

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<int>    history_iteration;
+    std::vector<scalar> history_max_torque;
+    std::vector<scalar> 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<scalar> force_max_abs_component_all;
 
-    // History of relevant quantities
-    std::map<std::string, std::vector<scalar>> history;
-
     //////////// General /////////////////////////////////////////////////////////
 
     // Systems the Solver will access

core/include/engine/Method_GNEB.hpp

@@ -75,8 +75,15 @@ class Method_GNEB : public Method_Solver<solver>
     std::vector<vectorfield> F_gradient;
     std::vector<vectorfield> F_spring;
     vectorfield f_shrink;
+
+    vectorfield F_translation_left;
+    vectorfield F_translation_right;
+
     // Last calculated tangents
     std::vector<vectorfield> tangents;
+
+    vectorfield tangent_endpoints_left;
+    vectorfield tangent_endpoints_right;
 };
 
 } // namespace Engine

core/include/engine/Solver_Depondt.hpp

@@ -16,7 +16,7 @@ inline void Method_Solver<Solver::Depondt>::Initialize()
         configurations_predictor[i] = std::shared_ptr<vectorfield>( 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<Solver::Depondt>::Iteration()
         // Get new spin conf n_new = R( (H+H')/2 ) * n
         Vectormath::rotate( conf, temp1, angle, conf );
     }
-};
+}
 
 template<>
 inline std::string Method_Solver<Solver::Depondt>::SolverName()
 {
     return "Depondt";
-};
+}
 
 template<>
 inline std::string Method_Solver<Solver::Depondt>::SolverFullName()
 {
     return "Depondt";
-};
+}

core/include/engine/Solver_Heun.hpp

@@ -16,7 +16,7 @@ inline void Method_Solver<Solver::Heun>::Initialize()
         configurations_predictor[i] = std::shared_ptr<vectorfield>( 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<Solver::Heun>::Iteration()
         // Copy out
         conf = conf_temp;
     }
-};
+}
 
 template<>
 inline std::string Method_Solver<Solver::Heun>::SolverName()
 {
     return "Heun";
-};
+}
 
 template<>
 inline std::string Method_Solver<Solver::Heun>::SolverFullName()
 {
     return "Heun";
-};
+}

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 );
 

core/include/engine/Solver_RK4.hpp

@@ -32,7 +32,7 @@ inline void Method_Solver<Solver::RungeKutta4>::Initialize()
         this->configurations_k4[i] = std::shared_ptr<vectorfield>( 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<Solver::RungeKutta4>::Iteration()
         // Copy out
         conf = conf_temp;
     }
-};
+}
 
 template<>
 inline std::string Method_Solver<Solver::RungeKutta4>::SolverName()
 {
     return "RK4";
-};
+}
 
 template<>
 inline std::string Method_Solver<Solver::RungeKutta4>::SolverFullName()
 {
     return "Runge Kutta (4th order)";
-};
+}