Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

added triplet interactions: ring-exchange and spin-chirality coupling #462

Open
wants to merge 3 commits into
base: develop
Choose a base branch
from
Open

added triplet interactions: ring-exchange and spin-chirality coupling #462

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
764d87c
added triplet interactions: ring-exchange and spin-chirality coupling
SpiritSuperUser Oct 26, 2018
da8afc9
Merge remote-tracking branch 'upstream/develop' into triplet
GPMueller May 30, 2021
263a7ec
Triplets: fixes from merge and added cuda backend.
GPMueller May 30, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 16 additions & 5 deletions core/docs/INPUT.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
```
Expand Down Expand Up @@ -506,4 +517,4 @@ inside the file.

---

[Home](Readme.md)
[Home](Readme.md)
15 changes: 13 additions & 2 deletions core/include/engine/Hamiltonian_Heisenberg.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<Data::Geometry> geometry,
intfield boundary_conditions
Expand All @@ -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<Data::Geometry> geometry,
intfield boundary_conditions
Expand Down Expand Up @@ -105,6 +107,10 @@ namespace Engine
scalarfield ddi_magnitudes;
vectorfield ddi_normals;

// ------------ Triplet Interactions ------------
tripletfield triplets;
GPMueller marked this conversation as resolved.
Show resolved Hide resolved
scalarfield triplet_magnitudes1, triplet_magnitudes2;

// ------------ Quadruplet Interactions ------------
quadrupletfield quadruplets;
scalarfield quadruplet_magnitudes;
Expand All @@ -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
Expand All @@ -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);

Expand Down Expand Up @@ -189,4 +200,4 @@ namespace Engine


}
#endif
#endif
4 changes: 3 additions & 1 deletion core/include/engine/Vectormath_Defines.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@ using Matrix3c = Eigen::Matrix<std::complex<scalar>, 3, 3>;
{
int i, j, k;
int d_j[3], d_k[3];
scalar n[3];
};
struct Quadruplet
{
Expand Down Expand Up @@ -73,6 +74,7 @@ using Matrix3c = Eigen::Matrix<std::complex<scalar>, 3, 3>;
{
int i, j, k;
std::array<int,3> d_j, d_k;
std::array<scalar, 3> n;
};
struct Quadruplet
{
Expand Down Expand Up @@ -100,4 +102,4 @@ using vectorfield = field<Vector3>;
using pairfield = field<Pair>;
using tripletfield = field<Triplet>;
using quadrupletfield = field<Quadruplet>;
using neighbourfield = field<Neighbour>;
using neighbourfield = field<Neighbour>;
3 changes: 3 additions & 0 deletions core/include/io/Dataparser.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<Data::Geometry> geometry, int& noq,
tripletfield& triplets, scalarfield& triplet_magnitudes, scalarfield& triplet_magnitudes2 );
void Quadruplets_from_File( const std::string quadrupletsFile,
const std::shared_ptr<Data::Geometry> geometry, int& noq,
quadrupletfield& quadruplets, scalarfield& quadruplet_magnitudes );
Expand Down
132 changes: 131 additions & 1 deletion core/src/engine/Hamiltonian_Heisenberg.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<Data::Geometry> geometry,
intfield boundary_conditions
Expand All @@ -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)
{
Expand All @@ -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<Data::Geometry> geometry,
intfield boundary_conditions
Expand All @@ -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)
{
Expand Down Expand Up @@ -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 )
{
Expand Down Expand Up @@ -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);
}
Expand Down Expand Up @@ -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<int, 3 > 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)
Expand Down Expand Up @@ -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)
{
Expand Down Expand Up @@ -615,6 +697,9 @@ namespace Engine
// DD
this->Gradient_DDI(spins, gradient);

// Triplets
this->Gradient_Triplet(spins, gradient);

// Quadruplets
this->Gradient_Quadruplet(spins, gradient);
}
Expand Down Expand Up @@ -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<int, 3 > 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)
{
Expand Down Expand Up @@ -1314,4 +1444,4 @@ namespace Engine
const std::string& Hamiltonian_Heisenberg::Name() { return name; }
}

#endif
#endif
Loading