functional.h

大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CHARMM并行计算。非常具有应用价值。

//
// functional.h --- definition of the dft functional
//
// Copyright (C) 1997 Limit Point Systems, Inc.
//
// Author: Curtis Janssen <[e-mail]>
// Maintainer: LPS
//
// This file is part of the SC Toolkit.
//
// The SC Toolkit is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
//
// The SC Toolkit is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//

#ifndef _chemistry_qc_dft_functional_h
#define _chemistry_qc_dft_functional_h

#ifdef __GNUC__
#pragma interface
#endif

#include <util/state/state.h>
#include <math/scmat/vector3.h>
#include <chemistry/qc/wfn/wfn.h>

namespace sc {

/** Contains data needed at each point by a DenFunctional. */
struct PointInputData {
    enum {X=0, Y=1, Z=2};
    enum {XX=0, YX=1, YY=2, ZX=3, ZY=4, ZZ=5};
    struct SpinData {
        double rho;
        // rho^(1/3)
        double rho_13;

        double del_rho[3];
        // gamma = (del rho).(del rho)
        double gamma;

        // hessian of rho
        double hes_rho[6];
        // del^2 rho
        double lap_rho;
    };
    SpinData a, b;

    // gamma_ab = (del rho_a).(del rho_b)
    double gamma_ab;

    const SCVector3 &r;

    // fill in derived quantities
    void compute_derived(int spin_polarized, int need_gradient);

    PointInputData(const SCVector3& r_): r(r_) {}
};

/** Contains data generated at each point by a DenFunctional. */
struct PointOutputData {
    // energy at r
    double energy;

    // derivative of functional wrt density
    double df_drho_a;
    double df_drho_b;

    // derivative of functional wrt density gradient
    double df_dgamma_aa;
    double df_dgamma_bb;
    double df_dgamma_ab;
 
  void zero(){energy=df_drho_a=df_drho_b=df_dgamma_aa=df_dgamma_bb=df_dgamma_ab=0.0;}

};

/** An abstract base class for density functionals. */
class DenFunctional: virtual public SavableState {
  protected:
    int spin_polarized_;
    int compute_potential_;
    double a0_;  // for ACM functionals

    void do_fd_point(PointInputData&id,double&in,double&out,
                     double lower_bound, double upper_bound);
  public:
    DenFunctional();
    DenFunctional(const Ref<KeyVal> &);
    DenFunctional(StateIn &);
    ~DenFunctional();
    void save_data_state(StateOut &);

    // Set to zero if dens_alpha == dens_beta everywhere.
    // The default is false.
    virtual void set_spin_polarized(int i);
    // Set to nonzero if the potential should be computed.
    // The default is false.
    virtual void set_compute_potential(int i);

    // Must return 1 if the density gradient must also be provided.
    // The default implementation returns 0.
    virtual int need_density_gradient();
    // Must return 1 if the density hessian must also be provided.
    // The default implementation returns 0.
    virtual int need_density_hessian();

    virtual void point(const PointInputData&, PointOutputData&) = 0;
    void gradient(const PointInputData&, PointOutputData&,
                  double *gradient, int acenter,
                  GaussianBasisSet *basis,
                  const double *dmat_a, const double *dmat_b,
                  int ncontrib_, const int *contrib_,
                  int ncontrib_bf_, const int *contrib_bf_,
                  const double *bs_values, const double *bsg_values,
                  const double *bsh_values);

    double a0() const { return a0_; }

    void fd_point(const PointInputData&, PointOutputData&);
    int test(const PointInputData &);
    int test();
};


/** The NElFunctional computes the number of electrons.
    It is primarily for testing the integrator. */
class NElFunctional: public DenFunctional {
  public:
    NElFunctional();
    NElFunctional(const Ref<KeyVal> &);
    NElFunctional(StateIn &);
    ~NElFunctional();
    void save_data_state(StateOut &);

    void point(const PointInputData&, PointOutputData&);
};

/** The SumDenFunctional computes energies and densities
    using the a sum of energy density functions method. */
class SumDenFunctional: public DenFunctional {
  protected:
    int n_;
    Ref<DenFunctional> *funcs_;
    double *coefs_;
  public:
    SumDenFunctional();
    SumDenFunctional(const Ref<KeyVal> &);
    SumDenFunctional(StateIn &);
    ~SumDenFunctional();
    void save_data_state(StateOut &);

    void set_spin_polarized(int);
    void set_compute_potential(int);
    int need_density_gradient();

    void point(const PointInputData&, PointOutputData&);

    void print(std::ostream& =ExEnv::out0()) const;
};

/** The StdDenFunctional class is used to construct the standard density
    functionals.

    The table below lists the functional names and the equivalent
    functionals in other packages.  The Name column gives the name as it is
    given in the input file (this is case sensitive).  Functional names
    with non-alpha-numeric names should be given in double quotes.  The
    description column gives the classes used to build up the functional
    and its coefficient, if it is other than one.  The G98 column lists the
    equivalent functional in Gaussian 98 A.6.  The NWChem column lists the
    equivalent functional in NWChem 3.3.1.

<table>
<tr><td>Name    <td> Description         <td> G98    <td> NWChem 
<tr><td>XALPHA  <td> XalphaFunctional    <td> XALPHA <td>        
<tr><td>HFS     <td> SlaterXFunctional   <td> HFS    <td> slater 
<tr><td>HFB     <td> Becke88XFunctional  <td> HFB    <td> becke88
<tr><td>HFG96   <td> G96XFunctional      <td>        <td>        
<tr><td>G96LYP  <td> G96XFunctional                      
                    +LYPCFunctional      <td> G96LYP <td>         
<tr><td>BLYP    <td> SlaterXFunctional                  
                    +Becke88XFunctional                
                    +LYPCFunctional           <td> BLYP  <td>        
<tr><td>SVWN1   <td> SlaterXFunctional
                    +VWN1LCFunctional         <td>       <td> slater vwn_1 
<tr><td>SVWN1RPA<td> SlaterXFunctional
                    +VWN1LCFunctional(1)      <td>       <td> slater vwn_1_rpa 
<tr><td>SVWN2   <td> SlaterXFunctional
                    +VWN2LCFunctional         <td>       <td> slater vwn_2 
<tr><td>SVWN3   <td> SlaterXFunctional
                    +VWN2LCFunctional         <td>       <td> slater vwn_3 
<tr><td>SVWN4   <td> SlaterXFunctional
                    +VWN4LCFunctional         <td>       <td> slater vwn_4 
<tr><td>SVWN5   <td> SlaterXFunctional
                    +VWN5LCFunctional         <td> SVWN5 <td> slater vwn_5 
<tr><td>SPZ81   <td> SlaterXFunctional
                    +PZ81LCFunctional         <td> SPL   <td>         
<tr><td>SPW92   <td> SlaterXFunctional
                    +PW92LCFunctional         <td>       <td> slater pw91lda 
<tr><td>BP86    <td> SlaterXFunctional
                    +Becke88XFunctional
                    +P86CFunctional
                    +PZ81LCFunctional         <td>       <td> becke88 perdue86 
<tr><td>B3LYP   <td> 0.2 HF-Exchange
                    + 0.8  SlaterXFunctional
                    + 0.72 Becke88XFunctional
                    + 0.19 VWN1LCFunctional(1)
                    + 0.81 LYPCFunctional     <td> B3LYP  <td>   b3lyp 
<tr><td>B3PW91  <td> 0.2 HF-Exchange
                    + 0.8  SlaterXFunctional
                    + 0.72 Becke88XFunctional
                    + 0.19 PW91CFunctional
                    + 0.81 PW92LCFunctional   <td> B3PW91 <td>        
<tr><td>B3P86   <td> 0.2 HF-Exchange
                   + 0.8 SlaterXFunctional
                   + 0.72 Becke88XFunctional
                   + 0.19 P86CFunctional
                   + 0.81 VWN1LCFunctional (1)<td>        <td>          
<tr><td>PBE     <td> PBEXFunctional
                    +PBECFunctional           <td>        <td> xpbe96 cpbe96 
<tr><td>PW91    <td> PW91XFunctional
                    +PW91CFunctional          <td>        <td>           
<tr><td>mPW(PW91)PW91
                <td> mPW91XFunctional(PW91)
                    +PW91CFunctional          <td>PW91PW91<td>       
<tr><td>mPWPW91 <td> mPW91XFunctional(mPW91)
                    +PW91CFunctional          <td>        <td>           
<tr><td>mPW1PW91<td> 0.16 HF-Exchange
                   + 0.84 mPW91XFunctional(mPW91)
                   +PW91CFunctional           <td>        <td>           

</table> */
class StdDenFunctional: public SumDenFunctional {
  protected:
    char *name_;
    void init_arrays(int n);
  public:
    StdDenFunctional();
    /** The "name" keyword is read from the input and is used to initialize
        the functional.  All other keywords will be ignored.
    */
    StdDenFunctional(const Ref<KeyVal> &);
    StdDenFunctional(StateIn &);
    ~StdDenFunctional();
    void save_data_state(StateOut &);

