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.

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
Show file tree
Hide file tree
Changes from all commits
Commits
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
24 changes: 18 additions & 6 deletions core/docs/Input.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
```
Expand Down Expand Up @@ -557,4 +569,4 @@ inside the file.

---

[Home](Readme.md)
[Home](Readme.md)
17 changes: 14 additions & 3 deletions core/include/engine/Hamiltonian_Heisenberg.hpp
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, bool ddi_pb_zero_padding, 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, bool ddi_pb_zero_padding, 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 @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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);

Expand All @@ -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
Expand All @@ -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);
Expand Down Expand Up @@ -203,4 +214,4 @@ namespace Engine


}
#endif
#endif
2 changes: 2 additions & 0 deletions core/include/engine/Vectormath_Defines.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,7 @@ using Vector2 = Eigen::Matrix<scalar, 2, 1>;
{
int i, j, k;
int d_j[3], d_k[3];
scalar n[3];
};
struct Quadruplet
{
Expand Down Expand Up @@ -77,6 +78,7 @@ using Vector2 = Eigen::Matrix<scalar, 2, 1>;
{
int i, j, k;
std::array<int,3> d_j, d_k;
std::array<scalar, 3> n;
};
struct Quadruplet
{
Expand Down
5 changes: 4 additions & 1 deletion core/include/io/Dataparser.hpp
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@

#pragma once
#ifndef SPIRIT_IO_DATAPARSER_HPP
#define SPIRIT_IO_DATAPARSER_HPP
Expand Down Expand Up @@ -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<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 ) noexcept;
Expand Down
133 changes: 131 additions & 2 deletions core/src/engine/Hamiltonian_Heisenberg.cpp
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, bool ddi_pb_zero_padding, 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_pb_zero_padding(ddi_pb_zero_padding) ,ddi_cutoff_radius(ddi_radius),
fft_plan_reverse(FFT::FFT_Plan()), fft_plan_spins(FFT::FFT_Plan())
Expand All @@ -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<Data::Geometry> geometry,
intfield boundary_conditions
Expand All @@ -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())
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 @@ -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<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 @@ -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)
{
}
}
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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<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 @@ -1392,4 +1521,4 @@ namespace Engine
const std::string& Hamiltonian_Heisenberg::Name() { return name; }
}

#endif
#endif
Loading