operator.h

Go to the documentation of this file.

/*
  This file is part of MADNESS.
 
  Copyright (C) 2007,2010 Oak Ridge National Laboratory
 
  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 
  For more information please contact:
 
  Robert J. Harrison
  Oak Ridge National Laboratory
  One Bethel Valley Road
  P.O. Box 2008, MS-6367
 
  email: harrisonrj@ornl.gov
  tel:   865-241-3937
  fax:   865-572-0680
*/
 
#ifndef MADNESS_MRA_OPERATOR_H__INCLUDED
#define MADNESS_MRA_OPERATOR_H__INCLUDED
 
/// \file mra/operator.h
/// \brief Implements most functionality of separated operators
 
/// \ingroup function
 
#include <type_traits>
#include <limits.h>
#include <madness/mra/adquad.h>
#include <madness/tensor/aligned.h>
#include <madness/tensor/tensor_lapack.h>
#include <madness/constants.h>
 
#include <madness/mra/simplecache.h>
#include <madness/mra/convolution1d.h>
#include <madness/mra/displacements.h>
#include <madness/mra/function_common_data.h>
#include <madness/mra/gfit.h>
#include <madness/mra/operatorinfo.h>
 
namespace madness {
 
    template<typename T, std::size_t NDIM>
    class Function;
 
    template<typename T, std::size_t NDIM>
    class SeparatedConvolution;
 
    template<typename T, std::size_t NDIM>
    class CCPairFunction;
 
    template <typename T, typename R, std::size_t NDIM, std::size_t KDIM>
    std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
    apply(const SeparatedConvolution<T,KDIM>& op, const std::vector< Function<R,NDIM> > f);
 
    template<typename T, std::size_t NDIM>
    std::vector<CCPairFunction<T,NDIM>> apply(const SeparatedConvolution<T,NDIM>& op, const std::vector<CCPairFunction<T,NDIM>>& argument);
 
    template<typename T, std::size_t NDIM>
    std::vector<CCPairFunction<T,NDIM>> apply(const SeparatedConvolution<T,NDIM/2>& op, const std::vector<CCPairFunction<T,NDIM>>& argument);
 
    template<typename T, std::size_t NDIM>
    CCPairFunction<T,NDIM> apply(const SeparatedConvolution<T,NDIM>& op, const CCPairFunction<T,NDIM>& argument);
 
    template<typename T, std::size_t NDIM>
    CCPairFunction<T,NDIM> apply(const SeparatedConvolution<T,NDIM/2>& op, const CCPairFunction<T,NDIM>& argument);
 
    /// SeparatedConvolutionInternal keeps data for 1 term and all dimensions and 1 displacement
    /// Why is this here?? Why don't you just use ConvolutionND in SeparatedConvolutionData??
    template <typename Q, std::size_t NDIM>
    struct SeparatedConvolutionInternal {
        double norm;
        const ConvolutionData1D<Q>* ops[NDIM];
    };
 
    /// SeparatedConvolutionData keeps data for all terms, all dimensions
 
    /// this struct is used to cache the data that are generated by
    template <typename Q, std::size_t NDIM>
    struct SeparatedConvolutionData {
        std::vector< SeparatedConvolutionInternal<Q,NDIM> > muops;
        double norm;
 
        SeparatedConvolutionData(int rank) : muops(rank), norm(0.0) {}
        SeparatedConvolutionData(const SeparatedConvolutionData<Q,NDIM>& q) {
            muops = q.muops;
            norm = q.norm;
        }
    };
 
 
    /// Convolutions in separated form (including Gaussian)
 
    /* this stuff is very confusing, poorly commented, and extremely poorly named!
 
    I think it works like this:
    We try to apply transition matrices to the compressed form of function coefficients.
    Most of the code is about caching these transition matrices. They are cached (key of the map is the displacement)
    in the SimpleCache "data", which is of type SeparatedConvolutionData, which keeps the matrices
    for all separated terms and dimensions. These SeparatedConvolutionData are constructed using
    ConvolutionND "ops", which is constructed at the construction of the SeparatedConvolution.
 
                        SeparatedConvolution (all terms, all dim, all displacements)
 
                    construction                                            storage
 
                                                                    SimpleCache<SeparatedConvolutionData>
                                                                        (all terms, all dim) / (all disp)
             vector<ConvolutionND>
             (1 term, all dim) / (all terms)
                                                                    vector<SeparatedConvolutionInternal>
                                                                        (1 term, all dim) / (all terms)
 
 
                                                                    vector<ConvolutionData1D>
                                                                        (1 term, 1 dim) / (all dim)
 
    ConvolutionND and SeparatedConvolutionInternal both point to the same data in ConvolutionData1D.
    Why we need SeparatedConvolutionInternal in the first place I have no idea. ConvolutionND has the global
    factor, and SeparatedConvolutionInternal has a norm.
 
 
    */
 
    /// The largest distance a kernel must represent along the axes that are not summed over an
    /// infinite lattice: within the cell, plus N cells along an axis summed over N images.
    /// This is the accuracy range of the finite axes for GFit::truncate_mixed_expansion.
    template <std::size_t NDIM>
    inline double lattice_finite_reach(const std::array<LatticeRange, NDIM>& lattice_ranges,
                                       const Tensor<double>& cell_width) {
        double reach2 = 0.0, diag2 = 0.0;
        for (std::size_t d = 0; d != NDIM; ++d) {
            const double w = cell_width(long(d));
            diag2 += w * w;
            if (lattice_ranges[d].infinite()) continue;
            const double r = (lattice_ranges[d].get_range() + 1) * w;
            reach2 += r * r;
        }
        return reach2 > 0.0 ? std::sqrt(reach2) : std::sqrt(diag2);
    }
 
    template <typename Q, std::size_t NDIM>
    class SeparatedConvolution : public WorldObject< SeparatedConvolution<Q,NDIM> > {
    public:
 
        typedef Q opT;  ///< The apply function uses this to infer resultT=opT*inputT
 
        OperatorInfo info;
 
        bool doleaves;  ///< If should be applied to leaf coefficients ... false by default
 
      private:
        array_of_bools<NDIM>
            lattice_summed_;    ///< If lattice_summed_[d] is true, sum 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
        array_of_bools<NDIM> func_domain_is_periodic_{false};    ///< If domain_is_periodic_[d]==false and lattice_summed_[d]==false,
                                                            ///< ignore periodicity of BC when applying this to function
        std::array<KernelRange, NDIM> range;  ///< kernel range is along axis d is limited by range[d] if it's nonnull
 
      public:
        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)
        bool print_timings=false;
 
        typedef Key<NDIM> keyT;
        const static size_t opdim=NDIM;
        Timer timer_full;
        Timer timer_low_transf;
        Timer timer_low_accumulate;
        Timer timer_stats_accumulate;
 
    private:
 
 
        mutable std::vector< ConvolutionND<Q,NDIM> > ops;   ///< ConvolutionND keeps data for 1 term, all dimensions, 1 displacement
        const int k;
        const FunctionCommonData<Q,NDIM>& cdata;
        int rank;
        const std::vector<long> vk;
        const std::vector<long> v2k;
        const std::vector<Slice> s0;
 
        // SeparatedConvolutionData keeps data for all terms and all dimensions and 1 displacement
        mutable SimpleCache< SeparatedConvolutionData<Q,NDIM>, NDIM > data; ///< cache for all terms, dims and displacements
        mutable SimpleCache< SeparatedConvolutionData<Q,NDIM>, 2*NDIM > mod_data; ///< cache for all terms, dims and displacements
 
    public:
 
        bool& modified() {return modified_;}
        const bool& modified() const {return modified_;}
 
        int& particle() {return particle_;}
        const int& particle() const {return particle_;}
        SeparatedConvolution<Q,NDIM>& set_particle(const int p) {
            if (p!=1 and p!=2) throw std::runtime_error("particle must be 1 or 2");
            particle_=p;
            return *this;
        }
 
        bool& destructive() {return destructive_;}
        const bool& destructive() const {return destructive_;}
 
        const double& gamma() const {return info.mu;}
        const double& mu() const {return info.mu;}
        int get_rank() const { return rank; }
        int get_k() const { return k; }
        const std::vector<ConvolutionND<Q,NDIM>>& get_ops() const { return ops; }
        const std::array<KernelRange, NDIM>& get_range() const { return range; }
        bool range_restricted() const { return std::any_of(range.begin(), range.end(), [](const auto& v) { return v.finite(); }); }
 
    private:
 
        /// laziness for calling lists: which terms to apply
        struct ApplyTerms {
            ApplyTerms() : r_term(false), t_term(false) {}
            bool r_term;
            bool t_term;
            bool any_terms() const {return r_term or t_term;}
        };
 
        /// too lazy for extended calling lists
        struct Transformation {
            long r;             // Effective rank of transformation
            const Q* U;         // Ptr to matrix
            const Q* VT;
        };
 
        static inline std::pair<Tensor<double>,Tensor<double>>
        make_coeff_for_operator(World& world, OperatorInfo& info,
                                const std::array<LatticeRange, NDIM>& lattice_ranges) {
 
          const Tensor<double> &cell_width =
              FunctionDefaults<NDIM>::get_cell_width();
          double hi = cell_width.normf(); // Diagonal width of cell
          // Extend kernel range for lattice summation
          // N.B. if have periodic boundaries, extend range just in case will be using periodic domain
          bool lattice_summed_any = std::any_of(
              lattice_ranges.begin(), lattice_ranges.end(), [](const LatticeRange& b){ return static_cast<bool>(b); });
          bool infinite_summed_any = false;
          for (size_t i = 0; i < NDIM; i++) {
            if (lattice_ranges[i].infinite() && info.range[i].infinite()) {
              infinite_summed_any = true;
              break;
            }
          }
          // A Gaussian's lattice sum is a gauge constant only if the kernel is summed over the
          // whole lattice. An axis with an infinite lattice range but a finite kernel range is
          // summed over the images the range reaches, so for the truncation it counts as a finite
          // sum over that many images -- consistent with infinite_summed_any above.
          std::array<LatticeRange, NDIM> summed_ranges = lattice_ranges;
          for (size_t i = 0; i < NDIM; i++) {
            if (lattice_ranges[i].infinite() && info.range[i].finite())
              summed_ranges[i] = LatticeRange((info.range[i].iextent_x2() + 1) / 2);
          }
          const double hi_fin = lattice_finite_reach(summed_ranges, cell_width);
          if (lattice_summed_any || FunctionDefaults<NDIM>::get_bc().is_periodic_any()) {
            hi *= 100;
          }
 
          info.hi = hi;
          GFit<double, NDIM> fit(info);
 
          Tensor<double> coeff = fit.coeffs();
          Tensor<double> expnt = fit.exponents();
 
          if (info.truncate_lowexp_gaussians.value_or(infinite_summed_any)) {
            fit.truncate_mixed_expansion(coeff, expnt, summed_ranges, cell_width, info.lo, hi_fin, info.thresh);
            info.truncate_lowexp_gaussians = true;
          }
 
          return std::make_pair(coeff, expnt);
        }
 
//        /// return the right block of the upsampled operator (modified NS only)
//
//        /// unlike the operator matrices on the natural level the upsampled operator
//        /// matrices are not Toeplitz, so we need more information than just the displacement
//        ///.@param[in]  source  the source key
//        /// @param[in]  disp    the displacement
//        /// @param[in]  upop    the unfiltered operator matrix from scale n-1
//        /// @return     (k,k) patch of the upop(2k,2k) matrix
//        static Tensor<Q> operator_patch(const Translation& source, const Translation& disp, const Tensor<Q>& upop) {
//
//            // which of the 4 upsampled matrices do we need?
//            Translation sx=source%2;              // source offset
//            Translation tx=(source+disp)%2;       // target offset
//
//            Tensor<Q> rij(k,k);
//            // those two are equivalent:
///*
//            if (sx==0 and tx==0) copy_2d_patch(rij.ptr(),             k, upop.ptr(), 2*k, k, k);
//            if (sx==1 and tx==0) copy_2d_patch(rij.ptr() + k,         k, upop.ptr(), 2*k, k, k);
//            if (sx==0 and tx==1) copy_2d_patch(rij.ptr() + 2*k*k,     k, upop.ptr(), 2*k, k, k);
//            if (sx==1 and tx==1) copy_2d_patch(rij.ptr() + 2*k*k + k, k, upop.ptr(), 2*k, k, k);
//*/
//            Slice s0(0,k-1), s1(k,2*k-1);
//            if (sx==0 and tx==0) rij=Rm(s0,s0);
//            if (sx==1 and tx==0) rij=Rm(s1,s0);
//            if (sx==0 and tx==1) rij=Rm(s0,s1);
//            if (sx==1 and tx==1) rij=Rm(s1,s1);
//
//            return rij;
//        }
 
 
 
        /// accumulate into result
        template <typename T, typename R>
        void apply_transformation(long dimk,
                                  const Transformation trans[NDIM],
                                  const Tensor<T>& f,
                                  Tensor<R>& work1,
                                  Tensor<R>& work2,
                                  const Q mufac,
                                  Tensor<R>& result) const {
 
            //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
            long size = 1;
            for (std::size_t i=0; i<NDIM; ++i) size *= dimk;
            long dimi = size/dimk;
 
            R* MADNESS_RESTRICT w1=work1.ptr();
            R* MADNESS_RESTRICT w2=work2.ptr();
 
            mTxmq(dimi, trans[0].r, dimk, w1, f.ptr(), trans[0].U, dimk);
 
            size = trans[0].r * size / dimk;
            dimi = size/dimk;
            for (std::size_t d=1; d<NDIM; ++d) {
                mTxmq(dimi, trans[d].r, dimk, w2, w1, trans[d].U, dimk);
                size = trans[d].r * size / dimk;
                dimi = size/dimk;
                std::swap(w1,w2);
            }
 
            // If all blocks are full rank we can skip the transposes
            bool doit = false;
            for (std::size_t d=0; d<NDIM; ++d) doit = doit || trans[d].VT;
 
            if (doit) {
                for (std::size_t d=0; d<NDIM; ++d) {
                    if (trans[d].VT) {
                        dimi = size/trans[d].r;
                        mTxmq(dimi, dimk, trans[d].r, w2, w1, trans[d].VT);
                        size = dimk*size/trans[d].r;
                    }
                    else {
                        fast_transpose(dimk, dimi, w1, w2);
                    }
                    std::swap(w1,w2);
                }
            }
            // Assuming here that result is contiguous and aligned
            aligned_axpy(size, result.ptr(), w1, mufac);
        }
 
 
        /// accumulate into result
        template <typename T, typename R>
        void apply_transformation3(const Tensor<T> trans2[NDIM],
                                  const Tensor<T>& f,
                                  const Q mufac,
                                  Tensor<R>& result) const {
 
            //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
 
            Tensor<R> result2=general_transform(f,trans2);
            result2.scale(mufac);
            result+=result2;
 
        }
 
 
 
        /// don't accumulate, since we want to do this at apply()
        template <typename T, typename R>
        void apply_transformation2(Level n, long dimk,  double tol,
                                  const Tensor<T> trans2[NDIM],
                                  const GenTensor<T>& f,
                                  GenTensor<R>& work1,
                                  GenTensor<R>& work2,
                                  const Q mufac,
                                  GenTensor<R>& result) const {
 
            //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
 
#if 1
            result=general_transform(f,trans2);
            result.scale(mufac);
 
#else
 
            long size = 1;
            for (std::size_t i=0; i<NDIM; ++i) size *= dimk;
            long dimi = size/dimk;
 
            R* MADNESS_RESTRICT w1=work1.ptr();
            R* MADNESS_RESTRICT w2=work2.ptr();
 
            mTxmq(dimi, trans[0].r, dimk, w1, f.ptr(), trans[0].U, dimk);
            size = trans[0].r * size / dimk;
            dimi = size/dimk;
            for (std::size_t d=1; d<NDIM; ++d) {
                mTxmq(dimi, trans[d].r, dimk, w2, w1, trans[d].U, dimk);
                size = trans[d].r * size / dimk;
                dimi = size/dimk;
                std::swap(w1,w2);
            }
 
            // If all blocks are full rank we can skip the transposes
            bool doit = false;
            for (std::size_t d=0; d<NDIM; ++d) doit = doit || trans[d].VT;
 
            if (doit) {
                for (std::size_t d=0; d<NDIM; ++d) {
                    if (trans[d].VT) {
                        dimi = size/trans[d].r;
                        mTxmq(dimi, dimk, trans[d].r, w2, w1, trans[d].VT);
                        size = dimk*size/trans[d].r;
                    }
                    else {
                        fast_transpose(dimk, dimi, w1, w2);
                    }
                    std::swap(w1,w2);
                }
            }
            // Assuming here that result is contiguous and aligned
            aligned_axpy(size, result.ptr(), w1, mufac);
            //    long one = 1;
            //daxpy_(&size, &mufac, w1, &one, result.ptr(), &one);
#endif
        }
 
 
        /// Apply one of the separated terms, accumulating into the result
        template <typename T>
        void muopxv_fast(ApplyTerms at,
                         const ConvolutionData1D<Q>* const ops_1d[NDIM],
                         const Tensor<T>& f, const Tensor<T>& f0,
                         Tensor<TENSOR_RESULT_TYPE(T,Q)>& result,
                         Tensor<TENSOR_RESULT_TYPE(T,Q)>& result0,
                         const double tol,
                         const Q mufac,
                         Tensor<TENSOR_RESULT_TYPE(T,Q)>& work1,
                         Tensor<TENSOR_RESULT_TYPE(T,Q)>& work2) const {
 
            //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
            Transformation trans[NDIM];
            Tensor<T> trans2[NDIM];
 
            double Rnorm = 1.0;
            for (std::size_t d=0; d<NDIM; ++d) Rnorm *= ops_1d[d]->Rnorm;
 
            if (at.r_term and (Rnorm > 1.e-20)) {
 
                const auto tol_Rs = tol/(Rnorm*NDIM);  // Errors are relative within here
 
                // Determine rank of SVD to use or if to use the full matrix
                long twok = 2*k;
                if (modified()) twok=k;
 
                long break_even;
                if (NDIM==1) break_even = long(0.5*twok);
                else if (NDIM==2) break_even = long(0.6*twok);
                else if (NDIM==3) break_even=long(0.65*twok);
                else break_even=long(0.7*twok);
                bool rank_is_zero = false;
                for (std::size_t d=0; d<NDIM; ++d) {
                    long r;
                    for (r=0; r<twok; ++r) {
                        if (ops_1d[d]->Rs[r] < tol_Rs) break;
                    }
                    if (r >= break_even) {
                        trans[d].r = twok;
                        trans[d].U = ops_1d[d]->R.ptr();
                        trans[d].VT = 0;
                    }
                    else {
 
#ifdef USE_GENTENSOR
                        r = std::max(2L,r+(r&1L)); // (needed for 6D == when GENTENSOR is on) NOLONGER NEED TO FORCE OPERATOR RANK TO BE EVEN
#endif
                        if (r == 0) {
                            rank_is_zero = true;
                            break;
                        }
                        trans[d].r = r;
                        trans[d].U = ops_1d[d]->RU.ptr();
                        trans[d].VT = ops_1d[d]->RVT.ptr();
                    }
                    trans2[d]=ops_1d[d]->R;
                }
 
                if (!rank_is_zero)
                    apply_transformation(twok, trans, f, work1, work2, mufac, result);
 
                //            apply_transformation2(n, twok, tol, trans2, f, work1, work2, mufac, result);
//                apply_transformation3(trans2, f, mufac, result);
            }
 
            double Tnorm = 1.0;
            for (std::size_t d=0; d<NDIM; ++d) Tnorm *= ops_1d[d]->Tnorm;
 
            if (at.t_term and (Tnorm>0.0)) {
                const auto tol_Ts = tol/(Tnorm*NDIM);  // Errors are relative within here
 
                long break_even;
                if (NDIM==1) break_even = long(0.5*k);
                else if (NDIM==2) break_even = long(0.6*k);
                else if (NDIM==3) break_even=long(0.65*k);
                else break_even=long(0.7*k);
                bool rank_is_zero = false;
                for (std::size_t d=0; d<NDIM; ++d) {
                    long r;
                    for (r=0; r<k; ++r) {
                        if (ops_1d[d]->Ts[r] < tol_Ts) break;
                    }
                    if (r >= break_even) {
                        trans[d].r = k;
                        trans[d].U = ops_1d[d]->T.ptr();
                        trans[d].VT = 0;
                    }
                    else {
 
#ifdef USE_GENTENSOR
                        r = std::max(2L,r+(r&1L)); // (needed for 6D == GENTENSOR is USED) NOLONGER NEED TO FORCE OPERATOR RANK TO BE EVEN
#endif
                        if (r == 0) {
                            rank_is_zero = true;
                            break;
                        }
                        trans[d].r = r;
                        trans[d].U = ops_1d[d]->TU.ptr();
                        trans[d].VT = ops_1d[d]->TVT.ptr();
                    }
                    trans2[d]=ops_1d[d]->T;
                }
                if (!rank_is_zero)
                    apply_transformation(k, trans, f0, work1, work2, -mufac, result0);
//                apply_transformation2(n, k, tol, trans2, f0, work1, work2, -mufac, result0);
//                apply_transformation3(trans2, f0, -mufac, result0);
            }
        }
 
 
        /// Apply one of the separated terms, accumulating into the result
        template <typename T>
        void muopxv_fast2(Level n,
                         const ConvolutionData1D<Q>* const ops_1d[NDIM],
                         const GenTensor<T>& f, const GenTensor<T>& f0,
                         GenTensor<TENSOR_RESULT_TYPE(T,Q)>& result,
                         GenTensor<TENSOR_RESULT_TYPE(T,Q)>& result0,
                         double tol,
                         const Q mufac,
                         GenTensor<TENSOR_RESULT_TYPE(T,Q)>& work1,
                         GenTensor<TENSOR_RESULT_TYPE(T,Q)>& work2) const {
 
            PROFILE_MEMBER_FUNC(SeparatedConvolution);
//            Transformation trans[NDIM];
            Tensor<T> trans2[NDIM];
//            MADNESS_EXCEPTION("no muopxv_fast2",1);
 
            double Rnorm = 1.0;
            for (std::size_t d=0; d<NDIM; ++d) Rnorm *= ops_1d[d]->Rnorm;
            if (Rnorm == 0.0) return;
 
            if (Rnorm > 1.e-20) {
 
                tol = tol/(Rnorm*NDIM);  // Errors are relative within here
 
                // Determine rank of SVD to use or if to use the full matrix
                long twok = 2*k;
                if (modified()) twok=k;
//              long break_even;
//              if (NDIM==1) break_even = long(0.5*twok);
//              else if (NDIM==2) break_even = long(0.6*twok);
//              else if (NDIM==3) break_even=long(0.65*twok);
//              else break_even=long(0.7*twok);
                for (std::size_t d=0; d<NDIM; ++d) {
                    // long r;
                    // for (r=0; r<twok; ++r) {
                    //  if (ops_1d[d]->Rs[r] < tol) break;
                    // }
//                  if (r >= break_even) {
//                      trans[d].r = twok;
//                      trans[d].U = ops_1d[d]->R.ptr();
//                      trans[d].VT = 0;
//                  }
//                  else {
//                      //r += std::max(2L,r&1L); // NOLONGER NEED TO FORCE OPERATOR RANK TO BE EVEN
//                      trans[d].r = r; 
//                      trans[d].U = ops_1d[d]->RU.ptr();
//                      trans[d].VT = ops_1d[d]->RVT.ptr();
//                  }
                    trans2[d]=ops_1d[d]->R;
                }
                apply_transformation2(n, twok, tol, trans2, f, work1, work2, mufac, result);
            }
 
            double Tnorm = 1.0;
            for (std::size_t d=0; d<NDIM; ++d) Tnorm *= ops_1d[d]->Tnorm;
 
            if (n > 0 and (Tnorm>1.e-20)) {
//              long break_even;
//
//                if (NDIM==1) break_even = long(0.5*k);
//                else if (NDIM==2) break_even = long(0.6*k);
//                else if (NDIM==3) break_even=long(0.65*k);
//                else break_even=long(0.7*k);
                for (std::size_t d=0; d<NDIM; ++d) {
                    // long r;
                    // for (r=0; r<k; ++r) {
                    //     if (ops_1d[d]->Ts[r] < tol) break;
                    // }
//                    if (r >= break_even) {
//                        trans[d].r = k;
//                        trans[d].U = ops_1d[d]->T.ptr();
//                        trans[d].VT = 0;
//                    }
//                    else {
//                        //r += std::max(2L,r&1L); // NOLONGER NEED TO FORCE OPERATOR RANK TO BE EVEN
//                        trans[d].r = r; 
//                        trans[d].U = ops_1d[d]->TU.ptr();
//                        trans[d].VT = ops_1d[d]->TVT.ptr();
//                    }
                    trans2[d]=ops_1d[d]->T;
                }
                apply_transformation2(n, k, tol, trans2, f0, work1, work2, -mufac, result0);
            }
        }
 
 
        /// Computes the Frobenius norm of one of the separated terms ... WITHOUT FACTOR INCLUDED
        /// compute for 1 term, all dim, 1 disp, essentially for SeparatedConvolutionInternal
        double munorm2(Level n, const ConvolutionData1D<Q>* ops[]) const {
            if (modified()) return munorm2_modified(n,ops);
            return munorm2_ns(n,ops);
        }
 
