---

functional.h

//
// functional.h --- definition of the dft functional
//
// Copyright (C) 1997 Limit Point Systems, Inc.
//
// Author: Curtis Janssen <cljanss@limitpt.com>
// Maintainer: LPS
//
// This file is part of the SC Toolkit.
//
// The SC Toolkit is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
//
// The SC Toolkit 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 Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//

#ifndef _chemistry_qc_dft_functional_h
#define _chemistry_qc_dft_functional_h

#ifdef __GNUC__
#pragma interface
#endif

#include <util/state/state.h>
#include <math/scmat/vector3.h>
#include <chemistry/qc/wfn/wfn.h>

struct PointInputData {
    enum {X=0, Y=1, Z=2};
    enum {XX=0, YX=1, YY=2, ZX=3, ZY=4, ZZ=5};
    struct SpinData {
        double rho;
        // rho^(1/3)
        double rho_13;

        double del_rho[3];
        // gamma = (del rho).(del rho)
        double gamma;

        // hessian of rho
        double hes_rho[6];
        // del^2 rho
        double lap_rho;
    };
    SpinData a, b;

    // gamma_ab = (del rho_a).(del rho_b)
    double gamma_ab;

    const SCVector3 &r;

    // fill in derived quantities
    void compute_derived(int spin_polarized, int need_gradient);

    PointInputData(const SCVector3& r_): r(r_) {}
};

struct PointOutputData {
    // energy at r
    double energy;

    // derivative of functional wrt density
    double df_drho_a;
    double df_drho_b;

    // derivative of functional wrt density gradient
    double df_dgamma_aa;
    double df_dgamma_bb;
    double df_dgamma_ab;
 
  void zero(){energy=df_drho_a=df_drho_b=df_dgamma_aa=df_dgamma_bb=df_dgamma_ab=0.0;}

};

class DenFunctional: virtual public SavableState {
  protected:
    int spin_polarized_;
    int compute_potential_;
    double a0_;  // for ACM functionals

    void do_fd_point(PointInputData&id,double&in,double&out,
                     double lower_bound, double upper_bound);
  public:
    DenFunctional();
    DenFunctional(const Ref<KeyVal> &);
    DenFunctional(StateIn &);
    ~DenFunctional();
    void save_data_state(StateOut &);

    // Set to zero if dens_alpha == dens_beta everywhere.
    // The default is false.
    virtual void set_spin_polarized(int i);
    // Set to nonzero if the potential should be computed.
    // The default is false.
    virtual void set_compute_potential(int i);

    // Must return 1 if the density gradient must also be provided.
    // The default implementation returns 0.
    virtual int need_density_gradient();
    // Must return 1 if the density hessian must also be provided.
    // The default implementation returns 0.
    virtual int need_density_hessian();

    virtual void point(const PointInputData&, PointOutputData&) = 0;
    void gradient(const PointInputData&, PointOutputData&,
                  double *gradient, int acenter,
                  GaussianBasisSet *basis,
                  const double *dmat_a, const double *dmat_b,
                  int ncontrib_, const int *contrib_,
                  int ncontrib_bf_, const int *contrib_bf_,
                  const double *bs_values, const double *bsg_values,
                  const double *bsh_values);

    double a0() const { return a0_; }

    void fd_point(const PointInputData&, PointOutputData&);
    int test(const PointInputData &);
    int test();
};


class NElFunctional: public DenFunctional {
  public:
    NElFunctional();
    NElFunctional(const Ref<KeyVal> &);
    NElFunctional(StateIn &);
    ~NElFunctional();
    void save_data_state(StateOut &);

    void point(const PointInputData&, PointOutputData&);
};

class SumDenFunctional: public DenFunctional {
  protected:
    int n_;
    Ref<DenFunctional> *funcs_;
    double *coefs_;
  public:
    SumDenFunctional();
    SumDenFunctional(const Ref<KeyVal> &);
    SumDenFunctional(StateIn &);
    ~SumDenFunctional();
    void save_data_state(StateOut &);

