Main Page Class Hierarchy Compound List File List Compound Members Related Pages

# 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
    00039 namespace sc {
    00040
    00042 struct PointInputData {
    00043     enum {X=0, Y=1, Z=2};
    00044     enum {XX=0, YX=1, YY=2, ZX=3, ZY=4, ZZ=5};
    00045     struct SpinData {
    00046         double rho;
    00047         // rho^(1/3)
    00048         double rho_13;
    00049
    00050         double del_rho[3];
    00051         // gamma = (del rho).(del rho)
    00052         double gamma;
    00053
    00054         // hessian of rho
    00055         double hes_rho[6];
    00056         // del^2 rho
    00057         double lap_rho;
    00058     };
    00059     SpinData a, b;
    00060
    00061     // gamma_ab = (del rho_a).(del rho_b)
    00062     double gamma_ab;
    00063
    00064     const SCVector3 &r;
    00065
    00066     // fill in derived quantities
    00067     void compute_derived(int spin_polarized, int need_gradient);
    00068
    00069     PointInputData(const SCVector3& r_): r(r_) {}
    00070 };
    00071
    00073 struct PointOutputData {
    00074     // energy at r
    00075     double energy;
    00076
    00077     // derivative of functional wrt density
    00078     double df_drho_a;
    00079     double df_drho_b;
    00080
    00081     // derivative of functional wrt density gradient
    00082     double df_dgamma_aa;
    00083     double df_dgamma_bb;
    00084     double df_dgamma_ab;
    00085
    00086   void zero(){energy=df_drho_a=df_drho_b=df_dgamma_aa=df_dgamma_bb=df_dgamma_ab=0.0;}
    00087
    00088 };
    00089
    00091 class DenFunctional: virtual public SavableState {
    00092   protected:
    00093     int spin_polarized_;
    00094     int compute_potential_;
    00095     double a0_;  // for ACM functionals
    00096
    00097     void do_fd_point(PointInputData&id,double&in,double&out,
    00098                      double lower_bound, double upper_bound);
    00099   public:
    00100     DenFunctional();
    00101     DenFunctional(const Ref<KeyVal> &);
    00102     DenFunctional(StateIn &);
    00103     ~DenFunctional();
    00104     void save_data_state(StateOut &);
    00105
    00106     // Set to zero if dens_alpha == dens_beta everywhere.
    00107     // The default is false.
    00108     virtual void set_spin_polarized(int i);
    00109     // Set to nonzero if the potential should be computed.
    00110     // The default is false.
    00111     virtual void set_compute_potential(int i);
    00112
    00113     // Must return 1 if the density gradient must also be provided.
    00114     // The default implementation returns 0.
    00115     virtual int need_density_gradient();
    00116     // Must return 1 if the density hessian must also be provided.
    00117     // The default implementation returns 0.
    00118     virtual int need_density_hessian();
    00119
    00120     virtual void point(const PointInputData&, PointOutputData&) = 0;
    00121     void gradient(const PointInputData&, PointOutputData&,
    00122                   double *gradient, int acenter,
    00123                   GaussianBasisSet *basis,
    00124                   const double *dmat_a, const double *dmat_b,
    00125                   int ncontrib_, const int *contrib_,
    00126                   int ncontrib_bf_, const int *contrib_bf_,
    00127                   const double *bs_values, const double *bsg_values,
    00128                   const double *bsh_values);
    00129
    00130     double a0() const { return a0_; }
    00131
    00132     void fd_point(const PointInputData&, PointOutputData&);
    00133     int test(const PointInputData &);
    00134     int test();
    00135 };
    00136
    00137
    00140 class NElFunctional: public DenFunctional {
    00141   public:
    00142     NElFunctional();
    00143     NElFunctional(const Ref<KeyVal> &);
    00144     NElFunctional(StateIn &);
    00145     ~NElFunctional();
    00146     void save_data_state(StateOut &);
    00147
    00148     void point(const PointInputData&, PointOutputData&);
    00149 };
    00150
    00153 class SumDenFunctional: public DenFunctional {
    00154   protected:
    00155     int n_;
    00156     Ref<DenFunctional> *funcs_;
    00157     double *coefs_;
    00158   public:
    00159     SumDenFunctional();
    00160     SumDenFunctional(const Ref<KeyVal> &);
    00161     SumDenFunctional(StateIn &);
    00162     ~SumDenFunctional();
    00163     void save_data_state(StateOut &);
    00164
    00165     void set_spin_polarized(int);
    00166     void set_compute_potential(int);
    00167     int need_density_gradient();
    00168
    00169     void point(const PointInputData&, PointOutputData&);
    00170
    00171     void print(std::ostream& =ExEnv::out0()) const;
    00172 };
    00173
    00246 class StdDenFunctional: public SumDenFunctional {
    00247   protected:
    00248     char *name_;
    00249     void init_arrays(int n);
    00250   public:
    00251     StdDenFunctional();
    00255     StdDenFunctional(const Ref<KeyVal> &);
    00256     StdDenFunctional(StateIn &);
    00257     ~StdDenFunctional();
    00258     void save_data_state(StateOut &);
    00259
    00260     void print(std::ostream& =ExEnv::out0()) const;
    00261 };
    00262
    00264 class LSDACFunctional: public DenFunctional {
    00265   protected:
    00266   public:
    00267     LSDACFunctional();
    00268     LSDACFunctional(const Ref<KeyVal> &);
    00269     LSDACFunctional(StateIn &);
    00270     ~LSDACFunctional();
    00271     void save_data_state(StateOut &);
    00272
    00273     void point(const PointInputData&, PointOutputData&);
    00274     virtual
    00275       void point_lc(const PointInputData&, PointOutputData&,
    00276                     double &ec_local, double &decrs, double &deczeta) = 0;
    00277
    00278 };
    00279
    00280
    00289 class PBECFunctional: public DenFunctional {
    00290   protected:
    00291     Ref<LSDACFunctional> local_;
    00292     double gamma;
    00293     double beta;
    00294     void init_constants();
    00295     double rho_deriv(double rho_a, double rho_b, double mdr,
    00296                      double ec_local, double ec_local_dra);
    00297     double gab_deriv(double rho, double phi, double mdr, double ec_local);
    00298   public:
    00299     PBECFunctional();
    00300     PBECFunctional(const Ref<KeyVal> &);
    00301     PBECFunctional(StateIn &);
    00302     ~PBECFunctional();
    00303     void save_data_state(StateOut &);
    00304     int need_density_gradient();
    00305     void point(const PointInputData&, PointOutputData&);
    00306     void set_spin_polarized(int);
    00307
    00308 };
    00309
    00320 class PW91CFunctional: public DenFunctional {
    00321   protected:
    00322     Ref<LSDACFunctional> local_;
    00323     double a;
    00324     double b;
    00325     double c;
    00326     double d;
    00327     double alpha;
    00328     double c_c0;
    00329     double c_x;
    00330     double nu;
    00331     void init_constants();
    00332     double limit_df_drhoa(double rhoa, double gamma,
    00333                           double ec, double decdrhoa);
    00334
    00335   public:
    00336     PW91CFunctional();
    00337     PW91CFunctional(const Ref<KeyVal> &);
    00338     PW91CFunctional(StateIn &);
    00339     ~PW91CFunctional();
    00340     void save_data_state(StateOut &);
    00341     int need_density_gradient();
    00342
    00343     void point(const PointInputData&, PointOutputData&);
    00344     void set_spin_polarized(int);
    00345
    00346 };
    00347
    00354 class P86CFunctional: public DenFunctional {
    00355   protected:
    00356     double a_;
    00357     double C1_;
    00358     double C2_;
    00359     double C3_;
    00360     double C4_;
    00361     double C5_;
    00362     double C6_;
    00363     double C7_;
    00364     void init_constants();
    00365   public:
    00366     P86CFunctional();
    00367     P86CFunctional(const Ref<KeyVal> &);
    00368     P86CFunctional(StateIn &);
    00369     ~P86CFunctional();
    00370     void save_data_state(StateOut &);
    00371     int need_density_gradient();
    00372     void point(const PointInputData&, PointOutputData&);
    00373
    00374 };
    00375
    00376
    00377 // The Perdew 1986 (P86) Correlation Functional computes energies and densities
    00378 //    using the designated local correlation functional.
    00379 class NewP86CFunctional: public DenFunctional {
    00380   protected:
    00381     double a_;
    00382     double C1_;
    00383     double C2_;
    00384     double C3_;
    00385     double C4_;
    00386     double C5_;
    00387     double C6_;
    00388     double C7_;
    00389     void init_constants();
    00390     double rho_deriv(double rho_a, double rho_b, double mdr);
    00391     double gab_deriv(double rho_a, double rho_b, double mdr);
    00392
    00393   public:
    00394     NewP86CFunctional();
    00395     NewP86CFunctional(const Ref<KeyVal> &);
    00396     NewP86CFunctional(StateIn &);
    00397     ~NewP86CFunctional();
    00398     void save_data_state(StateOut &);
    00399     int need_density_gradient();
    00400     void point(const PointInputData&, PointOutputData&);
    00401 };
    00402
    00406 class SlaterXFunctional: public DenFunctional {
    00407   protected:
    00408   public:
    00409     SlaterXFunctional();
    00410     SlaterXFunctional(const Ref<KeyVal> &);
    00411     SlaterXFunctional(StateIn &);
    00412     ~SlaterXFunctional();
    00413     void save_data_state(StateOut &);
    00414     void point(const PointInputData&, PointOutputData&);
    00415 };
    00416
    00424 class VWNLCFunctional: public LSDACFunctional {
    00425   protected:
    00426     double Ap_, Af_, A_alpha_;
    00427     double x0p_mc_, bp_mc_, cp_mc_, x0f_mc_, bf_mc_, cf_mc_;
    00428     double x0p_rpa_, bp_rpa_, cp_rpa_, x0f_rpa_, bf_rpa_, cf_rpa_;
    00429     double x0_alpha_mc_, b_alpha_mc_, c_alpha_mc_;
    00430     double x0_alpha_rpa_, b_alpha_rpa_, c_alpha_rpa_;
    00431     void init_constants();
    00432
    00433     double F(double x, double A, double x0, double b, double c);
    00434     double dFdr_s(double x, double A, double x0, double b, double c);
    00435   public:
    00436     VWNLCFunctional();
    00437     VWNLCFunctional(const Ref<KeyVal> &);
    00438     VWNLCFunctional(StateIn &);
    00439     ~VWNLCFunctional();
    00440     void save_data_state(StateOut &);
    00441
    00442     virtual
    00443       void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
    00444 };
    00445
    00448 class VWN1LCFunctional: public VWNLCFunctional {
    00449   protected:
    00450     double x0p_, bp_, cp_, x0f_, bf_, cf_;
    00451   public:
    00453     VWN1LCFunctional();
    00455     VWN1LCFunctional(int use_rpa);
    00461     VWN1LCFunctional(const Ref<KeyVal> &);
    00462     VWN1LCFunctional(StateIn &);
    00463     ~VWN1LCFunctional();
    00464     void save_data_state(StateOut &);
    00465
    00466     void point_lc(const PointInputData&, PointOutputData&,
    00467                   double &, double &, double &);
    00468 };
    00469
    00472 class VWN2LCFunctional: public VWNLCFunctional {
    00473   protected:
    00474   public:
    00476     VWN2LCFunctional();
    00478     VWN2LCFunctional(const Ref<KeyVal> &);
    00479     VWN2LCFunctional(StateIn &);
    00480     ~VWN2LCFunctional();
    00481     void save_data_state(StateOut &);
    00482
    00483     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
    00484 };
    00485
    00486
    00489 class VWN3LCFunctional: public VWNLCFunctional {
    00490   protected:
    00491     int monte_carlo_prefactor_;
    00492     int monte_carlo_e0_;
    00493   public:
    00494     VWN3LCFunctional(int mcp = 1, int mce0 = 1);
    00495     VWN3LCFunctional(const Ref<KeyVal> &);
    00496     VWN3LCFunctional(StateIn &);
    00497     ~VWN3LCFunctional();
    00498     void save_data_state(StateOut &);
    00499
    00500     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
    00501 };
    00502
    00505 class VWN4LCFunctional: public VWNLCFunctional {
    00506   protected:
    00507     int monte_carlo_prefactor_;
    00508   public:
    00509     VWN4LCFunctional();
    00510     VWN4LCFunctional(const Ref<KeyVal> &);
    00511     VWN4LCFunctional(StateIn &);
    00512     ~VWN4LCFunctional();
    00513     void save_data_state(StateOut &);
    00514
    00515     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
    00516 };
    00517
    00520 class VWN5LCFunctional: public VWNLCFunctional {
    00521   protected:
    00522   public:
    00523     VWN5LCFunctional();
    00524     VWN5LCFunctional(const Ref<KeyVal> &);
    00525     VWN5LCFunctional(StateIn &);
    00526     ~VWN5LCFunctional();
    00527     void save_data_state(StateOut &);
    00528
    00529     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
    00530 };
    00531
    00537 class PW92LCFunctional: public LSDACFunctional {
    00538   protected:
    00539     double F(double x, double A, double alpha_1, double beta_1, double beta_2,
    00540              double beta_3, double beta_4, double p);
    00541     double dFdr_s(double x, double A, double alpha_1, double beta_1, double beta_2,
    00542              double beta_3, double beta_4, double p);
    00543   public:
    00544     PW92LCFunctional();
    00545     PW92LCFunctional(const Ref<KeyVal> &);
    00546     PW92LCFunctional(StateIn &);
    00547     ~PW92LCFunctional();
    00548     void save_data_state(StateOut &);
    00549
    00550     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
    00551 };
    00552
    00558 class PZ81LCFunctional: public LSDACFunctional {
    00559   protected:
    00560     double Fec_rsgt1(double rs, double beta_1, double beta_2, double gamma);
    00561     double dFec_rsgt1_drho(double rs, double beta_1, double beta_2, double gamma,
    00562                            double &dec_drs);
    00563     double Fec_rslt1(double rs, double A, double B, double C, double D);
    00564     double dFec_rslt1_drho(double rs, double A, double B, double C, double D,
    00565                            double &dec_drs);
    00566   public:
    00567     PZ81LCFunctional();
    00568     PZ81LCFunctional(const Ref<KeyVal> &);
    00569     PZ81LCFunctional(StateIn &);
    00570     ~PZ81LCFunctional();
    00571     void save_data_state(StateOut &);
    00572
    00573     void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
    00574 };
    00575
    00577 class XalphaFunctional: public DenFunctional {
    00578   protected:
    00579     double alpha_;
    00580     double factor_;
    00581   public:
    00582     XalphaFunctional();
    00583     XalphaFunctional(const Ref<KeyVal> &);
    00584     XalphaFunctional(StateIn &);
    00585     ~XalphaFunctional();
    00586     void save_data_state(StateOut &);
    00587
    00588     void point(const PointInputData&, PointOutputData&);
    00589
    00590     void print(std::ostream& =ExEnv::out0()) const;
    00591 };
    00592
    00597 class Becke88XFunctional: public DenFunctional {
    00598   protected:
    00599     double beta_;
    00600     double beta6_;
    00601   public:
    00602     Becke88XFunctional();
    00603     Becke88XFunctional(const Ref<KeyVal> &);
    00604     Becke88XFunctional(StateIn &);
    00605     ~Becke88XFunctional();
    00606     void save_data_state(StateOut &);
    00607
    00608     int need_density_gradient();
    00609
    00610     void point(const PointInputData&, PointOutputData&);
    00611 };
    00612
    00621 class LYPCFunctional: public DenFunctional {
    00622   protected:
    00623     double a_;
    00624     double b_;
    00625     double c_;
    00626     double d_;
    00627     void init_constants();
    00628   public:
    00629     LYPCFunctional();
    00630     LYPCFunctional(const Ref<KeyVal> &);
    00631     LYPCFunctional(StateIn &);
    00632     ~LYPCFunctional();
    00633     void save_data_state(StateOut &);
    00634
    00635     int need_density_gradient();
    00636
    00637     void point(const PointInputData&, PointOutputData&);
    00638 };
    00639
    00644 class PW86XFunctional: public DenFunctional {
    00645   protected:
    00646     double a_;
    00647     double b_;
    00648     double c_;
    00649     double m_;
    00650     void init_constants();
    00651   public:
    00652     PW86XFunctional();
    00653     PW86XFunctional(const Ref<KeyVal> &);
    00654     PW86XFunctional(StateIn &);
    00655     ~PW86XFunctional();
    00656     void save_data_state(StateOut &);
    00657
    00658     int need_density_gradient();
    00659
    00660     void point(const PointInputData&, PointOutputData&);
    00661 };
    00662
    00678 class PBEXFunctional: public DenFunctional {
    00679   protected:
    00680     double mu;
    00681     double kappa;
    00682     void spin_contrib(const PointInputData::SpinData &,
    00683                       double &mpw, double &dmpw_dr, double &dmpw_dg);
    00684     void init_constants();
    00685   public:
    00686     PBEXFunctional();
    00687     PBEXFunctional(const Ref<KeyVal> &);
    00688     PBEXFunctional(StateIn &);
    00689     ~PBEXFunctional();
    00690     void save_data_state(StateOut &);
    00691
    00692     int need_density_gradient();
    00693
    00694     void point(const PointInputData&, PointOutputData&);
    00695 };
    00696
    00707 class PW91XFunctional: public DenFunctional {
    00708   protected:
    00709     double a;
    00710     double b;
    00711     double c;
    00712     double d;
    00713     double a_x;
    00714     void spin_contrib(const PointInputData::SpinData &,
    00715                       double &mpw, double &dmpw_dr, double &dmpw_dg);
    00716     void init_constants();
    00717   public:
    00718     PW91XFunctional();
    00719     PW91XFunctional(const Ref<KeyVal> &);
    00720     PW91XFunctional(StateIn &);
    00721     ~PW91XFunctional();
    00722     void save_data_state(StateOut &);
    00723
    00724     int need_density_gradient();
    00725
    00726     void point(const PointInputData&, PointOutputData&);
    00727 };
    00728
    00733 class mPW91XFunctional: public DenFunctional {
    00734   protected:
    00735     double b;
    00736     double beta;
    00737     double c;
    00738     double d;
    00739     double a_x;
    00740     double x_d_coef;
    00741
    00742     void spin_contrib(const PointInputData::SpinData &,
    00743                       double &mpw, double &dmpw_dr, double &dmpw_dg);
    00744   public:
    00745     enum Func { B88, PW91, mPW91 };
    00746
    00748     mPW91XFunctional();
    00751     mPW91XFunctional(Func variant);
    00770     mPW91XFunctional(const Ref<KeyVal> &);
    00771     mPW91XFunctional(StateIn &);
    00772     ~mPW91XFunctional();
    00773     void save_data_state(StateOut &);
    00774
    00775     int need_density_gradient();
    00776
    00777     void point(const PointInputData&, PointOutputData&);
    00778
    00779     void init_constants(Func);
    00780 };
    00781
    00786 class G96XFunctional: public DenFunctional {
    00787   protected:
    00788     double b_;
    00789     void init_constants();
    00790   public:
    00791     G96XFunctional();
    00792     G96XFunctional(const Ref<KeyVal> &);
    00793     G96XFunctional(StateIn &);
    00794     ~G96XFunctional();
    00795     void save_data_state(StateOut &);
    00796
    00797     int need_density_gradient();
    00798
    00799     void point(const PointInputData&, PointOutputData&);
    00800 };
    00801
    00802 }
    00803
    00804 #endif
    00805
    00806 // Local Variables:
    00807 // mode: c++
    00808 // c-file-style: "CLJ"
    00809 // End: