functional.cc

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

    << indent << "Standard Density Functional: " << n << endl;
  SumDenFunctional::print(o);
}

/////////////////////////////////////////////////////////////////////////////
// LSDACFunctional: All local correlation functionals inherit from this class.
// Coded by Matt Leininger
static ClassDesc LSDACFunctional_cd(
  typeid(LSDACFunctional),"LSDACFunctional",1,"public DenFunctional",
  0, 0, 0);

LSDACFunctional::LSDACFunctional(StateIn& s):
  SavableState(s),
  DenFunctional(s)
{
}

LSDACFunctional::LSDACFunctional()
{
}

LSDACFunctional::LSDACFunctional(const Ref<KeyVal>& keyval):
  DenFunctional(keyval)
{
}

LSDACFunctional::~LSDACFunctional()
{
}

void
LSDACFunctional::save_data_state(StateOut& s)
{
  DenFunctional::save_data_state(s);
}

void
LSDACFunctional::point(const PointInputData &id,
                       PointOutputData &od)
{
  double junk_1, junk_2, junk_3;
  point_lc(id, od, junk_1, junk_2, junk_3);
}

/////////////////////////////////////////////////////////////////////////////
// SlaterXFunctional

static ClassDesc SlaterXFunctional_cd(
  typeid(SlaterXFunctional),"SlaterXFunctional",1,"public DenFunctional",
  0, create<SlaterXFunctional>, create<SlaterXFunctional>);

SlaterXFunctional::SlaterXFunctional(StateIn& s):
  SavableState(s),
  DenFunctional(s)
{
}

SlaterXFunctional::SlaterXFunctional()
{
}

SlaterXFunctional::SlaterXFunctional(const Ref<KeyVal>& keyval):
  DenFunctional(keyval)
{
}

SlaterXFunctional::~SlaterXFunctional()
{
}

void
SlaterXFunctional::save_data_state(StateOut& s)
{
  DenFunctional::save_data_state(s);
}

void
SlaterXFunctional::point(const PointInputData &id,
                       PointOutputData &od)
{
  const double mcx2rthird = -0.9305257363491; // -1.5*(3/4pi)^1/3
  const double dmcx2rthird = -1.2407009817988; // 2*(3/4pi)^1/3
  od.zero();

  if (!spin_polarized_) {
      od.energy = mcx2rthird * 2.0 * id.a.rho * id.a.rho_13;
      if (compute_potential_) {
          od.df_drho_a = dmcx2rthird * id.a.rho_13;
          od.df_drho_b = od.df_drho_a;
          }
    }
  else {
      od.energy = mcx2rthird
                * (id.a.rho * id.a.rho_13 + id.b.rho * id.b.rho_13);
      if (compute_potential_) {
          od.df_drho_a = dmcx2rthird * id.a.rho_13;
          od.df_drho_b = dmcx2rthird * id.b.rho_13;
          }
    }
}

/////////////////////////////////////////////////////////////////////////////
// PW92LCFunctional
// Coded by Matt Leininger
static ClassDesc PW92LCFunctional_cd(
  typeid(PW92LCFunctional),"PW92LCFunctional",1,"public LSDACFunctional",
  0, create<PW92LCFunctional>, create<PW92LCFunctional>);

PW92LCFunctional::PW92LCFunctional(StateIn& s):
  SavableState(s),
  LSDACFunctional(s)
{
}

PW92LCFunctional::PW92LCFunctional()
{
}

PW92LCFunctional::PW92LCFunctional(const Ref<KeyVal>& keyval):
  LSDACFunctional(keyval)
{
}

PW92LCFunctional::~PW92LCFunctional()
{
}

void
PW92LCFunctional::save_data_state(StateOut& s)
{
  LSDACFunctional::save_data_state(s);
}

double
PW92LCFunctional::F(double x, double A, double alpha_1, double beta_1, double beta_2,
                    double beta_3, double beta_4, double p)
{
  double x2 = x*x; // r_s
  double denom = 2.*A*( beta_1 * x + beta_2 * x2 + beta_3 * x2*x + beta_4 * pow(x2,p+1.));
  double res = -2.*A*(1. + alpha_1*x2)*log(1.+ 1./denom);

  return res;
}

