diff --git a/core/docs/c-api/Quantities.md b/core/docs/c-api/Quantities.md index 7ea61fb7c..dee2ff93a 100644 --- a/core/docs/c-api/Quantities.md +++ b/core/docs/c-api/Quantities.md @@ -37,17 +37,3 @@ void Quantity_Get_Grad_Force_MinimumMode(State * state, float * gradient, float Minimum mode following information - - -### Quantity_Get_HTST_Prefactor - -```C -float Quantity_Get_HTST_Prefactor(State * state, int idx_image_minimum, int idx_image_sp, int idx_chain=-1) -``` - -HTST Prefactor for transition from minimum to saddle point. - -Note that the method assumes you gave it correct images, where the -gradient is zero and which correspond to a minimum and a saddle point -respectively. - diff --git a/core/include/Spirit/Quantities.h b/core/include/Spirit/Quantities.h index a3eade72d..e7b818b4d 100644 --- a/core/include/Spirit/Quantities.h +++ b/core/include/Spirit/Quantities.h @@ -22,13 +22,4 @@ PREFIX float Quantity_Get_Topological_Charge(State * state, int idx_image=-1, in // 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, int idx_chain=-1); -/* -HTST Prefactor for transition from minimum to saddle point. - -Note that the method assumes you gave it correct images, where the -gradient is zero and which correspond to a minimum and a saddle point -respectively. -*/ -PREFIX float Quantity_Get_HTST_Prefactor(State * state, int idx_image_minimum, int idx_image_sp, int idx_chain=-1); - #endif \ No newline at end of file diff --git a/core/include/engine/CMakeLists.txt b/core/include/engine/CMakeLists.txt index d754d4d4d..edcfaef46 100644 --- a/core/include/engine/CMakeLists.txt +++ b/core/include/engine/CMakeLists.txt @@ -5,7 +5,6 @@ set(HEADER_SPIRIT_ENGINE ${CMAKE_CURRENT_SOURCE_DIR}/Hamiltonian_Heisenberg.hpp ${CMAKE_CURRENT_SOURCE_DIR}/Hamiltonian_Gaussian.hpp ${CMAKE_CURRENT_SOURCE_DIR}/Eigenmodes.hpp - ${CMAKE_CURRENT_SOURCE_DIR}/HTST.hpp ${CMAKE_CURRENT_SOURCE_DIR}/Solver_SIB.hpp ${CMAKE_CURRENT_SOURCE_DIR}/Solver_Heun.hpp ${CMAKE_CURRENT_SOURCE_DIR}/Solver_Depondt.hpp diff --git a/core/include/engine/HTST.hpp b/core/include/engine/HTST.hpp deleted file mode 100644 index d9ece5fd4..000000000 --- a/core/include/engine/HTST.hpp +++ /dev/null @@ -1,34 +0,0 @@ -#pragma once -#ifndef HTST_H -#define HTST_H - -#include "Spirit_Defines.h" -#include -#include -#include - -namespace Engine -{ - namespace HTST - { - // Note the two images should correspond to one minimum and one saddle point - void Calculate_Prefactor(Data::HTST_Info & htst_info); - - // Calculate the 'a' component of the prefactor - void Calculate_Perpendicular_Velocity(const vectorfield & spins, const scalarfield & mu_s, const MatrixX & hessian, - const MatrixX & basis, const MatrixX & eigenbasis, VectorX & a); - - // Calculate the Velocity matrix - void Calculate_Dynamical_Matrix(const vectorfield & spins, const scalarfield & mu_s, const MatrixX & hessian, MatrixX & velocity); - - // Calculate the zero volume of a spin system - scalar Calculate_Zero_Volume(const std::shared_ptr system); - - // Generates the geodesic Hessian in 2N-representation and calculates it's eigenvalues and eigenvectors - void Geodesic_Eigen_Decomposition(const vectorfield & image, const vectorfield & gradient, const MatrixX & hessian, - MatrixX & hessian_geodesic_3N, MatrixX & hessian_geodesic_2N, VectorX & eigenvalues, MatrixX & eigenvectors); - }; -} - - -#endif \ No newline at end of file diff --git a/core/python/spirit/quantities.py b/core/python/spirit/quantities.py index 2268552f2..c7598d81d 100644 --- a/core/python/spirit/quantities.py +++ b/core/python/spirit/quantities.py @@ -27,14 +27,6 @@ def get_magnetization(p_state, idx_image=-1, idx_chain=-1): return [float(i) for i in magnetization] -_Get_HTST_Prefactor = _spirit.Quantity_Get_HTST_Prefactor -_Get_HTST_Prefactor.argtypes = [ctypes.c_void_p, ctypes.c_int, ctypes.c_int, ctypes.c_int] -_Get_HTST_Prefactor.restype = ctypes.c_float -def get_htst_prefactor(p_state, idx_image_minimum, idx_image_sp, idx_chain=-1): - """Calculates and returns the HTST rate prefactor.""" - return _Get_HTST_Prefactor(p_state, idx_image_minimum, idx_image_sp, idx_chain) - - ### Temporary info function for MMF _Get_MinimumMode = _spirit.Quantity_Get_Grad_Force_MinimumMode _Get_MinimumMode.argtypes = [ctypes.c_void_p, ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_float), diff --git a/core/src/Spirit/Quantities.cpp b/core/src/Spirit/Quantities.cpp index 4859c404b..2a579d99b 100644 --- a/core/src/Spirit/Quantities.cpp +++ b/core/src/Spirit/Quantities.cpp @@ -1,7 +1,6 @@ #include #include #include -#include #include #include #include @@ -70,29 +69,6 @@ float Quantity_Get_Topological_Charge(State * state, int idx_image, int idx_chai } -float Quantity_Get_HTST_Prefactor(State * state, int idx_image_minimum, int idx_image_sp, int idx_chain) -try -{ - std::shared_ptr image_minimum, image_sp; - std::shared_ptr chain; - from_indices(state, idx_image_minimum, idx_chain, image_minimum, chain); - from_indices(state, idx_image_sp, idx_chain, image_sp, chain); - - auto& info = chain->htst_info; - info.minimum = image_minimum; - info.saddle_point = image_sp; - - Engine::HTST::Calculate_Prefactor(chain->htst_info); - - return (float)info.prefactor; -} -catch( ... ) -{ - spirit_handle_exception_api(-1, 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) { using namespace Engine; diff --git a/core/src/engine/CMakeLists.txt b/core/src/engine/CMakeLists.txt index fb20c9b4d..5dac61362 100644 --- a/core/src/engine/CMakeLists.txt +++ b/core/src/engine/CMakeLists.txt @@ -6,7 +6,6 @@ set(SOURCE_SPIRIT_ENGINE ${CMAKE_CURRENT_SOURCE_DIR}/Hamiltonian_Heisenberg.cu ${CMAKE_CURRENT_SOURCE_DIR}/Hamiltonian_Gaussian.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Eigenmodes.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/HTST.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Method.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Method_LLG.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Method_GNEB.cpp diff --git a/core/src/engine/HTST.cpp b/core/src/engine/HTST.cpp deleted file mode 100644 index 3fc019234..000000000 --- a/core/src/engine/HTST.cpp +++ /dev/null @@ -1,519 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -// #include -#include -//#include -#include // Also includes -#include - -#include -#include - -namespace C = Utility::Constants; - -namespace Engine -{ - namespace HTST - { - // Note the two images should correspond to one minimum and one saddle point - // Non-extremal images may yield incorrect Hessians and thus incorrect results - void Calculate_Prefactor(Data::HTST_Info & htst_info) - { - Log(Utility::Log_Level::All, Utility::Log_Sender::HTST, "---- Prefactor calculation"); - - const scalar epsilon = 1e-4; - const scalar epsilon_force = 1e-8; - - auto& image_minimum = *htst_info.minimum->spins; - auto& image_sp = *htst_info.saddle_point->spins; - int nos = image_minimum.size(); - vectorfield force_tmp(nos, {0,0,0}); - std::vector block; - - // TODO - bool is_afm = false; - - // The gradient (unprojected) - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Evaluation of the gradient at the initial configuration..."); - vectorfield gradient_minimum(nos, {0,0,0}); - htst_info.minimum->hamiltonian->Gradient(image_minimum, gradient_minimum); - - // Check if the configuration is actually an extremum - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Checking if initial configuration is an extremum..."); - Vectormath::set_c_a(1, gradient_minimum, force_tmp); - Manifoldmath::project_tangential(force_tmp, image_minimum); - scalar fmax_minimum = Vectormath::max_abs_component(force_tmp); - if( fmax_minimum > epsilon_force ) - { - Log(Utility::Log_Level::Error, Utility::Log_Sender::All, fmt::format( - "HTST: the initial configuration is not a converged minimum, its max. force component is above the threshold ({} > {})!", fmax_minimum, epsilon_force )); - return; - } - - // The gradient (unprojected) - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Evaluation of the gradient at the transition configuration..."); - vectorfield gradient_sp(nos, {0,0,0}); - htst_info.saddle_point->hamiltonian->Gradient(image_sp, gradient_sp); - - // Check if the configuration is actually an extremum - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Checking if transition configuration is an extremum..."); - Vectormath::set_c_a(1, gradient_sp, force_tmp); - Manifoldmath::project_tangential(force_tmp, image_sp); - scalar fmax_sp = Vectormath::max_abs_component(force_tmp); - if( fmax_sp > epsilon_force ) - { - Log(Utility::Log_Level::Error, Utility::Log_Sender::All, fmt::format( - "HTST: the transition configuration is not a converged saddle point, its max. force component is above the threshold ({} > {})!", fmax_sp, epsilon_force )); - return; - } - - //////////////////////////////////////////////////////////////////////// - // Saddle point - { - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, "Calculation for the Saddle Point"); - - // Evaluation of the Hessian... - MatrixX hessian_sp = MatrixX::Zero(3*nos,3*nos); - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Evaluation of the Hessian..."); - htst_info.saddle_point->hamiltonian->Hessian(image_sp, hessian_sp); - - // Eigendecomposition - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Eigendecomposition..."); - MatrixX hessian_geodesic_sp_3N(3*nos, 3*nos); - MatrixX hessian_geodesic_sp_2N(2*nos, 2*nos); - VectorX eigenvalues_sp = VectorX::Zero(2*nos); - MatrixX eigenvectors_sp = MatrixX::Zero(2*nos, 2*nos); - Geodesic_Eigen_Decomposition(image_sp, gradient_sp, hessian_sp, - hessian_geodesic_sp_3N, hessian_geodesic_sp_2N, eigenvalues_sp, eigenvectors_sp); - - htst_info.eigenvalues_sp.assign(eigenvalues_sp.data(), eigenvalues_sp.data() + 2*nos); - // htst_info.eigenvectors_sp.assign(eigenvectors_sp.data(), eigenvectors_sp.data() + 2*3*nos); - - // Print some eigenvalues - block = std::vector{"10 lowest eigenvalues at saddle point:"}; - for( int i=0; i<10; ++i ) - block.push_back(fmt::format("ew[{}]={:^20e} ew[{}]={:^20e}", i, htst_info.eigenvalues_sp[i], i+2*nos-10, htst_info.eigenvalues_sp[i+2*nos-10])); - Log.SendBlock(Utility::Log_Level::Info, Utility::Log_Sender::HTST, block, -1, -1); - - // Check if lowest eigenvalue < 0 (else it's not a SP) - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, "Checking if actually a saddle point..."); - if( htst_info.eigenvalues_sp[0] > -epsilon ) - { - Log(Utility::Log_Level::Error, Utility::Log_Sender::All, fmt::format( - "HTST: the transition configuration is not a saddle point, its lowest eigenvalue is above the threshold ({} > {})!", htst_info.eigenvalues_sp[0], -epsilon )); - return; - } - - // Check if second-lowest eigenvalue < 0 (higher-order SP) - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, "Checking if higher order saddle point..."); - int n_negative = 0; - for( int i=0; i < htst_info.eigenvalues_sp.size(); ++i ) - { - if( htst_info.eigenvalues_sp[i] < -epsilon ) - ++n_negative; - } - if( n_negative > 1 ) - { - Log(Utility::Log_Level::Error, Utility::Log_Sender::All, fmt::format( - "HTST: the image you passed is a higher order saddle point (N={})!", n_negative )); - return; - } - - // Perpendicular velocity - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, "Calculating perpendicular velocity at saddle point ('a' factors)..."); - // Calculation of the 'a' parameters... - VectorX perpendicular_velocity_sp(2*nos); - MatrixX basis_sp(3*nos, 2*nos); - Manifoldmath::tangent_basis_spherical(image_sp, basis_sp); - // Manifoldmath::tangent_basis(image_sp, basis_sp); - // Calculate_Perpendicular_Velocity_2N(image_sp, hessian_geodesic_sp_2N, basis_sp, eigenvectors_sp, perpendicular_velocity_sp); - Calculate_Perpendicular_Velocity(image_sp, htst_info.saddle_point->geometry->mu_s, hessian_geodesic_sp_3N, basis_sp, eigenvectors_sp, perpendicular_velocity_sp); - htst_info.perpendicular_velocity.assign(perpendicular_velocity_sp.data(), perpendicular_velocity_sp.data() + 2*nos); - } - // End saddle point - //////////////////////////////////////////////////////////////////////// - - // Checking for zero modes at the saddle point... - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, "Checking for zero modes at the saddle point..."); - int n_zero_modes_sp = 0; - for( int i=0; i < htst_info.eigenvalues_sp.size(); ++i ) - { - if( std::abs(htst_info.eigenvalues_sp[i] ) <= epsilon) - ++n_zero_modes_sp; - } - - // Deal with zero modes if any (calculate volume) - htst_info.volume_sp = 1; - if( n_zero_modes_sp > 0 ) - { - Log(Utility::Log_Level::All, Utility::Log_Sender::HTST, fmt::format("ZERO MODES AT SADDLE POINT (N={})", n_zero_modes_sp)); - - if( is_afm ) - htst_info.volume_sp = Calculate_Zero_Volume(htst_info.saddle_point); - else - htst_info.volume_sp = Calculate_Zero_Volume(htst_info.saddle_point); - } - - // Calculate "s" - htst_info.s = 0; - for( int i = n_zero_modes_sp+1; i < 2*nos; ++i ) - htst_info.s += std::pow(htst_info.perpendicular_velocity[i], 2) / htst_info.eigenvalues_sp[i]; - htst_info.s = std::sqrt(htst_info.s); - - //////////////////////////////////////////////////////////////////////// - // Initial state minimum - { - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, "Calculation for the Minimum"); - - // Evaluation of the Hessian... - MatrixX hessian_minimum = MatrixX::Zero(3*nos,3*nos); - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Evaluation of the Hessian..."); - htst_info.minimum->hamiltonian->Hessian(image_minimum, hessian_minimum); - - // Eigendecomposition - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Eigendecomposition..."); - MatrixX hessian_geodesic_minimum_3N(3*nos, 3*nos); - MatrixX hessian_geodesic_minimum_2N(2*nos, 2*nos); - VectorX eigenvalues_minimum = VectorX::Zero(2*nos); - MatrixX eigenvectors_minimum = MatrixX::Zero(2*nos, 2*nos); - Geodesic_Eigen_Decomposition(image_minimum, gradient_minimum, hessian_minimum, - hessian_geodesic_minimum_3N, hessian_geodesic_minimum_2N, eigenvalues_minimum, eigenvectors_minimum); - - htst_info.eigenvalues_min.assign(eigenvalues_minimum.data(), eigenvalues_minimum.data() + 2*nos); - // htst_info.eigenvectors_min.assign(eigenvectors_minimum.data(), eigenvectors_minimum.data() + 2*3*nos); - - // Print some eigenvalues - block = std::vector{"10 lowest eigenvalues at minimum:"}; - for( int i=0; i<10; ++i ) - block.push_back(fmt::format("ew[{}]={:^20e} ew[{}]={:^20e}", i, htst_info.eigenvalues_min[i], i+2*nos-10, htst_info.eigenvalues_min[i+2*nos-10])); - Log.SendBlock(Utility::Log_Level::Info, Utility::Log_Sender::HTST, block, -1, -1); - - // Check for eigenvalues < 0 (i.e. not a minimum) - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, "Checking if actually a minimum..."); - if( htst_info.eigenvalues_min[0] < -epsilon ) - { - Log(Utility::Log_Level::Error, Utility::Log_Sender::All, fmt::format( - "HTST: the initial configuration is not a minimum, its lowest eigenvalue is below the threshold ({} < {})!", htst_info.eigenvalues_min[0], -epsilon )); - return; - } - } - // End initial state minimum - //////////////////////////////////////////////////////////////////////// - - // Checking for zero modes at the minimum... - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, "Checking for zero modes at the minimum..."); - int n_zero_modes_minimum = 0; - for( int i=0; i < htst_info.eigenvalues_min.size(); ++i ) - { - if( std::abs(htst_info.eigenvalues_min[i]) <= epsilon ) - ++n_zero_modes_minimum; - } - - // Deal with zero modes if any (calculate volume) - htst_info.volume_min = 1; - if( n_zero_modes_minimum > 0 ) - { - Log(Utility::Log_Level::All, Utility::Log_Sender::HTST, fmt::format("ZERO MODES AT MINIMUM (N={})", n_zero_modes_minimum)); - - if( is_afm ) - htst_info.volume_min = Calculate_Zero_Volume(htst_info.minimum); - else - htst_info.volume_min = Calculate_Zero_Volume(htst_info.minimum); - } - - //////////////////////////////////////////////////////////////////////// - // Calculation of the prefactor... - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, "Calculating prefactor..."); - - // Calculate the exponent for the temperature-dependence of the prefactor - // The exponent depends on the number of zero modes at the different states - htst_info.temperature_exponent = 0.5 * (n_zero_modes_minimum - n_zero_modes_sp); - - // Calculate "me" - htst_info.me = std::pow(2*C::Pi * C::k_B, htst_info.temperature_exponent); - - // Calculate Omega_0, i.e. the entropy contribution - htst_info.Omega_0 = 1; - if( n_zero_modes_minimum > n_zero_modes_sp+1 ) - { - for( int i = n_zero_modes_sp+1; i system) - { - int nos = system->geometry->nos; - auto& n_cells = system->geometry->n_cells; - auto& spins = *system->spins; - auto& spin_positions = system->geometry->positions; - auto& geometry = *system->geometry; - auto& bravais_vectors = system->geometry->bravais_vectors; - - // Dimensionality of the zero mode - int zero_mode_dimensionality = 0; - Vector3 zero_mode_length{0,0,0}; - for( int ibasis=0; ibasis<3; ++ibasis ) - { - // Only a periodical direction can be a true zero mode - if( system->hamiltonian->boundary_conditions[ibasis] && geometry.n_cells[ibasis] > 1 ) - { - // Vector3 shift_pos, test_pos; - vectorfield spins_shifted(nos, Vector3{0,0,0}); - - int shift = 0; - if( ibasis == 0 ) - shift = geometry.n_cell_atoms; - else if( ibasis == 1 ) - shift = geometry.n_cell_atoms * geometry.n_cells[0]; - else if( ibasis == 2 ) - shift = geometry.n_cell_atoms * geometry.n_cells[0] * geometry.n_cells[1]; - - for (int isite = 0; isite < nos; ++isite) - { - spins_shifted[(isite + shift) % nos] = spins[isite]; - } - - zero_mode_length[ibasis] = Manifoldmath::dist_geodesic(spins, spins_shifted); - - // Increment zero mode dimensionality - ++zero_mode_dimensionality; - } - } - - // Calculate the volume depending on the number of periodical boundaries - scalar zero_volume = 1; - if( zero_mode_dimensionality == 1 ) - { - zero_volume = zero_mode_length[0]; - } - else if( zero_mode_dimensionality == 2 ) - { - scalar area_factor = ( bravais_vectors[0].normalized().cross(bravais_vectors[1].normalized()) ).norm(); - zero_volume = zero_mode_length[0]*zero_mode_length[1] * area_factor; - } - else if( zero_mode_dimensionality == 3 ) - { - scalar volume_factor = std::abs( (bravais_vectors[0].normalized().cross(bravais_vectors[1].normalized()) ).dot( - bravais_vectors[2].normalized()) ); - zero_volume = zero_mode_length[0]*zero_mode_length[1]*zero_mode_length[2] * volume_factor; - } - - Log.SendBlock(Utility::Log_Level::Info, Utility::Log_Sender::HTST, - { - fmt::format("ZV zero mode dimensionality = {}", zero_mode_dimensionality), - fmt::format("ZV zero mode length = {}", zero_mode_length.transpose()), - fmt::format("ZV = {}", zero_volume) - }, -1, -1); - - // Return - return zero_volume; - } - - - void Calculate_Perpendicular_Velocity(const vectorfield & spins, const scalarfield & mu_s, const MatrixX & hessian, - const MatrixX & basis, const MatrixX & eigenbasis, VectorX & perpendicular_velocity) - { - int nos = spins.size(); - - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Calculate_Perpendicular_Velocity: calculate velocity matrix"); - - // Calculate the velocity matrix in the 3N-basis - MatrixX velocity(3*nos, 3*nos); - Calculate_Dynamical_Matrix(spins, mu_s, hessian, velocity); - - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Calculate_Perpendicular_Velocity: project velocity matrix"); - - // Project the velocity matrix into the 2N tangent space - MatrixX velocity_projected(2*nos, 2*nos); - velocity_projected = basis.transpose()*velocity*basis; - - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Calculate_Perpendicular_Velocity: calculate a"); - - // The velocity components orthogonal to the dividing surface - perpendicular_velocity = eigenbasis.col(0).transpose() * (velocity_projected*eigenbasis); - - // std::cerr << " Calculate_Perpendicular_Velocity: sorting" << std::endl; - // std::sort(perpendicular_velocity.data(),perpendicular_velocity.data()+perpendicular_velocity.size()); - - std::vector block(0); - for (int i=0; i<10; ++i) - block.push_back(fmt::format(" a[{}] = {}", i, perpendicular_velocity[i])); - Log.SendBlock(Utility::Log_Level::Info, Utility::Log_Sender::HTST, block, -1, -1); - - // std::cerr << "without units:" << std::endl; - // for (int i=0; i<10; ++i) - // std::cerr << " a[" << i << "] = " << perpendicular_velocity[i]/C::mu_B << std::endl; - } - - - void Calculate_Dynamical_Matrix(const vectorfield & spins, const scalarfield & mu_s, const MatrixX & hessian, MatrixX & velocity) - { - velocity.setZero(); - int nos = spins.size(); - - for (int i=0; i < nos; ++i) - { - Vector3 beff{0, 0, 0}; - - for (int j=0; j < nos; ++j) - { - velocity(3*i, 3*j) = spins[i][1]*hessian(3*i+2,3*j) - spins[i][2]*hessian(3*i+1,3*j); - velocity(3*i, 3*j+1) = spins[i][1]*hessian(3*i+2,3*j+1) - spins[i][2]*hessian(3*i+1,3*j+1); - velocity(3*i, 3*j+2) = spins[i][1]*hessian(3*i+2,3*j+2) - spins[i][2]*hessian(3*i+1,3*j+2); - - velocity(3*i+1, 3*j) = spins[i][2]*hessian(3*i,3*j) - spins[i][0]*hessian(3*i+2,3*j); - velocity(3*i+1, 3*j+1) = spins[i][2]*hessian(3*i,3*j+1) - spins[i][0]*hessian(3*i+2,3*j+1); - velocity(3*i+1, 3*j+2) = spins[i][2]*hessian(3*i,3*j+2) - spins[i][0]*hessian(3*i+2,3*j+2); - - velocity(3*i+2, 3*j) = spins[i][0]*hessian(3*i+1,3*j) - spins[i][1]*hessian(3*i,3*j); - velocity(3*i+2, 3*j+1) = spins[i][0]*hessian(3*i+1,3*j+1) - spins[i][1]*hessian(3*i,3*j+1); - velocity(3*i+2, 3*j+2) = spins[i][0]*hessian(3*i+1,3*j+2) - spins[i][1]*hessian(3*i,3*j+2); - - beff -= hessian.block<3, 3>(3*i, 3*j) * spins[j]; - } - - velocity(3*i, 3*i+1) -= beff[2]; - velocity(3*i, 3*i+2) += beff[1]; - velocity(3*i+1, 3*i) += beff[2]; - velocity(3*i+1, 3*i+2) -= beff[0]; - velocity(3*i+2, 3*i) -= beff[1]; - velocity(3*i+2, 3*i+1) += beff[0]; - - velocity.row(3*i) /= mu_s[i]; - velocity.row(3*i+1) /= mu_s[i]; - velocity.row(3*i+2) /= mu_s[i]; - } - } - - - void hessian_bordered_3N(const vectorfield & image, const vectorfield & gradient, const MatrixX & hessian, MatrixX & 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(); - hessian_out = hessian; - - VectorX lambda(nos); - for (int i=0; i matrix_solver(matrix); - - evalues = matrix_solver.eigenvalues().real(); - evectors = matrix_solver.eigenvectors().real(); - } - - void Eigen_Decomposition_Spectra(int nos, const MatrixX & matrix, VectorX & evalues, MatrixX & evectors, int n_decompose=1) - { - int n_steps = std::max(2, nos); - - // Create a Spectra solver - Spectra::DenseGenMatProd op(matrix); - Spectra::GenEigsSolver< scalar, Spectra::SMALLEST_REAL, Spectra::DenseGenMatProd > matrix_spectrum(&op, n_decompose, n_steps); - matrix_spectrum.init(); - - // Compute the specified spectrum - int nconv = matrix_spectrum.compute(); - - if (matrix_spectrum.info() == Spectra::SUCCESSFUL) - { - evalues = matrix_spectrum.eigenvalues().real(); - evectors = matrix_spectrum.eigenvectors().real(); - // Eigen::Ref evec = evectors.col(0); - } - else - { - Log(Utility::Log_Level::Error, Utility::Log_Sender::All, "Failed to calculate eigenvectors of the Matrix!"); - evalues.setZero(); - evectors.setZero(); - } - } - - void Geodesic_Eigen_Decomposition(const vectorfield & image, const vectorfield & gradient, const MatrixX & hessian, - MatrixX & hessian_geodesic_3N, MatrixX & hessian_geodesic_2N, VectorX & eigenvalues, MatrixX & eigenvectors) - { - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, "---------- Geodesic Eigen Decomposition"); - - int nos = image.size(); - - // Calculate geodesic Hessian in 3N-representation - hessian_geodesic_3N = MatrixX::Zero(3*nos, 3*nos); - hessian_bordered_3N(image, gradient, hessian, hessian_geodesic_3N); - - // Transform into geodesic Hessian - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Transforming Hessian into geodesic Hessian..."); - hessian_geodesic_2N = MatrixX::Zero(2*nos, 2*nos); - // Manifoldmath::hessian_bordered(image, gradient, hessian, hessian_geodesic_2N); - // Manifoldmath::hessian_projected(image, gradient, hessian, hessian_geodesic_2N); - // Manifoldmath::hessian_weingarten(image, gradient, hessian, hessian_geodesic_2N); - // Manifoldmath::hessian_spherical(image, gradient, hessian, hessian_geodesic_2N); - // Manifoldmath::hessian_covariant(image, gradient, hessian, hessian_geodesic_2N); - - // Do this manually - MatrixX basis = MatrixX::Zero(3*nos, 2*nos); - Manifoldmath::tangent_basis_spherical(image, basis); - // Manifoldmath::tangent_basis(image, basis); - hessian_geodesic_2N = basis.transpose() * hessian_geodesic_3N * basis; - - - // Calculate full eigenspectrum - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, " Calculation of full eigenspectrum..." ); - // std::cerr << hessian_geodesic_2N.cols() << " " << hessian_geodesic_2N.rows() << std::endl; - eigenvalues = VectorX::Zero(2*nos); - eigenvectors = MatrixX::Zero(2*nos, 2*nos); - Eigen_Decomposition(hessian_geodesic_2N, eigenvalues, eigenvectors); - // Eigen_Decomposition_Spectra(hessian_geodesic_2N, eigenvalues, eigenvectors); - - Log(Utility::Log_Level::Info, Utility::Log_Sender::HTST, "---------- Geodesic Eigen Decomposition Done"); - } - - - }// end namespace HTST -}// end namespace Engine \ No newline at end of file