<commit_before>/*
 * XCFun, an arbitrary order exchange-correlation library
 * Copyright (C) 2020 <NAME> and contributors.
 *
 * This file is part of XCFun.
 *
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
 *
 * For information on the complete list of contributors to the
 * XCFun library, see: <https://xcfun.readthedocs.io/>
 */

#include "config.hpp"
#include "constants.hpp"
#include "functional.hpp"
#include "pw92eps.hpp"

/*! Common code for SCAN and SCAN-like functionals
 *
 *  Implemented by <NAME>,
 *
 *  SCAN    - <NAME>, <NAME>, and <NAME>, Phys. Rev. Lett. 115, 036402
 * (2015) rSCAN   - <NAME> and <NAME>, J. Chem. Phys. 150, 161101 (2019)
 *  r2SCAN  - <NAME>, <NAME>, <NAME>, <NAME>, & J Sun; J. Chem. Phys. Lett.;
 * Accepted (DOI: 10.1021/acs.jpclett.0c02405)
 *
 *  The r++SCAN (aka "rppSCAN") and r4SCAN functionals are relatives of r2SCAN. They
 * are not recommended for general use.
 *
 *  r++SCAN, r4SCAN - <NAME>, <NAME>, <NAME>, <NAME>, & J Sun; in preparation.
 */

namespace SCAN_eps {
template <class num>
static num get_SCAN_Fx(const num,
                       const num,
                       const num,
                       const int,
                       const int,
                       const int);
template <class num>
static num SCAN_X_Fx(const num, const num, const parameter, const int, const int);
template <class num>
static num r2SCAN_C(const densvars<num> &, const int, const int, const int);
template <class num>
static num scan_ec0(const num,
                    const num,
                    const num,
                    const parameter,
                    const parameter,
                    const parameter);
template <class num>
static num lda_0(const num, const parameter, const parameter, const parameter);
template <class num>
static num scan_ec1(const num,
                    const num,
                    const num,
                    const parameter[8],
                    const parameter,
                    const parameter,
                    const parameter,
                    const parameter,
                    const int);
template <class num>
static void get_lsda1(const num, const num, const num, num &, num &);
template <class num>
static void gcor2(const parameter[6], const num, const num, num &, num &);

template <class num> static num fx_unif(const num & d) {
  return (-0.75 * pow(3 / PI, 1.0 / 3.0)) * pow(d, 4.0 / 3.0);
}

template <class num>
static num get_SCAN_Fx(const num d_n,
                       const num d_g,
                       const num d_tau,
                       const int IALPHA,
                       const int IINTERP,
                       const int IDELFX) {
  const parameter ETA = 1.0e-3;
  const parameter TAU_R = 1.0e-4;
  const parameter A_REG = 1.0e-3;

  num tauw = d_g / (8.0 * d_n);

  num tauUnif = 0.0;
  if (IALPHA == 1) {
    tauUnif = ((3.0 / 10.0) * pow(3 * PI2, 2.0 / 3.0) * pow(d_n, 5.0 / 3.0)) + TAU_R;
  } else {
    tauUnif = (0.3 * pow(3 * PI2, 2.0 / 3.0) * pow(d_n, 5.0 / 3.0));
  }

  num alpha = 0.0;
  if (IALPHA == 0) {
    // alpha (SCAN)
    if (abs(d_tau - tauw) > 1.0e-14) {
      alpha = (d_tau - tauw) / tauUnif;
    }

  } else if (IALPHA == 1) {
    // alpha' (rSCAN)
    alpha = (d_tau - tauw) / tauUnif;
    alpha = pow(alpha, 3) / (pow2(alpha) + A_REG);

  } else if (IALPHA == 2) {
    // \bar{alpha} (r2SCAN, r4SCAN)
    if (abs(d_tau - tauw) > 1.0e-14) {
      alpha = (d_tau - tauw) / (tauUnif + ETA * tauw);
    }

  } else {
    // Unknown IALPHA choice
    printf("ERROR: Unknown IALPHA %d\n", IALPHA);
  }

  num s = 0.0;
  if (abs(d_g) > 1.0e-16) {
    s = sqrt(d_g) /
        (2.0 * pow(3.0 * PI2, 1.0 / 3.0) * pow(d_n, 4.0 / 3.0)); // Eq. (4)
  };

  num Fx = SCAN_X_Fx(s, alpha, ETA, IINTERP, IDELFX);

  return Fx;
}

template <class num>
static num SCAN_X_Fx(const num s,
                     const num alpha,
                     const parameter ETA,
                     const int IINTERP,
                     const int IDELFX) {
  const parameter A1 = 4.9479;
  const parameter K1 = 0.065;
  const parameter K0 = 0.174;
  const parameter MU = 10.0 / 81.0;

  const parameter IE_PARAMS[8] = {1.0,
                                  -0.667,
                                  -0.4445555,
                                  -0.663086601049,
                                  1.451297044490,
                                  -0.887998041597,
                                  0.234528941479,
                                  -0.023185843322};
  const parameter CFX1 = 0.667;
  const parameter CFX2 = 0.8;
  const parameter CFDX1 = 1.24;

  const parameter D_DAMP2 = 0.361;
  const parameter DX_DAMP4_P = 0.232;
  const parameter DX_DAMP4_A = 0.232;
  const parameter B1 = 0.156632;
  const parameter B2 = 0.12083;
  const parameter B3 = 0.5;
  const parameter B4 = MU * MU / K1 - 0.112654;

  const parameter ALPHA_GE = 20.0 / 27.0 + ETA * 5.0 / 3.0;

  // Interpolation function
  num oma = 1.0 - alpha;
  num ief = 0.0;
  if (IINTERP == 0) {
    // SCAN
    if (alpha < 1.0) {
      ief = exp(-CFX1 * alpha / oma);
    } else {
      ief = -CFDX1 * exp(CFX2 / oma);
    }
  } else if (IINTERP == 1) {
    // rSCAN
    if (alpha < 1.0e-13) {
      ief = exp(-CFX1 * alpha / oma);
    } else if (alpha < 2.5) {
      for (int i = 0; i < 8; i++) {
        ief += IE_PARAMS[i] * pow(alpha, i);
      }
    } else {
      ief = -CFDX1 * exp(CFX2 / oma);
    }
  } else {
    printf("ERROR: Unknown IINTERP %d\n", IINTERP);
  }

  // Single orbital enhancement
  num h0x = 1.0 + K0;

  // Slowly varying enhancement
  num p = s * s;
  num del_f2 = 0.0;
  num C2 = 0.0;
  num h1x = 0.0;

  if (IDELFX == 0) {
    // 2nd and 4th order gradient expansion corrections for SCAN interpolation
    num wfac = B4 * pow2(p) * exp(-B4 * p / MU);
    num vfac = B1 * p + B2 * oma * exp(-B3 * pow2(oma));
    num yfac = MU * p + wfac + pow2(vfac);

    h1x = 1.0 + K1 - K1 / (1.0 + yfac / K1);
  } else if (IDELFX == 1 || IDELFX == 2) {
    // 2nd order GE corrections for rSCAN interpolation
    for (int i = 1; i < 8; i++) {
      del_f2 += i * IE_PARAMS[i];
    }
    C2 = -del_f2 * (1.0 - h0x);

    num damp = exp(-pow2(p) / pow(D_DAMP2, 4));
    h1x = 1.0 + K1 - K1 / (1.0 + p * (MU + ALPHA_GE * C2 * damp) / K1);
  } else {
    printf("ERROR: Unknown IDELFX %d\n", IDELFX);
  }

  // Scaling correction
  num gx = 1.0 - exp(-A1 / pow(p, 1.0 / 4.0));

  // 4th order gradient enhancement
  num del_fx = 0.0;
  if (IDELFX == 2) {<commit_msg>fixed grad reg for SCAN functionals, calculate p directly (#146)<commit_after>/*
 * XCFun, an arbitrary order exchange-correlation library
 * Copyright (C) 2020 <NAME> and contributors.
 *
 * This file is part of XCFun.
 *
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
 *
 * For information on the complete list of contributors to the
 * XCFun library, see: <https://xcfun.readthedocs.io/>
 */

#include "config.hpp"
#include "constants.hpp"
#include "functional.hpp"
#include "pw92eps.hpp"

/*! Common code for SCAN and SCAN-like functionals
 *
 *  Implemented by <NAME>,
 *
 *  SCAN    - <NAME>, <NAME>, and <NAME>, Phys. Rev. Lett. 115, 036402
 * (2015) rSCAN   - <NAME> and <NAME>, J. Chem. Phys. 150, 161101 (2019)
 *  r2SCAN  - <NAME>, <NAME>, <NAME>, <NAME>, & J Sun; J. Chem. Phys. Lett.;
 * Accepted (DOI: 10.1021/acs.jpclett.0c02405)
 *
 *  The r++SCAN (aka "rppSCAN") and r4SCAN functionals are relatives of r2SCAN. They
 * are not recommended for general use.
 *
 *  r++SCAN, r4SCAN - <NAME>, <NAME>, <NAME>, <NAME>, & J Sun; in preparation.
 */

namespace SCAN_eps {
template <class num>
static num get_SCAN_Fx(const num,
                       const num,
                       const num,
                       const int,
                       const int,
                       const int);
template <class num>
static num SCAN_X_Fx(const num, const num, const parameter, const int, const int);
template <class num>
static num r2SCAN_C(const densvars<num> &, const int, const int, const int);
template <class num>
static num scan_ec0(const num,
                    const num,
                    const num,
                    const parameter,
                    const parameter,
                    const parameter);
template <class num>
static num lda_0(const num, const parameter, const parameter, const parameter);
template <class num>
static num scan_ec1(const num,
                    const num,
                    const num,
                    const parameter[8],
                    const parameter,
                    const parameter,
                    const parameter,
                    const parameter,
                    const int);
template <class num>
static void get_lsda1(const num, const num, const num, num &, num &);
template <class num>
static void gcor2(const parameter[6], const num, const num, num &, num &);

template <class num> static num fx_unif(const num & d) {
  return (-0.75 * pow(3 / PI, 1.0 / 3.0)) * pow(d, 4.0 / 3.0);
}

template <class num>
static num get_SCAN_Fx(const num d_n,
                       const num d_g,
                       const num d_tau,
                       const int IALPHA,
                       const int IINTERP,
                       const int IDELFX) {
  const parameter ETA = 1.0e-3;
  const parameter TAU_R = 1.0e-4;
  const parameter A_REG = 1.0e-3;

  num tauw = d_g / (8.0 * d_n);

  num tauUnif = 0.0;
  if (IALPHA == 1) {
    tauUnif = ((3.0 / 10.0) * pow(3 * PI2, 2.0 / 3.0) * pow(d_n, 5.0 / 3.0)) + TAU_R;
  } else {
    tauUnif = (0.3 * pow(3 * PI2, 2.0 / 3.0) * pow(d_n, 5.0 / 3.0));
  }

  num alpha = 0.0;
  if (IALPHA == 0) {
    // alpha (SCAN)
    if (abs(d_tau - tauw) > 1.0e-14) {
      alpha = (d_tau - tauw) / tauUnif;
    }

  } else if (IALPHA == 1) {
    // alpha' (rSCAN)
    alpha = (d_tau - tauw) / tauUnif;
    alpha = pow(alpha, 3) / (pow2(alpha) + A_REG);

  } else if (IALPHA == 2) {
    // \bar{alpha} (r2SCAN, r4SCAN)
    if (abs(d_tau - tauw) > 1.0e-14) {
      alpha = (d_tau - tauw) / (tauUnif + ETA * tauw);
    }

  } else {
    // Unknown IALPHA choice
    printf("ERROR: Unknown IALPHA %d\n", IALPHA);
  }

  num p = 0.0;
  if (abs(d_g) > 1.0e-16) {
    p = d_g /
        (4.0 * pow(3.0 * PI2, 2.0 / 3.0) * pow(d_n, 8.0 / 3.0));
  } else {
    p = 1e-16/
        (4.0 * pow(3.0 * PI2, 2.0 / 3.0) * pow(d_n, 8.0 / 3.0));
  }

  num Fx = SCAN_X_Fx(p, alpha, ETA, IINTERP, IDELFX);

  return Fx;
}

template <class num>
static num SCAN_X_Fx(const num p,
                     const num alpha,
                     const parameter ETA,
                     const int IINTERP,
                     const int IDELFX) {
  const parameter A1 = 4.9479;
  const parameter K1 = 0.065;
  const parameter K0 = 0.174;
  const parameter MU = 10.0 / 81.0;

  const parameter IE_PARAMS[8] = {1.0,
                                  -0.667,
                                  -0.4445555,
                                  -0.663086601049,
                                  1.451297044490,
                                  -0.887998041597,
                                  0.234528941479,
                                  -0.023185843322};
  const parameter CFX1 = 0.667;
  const parameter CFX2 = 0.8;
  const parameter CFDX1 = 1.24;

  const parameter D_DAMP2 = 0.361;
  const parameter DX_DAMP4_P = 0.232;
  const parameter DX_DAMP4_A = 0.232;
  const parameter B1 = 0.156632;
  const parameter B2 = 0.12083;
  const parameter B3 = 0.5;
  const parameter B4 = MU * MU / K1 - 0.112654;

  const parameter ALPHA_GE = 20.0 / 27.0 + ETA * 5.0 / 3.0;

  // Interpolation function
  num oma = 1.0 - alpha;
  num ief = 0.0;
  if (IINTERP == 0) {
    // SCAN
    if (alpha < 1.0) {
      ief = exp(-CFX1 * alpha / oma);
    } else {
      ief = -CFDX1 * exp(CFX2 / oma);
    }
  } else if (IINTERP == 1) {
    // rSCAN
    if (alpha < 1.0e-13) {
      ief = exp(-CFX1 * alpha / oma);
    } else if (alpha < 2.5) {
      for (int i = 0; i < 8; i++) {
        ief += IE_PARAMS[i] * pow(alpha, i);
      }
    } else {
      ief = -CFDX1 * exp(CFX2 / oma);
    }
  } else {
    printf("ERROR: Unknown IINTERP %d\n", IINTERP);
  }

  // Single orbital enhancement
  num h0x = 1.0 + K0;

  // Slowly varying enhancement
  num del_f2 = 0.0;
  num C2 = 0.0;
  num h1x = 0.0;

  if (IDELFX == 0) {
    // 2nd and 4th order gradient expansion corrections for SCAN interpolation
    num wfac = B4 * pow2(p) * exp(-B4 * p / MU);
    num vfac = B1 * p + B2 * oma * exp(-B3 * pow2(oma));
    num yfac = MU * p + wfac + pow2(vfac);

    h1x = 1.0 + K1 - K1 / (1.0 + yfac / K1);
  } else if (IDELFX == 1 || IDELFX == 2) {
    // 2nd order GE corrections for rSCAN interpolation
    for (int i = 1; i < 8; i++) {
      del_f2 += i * IE_PARAMS[i];
    }
    C2 = -del_f2 * (1.0 - h0x);

    num damp = exp(-pow2(p) / pow(D_DAMP2, 4));
    h1x = 1.0 + K1 - K1 / (1.0 + p * (MU + ALPHA_GE * C2 * damp) / K1);
  } else {
    printf("ERROR: Unknown IDELFX %d\n", IDELFX);
  }

  // Scaling correction
  num gx = 1.0 - exp(-A1 / pow(p, 1.0 / 4.0));

  // 4th order gradient enhancement
  num del_fx = 0.0;
  if (IDELFX == 2) {
{'fraction_non_alphanumeric': 0.13441621339691529, 'fraction_numerical': 0.061545339393487816, 'mean_word_length': 1.7375076483785437, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 65, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2116656', 'n_tokens_mistral': 5474, 'n_tokens_neox': 4900, 'n_words': 1486}