function_interface.h

Go to the documentation of this file.

/*
  This file is part of MADNESS.

  Copyright (C) 2007,2010 Oak Ridge National Laboratory

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

  For more information please contact:

  Robert J. Harrison
  Oak Ridge National Laboratory
  One Bethel Valley Road
  P.O. Box 2008, MS-6367

  email: harrisonrj@ornl.gov
  tel:   865-241-3937
  fax:   865-572-0680

  $Id: function_factory_and_interface.h 3422 2014-03-24 09:16:15Z 3ru6ruWu $
*/


#ifndef MADNESS_MRA_FUNCTION_INTERFACE_H__INCLUDED
#define MADNESS_MRA_FUNCTION_INTERFACE_H__INCLUDED

#include <madness/tensor/tensor.h>
#include <madness/tensor/gentensor.h>
#include <madness/mra/key.h>

// needed for the TwoElectronInterface
#include <madness/mra/gfit.h>
#include <madness/mra/convolution1d.h>
#include <madness/mra/function_common_data.h>
namespace madness {

        // forward declaration needed for CompositeFunctorInterface
        template<typename T, std::size_t NDIM>
        class FunctionImpl;

    template<typename T, std::size_t NDIM>
    Tensor<T> fcube(const Key<NDIM>&, T (*f)(const Vector<double,NDIM>&), const Tensor<double>&);


        template<typename T, std::size_t NDIM>
        class FunctionFunctorInterface {
        public:

                typedef GenTensor<T> coeffT;
                typedef Key<NDIM> keyT;

                virtual bool screened(const Vector<double,NDIM>& c1, const Vector<double,NDIM>& c2) const {
                  return false;
                }
             
                virtual bool supports_vectorized() const {return false;}

                virtual void operator()(const Vector<double*,1>& xvals, T* fvals, int npts) const {
                    MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0);
                }