    void set_spin_polarized(int);
    void set_compute_potential(int);
    int need_density_gradient();

    void point(const PointInputData&, PointOutputData&);

    void print(std::ostream& =ExEnv::out()) const;
};

class StdDenFunctional: public SumDenFunctional {
  protected:
    char *name_;
    void init_arrays(int n);
  public:
    StdDenFunctional();
    StdDenFunctional(const Ref<KeyVal> &);
    StdDenFunctional(StateIn &);
    ~StdDenFunctional();
    void save_data_state(StateOut &);

    void print(std::ostream& =ExEnv::out()) const;
};

class LSDACFunctional: public DenFunctional {
  protected:
  public:
    LSDACFunctional();
    LSDACFunctional(const Ref<KeyVal> &);
    LSDACFunctional(StateIn &);
    ~LSDACFunctional();
    void save_data_state(StateOut &);

    void point(const PointInputData&, PointOutputData&);
    virtual 
      void point_lc(const PointInputData&, PointOutputData&, 
                    double &ec_local, double &decrs, double &deczeta) = 0;

};


class PBECFunctional: public DenFunctional {
  protected:
    Ref<LSDACFunctional> local_;
    double gamma;
    double beta;
    void init_constants();
    double rho_deriv(double rho_a, double rho_b, double mdr,
                     double ec_local, double ec_local_dra);
    double gab_deriv(double rho, double phi, double mdr, double ec_local);
  public:
    PBECFunctional();
    PBECFunctional(const Ref<KeyVal> &);
    PBECFunctional(StateIn &);
    ~PBECFunctional();
    void save_data_state(StateOut &);
    int need_density_gradient();
    void point(const PointInputData&, PointOutputData&);
    void set_spin_polarized(int);
  
};

class PW91CFunctional: public DenFunctional {
  protected:
    Ref<LSDACFunctional> local_;
    double a;
    double b;
    double c;
    double d;
    double alpha;
    double c_c0;
    double c_x;
    double nu;
    void init_constants();
    double limit_df_drhoa(double rhoa, double gamma,
                          double ec, double decdrhoa);

  public:
    PW91CFunctional();
    PW91CFunctional(const Ref<KeyVal> &);
    PW91CFunctional(StateIn &);
    ~PW91CFunctional();
    void save_data_state(StateOut &);
    int need_density_gradient();

    void point(const PointInputData&, PointOutputData&);
    void set_spin_polarized(int);
  
};

class P86CFunctional: public DenFunctional {
  protected:
    double a_;
    double C1_;
    double C2_;
    double C3_;
    double C4_;
    double C5_;
    double C6_;
    double C7_;
    void init_constants();
  public:
    P86CFunctional();
    P86CFunctional(const Ref<KeyVal> &);
    P86CFunctional(StateIn &);
    ~P86CFunctional();
    void save_data_state(StateOut &);
    int need_density_gradient();
    void point(const PointInputData&, PointOutputData&);
  
};


// The Perdew 1986 (P86) Correlation Functional computes energies and densities
//    using the designated local correlation functional.
class NewP86CFunctional: public DenFunctional {
  protected:
    double a_;
    double C1_;
    double C2_;
    double C3_;
    double C4_;
    double C5_;
    double C6_;
    double C7_;
    void init_constants();
    double rho_deriv(double rho_a, double rho_b, double mdr);
    double gab_deriv(double rho_a, double rho_b, double mdr);

  public:
    NewP86CFunctional();
    NewP86CFunctional(const Ref<KeyVal> &);
    NewP86CFunctional(StateIn &);
    ~NewP86CFunctional();
    void save_data_state(StateOut &);
    int need_density_gradient();
    void point(const PointInputData&, PointOutputData&);
};

class SlaterXFunctional: public DenFunctional {
  protected:
  public:
    SlaterXFunctional();
    SlaterXFunctional(const Ref<KeyVal> &);
    SlaterXFunctional(StateIn &);
    ~SlaterXFunctional();
    void save_data_state(StateOut &);
    void point(const PointInputData&, PointOutputData&);
};

