diff --git a/core/docs/Input.md b/core/docs/Input.md index cc2d11f9e..cdba042ee 100644 --- a/core/docs/Input.md +++ b/core/docs/Input.md @@ -261,6 +261,11 @@ i j da db dc Jij Dij Dijx Dijy Dijz 0 0 0 1 0 10.0 6.0 0.0 1.0 0.0 0 0 0 0 1 10.0 6.0 0.0 0.0 1.0 +### Triplets +n_interaction_triplets 1 +i j da_j db_j dc_j k da_k db_k dc_k na nb nc Q1 Q2 +0 0 1 0 0 0 0 1 0 0 0 1 3.0 4.0 + ### Quadruplets n_interaction_quadruplets 1 i j da_j db_j dc_j k da_k db_k dc_k l da_l db_l dc_l Q @@ -277,20 +282,27 @@ Note that instead of specifying the DM-vector as `Dijx Dijy Dijz`, you may speci `Dija Dijb Dijc` if you prefer. You may also specify the magnitude separately as a column `Dij`, but note that if you do, the vector (e.g. `Dijx Dijy Dijz`) will be normalized. -*Quadruplets:* Columns for these may also be placed in arbitrary order. +*Triplets:* +Columns for these may also be placed in arbitrary order. + +*Quadruplets:* +Columns for these may also be placed in arbitrary order. *Separate files:* -The anisotropy, pairs and quadruplets can be placed into separate files, -you can use `anisotropy_from_file`, `pairs_from_file` and `quadruplets_from_file`. +The anisotropy, pairs, triplets and quadruplets can be placed into separate files, +you can use `anisotropy_from_file`, `pairs_from_file`, `triplets_from_file` and `quadruplets_from_file`. -If the headers for anisotropies, pairs or quadruplets are at the top of the respective file, -it is not necessary to specify `n_anisotropy`, `n_interaction_pairs` or `n_interaction_quadruplets` +If the headers for anisotropies, pairs, triplets or quadruplets are at the top of the respective file, +it is not necessary to specify `n_anisotropy`, `n_interaction_pairs`, `n_interaction_triplets` or `n_interaction_quadruplets` respectively. ```Python ### Pairs interaction_pairs_file input/pairs.txt +### Triplets +interaction_triplets_file input/triplets.txt + ### Quadruplets interaction_quadruplets_file input/quadruplets.txt ``` @@ -557,4 +569,4 @@ inside the file. --- -[Home](Readme.md) \ No newline at end of file +[Home](Readme.md) diff --git a/core/include/engine/Hamiltonian_Heisenberg.hpp b/core/include/engine/Hamiltonian_Heisenberg.hpp index 59a71f51f..76aee1e97 100644 --- a/core/include/engine/Hamiltonian_Heisenberg.hpp +++ b/core/include/engine/Hamiltonian_Heisenberg.hpp @@ -36,6 +36,7 @@ namespace Engine pairfield exchange_pairs, scalarfield exchange_magnitudes, pairfield dmi_pairs, scalarfield dmi_magnitudes, vectorfield dmi_normals, DDI_Method ddi_method, intfield ddi_n_periodic_images, bool ddi_pb_zero_padding, scalar ddi_radius, + tripletfield triplets, scalarfield triplet_magnitudes1, scalarfield triplet_magnitudes2, quadrupletfield quadruplets, scalarfield quadruplet_magnitudes, std::shared_ptr geometry, intfield boundary_conditions @@ -47,6 +48,7 @@ namespace Engine scalarfield exchange_shell_magnitudes, scalarfield dmi_shell_magnitudes, int dm_chirality, DDI_Method ddi_method, intfield ddi_n_periodic_images, bool ddi_pb_zero_padding, scalar ddi_radius, + tripletfield triplets, scalarfield triplet_magnitudes1, scalarfield triplet_magnitudes2, quadrupletfield quadruplets, scalarfield quadruplet_magnitudes, std::shared_ptr geometry, intfield boundary_conditions @@ -70,7 +72,7 @@ namespace Engine // Hamiltonian name as string const std::string& Name() override; - + // ------------ Single Spin Interactions ------------ // External magnetic field across the sample scalar external_field_magnitude; @@ -111,6 +113,10 @@ namespace Engine scalarfield ddi_magnitudes; vectorfield ddi_normals; + // ------------ Triplet Interactions ------------ + tripletfield triplets; + scalarfield triplet_magnitudes1, triplet_magnitudes2; + // ------------ Quadruplet Interactions ------------ quadrupletfield quadruplets; scalarfield quadruplet_magnitudes; @@ -129,6 +135,8 @@ namespace Engine // Calculates the Dipole-Dipole contribution to the effective field of spin ispin within system s void Gradient_DDI(const vectorfield& spins, vectorfield & gradient); + // Triplet + void Gradient_Triplet(const vectorfield & spins, vectorfield & gradient); // Quadruplet void Gradient_Quadruplet(const vectorfield & spins, vectorfield & gradient); @@ -139,6 +147,7 @@ namespace Engine inline int Idx_Exchange() {return idx_exchange;}; inline int Idx_DMI() {return idx_dmi;}; inline int Idx_DDI() {return idx_ddi;}; + inline int Idx_Triplet() {return idx_triplet;}; inline int Idx_Quadruplet() {return idx_quadruplet;}; // Calculate the Zeeman energy of a Spin System @@ -151,11 +160,13 @@ namespace Engine void E_DMI(const vectorfield & spins, scalarfield & Energy); // Calculate the Dipole-Dipole energy void E_DDI(const vectorfield& spins, scalarfield & Energy); + // Calculate the Triplet energy + void E_Triplet(const vectorfield & spins, scalarfield & Energy); // Calculate the Quadruplet energy void E_Quadruplet(const vectorfield & spins, scalarfield & Energy); private: - int idx_zeeman, idx_anisotropy, idx_exchange, idx_dmi, idx_ddi, idx_quadruplet; + int idx_zeeman, idx_anisotropy, idx_exchange, idx_dmi, idx_ddi, idx_triplet, idx_quadruplet; void Gradient_DDI_Cutoff(const vectorfield& spins, vectorfield & gradient); void Gradient_DDI_Direct(const vectorfield& spins, vectorfield & gradient); void Gradient_DDI_FFT(const vectorfield& spins, vectorfield & gradient); @@ -203,4 +214,4 @@ namespace Engine } -#endif \ No newline at end of file +#endif diff --git a/core/include/engine/Vectormath_Defines.hpp b/core/include/engine/Vectormath_Defines.hpp index e8db08bf8..29faf7e79 100644 --- a/core/include/engine/Vectormath_Defines.hpp +++ b/core/include/engine/Vectormath_Defines.hpp @@ -50,6 +50,7 @@ using Vector2 = Eigen::Matrix; { int i, j, k; int d_j[3], d_k[3]; + scalar n[3]; }; struct Quadruplet { @@ -77,6 +78,7 @@ using Vector2 = Eigen::Matrix; { int i, j, k; std::array d_j, d_k; + std::array n; }; struct Quadruplet { diff --git a/core/include/io/Dataparser.hpp b/core/include/io/Dataparser.hpp index 460fe9375..992ed9454 100644 --- a/core/include/io/Dataparser.hpp +++ b/core/include/io/Dataparser.hpp @@ -1,4 +1,3 @@ - #pragma once #ifndef SPIRIT_IO_DATAPARSER_HPP #define SPIRIT_IO_DATAPARSER_HPP @@ -29,6 +28,10 @@ void Pairs_from_File( scalarfield & exchange_magnitudes, pairfield & dmi_pairs, scalarfield & dmi_magnitudes, vectorfield & dmi_normals ) noexcept; +void Triplets_from_File( + const std::string tripletsFile, const std::shared_ptr geometry, int & noq, tripletfield & triplets, + scalarfield & triplet_magnitudes, scalarfield & triplet_magnitudes2 ); + void Quadruplets_from_File( const std::string quadrupletsFile, const std::shared_ptr geometry, int & noq, quadrupletfield & quadruplets, scalarfield & quadruplet_magnitudes ) noexcept; diff --git a/core/src/engine/Hamiltonian_Heisenberg.cpp b/core/src/engine/Hamiltonian_Heisenberg.cpp index cc85a23ee..148ac1b6a 100644 --- a/core/src/engine/Hamiltonian_Heisenberg.cpp +++ b/core/src/engine/Hamiltonian_Heisenberg.cpp @@ -28,6 +28,7 @@ namespace Engine pairfield exchange_pairs, scalarfield exchange_magnitudes, pairfield dmi_pairs, scalarfield dmi_magnitudes, vectorfield dmi_normals, DDI_Method ddi_method, intfield ddi_n_periodic_images, bool ddi_pb_zero_padding, scalar ddi_radius, + tripletfield triplets, scalarfield triplet_magnitudes1, scalarfield triplet_magnitudes2, quadrupletfield quadruplets, scalarfield quadruplet_magnitudes, std::shared_ptr geometry, intfield boundary_conditions @@ -38,6 +39,7 @@ namespace Engine anisotropy_indices(anisotropy_indices), anisotropy_magnitudes(anisotropy_magnitudes), anisotropy_normals(anisotropy_normals), exchange_pairs_in(exchange_pairs), exchange_magnitudes_in(exchange_magnitudes), exchange_shell_magnitudes(0), dmi_pairs_in(dmi_pairs), dmi_magnitudes_in(dmi_magnitudes), dmi_normals_in(dmi_normals), dmi_shell_magnitudes(0), dmi_shell_chirality(0), + triplets(triplets), triplet_magnitudes1(triplet_magnitudes1), triplet_magnitudes2(triplet_magnitudes2), quadruplets(quadruplets), quadruplet_magnitudes(quadruplet_magnitudes), ddi_method(ddi_method), ddi_n_periodic_images(ddi_n_periodic_images), ddi_pb_zero_padding(ddi_pb_zero_padding) ,ddi_cutoff_radius(ddi_radius), fft_plan_reverse(FFT::FFT_Plan()), fft_plan_spins(FFT::FFT_Plan()) @@ -53,6 +55,7 @@ namespace Engine scalarfield exchange_shell_magnitudes, scalarfield dmi_shell_magnitudes, int dm_chirality, DDI_Method ddi_method, intfield ddi_n_periodic_images, bool ddi_pb_zero_padding, scalar ddi_radius, + tripletfield triplets, scalarfield triplet_magnitudes1, scalarfield triplet_magnitudes2, quadrupletfield quadruplets, scalarfield quadruplet_magnitudes, std::shared_ptr geometry, intfield boundary_conditions @@ -63,6 +66,7 @@ namespace Engine anisotropy_indices(anisotropy_indices), anisotropy_magnitudes(anisotropy_magnitudes), anisotropy_normals(anisotropy_normals), exchange_pairs_in(0), exchange_magnitudes_in(0), exchange_shell_magnitudes(exchange_shell_magnitudes), dmi_pairs_in(0), dmi_magnitudes_in(0), dmi_normals_in(0), dmi_shell_magnitudes(dmi_shell_magnitudes), dmi_shell_chirality(dm_chirality), + triplets(triplets), triplet_magnitudes1(triplet_magnitudes1), triplet_magnitudes2(triplet_magnitudes2), quadruplets(quadruplets), quadruplet_magnitudes(quadruplet_magnitudes), ddi_method(ddi_method), ddi_n_periodic_images(ddi_n_periodic_images), ddi_pb_zero_padding(ddi_pb_zero_padding), ddi_cutoff_radius(ddi_radius), fft_plan_reverse(FFT::FFT_Plan()), fft_plan_spins(FFT::FFT_Plan()) @@ -206,6 +210,13 @@ namespace Engine this->idx_ddi = this->energy_contributions_per_spin.size()-1; } else this->idx_ddi = -1; + // Triplets + if (this->triplets.size() > 0) + { + this->energy_contributions_per_spin.push_back({"triplets", scalarfield(0) }); + this->idx_triplet = this->energy_contributions_per_spin.size()-1; + } + else this->idx_triplet = -1; // Quadruplets if( this->quadruplets.size() > 0 ) { @@ -243,6 +254,8 @@ namespace Engine if( this->idx_dmi >=0 ) E_DMI(spins,contributions[idx_dmi].second); // DDI if( this->idx_ddi >=0 ) E_DDI(spins, contributions[idx_ddi].second); + // Triplets + if (this->idx_triplet >=0 ) E_Triplet(spins, contributions[idx_triplet].second); // Quadruplets if (this->idx_quadruplet >=0 ) E_Quadruplet(spins, contributions[idx_quadruplet].second); } @@ -420,6 +433,40 @@ namespace Engine } } + void Hamiltonian_Heisenberg::E_Triplet(const vectorfield & spins, scalarfield & Energy) + { + for (unsigned int itrip = 0; itrip < triplets.size(); ++itrip) + { + for (int da = 0; da < geometry->n_cells[0]; ++da) + { + for (int db = 0; db < geometry->n_cells[1]; ++db) + { + for (int dc = 0; dc < geometry->n_cells[2]; ++dc) + { + std::array translations = { da, db, dc }; + int ispin = triplets[itrip].i + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations); + int jspin = triplets[itrip].j + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, triplets[itrip].d_j); + int kspin = triplets[itrip].k + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, triplets[itrip].d_k); + Vector3 n = {triplets[itrip].n[0], triplets[itrip].n[1], triplets[itrip].n[2]}; + if ( check_atom_type(this->geometry->atom_types[ispin]) && check_atom_type(this->geometry->atom_types[jspin]) && + check_atom_type(this->geometry->atom_types[kspin])) + { + Energy[ispin] -= 1.0/3.0 * triplet_magnitudes1[itrip] * pow(spins[ispin].dot(spins[jspin].cross(spins[kspin])),2); + Energy[jspin] -= 1.0/3.0 * triplet_magnitudes1[itrip] * pow(spins[ispin].dot(spins[jspin].cross(spins[kspin])),2); + Energy[kspin] -= 1.0/3.0 * triplet_magnitudes1[itrip] * pow(spins[ispin].dot(spins[jspin].cross(spins[kspin])),2); + Energy[ispin] -= 1.0/3.0 * triplet_magnitudes2[itrip] * spins[ispin].dot(spins[jspin].cross(spins[kspin])) + * (n.dot(spins[ispin]+spins[jspin]+spins[kspin])); + Energy[jspin] -= 1.0/3.0 * triplet_magnitudes2[itrip] * spins[ispin].dot(spins[jspin].cross(spins[kspin])) + * (n.dot(spins[ispin]+spins[jspin]+spins[kspin])); + Energy[kspin] -= 1.0/3.0 * triplet_magnitudes2[itrip] * spins[ispin].dot(spins[jspin].cross(spins[kspin])) + * (n.dot(spins[ispin]+spins[jspin]+spins[kspin])); + } + } + } + } + } + } + void Hamiltonian_Heisenberg::E_Quadruplet(const vectorfield & spins, scalarfield & Energy) { for( unsigned int iquad = 0; iquad < quadruplets.size(); ++iquad ) @@ -532,8 +579,42 @@ namespace Engine } } + // Triplets + if (this->idx_triplet >= 0) + { + for (unsigned int itrip = 0; itrip < quadruplets.size(); ++itrip) + { + auto translations = Vectormath::translations_from_idx(geometry->n_cells, geometry->n_cell_atoms, icell); + int ispin = quadruplets[itrip].i + icell*geometry->n_cell_atoms; + int jspin = quadruplets[itrip].j + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, quadruplets[itrip].d_j); + int kspin = quadruplets[itrip].k + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, quadruplets[itrip].d_k); + Vector3 n = {triplets[itrip].n[0], triplets[itrip].n[1], triplets[itrip].n[2]}; + + if ( check_atom_type(this->geometry->atom_types[ispin]) && check_atom_type(this->geometry->atom_types[jspin]) && + check_atom_type(this->geometry->atom_types[kspin])) + { + Energy -= 1.0/3.0 * triplet_magnitudes1[itrip] * pow(spins[ispin].dot(spins[jspin].cross(spins[kspin])),2); + Energy -= 1.0/3.0 * triplet_magnitudes2[itrip] * spins[ispin].dot(spins[jspin].cross(spins[kspin])) + * (n.dot(spins[ispin]+spins[jspin]+spins[kspin])); + } + + #ifndef _OPENMP + // TODO: mirrored quadruplet when unique quadruplets are used + // jspin = quadruplets[iquad].j + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, quadruplets[iquad].d_j, true); + // kspin = quadruplets[iquad].k + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, quadruplets[iquad].d_k, true); + // lspin = quadruplets[iquad].l + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, quadruplets[iquad].d_l, true); + + // if ( check_atom_type(this->geometry->atom_types[ispin]) && check_atom_type(this->geometry->atom_types[jspin]) && + // check_atom_type(this->geometry->atom_types[kspin]) && check_atom_type(this->geometry->atom_types[lspin]) ) + // { + // Energy -= 0.25*quadruplet_magnitudes[iquad] * (spins[ispin].dot(spins[jspin])) * (spins[kspin].dot(spins[lspin])); + // } + #endif + } + } + // TODO: Quadruplets - if( this->idx_quadruplet >= 0 ) + if (this->idx_quadruplet >= 0) { } } @@ -566,6 +647,9 @@ namespace Engine if(idx_ddi >= 0) this->Gradient_DDI(spins, gradient); + // Triplets + this->Gradient_Triplet(spins, gradient); + // Quadruplets if(idx_quadruplet >= 0) this->Gradient_Quadruplet(spins, gradient); @@ -884,6 +968,51 @@ namespace Engine } } + void Hamiltonian_Heisenberg::Gradient_Triplet(const vectorfield & spins, vectorfield & gradient) + { + for (unsigned int itrip = 0; itrip < triplets.size(); ++itrip) + { + int i = triplets[itrip].i; + int j = triplets[itrip].j; + int k = triplets[itrip].k; + Vector3 n = {triplets[itrip].n[0], triplets[itrip].n[1], triplets[itrip].n[2]}; + for (int da = 0; da < geometry->n_cells[0]; ++da) + { + for (int db = 0; db < geometry->n_cells[1]; ++db) + { + for (int dc = 0; dc < geometry->n_cells[2]; ++dc) + { + std::array translations = { da, db, dc }; + int ispin = i + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations); + int jspin = j + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, triplets[itrip].d_j); + int kspin = k + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, triplets[itrip].d_k); + + if ( check_atom_type(this->geometry->atom_types[ispin]) && check_atom_type(this->geometry->atom_types[jspin]) && + check_atom_type(this->geometry->atom_types[kspin])) + { + gradient[ispin] -= 2.0 * triplet_magnitudes1[itrip] * spins[ispin].dot(spins[jspin].cross(spins[kspin])) + * spins[jspin].cross(spins[kspin]); + gradient[jspin] -= 2.0 * triplet_magnitudes1[itrip] * spins[ispin].dot(spins[jspin].cross(spins[kspin])) + * spins[kspin].cross(spins[ispin]); + gradient[kspin] -= 2.0 * triplet_magnitudes1[itrip] * spins[ispin].dot(spins[jspin].cross(spins[kspin])) + * spins[ispin].cross(spins[jspin]); + + gradient[ispin] -= triplet_magnitudes2[itrip] * (n.dot(spins[ispin]+spins[jspin]+spins[kspin]) + * spins[jspin].cross(spins[kspin]) + + spins[ispin].dot(spins[jspin].cross(spins[kspin])) * n); + gradient[jspin] -= triplet_magnitudes2[itrip] * (n.dot(spins[ispin]+spins[jspin]+spins[kspin]) + * spins[kspin].cross(spins[ispin]) + + spins[ispin].dot(spins[jspin].cross(spins[kspin])) * n); + gradient[kspin] -= triplet_magnitudes2[itrip] * (n.dot(spins[ispin]+spins[jspin]+spins[kspin]) + * spins[ispin].cross(spins[jspin]) + + spins[ispin].dot(spins[jspin].cross(spins[kspin])) * n); + std::cout << gradient[ispin] << " \n" << gradient[jspin] << " \n" << gradient[kspin] << "\n" << std::endl; + } + } + } + } + } + } void Hamiltonian_Heisenberg::Gradient_Quadruplet(const vectorfield & spins, vectorfield & gradient) { @@ -1392,4 +1521,4 @@ namespace Engine const std::string& Hamiltonian_Heisenberg::Name() { return name; } } -#endif \ No newline at end of file +#endif diff --git a/core/src/engine/Hamiltonian_Heisenberg.cu b/core/src/engine/Hamiltonian_Heisenberg.cu index ed3153fd0..3656149fa 100644 --- a/core/src/engine/Hamiltonian_Heisenberg.cu +++ b/core/src/engine/Hamiltonian_Heisenberg.cu @@ -31,6 +31,7 @@ namespace Engine pairfield exchange_pairs, scalarfield exchange_magnitudes, pairfield dmi_pairs, scalarfield dmi_magnitudes, vectorfield dmi_normals, DDI_Method ddi_method, intfield ddi_n_periodic_images, bool ddi_pb_zero_padding, scalar ddi_radius, + tripletfield triplets, scalarfield triplet_magnitudes1, scalarfield triplet_magnitudes2, quadrupletfield quadruplets, scalarfield quadruplet_magnitudes, std::shared_ptr geometry, intfield boundary_conditions @@ -41,6 +42,7 @@ namespace Engine anisotropy_indices(anisotropy_indices), anisotropy_magnitudes(anisotropy_magnitudes), anisotropy_normals(anisotropy_normals), exchange_pairs_in(exchange_pairs), exchange_magnitudes_in(exchange_magnitudes), exchange_shell_magnitudes(0), dmi_pairs_in(dmi_pairs), dmi_magnitudes_in(dmi_magnitudes), dmi_normals_in(dmi_normals), dmi_shell_magnitudes(0), dmi_shell_chirality(0), + triplets(triplets), triplet_magnitudes1(triplet_magnitudes1), triplet_magnitudes2(triplet_magnitudes2), quadruplets(quadruplets), quadruplet_magnitudes(quadruplet_magnitudes), ddi_method(ddi_method), ddi_n_periodic_images(ddi_n_periodic_images), ddi_pb_zero_padding(ddi_pb_zero_padding), ddi_cutoff_radius(ddi_radius), fft_plan_reverse(FFT::FFT_Plan()), fft_plan_spins(FFT::FFT_Plan()) @@ -56,6 +58,7 @@ namespace Engine scalarfield exchange_shell_magnitudes, scalarfield dmi_shell_magnitudes, int dmi_shell_chirality, DDI_Method ddi_method, intfield ddi_n_periodic_images, bool ddi_pb_zero_padding, scalar ddi_radius, + tripletfield triplets, scalarfield triplet_magnitudes1, scalarfield triplet_magnitudes2, quadrupletfield quadruplets, scalarfield quadruplet_magnitudes, std::shared_ptr geometry, intfield boundary_conditions @@ -66,6 +69,7 @@ namespace Engine anisotropy_indices(anisotropy_indices), anisotropy_magnitudes(anisotropy_magnitudes), anisotropy_normals(anisotropy_normals), exchange_pairs_in(0), exchange_magnitudes_in(0), exchange_shell_magnitudes(exchange_shell_magnitudes), dmi_pairs_in(0), dmi_magnitudes_in(0), dmi_normals_in(0), dmi_shell_magnitudes(dmi_shell_magnitudes), dmi_shell_chirality(dmi_shell_chirality), + triplets(triplets), triplet_magnitudes1(triplet_magnitudes1), triplet_magnitudes2(triplet_magnitudes2), quadruplets(quadruplets), quadruplet_magnitudes(quadruplet_magnitudes), ddi_method(ddi_method), ddi_n_periodic_images(ddi_n_periodic_images), ddi_pb_zero_padding(ddi_pb_zero_padding), ddi_cutoff_radius(ddi_radius), fft_plan_reverse(FFT::FFT_Plan()), fft_plan_spins(FFT::FFT_Plan()) @@ -205,6 +209,13 @@ namespace Engine this->idx_ddi = this->energy_contributions_per_spin.size()-1; } else this->idx_ddi = -1; + // Triplets + if (this->triplets.size() > 0) + { + this->energy_contributions_per_spin.push_back({"triplets", scalarfield(0) }); + this->idx_triplet = this->energy_contributions_per_spin.size()-1; + } + else this->idx_triplet = -1; // Quadruplets if( this->quadruplets.size() > 0 ) { @@ -242,7 +253,8 @@ namespace Engine if( this->idx_dmi >=0 ) E_DMI(spins,contributions[idx_dmi].second); // DDI if( this->idx_ddi >=0 ) E_DDI(spins, contributions[idx_ddi].second); - + // Triplets + if (this->idx_triplet >=0 ) E_Triplet(spins, contributions[idx_triplet].second); // Quadruplets if (this->idx_quadruplet >=0 ) E_Quadruplet(spins, contributions[idx_quadruplet].second); } @@ -511,6 +523,40 @@ namespace Engine } } + void Hamiltonian_Heisenberg::E_Triplet(const vectorfield & spins, scalarfield & Energy) + { + for (unsigned int itrip = 0; itrip < triplets.size(); ++itrip) + { + for (int da = 0; da < geometry->n_cells[0]; ++da) + { + for (int db = 0; db < geometry->n_cells[1]; ++db) + { + for (int dc = 0; dc < geometry->n_cells[2]; ++dc) + { + std::array translations = { da, db, dc }; + int ispin = triplets[itrip].i + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations); + int jspin = triplets[itrip].j + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, triplets[itrip].d_j); + int kspin = triplets[itrip].k + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, triplets[itrip].d_k); + Vector3 n = {triplets[itrip].n[0], triplets[itrip].n[1], triplets[itrip].n[2]}; + if ( check_atom_type(this->geometry->atom_types[ispin]) && check_atom_type(this->geometry->atom_types[jspin]) && + check_atom_type(this->geometry->atom_types[kspin])) + { + Energy[ispin] -= 1.0/3.0 * triplet_magnitudes1[itrip] * pow(spins[ispin].dot(spins[jspin].cross(spins[kspin])),2); + Energy[jspin] -= 1.0/3.0 * triplet_magnitudes1[itrip] * pow(spins[ispin].dot(spins[jspin].cross(spins[kspin])),2); + Energy[kspin] -= 1.0/3.0 * triplet_magnitudes1[itrip] * pow(spins[ispin].dot(spins[jspin].cross(spins[kspin])),2); + Energy[ispin] -= 1.0/3.0 * triplet_magnitudes2[itrip] * spins[ispin].dot(spins[jspin].cross(spins[kspin])) + * (n.dot(spins[ispin]+spins[jspin]+spins[kspin])); + Energy[jspin] -= 1.0/3.0 * triplet_magnitudes2[itrip] * spins[ispin].dot(spins[jspin].cross(spins[kspin])) + * (n.dot(spins[ispin]+spins[jspin]+spins[kspin])); + Energy[kspin] -= 1.0/3.0 * triplet_magnitudes2[itrip] * spins[ispin].dot(spins[jspin].cross(spins[kspin])) + * (n.dot(spins[ispin]+spins[jspin]+spins[kspin])); + } + } + } + } + } + } + void Hamiltonian_Heisenberg::E_Quadruplet(const vectorfield & spins, scalarfield & Energy) { for( unsigned int iquad = 0; iquad < quadruplets.size(); ++iquad ) @@ -608,6 +654,40 @@ namespace Engine } } + // Triplets + if (this->idx_triplet >= 0) + { + for (unsigned int itrip = 0; itrip < quadruplets.size(); ++itrip) + { + auto translations = Vectormath::translations_from_idx(geometry->n_cells, geometry->n_cell_atoms, icell); + int ispin = quadruplets[itrip].i + icell*geometry->n_cell_atoms; + int jspin = quadruplets[itrip].j + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, quadruplets[itrip].d_j); + int kspin = quadruplets[itrip].k + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, quadruplets[itrip].d_k); + Vector3 n = {triplets[itrip].n[0], triplets[itrip].n[1], triplets[itrip].n[2]}; + + if ( check_atom_type(this->geometry->atom_types[ispin]) && check_atom_type(this->geometry->atom_types[jspin]) && + check_atom_type(this->geometry->atom_types[kspin])) + { + Energy -= 1.0/3.0 * triplet_magnitudes1[itrip] * pow(spins[ispin].dot(spins[jspin].cross(spins[kspin])),2); + Energy -= 1.0/3.0 * triplet_magnitudes2[itrip] * spins[ispin].dot(spins[jspin].cross(spins[kspin])) + * (n.dot(spins[ispin]+spins[jspin]+spins[kspin])); + } + + #ifndef _OPENMP + // TODO: mirrored quadruplet when unique quadruplets are used + // jspin = quadruplets[iquad].j + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, quadruplets[iquad].d_j, true); + // kspin = quadruplets[iquad].k + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, quadruplets[iquad].d_k, true); + // lspin = quadruplets[iquad].l + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, quadruplets[iquad].d_l, true); + + // if ( check_atom_type(this->geometry->atom_types[ispin]) && check_atom_type(this->geometry->atom_types[jspin]) && + // check_atom_type(this->geometry->atom_types[kspin]) && check_atom_type(this->geometry->atom_types[lspin]) ) + // { + // Energy -= 0.25*quadruplet_magnitudes[iquad] * (spins[ispin].dot(spins[jspin])) * (spins[kspin].dot(spins[lspin])); + // } + #endif + } + } + // TODO: Quadruplets if (this->idx_quadruplet >= 0) { @@ -642,6 +722,9 @@ namespace Engine if(idx_ddi >= 0) this->Gradient_DDI(spins, gradient); + // Triplets + this->Gradient_Triplet(spins, gradient); + // Quadruplet if(idx_quadruplet) this->Gradient_Quadruplet(spins, gradient); @@ -928,6 +1011,52 @@ namespace Engine geometry->mu_s.data(), sublattice_size ); } // end Field_DipoleDipole + void Hamiltonian_Heisenberg::Gradient_Triplet(const vectorfield & spins, vectorfield & gradient) + { + for (unsigned int itrip = 0; itrip < triplets.size(); ++itrip) + { + int i = triplets[itrip].i; + int j = triplets[itrip].j; + int k = triplets[itrip].k; + Vector3 n = {triplets[itrip].n[0], triplets[itrip].n[1], triplets[itrip].n[2]}; + for (int da = 0; da < geometry->n_cells[0]; ++da) + { + for (int db = 0; db < geometry->n_cells[1]; ++db) + { + for (int dc = 0; dc < geometry->n_cells[2]; ++dc) + { + std::array translations = { da, db, dc }; + int ispin = i + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations); + int jspin = j + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, triplets[itrip].d_j); + int kspin = k + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations, triplets[itrip].d_k); + + if ( check_atom_type(this->geometry->atom_types[ispin]) && check_atom_type(this->geometry->atom_types[jspin]) && + check_atom_type(this->geometry->atom_types[kspin])) + { + gradient[ispin] -= 2.0 * triplet_magnitudes1[itrip] * spins[ispin].dot(spins[jspin].cross(spins[kspin])) + * spins[jspin].cross(spins[kspin]); + gradient[jspin] -= 2.0 * triplet_magnitudes1[itrip] * spins[ispin].dot(spins[jspin].cross(spins[kspin])) + * spins[kspin].cross(spins[ispin]); + gradient[kspin] -= 2.0 * triplet_magnitudes1[itrip] * spins[ispin].dot(spins[jspin].cross(spins[kspin])) + * spins[ispin].cross(spins[jspin]); + + gradient[ispin] -= triplet_magnitudes2[itrip] * (n.dot(spins[ispin]+spins[jspin]+spins[kspin]) + * spins[jspin].cross(spins[kspin]) + + spins[ispin].dot(spins[jspin].cross(spins[kspin])) * n); + gradient[jspin] -= triplet_magnitudes2[itrip] * (n.dot(spins[ispin]+spins[jspin]+spins[kspin]) + * spins[kspin].cross(spins[ispin]) + + spins[ispin].dot(spins[jspin].cross(spins[kspin])) * n); + gradient[kspin] -= triplet_magnitudes2[itrip] * (n.dot(spins[ispin]+spins[jspin]+spins[kspin]) + * spins[ispin].cross(spins[jspin]) + + spins[ispin].dot(spins[jspin].cross(spins[kspin])) * n); + std::cout << gradient[ispin] << " \n" << gradient[jspin] << " \n" << gradient[kspin] << "\n" << std::endl; + } + } + } + } + } + } + void Hamiltonian_Heisenberg::Gradient_Quadruplet(const vectorfield & spins, vectorfield & gradient) { for( unsigned int iquad = 0; iquad < quadruplets.size(); ++iquad ) diff --git a/core/src/io/Configparser.cpp b/core/src/io/Configparser.cpp index bca3b7a70..2f0f5aad8 100644 --- a/core/src/io/Configparser.cpp +++ b/core/src/io/Configparser.cpp @@ -653,7 +653,7 @@ Data::Pinning Pinning_from_Config( const std::string configFile, int n_cell_atom } Log( Log_Level::Parameter, Log_Sender::IO, "Pinning: read" ); return pinning; -#else // SPIRIT_ENABLE_PINNING +#else // SPIRIT_ENABLE_PINNING Log( Log_Level::Parameter, Log_Sender::IO, "Pinning is disabled" ); if( configFile != "" ) { @@ -1226,6 +1226,14 @@ std::unique_ptr Hamiltonian_Heisenberg_from_Conf scalar ddi_radius = 0.0; bool ddi_pb_zero_padding = true; + // ------------ Triplet Interactions ------------ + int n_triplets = 0; + std::string triplets_file = ""; + bool triplets_from_file = false; + tripletfield triplets( 0 ); + scalarfield triplet_magnitudes1( 0 ); + scalarfield triplet_magnitudes2( 0 ); + // ------------ Quadruplet Interactions ------------ int n_quadruplets = 0; std::string quadruplets_file = ""; @@ -1248,7 +1256,7 @@ std::unique_ptr Hamiltonian_Heisenberg_from_Conf boundary_conditions[0] = ( boundary_conditions_i[0] != 0 ); boundary_conditions[1] = ( boundary_conditions_i[1] != 0 ); boundary_conditions[2] = ( boundary_conditions_i[2] != 0 ); - } // end try + } catch( ... ) { spirit_handle_exception_core( @@ -1269,7 +1277,7 @@ std::unique_ptr Hamiltonian_Heisenberg_from_Conf Log( Log_Level::Warning, Log_Sender::IO, "Input for 'external_field_normal' had norm zero and has been set to (0,0,1)" ); } - } // end try + } catch( ... ) { spirit_handle_exception_core( @@ -1327,7 +1335,7 @@ std::unique_ptr Hamiltonian_Heisenberg_from_Conf anisotropy_normal = vectorfield( 0 ); } } - } // end try + } catch( ... ) { spirit_handle_exception_core( @@ -1359,7 +1367,7 @@ std::unique_ptr Hamiltonian_Heisenberg_from_Conf // been implemented yet."); throw Exception::System_not_Initialized; // // Not implemented! //} - } // end try + } catch( ... ) { spirit_handle_exception_core( @@ -1388,7 +1396,7 @@ std::unique_ptr Hamiltonian_Heisenberg_from_Conf "Hamiltonian_Heisenberg: Keyword 'jij' not found. Using Default: {}", exchange_magnitudes[0] ) ); } - } // end try + } catch( ... ) { spirit_handle_exception_core( @@ -1416,8 +1424,7 @@ std::unique_ptr Hamiltonian_Heisenberg_from_Conf dmi_magnitudes[0] ) ); } myfile.Read_Single( dm_chirality, "dm_chirality" ); - - } // end try + } catch( ... ) { spirit_handle_exception_core( @@ -1455,13 +1462,35 @@ std::unique_ptr Hamiltonian_Heisenberg_from_Conf // Dipole-dipole cutoff radius myfile.Read_Single( ddi_radius, "ddi_radius" ); - } // end try + } catch( ... ) { spirit_handle_exception_core( fmt::format( "Unable to read DDI radius from config file \"{}\"", configFile ) ); } + try + { + IO::Filter_File_Handle myfile( configFile ); + + // Interaction Triplets + if( myfile.Find( "n_interaction_triplets" ) ) + triplets_file = configFile; + else if( myfile.Find( "interaction_triplets_file" ) ) + myfile.iss >> triplets_file; + if( triplets_file.length() > 0 ) + { + // The file name should be valid so we try to read it + Triplets_from_File( + triplets_file, geometry, n_triplets, triplets, triplet_magnitudes1, triplet_magnitudes2 ); + } + } + catch( ... ) + { + spirit_handle_exception_core( + fmt::format( "Unable to read interaction triplets from config file \"{}\"", configFile ) ); + } + try { IO::Filter_File_Handle myfile( configFile ); @@ -1477,8 +1506,7 @@ std::unique_ptr Hamiltonian_Heisenberg_from_Conf // The file name should be valid so we try to read it Quadruplets_from_File( quadruplets_file, geometry, n_quadruplets, quadruplets, quadruplet_magnitudes ); } - - } // end try + } catch( ... ) { spirit_handle_exception_core( @@ -1524,15 +1552,17 @@ std::unique_ptr Hamiltonian_Heisenberg_from_Conf { hamiltonian = std::unique_ptr( new Engine::Hamiltonian_Heisenberg( B, B_normal, anisotropy_index, anisotropy_magnitude, anisotropy_normal, exchange_magnitudes, dmi_magnitudes, - dm_chirality, ddi_method, ddi_n_periodic_images, ddi_pb_zero_padding, ddi_radius, quadruplets, - quadruplet_magnitudes, geometry, boundary_conditions ) ); + dm_chirality, ddi_method, ddi_n_periodic_images, ddi_pb_zero_padding, ddi_radius, triplets, + triplet_magnitudes1, triplet_magnitudes2, quadruplets, quadruplet_magnitudes, geometry, + boundary_conditions ) ); } else { hamiltonian = std::unique_ptr( new Engine::Hamiltonian_Heisenberg( B, B_normal, anisotropy_index, anisotropy_magnitude, anisotropy_normal, exchange_pairs, exchange_magnitudes, dmi_pairs, dmi_magnitudes, dmi_normals, ddi_method, ddi_n_periodic_images, ddi_pb_zero_padding, ddi_radius, - quadruplets, quadruplet_magnitudes, geometry, boundary_conditions ) ); + triplets, triplet_magnitudes1, triplet_magnitudes2, quadruplets, quadruplet_magnitudes, geometry, + boundary_conditions ) ); } Log( Log_Level::Debug, Log_Sender::IO, "Hamiltonian_Heisenberg: built" ); return hamiltonian; diff --git a/core/src/io/Configwriter.cpp b/core/src/io/Configwriter.cpp index e7a8f454b..2ed83a4c8 100644 --- a/core/src/io/Configwriter.cpp +++ b/core/src/io/Configwriter.cpp @@ -314,6 +314,27 @@ void Hamiltonian_Heisenberg_to_Config( config += "### DDI cutoff radius (if cutoff is used)"; config += fmt::format( "ddi_radius {}\n", ham->ddi_cutoff_radius ); + // Triplets + config += "### Triplets:\n"; + config += fmt::format( "n_interaction_triplets {}\n", ham->quadruplets.size() ); + if( ham->quadruplets.size() > 0 ) + { + config += fmt::format( + "{:^3} {:^3} {:^3} {:^3} {:^3} {:^3} {:^3} {:^3} {:^3} {:^15} {:^15} {:^15} {:^15} {:^15}\n", + "i", "j", "k", "da_j", "db_j", "dc_j", "da_k", "db_k", "dc_k", "na", "nb", "nc", "q1", "q2" ); + for( unsigned int i = 0; i < ham->triplets.size(); ++i ) + { + config += fmt::format( + "{:^3} {:^3} {:^3} {:^3} {:^3} {:^3} {:^3} {:^3} {:^3} {:^15} {:^15} {:^15} {:^15.8f} " + "{:^15.8f}\n", + ham->triplets[i].i, ham->triplets[i].j, ham->triplets[i].k, ham->triplets[i].d_j[0], + ham->triplets[i].d_j[1], ham->triplets[i].d_j[2], ham->triplets[i].d_k[0], ham->triplets[i].d_k[1], + ham->triplets[i].d_k[2], ham->triplets[i].n[0], ham->triplets[i].n[1], ham->triplets[i].n[2], + ham->triplet_magnitudes1[i] ), + ham->triplet_magnitudes2[i]; + } + } + // Quadruplets config += "### Quadruplets:\n"; config += fmt::format( "n_interaction_quadruplets {}\n", ham->quadruplets.size() ); @@ -351,8 +372,8 @@ void Hamiltonian_Gaussian_to_Config( config += fmt::format( "{} {} {}\n", ham_gaussian->amplitude[i], ham_gaussian->width[i], ham_gaussian->center[i].transpose() ); } - } - Append_String_to_File( config, configFile ); +} +Append_String_to_File( config, configFile ); } } // namespace IO \ No newline at end of file diff --git a/core/src/io/Dataparser.cpp b/core/src/io/Dataparser.cpp index 45743dc80..213d8e099 100644 --- a/core/src/io/Dataparser.cpp +++ b/core/src/io/Dataparser.cpp @@ -468,6 +468,164 @@ catch( ... ) spirit_rethrow( fmt::format( "Could not read pairs file \"{}\"", pairsFile ) ); } +/* +Read from Triplet file +*/ +void Triplets_from_File( + const std::string tripletsFile, const std::shared_ptr, int & noq, tripletfield & triplets, + scalarfield & triplet_magnitudes1, scalarfield & triplet_magnitudes2 ) +{ + Log( Log_Level::Info, Log_Sender::IO, "Reading spin triplets from file " + tripletsFile ); + try + { + std::vector columns( 20 ); // at least: 4 (indices) + 3*3 (positions) + 1 (magnitude) + // column indices of pair indices and interactions + int col_i = -1; + int col_j = -1, col_da_j = -1, col_db_j = -1, col_dc_j = -1, periodicity_j = 0; + int col_k = -1, col_da_k = -1, col_db_k = -1, col_dc_k = -1, periodicity_k = 0; + int col_l = -1, col_da_l = -1, col_db_l = -1, col_dc_l = -1, periodicity_l = 0; + int col_Q1 = -1, col_Q2; + int col_na = -1, col_nb = -1, col_nc = -1; + bool Q1 = false; + bool Q2 = false; + int max_periods_a = 0, max_periods_b = 0, max_periods_c = 0; + int triplet_periodicity = 0; + int n_triplets = 0; + // Get column indices + Filter_File_Handle file( tripletsFile ); + if( file.Find( "n_interaction_triplets" ) ) + { + // Read n interaction triplets + file.iss >> n_triplets; + Log( Log_Level::Debug, Log_Sender::IO, + fmt::format( "File {} should have {} triplets", tripletsFile, n_triplets ) ); + } + else + { + // Read the whole file + n_triplets = (int)1e8; + // First line should contain the columns + file.ResetStream(); + Log( Log_Level::Info, Log_Sender::IO, "Trying to parse triplet columns from top of file " + tripletsFile ); + } + file.GetLine(); + for( unsigned int i = 0; i < columns.size(); ++i ) + { + file.iss >> columns[i]; + if( columns[i] == "i" ) + col_i = i; + else if( columns[i] == "j" ) + col_j = i; + else if( columns[i] == "da_j" ) + col_da_j = i; + else if( columns[i] == "db_j" ) + col_db_j = i; + else if( columns[i] == "dc_j" ) + col_dc_j = i; + else if( columns[i] == "k" ) + col_k = i; + else if( columns[i] == "da_k" ) + col_da_k = i; + else if( columns[i] == "db_k" ) + col_db_k = i; + else if( columns[i] == "dc_k" ) + col_dc_k = i; + else if( columns[i] == "q1" ) + { + col_Q1 = i; + Q1 = true; + } + else if( columns[i] == "na" ) + col_na = i; + else if( columns[i] == "nb" ) + col_nb = i; + else if( columns[i] == "nc" ) + col_nc = i; + else if( columns[i] == "q2" ) + { + col_Q2 = i; + Q2 = true; + } + } + // Check if interactions have been found in header + if( !Q1 ) + Log( Log_Level::Warning, Log_Sender::IO, + "No interactions could be found in header of triplets file " + tripletsFile ); + if( !Q2 ) + Log( Log_Level::Warning, Log_Sender::IO, + "No interactions could be found in header of triplets file " + tripletsFile ); + // triplet Indices + int q_i = 0; + int q_j = 0, q_da_j = 0, q_db_j = 0, q_dc_j = 0; + int q_k = 0, q_da_k = 0, q_db_k = 0, q_dc_k = 0; + scalar q_na = 0, q_nb = 0, q_nc = 0; + scalar q_Q1, q_Q2; + // Get actual triplets Data + int i_triplet = 0; + std::string sdump; + while( file.GetLine() && i_triplet < n_triplets ) + { + // Read a triplet from the File + for( unsigned int i = 0; i < columns.size(); ++i ) + { + // i + if( i == col_i ) + file.iss >> q_i; + // j + else if( i == col_j ) + file.iss >> q_j; + else if( i == col_da_j ) + file.iss >> q_da_j; + else if( i == col_db_j ) + file.iss >> q_db_j; + else if( i == col_dc_j ) + file.iss >> q_dc_j; + // k + else if( i == col_k ) + file.iss >> q_k; + else if( i == col_da_k ) + file.iss >> q_da_k; + else if( i == col_db_k ) + file.iss >> q_db_k; + else if( i == col_dc_k ) + file.iss >> q_dc_k; + // n + else if( i == col_na ) + file.iss >> q_na; + else if( i == col_nb ) + file.iss >> q_nb; + else if( i == col_nc ) + file.iss >> q_nc; + // triplet magnitude + else if( i == col_Q1 && Q1 ) + file.iss >> q_Q1; + else if( i == col_Q2 && Q2 ) + file.iss >> q_Q2; + // Otherwise dump the line + else + file.iss >> sdump; + } // end for columns + + // Add the indices and parameter to the corresponding list + if( q_Q1 != 0 || q_Q2 != 0 ) + { + triplets.push_back( + { q_i, q_j, q_k, { q_da_j, q_db_j, q_db_j }, { q_da_k, q_db_k, q_db_k }, { q_na, q_nb, q_nc } } ); + triplet_magnitudes1.push_back( q_Q1 ); + triplet_magnitudes2.push_back( q_Q2 ); + } + ++i_triplet; + } // end while GetLine + Log( Log_Level::Info, Log_Sender::IO, + fmt::format( "Done reading {} spin triplets from file {}", i_triplet, tripletsFile ) ); + noq = i_triplet; + } // end try + catch( ... ) + { + spirit_rethrow( fmt::format( "Could not read triplets from file \"{}\"", tripletsFile ) ); + } +} // End Triplets_from_File + // Read from Quadruplet file void Quadruplets_from_File( const std::string quadrupletsFile, const std::shared_ptr, int & noq, quadrupletfield & quadruplets,