functional.cc

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

VWNLCFunctional::VWNLCFunctional(const Ref<KeyVal>& keyval):
  LSDACFunctional(keyval)
{
  init_constants();
  Ap_  = keyval->doublevalue("Ap", KeyValValuedouble(Ap_));
  Af_  = keyval->doublevalue("Af", KeyValValuedouble(Af_));
  A_alpha_  = keyval->doublevalue("A_alpha", KeyValValuedouble(A_alpha_));
 
  x0p_mc_  = keyval->doublevalue("x0p_mc", KeyValValuedouble(x0p_mc_));
  bp_mc_   = keyval->doublevalue("bp_mc", KeyValValuedouble(bp_mc_));
  cp_mc_   = keyval->doublevalue("cp_mc", KeyValValuedouble(cp_mc_));
  x0f_mc_  = keyval->doublevalue("x0f_mc", KeyValValuedouble(x0f_mc_));
  bf_mc_   = keyval->doublevalue("bf_mc", KeyValValuedouble(bf_mc_));
  cf_mc_   = keyval->doublevalue("cf_mc", KeyValValuedouble(cf_mc_));
  x0_alpha_mc_ = keyval->doublevalue("x0_alpha_mc",
                                     KeyValValuedouble(x0_alpha_mc_));
  b_alpha_mc_  = keyval->doublevalue("b_alpha_mc",
                                     KeyValValuedouble(b_alpha_mc_));
  c_alpha_mc_  = keyval->doublevalue("c_alpha_mc",
                                     KeyValValuedouble(c_alpha_mc_));

  x0p_rpa_ = keyval->doublevalue("x0p_rpa", KeyValValuedouble(x0p_rpa_));
  bp_rpa_  = keyval->doublevalue("bp_rpa", KeyValValuedouble(bp_rpa_));
  cp_rpa_  = keyval->doublevalue("cp_rpa", KeyValValuedouble(cp_rpa_));
  x0f_rpa_ = keyval->doublevalue("x0f_rpa", KeyValValuedouble(x0f_rpa_));
  bf_rpa_  = keyval->doublevalue("bf_rpa", KeyValValuedouble(bf_rpa_));
  cf_rpa_  = keyval->doublevalue("cf_rpa", KeyValValuedouble(cf_rpa_));
  x0_alpha_rpa_ = keyval->doublevalue("x0_alpha_rpa",
                                      KeyValValuedouble(x0_alpha_rpa_));
  b_alpha_rpa_  = keyval->doublevalue("b_alpha_rpa",
                                      KeyValValuedouble(b_alpha_rpa_));
  c_alpha_rpa_  = keyval->doublevalue("c_alpha_rpa",
                                      KeyValValuedouble(c_alpha_rpa_));   
}

