MADNESS  0.10.1
function_interface.h
Go to the documentation of this file.
1 /*
2  This file is part of MADNESS.
3 
4  Copyright (C) 2007,2010 Oak Ridge National Laboratory
5 
6  This program is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2 of the License, or
9  (at your option) any later version.
10 
11  This program is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with this program; if not, write to the Free Software
18  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20  For more information please contact:
21 
22  Robert J. Harrison
23  Oak Ridge National Laboratory
24  One Bethel Valley Road
25  P.O. Box 2008, MS-6367
26 
27  email: harrisonrj@ornl.gov
28  tel: 865-241-3937
29  fax: 865-572-0680
30 
31  $Id: function_factory_and_interface.h 3422 2014-03-24 09:16:15Z 3ru6ruWu $
32 */
33 
34 
35 #ifndef MADNESS_MRA_FUNCTION_INTERFACE_H__INCLUDED
36 #define MADNESS_MRA_FUNCTION_INTERFACE_H__INCLUDED
37 
38 #include <madness/tensor/tensor.h>
40 #include <madness/mra/key.h>
41 
42 // needed for the TwoElectronInterface
44 #include <madness/mra/gfit.h>
47 
48 namespace madness {
49 
50  // forward declaration needed for CompositeFunctorInterface
51  template<typename T, std::size_t NDIM>
52  class FunctionImpl;
53 
54  template<typename T, std::size_t NDIM>
55  class Function;
56 
57  template<typename T, std::size_t NDIM>
58  Tensor<T> fcube(const Key<NDIM>&, T (*f)(const Vector<double,NDIM>&), const Tensor<double>&);
59 
60 // template <typename T, std::size_t NDIM>
61 // const std::vector<Function<T,NDIM>>& change_tree_state(const std::vector<Function<T,NDIM>>& v,
62 // const TreeState finalstate,
63 // const bool fence);
64 
65 
66  /// Abstract base class interface required for functors used as input to Functions
67  template<typename T, std::size_t NDIM>
69  public:
70 
72  typedef Key<NDIM> keyT;
73  typedef T value_type;
74 
76 
78 
79  /// adapt the special level to resolve the smallest length scale
80  void set_length_scale(double lo) {
82  double lo_sim=lo/Lmax; // lo in simulation coordinates;
83  special_level_=Level(-log2(lo_sim));
84  }
85 
86  /// Can we screen this function based on the bounding box information?
87  virtual bool screened(const Vector<double,NDIM>& c1, const Vector<double,NDIM>& c2) const {
88  return false;
89  }
90 
91  /// Does the interface support a vectorized operator()?
92  virtual bool supports_vectorized() const {return false;}
93 
94  virtual void operator()(const Vector<double*,1>& xvals, T* fvals, int npts) const {
95  MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0);
96  }
97 
98  virtual void operator()(const Vector<double*,2>& xvals, T* fvals, int npts) const {
99  MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0);
100  }
101 
102  virtual void operator()(const Vector<double*,3>& xvals, T* fvals, int npts) const {
103  MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0);
104  }
105 
106  virtual void operator()(const Vector<double*,4>& xvals, T* fvals, int npts) const {
107  MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0);
108  }
109 
110  virtual void operator()(const Vector<double*,5>& xvals, T* fvals, int npts) const {
111  MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0);
112  }
113 
114  virtual void operator()(const Vector<double*,6>& xvals, T* fvals, int npts) const {
115  MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0);
116  }
117 
118  /// You should implement this to return \c f(x)
119  virtual T operator()(const Vector<double, NDIM>& x) const = 0;
120 
121  /// Override this to return list of special points to be refined more deeply
122  virtual std::vector< Vector<double,NDIM> > special_points() const {
123  return std::vector< Vector<double,NDIM> >();
124  }
125 
126  /// Override this change level refinement for special points (default is 6)
127  virtual Level special_level() {return special_level_;}
128 
130 
131  virtual coeffT coeff(const keyT&) const {
132  MADNESS_EXCEPTION("implement coeff for FunctionFunctorInterface",0);
133  return coeffT();
134  }
135 
136  virtual coeffT values(const keyT& key, const Tensor<double>& tensor) const {
137  MADNESS_EXCEPTION("implement values for FunctionFunctorInterface",0);
138  return coeffT();
139  }
140 
141  /// does this functor directly provide sum coefficients? or only function values?
142  virtual bool provides_coeff() const {
143  return false;
144  }
145 
146  };
147 
148 
149 
150  ///forward declaration
151  template <typename T, std::size_t NDIM>
152  // void FunctionImpl<T,NDIM>::fcube(const keyT& key, const FunctionFunctorInterface<T,NDIM>& f, const Tensor<double>& qx, tensorT& fval) const {
153  void fcube(const Key<NDIM>& key, const FunctionFunctorInterface<T,NDIM>& f, const Tensor<double>& qx, Tensor<T>& fval);
154 
155  /// CompositeFunctorInterface implements a wrapper of holding several functions and functors
156 
157  /// Use this to "connect" several functions and/or functors and to return their coefficients
158  /// e.g. connect f1 and f2 with an addition, you can request the coefficients of any node
159  /// and they will be computed on the fly and returned. Mainly useful to connect a functor
160  /// with a function, if the functor is too large to be represented in MRA (e.g. 1/r12)
161  ///
162  /// as of now, the operation connecting the functions/functors is simply addition.
163  /// need to implement expression templates, if I only knew what that was...
164  template<typename T, std::size_t NDIM, std::size_t MDIM>
166 
167  typedef Vector<double, NDIM> coordT; ///< Type of vector holding coordinates
170  typedef std::shared_ptr<implT> pimplT;
171  typedef std::shared_ptr<implL> pimplL;
172 
174 
175  public:
176  /// various MRA functions of NDIM dimensionality
177  std::vector<std::shared_ptr<implT>> impl_ket_vector; ///< supposedly the pair function
178  std::shared_ptr<implT> impl_eri; ///< supposedly 1/r12
179 
180  /// various MRA functions of MDIM dimensionality (e.g. 3, if NDIM==6)
181  std::shared_ptr<implL> impl_m1; ///< supposedly 1/r1
182  std::shared_ptr<implL> impl_m2; ///< supposedly 1/r2
183  std::vector<std::shared_ptr<implL>> impl_p1_vector; ///< supposedly orbital 1
184  std::vector<std::shared_ptr<implL>> impl_p2_vector; ///< supposedly orbital 2
185 
186  public:
187 
188  /// constructor takes its Factory
189  CompositeFunctorInterface(World& world, std::vector<pimplT> ket, pimplT g12,
190  pimplL v1, pimplL v2, std::vector<pimplL> p1, std::vector<pimplL> p2)
191  : world(world), impl_ket_vector(ket), impl_eri(g12)
192  , impl_m1(v1), impl_m2(v2), impl_p1_vector(p1), impl_p2_vector(p2)
193  {
194 
195  // some consistency checks
196  // either a pair ket is provided, or two particles (tba)
197  MADNESS_CHECK_THROW(impl_p1_vector.size()==impl_p2_vector.size(), "CompositeFunctorInterface: p1 and p2 must have the same size");
198  MADNESS_CHECK_THROW(impl_ket_vector.size()>0 or (impl_p1_vector.size()>0),"CompositeFunctorInterface: either ket or p1 must be provided");
199 
200  }
201 
202  /// replicate low-dimensional functions over all ranks of this world
203  void replicate_low_dim_functions(const bool fence) {
204  if (impl_m1 and (not impl_m1->is_on_demand())) impl_m1->replicate(false);
205  if (impl_m2 and (not impl_m2->is_on_demand())) impl_m2->replicate(false);
206 
207  for (auto& p1 : impl_p1_vector) if (p1 and (not p1->is_on_demand())) p1->replicate(false);
208  for (auto& p2 : impl_p2_vector) if (p2 and (not p2->is_on_demand())) p2->replicate(false);
209  if (fence) world.gop.fence();
210  }
211 
212  void make_redundant(const bool fence) {
213  // prepare base functions that make this function
214  for (auto& k : impl_ket_vector) if (k and (not k->is_on_demand())) k->change_tree_state(redundant,false);
215  if (impl_eri) {
216  if (not impl_eri->is_on_demand()) impl_eri->change_tree_state(redundant,false);
217  }
218  if (impl_m1 and (not impl_m1->is_on_demand())) impl_m1->change_tree_state(redundant,false);
219  if (impl_m2 and (not impl_m2->is_on_demand())) impl_m2->change_tree_state(redundant,false);
220 
223 // for (auto& k : impl_p1_vector) if (k and (not k->is_on_demand())) k->change_tree_state(redundant,false);
224 // for (auto& k : impl_p2_vector) if (k and (not k->is_on_demand())) k->change_tree_state(redundant,false);
225 // if (impl_p1 and (not impl_p1->is_on_demand())) impl_p1->change_tree_state(redundant,false);
226 // if (impl_p2 and (not impl_p2->is_on_demand())) impl_p2->change_tree_state(redundant,false);
227  if (fence) world.gop.fence();
228  }
229 
230  /// return true if all constituent functions are in redundant tree state
231  bool check_redundant() const {
232  for (auto& k : impl_ket_vector) if (k and (not k->is_redundant())) return false;
233 // if (impl_ket and (not impl_ket->is_redundant())) return false;
234  if (impl_eri) MADNESS_ASSERT(impl_eri->is_on_demand());
235  if (impl_m1 and (not impl_m1->is_redundant())) return false;
236  if (impl_m2 and (not impl_m2->is_redundant())) return false;
237  for (auto& k : impl_p1_vector) if (k and (not k->is_redundant())) return false;
238  for (auto& k : impl_p2_vector) if (k and (not k->is_redundant())) return false;
239 // if (impl_p1 and (not impl_p1->is_redundant())) return false;
240 // if (impl_p2 and (not impl_p2->is_redundant())) return false;
241  return true;
242  }
243 
244  /// return value at point x; fairly inefficient
245  T operator()(const coordT& x) const {
246  print("there is no operator()(coordT&) in CompositeFunctorInterface, for good reason");
247  MADNESS_ASSERT(0);
248  return T(0);
249  };
250 
251  bool provides_coeff() const {
252  return false;
253  }
254 
255  };
256 
257 
258  /// ElementaryInterface (formerly FunctorInterfaceWrapper) interfaces a c-function
259 
260  /// hard-code your favorite function and interface it with this; Does only
261  /// provide function values, no MRA coefficients. Care must be taken if the
262  /// function we refer to is a singular function, and a on-demand function
263  /// at the same time, since direct computation of coefficients via mraimpl::project
264  /// might suffer from inaccurate quadrature.
265  template<typename T, std::size_t NDIM>
267 
268  public:
269  typedef Vector<double, NDIM> coordT; ///< Type of vector holding coordinates
271 
272  T (*f)(const coordT&);
273 
274  ElementaryInterface(T (*f)(const coordT&)) : f(f) {}
275 
276  T operator()(const coordT& x) const {return f(x);}
277 
278  coeffT values(const Key<NDIM>& key, const Tensor<double>& quad_x) const {
279  typedef Tensor<T> tensorT;
280  tensorT fval=madness::fcube(key,f,quad_x);
282  }
283  };
284 
285  /// FunctorInterface interfaces a class or struct with an operator()()
286  template<typename T, std::size_t NDIM, typename opT>
288 
289  public:
290  typedef Vector<double, NDIM> coordT; ///< Type of vector holding coordinates
292 
293  opT op;
294 
295  FunctorInterface(const opT& op) : op(op) {}
296 
297  T operator()(const coordT& x) const {return op(x);}
298  };
299 
300  /// FunctionInterface implements a wrapper around any class with the operator()()
301  template<typename T, size_t NDIM, typename opT>
303 
305  typedef Vector<double, NDIM> coordT; ///< Type of vector holding coordinates
306 
307  const opT op;
308 
309  public:
310  FunctionInterface(const opT& op) : op(op) {}
311 
312  T operator()(const coordT& coord) const {return op(coord);}
313 
314  bool provides_coeff() const {return false;}
315 
316  };
317 
318  /// base class to compute the wavelet coefficients for an isotropic 2e-operator
319 
320  /// all classes that derive from this base class use the Gaussian fitting
321  /// procedure that has been developed for the BSH operator. We simply
322  /// reuse the wavelet coefficients that we get from there to avoid
323  /// evaluating the functions themselves, since the quadrature of singular
324  /// functions is imprecise and slow.
325  template<typename T, std::size_t NDIM>
327  protected:
328  static constexpr std::size_t LDIM=NDIM/2;
329  public:
330 
332 
333  /// constructor: cf the Coulomb kernel
334 
335  /// @param[in] lo the smallest length scale to be resolved
336  /// @param[in] eps the accuracy threshold
337  TwoElectronInterface(double lo, double eps,
340  :rank(), k(kk), lo(lo), hi(1.0) {
341 
342  // Presently we must have periodic or non-periodic in all dimensions.
343  for (std::size_t d=1; d<NDIM; ++d) {MADNESS_ASSERT(bc(d,0)==bc(0,0));}
344 
346  hi = width.normf(); // Diagonal width of cell
347  if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
348 
349  }
350 
351  bool provides_coeff() const {
352  return true;
353  }
354 
355  /// return the coefficients of the function in 6D (x1,y1,z1, x2,y2,z2)
356  coeffT coeff(const Key<NDIM>& key) const {
359  }
360 
362  print("there is no operator()(coordT&) in TwoElectronInterface, for good reason");
363  MADNESS_ASSERT(0);
364  return T(0);
365  }
366 
367  protected:
368 
369  /// make the coefficients from the 1d convolution
370  Tensor<double> make_coeff(const Key<NDIM>& key) const {
371  const Level n=key.level();
373 
374  // get the displacements for all 3 dimensions: x12, y12, z12
375  Translation l0, l1, l2;
376  if (NDIM==2) {
377  l0=l[0]-l[1];
378  } else if (NDIM==4) {
379  l0=(l[0]-l[2]);
380  l1=(l[1]-l[3]);
381  } else if (NDIM==6) {
382  l0=(l[0]-l[3]);
383  l1=(l[1]-l[4]);
384  l2=(l[2]-l[5]);
385  } else {
386  MADNESS_EXCEPTION("TwoElectronInterface: NDIM must be 2, 4, or 6",1);
387  }
388 
389  Tensor<double> scr1(rank,k*k), scr2(rank,k*k,k*k), scr3(rank,k*k);
390 
391  // lump all the terms together
392  for (long mu=0; mu<rank; mu++) {
393  Tensor<double> r0, r1, r2;
394  if (NDIM>=2) r0=(ops[mu].getop(0)->rnlij(n,l0)).reshape(k*k);
395  if (NDIM>=4) r1=(ops[mu].getop(1)->rnlij(n,l1)).reshape(k*k);
396  if (NDIM>=6) r2=(ops[mu].getop(2)->rnlij(n,l2)).reshape(k*k);
397 
398  // include weights in first vector
399  scr1(mu,Slice(_))=r0*ops[mu].getfac();
400 
401  if (NDIM==2) {
402  ;
403  } else if (NDIM==4) {
404  scr3(mu,Slice(_))=r1;
405  } else if (NDIM==6) {
406  // merge second and third vector to scr(r,k1,k2)
407  scr2(mu,Slice(_),Slice(_))=outer(r1,r2);
408  } else {
409  MADNESS_EXCEPTION("TwoElectronInterface: NDIM must be 2, 4, or 6",1);
410  }
411  }
412 
413  if (NDIM==2) {
414  // perform sum over the rank
415  Tensor<double> result(scr1.dim(1));
416  for (long mu=0; mu<rank; ++mu) result(_)+= scr1(mu,_);
417  return result;
418  }
419  else if (NDIM==4) return inner(scr1,scr3,0,0);
420  else if (NDIM==6) return inner(scr1,scr2,0,0);
421  else {
422  MADNESS_EXCEPTION("TwoElectronInterface: NDIM must be 2, 4, or 6",1);
423  return Tensor<double>();
424  }
425  }
426 
427  /// the dimensions are a bit confused (x1,x2, y1,y2, z1,z2) -> (x1,y1,z1, x2,y2,z2)
429  std::vector<long> map(NDIM);
430  if (NDIM==2) {
431  map[0]=0; map[1]=1;
432  return copy(c.reshape(k,k));
433  } else if (NDIM==4) {
434  map[0]=0; map[1]=2;
435  map[2]=1; map[3]=3;
436  return copy(c.reshape(k,k,k,k).mapdim(map));
437  } else if (NDIM==6) {
438  map[0]=0; map[1]=3; map[2]=1;
439  map[3]=4; map[4]=2; map[5]=5;
440  return copy(c.reshape(k,k,k,k,k,k).mapdim(map));
441  }
442  return Tensor<double>();
443  }
444 
445  /// initialize the Gaussian fit; uses the virtual function fit() to fit
446  void initialize(const double eps) {
447  GFit<double,LDIM> fit=this->fit(eps);
448  Tensor<double> coeff=fit.coeffs();
449  Tensor<double> expnt=fit.exponents();
450 
451  // set some parameters
452  rank=coeff.dim(0);
453  ops.resize(rank);
455 
456  // construct all the terms
457  for (int mu=0; mu<rank; ++mu) {
458  // double c = std::pow(sqrt(expnt(mu)/pi),static_cast<int>(NDIM)); // Normalization coeff
459  double c = std::pow(sqrt(expnt(mu)/constants::pi),LDIM); // Normalization coeff
460 
461  // We cache the normalized operator so the factor is the value we must multiply
462  // by to recover the coeff we want.
463  ops[mu].setfac(coeff(mu)/c);
464 
465  // only 3 dimensions here!
466  for (std::size_t d=0; d<LDIM; ++d) {
467  ops[mu].setop(d,GaussianConvolution1DCache<double>::get(k, expnt(mu)*width[d]*width[d], 0, false));
468  }
469  }
470  }
471 
472  /// derived classes must implement this -- cf GFit.h
473  virtual GFit<double,LDIM> fit(const double eps) const = 0;
474 
475  /// storing the coefficients
476  mutable std::vector< ConvolutionND<double,NDIM> > ops;
477 
478  /// the number of terms in the Gaussian quadrature
479  int rank;
480 
481  /// the wavelet order
482  int k;
483 
484  /// the smallest length scale that needs to be represented
485  double lo;
486 
487  /// the largest length scale that needs to be represented
488  double hi;
489 
490  };
491 
492 
493  /// a function like f(x)=1/x
494  template<typename T, std::size_t NDIM>
496  public:
497 
498  /// constructor: cf the Coulomb kernel
499 
500  /// @param[in] lo the smallest length scale to be resolved
501  /// @param[in] eps the accuracy threshold
506 
507  if (info.hi<0) {
509  if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
510  this->info.hi=hi;
511  }
512  this->initialize(info.thresh);
513  }
514 
515  private:
517  static constexpr std::size_t LDIM=NDIM/2;
518 
519  GFit<double,LDIM> fit(const double eps) const {
520  return GFit<double,LDIM>(info);
521  }
522  };
523 
524  /// a function like f(x)=1/x
525  template<typename T, std::size_t NDIM>
526  class ElectronRepulsionInterface : public TwoElectronInterface<double,NDIM> {
527 
528  public:
529 
530  /// constructor: cf the Coulomb kernel
531 
532  /// @param[in] lo the smallest length scale to be resolved
533  /// @param[in] eps the accuracy threshold
534  ElectronRepulsionInterface(double lo,double eps,
537  : TwoElectronInterface<double,NDIM>(lo,eps,bc,kk) {
538 
539  this->initialize(eps);
540  }
541 
542  private:
543  static constexpr std::size_t LDIM=NDIM/2;
544 
545  GFit<double,LDIM> fit(const double eps) const {
546  return GFit<double,LDIM>::CoulombFit(this->lo,this->hi,eps,false);
547  }
548  };
549 
550  /// a function like f(x) = exp(-mu x)/x
551  class BSHFunctionInterface : public TwoElectronInterface<double,6> {
552  public:
553 
554  /// constructor: cf the Coulomb kernel
555 
556  /// @param[in] mu the exponent of the BSH/inverse Laplacian
557  /// @param[in] lo the smallest length scale to be resolved
558  /// @param[in] eps the accuracy threshold
559  BSHFunctionInterface(double mu, double lo, double eps,
562  : TwoElectronInterface<double,6>(lo,eps,bc,kk), mu(mu) {
563 
564  initialize(eps);
565  }
566 
567  private:
568 
569  double mu;
570 
571  GFit<double,3> fit(const double eps) const {
572  return GFit<double,3>::BSHFit(mu,lo,hi,eps,false);
573  }
574  };
575 
576  /// a function like f(x)=exp(-mu x)
578  public:
579 
580  /// constructor: cf the Coulomb kernel
581 
582  /// @param[in] mu the exponent of the Slater function
583  /// @param[in] lo the smallest length scale to be resolved
584  /// @param[in] eps the accuracy threshold
585  SlaterFunctionInterface(double mu, double lo, double eps,
588  : TwoElectronInterface<double,6>(lo,eps,bc,kk), mu(mu) {
589  initialize(eps);
590  }
591 
592  private:
593 
594  double mu;
595 
596  GFit<double,3> fit(const double eps) const {
597  return GFit<double,3>::SlaterFit(mu,lo,hi,eps,false);
598  }
599  };
600 
601  /// a function like f(x) = (1 - exp(-mu x))/(2 gamma)
602  class SlaterF12Interface : public TwoElectronInterface<double,6> {
603  public:
604 
605  /// constructor: cf the Coulomb kernel
606 
607  /// @param[in] mu the exponent of the Slater function
608  /// @param[in] lo the smallest length scale to be resolved
609  /// @param[in] eps the accuracy threshold
610  SlaterF12Interface(double mu, double lo, double eps,
613  : TwoElectronInterface<double,6>(lo,eps,bc,kk), mu(mu) {
614 
615  initialize(eps);
616  }
617 
618 // /// overload the function of the base class
619 // coeffT coeff(const Key<6>& key) const {
620 //
621 // Tensor<double> c=make_coeff(key);
622 //
623 // // subtract 1 from the (0,0,..,0) element of the tensor,
624 // // which is the 0th order polynomial coefficient
625 // double one_coeff1=1.0*sqrt(FunctionDefaults<6>::get_cell_volume())
626 // *pow(0.5,0.5*6*key.level());
627 // std::vector<long> v0(6,0L);
628 // c(v0)-=one_coeff1;
629 //
630 // c.scale(-0.5/mu);
631 // return coeffT(map_coeff(c),FunctionDefaults<6>::get_thresh(),TT_FULL);
632 // }
633 
634  private:
635 
636  double mu;
637 
638  GFit<double,3> fit(const double eps) const {
639  return GFit<double,3>::F12Fit(mu,lo,hi,eps,false);
640  }
641  };
642 
643 // Not right
644 // /// a function like f(x) = (1 - exp(-mu x))/x
645 // class FGInterface : public TwoElectronInterface<double,6> {
646 // public:
647 //
648 // /// constructor: cf the Coulomb kernel
649 //
650 // /// @param[in] mu the exponent of the Slater function
651 // /// @param[in] lo the smallest length scale to be resolved
652 // /// @param[in] eps the accuracy threshold
653 // FGInterface(double mu, double lo, double eps,
654 // const BoundaryConditions<6>& bc=FunctionDefaults<6>::get_bc(),
655 // int kk=FunctionDefaults<6>::get_k())
656 // : TwoElectronInterface<double,6>(lo,eps,bc,kk), mu(mu) {
657 //
658 // initialize(eps);
659 // }
660 //
661 // private:
662 //
663 // double mu;
664 //
665 // GFit<double,3> fit(const double eps) const {
666 // return GFit<double,3>::SlaterFit(mu,lo,hi,eps,false);
667 // }
668 // };
669 
670 
671 #if 0
672 
673  /// ElectronRepulsionInterface implements the electron repulsion term 1/r12
674 
675  /// this is essentially just a wrapper around ElectronRepulsion
676  template<typename T, std::size_t NDIM>
677  class ElectronRepulsionInterface : public FunctionFunctorInterface<T,NDIM> {
678 
679  typedef GenTensor<T> coeffT;
680  typedef Vector<double, NDIM> coordT; ///< Type of vector holding coordinates
681 
682  /// the class computing the coefficients
683  ElectronRepulsion eri;
684 
685  public:
686 
687  /// constructor takes the same parameters as the Coulomb operator
688  /// which it uses to compute the coefficients
689  ElectronRepulsionInterface(World& world,double lo,double eps,
690  const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(),
692  : eri(ElectronRepulsion(eps,eps,bc,k)) {
693  }
694 
695 
696  /// return value at point x; fairly inefficient
697  T operator()(const coordT& x) const {
698  print("there is no operator()(coordT&) in ElectronRepulsionInterface, for good reason");
699  MADNESS_ASSERT(0);
700  return T(0);
701  };
702 
703 
704  /// return sum coefficients for imagined node at key
705  coeffT coeff(const Key<NDIM>& key) const {
706  return coeffT(this->eri.coeff(key),FunctionDefaults<NDIM>::get_thresh(),
707  TT_FULL);
708  }
709 
710  };
711 
712  /// FGIntegralInterface implements the two-electron integral (1-exp(-gamma*r12))/r12
713 
714  /// this is essentially just a wrapper around ElectronRepulsion
715  /// The integral expressed as: 1/r12 - exp(-gamma*r12)/r12
716  /// which can be expressed with an eri and a bsh
717  template<typename T, std::size_t NDIM>
718  class FGIntegralInterface : public FunctionFunctorInterface<T,NDIM> {
719 
720  typedef GenTensor<T> coeffT;
721  typedef Vector<double, NDIM> coordT; ///< Type of vector holding coordinates
722 
723  /// the class computing the coefficients
724  ElectronRepulsion eri;
725  BSHFunction bsh;
726 
727  public:
728 
729  /// constructor takes the same parameters as the Coulomb operator
730  /// which it uses to compute the coefficients
731  FGIntegralInterface(World& world, double lo, double eps, double gamma,
732  const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(),
734  : eri(ElectronRepulsion(eps,eps,0.0,bc,k))
735  , bsh(BSHFunction(eps,eps,gamma,bc,k)) {
736  }
737 
738  bool provides_coeff() const {
739  return true;
740  }
741 
742  /// return value at point x; fairly inefficient
743  T operator()(const coordT& x) const {
744  print("there is no operator()(coordT&) in FGIntegralInterface, for good reason");
745  MADNESS_ASSERT(0);
746  return T(0);
747  };
748 
749  /// return sum coefficients for imagined node at key
750  coeffT coeff(const Key<NDIM>& key) const {
751  typedef Tensor<T> tensorT;
752  tensorT e_b=eri.coeff(key)-bsh.coeff(key);
753  return coeffT(e_b,FunctionDefaults<NDIM>::get_thresh(),TT_FULL);
754  }
755 
756  };
757 
758 #endif
759 
760 }
761 
762 #endif // MADNESS_MRA_FUNCTION_INTERFACE_H__INCLUDED
a function like f(x) = exp(-mu x)/x
Definition: function_interface.h:551
GFit< double, 3 > fit(const double eps) const
derived classes must implement this – cf GFit.h
Definition: function_interface.h:571
double mu
Definition: function_interface.h:569
BSHFunctionInterface(double mu, double lo, double eps, const BoundaryConditions< 6 > &bc=FunctionDefaults< 6 >::get_bc(), int kk=FunctionDefaults< 6 >::get_k())
constructor: cf the Coulomb kernel
Definition: function_interface.h:559
This class is used to specify boundary conditions for all operators.
Definition: funcdefaults.h:101
CompositeFunctorInterface implements a wrapper of holding several functions and functors.
Definition: function_interface.h:165
FunctionImpl< T, NDIM > implT
Definition: function_interface.h:168
CompositeFunctorInterface(World &world, std::vector< pimplT > ket, pimplT g12, pimplL v1, pimplL v2, std::vector< pimplL > p1, std::vector< pimplL > p2)
constructor takes its Factory
Definition: function_interface.h:189
std::vector< std::shared_ptr< implT > > impl_ket_vector
various MRA functions of NDIM dimensionality
Definition: function_interface.h:177
std::shared_ptr< implT > pimplT
Definition: function_interface.h:170
bool provides_coeff() const
does this functor directly provide sum coefficients? or only function values?
Definition: function_interface.h:251
T operator()(const coordT &x) const
return value at point x; fairly inefficient
Definition: function_interface.h:245
World & world
Definition: function_interface.h:173
std::vector< std::shared_ptr< implL > > impl_p1_vector
supposedly orbital 1
Definition: function_interface.h:183
std::shared_ptr< implL > pimplL
Definition: function_interface.h:171
std::shared_ptr< implL > impl_m1
various MRA functions of MDIM dimensionality (e.g. 3, if NDIM==6)
Definition: function_interface.h:181
FunctionImpl< T, MDIM > implL
Definition: function_interface.h:169
std::shared_ptr< implT > impl_eri
supposedly 1/r12
Definition: function_interface.h:178
std::shared_ptr< implL > impl_m2
supposedly 1/r2
Definition: function_interface.h:182
void make_redundant(const bool fence)
Definition: function_interface.h:212
bool check_redundant() const
return true if all constituent functions are in redundant tree state
Definition: function_interface.h:231
std::vector< std::shared_ptr< implL > > impl_p2_vector
supposedly orbital 2
Definition: function_interface.h:184
Vector< double, NDIM > coordT
Type of vector holding coordinates.
Definition: function_interface.h:167
void replicate_low_dim_functions(const bool fence)
replicate low-dimensional functions over all ranks of this world
Definition: function_interface.h:203
a function like f(x)=1/x
Definition: function_interface.h:526
static constexpr std::size_t LDIM
Definition: function_interface.h:543
GFit< double, LDIM > fit(const double eps) const
derived classes must implement this – cf GFit.h
Definition: function_interface.h:545
ElectronRepulsionInterface(double lo, double eps, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), int kk=FunctionDefaults< NDIM >::get_k())
constructor: cf the Coulomb kernel
Definition: function_interface.h:534
ElementaryInterface (formerly FunctorInterfaceWrapper) interfaces a c-function.
Definition: function_interface.h:266
ElementaryInterface(T(*f)(const coordT &))
Definition: function_interface.h:274
T operator()(const coordT &x) const
You should implement this to return f(x)
Definition: function_interface.h:276
T(* f)(const coordT &)
Definition: function_interface.h:272
coeffT values(const Key< NDIM > &key, const Tensor< double > &quad_x) const
Definition: function_interface.h:278
GenTensor< T > coeffT
Definition: function_interface.h:270
Vector< double, NDIM > coordT
Type of vector holding coordinates.
Definition: function_interface.h:269
FunctionDefaults holds default paramaters as static class members.
Definition: funcdefaults.h:204
static int get_k()
Returns the default wavelet order.
Definition: funcdefaults.h:266
static const double & get_thresh()
Returns the default threshold.
Definition: funcdefaults.h:279
static const BoundaryConditions< NDIM > & get_bc()
Returns the default boundary conditions.
Definition: funcdefaults.h:413
static const Tensor< double > & get_cell_width()
Returns the width of each user cell dimension.
Definition: funcdefaults.h:468
Abstract base class interface required for functors used as input to Functions.
Definition: function_interface.h:68
virtual ~FunctionFunctorInterface()
Definition: function_interface.h:129
virtual T operator()(const Vector< double, NDIM > &x) const =0
You should implement this to return f(x)
virtual bool provides_coeff() const
does this functor directly provide sum coefficients? or only function values?
Definition: function_interface.h:142
virtual void operator()(const Vector< double *, 4 > &xvals, T *fvals, int npts) const
Definition: function_interface.h:106
virtual bool supports_vectorized() const
Does the interface support a vectorized operator()?
Definition: function_interface.h:92
virtual std::vector< Vector< double, NDIM > > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: function_interface.h:122
GenTensor< T > coeffT
Definition: function_interface.h:71
virtual coeffT coeff(const keyT &) const
Definition: function_interface.h:131
virtual coeffT values(const keyT &key, const Tensor< double > &tensor) const
Definition: function_interface.h:136
virtual bool screened(const Vector< double, NDIM > &c1, const Vector< double, NDIM > &c2) const
Can we screen this function based on the bounding box information?
Definition: function_interface.h:87
T value_type
Definition: function_interface.h:73
virtual void operator()(const Vector< double *, 1 > &xvals, T *fvals, int npts) const
Definition: function_interface.h:94
Level special_level_
Definition: function_interface.h:75
void set_length_scale(double lo)
adapt the special level to resolve the smallest length scale
Definition: function_interface.h:80
virtual void operator()(const Vector< double *, 2 > &xvals, T *fvals, int npts) const
Definition: function_interface.h:98
virtual void operator()(const Vector< double *, 5 > &xvals, T *fvals, int npts) const
Definition: function_interface.h:110
virtual Level special_level()
Override this change level refinement for special points (default is 6)
Definition: function_interface.h:127
virtual void operator()(const Vector< double *, 3 > &xvals, T *fvals, int npts) const
Definition: function_interface.h:102
FunctionFunctorInterface()
Definition: function_interface.h:77
Key< NDIM > keyT
Definition: function_interface.h:72
virtual void operator()(const Vector< double *, 6 > &xvals, T *fvals, int npts) const
Definition: function_interface.h:114
FunctionImpl holds all Function state to facilitate shallow copy semantics.
Definition: funcimpl.h:941
FunctionInterface implements a wrapper around any class with the operator()()
Definition: function_interface.h:302
const opT op
Definition: function_interface.h:307
bool provides_coeff() const
does this functor directly provide sum coefficients? or only function values?
Definition: function_interface.h:314
Vector< double, NDIM > coordT
Type of vector holding coordinates.
Definition: function_interface.h:305
GenTensor< T > coeffT
Definition: function_interface.h:304
FunctionInterface(const opT &op)
Definition: function_interface.h:310
T operator()(const coordT &coord) const
You should implement this to return f(x)
Definition: function_interface.h:312
FunctorInterface interfaces a class or struct with an operator()()
Definition: function_interface.h:287
GenTensor< T > coeffT
Definition: function_interface.h:291
Vector< double, NDIM > coordT
Type of vector holding coordinates.
Definition: function_interface.h:290
FunctorInterface(const opT &op)
Definition: function_interface.h:295
opT op
Definition: function_interface.h:293
T operator()(const coordT &x) const
You should implement this to return f(x)
Definition: function_interface.h:297
Definition: gfit.h:57
static GFit BSHFit(double mu, double lo, double hi, double eps, bool prnt=false)
return a fit for the bound-state Helmholtz function
Definition: gfit.h:117
static GFit F12Fit(double gamma, double lo, double hi, double eps, bool prnt=false)
return a fit for the F12 correlation factor
Definition: gfit.h:183
static GFit SlaterFit(double gamma, double lo, double hi, double eps, bool prnt=false)
return a fit for the Slater function
Definition: gfit.h:140
static GFit CoulombFit(double lo, double hi, double eps, bool prnt=false)
return a fit for the Coulomb function
Definition: gfit.h:102
Definition: lowranktensor.h:59
long dim(const int i) const
return the number of entries in dimension i
Definition: lowranktensor.h:391
a function like f(x)=1/x
Definition: function_interface.h:495
GFit< double, LDIM > fit(const double eps) const
derived classes must implement this – cf GFit.h
Definition: function_interface.h:519
static constexpr std::size_t LDIM
Definition: function_interface.h:517
GeneralTwoElectronInterface(OperatorInfo info, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), int kk=FunctionDefaults< NDIM >::get_k())
constructor: cf the Coulomb kernel
Definition: function_interface.h:502
OperatorInfo info
Definition: function_interface.h:516
Key is the index for a node of the 2^NDIM-tree.
Definition: key.h:66
Level level() const
Definition: key.h:159
const Vector< Translation, NDIM > & translation() const
Definition: key.h:164
a function like f(x) = (1 - exp(-mu x))/(2 gamma)
Definition: function_interface.h:602
double mu
Definition: function_interface.h:636
SlaterF12Interface(double mu, double lo, double eps, const BoundaryConditions< 6 > &bc=FunctionDefaults< 6 >::get_bc(), int kk=FunctionDefaults< 6 >::get_k())
constructor: cf the Coulomb kernel
Definition: function_interface.h:610
GFit< double, 3 > fit(const double eps) const
derived classes must implement this – cf GFit.h
Definition: function_interface.h:638
a function like f(x)=exp(-mu x)
Definition: function_interface.h:577
GFit< double, 3 > fit(const double eps) const
derived classes must implement this – cf GFit.h
Definition: function_interface.h:596
SlaterFunctionInterface(double mu, double lo, double eps, const BoundaryConditions< 6 > &bc=FunctionDefaults< 6 >::get_bc(), int kk=FunctionDefaults< 6 >::get_k())
constructor: cf the Coulomb kernel
Definition: function_interface.h:585
double mu
Definition: function_interface.h:594
A slice defines a sub-range or patch of a dimension.
Definition: slice.h:103
float_scalar_type normf() const
Returns the Frobenius norm of the tensor.
Definition: tensor.h:1726
T max(long *ind=0) const
Return the maximum value (and if ind is non-null, its index) in the Tensor.
Definition: tensor.h:1703
base class to compute the wavelet coefficients for an isotropic 2e-operator
Definition: function_interface.h:326
virtual GFit< double, LDIM > fit(const double eps) const =0
derived classes must implement this – cf GFit.h
double hi
the largest length scale that needs to be represented
Definition: function_interface.h:488
coeffT coeff(const Key< NDIM > &key) const
return the coefficients of the function in 6D (x1,y1,z1, x2,y2,z2)
Definition: function_interface.h:356
TwoElectronInterface(double lo, double eps, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), int kk=FunctionDefaults< NDIM >::get_k())
constructor: cf the Coulomb kernel
Definition: function_interface.h:337
int rank
the number of terms in the Gaussian quadrature
Definition: function_interface.h:479
T operator()(const Vector< double, NDIM > &x) const
You should implement this to return f(x)
Definition: function_interface.h:361
std::vector< ConvolutionND< double, NDIM > > ops
storing the coefficients
Definition: function_interface.h:476
bool provides_coeff() const
does this functor directly provide sum coefficients? or only function values?
Definition: function_interface.h:351
Tensor< double > make_coeff(const Key< NDIM > &key) const
make the coefficients from the 1d convolution
Definition: function_interface.h:370
double lo
the smallest length scale that needs to be represented
Definition: function_interface.h:485
static constexpr std::size_t LDIM
Definition: function_interface.h:328
void initialize(const double eps)
initialize the Gaussian fit; uses the virtual function fit() to fit
Definition: function_interface.h:446
int k
the wavelet order
Definition: function_interface.h:482
Tensor< double > map_coeff(const Tensor< double > &c) const
the dimensions are a bit confused (x1,x2, y1,y2, z1,z2) -> (x1,y1,z1, x2,y2,z2)
Definition: function_interface.h:428
GenTensor< T > coeffT
Definition: function_interface.h:331
void fence(bool debug=false)
Synchronizes all processes in communicator AND globally ensures no pending AM or tasks.
Definition: worldgop.cc:161
A parallel world class.
Definition: world.h:132
WorldGopInterface & gop
Global operations.
Definition: world.h:205
Compuates most matrix elements over 1D operators (including Gaussians)
static double lo
Definition: dirac-hatom.cc:23
fit isotropic functions to a set of Gaussians with controlled precision
auto T(World &world, response_space &f) -> response_space
Definition: global_functions.cc:34
Multidimension Key for MRA tree and associated iterators.
static double pow(const double *a, const double *b)
Definition: lda.h:74
#define MADNESS_EXCEPTION(msg, value)
Macro for throwing a MADNESS exception.
Definition: madness_exception.h:119
#define MADNESS_ASSERT(condition)
Assert a condition that should be free of side-effects since in release builds this might be a no-op.
Definition: madness_exception.h:134
#define MADNESS_CHECK_THROW(condition, msg)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition: madness_exception.h:210
Vector< double, 3 > coordT
Definition: mcpfit.cc:48
const double pi
Mathematical constant .
Definition: constants.h:48
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
@ BC_PERIODIC
Definition: funcdefaults.h:56
Vector< double, 3 > coordT
Definition: corepotential.cc:54
const std::vector< Function< T, NDIM > > & change_tree_state(const std::vector< Function< T, NDIM >> &v, const TreeState finalstate, const bool fence=true)
change tree state of the functions
Definition: vmra.h:277
@ redundant
s coeffs everywhere
Definition: funcdefaults.h:64
static double r2(const coord_3d &x)
Definition: smooth.h:45
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition: mra.h:2002
Tensor< double > tensorT
Definition: distpm.cc:21
int64_t Translation
Definition: key.h:54
static const Slice _(0,-1, 1)
int Level
Definition: key.h:55
void print(const T &t, const Ts &... ts)
Print items to std::cout (items separated by spaces) and terminate with a new line.
Definition: print.h:225
Tensor< T > fcube(const Key< NDIM > &, T(*f)(const Vector< double, NDIM > &), const Tensor< double > &)
Definition: mraimpl.h:2163
@ TT_FULL
Definition: gentensor.h:120
NDIM & f
Definition: mra.h:2416
double inner(response_space &a, response_space &b)
Definition: response_functions.h:442
std::vector< Function< T, NDIM > > impl2function(const std::vector< std::shared_ptr< FunctionImpl< T, NDIM >>> vimpl)
Definition: vmra.h:676
std::enable_if< std::is_base_of< ProjectorBase, projT >::value, OuterProjector< projT, projQ > >::type outer(const projT &p0, const projQ &p1)
Definition: projector.h:457
const double mu
Definition: navstokes_cosines.cc:95
static const double c
Definition: relops.cc:10
static const double thresh
Definition: rk.cc:45
static const long k
Definition: rk.cc:44
Definition: test_dc.cc:47
Definition: convolution1d.h:849
Definition: operatorinfo.h:58
double hi
Definition: operatorinfo.h:65
double thresh
Definition: operatorinfo.h:63
Defines and implements most of Tensor.
void d()
Definition: test_sig.cc:79
static const std::size_t NDIM
Definition: testpdiff.cc:42