Merge branch 'develop' into feature-imgui

core/include/Spirit/System.h

@@ -53,8 +53,13 @@ PREFIX float System_Get_Rx(State * state, int idx_image=-1, int idx_chain=-1) SU
 // Returns the energy of a spin system.
 PREFIX float System_Get_Energy(State * state, int idx_image=-1, int idx_chain=-1) SUFFIX;
 
+// Retrieves the names of the energy contributions, represented as a single string and separated by "|". E.g "Zeeman|Exchange|DMI"
+// If 'names' is a nullptr, the required length of the char array is returned.
+PREFIX int System_Get_Energy_Array_Names(State * state, char* names, int idx_image=-1, int idx_chain=-1) SUFFIX;
+
 // Retrieves the energy contributions of a spin system.
-PREFIX void System_Get_Energy_Array(State * state, float * energies, int idx_image=-1, int idx_chain=-1) SUFFIX;
+// If 'energies' is a nullptr, the required length of the energies array is returned.
+PREFIX int System_Get_Energy_Array(State * state, float * energies, bool divide_by_nspins=true, int idx_image=-1, int idx_chain=-1) SUFFIX;
 
 // Retrieves the eigenvalues of a spin system
 PREFIX void System_Get_Eigenvalues(State * state, float * eigenvalues, int idx_image=-1, int idx_chain=-1) SUFFIX;

core/include/engine/FFT.hpp

@@ -3,6 +3,7 @@
 #define FFT_H
 #include <Eigen/Core>
 #include <data/Geometry.hpp>
+#include <utility/Logging.hpp>
 #include <engine/Vectormath_Defines.hpp>
 #include <iostream>
 #include <complex>
@@ -25,7 +26,7 @@
 #include <cufft.h>
 #endif
 
-namespace Engine 
+namespace Engine
 {
     namespace FFT
     {
@@ -84,7 +85,7 @@ namespace Engine
             FFT_cpx_type res;
             res[0] = d1[0] * s1[0] + d2[0] * s2[0] + d3[0] * s3[0] - d1[1] * s1[1] - d2[1] * s2[1] - d3[1] * s3[1];
             res[1] = d1[0] * s1[1] + d2[0] * s2[1] + d3[0] * s3[1] + d1[1] * s1[0] + d2[1] * s2[0] + d3[1] * s3[0];
-            return res;     
+            return res;
         }
 
         inline void addTo(FFT_cpx_type & a, const FFT_cpx_type & b, bool overwrite)
@@ -112,7 +113,7 @@ namespace Engine
             FFT_cpx_type res;
             res.x = d1.x * s1.x + d2.x * s2.x + d3.x * s3.x - d1.y * s1.y - d2.y * s2.y - d3.y * s3.y;
             res.y = d1.x * s1.y + d2.x * s2.y + d3.x * s3.y + d1.y * s1.x + d2.y * s2.x + d3.y * s3.x;
-            return res;   
+            return res;
         }
 
         inline __device__ void addTo(FFT_cpx_type & a, const FFT_cpx_type & b, bool overwrite = false)
@@ -135,7 +136,7 @@ namespace Engine
                 fftw_plan_with_nthreads(omp_get_max_threads());
             #endif
         }
-        
+
         inline void get_strides(field<int*> & strides, const field<int> & maxVal)
         {
             strides.resize(maxVal.size());
@@ -173,7 +174,7 @@ namespace Engine
 
             //Constructor delegation
             FFT_Plan() : FFT_Plan({2,2,2}, true, 1, 8)
-            {} 
+            {}
 
             FFT_Plan(std::vector<int> dims, bool inverse, int n_transforms, int len) :
                 dims(dims),
@@ -187,7 +188,7 @@ namespace Engine
             }
 
             //copy constructor
-            FFT_Plan(FFT_Plan const & other) 
+            FFT_Plan(FFT_Plan const & other)
             {
                 this->dims     = other.dims;
                 this->inverse  = other.inverse;
@@ -219,7 +220,7 @@ namespace Engine
                 }
                 return *this;
             }
-            
+
             //move assignment operator
             FFT_Plan& operator=(FFT_Plan const && other)
             {

core/include/engine/Sparse_HTST.hpp

@@ -52,7 +52,7 @@ namespace Engine
             SpMatrixX & hessian_geodesic_3N, SpMatrixX & hessian_geodesic_2N, SpMatrixX & tangent_basis, VectorX & eigenvalues, MatrixX & eigenvectors);
 
         void sparse_hessian_bordered_3N(const vectorfield & image, const vectorfield & gradient, const SpMatrixX & hessian, SpMatrixX & hessian_out);
-        
+
     };
 }
 

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/python/spirit/system.py

@@ -98,16 +98,37 @@ def get_eigenvalues(p_state, idx_image=-1, idx_chain=-1):
     _Get_Eigenvalues(ctypes.c_void_p(p_state), eigenvalues, ctypes.c_int(idx_image), ctypes.c_int(idx_chain))
     return eigenvalues
 