class VWNLCFunctional: public LSDACFunctional {
  protected:
    double Ap_, Af_, A_alpha_;
    double x0p_mc_, bp_mc_, cp_mc_, x0f_mc_, bf_mc_, cf_mc_;
    double x0p_rpa_, bp_rpa_, cp_rpa_, x0f_rpa_, bf_rpa_, cf_rpa_;
    double x0_alpha_mc_, b_alpha_mc_, c_alpha_mc_;
    double x0_alpha_rpa_, b_alpha_rpa_, c_alpha_rpa_;
    void init_constants();

    double F(double x, double A, double x0, double b, double c);
    double dFdr_s(double x, double A, double x0, double b, double c);
  public:
    VWNLCFunctional();
    VWNLCFunctional(const Ref<KeyVal> &);
    VWNLCFunctional(StateIn &);
    ~VWNLCFunctional();
    void save_data_state(StateOut &);

    virtual
      void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
};
    
class VWN1LCFunctional: public VWNLCFunctional {
  protected:
    double x0p_, bp_, cp_, x0f_, bf_, cf_;
  public:
    VWN1LCFunctional();
    VWN1LCFunctional(int use_rpa);
    VWN1LCFunctional(const Ref<KeyVal> &);
    VWN1LCFunctional(StateIn &);
    ~VWN1LCFunctional();
    void save_data_state(StateOut &);

    void point_lc(const PointInputData&, PointOutputData&,
                  double &, double &, double &);
};

class VWN2LCFunctional: public VWNLCFunctional {
  protected:
  public:
    VWN2LCFunctional();
    VWN2LCFunctional(const Ref<KeyVal> &);
    VWN2LCFunctional(StateIn &);
    ~VWN2LCFunctional();
    void save_data_state(StateOut &);

    void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
};


class VWN3LCFunctional: public VWNLCFunctional {
  protected:
    int monte_carlo_prefactor_;
    int monte_carlo_e0_;
  public:
    VWN3LCFunctional(int mcp = 1, int mce0 = 1);
    VWN3LCFunctional(const Ref<KeyVal> &);
    VWN3LCFunctional(StateIn &);
    ~VWN3LCFunctional();
    void save_data_state(StateOut &);

    void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
};

class VWN4LCFunctional: public VWNLCFunctional {
  protected:
    int monte_carlo_prefactor_;
  public:
    VWN4LCFunctional();
    VWN4LCFunctional(const Ref<KeyVal> &);
    VWN4LCFunctional(StateIn &);
    ~VWN4LCFunctional();
    void save_data_state(StateOut &);

    void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
};

class VWN5LCFunctional: public VWNLCFunctional {
  protected:
  public:
    VWN5LCFunctional();
    VWN5LCFunctional(const Ref<KeyVal> &);
    VWN5LCFunctional(StateIn &);
    ~VWN5LCFunctional();
    void save_data_state(StateOut &);

    void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
};

class PW92LCFunctional: public LSDACFunctional {
  protected:
    double F(double x, double A, double alpha_1, double beta_1, double beta_2, 
             double beta_3, double beta_4, double p);
    double dFdr_s(double x, double A, double alpha_1, double beta_1, double beta_2, 
             double beta_3, double beta_4, double p);
  public:
    PW92LCFunctional();
    PW92LCFunctional(const Ref<KeyVal> &);
    PW92LCFunctional(StateIn &);
    ~PW92LCFunctional();
    void save_data_state(StateOut &);

    void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
};

class PZ81LCFunctional: public LSDACFunctional {
  protected:
    double Fec_rsgt1(double rs, double beta_1, double beta_2, double gamma);
    double dFec_rsgt1_drho(double rs, double beta_1, double beta_2, double gamma,
                           double &dec_drs);
    double Fec_rslt1(double rs, double A, double B, double C, double D);
    double dFec_rslt1_drho(double rs, double A, double B, double C, double D,
                           double &dec_drs);
  public:
    PZ81LCFunctional();
    PZ81LCFunctional(const Ref<KeyVal> &);
    PZ81LCFunctional(StateIn &);
    ~PZ81LCFunctional();
    void save_data_state(StateOut &);

