functional.cc

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

}

// based on the equations given on a NIST WWW site
void
VWN3LCFunctional::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;
  
  // Monte Carlo fitting parameters 
  double epc0_mc    = F(x, Ap_, x0p_mc_, bp_mc_, cp_mc_);
  double efc1_mc    = F(x, Af_, x0f_mc_, bf_mc_, cf_mc_);
  // double alphac_mc = F(x, A_alpha_, x0_alpha_mc_, b_alpha_mc_, c_alpha_mc_);
  // RPA fitting parameters
  double epc0_rpa    = F(x, Ap_, x0p_rpa_, bp_rpa_, cp_rpa_);
  double efc1_rpa    = F(x, Af_, x0f_rpa_, bf_rpa_, cf_rpa_);
  double alphac_rpa = F(x, A_alpha_, x0_alpha_rpa_, b_alpha_rpa_, c_alpha_rpa_);
    
  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_e_rpa = efc1_rpa - epc0_rpa;
  double delta_e_mc  = efc1_mc  - epc0_mc;
  double delta_erpa_rszeta = alphac_rpa * f / fpp0 * (1.-zeta4) + f * zeta4 * delta_e_rpa;
  double delta_ec;
  if (!monte_carlo_prefactor_) delta_ec = delta_erpa_rszeta;
  else delta_ec = delta_e_mc/delta_e_rpa * delta_erpa_rszeta;

  double ec;
  if (monte_carlo_e0_) ec = epc0_mc;
  else ec = epc0_rpa;
  ec += delta_ec;
  
  od.energy = ec * rho;
  ec_local = ec;

  if (compute_potential_) {
      // Monte Carlo fitting parameters
      double depc_dr_s0_mc = dFdr_s(x, Ap_, x0p_mc_, bp_mc_, cp_mc_);
      double defc_dr_s1_mc = dFdr_s(x, Af_, x0f_mc_, bf_mc_, cf_mc_);
      // double dalphac_dr_s_mc = dFdr_s(x, A_alpha_, x0_alpha_mc_, b_alpha_mc_, c_alpha_mc_);
      // RPA fitting parameters
      double depc_dr_s0_rpa = dFdr_s(x, Ap_, x0p_rpa_, bp_rpa_, cp_rpa_);
      double defc_dr_s1_rpa = dFdr_s(x, Af_, x0f_rpa_, bf_rpa_, cf_rpa_);
      double dalphac_dr_s_rpa = dFdr_s(x, A_alpha_, x0_alpha_rpa_, b_alpha_rpa_, c_alpha_rpa_);

      double fp = two_thirds * (pow((1+zeta),one_third)
                  - pow((1-zeta),one_third))/(pow(2.,one_third)-1);
           
      // RPA fitting parameters
      double ddelta_e_rpa = defc_dr_s1_rpa - depc_dr_s0_rpa;
      double ddeltae_rpa_dr_s = f / fpp0 * (1 - zeta4)* dalphac_dr_s_rpa 
                    + f * zeta4 * ddelta_e_rpa;
      double ddeltae_rpa_dzeta = alphac_rpa / fpp0 * 
         ( fp * (1.-zeta4) - 4.* f * zeta*zeta2) 
      + delta_e_rpa * ( fp*zeta4 + 4.*f*zeta*zeta2 );

      // Monte Carlo fitting parameters
      double ddelta_e_mc = defc_dr_s1_mc - depc_dr_s0_mc;
      
      double dec_dzeta, dec_dr_s;
      if (!monte_carlo_prefactor_) {
        dec_dzeta = ddeltae_rpa_dzeta;
        if (monte_carlo_e0_) dec_dr_s = depc_dr_s0_mc;
        else dec_dr_s = depc_dr_s0_rpa;
        dec_dr_s += ddeltae_rpa_dr_s;
        }
      else {
        dec_dzeta = delta_e_mc / delta_e_rpa * ddeltae_rpa_dzeta;
        if (monte_carlo_e0_) dec_dr_s = depc_dr_s0_mc;
        else dec_dr_s = depc_dr_s0_rpa;
        dec_dr_s += delta_erpa_rszeta/delta_e_rpa * ddelta_e_mc
                   + delta_e_mc/delta_e_rpa * ddeltae_rpa_dr_s 
                   - delta_erpa_rszeta*delta_e_mc/(delta_e_rpa*delta_e_rpa)
                   * ddelta_e_rpa;
        }        
        
      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;
    }
}

