# 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
    00141     virtual double a0() const;
    00142
    00143     void fd_point(const PointInputData&, PointOutputData&);
    00144     int test(const PointInputData &);
    00145     int test();
    00146 };
    00147
    00148
    00151 class NElFunctional: public DenFunctional {
    00152   public:
    00153     NElFunctional();
    00154     NElFunctional(const Ref<KeyVal> &);
    00155     NElFunctional(StateIn &);
    00156     ~NElFunctional();
    00157     void save_data_state(StateOut &);
    00158
    00159     void point(const PointInputData&, PointOutputData&);
    00160 };
    00161
    00164 class SumDenFunctional: public DenFunctional {
    00165   protected:
    00166     int n_;
    00167     Ref<DenFunctional> *funcs_;
    00168     double *coefs_;
    00169   public:
    00170     SumDenFunctional();
    00198     SumDenFunctional(const Ref<KeyVal> &);
    00199     SumDenFunctional(StateIn &);
    00200     ~SumDenFunctional();
    00201     void save_data_state(StateOut &);
    00202
    00203     void set_spin_polarized(int);
    00204     void set_compute_potential(int);
    00205     int need_density_gradient();
    00206
    00207     void point(const PointInputData&, PointOutputData&);
    00208
    00209     void print(std::ostream& =ExEnv::out0()) const;
    00210
    00213     double a0() const;
    00214 };
    00215
    00288 class StdDenFunctional: public SumDenFunctional {
    00289   protected:
    00290     char *name_;
    00291     void init_arrays(int n);
    00292   public:
    00293     StdDenFunctional();
    00298     StdDenFunctional(const Ref<KeyVal> &);
    00299     StdDenFunctional(StateIn &);
    00300     ~StdDenFunctional();
    00301     void save_data_state(StateOut &);
    00302
    00303     void print(std::ostream& =ExEnv::out0()) const;
    00304 };
    00305
    00307 class LSDACFunctional: public DenFunctional {
    00308   protected:
    00309   public:
    00310     LSDACFunctional();
    00311     LSDACFunctional(const Ref<KeyVal> &);
    00312     LSDACFunctional(StateIn &);
    00313     ~LSDACFunctional();
    00314     void save_data_state(StateOut &);
    00315
    00316     void point(const PointInputData&, PointOutputData&);
    00317     virtual
    00318       void point_lc(const PointInputData&, PointOutputData&,
    00319                     double &ec_local, double &decrs, double &deczeta) = 0;
    00320
    00321 };
    00322
    00323
    00332 class PBECFunctional: public DenFunctional {
    00333   protected:
    00334     Ref<LSDACFunctional> local_;
    00335     double gamma;
    00336     double beta;
    00337     void init_constants();
    00338     double rho_deriv(double rho_a, double rho_b, double mdr,
    00339                      double ec_local, double ec_local_dra);
    00340     double gab_deriv(double rho, double phi, double mdr, double ec_local);
    00341   public:
    00342     PBECFunctional();
    00343     PBECFunctional(const Ref<KeyVal> &);
    00344     PBECFunctional(StateIn &);
    00345     ~PBECFunctional();
    00346     void save_data_state(StateOut &);
    00347     int need_density_gradient();
    00348     void point(const PointInputData&, PointOutputData&);
    00349     void set_spin_polarized(int);
    00350
    00351 };
    00352
    00363 class PW91CFunctional: public DenFunctional {
    00364   protected:
    00365     Ref<LSDACFunctional> local_;
    00366     double a;
    00367     double b;
    00368     double c;
    00369     double d;
    00370     double alpha;
    00371     double c_c0;
    00372     double c_x;
    00373     double nu;
    00374     void init_constants();
    00375     double limit_df_drhoa(double rhoa, double gamma,
    00376                           double ec, double decdrhoa);
    00377
    00378   public:
    00379     PW91CFunctional();
    00380     PW91CFunctional(const Ref<KeyVal> &);
    00381     PW91CFunctional(StateIn &);
    00382     ~PW91CFunctional();
    00383     void save_data_state(StateOut &);
    00384     int need_density_gradient();
    00385
    00386     void point(const PointInputData&, PointOutputData&);
    00387     void set_spin_polarized(int);
    00388
    00389 };
    00390
    00397 class P86CFunctional: public DenFunctional {
    00398   protected:
    00399     double a_;
    00400     double C1_;
    00401     double C2_;
    00402     double C3_;
    00403     double C4_;
    00404     double C5_;
    00405     double C6_;
    00406     double C7_;
    00407     void init_constants();
    00408   public:
    00409     P86CFunctional();
    00410     P86CFunctional(const Ref<KeyVal> &);
    00411     P86CFunctional(StateIn &);
    00412     ~P86CFunctional();
    00413     void save_data_state(StateOut &);
    00414     int need_density_gradient();
    00415     void point(const PointInputData&, PointOutputData&);
    00416
    00417 };
    00418
    00419
    00420 // The Perdew 1986 (P86) Correlation Functional computes energies and densities
    00421 //    using the designated local correlation functional.
    00422 class NewP86CFunctional: public DenFunctional {
    00423   protected:
    00424     double a_;
    00425     double C1_;
    00426     double C2_;
    00427     double C3_;
    00428     double C4_;
    00429     double C5_;
    00430     double C6_;
    00431     double C7_;
    00432     void init_constants();
    00433     double rho_deriv(double rho_a, double rho_b, double mdr);
    00434     double gab_deriv(double rho_a, double rho_b, double mdr);
    00435
    00436   public:
    00437     NewP86CFunctional();
    00438     NewP86CFunctional(const Ref<KeyVal> &);
    00439     NewP86CFunctional(StateIn &);
    00440     ~NewP86CFunctional();
    00441     void save_data_state(StateOut &);
    00442     int need_density_gradient();
    00443     void point(const PointInputData&, PointOutputData&);
    00444 };
    00445
    00449 class SlaterXFunctional: public DenFunctional {
    00450   protected:
    00451   public:
    00452     SlaterXFunctional();
    00453     SlaterXFunctional(const Ref<KeyVal> &);
    00454     SlaterXFunctional(StateIn &);
    00455     ~SlaterXFunctional();
    00456     void save_data_state(StateOut &);
    00457     void point(const PointInputData&, PointOutputData&);
    00458 };
    00459
    00467 class VWNLCFunctional: public LSDACFunctional {
    00468   protected:
    00469     double Ap_, Af_, A_alpha_;
    00470     double x0p_mc_, bp_mc_, cp_mc_, x0f_mc_, bf_mc_, cf_mc_;
    00471     double x0p_rpa_, bp_rpa_, cp_rpa_, x0f_rpa_, bf_rpa_, cf_rpa_;
    00472     double x0_alpha_mc_, b_alpha_mc_, c_alpha_mc_;
    00473     double x0_alpha_rpa_, b_alpha_rpa_, c_alpha_rpa_;
    00474     void init_constants();
    00475
    00476     double F(double x, double A, double x0, double b, double c);
    00477     double dFdr_s(double x, double A, double x0, double b, double c);
    00478   public:
    00479     VWNLCFunctional();
    00480     VWNLCFunctional(const Ref<KeyVal> &);
    00481     VWNLCFunctional(StateIn &);
    00482     ~VWNLCFunctional();
    00483     void save_data_state(StateOut &);
    00484
    00485     virtual
    00486       void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
    00487 };
    00488
    00491 class VWN1LCFunctional: public VWNLCFunctional {
    00492   protected:
    00493     double x0p_, bp_, cp_, x0f_, bf_, cf_;
    00494   public:
    00496     VWN1LCFunctional();
    00498     VWN1LCFunctional(int use_rpa);
    00504     VWN1LCFunctional(const Ref<KeyVal> &);
    00505     VWN1LCFunctional(StateIn &);
    00506     ~VWN1LCFunctional();
    00507     void save_data_state(StateOut &);
    00508
    00509     void point_lc(const PointInputData&, PointOutputData&,
    00510                   double &, double &, double &);
    00511 };
    00512
    00515 class VWN2LCFunctional: public VWNLCFunctional {
    00516   protected:
    00517   public:
    00519     VWN2LCFunctional();
    00521     VWN2LCFunctional(const Ref<KeyVal> &);
    00522     VWN2LCFunctional(StateIn &);
    00523     ~VWN2LCFunctional();
    00524     void save_data_state(StateOut &);
    00525
    00526     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
    00527 };
    00528
    00529
    00532 class VWN3LCFunctional: public VWNLCFunctional {
    00533   protected:
    00534     int monte_carlo_prefactor_;
    00535     int monte_carlo_e0_;
    00536   public:
    00537     VWN3LCFunctional(int mcp = 1, int mce0 = 1);
    00538     VWN3LCFunctional(const Ref<KeyVal> &);
    00539     VWN3LCFunctional(StateIn &);
    00540     ~VWN3LCFunctional();
    00541     void save_data_state(StateOut &);
    00542
    00543     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
    00544 };
    00545
    00548 class VWN4LCFunctional: public VWNLCFunctional {
    00549   protected:
    00550     int monte_carlo_prefactor_;
    00551   public:
    00552     VWN4LCFunctional();
    00553     VWN4LCFunctional(const Ref<KeyVal> &);
    00554     VWN4LCFunctional(StateIn &);
    00555     ~VWN4LCFunctional();
    00556     void save_data_state(StateOut &);
    00557
    00558     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
    00559 };
    00560
    00563 class VWN5LCFunctional: public VWNLCFunctional {
    00564   protected:
    00565   public:
    00566     VWN5LCFunctional();
    00567     VWN5LCFunctional(const Ref<KeyVal> &);
    00568     VWN5LCFunctional(StateIn &);
    00569     ~VWN5LCFunctional();
    00570     void save_data_state(StateOut &);
    00571
    00572     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
    00573 };
    00574
    00580 class PW92LCFunctional: public LSDACFunctional {
    00581   protected:
    00582     double F(double x, double A, double alpha_1, double beta_1, double beta_2,
    00583              double beta_3, double beta_4, double p);
    00584     double dFdr_s(double x, double A, double alpha_1, double beta_1, double beta_2,
    00585              double beta_3, double beta_4, double p);
    00586   public:
    00587     PW92LCFunctional();
    00588     PW92LCFunctional(const Ref<KeyVal> &);
    00589     PW92LCFunctional(StateIn &);
    00590     ~PW92LCFunctional();
    00591     void save_data_state(StateOut &);
    00592
    00593     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
    00594 };
    00595
    00601 class PZ81LCFunctional: public LSDACFunctional {
    00602   protected:
    00603     double Fec_rsgt1(double rs, double beta_1, double beta_2, double gamma);
    00604     double dFec_rsgt1_drho(double rs, double beta_1, double beta_2, double gamma,
    00605                            double &dec_drs);
    00606     double Fec_rslt1(double rs, double A, double B, double C, double D);
    00607     double dFec_rslt1_drho(double rs, double A, double B, double C, double D,
    00608                            double &dec_drs);
    00609   public:
    00610     PZ81LCFunctional();
    00611     PZ81LCFunctional(const Ref<KeyVal> &);
    00612     PZ81LCFunctional(StateIn &);
    00613     ~PZ81LCFunctional();
    00614     void save_data_state(StateOut &);
    00615
    00616     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
    00617 };
    00618
    00620 class XalphaFunctional: public DenFunctional {
    00621   protected:
    00622     double alpha_;
    00623     double factor_;
    00624   public:
    00625     XalphaFunctional();
    00626     XalphaFunctional(const Ref<KeyVal> &);
    00627     XalphaFunctional(StateIn &);
    00628     ~XalphaFunctional();
    00629     void save_data_state(StateOut &);
    00630
    00631     void point(const PointInputData&, PointOutputData&);
    00632
    00633     void print(std::ostream& =ExEnv::out0()) const;
    00634 };
    00635
    00640 class Becke88XFunctional: public DenFunctional {
    00641   protected:
    00642     double beta_;
    00643     double beta6_;
    00644   public:
    00645     Becke88XFunctional();
    00646     Becke88XFunctional(const Ref<KeyVal> &);
    00647     Becke88XFunctional(StateIn &);
    00648     ~Becke88XFunctional();
    00649     void save_data_state(StateOut &);
    00650
    00651     int need_density_gradient();
    00652
    00653     void point(const PointInputData&, PointOutputData&);
    00654 };
    00655
    00664 class LYPCFunctional: public DenFunctional {
    00665   protected:
    00666     double a_;
    00667     double b_;
    00668     double c_;
    00669     double d_;
    00670     void init_constants();
    00671   public:
    00672     LYPCFunctional();
    00673     LYPCFunctional(const Ref<KeyVal> &);
    00674     LYPCFunctional(StateIn &);
    00675     ~LYPCFunctional();
    00676     void save_data_state(StateOut &);
    00677
    00678     int need_density_gradient();
    00679
    00680     void point(const PointInputData&, PointOutputData&);
    00681 };
    00682
    00687 class PW86XFunctional: public DenFunctional {
    00688   protected:
    00689     double a_;
    00690     double b_;
    00691     double c_;
    00692     double m_;
    00693     void init_constants();
    00694   public:
    00695     PW86XFunctional();
    00696     PW86XFunctional(const Ref<KeyVal> &);
    00697     PW86XFunctional(StateIn &);
    00698     ~PW86XFunctional();
    00699     void save_data_state(StateOut &);
    00700
    00701     int need_density_gradient();
    00702
    00703     void point(const PointInputData&, PointOutputData&);
    00704 };
    00705
    00721 class PBEXFunctional: public DenFunctional {
    00722   protected:
    00723     double mu;
    00724     double kappa;
    00725     void spin_contrib(const PointInputData::SpinData &,
    00726                       double &mpw, double &dmpw_dr, double &dmpw_dg);
    00727     void init_constants();
    00728   public:
    00729     PBEXFunctional();
    00730     PBEXFunctional(const Ref<KeyVal> &);
    00731     PBEXFunctional(StateIn &);
    00732     ~PBEXFunctional();
    00733     void save_data_state(StateOut &);
    00734
    00735     int need_density_gradient();
    00736
    00737     void point(const PointInputData&, PointOutputData&);
    00738 };
    00739
    00750 class PW91XFunctional: public DenFunctional {
    00751   protected:
    00752     double a;
    00753     double b;
    00754     double c;
    00755     double d;
    00756     double a_x;
    00757     void spin_contrib(const PointInputData::SpinData &,
    00758                       double &mpw, double &dmpw_dr, double &dmpw_dg);
    00759     void init_constants();
    00760   public:
    00761     PW91XFunctional();
    00762     PW91XFunctional(const Ref<KeyVal> &);
    00763     PW91XFunctional(StateIn &);
    00764     ~PW91XFunctional();
    00765     void save_data_state(StateOut &);
    00766
    00767     int need_density_gradient();
    00768
    00769     void point(const PointInputData&, PointOutputData&);
    00770 };
    00771
    00776 class mPW91XFunctional: public DenFunctional {
    00777   protected:
    00778     double b;
    00779     double beta;
    00780     double c;
    00781     double d;
    00782     double a_x;
    00783     double x_d_coef;
    00784
    00785     void spin_contrib(const PointInputData::SpinData &,
    00786                       double &mpw, double &dmpw_dr, double &dmpw_dg);
    00787   public:
    00788     enum Func { B88, PW91, mPW91 };
    00789
    00791     mPW91XFunctional();
    00794     mPW91XFunctional(Func variant);
    00813     mPW91XFunctional(const Ref<KeyVal> &);
    00814     mPW91XFunctional(StateIn &);
    00815     ~mPW91XFunctional();
    00816     void save_data_state(StateOut &);
    00817
    00818     int need_density_gradient();
    00819
    00820     void point(const PointInputData&, PointOutputData&);
    00821
    00822     void init_constants(Func);
    00823 };
    00824
    00829 class G96XFunctional: public DenFunctional {
    00830   protected:
    00831     double b_;
    00832     void init_constants();
    00833   public:
    00834     G96XFunctional();
    00835     G96XFunctional(const Ref<KeyVal> &);
    00836     G96XFunctional(StateIn &);
    00837     ~G96XFunctional();
    00838     void save_data_state(StateOut &);
    00839
    00840     int need_density_gradient();
    00841
    00842     void point(const PointInputData&, PointOutputData&);
    00843 };
    00844
    00845 }
    00846
    00847 #endif
    00848
    00849 // Local Variables:
    00850 // mode: c++
    00851 // c-file-style: "CLJ"
    00852 // End: