functional.h

00001 //
00002 // functional.h --- definition of the dft functional
00003 //
00004 // Copyright (C) 1997 Limit Point Systems, Inc.
00005 //
00006 // Author: Curtis Janssen <cljanss@limitpt.com>
00007 // Maintainer: LPS
00008 //
00009 // This file is part of the SC Toolkit.
00010 //
00011 // The SC Toolkit is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Library General Public License as published by
00013 // the Free Software Foundation; either version 2, or (at your option)
00014 // any later version.
00015 //
00016 // The SC Toolkit is distributed in the hope that it will be useful,
00017 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00019 // GNU Library General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Library General Public License
00022 // along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
00023 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
00024 //
00025 // The U.S. Government is granted a limited license as per AL 91-7.
00026 //
00027 
00028 #ifndef _chemistry_qc_dft_functional_h
00029 #define _chemistry_qc_dft_functional_h
00030 
00031 #ifdef __GNUC__
00032 #pragma interface
00033 #endif
00034 
00035 #include <util/state/state.h>
00036 #include <math/scmat/vector3.h>
00037 #include <chemistry/qc/wfn/wfn.h>
00038 #include <chemistry/qc/wfn/density.h>
00039 
00040 namespace sc {
00041 
00043 struct PointInputData {
00044     enum {X=BatchElectronDensity::X,
00045           Y=BatchElectronDensity::Y,
00046           Z=BatchElectronDensity::Z};
00047     enum {XX=BatchElectronDensity::XX,
00048           YX=BatchElectronDensity::YX,
00049           YY=BatchElectronDensity::YY,
00050           ZX=BatchElectronDensity::ZX,
00051           ZY=BatchElectronDensity::ZY,
00052           ZZ=BatchElectronDensity::ZZ};
00053     struct SpinData {
00054         double rho;
00055         // rho^(1/3)
00056         double rho_13;
00057 
00058         double del_rho[3];
00059         // gamma = (del rho).(del rho)
00060         double gamma;
00061 
00062         // hessian of rho
00063         double hes_rho[6];
00064         // del^2 rho
00065         double lap_rho;
00066     };
00067     SpinData a, b;
00068 
00069     // gamma_ab = (del rho_a).(del rho_b)
00070     double gamma_ab;
00071 
00072     const SCVector3 &r;
00073 
00074     // fill in derived quantities
00075     void compute_derived(int spin_polarized,
00076                          int need_gradient,
00077                          int need_hessian);
00078 
00079     PointInputData(const SCVector3& r_): r(r_) {}
00080 };
00081 
00083 struct PointOutputData {
00084     // energy at r
00085     double energy;
00086 
00087     // derivative of functional wrt density
00088     double df_drho_a;
00089     double df_drho_b;
00090 
00091     // derivative of functional wrt density gradient
00092     double df_dgamma_aa;
00093     double df_dgamma_bb;
00094     double df_dgamma_ab;
00095  
00096   void zero(){energy=df_drho_a=df_drho_b=df_dgamma_aa=df_dgamma_bb=df_dgamma_ab=0.0;}
00097 
00098 };
00099 
00101 class DenFunctional: virtual public SavableState {
00102   protected:
00103     int spin_polarized_;
00104     int compute_potential_;
00105     double a0_;  // for ACM functionals
00106 
00107     void do_fd_point(PointInputData&id,double&in,double&out,
00108                      double lower_bound, double upper_bound);
00109   public:
00110     DenFunctional();
00111     DenFunctional(const Ref<KeyVal> &);
00112     DenFunctional(StateIn &);
00113     ~DenFunctional();
00114     void save_data_state(StateOut &);
00115 
00116     // Set to zero if dens_alpha == dens_beta everywhere.
00117     // The default is false.
00118     virtual void set_spin_polarized(int i);
00119     // Set to nonzero if the potential should be computed.
00120     // The default is false.
00121     virtual void set_compute_potential(int i);
00122 
00123     // Must return 1 if the density gradient must also be provided.
00124     // The default implementation returns 0.
00125     virtual int need_density_gradient();
00126     // Must return 1 if the density hessian must also be provided.
00127     // The default implementation returns 0.
00128     virtual int need_density_hessian();
00129 
00130     virtual void point(const PointInputData&, PointOutputData&) = 0;
00131     void gradient(const PointInputData&, PointOutputData&,
00132                   double *gradient, int acenter,
00133                   GaussianBasisSet *basis,
00134                   const double *dmat_a, const double *dmat_b,
00135                   int ncontrib, const int *contrib,
00136                   int ncontrib_bf, const int *contrib_bf,
00137                   const double *bs_values, const double *bsg_values,
00138                   const double *bsh_values);
00139 
00140     double a0() const { return a0_; }
00141 
00142     void fd_point(const PointInputData&, PointOutputData&);
00143     int test(const PointInputData &);
00144     int test();
00145 };
00146 
00147 
00150 class NElFunctional: public DenFunctional {
00151   public:
00152     NElFunctional();
00153     NElFunctional(const Ref<KeyVal> &);
00154     NElFunctional(StateIn &);
00155     ~NElFunctional();
00156     void save_data_state(StateOut &);
00157 
00158     void point(const PointInputData&, PointOutputData&);
00159 };
00160 
00163 class SumDenFunctional: public DenFunctional {
00164   protected:
00165     int n_;
00166     Ref<DenFunctional> *funcs_;
00167     double *coefs_;
00168   public:
00169     SumDenFunctional();
00197     SumDenFunctional(const Ref<KeyVal> &);
00198     SumDenFunctional(StateIn &);
00199     ~SumDenFunctional();
00200     void save_data_state(StateOut &);
00201 
00202     void set_spin_polarized(int);
00203     void set_compute_potential(int);
00204     int need_density_gradient();
00205 
00206     void point(const PointInputData&, PointOutputData&);
00207 
00208     void print(std::ostream& =ExEnv::out0()) const;
00209 };
00210 
00283 class StdDenFunctional: public SumDenFunctional {
00284   protected:
00285     char *name_;
00286     void init_arrays(int n);
00287   public:
00288     StdDenFunctional();
00293     StdDenFunctional(const Ref<KeyVal> &);
00294     StdDenFunctional(StateIn &);
00295     ~StdDenFunctional();
00296     void save_data_state(StateOut &);
00297 
00298     void print(std::ostream& =ExEnv::out0()) const;
00299 };
00300 
00302 class LSDACFunctional: public DenFunctional {
00303   protected:
00304   public:
00305     LSDACFunctional();
00306     LSDACFunctional(const Ref<KeyVal> &);
00307     LSDACFunctional(StateIn &);
00308     ~LSDACFunctional();
00309     void save_data_state(StateOut &);
00310 
00311     void point(const PointInputData&, PointOutputData&);
00312     virtual 
00313       void point_lc(const PointInputData&, PointOutputData&, 
00314                     double &ec_local, double &decrs, double &deczeta) = 0;
00315 
00316 };
00317 
00318 
00327 class PBECFunctional: public DenFunctional {
00328   protected:
00329     Ref<LSDACFunctional> local_;
00330     double gamma;
00331     double beta;
00332     void init_constants();
00333     double rho_deriv(double rho_a, double rho_b, double mdr,
00334                      double ec_local, double ec_local_dra);
00335     double gab_deriv(double rho, double phi, double mdr, double ec_local);
00336   public:
00337     PBECFunctional();
00338     PBECFunctional(const Ref<KeyVal> &);
00339     PBECFunctional(StateIn &);
00340     ~PBECFunctional();
00341     void save_data_state(StateOut &);
00342     int need_density_gradient();
00343     void point(const PointInputData&, PointOutputData&);
00344     void set_spin_polarized(int);
00345   
00346 };
00347 
00358 class PW91CFunctional: public DenFunctional {
00359   protected:
00360     Ref<LSDACFunctional> local_;
00361     double a;
00362     double b;
00363     double c;
00364     double d;
00365     double alpha;
00366     double c_c0;
00367     double c_x;
00368     double nu;
00369     void init_constants();
00370     double limit_df_drhoa(double rhoa, double gamma,
00371                           double ec, double decdrhoa);
00372 
00373   public:
00374     PW91CFunctional();
00375     PW91CFunctional(const Ref<KeyVal> &);
00376     PW91CFunctional(StateIn &);
00377     ~PW91CFunctional();
00378     void save_data_state(StateOut &);
00379     int need_density_gradient();
00380 
00381     void point(const PointInputData&, PointOutputData&);
00382     void set_spin_polarized(int);
00383   
00384 };
00385 
00392 class P86CFunctional: public DenFunctional {
00393   protected:
00394     double a_;
00395     double C1_;
00396     double C2_;
00397     double C3_;
00398     double C4_;
00399     double C5_;
00400     double C6_;
00401     double C7_;
00402     void init_constants();
00403   public:
00404     P86CFunctional();
00405     P86CFunctional(const Ref<KeyVal> &);
00406     P86CFunctional(StateIn &);
00407     ~P86CFunctional();
00408     void save_data_state(StateOut &);
00409     int need_density_gradient();
00410     void point(const PointInputData&, PointOutputData&);
00411   
00412 };
00413 
00414 
00415 // The Perdew 1986 (P86) Correlation Functional computes energies and densities
00416 //    using the designated local correlation functional.
00417 class NewP86CFunctional: public DenFunctional {
00418   protected:
00419     double a_;
00420     double C1_;
00421     double C2_;
00422     double C3_;
00423     double C4_;
00424     double C5_;
00425     double C6_;
00426     double C7_;
00427     void init_constants();
00428     double rho_deriv(double rho_a, double rho_b, double mdr);
00429     double gab_deriv(double rho_a, double rho_b, double mdr);
00430 
00431   public:
00432     NewP86CFunctional();
00433     NewP86CFunctional(const Ref<KeyVal> &);
00434     NewP86CFunctional(StateIn &);
00435     ~NewP86CFunctional();
00436     void save_data_state(StateOut &);
00437     int need_density_gradient();
00438     void point(const PointInputData&, PointOutputData&);
00439 };
00440 
00444 class SlaterXFunctional: public DenFunctional {
00445   protected:
00446   public:
00447     SlaterXFunctional();
00448     SlaterXFunctional(const Ref<KeyVal> &);
00449     SlaterXFunctional(StateIn &);
00450     ~SlaterXFunctional();
00451     void save_data_state(StateOut &);
00452     void point(const PointInputData&, PointOutputData&);
00453 };
00454 
00462 class VWNLCFunctional: public LSDACFunctional {
00463   protected:
00464     double Ap_, Af_, A_alpha_;
00465     double x0p_mc_, bp_mc_, cp_mc_, x0f_mc_, bf_mc_, cf_mc_;
00466     double x0p_rpa_, bp_rpa_, cp_rpa_, x0f_rpa_, bf_rpa_, cf_rpa_;
00467     double x0_alpha_mc_, b_alpha_mc_, c_alpha_mc_;
00468     double x0_alpha_rpa_, b_alpha_rpa_, c_alpha_rpa_;
00469     void init_constants();
00470 
00471     double F(double x, double A, double x0, double b, double c);
00472     double dFdr_s(double x, double A, double x0, double b, double c);
00473   public:
00474     VWNLCFunctional();
00475     VWNLCFunctional(const Ref<KeyVal> &);
00476     VWNLCFunctional(StateIn &);
00477     ~VWNLCFunctional();
00478     void save_data_state(StateOut &);
00479 
00480     virtual
00481       void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
00482 };
00483     
00486 class VWN1LCFunctional: public VWNLCFunctional {
00487   protected:
00488     double x0p_, bp_, cp_, x0f_, bf_, cf_;
00489   public:
00491     VWN1LCFunctional();
00493     VWN1LCFunctional(int use_rpa);
00499     VWN1LCFunctional(const Ref<KeyVal> &);
00500     VWN1LCFunctional(StateIn &);
00501     ~VWN1LCFunctional();
00502     void save_data_state(StateOut &);
00503 
00504     void point_lc(const PointInputData&, PointOutputData&,
00505                   double &, double &, double &);
00506 };
00507 
00510 class VWN2LCFunctional: public VWNLCFunctional {
00511   protected:
00512   public:
00514     VWN2LCFunctional();
00516     VWN2LCFunctional(const Ref<KeyVal> &);
00517     VWN2LCFunctional(StateIn &);
00518     ~VWN2LCFunctional();
00519     void save_data_state(StateOut &);
00520 
00521     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
00522 };
00523 
00524 
00527 class VWN3LCFunctional: public VWNLCFunctional {
00528   protected:
00529     int monte_carlo_prefactor_;
00530     int monte_carlo_e0_;
00531   public:
00532     VWN3LCFunctional(int mcp = 1, int mce0 = 1);
00533     VWN3LCFunctional(const Ref<KeyVal> &);
00534     VWN3LCFunctional(StateIn &);
00535     ~VWN3LCFunctional();
00536     void save_data_state(StateOut &);
00537 
00538     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
00539 };
00540 
00543 class VWN4LCFunctional: public VWNLCFunctional {
00544   protected:
00545     int monte_carlo_prefactor_;
00546   public:
00547     VWN4LCFunctional();
00548     VWN4LCFunctional(const Ref<KeyVal> &);
00549     VWN4LCFunctional(StateIn &);
00550     ~VWN4LCFunctional();
00551     void save_data_state(StateOut &);
00552 
00553     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
00554 };
00555 
00558 class VWN5LCFunctional: public VWNLCFunctional {
00559   protected:
00560   public:
00561     VWN5LCFunctional();
00562     VWN5LCFunctional(const Ref<KeyVal> &);
00563     VWN5LCFunctional(StateIn &);
00564     ~VWN5LCFunctional();
00565     void save_data_state(StateOut &);
00566 
00567     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
00568 };
00569 
00575 class PW92LCFunctional: public LSDACFunctional {
00576   protected:
00577     double F(double x, double A, double alpha_1, double beta_1, double beta_2, 
00578              double beta_3, double beta_4, double p);
00579     double dFdr_s(double x, double A, double alpha_1, double beta_1, double beta_2, 
00580              double beta_3, double beta_4, double p);
00581   public:
00582     PW92LCFunctional();
00583     PW92LCFunctional(const Ref<KeyVal> &);
00584     PW92LCFunctional(StateIn &);
00585     ~PW92LCFunctional();
00586     void save_data_state(StateOut &);
00587 
00588     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
00589 };
00590 
00596 class PZ81LCFunctional: public LSDACFunctional {
00597   protected:
00598     double Fec_rsgt1(double rs, double beta_1, double beta_2, double gamma);
00599     double dFec_rsgt1_drho(double rs, double beta_1, double beta_2, double gamma,
00600                            double &dec_drs);
00601     double Fec_rslt1(double rs, double A, double B, double C, double D);
00602     double dFec_rslt1_drho(double rs, double A, double B, double C, double D,
00603                            double &dec_drs);
00604   public:
00605     PZ81LCFunctional();
00606     PZ81LCFunctional(const Ref<KeyVal> &);
00607     PZ81LCFunctional(StateIn &);
00608     ~PZ81LCFunctional();
00609     void save_data_state(StateOut &);
00610 
00611     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
00612 };
00613 
00615 class XalphaFunctional: public DenFunctional {
00616   protected:
00617     double alpha_;
00618     double factor_;
00619   public:
00620     XalphaFunctional();
00621     XalphaFunctional(const Ref<KeyVal> &);
00622     XalphaFunctional(StateIn &);
00623     ~XalphaFunctional();
00624     void save_data_state(StateOut &);
00625 
00626     void point(const PointInputData&, PointOutputData&);
00627 
00628     void print(std::ostream& =ExEnv::out0()) const;
00629 };
00630 
00635 class Becke88XFunctional: public DenFunctional {
00636   protected:
00637     double beta_;
00638     double beta6_;
00639   public:
00640     Becke88XFunctional();
00641     Becke88XFunctional(const Ref<KeyVal> &);
00642     Becke88XFunctional(StateIn &);
00643     ~Becke88XFunctional();
00644     void save_data_state(StateOut &);
00645 
00646     int need_density_gradient();
00647 
00648     void point(const PointInputData&, PointOutputData&);
00649 };
00650 
00659 class LYPCFunctional: public DenFunctional {
00660   protected:
00661     double a_;
00662     double b_;
00663     double c_;
00664     double d_;
00665     void init_constants();
00666   public:
00667     LYPCFunctional();
00668     LYPCFunctional(const Ref<KeyVal> &);
00669     LYPCFunctional(StateIn &);
00670     ~LYPCFunctional();
00671     void save_data_state(StateOut &);
00672 
00673     int need_density_gradient();
00674 
00675     void point(const PointInputData&, PointOutputData&);
00676 };
00677 
00682 class PW86XFunctional: public DenFunctional {
00683   protected:
00684     double a_;
00685     double b_;
00686     double c_;
00687     double m_;
00688     void init_constants();
00689   public:
00690     PW86XFunctional();
00691     PW86XFunctional(const Ref<KeyVal> &);
00692     PW86XFunctional(StateIn &);
00693     ~PW86XFunctional();
00694     void save_data_state(StateOut &);
00695 
00696     int need_density_gradient();
00697 
00698     void point(const PointInputData&, PointOutputData&);
00699 };
00700 
00716 class PBEXFunctional: public DenFunctional {
00717   protected:
00718     double mu;
00719     double kappa;
00720     void spin_contrib(const PointInputData::SpinData &,
00721                       double &mpw, double &dmpw_dr, double &dmpw_dg);
00722     void init_constants();
00723   public:
00724     PBEXFunctional();
00725     PBEXFunctional(const Ref<KeyVal> &);
00726     PBEXFunctional(StateIn &);
00727     ~PBEXFunctional();
00728     void save_data_state(StateOut &);
00729 
00730     int need_density_gradient();
00731 
00732     void point(const PointInputData&, PointOutputData&);
00733 };
00734 
00745 class PW91XFunctional: public DenFunctional {
00746   protected:
00747     double a;
00748     double b;
00749     double c;
00750     double d;
00751     double a_x;
00752     void spin_contrib(const PointInputData::SpinData &,
00753                       double &mpw, double &dmpw_dr, double &dmpw_dg);
00754     void init_constants();
00755   public:
00756     PW91XFunctional();
00757     PW91XFunctional(const Ref<KeyVal> &);
00758     PW91XFunctional(StateIn &);
00759     ~PW91XFunctional();
00760     void save_data_state(StateOut &);
00761 
00762     int need_density_gradient();
00763 
00764     void point(const PointInputData&, PointOutputData&);
00765 };
00766 
00771 class mPW91XFunctional: public DenFunctional {
00772   protected:
00773     double b;
00774     double beta;
00775     double c;
00776     double d;
00777     double a_x;
00778     double x_d_coef;
00779 
00780     void spin_contrib(const PointInputData::SpinData &,
00781                       double &mpw, double &dmpw_dr, double &dmpw_dg);
00782   public:
00783     enum Func { B88, PW91, mPW91 };
00784 
00786     mPW91XFunctional();
00789     mPW91XFunctional(Func variant);
00808     mPW91XFunctional(const Ref<KeyVal> &);
00809     mPW91XFunctional(StateIn &);
00810     ~mPW91XFunctional();
00811     void save_data_state(StateOut &);
00812 
00813     int need_density_gradient();
00814 
00815     void point(const PointInputData&, PointOutputData&);
00816 
00817     void init_constants(Func);
00818 };
00819 
00824 class G96XFunctional: public DenFunctional {
00825   protected:
00826     double b_;
00827     void init_constants();
00828   public:
00829     G96XFunctional();
00830     G96XFunctional(const Ref<KeyVal> &);
00831     G96XFunctional(StateIn &);
00832     ~G96XFunctional();
00833     void save_data_state(StateOut &);
00834 
00835     int need_density_gradient();
00836 
00837     void point(const PointInputData&, PointOutputData&);
00838 };
00839 
00840 }
00841 
00842 #endif
00843 
00844 // Local Variables:
00845 // mode: c++
00846 // c-file-style: "CLJ"
00847 // End:

Generated at Sun Dec 10 18:20:43 2006 for MPQC 2.3.0 using the documentation package Doxygen 1.4.7.