Hamiltonian_Heisenberg: DDI FFT optional skipping of zero-padding for PB

- added input keyword to trigger skipping of zero-padding for periodical directions
- also added it to the API and GUI
- fixed get_ddi in python API

core/docs/Input.md

@@ -185,6 +185,8 @@ ddi_n_periodic_images    4 4 4
 
 ### DDI cutoff radius (if cutoff is used)
 ddi_radius               0.0
+
+ddi_pb_zero_padding      1.0
 ```
 
 *Anisotropy:*
@@ -197,7 +199,7 @@ a magnitude `K`.
 Via the keyword `ddi_method` the method employed to calculate the dipole-dipole interactions is specified.
 
       `none`   -  Dipole-Dipole interactions are neglected
-      `fft`    -  Uses a fast convolution method to accelerate the calculation
+      `fft`    -  Uses a fast convolution method to accelerate the calculation (RECOMMENDED)
       `cutoff` -  Lets only spins within a maximal distance of 'ddi_radius' interact
       `fmm`    -  Uses the Fast-Multipole-Method (NOT YET IMPLEMENTED!)
 
@@ -208,6 +210,10 @@ If the boundary conditions are periodic `ddi_n_periodic_images` specifies how ma
 *Note:* The images are appended on both sides (the edges get filled too)
 i.e. 1 0 0 -> one image in +a direction and one image in -a direction
 
+If the boundary conditions are open in a lattice direction and sufficiently many periodic images are chosen, zero-padding in that direction can be skipped. 
+This improves the speed and memory footprint of the calculation, but comes at the cost of a very slight asymmetry in the interactions (decreasing with increasing periodic images). 
+If `ddi_pb_zero_padding` is set to 1, zero-padding is performed - even if the boundary condition is periodic in a direction. If it is set to 0, zero-padding is skipped.
+
 **Neighbour shells**:
 
 Using `hamiltonian heisenberg_neighbours`, pair-wise interactions are handled in terms of

core/include/Spirit/Hamiltonian.h

@@ -85,8 +85,9 @@ Configure the dipole-dipole interaction
   boundary conditions are used
 - `cutoff_radius`: the distance at which to stop the direct summation,
   if used
+- `pb_zero_padding`: if `True` zero padding is used even for periodical directions
 */
-PREFIX void Hamiltonian_Set_DDI(State *state, int ddi_method, int n_periodic_images[3], float cutoff_radius=0, int idx_image=-1, int idx_chain=-1) SUFFIX;
+PREFIX void Hamiltonian_Set_DDI(State *state, int ddi_method, int n_periodic_images[3], float cutoff_radius=0, bool pb_zero_padding=true, int idx_image=-1, int idx_chain=-1) SUFFIX;
 
 /*
 Getters
@@ -133,13 +134,12 @@ PREFIX int  Hamiltonian_Get_DMI_N_Pairs(State *state, int idx_image=-1, int idx_
 Retrieves the dipole-dipole interaction configuration.
 
 - `ddi_method`: see integers defined above
-- `n_periodic_images`: how many repetition of the spin configuration to
-  append along the translation directions [a, b, c], if periodical
-  boundary conditions are used
-- `cutoff_radius`: the distance at which to stop the direct summation,
-  if used
+- `n_periodic_images`: how many repetitions of the spin configuration to
+  append along the translation directions [a, b, c], if periodical boundary conditions are used
+- `cutoff_radius`: the distance at which to stop the direct summation, if method_cutoff is used
+- `pb_zero_padding`: if `True` zero padding is used even for periodical directions
 */
-PREFIX void Hamiltonian_Get_DDI(State *state, int * ddi_method, int n_periodic_images[3], float * cutoff_radius, int idx_image=-1, int idx_chain=-1) SUFFIX;
+PREFIX void Hamiltonian_Get_DDI(State *state, int * ddi_method, int n_periodic_images[3], float * cutoff_radius,  bool * pb_zero_padding, int idx_image=-1, int idx_chain=-1) SUFFIX;
 
 #include "DLL_Undefine_Export.h"
 #endif

core/include/engine/Hamiltonian_Heisenberg.hpp

@@ -35,7 +35,7 @@ namespace Engine
             intfield anisotropy_indices, scalarfield anisotropy_magnitudes, vectorfield anisotropy_normals,
             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,
+            DDI_Method ddi_method, intfield ddi_n_periodic_images, bool ddi_pb_zero_padding, scalar ddi_radius,
             quadrupletfield quadruplets, scalarfield quadruplet_magnitudes,
             std::shared_ptr<Data::Geometry> geometry,
             intfield boundary_conditions
@@ -46,7 +46,7 @@ namespace Engine
             intfield anisotropy_indices, scalarfield anisotropy_magnitudes, vectorfield anisotropy_normals,
             scalarfield exchange_shell_magnitudes,
             scalarfield dmi_shell_magnitudes, int dm_chirality,
-            DDI_Method ddi_method, intfield ddi_n_periodic_images, scalar ddi_radius,
+            DDI_Method ddi_method, intfield ddi_n_periodic_images, bool ddi_pb_zero_padding, scalar ddi_radius,
             quadrupletfield quadruplets, scalarfield quadruplet_magnitudes,
             std::shared_ptr<Data::Geometry> geometry,
             intfield boundary_conditions
@@ -100,6 +100,7 @@ namespace Engine
         // Dipole Dipole interaction
         DDI_Method  ddi_method;
         intfield    ddi_n_periodic_images;
+        bool        ddi_pb_zero_padding;
         //      ddi cutoff variables
         scalar      ddi_cutoff_radius;
         pairfield   ddi_pairs;
@@ -177,12 +178,12 @@ namespace Engine
         FFT::StrideContainer spin_stride;
         FFT::StrideContainer dipole_stride;
 
-        //Calculate the FT of the padded D matriess
+        // Calculate the FT of the padded D-matrics
         void FFT_Dipole_Matrices(FFT::FFT_Plan & fft_plan_dipole, int img_a, int img_b, int img_c);
-        //Calculate the FT of the padded spins
+        // Calculate the FT of the padded spins
         void FFT_Spins(const vectorfield & spins);
 
-        //Bounds for nested for loops. Only important for the CUDA version
+        // Bounds for nested for loops. Only important for the CUDA version
         field<int> it_bounds_pointwise_mult;
         field<int> it_bounds_write_gradients;
         field<int> it_bounds_write_spins;

core/python/spirit/hamiltonian.py

@@ -93,19 +93,20 @@ def set_dmi(p_state, n_shells, D_ij, chirality=CHIRALITY_BLOCH, idx_image=-1, id
              ctypes.c_int(chirality), ctypes.c_int(idx_image), ctypes.c_int(idx_chain))
 
 _Set_DDI             = _spirit.Hamiltonian_Set_DDI
-_Set_DDI.argtypes    = [ctypes.c_void_p, ctypes.c_int, ctypes.POINTER(ctypes.c_int), ctypes.c_float,
+_Set_DDI.argtypes    = [ctypes.c_void_p, ctypes.c_int, ctypes.POINTER(ctypes.c_int), ctypes.c_float, ctypes.c_bool,
                         ctypes.c_int, ctypes.c_int]
 _Set_DDI.restype     = None
-def set_ddi(p_state, ddi_method, n_periodic_images=[4,4,4], radius=0.0, idx_image=-1, idx_chain=-1):
+def set_ddi(p_state, ddi_method, n_periodic_images=[4,4,4], radius=0.0, pb_zero_padding=True, idx_image=-1, idx_chain=-1):
     """Set the dipolar interaction calculation method.
 
     - `ddi_method`: one of the integers defined above
     - `n_periodic_images`: the number of periodical images in the three translation directions,
       taken into account when boundaries in the corresponding direction are periodical
     - `radius`: the cutoff radius for the direct summation method
+    - `pb_zero_padding`: if `True` zero padding is used for periodical directions
     """
     vec3 = ctypes.c_int * 3
-    _Set_DDI(ctypes.c_void_p(p_state), ctypes.c_int(ddi_method) , vec3(*n_periodic_images), ctypes.c_float(radius),
+    _Set_DDI(ctypes.c_void_p(p_state), ctypes.c_int(ddi_method) , vec3(*n_periodic_images), ctypes.c_float(radius), ctypes.c_bool(pb_zero_padding),
              ctypes.c_int(idx_image), ctypes.c_int(idx_chain))
 
 ### ---------------------------------- Get ----------------------------------
@@ -145,9 +146,16 @@ def get_field(p_state, idx_image=-1, idx_chain=-1):
     return float(magnitude), [n for n in normal]
 
 _Get_DDI          = _spirit.Hamiltonian_Get_DDI
-_Get_DDI.argtypes = [ctypes.c_void_p, ctypes.c_int, ctypes.c_int]
-_Get_DDI.restype  = ctypes.c_float
+_Get_DDI.argtypes = [ctypes.c_void_p, ctypes.POINTER(ctypes.c_int), ctypes.POINTER(ctypes.c_int), ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_bool), ctypes.c_int, ctypes.c_int]
+_Get_DDI.restype  = None
 def get_ddi(p_state, idx_image=-1, idx_chain=-1):
-    """Returns the cutoff radius of the DDI."""
-    return float(_Get_DDI(ctypes.c_void_p(p_state), ctypes.c_int(idx_image),
-                          ctypes.c_int(idx_chain)))
+    """Returns a dictionary, containing information about the used ddi settings"""
+    ddi_method = ctypes.c_int()
+    n_periodic_images = (3*ctypes.c_int)()
+    cutoff_radius = ctypes.c_float()
+    pb_zero_padding = ctypes.c_bool()
+    _Get_DDI(ctypes.c_void_p(p_state), ctypes.byref(ddi_method), n_periodic_images, ctypes.byref(cutoff_radius), ctypes.byref(pb_zero_padding), ctypes.c_int(idx_image), ctypes.c_int(idx_chain))
+    return { "ddi_method" : ddi_method.value,
+            "n_periodic_images" : [ i for i in n_periodic_images ],
+            "cutoff_radius" : cutoff_radius.value,
+            "pb_zero_padding" : pb_zero_padding.value }

core/src/Spirit/Hamiltonian.cpp

@@ -267,7 +267,7 @@ void Hamiltonian_Set_DMI(State *state, int n_shells, const float * dij, int chir
     }
 }
 
-void Hamiltonian_Set_DDI(State *state, int ddi_method, int n_periodic_images[3], float cutoff_radius, int idx_image, int idx_chain) noexcept
+void Hamiltonian_Set_DDI(State *state, int ddi_method, int n_periodic_images[3], float cutoff_radius, bool pb_zero_padding, int idx_image, int idx_chain) noexcept
 {
     try
     {
@@ -289,11 +289,12 @@ void Hamiltonian_Set_DDI(State *state, int ddi_method, int n_periodic_images[3],
                 ham->ddi_n_periodic_images[1] = n_periodic_images[1];
                 ham->ddi_n_periodic_images[2] = n_periodic_images[2];
                 ham->ddi_cutoff_radius = cutoff_radius;
+                ham->ddi_pb_zero_padding = pb_zero_padding;
                 ham->Update_Interactions();
 
                 Log( Utility::Log_Level::Info, Utility::Log_Sender::API, fmt::format(
-                    "Set ddi to method {}, periodic images {} {} {} and cutoff radius {}",
-                    ddi_method, n_periodic_images[0], n_periodic_images[1], n_periodic_images[2], cutoff_radius), idx_image, idx_chain );
+                    "Set ddi to method {}, periodic images {} {} {}, cutoff radius {} and pb_zero_padding {}",
+                    ddi_method, n_periodic_images[0], n_periodic_images[1], n_periodic_images[2], cutoff_radius, pb_zero_padding), idx_image, idx_chain );
             }
             else
                 Log( Utility::Log_Level::Warning, Utility::Log_Sender::API, "DDI cannot be set on " + 
@@ -552,7 +553,7 @@ int  Hamiltonian_Get_DMI_N_Pairs(State *state, int idx_image, int idx_chain) noe
     }
 }
 
-void Hamiltonian_Get_DDI(State *state, int * ddi_method, int n_periodic_images[3], float * cutoff_radius, int idx_image, int idx_chain) noexcept
+void Hamiltonian_Get_DDI(State *state, int * ddi_method, int n_periodic_images[3], float * cutoff_radius, bool * pb_zero_padding, int idx_image, int idx_chain) noexcept
 {
     try
     {
@@ -571,6 +572,7 @@ void Hamiltonian_Get_DDI(State *state, int * ddi_method, int n_periodic_images[3
             n_periodic_images[1] = (int)ham->ddi_n_periodic_images[1];
             n_periodic_images[2] = (int)ham->ddi_n_periodic_images[2];
             *cutoff_radius       = (float)ham->ddi_cutoff_radius;
+            *pb_zero_padding     = ham->ddi_pb_zero_padding;
         }
     }
     catch( ... )

core/src/engine/Hamiltonian_Heisenberg.cpp

@@ -27,7 +27,7 @@ namespace Engine
         intfield anisotropy_indices, scalarfield anisotropy_magnitudes, vectorfield anisotropy_normals,
         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,
+        DDI_Method ddi_method, intfield ddi_n_periodic_images, bool ddi_pb_zero_padding, scalar ddi_radius,
         quadrupletfield quadruplets, scalarfield quadruplet_magnitudes,
         std::shared_ptr<Data::Geometry> geometry,
         intfield boundary_conditions
@@ -39,7 +39,7 @@ namespace Engine
         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),
         quadruplets(quadruplets), quadruplet_magnitudes(quadruplet_magnitudes),
-        ddi_method(ddi_method), ddi_n_periodic_images(ddi_n_periodic_images), ddi_cutoff_radius(ddi_radius),
+        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())
     {
         // Generate interaction pairs, constants etc.
@@ -52,7 +52,7 @@ namespace Engine
         intfield anisotropy_indices, scalarfield anisotropy_magnitudes, vectorfield anisotropy_normals,
         scalarfield exchange_shell_magnitudes,
         scalarfield dmi_shell_magnitudes, int dm_chirality,
-        DDI_Method ddi_method, intfield ddi_n_periodic_images, scalar ddi_radius,
+        DDI_Method ddi_method, intfield ddi_n_periodic_images, bool ddi_pb_zero_padding, scalar ddi_radius,
         quadrupletfield quadruplets, scalarfield quadruplet_magnitudes,
         std::shared_ptr<Data::Geometry> geometry,
         intfield boundary_conditions
@@ -64,7 +64,7 @@ namespace Engine
         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),
         quadruplets(quadruplets), quadruplet_magnitudes(quadruplet_magnitudes),
-        ddi_method(ddi_method), ddi_n_periodic_images(ddi_n_periodic_images), ddi_cutoff_radius(ddi_radius),
+        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())
 
     {
@@ -1165,15 +1165,23 @@ namespace Engine
         if( ddi_method != DDI_Method::FFT )
             return;
 
+        // We perform zero-padding in a lattice direction if the dimension of the system is greater than 1 *and*
+        //  - the boundary conditions are open, or
+        //  - the boundary conditions are periodic and zero-padding is explicitly requested
         n_cells_padded.resize(3);
-        n_cells_padded[0] = (geometry->n_cells[0] > 1) ? 2 * geometry->n_cells[0] : 1;
-        n_cells_padded[1] = (geometry->n_cells[1] > 1) ? 2 * geometry->n_cells[1] : 1;
-        n_cells_padded[2] = (geometry->n_cells[2] > 1) ? 2 * geometry->n_cells[2] : 1;
+        for(int i=0; i<3; i++)
+        {
+            n_cells_padded[i] = geometry->n_cells[i];
+            bool perform_zero_padding = geometry->n_cells[i] > 1 && (boundary_conditions[i] == 0 || ddi_pb_zero_padding);
+            if(perform_zero_padding)
+                n_cells_padded[i] *= 2;
+        }
+        sublattice_size = n_cells_padded[0] * n_cells_padded[1] * n_cells_padded[2];
 
         FFT::FFT_Init();
 
-        //workaround for bug in kissfft
-        //kissfft_ndr does not perform one-dimensional ffts properly
+        // Workaround for bug in kissfft
+        // kissfft_ndr does not perform one-dimensional FFTs properly
         #ifndef SPIRIT_USE_FFTW
         int number_of_one_dims = 0;
         for( int i=0; i<3; i++ )
@@ -1185,15 +1193,15 @@ namespace Engine
 
         inter_sublattice_lookup.resize(geometry->n_cell_atoms * geometry->n_cell_atoms);
 
-        //we dont need to transform over length 1 dims
+        // We dont need to transform over length 1 dims
         std::vector<int> fft_dims;
         for( int i = 2; i >= 0; i-- ) //notice that reverse order is important!
         {
             if(n_cells_padded[i] > 1)
                 fft_dims.push_back(n_cells_padded[i]);
         }
 
-        //Count how many distinct inter-lattice contributions we need to store
+        // Count how many distinct inter-lattice contributions we need to store
         n_inter_sublattice = 0;
         for( int i = 0; i < geometry->n_cell_atoms; i++ )
         {
@@ -1204,7 +1212,7 @@ namespace Engine
             }
         }
 
-        //Create fft plans.
+        // Create FFT plans
         FFT::FFT_Plan fft_plan_dipole  = FFT::FFT_Plan(fft_dims, false, 6 * n_inter_sublattice, sublattice_size);
         fft_plan_spins   = FFT::FFT_Plan(fft_dims, false, 3 * geometry->n_cell_atoms, sublattice_size);
         fft_plan_reverse = FFT::FFT_Plan(fft_dims, true, 3 * geometry->n_cell_atoms, sublattice_size);
@@ -1230,7 +1238,7 @@ namespace Engine
             (it_bounds_pointwise_mult[fft_dims.size() - 1] /= 2 )++;
         #endif
 
-        //perform FFT of dipole matrices
+        // Perform FFT of dipole matrices
         int img_a = boundary_conditions[0] == 0 ? 0 : ddi_n_periodic_images[0];
         int img_b = boundary_conditions[1] == 0 ? 0 : ddi_n_periodic_images[1];
         int img_c = boundary_conditions[2] == 0 ? 0 : ddi_n_periodic_images[2];