<commit_before>#include "functional.hpp"

template<class num>
static num becke_alpha(const num &na, const num &gaa)
{
  const parameter c = pow(81/(4*M_PI),1.0/3.0)/2;
  const parameter d = 0.0042;
  num na43 = pow(na,4.0/3.0);
  num lda = -c*na43;
  num chi2 = gaa*pow(na,-8.0/3.0);
  num b88 = -(d*na43*chi2)/(1+6*d*sqrtx_asinh_sqrtx(chi2));
  return lda + b88;
}

template<class num>
static num becke_corr(const num &na, const num &gaa)
{
  const parameter d = 0.0042;
  num na43 = pow(na,4.0/3.0);
  num chi2 = gaa*pow(na,-8.0/3.0);
  return -(d*na43*chi2)/(1+6*d*sqrtx_asinh_sqrtx(chi2));
}

// Short range becke exchange as used in camb3lyp If mu=0 this reduces
// to the standard beckex for which must be used instead.
// FIXME: The erf + something is basically the erf minus
// its asymptotic expansion. This is horribel for numerics,
// will have to code a special function. 
// As coded here the code will fail if mu = 0, in which case the
// regular beckex should be used.

template<class num>
static num becke_sr(parameter mu, const num &na, const num &gaa)
{
  const parameter cparam = pow(81/(4*M_PI),1.0/3.0)/2;
  const parameter d = 0.0042;
  num na43 = pow(na,4.0/3.0);
  num chi2 = gaa*pow(na,-8.0/3.0);
  num K = 2*(cparam + (d*chi2)/(1+6*d*sqrtx_asinh_sqrtx(chi2)));
  num a = mu*sqrt(K)/(6*M_PI*pow(na,1.0/3.0));
  num b = expm1(-1/(4*a*a));
  num c = 2*a*a*b + 0.5;
  return -0.5*na43*K*(1-8.0/3.0*a*(sqrt(M_PI)*erf(1/(2*a))+2*a*(b-c)));
}

template<class num>
static num beckex(const densvars<num> &d)
{
  return becke_alpha(d.a,d.gaa) + becke_alpha(d.b,d.gbb);
}

template<class num>
static num beckexcorr(const densvars<num> &d)
{
  return becke_corr(d.a,d.gaa) + becke_corr(d.b,d.gbb);
}

template<class num>
static num beckesrx(const densvars<num> &d)
{
  parameter mu = d.get_param(XC_RANGESEP_MU);
  return becke_sr(mu,d.a,d.gaa) + becke_sr(mu,d.b,d.gbb);
}

FUNCTIONAL(XC_BECKEX) = {
  "Becke 88 exchange",
  "Becke 88 exchange including Slater part\n"
  "A.D. Becke, Density-functional exchange-energy approximation\n"
  "with correct asymptotic behaviour, Phys. Rev. A38 (1988) 3098-3100.\n"
  "Implemented by <NAME>\n"<commit_msg>Bugfix for short range becke contributed by <NAME>.<commit_after>#include "functional.hpp"

template<class num>
static num becke_alpha(const num &na, const num &gaa)
{
  const parameter c = pow(81/(4*M_PI),1.0/3.0)/2;
  const parameter d = 0.0042;
  num na43 = pow(na,4.0/3.0);
  num lda = -c*na43;
  num chi2 = gaa*pow(na,-8.0/3.0);
  num b88 = -(d*na43*chi2)/(1+6*d*sqrtx_asinh_sqrtx(chi2));
  return lda + b88;
}

template<class num>
static num becke_corr(const num &na, const num &gaa)
{
  const parameter d = 0.0042;
  num na43 = pow(na,4.0/3.0);
  num chi2 = gaa*pow(na,-8.0/3.0);
  return -(d*na43*chi2)/(1+6*d*sqrtx_asinh_sqrtx(chi2));
}

// Short range becke exchange as used in camb3lyp If mu=0 this reduces
// to the standard beckex for which must be used instead.
// FIXME: The erf + something is basically the erf minus
// its asymptotic expansion. This is horribel for numerics,
// will have to code a special function. 
// As coded here the code will fail if mu = 0, in which case the
// regular beckex should be used.

template<class num>
static num becke_sr(parameter mu, const num &na, const num &gaa)
{
  const parameter cparam = pow(81/(4*M_PI),1.0/3.0)/2;
  const parameter d = 0.0042;
  num na43 = pow(na,4.0/3.0);
  num chi2 = gaa*pow(na,-8.0/3.0);
  num K = 2*(cparam + (d*chi2)/(1+6*d*sqrtx_asinh_sqrtx(chi2)));
  num a = mu*sqrt(K)/(6*sqrt(M_PI)*pow(na,1.0/3.0));
  num b = expm1(-1/(4*a*a));
  num c = 2*a*a*b + 0.5;
  return -0.5*na43*K*(1-8.0/3.0*a*(sqrt(M_PI)*erf(1/(2*a))+2*a*(b-c)));
}

template<class num>
static num beckex(const densvars<num> &d)
{
  return becke_alpha(d.a,d.gaa) + becke_alpha(d.b,d.gbb);
}

template<class num>
static num beckexcorr(const densvars<num> &d)
{
  return becke_corr(d.a,d.gaa) + becke_corr(d.b,d.gbb);
}

template<class num>
static num beckesrx(const densvars<num> &d)
{
  parameter mu = d.get_param(XC_RANGESEP_MU);
  return becke_sr(mu,d.a,d.gaa) + becke_sr(mu,d.b,d.gbb);
}

FUNCTIONAL(XC_BECKEX) = {
  "Becke 88 exchange",
  "Becke 88 exchange including Slater part\n"
  "A.D. Becke, Density-functional exchange-energy approximation\n"
  "with correct asymptotic behaviour, Phys. Rev. A38 (1988) 3098-3100.\n"
  "Implemented by <NAME>rom\n"
{'fraction_non_alphanumeric': 0.16899153898925223, 'fraction_numerical': 0.05991310313286074, 'mean_word_length': 3.9648127128263337, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3920201', 'n_tokens_mistral': 2021, 'n_tokens_neox': 1915, 'n_words': 454}