Hamiltonian_Heisenberg:

- Exchange/DMI: cu_idx_from_pair was erroneously given the cell index as first argument (it should be the spin index)
- Zeeman: The index computation was missing a factor of n_cell_atoms

Vectormath:
- cu_idx_from_pair was different from idx_from_pair and missing a division by N (aka n_cell_atoms)

core/include/engine/Vectormath.hpp

@@ -272,7 +272,7 @@ namespace Engine
             // Translations (cell) of spin i
             int nic = ispin / (N*Na*Nb);
             int nib = (ispin - nic*N*Na*Nb) / (N*Na);
-            int nia = ispin - nic*N*Na*Nb - nib*N*Na;
+            int nia = (ispin - nic*N*Na*Nb - nib*N*Na) / N;
 
             // Translations (cell) of spin j (possibly outside of non-periodical domain)
             int pm = 1;
@@ -336,8 +336,7 @@ namespace Engine
             // Calculate the index of spin j according to it's translations
             int jspin = pair.j + (nja)*N + (njb)*N*Na + (njc)*N*Na*Nb;
 
-            // Invalid index if atom type of spin j is not correct
-            if ( pair.j != jspin%N || !cu_check_atom_type(atom_types[jspin]) )
+            if ( !cu_check_atom_type(atom_types[jspin]) )
                 return -1;
 
             // Return a valid index

core/src/engine/Hamiltonian_Heisenberg.cu

@@ -256,7 +256,7 @@ namespace Engine
         {
             for (int ibasis=0; ibasis<n_cell_atoms; ++ibasis)
             {
-                int ispin = icell + ibasis;
+                int ispin = n_cell_atoms * icell + ibasis;
                 if ( cu_check_atom_type(atom_types[ispin]) )
                     Energy[ispin] -= mu_s[ispin] * external_field_magnitude * external_field_normal.dot(spins[ispin]);
             }
@@ -305,7 +305,7 @@ namespace Engine
             for(auto ipair = 0; ipair < n_pairs; ++ipair)
             {
                 int ispin = pairs[ipair].i + icell*n_cell_atoms;
-                int jspin = cu_idx_from_pair(icell, bc, nc, n_cell_atoms, atom_types, pairs[ipair]);
+                int jspin = cu_idx_from_pair(ispin, bc, nc, n_cell_atoms, atom_types, pairs[ipair]);
                 if (jspin >= 0)
                 {
                     Energy[ispin] -= 0.5 * magnitudes[ipair] * spins[ispin].dot(spins[jspin]);
@@ -335,7 +335,7 @@ namespace Engine
             for(auto ipair = 0; ipair < n_pairs; ++ipair)
             {
                 int ispin = pairs[ipair].i + icell*n_cell_atoms;
-                int jspin = cu_idx_from_pair(icell, bc, nc, n_cell_atoms, atom_types, pairs[ipair]);
+                int jspin = cu_idx_from_pair(ispin, bc, nc, n_cell_atoms, atom_types, pairs[ipair]);
                 if (jspin >= 0)
                 {
                     Energy[ispin] -= 0.5 * magnitudes[ipair] * normals[ipair].dot(spins[ispin].cross(spins[jspin]));
@@ -708,7 +708,7 @@ namespace Engine
         {
             for (int ibasis=0; ibasis<n_cell_atoms; ++ibasis)
             {
-                int ispin = icell + ibasis;
+                int ispin = n_cell_atoms * icell + ibasis;
                 if ( cu_check_atom_type(atom_types[ispin]) )
                     gradient[ispin] -= mu_s[ispin] * external_field_magnitude*external_field_normal;
             }
@@ -760,7 +760,7 @@ namespace Engine
             for(auto ipair = 0; ipair < n_pairs; ++ipair)
             {
                 int ispin = pairs[ipair].i + icell*n_cell_atoms;
-                int jspin = cu_idx_from_pair(icell, bc, nc, n_cell_atoms, atom_types, pairs[ipair]);
+                int jspin = cu_idx_from_pair(ispin, bc, nc, n_cell_atoms, atom_types, pairs[ipair]);
                 if (jspin >= 0)
                 {
                     gradient[ispin] -= magnitudes[ipair]*spins[jspin];
@@ -790,7 +790,7 @@ namespace Engine
             for(auto ipair = 0; ipair < n_pairs; ++ipair)
             {
                 int ispin = pairs[ipair].i + icell*n_cell_atoms;
-                int jspin = cu_idx_from_pair(icell, bc, nc, n_cell_atoms, atom_types, pairs[ipair]);
+                int jspin = cu_idx_from_pair(ispin, bc, nc, n_cell_atoms, atom_types, pairs[ipair]);
                 if (jspin >= 0)
                 {
                     gradient[ispin] -= magnitudes[ipair]*spins[jspin].cross(normals[ipair]);