Rough implementation of transition plane sampling

core/include/engine/CMakeLists.txt

@@ -21,6 +21,7 @@ set(HEADER_SPIRIT_ENGINE
     ${CMAKE_CURRENT_SOURCE_DIR}/Method_MC.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/Method_MMF.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/Method_EMA.hpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/Method_TS_Sampling.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/Vectormath_Defines.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/Vectormath.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/Manifoldmath.hpp

core/include/engine/Manifoldmath.hpp

@@ -32,6 +32,11 @@ void project_parallel( vectorfield & vf1, const vectorfield & vf2 );
 // Project v1 to be orthogonal to v2
 //      Note: this assumes normalized vectorfields
 void project_orthogonal( vectorfield & vf1, const vectorfield & vf2 );
+
+// Project v1 to be orthogonal to v2
+//      Note: this assumes normalized vectorfields
+void project_orthogonal( vectorfield & vf1, const vectorfield & vf2, scalar dot_product );
+
 // Invert v1's component parallel to v2
 //      Note: this assumes normalized vectorfields
 void invert_parallel( vectorfield & vf1, const vectorfield & vf2 );

core/include/engine/Method_MC.hpp

@@ -25,11 +25,18 @@ class Method_MC : public Method
     // Method name as string
     std::string Name() override;
 
-private:
-    // Solver_Iteration represents one iteration of a certain Solver
-    void Iteration() override;
+protected:
+    virtual void Displace_Spin(int ispin, vectorfield & spins_new, std::uniform_real_distribution<scalar> & distribution, std::vector<int> & changed_indices, vectorfield & old_spins);
+    virtual void Spin_Trial(int spin, vectorfield & spins_new, const vectorfield & spins_old, std::uniform_real_distribution<scalar> & distribution );
+    virtual scalar Compute_Energy_Diff(const std::vector<int> & changed_indices, vectorfield & spins_new, const vectorfield & spins_old);
+    virtual void Reject(const std::vector<int> & rejected_indices, const vectorfield & spins_old);
+
+    std::shared_ptr<Data::Parameters_Method_MC> parameters_mc;
+    // Vector to save the previous spin directions
+    vectorfield spins_new;
 
-    void Displace_Spin(int ispin, vectorfield & spins_new, const vectorfield & spins_old, std::uniform_real_distribution<scalar> & distribution);
+    // Solver_Iteration represents one iteration of a certain Solver
+    virtual void Iteration() override;
 
     // Metropolis iteration with adaptive cone radius
     void Metropolis( const vectorfield & spins_old, vectorfield & spins_new );
@@ -51,7 +58,6 @@ class Method_MC : public Method
     void Message_Step() override;
     void Message_End() override;
 
-    std::shared_ptr<Data::Parameters_Method_MC> parameters_mc;
 
     // Cosine of current cone angle
     scalar cone_angle;
@@ -61,9 +67,6 @@ class Method_MC : public Method
 
     // Random vector array
     vectorfield xi;
-
-    // Vector to save the previous spin directions
-    vectorfield spins_new;
 };
 
 } // namespace Engine

core/include/engine/Method_TS_Sampling.hpp

@@ -0,0 +1,34 @@
+#pragma once
+#ifndef SPIRIT_CORE_ENGINE_METHOD_TS_SAMPLING_HPP
+#define SPIRIT_CORE_ENGINE_METHOD_TS_SAMPLING_HPP
+
+#include <engine/Method_MC.hpp>
+#include <deque>
+
+
+namespace Engine
+{
+/*
+    Transition plane sampling
+*/
+class Method_TS_Sampling : public Method_MC
+{
+    public:
+    Method_TS_Sampling( std::shared_ptr<Data::Spin_System> system, int idx_img, int idx_chain );
+
+    void Set_Transition_Plane_Normal( vectorfield & spins_minimum, vectorfield & unstable_mode );
+    void Compute_Transition_Plane_Normal( );
+
+    std::string Name() override;
+
+    private:
+    void Iteration() override;
+    // void Displace_Spin(int ispin, vectorfield & spins_new, std::uniform_real_distribution<scalar> & distribution, std::vector<int> & changed_indices, vectorfield & old_spins) override;
+
+    vectorfield transition_plane_normal;
+    std::uniform_int_distribution<> distribution_idx;
+    std::deque<int> rejected;
+};
+
+}
+#endif

