message: Core: attempt to fix DDI in both Heisenberg Hamiltonians.
  Added mu_s factors and 1/2 to energies and gradients, which should hopefully fix the bad DDI calculation.
diff: @@ -269,13 +269,15 @@ namespace Engine
  for (int dc = 0; dc < geometry->n_cells[2]; ++dc)
  {
  std::array<int, 3 > translations = { da, db, dc };
  - int ispin = ddi_pairs[i_pair].i + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations);
  + int i = ddi_pairs[i_pair].i;
  + int j = ddi_pairs[i_pair].j;
  + int ispin = i + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations);
  int jspin = idx_from_pair(ispin, boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, ddi_pairs[i_pair]);
  if (jspin >= 0)
  {
  - Energy[ispin] -= mult / std::pow(ddi_magnitudes[i_pair], 3.0) *
  + Energy[ispin] -= 0.5 * this->mu_s[i] * this->mu_s[j] * mult / std::pow(ddi_magnitudes[i_pair], 3.0) *
  (3 * spins[ispin].dot(ddi_normals[i_pair]) * spins[ispin].dot(ddi_normals[i_pair]) - spins[ispin].dot(spins[ispin]));
  - Energy[ispin] -= mult / std::pow(ddi_magnitudes[i_pair], 3.0) *
  + Energy[jspin] -= 0.5 * this->mu_s[i] * this->mu_s[j] * mult / std::pow(ddi_magnitudes[i_pair], 3.0) *
  (3 * spins[ispin].dot(ddi_normals[i_pair]) * spins[ispin].dot(ddi_normals[i_pair]) - spins[ispin].dot(spins[ispin]));
  }
  }
  @@ -474,12 +476,14 @@ namespace Engine
  scalar skalar_contrib = mult / std::pow(ddi_magnitudes[i_pair], 3.0);
  std::array<int, 3 > translations = { da, db, dc };
  - int ispin = ddi_pairs[i_pair].i + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations);
  + int i = ddi_pairs[i_pair].i;
  + int j = ddi_pairs[i_pair].j;
  + int ispin = i + Vectormath::idx_from_translations(geometry->n_cells, geometry->n_cell_atoms, translations);
  int jspin = idx_from_pair(ispin, boundary_conditions, geometry->n_cells, geometry->n_cell_atoms, geometry->atom_types, ddi_pairs[i_pair]);
  if (jspin >= 0)
  {
  - gradient[ispin] -= skalar_contrib * (3 * ddi_normals[i_pair] * spins[jspin].dot(ddi_normals[i_pair]) - spins[jspin]);
  - gradient[jspin] -= skalar_contrib * (3 * ddi_normals[i_pair] * spins[ispin].dot(ddi_normals[i_pair]) - spins[ispin]);
  + gradient[ispin] -= this->mu_s[j] * skalar_contrib * (3 * ddi_normals[i_pair] * spins[jspin].dot(ddi_normals[i_pair]) - spins[jspin]);
  + gradient[jspin] -= this->mu_s[i] * skalar_contrib * (3 * ddi_normals[i_pair] * spins[ispin].dot(ddi_normals[i_pair]) - spins[ispin]);
  }
  }
  }
  