double
PW92LCFunctional::dFdr_s(double x, double A, double alpha_1, double beta_1, double beta_2,
                    double beta_3, double beta_4, double p)
{
  double x2 = x*x; // r_s
  double Q_0 = -2.*A*(1. + alpha_1*x2);
  double Q_1 =  2.*A*(beta_1 * x + beta_2 * x2 + beta_3*x*x2 + beta_4 * pow(x2,p+1.));
  double Q_1prime = A * 
           ( beta_1 * 1./x + 2.*beta_2 + 3.*beta_3*x + 2.*(p+1.)*beta_4*pow(x2,p));
  double res = -2.*A*alpha_1*log(1. + 1./Q_1) - Q_0*Q_1prime/(Q_1*Q_1 + Q_1);

  return res;
}

void
PW92LCFunctional::point_lc(const PointInputData &id, PointOutputData &od, 
                           double &ec_local, double &decrs, double &deczeta)
{
  od.zero();
  const double fpp0 = 4./9. * 1./(pow(2., (1./3.)) - 1.);
  const double sixth = 1./6.;
  const double four_thirds = 4./3.;
  const double one_third = 1./3.;
  const double two_thirds = 2./3.;

  double rho = id.a.rho + id.b.rho;
  double zeta = (id.a.rho - id.b.rho)/rho;
  double x = pow(3./(4.*M_PI*rho), sixth);
  double rs = x*x;
  
  double epc    = F(x, 0.0310907,  0.21370, 7.5957,  3.5876, 1.6382,  0.49294, 1.00);
  double efc    = F(x, 0.01554535, 0.20548, 14.1189, 6.1977, 3.3662,  0.62517, 1.00);
  double alphac = F(x, 0.0168869, 0.11125, 10.357,  3.6231, 0.88026, 0.49671, 1.00);
     
  double f = 9./8.*fpp0*(pow(1.+zeta, four_thirds)+pow(1.-zeta, four_thirds)-2.);
  double zeta2 = zeta*zeta;
  double zeta4 = zeta2*zeta2;
  double delta_ec = -alphac * f / fpp0 * (1. - zeta4) + (efc - epc) * f * zeta4;
  double ec = epc + delta_ec;

  od.energy = ec * rho;
  ec_local = ec;

  if (compute_potential_) {
    if (!spin_polarized_) {
      double depc_dr_s0 = 
             dFdr_s(x, 0.0310907, 0.21370, 7.5957, 3.5876, 1.6382, 0.49294, 1.00);
      double dec_dr_s = depc_dr_s0;
      od.df_drho_a = od.df_drho_b = ec - (rs/3.)*dec_dr_s;
      decrs = dec_dr_s;
      deczeta = 0.;
    }
    else {
      double zeta3 = zeta2*zeta;
      double depc_dr_s0 = dFdr_s(x, 0.0310907, 0.21370, 7.5957,  
                                    3.5876, 1.6382,  0.49294, 1.00);
      double defc_dr_s1 = dFdr_s(x, 0.01554535, 0.20548, 14.1189, 
                                    6.1977, 3.3662,  0.62517, 1.00);
      double dalphac_dr_s = dFdr_s(x, 0.0168869, 0.11125, 10.357,  
                                    3.6231, 0.88026, 0.49671, 1.00);
      double dec_dr_s = depc_dr_s0*(1 - f*zeta4) + defc_dr_s1 * f * zeta4
                        + -dalphac_dr_s * f / fpp0 * (1 - zeta4);
      double fp = two_thirds * (pow((1+zeta),one_third) 
            - pow((1-zeta),one_third))/(pow(2.,one_third)-1);
      double dec_dzeta = 4.* zeta3 * f * (efc - epc - (-alphac/fpp0))
              + fp * (zeta4 * (efc - epc) + (1-zeta4)*(-alphac/fpp0));
      od.df_drho_a = ec - (rs/3.)*dec_dr_s - (zeta-1)*dec_dzeta;
      od.df_drho_b = ec - (rs/3.)*dec_dr_s - (zeta+1)*dec_dzeta;
      decrs = dec_dr_s;
      deczeta = dec_dzeta;
      } 
  }
}

