source: pcp/src/excor.c@ f70c2a

/** \file excor.c
 * Exchange correlation energy.
 * 
 * The Exchange correlation energy is calculated via a local spin-density-approximation in
 * a parametrized manner. It is the sum of the exchange and the correlation, per electron
 * (thus, (discretely) integrated over the density). 
 * Due to the parametrization there is a huge number of functions evaluating small
 * expressions or auxiliary variables with its first and second derivatives
 * Calcf(), CalcDf(), CalcD2f() and CalcZeta(), CalcDZeta(), CalcD2Zeta(),
 * which then help to evaluate the energy CalcECr(), CalcEXrUP(), in the summing
 * CalcSECr(), CalcSEXr(), potential CalcVCr(), CalcVXrUP(), in its summation 
 * CalcSVCr(), CalcSVXr()
 * and derivatives CalcVVCr(), CalcVVXrUP(), in its summation CalcSVVCr(), CalcSVVXr().
 * All these are needed by the super functions evaluating exchange correlatiob
 * potential CalculateXCPotentialNoRT(), exchange correlation energy without RiemannTensor
 * CalculateXCEnergyNoRT() and with RiemannTensor CalculateXCEnergyUseRT(),
 * and the second derivative CalculateXCddEddt0NoRT() (needed for conjugate gradient).
 * During initilization the ExCor structure's entries is set to
 * these fixed parameters InitExchangeCorrelationEnergy().
 * 
 * 
  Project: ParallelCarParrinello
  \author Jan Hamaekers
 \date 2000

  File: excor.c
  $Id: excor.c,v 1.47 2007-03-29 13:37:25 foo Exp $
*/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdlib.h>
#include <stdio.h>
#include <stddef.h>
#include <math.h>

#include "data.h"
#include "mymath.h"
#include "run.h"
#include "myfft.h"
#include "errors.h"
#include "excor.h"
#include "helpers.h"

/** Calculate Wigner-Seitz-Radius \f$r_s\f$
 * \param *EC ExCor exchange correlation structure
 * \param p electron density
 * \return \f$(\frac{3}{4\pi \cdot p})^{1/3}   \qquad (2.30)\f$
 */
