MPQC 2.3.1
|
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: