chem/xcfunctional.h

Go to the documentation of this file.

#ifndef MADNESS_CHEM_XCFUNCTIONAL_H__INCLUDED
#define MADNESS_CHEM_XCFUNCTIONAL_H__INCLUDED

#include <madness/madness_config.h>


#include <madness/tensor/tensor.h>
#include <vector>
#include <algorithm>
#include <utility>
#include <madness/mra/key.h>

#ifdef MADNESS_HAS_LIBXC
#include <xc.h>
#endif

namespace madness {
struct xc_lda_potential {
    xc_lda_potential() {}

    void operator()(const Key<3> & key, Tensor<double>& t) const
    {
        int x_rks_s__(const double *r__, double *f, double * dfdra);
        int c_rks_vwn5__(const double *r__, double *f, double * dfdra);
        double* rho = t.ptr();
        for (int i=0; i<t.size(); i++) {
            double r = std::max(rho[i],1e-12);
            double q, dq1, dq2;
            x_rks_s__(&r, &q, &dq1);
            c_rks_vwn5__(&r, &q, &dq2);
            rho[i] = dq1 + dq2;
        }
    }
};

class XCfunctional {
protected:
    bool spin_polarized;        
    double hf_coeff;            
    double rhomin, rhotol, sigmin, sigtol; // See initialize and munge*

#ifdef MADNESS_HAS_LIBXC
    std::vector< std::pair<xc_func_type*,double> > funcs;
    void make_libxc_args(const std::vector< madness::Tensor<double> >& t,
                         madness::Tensor<double>& rho,
                         madness::Tensor<double>& sigma) const;
    int nderiv;
#endif



    static void polyn(const double x, double& p, double& dpdx) {
        // All of the static const stuff is evaluated at compile time

        static const double xmin = 1e-10; // <<<< MINIMUM VALUE OF DENSITY
        static const double xmax = 1e-8;  // <<<< DENSITY SMOOTHLY MODIFIED BELOW THIS VALUE

        static const double xmax2 = xmax*xmax;
        static const double xmax3 = xmax2*xmax;
        static const double xmin2 = xmin*xmin;
        static const double xmin3 = xmin2*xmin;
        static const double r = 1.0/((xmax-xmin)*(-xmin3+(3.0*xmin2+(-3.0*xmin+xmax)*xmax)*xmax));
        static const double a0 = xmax3*xmin*(xmax-4.0*xmin)*r;
        static const double a = xmin2*(xmin2+(-4.0*xmin+18.0*xmax)*xmax)*r;
        static const double b = -6.0*xmin*xmax*(3.0*xmax+2.0*xmin)*r;
        static const double c = (4.0*xmin2+(20.0*xmin+6.0*xmax)*xmax)*r;
        static const double d = -(8.0*xmax+7.0*xmin)*r;
        static const double e = 3.0*r;

        if (x > xmax) {
            p = x;
            dpdx = 1.0;
        }
        else if (x < xmin) {
            p = xmin;
            dpdx = 0.0;
        }
        else {
            p = a0+(a+(b+(c+(d+e*x)*x)*x)*x)*x;
            dpdx = a+(2.0*b+(3.0*c+(4.0*d+5.0*e*x)*x)*x)*x;
        }
    }

    static double munge_old(double rho) {
        double p, dpdx;
        polyn(rho, p, dpdx);
        return p;
    }


    double munge(double rho) const {
        if (rho <= rhotol) rho=rhomin;
        return rho;
    }

    void munge2(double& rho, double& sigma) const {
        if (rho < rhotol) rho=rhomin;
        if (rho < rhotol || sigma < sigtol) sigma=sigmin;
    }

    void munge5(double& rhoa, double& rhob, double& saa, double& sab, double& sbb) const {
        if (rhoa < rhotol || rhob < rhotol || sab < sigtol) sab=sigmin; // ??????????

        if (rhoa < rhotol) rhoa=rhomin;
        if (rhoa < rhotol || saa < sigtol) saa=sigmin;

        if (rhob < rhotol) rhob=rhomin;
        if (rhob < rhotol || sbb < sigtol) sbb=sigmin;
    }

public:
    XCfunctional();


    void initialize(const std::string& input_line, bool polarized);

    ~XCfunctional();

    bool is_lda() const;

    bool is_gga() const;

    bool is_meta() const;

    bool is_dft() const;

    bool is_spin_polarized() const
    {
        return spin_polarized;
    }

    bool has_fxc() const;

    bool has_kxc() const;

    double hf_exchange_coefficient() const
    {
        return hf_coeff;
    }


    madness::Tensor<double> exc(const std::vector< madness::Tensor<double> >& t , const int ispin) const;



    madness::Tensor<double> vxc(const std::vector< madness::Tensor<double> >& t, const int ispin, const int what) const;

    madness::Tensor<double> fxc(const std::vector< madness::Tensor<double> >& t, const int ispin, const int what) const;

    void plot() const {
        long npt = 1001;
        double lo=1e-6, hi=1e+1, s=std::pow(hi/lo, 1.0/(npt-1));

        madness::Tensor<double> rho(npt);
        for (int i=0; i<npt; i++) {
            rho[i] = lo;
            lo *= s;
        }
        std::vector< madness::Tensor<double> > t;
        t.push_back(rho);
        if (is_spin_polarized()) t.push_back(rho);
        madness::Tensor<double> f  = exc(t,0); //pending UGHHHHH
        madness::Tensor<double> va = vxc(t,0,0);
        madness::Tensor<double> vb = vxc(t,0,1);
        for (long i=0; i<npt; i++) {
            printf("%.3e %.3e %.3e %.3e\n", rho[i], f[i], va[i], vb[i]);
        }
    }
};

struct xc_functional {
    const XCfunctional* xc;
    const int ispin;

    xc_functional(const XCfunctional& xc, int ispin)
        : xc(&xc), ispin(ispin)
    {}

    madness::Tensor<double> operator()(const madness::Key<3> & key, const std::vector< madness::Tensor<double> >& t) const
    {
        MADNESS_ASSERT(xc);
        return xc->exc(t,ispin);
    }
};

struct xc_potential {
    const XCfunctional* xc;
    const int what;
    const int ispin;

    xc_potential(const XCfunctional& xc, int ispin,int what)
        : xc(&xc), what(what), ispin(ispin)
    {}

    madness::Tensor<double> operator()(const madness::Key<3> & key, const std::vector< madness::Tensor<double> >& t) const
    {
        MADNESS_ASSERT(xc);
        madness::Tensor<double> r = xc->vxc(t, ispin, what);
        return r;
    }
};

struct xc_kernel {
    const XCfunctional* xc;
    const int what;
    const int ispin;

    xc_kernel(const XCfunctional& xc, int ispin,int what)
        : xc(&xc), what(what), ispin(ispin)
    {}

    madness::Tensor<double> operator()(const madness::Key<3> & key, const std::vector< madness::Tensor<double> >& t) const
    {
        MADNESS_ASSERT(xc);
        madness::Tensor<double> r = xc->fxc(t, ispin, what);
        return r;
    }
};
}

#endif