From 51396c63fadfe64ac382c4d94f3313888a2569da Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Wed, 17 Nov 2021 10:42:02 +0000 Subject: [PATCH 01/41] Method_EMA: WIP implementation for sparse matrices --- core/include/engine/Eigenmodes.hpp | 9 ++- core/include/engine/Manifoldmath.hpp | 3 + core/src/engine/Eigenmodes.cpp | 85 ++++++++++++++++++++++++---- core/src/engine/Manifoldmath.cpp | 31 ++++++++++ 4 files changed, 117 insertions(+), 11 deletions(-) diff --git a/core/include/engine/Eigenmodes.hpp b/core/include/engine/Eigenmodes.hpp index 5fbf0c570..cf033447f 100644 --- a/core/include/engine/Eigenmodes.hpp +++ b/core/include/engine/Eigenmodes.hpp @@ -33,7 +33,14 @@ bool Hessian_Full_Spectrum( bool Hessian_Partial_Spectrum( const std::shared_ptr parameters, const vectorfield & spins, const vectorfield & gradient, const MatrixX & hessian, int n_modes, MatrixX & tangent_basis, MatrixX & hessian_constrained, VectorX & eigenvalues, - MatrixX & eigenvectors ); + MatrixX & eigenvcetors ); + +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 #endif \ No newline at end of file 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/src/engine/Eigenmodes.cpp b/core/src/engine/Eigenmodes.cpp index 7b7ed8c09..5e081f1a3 100644 --- a/core/src/engine/Eigenmodes.cpp +++ b/core/src/engine/Eigenmodes.cpp @@ -4,6 +4,8 @@ // #include #include // Also includes +#include // Also includes + #include #include @@ -65,7 +67,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 ); @@ -75,18 +76,41 @@ void Calculate_Eigenmodes( std::shared_ptr system, int idx_im // g[idx] = mask[idx]*g[idx]; // }); 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); + + // TODO: do this properly + bool sparse = true; + 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 ) { @@ -216,5 +240,46 @@ 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) + + // Create the Spectra Matrix product operation + Spectra::SparseSymMatProd op( hessian_constrained ); + // Create and initialize a Spectra solver + Spectra::SymEigsSolver> hessian_spectrum( + &op, n_modes, 2 * nos ); + hessian_spectrum.init(); + + // Compute the specified spectrum, sorted by smallest real eigenvalue + int nconv = hessian_spectrum.compute( 1000, 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/Manifoldmath.cpp b/core/src/engine/Manifoldmath.cpp index e9c55f934..18447c7a3 100644 --- a/core/src/engine/Manifoldmath.cpp +++ b/core/src/engine/Manifoldmath.cpp @@ -469,6 +469,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 +529,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 ) From f87ec6badb9cbd63b5ac4a4604657e6182ccb1e6 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Wed, 17 Nov 2021 11:29:43 +0000 Subject: [PATCH 02/41] EMA: Can now control `sparse` parameter via API. Also added function to clear the modes. --- core/include/Spirit/Parameters_EMA.h | 9 +++ core/include/data/Parameters_Method_EMA.hpp | 1 + core/python/spirit/parameters/ema.py | 79 ++++++++++++++++----- core/src/Spirit/Parameters_EMA.cpp | 56 +++++++++++++++ core/src/engine/Eigenmodes.cpp | 26 ++++--- 5 files changed, 141 insertions(+), 30 deletions(-) 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/data/Parameters_Method_EMA.hpp b/core/include/data/Parameters_Method_EMA.hpp index 77612a914..f8bd907f6 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 = true; // ----------------- Output -------------- // Energy output settings diff --git a/core/python/spirit/parameters/ema.py b/core/python/spirit/parameters/ema.py index 6996d3ea9..f6d229fc8 100644 --- a/core/python/spirit/parameters/ema.py +++ b/core/python/spirit/parameters/ema.py @@ -9,41 +9,88 @@ 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 -_EMA_Set_N_Modes = _spirit.Parameters_EMA_Set_N_Modes + +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, - ctypes.c_int, ctypes.c_int] -_EMA_Set_N_Modes.restype = None + ctypes.c_int, ctypes.c_int] +_EMA_Set_N_Modes.restype = None + + def set_n_modes(p_state, n_modes, idx_image=-1, idx_chain=-1): """Set the number of modes to calculate or use.""" _EMA_Set_N_Modes(ctypes.c_void_p(p_state), ctypes.c_int(n_modes), - ctypes.c_int(idx_image), ctypes.c_int(idx_chain)) + ctypes.c_int(idx_image), ctypes.c_int(idx_chain)) -_EMA_Set_N_Mode_Follow = _spirit.Parameters_EMA_Set_N_Mode_Follow + +_EMA_Set_N_Mode_Follow = _spirit.Parameters_EMA_Set_N_Mode_Follow _EMA_Set_N_Mode_Follow.argtypes = [ctypes.c_void_p, ctypes.c_int, - ctypes.c_int, ctypes.c_int] -_EMA_Set_N_Mode_Follow.restype = None + ctypes.c_int, ctypes.c_int] +_EMA_Set_N_Mode_Follow.restype = None + + def set_n_mode_follow(p_state, n_mode, idx_image=-1, idx_chain=-1): """Set the index of the mode to use.""" _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)) + ctypes.c_int(idx_image), ctypes.c_int(idx_chain)) + + +_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 -## ---------------------------------- Get ---------------------------------- -_EMA_Get_N_Modes = _spirit.Parameters_EMA_Get_N_Modes +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] -_EMA_Get_N_Modes.restype = ctypes.c_int +_EMA_Get_N_Modes.restype = ctypes.c_int + + def get_n_modes(p_state, idx_image=-1, idx_chain=-1): """Returns the number of modes to calculate or use.""" return int(_EMA_Get_N_Modes(p_state, ctypes.c_int(idx_image), ctypes.c_int(idx_chain))) -_EMA_Get_N_Mode_Follow = _spirit.Parameters_EMA_Get_N_Mode_Follow + +_EMA_Get_N_Mode_Follow = _spirit.Parameters_EMA_Get_N_Mode_Follow _EMA_Get_N_Mode_Follow.argtypes = [ctypes.c_void_p, ctypes.c_int, ctypes.c_int] -_EMA_Get_N_Mode_Follow.restype = ctypes.c_int +_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/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/engine/Eigenmodes.cpp b/core/src/engine/Eigenmodes.cpp index 5e081f1a3..b0017cbbd 100644 --- a/core/src/engine/Eigenmodes.cpp +++ b/core/src/engine/Eigenmodes.cpp @@ -3,8 +3,8 @@ #include // #include -#include // Also includes #include // Also includes +#include // Also includes #include #include @@ -76,16 +76,14 @@ void Calculate_Eigenmodes( std::shared_ptr system, int idx_im // g[idx] = mask[idx]*g[idx]; // }); Vectormath::set_c_a( 1, gradient, gradient, system->geometry->mask_unpinned ); - VectorX eigenvalues; MatrixX eigenvectors; - SpMatrixX tangent_basis = SpMatrixX(3*nos, 2*nos); + SpMatrixX tangent_basis = SpMatrixX( 3 * nos, 2 * nos ); - // TODO: do this properly - bool sparse = true; + bool sparse = system->ema_parameters->sparse; bool successful; - if(sparse) + if( sparse ) { // The Hessian (unprojected) SpMatrixX hessian( 3 * nos, 3 * nos ); @@ -96,14 +94,15 @@ void Calculate_Eigenmodes( std::shared_ptr system, int idx_im successful = Eigenmodes::Sparse_Hessian_Partial_Spectrum( system->ema_parameters, spins_initial, gradient, hessian, n_modes, tangent_basis, hessian_constrained, eigenvalues, eigenvectors ); - - } else { + } + 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); + MatrixX _tangent_basis = MatrixX( tangent_basis ); successful = Eigenmodes::Hessian_Partial_Spectrum( system->ema_parameters, spins_initial, gradient, hessian, n_modes, _tangent_basis, hessian_constrained, @@ -242,8 +241,8 @@ bool Hessian_Partial_Spectrum( 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 ) + const SpMatrixX & hessian, int n_modes, SpMatrixX & tangent_basis, SpMatrixX & hessian_constrained, + VectorX & eigenvalues, MatrixX & eigenvectors ) { int nos = spins.size(); @@ -251,9 +250,9 @@ bool Sparse_Hessian_Partial_Spectrum( 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); + Manifoldmath::sparse_tangent_basis_spherical( spins, tangent_basis ); - SpMatrixX hessian_constrained_3N = SpMatrixX(3*nos, 3*nos); + 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; @@ -280,6 +279,5 @@ bool Sparse_Hessian_Partial_Spectrum( return ( hessian_spectrum.info() == Spectra::SUCCESSFUL ) && ( nconv > 0 ); } - } // namespace Eigenmodes } // namespace Engine From 3c8a7e087126843cc4eabc4c8d767ebb790c304d Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Thu, 18 Nov 2021 16:23:10 +0000 Subject: [PATCH 03/41] Eigenmodes: Use a better value of ncv for Sparse_Hessian_Partial_Spectrum --- core/src/engine/Eigenmodes.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/core/src/engine/Eigenmodes.cpp b/core/src/engine/Eigenmodes.cpp index b0017cbbd..757f28053 100644 --- a/core/src/engine/Eigenmodes.cpp +++ b/core/src/engine/Eigenmodes.cpp @@ -259,15 +259,19 @@ bool Sparse_Hessian_Partial_Spectrum( // 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, 2 * nos ); + &op, n_modes, ncv); hessian_spectrum.init(); // Compute the specified spectrum, sorted by smallest real eigenvalue - int nconv = hessian_spectrum.compute( 1000, 1e-10, int( Spectra::SMALLEST_ALGE ) ); + int nconv = hessian_spectrum.compute( max_iter, 1e-10, int( Spectra::SMALLEST_ALGE ) ); // Extract real eigenvalues eigenvalues = hessian_spectrum.eigenvalues().real(); From 9b998cb628e0bbb33b030bf8849974965a909786 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Thu, 25 Nov 2021 19:42:19 +0000 Subject: [PATCH 04/41] Fixed an issue when reading in eigenmodes from the GUI --- core/src/Spirit/IO.cpp | 2 +- ui-cpp/ui-qt/src/ParametersWidget.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) 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/ui-cpp/ui-qt/src/ParametersWidget.cpp b/ui-cpp/ui-qt/src/ParametersWidget.cpp index 9c8843698..9495f01b3 100644 --- a/ui-cpp/ui-qt/src/ParametersWidget.cpp +++ b/ui-cpp/ui-qt/src/ParametersWidget.cpp @@ -580,7 +580,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() ) ); From ece216131bccf0f906a3eb86eb4713756cd7ec3a Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Fri, 26 Nov 2021 11:28:06 +0000 Subject: [PATCH 05/41] Rough implementation of moving endpoints gneb --- core/include/data/Parameters_Method_GNEB.hpp | 2 + core/src/engine/Method_GNEB.cpp | 40 +++++++++++++++++++- 2 files changed, 41 insertions(+), 1 deletion(-) diff --git a/core/include/data/Parameters_Method_GNEB.hpp b/core/include/data/Parameters_Method_GNEB.hpp index b80a06883..c98cb14f2 100644 --- a/core/include/data/Parameters_Method_GNEB.hpp +++ b/core/include/data/Parameters_Method_GNEB.hpp @@ -33,6 +33,8 @@ struct Parameters_Method_GNEB : Parameters_Method_Solver // Mersenne twister PRNG std::mt19937 prng = std::mt19937( rng_seed ); + bool moving_endpoints = false; + // ----------------- Output -------------- bool output_energies_step = false; bool output_energies_divide_by_nspins = true; diff --git a/core/src/engine/Method_GNEB.cpp b/core/src/engine/Method_GNEB.cpp index edac09353..5f7b6dd2c 100644 --- a/core/src/engine/Method_GNEB.cpp +++ b/core/src/engine/Method_GNEB.cpp @@ -252,6 +252,44 @@ 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 ) + if( true ) + { + for(int img : {0,chain->noi-1}) + { + int nos = chain->images[img]->nos; + scalar _Rx0 = 1e-3 * nos; + + auto _F_total = F_total[img].data(); + auto _forces = forces[img].data(); + auto _F_gradient = F_gradient[img].data(); + auto _tangents = tangents[img].data(); + + scalar _Rx = (img==0 ? Rx[1]-Rx[0] : Rx[chain->noi-1] - Rx[chain->noi-2]); + // scalar _Rx = (img==0 ? lengths[1] : lengths[chain->noi-1]); + + auto dE_dRx = (img==0 ? energies[1] - energies[0] : energies[chain->noi-1] - energies[chain->noi - 2]) / _Rx; + + auto spring_constant = dE_dRx * this->chain->gneb_parameters->spring_constant; + + std::cout << "image" << img << "\n"; + std::cout << _tangents[0].transpose() << "\n"; + std::cout << _Rx << "\n"; + std::cout << _F_gradient[0].transpose() << "\n--\n"; + + + Backend::par::apply( + nos, + [_forces, _F_gradient, _Rx, _Rx0, _tangents, spring_constant ] SPIRIT_LAMBDA( int idx ) + { + _forces[idx] = _F_gradient[idx] - _F_gradient[idx].dot(_tangents[idx]) * _tangents[idx] + 1 * spring_constant * ( _Rx - _Rx0 ) * _tangents[idx]; + } + ); + } + } + } // end Calculate template @@ -262,7 +300,7 @@ void Method_GNEB::Calculate_Force_Virtual( using namespace Utility; // Calculate the cross product with the spin configuration to get direct minimization - for( unsigned int i = 1; i < configurations.size() - 1; ++i ) + for( unsigned int i = 0; i < configurations.size(); ++i ) { auto & image = *configurations[i]; auto & force = forces[i]; From 8bb8a71967c087d77965c540aa3195da3a97c815 Mon Sep 17 00:00:00 2001 From: "G. P. Mueller" Date: Thu, 9 Dec 2021 17:50:38 +0100 Subject: [PATCH 06/41] Core: fix clang-tidy linting in Exception.hpp --- core/include/utility/Exception.hpp | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/core/include/utility/Exception.hpp b/core/include/utility/Exception.hpp index b59c4acd6..d1f70f713 100644 --- a/core/include/utility/Exception.hpp +++ b/core/include/utility/Exception.hpp @@ -81,22 +81,24 @@ void Handle_Exception_API( void Handle_Exception_Core( const std::string & message, const char * file, unsigned int line, const char * function ); // 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 From 06d5bc6499af83e7298588ceb1326534ad6484c8 Mon Sep 17 00:00:00 2001 From: Moritz Date: Tue, 14 Dec 2021 15:15:08 +0100 Subject: [PATCH 07/41] Quantities: Get_Magnetization now respects mu_s Previously, the values of mu_s were ignored, so that the computed quantity was only the average of the spin directions. Unfortunately computing the magnetization now requires a temporary field of Vector3 (especially when using CUDA), which has negative performance implications. Therefore the new function `Quantity_Get_Average_Spin` has been introduced which is now used to display the average of the spin direction in the QT GUI. Also the magnetization has been removed from the `history` of Method_LLG. Some unit tests had to be adjusted, because they were relying on the old (wrong) behaviour of Quantity_Get_Magnetization --- core/include/Spirit/Quantities.h | 3 +++ core/include/engine/Vectormath.hpp | 2 +- core/python/test/quantities.py | 6 ++++-- core/src/Spirit/Quantities.cpp | 23 ++++++++++++++++++++++- core/src/engine/Method_LLG.cpp | 5 +++-- core/src/engine/Vectormath.cpp | 6 ++++-- core/test/test_solvers.cpp | 4 ++-- core/test/test_vectormath.cpp | 6 ++++-- ui-cpp/ui-qt/src/MainWindow.cpp | 4 ++-- 9 files changed, 45 insertions(+), 14 deletions(-) diff --git a/core/include/Spirit/Quantities.h b/core/include/Spirit/Quantities.h index 6826d49f8..9252a6fb3 100644 --- a/core/include/Spirit/Quantities.h +++ b/core/include/Spirit/Quantities.h @@ -14,6 +14,9 @@ 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 ); diff --git a/core/include/engine/Vectormath.hpp b/core/include/engine/Vectormath.hpp index e58979dc9..cf45228f0 100644 --- a/core/include/engine/Vectormath.hpp +++ b/core/include/engine/Vectormath.hpp @@ -531,7 +531,7 @@ 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 inside a vectorfield scalar TopologicalCharge( const vectorfield & vf, const Data::Geometry & geom, const intfield & boundary_conditions ); diff --git a/core/python/test/quantities.py b/core/python/test/quantities.py index 0eff6bfe6..9f09de090 100644 --- a/core/python/test/quantities.py +++ b/core/python/test/quantities.py @@ -5,7 +5,7 @@ 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 @@ -25,10 +25,12 @@ 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) ######### diff --git a/core/src/Spirit/Quantities.cpp b/core/src/Spirit/Quantities.cpp index 8642e6f53..9bd0fd780 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(); diff --git a/core/src/engine/Method_LLG.cpp b/core/src/engine/Method_LLG.cpp index a3a459ebb..3f217f0f9 100644 --- a/core/src/engine/Method_LLG.cpp +++ b/core/src/engine/Method_LLG.cpp @@ -307,8 +307,9 @@ void Method_LLG::Save_Current( std::string starttime, int iteration, boo 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] ); + // 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 ) diff --git a/core/src/engine/Vectormath.cpp b/core/src/engine/Vectormath.cpp index f81094ad4..374e1201f 100644 --- a/core/src/engine/Vectormath.cpp +++ b/core/src/engine/Vectormath.cpp @@ -492,9 +492,11 @@ 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; } 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-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; From 91de35141cf6f22dea8c8c6f0aacd8d4ed9f655c Mon Sep 17 00:00:00 2001 From: Moritz Date: Tue, 14 Dec 2021 16:42:43 +0100 Subject: [PATCH 08/41] Quantities: Implemented a new function Quantity_Get_Toplogical_Charge_Density Gives the spatial distribution of the topological charge, as well as the indices of the used triangulation --- core/include/Spirit/Quantities.h | 4 +++ core/include/engine/Vectormath.hpp | 6 ++++ core/python/spirit/quantities.py | 28 ++++++++++++++++++- core/src/Spirit/Quantities.cpp | 45 ++++++++++++++++++++++++++++++ core/src/engine/Vectormath.cpp | 27 +++++++++++++++--- ui-cpp/ui-qt/src/InfoWidget.cpp | 8 +++--- 6 files changed, 109 insertions(+), 9 deletions(-) diff --git a/core/include/Spirit/Quantities.h b/core/include/Spirit/Quantities.h index 9252a6fb3..21e14688c 100644 --- a/core/include/Spirit/Quantities.h +++ b/core/include/Spirit/Quantities.h @@ -23,6 +23,10 @@ PREFIX void Quantity_Get_Magnetization( State * state, float m[3], int idx_image // 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/engine/Vectormath.hpp b/core/include/engine/Vectormath.hpp index cf45228f0..0a5b403a8 100644 --- a/core/include/engine/Vectormath.hpp +++ b/core/include/engine/Vectormath.hpp @@ -532,6 +532,12 @@ inline int idx_from_pair( // Calculate the mean of a vectorfield 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/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/src/Spirit/Quantities.cpp b/core/src/Spirit/Quantities.cpp index 9bd0fd780..b621104e7 100644 --- a/core/src/Spirit/Quantities.cpp +++ b/core/src/Spirit/Quantities.cpp @@ -86,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/engine/Vectormath.cpp b/core/src/engine/Vectormath.cpp index 374e1201f..623b2ee2a 100644 --- a/core/src/engine/Vectormath.cpp +++ b/core/src/engine/Vectormath.cpp @@ -501,9 +501,12 @@ std::array Magnetization( const vectorfield & vf, const scalarfield & 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 @@ -569,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] ); @@ -602,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/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() ); From 8c554edc1fb32ab8491683b41f48ff6964684f07 Mon Sep 17 00:00:00 2001 From: Moritz Date: Wed, 15 Dec 2021 12:47:27 +0100 Subject: [PATCH 09/41] Quantities: Added topological charge to unit test Also changed the `api.cfg` input file in the quantities unit test to be more general and use an additional basis atom instead of a monoatomic basis. The quantities unit test now uses the `api.cfg` input file instead of `fd_neighbours.cfg`, because the system size should be actually be able to fit an entire skyrmion. --- core/python/test/quantities.py | 12 ++++++++++-- core/test/input/api.cfg | 8 ++++++-- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/core/python/test/quantities.py b/core/python/test/quantities.py index 9f09de090..3bc3c6691 100644 --- a/core/python/test/quantities.py +++ b/core/python/test/quantities.py @@ -11,7 +11,7 @@ ########## -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 @@ -31,7 +31,15 @@ def test_magnetization(self): self.assertAlmostEqual(M[0], 0) self.assertAlmostEqual(M[1], 0) 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/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 From 0706699ecfb9830993e87b337aca912c8b5625ff Mon Sep 17 00:00:00 2001 From: Moritz Date: Thu, 16 Dec 2021 11:50:28 +0100 Subject: [PATCH 10/41] Configparser: Fixed a small error where the log of individual pinned sites would just repeat the direction of the first fixed spin. --- core/src/io/Configparser.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/core/src/io/Configparser.cpp b/core/src/io/Configparser.cpp index 3793e1562..b09ad94ba 100644 --- a/core/src/io/Configparser.cpp +++ b/core/src/io/Configparser.cpp @@ -528,7 +528,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() ) ); } } } From 41403b4e1148172bb642931ba585a912774a2dce Mon Sep 17 00:00:00 2001 From: Moritz Date: Thu, 16 Dec 2021 12:03:10 +0100 Subject: [PATCH 11/41] Hamiltonian_Heisenberg: Fixed race condition in DDI In the pointwise multiplication the parallelization was collapsing one loop too many, and thus parallelized the summation of the inter-sublattice contributions. This lead to a racecondition when n_cell_atoms > 1. The effects of this racecondition were only really seen when using a very high number of threads compared to the system size. --- core/src/engine/Hamiltonian_Heisenberg.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/core/src/engine/Hamiltonian_Heisenberg.cpp b/core/src/engine/Hamiltonian_Heisenberg.cpp index 9d99f1603..d84a9f0ce 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 ) From 7ec8ee469dec0069c9d557ff95c4e16023840ab3 Mon Sep 17 00:00:00 2001 From: Moritz Date: Thu, 16 Dec 2021 15:46:34 +0100 Subject: [PATCH 12/41] API: implemented force and energy history in Simulation_Run_Info --- core/include/Spirit/Simulation.h | 9 ++++ core/include/engine/Method.hpp | 8 +-- core/python/spirit/simulation.py | 20 +++++++- core/src/Spirit/Simulation.cpp | 30 +++++++++++ core/src/engine/Method.cpp | 5 +- core/src/engine/Method_GNEB.cpp | 6 ++- core/src/engine/Method_LLG.cpp | 85 +++++++++++++++----------------- core/src/engine/Method_MC.cpp | 6 +-- core/src/engine/Method_MMF.cpp | 5 +- 9 files changed, 117 insertions(+), 57 deletions(-) diff --git a/core/include/Spirit/Simulation.h b/core/include/Spirit/Simulation.h index a62e6d7de..0a807be03 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 = NULL; + int n_history_max_torque = 0; + float * history_max_torque = NULL; + int n_history_energy = 0; + float * history_energy = NULL; }; +PREFIX void free_run_info(Simulation_Run_Info info) SUFFIX; + /* Start or stop a simulation -------------------------------------------------------------------- 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/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/src/Spirit/Simulation.cpp b/core/src/Spirit/Simulation.cpp index 23bbdc67f..18012e279 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); + } } } } 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..d416197b0 100644 --- a/core/src/engine/Method_GNEB.cpp +++ b/core/src/engine/Method_GNEB.cpp @@ -55,7 +55,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(); @@ -491,7 +491,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 ) diff --git a/core/src/engine/Method_LLG.cpp b/core/src/engine/Method_LLG.cpp index 3f217f0f9..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,9 +299,13 @@ 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 ); + 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] ); @@ -335,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 ) From 6b0ae0bbffbe9dde28ceb4d3583315f10f0b294d Mon Sep 17 00:00:00 2001 From: Moritz Date: Sat, 12 Feb 2022 11:48:23 +0100 Subject: [PATCH 13/41] Python Package: Fixed name of shared library in `setup.py` for Windows. Was erronously looking for `libSpirit.dll` instead of `Spirit.dll`. --- core/python/setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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', From c95766c4e0f34a2b3c66eac77c4e7cda29a95042 Mon Sep 17 00:00:00 2001 From: Moritz Date: Sat, 12 Feb 2022 11:52:02 +0100 Subject: [PATCH 14/41] C-API: Changed `NULL` to `nullptr` in `Simulation.h` --- core/include/Spirit/Simulation.h | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/core/include/Spirit/Simulation.h b/core/include/Spirit/Simulation.h index 0a807be03..bcfd84bb8 100644 --- a/core/include/Spirit/Simulation.h +++ b/core/include/Spirit/Simulation.h @@ -63,11 +63,11 @@ struct Simulation_Run_Info float max_torque = 0; int n_history_iteration = 0; - int * history_iteration = NULL; + int * history_iteration = nullptr; int n_history_max_torque = 0; - float * history_max_torque = NULL; + float * history_max_torque = nullptr; int n_history_energy = 0; - float * history_energy = NULL; + float * history_energy = nullptr; }; PREFIX void free_run_info(Simulation_Run_Info info) SUFFIX; @@ -80,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 From a30dc84170699461086ef10565e6301b99daf8a8 Mon Sep 17 00:00:00 2001 From: Moritz Date: Sat, 12 Feb 2022 12:35:00 +0100 Subject: [PATCH 15/41] GNEB: Added get and set functions for moving endpoints to API --- core/include/Spirit/Parameters_GNEB.h | 6 ++++ core/python/spirit/parameters/gneb.py | 15 +++++++++ core/src/Spirit/Parameters_GNEB.cpp | 43 ++++++++++++++++++++++++ core/src/engine/Method_GNEB.cpp | 47 ++++++++++++--------------- 4 files changed, 84 insertions(+), 27 deletions(-) diff --git a/core/include/Spirit/Parameters_GNEB.h b/core/include/Spirit/Parameters_GNEB.h index 9b930e7cc..2dc7e8cbf 100644 --- a/core/include/Spirit/Parameters_GNEB.h +++ b/core/include/Spirit/Parameters_GNEB.h @@ -111,6 +111,9 @@ 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 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 +172,9 @@ 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; + /* Returns the integer of whether an image is regular, climbing, falling, or stationary. diff --git a/core/python/spirit/parameters/gneb.py b/core/python/spirit/parameters/gneb.py index bbe549498..0d6ff4612 100644 --- a/core/python/spirit/parameters/gneb.py +++ b/core/python/spirit/parameters/gneb.py @@ -135,6 +135,13 @@ 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_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 +205,14 @@ 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_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/src/Spirit/Parameters_GNEB.cpp b/core/src/Spirit/Parameters_GNEB.cpp index d2e22118b..dffd7ecaa 100644 --- a/core/src/Spirit/Parameters_GNEB.cpp +++ b/core/src/Spirit/Parameters_GNEB.cpp @@ -237,6 +237,30 @@ 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 ); +} + void Parameters_GNEB_Set_Climbing_Falling( State * state, int image_type, int idx_image, int idx_chain ) noexcept try { @@ -507,6 +531,25 @@ 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; +} + int Parameters_GNEB_Get_Climbing_Falling( State * state, int idx_image, int idx_chain ) noexcept try { diff --git a/core/src/engine/Method_GNEB.cpp b/core/src/engine/Method_GNEB.cpp index 1ce01f4be..39f25a2f7 100644 --- a/core/src/engine/Method_GNEB.cpp +++ b/core/src/engine/Method_GNEB.cpp @@ -100,13 +100,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 ) { @@ -253,23 +250,23 @@ void Method_GNEB::Calculate_Force( } // end for img=1..noi-1 // Moving endpoints - // if( chain->gneb_parameters->moving_endpoints ) - if( true ) + if( chain->gneb_parameters->moving_endpoints ) { - for(int img : {0,chain->noi-1}) + for( int img : { 0, chain->noi - 1 } ) { - int nos = chain->images[img]->nos; + int nos = chain->images[img]->nos; scalar _Rx0 = 1e-3 * nos; - auto _F_total = F_total[img].data(); - auto _forces = forces[img].data(); + auto _F_total = F_total[img].data(); + auto _forces = forces[img].data(); auto _F_gradient = F_gradient[img].data(); - auto _tangents = tangents[img].data(); + auto _tangents = tangents[img].data(); - scalar _Rx = (img==0 ? Rx[1]-Rx[0] : Rx[chain->noi-1] - Rx[chain->noi-2]); + scalar _Rx = ( img == 0 ? Rx[1] - Rx[0] : Rx[chain->noi - 1] - Rx[chain->noi - 2] ); // scalar _Rx = (img==0 ? lengths[1] : lengths[chain->noi-1]); - auto dE_dRx = (img==0 ? energies[1] - energies[0] : energies[chain->noi-1] - energies[chain->noi - 2]) / _Rx; + auto dE_dRx + = ( img == 0 ? energies[1] - energies[0] : energies[chain->noi - 1] - energies[chain->noi - 2] ) / _Rx; auto spring_constant = dE_dRx * this->chain->gneb_parameters->spring_constant; @@ -278,14 +275,11 @@ void Method_GNEB::Calculate_Force( std::cout << _Rx << "\n"; std::cout << _F_gradient[0].transpose() << "\n--\n"; - Backend::par::apply( - nos, - [_forces, _F_gradient, _Rx, _Rx0, _tangents, spring_constant ] SPIRIT_LAMBDA( int idx ) - { - _forces[idx] = _F_gradient[idx] - _F_gradient[idx].dot(_tangents[idx]) * _tangents[idx] + 1 * spring_constant * ( _Rx - _Rx0 ) * _tangents[idx]; - } - ); + nos, [_forces, _F_gradient, _Rx, _Rx0, _tangents, spring_constant] SPIRIT_LAMBDA( int idx ) { + _forces[idx] = _F_gradient[idx] - _F_gradient[idx].dot( _tangents[idx] ) * _tangents[idx] + + 1 * spring_constant * ( _Rx - _Rx0 ) * _tangents[idx]; + } ); } } @@ -554,8 +548,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 @@ -599,8 +593,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; From 20981f364b25cfcf033d1b724daa165a322eac1d Mon Sep 17 00:00:00 2001 From: Moritz Date: Mon, 21 Feb 2022 10:11:08 +0100 Subject: [PATCH 16/41] Method_GNEB: fixed moving endpoints implementation --- core/src/engine/Method_GNEB.cpp | 24 +++++++----------------- 1 file changed, 7 insertions(+), 17 deletions(-) diff --git a/core/src/engine/Method_GNEB.cpp b/core/src/engine/Method_GNEB.cpp index 39f25a2f7..e87b0f432 100644 --- a/core/src/engine/Method_GNEB.cpp +++ b/core/src/engine/Method_GNEB.cpp @@ -254,31 +254,21 @@ void Method_GNEB::Calculate_Force( { for( int img : { 0, chain->noi - 1 } ) { - int nos = chain->images[img]->nos; - scalar _Rx0 = 1e-3 * nos; + scalar delta_Rx0 = 0.1; auto _F_total = F_total[img].data(); auto _forces = forces[img].data(); auto _F_gradient = F_gradient[img].data(); auto _tangents = tangents[img].data(); - scalar _Rx = ( img == 0 ? Rx[1] - Rx[0] : Rx[chain->noi - 1] - Rx[chain->noi - 2] ); - // scalar _Rx = (img==0 ? lengths[1] : lengths[chain->noi-1]); - - auto dE_dRx - = ( img == 0 ? energies[1] - energies[0] : energies[chain->noi - 1] - energies[chain->noi - 2] ) / _Rx; - - auto spring_constant = dE_dRx * this->chain->gneb_parameters->spring_constant; - - std::cout << "image" << img << "\n"; - std::cout << _tangents[0].transpose() << "\n"; - std::cout << _Rx << "\n"; - std::cout << _F_gradient[0].transpose() << "\n--\n"; + 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] ); Backend::par::apply( - nos, [_forces, _F_gradient, _Rx, _Rx0, _tangents, spring_constant] SPIRIT_LAMBDA( int idx ) { - _forces[idx] = _F_gradient[idx] - _F_gradient[idx].dot( _tangents[idx] ) * _tangents[idx] - + 1 * spring_constant * ( _Rx - _Rx0 ) * _tangents[idx]; + nos, [_forces, _F_total, _F_gradient, delta_Rx, delta_Rx0, _tangents, spring_constant, projection] SPIRIT_LAMBDA( int idx ) { + _forces[idx] = _F_gradient[idx] - projection * _tangents[idx] + spring_constant * ( delta_Rx - delta_Rx0 ) * _tangents[idx]; + _F_total[idx] = _forces[idx]; } ); } } From cfa686163060022875b1947a24e75c72d6533aef Mon Sep 17 00:00:00 2001 From: Moritz Date: Mon, 21 Feb 2022 11:14:24 +0100 Subject: [PATCH 17/41] Configparser: Added gneb_moving_endpoints to --- core/src/io/Configparser.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/core/src/io/Configparser.cpp b/core/src/io/Configparser.cpp index b09ad94ba..8cce9583c 100644 --- a/core/src/io/Configparser.cpp +++ b/core/src/io/Configparser.cpp @@ -1052,6 +1052,7 @@ 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" ); } catch( ... ) { @@ -1074,6 +1075,7 @@ std::unique_ptr Parameters_Method_GNEB_from_Config 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} = {}", "moving_endpoints", parameters->moving_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 ) ); From 614d334628c2eff0638511658855015ee68bf14e Mon Sep 17 00:00:00 2001 From: Moritz Date: Mon, 21 Feb 2022 11:16:21 +0100 Subject: [PATCH 18/41] GNEB: more fixes for moving endpoints --- core/src/engine/Method_GNEB.cpp | 30 +++++++++++++++++++++++------- 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/core/src/engine/Method_GNEB.cpp b/core/src/engine/Method_GNEB.cpp index e87b0f432..c383ee939 100644 --- a/core/src/engine/Method_GNEB.cpp +++ b/core/src/engine/Method_GNEB.cpp @@ -256,18 +256,28 @@ void Method_GNEB::Calculate_Force( { scalar delta_Rx0 = 0.1; + auto & image = *configurations[img]; + Manifoldmath::project_tangential( F_gradient[img], image ); + auto _F_total = F_total[img].data(); auto _forces = forces[img].data(); auto _F_gradient = F_gradient[img].data(); auto _tangents = tangents[img].data(); - 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] ); + 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] ); + + // std::cout << " === " << img << " ===\n"; + // fmt::print( "{}\n", _F_total[0].transpose() ); + // fmt::print( "{}\n", _F_gradient[0].transpose() ); + // fmt::print( "{}\n", delta_Rx ); Backend::par::apply( - nos, [_forces, _F_total, _F_gradient, delta_Rx, delta_Rx0, _tangents, spring_constant, projection] SPIRIT_LAMBDA( int idx ) { - _forces[idx] = _F_gradient[idx] - projection * _tangents[idx] + spring_constant * ( delta_Rx - delta_Rx0 ) * _tangents[idx]; + nos, [_forces, _F_total, _F_gradient, delta_Rx, delta_Rx0, _tangents, spring_constant, + projection] SPIRIT_LAMBDA( int idx ) { + _forces[idx] = _F_gradient[idx] - projection * _tangents[idx] + + spring_constant * ( delta_Rx - delta_Rx0 ) * _tangents[idx]; _F_total[idx] = _forces[idx]; } ); } @@ -283,8 +293,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]; @@ -327,7 +343,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 From c7e84313c24a13655f29677bb056ea011fa07fc7 Mon Sep 17 00:00:00 2001 From: Moritz Date: Wed, 2 Mar 2022 22:41:42 +0100 Subject: [PATCH 19/41] GNEB: equilibrium delta_Rx for moving endpoints can now be controlled via the API --- core/include/Spirit/Parameters_GNEB.h | 3 +++ core/include/data/Parameters_Method_GNEB.hpp | 2 ++ core/python/spirit/parameters/gneb.py | 7 ++++++ core/src/Spirit/Parameters_GNEB.cpp | 26 ++++++++++++++++++++ core/src/engine/Method_GNEB.cpp | 4 +-- 5 files changed, 40 insertions(+), 2 deletions(-) diff --git a/core/include/Spirit/Parameters_GNEB.h b/core/include/Spirit/Parameters_GNEB.h index 2dc7e8cbf..3bd812d83 100644 --- a/core/include/Spirit/Parameters_GNEB.h +++ b/core/include/Spirit/Parameters_GNEB.h @@ -114,6 +114,9 @@ PREFIX void Parameters_GNEB_Set_Path_Shortening_Constant( // 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 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; diff --git a/core/include/data/Parameters_Method_GNEB.hpp b/core/include/data/Parameters_Method_GNEB.hpp index c98cb14f2..213a42aa5 100644 --- a/core/include/data/Parameters_Method_GNEB.hpp +++ b/core/include/data/Parameters_Method_GNEB.hpp @@ -34,6 +34,8 @@ struct Parameters_Method_GNEB : Parameters_Method_Solver std::mt19937 prng = std::mt19937( rng_seed ); bool moving_endpoints = false; + scalar equilibrium_delta_Rx_left = 1.0; + scalar equilibrium_delta_Rx_right = 1.0; // ----------------- Output -------------- bool output_energies_step = false; diff --git a/core/python/spirit/parameters/gneb.py b/core/python/spirit/parameters/gneb.py index 0d6ff4612..5de4897a8 100644 --- a/core/python/spirit/parameters/gneb.py +++ b/core/python/spirit/parameters/gneb.py @@ -142,6 +142,13 @@ 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_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 diff --git a/core/src/Spirit/Parameters_GNEB.cpp b/core/src/Spirit/Parameters_GNEB.cpp index dffd7ecaa..896bc8817 100644 --- a/core/src/Spirit/Parameters_GNEB.cpp +++ b/core/src/Spirit/Parameters_GNEB.cpp @@ -261,6 +261,32 @@ 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 { diff --git a/core/src/engine/Method_GNEB.cpp b/core/src/engine/Method_GNEB.cpp index c383ee939..f0fc36b34 100644 --- a/core/src/engine/Method_GNEB.cpp +++ b/core/src/engine/Method_GNEB.cpp @@ -254,8 +254,6 @@ void Method_GNEB::Calculate_Force( { for( int img : { 0, chain->noi - 1 } ) { - scalar delta_Rx0 = 0.1; - auto & image = *configurations[img]; Manifoldmath::project_tangential( F_gradient[img], image ); @@ -264,7 +262,9 @@ void Method_GNEB::Calculate_Force( auto _F_gradient = F_gradient[img].data(); auto _tangents = tangents[img].data(); + 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] ); From cf5efbef40bc5c4ae5544a36a8690a4d2ab1ccdd Mon Sep 17 00:00:00 2001 From: Moritz Date: Wed, 2 Mar 2022 23:19:35 +0100 Subject: [PATCH 20/41] Method_GNEB: now also updates energy of 0th image (useful in case of moving endpoints) --- core/src/engine/Method_GNEB.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/core/src/engine/Method_GNEB.cpp b/core/src/engine/Method_GNEB.cpp index f0fc36b34..3e1c73ba3 100644 --- a/core/src/engine/Method_GNEB.cpp +++ b/core/src/engine/Method_GNEB.cpp @@ -376,7 +376,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]; From 2d4df714807d873ba39a1283b765cb89c84a9cd4 Mon Sep 17 00:00:00 2001 From: Moritz Date: Mon, 7 Mar 2022 13:02:17 +0100 Subject: [PATCH 21/41] Method_GNEB: Can now use attracting endpoints --- core/include/Spirit/Parameters_GNEB.h | 3 ++ core/include/data/Parameters_Method_GNEB.hpp | 2 + core/include/engine/Method_GNEB.hpp | 2 + core/python/spirit/parameters/gneb.py | 7 ++++ core/src/Spirit/Parameters_GNEB.cpp | 24 +++++++++++ core/src/engine/Method_GNEB.cpp | 43 +++++++++++++++++--- 6 files changed, 76 insertions(+), 5 deletions(-) diff --git a/core/include/Spirit/Parameters_GNEB.h b/core/include/Spirit/Parameters_GNEB.h index 3bd812d83..3c5c0a904 100644 --- a/core/include/Spirit/Parameters_GNEB.h +++ b/core/include/Spirit/Parameters_GNEB.h @@ -114,6 +114,9 @@ PREFIX void Parameters_GNEB_Set_Path_Shortening_Constant( // 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_Attracting_Endpoints( State * state, bool attracting, 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; diff --git a/core/include/data/Parameters_Method_GNEB.hpp b/core/include/data/Parameters_Method_GNEB.hpp index 213a42aa5..e358cec32 100644 --- a/core/include/data/Parameters_Method_GNEB.hpp +++ b/core/include/data/Parameters_Method_GNEB.hpp @@ -34,6 +34,8 @@ struct Parameters_Method_GNEB : Parameters_Method_Solver std::mt19937 prng = std::mt19937( rng_seed ); bool moving_endpoints = false; + bool attracting_endpoints = false; + scalar equilibrium_delta_Rx_left = 1.0; scalar equilibrium_delta_Rx_right = 1.0; diff --git a/core/include/engine/Method_GNEB.hpp b/core/include/engine/Method_GNEB.hpp index 5509b2c7d..40e08b743 100644 --- a/core/include/engine/Method_GNEB.hpp +++ b/core/include/engine/Method_GNEB.hpp @@ -77,6 +77,8 @@ class Method_GNEB : public Method_Solver vectorfield f_shrink; // Last calculated tangents std::vector tangents; + vectorfield tangent_endpoints_left; + vectorfield tangent_endpoints_right; }; } // namespace Engine diff --git a/core/python/spirit/parameters/gneb.py b/core/python/spirit/parameters/gneb.py index 5de4897a8..33d9ea70a 100644 --- a/core/python/spirit/parameters/gneb.py +++ b/core/python/spirit/parameters/gneb.py @@ -142,6 +142,13 @@ 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_Attracting_Endpoints = _spirit.Parameters_GNEB_Set_Attracting_Endpoints +_GNEB_Set_Attracting_Endpoints.argtypes = [ctypes.c_void_p, ctypes.c_bool, ctypes.c_int] +_GNEB_Set_Attracting_Endpoints.restype = None +def set_attracting_endpoints(p_state, attracting_endpoints, idx_chain=-1): + """Set if attracting endpoints should be used.""" + _GNEB_Set_Attracting_Endpoints(ctypes.c_void_p(p_state), ctypes.c_bool(attracting_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 diff --git a/core/src/Spirit/Parameters_GNEB.cpp b/core/src/Spirit/Parameters_GNEB.cpp index 896bc8817..7b9c4ec79 100644 --- a/core/src/Spirit/Parameters_GNEB.cpp +++ b/core/src/Spirit/Parameters_GNEB.cpp @@ -261,6 +261,30 @@ catch( ... ) spirit_handle_exception_api( -1, idx_chain ); } +// Set if moving endpoints should be used +void Parameters_GNEB_Set_Attracting_Endpoints( State * state, bool attracting_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->attracting_endpoints = attracting_endpoints; + chain->Unlock(); + + Log( Utility::Log_Level::Parameter, Utility::Log_Sender::API, + fmt::format( "Set GNEB attracting endpoints = {}", attracting_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 { diff --git a/core/src/engine/Method_GNEB.cpp b/core/src/engine/Method_GNEB.cpp index 3e1c73ba3..fb86def3e 100644 --- a/core/src/engine/Method_GNEB.cpp +++ b/core/src/engine/Method_GNEB.cpp @@ -40,10 +40,12 @@ 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 } ); // 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; @@ -252,6 +254,33 @@ void Method_GNEB::Calculate_Force( // Moving endpoints if( chain->gneb_parameters->moving_endpoints ) { + + scalar delta_Rx_endpoints = Manifoldmath::dist_geodesic( *this->chain->images[0]->spins, *this->chain->images[chain->noi - 1]->spins); + scalar delta_Rx_endpoints0 = chain->gneb_parameters->equilibrium_delta_Rx_left + chain->gneb_parameters->equilibrium_delta_Rx_right; + + auto _tangent_endpoints_left = tangent_endpoints_left.data(); + auto _tangent_endpoints_right = tangent_endpoints_right.data(); + auto _spins_left = this->chain->images[0]->spins->data(); + auto _spins_right = this->chain->images[chain->noi - 1]->spins->data(); + + if(chain->gneb_parameters->attracting_endpoints) + { + Backend::par::apply(nos, + [ + _tangent_endpoints_left, + _tangent_endpoints_right, + _spins_left, + _spins_right + ] + SPIRIT_LAMBDA( int idx ) + { + const Vector3 ds = _spins_right[idx] - _spins_left[idx]; + _tangent_endpoints_left[idx] = ds - ds.dot( _spins_left[idx] ) * ds; + _tangent_endpoints_right[idx] = ds - ds.dot( _spins_right[idx] ) * ds; + } + ); + } + for( int img : { 0, chain->noi - 1 } ) { auto & image = *configurations[img]; @@ -262,8 +291,11 @@ void Method_GNEB::Calculate_Force( auto _F_gradient = F_gradient[img].data(); auto _tangents = tangents[img].data(); - 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] ); + 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 _tangent_endpoints = ( img==0 ) ? _tangent_endpoints_left : _tangent_endpoints_right; auto spring_constant = ( ( img == 0 ) ? 1.0 : -1.0 ) * this->chain->gneb_parameters->spring_constant; auto projection = Vectormath::dot( F_gradient[img], tangents[img] ); @@ -274,10 +306,11 @@ void Method_GNEB::Calculate_Force( // fmt::print( "{}\n", delta_Rx ); Backend::par::apply( - nos, [_forces, _F_total, _F_gradient, delta_Rx, delta_Rx0, _tangents, spring_constant, + nos, [_forces, _F_total, _F_gradient, delta_Rx, delta_Rx0, delta_Rx_endpoints, delta_Rx_endpoints0, _tangent_endpoints, _tangents, spring_constant, projection] SPIRIT_LAMBDA( int idx ) { _forces[idx] = _F_gradient[idx] - projection * _tangents[idx] - + spring_constant * ( delta_Rx - delta_Rx0 ) * _tangents[idx]; + + spring_constant * ( delta_Rx - delta_Rx0 ) * _tangents[idx] + + spring_constant * ( delta_Rx_endpoints - delta_Rx_endpoints0 ) * _tangent_endpoints[idx]; _F_total[idx] = _forces[idx]; } ); } From 8bce187c0df3076f37f92a29da3b2ec0208eee2e Mon Sep 17 00:00:00 2001 From: "M. Sallermann" Date: Tue, 8 Mar 2022 13:55:16 +0100 Subject: [PATCH 22/41] Method_GNEB: can now set equilibrium_delta_Rx_left/right via input file --- core/src/io/Configparser.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/core/src/io/Configparser.cpp b/core/src/io/Configparser.cpp index 8cce9583c..c024257db 100644 --- a/core/src/io/Configparser.cpp +++ b/core/src/io/Configparser.cpp @@ -1053,6 +1053,9 @@ std::unique_ptr Parameters_Method_GNEB_from_Config 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->attracting_endpoints, "gneb_attracting_endpoints" ); } catch( ... ) { @@ -1076,6 +1079,9 @@ std::unique_ptr Parameters_Method_GNEB_from_Config 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} = {}", "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} = {}", "attracting_endpoints", parameters->attracting_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 ) ); From 5b7f195452ea9e571323fdec63191c9148a3238d Mon Sep 17 00:00:00 2001 From: Moritz Date: Wed, 16 Mar 2022 16:38:59 +0100 Subject: [PATCH 23/41] Ui-CPP: Added two options to the plot widget 1) divide by nos and 2) renormalize Rx. Both are useful when observing moving_endpoints calculations --- ui-cpp/ui-qt/include/PlotWidget.hpp | 2 + ui-cpp/ui-qt/src/PlotWidget.cpp | 47 +++++++++++++++++--- ui-cpp/ui-qt/src/PlotsWidget.cpp | 6 +++ ui-cpp/ui-qt/ui/PlotsWidget.ui | 68 +++++++++++++++++------------ 4 files changed, 91 insertions(+), 32 deletions(-) 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/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/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 + + + From 17771debc7b348fe230ef4b311c681a7d44184df Mon Sep 17 00:00:00 2001 From: Moritz Date: Wed, 16 Mar 2022 16:49:35 +0100 Subject: [PATCH 24/41] Method_GNEB: Replaced attracting endpoints with translating endpoints --- core/include/data/Parameters_Method_GNEB.hpp | 2 +- core/include/engine/Method_GNEB.hpp | 5 + core/src/Spirit/Parameters_GNEB.cpp | 8 +- core/src/engine/Method_GNEB.cpp | 104 +++++++++++-------- core/src/io/Configparser.cpp | 4 +- 5 files changed, 74 insertions(+), 49 deletions(-) diff --git a/core/include/data/Parameters_Method_GNEB.hpp b/core/include/data/Parameters_Method_GNEB.hpp index e358cec32..d52818d95 100644 --- a/core/include/data/Parameters_Method_GNEB.hpp +++ b/core/include/data/Parameters_Method_GNEB.hpp @@ -34,7 +34,7 @@ struct Parameters_Method_GNEB : Parameters_Method_Solver std::mt19937 prng = std::mt19937( rng_seed ); bool moving_endpoints = false; - bool attracting_endpoints = false; + bool translating_endpoints = false; scalar equilibrium_delta_Rx_left = 1.0; scalar equilibrium_delta_Rx_right = 1.0; diff --git a/core/include/engine/Method_GNEB.hpp b/core/include/engine/Method_GNEB.hpp index 40e08b743..572fe383b 100644 --- a/core/include/engine/Method_GNEB.hpp +++ b/core/include/engine/Method_GNEB.hpp @@ -75,8 +75,13 @@ 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; }; diff --git a/core/src/Spirit/Parameters_GNEB.cpp b/core/src/Spirit/Parameters_GNEB.cpp index 7b9c4ec79..c25ab71cc 100644 --- a/core/src/Spirit/Parameters_GNEB.cpp +++ b/core/src/Spirit/Parameters_GNEB.cpp @@ -261,8 +261,8 @@ catch( ... ) spirit_handle_exception_api( -1, idx_chain ); } -// Set if moving endpoints should be used -void Parameters_GNEB_Set_Attracting_Endpoints( State * state, bool attracting_endpoints, int idx_chain ) noexcept +// 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; @@ -274,11 +274,11 @@ try chain->Lock(); auto p = chain->gneb_parameters; - p->attracting_endpoints = attracting_endpoints; + p->moving_endpoints = translating_endpoints; chain->Unlock(); Log( Utility::Log_Level::Parameter, Utility::Log_Sender::API, - fmt::format( "Set GNEB attracting endpoints = {}", attracting_endpoints ), idx_image, idx_chain ); + fmt::format( "Set GNEB translating endpoints = {}", translating_endpoints ), idx_image, idx_chain ); } catch( ... ) { diff --git a/core/src/engine/Method_GNEB.cpp b/core/src/engine/Method_GNEB.cpp index fb86def3e..912a84bc0 100644 --- a/core/src/engine/Method_GNEB.cpp +++ b/core/src/engine/Method_GNEB.cpp @@ -9,7 +9,7 @@ #include #include #include - +#include #include #include @@ -42,6 +42,10 @@ Method_GNEB::Method_GNEB( std::shared_ptr chain this->f_shrink = 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] @@ -254,65 +258,81 @@ void Method_GNEB::Calculate_Force( // Moving endpoints if( chain->gneb_parameters->moving_endpoints ) { - - scalar delta_Rx_endpoints = Manifoldmath::dist_geodesic( *this->chain->images[0]->spins, *this->chain->images[chain->noi - 1]->spins); - scalar delta_Rx_endpoints0 = chain->gneb_parameters->equilibrium_delta_Rx_left + chain->gneb_parameters->equilibrium_delta_Rx_right; + // Overall translational force + int noi = chain->noi; + Manifoldmath::project_tangential( F_gradient[0], *configurations[0] ); + Manifoldmath::project_tangential( F_gradient[noi-1], *configurations[noi-1] ); - auto _tangent_endpoints_left = tangent_endpoints_left.data(); - auto _tangent_endpoints_right = tangent_endpoints_right.data(); - auto _spins_left = this->chain->images[0]->spins->data(); - auto _spins_right = this->chain->images[chain->noi - 1]->spins->data(); - - if(chain->gneb_parameters->attracting_endpoints) + if(chain->gneb_parameters->translating_endpoints) { - Backend::par::apply(nos, + Backend::par::apply(nos, [ - _tangent_endpoints_left, - _tangent_endpoints_right, - _spins_left, - _spins_right - ] - SPIRIT_LAMBDA( int idx ) + 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 ds = _spins_right[idx] - _spins_left[idx]; - _tangent_endpoints_left[idx] = ds - ds.dot( _spins_left[idx] ) * ds; - _tangent_endpoints_right[idx] = ds - ds.dot( _spins_right[idx] ) * ds; + Vector3 axis = spins_left[idx].cross(spins_right[idx]); + 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 (std::abs(angle) < 1e-6 || std::isnan(angle)) // Angle can become nan for collinear spins + rotation_matrix = Matrix3::Identity(); + + 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); + + 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]); } ); + Manifoldmath::project_parallel( F_translation_left, tangents[0] ); + Manifoldmath::project_parallel( F_translation_right, tangents[noi-1]); } for( int img : { 0, chain->noi - 1 } ) { - auto & image = *configurations[img]; - Manifoldmath::project_tangential( F_gradient[img], image ); - - auto _F_total = F_total[img].data(); - auto _forces = forces[img].data(); - auto _F_gradient = F_gradient[img].data(); - auto _tangents = tangents[img].data(); - 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 _tangent_endpoints = ( img==0 ) ? _tangent_endpoints_left : _tangent_endpoints_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( "{}\n", _F_total[0].transpose() ); - // fmt::print( "{}\n", _F_gradient[0].transpose() ); - // fmt::print( "{}\n", delta_Rx ); + // 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 ); Backend::par::apply( - nos, [_forces, _F_total, _F_gradient, delta_Rx, delta_Rx0, delta_Rx_endpoints, delta_Rx_endpoints0, _tangent_endpoints, _tangents, spring_constant, - projection] SPIRIT_LAMBDA( int idx ) { - _forces[idx] = _F_gradient[idx] - projection * _tangents[idx] - + spring_constant * ( delta_Rx - delta_Rx0 ) * _tangents[idx] - + spring_constant * ( delta_Rx_endpoints - delta_Rx_endpoints0 ) * _tangent_endpoints[idx]; - _F_total[idx] = _forces[idx]; - } ); + 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 * std::tanh(delta_Rx - delta_Rx0), + F_translation, + projection + ] SPIRIT_LAMBDA ( int idx ) + { + forces[idx] = 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]; + } + ); + + Manifoldmath::project_tangential( F_gradient[0], *configurations[0] ); + Manifoldmath::project_tangential( F_gradient[noi-1], *configurations[noi-1] ); } } diff --git a/core/src/io/Configparser.cpp b/core/src/io/Configparser.cpp index c024257db..588749b33 100644 --- a/core/src/io/Configparser.cpp +++ b/core/src/io/Configparser.cpp @@ -1055,7 +1055,7 @@ std::unique_ptr Parameters_Method_GNEB_from_Config 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->attracting_endpoints, "gneb_attracting_endpoints" ); + config_file_handle.Read_Single( parameters->translating_endpoints, "gneb_translating_endpoints" ); } catch( ... ) { @@ -1081,7 +1081,7 @@ std::unique_ptr Parameters_Method_GNEB_from_Config 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} = {}", "attracting_endpoints", parameters->attracting_endpoints ) ); + 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 ) ); From bd830154337e1a4cc054e6e145497996482d2de2 Mon Sep 17 00:00:00 2001 From: "M. Sallermann" Date: Thu, 17 Mar 2022 14:33:07 +0100 Subject: [PATCH 25/41] formatted Method_GNEB.cpp --- core/src/engine/Method_GNEB.cpp | 41 ++++++++++++++++++--------------- 1 file changed, 23 insertions(+), 18 deletions(-) diff --git a/core/src/engine/Method_GNEB.cpp b/core/src/engine/Method_GNEB.cpp index 912a84bc0..8494f52d3 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,12 +8,10 @@ #include #include #include +#include #include #include #include -#include -#include -#include #include @@ -45,11 +45,10 @@ Method_GNEB::Method_GNEB( std::shared_ptr chain 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] + 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; @@ -148,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 ) { @@ -261,10 +260,11 @@ void Method_GNEB::Calculate_Force( // Overall translational force int noi = chain->noi; Manifoldmath::project_tangential( F_gradient[0], *configurations[0] ); - Manifoldmath::project_tangential( F_gradient[noi-1], *configurations[noi-1] ); + Manifoldmath::project_tangential( F_gradient[noi - 1], *configurations[noi - 1] ); - if(chain->gneb_parameters->translating_endpoints) + if( chain->gneb_parameters->translating_endpoints ) { + // clang-format off Backend::par::apply(nos, [ F_translation_left = F_translation_left.data(), @@ -277,7 +277,7 @@ void Method_GNEB::Calculate_Force( { Vector3 axis = spins_left[idx].cross(spins_right[idx]); 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(); @@ -291,14 +291,17 @@ void Method_GNEB::Calculate_Force( F_translation_right[idx] = -0.5 * (F_gradient_left_rotated + F_gradient_right[idx]); } ); - Manifoldmath::project_parallel( F_translation_left, tangents[0] ); - Manifoldmath::project_parallel( F_translation_right, tangents[noi-1]); + // clang-format on + + Manifoldmath::project_parallel( F_translation_left, tangents[0] ); + Manifoldmath::project_parallel( F_translation_right, tangents[noi - 1] ); } 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]; + 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] ); @@ -310,8 +313,9 @@ void Method_GNEB::Calculate_Force( // fmt::print( "delta_Rx {}\n", delta_Rx ); // fmt::print( "delta_Rx0 {}\n", delta_Rx0 ); + // clang-format off Backend::par::apply( - nos, + nos, [ F_total = F_total[img].data(), F_gradient = F_gradient[img].data(), @@ -322,7 +326,7 @@ void Method_GNEB::Calculate_Force( projection ] SPIRIT_LAMBDA ( int idx ) { - forces[idx] = F_gradient[idx] - projection * tangents[idx] + forces[idx] = F_gradient[idx] - projection * tangents[idx] + tangent_coeff * tangents[idx] + F_translation[idx]; @@ -330,9 +334,10 @@ void Method_GNEB::Calculate_Force( F_total[idx] = forces[idx]; } ); + // clang-format on Manifoldmath::project_tangential( F_gradient[0], *configurations[0] ); - Manifoldmath::project_tangential( F_gradient[noi-1], *configurations[noi-1] ); + Manifoldmath::project_tangential( F_gradient[noi - 1], *configurations[noi - 1] ); } } From 733d9821c1fb47d9131daebdc2d82e95b3dd23d7 Mon Sep 17 00:00:00 2001 From: "M. Sallermann" Date: Thu, 17 Mar 2022 14:42:59 +0100 Subject: [PATCH 26/41] Method_GNEB: Adjusted API functions from attracting_endpoints to translating_endpoints --- core/include/Spirit/Parameters_GNEB.h | 2 +- core/python/spirit/parameters/gneb.py | 10 +++++----- core/src/Spirit/Parameters_GNEB.cpp | 4 ++-- core/src/engine/Method_GNEB.cpp | 2 +- 4 files changed, 9 insertions(+), 9 deletions(-) diff --git a/core/include/Spirit/Parameters_GNEB.h b/core/include/Spirit/Parameters_GNEB.h index 3c5c0a904..316b767f6 100644 --- a/core/include/Spirit/Parameters_GNEB.h +++ b/core/include/Spirit/Parameters_GNEB.h @@ -115,7 +115,7 @@ PREFIX void Parameters_GNEB_Set_Path_Shortening_Constant( 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_Attracting_Endpoints( State * state, bool attracting, int idx_chain = -1 ) SUFFIX; +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; diff --git a/core/python/spirit/parameters/gneb.py b/core/python/spirit/parameters/gneb.py index 33d9ea70a..52c3b95ea 100644 --- a/core/python/spirit/parameters/gneb.py +++ b/core/python/spirit/parameters/gneb.py @@ -142,12 +142,12 @@ 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_Attracting_Endpoints = _spirit.Parameters_GNEB_Set_Attracting_Endpoints -_GNEB_Set_Attracting_Endpoints.argtypes = [ctypes.c_void_p, ctypes.c_bool, ctypes.c_int] -_GNEB_Set_Attracting_Endpoints.restype = None -def set_attracting_endpoints(p_state, attracting_endpoints, idx_chain=-1): +_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_Attracting_Endpoints(ctypes.c_void_p(p_state), ctypes.c_bool(attracting_endpoints), ctypes.c_int(idx_chain)) + _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] diff --git a/core/src/Spirit/Parameters_GNEB.cpp b/core/src/Spirit/Parameters_GNEB.cpp index c25ab71cc..6fb9b3862 100644 --- a/core/src/Spirit/Parameters_GNEB.cpp +++ b/core/src/Spirit/Parameters_GNEB.cpp @@ -273,8 +273,8 @@ try from_indices( state, idx_image, idx_chain, image, chain ); chain->Lock(); - auto p = chain->gneb_parameters; - p->moving_endpoints = translating_endpoints; + auto p = chain->gneb_parameters; + p->translating_endpoints = translating_endpoints; chain->Unlock(); Log( Utility::Log_Level::Parameter, Utility::Log_Sender::API, diff --git a/core/src/engine/Method_GNEB.cpp b/core/src/engine/Method_GNEB.cpp index 8494f52d3..a6fa490ca 100644 --- a/core/src/engine/Method_GNEB.cpp +++ b/core/src/engine/Method_GNEB.cpp @@ -257,11 +257,11 @@ void Method_GNEB::Calculate_Force( // Moving endpoints if( chain->gneb_parameters->moving_endpoints ) { - // Overall translational force 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 From adfacce39c621cb57da9241e0ec07a9b0497938d Mon Sep 17 00:00:00 2001 From: "M. Sallermann" Date: Thu, 17 Mar 2022 14:43:31 +0100 Subject: [PATCH 27/41] Simulation.cpp: Can now start GNEB calculations with 2 images. --- core/src/Spirit/Simulation.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/core/src/Spirit/Simulation.cpp b/core/src/Spirit/Simulation.cpp index 18012e279..c6e68d7c6 100644 --- a/core/src/Spirit/Simulation.cpp +++ b/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 From 6b0c6129f2c2989ba080741169754e8af320eebb Mon Sep 17 00:00:00 2001 From: "M. Sallermann" Date: Thu, 17 Mar 2022 14:51:00 +0100 Subject: [PATCH 28/41] Parameters_GNEB: implemented API getters for translating_endpoints --- core/include/Spirit/Parameters_GNEB.h | 3 +++ core/python/spirit/parameters/gneb.py | 10 +++++++++- core/src/Spirit/Parameters_GNEB.cpp | 19 +++++++++++++++++++ 3 files changed, 31 insertions(+), 1 deletion(-) diff --git a/core/include/Spirit/Parameters_GNEB.h b/core/include/Spirit/Parameters_GNEB.h index 316b767f6..b1e7ae4f8 100644 --- a/core/include/Spirit/Parameters_GNEB.h +++ b/core/include/Spirit/Parameters_GNEB.h @@ -181,6 +181,9 @@ PREFIX float Parameters_GNEB_Get_Path_Shortening_Constant( State * state, int id // 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; + /* Returns the integer of whether an image is regular, climbing, falling, or stationary. diff --git a/core/python/spirit/parameters/gneb.py b/core/python/spirit/parameters/gneb.py index 52c3b95ea..39637e7c5 100644 --- a/core/python/spirit/parameters/gneb.py +++ b/core/python/spirit/parameters/gneb.py @@ -225,7 +225,15 @@ def get_path_shortening_constant(p_state, idx_image=-1, idx_chain=-1): 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))) + 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_Climbing_Falling = _spirit.Parameters_GNEB_Get_Climbing_Falling _GNEB_Get_Climbing_Falling.argtypes = [ctypes.c_void_p, ctypes.c_int, ctypes.c_int] diff --git a/core/src/Spirit/Parameters_GNEB.cpp b/core/src/Spirit/Parameters_GNEB.cpp index 6fb9b3862..2dee7eca0 100644 --- a/core/src/Spirit/Parameters_GNEB.cpp +++ b/core/src/Spirit/Parameters_GNEB.cpp @@ -600,6 +600,25 @@ catch( ... ) 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; +} + int Parameters_GNEB_Get_Climbing_Falling( State * state, int idx_image, int idx_chain ) noexcept try { From cfe8fd861fade14518ac6168d3abe569491c58b0 Mon Sep 17 00:00:00 2001 From: "M. Sallermann" Date: Thu, 17 Mar 2022 15:02:26 +0100 Subject: [PATCH 29/41] GNEB: API getter for equilibrium_delta_Rx --- core/include/Spirit/Parameters_GNEB.h | 3 +++ core/python/spirit/parameters/gneb.py | 12 ++++++++++++ core/src/Spirit/Parameters_GNEB.cpp | 19 +++++++++++++++++++ 3 files changed, 34 insertions(+) diff --git a/core/include/Spirit/Parameters_GNEB.h b/core/include/Spirit/Parameters_GNEB.h index b1e7ae4f8..ff25fd0a2 100644 --- a/core/include/Spirit/Parameters_GNEB.h +++ b/core/include/Spirit/Parameters_GNEB.h @@ -184,6 +184,9 @@ PREFIX bool Parameters_GNEB_Get_Moving_Endpoints( State * state, int idx_chain = // 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/python/spirit/parameters/gneb.py b/core/python/spirit/parameters/gneb.py index 39637e7c5..1c46dc60d 100644 --- a/core/python/spirit/parameters/gneb.py +++ b/core/python/spirit/parameters/gneb.py @@ -235,6 +235,18 @@ def get_translating_endpoints(p_state, idx_chain=-1): 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/src/Spirit/Parameters_GNEB.cpp b/core/src/Spirit/Parameters_GNEB.cpp index 2dee7eca0..27e9a500b 100644 --- a/core/src/Spirit/Parameters_GNEB.cpp +++ b/core/src/Spirit/Parameters_GNEB.cpp @@ -619,6 +619,25 @@ catch( ... ) 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 { From 5dee3d9dec787bf994567aac58bdfe2a2b8bf008 Mon Sep 17 00:00:00 2001 From: "M. Sallermann" Date: Thu, 17 Mar 2022 15:37:09 +0100 Subject: [PATCH 30/41] Ui-QT: added moving endpoints to parameters widget --- ui-cpp/ui-qt/src/ParametersWidget.cpp | 27 + ui-cpp/ui-qt/ui/ParametersWidget.ui | 801 ++++++++++++++------------ 2 files changed, 464 insertions(+), 364 deletions(-) diff --git a/ui-cpp/ui-qt/src/ParametersWidget.cpp b/ui-cpp/ui-qt/src/ParametersWidget.cpp index 9c8843698..f45f89c47 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() @@ -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/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 @@ - + From b96c765ecbd7f9c00e0ed59c197a721d662db392 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Mon, 21 Mar 2022 10:33:27 +0100 Subject: [PATCH 31/41] Method_GNEB: added escape_first parameters, which forbids rotations in the convex region --- core/include/data/Parameters_Method_GNEB.hpp | 2 ++ core/src/engine/Method_GNEB.cpp | 22 ++++++++++++++++---- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/core/include/data/Parameters_Method_GNEB.hpp b/core/include/data/Parameters_Method_GNEB.hpp index d52818d95..6b4b0165b 100644 --- a/core/include/data/Parameters_Method_GNEB.hpp +++ b/core/include/data/Parameters_Method_GNEB.hpp @@ -39,6 +39,8 @@ struct Parameters_Method_GNEB : Parameters_Method_Solver 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/src/engine/Method_GNEB.cpp b/core/src/engine/Method_GNEB.cpp index a6fa490ca..4d6c9dfe7 100644 --- a/core/src/engine/Method_GNEB.cpp +++ b/core/src/engine/Method_GNEB.cpp @@ -294,7 +294,19 @@ void Method_GNEB::Calculate_Force( // clang-format on Manifoldmath::project_parallel( F_translation_left, tangents[0] ); - Manifoldmath::project_parallel( F_translation_right, tangents[noi - 1] ); + 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 } ) @@ -326,7 +338,7 @@ void Method_GNEB::Calculate_Force( projection ] SPIRIT_LAMBDA ( int idx ) { - forces[idx] = F_gradient[idx] - projection * tangents[idx] + forces[idx] = rotational_coeff * (F_gradient[idx] - projection * tangents[idx]) + tangent_coeff * tangents[idx] + F_translation[idx]; @@ -336,8 +348,10 @@ void Method_GNEB::Calculate_Force( ); // clang-format on - Manifoldmath::project_tangential( F_gradient[0], *configurations[0] ); - Manifoldmath::project_tangential( F_gradient[noi - 1], *configurations[noi - 1] ); + + // Manifoldmath::project_tangential( F_gradient[0], *configurations[0] ); + // Manifoldmath::project_tangential( F_gradient[noi - 1], *configurations[noi - 1] ); + } } From 299f885d97c7a8a13e87ac458f477102ec07dff6 Mon Sep 17 00:00:00 2001 From: "M. Sallermann" Date: Sun, 27 Mar 2022 18:34:20 +0200 Subject: [PATCH 32/41] Ui-Imgui: Added moving endpoints settings to parameters widget --- ui-cpp/ui-imgui/include/parameters_widget.hpp | 6 +++ ui-cpp/ui-imgui/src/parameters_widget.cpp | 45 +++++++++++++++++++ 2 files changed, 51 insertions(+) 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..2ac6c3253 100644 --- a/ui-cpp/ui-imgui/src/parameters_widget.cpp +++ b/ui-cpp/ui-imgui/src/parameters_widget.cpp @@ -550,6 +550,51 @@ 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 ) { From 4827222ddc761307ef76cea6a417e241a54efc37 Mon Sep 17 00:00:00 2001 From: "M. Sallermann" Date: Sun, 27 Mar 2022 18:49:26 +0200 Subject: [PATCH 33/41] Method_GNEB: fixed a bug in a lambda capture --- core/src/engine/Method_GNEB.cpp | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/core/src/engine/Method_GNEB.cpp b/core/src/engine/Method_GNEB.cpp index 4d6c9dfe7..3e1476cfd 100644 --- a/core/src/engine/Method_GNEB.cpp +++ b/core/src/engine/Method_GNEB.cpp @@ -335,7 +335,8 @@ void Method_GNEB::Calculate_Force( tangents = tangents[img].data(), tangent_coeff = spring_constant * std::tanh(delta_Rx - delta_Rx0), F_translation, - projection + projection, + rotational_coeff ] SPIRIT_LAMBDA ( int idx ) { forces[idx] = rotational_coeff * (F_gradient[idx] - projection * tangents[idx]) @@ -347,11 +348,6 @@ void Method_GNEB::Calculate_Force( } ); // clang-format on - - - // Manifoldmath::project_tangential( F_gradient[0], *configurations[0] ); - // Manifoldmath::project_tangential( F_gradient[noi - 1], *configurations[noi - 1] ); - } } From be9bceac5ef909c65d866fa81fb8d27dece220b9 Mon Sep 17 00:00:00 2001 From: "M. Sallermann" Date: Tue, 29 Mar 2022 13:45:23 +0200 Subject: [PATCH 34/41] Manifoldmath: Added a new function Geodesic_Tangent that computes the tangents more accurately. Use it to compute the tangents at the endpoints --- core/src/engine/Manifoldmath.cpp | 55 +++++++++++++++++++++++++------- 1 file changed, 44 insertions(+), 11 deletions(-) diff --git a/core/src/engine/Manifoldmath.cpp b/core/src/engine/Manifoldmath.cpp index e9c55f934..0fc349d5d 100644 --- a/core/src/engine/Manifoldmath.cpp +++ b/core/src/engine/Manifoldmath.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include @@ -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 From 9f02f90813f2f671dd8f1b2e9e3e4b4abef3f9d6 Mon Sep 17 00:00:00 2001 From: "M. Sallermann" Date: Tue, 29 Mar 2022 13:47:49 +0200 Subject: [PATCH 35/41] Method_GNEB: After some testing removed tanh function in tangent_coeff --- core/src/engine/Method_GNEB.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/core/src/engine/Method_GNEB.cpp b/core/src/engine/Method_GNEB.cpp index 3e1476cfd..d2d17291e 100644 --- a/core/src/engine/Method_GNEB.cpp +++ b/core/src/engine/Method_GNEB.cpp @@ -333,7 +333,7 @@ void Method_GNEB::Calculate_Force( F_gradient = F_gradient[img].data(), forces = forces[img].data(), tangents = tangents[img].data(), - tangent_coeff = spring_constant * std::tanh(delta_Rx - delta_Rx0), + tangent_coeff = spring_constant * (delta_Rx - delta_Rx0), F_translation, projection, rotational_coeff From d07ce032fb6a0ee8c82475e5fcfd3e6375c659f8 Mon Sep 17 00:00:00 2001 From: "M. Sallermann" Date: Tue, 29 Mar 2022 13:56:21 +0200 Subject: [PATCH 36/41] Ui-Imgui: added correct getters for moving endpoints data and formatted --- ui-cpp/ui-imgui/src/parameters_widget.cpp | 37 ++++++++++------------- 1 file changed, 16 insertions(+), 21 deletions(-) diff --git a/ui-cpp/ui-imgui/src/parameters_widget.cpp b/ui-cpp/ui-imgui/src/parameters_widget.cpp index 2ac6c3253..33203fe30 100644 --- a/ui-cpp/ui-imgui/src/parameters_widget.cpp +++ b/ui-cpp/ui-imgui/src/parameters_widget.cpp @@ -551,48 +551,38 @@ void ParametersWidget::show_content() Parameters_GNEB_Set_Path_Shortening_Constant( state.get(), parameters_gneb.path_shortening_constant ); } - ImGui::BeginTable("MovingEndpointsTable", 2); + ImGui::BeginTable( "MovingEndpointsTable", 2 ); ImGui::TableNextColumn(); ImGui::TextUnformatted( "Moving endpoints" ); ImGui::TableNextColumn(); - if( ImGui::Checkbox( - "##gneb_moving_endpoints", - ¶meters_gneb.moving_endpoints ) ) + if( ImGui::Checkbox( "##gneb_moving_endpoints", ¶meters_gneb.moving_endpoints ) ) { - Parameters_GNEB_Set_Moving_Endpoints( - state.get(), parameters_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 ) ) + if( ImGui::Checkbox( "##gneb_translating_endpoints", ¶meters_gneb.translating_endpoints ) ) { - Parameters_GNEB_Set_Translating_Endpoints( - state.get(), parameters_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 ) ) + 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 ); + 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 ) ) + 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 ); + Parameters_GNEB_Set_Equilibrium_Delta_Rx( + state.get(), parameters_gneb.delta_Rx_left, parameters_gneb.delta_Rx_right ); } ImGui::EndTable(); } @@ -813,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 From 21eba6fc82d39fe8826a71e04f45d691a9a2fbc9 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Sun, 10 Apr 2022 21:02:24 +0200 Subject: [PATCH 37/41] Sparse_HTST: Use sparse lanczos for computation of unstable mode and zero modes --- core/src/engine/Sparse_HTST.cpp | 58 +++++++++++++++++++++------------ 1 file changed, 38 insertions(+), 20 deletions(-) diff --git a/core/src/engine/Sparse_HTST.cpp b/core/src/engine/Sparse_HTST.cpp index 6752a22f9..79417e974 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]; @@ -430,7 +444,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 ..." ); From aca472f79a0ddd7e4ed67fdf9319c2e51ea9ee14 Mon Sep 17 00:00:00 2001 From: Moritz Date: Wed, 20 Apr 2022 11:19:54 +0200 Subject: [PATCH 38/41] Manifoldmath: more performant implementation of project_parallel --- core/src/engine/Manifoldmath.cpp | 12 ++++++------ core/src/engine/Manifoldmath.cu | 10 +++++++--- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/core/src/engine/Manifoldmath.cpp b/core/src/engine/Manifoldmath.cpp index 0fc349d5d..cd4affd54 100644 --- a/core/src/engine/Manifoldmath.cpp +++ b/core/src/engine/Manifoldmath.cpp @@ -22,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 ) 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) From c4897914da0d4ed1af9e87f359f4c222a0da931d Mon Sep 17 00:00:00 2001 From: Moritz Date: Wed, 20 Apr 2022 11:20:46 +0200 Subject: [PATCH 39/41] Method_GNEB: slight improvements in calculation of translational force --- core/src/engine/Method_GNEB.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/core/src/engine/Method_GNEB.cpp b/core/src/engine/Method_GNEB.cpp index d2d17291e..5d410bc09 100644 --- a/core/src/engine/Method_GNEB.cpp +++ b/core/src/engine/Method_GNEB.cpp @@ -275,19 +275,19 @@ void Method_GNEB::Calculate_Force( spins_right = this->chain->images[noi-1]->spins->data() ] SPIRIT_LAMBDA ( int idx) { - Vector3 axis = spins_left[idx].cross(spins_right[idx]); - scalar angle = acos(spins_left[idx].dot(spins_right[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 (std::abs(angle) < 1e-6 || std::isnan(angle)) // Angle can become nan for collinear spins + if ( abs(spins_left[idx].dot(spins_right[idx])) >= 1.0 ) // Angle can become nan for collinear spins rotation_matrix = Matrix3::Identity(); - Vector3 F_gradient_right_rotated = rotation_matrix * F_gradient_right[idx]; + 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); - Vector3 F_gradient_left_rotated = rotation_matrix.transpose() * F_gradient_left[idx]; + 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]); } ); From a4b451785dd30a4299a0f2bedbccf7669ca68e93 Mon Sep 17 00:00:00 2001 From: Moritz Date: Mon, 25 Apr 2022 11:05:24 +0200 Subject: [PATCH 40/41] Timing: Fixed issues in DurationFromString Previously this function was implemented wrongly, which caused problems with the max_walltime of simulations. --- core/src/utility/Timing.cpp | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) 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 ); } From f49c11432087e9b4d3e481acacbbe580717f2a9e Mon Sep 17 00:00:00 2001 From: "G. P. Mueller" Date: Mon, 25 Apr 2022 21:54:38 +0200 Subject: [PATCH 41/41] Core: minor, nonintrusive improvements - remove some unnecessary semicolons - added missing `Exception` copy and move constructors and operators definitions, add `noexcept` specifiers - fixed unnecessary copying in `Geometry` constructor --- core/include/data/Geometry.hpp | 5 +- core/include/engine/Hamiltonian.hpp | 2 +- core/include/engine/Hamiltonian_Gaussian.hpp | 2 +- .../include/engine/Hamiltonian_Heisenberg.hpp | 2 +- core/include/engine/Solver_Depondt.hpp | 8 +- core/include/engine/Solver_Heun.hpp | 8 +- core/include/engine/Solver_Kernels.hpp | 1 + core/include/engine/Solver_RK4.hpp | 8 +- core/include/engine/Solver_SIB.hpp | 8 +- core/include/engine/Solver_VP.hpp | 8 +- core/include/engine/Solver_VP_OSO.hpp | 6 +- core/include/utility/Exception.hpp | 51 +++++++++---- core/src/data/Geometry.cpp | 5 +- core/src/engine/Hamiltonian.cpp | 8 -- core/src/engine/Hamiltonian_Gaussian.cpp | 2 +- core/src/engine/Hamiltonian_Heisenberg.cpp | 2 +- core/src/engine/Hamiltonian_Heisenberg.cu | 5 +- core/src/utility/Exception.cpp | 75 +++++-------------- 18 files changed, 94 insertions(+), 112 deletions(-) diff --git a/core/include/data/Geometry.hpp b/core/include/data/Geometry.hpp index 53b0e80a1..8e33f023d 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/engine/Hamiltonian.hpp b/core/include/engine/Hamiltonian.hpp index d1c95005b..0a5519c4e 100644 --- a/core/include/engine/Hamiltonian.hpp +++ b/core/include/engine/Hamiltonian.hpp @@ -86,7 +86,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 db3cfe275..056ee5aea 100644 --- a/core/include/engine/Hamiltonian_Heisenberg.hpp +++ b/core/include/engine/Hamiltonian_Heisenberg.hpp @@ -64,7 +64,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/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/utility/Exception.hpp b/core/include/utility/Exception.hpp index d1f70f713..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,25 +72,38 @@ 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 #define spirit_throw( classifier, level, message ) \ diff --git a/core/src/data/Geometry.cpp b/core/src/data/Geometry.cpp index 288e5c9d8..5b9cf2834 100644 --- a/core/src/data/Geometry.cpp +++ b/core/src/data/Geometry.cpp @@ -20,8 +20,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/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 d84a9f0ce..2419163bc 100644 --- a/core/src/engine/Hamiltonian_Heisenberg.cpp +++ b/core/src/engine/Hamiltonian_Heisenberg.cpp @@ -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/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