        /// Computes the Frobenius norm of one of the separated terms for the NS form
        ///       ... WITHOUT FACTOR INCLUDED
        /// compute for 1 term, all dim, 1 disp, essentially for SeparatedConvolutionInternal
        double munorm2_ns(Level n, const ConvolutionData1D<Q>* ops[]) const {
            //PROFILE_MEMBER_FUNC(SeparatedConvolution);
            
            double prod=1.0, sum=0.0;
            for (std::size_t d=0; d<NDIM; ++d) {
                double a = ops[d]->NSnormf;
                double b = ops[d]->Tnormf;
                double aa = std::min(a,b);
                double bb = std::max(a,b);
                prod *= bb;
                if (bb > 0.0) sum +=(aa/bb);
            }
            if (n) prod *= sum;
 
            return prod;
        }
 
 
        /// Computes the operator norm of one of the separated terms of the modified NS form
        ///    ... WITHOUT FACTOR INCLUDED
        /// compute for 1 term, all dim, 1 disp, essentially for SeparatedConvolutionInternal
        double munorm2_modified(Level n, const ConvolutionData1D<Q>* ops_1d[]) const {
            PROFILE_MEMBER_FUNC(SeparatedConvolution);
 
            // follows Eq. (21) ff of Beylkin 2008 (Beylkin Appl. Comput. Harmon. Anal. 24, pp 354)
 
            // we have all combinations of difference, upsampled, F terms (d, u, f),
            // with the constraint that d is in each term exactly once. In the mixed terms (udf)
            // we just get all possible combinations, in the pure terms (dff, duu) we have
            // to multiply each term (dff, fdf, ffd) with (NDIM-1)!, to get the right number.
 
            double dff = 0.0;
            double duu = 0.0;
            double udf = 0.0;
 
            // loop over d shifting over the dimensions dxx, xdx, xxd,
            for (size_t d=0; d<NDIM; ++d) {
                double dff_tmp = ops_1d[d]->N_diff;
                double duu_tmp = ops_1d[d]->N_diff;
                double udf_tmp = ops_1d[d]->N_diff;
 
 
                for (size_t dd=0; dd<NDIM; ++dd) {
                    if (dd!=d) {
                        dff_tmp *= ops_1d[dd]->N_F;
                        duu_tmp *= ops_1d[dd]->N_up;
 
                        udf_tmp *= ops_1d[dd]->N_F;
                        for (size_t ddd=0; ddd<NDIM; ++ddd) {
                            if (ddd!=dd) udf += udf_tmp * ops_1d[ddd]->N_up;
                        }
                    }
                }
 
                dff+=dff_tmp;
                duu+=duu_tmp;
            }
 
            // finalize with the factorial
            double factorial=1.0;
            for (int i=1; i<static_cast<int>(NDIM)-1; ++i) factorial*=double(i);
            dff*=factorial;
            duu*=factorial;
 
            // Eq. (23) of Beylkin 2008, for one separated term WITHOUT the factor
            double norm=(dff + udf + duu) /(factorial * double(NDIM));
 
//            // double check
//            if (NDIM==3) {
//                Tensor<Q> R_full=outer(ops_1d[0]->R,outer(ops_1d[1]->R,ops_1d[2]->R));
//                Tensor<Q> T_full=outer(ops_1d[0]->T,outer(ops_1d[1]->T,ops_1d[2]->T));
//                double n2=(R_full-T_full).normf();
////                print("norm estimate, norm",norm, n2, norm<n2);
//                norm=n2;
//            }
 
            return norm;
            
        }
 
 
        /// get the transformation matrices for 1 term and all dimensions and one displacement
 
        /// use ConvolutionND, which uses ConvolutionData1D to collect the transformation matrices
        const SeparatedConvolutionInternal<Q,NDIM> getmuop(int mu, Level n, const Key<NDIM>& disp) const {
            //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
            SeparatedConvolutionInternal<Q,NDIM> op;
            for (std::size_t d=0; d<NDIM; ++d) {
                op.ops[d] = ops[mu].getop(d)->nonstandard(n, disp.translation()[d]);
            }
            op.norm = munorm2(n, op.ops)*std::abs(ops[mu].getfac());
 
//             double newnorm = munorm2(n, op.ops);
//             // This rescaling empirically based upon BSH separated expansion
//             // ... needs more testing.  OK also for TDSE.
//             // All is good except for some 000 blocks which are up to sqrt(k^d) off.
//             for (int d=0; d<NDIM; ++d)  {
//                 if (disp[d] == 0) newnorm *= 0.5;
//                 else if (std::abs(disp[d]) == 1) newnorm *= 0.8;
//             }
//            double oldnorm = munorm(n, op.ops);
//             if (oldnorm > 1e-13 && (newnorm < 0.5*oldnorm || newnorm > 2.0*oldnorm) )
//                 print("munorm", n, disp, mu, newnorm, oldnorm, newnorm/oldnorm);
 
            return op;
        }
 
 
 
        /// get the transformation matrices for 1 term and all dimensions and one displacement
 
        /// use ConvolutionND, which uses ConvolutionData1D to collect the transformation matrices
        const SeparatedConvolutionInternal<Q,NDIM>
        getmuop_modified(int mu, Level n, const Key<NDIM>& disp, const Key<NDIM>& source) const {
            //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
 
 
            // SeparatedConvolutionInternal keeps data for 1 term and all dimensions
            SeparatedConvolutionInternal<Q,NDIM> op;
 
            // in the modified NS form we need not only the displacement, but also the source Translation
            // for correctly constructing the operator, b/c the operator is not Toeplitz
 
            // op.ops is of type ConvolutionData1D (1 term, 1 dim, 1 disp)
            // ops[mu] is of type ConvolutionND (1 term, all dim, 1 disp)
            for (std::size_t d=0; d<NDIM; ++d) {
                Translation sx=source.translation()[d];                          // source translation
                Translation tx=source.translation()[d]+disp.translation()[d];    // target translation
 
                Key<2> op_key(n,Vector<Translation,2>{sx,tx});
                op.ops[d] = ops[mu].getop(d)->mod_nonstandard(op_key);
            }
 
            // works for both modified and not modified NS form
            op.norm = munorm2(n, op.ops)*std::abs(ops[mu].getfac());
//            op.norm=1.0;
            return op;
        }
 
        /// get the data for all terms and all dimensions for one displacement
        const SeparatedConvolutionData<Q,NDIM>* getop(Level n, const Key<NDIM>& d, const Key<NDIM>& source) const {
 
            // in the NS form the operator depends only on the displacement
            if (not modified()) return getop_ns(n,d);
            return getop_modified(n, d, source);
        }
 
 
        /// get the data for all terms and all dimensions for one displacement
 