#ifdef HAVE_INLINE
inline double Calcrs(struct ExCor *EC, double p) {
#else
double Calcrs(struct ExCor *EC, double p) {
#endif
  if (fabs(p) < EC->epsilon0) return(0);
  //return (pow(3./(4.*PI*p),1./3.));
  return (cbrt(3./(4.*PI*p)));
}

/** Calculates exchange potential \f$V_{x}\f$.
 * \f[
 *    V_{x} = {\cal E}_{x} + \Bigr ( \frac{d{\cal E}_{x}}{dn} \Bigl )_{n=n(r_s)} n(r_s) \qquad (2.32)
 * \f]
 * \param *EC ExCor exchange correlation structure
 * \param rs Wigner-Seitz-Radius \f$r_s \f$, see Calcrs()
 * \return \f$-\frac{1}{2} (\frac{6}{\pi})^{2/3} \cdot \frac{1}{rs}\f$
 * \note The formula from CalcEXrUP() was used for the derivation.
 */
#ifdef HAVE_INLINE
inline static double CalcVXrUP(struct ExCor *EC, double rs) {
#else
static double CalcVXrUP(struct ExCor *EC, double rs) {
#endif
  if (fabs(rs) < EC->epsilon0) return(0); 
  return(-0.5*EC->fac6PI23/rs); // formula checked (16.5.06)
}

/** Calculates first derivative of exchange potential \f$\frac{\delta V_{x}}{\delta n}\f$.
 * \param *EC ExCor exchange correlation structure
 * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
 * \return \f$-\frac{2}{9\pi} (6\pi^2)^2/3 \cdot (rs)^2\f$
 * \sa CalcVXrUP()
 */
#ifdef HAVE_INLINE
inline static double CalcVVXrUP(struct ExCor *EC, double rs) {
#else
static double CalcVVXrUP(struct ExCor *EC, double rs) {
#endif
  return (-2./(9*PI)*EC->fac6PIPI23*rs*rs); // formula checked (16.5.06)
}

/** Calculates second derivative of exchange potential \f$\frac{\delta^2 V_{x}}{\delta n^2}\f$.
 * \param *EC ExCor exchange correlation structure
 * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
 * \return \f$\frac{16}{81} (6\pi^2)^{2/3} \cdot (rs)^5\f$
 * \sa CalcVXrUP(), CalcVVXrUP()
 */
#ifdef HAVE_INLINE
inline static double CalcVVVXrUP(struct ExCor *EC, double rs) {
#else
static double CalcVVVXrUP(struct ExCor *EC, double rs) {
#endif
  return (16./81.*EC->fac6PIPI23*rs*rs*rs*rs*rs); // formula checked (16.5.06)
}

/** Calculates exchange energy over density for homogeneous charge within given Wigner-Seitz-Radius, \f${\cal E}_x\f$.
 * The formula specified on the return statement is derived if into
 * \f[ \forall 0 \geq \zeta \geq 1: \quad 
 *      E_x = {\cal E}_x n = -\frac{3}{8} \Bigr (\frac{3(n_{up}+n_{down})}{\pi} \Bigl)^{1/3} 
 *        \bigr[ (1+\zeta)^{4/3}+ (1-\zeta)^{4/3} \bigl ] \cdot (n_{up}+n_{down}) \qquad (2.52)
 * \f] 
 * the definition of spin polarisation \f$\frac{n_{up} - n_{down}}{n_{up} + n_{down}}\f$
 * is inserted and the expression evaluated.
 * \param *EC ExCor exchange correlation structure
 * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
 * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
 * \return \f$-\frac{3}{8} (\frac{3}{2\pi})^{2/3} \Bigl [ (1+\zeta)^{\frac{4}{3}} + (1-\zeta)^{\frac{4}{3}}\Bigr] \cdot \frac{1}{rs}\f$
 */
#ifdef HAVE_INLINE
inline double CalcEXr(struct ExCor *EC, double rs, double zeta) {
#else
double CalcEXr(struct ExCor *EC, double rs, double zeta) {
#endif
  if (fabs(rs) < EC->epsilon0) return(0);
  //return (-3./8.*pow(3./(2.*PI),2./3.)/rs *(pow(1.+zeta,4./3.)+pow(1.-zeta,4./3.)));  // formula checked
  return (-3./8.*cbrt((3./(2.*PI))*(3./(2.*PI)))/rs *((1.+zeta)*cbrt(1.+zeta)+(1.-zeta)*cbrt(1.-zeta)));
}

/** Calculates exchange energy over density for homogeneous charge within given Wigner-Seitz-Radius, \f${\cal E}_x\f$.
 * The formula specified on the return statement is derived if into
 * \f[ \forall 0 \geq \zeta \geq 1: \quad 
 *      E_x = {\cal E}_x n = -\frac{3}{8} \Bigr (\frac{3(n_{up}+n_{down})}{\pi} \Bigl)^{1/3} 
 *        \bigr[ (1+\zeta)^{4/3}+ (1-\zeta)^{4/3} \bigl ] \cdot (n_{up}+n_{down}) \qquad (2.52)
 * \f] 
 * the definition of spin polarisation \f$\frac{n_{up} - n_{down}}{n_{up} + n_{down}}\f$
 * is inserted and the expression evaluated.
 * \param *EC ExCor exchange correlation structure
 * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
 * \return \f$-\frac{3}{8} (\frac{6}{\pi})^{2/3} \cdot \frac{1}{rs}\f$
 */
#ifdef HAVE_INLINE
inline static double CalcEXrUP(struct ExCor *EC, double rs) {
#else
static double CalcEXrUP(struct ExCor *EC, double rs) {
#endif
  if (fabs(rs) < EC->epsilon0) return(0); 
  return (-3./8.*EC->fac6PI23/rs);  // formula checked
}

/** Calculates correlation energy \f${\cal E}_c\f$ per electron for un-/polarised case.
 * Formula derived from Monte-Carlo-simulations by Ceperley and Adler [CA80],
 * parametrized by Perdew and Zunger [PZ81] for the polarised (\f$\zeta = 1\f$) or
 * unpolarised (\f$\zeta = 0\f$) case.
 * \param *EC ExCor exchange correlation structure
 * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
 * \param up UnPolarised: whether it's the polarised or unpolarised case
 * \return In the unpolarised case: <br>
 *         \f$r_s\f$>=1: \f$\gamma_{up} (1+\beta_{1,up}\sqrt{r_s}+\beta_{2,up}r_s)^{-1} = -0.1423\cdot(1+1.0529\sqrt{r_s}+0.3334r_s)^{-1} \qquad (2.31a)\f$<br>
 *         \f$r_s\f$<1:  \f$-0.0480+0.0311\ln(r_s)-0.0116r_s+0.0020r_s\ln(r_s) \qquad (2.31b)\f$<br>
 *         In the polarised case: <br>
 *         \f$r_s\f$>=1: \f$-0.0843\cdot(1+1.3981\sqrt{r_s}+0.2611r_s)^{-1}\f$<br>
 *         \f$r_s\f$<1:  \f$-0.0269+0.01555\ln(r_s)-0.0048r_s+0.0007r_s\ln(r_s)\f$<br>
 */
#ifdef HAVE_INLINE
inline double CalcECr(struct ExCor *EC, double rs, enum UnPolarised up) {
#else
double CalcECr(struct ExCor *EC, double rs, enum UnPolarised up) {
#endif
  double lrs;
  double res=0;
  if (rs >= 1) {
    res= (EC->gamma[up]/(1.+EC->beta_1[up]*sqrt(rs)+EC->beta_2[up]*rs));
    return (res);
  }
  if (rs <= EC->epsilon0) {
    res = 0;
    return(res);
  }
  lrs = log(rs);
  res = (EC->A[up]*lrs+EC->B[up]+EC->C[up]*rs*lrs+EC->D[up]*rs);
  return (res);
}

/** Calculates correlation potential \f$V_{c}\f$.
 * \f[
 *    V_{c} = {\cal E}_{c} + \Bigr ( \frac{d{\cal E}_{c}}{dn} \Bigl )_{n=n(r)} n(r) \qquad (2.32)
 * \f]
 * \param *EC ExCor exchange correlation structure
 * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
 * \param up UnPolarised
 * \return in the un-/polarised case: <br>
 *         \f$r_s\f$>=1: \f$\frac{{\cal E}_c}{1+\beta_1[up]\sqrt{r_s}+\beta_2[up]r_s} \cdot (1+\frac{7}{6}\beta_1[up]\sqrt{r_s} + \frac{4}{3}\beta_2[up]r_s)\f$<br>
 *         \f$r_s\f$=0: 0 <br>
 *         else:  \f$A[up]\cdot\log(rs) + (B[up] - \frac{1}{3} A[up]) + \frac{2}{3} C[up]\cdot rs \log(rs)+ \frac{1}{3} (2 D[up]-C[up]) \cdot rs\f$
 */
#ifdef HAVE_INLINE
inline static double CalcVCr(struct ExCor *EC, double rs, enum UnPolarised up) {
#else
static double CalcVCr(struct ExCor *EC, double rs, enum UnPolarised up) {
#endif
  double srs,lrs;
  if (rs >= 1) {
    srs = sqrt(rs);
    return ((EC->gamma[up]/(1.+EC->beta_1[up]*srs+EC->beta_2[up]*rs))*(1.+(7./6.)*EC->beta_1[up]*srs+(4./3.)*EC->beta_2[up]*rs)/(1.+EC->beta_1[up]*srs+EC->beta_2[up]*rs)); // formula checked  (15.5.06)
  }
  if (rs <= EC->epsilon0) return (0);
  lrs = log(rs);
  return (EC->A[up]*lrs+(EC->B[up]-(1./3.)*EC->A[up])+(2./3.)*EC->C[up]*rs*lrs+(1./3.)*(2.*EC->D[up]-EC->C[up])*rs); // formula checked  (15.5.06)
}

/** Calculates first derivative of correlation potential \f$\frac{\delta V_{c}}{\delta n}\f$.
 * \f[
 *  \frac{\delta V_{c}}{\delta n} = \frac{\delta^2 {\cal E}_{c}}{\delta n^2}
 * \f]
 * \param *EC ExCor exchange correlation structure
 * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
 * \param up UnPolarised
 * \return In the un-/polarised case: <br>
 *        \f$r_s\f$>=1: \f$\frac{{\cal E}_c[up] \pi}{27\cdot(1+\beta_1[up]\sqrt{r_s}+\beta_2[up]r_s)} \Bigr ( 5\beta_1[up]r_s^{7/2} + 8\beta_2[up]r_s^4
 *            + \frac{1}{1+\beta_1[up]\sqrt{r_s}+\beta_2[up]r_s} (8\beta_1[up]\beta_2[up]r_s^{9/2} + 2\beta_1^2[up]r_s^4 + 8\beta_2^2[up]r_s^5) \Bigl)\f$<br>
 *        \f$r_s\f$< 1: \f$-\frac{4\pi r_s^3}{9} \bigr (A[up] + \frac{r_s}{3} \cdot(C[up]+2C[up]\log(r_s)+2D[up]) \bigl)\f$<br>
 * \sa CalcVCr()
 */
#ifdef HAVE_INLINE
inline static double CalcVVCr(struct ExCor *EC, double rs, enum UnPolarised up) {
#else
static double CalcVVCr(struct ExCor *EC, double rs, enum UnPolarised up) {
#endif
  double eps,deps;
  if (rs <= EC->epsilon0) return (0);
  if (rs >= 1) {
    deps = 1.+EC->beta_1[up]*sqrt(rs)+EC->beta_2[up]*rs;
    eps = EC->gamma[up]/deps;
    //return (eps*PI/(27.*deps)*(8.*EC->beta_2[up]*rs*rs*rs*rs+5.*EC->beta_1[up]*pow(rs,7./2.)+(1./deps)*(8.*EC->beta_1[up]*EC->beta_2[up]*pow(rs,9./2.)+2.*EC->beta_1[up]*EC->beta_1[up]*rs*rs*rs*rs+8.*EC->beta_2[up]*EC->beta_2[up]*rs*rs*rs*rs*rs)));
    //return (eps*PI/(27.*deps*deps)*(8.*EC->beta_2[up]*rs*rs*rs*rs+5.*EC->beta_1[up]*pow(rs,7./2.)+21.*EC->beta_1[up]*EC->beta_2[up]*pow(rs,9./2.)+7.*EC->beta_1[up]*EC->beta_1[up]*rs*rs*rs*rs+16.*EC->beta_2[up]*EC->beta_2[up]*rs*rs*rs*rs*rs)); //formula checked (15.05.06)
    return (eps*PI/(27.*deps*deps)*(8.*EC->beta_2[up]*rs*rs*rs*rs+5.*EC->beta_1[up]*(sqrt(rs)*rs*rs*rs)+21.*EC->beta_1[up]*EC->beta_2[up]*(sqrt(rs)*rs*rs*rs*rs)+7.*EC->beta_1[up]*EC->beta_1[up]*rs*rs*rs*rs+16.*EC->beta_2[up]*EC->beta_2[up]*rs*rs*rs*rs*rs)); //formula checked (15.05.06)
  }
  return (-4.*PI*rs*rs*rs/9.*(EC->A[up]+rs/3.*(EC->C[up]+2.*EC->C[up]*log(rs)+2.*EC->D[up]))); //formula checked (23.05.06)
}

/** Calculates second derivative of correlation potential \f$\frac{\delta^2 V_{c}}{\delta n^2}\f$.
 * \f[
 *  \frac{\delta^2 V_{c}}{\delta n^2} = \frac{\delta^3 {\cal E}_{c}}{\delta n^3}
 * \f]
 * \param *EC ExCor exchange correlation structure
 * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
 * \param up UnPolarised
 * \return In the un-/polarised case: <br>
 *        \f$r_s\f$>=1: \f$-\frac{2\pi^2 {\cal E}_c[up]}{243\cdot(1+\beta_1[up]\sqrt{r_s}+\beta_2[up]r_s)^4} 
 *            \Bigr ( 35\beta_1[up]r_s^{13/2} + r_s^7(64\beta_2[up]+76\beta_1^2[up])
 *              + r_s^{15/2} (35\beta^3_1[up]+243\beta_1[up]\beta_2[up])
 *              + r_s^8 (176\beta^2_2[up]+140\beta^2_1[up]\beta_2[up])
 *              + r_s^{17/2} (175\beta_1[up]\beta^2_2[up])
 *              + r_s^9 (64\beta^3_2[up])      
 *            \Bigl)\f$<br>
 *        \f$r_s\f$< 1: \f$\frac{16\pi^2 r_s^6}{243} \bigr (9A[up] + r_s\cdot(6C[up]+8C[up]\log(r_s)+8D[up]) \bigl)\f$<br>
 * \sa CalcVCr()
 */
#ifdef HAVE_INLINE
inline static double CalcVVVCr(struct ExCor *EC, double rs, enum UnPolarised up) {
#else
static double CalcVVVCr(struct ExCor *EC, double rs, enum UnPolarised up) {
#endif
  double eps,deps;
  double beta11,beta12,beta13,beta21,beta22,beta23;
  if (rs <= EC->epsilon0) return (0);
  if (rs >= 1) {
    deps = 1.+EC->beta_1[up]*sqrt(rs)+EC->beta_2[up]*rs;
    eps=EC->gamma[up]/deps;
    beta11 = EC->beta_1[up];
    beta12 = beta11*EC->beta_1[up];
    beta13 = beta12*EC->beta_1[up];
    beta21 = EC->beta_2[up];
    beta22 = beta21*EC->beta_2[up];
    beta23 = beta22*EC->beta_2[up];
//    return (-2.*eps*PI*PI/(243.*deps*deps*deps)*(
//            pow(rs,13./2.) * (35.*beta11)
//          + pow(rs,7.) * (64.*beta21 + 76.*beta12)
//          + pow(rs,15./2.) * (35.*beta13 + 234.*beta11*beta21)
//          + pow(rs,8.) * (176.*beta22 + 140.*beta12*beta21)
//          + pow(rs,17./2.) * (175.*beta11*beta22)
//          + pow(rs,9.) * (64.*beta23) ));
    return (-2.*eps*PI*PI/(243.*deps*deps*deps)*(
            rs*rs*rs*rs*rs*rs*(
            sqrt(rs) *((35.*beta11)  // fractial ones
          + rs *((35.*beta13 + 234.*beta11*beta21)
          + rs * (175.*beta11*beta22))
            )
          + rs * ((64.*beta21 + 76.*beta12) // integer ones
          + rs * ((176.*beta22 + 140.*beta12*beta21)
          + rs * ((64.*beta23) ))))));
  }
  //return (16.*PI*PI*pow(rs,6)/243. * (9.*EC->A[up] + rs*(8.*EC->C[up]*log(rs) + 6.*EC->C[up] + 8.*EC->D[up])));
  return (16.*PI*PI*rs*rs*rs*rs*rs*rs/243. * (9.*EC->A[up] + rs*(8.*EC->C[up]*log(rs) + 6.*EC->C[up] + 8.*EC->D[up])));
}

/** Calculates spin polarisation \f$\zeta\f$
 * \param *EC ExCor exchange correlation structure
 * \param pUp SpinUp electron density
 * \param pDown SpinDown electron density
 * \return \f$\frac{pUp - pDown}{pUp + pDown}\f$
 * \note zeta is never less than ExCor::epsilon0
 */
#ifdef HAVE_INLINE
inline double CalcZeta(struct ExCor *EC, double pUp, double pDown) {
#else
double CalcZeta(struct ExCor *EC, double pUp, double pDown) {
#endif
  if (fabs(pUp+pDown) < EC->epsilon0) return(0);
  return ((pUp - pDown)/(pUp + pDown));
}

/** Calculates derivative \f$\frac{\delta \zeta}{\delta p}\f$
 * \param *EC ExCor exchange correlation structure
 * \param ST SpinType
 * \param pUp SpinUp density
 * \param pDown SpinDown density
 * \return p is pUp or pDown: \f$\pm 2 \frac{p}{ (pUp+pDown)^2 }\f$
 * \sa CalcZeta()
 */
#ifdef HAVE_INLINE
inline static double CalcDzeta(struct ExCor *EC, enum SpinType ST, double pUp, double pDown) {
#else
static double CalcDzeta(struct ExCor *EC, enum SpinType ST, double pUp, double pDown) {
#endif
  double res = 0;
  if (fabs(pUp+pDown) < EC->epsilon0) return(0);
  //EC = EC;
  switch(ST) {
  case SpinUp:
    res = 2.*pDown/(pUp+pDown)/(pUp+pDown); // formula checked
    break;
  case SpinDown:
    res = -2.*pUp/(pUp+pDown)/(pUp+pDown); // formula checked
    break;
  case SpinDouble:
    res = 0;
    break;
  }
  return (res);
}

/** Calculates second derivative \f$\frac{\delta^2 \zeta}{\delta p^2}\f$
 * \param *EC ExCor exchange correlation structure
 * \param ST SpinType
 * \param pUp SpinUp density
 * \param pDown SpinDown density
 * \return p is pUp or pDown: \f$\pm 4 \frac{p}{ (pUp+pDown)^3 }\f$
 * \sa CalcZeta(), CalcDZeta()
 */
#ifdef HAVE_INLINE
inline static double CalcD2zeta(struct ExCor *EC, enum SpinType ST, double pUp, double pDown) {
#else
static double CalcD2zeta(struct ExCor *EC, enum SpinType ST, double pUp, double pDown) {
#endif
  double res = 0;
  if (fabs(pUp+pDown) < EC->epsilon0) return(0);
  //EC = EC;
  switch(ST) {
  case SpinUp:
    res = -4.*pDown/(pUp+pDown)/(pUp+pDown)/(pUp+pDown);
    break;
  case SpinDown:
    res = 4.*pUp/(pUp+pDown)/(pUp+pDown)/(pUp+pDown);
    break;
  case SpinDouble:
    res = 0;
    break;
  }
  return (res);
}

/** Calculates third derivative \f$\frac{\delta^3 \zeta}{\delta p^3}\f$
 * \param *EC ExCor exchange correlation structure
 * \param ST SpinType
 * \param pUp SpinUp density
 * \param pDown SpinDown density
 * \return p is pUp or pDown: \f$\pm -12 \frac{p}{ (pUp+pDown)^4 }\f$
 * \sa CalcZeta(), CalcDZeta(), CalcD2Zeta()
 */
#ifdef HAVE_INLINE
inline static double CalcD3zeta(struct ExCor *EC, enum SpinType ST, double pUp, double pDown) {
#else
static double CalcD3zeta(struct ExCor *EC, enum SpinType ST, double pUp, double pDown) {
#endif
  double res = 0;
  if (fabs(pUp+pDown) < EC->epsilon0) return(0);
  //EC = EC;
  switch(ST) {
  case SpinUp:
    res = +12.*pDown/(pUp+pDown)/(pUp+pDown)/(pUp+pDown)/(pUp+pDown);
    break;
  case SpinDown:
    res = -12.*pUp/(pUp+pDown)/(pUp+pDown)/(pUp+pDown)/(pUp+pDown);
    break;
  case SpinDouble:
    res = 0;
    break;
  }
  return (res);
}

/** Calculates auxiliary factor \f$f(\zeta)\f$.
 * This auxiliary factor is needed in the parametrized evaluation of the
 * full spin-polarised correlation energy \f${\cal E}_c\f$ (see section 2.3.6)
 * \param *EC ExCor exchange correlation structure
 * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
 * \return \f$\frac{(1+\zeta)^{4/3} + (1-\zeta)^{4/3} - 2} {2^{4/3}-2}\f$
 */
#ifdef HAVE_INLINE
inline static double Calcf(struct ExCor *EC, double zeta) {
#else
static double Calcf(struct ExCor *EC, double zeta) {
#endif
  double res =0;
  EC = EC;
  if (fabs(zeta) < EC->epsilon0) return(0.); // unpolarised case
  //res = ((pow(1.+zeta,4./3.)+pow(1.-zeta,4./3.)-2.)/(EC->fac243-2.));
  res = ((cbrt(1.+zeta)*(1.+zeta)+cbrt(1.-zeta)*(1.-zeta)-2.)/(EC->fac243-2.));
  return (res);
}

/** Calculates the derivative \f$\frac{\delta f}{\delta \zeta}\f$.
 * \param *EC ExCor exchange correlation structure
 * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
 * \return \f$\frac{4}{3} \frac{(\zeta+1)^{1/3} - (1-\zeta)^{1/3}} {2^{4/3}-2}\f$
 * \sa Calcf()
 */
#ifdef HAVE_INLINE
inline static double CalcDf(struct ExCor *EC, double zeta) {
#else
static double CalcDf(struct ExCor *EC, double zeta) {
#endif
  if (fabs(zeta) < EC->epsilon0) return(0.); // unpolarised case
  //return((4./3.)*(pow(zeta+1.,1./3.)-pow(1.-zeta,1./3.))/(EC->fac243-2.)); // formula checked (16.5.06)
  return((4./3.)*(cbrt(zeta+1.)-cbrt(1.-zeta))/(EC->fac243-2.)); // formula checked (16.5.06)
}

/** Calculates second derivative \f$\frac{\delta^2 f}{\delta \zeta^2}\f$.
 * \param *EC ExCor exchange correlation structure
 * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
 * \return \f$\frac{4}{9} \frac{(\zeta+1)^{-2/3}-(1-\zeta)^{-2/3} }{2^{4/3}-2}\f$
 * \sa Calcf(), CalcDf()
 */
#ifdef HAVE_INLINE
inline static double CalcD2f(struct ExCor *EC, double zeta) {
#else
static double CalcD2f(struct ExCor *EC, double zeta) {
#endif
  if (fabs(zeta) < EC->epsilon0) return(0.); // unpolarised case
  if (1.-fabs(zeta) < EC->epsilon0) return (0.);  // second derivative not defined for these cases (which are needed!) 
  //if (1-fabs(zeta) < EC->epsilon0) { fprintf(stderr,"CalcD2f: zeta = %lg\n",zeta); return(0); }
  //return((4./9.)*(pow(zeta+1.,-2./3.)+pow(1.-zeta,-2./3.))/(EC->fac243-2.));  // formula checked (16.5.06)
  return((4./9.)*(1./cbrt((zeta+1.)*(zeta+1.))+1./cbrt((1.-zeta)*(1.-zeta)))/(EC->fac243-2.));  // formula checked (16.5.06)
}

/** Calculates third derivative \f$\frac{\delta^3 f}{\delta \zeta^3}\f$.
 * \param *EC ExCor exchange correlation structure
 * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
 * \return \f$\frac{8}{27} \frac{-(\zeta+1)^{-5/3}+(1-\zeta)^{-5/3} }{2^{4/3}-2}\f$
 * \sa Calcf(), CalcDf()
 */
#ifdef HAVE_INLINE
inline static double CalcD3f(struct ExCor *EC, double zeta) {
#else
static double CalcD3f(struct ExCor *EC, double zeta) {
#endif
  if (fabs(zeta) < EC->epsilon0) return(0.); // unpolarised case
  //if (1-fabs(zeta) < EC->epsilon0) return (0.);  // second derivative not defined for these cases (which are needed!) 
  if (1-fabs(zeta) < EC->epsilon0) { fprintf(stderr,"CalcD2f: zeta = %lg\n",zeta); return(0.); }
  //return((-8./27.)*(pow(zeta+1.,-5./3.)-pow(1.-zeta,-5./3.))/(EC->fac243-2.));  // formula checked (16.5.06)
  return((-8./27.)*(1/((zeta+1)*cbrt((zeta+1.)*(zeta+1.)))-1/((zeta+1)*cbrt((1.-zeta)*(1.-zeta))))/(EC->fac243-2.));  // formula checked (16.5.06)
}
  
/** Calculates correlation energy with free spin-polarisation \f$E_c (n,\zeta)\f$.
 * \param *EC exchange correlation structure
 * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
 * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
 * \param p electron density \f$n\f$
 * \return \f$0 \leq \zeta \leq 1:\quad E_c (n,\zeta) = {\cal E}_c (n,\zeta)\cdot n = \bigr ({\cal E}_c(n,0) + [{\cal E}_c(n,1) - {\cal E}_c(n,0)] f(\zeta)\bigl ) \cdot n \qquad (2.3.6)\f$
 * \sa CalcECr()
 */
#ifdef HAVE_INLINE
inline double CalcSECr(struct ExCor *EC, double rs, double zeta, double p) {
#else
double CalcSECr(struct ExCor *EC, double rs, double zeta, double p) {
#endif
  double res =0;
  double unpol, pol;
  unpol = CalcECr(EC,rs,unpolarised); 
  pol = CalcECr(EC,rs,polarised); 
  res = (unpol+Calcf(EC,zeta)*(pol-unpol))*p;
  return (res);
}

/** Calculates SpinType-dependently correlation potential with free spin-polarisation \f$V_c (n,\zeta)\f$.
 * \f[
 *  V_c = {\cal E}_c + \frac{\delta{\ E}_c}{\delta n}n = V_c(n,0) + [V_c(n,1)-V_c(n,0)]f(\zeta) 
 *      + [{\cal E}_c(n,1)-{\cal E}_c(N,0)]\frac{\delta f(\zeta)}{\delta \zeta}\frac{\delta\zeta}{\delta n}
 * \f]
 * \param *EC ExCor exchange correlation structure
 * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
 * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
 * \param ST SpinType
 * \return \f$V_c(n,0) + [V_c(n,1)-V_c(n,0)]f(\zeta) + [{\cal E}_c(n,1)-{\cal E}_c(N,0)](\pm1-\zeta)\frac{\delta f(\zeta)}{\delta \zeta}\f$
 */
#ifdef HAVE_INLINE
inline static double CalcSVCr(struct ExCor *EC, double rs, double zeta, enum SpinType ST) {
#else
static double CalcSVCr(struct ExCor *EC, double rs, double zeta, enum SpinType ST) {
#endif
  double res = 0;
  double VCR_unpol = CalcVCr(EC,rs,unpolarised);
  switch (ST) {
  case SpinUp:
    res = VCR_unpol + Calcf(EC,zeta)*(CalcVCr(EC,rs,polarised)-VCR_unpol) + (CalcECr(EC,rs,polarised)-CalcECr(EC,rs,unpolarised))*(1.-zeta)*CalcDf(EC,zeta); // formula checked (15.5.06)
    break;
  case SpinDown:
    res = VCR_unpol + Calcf(EC,zeta)*(CalcVCr(EC,rs,polarised)-VCR_unpol) + (CalcECr(EC,rs,polarised)-CalcECr(EC,rs,unpolarised))*(-1.-zeta)*CalcDf(EC,zeta); // formula checked (15.5.06)
    break;
  case SpinDouble:
    res = VCR_unpol; /*EC->fac1213**/
    break;
  }
  return (res);
}

/** Calculates complete exchange energy \f$E_x (n)\f$.
 * \param *EC exchange correlation structure
 * \param rsUp Wigner-Seitz-Radius \f$r_s\f$ for pUp
 * \param rsDown Wigner-Seitz-Radius \f$r_s\f$ for pDown
 * \param pUp SpinUp electron density
 * \param pDown SpinDown electron density
 * \return \f$E_x = {\cal E}_x \cdot n\f$
 * \sa CalcEXrUP(), CalculateXCEnergyNoRT()
 */
#ifdef HAVE_INLINE
inline double CalcSEXr(struct ExCor *EC, double rsUp, double rsDown, double pUp, double pDown) {
#else
double CalcSEXr(struct ExCor *EC, double rsUp, double rsDown, double pUp, double pDown) {
#endif
  return(CalcEXrUP(EC,rsUp)*pUp+CalcEXrUP(EC,rsDown)*pDown);
}

/** Calculates SpinType-dependently exchange potential \f$V_x (n)\f$.
 * Essentially, just CalcVXrUP() is called, even in SpinDouble case
 * \param *EC ExCor exchange correlation structure
 * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
 * \param ST SpinType
 * \return CalcVXrUP()
 */
#ifdef HAVE_INLINE
inline static double CalcSVXr(struct ExCor *EC, double rs, enum SpinType ST) {
#else
static double CalcSVXr(struct ExCor *EC, double rs, enum SpinType ST) {
#endif
  double res = 0;
  switch (ST) {
  case SpinUp:
  case SpinDown:
    res = CalcVXrUP(EC,rs);
    break;
  case SpinDouble:
    res = EC->fac1213*CalcVXrUP(EC,rs);
    break;
  }
  return (res);
}

/** Calculates SpinType-dependently derivative of exchange potential \f$\frac{\delta V_x}{\delta n}\f$.
 * Essentially, just CalcVVXrUP() is called, even in SpinDouble case 
 * \param *EC ExCor exchange correlation structure
 * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
 * \param ST SpinType
 * \return CalcVVXrUP()
 */
#ifdef HAVE_INLINE
inline double CalcSVVXr(struct ExCor *EC, double rs, enum SpinType ST) {
#else
double CalcSVVXr(struct ExCor *EC, double rs, enum SpinType ST) {
#endif
  double res = 0;
  switch (ST) {
  case SpinUp:
  case SpinDown:
    res = CalcVVXrUP(EC,rs);
    break;
  case SpinDouble:
    res = EC->fac1213*CalcVVXrUP(EC,rs);
    break;
  }
  return (res);
}

/** Calculates SpinType-dependently second derivative of exchange potential \f$\frac{\delta^2 V_x}{\delta n^2}\f$.
 * Essentially, just CalcVVVXrUP() is called, even in SpinDouble case 
 * \param *EC ExCor exchange correlation structure
 * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
 * \param ST SpinType
 * \return CalcVVVXrUP()
 */
#ifdef HAVE_INLINE
inline double CalcSVVVXr(struct ExCor *EC, double rs, enum SpinType ST) {
#else
double CalcSVVVXr(struct ExCor *EC, double rs, enum SpinType ST) {
#endif
  double res = 0;
  switch (ST) {
  case SpinUp:
  case SpinDown:
    res = CalcVVVXrUP(EC,rs);
    break;
  case SpinDouble:
    res = EC->fac1213*CalcVVVXrUP(EC,rs);
    break;
  }
  return (res);
}

/** Calculates SpinType-dependently first derivative of correlation potential with free spin-polarisation \f$\frac{\delta V_c}{\delta n} (n,\zeta)\f$.
 * \param *EC ExCor exchange correlation structure
 * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
 * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
 * \param ST SpinType
 * \param pUp SpinUp density
 * \param pDown SpinDown density
 * \return \f$\frac{\delta V_c}{\delta n} (n,0) + f(\zeta) \bigr ( \frac{\delta V_c}{\delta n} (n,1) - \frac{\delta V_c}{\delta n} (n,0)\bigl )
 *  + 2 \cdot \frac{\delta f}{\delta \zeta} \frac{\delta \zeta}{\delta n} \bigr ( V_c (n,1) - V_c (n,0)\bigl ) + \ldots\f$<br>
 *  \f$\ldots + ({\cal E}_c (n,1) - {\cal E}_c (n,0) ) n (\frac{\delta^2 \zeta}{\delta n^2}\frac{\delta f}{\delta \zeta} + \bigr(\frac{\delta \zeta}{\delta n}\bigl)^2
 * \frac{\delta f}{\delta \zeta})\f$
 * \sa CalcSVCr()
 */
#ifdef HAVE_INLINE
inline double CalcSVVCr(struct ExCor *EC, double rs, double zeta, enum SpinType ST, double pUp, double pDown) {
#else
double CalcSVVCr(struct ExCor *EC, double rs, double zeta, enum SpinType ST, double pUp, double pDown) {
#endif
  double res = 0;
  double VVCR_pol, VVCR_unpol, VCR_pol, VCR_unpol, ECr_pol, ECr_unpol, f, Df, D2f, DZeta, D2Zeta;
  switch (ST) {
  case SpinUp:
  case SpinDown:
    // store some function calls for faster and easierly understandable operation
    VVCR_unpol = CalcVVCr(EC,rs,unpolarised);
    VVCR_pol = CalcVVCr(EC,rs,polarised);
    ECr_pol = CalcECr(EC,rs,polarised);
    ECr_unpol = CalcECr(EC,rs,unpolarised);
    f = Calcf(EC,zeta);
    Df = CalcDf(EC,zeta);
    D2Zeta = CalcD2zeta(EC,ST,pUp,pDown);
    if (CalcDzeta(EC,ST,pUp,pDown) != 0.0) {
      VCR_unpol = CalcVCr(EC,rs,unpolarised);
      VCR_pol = CalcVCr(EC,rs,polarised);
      D2f = CalcD2f(EC,zeta);
      DZeta = CalcDzeta(EC,ST,pUp,pDown);
      res = VVCR_unpol + f*(VVCR_pol-VVCR_unpol) + 2.*DZeta*Df*(VCR_pol-VCR_unpol) + (ECr_pol-ECr_unpol)*(pUp+pDown)*(D2Zeta*Df + DZeta*DZeta*D2f); //formula checked (15.5.06)
    } else {
      res = VVCR_unpol + f*(VVCR_pol-VVCR_unpol) + (ECr_pol-ECr_unpol)*(pUp+pDown)*(D2Zeta*Df);  //formula checked (15.5.06)
    }
    break;
  case SpinDouble:
    res = CalcVVCr(EC,rs,unpolarised);
    break;
  }
  return (res);
}

/** Calculates SpinType-dependently second derivative of correlation potential with free spin-polarisation \f$\frac{\delta V_c}{\delta n} (n,\zeta)\f$.
 * \param *EC ExCor exchange correlation structure
 * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
 * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
 * \param ST SpinType
 * \param pUp SpinUp density
 * \param pDown SpinDown density
 * \return \f$\frac{\delta^2 V_c}{\delta n^2} (n,0) + f(\zeta) \bigr ( \frac{\delta^2 V_c}{\delta n^2} (n,1) - \frac{\delta^2 V_c}{\delta n^2} (n,0)\bigl )
 *  + 3 \cdot \frac{\delta f}{\delta \zeta} \frac{\delta \zeta}{\delta n} \bigr ( \frac{\delta V_c}{\delta n} (n,1) - \frac{\delta V_c}{\delta n} (n,0) \bigl ) + \ldots\f$<br>
 *  \f$\ldots + 3 \cdot (\frac{\delta^2 f}{\delta \zeta^2} (\frac{\delta \zeta}{\delta n})^2 + \frac{\delta f}{\delta \zeta} \frac{\delta^2 \zeta}{\delta n^2} ) \bigr ( V_c (n,1) - V_c (n,0)\bigl ) +
 *  + ({\cal E}_c (n,1) - {\cal E}_c (n,0) ) n \bigr (\frac{\delta^3 \zeta}{\delta n^3}(\frac{\delta f}{\delta \zeta})^3 + 3\cdot \frac{\delta^2 f}{\delta \zeta^2}\frac{\delta \zeta}{\delta n}\frac{\delta^2 \zeta}{\delta n^2} + \frac{\delta f}{\delta \zeta}\frac{\delta^3 \zeta}{\delta n^3}\bigl)^2\f$
 * \sa CalcSVCr(), CalcSVVCr()
 */
#ifdef HAVE_INLINE
inline double CalcSVVVCr(struct ExCor *EC, double rs, double zeta, enum SpinType ST, double pUp, double pDown) {
#else
double CalcSVVVCr(struct ExCor *EC, double rs, double zeta, enum SpinType ST, double pUp, double pDown) {
#endif
  double res = 0;
  double VVVCR_pol, VVVCR_unpol, VVCR_pol, VVCR_unpol, VCR_pol, VCR_unpol, ECr_pol, ECr_unpol, f, Df, D2f, D3f, DZeta, D2Zeta, D3Zeta;
  switch (ST) {
  case SpinUp:
  case SpinDown:
    // store some function calls for faster and easierly understandable operation
    VVVCR_unpol = CalcVVVCr(EC,rs,unpolarised);
    VVVCR_pol = CalcVVVCr(EC,rs,polarised);
    VCR_unpol = CalcVCr(EC,rs,unpolarised);
    VCR_pol = CalcVCr(EC,rs,polarised);
    ECr_pol = CalcECr(EC,rs,polarised);
    ECr_unpol = CalcECr(EC,rs,unpolarised);
    f = Calcf(EC,zeta);
    Df = CalcDf(EC,zeta);
    D2Zeta = CalcD2zeta(EC,ST,pUp,pDown);
    D3Zeta = CalcD3zeta(EC,ST,pUp,pDown);
    if (CalcDzeta(EC,ST,pUp,pDown) != 0.0) {
      VVCR_unpol = CalcVVCr(EC,rs,unpolarised);
      VVCR_pol = CalcVVCr(EC,rs,polarised);
      D2f = CalcD2f(EC,zeta);
      D3f = CalcD3f(EC,zeta);
      DZeta = CalcDzeta(EC,ST,pUp,pDown);
      res = VVVCR_unpol + f*(VVVCR_pol-VVVCR_unpol) + 3.*DZeta*Df*(VVCR_pol-VVCR_unpol) + 3.*(DZeta*DZeta*D2f + Df*D2Zeta)*(VCR_pol-VCR_unpol) + (ECr_pol-ECr_unpol)*(pUp+pDown)*(D3Zeta*Df + 3*D2f*DZeta*D2Zeta + DZeta*DZeta*DZeta*D3f); //formula checked (16.5.06)
    } else {
      res = VVVCR_unpol + f*(VVVCR_pol-VVVCR_unpol) + 3.*Df*D2Zeta*(VCR_pol-VCR_unpol) + (ECr_pol-ECr_unpol)*(pUp+pDown)*(D3Zeta*Df);  //formula checked (16.5.06)
    }
    break;
  case SpinDouble:
    res = CalcVVVCr(EC,rs,unpolarised);
    break;
  }
  return (res);
}

/** SpinType-dependent calculation of exchange and correlation energy with free spin-polarisation without riemann tensor.
 * Discretely does the integration of \f${\cal E}_{xc} (n,\zeta) \cdot n\f$ on
 * the radial mesh of the respective real densities.<br>
 * Takes either whole density or separated for each SpinType and calls
 * Calcrs(), then summing each CalcSEXr() for each Density::LocalSizeR 
 * and times factor RunStruct::XCEnergyFactor / LatticeLevel::MaxN
 * \param *P Problem at hand
 */
void CalculateXCEnergyNoRT(struct Problem *P)
{
  struct Lattice *Lat = &P->Lat;
  struct Energy *E = Lat->E;
  struct PseudoPot *PP = &P->PP;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct Psis *Psi = &Lat->Psi;
  struct Density *Dens = Lev0->Dens;
  struct ExCor *EC = &P->ExCo;
  double SumEc = 0;
  //double SumEx_GC = 0.; // Gradient correction part according to Becke '92
  double rs, p = 0.0, pUp = 0.0, pDown = 0.0, zeta, rsUp, rsDown, SumEx=0.0;
  double Factor = R->XCEnergyFactor/Lev0->MaxN;
  //double Dp, DpUp, DpDown;
  int i;
  
  for (i = 0; i < Dens->LocalSizeR; i++) {  // for each node in radial mesh
    // put (corecorrected) densities in p, pUp, pDown  
    p = Dens->DensityArray[TotalDensity][i];
    //Dp = DensityGradient(Dens->DensityArray[TotalDensity], i, Lev0, Lat);
    if (R->CurrentMin > UnOccupied)
    if (PP->corecorr == CoreCorrected)
      p += Dens->DensityArray[CoreWaveDensity][i];
    switch (Psi->PsiST) {
    default:
    case SpinDouble:
      pUp = 0.5*p; 
      pDown = 0.5*p;
      //DpUp = 0.5*Dp; 
      //DpDown = 0.5*Dp;
      break;
    case SpinUp:  
    case SpinDown:
      pUp = Dens->DensityArray[TotalUpDensity][i];
      pDown = Dens->DensityArray[TotalDownDensity][i];
      //DpUp = DensityGradient(Dens->DensityArray[TotalUpDensity], i, Lev0, Lat);
      //DpDown = DensityGradient(Dens->DensityArray[TotalDownDensity], i, Lev0, Lat);
      if (PP->corecorr == CoreCorrected) {
                                pUp += 0.5*Dens->DensityArray[CoreWaveDensity][i];
                                pDown += 0.5*Dens->DensityArray[CoreWaveDensity][i];
      }
      break;
    }
    // set all to zero if one of them is negative
    if ((p < 0) || (pUp < 0) || (pDown < 0)) {
      /*fprintf(stderr,"index %i pc %g p %g\n",i,Dens->DensityArray[CoreWaveDensity][i],p);*/
      p = 0.0; 
      pUp = 0.0;
      pDown = 0.0;
    }
    // Calculation with the densities and summation
    rs = Calcrs(EC,p);
    rsUp = Calcrs(EC,pUp);
    rsDown = Calcrs(EC,pDown);
    SumEx += CalcSEXr(EC, rsUp, rsDown, pUp, pDown);
    zeta = CalcZeta(EC,pUp,pDown);
    SumEc += CalcSECr(EC, rs, zeta, p);
    //SumEx_GC += CalcSE_GC(EC, pUp, DpUp);
    //SumEx_GC += CalcSE_GC(EC, pDown, DpDown);
  }
  E->AllLocalDensityEnergy[CorrelationEnergy] = Factor*SumEc;
  E->AllLocalDensityEnergy[ExchangeEnergy] = Factor*(SumEx); // - SumEx_GC);
}

/** Returns \f$y = \sinh^{-1}(x)\f$.
 * Recall that \f$ x = sinh (y) =  \frac{\exp(y)-\exp(-y)}{2}\f$. Then, solve the in \f$\exp(y)\f$ quadratic
 * equation: \f$\exp(2y) - 1 - 2x\exp(y) = 0\f$ by pq-formula.
 * \param arg argument \f$x\f$
 * \return \f$\sinh^{-1}(x) = \ln(x+\sqrt{x^2+1})\f$
 */
inline static double sinh_inverse(double arg)
{
  return log(arg+sqrt(arg*arg+1));
}

/** Gradient correction according to Becke '92.
 * \param EC ExCor structure
 * \param p \f$\rho\f$ local density
 * \param Dp \f$|\nabla \rho|\f$ local magnitude of density gradient
 * \return \f$b \rho^{4/3} \frac{x^2}{(1+6bx \sinh^-1 x}\f$ with \f$x = \frac{|\nabla \rho|}{\rho^{4/3}}\f$
 */
double CalcSE_GC(struct ExCor *EC, double p, double Dp)
{
  double res = 0.;
  double b = 0.0042;  // in atomic units, empirical constant fitted to noble gas atoms
  double x = Dp/pow(p, 4./3.);
  res = b * pow(p, 4./3.) * x/(1+6.*b*x*sinh_inverse(x));
  
  return res;
}

/** Evaluates magnitude of gradient by simple 3-star.
 * \param *density density array
 * \param i current node
 * \param *Lev LatticeLevel structure for number of nodes per axis
 * \param *Lat Lattice structure for axis lengths
 * \return \f$|\nabla\rho|\f$
 */
double DensityGradient(fftw_real *density, int i, struct LatticeLevel *Lev, struct Lattice *Lat)
{
  double res=0., right_diff, left_diff;
  int neighbour[NDIM], nodes[NDIM];
  nodes[0] = Lev->Plan0.plan->local_nx;
  nodes[1] = Lev->Plan0.plan->N[1];
  nodes[2] = Lev->Plan0.plan->N[2];  
  neighbour[0] = Lev->Plan0.plan->N[1] * Lev->Plan0.plan->N[2];
  neighbour[1] = Lev->Plan0.plan->N[2];
  neighbour[2] = 1.;
  int k, l, nr, i_check = i;
  double h;

  //iS = n[2] + N[2]*(n[1] + N[1]*n0);  // howto access the array ...

  for (k=NDIM-1;k>=0;k--) { // for each axis
    h = 0.;
    for (l=0;l<NDIM;l++)
      h += Lat->RealBasis[k*NDIM+l]/Lev->Plan0.plan->N[k];  // finite distance
    h = sqrt(h);
    // check which limit exists: right, left, both?
    right_diff = 0.;
    left_diff = 0.;
    nr = 0; // neighbour counter
    if (i_check % nodes[k] != nodes[k]-1) {// right neighbour?
      right_diff = density[i] - density[i+neighbour[k]];
      nr++;
    }
    if (i_check % nodes[k] != 0) { // left neighbour?
      left_diff = density[i] - density[i-neighbour[k]];
      nr++;
    }
    res += (left_diff + right_diff)/(h*nr) * (left_diff + right_diff)/(h*nr);
    i_check /= nodes[k];  // remove axis from i_check
  }
  
  return sqrt(res);
}

/** SpinType-dependent calculation of exchange and correlation with free spin-polarisation energy with riemann tensor.
 * Discretely does the integration of \f${\cal E}_{xc} (n,\zeta) \cdot n\f$ on
 * the radial mesh of the respective real densities.<br>
 * 
 * Like CalculateXCEnergyNoRT(), only CalcSEXr() and CalcSECr() are
 * divided by RTFactor Lattice->RT.DensityR
 * \param *P Problem at hand
 */
void CalculateXCEnergyUseRT(struct Problem *P)
{
  struct Lattice *Lat = &P->Lat;
  struct Energy *E = Lat->E;
  struct PseudoPot *PP = &P->PP;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct Psis *Psi = &Lat->Psi;
  struct Density *Dens = Lev0->Dens;
  struct ExCor *EC = &P->ExCo;
  double SumEc = 0;
  double rs, p = 0.0, pUp = 0.0, pDown = 0.0, zeta, rsUp, rsDown, SumEx = 0.0;
  double Factor = R->XCEnergyFactor/Lev0->MaxN;
  int i;
  fftw_real *RTFactor = Lat->RT.DensityR[RTADetPreRT];
  
  for (i = 0; i < Dens->LocalSizeR; i++) {  // for each point of radial mesh take density p[i]
    p = Dens->DensityArray[TotalDensity][i];
    if (PP->corecorr == CoreCorrected)
      p += Dens->DensityArray[CoreWaveDensity][i];
    switch (Psi->PsiST) {
    case SpinDouble:
      pUp = 0.5*p; 
      pDown = 0.5*p;
      break;
    case SpinUp:  
    case SpinDown:
      pUp = Dens->DensityArray[TotalUpDensity][i];
      pDown = Dens->DensityArray[TotalDownDensity][i];
      if (PP->corecorr == CoreCorrected) {
        pUp += 0.5*Dens->DensityArray[CoreWaveDensity][i];
          pDown += 0.5*Dens->DensityArray[CoreWaveDensity][i];
      }
      break;
    default:
      ;
    }
    if ((p < 0) || (pUp < 0) || (pDown < 0)) {
      /*fprintf(stderr,"index %i pc %g p %g\n",i,Dens->DensityArray[CoreWaveDensity][i],p);*/
      p = 0.0; 
      pUp = 0.0;
      pDown = 0.0;
    }
    // .. calculate Ec and Ex and sum them up ...
    rs = Calcrs(EC,p);
    rsUp = Calcrs(EC,pUp);  
    rsDown = Calcrs(EC,pDown);
    SumEx += CalcSEXr(EC, rsUp, rsDown, pUp, pDown)/fabs(RTFactor[i]);
    zeta = CalcZeta(EC,pUp,pDown);
    SumEc += CalcSECr(EC, rs, zeta, p)/fabs(RTFactor[i]);
  }
  // ... and factorise with discrete integration width
  E->AllLocalDensityEnergy[CorrelationEnergy] = Factor*SumEc;
  E->AllLocalDensityEnergy[ExchangeEnergy] = Factor*SumEx;
}

/** Calculates SpinType-dependently exchange correlation potential with free spin-polarisation without Riemann tensor.
 * Goes through all possible nodes, calculates the potential and stores it in \a *HGR
 * \param *P Problem at hand
 * \param *HGR pointer storage array for (added!) result: \f$V_c (n,\zeta)(R) + V_x (n,\zeta)(R)\f$
 * \sa CalcSVXr(), CalcSVCr()
 */
void CalculateXCPotentialNoRT(struct Problem *P, fftw_real *HGR)
{
  struct Lattice *Lat = &P->Lat;
  struct PseudoPot *PP = &P->PP;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct LatticeLevel *LevS = R->LevS;
  struct Psis *Psi = &Lat->Psi;
  struct Density *Dens = Lev0->Dens;
  struct ExCor *EC = &P->ExCo;
  double rsX, rsC, zeta, p = 0.0, pUp = 0.0, pDown = 0.0;
  int nx,ny,nz,i;
  const int Nx = LevS->Plan0.plan->local_nx;
  const int Ny = LevS->Plan0.plan->N[1];
  const int Nz = LevS->Plan0.plan->N[2];
  const int NUpx = LevS->NUp[0];  // factors due to density being calculated on a finer grid
  const int NUpy = LevS->NUp[1];
  const int NUpz = LevS->NUp[2];
  double tmp;
  for (nx=0;nx<Nx;nx++)  
    for (ny=0;ny<Ny;ny++) 
      for (nz=0;nz<Nz;nz++) {
        i = nz*NUpz+Nz*NUpz*(ny*NUpy+Ny*NUpy*nx*NUpx);
        p = Dens->DensityArray[TotalDensity][i];
        //if (isnan(p)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): p_%i = NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
        if (PP->corecorr == CoreCorrected)
          p += Dens->DensityArray[CoreWaveDensity][i];
        switch (Psi->PsiST) {
        case SpinDouble:
          pUp = 0.5*p; 
          pDown = 0.5*p;
          break;
        case SpinUp:  
        case SpinDown:
          pUp = Dens->DensityArray[TotalUpDensity][i];
          pDown = Dens->DensityArray[TotalDownDensity][i];
          if (PP->corecorr == CoreCorrected) {   // additional factor due to PseudoPot'entials
            pUp += 0.5*Dens->DensityArray[CoreWaveDensity][i];
            pDown += 0.5*Dens->DensityArray[CoreWaveDensity][i];
          }
          break;
        default:
          ;
        }
        //if (isnan(pUp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): pUp_%i = NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
        //if (isnan(pDown)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): pDown_%i = NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
        if ((p < 0) || (pUp < 0) || (pDown < 0)) {
          p = 0; 
          pUp = 0;
          pDown = 0;
        }
        switch (Psi->PsiST) {
        case SpinUp:
          rsX = Calcrs(EC, pUp);
          break;
        case SpinDown:
          rsX = Calcrs(EC, pDown);
          break;
        case SpinDouble:
          //rsX = Calcrs(EC,p/2.);  // why half of it???
          rsX = Calcrs(EC, p);
          break;
        default:
          rsX = 0.0;
        }
        rsC = Calcrs(EC,p);
        zeta = CalcZeta(EC,pUp,pDown);
        //if (isnan(rsC)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): rsC%i = NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
        //if (isnan(zeta)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): zeta%i = NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
        switch (R->CurrentMin) {
          case UnOccupied:  // here epsilon appears instead of the potential in the integrand due to different variation
            tmp = CalcSEXr(EC,Calcrs(EC, pUp),Calcrs(EC, pDown),pUp,pDown);     // \todo ../p was here before test!
            //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i = NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
            tmp += CalcSECr(EC, rsX,zeta, 1);  
            //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i += NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
            break;
          default:
            tmp = CalcSVXr(EC,rsX,Psi->PsiST);
            //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i def= NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
            tmp += CalcSVCr(EC, rsC,zeta,Psi->PsiST);
            //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i def+= NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
            break;
        }
        //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i := NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
        HGR[i] += tmp;
        //if (isnan(HGR[i])) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): HGR[%i] = NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
      }
}

/** Calculates SpinType-dependently derivative of exchange correlation potential with free spin-polarisation without Riemann tensor.
 * \f[
 *    A^{XC} = \sum_R \frac{\delta V^{XC}}{\delta n} \rho^2 \qquad (\textnormal{section 5.1, line search})
 * \f]
 * With the density from Dens::DensityArray calculates discretely the integral over the summed
 * derivative, SpinType is taken from Psi::PsiST. Due to the principle of the theory the wave
 * functions themselves are not needed explicitely, see also CalculateXCEnergyNoRT()
 * \param *P Problem at hand
 * \param *PsiCD \f$\rho (r)\f$
 * \return \f$A^{XC}\f$
 * \sa CalcSVVXr(), CalcSVVCr() - derivatives of exchange and correlation potential
 */
double CalculateXCddEddt0NoRT(struct Problem *P, fftw_real *PsiCD)
{
  struct Lattice *Lat = &P->Lat;
  struct PseudoPot *PP = &P->PP;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct Psis *Psi = &Lat->Psi;
  struct Density *Dens = Lev0->Dens;
  struct ExCor *EC = &P->ExCo;
  double SumExc = 0.;
  double rsX=0.0, rsC, p = 0.0, pUp = 0.0, pDown = 0.0, zeta;
  double Factor = R->XCEnergyFactor/Lev0->MaxN;
  int i;
  
  for (i = 0; i < Dens->LocalSizeR; i++) {
    p = Dens->DensityArray[TotalDensity][i];
    if (PP->corecorr == CoreCorrected)
      p += Dens->DensityArray[CoreWaveDensity][i];
    switch (Psi->PsiST) {
      case SpinDouble:
        pUp = 0.5*p; 
        pDown = 0.5*p;
        break;
      case SpinUp:  
      case SpinDown:
        pUp = Dens->DensityArray[TotalUpDensity][i];
        pDown = Dens->DensityArray[TotalDownDensity][i];
        if (PP->corecorr == CoreCorrected) {
                pUp += 0.5*Dens->DensityArray[CoreWaveDensity][i];
                pDown += 0.5*Dens->DensityArray[CoreWaveDensity][i];
        }
      break;
    default:
      ;
    }
    if ((p < 0) || (pUp < 0) || (pDown < 0)) {
      /*fprintf(stderr,"index %i pc %g p %g\n",i,Dens->DensityArray[CoreWaveDensity][i],p);*/
      p = 0.0; 
      pUp = 0.0;
      pDown = 0.0;
    }
    switch (Psi->PsiST) {
    case SpinUp:
      rsX = Calcrs(EC,pUp);
      break;
    case SpinDown:
      rsX = Calcrs(EC,pDown);
      break;
    case SpinDouble:
      //rsX = Calcrs(EC,p/2.);  // why half of it???
      rsX = Calcrs(EC,p);
      break;
    }
    rsC = Calcrs(EC,p);
    zeta = CalcZeta(EC,pUp,pDown);
    SumExc += (CalcSVVXr(EC,rsX,Psi->PsiST) + CalcSVVCr(EC, rsC,zeta,Psi->PsiST,pUp,pDown))*PsiCD[i]*PsiCD[i];
  }
  return (SumExc*Factor);
}

/** Initialises ExCor structure with parametrization values.
 * All of the entries in the structure are set to parametrization values, see [CA80].
 * \param *P Problem at hand
 * \param *EC Exchange correlation ExCor structure
 */
void InitExchangeCorrelationEnergy(struct Problem *P, struct ExCor *EC) 
{
  EC->gamma[unpolarised] = -0.1423;
  EC->beta_1[unpolarised] = 1.0529;
  EC->beta_2[unpolarised] = 0.3334;
  EC->C[unpolarised] = 0.0020;
  EC->D[unpolarised] = -0.0116;
  EC->A[unpolarised] = 0.0311;
  EC->B[unpolarised] = -0.048;
  EC->gamma[polarised] = -0.0843;
  EC->beta_1[polarised] = 1.3981;
  EC->beta_2[polarised] = 0.2611;
  EC->C[polarised] = 0.0007;
  EC->D[polarised] = -0.0048;
  EC->A[polarised] = 0.01555;
  EC->B[polarised] = -0.0269;
  EC->fac34pi = -3./(4.*PI);
  EC->facexrs = EC->fac34pi*pow(9.*PI/4.,1./3.); 
  EC->epsilon0 = MYEPSILON;
  EC->fac6PI23 = pow(6./PI,2./3.);
  EC->facPI213 = pow(PI/2.,1./3.);
  EC->fac3PI23 = pow(3.*PI,2./3.);
  EC->fac6PIPI23 = pow(6.*PI*PI,2./3.);
  EC->fac243 = pow(2.,4./3.);
  EC->fac1213 = pow(1./2.,1./3.);
}