diff --git a/core/include/engine/Vectormath_Defines.hpp b/core/include/engine/Vectormath_Defines.hpp index a5687913b..b7be1d51b 100644 --- a/core/include/engine/Vectormath_Defines.hpp +++ b/core/include/engine/Vectormath_Defines.hpp @@ -144,6 +144,7 @@ struct Neighbour : Pair using intfield = field; using scalarfield = field; using vectorfield = field; +using matrixfield = field; // Additional fields using pairfield = field; diff --git a/core/include/engine/spin/Hamiltonian.hpp b/core/include/engine/spin/Hamiltonian.hpp index 8c92fc365..01a24b7fa 100644 --- a/core/include/engine/spin/Hamiltonian.hpp +++ b/core/include/engine/spin/Hamiltonian.hpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include @@ -102,8 +103,8 @@ struct HamiltonianVariantTypes using Gaussian = Hamiltonian; using Heisenberg = Hamiltonian< state_t, AdaptorType, Interaction::Zeeman, Interaction::Anisotropy, Interaction::Biaxial_Anisotropy, - Interaction::Cubic_Anisotropy, Interaction::Exchange, Interaction::DMI, Interaction::Quadruplet, - Interaction::DDI>; + Interaction::Cubic_Anisotropy, Interaction::Exchange, Interaction::DMI, Interaction::Two_Site_Anisotropy, + Interaction::Quadruplet, Interaction::DDI>; using Variant = std::variant; }; diff --git a/core/include/engine/spin/interaction/Two_Site_Anisotropy.hpp b/core/include/engine/spin/interaction/Two_Site_Anisotropy.hpp new file mode 100644 index 000000000..239cea6df --- /dev/null +++ b/core/include/engine/spin/interaction/Two_Site_Anisotropy.hpp @@ -0,0 +1,176 @@ +#pragma once +#ifndef SPIRIT_CORE_ENGINE_INTERACTION_TWO_ION_ANISOTROPY_HPP +#define SPIRIT_CORE_ENGINE_INTERACTION_TWO_ION_ANISOTROPY_HPP + +#include +#include +#include +#include + +#include + +namespace Engine +{ + +namespace Spin +{ + +namespace Interaction +{ + +struct Two_Site_Anisotropy +{ + using state_t = StateType; + + struct Data + { + pairfield pairs{}; + matrixfield matrices{}; + + Data() = default; + Data( pairfield pairs, matrixfield matrices ) + : pairs( std::move( pairs ) ), matrices( std::move( matrices ) ) {}; + }; + + static bool valid_data( const Data & data ) + { + if( !data.pairs.empty() && ( data.pairs.size() != data.matrices.size() ) ) + return false; + + return true; + }; + + struct Cache + { + pairfield pairs{}; + matrixfield matrices{}; + }; + + static bool is_contributing( const Data &, const Cache & cache ) + { + return !cache.pairs.empty(); + } + + struct IndexType + { + int ispin, jspin, ipair; + }; + + using Index = Engine::Span; + using IndexStorage = Backend::vector; + + using Energy = Functor::Local::Energy_Functor>; + using Gradient = Functor::Local::Gradient_Functor>; + using Hessian = Functor::Local::Hessian_Functor>; + + static std::size_t Sparse_Hessian_Size_per_Cell( const Data &, const Cache & cache ) + { + return cache.pairs.size() * 9; + }; + + // Calculate the total energy for a single spin to be used in Monte Carlo. + // Note: therefore the energy of pairs is weighted x2 and of quadruplets x4. + using Energy_Single_Spin = Functor::Local::Energy_Single_Spin_Functor; + + // Interaction name as string + static constexpr std::string_view name = "Two Site Anisotropy"; + + template + static void applyGeometry( + const ::Data::Geometry & geometry, const intfield & boundary_conditions, const Data & data, Cache & cache, + IndexStorageVector & indices ) + { + using Indexing::idx_from_pair; + + // redundant neighbours are captured when expanding pairs below + static constexpr bool use_redundant_neighbours = false; + + // Use direct list of pairs + // TODO: find a more convenient approach for implementing these in terms of neighbours. + cache.pairs = data.pairs; + cache.matrices = data.matrices; + + for( unsigned int icell = 0; icell < geometry.n_cells_total; ++icell ) + { + for( unsigned int i_pair = 0; i_pair < cache.pairs.size(); ++i_pair ) + { + int ispin = cache.pairs[i_pair].i + icell * geometry.n_cell_atoms; + int jspin = idx_from_pair( + ispin, boundary_conditions, geometry.n_cells, geometry.n_cell_atoms, geometry.atom_types, + cache.pairs[i_pair] ); + if( jspin >= 0 ) + { + Backend::get( indices[ispin] ).push_back( IndexType{ ispin, jspin, (int)i_pair } ); + Backend::get( indices[jspin] ).push_back( IndexType{ jspin, ispin, (int)i_pair } ); + } + }; + } + } +}; + +template<> +struct Functor::Local::DataRef +{ + using Interaction = Two_Site_Anisotropy; + using Data = typename Interaction::Data; + using Cache = typename Interaction::Cache; + + DataRef( const Data & data, const Cache & cache ) noexcept + : is_contributing( Interaction::is_contributing( data, cache ) ), matrices( cache.matrices.data() ) + { + } + + const bool is_contributing; + +protected: + const Matrix3 * matrices; +}; + +template<> +inline scalar Two_Site_Anisotropy::Energy::operator()( const Index & index, quantity state ) const +{ + // don't need to check for `is_contributing` here, because `transform_reduce` will short circuit properly + return Backend::transform_reduce( + index.begin(), index.end(), scalar( 0.0 ), Backend::plus{}, + [this, state] SPIRIT_LAMBDA( const Two_Site_Anisotropy::IndexType & idx ) -> scalar + { return state.spin[idx.ispin].dot( matrices[idx.ipair] * state.spin[idx.jspin] ); } ); +} + +template<> +inline Vector3 Two_Site_Anisotropy::Gradient::operator()( const Index & index, quantity state ) const +{ + // don't need to check for `is_contributing` here, because `transform_reduce` will short circuit properly + return Backend::transform_reduce( + index.begin(), index.end(), Vector3{ 0.0, 0.0, 0.0 }, Backend::plus{}, + [this, state] SPIRIT_LAMBDA( const Two_Site_Anisotropy::IndexType & idx ) -> Vector3 + { return matrices[idx.ipair] * state.spin[idx.jspin]; } ); +} + +template<> +template +void Two_Site_Anisotropy::Hessian::operator()( const Index & index, const StateType &, Callable & hessian ) const +{ + // don't need to check for `is_contributing` here, because `for_each` will short circuit properly + Backend::cpu::for_each( + index.begin(), index.end(), + [this, &index, &hessian]( const Two_Site_Anisotropy::IndexType & idx ) + { + const int i = 3 * idx.ispin; + const int j = 3 * idx.jspin; + + for( int alpha = 0; alpha < 3; ++alpha ) + { + for( int beta = 0; beta < 3; ++beta ) + { + hessian( i + alpha, j + beta, matrices[idx.ipair]( alpha, beta ) ); + } + } + } ); +}; + +} // namespace Interaction + +} // namespace Spin + +} // namespace Engine +#endif diff --git a/core/include/io/hamiltonian/Interactions.hpp b/core/include/io/hamiltonian/Interactions.hpp index 74118f784..9aceb941e 100644 --- a/core/include/io/hamiltonian/Interactions.hpp +++ b/core/include/io/hamiltonian/Interactions.hpp @@ -30,7 +30,7 @@ void Biaxial_Anisotropy_from_Config( void Pair_Interactions_from_Pairs_from_Config( const std::string & config_file_name, const Data::Geometry & geometry, std::vector & parameter_log, pairfield & exchange_pairs, scalarfield & exchange_magnitudes, pairfield & dmi_pairs, scalarfield & dmi_magnitudes, - vectorfield & dmi_normals ); + vectorfield & dmi_normals, Engine::Spin::Interaction::Two_Site_Anisotropy::Data & anisotropy ); void Pair_Interactions_from_Shells_from_Config( const std::string & config_file_name, const Data::Geometry & geometry, std::vector & parameter_log, diff --git a/core/src/io/hamiltonian/Hamiltonian.cpp b/core/src/io/hamiltonian/Hamiltonian.cpp index 1339893c6..6b7458ea5 100644 --- a/core/src/io/hamiltonian/Hamiltonian.cpp +++ b/core/src/io/hamiltonian/Hamiltonian.cpp @@ -86,6 +86,10 @@ std::unique_ptr Hamiltonian_Heisenberg_from_Co scalar ddi_radius = 0.0; bool ddi_pb_zero_padding = false; + Engine::Spin::Interaction::Two_Site_Anisotropy::Data two_site_anisotropy{}; + two_site_anisotropy.pairs = pairfield( 0 ); + two_site_anisotropy.matrices = matrixfield( 0 ); + // ------------ Quadruplet Interactions ------------ auto quadruplets = quadrupletfield( 0 ); auto quadruplet_magnitudes = scalarfield( 0 ); @@ -111,7 +115,7 @@ std::unique_ptr Hamiltonian_Heisenberg_from_Co else Pair_Interactions_from_Pairs_from_Config( config_file_name, geometry, parameter_log, exchange_pairs, exchange_magnitudes, dmi_pairs, - dmi_magnitudes, dmi_normals ); + dmi_magnitudes, dmi_normals, two_site_anisotropy ); DDI_from_Config( config_file_name, geometry, parameter_log, ddi_method, ddi_n_periodic_images, ddi_pb_zero_padding, @@ -154,7 +158,7 @@ std::unique_ptr Hamiltonian_Heisenberg_from_Co auto hamiltonian = std::make_unique( HamiltonianVariant::Heisenberg( std::move( geometry ), std::move( boundary_conditions ), zeeman, anisotropy, biaxial_anisotropy, - cubic_anisotropy, exchange, dmi, quadruplet, ddi ) ); + cubic_anisotropy, exchange, dmi, two_site_anisotropy, quadruplet, ddi ) ); assert( hamiltonian->Name() == "Heisenberg" ); Log( Log_Level::Debug, Log_Sender::IO, fmt::format( "Hamiltonian_{}: built", hamiltonian->Name() ) ); diff --git a/core/src/io/hamiltonian/Pairs.cpp b/core/src/io/hamiltonian/Pairs.cpp index df9677f6f..7fae6bc8b 100644 --- a/core/src/io/hamiltonian/Pairs.cpp +++ b/core/src/io/hamiltonian/Pairs.cpp @@ -3,6 +3,7 @@ #include #include +#include #include using Utility::Log_Level, Utility::Log_Sender; @@ -16,54 +17,107 @@ namespace // Read from Pairs file by Markus & Bernd void Pairs_from_File( const std::string & pairs_file, const Data::Geometry & geometry, int & nop, pairfield & exchange_pairs, - scalarfield & exchange_magnitudes, pairfield & dmi_pairs, scalarfield & dmi_magnitudes, - vectorfield & dmi_normals ) noexcept + scalarfield & exchange_magnitudes, pairfield & dmi_pairs, scalarfield & dmi_magnitudes, vectorfield & dmi_normals, + Engine::Spin::Interaction::Two_Site_Anisotropy::Data & anisotropy ) noexcept try { Log( Log_Level::Debug, Log_Sender::IO, fmt::format( "Reading spin pairs from file \"{}\"", pairs_file ) ); - using PairTableParser = TableParserInit, std::array>; - const PairTableParser parser( - { "i", "j", "da", "db", "dc", "dij", "dijx", "dijy", "dijz", "dija", "dijb", "dijc", "jij" } ); + using PairTableParser = TableParserInit, std::array>; + const PairTableParser parser( { "i", "j", "da", "db", "dc", "dij", "dijx", "dijy", "dijz", "dija", "dijb", "dijc", + "jij", "kijxx", "kijyy", "kijzz", "kijxy", "kijxz", "kijyz" } ); auto transform_factory = [&pairs_file, &geometry]( const std::map & idx ) { - bool DMI_xyz = false, DMI_abc = false, DMI_magnitude = false; + const bool DMI_xyz = idx.at( "dijx" ) >= 0 && idx.at( "dijy" ) >= 0 && idx.at( "dijz" ) >= 0; + const bool DMI_abc = idx.at( "dija" ) >= 0 && idx.at( "dijb" ) >= 0 && idx.at( "dijc" ) >= 0; + const bool DMI_magnitude = idx.at( "dij" ) >= 0; - if( idx.at( "dijx" ) >= 0 && idx.at( "dijy" ) >= 0 && idx.at( "dijz" ) >= 0 ) - DMI_xyz = true; - if( idx.at( "dija" ) >= 0 && idx.at( "dijb" ) >= 0 && idx.at( "dijc" ) >= 0 ) - DMI_abc = true; - if( idx.at( "dij" ) >= 0 ) - DMI_magnitude = true; - - if( idx.at( "j" ) < 0 && !DMI_xyz && !DMI_abc ) + const std::bitset<3> Kdiag = [&idx] + { + std::bitset<3> flags{}; + flags[0] = idx.at( "kijxx" ) >= 0; + flags[1] = idx.at( "kijyy" ) >= 0; + flags[2] = idx.at( "kijzz" ) >= 0; + return flags; + }(); + + if( idx.at( "j" ) < 0 && !DMI_xyz && !DMI_abc && Kdiag.count() == 1 ) Log( Log_Level::Warning, Log_Sender::IO, fmt::format( "No interactions could be found in pairs file \"{}\"", pairs_file ) ); - return [DMI_xyz, DMI_abc, DMI_magnitude, - &geometry]( const PairTableParser::read_row_t & row ) -> std::tuple + if( Kdiag.count() == 3 ) + Log( Log_Level::Warning, Log_Sender::IO, + fmt::format( + "Found three diagonal elements for two-site anisotropy, adjusting Heisenberg exchange!" + "file: \"{}\"! ", + pairs_file ) ); + + return [DMI_xyz, DMI_abc, DMI_magnitude, Kdiag, &geometry]( const PairTableParser::read_row_t & row ) + -> std::tuple> { - auto [i, j, da, db, dc, Dij, Dijx, Dijy, Dijz, Dija, Dijb, Dijc, Jij] = row; + struct RowData + { + int i, j, da, db, dc; + scalar Dij, Dijx, Dijy, Dijz, Dija, Dijb, Dijc, Jij; + scalar Kijxx, Kijyy, Kijzz, Kijxy, Kijxz, Kijyz; + }; + + auto data = IO::make_from_tuple( row ); Vector3 D_temp = Vector3::Zero(); if( DMI_xyz ) - D_temp = { Dijx, Dijy, Dijz }; + D_temp = { data.Dijx, data.Dijy, data.Dijz }; // Anisotropy vector orientation if( DMI_abc ) { - D_temp = { Dija, Dijb, Dijc }; + D_temp = { data.Dija, data.Dijb, data.Dijc }; D_temp = { D_temp.dot( geometry.lattice_constant * geometry.bravais_vectors[0] ), D_temp.dot( geometry.lattice_constant * geometry.bravais_vectors[1] ), D_temp.dot( geometry.lattice_constant * geometry.bravais_vectors[2] ) }; } if( !DMI_magnitude ) - Dij = D_temp.norm(); + data.Dij = D_temp.norm(); D_temp.normalize(); - return std::make_tuple( Pair{ i, j, { da, db, dc } }, Jij, D_temp, Dij ); + const auto Kij = [&d = data, &Kdiag] + { + std::optional Kij = std::nullopt; + if( Kdiag.count() < 2 ) + return Kij; + else if( d.Kijxx == 0 && d.Kijxy == 0 && d.Kijxz == 0 && d.Kijyy == 0 && d.Kijyz == 0 && d.Kijzz == 0 ) + return Kij; + else if( Kdiag.count() == 2 ) + { + if( !Kdiag[0] ) + d.Kijxx = -( d.Kijyy + d.Kijzz ); + else if( !Kdiag[1] ) + d.Kijyy = -( d.Kijxx + d.Kijzz ); + else + d.Kijzz = -( d.Kijxx + d.Kijyy ); + } + else + { + const scalar tr = ( d.Kijxx + d.Kijyy + d.Kijzz ) / 3.0; + d.Kijxx -= tr; + d.Kijyy -= tr; + d.Kijyy -= tr; + d.Jij += tr; + } + + Kij.emplace( Matrix3{} ); + // clang-format off + (*Kij) << d.Kijxx, d.Kijxy, d.Kijxz, + d.Kijxy, d.Kijyy, d.Kijyz, + d.Kijxz, d.Kijyz, d.Kijzz; + // clang-format on + return Kij; + }(); + + return std::make_tuple( + Pair{ data.i, data.j, { data.da, data.db, data.dc } }, data.Jij, D_temp, data.Dij, Kij ); }; }; @@ -85,7 +139,7 @@ try }; // Add the indices and parameters to the corresponding lists and deduplicate entries - for( const auto & [pair, Jij, D_vec, Dij] : data ) + for( const auto & [pair, Jij, D_vec, Dij, Kij] : data ) { if( Jij != 0 ) { @@ -141,13 +195,37 @@ try dmi_normals.push_back( D_vec ); } } + if( Kij.has_value() ) + { + bool already_in{ false }; + int atposition = -1; + for( std::size_t icheck = 0; icheck < anisotropy.pairs.size(); ++icheck ) + { + if( predicate( pair, anisotropy.pairs[icheck] ) == 0 ) + continue; + + already_in = true; + atposition = icheck; + break; + } + if( already_in ) + { + anisotropy.matrices[atposition] += *Kij; + } + else + { + anisotropy.pairs.push_back( pair ); + anisotropy.matrices.push_back( *Kij ); + } + } } } Log( Log_Level::Parameter, Log_Sender::IO, fmt::format( - "Done reading {} spin pairs from file \"{}\", giving {} exchange and {} DM (symmetry-reduced) pairs.", nop, - pairs_file, exchange_pairs.size(), dmi_pairs.size() ) ); + "Done reading {} spin pairs from file \"{}\", giving {} exchange, {} DM (symmetry-reduced), and {} " + "two-site anisotropy pairs.", + nop, pairs_file, exchange_pairs.size(), dmi_pairs.size(), anisotropy.pairs.size() ) ); } catch( ... ) { @@ -159,7 +237,7 @@ catch( ... ) void Pair_Interactions_from_Pairs_from_Config( const std::string & config_file_name, const Data::Geometry & geometry, std::vector & parameter_log, pairfield & exchange_pairs, scalarfield & exchange_magnitudes, pairfield & dmi_pairs, scalarfield & dmi_magnitudes, - vectorfield & dmi_normals ) + vectorfield & dmi_normals, Engine::Spin::Interaction::Two_Site_Anisotropy::Data & anisotropy ) { std::string interaction_pairs_file{}; int n_pairs = 0; @@ -179,7 +257,7 @@ void Pair_Interactions_from_Pairs_from_Config( // The file name should be valid so we try to read it Pairs_from_File( interaction_pairs_file, geometry, n_pairs, exchange_pairs, exchange_magnitudes, dmi_pairs, - dmi_magnitudes, dmi_normals ); + dmi_magnitudes, dmi_normals, anisotropy ); } // else //{ @@ -257,6 +335,9 @@ void Pair_Interactions_from_Shells_from_Config( fmt::format( "Failed to read DMI parameters from config file \"{}\"", config_file_name ) ); } + Log( Log_Level::Warning, Log_Sender::IO, + "Hamiltonian_Heisenberg: two-site anisotropy not supported when using neighbours." ); + parameter_log.emplace_back( fmt::format( " {:<21} = {}", "n_shells_exchange", n_shells_exchange ) ); if( n_shells_exchange > 0 ) parameter_log.emplace_back( fmt::format( " {:<21} = {}", "J_ij[0]", exchange_magnitudes[0] ) );