diff --git a/src/madness/mra/bc.h b/src/madness/mra/bc.h
index 3259e14d725..5fdc62393f2 100644
--- a/src/madness/mra/bc.h
+++ b/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,11 +126,14 @@ 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;
   }
 
@@ -138,11 +141,25 @@ template <std::size_t NDIM> class BoundaryConditions {
 
   /// @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
diff --git a/src/madness/mra/funcimpl.h b/src/madness/mra/funcimpl.h
index c57261e35cc..2985c44f6cb 100644
--- a/src/madness/mra/funcimpl.h
+++ b/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);
diff --git a/src/madness/mra/key.h b/src/madness/mra/key.h
index 47e2126c207..b2b730d2597 100644
--- a/src/madness/mra/key.h
+++ b/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)
diff --git a/src/madness/mra/leafop.h b/src/madness/mra/leafop.h
index ee9e63f9e2e..24cb439f8ff 100644
--- a/src/madness/mra/leafop.h
+++ b/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);
diff --git a/src/madness/mra/mraimpl.h b/src/madness/mra/mraimpl.h
index 52a739683e3..385470aafdd 100644
--- a/src/madness/mra/mraimpl.h
+++ b/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) {
diff --git a/src/madness/mra/operator.h b/src/madness/mra/operator.h
index 8f659d03b4c..e2ca2d6d3cd 100644
--- a/src/madness/mra/operator.h
+++ b/src/madness/mra/operator.h
@@ -141,8 +141,8 @@ namespace madness {
         OperatorInfo info;
 
         bool doleaves;  ///< If should be applied to leaf coefficients ... false by default
-        bool isperiodicsum;///< If true the operator 1D kernels have been summed over lattice translations
-                           ///< and may be non-zero at both ends of the unit cell
+        std::array<bool, NDIM> lattice_sum;  ///< If lattice_sum[d] is true summation over lattice translations along axis d
+                                             ///< N.B. the resulting kernel can be non-zero at both ends of the simulation cell along that axis
         bool modified_=false;     ///< use modified NS form
         int particle_=1;        ///< must only be 1 or 2
         bool destructive_=false;	///< destroy the argument or restore it (expensive for 6d functions)
@@ -159,7 +159,6 @@ namespace madness {
 
 
         mutable std::vector< ConvolutionND<Q,NDIM> > ops;   ///< ConvolutionND keeps data for 1 term, all dimensions, 1 displacement
-        const BoundaryConditions<NDIM> bc;
         const int k;
         const FunctionCommonData<Q,NDIM>& cdata;
         int rank;
@@ -209,10 +208,10 @@ namespace madness {
 
         static inline std::pair<Tensor<Q>,Tensor<Q>>
         make_coeff_for_operator(World& world, double mu, double lo, double eps, OpType type,
-                                const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc()) {
+                                const std::array<bool, NDIM>& lattice_sum) {
 
             OperatorInfo info(mu,lo,eps,type);
-            return make_coeff_for_operator(world, info, bc);
+            return make_coeff_for_operator(world, info, lattice_sum);
 //            const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
 //            double hi = cell_width.normf(); // Diagonal width of cell
 //            if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
@@ -233,26 +232,32 @@ namespace madness {
 
         static inline std::pair<Tensor<double>,Tensor<double>>
         make_coeff_for_operator(World& world, OperatorInfo& info,
-                                const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc()) {
-
-            const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
-            double hi = cell_width.normf(); // Diagonal width of cell
-            if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
-
-            info.hi=hi;
-            GFit<Q,NDIM> fit(info);
-
-            Tensor<Q> coeff=fit.coeffs();
-            Tensor<Q> expnt=fit.exponents();
-
-            // WARNING! More fine-grained control over the last argument is needed.
-            // This is a hotfix.
-            if (info.truncate_lowexp_gaussians.value_or(bc(0,0) == BC_PERIODIC)) {
-                fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), true);
-                info.truncate_lowexp_gaussians = true;
-            }
-
-            return std::make_pair(coeff,expnt);
+                                const std::array<bool, NDIM>& lattice_sum) {
+
+          const Tensor<double> &cell_width =
+              FunctionDefaults<3>::get_cell_width();
+          double hi = cell_width.normf(); // Diagonal width of cell
+          // Extend kernel range for lattice summation
+          const auto any_lattice_sum = std::accumulate(
+              lattice_sum.begin(), lattice_sum.end(), false, std::logical_or{});
+          if (any_lattice_sum) {
+            hi *= 100;
+          }
+
+          info.hi = hi;
+          GFit<Q, NDIM> fit(info);
+
+          Tensor<Q> coeff = fit.coeffs();
+          Tensor<Q> expnt = fit.exponents();
+
+          // WARNING! More fine-grained control over the last argument is needed. This is a hotfix.
+          if (info.truncate_lowexp_gaussians.value_or(any_lattice_sum)) {
+            fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(),
+                                            true);
+            info.truncate_lowexp_gaussians = true;
+          }
+
+          return std::make_pair(coeff, expnt);
         }
 
 //        /// return the right block of the upsampled operator (modified NS only)
@@ -947,17 +952,16 @@ namespace madness {
         // For separated convolutions with same operator in each direction (isotropic)
         SeparatedConvolution(World& world,
                              std::vector< std::shared_ptr< Convolution1D<Q> > >& argops,
-                             const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
+                             const std::array<bool, NDIM>& lattice_sum = FunctionDefaults<NDIM>::get_bc().is_periodic(),
                              long k = FunctionDefaults<NDIM>::get_k(),
                              bool doleaves = false)
                 : WorldObject< SeparatedConvolution<Q,NDIM> >(world)
                 , info()
                 , doleaves(doleaves)
-                , isperiodicsum(bc(0,0)==BC_PERIODIC)
+                , lattice_sum(lattice_sum)
                 , modified_(false)
                 , particle_(1)
                 , destructive_(false)
-                , bc(bc)
                 , k(k)
                 , cdata(FunctionCommonData<Q,NDIM>::get(k))
                 , rank(argops.size())
@@ -965,10 +969,6 @@ namespace madness {
                 , v2k(NDIM,2*k)
                 , s0(std::max<std::size_t>(2,NDIM),Slice(0,k-1))
         {
-            // Presently we must have periodic or non-periodic in all dimensions.
-            for (std::size_t d=1; d<NDIM; ++d) {
-                MADNESS_ASSERT(bc(d,0)==bc(0,0));
-            }
             check_cubic();
 
             for (unsigned int mu=0; mu < argops.size(); ++mu) {
@@ -981,18 +981,17 @@ namespace madness {
         // For general convolutions
         SeparatedConvolution(World& world,
                              std::vector< ConvolutionND<Q,NDIM> >& argops,
-                             const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
+                             const std::array<bool, NDIM>& lattice_sum = FunctionDefaults<NDIM>::get_bc().is_periodic(),
                              long k = FunctionDefaults<NDIM>::get_k(),
                              bool doleaves = false)
                 : WorldObject< SeparatedConvolution<Q,NDIM> >(world)
                 , info()
                 , doleaves(doleaves)
-                , isperiodicsum(bc(0,0)==BC_PERIODIC)
+                , lattice_sum(lattice_sum)
                 , modified_(false)
                 , particle_(1)
                 , destructive_(false)
                 , ops(argops)
-                , bc(bc)
                 , k(k)
                 , cdata(FunctionCommonData<Q,NDIM>::get(k))
                 , rank(argops.size())
@@ -1000,22 +999,18 @@ namespace madness {
                 , v2k(NDIM,2*k)
                 , s0(std::max<std::size_t>(2,NDIM),Slice(0,k-1))
         {
-            // Presently we must have periodic or non-periodic in all dimensions.
-            for (std::size_t d=1; d<NDIM; ++d) {
-                MADNESS_ASSERT(bc(d,0)==bc(0,0));
-            }
             this->process_pending();
         }
 
         /// Constructor for Gaussian Convolutions (mostly for backward compatability)
         SeparatedConvolution(World& world, const OperatorInfo info1,
-                             const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
+                             const std::array<bool, NDIM>& lattice_sum = FunctionDefaults<NDIM>::get_bc().is_periodic(),
                              int k=FunctionDefaults<NDIM>::get_k(),
                              bool doleaves = false)
-               : SeparatedConvolution(world,Tensor<double>(0l),Tensor<double>(0l),info1.lo,info1.thresh,bc,k,doleaves,info1.mu) {
+               : SeparatedConvolution(world,Tensor<double>(0l),Tensor<double>(0l),info1.lo,info1.thresh,lattice_sum,k,doleaves,info1.mu) {
             info.type=info1.type;
             info.truncate_lowexp_gaussians = info1.truncate_lowexp_gaussians;
-            auto [coeff, expnt] = make_coeff_for_operator(world, info, bc);
+            auto [coeff, expnt] = make_coeff_for_operator(world, info, lattice_sum);
             rank=coeff.dim(0);
             ops.resize(rank);
             initialize(coeff,expnt);
@@ -1025,16 +1020,15 @@ namespace madness {
         SeparatedConvolution(World& world,
                              const Tensor<Q>& coeff, const Tensor<double>& expnt,
                              double lo, double thresh,
-                             const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
+                             const std::array<bool, NDIM>& lattice_sum = FunctionDefaults<NDIM>::get_bc().is_periodic(),
                              int k=FunctionDefaults<NDIM>::get_k(),
                              bool doleaves = false,
                              double mu=0.0)
                 : WorldObject< SeparatedConvolution<Q,NDIM> >(world)
                 , info(mu,lo,thresh,OT_UNDEFINED)
                 , doleaves(doleaves)
-                , isperiodicsum(bc(0,0)==BC_PERIODIC)
+                , lattice_sum(lattice_sum)
                 , ops(coeff.dim(0))
-                , bc(bc)
                 , k(k)
                 , cdata(FunctionCommonData<Q,NDIM>::get(k))
                 , rank(coeff.dim(0))
@@ -1045,11 +1039,6 @@ namespace madness {
         }
 
         void initialize(const Tensor<Q>& coeff, const Tensor<double>& expnt) {
-            // Presently we must have periodic or non-periodic in all dimensions.
-            for (std::size_t d=1; d<NDIM; ++d) {
-                MADNESS_ASSERT(bc(d,0)==bc(0,0));
-            }
-
             const Tensor<double>& width = FunctionDefaults<NDIM>::get_cell_width();
             const double pi = constants::pi;
 
@@ -1061,7 +1050,7 @@ namespace madness {
                 ops[mu].setfac(coeff(mu)/c);
 
                 for (std::size_t d=0; d<NDIM; ++d) {
-                  ops[mu].setop(d,GaussianConvolution1DCache<Q>::get(k, expnt(mu)*width[d]*width[d], 0, isperiodicsum));
+                  ops[mu].setop(d,GaussianConvolution1DCache<Q>::get(k, expnt(mu)*width[d]*width[d], 0, lattice_sum[d]));
                 }
             }
         }
@@ -1070,18 +1059,17 @@ namespace madness {
         SeparatedConvolution(World& world,
                              Vector<double,NDIM> args,
                              const Tensor<Q>& coeff, const Tensor<double>& expnt,
-                             const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
+                             const std::array<bool, NDIM>& lattice_sum = FunctionDefaults<NDIM>::get_bc().is_periodic(),
                              int k=FunctionDefaults<NDIM>::get_k(),
                              bool doleaves=false)
                 : WorldObject< SeparatedConvolution<Q,NDIM> >(world)
                 , info(0.0,0.0,0.0,OT_UNDEFINED)
                 , doleaves(doleaves)
-                , isperiodicsum(bc(0,0)==BC_PERIODIC)
+                , lattice_sum(lattice_sum)
                 , modified_(false)
                 , particle_(1)
                 , destructive_(false)
                 , ops(coeff.dim(0))
-                , bc(bc)
                 , k(k)
                 , cdata(FunctionCommonData<Q,NDIM>::get(k))
                 , rank(coeff.dim(0))
@@ -1089,11 +1077,6 @@ namespace madness {
                 , v2k(NDIM,2*k)
                 , s0(std::max<std::size_t>(2,NDIM),Slice(0,k-1))
         {
-            // Presently we must have periodic or non-periodic in all dimensions.
-            for (std::size_t d=1; d<NDIM; ++d) {
-                MADNESS_ASSERT(bc(d,0)==bc(0,0));
-            }
-
             const Tensor<double>& width = FunctionDefaults<NDIM>::get_cell_width();
 
             for (int mu=0; mu<rank; ++mu) {
@@ -1103,7 +1086,7 @@ namespace madness {
                   double c2 = sqrt(expnt[mu]*width[d]*width[d]/madness::constants::pi);
                   std::shared_ptr<GaussianConvolution1D<double_complex> >
                       gcptr(new GaussianConvolution1D<double_complex>(k, c2, 
-                            expnt(mu)*width[d]*width[d], 0, isperiodicsum, args[d]));
+                            expnt(mu)*width[d]*width[d], 0, lattice_sum[d], args[d]));
                   ops[mu].setop(d,gcptr);
                 }
             }
@@ -1127,13 +1110,11 @@ namespace madness {
         	}
         }
 
-        const BoundaryConditions<NDIM>& get_bc() const {return bc;}
-
         const std::vector< Key<NDIM> >& get_disp(Level n) const {
-            return Displacements<NDIM>().get_disp(n, isperiodicsum);
+            return Displacements<NDIM>().get_disp(n, false);
         }
 
-        bool get_isperiodicsum() const {return isperiodicsum;}
+        const std::array<bool, NDIM>& includes_lattice_sum() const { return lattice_sum; }
 
         /// return the operator norm for all terms, all dimensions and 1 displacement
         double norm(Level n, const Key<NDIM>& d, const Key<NDIM>& source_key) const {
@@ -1697,9 +1678,10 @@ namespace madness {
                                                     const SeparatedConvolution<Q,NDIM>& right) {
             MADNESS_CHECK(can_combine(left,right));
             MADNESS_CHECK(left.get_world().id()==right.get_world().id());
+            MADNESS_CHECK(left.includes_lattice_sum() == right.includes_lattice_sum());
 
             auto info=combine_OT(left,right);
-            return SeparatedConvolution<Q,NDIM>(left.get_world(),info,left.bc,left.k);
+            return SeparatedConvolution<Q,NDIM>(left.get_world(),info,left.includes_lattice_sum(),left.k);
         }
 
         /// combine 2 convolution operators to one
@@ -1728,24 +1710,28 @@ namespace madness {
                                                    Vector<double,3> args,
                                                    double lo,
                                                    double eps,
-                                                   const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
-                                                   int k=FunctionDefaults<3>::get_k())
-    {
-        const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
-        double hi = cell_width.normf(); // Diagonal width of cell
-        if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
-
-        GFit<double,3> fit=GFit<double,3>::CoulombFit(lo,hi,eps,false);
-		Tensor<double> coeff=fit.coeffs();
-		Tensor<double> expnt=fit.exponents();
-
-		if (bc(0,0) == BC_PERIODIC) {
-        	fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), true);
-        }
-
-        return SeparatedConvolution<double_complex,3>(world, args, coeff, expnt, bc, k, false);
-//        return SeparatedConvolution<double_complex,3>(world, coeff, expnt, bc, k);
-
+                                                   const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
+                                                   int k=FunctionDefaults<3>::get_k()) {
+      const Tensor<double> &cell_width = FunctionDefaults<3>::get_cell_width();
+      double hi = cell_width.normf(); // Diagonal width of cell
+
+      // Extend kernel range for lattice summation
+      const auto any_lattice_sum = std::accumulate(
+          lattice_sum.begin(), lattice_sum.end(), false, std::logical_or{});
+      if (any_lattice_sum) {
+        hi *= 100;
+      }
+
+      GFit<double, 3> fit = GFit<double, 3>::CoulombFit(lo, hi, eps, false);
+      Tensor<double> coeff = fit.coeffs();
+      Tensor<double> expnt = fit.exponents();
+
+      if (any_lattice_sum) {
+        fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), true);
+      }
+
+      return SeparatedConvolution<double_complex, 3>(world, args, coeff, expnt,
+                                                     lattice_sum, k, false);
     }
 
     /// Factory function generating separated kernel for convolution with 1/r in 3D.
@@ -1754,11 +1740,10 @@ namespace madness {
     SeparatedConvolution<double,3> CoulombOperator(World& world,
                                                    double lo,
                                                    double eps,
-                                                   const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
-
-												   int k=FunctionDefaults<3>::get_k())
+                                                   const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
+						   int k=FunctionDefaults<3>::get_k())
     {
-        return SeparatedConvolution<double,3>(world,OperatorInfo(0.0,lo,eps,OT_G12),bc,k);
+        return SeparatedConvolution<double,3>(world,OperatorInfo(0.0,lo,eps,OT_G12),lattice_sum,k);
     }
 
 
@@ -1768,10 +1753,10 @@ namespace madness {
     SeparatedConvolution<double,3>* CoulombOperatorPtr(World& world,
                                                        double lo,
                                                        double eps,
-                                                       const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
+                                                       const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
                                                        int k=FunctionDefaults<3>::get_k())
     {
-        return new SeparatedConvolution<double,3>(world,OperatorInfo(0.0,lo,eps,OT_G12),bc,k);
+        return new SeparatedConvolution<double,3>(world,OperatorInfo(0.0,lo,eps,OT_G12),lattice_sum,k);
     }
 
 
@@ -1780,13 +1765,13 @@ namespace madness {
     static inline
     SeparatedConvolution<double,NDIM>
     BSHOperator(World& world, double mu, double lo, double eps,
-                const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(),
+                const std::array<bool, NDIM>& lattice_sum = FunctionDefaults<NDIM>::get_bc().is_periodic(),
                 int k=FunctionDefaults<NDIM>::get_k()) {
     	if (eps>1.e-4) {
     		if (world.rank()==0) print("the accuracy in BSHOperator is too small, tighten the threshold",eps);
     		MADNESS_EXCEPTION("0",1);
     	}
-        return SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_BSH),bc,k);
+        return SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_BSH),lattice_sum,k);
     }
 
     /// Factory function generating separated kernel for convolution with BSH kernel in general NDIM
@@ -1794,22 +1779,22 @@ namespace madness {
     static inline
     SeparatedConvolution<double,NDIM>*
     BSHOperatorPtr(World& world, double mu, double lo, double eps,
-                   const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(),
+                   const std::array<bool, NDIM>& lattice_sum = FunctionDefaults<NDIM>::get_bc().is_periodic(),
                    int k=FunctionDefaults<NDIM>::get_k()) {
         if (eps>1.e-4) {
             if (world.rank()==0) print("the accuracy in BSHOperator is too small, tighten the threshold",eps);
             MADNESS_EXCEPTION("0",1);
         }
-        return new SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_BSH),bc,k);
+        return new SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_BSH),lattice_sum,k);
     }
 
 
     /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D
     static inline SeparatedConvolution<double,3>
     BSHOperator3D(World& world, double mu, double lo, double eps,
-                  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
+                  const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
                   int k=FunctionDefaults<3>::get_k()) {
-        return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_BSH),bc,k);
+        return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_BSH),lattice_sum,k);
     }
 
     /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D