core/src/Spirit/Simulation.cpp

@@ -7,6 +7,7 @@
 #include <engine/Method_GNEB.hpp>
 #include <engine/Method_LLG.hpp>
 #include <engine/Method_MC.hpp>
+#include <engine/Method_TS_Sampling.hpp>
 #include <engine/Method_MMF.hpp>
 #include <utility/Exception.hpp>
 #include <utility/Logging.hpp>
@@ -117,6 +118,7 @@ try
             image->mc_parameters->n_iterations_log = n_iterations_log;
 
         auto method = std::shared_ptr<Engine::Method>( new Engine::Method_MC( image, idx_image, idx_chain ) );
+        // auto method = std::shared_ptr<Engine::Method>( new Engine::Method_TS_Sampling( image, idx_image, idx_chain ) );
 
         image->Unlock();
 

core/src/engine/CMakeLists.txt

@@ -14,6 +14,7 @@ set(SOURCE_SPIRIT_ENGINE
     ${CMAKE_CURRENT_SOURCE_DIR}/Method_MC.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/Method_MMF.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/Method_EMA.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/Method_TS_Sampling.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/Vectormath.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/Vectormath.cu
     ${CMAKE_CURRENT_SOURCE_DIR}/Manifoldmath.cpp

core/src/engine/Manifoldmath.cpp

@@ -26,17 +26,22 @@ void project_parallel( vectorfield & vf1, const vectorfield & 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 )
+void project_orthogonal( vectorfield & vf1, const vectorfield & vf2, scalar proj )
 {
-    scalar x = Vectormath::dot( vf1, vf2 );
 // TODO: replace the loop with Vectormath Kernel
 #pragma omp parallel for
     for( unsigned int i = 0; i < vf1.size(); ++i )
-        vf1[i] -= x * vf2[i];
+        vf1[i] -= proj * vf2[i];
+}
+
+void project_orthogonal( vectorfield & vf1, const vectorfield & vf2 )
+{
+    scalar x = Vectormath::dot( vf1, vf2 );
+    project_orthogonal(vf1, vf2, x);
 }
 
 void invert_parallel( vectorfield & vf1, const vectorfield & vf2 )

core/src/engine/Method_MC.cpp

@@ -66,7 +66,7 @@ void Method_MC::Iteration()
     Vectormath::set_c_a( 1, this->spins_new, spins_old );
 }
 
