SeparatedConvolution can deal with mixed boundary conditions (periodi…

…c in some directions, nonperiodic in others) ... ?

periodicity flags use std::array, which already revealed some sketchy uses of is_periodic

src/madness/mra/bc.h

@@ -39,9 +39,9 @@
 
 #include <madness/world/madness_exception.h>
 
+#include <array>
 #include <cstddef>
 #include <iostream>
-#include <vector>
 
 namespace madness {
 
@@ -126,23 +126,40 @@ template <std::size_t NDIM> class BoundaryConditions {
 
   /// Convenience for application of integral operators
 
-  /// @return Returns a vector indicating if each dimension is periodic
-  std::vector<bool> is_periodic() const {
-    std::vector<bool> v(NDIM);
-    for (std::size_t d = 0; d < NDIM; ++d)
+  /// @return Returns a vector indicating if dimensions [0, ND) are periodic
+  template <std::size_t ND = NDIM>
+  std::enable_if_t<ND <= NDIM,std::array<bool, ND>> is_periodic() const {
+    std::array<bool, ND> v;
+    for (std::size_t d = 0; d < ND; ++d) {
+      MADNESS_ASSERT(bc[2 * d + 1] == bc[2 * d]);
       v[d] = (bc[2 * d] == BC_PERIODIC);
+    }
     return v;
   }
 
   /// Checks whether the boundary condition along any axis is periodic
 
   /// @return Returns true if any dimension is periodic
   bool is_periodic_any() const {
-    for (std::size_t d = 0; d < NDIM; ++d)
+    for (std::size_t d = 0; d < NDIM; ++d) {
+      MADNESS_ASSERT(bc[2 * d + 1] == bc[2 * d]);
       if (bc[2 * d] == BC_PERIODIC)
         return true;
+    }
     return false;
   }
+
+  /// Checks whether the boundary condition along all axes is periodic
+
+  /// @return Returns true if every dimension is periodic
+  bool is_periodic_all() const {
+    for (std::size_t d = 0; d < NDIM; ++d) {
+      MADNESS_ASSERT(bc[2 * d + 1] == bc[2 * d]);
+      if (bc[2 * d] != BC_PERIODIC)
+        return false;
+    }
+    return true;
+  }
 };
 
 template <std::size_t NDIM>
@@ -159,6 +176,13 @@ static inline std::ostream &operator<<(std::ostream &s,
   return s;
 }
 
+template <std::size_t NDIM>
+std::array<bool, NDIM> no_lattice_sum() {
+  std::array<bool, NDIM> result;
+  result.fill(false);
+  return result;
+}
+
 }  // namespace madness
 
 #endif // MADNESS_MRA_BC_H__INCLUDED

src/madness/mra/funcimpl.h

@@ -35,8 +35,6 @@
 /// \file funcimpl.h
 /// \brief Provides FunctionCommonData, FunctionImpl and FunctionFactory
 
-#include <iostream>
-#include <type_traits>
 #include <madness/world/MADworld.h>
 #include <madness/world/print.h>
 #include <madness/misc/misc.h>
@@ -50,7 +48,11 @@
 #include <madness/mra/function_factory.h>
 #include <madness/mra/displacements.h>
 
-#include "leafop.h"
+#include <madness/mra/leafop.h>
+
+#include <array>
+#include <iostream>
+#include <type_traits>
 
 namespace madness {
     template <typename T, std::size_t NDIM>
@@ -3820,7 +3822,7 @@ template<size_t NDIM>
         /// Out of volume keys are mapped to enforce the BC as follows.
         ///   * Periodic BC map back into the volume and return the correct key
         ///   * non-periodic BC - returns invalid() to indicate out of volume
-        keyT neighbor(const keyT& key, const keyT& disp, const std::vector<bool>& is_periodic) const;
+        keyT neighbor(const keyT& key, const keyT& disp, const std::array<bool, NDIM>& is_periodic) const;
 
         /// Returns key of general neighbor that resides in-volume
 
@@ -4524,7 +4526,7 @@ template<size_t NDIM>
         void zero_norm_tree();
 
         // Broaden tree
-        void broaden(std::vector<bool> is_periodic, bool fence);
+        void broaden(const std::array<bool, NDIM>& is_periodic, bool fence);
 
         /// sum all the contributions from all scales after applying an operator in mod-NS form
         void trickle_down(bool fence);
@@ -4810,10 +4812,25 @@ template<size_t NDIM>
 
             double cnorm = c.normf();
 
-            const std::vector<opkeyT>& disp = Displacements<opdim>().get_disp(key.level(), /* periodic = */ false); // list of displacements sorted in order of increasing distance; N.B. periodic displacements skipped since operator already includes lattice sum
+            // BC handling:
+            // - if operator is lattice-summed then treat this as nonperiodic (i.e. tell neighbor() to stay in simulation cell)
+            // - if operator is NOT lattice-summed then obey BC (i.e. tell neighbor() to go outside the simulation cell along periodic dimensions)
+            const auto bc_is_periodic = FunctionDefaults<NDIM>::get_bc().is_periodic();
+            const auto lattice_sum_in_op = op->includes_lattice_sum();
+            std::array<bool, NDIM> treat_this_as_periodic;
+            for(auto axis=0; axis!=NDIM; ++axis) {
+              treat_this_as_periodic[axis] =
+                  bc_is_periodic[axis] && !lattice_sum_in_op[axis];
+            }
+            // if this is treated as periodic in ANY direction need to use displacements that are periodic in EVERY direction
+            // because these are the only flavor of periodic displacements that we have
+            // don't worry, neighbor() will still treat nonperiodic axes as nonperiodic
+            const bool treat_this_as_periodic_in_any_direction = std::accumulate(treat_this_as_periodic.begin(), treat_this_as_periodic.end(), false, std::logical_or{});
+            const std::vector<opkeyT>& disp = Displacements<opdim>().get_disp(key.level(), /* periodic = */ treat_this_as_periodic_in_any_direction); // list of displacements sorted in order of increasing distance
+
             int nvalid=1;       // Counts #valid at each distance
             int nused=1;	// Counts #used at each distance
-	    uint64_t distsq = 99999999999999;
+            uint64_t distsq = 99999999999999;
             for (typename std::vector<opkeyT>::const_iterator it=disp.begin(); it != disp.end(); ++it) {
 	        keyT d;
                 Key<NDIM-opdim> nullkey(key.level());
@@ -4834,7 +4851,7 @@ template<size_t NDIM>
 		    distsq = dsq;
 		}
 
-                keyT dest = neighbor_in_volume(key, d);
+                keyT dest = neighbor(key, d, treat_this_as_periodic);
                 if (dest.is_valid()) {
                     nvalid++;
                     double opnorm = op->norm(key.level(), *it, source);

src/madness/mra/key.h

@@ -247,7 +247,7 @@ namespace madness {
 
         /// Assumes key and this are at the same level
         bool
-        is_neighbor_of(const Key& key, const std::vector<bool>& bperiodic) const {
+        is_neighbor_of(const Key& key, const std::array<bool, NDIM>& bperiodic) const {
           Translation dist = 0;
           Translation TWON1 = (Translation(1)<<n) - 1;
         	for (std::size_t i=0; i<NDIM; ++i)

src/madness/mra/leafop.h

@@ -84,7 +84,7 @@ struct Specialbox_op {
         if (key.level() > 1 and box_is_at_boundary(key)) return false;
         // check for all special points if they are neighbours of the current box
         BoundaryConditions<NDIM> bc = FunctionDefaults<NDIM>::get_bc();
-        std::vector<bool> bperiodic = bc.is_periodic();
+        const auto bperiodic = bc.is_periodic();
         for (size_t i = 0; i < special_points.size(); ++i) {
             Vector<double, NDIM> simpt;
             user_to_sim(special_points[i], simpt);
@@ -151,7 +151,7 @@ struct ElectronCuspyBox_op : public Specialbox_op<T, NDIM> {
     std::string name() const { return "Cuspybox_op"; }
 
     /// Operator which decides if the key belongs to a special box
-    /// The key is broken appart in two lower dimensional keys (for electron cusps this is 6D -> 2x3D)
+    /// The key is broken apart in two lower dimensional keys (for electron cusps this is 6D -> 2x3D)
     /// if the keys are neighbours then refinement up to the special level is enforced (means the 6D box is close to the cusp or contains it)
     /// if the refinement level is already beyond half of the special_level then refinement is only enforded if the broken keys are the same (6D box contains cusp)
     bool operator()(const Key<NDIM>& key, const FunctionImpl<T, NDIM> *const f) const {
@@ -160,21 +160,27 @@ struct ElectronCuspyBox_op : public Specialbox_op<T, NDIM> {
             if (this->box_is_at_boundary(key)) return false;
         }
 
-        if (NDIM % 2 != 0) MADNESS_EXCEPTION("Cuspybox_op only valid for even dimensions",
-                                             1);     // if uneven dims are needed just make a new class with NDIM+1/2 and NDIM-LDIM
-        BoundaryConditions<NDIM> bc = FunctionDefaults<NDIM>::get_bc();
-        std::vector<bool> bperiodic = bc.is_periodic();
-        Key<NDIM / 2> key1;
-        Key<NDIM / 2> key2;
-        key.break_apart(key1, key2);
-        int ll = this->get_half_of_special_level();
-        if (ll < f->get_initial_level()) ll = f->get_initial_level();
-        if (key.level() > ll) {
-            if (key1 == key2) return true;
-            else return false;
-        } else {
-            if (key1.is_neighbor_of(key2, bperiodic)) return true;
-            else return false;
+        // Cuspybox_op only valid for even dimensions, if uneven dims are needed just make a new class with NDIM+1/2 and NDIM-LDIM
+        if constexpr (NDIM % 2 == 0) {
+          BoundaryConditions<NDIM> bc = FunctionDefaults<NDIM>::get_bc();
+          const auto bperiodic = bc.template is_periodic<NDIM / 2>();
+          Key<NDIM / 2> key1;
+          Key<NDIM / 2> key2;
+          key.break_apart(key1, key2);
+          int ll = this->get_half_of_special_level();
+          if (ll < f->get_initial_level())
+            ll = f->get_initial_level();
+          if (key.level() > ll) {
+            if (key1 == key2)
+              return true;
+            else
+              return false;
+          } else {
+            if (key1.is_neighbor_of(key2, bperiodic))
+              return true;
+            else
+              return false;
+          }
         }
         MADNESS_EXCEPTION("We should not end up here (check further of cuspy box)", 1);
         return false;
@@ -239,7 +245,7 @@ struct NuclearCuspyBox_op : public Specialbox_op<T, NDIM> {
 
         // now break the key appart and check if one if the results is in the neighbourhood of a special point
         BoundaryConditions<NDIM / 2> bc = FunctionDefaults<NDIM / 2>::get_bc();
-        std::vector<bool> bperiodic = bc.is_periodic();
+        const auto bperiodic = bc.is_periodic();
         Key<NDIM / 2> key1;
         Key<NDIM / 2> key2;
 
@@ -346,7 +352,7 @@ class Leaf_op {
         sanity();
         const double cnorm = coeff.normf();
         BoundaryConditions<NDIM> bc = FunctionDefaults<NDIM>::get_bc();
-        std::vector<bool> bperiodic = bc.is_periodic();
+        const auto bperiodic = bc.is_periodic();
 
         typedef Key<opT::opdim> opkeyT;
         const opkeyT source = op->get_source_key(key);

src/madness/mra/mraimpl.h

@@ -1279,7 +1279,7 @@ namespace madness {
 
     // Broaden tree
     template <typename T, std::size_t NDIM>
-    void FunctionImpl<T,NDIM>::broaden(std::vector<bool> is_periodic, bool fence) {
+    void FunctionImpl<T,NDIM>::broaden(const std::array<bool, NDIM>& is_periodic, bool fence) {
         typename dcT::iterator end = coeffs.end();
         for (typename dcT::iterator it=coeffs.begin(); it!=end; ++it) {
             const keyT& key = it->first;
@@ -2491,7 +2491,7 @@ namespace madness {
             std::vector<Vector<double,NDIM> > newspecialpts;
             if (key.level() < functor->special_level() && specialpts.size() > 0) {
                 BoundaryConditions<NDIM> bc = FunctionDefaults<NDIM>::get_bc();
-                std::vector<bool> bperiodic = bc.is_periodic();
+                const auto bperiodic = bc.is_periodic();
                 for (unsigned int i = 0; i < specialpts.size(); ++i) {
                     coordT simpt;
                     user_to_sim(specialpts[i], simpt);
@@ -3239,7 +3239,7 @@ template <typename T, std::size_t NDIM>
     }
 
     template <typename T, std::size_t NDIM>
-    Key<NDIM> FunctionImpl<T,NDIM>::neighbor(const keyT& key, const Key<NDIM>& disp, const std::vector<bool>& is_periodic) const {
+    Key<NDIM> FunctionImpl<T,NDIM>::neighbor(const keyT& key, const Key<NDIM>& disp, const std::array<bool, NDIM>& is_periodic) const {
         Vector<Translation,NDIM> l = key.translation();
 
         for (std::size_t axis=0; axis<NDIM; ++axis) {