/////////////////////////////////////////////////////////////////////////////
// VWN4LCFunctional
// Coded by Matt Leininger
static ClassDesc VWN4LCFunctional_cd(
  typeid(VWN4LCFunctional),"VWN4LCFunctional",1,"public VWNLCFunctional",
  0, create<VWN4LCFunctional>, create<VWN4LCFunctional>);

VWN4LCFunctional::VWN4LCFunctional(StateIn& s):
  SavableState(s),
  VWNLCFunctional(s)
{
  s.get(monte_carlo_prefactor_);
}

VWN4LCFunctional::VWN4LCFunctional()
{
  monte_carlo_prefactor_ = 0;
}

VWN4LCFunctional::VWN4LCFunctional(const Ref<KeyVal>& keyval):
  VWNLCFunctional(keyval)
{
  monte_carlo_prefactor_ = keyval->booleanvalue("monte_carlo_prefactor",
                                                KeyValValueboolean(0));
}

VWN4LCFunctional::~VWN4LCFunctional()
{
}

void
VWN4LCFunctional::save_data_state(StateOut& s)
{
  VWNLCFunctional::save_data_state(s);
  s.put(monte_carlo_prefactor_);
}

// based on the equations given on a NIST WWW site
void
VWN4LCFunctional::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;
  
  // Monte Carlo fitting parameters 
  double epc_mc    = F(x, Ap_, x0p_mc_, bp_mc_, cp_mc_);
  double efc_mc    = F(x, Af_, x0f_mc_, bf_mc_, cf_mc_);
  // double alphac_mc = F(x, A_alpha_, x0_alpha_mc_, b_alpha_mc_, c_alpha_mc_);
  // RPA fitting parameters
  double epc_rpa    = F(x, Ap_, x0p_rpa_, bp_rpa_, cp_rpa_);
  double efc_rpa    = F(x, Af_, x0f_rpa_, bf_rpa_, cf_rpa_);
  double alphac_rpa = F(x, A_alpha_, x0_alpha_rpa_, b_alpha_rpa_, c_alpha_rpa_);

  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_e_rpa = efc_rpa - epc_rpa;
  double delta_e_mc  = efc_mc - epc_mc;
  // double beta = fpp0 * delta_e_mc / alphac_rpa - 1.;
  // double delta_erpa_rszeta = alphac_rpa * f / fpp0 * (1. + beta * zeta4);
  // use delta_e_rpa here instead of delta_e_mc to get VWN3
  double delta_erpa_rszeta = alphac_rpa * f / fpp0 * (1. - zeta4)
                           + f * delta_e_mc * zeta4;
  double delta_ec;
  if (!monte_carlo_prefactor_) delta_ec = delta_erpa_rszeta;
  else delta_ec = delta_e_mc/delta_e_rpa * delta_erpa_rszeta;
  // double ec = epc_rpa + delta_ec;
  double ec = epc_mc + delta_ec;
  
  od.energy = ec * rho;
  ec_local = ec;

  if (compute_potential_) {
      double zeta3 = zeta2*zeta;
      // Monte Carlo fitting parameters
      double depc_dr_s0_mc = dFdr_s(x, Ap_, x0p_mc_, bp_mc_, cp_mc_);
      double defc_dr_s1_mc = dFdr_s(x, Af_, x0f_mc_, bf_mc_, cf_mc_);
      // double dalphac_dr_s_mc = dFdr_s(x, A_alpha_, x0_alpha_mc_, b_alpha_mc_, c_alpha_mc_);
      // RPA fitting parameters
      double depc_dr_s0_rpa = dFdr_s(x, Ap_, x0p_rpa_, bp_rpa_, cp_rpa_);
      double defc_dr_s1_rpa = dFdr_s(x, Af_, x0f_rpa_, bf_rpa_, cf_rpa_);
      double dalphac_dr_s_rpa = dFdr_s(x, A_alpha_, x0_alpha_rpa_, b_alpha_rpa_,
                                       c_alpha_rpa_);
    
      double fp = two_thirds * (pow((1+zeta),one_third)
             - pow((1-zeta),one_third))/(pow(2.,one_third)-1);
      
      double ddelta_e_rpa = defc_dr_s1_rpa - depc_dr_s0_rpa;
      double ddelta_e_mc = defc_dr_s1_mc - depc_dr_s0_mc;
      // use ddelta_e_rpa here instead of ddelta_e_mc to get VWN3
      double ddeltae_rpa_dr_s = f / fpp0 * (1 - zeta4)* dalphac_dr_s_rpa 
                    + f * zeta4 * ddelta_e_mc;
      // use ddelta_e_rpa here instead of ddelta_e_mc to get VWN3
      double ddeltae_rpa_dzeta = alphac_rpa / fpp0 * 
        ( fp * (1.-zeta4) - 4.* f * zeta3) 
        + delta_e_mc * ( fp*zeta4 + 4.*f*zeta3);

      double dec_dzeta, dec_dr_s;
      if (!monte_carlo_prefactor_) {
          // dec_dr_s = depc_dr_s0_rpa + ddeltae_rpa_dr_s;
          dec_dr_s = depc_dr_s0_mc + ddeltae_rpa_dr_s;
          dec_dzeta = ddeltae_rpa_dzeta;
        }
      else {
          dec_dr_s = depc_dr_s0_rpa + delta_erpa_rszeta/delta_e_rpa * ddelta_e_mc
                   + delta_e_mc/delta_e_rpa * ddeltae_rpa_dr_s 
                   - delta_erpa_rszeta*delta_e_mc/(delta_e_rpa*delta_e_rpa)* ddelta_e_rpa;
          dec_dzeta = delta_e_mc / delta_e_rpa * ddeltae_rpa_dzeta;      
        }
      
      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;
    }
}

