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> 39 #include <madness/tensor/gentensor.h> 40 #include <madness/mra/key.h> 41 42 // needed for the TwoElectronInterface 43 #include <madness/mra/gfit.h> 44 #include <madness/mra/convolution1d.h> 45 #include <madness/mra/function_common_data.h> 46 namespace madness { 47 48 // forward declaration needed for CompositeFunctorInterface 49 template<typename T, std::size_t NDIM> 50 class FunctionImpl; 51 52 template<typename T, std::size_t NDIM> 53 Tensor<T> fcube(const Key<NDIM>&, T (*f)(const Vector<double,NDIM>&), const Tensor<double>&); 54 55 56 57 /// Abstract base class interface required for functors used as input to Functions 58 template<typename T, std::size_t NDIM> 59 class FunctionFunctorInterface { 60 public: 61 62 typedef GenTensor<T> coeffT; 63 typedef Key<NDIM> keyT; 64 typedef T value_type; 65 66 Level special_level_; 67 FunctionFunctorInterface()68 FunctionFunctorInterface() : special_level_(6) {} 69 70 /// adapt the special level to resolve the smallest length scale set_length_scale(double lo)71 void set_length_scale(double lo) { 72 double Lmax=FunctionDefaults<NDIM>::get_cell_width().max(); 73 double lo_sim=lo/Lmax; // lo in simulation coordinates; 74 special_level_=Level(-log2(lo_sim)); 75 } 76 77 /// Can we screen this function based on the bounding box information? screened(const Vector<double,NDIM> & c1,const Vector<double,NDIM> & c2)78 virtual bool screened(const Vector<double,NDIM>& c1, const Vector<double,NDIM>& c2) const { 79 return false; 80 } 81 82 /// Does the interface support a vectorized operator()? supports_vectorized()83 virtual bool supports_vectorized() const {return false;} 84 operator()85 virtual void operator()(const Vector<double*,1>& xvals, T* fvals, int npts) const { 86 MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0); 87 } 88 operator()89 virtual void operator()(const Vector<double*,2>& xvals, T* fvals, int npts) const { 90 MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0); 91 } 92 operator()93 virtual void operator()(const Vector<double*,3>& xvals, T* fvals, int npts) const { 94 MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0); 95 } 96 operator()97 virtual void operator()(const Vector<double*,4>& xvals, T* fvals, int npts) const { 98 MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0); 99 } 100 operator()101 virtual void operator()(const Vector<double*,5>& xvals, T* fvals, int npts) const { 102 MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0); 103 } 104 operator()105 virtual void operator()(const Vector<double*,6>& xvals, T* fvals, int npts) const { 106 MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0); 107 } 108 109 /// You should implement this to return \c f(x) 110 virtual T operator()(const Vector<double, NDIM>& x) const = 0; 111 112 /// Override this to return list of special points to be refined more deeply special_points()113 virtual std::vector< Vector<double,NDIM> > special_points() const { 114 return std::vector< Vector<double,NDIM> >(); 115 } 116 117 /// Override this change level refinement for special points (default is 6) special_level()118 virtual Level special_level() {return special_level_;} 119 ~FunctionFunctorInterface()120 virtual ~FunctionFunctorInterface() {} 121 coeff(const keyT &)122 virtual coeffT coeff(const keyT&) const { 123 MADNESS_EXCEPTION("implement coeff for FunctionFunctorInterface",0); 124 return coeffT(); 125 } 126 values(const keyT & key,const Tensor<double> & tensor)127 virtual coeffT values(const keyT& key, const Tensor<double>& tensor) const { 128 MADNESS_EXCEPTION("implement values for FunctionFunctorInterface",0); 129 return coeffT(); 130 } 131 132 /// does this functor directly provide sum coefficients? or only function values? provides_coeff()133 virtual bool provides_coeff() const { 134 return false; 135 } 136 137 }; 138 139 140 141 ///forward declaration 142 template <typename T, std::size_t NDIM> 143 // void FunctionImpl<T,NDIM>::fcube(const keyT& key, const FunctionFunctorInterface<T,NDIM>& f, const Tensor<double>& qx, tensorT& fval) const { 144 void fcube(const Key<NDIM>& key, const FunctionFunctorInterface<T,NDIM>& f, const Tensor<double>& qx, Tensor<T>& fval); 145 146 /// CompositeFunctorInterface implements a wrapper of holding several functions and functors 147 148 /// Use this to "connect" several functions and/or functors and to return their coefficients 149 /// e.g. connect f1 and f2 with an addition, you can request the coefficients of any node 150 /// and they will be computed on the fly and returned. Mainly useful to connect a functor 151 /// with a function, if the functor is too large to be represented in MRA (e.g. 1/r12) 152 /// 153 /// as of now, the operation connecting the functions/functors is simply addition. 154 /// need to implement expression templates, if I only knew what that was... 155 template<typename T, std::size_t NDIM, std::size_t MDIM> 156 class CompositeFunctorInterface : public FunctionFunctorInterface<T,NDIM> { 157 158 typedef Vector<double, NDIM> coordT; ///< Type of vector holding coordinates 159 typedef FunctionImpl<T,NDIM> implT; 160 typedef FunctionImpl<T,MDIM> implL; 161 typedef std::shared_ptr<implT> pimplT; 162 typedef std::shared_ptr<implL> pimplL; 163 164 World& world; 165 166 public: 167 /// various MRA functions of NDIM dimensionality 168 std::shared_ptr<implT> impl_ket; ///< supposedly the pair function 169 std::shared_ptr<implT> impl_eri; ///< supposedly 1/r12 170 171 /// various MRA functions of MDIM dimensionality (e.g. 3, if NDIM==6) 172 std::shared_ptr<implL> impl_m1; ///< supposedly 1/r1 173 std::shared_ptr<implL> impl_m2; ///< supposedly 1/r2 174 std::shared_ptr<implL> impl_p1; ///< supposedly orbital 1 175 std::shared_ptr<implL> impl_p2; ///< supposedly orbital 2 176 177 public: 178 179 /// constructor takes its Factory CompositeFunctorInterface(World & world,pimplT ket,pimplT g12,pimplL v1,pimplL v2,pimplL p1,pimplL p2)180 CompositeFunctorInterface(World& world, pimplT ket, pimplT g12, 181 pimplL v1, pimplL v2, pimplL p1, pimplL p2) 182 : world(world), impl_ket(ket), impl_eri(g12) 183 , impl_m1(v1), impl_m2(v2), impl_p1(p1), impl_p2(p2) 184 { 185 186 // some consistency checks 187 // either a pair ket is provided, or two particles (tba) 188 MADNESS_ASSERT(impl_ket or (impl_p1 and impl_p2)); 189 190 // prepare base functions that make this function 191 if (impl_ket and (not impl_ket->is_on_demand())) impl_ket->make_redundant(false); 192 if (impl_eri) { 193 if (not impl_eri->is_on_demand()) impl_eri->make_redundant(false); 194 } 195 if (impl_m1 and (not impl_m1->is_on_demand())) impl_m1->make_redundant(false); 196 if (impl_m2 and (not impl_m2->is_on_demand())) impl_m2->make_redundant(false); 197 198 if (impl_p1 and (not impl_p1->is_on_demand())) impl_p1->make_redundant(false); 199 if (impl_p2 and (not impl_p2->is_on_demand())) impl_p2->make_redundant(false); 200 world.gop.fence(); 201 202 } 203 204 /// return value at point x; fairly inefficient operator()205 T operator()(const coordT& x) const { 206 print("there is no operator()(coordT&) in CompositeFunctorInterface, for good reason"); 207 MADNESS_ASSERT(0); 208 return T(0); 209 }; 210 provides_coeff()211 bool provides_coeff() const { 212 return false; 213 } 214 215 }; 216 217 218 /// ElementaryInterface (formerly FunctorInterfaceWrapper) interfaces a c-function 219 220 /// hard-code your favorite function and interface it with this; Does only 221 /// provide function values, no MRA coefficients. Care must be taken if the 222 /// function we refer to is a singular function, and a on-demand function 223 /// at the same time, since direct computation of coefficients via mraimpl::project 224 /// might suffer from inaccurate quadrature. 225 template<typename T, std::size_t NDIM> 226 class ElementaryInterface : public FunctionFunctorInterface<T,NDIM> { 227 228 public: 229 typedef Vector<double, NDIM> coordT; ///< Type of vector holding coordinates 230 typedef GenTensor<T> coeffT; 231 232 T (*f)(const coordT&); 233 ElementaryInterface(T (* f)(const coordT &))234 ElementaryInterface(T (*f)(const coordT&)) : f(f) {} 235 operator()236 T operator()(const coordT& x) const {return f(x);} 237 values(const Key<NDIM> & key,const Tensor<double> & quad_x)238 coeffT values(const Key<NDIM>& key, const Tensor<double>& quad_x) const { 239 typedef Tensor<T> tensorT; 240 tensorT fval=madness::fcube(key,f,quad_x); 241 return coeffT(fval,FunctionDefaults<NDIM>::get_thresh(),TT_FULL); 242 } 243 }; 244 245 /// FunctorInterface interfaces a class or struct with an operator()() 246 template<typename T, std::size_t NDIM, typename opT> 247 class FunctorInterface : public FunctionFunctorInterface<T,NDIM> { 248 249 public: 250 typedef Vector<double, NDIM> coordT; ///< Type of vector holding coordinates 251 typedef GenTensor<T> coeffT; 252 253 opT op; 254 FunctorInterface(const opT & op)255 FunctorInterface(const opT& op) : op(op) {} 256 operator()257 T operator()(const coordT& x) const {return op(x);} 258 }; 259 260 /// FunctionInterface implements a wrapper around any class with the operator()() 261 template<typename T, size_t NDIM, typename opT> 262 class FunctionInterface : public FunctionFunctorInterface<T,NDIM> { 263 264 typedef GenTensor<T> coeffT; 265 typedef Vector<double, NDIM> coordT; ///< Type of vector holding coordinates 266 267 const opT op; 268 269 public: FunctionInterface(const opT & op)270 FunctionInterface(const opT& op) : op(op) {} 271 operator()272 T operator()(const coordT& coord) const {return op(coord);} 273 provides_coeff()274 bool provides_coeff() const {return false;} 275 276 }; 277 278 /// base class to compute the wavelet coefficients for an isotropic 2e-operator 279 280 /// all classes that derive from this base class use the Gaussian fitting 281 /// procedure that has been developed for the BSH operator. We simply 282 /// reuse the wavelet coefficients that we get from there to avoid 283 /// evaluating the functions themselves, since the quadrature of singular 284 /// functions is imprecise and slow. 285 template<typename T, std::size_t NDIM> 286 class TwoElectronInterface : public FunctionFunctorInterface<T,NDIM> { 287 public: 288 289 typedef GenTensor<T> coeffT; 290 291 /// constructor: cf the Coulomb kernel 292 293 /// @param[in] lo the smallest length scale to be resolved 294 /// @param[in] eps the accuracy threshold 295 TwoElectronInterface(double lo, double eps, 296 const BoundaryConditions<6>& bc=FunctionDefaults<6>::get_bc(), 297 int kk=FunctionDefaults<6>::get_k()) rank()298 :rank(), k(kk), lo(lo), hi(1.0) { 299 300 // Presently we must have periodic or non-periodic in all dimensions. 301 for (std::size_t d=1; d<6; ++d) {MADNESS_ASSERT(bc(d,0)==bc(0,0));} 302 303 const Tensor<double>& width = FunctionDefaults<6>::get_cell_width(); 304 hi = width.normf(); // Diagonal width of cell 305 if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation 306 307 } 308 provides_coeff()309 bool provides_coeff() const { 310 return true; 311 } 312 313 /// return the coefficients of the function in 6D (x1,y1,z1, x2,y2,z2) coeff(const Key<NDIM> & key)314 coeffT coeff(const Key<NDIM>& key) const { 315 Tensor<double> c=make_coeff(key); 316 return coeffT(map_coeff(c),FunctionDefaults<6>::get_thresh(),TT_FULL); 317 } 318 operator()319 T operator()(const Vector<double, NDIM>& x) const { 320 print("there is no operator()(coordT&) in TwoElectronInterface, for good reason"); 321 MADNESS_ASSERT(0); 322 return T(0); 323 } 324 325 protected: 326 327 /// make the coefficients from the 1d convolution make_coeff(const Key<6> & key)328 Tensor<double> make_coeff(const Key<6>& key) const { 329 const Level n=key.level(); 330 const Vector<Translation,6> l=key.translation(); 331 332 // get the displacements for all 3 dimensions: x12, y12, z12 333 const Translation l0=(l[0]-l[3]); 334 const Translation l1=(l[1]-l[4]); 335 const Translation l2=(l[2]-l[5]); 336 337 Tensor<double> scr1(rank,k*k), scr2(rank,k*k,k*k); 338 339 // lump all the terms together 340 for (long mu=0; mu<rank; mu++) { 341 const Tensor<double> r0=(ops[mu].getop(0)->rnlij(n,l0)).reshape(k*k); 342 const Tensor<double> r1=(ops[mu].getop(1)->rnlij(n,l1)).reshape(k*k); 343 const Tensor<double> r2=(ops[mu].getop(2)->rnlij(n,l2)).reshape(k*k); 344 345 // include weights in first vector 346 scr1(mu,Slice(_))=r0*ops[mu].getfac(); 347 348 // merge second and third vector to scr(r,k1,k2) 349 scr2(mu,Slice(_),Slice(_))=outer(r1,r2); 350 } 351 352 Tensor<double> c=inner(scr1,scr2,0,0); 353 return c; 354 } 355 356 /// the dimensions are a bit confused (x1,x2, y1,y2, z1,z2) -> (x1,y1,z1, x2,y2,z2) map_coeff(const Tensor<double> & c)357 Tensor<double> map_coeff(const Tensor<double>& c) const { 358 std::vector<long> map(6); 359 map[0]=0; map[1]=3; map[2]=1; 360 map[3]=4; map[4]=2; map[5]=5; 361 return copy(c.reshape(k,k,k,k,k,k).mapdim(map)); 362 } 363 364 /// initialize the Gaussian fit; uses the virtual function fit() to fit initialize(const double eps)365 void initialize(const double eps) { 366 GFit<double,3> fit=this->fit(eps); 367 Tensor<double> coeff=fit.coeffs(); 368 Tensor<double> expnt=fit.exponents(); 369 370 // set some parameters 371 rank=coeff.dim(0); 372 ops.resize(rank); 373 const Tensor<double>& width = FunctionDefaults<6>::get_cell_width(); 374 375 // construct all the terms 376 for (int mu=0; mu<rank; ++mu) { 377 // double c = std::pow(sqrt(expnt(mu)/pi),static_cast<int>(NDIM)); // Normalization coeff 378 double c = std::pow(sqrt(expnt(mu)/constants::pi),3); // Normalization coeff 379 380 // We cache the normalized operator so the factor is the value we must multiply 381 // by to recover the coeff we want. 382 ops[mu].setfac(coeff(mu)/c); 383 384 // only 3 dimensions here! 385 for (std::size_t d=0; d<3; ++d) { 386 ops[mu].setop(d,GaussianConvolution1DCache<double>::get(k, expnt(mu)*width[d]*width[d], 0, false)); 387 } 388 } 389 } 390 391 /// derived classes must implement this -- cf GFit.h 392 virtual GFit<double,3> fit(const double eps) const = 0; 393 394 /// storing the coefficients 395 mutable std::vector< ConvolutionND<double,6> > ops; 396 397 /// the number of terms in the Gaussian quadrature 398 int rank; 399 400 /// the wavelet order 401 int k; 402 403 /// the smallest length scale that needs to be represented 404 double lo; 405 406 /// the largest length scale that needs to be represented 407 double hi; 408 409 }; 410 411 /// a function like f(x)=1/x 412 class ElectronRepulsionInterface : public TwoElectronInterface<double,6> { 413 public: 414 415 /// constructor: cf the Coulomb kernel 416 417 /// @param[in] lo the smallest length scale to be resolved 418 /// @param[in] eps the accuracy threshold 419 ElectronRepulsionInterface(double lo,double eps, 420 const BoundaryConditions<6>& bc=FunctionDefaults<6>::get_bc(), 421 int kk=FunctionDefaults<6>::get_k()) 422 : TwoElectronInterface<double,6>(lo,eps,bc,kk) { 423 424 initialize(eps); 425 } 426 427 private: 428 fit(const double eps)429 GFit<double,3> fit(const double eps) const { 430 return GFit<double,3>::CoulombFit(lo,hi,eps,false); 431 } 432 }; 433 434 /// a function like f(x) = exp(-mu x)/x 435 class BSHFunctionInterface : public TwoElectronInterface<double,6> { 436 public: 437 438 /// constructor: cf the Coulomb kernel 439 440 /// @param[in] mu the exponent of the BSH/inverse Laplacian 441 /// @param[in] lo the smallest length scale to be resolved 442 /// @param[in] eps the accuracy threshold 443 BSHFunctionInterface(double mu, double lo, double eps, 444 const BoundaryConditions<6>& bc=FunctionDefaults<6>::get_bc(), 445 int kk=FunctionDefaults<6>::get_k()) 446 : TwoElectronInterface<double,6>(lo,eps,bc,kk), mu(mu) { 447 448 initialize(eps); 449 } 450 451 private: 452 453 double mu; 454 fit(const double eps)455 GFit<double,3> fit(const double eps) const { 456 return GFit<double,3>::BSHFit(mu,lo,hi,eps,false); 457 } 458 }; 459 460 /// a function like f(x)=exp(-mu x) 461 class SlaterFunctionInterface : public TwoElectronInterface<double,6> { 462 public: 463 464 /// constructor: cf the Coulomb kernel 465 466 /// @param[in] mu the exponent of the Slater function 467 /// @param[in] lo the smallest length scale to be resolved 468 /// @param[in] eps the accuracy threshold 469 SlaterFunctionInterface(double mu, double lo, double eps, 470 const BoundaryConditions<6>& bc=FunctionDefaults<6>::get_bc(), 471 int kk=FunctionDefaults<6>::get_k()) 472 : TwoElectronInterface<double,6>(lo,eps,bc,kk), mu(mu) { 473 initialize(eps); 474 } 475 476 private: 477 478 double mu; 479 fit(const double eps)480 GFit<double,3> fit(const double eps) const { 481 return GFit<double,3>::SlaterFit(mu,lo,hi,eps,false); 482 } 483 }; 484 485 /// a function like f(x) = (1 - exp(-mu x))/(2 gamma) 486 class SlaterF12Interface : public TwoElectronInterface<double,6> { 487 public: 488 489 /// constructor: cf the Coulomb kernel 490 491 /// @param[in] mu the exponent of the Slater function 492 /// @param[in] lo the smallest length scale to be resolved 493 /// @param[in] eps the accuracy threshold 494 SlaterF12Interface(double mu, double lo, double eps, 495 const BoundaryConditions<6>& bc=FunctionDefaults<6>::get_bc(), 496 int kk=FunctionDefaults<6>::get_k()) 497 : TwoElectronInterface<double,6>(lo,eps,bc,kk), mu(mu) { 498 499 initialize(eps); 500 } 501 502 /// overload the function of the base class coeff(const Key<6> & key)503 coeffT coeff(const Key<6>& key) const { 504 505 Tensor<double> c=make_coeff(key); 506 507 // subtract 1 from the (0,0,..,0) element of the tensor, 508 // which is the 0th order polynomial coefficient 509 double one_coeff1=1.0*sqrt(FunctionDefaults<6>::get_cell_volume()) 510 *pow(0.5,0.5*6*key.level()); 511 std::vector<long> v0(6,0L); 512 c(v0)-=one_coeff1; 513 514 c.scale(-0.5/mu); 515 return coeffT(map_coeff(c),FunctionDefaults<6>::get_thresh(),TT_FULL); 516 } 517 518 private: 519 520 double mu; 521 fit(const double eps)522 GFit<double,3> fit(const double eps) const { 523 return GFit<double,3>::SlaterFit(mu,lo,hi,eps,false); 524 } 525 }; 526 527 // Not right 528 // /// a function like f(x) = (1 - exp(-mu x))/x 529 // class FGInterface : public TwoElectronInterface<double,6> { 530 // public: 531 // 532 // /// constructor: cf the Coulomb kernel 533 // 534 // /// @param[in] mu the exponent of the Slater function 535 // /// @param[in] lo the smallest length scale to be resolved 536 // /// @param[in] eps the accuracy threshold 537 // FGInterface(double mu, double lo, double eps, 538 // const BoundaryConditions<6>& bc=FunctionDefaults<6>::get_bc(), 539 // int kk=FunctionDefaults<6>::get_k()) 540 // : TwoElectronInterface<double,6>(lo,eps,bc,kk), mu(mu) { 541 // 542 // initialize(eps); 543 // } 544 // 545 // private: 546 // 547 // double mu; 548 // 549 // GFit<double,3> fit(const double eps) const { 550 // return GFit<double,3>::SlaterFit(mu,lo,hi,eps,false); 551 // } 552 // }; 553 554 555 #if 0 556 557 /// ElectronRepulsionInterface implements the electron repulsion term 1/r12 558 559 /// this is essentially just a wrapper around ElectronRepulsion 560 template<typename T, std::size_t NDIM> 561 class ElectronRepulsionInterface : public FunctionFunctorInterface<T,NDIM> { 562 563 typedef GenTensor<T> coeffT; 564 typedef Vector<double, NDIM> coordT; ///< Type of vector holding coordinates 565 566 /// the class computing the coefficients 567 ElectronRepulsion eri; 568 569 public: 570 571 /// constructor takes the same parameters as the Coulomb operator 572 /// which it uses to compute the coefficients 573 ElectronRepulsionInterface(World& world,double lo,double eps, 574 const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(), 575 int k=FunctionDefaults<NDIM>::get_k()) 576 : eri(ElectronRepulsion(eps,eps,bc,k)) { 577 } 578 579 580 /// return value at point x; fairly inefficient 581 T operator()(const coordT& x) const { 582 print("there is no operator()(coordT&) in ElectronRepulsionInterface, for good reason"); 583 MADNESS_ASSERT(0); 584 return T(0); 585 }; 586 587 588 /// return sum coefficients for imagined node at key 589 coeffT coeff(const Key<NDIM>& key) const { 590 return coeffT(this->eri.coeff(key),FunctionDefaults<NDIM>::get_thresh(), 591 TT_FULL); 592 } 593 594 }; 595 596 /// FGIntegralInterface implements the two-electron integral (1-exp(-gamma*r12))/r12 597 598 /// this is essentially just a wrapper around ElectronRepulsion 599 /// The integral expressed as: 1/r12 - exp(-gamma*r12)/r12 600 /// which can be expressed with an eri and a bsh 601 template<typename T, std::size_t NDIM> 602 class FGIntegralInterface : public FunctionFunctorInterface<T,NDIM> { 603 604 typedef GenTensor<T> coeffT; 605 typedef Vector<double, NDIM> coordT; ///< Type of vector holding coordinates 606 607 /// the class computing the coefficients 608 ElectronRepulsion eri; 609 BSHFunction bsh; 610 611 public: 612 613 /// constructor takes the same parameters as the Coulomb operator 614 /// which it uses to compute the coefficients 615 FGIntegralInterface(World& world, double lo, double eps, double gamma, 616 const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(), 617 int k=FunctionDefaults<NDIM>::get_k()) 618 : eri(ElectronRepulsion(eps,eps,0.0,bc,k)) 619 , bsh(BSHFunction(eps,eps,gamma,bc,k)) { 620 } 621 622 bool provides_coeff() const { 623 return true; 624 } 625 626 /// return value at point x; fairly inefficient 627 T operator()(const coordT& x) const { 628 print("there is no operator()(coordT&) in FGIntegralInterface, for good reason"); 629 MADNESS_ASSERT(0); 630 return T(0); 631 }; 632 633 /// return sum coefficients for imagined node at key 634 coeffT coeff(const Key<NDIM>& key) const { 635 typedef Tensor<T> tensorT; 636 tensorT e_b=eri.coeff(key)-bsh.coeff(key); 637 return coeffT(e_b,FunctionDefaults<NDIM>::get_thresh(),TT_FULL); 638 } 639 640 }; 641 642 #endif 643 644 } 645 646 #endif // MADNESS_MRA_FUNCTION_INTERFACE_H__INCLUDED 647