<commit_before>
/* Needs some work to work with ctaylor */
#include "functional.hpp"
#include "slater.hpp"
#include "taylor.hpp"

// This is the function we want to find roots for
template<class T>
T BR_y(const T &x)
{
  return x*exp(-2.0/3.0*x)/(x-2);
}


// Determine x as a function of y = BR_y(x), using first 
// a guess and then root polishing.
// We need about 3-5 NR iterations depending on the y value.
ireal_t BR(ireal_t y)
{
  ireal_t d;
  taylor<ireal_t,1,2> x(0,0),fx;
  // More or less clever starting guesses
  if (y < 0)
    {
      if (y > -0.65)
	x[0] = -2*y;
      else
	x[0] = (6*exp(4.0/3.0)*y + 8)/(3*exp(4.0/3.0)*y+1);
    }
  else
    {
      if (y > 0.1)
	{
	  x[0] = (6*exp(4.0/3.0)*y + 8)/(3*exp(4.0/3.0)*y+1); 
	  if (y < 0.7)
	    x[0] -= 0.7*exp(-5*y); // Pragmatic correction
	}
      else
	{
	  x[0] = -3.0/2.0*log(y) + y*(47 + y*(-850 + y*4800));
	}
    }
  //Polish root with Newton-Raphson (or Halley's method)
  int niter = 0;
  do
    {
      fx = BR_y(x);
      // Newton:
      // d = (y-fx[0])/fx[1];
      // Halley:
      ireal_t yn = fx[0] - y;
      d = -yn*fx[1]/(fx[1]*fx[1] - yn*fx[2]);
      x[0] += d;
      if (++niter > 100)
	{
	  assert(0 && "Newton not converging in BR"); 
	  return 0;
	}
    }
  while (fabs(d)>1e-12);
  return x[0];
}


// Obtain the Taylor expansion of x(y), which is the
// inverse of BR_y. Use linear method for simplicity.
template<class T, int Ndeg>
void BR_taylor(const T &y0, taylor<T,1,Ndeg> &t)
{
  taylor<T,1,Ndeg> f,d;
  t = 0;
  t[0] = BR(y0);
  t[1] = 1;
  f = BR_y(t);
  t[1] = 1/f[1];
  // Linear method, for quadratic see i.e. Brent & Kung ~197x
  for (int i=2;i<=Ndeg;i++) 
    {
      f = BR_y(t);
      t[i] = -f[i]*t[1];
    }
}

/* This is a fully differentiable solver for Eq.(21) in 
   Becke and Roussel, PRA 39, 1989. t is the right hand
   side value, x is returned.
 */ 
template<class T,int Nvar>
static ctaylor<T,Nvar> BR(const ctaylor<T,Nvar> &t)
{
  taylor<T,1,Nvar> tmp;
  BR_taylor(t.c[0],tmp);
  printf("here x is %.15e\n",tmp[0]);

  ctaylor<T,Nvar> res = tmp[0];
  for (int i=1;i<=Nvar;i++)
    res += tmp[i]*pow(t-t.c[0],i);
  return res;
}

template<class num>
static num polarized(const num &na,
		     const num &gaa,
		     const num &lapa,
		     const num &taua,
		     const num &jpaa) // Becke tau here, no factor 1/2
{
  //  const parameter gam = 1.00; // like in the 96 paper
  //  const parameter gam = 0.80;
  printf("na %e, gaa %e, lapa %e, taua %e",na.c[0],gaa.c[0],lapa.c[0],taua.c[0]);
  num Q = (lapa - 2*taua + (0.5*gaa + 2*jpaa)/na)/6.0;
  printf("Q[0] is %.15e\n",Q.c[0]);
  num x = BR(2.0/3.0*pow(M_PI,2.0/3.0)*pow(na,5.0/3.0)/Q);
  printf("x[0] is %.15e\n",x.c[0]);
  num BRarg = 2.0/3.0*pow(M_PI,2.0/3.0)*pow(na,5.0/3.0)/Q;
  printf("br arg is %.15e\n",BRarg.c[0]);
  num b = cbrt(pow3(x)*exp(-x)/(8*M_PI*na));
  return -(1-(1+0.5*x)*exp(-x))/b; //FIXME: use expm1
}