@@ -1820,22 +1805,28 @@ namespace madness {
                                                          double mu,
                                                          double lo,
                                                          double eps,
-                                                         const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
+                                                         const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
                                                          int k=FunctionDefaults<3>::get_k())
 
     {
-        const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
-        double hi = cell_width.normf(); // Diagonal width of cell
-        if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
-
-        GFit<double,3> fit=GFit<double,3>::BSHFit(mu,lo,hi,eps,false);
-		Tensor<double> coeff=fit.coeffs();
-		Tensor<double> expnt=fit.exponents();
-
-        if (bc(0,0) == BC_PERIODIC) {
-            fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false);
-        }
-        return SeparatedConvolution<double_complex,3>(world, args, coeff, expnt, bc, k);
+      const Tensor<double> &cell_width = FunctionDefaults<3>::get_cell_width();
+      double hi = cell_width.normf(); // Diagonal width of cell
+      // Extend kernel range for lattice summation
+      const auto any_lattice_sum = std::accumulate(
+          lattice_sum.begin(), lattice_sum.end(), false, std::logical_or{});
+      if (any_lattice_sum) {
+        hi *= 100;
+      }
+
+      GFit<double, 3> fit = GFit<double, 3>::BSHFit(mu, lo, hi, eps, false);
+      Tensor<double> coeff = fit.coeffs();
+      Tensor<double> expnt = fit.exponents();
+
+      if (any_lattice_sum) {
+        fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false);
+      }
+      return SeparatedConvolution<double_complex, 3>(world, args, coeff, expnt,
+                                                     lattice_sum, k);
     }
 
     /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D
