Merge branch 'moving_endpoints_gneb' of github.com:MSallermann/spirit…

… into develop

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/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/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/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

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<Data::Spin_System> image;
+    std::shared_ptr<Data::Spin_System_Chain> 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<Data::Spin_System> image;
+    std::shared_ptr<Data::Spin_System_Chain> 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<Data::Spin_System> image;
+    std::shared_ptr<Data::Spin_System_Chain> 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<Data::Spin_System> image;
+    std::shared_ptr<Data::Spin_System_Chain> chain;
+
+    // Fetch correct indices and pointers
+    from_indices( state, idx_image, idx_chain, image, chain );
+
+    auto p = chain->gneb_parameters;
+    return static_cast<bool>( 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<Data::Spin_System> image;
+    std::shared_ptr<Data::Spin_System_Chain> chain;
+
+    // Fetch correct indices and pointers
+    from_indices( state, idx_image, idx_chain, image, chain );
+
+    auto p = chain->gneb_parameters;
+    return static_cast<bool>( 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<Data::Spin_System> image;
+    std::shared_ptr<Data::Spin_System_Chain> 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
 {

core/src/Spirit/Simulation.cpp

@@ -255,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

core/src/engine/Manifoldmath.cpp

@@ -3,6 +3,7 @@
 #include <utility/Constants.hpp>
 #include <utility/Exception.hpp>
 #include <utility/Logging.hpp>
+#include <engine/Backend_par.hpp>
 
 #include <Eigen/Dense>
 
@@ -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