template<class num>
static num brx(const densvars<num> &d)
{
  return 0.5*(d.a*polarized(d.a,d.gaa,d.lapa,2*d.taua,d.jpaa) +
	      d.a*polarized(d.b,d.gbb,d.lapb,2*d.taub,d.jpbb));
}

template<class num>
static num brc(const densvars<num> &d)
{
  parameter cab = 0.63, caa = 0.88;
  num UXa = polarized(d.a,d.gaa,d.lapa,2*d.taua,d.jpaa);
  num UXb = polarized(d.b,d.gbb,d.lapb,2*d.taub,d.jpbb);
  num zaa = abs(caa*(2.0/UXa));
  num zbb = abs(caa*(2.0/UXb));
  num zab = abs(cab*(1.0/UXa + 1.0/UXb));
  num ECopp = -0.8*d.a*d.b*zab*zab*(1-log(1+zab)/zab);
  num ECaa = -0.01*d.a*(d.taua - (0.25*d.gaa + d.jpaa)/d.a)*pow(zaa,4)*(1-2/zaa*log(1+zaa/2));
  num ECbb = -0.01*d.b*(d.taub - (0.25*d.gbb + d.jpbb)/d.b)*pow(zbb,4)*(1-2/zbb*log(1+zbb/2));
  return ECopp + ECaa + ECbb;
}


FUNCTIONAL(XC_BRX) = {
  "Becke-Roussells exchange with jp dependence",
  "Add info here"
  "Implemented by Ulf Ekstrom\n",
  XC_DENSITY | XC_GRADIENT | XC_KINETIC | XC_LAPLACIAN | XC_JP,
  ENERGY_FUNCTION(brx)
};

FUNCTIONAL(XC_BRC) = {
  "Becke-Roussells correlation with jp dependence",
  "Add info here"
  "Implemented by Ulf Ekstrom\n",
  XC_DENSITY | XC_GRADIENT | XC_KINETIC | XC_LAPLACIAN | XC_JP,
  ENERGY_FUNCTION(brc)
};
<commit_msg>BRX now verified against the implementation by Tellgren.<commit_after>
/* Needs some work to work with ctaylor */
#include "functional.hpp"
#include "slater.hpp"
#include "taylor.hpp"

// This is the function we want to find roots for
template<class T>
T BR_y(const T &x)
{
  return x*exp(-2.0/3.0*x)/(x-2);
}


// Determine x as a function of y = BR_y(x), using first 
// a guess and then root polishing.
// We need about 3-5 NR iterations depending on the y value.
ireal_t BR(ireal_t y)
{
  ireal_t d;
  taylor<ireal_t,1,2> x(0,0),fx;
  // More or less clever starting guesses
  if (y < 0)
    {
      if (y > -0.65)
	x[0] = -2*y;
      else
	x[0] = (6*exp(4.0/3.0)*y + 8)/(3*exp(4.0/3.0)*y+1);
    }
  else
    {
      if (y > 0.1)
	{
	  x[0] = (6*exp(4.0/3.0)*y + 8)/(3*exp(4.0/3.0)*y+1); 
	  if (y < 0.7)
	    x[0] -= 0.7*exp(-5*y); // Pragmatic correction
	}
      else
	{
	  x[0] = -3.0/2.0*log(y) + y*(47 + y*(-850 + y*4800));
	}
    }
  //Polish root with Newton-Raphson (or Halley's method)
  int niter = 0;
  do
    {
      fx = BR_y(x);
      // Newton:
      // d = (y-fx[0])/fx[1];
      // Halley:
      ireal_t yn = fx[0] - y;
      d = -yn*fx[1]/(fx[1]*fx[1] - yn*fx[2]);
      x[0] += d;
      if (++niter > 100)
	{
	  assert(0 && "Newton not converging in BR"); 
	  return 0;
	}
    }
  while (fabs(d)>1e-12);
  return x[0];
}