@@ -1845,52 +1836,59 @@ namespace madness {
                                                          double mu,
                                                          double lo,
                                                          double eps,
-                                                         const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
+                                                         const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
                                                          int k=FunctionDefaults<3>::get_k())
 
     {
-        const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
-        double hi = cell_width.normf(); // Diagonal width of cell
-        if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
-
-        GFit<double,3> fit=GFit<double,3>::BSHFit(mu,lo,hi,eps,false);
-		Tensor<double> coeff=fit.coeffs();
-		Tensor<double> expnt=fit.exponents();
-
-        if (bc(0,0) == BC_PERIODIC) {
-            fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false);
-        }
-        return new SeparatedConvolution<double_complex,3>(world, args, coeff, expnt, bc, k);
+      const Tensor<double> &cell_width = FunctionDefaults<3>::get_cell_width();
+      double hi = cell_width.normf(); // Diagonal width of cell
+      // Extend kernel range for lattice summation
+      const auto any_lattice_sum = std::accumulate(
+          lattice_sum.begin(), lattice_sum.end(), false, std::logical_or{});
+      if (any_lattice_sum) {
+        hi *= 100;
+      }
+
+      GFit<double, 3> fit = GFit<double, 3>::BSHFit(mu, lo, hi, eps, false);
+      Tensor<double> coeff = fit.coeffs();
+      Tensor<double> expnt = fit.exponents();
+
+      if (any_lattice_sum) {
+        fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false);
+      }
+      return new SeparatedConvolution<double_complex, 3>(world, args, coeff,
+                                                         expnt, lattice_sum, k);
     }
 
 
     static inline SeparatedConvolution<double,3>
     SlaterF12Operator(World& world, double mu, double lo, double eps,
-                      const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), int k=FunctionDefaults<3>::get_k()) {
-        return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F12),bc,k);
+                      const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
+                      int k=FunctionDefaults<3>::get_k()) {
+        return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F12),lattice_sum,k);
     }
 
     static inline SeparatedConvolution<double,3> SlaterF12sqOperator(World& world,
                                                                    double mu, double lo, double eps,
-                                                                   const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
+                                                                   const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
                                                                    int k=FunctionDefaults<3>::get_k()) {
-        return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F212),bc,k);
+        return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F212),lattice_sum,k);
     }
 
     static inline SeparatedConvolution<double,3>* SlaterF12sqOperatorPtr(World& world,
                                                                        double mu, double lo, double eps,
-                                                                       const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
+                                                                       const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
                                                                        int k=FunctionDefaults<3>::get_k()) {
-        return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F212),bc,k);
+        return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F212),lattice_sum,k);
     }
 
     /// Factory function generating separated kernel for convolution with exp(-mu*r) in 3D
     template<std::size_t NDIM=3>
     static inline SeparatedConvolution<double,NDIM> SlaterOperator(World& world,
     		double mu, double lo, double eps,
-    		const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(),
+                const std::array<bool, NDIM>& lattice_sum = FunctionDefaults<NDIM>::get_bc().is_periodic(),
     		int k=FunctionDefaults<NDIM>::get_k()) {
-        return SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_SLATER),bc,k);
+        return SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_SLATER),lattice_sum,k);
     }
 
     /// Factory function generating separated kernel for convolution with exp(-mu*r*r)