                virtual void operator()(const Vector<double*,2>& xvals, T* fvals, int npts) const {
                    MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0);
                }

                virtual void operator()(const Vector<double*,3>& xvals, T* fvals, int npts) const {
                    MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0);
                }

                virtual void operator()(const Vector<double*,4>& xvals, T* fvals, int npts) const {
                    MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0);
                }

                virtual void operator()(const Vector<double*,5>& xvals, T* fvals, int npts) const {
                    MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0);
                }

                virtual void operator()(const Vector<double*,6>& xvals, T* fvals, int npts) const {
                    MADNESS_EXCEPTION("FunctionFunctorInterface: This function should not be called!", 0);
                }

                virtual T operator()(const Vector<double, NDIM>& x) const = 0;

                virtual std::vector< Vector<double,NDIM> > special_points() const {
                        return std::vector< Vector<double,NDIM> >();
                }

                virtual Level special_level() {return 6;}

                virtual ~FunctionFunctorInterface() {}

                virtual coeffT coeff(const keyT&) const {
                        MADNESS_EXCEPTION("implement coeff for FunctionFunctorInterface",0);
                        return coeffT();
                }

        virtual coeffT values(const keyT& key, const Tensor<double>& tensor) const {
            MADNESS_EXCEPTION("implement values for FunctionFunctorInterface",0);
            return coeffT();
        }

                virtual bool provides_coeff() const {
                        return false;
                }

        };



        template<typename T, std::size_t NDIM, std::size_t MDIM>
        class CompositeFunctorInterface : public FunctionFunctorInterface<T,NDIM> {

                typedef Vector<double, NDIM> coordT; 
                typedef FunctionImpl<T,NDIM> implT;
                typedef FunctionImpl<T,MDIM> implL;
                typedef std::shared_ptr<implT> pimplT;
                typedef std::shared_ptr<implL> pimplL;

                World& world;

        public:
                std::shared_ptr<implT> impl_ket;        
                std::shared_ptr<implT> impl_eri;        

                std::shared_ptr<implL> impl_m1; 
                std::shared_ptr<implL> impl_m2; 
                std::shared_ptr<implL> impl_p1; 
                std::shared_ptr<implL> impl_p2; 

        public:

                CompositeFunctorInterface(World& world, pimplT ket, pimplT g12,
                                pimplL v1, pimplL v2, pimplL p1, pimplL p2)
                        : world(world), impl_ket(ket), impl_eri(g12)
                        , impl_m1(v1), impl_m2(v2), impl_p1(p1), impl_p2(p2)
                {

                        // some consistency checks
                        // either a pair ket is provided, or two particles (tba)
                        MADNESS_ASSERT(impl_ket or (impl_p1 and impl_p2));

                        // prepare base functions that make this function
                        if (impl_ket and (not impl_ket->is_on_demand())) impl_ket->make_redundant(false);
                        if (impl_eri) {
                                if (not impl_eri->is_on_demand()) impl_eri->make_redundant(false);
                        }
                        if (impl_m1 and (not impl_m1->is_on_demand())) impl_m1->make_redundant(false);
                        if (impl_m2 and (not impl_m2->is_on_demand())) impl_m2->make_redundant(false);

                        if (impl_p1 and (not impl_p1->is_on_demand())) impl_p1->make_redundant(false);
                        if (impl_p2 and (not impl_p2->is_on_demand())) impl_p2->make_redundant(false);
                        world.gop.fence();

                }

                T operator()(const coordT& x) const {
                        print("there is no operator()(coordT&) in CompositeFunctorInterface, for good reason");
                        MADNESS_ASSERT(0);
                        return T(0);
                };

                bool provides_coeff() const {
                        return false;
                }

        };



        template<typename T, std::size_t NDIM>
        class ElementaryInterface : public FunctionFunctorInterface<T,NDIM> {

        public:
                typedef Vector<double, NDIM> coordT; 
        typedef GenTensor<T> coeffT;

                T (*f)(const coordT&);

                ElementaryInterface(T (*f)(const coordT&)) : f(f) {}

                T operator()(const coordT& x) const {return f(x);}

                coeffT values(const Key<NDIM>& key, const Tensor<double>& quad_x) const {
                typedef Tensor<T> tensorT;
            tensorT fval=madness::fcube(key,f,quad_x);
            return coeffT(fval,FunctionDefaults<NDIM>::get_thresh(),TT_FULL);
                }
        };

    template<typename T, std::size_t NDIM, typename opT>
    class FunctorInterface : public FunctionFunctorInterface<T,NDIM> {

    public:
        typedef Vector<double, NDIM> coordT; 
        typedef GenTensor<T> coeffT;

        opT op;

        FunctorInterface(const opT& op) : op(op) {}

        T operator()(const coordT& x) const {return op(x);}
    };

        template<typename T, size_t NDIM, typename opT>
        class FunctionInterface : public FunctionFunctorInterface<T,NDIM> {

            typedef GenTensor<T> coeffT;
        typedef Vector<double, NDIM> coordT; 

        const opT op;

    public:
        FunctionInterface(const opT& op) : op(op) {}

        T operator()(const coordT& coord) const {return op(coord);}

        bool provides_coeff() const {return false;}

        };


        template<typename T, std::size_t NDIM>
    class TwoElectronInterface : public FunctionFunctorInterface<T,NDIM> {
        public:

                typedef GenTensor<T> coeffT;


                TwoElectronInterface(double lo,double eps,
                                const BoundaryConditions<6>& bc=FunctionDefaults<6>::get_bc(),
                                int kk=FunctionDefaults<6>::get_k())
                                :rank(), k(kk), lo(lo), hi(1.0) {

                        // Presently we must have periodic or non-periodic in all dimensions.
                        for (std::size_t d=1; d<6; ++d) {MADNESS_ASSERT(bc(d,0)==bc(0,0));}

                        const Tensor<double>& width = FunctionDefaults<6>::get_cell_width();
                        hi = width.normf(); // Diagonal width of cell
                        if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation

                }

                bool provides_coeff() const {
                        return true;
                }

                coeffT coeff(const Key<NDIM>& key) const {
                        Tensor<double> c=make_coeff(key);
            return coeffT(map_coeff(c),FunctionDefaults<6>::get_thresh(),TT_FULL);
                }

                T operator()(const Vector<double, NDIM>& x) const {
                        print("there is no operator()(coordT&) in TwoElectronInterface, for good reason");
                        MADNESS_ASSERT(0);
                        return T(0);
                }

        protected:

                Tensor<double> make_coeff(const Key<6>& key) const {
                        const Level n=key.level();
                        const Vector<Translation,6> l=key.translation();

                        // get the displacements for all 3 dimensions: x12, y12, z12
                        const Translation l0=(l[0]-l[3]);
                        const Translation l1=(l[1]-l[4]);
                        const Translation l2=(l[2]-l[5]);

                        Tensor<double> scr1(rank,k*k), scr2(rank,k*k,k*k);

                        // lump all the terms together
                        for (long mu=0; mu<rank; mu++) {
                                const Tensor<double> r0=(ops[mu].getop(0)->rnlij(n,l0)).reshape(k*k);
                                const Tensor<double> r1=(ops[mu].getop(1)->rnlij(n,l1)).reshape(k*k);
                                const Tensor<double> r2=(ops[mu].getop(2)->rnlij(n,l2)).reshape(k*k);

                                // include weights in first vector
                                scr1(mu,Slice(_))=r0*ops[mu].getfac();

                                // merge second and third vector to scr(r,k1,k2)
                                scr2(mu,Slice(_),Slice(_))=outer(r1,r2);
                        }

                        Tensor<double> c=inner(scr1,scr2,0,0);
                        return c;
                }

                Tensor<double> map_coeff(const Tensor<double>& c) const {
                        std::vector<long> map(6);
                        map[0]=0;       map[1]=3;       map[2]=1;
                        map[3]=4;       map[4]=2;       map[5]=5;
                        return copy(c.reshape(k,k,k,k,k,k).mapdim(map));
                }

                void initialize(const double eps) {
                        GFit<double,3> fit=this->fit(eps);
                        Tensor<double> coeff=fit.coeffs();
                        Tensor<double> expnt=fit.exponents();

                        // set some parameters
                        rank=coeff.dim(0);
                        ops.resize(rank);
                        const Tensor<double>& width = FunctionDefaults<6>::get_cell_width();

                        // construct all the terms
                        for (int mu=0; mu<rank; ++mu) {
                                //                double c = std::pow(sqrt(expnt(mu)/pi),static_cast<int>(NDIM)); // Normalization coeff
                                double c = std::pow(sqrt(expnt(mu)/constants::pi),3); // Normalization coeff

                                // We cache the normalized operator so the factor is the value we must multiply
                                // by to recover the coeff we want.
                                ops[mu].setfac(coeff(mu)/c);

                                // only 3 dimensions here!
                                for (std::size_t d=0; d<3; ++d) {
                                        ops[mu].setop(d,GaussianConvolution1DCache<double>::get(k, expnt(mu)*width[d]*width[d], 0, false));
                                }
                        }
                }

                virtual GFit<double,3> fit(const double eps) const = 0;

                mutable std::vector< ConvolutionND<double,6> > ops;

                int rank;

                int k;

                double lo;

                double hi;

        };

        class ElectronRepulsionInterface : public TwoElectronInterface<double,6> {
        public:


                ElectronRepulsionInterface(double lo,double eps,
                                const BoundaryConditions<6>& bc=FunctionDefaults<6>::get_bc(),
                                int kk=FunctionDefaults<6>::get_k())
                  : TwoElectronInterface<double,6>(lo,eps,bc,kk) {

                        initialize(eps);
                }

        private:

                GFit<double,3> fit(const double eps) const {
                        return GFit<double,3>::CoulombFit(lo,hi,eps,false);
                }
        };

        class BSHFunctionInterface : public TwoElectronInterface<double,6> {
        public:


                BSHFunctionInterface(double mu, double lo,double eps,
                                const BoundaryConditions<6>& bc=FunctionDefaults<6>::get_bc(),
                                int kk=FunctionDefaults<6>::get_k())
                  : TwoElectronInterface<double,6>(lo,eps,bc,kk), mu(mu) {

                        initialize(eps);
                }

        private:

                double mu;

                GFit<double,3> fit(const double eps) const {
                        return GFit<double,3>::BSHFit(mu,lo,hi,eps,false);
                }
        };

        class SlaterFunctionInterface : public TwoElectronInterface<double,6> {
        public:


                SlaterFunctionInterface(double mu, double lo,double eps,
                                const BoundaryConditions<6>& bc=FunctionDefaults<6>::get_bc(),
                                int kk=FunctionDefaults<6>::get_k())
                  : TwoElectronInterface<double,6>(lo,eps,bc,kk), mu(mu) {

                        initialize(eps);
                }

        private:

                double mu;

                GFit<double,3> fit(const double eps) const {
                        return GFit<double,3>::SlaterFit(mu,lo,hi,eps,false);
                }
        };

        class SlaterF12Interface : public TwoElectronInterface<double,6> {
        public:


                SlaterF12Interface(double mu, double lo,double eps,
                                const BoundaryConditions<6>& bc=FunctionDefaults<6>::get_bc(),
                                int kk=FunctionDefaults<6>::get_k())
                  : TwoElectronInterface<double,6>(lo,eps,bc,kk), mu(mu) {

                        initialize(eps);
                }

                coeffT coeff(const Key<6>& key) const {

                        Tensor<double> c=make_coeff(key);

                        // subtract 1 from the (0,0,..,0) element of the tensor,
                        // which is the 0th order polynomial coefficient
                double one_coeff1=1.0*sqrt(FunctionDefaults<6>::get_cell_volume())
                                *pow(0.5,0.5*6*key.level());
            std::vector<long> v0(6,0L);
            c(v0)-=one_coeff1;

                        c.scale(-0.5/mu);
            return coeffT(map_coeff(c),FunctionDefaults<6>::get_thresh(),TT_FULL);
                }

        private:

                double mu;

                GFit<double,3> fit(const double eps) const {
                        return GFit<double,3>::SlaterFit(mu,lo,hi,eps,false);
                }
        };

        class FGInterface : public TwoElectronInterface<double,6> {
        public:


                FGInterface(double mu, double lo,double eps,
                                const BoundaryConditions<6>& bc=FunctionDefaults<6>::get_bc(),
                                int kk=FunctionDefaults<6>::get_k())
                  : TwoElectronInterface<double,6>(lo,eps,bc,kk), mu(mu) {

                        initialize(eps);
                }

        private:

                double mu;

                GFit<double,3> fit(const double eps) const {
                        return GFit<double,3>::SlaterFit(mu,lo,hi,eps,false);
                }
        };