        /// uses SeparatedConvolutionInternal (ConvolutionND, ConvolutionData1D) to construct
        /// the transformation matrices.
        /// @param[in]  d   displacement
        /// @return pointer to cached operator
        const SeparatedConvolutionData<Q,NDIM>* getop_ns(Level n, const Key<NDIM>& d) const {
            //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
            const SeparatedConvolutionData<Q,NDIM>* p = data.getptr(n,d);
            if (p) return p;
 
            // get the data for each term
            SeparatedConvolutionData<Q,NDIM> op(rank);
            for (int mu=0; mu<rank; ++mu) {
                // op.muops is of type SeparatedConvolutionInternal (1 term, all dim, 1 disp)
                // getmuop uses ConvolutionND
                op.muops[mu] = getmuop(mu, n, d);
            }
 
            double norm = 0.0;
            for (int mu=0; mu<rank; ++mu) {
                const double munorm = op.muops[mu].norm;
                norm += munorm*munorm;
            }
        //print("getop", n, d, norm);
            op.norm = sqrt(norm);
            return data.set(n, d, op);
        }
 
 
 
        /// get the data for all terms and all dimensions for one displacement (modified NS form)
 
        /// remember that the operator in the modified NS form is not Toeplitz, so we need
        /// information about the displacement and the source key
        /// @param[in]  n       level (=scale) (actually redundant, since included in source)
        /// @param[in]  disp    displacement key
        /// @param[in]  source  source key
        /// @return pointer to cached operator
        const SeparatedConvolutionData<Q,NDIM>* getop_modified(Level n, const Key<NDIM>& disp, const Key<NDIM>& source) const {
            //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
 
            // in the modified NS form the upsampled part of the operator depends on the modulus of the source
            Vector<Translation,NDIM> t=source.translation();
            for (size_t i=0; i<NDIM; ++i) t[i]=t[i]%2;
            Key<2*NDIM> key=disp.merge_with(Key<NDIM>(source.level(),t));
 
            const SeparatedConvolutionData<Q,NDIM>* p = mod_data.getptr(n,key);
            if (p) return p;
 
            // get the data for each term
            // op.muops is of type SeparatedConvolutionInternal (1 term, all dim, 1 disp)
            // getmuop uses ConvolutionND
            SeparatedConvolutionData<Q,NDIM> op(rank);
            for (int mu=0; mu<rank; ++mu) op.muops[mu] = getmuop_modified(mu, n, disp, source);
 
            double norm = 0.0;
            for (int mu=0; mu<rank; ++mu) {
                const double munorm = op.muops[mu].norm;
                norm += munorm*munorm;
            }
 
            op.norm = sqrt(norm);
            return mod_data.set(n, key, op);
        }
 
 
        void check_cubic() {
            // !!! NB ... cell volume obtained from global defaults
            const Tensor<double>& cell_width = FunctionDefaults<NDIM>::get_cell_width();
            // Check that the cell is cubic since currently is assumed
            for (std::size_t d=1; d<NDIM; ++d) {
                MADNESS_CHECK(fabs(cell_width(d)-cell_width(0L)) < 1e-14*cell_width(0L));
            }
        }
 
 
        /// upsample some of the dimensions of coeff to its child indicated by key
 
        /// @param[in]  coeff   the coeffs of dim 2*NDIM that will be upsampled
        /// @param[in]  key     the key indicating the child -- only some dimensions will be "reproductive"
        /// @param[in]  particle    if 0: upsample dimensions 0-2
        ///                         if 1: upsample dimensions 3-5
        /// @return     a partially upsampled coefficient tensor
        template<typename T, size_t FDIM>
        GenTensor<T> partial_upsample(const Key<FDIM>& key, const GenTensor<T>& coeff, const int particle) const {
 
            if (coeff.rank()==0) return GenTensor<T>();
            MADNESS_ASSERT(coeff.dim(0)==k);
            if (NDIM==coeff.ndim()) {
                MADNESS_ASSERT(particle==1);    // other particle, leave this particle unchanged
                return coeff;
            }
 
            MADNESS_ASSERT(coeff.ndim()==FDIM);
            MADNESS_ASSERT(particle==0 or (2*NDIM==FDIM));
 
            // the twoscale coefficients: for upsampling use h0/h1; see Alpert Eq (3.35a/b)
            // handle the spectator dimensions with the identity matrix
            const Tensor<T> h[2] = {cdata.h0, cdata.h1};
            Tensor<T> identity(k,k);
            for (int i=0; i<k; ++i) identity(i,i)=1.0;
            Tensor<T> matrices[2*NDIM];
 
            // get the appropriate twoscale coefficients for each dimension
            if (particle==0) {
                for (size_t ii=0; ii<NDIM; ++ii) matrices[ii]=h[key.translation()[ii]%2];
                for (size_t ii=0; ii<NDIM; ++ii) matrices[ii+NDIM]=identity;
            } else if (particle==1) {
                for (size_t ii=0; ii<NDIM; ++ii) matrices[ii]=identity;
                for (size_t ii=0; ii<NDIM; ++ii) matrices[ii+NDIM]=h[key.translation()[ii+NDIM]%2];
            } else {
                MADNESS_EXCEPTION("unknown particle",1);
            }
 
            // transform and accumulate on the result
            const GenTensor<T> result=general_transform(coeff,matrices);
            return result;
        }
 
 
        /// upsample the sum coefficients of level 1 to sum coeffs on level n+1
 
        /// specialization of the unfilter method, will transform only the sum coefficients
        /// @param[in]  key     key of level n+1
        /// @param[in]  coeff   sum coefficients of level n (does NOT belong to key!!)
        /// @return     sum     coefficients on level n+1
        template<typename T, size_t FDIM>
        GenTensor<T> upsample(const Key<FDIM>& key, const GenTensor<T>& coeff) const {
 
            // the twoscale coefficients: for upsampling use h0/h1; see Alpert Eq (3.35a/b)
            // note there are no difference coefficients; if you want that use unfilter
            const Tensor<T> h[2] = {cdata.h0, cdata.h1};
            Tensor<T> matrices[FDIM];
 
            // get the appropriate twoscale coefficients for each dimension
            for (size_t ii=0; ii<FDIM; ++ii) matrices[ii]=h[key.translation()[ii]%2];
 
            // transform and accumulate on the result
            const GenTensor<T> result=general_transform(coeff,matrices);
            return result;
        }
 
        /// initializes range using range of ops[0]
        /// @pre `ops[i].range == ops[0].range`
        void init_range() {
          if (!ops.empty()) {
            for (int d = 0; d != NDIM; ++d) {
              for(const auto & op: ops) {
                MADNESS_ASSERT(op.getop(d)->range == ops[0].getop(d)->range);
              }
              range[d] = ops[0].getop(d)->range;
            }
          }
        }
 
        /// initializes lattice_sum using `ops[0].lattice_summed()`
        /// @pre `ops[i].lattice_summed() == ops[0].lattice_summed()`
        void init_lattice_summed() {
          if (!ops.empty()) {
            for (int d = 0; d != NDIM; ++d) {
              for (const auto &op : ops) {
                MADNESS_ASSERT(op.lattice_summed() ==
                               ops[0].lattice_summed());
              }
              lattice_summed_ = ops[0].lattice_summed();
            }
          }
        }
 
    public:
 
        // For separated convolutions with same operator in each direction (isotropic)
        SeparatedConvolution(World& world,
                             const std::vector< std::shared_ptr< Convolution1D<Q> > >& argops,
                             long k = FunctionDefaults<NDIM>::get_k(),
                             bool doleaves = false)
                : WorldObject< SeparatedConvolution<Q,NDIM> >(world)
                , info()
                , doleaves(doleaves)
                , lattice_summed_(false)  // this will be overridden by init_lattice_summed below
                , modified_(false)
                , particle_(1)
                , destructive_(false)
                , k(k)
                , cdata(FunctionCommonData<Q,NDIM>::get(k))
                , rank(argops.size())
                , vk(NDIM,k)
                , v2k(NDIM,2*k)
                , s0(std::max<std::size_t>(2,NDIM),Slice(0,k-1))
        {
 
            for (unsigned int mu=0; mu < argops.size(); ++mu) {
              this->ops.push_back(ConvolutionND<Q,NDIM>(argops[mu]));
            }
            init_range();
            init_lattice_summed();
 
            this->process_pending();
        }
 
        // For general convolutions
        SeparatedConvolution(World& world,
                             const std::vector< ConvolutionND<Q,NDIM> >& argops,
                             long k = FunctionDefaults<NDIM>::get_k(),
                             bool doleaves = false)
                : WorldObject< SeparatedConvolution<Q,NDIM> >(world)
                , info()
                , doleaves(doleaves)
                , lattice_summed_(false)  // this will be overridden by init_lattice_summed below
                , modified_(false)
                , particle_(1)
                , destructive_(false)
                , ops(argops)
                , k(k)
                , cdata(FunctionCommonData<Q,NDIM>::get(k))
                , rank(argops.size())
                , vk(NDIM,k)
                , v2k(NDIM,2*k)
                , s0(std::max<std::size_t>(2,NDIM),Slice(0,k-1))
        {
            init_range();
            init_lattice_summed();
            this->process_pending();
        }
 
        /// Constructor for Gaussian Convolutions
        /// WARNING! bloch_k should only ever be nonzero if `Q` is a complex type.
        SeparatedConvolution(World& world, const OperatorInfo info1,
                             const std::array<LatticeRange, NDIM>& lattice_ranges = FunctionDefaults<NDIM>::get_bc().lattice_range(),
                             int k=FunctionDefaults<NDIM>::get_k(),
                             bool doleaves = false,
                             const Vector<double, NDIM>& bloch_k = Vector<double, NDIM>(0.0))
               : SeparatedConvolution(world,Tensor<double>(0l),Tensor<double>(0l),info1.lo,info1.thresh,lattice_ranges,k,doleaves,info1.mu) {
            info.type=info1.type;
            info.truncate_lowexp_gaussians = info1.truncate_lowexp_gaussians;
            info.range = info1.range;
            auto [coeff, expnt] = make_coeff_for_operator(world, info, lattice_ranges);
            rank=coeff.dim(0);
            range = info.template range_as_array<NDIM>();
            ops.resize(rank);
            initialize(coeff,expnt,lattice_ranges,range,bloch_k);
            init_lattice_summed();
        }
 
        /// Constructor for Gaussian Convolutions (mostly for backward compatability)
        SeparatedConvolution(World& world,
                             const Tensor<Q>& coeff, const Tensor<double>& expnt,
                             double lo, double thresh,
                             const std::array<LatticeRange, NDIM>& lattice_ranges = FunctionDefaults<NDIM>::get_bc().lattice_range(),
                             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)
                ,
              lattice_summed_(false) // will be updated by init_lattice_summed
                , ops(coeff.dim(0))
                , k(k)
                , cdata(FunctionCommonData<Q,NDIM>::get(k))
                , rank(coeff.dim(0))
                , vk(NDIM,k)
                , v2k(NDIM,2*k)
                , s0(std::max<std::size_t>(2,NDIM),Slice(0,k-1)) {
            initialize(coeff,expnt,lattice_ranges);
            init_range();
            init_lattice_summed();
        }
 