@@ -1899,9 +1897,9 @@ namespace madness {
     template<std::size_t NDIM>
     static inline SeparatedConvolution<double,NDIM> GaussOperator(World& world,
                                                                 double mu, double lo=0.0, double eps=0.0,
-                                                                const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(),
+                                                                const std::array<bool, NDIM>& lattice_sum = FunctionDefaults<NDIM>::get_bc().is_periodic(),
                                                                 int k=FunctionDefaults<NDIM>::get_k()) {
-        return SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_GAUSS),bc,k);
+        return SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_GAUSS),lattice_sum,k);
     }
 
     /// Factory function generating separated kernel for convolution with exp(-mu*r*r) in 3D
@@ -1910,9 +1908,9 @@ namespace madness {
     template<std::size_t NDIM>
     static inline SeparatedConvolution<double, NDIM>* GaussOperatorPtr(World& world,
                                                                      double mu, double lo=0.0, double eps=0.0,
-                                                                     const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
+                                                                     const std::array<bool, NDIM>& lattice_sum = FunctionDefaults<NDIM>::get_bc().is_periodic(),
                                                                      int k = FunctionDefaults<NDIM>::get_k()) {
-        return new SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_GAUSS),bc,k);
+        return new SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_GAUSS),lattice_sum,k);
     }
 
 