/////////////////////////////////////////////////////////////////////////////
// PZ81LCFunctional
// Coded by Matt Leininger
// Used in P86 correlation functional
// J. P. Perdew and A. Zunger, Phys. Rev. B, 23, 5048, 1981.
// C. W. Murray, N. C. Handy, G. J. Laming, Mol. Phys., 78, 997, 1993.
//
static ClassDesc PZ81LCFunctional_cd(
  typeid(PZ81LCFunctional),"PZ81LCFunctional",1,"public LSDACFunctional",
  0, create<PZ81LCFunctional>, create<PZ81LCFunctional>);

PZ81LCFunctional::PZ81LCFunctional(StateIn& s):
  SavableState(s),
  LSDACFunctional(s)
{
}

PZ81LCFunctional::PZ81LCFunctional()
{
}

PZ81LCFunctional::PZ81LCFunctional(const Ref<KeyVal>& keyval):
  LSDACFunctional(keyval)
{
}

PZ81LCFunctional::~PZ81LCFunctional()
{
}

void
PZ81LCFunctional::save_data_state(StateOut& s)
{
  LSDACFunctional::save_data_state(s);
}

double
PZ81LCFunctional::Fec_rsgt1(double rs, double beta1, double beta2, double gamma)
{
  double sqrt_rs = sqrt(rs);
  double res = gamma / (1. + beta1*sqrt_rs + beta2*rs);
      
  return res;
}

double
PZ81LCFunctional::dFec_rsgt1_drho(double rs, double beta1, double beta2, double gamma,
                                  double &dec_drs)
{
  double ec = Fec_rsgt1(rs, beta1, beta2, gamma);
  double sqrt_rs = sqrt(rs);
  // double numer = 1.+ 7./6.*beta1*sqrt_rs + 4./3.*beta2*rs;
  double denom = 1. + beta1*sqrt_rs + beta2*rs;
  dec_drs = -ec/denom * (beta1/(2.*sqrt_rs) + beta2);
  // double res = ec * numer / denom;
  double res = (ec - rs/3.*dec_drs);
  
  return res;
}

double
PZ81LCFunctional::Fec_rslt1(double rs, double A, double B, double C, double D)
{
  double lnrs = log(rs);
  double res = A*lnrs + B + C*rs*lnrs + D*rs;
      
  return res;
}

double
PZ81LCFunctional::dFec_rslt1_drho(double rs, double A, double B, double C, double D,
                                  double &dec_drs)
{
  double lnrs = log(rs);
  double res = A*lnrs + B - A/3. + 2./3.*C*rs*lnrs + 1./3.*(2.*D - C)*rs;
  dec_drs = A/rs + C*lnrs + C + D;
  return res;
}


void
PZ81LCFunctional::point_lc(const PointInputData &id, PointOutputData &od, 
                          double &ec_local, double &decrs, double &deczeta)
{
  od.zero();

  const double Au = 0.0311;
  const double Ap = 0.01555;
  const double Bu = -0.048;
  const double Bp = -0.0269;
  const double Cu = 0.0020;
  const double Cp = 0.0007;
  const double Du = -0.0116;
  const double Dp = -0.0048;
  const double beta1u = 1.0529;
  const double beta1p = 1.3981;
  const double beta2u = 0.3334;
  const double beta2p = 0.2611;
  const double gammau = -0.1423;
  const double gammap = -0.0843;
  double rho = id.a.rho + id.b.rho;
  double zeta = (id.a.rho - id.b.rho)/rho;
  double rs = pow(3./(4.*M_PI*rho), (1./3.) );
  double fzeta = ( pow((1.+zeta), (4./3.)) + pow((1.-zeta), (4./3.)) - 2.)
                 / ( pow(2., (4./3.)) - 2. );

  double euc, epc;
  if (rs >= 1.) {
    // Ceperley U
    // euc = Fec_rsgt1(rs, 1.1581, 0.3446, -0.1471);
    // Ceperley P
    // epc = Fec_rsgt1(rs, 1.2520, 0.2567, -0.0790);
    // Ceperley-Adler U
    euc = Fec_rsgt1(rs, beta1u, beta2u, gammau);
    // Ceperley-Adler P
    epc = Fec_rsgt1(rs, beta1p, beta2p, gammap);
    }
  else { // rs < 1.
      // Ceperley U with A_u and B_u
      // euc = Fec_rslt1(rs, 0.0311, -0.048, 0.0014, -0.0108);
      // Ceperley P with A_p and B_p
      // epc = Fec_rslt1(rs, 0.01555, -0.0269, 0.0001, -0.0046);
      // Ceperley-Adler U with A_u and B_u
      euc = Fec_rslt1(rs, Au, Bu, Cu, Du);
      // Ceperley-Adler P with A_p and B_p
      epc = Fec_rslt1(rs, Ap, Bp, Cp, Dp);
    }
  double ec = euc + fzeta*(epc-euc);
  ec_local = ec;
  od.energy = ec * rho;

  if (compute_potential_) {
      double deuc_drs = 0.;
      double depc_drs = 0.;
      double deuc_drho, depc_drho;
      if (rs > 1.) {
          // Ceperley U
          // deuc_drho = dFec_rsgt1_drho(rs, 1.1581, 0.3446, -0.1471, deuc_drs);
          // Ceperley P
          // depc_drho = dFec_rsgt1_drho(rs, 1.2520, 0.2567, -0.0790, depc_drs);
          // Ceperley-Adler U
          deuc_drho = dFec_rsgt1_drho(rs, beta1u, beta2u, gammau, deuc_drs);
          // Ceperley-Adler P
          depc_drho = dFec_rsgt1_drho(rs, beta1p, beta2p, gammap, depc_drs);
         }
      else { // rs < 1.
         // Ceperley U with A_u and B_u
         // deuc_drho = dFec_rslt1_drho(rs, 0.0311, -0.048, 0.0014, -0.0108, deuc_drs);
         // Ceperley P with A_p and B_p
         // depc_drho = dFec_rslt1_drho(rs, 0.01555, -0.0269, 0.0001, -0.0046, depc_drs);
         // Ceperley-Adler U with A_u and B_u
         deuc_drho = dFec_rslt1_drho(rs, Au, Bu, Cu, Du, deuc_drs);
         // Ceperley-Adler P with A_p and B_p
         depc_drho = dFec_rslt1_drho(rs, Ap, Bp, Cp, Dp, depc_drs);
         }
      double dfzeta_dzeta = 4./3.*( pow((1.+zeta), (1./3.)) - pow((1.-zeta), (1./3.)) )
                      / (pow(2., (4./3.)) - 2.);
      decrs = deuc_drs + fzeta*(depc_drs - deuc_drs);
      deczeta = dfzeta_dzeta*(epc - euc);
      od.df_drho_a = deuc_drho + fzeta*(depc_drho - deuc_drho) +
                     (epc - euc)*(1.-zeta)*dfzeta_dzeta;
      od.df_drho_b = deuc_drho + fzeta*(depc_drho - deuc_drho) +
                     (epc - euc)*(-1.-zeta)*dfzeta_dzeta;
           
      } 
}

/////////////////////////////////////////////////////////////////////////////
// VWNLCFunctional
// Coded by Matt Leininger

static ClassDesc VWNLCFunctional_cd(
  typeid(VWNLCFunctional),"VWNLCFunctional",1,"public LSDACFunctional",
  0, create<VWNLCFunctional>, create<VWNLCFunctional>);

VWNLCFunctional::VWNLCFunctional(StateIn& s):
  SavableState(s), LSDACFunctional(s)
{
  s.get(Ap_);
  s.get(Af_);
  s.get(A_alpha_);

  s.get(x0p_mc_);
  s.get(bp_mc_);
  s.get(cp_mc_);
  s.get(x0f_mc_);
  s.get(bf_mc_);
  s.get(cf_mc_);
  s.get(x0_alpha_mc_);
  s.get(b_alpha_mc_);
  s.get(c_alpha_mc_);
  
  s.get(x0p_rpa_);
  s.get(bp_rpa_);
  s.get(cp_rpa_);
  s.get(x0f_rpa_);
  s.get(bf_rpa_);
  s.get(cf_rpa_);
  s.get(x0_alpha_rpa_);
  s.get(b_alpha_rpa_);
  s.get(c_alpha_rpa_);
}

VWNLCFunctional::VWNLCFunctional()
{
  init_constants();
}

VWNLCFunctional::~VWNLCFunctional()
{
}