/** \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 $ */ #include #include #include #include #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:
* \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$
* \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$
* In the polarised case:
* \f$r_s\f$>=1: \f$-0.0843\cdot(1+1.3981\sqrt{r_s}+0.2611r_s)^{-1}\f$
* \f$r_s\f$<1: \f$-0.0269+0.01555\ln(r_s)-0.0048r_s+0.0007r_s\ln(r_s)\f$
*/ #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:
* \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$
* \f$r_s\f$=0: 0
* 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:
* \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$
* \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$
* \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:
* \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$
* \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$
* \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 = /*EC->fac1213**/VCR_unpol; 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$
* \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$
* \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.
* 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 rs, p = 0.0, pUp = 0.0, pDown = 0.0, zeta, rsUp, rsDown, SumEx=0.0; double Factor = R->XCEnergyFactor/Lev0->MaxN; 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]; if (R->CurrentMin > UnOccupied) 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: ; } // 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); } E->AllLocalDensityEnergy[CorrelationEnergy] = Factor*SumEc; E->AllLocalDensityEnergy[ExchangeEnergy] = Factor*SumEx; } /** 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.
* * 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;nxDensityArray[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) { // additional factor due to PseudoPot'entials pUp += 0.5*Dens->DensityArray[CoreWaveDensity][i]; pDown += 0.5*Dens->DensityArray[CoreWaveDensity][i]; } break; default: ; } 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); 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)/p; //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i un= NaN!\n", i); Error(SomeError, "NaN-Fehler!"); } tmp += CalcSECr(EC, rsX,zeta, 1); //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i un+= 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.); }