@@ -1921,18 +1919,18 @@ namespace madness {
     template<std::size_t NDIM>
     static inline SeparatedConvolution<double, NDIM>* SlaterOperatorPtr_ND(World& world,
                                                                      double mu, double lo, double eps,
-                                                                     const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
+                                                                     const std::array<bool, NDIM>& lattice_sum = FunctionDefaults<NDIM>::get_bc().is_periodic(),
                                                                      int k = FunctionDefaults<NDIM>::get_k()) {
-        return new SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_SLATER),bc,k);
+        return new SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_SLATER),lattice_sum,k);
     }
 
     /// Factory function generating separated kernel for convolution with exp(-mu*r) in 3D
     /// Note that the 1/(2mu) factor of SlaterF12Operator is not included, this is just the exponential function
     static inline SeparatedConvolution<double, 3>* SlaterOperatorPtr(World& world,
                                                                  double mu, double lo, double eps,
-                                                                 const BoundaryConditions<3>& bc = FunctionDefaults<3>::get_bc(),
+                                                                 const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
                                                                  int k = FunctionDefaults<3>::get_k()) {
-        return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_SLATER),bc,k);
+        return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_SLATER),lattice_sum,k);
     }
 
     /// Factory function generating separated kernel for convolution with (1 - exp(-mu*r))/(2 mu) in 3D