-void Method_MC::Displace_Spin(int ispin, vectorfield & spins_new, const vectorfield & spins_old, std::uniform_real_distribution<scalar> & distribution)
+void Method_MC::Displace_Spin(int ispin, vectorfield & spins_new, std::uniform_real_distribution<scalar> & distribution, std::vector<int> & changed_indices, vectorfield & old_spins)
 {
     // One Metropolis step for each spin
     const Vector3 e_z{ 0, 0, 1 };
@@ -78,9 +78,9 @@ void Method_MC::Displace_Spin(int ispin, vectorfield & spins_new, const vectorfi
     if( this->parameters_mc->metropolis_step_cone )
     {
         // Calculate local basis for the spin
-        if( spins_old[ispin].z() < 1 - 1e-10 )
+        if( spins_new[ispin].z() < 1 - 1e-10 )
         {
-            local_basis.col( 2 ) = spins_old[ispin];
+            local_basis.col( 2 ) = spins_new[ispin];
             local_basis.col( 0 ) = ( local_basis.col( 2 ).cross( e_z ) ).normalized();
             local_basis.col( 1 ) = local_basis.col( 2 ).cross( local_basis.col( 0 ) );
         }
@@ -119,12 +119,76 @@ void Method_MC::Displace_Spin(int ispin, vectorfield & spins_new, const vectorfi
     }
 };
 
+scalar Method_MC::Compute_Energy_Diff(const std::vector<int> & changed_indices, vectorfield & spins_new, const vectorfield & spins_old)
+{
+    // Energy difference of configurations with and without displacement
+    scalar Eold = 0;
+    for(auto ispin : changed_indices)
+        Eold += this->systems[0]->hamiltonian->Energy_Single_Spin( ispin, spins_old );
+
+    scalar Enew = 0;
+    for(auto ispin : changed_indices)
+        Enew  += this->systems[0]->hamiltonian->Energy_Single_Spin( ispin, spins_new );
+
+    scalar Ediff = Enew - Eold;
+    return Ediff;
+}
+
+void Method_MC::Reject(const std::vector<int> & rejected_indices, const vectorfield & spins_old)
+{
+    // Restore the spin
+    for(int i=0; i<rejected_indices.size(); i++)
+    {
+        spins_new[rejected_indices[i]] = spins_old[i];
+    }
+    // Counter for the number of rejections
+    ++this->n_rejected;
+}
+
+void Method_MC::Spin_Trial(int ispin, vectorfield & spins_new, const vectorfield & spins_old, std::uniform_real_distribution<scalar> & distribution )
+{
+    const scalar kB_T = Constants::k_B * this->parameters_mc->temperature;
+
+    auto indices   = std::vector<int>{ispin};
+    auto old_spins = std::vector<Vector3>{spins_new[ispin]};
+
+    Displace_Spin( ispin, spins_new, distribution, indices, old_spins);
+
+    // Energy difference of configurations with and without displacement
+    scalar Ediff = Compute_Energy_Diff(indices, spins_new, spins_old);
+
+    bool reject = false;
+
+    // Metropolis criterion: reject the step if energy rose
+    if( Ediff > 1e-14 )
+    {
+        if( this->parameters_mc->temperature < 1e-12 )
+        {
+            reject = true;
+        }
+        else
+        {
+            // Exponential factor
+            scalar exp_ediff = std::exp( -Ediff / kB_T );
+            // Metropolis random number
+            scalar x_metropolis = distribution( this->parameters_mc->prng );
+            // Only reject if random number is larger than exponential
+            if( exp_ediff < x_metropolis )
+            {
+                reject = true;
+            }
+        }
+    }
+
+    if (reject)
+        Reject(indices, old_spins);
+}
+
 // Simple metropolis step
 void Method_MC::Metropolis( const vectorfield & spins_old, vectorfield & spins_new )
 {
     auto distribution     = std::uniform_real_distribution<scalar>( 0, 1 );
     auto distribution_idx = std::uniform_int_distribution<>( 0, this->nos - 1 );
-    scalar kB_T           = Constants::k_B * this->parameters_mc->temperature;
 
     scalar diff = 0.01;
 
@@ -158,41 +222,7 @@ void Method_MC::Metropolis( const vectorfield & spins_old, vectorfield & spins_n
 
         if( Vectormath::check_atom_type( this->systems[0]->geometry->atom_types[ispin] ) )
         {
-
-            Displace_Spin( ispin, spins_new, spins_old, distribution );
-
-            // Energy difference of configurations with and without displacement
-            scalar Eold  = this->systems[0]->hamiltonian->Energy_Single_Spin( ispin, spins_old );
-            scalar Enew  = this->systems[0]->hamiltonian->Energy_Single_Spin( ispin, spins_new );
-            scalar Ediff = Enew - Eold;
-
-            // Metropolis criterion: reject the step if energy rose
-            if( Ediff > 1e-14 )
-            {
-                if( this->parameters_mc->temperature < 1e-12 )
-                {
-                    // Restore the spin
-                    spins_new[ispin] = spins_old[ispin];
-                    // Counter for the number of rejections
-                    ++this->n_rejected;
-                }
-                else
-                {
-                    // Exponential factor
-                    scalar exp_ediff = std::exp( -Ediff / kB_T );
-                    // Metropolis random number
-                    scalar x_metropolis = distribution( this->parameters_mc->prng );
-
-                    // Only reject if random number is larger than exponential
-                    if( exp_ediff < x_metropolis )
-                    {
-                        // Restore the spin
-                        spins_new[ispin] = spins_old[ispin];
-                        // Counter for the number of rejections
-                        ++this->n_rejected;
-                    }
-                }
-            }
+            Spin_Trial(ispin, spins_new, spins_old, distribution);
         }
     }
 }