#if 0


        template<typename T, std::size_t NDIM>
        class ElectronRepulsionInterface : public FunctionFunctorInterface<T,NDIM> {

                typedef GenTensor<T> coeffT;
                typedef Vector<double, NDIM> coordT; 

                ElectronRepulsion eri;

        public:

                ElectronRepulsionInterface(World& world,double lo,double eps,
                const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(),
                int k=FunctionDefaults<NDIM>::get_k())
                        : eri(ElectronRepulsion(eps,eps,bc,k)) {
                }


                T operator()(const coordT& x) const {
                        print("there is no operator()(coordT&) in ElectronRepulsionInterface, for good reason");
                        MADNESS_ASSERT(0);
                        return T(0);
                };


                coeffT coeff(const Key<NDIM>& key) const {
            return coeffT(this->eri.coeff(key),FunctionDefaults<NDIM>::get_thresh(),
                    TT_FULL);
                }

        };


        template<typename T, std::size_t NDIM>
        class FGIntegralInterface : public FunctionFunctorInterface<T,NDIM> {

                typedef GenTensor<T> coeffT;
                typedef Vector<double, NDIM> coordT; 

                ElectronRepulsion eri;
                BSHFunction bsh;

        public:

                FGIntegralInterface(World& world, double lo, double eps, double gamma,
                const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(),
                int k=FunctionDefaults<NDIM>::get_k())
                        : eri(ElectronRepulsion(eps,eps,0.0,bc,k))
                        , bsh(BSHFunction(eps,eps,gamma,bc,k)) {
                }

                bool provides_coeff() const {
                        return true;
                }

                T operator()(const coordT& x) const {
                        print("there is no operator()(coordT&) in FGIntegralInterface, for good reason");
                        MADNESS_ASSERT(0);
                        return T(0);
                };

                coeffT coeff(const Key<NDIM>& key) const {
                typedef Tensor<T> tensorT;
                        tensorT e_b=eri.coeff(key)-bsh.coeff(key);
            return coeffT(e_b,FunctionDefaults<NDIM>::get_thresh(),TT_FULL);
                }

        };

#endif

}

#endif // MADNESS_MRA_FUNCTION_INTERFACE_H__INCLUDED