@@ -1940,9 +1938,9 @@ namespace madness {
     /// includes the factor 1/(2 mu)
     static inline SeparatedConvolution<double,3>* SlaterF12OperatorPtr(World& world,
     		double mu, double lo, double eps,
-    		const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
+                const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
     		int k=FunctionDefaults<3>::get_k()) {
-        return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F12),bc,k);
+        return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F12),lattice_sum,k);
     }
 
 
@@ -1952,9 +1950,9 @@ namespace madness {
     /// includes the factor 1/(2 mu)
     static inline SeparatedConvolution<double,3>
     FGOperator(World& world, double mu, double lo, double eps,
-               const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
+               const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
                int k=FunctionDefaults<3>::get_k()) {
-        return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_FG12),bc,k);
+        return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_FG12),lattice_sum,k);
     }
 
     /// Factory function generating separated kernel for convolution with 1/(2 mu)*(1 - exp(-mu*r))/r in 3D
@@ -1963,9 +1961,9 @@ namespace madness {
     /// includes the factor 1/(2 mu)
     static inline SeparatedConvolution<double,3>*
     FGOperatorPtr(World& world, double mu, double lo, double eps,
-                  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
+                  const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
                   int k=FunctionDefaults<3>::get_k()) {
-        return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_FG12),bc,k);
+        return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_FG12),lattice_sum,k);
     }
 
     /// Factory function generating separated kernel for convolution with (1/(2 mu)*(1 - exp(-mu*r)))^2/r in 3D
@@ -1975,9 +1973,9 @@ namespace madness {
     /// includes the factor 1/(2 mu)^2
     static inline SeparatedConvolution<double,3>*
     F2GOperatorPtr(World& world, double mu, double lo, double eps,
-                   const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
+                   const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
                    int k=FunctionDefaults<3>::get_k()) {
-        return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F2G12),bc,k);
+        return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F2G12),lattice_sum,k);
     }
 
     /// Factory function generating separated kernel for convolution with (1/(2 mu)*(1 - exp(-mu*r)))^2/r in 3D
@@ -1987,9 +1985,9 @@ namespace madness {
     /// includes the factor 1/(2 mu)^2
     static inline SeparatedConvolution<double,3>
     F2GOperator(World& world, double mu, double lo, double eps,
-                const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
+                const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
                 int k=FunctionDefaults<3>::get_k()) {
-        return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F2G12),bc,k);
+        return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F2G12),lattice_sum,k);
     }
 
 
@@ -1997,15 +1995,14 @@ namespace madness {
     /// Gaussian (aka a widened delta function)
     static inline SeparatedConvolution<double,3> SmoothingOperator3D(World& world,
             double eps,
-            const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
+            const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
             int k=FunctionDefaults<3>::get_k()) {
 
         double exponent = 1.0/(2.0*eps);
         Tensor<double> coeffs(1), exponents(1);
         exponents(0L) =  exponent;
         coeffs(0L)=pow(exponent/M_PI,0.5*3.0);  // norm of the gaussian
-        return SeparatedConvolution<double,3>(world, coeffs, exponents, 1.e-8, eps);
-
+        return SeparatedConvolution<double,3>(world, coeffs, exponents, 1.e-8, eps, lattice_sum, k);
     }
 
     /// Factory function generating separated kernel for convolution a normalized
@@ -2013,14 +2010,14 @@ namespace madness {
     template<std::size_t NDIM>
     static inline SeparatedConvolution<double,NDIM> SmoothingOperator(World& world,
             double eps,
-            const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(),
+            const std::array<bool, NDIM>& lattice_sum = FunctionDefaults<NDIM>::get_bc().is_periodic(),
             int k=FunctionDefaults<NDIM>::get_k()) {
 
         double exponent = 1.0/(2.0*eps);
         Tensor<double> coeffs(1), exponents(1);
         exponents(0L) =  exponent;
         coeffs(0L)=pow(exponent/M_PI,0.5*NDIM);  // norm of the gaussian
-        return SeparatedConvolution<double,NDIM>(world, coeffs, exponents, 1.e-8, eps);
+        return SeparatedConvolution<double,NDIM>(world, coeffs, exponents, 1.e-8, eps, lattice_sum, k);
     }
 
     /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D
@@ -2030,21 +2027,26 @@ namespace madness {
                                                      double mu,
                                                      double lo,
                                                      double eps,
-                                                     const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
-                                                     int k=FunctionDefaults<3>::get_k())
-    {
-        const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
-        double hi = cell_width.normf(); // Diagonal width of cell
-        if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
-
-        GFit<double,3> fit=GFit<double,3>::BSHFit(mu,lo,hi,eps,false);
-		Tensor<double> coeff=fit.coeffs();
-		Tensor<double> expnt=fit.exponents();
-
-        if (bc(0,0) == BC_PERIODIC) {
-            fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false);
-        }
-        return new SeparatedConvolution<double,3>(world, coeff, expnt, lo, eps, bc, k);
+                                                     const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
+                                                     int k=FunctionDefaults<3>::get_k()) {
+      const Tensor<double> &cell_width = FunctionDefaults<3>::get_cell_width();
+      double hi = cell_width.normf(); // Diagonal width of cell
+      // Extend kernel range for lattice summation
+      const auto any_lattice_sum = std::accumulate(
+          lattice_sum.begin(), lattice_sum.end(), false, std::logical_or{});
+      if (any_lattice_sum) {
+        hi *= 100;
+      }
+
+      GFit<double, 3> fit = GFit<double, 3>::BSHFit(mu, lo, hi, eps, false);
+      Tensor<double> coeff = fit.coeffs();
+      Tensor<double> expnt = fit.exponents();
+
+      if (any_lattice_sum) {
+        fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false);
+      }
+      return new SeparatedConvolution<double, 3>(world, coeff, expnt, lo, eps,
+                                                 lattice_sum, k);
     }
 
 