void
VWNLCFunctional::save_data_state(StateOut& s)
{
  LSDACFunctional::save_data_state(s);

  s.put(Ap_);
  s.put(Af_);
  s.put(A_alpha_);

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

void
VWNLCFunctional::init_constants()
{
  Ap_ = 0.0310907;
  Af_ = 0.01554535;
  A_alpha_ = -1./(6.*M_PI*M_PI);

  x0p_mc_ = -0.10498;
  bp_mc_  = 3.72744;
  cp_mc_  = 12.9352;
  x0f_mc_ = -0.32500;
  bf_mc_  = 7.06042;
  cf_mc_  = 18.0578;
  x0_alpha_mc_ = -0.00475840;
  b_alpha_mc_  = 1.13107;
  c_alpha_mc_  = 13.0045;
  
  x0p_rpa_ = -0.409286;
  bp_rpa_  = 13.0720;
  cp_rpa_  = 42.7198;
  x0f_rpa_ = -0.743294;
  bf_rpa_  = 20.1231;
  cf_rpa_  = 101.578;
  x0_alpha_rpa_ = -0.228344;
  b_alpha_rpa_  = 1.06835;
  c_alpha_rpa_  = 11.4813;
}

double
VWNLCFunctional::F(double x, double A, double x0, double b, double c)
{
  double x2 = x*x;
  double x02 = x0*x0;
  double Xx = x2 + b*x + c;
  double Xx0 = x02 + b*x0 + c;
  double Q = sqrt(4.*c-b*b);
  double res
      = A * ( log(x2/Xx)
              + 2.*b/Q * atan(Q/(2.*x+b))
              - b*x0/Xx0 * ( log((x-x0)*(x-x0)/Xx)
                           + 2.*(b+2.*x0)/Q * atan(Q/(2.*x+b)) ) );
  return res;
}

double
VWNLCFunctional::dFdr_s(double x, double A, double x0, double b, double c)
{
  double x2 = x*x;
  double x02 = x0*x0;
  double Xx = x2 + b*x +c;
  double Xx0 = x02 + b*x0 + c;
  double Q = sqrt(4.*c-b*b);
  double res
      = A * ( 1./x2 - 1./Xx - b/(2.*Xx*x) 
          + ((x0*(2.*x0+b))/Xx0 - 1) * (2.*b)/(x*(Q*Q+(2.*x+b)*(2.*x+b))) 
          - (b*x0)/(x*(x-x0)*Xx0) + (b*x0*(1+(b/(2.*x))))/(Xx0*Xx) );
  return res;
}

void
VWNLCFunctional::point_lc(const PointInputData &id, PointOutputData &od, 
                          double &ec_local, double &decrs, double &deczeta)
{
}

/////////////////////////////////////////////////////////////////////////////
// VWN1LCFunctional
// Coded by Matt Leininger

static ClassDesc VWN1LCFunctional_cd(
  typeid(VWN1LCFunctional),"VWN1LCFunctional",1,"public LSDACFunctional",
  0, create<VWN1LCFunctional>, create<VWN1LCFunctional>);

VWN1LCFunctional::VWN1LCFunctional(StateIn& s):
  SavableState(s),
  VWNLCFunctional(s)
{
  s.get(x0p_);
  s.get(bp_);
  s.get(cp_);
  s.get(x0f_);
  s.get(bf_);
  s.get(cf_);
}

VWN1LCFunctional::VWN1LCFunctional()
{
  x0p_ = x0p_mc_;
  bp_  = bp_mc_;
  cp_  = cp_mc_;
  x0f_ = x0f_mc_;
  bf_  = bf_mc_;
  cf_  = cf_mc_;
}

VWN1LCFunctional::VWN1LCFunctional(int use_rpa)
{
  if (use_rpa) {
      x0p_ = x0p_rpa_;
      bp_  = bp_rpa_;
      cp_  = cp_rpa_;
      x0f_ = x0f_rpa_;
      bf_  = bf_rpa_;
      cf_  = cf_rpa_;
    }
  else {
      x0p_ = x0p_mc_;
      bp_  = bp_mc_;
      cp_  = cp_mc_;
      x0f_ = x0f_mc_;
      bf_  = bf_mc_;
      cf_  = cf_mc_;
    }
  
}

VWN1LCFunctional::VWN1LCFunctional(const Ref<KeyVal>& keyval):
  VWNLCFunctional(keyval)
{
  int vwn1rpa = keyval->booleanvalue("rpa", KeyValValueboolean(0));
  if (vwn1rpa) {
      x0p_ = x0p_rpa_;
      bp_  = bp_rpa_;
      cp_  = cp_rpa_;
      x0f_ = x0f_rpa_;
      bf_  = bf_rpa_;
      cf_  = cf_rpa_;
    }
  else {
      x0p_ = x0p_mc_;
      bp_  = bp_mc_;
      cp_  = cp_mc_;
      x0f_ = x0f_mc_;
      bf_  = bf_mc_;
      cf_  = cf_mc_;
    }

  x0p_ = keyval->doublevalue("x0p", KeyValValuedouble(x0p_));
  bp_ = keyval->doublevalue("bp", KeyValValuedouble(bp_));
  cp_ = keyval->doublevalue("cp", KeyValValuedouble(cp_));
  x0f_ = keyval->doublevalue("x0f", KeyValValuedouble(x0f_));
  bf_ = keyval->doublevalue("bf", KeyValValuedouble(bf_));
  cf_ = keyval->doublevalue("cf", KeyValValuedouble(cf_));
}

VWN1LCFunctional::~VWN1LCFunctional()
{
}

void
VWN1LCFunctional::save_data_state(StateOut& s)
{
  VWNLCFunctional::save_data_state(s);
  s.put(x0p_);
  s.put(bp_);
  s.put(cp_);
  s.put(x0f_);
  s.put(bf_);
  s.put(cf_);
}

// Based on the VWN1 functional in Vosko, Wilk, and Nusair, Can. J. Phys.
// 58, 1200, (1980).
void
VWN1LCFunctional::point_lc(const PointInputData &id, PointOutputData &od, 
                          double &ec_local, double &decrs, double &deczeta)
{
  od.zero();

  double four_thirds = 4./3.;
  double one_third = 1./3.;
  double two_thirds = 2./3.;
  double sixth = 1./6.;
  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_, bp_, cp_);
  double efc    = F(x, Af_, x0f_, bf_, cf_);

  const double fpp0 = 4./9. * 1./(pow(2., (1./3.)) - 1.);
  double f = 9./8.*fpp0*(pow(1.+zeta, four_thirds)+pow(1.-zeta, four_thirds)-2.);
  double delta_ec = f * (efc - epc);
  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_, bp_, cp_);
      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 depc_dr_s0 = dFdr_s(x, Ap_, x0p_, bp_, cp_);
      double defc_dr_s1 = dFdr_s(x, Af_, x0f_, bf_, cf_);
      double dec_dr_s = depc_dr_s0 + f * (defc_dr_s1 - depc_dr_s0);
      double fp = two_thirds * (pow((1+zeta),one_third) 
            - pow((1-zeta),one_third))/(pow(2.,one_third)-1);
      // double dec_dzeta = depc_dr_s0 + fp * (efc - epc) + f * (defc_dr_s1 - depc_dr_s0);
      double dec_dzeta = fp * (efc - epc);       
      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;
      } 
  }
}

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

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