        void initialize(const Tensor<Q>& coeff, const Tensor<double>& expnt, std::array<LatticeRange, NDIM> lattice_range, std::array<KernelRange, NDIM> range = {}, const Vector<double, NDIM>& bloch_k = Vector<double, NDIM>(0.0)) {
            const Tensor<double>& width = FunctionDefaults<NDIM>::get_cell_width();
            const double pi = constants::pi;
 
            for (int mu=0; mu<rank; ++mu) {
                Q c = std::pow(sqrt(expnt(mu)/pi),static_cast<int>(NDIM)); // Normalization coeff
 
                // We cache the normalized operator so the factor is the value we must multiply
                // by to recover the coeff we want.
                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,
                                       lattice_range[d], bloch_k[d], range[d]));
                }
            }
        }
 
        virtual ~SeparatedConvolution() { }
 
        void print_timer() const {
            if (this->get_world().rank()==0) {
                timer_full.print("op full tensor       ");
                timer_low_transf.print("op low rank transform");
                timer_low_accumulate.print("op low rank addition ");
            }
        }
 
        void reset_timer() const {
            if (this->get_world().rank()==0) {
                timer_full.reset();
                timer_low_transf.reset();
                timer_low_accumulate.reset();
            }
        }
 
        const std::vector< Key<NDIM> >& get_disp(Level n) const {
            return Displacements<NDIM>().get_disp(n, lattice_summed());
        }
 
        /// @return flag for each axis indicating whether lattice summation is performed in that direction
        const array_of_bools<NDIM>& lattice_summed() const { return lattice_summed_; }
        /// @return flag for each axis indicating whether functions that this op acts on are periodic on the box in that direction (false by default)
        ///         it is normally preferred to lattice sum, obtaining an equivalent problem with a lattice-summed operator on a non-periodic function
        const array_of_bools<NDIM>& func_domain_is_periodic() const { return func_domain_is_periodic_; }
        /// changes domain periodicity
        /// \param domain_is_periodic
        void set_domain_periodicity(const array_of_bools<NDIM>& domain_is_periodic) { func_domain_is_periodic_ = domain_is_periodic;}
 
        /// 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 {
            // SeparatedConvolutionData keeps data for all terms and all dimensions and 1 displacement