@@ -2058,47 +2060,54 @@ namespace madness {
     GradCoulombOperator(World& world,
                         double lo,
                         double eps,
-                        const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
-                        int k=FunctionDefaults<3>::get_k())
-    {
-        typedef SeparatedConvolution<double,3> real_convolution_3d;
-        typedef std::shared_ptr<real_convolution_3d> real_convolution_3d_ptr;
-        const double pi = constants::pi;
-        const Tensor<double> width = FunctionDefaults<3>::get_cell_width();
-        double hi = width.normf(); // Diagonal width of cell
-        const bool isperiodicsum = (bc(0,0)==BC_PERIODIC);
-        if (isperiodicsum) hi *= 100; // Extend range for periodic summation
-
-        GFit<double,3> fit=GFit<double,3>::CoulombFit(lo,hi,eps,false);
-		Tensor<double> coeff=fit.coeffs();
-		Tensor<double> expnt=fit.exponents();
-
-        if (bc(0,0) == BC_PERIODIC) {
-            fit.truncate_periodic_expansion(coeff, expnt, width.max(), true);
+                        const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
+                        int k=FunctionDefaults<3>::get_k()) {
+      typedef SeparatedConvolution<double, 3> real_convolution_3d;
+      typedef std::shared_ptr<real_convolution_3d> real_convolution_3d_ptr;
+      const double pi = constants::pi;
+      const Tensor<double> width = FunctionDefaults<3>::get_cell_width();
+      double hi = width.normf(); // Diagonal width of cell
+      // Extend kernel range for lattice summation
+      const auto any_lattice_sum = std::accumulate(
+          lattice_sum.begin(), lattice_sum.end(), false, std::logical_or{});
+      if (any_lattice_sum) {
+        hi *= 100;
+      }
+
+      GFit<double, 3> fit = GFit<double, 3>::CoulombFit(lo, hi, eps, false);
+      Tensor<double> coeff = fit.coeffs();
+      Tensor<double> expnt = fit.exponents();
+
+      if (any_lattice_sum) {
+        fit.truncate_periodic_expansion(coeff, expnt, width.max(), true);
+      }
+
+      int rank = coeff.dim(0);
+
+      std::vector<real_convolution_3d_ptr> gradG(3);
+
+      for (int dir = 0; dir < 3; dir++) {
+        std::vector<ConvolutionND<double, 3>> ops(rank);
+        for (int mu = 0; mu < rank; mu++) {
+          // We cache the normalized operator so the factor is the value we must multiply by to recover the coeff we want.
+          double c = std::pow(sqrt(expnt(mu) / pi), 3); // Normalization coeff
+          ops[mu].setfac(coeff(mu) / c / width[dir]);
+
+          for (int d = 0; d < 3; d++) {
+            if (d != dir)
+              ops[mu].setop(d, GaussianConvolution1DCache<double>::get(
+                                   k, expnt(mu) * width[d] * width[d], 0,
+                                   lattice_sum[d]));
+          }
+          ops[mu].setop(dir, GaussianConvolution1DCache<double>::get(
+                                 k, expnt(mu) * width[dir] * width[dir], 1,
+                                 lattice_sum[dir]));
         }
+        gradG[dir] = real_convolution_3d_ptr(
+            new SeparatedConvolution<double, 3>(world, ops));
+      }
 
-        int rank = coeff.dim(0);
-
-        std::vector<real_convolution_3d_ptr> gradG(3);
-
-        for (int dir=0; dir<3; dir++) {
-            std::vector< ConvolutionND<double,3> > ops(rank);
-            for (int mu=0; mu<rank; mu++) {
-                // We cache the normalized operator so the factor is the value we must multiply
-                // by to recover the coeff we want.
-                double c = std::pow(sqrt(expnt(mu)/pi),3); // Normalization coeff
-                ops[mu].setfac(coeff(mu)/c/width[dir]);
-
-                for (int d=0; d<3; d++) {
-                    if (d != dir)
-                        ops[mu].setop(d,GaussianConvolution1DCache<double>::get(k, expnt(mu)*width[d]*width[d], 0, isperiodicsum));
-                }
-                ops[mu].setop(dir,GaussianConvolution1DCache<double>::get(k, expnt(mu)*width[dir]*width[dir], 1, isperiodicsum));
-            }
-            gradG[dir] = real_convolution_3d_ptr(new SeparatedConvolution<double,3>(world, ops));
-        }
-
-        return gradG;
+      return gradG;
     }
 
     /// Factory function generating operator for convolution with grad(bsh) in 3D
@@ -2112,47 +2121,54 @@ namespace madness {
                         double mu,
                         double lo,
                         double eps,
-                        const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
-                        int k=FunctionDefaults<3>::get_k())
-    {
-        typedef SeparatedConvolution<double,3> real_convolution_3d;
-        typedef std::shared_ptr<real_convolution_3d> real_convolution_3d_ptr;
-        const double pi = constants::pi;
-        const Tensor<double> width = FunctionDefaults<3>::get_cell_width();
-        double hi = width.normf(); // Diagonal width of cell
-        const bool isperiodicsum = (bc(0,0)==BC_PERIODIC);
-        if (isperiodicsum) hi *= 100; // Extend range for periodic summation
-
-        GFit<double,3> fit=GFit<double,3>::BSHFit(mu,lo,hi,eps,false);
-        Tensor<double> coeff=fit.coeffs();
-        Tensor<double> expnt=fit.exponents();
-
-        if (bc(0,0) == BC_PERIODIC) {
-            fit.truncate_periodic_expansion(coeff, expnt, width.max(), true);
-        }
-
-        int rank = coeff.dim(0);
-
-        std::vector<real_convolution_3d_ptr> gradG(3);
-
-        for (int dir=0; dir<3; dir++) {
-            std::vector< ConvolutionND<double,3> > ops(rank);
-            for (int mu=0; mu<rank; mu++) {
-                // We cache the normalized operator so the factor is the value we must multiply
-                // by to recover the coeff we want.
-                double c = std::pow(sqrt(expnt(mu)/pi),3); // Normalization coeff
-                ops[mu].setfac(coeff(mu)/c/width[dir]);
-
-                for (int d=0; d<3; d++) {
-                    if (d != dir)
-                        ops[mu].setop(d,GaussianConvolution1DCache<double>::get(k, expnt(mu)*width[d]*width[d], 0, isperiodicsum));
-                }
-                ops[mu].setop(dir,GaussianConvolution1DCache<double>::get(k, expnt(mu)*width[dir]*width[dir], 1, isperiodicsum));
-            }
-            gradG[dir] = real_convolution_3d_ptr(new SeparatedConvolution<double,3>(world, ops));
+                        const std::array<bool, 3>& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(),
+                        int k=FunctionDefaults<3>::get_k()) {
+      typedef SeparatedConvolution<double, 3> real_convolution_3d;
+      typedef std::shared_ptr<real_convolution_3d> real_convolution_3d_ptr;
+      const double pi = constants::pi;
+      const Tensor<double> width = FunctionDefaults<3>::get_cell_width();
+      double hi = width.normf(); // Diagonal width of cell
+      // Extend kernel range for lattice summation
+      const auto any_lattice_sum = std::accumulate(
+          lattice_sum.begin(), lattice_sum.end(), false, std::logical_or{});
+      if (any_lattice_sum) {
+        hi *= 100;
+      }
+
+      GFit<double, 3> fit = GFit<double, 3>::BSHFit(mu, lo, hi, eps, false);
+      Tensor<double> coeff = fit.coeffs();
+      Tensor<double> expnt = fit.exponents();
+
+      if (any_lattice_sum) {
+        fit.truncate_periodic_expansion(coeff, expnt, width.max(), true);
+      }
+
+      int rank = coeff.dim(0);
+
+      std::vector<real_convolution_3d_ptr> gradG(3);
+
+      for (int dir = 0; dir < 3; dir++) {
+        std::vector<ConvolutionND<double, 3>> ops(rank);
+        for (int mu = 0; mu < rank; mu++) {
+          // We cache the normalized operator so the factor is the value we must multiply by to recover the coeff we want.
+          double c = std::pow(sqrt(expnt(mu) / pi), 3); // Normalization coeff
+          ops[mu].setfac(coeff(mu) / c / width[dir]);
+
+          for (int d = 0; d < 3; d++) {
+            if (d != dir)
+              ops[mu].setop(d, GaussianConvolution1DCache<double>::get(
+                                   k, expnt(mu) * width[d] * width[d], 0,
+                                   lattice_sum[d]));
+          }
+          ops[mu].setop(dir, GaussianConvolution1DCache<double>::get(
+                                 k, expnt(mu) * width[dir] * width[dir], 1,
+                                 lattice_sum[dir]));
         }
+        gradG[dir] = real_convolution_3d_ptr(
+            new SeparatedConvolution<double, 3>(world, ops));
+      }
 
-        return gradG;
+      return gradG;
     }
 
 
diff --git a/src/madness/mra/qmprop.cc b/src/madness/mra/qmprop.cc
index 7108480f906..f4a61536b4a 100644
--- a/src/madness/mra/qmprop.cc
+++ b/src/madness/mra/qmprop.cc
@@ -172,7 +172,7 @@ namespace madness {
         double width = FunctionDefaults<NDIM>::get_cell_min_width(); // Assuming cubic so all dim equal
         std::vector< std::shared_ptr< Convolution1D<double_complex> > > q(1);
         q[0].reset(qm_1d_free_particle_propagator(k, bandlimit, timestep, width));
-        return SeparatedConvolution<double_complex,NDIM>(world, q, BoundaryConditions<NDIM>(BC_FREE), k, true);
+        return SeparatedConvolution<double_complex,NDIM>(world, q, no_lattice_sum<NDIM>(), k, true);
     }
 
     template <std::size_t NDIM>
@@ -181,7 +181,7 @@ namespace madness {
         double width = FunctionDefaults<NDIM>::get_cell_min_width(); // Assuming cubic so all dim equal
         std::vector< std::shared_ptr< Convolution1D<double_complex> > > q(1);
         q[0].reset(qm_1d_free_particle_propagator(k, bandlimit, timestep, width));
-        return new SeparatedConvolution<double_complex,NDIM>(world, q, BoundaryConditions<NDIM>(BC_FREE), k, true);
+        return new SeparatedConvolution<double_complex,NDIM>(world, q, no_lattice_sum<NDIM>(), k, true);
     }
 
 #ifdef FUNCTION_INSTANTIATE_1