VWN2LCFunctional::VWN2LCFunctional()
{
}

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

VWN2LCFunctional::~VWN2LCFunctional()
{
}

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

void
VWN2LCFunctional::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_rpa / alphac_rpa) - 1.;
  double delta_erpa_rszeta = alphac_rpa * f / fpp0 * (1. + beta * zeta4);
  double delta_ec =  delta_erpa_rszeta + f*(delta_e_mc - delta_e_rpa);
  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);
      
      // 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_rpa - depc_dr_s0_mc;
      // double dec_dr_s_mc = depc_dr_s0*(1 - f*zeta4) + defc_dr_s1 * f * zeta4
      //                   + dalphac_dr_s * f / fpp0 * (1 - zeta4);
      // double dec_dzeta_mc = 4.* zeta3 * f * (efc_mc - epc_mc - (alphac_mc/fpp0))
      //        + fp * (zeta4 * (efc_mc - epc_mc) + (1-zeta4)*(alphac_mc/fpp0));
      
      double dec_dr_s = depc_dr_s0_mc + ddeltae_rpa_dr_s + f * (ddelta_e_mc - ddelta_e_rpa);
      
      double dec_dzeta = ddeltae_rpa_dzeta + fp * (delta_e_mc - delta_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;
    }
}

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

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

VWN3LCFunctional::VWN3LCFunctional(int mcp, int mce0)
{
  monte_carlo_prefactor_ = mcp;
  monte_carlo_e0_ = mce0;
}

VWN3LCFunctional::VWN3LCFunctional(const Ref<KeyVal>& keyval):
  VWNLCFunctional(keyval)
{
    monte_carlo_prefactor_ = keyval->booleanvalue("monte_carlo_prefactor",
                                                  KeyValValueboolean(1));
    monte_carlo_e0_ = keyval->booleanvalue("monte_carlo_e0",
                                           KeyValValueboolean(1));
}

VWN3LCFunctional::~VWN3LCFunctional()
{
}

void
VWN3LCFunctional::save_data_state(StateOut& s)
{
  VWNLCFunctional::save_data_state(s);
  s.put(monte_carlo_prefactor_);
  s.put(monte_carlo_e0_);