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>
#include <chemistry/qc/wfn/density.h>
namespace sc {
struct PointInputData {
enum {X=BatchElectronDensity::X,
Y=BatchElectronDensity::Y,
Z=BatchElectronDensity::Z};
enum {XX=BatchElectronDensity::XX,
YX=BatchElectronDensity::YX,
YY=BatchElectronDensity::YY,
ZX=BatchElectronDensity::ZX,
ZY=BatchElectronDensity::ZY,
ZZ=BatchElectronDensity::ZZ};
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,
int need_hessian);
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);
virtual double a0() const;
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::out0()) const;
double a0() 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::out0()) 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::out0()) 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: