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:

Generated at Mon Dec 3 23:23:36 2007 for MPQC 2.3.1 using the documentation package Doxygen 1.5.2.