/////////////////////////////////////////////////////////////////////////////
// VWN5LCFunctional
// Coded by Matt Leininger

static ClassDesc VWN5LCFunctional_cd(
  typeid(VWN5LCFunctional),"VWN5LCFunctional",1,"public VWNLCFunctional",
  0, create<VWN5LCFunctional>, create<VWN5LCFunctional>);

VWN5LCFunctional::VWN5LCFunctional(StateIn& s):
  SavableState(s),
  VWNLCFunctional(s)
{
}

VWN5LCFunctional::VWN5LCFunctional()
{
}

VWN5LCFunctional::VWN5LCFunctional(const Ref<KeyVal>& keyval):
  VWNLCFunctional(keyval)
{
}

VWN5LCFunctional::~VWN5LCFunctional()
{
}

void
VWN5LCFunctional::save_data_state(StateOut& s)
{
  VWNLCFunctional::save_data_state(s);
}

// based on the equations given on a NIST WWW site
// based on the VWN5 functional in Vosko, Wilk, and Nusair, Can. J. Phys.
// 58, 1200, (1980). and Perdew and Wang, Phys. Rev. B, 45, 13244, (1992).
void
VWN5LCFunctional::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, Ap_, x0p_mc_, bp_mc_, cp_mc_);
  double efc    = F(x, Af_, x0f_mc_, bf_mc_, cf_mc_);
  double alphac = F(x, A_alpha_, x0_alpha_mc_, b_alpha_mc_, c_alpha_mc_);
     
  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 beta = fpp0 * (efc - epc) / alphac - 1.;
  double delta_ec = alphac * f / fpp0 * (1. + beta * zeta4);
  double ec = epc + delta_ec;

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

  if (compute_potential_) {
    if (!spin_polarized_) {
      double depc_dr_s0 = dFdr_s(x, Ap_, x0p_mc_, bp_mc_, cp_mc_); 
      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, Ap_, x0p_mc_, bp_mc_, cp_mc_);
      double defc_dr_s1 = dFdr_s(x, Af_, x0f_mc_, bf_mc_, cf_mc_);
      double dalphac_dr_s = dFdr_s(x, A_alpha_, x0_alpha_mc_, b_alpha_mc_, c_alpha_mc_);
      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;
      } 
  }
}