// Obtain the Taylor expansion of x(y), which is the
// inverse of BR_y. Use linear method for simplicity.
template<class T, int Ndeg>
void BR_taylor(const T &y0, taylor<T,1,Ndeg> &t)
{
  taylor<T,1,Ndeg> f,d;
  t = 0;
  t[0] = BR(y0);
  t[1] = 1;
  f = BR_y(t);
  t[1] = 1/f[1];
  // Linear method, for quadratic see i.e. Brent & Kung ~197x
  for (int i=2;i<=Ndeg;i++) 
    {
      f = BR_y(t);
      t[i] = -f[i]*t[1];
    }
}

/* This is a fully differentiable solver for Eq.(21) in 
   Becke and Roussel, PRA 39, 1989. t is the right hand
   side value, x is returned.
 */ 
template<class T,int Nvar>
static ctaylor<T,Nvar> BR(const ctaylor<T,Nvar> &t)
{
  taylor<T,1,Nvar> tmp;
  BR_taylor(t.c[0],tmp);

  ctaylor<T,Nvar> res = tmp[0];
  for (int i=1;i<=Nvar;i++)
    res += tmp[i]*pow(t-t.c[0],i);
  return res;
}

template<class num>
static num polarized(const num &na,
		     const num &gaa,
		     const num &lapa,
		     const num &taua,
		     const num &jpaa) // Becke tau here, no factor 1/2
{
  // The original BR article has a gamma constant, here this is put to 1.0
  num Q = (lapa - 2*taua + (0.5*gaa + 2*jpaa)/na)/6.0;
  num x = BR(2.0/3.0*pow(M_PI,2.0/3.0)*pow(na,5.0/3.0)/Q);
  num BRarg = 2.0/3.0*pow(M_PI,2.0/3.0)*pow(na,5.0/3.0)/Q;
  num b = cbrt(pow3(x)*exp(-x)/(8*M_PI*na));
  return -(1-(1+0.5*x)*exp(-x))/b; //FIXME: use expm1
}

template<class num>
static num brx(const densvars<num> &d)
{
  return 0.5*(d.a*polarized(d.a,d.gaa,d.lapa,2*d.taua,d.jpaa) +
	      d.b*polarized(d.b,d.gbb,d.lapb,2*d.taub,d.jpbb));
}

template<class num>
static num brc(const densvars<num> &d)
{
  parameter cab = 0.63, caa = 0.88;
  num UXa = polarized(d.a,d.gaa,d.lapa,2*d.taua,d.jpaa);
  num UXb = polarized(d.b,d.gbb,d.lapb,2*d.taub,d.jpbb);
  num zaa = abs(caa*(2.0/UXa));
  num zbb = abs(caa*(2.0/UXb));
  num zab = abs(cab*(1.0/UXa + 1.0/UXb));
  num ECopp = -0.8*d.a*d.b*zab*zab*(1-log(1+zab)/zab);
  num ECaa = -0.01*d.a*(d.taua - (0.25*d.gaa + d.jpaa)/d.a)*pow(zaa,4)*(1-2/zaa*log(1+zaa/2));
  num ECbb = -0.01*d.b*(d.taub - (0.25*d.gbb + d.jpbb)/d.b)*pow(zbb,4)*(1-2/zbb*log(1+zbb/2));
  return ECopp + ECaa + ECbb;
}


FUNCTIONAL(XC_BRX) = {
  "Becke-Roussells exchange with jp dependence",
  "See Becke, Canadian Journal of Chemistry, 1996, 74(6): 995-997"
  "Implemented by Ulf Ekstrom\n",
  XC_DENSITY | XC_GRADIENT | XC_KINETIC | XC_LAPLACIAN | XC_JP,
  ENERGY_FUNCTION(brx)
};

FUNCTIONAL(XC_BRC) = {
  "Becke-Roussells correlation with jp dependence",
  "See Becke, Canadian Journal of Chemistry, 1996, 74(6): 995-997"
  "Implemented by Ulf Ekstrom\n",
  XC_DENSITY | XC_GRADIENT | XC_KINETIC | XC_LAPLACIAN | XC_JP,
  ENERGY_FUNCTION(brc)
};

{'fraction_non_alphanumeric': 0.20584622942765904, 'fraction_numerical': 0.05502333578973225, 'mean_word_length': 2.8122659176029963, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 41, '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': '6683589', 'n_tokens_mistral': 4254, 'n_tokens_neox': 3962, 'n_words': 836}