From 764d87cceddd20f767b14567cb2bb65fbdb19034 Mon Sep 17 00:00:00 2001 From: Markus Hoffmann Date: Fri, 26 Oct 2018 16:18:52 +0200 Subject: [PATCH 1/2] added triplet interactions: ring-exchange and spin-chirality coupling --- core/docs/INPUT.md | 21 ++- .../include/engine/Hamiltonian_Heisenberg.hpp | 15 +- core/include/engine/Vectormath_Defines.hpp | 4 +- core/include/io/Dataparser.hpp | 3 + core/src/engine/Hamiltonian_Heisenberg.cpp | 132 +++++++++++++++++- core/src/io/Configparser.cpp | 30 ++++ core/src/io/Configwriter.cpp | 20 ++- core/src/io/Dataparser.cpp | 132 ++++++++++++++++++ 8 files changed, 347 insertions(+), 10 deletions(-) diff --git a/core/docs/INPUT.md b/core/docs/INPUT.md index e0e1f92db..2d9e1a915 100644 --- a/core/docs/INPUT.md +++ b/core/docs/INPUT.md @@ -237,6 +237,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 @@ -253,21 +258,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. +*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 ``` @@ -506,4 +517,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 982c4136e..46d20b1b5 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, 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, scalar ddi_radius, + tripletfield triplets, scalarfield triplet_magnitudes1, scalarfield triplet_magnitudes2, quadrupletfield quadruplets, scalarfield quadruplet_magnitudes, std::shared_ptr geometry, intfield boundary_conditions @@ -105,6 +107,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; @@ -127,12 +133,14 @@ namespace Engine void Gradient_DDI_Direct(const vectorfield& spins, vectorfield & gradient); void Gradient_DDI_FFT(const vectorfield& spins, vectorfield & gradient); + // Triplet + void Gradient_Triplet(const vectorfield & spins, vectorfield & gradient); // Quadruplet void Gradient_Quadruplet(const vectorfield & spins, vectorfield & gradient); // ------------ Energy Functions ------------ // Indices for Energy vector - 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; // Calculate the Zeeman energy of a Spin System void E_Zeeman(const vectorfield & spins, scalarfield & Energy); // Calculate the Anisotropy energy of a Spin System @@ -147,6 +155,9 @@ namespace Engine void E_DDI_Cutoff(const vectorfield& spins, scalarfield & Energy); void E_DDI_FFT(const vectorfield& spins, scalarfield & Energy); + // Triplet + void E_Triplet(const vectorfield & spins, scalarfield & Energy); + // Quadruplet void E_Quadruplet(const vectorfield & spins, scalarfield & Energy); @@ -189,4 +200,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 092077ed5..8c404ec14 100644 --- a/core/include/engine/Vectormath_Defines.hpp +++ b/core/include/engine/Vectormath_Defines.hpp @@ -46,6 +46,7 @@ using Matrix3c = Eigen::Matrix, 3, 3>; { int i, j, k; int d_j[3], d_k[3]; + scalar n[3]; }; struct Quadruplet { @@ -73,6 +74,7 @@ using Matrix3c = Eigen::Matrix, 3, 3>; { int i, j, k; std::array d_j, d_k; + std::array n; }; struct Quadruplet { @@ -100,4 +102,4 @@ using vectorfield = field; using pairfield = field; using tripletfield = field; using quadrupletfield = field; -using neighbourfield = field; \ No newline at end of file +using neighbourfield = field; diff --git a/core/include/io/Dataparser.hpp b/core/include/io/Dataparser.hpp index c067f4c20..2c0820836 100644 --- a/core/include/io/Dataparser.hpp +++ b/core/include/io/Dataparser.hpp @@ -29,6 +29,9 @@ namespace IO pairfield& exchange_pairs, scalarfield& exchange_magnitudes, pairfield& dmi_pairs, scalarfield& dmi_magnitudes, vectorfield& dmi_normals ); + 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 ); diff --git a/core/src/engine/Hamiltonian_Heisenberg.cpp b/core/src/engine/Hamiltonian_Heisenberg.cpp index f51289d4d..af5595418 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, 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_cutoff_radius(ddi_radius) { @@ -52,6 +54,7 @@ namespace Engine scalarfield exchange_shell_magnitudes, scalarfield dmi_shell_magnitudes, int dm_chirality, DDI_Method ddi_method, intfield ddi_n_periodic_images, scalar ddi_radius, + tripletfield triplets, scalarfield triplet_magnitudes1, scalarfield triplet_magnitudes2, quadrupletfield quadruplets, scalarfield quadruplet_magnitudes, std::shared_ptr geometry, intfield boundary_conditions @@ -62,6 +65,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_cutoff_radius(ddi_radius) { @@ -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); } @@ -422,6 +435,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) @@ -559,6 +606,41 @@ 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 + } + } + // Quadruplets if (this->idx_quadruplet >= 0) { @@ -615,6 +697,9 @@ namespace Engine // DD this->Gradient_DDI(spins, gradient); + // Triplets + this->Gradient_Triplet(spins, gradient); + // Quadruplets this->Gradient_Quadruplet(spins, gradient); } @@ -873,6 +958,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) { @@ -1314,4 +1444,4 @@ namespace Engine const std::string& Hamiltonian_Heisenberg::Name() { return name; } } -#endif \ No newline at end of file +#endif diff --git a/core/src/io/Configparser.cpp b/core/src/io/Configparser.cpp index ffbaaf6ec..de9a395c3 100644 --- a/core/src/io/Configparser.cpp +++ b/core/src/io/Configparser.cpp @@ -1247,6 +1247,13 @@ namespace IO intfield ddi_n_periodic_images = { 4, 4, 4 }; scalar ddi_radius = 0.0; + // ------------ 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 = ""; @@ -1468,6 +1475,27 @@ namespace IO 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); @@ -1526,6 +1554,7 @@ namespace IO exchange_magnitudes, dmi_magnitudes, dm_chirality, ddi_method, ddi_n_periodic_images, ddi_radius, + triplets, triplet_magnitudes1, triplet_magnitudes2, quadruplets, quadruplet_magnitudes, geometry, boundary_conditions @@ -1539,6 +1568,7 @@ namespace IO exchange_pairs, exchange_magnitudes, dmi_pairs, dmi_magnitudes, dmi_normals, ddi_method, ddi_n_periodic_images, ddi_radius, + triplets, triplet_magnitudes1, triplet_magnitudes2, quadruplets, quadruplet_magnitudes, geometry, boundary_conditions diff --git a/core/src/io/Configwriter.cpp b/core/src/io/Configwriter.cpp index f428c37f1..d2554bc89 100644 --- a/core/src/io/Configwriter.cpp +++ b/core/src/io/Configwriter.cpp @@ -281,6 +281,24 @@ namespace IO 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; itriplets.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()); @@ -317,4 +335,4 @@ namespace IO } Append_String_to_File(config, configFile); }// end Hamiltonian_Gaussian_to_Config -}// end namespace IO \ No newline at end of file +}// end namespace IO diff --git a/core/src/io/Dataparser.cpp b/core/src/io/Dataparser.cpp index 5761447d0..c4d99f9b6 100644 --- a/core/src/io/Dataparser.cpp +++ b/core/src/io/Dataparser.cpp @@ -492,6 +492,138 @@ namespace IO } } + /* + 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 /* From 263a7ec3ac6f790a4277c38253c8882a70926efb Mon Sep 17 00:00:00 2001 From: "G. P. Mueller" Date: Sun, 30 May 2021 20:16:41 +0200 Subject: [PATCH 2/2] Triplets: fixes from merge and added cuda backend. --- core/src/engine/Hamiltonian_Heisenberg.cpp | 6 +- core/src/engine/Hamiltonian_Heisenberg.cu | 131 ++++++++++++++++++++- core/src/io/Configparser.cpp | 8 +- 3 files changed, 138 insertions(+), 7 deletions(-) diff --git a/core/src/engine/Hamiltonian_Heisenberg.cpp b/core/src/engine/Hamiltonian_Heisenberg.cpp index 650ca4d24..148ac1b6a 100644 --- a/core/src/engine/Hamiltonian_Heisenberg.cpp +++ b/core/src/engine/Hamiltonian_Heisenberg.cpp @@ -210,7 +210,7 @@ namespace Engine this->idx_ddi = this->energy_contributions_per_spin.size()-1; } else this->idx_ddi = -1; - // Triplets + // Triplets if (this->triplets.size() > 0) { this->energy_contributions_per_spin.push_back({"triplets", scalarfield(0) }); @@ -433,7 +433,7 @@ namespace Engine } } -void Hamiltonian_Heisenberg::E_Triplet(const vectorfield & spins, scalarfield & Energy) + void Hamiltonian_Heisenberg::E_Triplet(const vectorfield & spins, scalarfield & Energy) { for (unsigned int itrip = 0; itrip < triplets.size(); ++itrip) { @@ -968,7 +968,7 @@ void Hamiltonian_Heisenberg::E_Triplet(const vectorfield & spins, scalarfield & } } -void Hamiltonian_Heisenberg::Gradient_Triplet(const vectorfield & spins, vectorfield & gradient) + void Hamiltonian_Heisenberg::Gradient_Triplet(const vectorfield & spins, vectorfield & gradient) { for (unsigned int itrip = 0; itrip < triplets.size(); ++itrip) { 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 e22cb7e40..2f0f5aad8 100644 --- a/core/src/io/Configparser.cpp +++ b/core/src/io/Configparser.cpp @@ -1552,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;