    void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
};

class XalphaFunctional: public DenFunctional {
  protected:
    double alpha_;
    double factor_;
  public:
    XalphaFunctional();
    XalphaFunctional(const Ref<KeyVal> &);
    XalphaFunctional(StateIn &);
    ~XalphaFunctional();
    void save_data_state(StateOut &);

    void point(const PointInputData&, PointOutputData&);

    void print(std::ostream& =ExEnv::out()) const;
};

class Becke88XFunctional: public DenFunctional {
  protected:
    double beta_;
    double beta6_;
  public:
    Becke88XFunctional();
    Becke88XFunctional(const Ref<KeyVal> &);
    Becke88XFunctional(StateIn &);
    ~Becke88XFunctional();
    void save_data_state(StateOut &);

    int need_density_gradient();

    void point(const PointInputData&, PointOutputData&);
};

class LYPCFunctional: public DenFunctional {
  protected:
    double a_;
    double b_;
    double c_;
    double d_;
    void init_constants();
  public:
    LYPCFunctional();
    LYPCFunctional(const Ref<KeyVal> &);
    LYPCFunctional(StateIn &);
    ~LYPCFunctional();
    void save_data_state(StateOut &);

    int need_density_gradient();

    void point(const PointInputData&, PointOutputData&);
};

class PW86XFunctional: public DenFunctional {
  protected:
    double a_;
    double b_;
    double c_;
    double m_;
    void init_constants();
  public:
    PW86XFunctional();
    PW86XFunctional(const Ref<KeyVal> &);
    PW86XFunctional(StateIn &);
    ~PW86XFunctional();
    void save_data_state(StateOut &);

    int need_density_gradient();

    void point(const PointInputData&, PointOutputData&);
};

class PBEXFunctional: public DenFunctional {
  protected:
    double mu;
    double kappa;
    void spin_contrib(const PointInputData::SpinData &,
                      double &mpw, double &dmpw_dr, double &dmpw_dg);
    void init_constants();
  public:
    PBEXFunctional();
    PBEXFunctional(const Ref<KeyVal> &);
    PBEXFunctional(StateIn &);
    ~PBEXFunctional();
    void save_data_state(StateOut &);

    int need_density_gradient();

    void point(const PointInputData&, PointOutputData&);
};

class PW91XFunctional: public DenFunctional {
  protected:
    double a;
    double b;
    double c;
    double d;
    double a_x;
    void spin_contrib(const PointInputData::SpinData &,
                      double &mpw, double &dmpw_dr, double &dmpw_dg);
    void init_constants();
  public:
    PW91XFunctional();
    PW91XFunctional(const Ref<KeyVal> &);
    PW91XFunctional(StateIn &);
    ~PW91XFunctional();
    void save_data_state(StateOut &);

    int need_density_gradient();

    void point(const PointInputData&, PointOutputData&);
};

class mPW91XFunctional: public DenFunctional {
  protected:
    double b;
    double beta;
    double c;
    double d;
    double a_x;
    double x_d_coef;

    void spin_contrib(const PointInputData::SpinData &,
                      double &mpw, double &dmpw_dr, double &dmpw_dg);
  public:
    enum Func { B88, PW91, mPW91 };

    mPW91XFunctional();
    mPW91XFunctional(Func variant);
    mPW91XFunctional(const Ref<KeyVal> &);
    mPW91XFunctional(StateIn &);
    ~mPW91XFunctional();
    void save_data_state(StateOut &);

    int need_density_gradient();

    void point(const PointInputData&, PointOutputData&);

    void init_constants(Func);
};

class G96XFunctional: public DenFunctional {
  protected:
    double b_;
    void init_constants();
  public:
    G96XFunctional();
    G96XFunctional(const Ref<KeyVal> &);
    G96XFunctional(StateIn &);
    ~G96XFunctional();
    void save_data_state(StateOut &);

    int need_density_gradient();

    void point(const PointInputData&, PointOutputData&);
};

#endif

// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End:

---