#include "param.h"

/*********************************************/
/* Austausch-Korrelations-Wechselwirkung: 3D */
/*********************************************/

void evxc3d (double dr,      /* Schrittweite                      */
             int nr,         /* Anzahl Stuetzstellen              */
             double rf[],    /* Ortsraumgitter                    */
             double rho[],   /* Ladungsdichte                     */
             double pot[],   /* Potential                         */
             double *vxctot, /* Erwartungswert <Vxc>              */
             double *exctot, /* Erwartungswert <Exc>              */
             int modus)      /* Typ Austausch-Korrelationspot.    */
{
  int n;
  double vxc, fak, exc;
  *vxctot = 0.0;
  *exctot = 0.0;
  for (n = 1; n <= nr; n++) {
    if (modus == 1)
      vxc = vxc_hl (rho[n], &exc);
    else if (modus == 2)
      vxc = vxc_pz (rho[n], &exc);
    else {
      printf ("modus = %i undefined", modus);
      exit (1);
    }
    pot[n] += vxc;

    fak = 4.0 * PI * rf[n] * rf[n] * dr;
    /* Erwartungswert: Austausch-Korrelations-Potential */
    *vxctot += fak * vxc * rho[n];
    /* Erwartungswert: Austausch-Korrelations-Energie */
    *exctot += fak * exc;
  }
}

/**********************************************/
/* Austausch-Korrelations-Wechselwirkung      */
/* Hedin, Lundquist, J. Phys. C4, 2046 (1971) */
/**********************************************/

#define exfak1 0.0290898598909931995857496
/*             = 1/7 * (1 / (12 * pi **2)) **1/3  */
#define exfak2 33.8518310343458624550564
/*             = 21 * (4/3 * pi) **1/3  */
#define exfak3 0.7734
/*             = phaenomenolog. Parameter */
#define z1dz3  0.333333333333333333333333

double vxc_hl (double rho, double *exc)
{
  double arho = fabs (rho);

  if (arho > 1e-13) {
    double fak1 = exfak2 * pow (arho, z1dz3);
    double fak2 = 1.0 / fak1;
    *exc = - arho * exfak1
      * (0.75 * fak1 + exfak3 * ((1.0 + fak2*fak2*fak2) * log (1.0 +fak1)
                                 + 0.5 * fak2 - fak2 * fak2 - z1dz3));
    return (- exfak1 * (fak1 + exfak3 * log (1.0 + fak1)));
  }
  else {
    *exc = 0.0;
    return (0.0);
  }
}

#undef exfak1
#undef exfak2
#undef exfak3
#undef z1dz3

/**************************************************************/
/* Exchange-Correlation Potential:                            */
/* The correlation energy is due to Ceperley & Alder,         */
/* as param. by Perdew & Zunger,  Phys. Rev. B 23,5048 (1981) */
/**************************************************************/

double vxc_pz (double rho, double *exc)
{
  double au,bu,cu,du,aph,bt1,bt2,aa,bu1,cu1,du1,p13,pi4,c76,c43;
  double bpr, ec, ex, rrs, rs, rsl, rss, uc, ux, ax;
  au  =  0.0311; bu = -0.048; cu = 0.002;
  du  = -0.0116;
  aph = -0.1423; bt1 = 1.0529; bt2 = 0.3334;
  aa  = -3.0 / (4.0 * PI); bu1 = bu - au / 3.0;
  cu1 = 2.0 * cu / 3.0; du1 =  (2.0 * du - cu) / 3.0;
  p13 = 1.0 / 3.0; pi4 = 4.0 * PI; c76 = 7.0 / 6.0;
  c43 = 4.0 / 3.0;

  ax  = aa * pow (9.0 * PI / 4.0, p13);
  if (rho < 1.0e13) {
    *exc = 0.0;
    return (0.0);
  }
  else {
    rs  = pow (3.0 / (pi4 * rho), p13);
    rrs = sqrt (rs);
    rss = rs * rs;

    /* Exchange */
    ex = ax / rs;
    ux = c43 * ex;

    /* Correlation */
    if (rs >= 1.0) {
      bpr = bt1 * rrs + 1.0 + bt2 * rs;
      ec  = aph / bpr;
      uc  = (c76 * bt1 * rrs + 1.0 + c43 * bt2 * rs) * ec / bpr;
    }
    else {
      rsl = log (rs);
      ec  = du  * rs + bu  + au * rsl + cu  * rs * rsl;
      uc  = du1 * rs + bu1 + au * rsl + cu1 * rs * rsl;
    }
    *exc = rho * (ex + ec);
    return (ux + uc);
  }
}