//            return 1.0;
            return getop(n, d, source_key)->norm;
        }
 
        /// return that part of a hi-dim key that serves as the base for displacements of this operator
 
        /// if the function and the operator have the same dimension return key
        /// if the function has a higher dimension than the operator (e.g. in the exchange operator)
        /// return only that part of key that corresponds to the particle this operator works on
        /// @param[in]  key hi-dim key
        /// @return     a lo-dim part of key; typically first or second half
        template<size_t FDIM>
        typename std::enable_if<FDIM!=NDIM, Key<NDIM> >::type
        get_source_key(const Key<FDIM> key) const {
            Key<NDIM> source;
            Key<FDIM-NDIM> dummykey;
            if (particle()==1) key.break_apart(source,dummykey);
            if (particle()==2) key.break_apart(dummykey,source);
            return source;
        }
 
        /// return that part of a hi-dim key that serves as the base for displacements of this operator
 
        /// if the function and the operator have the same dimension return key
        /// if the function has a higher dimension than the operator (e.g. in the exchange operator)
        /// return only that part of key that corresponds to the particle this operator works on
        /// @param[in]  key hi-dim key
        /// @return     a lo-dim part of key; typically first or second half
        template<size_t FDIM>
        typename std::enable_if<FDIM==NDIM, Key<NDIM> >::type
        get_source_key(const Key<FDIM> key) const {
            return key;
        }
 
        /// apply this operator on a function f
 
        /// the operator does not need to have the same dimension as the function, e,g,
        /// the Poisson kernel for the exchange operator acts only on 1 electron of a
        /// given (pair) function.
        /// @param[in]  f   a function of same or different dimension as this operator
        /// @return     the result function of the same dimensionality as the input function f
        template <typename T, size_t FDIM>
        Function<TENSOR_RESULT_TYPE(T,Q),FDIM> operator()(const Function<T,FDIM>& f) const {
            return madness::apply(*this, f);
        }
 
        /// apply this on a vector of functions
        template <typename T, size_t FDIM>
        std::vector<Function<TENSOR_RESULT_TYPE(T,Q),FDIM>> operator()(const std::vector<Function<T,FDIM>>& f) const {
            return madness::apply(*this, f);
        }
 
        /// apply this operator on a separable function f(1,2) = f(1) f(2)
 
        /// @param[in]  f1   a function of dim LDIM
        /// @param[in]  f2   a function of dim LDIM
        /// @return     the result function of dim NDIM=2*LDIM: g(1,2) = G(1,1',2,2') f(1',2')
        template <typename T, size_t LDIM>
        Function<TENSOR_RESULT_TYPE(T,Q),LDIM+LDIM>
        operator()(const Function<T,LDIM>& f1, const Function<Q,LDIM>& f2) const {
            return madness::apply(*this, std::vector<Function<Q,LDIM>>({f1}),
                                  std::vector<Function<Q,LDIM>>({f2}));
        }
 
        /// apply this operator on a sum of separable functions f(1,2) = \sum_i f_i(1) f_i(2)
 
        /// @param[in]  f1   a function of dim LDIM
        /// @param[in]  f2   a function of dim LDIM
        /// @return     the result function of dim NDIM=2*LDIM: g(1,2) = G(1,1',2,2') f(1',2')
        template <typename T, size_t LDIM>
        Function<TENSOR_RESULT_TYPE(T,Q),LDIM+LDIM>
        operator()(const std::vector<Function<T,LDIM>>& f1, const std::vector<Function<Q,LDIM>>& f2) const {
            return madness::apply(*this, f1, f2);
        }
 
        /// apply this onto another suitable argument, returning the same type
 
        /// argT must implement argT::apply(const SeparatedConvolution& op, const argT& arg)
        template<typename argT>
        argT operator()(const argT& argument) const {
            return madness::apply(*this,argument);
        }
 
 
        /// apply this operator on coefficients in full rank form
 
        /// @param[in]  coeff   source coeffs in full rank
        /// @param[in]  source  the source key
        /// @param[in]  shift   the displacement, where the source coeffs come from
        /// @param[in]  tol     thresh/#neigh*cnorm
        /// @return     a tensor of full rank with the result op(coeff)
        template <typename T>
        Tensor<TENSOR_RESULT_TYPE(T,Q)> apply(const Key<NDIM>& source,
                                              const Key<NDIM>& shift,
                                              const Tensor<T>& coeff,
                                              double tol) const {
            //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
            MADNESS_ASSERT(coeff.ndim()==NDIM);
 
            double cpu0=cpu_time();
 
            typedef TENSOR_RESULT_TYPE(T,Q) resultT;
            const Tensor<T>* input = &coeff;
            Tensor<T> dummy;
 
            if (not modified()) {
                if (coeff.dim(0) == k) {
                    // This processes leaf nodes with only scaling
                    // coefficients ... FuncImpl::apply by default does not
                    // apply the operator to these since for smoothing operators
                    // it is not necessary.  It is necessary for operators such
                    // as differentiation and time evolution and will also occur
                    // if the application of the operator widens the tree.
                    dummy = Tensor<T>(v2k);
                    dummy(s0) = coeff;
                    input = &dummy;
                }
                else {
                    MADNESS_ASSERT(coeff.dim(0)==2*k);
                }
            }
 
            tol = 0.01*tol/rank; // Error is per separated term
            ApplyTerms at;
            at.r_term=true;
            at.t_term=(source.level()>0);
 
            /// SeparatedConvolutionData keeps data for all terms and all dimensions and 1 displacement
            const SeparatedConvolutionData<Q,NDIM>* op = getop(source.level(), shift, source);
 
            //print("sepop",source,shift,op->norm,tol);
 
            Tensor<resultT> r(v2k), r0(vk);
            Tensor<resultT> work1(v2k,false), work2(v2k,false);
 
            if (modified()) {
                   r=Tensor<resultT>(vk);
                   work1=Tensor<resultT>(vk,false);
                   work2=Tensor<resultT>(vk,false);
            }
 
            const Tensor<T> f0 = copy(coeff(s0));
            for (int mu=0; mu<rank; ++mu) {
                // SeparatedConvolutionInternal keeps data for 1 term and all dimensions and 1 displacement
                const SeparatedConvolutionInternal<Q,NDIM>& muop =  op->muops[mu];
                if (muop.norm > tol) {
                    // ops is of ConvolutionND, returns data for 1 term and all dimensions
                    Q fac = ops[mu].getfac();
                    muopxv_fast(at, muop.ops, *input, f0, r, r0, tol/std::abs(fac), fac,
                                work1, work2);
                }
            }
 
            r(s0).gaxpy(1.0,r0,1.0);
            double cpu1=cpu_time();
            timer_full.accumulate(cpu1-cpu0);
 
            return r;
        }
 
 
        /// apply this operator on only 1 particle of the coefficients in low rank form
 
        /// note the unfortunate mess with NDIM: here NDIM is the operator dimension, and FDIM is the
        /// function's dimension, whereas in the function we have OPDIM for the operator and NDIM for
        /// the function
        /// @tparam T   the dimension of the function this operator is applied on. \todo MGR: Make sure info on T is correct. Was previously labeled FDIM.
        /// @param[in]  coeff   source coeffs in SVD (=optimal!) form, in high dimensionality (FDIM)
        /// @param[in]  source  the source key in low dimensionality (NDIM)
        /// @param[in]  shift   the displacement in low dimensionality (NDIM)
        /// @param[in]  tol     thresh/(#neigh*cnorm)
        /// @param[in]  tol2    thresh/#neigh
        /// @return     coeff result
        template<typename T>
        GenTensor<TENSOR_RESULT_TYPE(T,Q)> apply2_lowdim(const Key<NDIM>& source,
                const Key<NDIM>& shift, const GenTensor<T>& coeff, double tol, double tol2) const {
 
            typedef TENSOR_RESULT_TYPE(T,Q) resultT;
 
            // prepare access to the singular vectors
            const SVDTensor<T>& svdcoeff=coeff.get_svdtensor();
//            std::vector<Slice> s(coeff.config().dim_per_vector()+1,_);
            std::vector<Slice> s(svdcoeff.dim_per_vector(particle()-1)+1,_);
            // can't use predefined slices and vectors -- they have the wrong dimension
            const std::vector<Slice> s00(coeff.ndim(),Slice(0,k-1));
 
            // some checks
            MADNESS_ASSERT(coeff.is_svd_tensor());           // for now
            MADNESS_ASSERT(not modified());
            MADNESS_ASSERT(not doleaves);
            MADNESS_ASSERT(coeff.dim(0)==2*k);
            MADNESS_ASSERT(2*NDIM==coeff.ndim());
 
            double cpu0=cpu_time();
            const SeparatedConvolutionData<Q,NDIM>* op = getop(source.level(), shift, source);
 
            // some workspace
            Tensor<resultT> work1(v2k,false), work2(v2k,false);
 
            // sliced input and final result
            const GenTensor<T> f0 = copy(coeff(s00));
            GenTensor<resultT> final=copy(coeff);
            GenTensor<resultT> final0=copy(f0);
 
            tol = tol/rank*0.01; // Error is per separated term
            tol2= tol2/rank;
 
            // the operator norm is missing the identity working on the other particle
            // use as (muop.norm*exchange_norm < tol)
            // for some reason the screening is not working at all..
//            double exchange_norm=std::pow(2.0*k,1.5);
 
            for (int r=0; r<coeff.rank(); ++r) {
 
                // get the appropriate singular vector (left or right depends on particle)
                // and apply the full tensor muopxv_fast on it, term by term
                s[0]=Slice(r,r);
                const Tensor<T> chunk=svdcoeff.ref_vector(particle()-1)(s).reshape(v2k);
                const Tensor<T> chunk0=f0.get_svdtensor().ref_vector(particle()-1)(s).reshape(vk);
//                const double weight=std::abs(coeff.config().weights(r));
 
                // accumulate all terms of the operator for a specific term of the function
                Tensor<resultT> result(v2k), result0(vk);
 
                ApplyTerms at;
                at.r_term=true;
                at.t_term=source.level()>0;
 
                // this loop will return on result and result0 the terms [(P+Q) G (P+Q)]_1,
                // and [P G P]_1, respectively
                for (int mu=0; mu<rank; ++mu) {
                    const SeparatedConvolutionInternal<Q,NDIM>& muop =  op->muops[mu];
                    Q fac = ops[mu].getfac();
                    muopxv_fast(at, muop.ops, chunk, chunk0, result, result0,
                            tol/std::abs(fac), fac, work1, work2);
                }
 
                // reinsert the transformed terms into result, leaving the other particle unchanged
                MADNESS_ASSERT(final.get_svdtensor().has_structure());
                final.get_svdtensor().ref_vector(particle()-1)(s)=result;
 
                if (source.level()>0) {
                    final0.get_svdtensor().ref_vector(particle()-1)(s)=result0;
                } else {
                    final0.get_svdtensor().ref_vector(0)(s)=0.0;
                    final0.get_svdtensor().ref_vector(1)(s)=0.0;
                }
 
            }
            double cpu1=cpu_time();
            timer_low_transf.accumulate(cpu1-cpu0);
 
            double cpu00=cpu_time();
 
            final.reduce_rank(tol2*0.5);
            final0.reduce_rank(tol2*0.5);
            final(s00)+=final0;
            final.reduce_rank(tol2);
 
            double cpu11=cpu_time();
            timer_low_accumulate.accumulate(cpu11-cpu00);
            return final;
        }
 
        /// apply this operator on coefficients in low rank form
 
        /// @param[in]  coeff   source coeffs in SVD (=optimal!) form
        /// @param[in]  tol     thresh/#neigh*cnorm
        /// @param[in]  tol2    thresh/#neigh
        template <typename T>
        GenTensor<TENSOR_RESULT_TYPE(T,Q)> apply2(const Key<NDIM>& source,
                                              const Key<NDIM>& shift,
                                              const GenTensor<T>& coeff,
                                              double tol, double tol2) const {
            PROFILE_MEMBER_FUNC(SeparatedConvolution);
            typedef TENSOR_RESULT_TYPE(T,Q) resultT;
 
            MADNESS_ASSERT(coeff.ndim()==NDIM);
            MADNESS_ASSERT(coeff.is_svd_tensor());  // we use the rank below
//            MADNESS_EXCEPTION("no apply2",1);
            const TensorType tt=TT_2D;
 
            const GenTensor<T>* input = &coeff;
            GenTensor<T> dummy;
 
            if (not modified()) {
                if (coeff.dim(0) == k) {
                    // This processes leaf nodes with only scaling
                    // coefficients ... FuncImpl::apply by default does not
                    // apply the operator to these since for smoothing operators
                    // it is not necessary.  It is necessary for operators such
                    // as differentiation and time evolution and will also occur
                    // if the application of the operator widens the tree.
                    dummy = GenTensor<T>(v2k,TT_2D);
                    dummy(s0) += coeff;
                    input = &dummy;
                }
                else {
                    MADNESS_ASSERT(coeff.dim(0)==2*k);
                }
            }
 
            tol = tol/rank; // Error is per separated term
            tol2= tol2/rank;
 
            const SeparatedConvolutionData<Q,NDIM>* op = getop(source.level(), shift, source);
 
            GenTensor<resultT> r, r0, result, result0;
            GenTensor<resultT> work1(v2k,tt), work2(v2k,tt);
 
            if (modified()) {
                r=GenTensor<resultT>(vk,tt);
                work1=GenTensor<resultT>(vk,tt);
                work2=GenTensor<resultT>(vk,tt);
            }
 
            // collect the results of the individual operator terms
            std::list<GenTensor<T> > r_list;
            std::list<GenTensor<T> > r0_list;
 
//            const GenTensor<T> f0 = copy(coeff(s0));
            const GenTensor<T> f0 = copy((*input)(s0));
            for (int mu=0; mu<rank; ++mu) {
                const SeparatedConvolutionInternal<Q,NDIM>& muop =  op->muops[mu];
                //print("muop",source, shift, mu, muop.norm);
 
                // delta(g)  <  delta(T) * || f ||
                if (muop.norm > tol) {
 
                    // get maximum rank of coeff to contribute:
                    //  delta(g)  <  eps  <  || T || * delta(f)
                    //  delta(coeff) * || T || < tol2
                    const int r_max=SRConf<T>::max_sigma(tol2/muop.norm,coeff.rank(),coeff.get_svdtensor().weights_);
                    //                  print("r_max",coeff.config().weights(r_max));
 
                    // note that max_sigma is inclusive!
                    if (r_max>=0) {
                        const GenTensor<resultT> chunk=SVDTensor<resultT>(input->get_svdtensor().get_configs(0,r_max));
                        const GenTensor<resultT> chunk0=SVDTensor<resultT>(f0.get_svdtensor().get_configs(0,r_max));
 
                        double cpu0=cpu_time();
 
                        Q fac = ops[mu].getfac();
                        muopxv_fast2(source.level(), muop.ops, chunk, chunk0, r, r0,
                                tol/std::abs(fac), fac, work1, work2);
                        double cpu1=cpu_time();
                        timer_low_transf.accumulate(cpu1-cpu0);
 
                        r_list.push_back(r);
                        r0_list.push_back(r0);
                    }
                }
            }
 
            // finally accumulate all the resultant terms into one tensor
            double cpu0=cpu_time();
 
            result0=reduce(r0_list,tol2*rank);
            if (r_list.size()>0) r_list.front()(s0)+=result0;
            result=reduce(r_list,tol2*rank);
//            result.reduce_rank(tol2*rank);
 
            double cpu1=cpu_time();
            timer_low_accumulate.accumulate(cpu1-cpu0);
            timer_stats_accumulate.accumulate(result.rank());
            return result;
        }
 
        /// estimate the ratio of cost of full rank versus low rank
 
        /// @param[in]  source  source key
        /// @param[in]  shift   displacement
        /// @param[in]  tol     thresh/#neigh/cnorm
        /// @param[in]  tol2    thresh/#neigh
        /// @return cost_ratio  r=-1:   no terms left
        ///                     0<r<1:  better to do full rank
        ///                     1<r:    better to do low rank
        template<typename T>
        double estimate_costs(const Key<NDIM>& source,
                const Key<NDIM>& shift,
                const GenTensor<T>& coeff,
                double tol, double tol2) const {
 
            if (coeff.is_full_tensor()) return 0.5;
            if (2*NDIM==coeff.ndim()) return 1.5;
            MADNESS_ASSERT(NDIM==coeff.ndim());
            MADNESS_ASSERT(coeff.is_svd_tensor());
 
            const SeparatedConvolutionData<Q,NDIM>* op = getop(source.level(), shift, source);
 
            tol = tol/rank; // Error is per separated term
            tol2= tol2/rank;
 
            const double full_operator_cost=pow(coeff.dim(0),NDIM+1);
            const double low_operator_cost=pow(coeff.dim(0),NDIM/2+1);
            const double low_reduction_cost=pow(coeff.dim(0),NDIM/2);
 
            double full_cost=0.0;
            double low_cost=0.0;
 
            long initial_rank=0;
            long final_rank=sqrt(coeff.size())*0.05;    // size=ncol*nrow; final rank is 5% of max rank
 
            for (int mu=0; mu<rank; ++mu) {
                const SeparatedConvolutionInternal<Q,NDIM>& muop =  op->muops[mu];
 
                // delta(g)  <  delta(T) * || f ||
                if (muop.norm > tol) {
                    // note that max_sigma is inclusive: it returns a slice w(Slice(0,i))
                    long nterms=SRConf<T>::max_sigma(tol2/muop.norm,coeff.rank(),coeff.get_svdtensor().weights_)+1;
 
                    // take only the first overlap computation of rank reduction into account
//                    low_cost+=nterms*low_operator_cost + 2.0*nterms*nterms*low_reduction_cost;
                    initial_rank+=nterms;
 
                    full_cost+=full_operator_cost;
                }
            }
            low_cost=initial_rank*low_operator_cost + initial_rank*final_rank*low_reduction_cost;
 
            // include random empirical factor of 2
            double ratio=-1.0;
            if (low_cost>0.0) ratio=full_cost/low_cost;
//            print("nterms, full, low, full/low", full_cost, low_cost,shift.distsq(), ratio);
            return ratio;
 
        }
 
        /// construct the tensortrain representation of the operator
 
        /// @param[in]  source  source coefficient box
        /// @param[in]  shift   displacement
        /// @param[in]  tol     threshold for the TT truncation
        /// @param[in]  do_R    compute the R term of the operator (2k^d)
        /// @param[in]  do_T    compute the T term of the operator (k^d), including factor -1
        /// Both do_R and do_T may be used simultaneously, then the final
        /// operator will have dimensions (2k^d)
        TensorTrain<double> make_tt_representation(const Key<NDIM>& source,
                const Key<NDIM>& shift, double tol, bool do_R, bool do_T) const {
 
            if (not (do_R or do_T)) {
                print("no operator requested in make_tt_representation??");
                MADNESS_EXCEPTION("you're sure you know what you're doing?",1);
            }
 
 
            const SeparatedConvolutionData<Q,NDIM>* op = getop(source.level(), shift, source);
 
            // check for significant ranks since the R/T matrices' construction
            // might have been omitted. Tnorm is always smaller than Rnorm
            long lo=0,hi=rank;
            for (int mu=0; mu<rank; ++mu) {
                double Rnorm=1.0;
                for (std::size_t d=0; d<NDIM; ++d) Rnorm *= op->muops[mu].ops[d]->Rnorm;
                if (Rnorm>1.e-20) hi=mu;
                if ((Rnorm<1.e-20) and (mu<hi)) lo=mu;
            }
            hi++;lo++;
 
            // think about dimensions
            long rank_eff=(hi-lo);    // R or T matrices
            long step=1;
            if (do_R and do_T) {        // R and T matrices
                rank_eff*=2;
                step*=2;
            }
 
            long k2k=k;             // T matrices
            if (do_R) k2k=2*k;      // R matrices
 
 
            // construct empty TT cores and fill them with the significant R/T matrices
            std::vector<Tensor<double> > cores(NDIM,Tensor<double>(rank_eff,k2k,k2k,rank_eff));
            cores[0]=Tensor<double>(k2k,k2k,rank_eff);
            cores[NDIM-1]=Tensor<double>(rank_eff,k2k,k2k);
 
 
            for (int mu=lo, r=0; mu<hi; ++mu, ++r) {
                const SeparatedConvolutionInternal<Q,NDIM>& muop =  op->muops[mu];
                const Q fac = ops[mu].getfac();
                const Slice sr0(step*r,  step*r,  0);
                const Slice sr1(step*r+step-1,step*r+step-1,0);
                const Slice s00(0,k-1,1);
 
                if (do_R) {
                    cores[0](_,  _  ,sr0)=muop.ops[0]->R;
                    for (std::size_t idim=1; idim<NDIM-1; ++idim) {
                          cores[idim](sr0,_  ,_  ,sr0)=muop.ops[idim]->R;
                    }
                    cores[NDIM-1](sr0,_  ,_  )=muop.ops[NDIM-1]->R*fac;
                }
 
                if (do_T) {
                    cores[0](s00,s00,sr1)=muop.ops[0]->T;
                    for (std::size_t idim=1; idim<NDIM-1; ++idim) {
                        cores[idim](sr1,s00,s00,sr1)=muop.ops[idim]->T;
                    }
                    cores[NDIM-1](sr1,s00,s00)=muop.ops[NDIM-1]->T*(-fac);
                }
            }
 
            // construct TT representation
            TensorTrain<double> tt(cores);
 
            // need to reshape for the TT truncation
            tt.make_tensor();
            tt.truncate(tol*GenTensor<double>::fac_reduce());
            tt.make_operator();
 
            return tt;
        }
 
 
        static bool can_combine(const SeparatedConvolution<Q,NDIM>& left, const SeparatedConvolution<Q,NDIM>& right) {
            return (combine_OT(left,right).type!=OT_UNDEFINED);
        }
 
        /// return operator type and other info of the combined operator (e.g. fg = f(1,2)* g(1,2)
        static OperatorInfo combine_OT(const SeparatedConvolution<Q,NDIM>& left, const SeparatedConvolution<Q,NDIM>& right) {
            OperatorInfo info=left.info;
            if ((left.info.type==OT_F12) and (right.info.type==OT_G12)) {
                info.type=OT_FG12;
            } else if ((left.info.type==OT_GAUSS) and (right.info.type==OT_GAUSS)) {
                info=right.info;
                info.type=OT_GAUSS;
                info.mu=2.0*right.info.mu;
            } else if ((left.info.type==OT_SLATER) and (right.info.type==OT_SLATER)) {
                info=right.info;
                info.type=OT_SLATER;
                info.mu=2.0*right.info.mu;
            } else if ((left.info.type==OT_G12) and (right.info.type==OT_F12)) {
                info=right.info;
                info.type=OT_FG12;
            } else if ((left.info.type==OT_G12) and (right.info.type==OT_F212)) {
                info=right.info;
                info.type=OT_F2G12;
            } else if (((left.info.type==OT_F212) and (right.info.type==OT_G12)) or
                ((left.info.type==OT_F12) and (right.info.type==OT_FG12)) or
                ((left.info.type==OT_FG12) and (right.info.type==OT_F12))) {
                info=right.info;
                info.type=OT_F2G12;
                if (right.info.type!=OT_G12) MADNESS_CHECK(right.info.mu == left.info.mu);
            } else if ((left.info.type==OT_F12) and (right.info.type==OT_F12)) {
                info.type=OT_F212;
                // keep the original gamma
                // (f12)^2 = (1- slater12)^2  = 1/(4 gamma) (1 - 2 exp(-gamma) + exp(-2 gamma))
                MADNESS_CHECK(right.info.mu == left.info.mu);
            } else {
                MADNESS_EXCEPTION("unknown combination of SeparatedConvolutions: feel free to extend in operator.h",1);
            }
            return info;
        }
 
 
        /// combine 2 convolution operators to one
        /// There's some information loss with this constructor - not recommended
        static SeparatedConvolution<Q,NDIM> combine(const SeparatedConvolution<Q,NDIM>& left,
                                                    const SeparatedConvolution<Q,NDIM>& right) {
            MADNESS_CHECK(can_combine(left,right));
            MADNESS_CHECK(left.get_world().id()==right.get_world().id());
            MADNESS_CHECK(left.lattice_summed() == right.lattice_summed());
            std::array<LatticeRange, NDIM> lattice_summed;
            for (std::size_t i = 0; i < NDIM; ++i) {
              if (left.lattice_summed()[i]) lattice_summed[i].set_infinite();
            }
 
            auto info=combine_OT(left,right);
            return SeparatedConvolution<Q,NDIM>(left.get_world(),info,lattice_summed,left.k);
        }
 
        /// combine 2 convolution operators to one
        friend SeparatedConvolution<Q,NDIM> combine(const std::shared_ptr<SeparatedConvolution<Q,NDIM>> left,
                                                    const std::shared_ptr<SeparatedConvolution<Q,NDIM>> right) {
            SeparatedConvolution<Q,NDIM> result;
            if (left and right) {
                return combine(*left, *right);
            } else if (left) {
                return *left;
            } else if (right) {
                return *right;
            } else {
                MADNESS_EXCEPTION("can't combine empty SeparatedConvolutions",1);
            }
            return result;
        }
    };
 
 
 
    /// Factory function generating separated kernel for convolution with 1/r in 3D.
    /// N.B. bloch_k species how the function's phase changes as the **operator** translates,
    /// which (by translational invariance) is the additive inverse of how the function's
    /// phase changes as the **cell** shifts
    static
    inline
    SeparatedConvolution<double_complex,3> PeriodicHFExchangeOperator(World& world,
                                                   Vector<double,3> bloch_k,
                                                   double lo,
                                                   double eps,
                                                   const std::array<KernelRange, 3>& kernel_ranges = std::array<KernelRange, 3>(),
                                                   const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                                                   int k=FunctionDefaults<3>::get_k()) {
      return SeparatedConvolution<double_complex, 3>(world, OperatorInfo(0.0, lo, eps, OT_G12, std::nullopt, kernel_ranges),
                                                     lattice_ranges, k, false, bloch_k);
    }
 
    /// Factory function generating separated kernel for convolution with 1/r in 3D.
    /// N.B. bloch_k species how the function's phase changes as the **operator** translates,
    /// which (by translational invariance) is the additive inverse of how the function's
    /// phase changes as the **cell** shifts
    static
    inline
    SeparatedConvolution<double_complex,3>* PeriodicHFExchangeOperatorPtr(World& world,
                                                   Vector<double,3> bloch_k,
                                                   double lo,
                                                   double eps,
                                                   const std::array<KernelRange, 3>& kernel_ranges = std::array<KernelRange, 3>(),
                                                   const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                                                   int k=FunctionDefaults<3>::get_k()) {
      return new SeparatedConvolution<double_complex, 3>(world, OperatorInfo(0.0, lo, eps, OT_G12, std::nullopt, kernel_ranges),
                                                         lattice_ranges, k, false, bloch_k);
    }
 
    /// Factory function generating separated kernel for convolution with 1/r in 3D.
    static
    inline
    SeparatedConvolution<double,3> CoulombOperator(World& world,
                                                   double lo,
                                                   double eps,
                                                   const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                           int k=FunctionDefaults<3>::get_k())
    {
        return SeparatedConvolution<double,3>(world,OperatorInfo(0.0,lo,eps,OT_G12),lattice_ranges,k);
    }
 
 
    /// Factory function generating separated kernel for convolution with 1/r in 3D.
    static
    inline
    SeparatedConvolution<double,3>* CoulombOperatorPtr(World& world,
                                                       double lo,
                                                       double eps,
                                                       const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                                                       int k=FunctionDefaults<3>::get_k())
    {
        return new SeparatedConvolution<double,3>(world,OperatorInfo(0.0,lo,eps,OT_G12),lattice_ranges,k);
    }
 
 
    /// Factory function generating separated kernel for convolution with BSH kernel in general NDIM
    template <std::size_t NDIM>
    static inline
    SeparatedConvolution<double,NDIM>
    BSHOperator(World& world, double mu, double lo, double eps,
    const std::array<LatticeRange, NDIM>& lattice_ranges = FunctionDefaults<NDIM>::get_bc().lattice_range(),
                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),lattice_ranges,k);
    }
 
    /// Factory function generating separated kernel for convolution with BSH kernel in general NDIM
    template <std::size_t NDIM>
    static inline
    SeparatedConvolution<double,NDIM>*
    BSHOperatorPtr(World& world, double mu, double lo, double eps,
      const std::array<LatticeRange, NDIM>& lattice_ranges = FunctionDefaults<NDIM>::get_bc().lattice_range(),
                   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),lattice_ranges,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 std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                  int k=FunctionDefaults<3>::get_k()) {
        return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_BSH),lattice_ranges,k);
    }
 
    /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D
    static
    inline
    SeparatedConvolution<double_complex,3> PeriodicBSHOperator3D(World& world,
                                                         Vector<double,3> bloch_k,
                                                         double mu,
                                                         double lo,
                                                         double eps,
                                                         const std::array<KernelRange, 3>& kernel_ranges = std::array<KernelRange, 3>(),
                                                         const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                                                         int k=FunctionDefaults<3>::get_k())
 
    {
      return SeparatedConvolution<double_complex, 3>(world, OperatorInfo(mu, lo, eps, OT_BSH, std::nullopt, kernel_ranges),
                                                     lattice_ranges, k, false, bloch_k);
    }
 
    /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D
    static inline
    SeparatedConvolution<double_complex,3>* PeriodicBSHOperatorPtr3D(World& world,
                                                         Vector<double,3> bloch_k,
                                                         double mu,
                                                         double lo,
                                                         double eps,
                                                         const std::array<KernelRange, 3>& kernel_ranges = std::array<KernelRange, 3>(),
                                                         const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                                                         int k=FunctionDefaults<3>::get_k())
 
    {
      return new SeparatedConvolution<double_complex, 3>(world, OperatorInfo(mu, lo, eps, OT_BSH, std::nullopt, kernel_ranges),
                                                     lattice_ranges, k, false, bloch_k);
    }
 
 
    static inline SeparatedConvolution<double,3>
    SlaterF12Operator(World& world, double mu, double lo, double eps,
                      const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                      int k=FunctionDefaults<3>::get_k()) {
        return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F12),lattice_ranges,k);
    }
 
    static inline SeparatedConvolution<double,3> SlaterF12sqOperator(World& world,
                                                                   double mu, double lo, double eps,
                                                                   const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                                                                   int k=FunctionDefaults<3>::get_k()) {
        return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F212),lattice_ranges,k);
    }
 
    static inline SeparatedConvolution<double,3>* SlaterF12sqOperatorPtr(World& world,
                                                                       double mu, double lo, double eps,
                                                                       const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                                                                       int k=FunctionDefaults<3>::get_k()) {
        return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F212),lattice_ranges,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 std::array<LatticeRange, NDIM>& lattice_ranges = FunctionDefaults<NDIM>::get_bc().lattice_range(),
            int k=FunctionDefaults<NDIM>::get_k()) {
        return SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_SLATER),lattice_ranges,k);
    }
 
    /// Factory function generating separated kernel for convolution with exp(-mu*r*r)
 
    /// lo and eps are not used here
    template<std::size_t NDIM>
    static inline SeparatedConvolution<double,NDIM> GaussOperator(World& world,
                                                                double mu, double lo=0.0, double eps=0.0,
                                                                const std::array<LatticeRange, NDIM>& lattice_ranges = FunctionDefaults<NDIM>::get_bc().lattice_range(),
                                                                int k=FunctionDefaults<NDIM>::get_k()) {
        return SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_GAUSS),lattice_ranges,k);
    }
 
    /// Factory function generating separated kernel for convolution with exp(-mu*r*r) in 3D
 
    /// lo and eps are not used here
    template<std::size_t NDIM>
    static inline SeparatedConvolution<double, NDIM>* GaussOperatorPtr(World& world,
                                                                     double mu, double lo=0.0, double eps=0.0,
                                                                     const std::array<LatticeRange, NDIM>& lattice_ranges = FunctionDefaults<NDIM>::get_bc().lattice_range(),
                                                                     int k = FunctionDefaults<NDIM>::get_k()) {
        return new SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_GAUSS),lattice_ranges,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
    template<std::size_t NDIM>
    static inline SeparatedConvolution<double, NDIM>* SlaterOperatorPtr_ND(World& world,
                                                                     double mu, double lo, double eps,
                                                                     const std::array<LatticeRange, NDIM>& lattice_ranges = FunctionDefaults<NDIM>::get_bc().lattice_range(),
                                                                     int k = FunctionDefaults<NDIM>::get_k()) {
        return new SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_SLATER),lattice_ranges,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 std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                                                                 int k = FunctionDefaults<3>::get_k()) {
        return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_SLATER),lattice_ranges,k);
    }
 
    /// Factory function generating separated kernel for convolution with (1 - exp(-mu*r))/(2 mu) in 3D
 
    /// includes the factor 1/(2 mu)
    static inline SeparatedConvolution<double,3>* SlaterF12OperatorPtr(World& world,
            double mu, double lo, double eps,
            const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
            int k=FunctionDefaults<3>::get_k()) {
        return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F12),lattice_ranges,k);
    }
 
 
    /// Factory function generating separated kernel for convolution with 1/(2 mu)*(1 - exp(-mu*r))/r in 3D
 
    /// fg = (1 - exp(-gamma r12))  / r12 = 1/r12 - exp(-gamma r12)/r12 = coulomb - bsh
    /// includes the factor 1/(2 mu)
    static inline SeparatedConvolution<double,3>
    FGOperator(World& world, double mu, double lo, double eps,
    const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
               int k=FunctionDefaults<3>::get_k()) {
        return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_FG12),lattice_ranges,k);
    }
 
    /// Factory function generating separated kernel for convolution with 1/(2 mu)*(1 - exp(-mu*r))/r in 3D
 
    /// fg = (1 - exp(-gamma r12))  / r12 = 1/r12 - exp(-gamma r12)/r12 = coulomb - bsh
    /// includes the factor 1/(2 mu)
    static inline SeparatedConvolution<double,3>*
    FGOperatorPtr(World& world, double mu, double lo, double eps,
    const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                  int k=FunctionDefaults<3>::get_k()) {
        return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_FG12),lattice_ranges,k);
    }
 
    /// Factory function generating separated kernel for convolution with (1/(2 mu)*(1 - exp(-mu*r)))^2/r in 3D
 
    /// f2g = (1/(2 gamma) (1 - exp(-gamma r12)))^2  / r12
    ///     = 1/(4 gamma) * [ 1/r12 - 2 exp(-gamma r12)/r12 + exp(-2 gamma r12)/r12 ]
    /// includes the factor 1/(2 mu)^2
    static inline SeparatedConvolution<double,3>*
    F2GOperatorPtr(World& world, double mu, double lo, double eps,
    const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                   int k=FunctionDefaults<3>::get_k()) {
        return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F2G12),lattice_ranges,k);
    }
 
    /// Factory function generating separated kernel for convolution with (1/(2 mu)*(1 - exp(-mu*r)))^2/r in 3D
 
    /// f2g = (1/(2 gamma) (1 - exp(-gamma r12)))^2  / r12
    ///     = 1/(4 gamma) * [ 1/r12 - 2 exp(-gamma r12)/r12 + exp(-2 gamma r12)/r12 ]
    /// includes the factor 1/(2 mu)^2
    static inline SeparatedConvolution<double,3>
    F2GOperator(World& world, double mu, double lo, double eps,
    const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                int k=FunctionDefaults<3>::get_k()) {
        return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F2G12),lattice_ranges,k);
    }
 
 
    /// Factory function generating separated kernel for convolution a normalized
    /// Gaussian (aka a widened delta function)
    static inline SeparatedConvolution<double,3> SmoothingOperator3D(World& world,
            double eps,
            const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
            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, lattice_ranges, k);
    }
 
    /// Factory function generating separated kernel for convolution a normalized
    /// Gaussian (aka a widened delta function)
    template<std::size_t NDIM>
    static inline SeparatedConvolution<double,NDIM> SmoothingOperator(World& world,
            double eps,
            const std::array<LatticeRange, NDIM>& lattice_ranges = FunctionDefaults<NDIM>::get_bc().lattice_range(),
            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, lattice_ranges, k);
    }
 
    /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D
    static
    inline
    SeparatedConvolution<double,3>* BSHOperatorPtr3D(World& world,
                                                     double mu,
                                                     double lo,
                                                     double eps,
                                                     const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                                                     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
      // N.B. if have periodic boundaries, extend range just in case will be using periodic domain
      const auto lattice_summed_any = std::any_of(lattice_ranges.begin(), lattice_ranges.end(), [](const auto& b) { return static_cast<bool>(b);});
      const auto infinite_any = std::any_of(lattice_ranges.begin(), lattice_ranges.end(), [](const auto& b) { return b.infinite();});
      const double hi_fin = lattice_finite_reach(lattice_ranges, cell_width);
      if (lattice_summed_any) {
        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 (infinite_any) {
        fit.truncate_mixed_expansion(coeff, expnt, lattice_ranges, cell_width, lo, hi_fin, eps);
      }
      return new SeparatedConvolution<double, 3>(world, coeff, expnt, lo, eps,
                                                 lattice_ranges, k);
    }
 
 
    /// Factory function generating operator for convolution with grad(1/r) in 3D
 
    /// Returns a 3-vector containing the convolution operator for the
    /// x, y, and z components of grad(1/r)
    static
    inline
    std::vector< std::shared_ptr< SeparatedConvolution<double,3> > >
    GradCoulombOperator(World& world,
                        double lo,
                        double eps,
                        const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                        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 lattice_sum_any = std::any_of(lattice_ranges.begin(), lattice_ranges.end(), [](const LatticeRange& b){ return static_cast<bool>(b); });
      const auto infinite_any = std::any_of(lattice_ranges.begin(), lattice_ranges.end(), [](const LatticeRange& b){ return b.infinite(); });
      if (lattice_sum_any) {
        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 (infinite_any) {
        fit.truncate_mixed_expansion(coeff, expnt, lattice_ranges, width, lo, lattice_finite_reach(lattice_ranges, width), eps);
      }
 
      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_ranges[d]));
          }
          ops[mu].setop(dir, GaussianConvolution1DCache<double>::get(
                                 k, expnt(mu) * width[dir] * width[dir], 1,
                                 lattice_ranges[dir]));
        }
        gradG[dir] = real_convolution_3d_ptr(
            new SeparatedConvolution<double, 3>(world, ops));
      }
 
      return gradG;
    }
 
    /// Factory function generating operator for convolution with grad(bsh) in 3D
 
    /// Returns a 3-vector containing the convolution operator for the
    /// x, y, and z components of grad(bsh)
    static
    inline
    std::vector< std::shared_ptr< SeparatedConvolution<double,3> > >
    GradBSHOperator(World& world,
                        double mu,
                        double lo,
                        double eps,
                        const std::array<LatticeRange, 3>& lattice_ranges = FunctionDefaults<3>::get_bc().lattice_range(),
                        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
      bool lattice_sum_any = std::any_of(lattice_ranges.begin(), lattice_ranges.end(), [](const LatticeRange& b){ return b.get_range(); });
      const auto infinite_any = std::any_of(lattice_ranges.begin(), lattice_ranges.end(), [](const LatticeRange& b){ return b.infinite(); });
      if (lattice_sum_any) {
        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 (infinite_any) {
        fit.truncate_mixed_expansion(coeff, expnt, lattice_ranges, width, lo, lattice_finite_reach(lattice_ranges, width), eps);
      }
 
      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_ranges[d]));
          }
          ops[mu].setop(dir, GaussianConvolution1DCache<double>::get(
                                 k, expnt(mu) * width[dir] * width[dir], 1,
                                 lattice_ranges[dir]));
        }
        gradG[dir] = real_convolution_3d_ptr(
            new SeparatedConvolution<double, 3>(world, ops));
      }
 
      return gradG;
    }
 
 
 
    namespace archive {
        template <class Archive, class T, std::size_t NDIM>
        struct ArchiveLoadImpl<Archive,const SeparatedConvolution<T,NDIM>*> {
            static inline void load(const Archive& ar, const SeparatedConvolution<T,NDIM>*& ptr) {
                WorldObject< SeparatedConvolution<T,NDIM> >* p = nullptr;
                ar & p;
                ptr = static_cast< const SeparatedConvolution<T,NDIM>* >(p);
            }
        };
 
        template <class Archive, class T, std::size_t NDIM>
        struct ArchiveStoreImpl<Archive,const SeparatedConvolution<T,NDIM>*> {
            static inline void store(const Archive& ar, const SeparatedConvolution<T,NDIM>*const& ptr) {
                ar & static_cast< const WorldObject< SeparatedConvolution<T,NDIM> >* >(ptr);
            }
        };
    }
 
}
 
 
 
 
#endif // MADNESS_MRA_OPERATOR_H__INCLUDED