    void print(std::ostream& =ExEnv::out0()) const;
};

/** An abstract base class for local correlation functionals. */
class LSDACFunctional: public DenFunctional {
  protected:
  public:
    LSDACFunctional();
    LSDACFunctional(const Ref<KeyVal> &);
    LSDACFunctional(StateIn &);
    ~LSDACFunctional();
    void save_data_state(StateOut &);

    void point(const PointInputData&, PointOutputData&);
    virtual 
      void point_lc(const PointInputData&, PointOutputData&, 
                    double &ec_local, double &decrs, double &deczeta) = 0;

};


/** Implements the Perdew-Burke-Ernzerhof (PBE) correlation functional.

    John P. Perdew, Kieron Burke, and Yue Wang, Phys. Rev. B, 54(23),
    pp. 16533-16539, 1996.

    John P. Perdew, Kieron Burke, and Matthias Ernzerhof, Phys. Rev. Lett.,
    77(18), pp. 3865-3868, 1996.
*/
class PBECFunctional: public DenFunctional {
  protected:
    Ref<LSDACFunctional> local_;
    double gamma;
    double beta;
    void init_constants();
    double rho_deriv(double rho_a, double rho_b, double mdr,
                     double ec_local, double ec_local_dra);
    double gab_deriv(double rho, double phi, double mdr, double ec_local);
  public:
    PBECFunctional();
    PBECFunctional(const Ref<KeyVal> &);
    PBECFunctional(StateIn &);
    ~PBECFunctional();
    void save_data_state(StateOut &);
    int need_density_gradient();
    void point(const PointInputData&, PointOutputData&);
    void set_spin_polarized(int);
  
};

/** The Perdew-Wang 1991 correlation functional computes energies and
    densities using the designated local correlation functional.

    J. P. Perdew, Proceedings of the 75. WE-Heraeus-Seminar and 21st Annual
    International Symposium on Electronic Structure of Solids held in
    Gaussig (Germany), March 11-15, 1991, P. Ziesche and H. Eschrig, eds.,
    pp. 11-20.

    J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson,
    and D. J. Singh, Phys. Rev. B, 46, 6671, 1992.  */
class PW91CFunctional: public DenFunctional {
  protected:
    Ref<LSDACFunctional> local_;
    double a;
    double b;
    double c;
    double d;
    double alpha;
    double c_c0;
    double c_x;
    double nu;
    void init_constants();
    double limit_df_drhoa(double rhoa, double gamma,
                          double ec, double decdrhoa);

  public:
    PW91CFunctional();
    PW91CFunctional(const Ref<KeyVal> &);
    PW91CFunctional(StateIn &);
    ~PW91CFunctional();
    void save_data_state(StateOut &);
    int need_density_gradient();

    void point(const PointInputData&, PointOutputData&);
    void set_spin_polarized(int);
  
};

/** Implements the Perdew 1986 (P86) correlation functional.

    J. P. Perdew, Phys. Rev. B, 33(12), pp. 8822-8824.

    J. P. Perdew, Phys. Rev. B. 34(10), pp. 7406.
 */
class P86CFunctional: public DenFunctional {
  protected:
    double a_;
    double C1_;
    double C2_;
    double C3_;
    double C4_;
    double C5_;
    double C6_;
    double C7_;
    void init_constants();
  public:
    P86CFunctional();
    P86CFunctional(const Ref<KeyVal> &);
    P86CFunctional(StateIn &);
    ~P86CFunctional();
    void save_data_state(StateOut &);
    int need_density_gradient();
    void point(const PointInputData&, PointOutputData&);
  
};


// The Perdew 1986 (P86) Correlation Functional computes energies and densities
//    using the designated local correlation functional.
class NewP86CFunctional: public DenFunctional {
  protected:
    double a_;
    double C1_;
    double C2_;
    double C3_;
    double C4_;
    double C5_;
    double C6_;
    double C7_;
    void init_constants();
    double rho_deriv(double rho_a, double rho_b, double mdr);
    double gab_deriv(double rho_a, double rho_b, double mdr);

  public:
    NewP86CFunctional();
    NewP86CFunctional(const Ref<KeyVal> &);
    NewP86CFunctional(StateIn &);
    ~NewP86CFunctional();
    void save_data_state(StateOut &);
    int need_density_gradient();
    void point(const PointInputData&, PointOutputData&);
};

/**
   Implements the Slater exchange functional.
*/