/////////////////////////////////////////////////////////////////////////////
// XalphaFunctional

static ClassDesc XalphaFunctional_cd(
  typeid(XalphaFunctional),"XalphaFunctional",1,"public DenFunctional",
  0, create<XalphaFunctional>, create<XalphaFunctional>);

XalphaFunctional::XalphaFunctional(StateIn& s):
  SavableState(s),
  DenFunctional(s)
{
  s.get(alpha_);
  factor_ = alpha_ * 2.25 * pow(3.0/(4.*M_PI), 1.0/3.0);
}

XalphaFunctional::XalphaFunctional()
{
  alpha_ = 0.70;
  factor_ = alpha_ * 2.25 * pow(3.0/(4.*M_PI), 1.0/3.0);
}

XalphaFunctional::XalphaFunctional(const Ref<KeyVal>& keyval):
  DenFunctional(keyval)
{
  alpha_ = keyval->doublevalue("alpha", KeyValValuedouble(0.70));
  factor_ = alpha_ * 2.25 * pow(3.0/(4.*M_PI), 1.0/3.0);
}

XalphaFunctional::~XalphaFunctional()
{
}

void
XalphaFunctional::save_data_state(StateOut& s)
{
  DenFunctional::save_data_state(s);
  s.put(alpha_);
}

void
XalphaFunctional::point(const PointInputData &id,
                        PointOutputData &od)
{
  od.zero();
  const double four_thirds = 4./3.;
  
  if (!spin_polarized_) {
      od.energy = - 2.0 * factor_ * id.a.rho * id.a.rho_13;
      if (compute_potential_) {
          od.df_drho_a = -four_thirds * factor_ * id.a.rho_13;
          od.df_drho_b = od.df_drho_a;
              }
    }
  else {
      double rhoa43 = id.a.rho * id.a.rho_13;
      double rhob43 = id.b.rho * id.b.rho_13;
      od.energy = - factor_ * (rhoa43 + rhob43);
      if (compute_potential_) {
          od.df_drho_a = -four_thirds * factor_ * id.a.rho_13;
          od.df_drho_b = -four_thirds * factor_ * id.b.rho_13;
              }
    }
}

void
XalphaFunctional::print(ostream& o) const
{
  o
    << indent << scprintf("XalphaFunctional: alpha = %12.8f", alpha_) << endl;
}

/////////////////////////////////////////////////////////////////////////////
// Becke88XFunctional

static ClassDesc Becke88XFunctional_cd(
  typeid(Becke88XFunctional),"Becke88XFunctional",1,"public DenFunctional",
  0, create<Becke88XFunctional>, create<Becke88XFunctional>);

Becke88XFunctional::Becke88XFunctional(StateIn& s):
  SavableState(s),
  DenFunctional(s)
{
  s.get(beta_);
  beta6_ = 6. * beta_;
}

Becke88XFunctional::Becke88XFunctional()
{
  beta_ = 0.0042;
  beta6_ = 6. * beta_;
}

Becke88XFunctional::Becke88XFunctional(const Ref<KeyVal>& keyval):
  DenFunctional(keyval)
{
  beta_ = keyval->doublevalue("beta", KeyValValuedouble(0.0042));
  beta6_ = 6. * beta_;
}

Becke88XFunctional::~Becke88XFunctional()
{
}

void
Becke88XFunctional::save_data_state(StateOut& s)
{
  DenFunctional::save_data_state(s);
  s.put(beta_);
}

in