-# NOTE: excluded since there is no clean way to get the C++ pairs
-### Get Energy array
-# _Get_Energy_Array          = _spirit.System_Get_Energy_Array
-# _Get_Energy_Array.argtypes = [ctypes.c_void_p, ctypes.POINTER(ctypes.c_float),
-#                               ctypes.c_int, ctypes.c_int]
-# _Get_Energy_Array.restype  = None
-# def Get_Energy_Array(p_state, idx_image=-1, idx_chain=-1):
-#     Energies
-#     _Get_Energy_Array(ctypes.c_void_p(p_state), energies,
-#                       ctypes.c_int(idx_image), ctypes.c_int(idx_chain))
+### Get Energy Contributions
+### The result is a dictionary with strings as keys and floats as values
+### The keys are the names of the energy contributions, the values the energy_contribution in meV
+_Get_Energy_Array       = _spirit.System_Get_Energy_Array
+_Get_Energy_Array.argtypes = [ctypes.c_void_p, ctypes.POINTER(ctypes.c_float), ctypes.c_bool,
+                              ctypes.c_int, ctypes.c_int]
+_Get_Energy_Array.restype  = None
+
+_Get_Energy_Array_Names = _spirit.System_Get_Energy_Array_Names
+_Get_Energy_Array_Names.argtypes = [ctypes.c_void_p, ctypes.POINTER(ctypes.c_char),
+                                    ctypes.c_int, ctypes.c_int]
+_Get_Energy_Array_Names.restype  = ctypes.c_int
+def get_energy_contributions(p_state, divide_by_nspins = True, idx_image=-1, idx_chain=-1):
+    NULL = ctypes.POINTER(ctypes.c_char)()
+
+    n_char_array = _Get_Energy_Array_Names(ctypes.c_void_p(p_state), NULL,
+                    ctypes.c_int(idx_image), ctypes.c_int(idx_chain))
+
+    energy_array_names = (n_char_array*ctypes.c_char)()
+
+    _Get_Energy_Array_Names(ctypes.c_void_p(p_state), energy_array_names,
+                    ctypes.c_int(idx_image), ctypes.c_int(idx_chain))
+
+    contrib_names = str(energy_array_names[:].decode("utf-8")).split("|")
+    n_contribs = len(contrib_names)
+    energies = (n_contribs*ctypes.c_float)()
+
+    _Get_Energy_Array(ctypes.c_void_p(p_state), energies, divide_by_nspins,
+                      ctypes.c_int(idx_image), ctypes.c_int(idx_chain))
+
+    return dict(zip(contrib_names, energies))
 
 ### Get Chain number of images
 _Update_Data            = _spirit.System_Update_Data

core/python/test/system.py

@@ -5,7 +5,7 @@
 spirit_py_dir = os.path.abspath(os.path.join(os.path.dirname( __file__ ), ".."))
 sys.path.insert(0, spirit_py_dir)
 
-from spirit import state, system, configuration
+from spirit import state, system, configuration, hamiltonian
 
 import unittest
 
@@ -42,8 +42,16 @@ def test_get_spin_directions(self):
     def test_get_energy(self):
         # NOTE: that test is trivial
         E = system.get_energy(self.p_state)
-
-
+
+    def test_get_energy_contributions(self):
+        configuration.plus_z(self.p_state)
+        configuration.domain(self.p_state, [0,0,-1], border_cylindrical=2)
+        system.update_data(self.p_state)
+        E_contribs = system.get_energy_contributions(self.p_state, divide_by_nspins=False)
+        E = system.get_energy(self.p_state)
+        system.print_energy_array(p_state)
+        self.assertEqual( len(E_contribs.values()), 3 ) # There should be 3 contributions
+        self.assertAlmostEqual( sum(E_contribs.values()), E, places=5 ) #TODO: Apparently we can not go higher with the number of decimal places, because the order of summation differs. This Should be invesitgated.
     # NOTE: there is no way to test the system.Update_Data() and system.Print_Energy_Array()
 
 #########

core/src/Spirit/System.cpp

@@ -133,7 +133,7 @@ catch( ... )
     return 0;
 }
 
-void System_Get_Energy_Array(State * state, float * energies, int idx_image, int idx_chain) noexcept
+int System_Get_Energy_Array_Names(State * state, char* names, int idx_image, int idx_chain) noexcept
 try
 {
     std::shared_ptr<Data::Spin_System> image;
@@ -142,14 +142,63 @@ try
     // Fetch correct indices and pointers
     from_indices( state, idx_image, idx_chain, image, chain );
 
+    int n_char_array = -1; // Start of with offset -1, because the last contributions gets no "|" delimiter
     for (unsigned int i=0; i<image->E_array.size(); ++i)
     {
-        energies[i] = (float)image->E_array[i].second;
+        n_char_array += image->E_array[i].first.size() + 1; // Add +1 because we separate the contribution names with the character "|"
+    }
+
+    // If 'names' is a nullptr, we return the required length of the names array
+    if(names==nullptr)
+    {
+        return n_char_array;
+    } else { // Else we try to fill the provided char array
+        int idx=0;
+        for (unsigned int i=0; i<image->E_array.size(); ++i)
+        {
+            for(const char & cur_char : (image->E_array[i]).first)
+            {
+                names[idx++] = cur_char;
+            }
+            if(i != image->E_array.size()-1)
+                names[idx++] = '|';
+        }
+        return -1;
     }
 }
 catch( ... )
 {
     spirit_handle_exception_api(idx_image, idx_chain);
+    return -1;
+}
+
+
+int System_Get_Energy_Array(State * state, float * energies, bool divide_by_nspins, int idx_image, int idx_chain) noexcept
+try
+{
+    std::shared_ptr<Data::Spin_System> image;
+    std::shared_ptr<Data::Spin_System_Chain> chain;
+
+    // Fetch correct indices and pointers
+    from_indices( state, idx_image, idx_chain, image, chain );
+
+    scalar nd = divide_by_nspins ? 1/(scalar)image->nos : 1;
+
+    if(energies == nullptr)
+    {
+        return image->E_array.size();
+    } else {
+        for (unsigned int i=0; i<image->E_array.size(); ++i)
+        {
+            energies[i] = nd * (float)image->E_array[i].second;
+        }
+        return -1;
+    }
+}
+catch( ... )
+{
+    spirit_handle_exception_api(idx_image, idx_chain);
+    return -1;
 }
 
 void System_Get_Eigenvalues(State * state, float * eigenvalues, int idx_image, int idx_chain) noexcept

core/src/engine/FFT.cu

@@ -13,20 +13,29 @@ namespace Engine
         {
             std::cerr << "NOT IMPLEMENTED FOR cuFFT" << std::endl;
         }
+
         void iFour_3D(FFT_cfg cfg, FFT_cpx_type * in, FFT_real_type * out)
         {
             std::cerr << "NOT IMPLEMENTED FOR cuFFT" << std::endl;
         }
 
         void batch_Four_3D(FFT_Plan & plan)
         {
-            cufftExecR2C(plan.cfg, plan.real_ptr.data(), plan.cpx_ptr.data());
+            auto res = cufftExecR2C(plan.cfg, plan.real_ptr.data(), plan.cpx_ptr.data());
+            if(res != CUFFT_SUCCESS)
+            {
+                Log(Utility::Log_Level::Error, Utility::Log_Sender::All, fmt::format("cufftExecR2C failed with error: {}", res));
+            }
             cudaDeviceSynchronize();
         }
 
         void batch_iFour_3D(FFT_Plan & plan)
         {
-            cufftExecC2R(plan.cfg, plan.cpx_ptr.data(), plan.real_ptr.data());
+            auto res = cufftExecC2R(plan.cfg, plan.cpx_ptr.data(), plan.real_ptr.data());
+            if(res != CUFFT_SUCCESS)
+            {
+                Log(Utility::Log_Level::Error, Utility::Log_Sender::All, fmt::format("cufftExecC2R failed with error: {}", res));
+            }
             cudaDeviceSynchronize();
         }
 
@@ -47,16 +56,28 @@ namespace Engine
 
             if(this->inverse == false)
             {
-                cufftPlanMany(&this->cfg, rank, n, inembed, istride, idist, onembed, ostride, odist, CUFFT_R2C, n_transforms);
+                auto res = cufftPlanMany(&this->cfg, rank, n, inembed, istride, idist, onembed, ostride, odist, CUFFT_R2C, n_transforms);
+                if(res != CUFFT_SUCCESS)
+                {
+                    Log(Utility::Log_Level::Error, Utility::Log_Sender::All, fmt::format("cufftPlanMany failed with error: {}", res));
+                }
             } else
             {
-                cufftPlanMany(&this->cfg, rank, n, inembed, istride, idist, onembed, ostride, odist, CUFFT_C2R, n_transforms);
+                auto res = cufftPlanMany(&this->cfg, rank, n, inembed, istride, idist, onembed, ostride, odist, CUFFT_C2R, n_transforms);
+                if(res != CUFFT_SUCCESS)
+                {
+                    Log(Utility::Log_Level::Error, Utility::Log_Sender::All, fmt::format("cufftPlanMany failed with error: {}", res));
+                }
             }
         }
 
         void FFT_Plan::Free_Configuration()
         {
-            cufftDestroy(this->cfg);
+            auto res = cufftDestroy(this->cfg);
+            if(res != CUFFT_SUCCESS)
+            {
+                Log(Utility::Log_Level::Error, Utility::Log_Sender::All, fmt::format("cufftDestroy failed with error: {}", res));
+            }
         }
     }
 }