DFcode/xcfunctional.h

Go to the documentation of this file.

#ifndef MOLDFT_XCMOLDFT_H
#define MOLDFT_XCMOLDFT_H

#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

struct xc_lda_potential {
    xc_lda_potential() {}

    void operator()(const madness::Key<3> & key, madness::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;
        }
    }
};

#ifdef MADNESS_HAS_MADXC
//#error
struct xc_func_type {
  int   number;   /* indentifier number */                                                                                                   
  bool nspin;      /* XC_UNPOLARIZED or XC_POLARIZED  */
  int  kind;     /* XC_EXCHANGE, XC_CORRELATION, or XC_EXCHANGE_CORRELATION */                                                              
  std::string name;     /* name of the functional, e.g. "PBE" */                                                                             
  int   family;   /* type of the functional, e.g. XC_FAMILY_GGA */                                                                           
  double   coef;   /* coefficient */
};  
#endif

class XCfunctional {
protected:
    bool spin_polarized;        
    double hf_coeff;            
#if defined(MADNESS_HAS_MADXC) ||defined( MADNESS_HAS_LIBXC)
    std::vector< std::pair<xc_func_type*,double> > funcs;
    int nderiv;
#endif

#ifdef MADNESS_HAS_LIBXC
    void make_libxc_args(const std::vector< madness::Tensor<double> >& t,
                         madness::Tensor<double>& rho,
                         madness::Tensor<double>& sigma,
                         madness::Tensor<double>& delrho) const;
#endif

#ifdef MADNESS_HAS_MADXC
    void make_xc_args(const std::vector< madness::Tensor<double> >& t,
                         std::vector< madness::Tensor<double> >& rho,
                         std::vector< madness::Tensor<double> >& sigma,
                         std::vector< madness::Tensor<double> >& delrho) const;
#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(double rho) {
        double p, dpdx;
        polyn(rho, p, dpdx);
        return p;
    }

    static void munge2(double& rho, double& sigma) {
        // rho(x) --> p(rho(x))
        // d/dx p(rho(x)) --> dp/drho * drho/dx
        if (sigma < 0.0) sigma = 0.0;
        double p, dpdx;
        polyn(rho, p, dpdx);
        rho = p;
        sigma *= dpdx*dpdx;
    }
    static void munge_der(double& rhoa, double& sigma, double& drx, double& dry, double& drz) {
        // rho(x) --> p(rho(x))
        // d/dx p(rho(x)) --> dp/drho * drho/dx
        if (sigma < 0.0) sigma = 0.0;
        double p, dpdx;
        polyn(rhoa, p, dpdx);
        rhoa = p;
        sigma *= dpdx*dpdx;
        drx *=dpdx;
        dry *=dpdx;
        drz *=dpdx;
    }
    static void munge5(double& rhoa, double& rhob, double& saa, double& sab, double& sbb) {
        // rho(x) --> p(rho(x))
        // d/dx p(rho(x)) --> dp/drho * drho/dx
        if (saa < 0.0) saa = 0.0;
        if (sab < 0.0) sab = 0.0;
        if (sbb < 0.0) sbb = 0.0;
        double pa, pb, dpadx, dpbdx;
        polyn(rhoa, pa, dpadx);
        polyn(rhob, pb, dpbdx);
        rhoa = pa;
        rhob = pb;
        saa *= dpadx*dpadx;
        sab *= dpadx*dpbdx;
        sbb *= dpbdx*dpbdx;
    }
    static void munge5_der(double& rhoa, double& rhob, double& saa, double& sab, double& sbb,
                           double& drax, double& dray, double& draz,
                           double& drbx, double& drby, double& drbz) {
        // rho(x) --> p(rho(x))
        // d/dx p(rho(x)) --> dp/drho * drho/dx
        if (saa < 0.0) saa = 0.0;
        if (sab < 0.0) sab = 0.0;
        if (sbb < 0.0) sbb = 0.0;
        double pa, pb, dpadx, dpbdx;
        polyn(rhoa, pa, dpadx);
        polyn(rhob, pb, dpbdx);
        rhoa = pa;
        rhob = pb;
        saa *= dpadx*dpadx;
        sab *= dpadx*dpbdx;
        sbb *= dpbdx*dpbdx;

        drax *=dpadx;
        dray *=dpadx;
        draz *=dpadx;

        drbx *=dpbdx;
        drby *=dpbdx;
        drbz *=dpbdx;
    }

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;

    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;
    }
};

#endif