diff --git a/src/examples/periodic/testfilter.cc b/src/examples/periodic/testfilter.cc index ee63e4d7b9d..1295131fa70 100644 --- a/src/examples/periodic/testfilter.cc +++ b/src/examples/periodic/testfilter.cc @@ -100,7 +100,7 @@ int main(int argc, char**argv) { rhoN.truncate(); BoundaryConditions<3> bc_periodic(BC_PERIODIC); - SeparatedConvolution pop = CoulombOperator(world, 1e-4, eps, bc_periodic); + SeparatedConvolution pop = CoulombOperator(world, 1e-4, eps, bc_periodic.is_periodic()); printf("applying periodic operator ...\n\n"); V_periodic = apply(pop, rho); V_periodic.truncate(); diff --git a/src/madness/mra/bc.h b/src/madness/mra/bc.h index 3259e14d725..5fdc62393f2 100644 --- a/src/madness/mra/bc.h +++ b/src/madness/mra/bc.h @@ -39,9 +39,9 @@ #include +#include #include #include -#include namespace madness { @@ -126,11 +126,14 @@ template class BoundaryConditions { /// Convenience for application of integral operators - /// @return Returns a vector indicating if each dimension is periodic - std::vector is_periodic() const { - std::vector v(NDIM); - for (std::size_t d = 0; d < NDIM; ++d) + /// @return Returns a vector indicating if dimensions [0, ND) are periodic + template + std::enable_if_t> is_periodic() const { + std::array v; + for (std::size_t d = 0; d < ND; ++d) { + MADNESS_ASSERT(bc[2 * d + 1] == bc[2 * d]); v[d] = (bc[2 * d] == BC_PERIODIC); + } return v; } @@ -138,11 +141,25 @@ template class BoundaryConditions { /// @return Returns true if any dimension is periodic bool is_periodic_any() const { - for (std::size_t d = 0; d < NDIM; ++d) + for (std::size_t d = 0; d < NDIM; ++d) { + MADNESS_ASSERT(bc[2 * d + 1] == bc[2 * d]); if (bc[2 * d] == BC_PERIODIC) return true; + } return false; } + + /// Checks whether the boundary condition along all axes is periodic + + /// @return Returns true if every dimension is periodic + bool is_periodic_all() const { + for (std::size_t d = 0; d < NDIM; ++d) { + MADNESS_ASSERT(bc[2 * d + 1] == bc[2 * d]); + if (bc[2 * d] != BC_PERIODIC) + return false; + } + return true; + } }; template @@ -159,6 +176,13 @@ static inline std::ostream &operator<<(std::ostream &s, return s; } +template +std::array no_lattice_sum() { + std::array result; + result.fill(false); + return result; +} + } // namespace madness #endif // MADNESS_MRA_BC_H__INCLUDED diff --git a/src/madness/mra/funcimpl.h b/src/madness/mra/funcimpl.h index c57261e35cc..2985c44f6cb 100644 --- a/src/madness/mra/funcimpl.h +++ b/src/madness/mra/funcimpl.h @@ -35,8 +35,6 @@ /// \file funcimpl.h /// \brief Provides FunctionCommonData, FunctionImpl and FunctionFactory -#include -#include #include #include #include @@ -50,7 +48,11 @@ #include #include -#include "leafop.h" +#include + +#include +#include +#include namespace madness { template @@ -3820,7 +3822,7 @@ template /// Out of volume keys are mapped to enforce the BC as follows. /// * Periodic BC map back into the volume and return the correct key /// * non-periodic BC - returns invalid() to indicate out of volume - keyT neighbor(const keyT& key, const keyT& disp, const std::vector& is_periodic) const; + keyT neighbor(const keyT& key, const keyT& disp, const std::array& is_periodic) const; /// Returns key of general neighbor that resides in-volume @@ -4524,7 +4526,7 @@ template void zero_norm_tree(); // Broaden tree - void broaden(std::vector is_periodic, bool fence); + void broaden(const std::array& is_periodic, bool fence); /// sum all the contributions from all scales after applying an operator in mod-NS form void trickle_down(bool fence); @@ -4810,10 +4812,25 @@ template double cnorm = c.normf(); - const std::vector& disp = Displacements().get_disp(key.level(), /* periodic = */ false); // list of displacements sorted in order of increasing distance; N.B. periodic displacements skipped since operator already includes lattice sum + // BC handling: + // - if operator is lattice-summed then treat this as nonperiodic (i.e. tell neighbor() to stay in simulation cell) + // - if operator is NOT lattice-summed then obey BC (i.e. tell neighbor() to go outside the simulation cell along periodic dimensions) + const auto bc_is_periodic = FunctionDefaults::get_bc().is_periodic(); + const auto lattice_sum_in_op = op->includes_lattice_sum(); + std::array treat_this_as_periodic; + for(auto axis=0; axis!=NDIM; ++axis) { + treat_this_as_periodic[axis] = + bc_is_periodic[axis] && !lattice_sum_in_op[axis]; + } + // if this is treated as periodic in ANY direction need to use displacements that are periodic in EVERY direction + // because these are the only flavor of periodic displacements that we have + // don't worry, neighbor() will still treat nonperiodic axes as nonperiodic + const bool treat_this_as_periodic_in_any_direction = std::accumulate(treat_this_as_periodic.begin(), treat_this_as_periodic.end(), false, std::logical_or{}); + const std::vector& disp = Displacements().get_disp(key.level(), /* periodic = */ treat_this_as_periodic_in_any_direction); // list of displacements sorted in order of increasing distance + int nvalid=1; // Counts #valid at each distance int nused=1; // Counts #used at each distance - uint64_t distsq = 99999999999999; + uint64_t distsq = 99999999999999; for (typename std::vector::const_iterator it=disp.begin(); it != disp.end(); ++it) { keyT d; Key nullkey(key.level()); @@ -4834,7 +4851,7 @@ template distsq = dsq; } - keyT dest = neighbor_in_volume(key, d); + keyT dest = neighbor(key, d, treat_this_as_periodic); if (dest.is_valid()) { nvalid++; double opnorm = op->norm(key.level(), *it, source); diff --git a/src/madness/mra/key.h b/src/madness/mra/key.h index 47e2126c207..b2b730d2597 100644 --- a/src/madness/mra/key.h +++ b/src/madness/mra/key.h @@ -247,7 +247,7 @@ namespace madness { /// Assumes key and this are at the same level bool - is_neighbor_of(const Key& key, const std::vector& bperiodic) const { + is_neighbor_of(const Key& key, const std::array& bperiodic) const { Translation dist = 0; Translation TWON1 = (Translation(1)< 1 and box_is_at_boundary(key)) return false; // check for all special points if they are neighbours of the current box BoundaryConditions bc = FunctionDefaults::get_bc(); - std::vector bperiodic = bc.is_periodic(); + const auto bperiodic = bc.is_periodic(); for (size_t i = 0; i < special_points.size(); ++i) { Vector simpt; user_to_sim(special_points[i], simpt); @@ -151,7 +151,7 @@ struct ElectronCuspyBox_op : public Specialbox_op { std::string name() const { return "Cuspybox_op"; } /// Operator which decides if the key belongs to a special box - /// The key is broken appart in two lower dimensional keys (for electron cusps this is 6D -> 2x3D) + /// The key is broken apart in two lower dimensional keys (for electron cusps this is 6D -> 2x3D) /// if the keys are neighbours then refinement up to the special level is enforced (means the 6D box is close to the cusp or contains it) /// if the refinement level is already beyond half of the special_level then refinement is only enforded if the broken keys are the same (6D box contains cusp) bool operator()(const Key& key, const FunctionImpl *const f) const { @@ -160,21 +160,27 @@ struct ElectronCuspyBox_op : public Specialbox_op { if (this->box_is_at_boundary(key)) return false; } - if (NDIM % 2 != 0) MADNESS_EXCEPTION("Cuspybox_op only valid for even dimensions", - 1); // if uneven dims are needed just make a new class with NDIM+1/2 and NDIM-LDIM - BoundaryConditions bc = FunctionDefaults::get_bc(); - std::vector bperiodic = bc.is_periodic(); - Key key1; - Key key2; - key.break_apart(key1, key2); - int ll = this->get_half_of_special_level(); - if (ll < f->get_initial_level()) ll = f->get_initial_level(); - if (key.level() > ll) { - if (key1 == key2) return true; - else return false; - } else { - if (key1.is_neighbor_of(key2, bperiodic)) return true; - else return false; + // Cuspybox_op only valid for even dimensions, if uneven dims are needed just make a new class with NDIM+1/2 and NDIM-LDIM + if constexpr (NDIM % 2 == 0) { + BoundaryConditions bc = FunctionDefaults::get_bc(); + const auto bperiodic = bc.template is_periodic(); + Key key1; + Key key2; + key.break_apart(key1, key2); + int ll = this->get_half_of_special_level(); + if (ll < f->get_initial_level()) + ll = f->get_initial_level(); + if (key.level() > ll) { + if (key1 == key2) + return true; + else + return false; + } else { + if (key1.is_neighbor_of(key2, bperiodic)) + return true; + else + return false; + } } MADNESS_EXCEPTION("We should not end up here (check further of cuspy box)", 1); return false; @@ -239,7 +245,7 @@ struct NuclearCuspyBox_op : public Specialbox_op { // now break the key appart and check if one if the results is in the neighbourhood of a special point BoundaryConditions bc = FunctionDefaults::get_bc(); - std::vector bperiodic = bc.is_periodic(); + const auto bperiodic = bc.is_periodic(); Key key1; Key key2; @@ -346,7 +352,7 @@ class Leaf_op { sanity(); const double cnorm = coeff.normf(); BoundaryConditions bc = FunctionDefaults::get_bc(); - std::vector bperiodic = bc.is_periodic(); + const auto bperiodic = bc.is_periodic(); typedef Key opkeyT; const opkeyT source = op->get_source_key(key); diff --git a/src/madness/mra/mraimpl.h b/src/madness/mra/mraimpl.h index 52a739683e3..385470aafdd 100644 --- a/src/madness/mra/mraimpl.h +++ b/src/madness/mra/mraimpl.h @@ -1279,7 +1279,7 @@ namespace madness { // Broaden tree template - void FunctionImpl::broaden(std::vector is_periodic, bool fence) { + void FunctionImpl::broaden(const std::array& is_periodic, bool fence) { typename dcT::iterator end = coeffs.end(); for (typename dcT::iterator it=coeffs.begin(); it!=end; ++it) { const keyT& key = it->first; @@ -2491,7 +2491,7 @@ namespace madness { std::vector > newspecialpts; if (key.level() < functor->special_level() && specialpts.size() > 0) { BoundaryConditions bc = FunctionDefaults::get_bc(); - std::vector bperiodic = bc.is_periodic(); + const auto bperiodic = bc.is_periodic(); for (unsigned int i = 0; i < specialpts.size(); ++i) { coordT simpt; user_to_sim(specialpts[i], simpt); @@ -3239,7 +3239,7 @@ template } template - Key FunctionImpl::neighbor(const keyT& key, const Key& disp, const std::vector& is_periodic) const { + Key FunctionImpl::neighbor(const keyT& key, const Key& disp, const std::array& is_periodic) const { Vector l = key.translation(); for (std::size_t axis=0; axis lattice_sum; ///< If lattice_sum[d] is true summation over lattice translations along axis d + ///< N.B. the resulting kernel can be non-zero at both ends of the simulation cell along that axis bool modified_=false; ///< use modified NS form int particle_=1; ///< must only be 1 or 2 bool destructive_=false; ///< destroy the argument or restore it (expensive for 6d functions) @@ -159,7 +159,6 @@ namespace madness { mutable std::vector< ConvolutionND > ops; ///< ConvolutionND keeps data for 1 term, all dimensions, 1 displacement - const BoundaryConditions bc; const int k; const FunctionCommonData& cdata; int rank; @@ -209,10 +208,10 @@ namespace madness { static inline std::pair,Tensor> make_coeff_for_operator(World& world, double mu, double lo, double eps, OpType type, - const BoundaryConditions& bc=FunctionDefaults::get_bc()) { + const std::array& lattice_sum) { OperatorInfo info(mu,lo,eps,type); - return make_coeff_for_operator(world, info, bc); + return make_coeff_for_operator(world, info, lattice_sum); // const Tensor& cell_width = FunctionDefaults<3>::get_cell_width(); // double hi = cell_width.normf(); // Diagonal width of cell // if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation @@ -233,26 +232,32 @@ namespace madness { static inline std::pair,Tensor> make_coeff_for_operator(World& world, OperatorInfo& info, - const BoundaryConditions& bc=FunctionDefaults::get_bc()) { - - const Tensor& cell_width = FunctionDefaults<3>::get_cell_width(); - double hi = cell_width.normf(); // Diagonal width of cell - if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation - - info.hi=hi; - GFit fit(info); - - Tensor coeff=fit.coeffs(); - Tensor expnt=fit.exponents(); - - // WARNING! More fine-grained control over the last argument is needed. - // This is a hotfix. - if (info.truncate_lowexp_gaussians.value_or(bc(0,0) == BC_PERIODIC)) { - fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), true); - info.truncate_lowexp_gaussians = true; - } - - return std::make_pair(coeff,expnt); + const std::array& lattice_sum) { + + const Tensor &cell_width = + FunctionDefaults<3>::get_cell_width(); + double hi = cell_width.normf(); // Diagonal width of cell + // Extend kernel range for lattice summation + const auto any_lattice_sum = std::accumulate( + lattice_sum.begin(), lattice_sum.end(), false, std::logical_or{}); + if (any_lattice_sum) { + hi *= 100; + } + + info.hi = hi; + GFit fit(info); + + Tensor coeff = fit.coeffs(); + Tensor expnt = fit.exponents(); + + // WARNING! More fine-grained control over the last argument is needed. This is a hotfix. + if (info.truncate_lowexp_gaussians.value_or(any_lattice_sum)) { + fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), + true); + info.truncate_lowexp_gaussians = true; + } + + return std::make_pair(coeff, expnt); } // /// return the right block of the upsampled operator (modified NS only) @@ -947,17 +952,16 @@ namespace madness { // For separated convolutions with same operator in each direction (isotropic) SeparatedConvolution(World& world, std::vector< std::shared_ptr< Convolution1D > >& argops, - const BoundaryConditions& bc = FunctionDefaults::get_bc(), + const std::array& lattice_sum = FunctionDefaults::get_bc().is_periodic(), long k = FunctionDefaults::get_k(), bool doleaves = false) : WorldObject< SeparatedConvolution >(world) , info() , doleaves(doleaves) - , isperiodicsum(bc(0,0)==BC_PERIODIC) + , lattice_sum(lattice_sum) , modified_(false) , particle_(1) , destructive_(false) - , bc(bc) , k(k) , cdata(FunctionCommonData::get(k)) , rank(argops.size()) @@ -965,10 +969,6 @@ namespace madness { , v2k(NDIM,2*k) , s0(std::max(2,NDIM),Slice(0,k-1)) { - // Presently we must have periodic or non-periodic in all dimensions. - for (std::size_t d=1; d >& argops, - const BoundaryConditions& bc = FunctionDefaults::get_bc(), + const std::array& lattice_sum = FunctionDefaults::get_bc().is_periodic(), long k = FunctionDefaults::get_k(), bool doleaves = false) : WorldObject< SeparatedConvolution >(world) , info() , doleaves(doleaves) - , isperiodicsum(bc(0,0)==BC_PERIODIC) + , lattice_sum(lattice_sum) , modified_(false) , particle_(1) , destructive_(false) , ops(argops) - , bc(bc) , k(k) , cdata(FunctionCommonData::get(k)) , rank(argops.size()) @@ -1000,22 +999,18 @@ namespace madness { , v2k(NDIM,2*k) , s0(std::max(2,NDIM),Slice(0,k-1)) { - // Presently we must have periodic or non-periodic in all dimensions. - for (std::size_t d=1; dprocess_pending(); } /// Constructor for Gaussian Convolutions (mostly for backward compatability) SeparatedConvolution(World& world, const OperatorInfo info1, - const BoundaryConditions& bc = FunctionDefaults::get_bc(), + const std::array& lattice_sum = FunctionDefaults::get_bc().is_periodic(), int k=FunctionDefaults::get_k(), bool doleaves = false) - : SeparatedConvolution(world,Tensor(0l),Tensor(0l),info1.lo,info1.thresh,bc,k,doleaves,info1.mu) { + : SeparatedConvolution(world,Tensor(0l),Tensor(0l),info1.lo,info1.thresh,lattice_sum,k,doleaves,info1.mu) { info.type=info1.type; info.truncate_lowexp_gaussians = info1.truncate_lowexp_gaussians; - auto [coeff, expnt] = make_coeff_for_operator(world, info, bc); + auto [coeff, expnt] = make_coeff_for_operator(world, info, lattice_sum); rank=coeff.dim(0); ops.resize(rank); initialize(coeff,expnt); @@ -1025,16 +1020,15 @@ namespace madness { SeparatedConvolution(World& world, const Tensor& coeff, const Tensor& expnt, double lo, double thresh, - const BoundaryConditions& bc = FunctionDefaults::get_bc(), + const std::array& lattice_sum = FunctionDefaults::get_bc().is_periodic(), int k=FunctionDefaults::get_k(), bool doleaves = false, double mu=0.0) : WorldObject< SeparatedConvolution >(world) , info(mu,lo,thresh,OT_UNDEFINED) , doleaves(doleaves) - , isperiodicsum(bc(0,0)==BC_PERIODIC) + , lattice_sum(lattice_sum) , ops(coeff.dim(0)) - , bc(bc) , k(k) , cdata(FunctionCommonData::get(k)) , rank(coeff.dim(0)) @@ -1045,11 +1039,6 @@ namespace madness { } void initialize(const Tensor& coeff, const Tensor& expnt) { - // Presently we must have periodic or non-periodic in all dimensions. - for (std::size_t d=1; d& width = FunctionDefaults::get_cell_width(); const double pi = constants::pi; @@ -1061,7 +1050,7 @@ namespace madness { ops[mu].setfac(coeff(mu)/c); for (std::size_t d=0; d::get(k, expnt(mu)*width[d]*width[d], 0, isperiodicsum)); + ops[mu].setop(d,GaussianConvolution1DCache::get(k, expnt(mu)*width[d]*width[d], 0, lattice_sum[d])); } } } @@ -1070,18 +1059,17 @@ namespace madness { SeparatedConvolution(World& world, Vector args, const Tensor& coeff, const Tensor& expnt, - const BoundaryConditions& bc = FunctionDefaults::get_bc(), + const std::array& lattice_sum = FunctionDefaults::get_bc().is_periodic(), int k=FunctionDefaults::get_k(), bool doleaves=false) : WorldObject< SeparatedConvolution >(world) , info(0.0,0.0,0.0,OT_UNDEFINED) , doleaves(doleaves) - , isperiodicsum(bc(0,0)==BC_PERIODIC) + , lattice_sum(lattice_sum) , modified_(false) , particle_(1) , destructive_(false) , ops(coeff.dim(0)) - , bc(bc) , k(k) , cdata(FunctionCommonData::get(k)) , rank(coeff.dim(0)) @@ -1089,11 +1077,6 @@ namespace madness { , v2k(NDIM,2*k) , s0(std::max(2,NDIM),Slice(0,k-1)) { - // Presently we must have periodic or non-periodic in all dimensions. - for (std::size_t d=1; d& width = FunctionDefaults::get_cell_width(); for (int mu=0; mu > gcptr(new GaussianConvolution1D(k, c2, - expnt(mu)*width[d]*width[d], 0, isperiodicsum, args[d])); + expnt(mu)*width[d]*width[d], 0, lattice_sum[d], args[d])); ops[mu].setop(d,gcptr); } } @@ -1127,13 +1110,11 @@ namespace madness { } } - const BoundaryConditions& get_bc() const {return bc;} - const std::vector< Key >& get_disp(Level n) const { - return Displacements().get_disp(n, isperiodicsum); + return Displacements().get_disp(n, false); } - bool get_isperiodicsum() const {return isperiodicsum;} + const std::array& includes_lattice_sum() const { return lattice_sum; } /// return the operator norm for all terms, all dimensions and 1 displacement double norm(Level n, const Key& d, const Key& source_key) const { @@ -1697,9 +1678,10 @@ namespace madness { const SeparatedConvolution& right) { MADNESS_CHECK(can_combine(left,right)); MADNESS_CHECK(left.get_world().id()==right.get_world().id()); + MADNESS_CHECK(left.includes_lattice_sum() == right.includes_lattice_sum()); auto info=combine_OT(left,right); - return SeparatedConvolution(left.get_world(),info,left.bc,left.k); + return SeparatedConvolution(left.get_world(),info,left.includes_lattice_sum(),left.k); } /// combine 2 convolution operators to one @@ -1728,24 +1710,28 @@ namespace madness { Vector args, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), - int k=FunctionDefaults<3>::get_k()) - { - const Tensor& cell_width = FunctionDefaults<3>::get_cell_width(); - double hi = cell_width.normf(); // Diagonal width of cell - if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation - - GFit fit=GFit::CoulombFit(lo,hi,eps,false); - Tensor coeff=fit.coeffs(); - Tensor expnt=fit.exponents(); - - if (bc(0,0) == BC_PERIODIC) { - fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), true); - } - - return SeparatedConvolution(world, args, coeff, expnt, bc, k, false); -// return SeparatedConvolution(world, coeff, expnt, bc, k); - + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), + int k=FunctionDefaults<3>::get_k()) { + const Tensor &cell_width = FunctionDefaults<3>::get_cell_width(); + double hi = cell_width.normf(); // Diagonal width of cell + + // Extend kernel range for lattice summation + const auto any_lattice_sum = std::accumulate( + lattice_sum.begin(), lattice_sum.end(), false, std::logical_or{}); + if (any_lattice_sum) { + hi *= 100; + } + + GFit fit = GFit::CoulombFit(lo, hi, eps, false); + Tensor coeff = fit.coeffs(); + Tensor expnt = fit.exponents(); + + if (any_lattice_sum) { + fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), true); + } + + return SeparatedConvolution(world, args, coeff, expnt, + lattice_sum, k, false); } /// Factory function generating separated kernel for convolution with 1/r in 3D. @@ -1754,11 +1740,10 @@ namespace madness { SeparatedConvolution CoulombOperator(World& world, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), - - int k=FunctionDefaults<3>::get_k()) + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), + int k=FunctionDefaults<3>::get_k()) { - return SeparatedConvolution(world,OperatorInfo(0.0,lo,eps,OT_G12),bc,k); + return SeparatedConvolution(world,OperatorInfo(0.0,lo,eps,OT_G12),lattice_sum,k); } @@ -1768,10 +1753,10 @@ namespace madness { SeparatedConvolution* CoulombOperatorPtr(World& world, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), int k=FunctionDefaults<3>::get_k()) { - return new SeparatedConvolution(world,OperatorInfo(0.0,lo,eps,OT_G12),bc,k); + return new SeparatedConvolution(world,OperatorInfo(0.0,lo,eps,OT_G12),lattice_sum,k); } @@ -1780,13 +1765,13 @@ namespace madness { static inline SeparatedConvolution BSHOperator(World& world, double mu, double lo, double eps, - const BoundaryConditions& bc=FunctionDefaults::get_bc(), + const std::array& lattice_sum = FunctionDefaults::get_bc().is_periodic(), int k=FunctionDefaults::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(world,OperatorInfo(mu,lo,eps,OT_BSH),bc,k); + return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_BSH),lattice_sum,k); } /// Factory function generating separated kernel for convolution with BSH kernel in general NDIM @@ -1794,22 +1779,22 @@ namespace madness { static inline SeparatedConvolution* BSHOperatorPtr(World& world, double mu, double lo, double eps, - const BoundaryConditions& bc=FunctionDefaults::get_bc(), + const std::array& lattice_sum = FunctionDefaults::get_bc().is_periodic(), int k=FunctionDefaults::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(world,OperatorInfo(mu,lo,eps,OT_BSH),bc,k); + return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_BSH),lattice_sum,k); } /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D static inline SeparatedConvolution BSHOperator3D(World& world, double mu, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), int k=FunctionDefaults<3>::get_k()) { - return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_BSH),bc,k); + return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_BSH),lattice_sum,k); } /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D @@ -1820,22 +1805,28 @@ namespace madness { double mu, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), int k=FunctionDefaults<3>::get_k()) { - const Tensor& cell_width = FunctionDefaults<3>::get_cell_width(); - double hi = cell_width.normf(); // Diagonal width of cell - if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation - - GFit fit=GFit::BSHFit(mu,lo,hi,eps,false); - Tensor coeff=fit.coeffs(); - Tensor expnt=fit.exponents(); - - if (bc(0,0) == BC_PERIODIC) { - fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false); - } - return SeparatedConvolution(world, args, coeff, expnt, bc, k); + const Tensor &cell_width = FunctionDefaults<3>::get_cell_width(); + double hi = cell_width.normf(); // Diagonal width of cell + // Extend kernel range for lattice summation + const auto any_lattice_sum = std::accumulate( + lattice_sum.begin(), lattice_sum.end(), false, std::logical_or{}); + if (any_lattice_sum) { + hi *= 100; + } + + GFit fit = GFit::BSHFit(mu, lo, hi, eps, false); + Tensor coeff = fit.coeffs(); + Tensor expnt = fit.exponents(); + + if (any_lattice_sum) { + fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false); + } + return SeparatedConvolution(world, args, coeff, expnt, + lattice_sum, k); } /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D @@ -1845,52 +1836,59 @@ namespace madness { double mu, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), int k=FunctionDefaults<3>::get_k()) { - const Tensor& cell_width = FunctionDefaults<3>::get_cell_width(); - double hi = cell_width.normf(); // Diagonal width of cell - if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation - - GFit fit=GFit::BSHFit(mu,lo,hi,eps,false); - Tensor coeff=fit.coeffs(); - Tensor expnt=fit.exponents(); - - if (bc(0,0) == BC_PERIODIC) { - fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false); - } - return new SeparatedConvolution(world, args, coeff, expnt, bc, k); + const Tensor &cell_width = FunctionDefaults<3>::get_cell_width(); + double hi = cell_width.normf(); // Diagonal width of cell + // Extend kernel range for lattice summation + const auto any_lattice_sum = std::accumulate( + lattice_sum.begin(), lattice_sum.end(), false, std::logical_or{}); + if (any_lattice_sum) { + hi *= 100; + } + + GFit fit = GFit::BSHFit(mu, lo, hi, eps, false); + Tensor coeff = fit.coeffs(); + Tensor expnt = fit.exponents(); + + if (any_lattice_sum) { + fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false); + } + return new SeparatedConvolution(world, args, coeff, + expnt, lattice_sum, k); } static inline SeparatedConvolution SlaterF12Operator(World& world, double mu, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), int k=FunctionDefaults<3>::get_k()) { - return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F12),bc,k); + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), + int k=FunctionDefaults<3>::get_k()) { + return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F12),lattice_sum,k); } static inline SeparatedConvolution SlaterF12sqOperator(World& world, double mu, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), int k=FunctionDefaults<3>::get_k()) { - return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F212),bc,k); + return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F212),lattice_sum,k); } static inline SeparatedConvolution* SlaterF12sqOperatorPtr(World& world, double mu, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), int k=FunctionDefaults<3>::get_k()) { - return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F212),bc,k); + return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F212),lattice_sum,k); } /// Factory function generating separated kernel for convolution with exp(-mu*r) in 3D template static inline SeparatedConvolution SlaterOperator(World& world, double mu, double lo, double eps, - const BoundaryConditions& bc=FunctionDefaults::get_bc(), + const std::array& lattice_sum = FunctionDefaults::get_bc().is_periodic(), int k=FunctionDefaults::get_k()) { - return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_SLATER),bc,k); + return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_SLATER),lattice_sum,k); } /// Factory function generating separated kernel for convolution with exp(-mu*r*r) @@ -1899,9 +1897,9 @@ namespace madness { template static inline SeparatedConvolution GaussOperator(World& world, double mu, double lo=0.0, double eps=0.0, - const BoundaryConditions& bc=FunctionDefaults::get_bc(), + const std::array& lattice_sum = FunctionDefaults::get_bc().is_periodic(), int k=FunctionDefaults::get_k()) { - return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_GAUSS),bc,k); + return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_GAUSS),lattice_sum,k); } /// Factory function generating separated kernel for convolution with exp(-mu*r*r) in 3D @@ -1910,9 +1908,9 @@ namespace madness { template static inline SeparatedConvolution* GaussOperatorPtr(World& world, double mu, double lo=0.0, double eps=0.0, - const BoundaryConditions& bc = FunctionDefaults::get_bc(), + const std::array& lattice_sum = FunctionDefaults::get_bc().is_periodic(), int k = FunctionDefaults::get_k()) { - return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_GAUSS),bc,k); + return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_GAUSS),lattice_sum,k); } @@ -1921,18 +1919,18 @@ namespace madness { template static inline SeparatedConvolution* SlaterOperatorPtr_ND(World& world, double mu, double lo, double eps, - const BoundaryConditions& bc = FunctionDefaults::get_bc(), + const std::array& lattice_sum = FunctionDefaults::get_bc().is_periodic(), int k = FunctionDefaults::get_k()) { - return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_SLATER),bc,k); + return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_SLATER),lattice_sum,k); } /// Factory function generating separated kernel for convolution with exp(-mu*r) in 3D /// Note that the 1/(2mu) factor of SlaterF12Operator is not included, this is just the exponential function static inline SeparatedConvolution* SlaterOperatorPtr(World& world, double mu, double lo, double eps, - const BoundaryConditions<3>& bc = FunctionDefaults<3>::get_bc(), + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), int k = FunctionDefaults<3>::get_k()) { - return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_SLATER),bc,k); + return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_SLATER),lattice_sum,k); } /// Factory function generating separated kernel for convolution with (1 - exp(-mu*r))/(2 mu) in 3D @@ -1940,9 +1938,9 @@ namespace madness { /// includes the factor 1/(2 mu) static inline SeparatedConvolution* SlaterF12OperatorPtr(World& world, double mu, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), int k=FunctionDefaults<3>::get_k()) { - return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F12),bc,k); + return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F12),lattice_sum,k); } @@ -1952,9 +1950,9 @@ namespace madness { /// includes the factor 1/(2 mu) static inline SeparatedConvolution FGOperator(World& world, double mu, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), int k=FunctionDefaults<3>::get_k()) { - return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_FG12),bc,k); + return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_FG12),lattice_sum,k); } /// Factory function generating separated kernel for convolution with 1/(2 mu)*(1 - exp(-mu*r))/r in 3D @@ -1963,9 +1961,9 @@ namespace madness { /// includes the factor 1/(2 mu) static inline SeparatedConvolution* FGOperatorPtr(World& world, double mu, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), int k=FunctionDefaults<3>::get_k()) { - return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_FG12),bc,k); + return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_FG12),lattice_sum,k); } /// Factory function generating separated kernel for convolution with (1/(2 mu)*(1 - exp(-mu*r)))^2/r in 3D @@ -1975,9 +1973,9 @@ namespace madness { /// includes the factor 1/(2 mu)^2 static inline SeparatedConvolution* F2GOperatorPtr(World& world, double mu, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), int k=FunctionDefaults<3>::get_k()) { - return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F2G12),bc,k); + return new SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F2G12),lattice_sum,k); } /// Factory function generating separated kernel for convolution with (1/(2 mu)*(1 - exp(-mu*r)))^2/r in 3D @@ -1987,9 +1985,9 @@ namespace madness { /// includes the factor 1/(2 mu)^2 static inline SeparatedConvolution F2GOperator(World& world, double mu, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), int k=FunctionDefaults<3>::get_k()) { - return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F2G12),bc,k); + return SeparatedConvolution(world,OperatorInfo(mu,lo,eps,OT_F2G12),lattice_sum,k); } @@ -1997,15 +1995,14 @@ namespace madness { /// Gaussian (aka a widened delta function) static inline SeparatedConvolution SmoothingOperator3D(World& world, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), int k=FunctionDefaults<3>::get_k()) { double exponent = 1.0/(2.0*eps); Tensor coeffs(1), exponents(1); exponents(0L) = exponent; coeffs(0L)=pow(exponent/M_PI,0.5*3.0); // norm of the gaussian - return SeparatedConvolution(world, coeffs, exponents, 1.e-8, eps); - + return SeparatedConvolution(world, coeffs, exponents, 1.e-8, eps, lattice_sum, k); } /// Factory function generating separated kernel for convolution a normalized @@ -2013,14 +2010,14 @@ namespace madness { template static inline SeparatedConvolution SmoothingOperator(World& world, double eps, - const BoundaryConditions& bc=FunctionDefaults::get_bc(), + const std::array& lattice_sum = FunctionDefaults::get_bc().is_periodic(), int k=FunctionDefaults::get_k()) { double exponent = 1.0/(2.0*eps); Tensor coeffs(1), exponents(1); exponents(0L) = exponent; coeffs(0L)=pow(exponent/M_PI,0.5*NDIM); // norm of the gaussian - return SeparatedConvolution(world, coeffs, exponents, 1.e-8, eps); + return SeparatedConvolution(world, coeffs, exponents, 1.e-8, eps, lattice_sum, k); } /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D @@ -2030,21 +2027,26 @@ namespace madness { double mu, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), - int k=FunctionDefaults<3>::get_k()) - { - const Tensor& cell_width = FunctionDefaults<3>::get_cell_width(); - double hi = cell_width.normf(); // Diagonal width of cell - if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation - - GFit fit=GFit::BSHFit(mu,lo,hi,eps,false); - Tensor coeff=fit.coeffs(); - Tensor expnt=fit.exponents(); - - if (bc(0,0) == BC_PERIODIC) { - fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false); - } - return new SeparatedConvolution(world, coeff, expnt, lo, eps, bc, k); + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), + int k=FunctionDefaults<3>::get_k()) { + const Tensor &cell_width = FunctionDefaults<3>::get_cell_width(); + double hi = cell_width.normf(); // Diagonal width of cell + // Extend kernel range for lattice summation + const auto any_lattice_sum = std::accumulate( + lattice_sum.begin(), lattice_sum.end(), false, std::logical_or{}); + if (any_lattice_sum) { + hi *= 100; + } + + GFit fit = GFit::BSHFit(mu, lo, hi, eps, false); + Tensor coeff = fit.coeffs(); + Tensor expnt = fit.exponents(); + + if (any_lattice_sum) { + fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false); + } + return new SeparatedConvolution(world, coeff, expnt, lo, eps, + lattice_sum, k); } @@ -2058,47 +2060,54 @@ namespace madness { GradCoulombOperator(World& world, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), - int k=FunctionDefaults<3>::get_k()) - { - typedef SeparatedConvolution real_convolution_3d; - typedef std::shared_ptr real_convolution_3d_ptr; - const double pi = constants::pi; - const Tensor width = FunctionDefaults<3>::get_cell_width(); - double hi = width.normf(); // Diagonal width of cell - const bool isperiodicsum = (bc(0,0)==BC_PERIODIC); - if (isperiodicsum) hi *= 100; // Extend range for periodic summation - - GFit fit=GFit::CoulombFit(lo,hi,eps,false); - Tensor coeff=fit.coeffs(); - Tensor expnt=fit.exponents(); - - if (bc(0,0) == BC_PERIODIC) { - fit.truncate_periodic_expansion(coeff, expnt, width.max(), true); + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), + int k=FunctionDefaults<3>::get_k()) { + typedef SeparatedConvolution real_convolution_3d; + typedef std::shared_ptr real_convolution_3d_ptr; + const double pi = constants::pi; + const Tensor width = FunctionDefaults<3>::get_cell_width(); + double hi = width.normf(); // Diagonal width of cell + // Extend kernel range for lattice summation + const auto any_lattice_sum = std::accumulate( + lattice_sum.begin(), lattice_sum.end(), false, std::logical_or{}); + if (any_lattice_sum) { + hi *= 100; + } + + GFit fit = GFit::CoulombFit(lo, hi, eps, false); + Tensor coeff = fit.coeffs(); + Tensor expnt = fit.exponents(); + + if (any_lattice_sum) { + fit.truncate_periodic_expansion(coeff, expnt, width.max(), true); + } + + int rank = coeff.dim(0); + + std::vector gradG(3); + + for (int dir = 0; dir < 3; dir++) { + std::vector> 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::get( + k, expnt(mu) * width[d] * width[d], 0, + lattice_sum[d])); + } + ops[mu].setop(dir, GaussianConvolution1DCache::get( + k, expnt(mu) * width[dir] * width[dir], 1, + lattice_sum[dir])); } + gradG[dir] = real_convolution_3d_ptr( + new SeparatedConvolution(world, ops)); + } - int rank = coeff.dim(0); - - std::vector gradG(3); - - for (int dir=0; dir<3; dir++) { - std::vector< ConvolutionND > ops(rank); - for (int mu=0; mu::get(k, expnt(mu)*width[d]*width[d], 0, isperiodicsum)); - } - ops[mu].setop(dir,GaussianConvolution1DCache::get(k, expnt(mu)*width[dir]*width[dir], 1, isperiodicsum)); - } - gradG[dir] = real_convolution_3d_ptr(new SeparatedConvolution(world, ops)); - } - - return gradG; + return gradG; } /// Factory function generating operator for convolution with grad(bsh) in 3D @@ -2112,47 +2121,54 @@ namespace madness { double mu, double lo, double eps, - const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), - int k=FunctionDefaults<3>::get_k()) - { - typedef SeparatedConvolution real_convolution_3d; - typedef std::shared_ptr real_convolution_3d_ptr; - const double pi = constants::pi; - const Tensor width = FunctionDefaults<3>::get_cell_width(); - double hi = width.normf(); // Diagonal width of cell - const bool isperiodicsum = (bc(0,0)==BC_PERIODIC); - if (isperiodicsum) hi *= 100; // Extend range for periodic summation - - GFit fit=GFit::BSHFit(mu,lo,hi,eps,false); - Tensor coeff=fit.coeffs(); - Tensor expnt=fit.exponents(); - - if (bc(0,0) == BC_PERIODIC) { - fit.truncate_periodic_expansion(coeff, expnt, width.max(), true); - } - - int rank = coeff.dim(0); - - std::vector gradG(3); - - for (int dir=0; dir<3; dir++) { - std::vector< ConvolutionND > ops(rank); - for (int mu=0; mu::get(k, expnt(mu)*width[d]*width[d], 0, isperiodicsum)); - } - ops[mu].setop(dir,GaussianConvolution1DCache::get(k, expnt(mu)*width[dir]*width[dir], 1, isperiodicsum)); - } - gradG[dir] = real_convolution_3d_ptr(new SeparatedConvolution(world, ops)); + const std::array& lattice_sum = FunctionDefaults<3>::get_bc().is_periodic(), + int k=FunctionDefaults<3>::get_k()) { + typedef SeparatedConvolution real_convolution_3d; + typedef std::shared_ptr real_convolution_3d_ptr; + const double pi = constants::pi; + const Tensor width = FunctionDefaults<3>::get_cell_width(); + double hi = width.normf(); // Diagonal width of cell + // Extend kernel range for lattice summation + const auto any_lattice_sum = std::accumulate( + lattice_sum.begin(), lattice_sum.end(), false, std::logical_or{}); + if (any_lattice_sum) { + hi *= 100; + } + + GFit fit = GFit::BSHFit(mu, lo, hi, eps, false); + Tensor coeff = fit.coeffs(); + Tensor expnt = fit.exponents(); + + if (any_lattice_sum) { + fit.truncate_periodic_expansion(coeff, expnt, width.max(), true); + } + + int rank = coeff.dim(0); + + std::vector gradG(3); + + for (int dir = 0; dir < 3; dir++) { + std::vector> 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::get( + k, expnt(mu) * width[d] * width[d], 0, + lattice_sum[d])); + } + ops[mu].setop(dir, GaussianConvolution1DCache::get( + k, expnt(mu) * width[dir] * width[dir], 1, + lattice_sum[dir])); } + gradG[dir] = real_convolution_3d_ptr( + new SeparatedConvolution(world, ops)); + } - return gradG; + return gradG; } diff --git a/src/madness/mra/qmprop.cc b/src/madness/mra/qmprop.cc index 7108480f906..f4a61536b4a 100644 --- a/src/madness/mra/qmprop.cc +++ b/src/madness/mra/qmprop.cc @@ -172,7 +172,7 @@ namespace madness { double width = FunctionDefaults::get_cell_min_width(); // Assuming cubic so all dim equal std::vector< std::shared_ptr< Convolution1D > > q(1); q[0].reset(qm_1d_free_particle_propagator(k, bandlimit, timestep, width)); - return SeparatedConvolution(world, q, BoundaryConditions(BC_FREE), k, true); + return SeparatedConvolution(world, q, no_lattice_sum(), k, true); } template @@ -181,7 +181,7 @@ namespace madness { double width = FunctionDefaults::get_cell_min_width(); // Assuming cubic so all dim equal std::vector< std::shared_ptr< Convolution1D > > q(1); q[0].reset(qm_1d_free_particle_propagator(k, bandlimit, timestep, width)); - return new SeparatedConvolution(world, q, BoundaryConditions(BC_FREE), k, true); + return new SeparatedConvolution(world, q, no_lattice_sum(), k, true); } #ifdef FUNCTION_INSTANTIATE_1