[a0bcf1] | 1 | /** \file excor.c
|
---|
| 2 | * Exchange correlation energy.
|
---|
| 3 | *
|
---|
| 4 | * The Exchange correlation energy is calculated via a local spin-density-approximation in
|
---|
| 5 | * a parametrized manner. It is the sum of the exchange and the correlation, per electron
|
---|
| 6 | * (thus, (discretely) integrated over the density).
|
---|
| 7 | * Due to the parametrization there is a huge number of functions evaluating small
|
---|
| 8 | * expressions or auxiliary variables with its first and second derivatives
|
---|
| 9 | * Calcf(), CalcDf(), CalcD2f() and CalcZeta(), CalcDZeta(), CalcD2Zeta(),
|
---|
| 10 | * which then help to evaluate the energy CalcECr(), CalcEXrUP(), in the summing
|
---|
| 11 | * CalcSECr(), CalcSEXr(), potential CalcVCr(), CalcVXrUP(), in its summation
|
---|
| 12 | * CalcSVCr(), CalcSVXr()
|
---|
| 13 | * and derivatives CalcVVCr(), CalcVVXrUP(), in its summation CalcSVVCr(), CalcSVVXr().
|
---|
| 14 | * All these are needed by the super functions evaluating exchange correlatiob
|
---|
| 15 | * potential CalculateXCPotentialNoRT(), exchange correlation energy without RiemannTensor
|
---|
| 16 | * CalculateXCEnergyNoRT() and with RiemannTensor CalculateXCEnergyUseRT(),
|
---|
| 17 | * and the second derivative CalculateXCddEddt0NoRT() (needed for conjugate gradient).
|
---|
| 18 | * During initilization the ExCor structure's entries is set to
|
---|
| 19 | * these fixed parameters InitExchangeCorrelationEnergy().
|
---|
| 20 | *
|
---|
| 21 | *
|
---|
| 22 | Project: ParallelCarParrinello
|
---|
| 23 | \author Jan Hamaekers
|
---|
| 24 | \date 2000
|
---|
| 25 |
|
---|
| 26 | File: excor.c
|
---|
| 27 | $Id: excor.c,v 1.47 2007-03-29 13:37:25 foo Exp $
|
---|
| 28 | */
|
---|
| 29 |
|
---|
| 30 | #include <stdlib.h>
|
---|
| 31 | #include <stdio.h>
|
---|
| 32 | #include <stddef.h>
|
---|
| 33 | #include <math.h>
|
---|
| 34 |
|
---|
| 35 | #include "data.h"
|
---|
| 36 | #include "mymath.h"
|
---|
| 37 | #include "run.h"
|
---|
| 38 | #include "myfft.h"
|
---|
| 39 | #include "errors.h"
|
---|
| 40 | #include "excor.h"
|
---|
| 41 | #include "helpers.h"
|
---|
| 42 |
|
---|
| 43 | /** Calculate Wigner-Seitz-Radius \f$r_s\f$
|
---|
| 44 | * \param *EC ExCor exchange correlation structure
|
---|
| 45 | * \param p electron density
|
---|
| 46 | * \return \f$(\frac{3}{4\pi \cdot p})^{1/3} \qquad (2.30)\f$
|
---|
| 47 | */
|
---|
| 48 | #ifdef HAVE_INLINE
|
---|
| 49 | inline double Calcrs(struct ExCor *EC, double p) {
|
---|
| 50 | #else
|
---|
| 51 | double Calcrs(struct ExCor *EC, double p) {
|
---|
| 52 | #endif
|
---|
| 53 | if (fabs(p) < EC->epsilon0) return(0);
|
---|
| 54 | //return (pow(3./(4.*PI*p),1./3.));
|
---|
| 55 | return (cbrt(3./(4.*PI*p)));
|
---|
| 56 | }
|
---|
| 57 |
|
---|
| 58 | /** Calculates exchange potential \f$V_{x}\f$.
|
---|
| 59 | * \f[
|
---|
| 60 | * V_{x} = {\cal E}_{x} + \Bigr ( \frac{d{\cal E}_{x}}{dn} \Bigl )_{n=n(r_s)} n(r_s) \qquad (2.32)
|
---|
| 61 | * \f]
|
---|
| 62 | * \param *EC ExCor exchange correlation structure
|
---|
| 63 | * \param rs Wigner-Seitz-Radius \f$r_s \f$, see Calcrs()
|
---|
| 64 | * \return \f$-\frac{1}{2} (\frac{6}{\pi})^{2/3} \cdot \frac{1}{rs}\f$
|
---|
| 65 | * \note The formula from CalcEXrUP() was used for the derivation.
|
---|
| 66 | */
|
---|
| 67 | #ifdef HAVE_INLINE
|
---|
| 68 | inline static double CalcVXrUP(struct ExCor *EC, double rs) {
|
---|
| 69 | #else
|
---|
| 70 | static double CalcVXrUP(struct ExCor *EC, double rs) {
|
---|
| 71 | #endif
|
---|
| 72 | if (fabs(rs) < EC->epsilon0) return(0);
|
---|
| 73 | return(-0.5*EC->fac6PI23/rs); // formula checked (16.5.06)
|
---|
| 74 | }
|
---|
| 75 |
|
---|
| 76 | /** Calculates first derivative of exchange potential \f$\frac{\delta V_{x}}{\delta n}\f$.
|
---|
| 77 | * \param *EC ExCor exchange correlation structure
|
---|
| 78 | * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
|
---|
| 79 | * \return \f$-\frac{2}{9\pi} (6\pi^2)^2/3 \cdot (rs)^2\f$
|
---|
| 80 | * \sa CalcVXrUP()
|
---|
| 81 | */
|
---|
| 82 | #ifdef HAVE_INLINE
|
---|
| 83 | inline static double CalcVVXrUP(struct ExCor *EC, double rs) {
|
---|
| 84 | #else
|
---|
| 85 | static double CalcVVXrUP(struct ExCor *EC, double rs) {
|
---|
| 86 | #endif
|
---|
| 87 | return (-2./(9*PI)*EC->fac6PIPI23*rs*rs); // formula checked (16.5.06)
|
---|
| 88 | }
|
---|
| 89 |
|
---|
| 90 | /** Calculates second derivative of exchange potential \f$\frac{\delta^2 V_{x}}{\delta n^2}\f$.
|
---|
| 91 | * \param *EC ExCor exchange correlation structure
|
---|
| 92 | * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
|
---|
| 93 | * \return \f$\frac{16}{81} (6\pi^2)^{2/3} \cdot (rs)^5\f$
|
---|
| 94 | * \sa CalcVXrUP(), CalcVVXrUP()
|
---|
| 95 | */
|
---|
| 96 | #ifdef HAVE_INLINE
|
---|
| 97 | inline static double CalcVVVXrUP(struct ExCor *EC, double rs) {
|
---|
| 98 | #else
|
---|
| 99 | static double CalcVVVXrUP(struct ExCor *EC, double rs) {
|
---|
| 100 | #endif
|
---|
| 101 | return (16./81.*EC->fac6PIPI23*rs*rs*rs*rs*rs); // formula checked (16.5.06)
|
---|
| 102 | }
|
---|
| 103 |
|
---|
| 104 | /** Calculates exchange energy over density for homogeneous charge within given Wigner-Seitz-Radius, \f${\cal E}_x\f$.
|
---|
| 105 | * The formula specified on the return statement is derived if into
|
---|
| 106 | * \f[ \forall 0 \geq \zeta \geq 1: \quad
|
---|
| 107 | * E_x = {\cal E}_x n = -\frac{3}{8} \Bigr (\frac{3(n_{up}+n_{down})}{\pi} \Bigl)^{1/3}
|
---|
| 108 | * \bigr[ (1+\zeta)^{4/3}+ (1-\zeta)^{4/3} \bigl ] \cdot (n_{up}+n_{down}) \qquad (2.52)
|
---|
| 109 | * \f]
|
---|
| 110 | * the definition of spin polarisation \f$\frac{n_{up} - n_{down}}{n_{up} + n_{down}}\f$
|
---|
| 111 | * is inserted and the expression evaluated.
|
---|
| 112 | * \param *EC ExCor exchange correlation structure
|
---|
| 113 | * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
|
---|
| 114 | * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
|
---|
| 115 | * \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$
|
---|
| 116 | */
|
---|
| 117 | #ifdef HAVE_INLINE
|
---|
| 118 | inline double CalcEXr(struct ExCor *EC, double rs, double zeta) {
|
---|
| 119 | #else
|
---|
| 120 | double CalcEXr(struct ExCor *EC, double rs, double zeta) {
|
---|
| 121 | #endif
|
---|
| 122 | if (fabs(rs) < EC->epsilon0) return(0);
|
---|
| 123 | //return (-3./8.*pow(3./(2.*PI),2./3.)/rs *(pow(1.+zeta,4./3.)+pow(1.-zeta,4./3.))); // formula checked
|
---|
| 124 | return (-3./8.*cbrt((3./(2.*PI))*(3./(2.*PI)))/rs *((1.+zeta)*cbrt(1.+zeta)+(1.-zeta)*cbrt(1.-zeta)));
|
---|
| 125 | }
|
---|
| 126 |
|
---|
| 127 | /** Calculates exchange energy over density for homogeneous charge within given Wigner-Seitz-Radius, \f${\cal E}_x\f$.
|
---|
| 128 | * The formula specified on the return statement is derived if into
|
---|
| 129 | * \f[ \forall 0 \geq \zeta \geq 1: \quad
|
---|
| 130 | * E_x = {\cal E}_x n = -\frac{3}{8} \Bigr (\frac{3(n_{up}+n_{down})}{\pi} \Bigl)^{1/3}
|
---|
| 131 | * \bigr[ (1+\zeta)^{4/3}+ (1-\zeta)^{4/3} \bigl ] \cdot (n_{up}+n_{down}) \qquad (2.52)
|
---|
| 132 | * \f]
|
---|
| 133 | * the definition of spin polarisation \f$\frac{n_{up} - n_{down}}{n_{up} + n_{down}}\f$
|
---|
| 134 | * is inserted and the expression evaluated.
|
---|
| 135 | * \param *EC ExCor exchange correlation structure
|
---|
| 136 | * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
|
---|
| 137 | * \return \f$-\frac{3}{8} (\frac{6}{\pi})^{2/3} \cdot \frac{1}{rs}\f$
|
---|
| 138 | */
|
---|
| 139 | #ifdef HAVE_INLINE
|
---|
| 140 | inline static double CalcEXrUP(struct ExCor *EC, double rs) {
|
---|
| 141 | #else
|
---|
| 142 | static double CalcEXrUP(struct ExCor *EC, double rs) {
|
---|
| 143 | #endif
|
---|
| 144 | if (fabs(rs) < EC->epsilon0) return(0);
|
---|
| 145 | return (-3./8.*EC->fac6PI23/rs); // formula checked
|
---|
| 146 | }
|
---|
| 147 |
|
---|
| 148 | /** Calculates correlation energy \f${\cal E}_c\f$ per electron for un-/polarised case.
|
---|
| 149 | * Formula derived from Monte-Carlo-simulations by Ceperley and Adler [CA80],
|
---|
| 150 | * parametrized by Perdew and Zunger [PZ81] for the polarised (\f$\zeta = 1\f$) or
|
---|
| 151 | * unpolarised (\f$\zeta = 0\f$) case.
|
---|
| 152 | * \param *EC ExCor exchange correlation structure
|
---|
| 153 | * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
|
---|
| 154 | * \param up UnPolarised: whether it's the polarised or unpolarised case
|
---|
| 155 | * \return In the unpolarised case: <br>
|
---|
| 156 | * \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>
|
---|
| 157 | * \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>
|
---|
| 158 | * In the polarised case: <br>
|
---|
| 159 | * \f$r_s\f$>=1: \f$-0.0843\cdot(1+1.3981\sqrt{r_s}+0.2611r_s)^{-1}\f$<br>
|
---|
| 160 | * \f$r_s\f$<1: \f$-0.0269+0.01555\ln(r_s)-0.0048r_s+0.0007r_s\ln(r_s)\f$<br>
|
---|
| 161 | */
|
---|
| 162 | #ifdef HAVE_INLINE
|
---|
| 163 | inline double CalcECr(struct ExCor *EC, double rs, enum UnPolarised up) {
|
---|
| 164 | #else
|
---|
| 165 | double CalcECr(struct ExCor *EC, double rs, enum UnPolarised up) {
|
---|
| 166 | #endif
|
---|
| 167 | double lrs;
|
---|
| 168 | double res=0;
|
---|
| 169 | if (rs >= 1) {
|
---|
| 170 | res= (EC->gamma[up]/(1.+EC->beta_1[up]*sqrt(rs)+EC->beta_2[up]*rs));
|
---|
| 171 | return (res);
|
---|
| 172 | }
|
---|
| 173 | if (rs <= EC->epsilon0) {
|
---|
| 174 | res = 0;
|
---|
| 175 | return(res);
|
---|
| 176 | }
|
---|
| 177 | lrs = log(rs);
|
---|
| 178 | res = (EC->A[up]*lrs+EC->B[up]+EC->C[up]*rs*lrs+EC->D[up]*rs);
|
---|
| 179 | return (res);
|
---|
| 180 | }
|
---|
| 181 |
|
---|
| 182 | /** Calculates correlation potential \f$V_{c}\f$.
|
---|
| 183 | * \f[
|
---|
| 184 | * V_{c} = {\cal E}_{c} + \Bigr ( \frac{d{\cal E}_{c}}{dn} \Bigl )_{n=n(r)} n(r) \qquad (2.32)
|
---|
| 185 | * \f]
|
---|
| 186 | * \param *EC ExCor exchange correlation structure
|
---|
| 187 | * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
|
---|
| 188 | * \param up UnPolarised
|
---|
| 189 | * \return in the un-/polarised case: <br>
|
---|
| 190 | * \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>
|
---|
| 191 | * \f$r_s\f$=0: 0 <br>
|
---|
| 192 | * 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$
|
---|
| 193 | */
|
---|
| 194 | #ifdef HAVE_INLINE
|
---|
| 195 | inline static double CalcVCr(struct ExCor *EC, double rs, enum UnPolarised up) {
|
---|
| 196 | #else
|
---|
| 197 | static double CalcVCr(struct ExCor *EC, double rs, enum UnPolarised up) {
|
---|
| 198 | #endif
|
---|
| 199 | double srs,lrs;
|
---|
| 200 | if (rs >= 1) {
|
---|
| 201 | srs = sqrt(rs);
|
---|
| 202 | 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)
|
---|
| 203 | }
|
---|
| 204 | if (rs <= EC->epsilon0) return (0);
|
---|
| 205 | lrs = log(rs);
|
---|
| 206 | 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)
|
---|
| 207 | }
|
---|
| 208 |
|
---|
| 209 | /** Calculates first derivative of correlation potential \f$\frac{\delta V_{c}}{\delta n}\f$.
|
---|
| 210 | * \f[
|
---|
| 211 | * \frac{\delta V_{c}}{\delta n} = \frac{\delta^2 {\cal E}_{c}}{\delta n^2}
|
---|
| 212 | * \f]
|
---|
| 213 | * \param *EC ExCor exchange correlation structure
|
---|
| 214 | * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
|
---|
| 215 | * \param up UnPolarised
|
---|
| 216 | * \return In the un-/polarised case: <br>
|
---|
| 217 | * \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
|
---|
| 218 | * + \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>
|
---|
| 219 | * \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>
|
---|
| 220 | * \sa CalcVCr()
|
---|
| 221 | */
|
---|
| 222 | #ifdef HAVE_INLINE
|
---|
| 223 | inline static double CalcVVCr(struct ExCor *EC, double rs, enum UnPolarised up) {
|
---|
| 224 | #else
|
---|
| 225 | static double CalcVVCr(struct ExCor *EC, double rs, enum UnPolarised up) {
|
---|
| 226 | #endif
|
---|
| 227 | double eps,deps;
|
---|
| 228 | if (rs <= EC->epsilon0) return (0);
|
---|
| 229 | if (rs >= 1) {
|
---|
| 230 | deps = 1.+EC->beta_1[up]*sqrt(rs)+EC->beta_2[up]*rs;
|
---|
| 231 | eps = EC->gamma[up]/deps;
|
---|
| 232 | //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)));
|
---|
| 233 | //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)
|
---|
| 234 | 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)
|
---|
| 235 | }
|
---|
| 236 | 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)
|
---|
| 237 | }
|
---|
| 238 |
|
---|
| 239 | /** Calculates second derivative of correlation potential \f$\frac{\delta^2 V_{c}}{\delta n^2}\f$.
|
---|
| 240 | * \f[
|
---|
| 241 | * \frac{\delta^2 V_{c}}{\delta n^2} = \frac{\delta^3 {\cal E}_{c}}{\delta n^3}
|
---|
| 242 | * \f]
|
---|
| 243 | * \param *EC ExCor exchange correlation structure
|
---|
| 244 | * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
|
---|
| 245 | * \param up UnPolarised
|
---|
| 246 | * \return In the un-/polarised case: <br>
|
---|
| 247 | * \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}
|
---|
| 248 | * \Bigr ( 35\beta_1[up]r_s^{13/2} + r_s^7(64\beta_2[up]+76\beta_1^2[up])
|
---|
| 249 | * + r_s^{15/2} (35\beta^3_1[up]+243\beta_1[up]\beta_2[up])
|
---|
| 250 | * + r_s^8 (176\beta^2_2[up]+140\beta^2_1[up]\beta_2[up])
|
---|
| 251 | * + r_s^{17/2} (175\beta_1[up]\beta^2_2[up])
|
---|
| 252 | * + r_s^9 (64\beta^3_2[up])
|
---|
| 253 | * \Bigl)\f$<br>
|
---|
| 254 | * \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>
|
---|
| 255 | * \sa CalcVCr()
|
---|
| 256 | */
|
---|
| 257 | #ifdef HAVE_INLINE
|
---|
| 258 | inline static double CalcVVVCr(struct ExCor *EC, double rs, enum UnPolarised up) {
|
---|
| 259 | #else
|
---|
| 260 | static double CalcVVVCr(struct ExCor *EC, double rs, enum UnPolarised up) {
|
---|
| 261 | #endif
|
---|
| 262 | double eps,deps;
|
---|
| 263 | double beta11,beta12,beta13,beta21,beta22,beta23;
|
---|
| 264 | if (rs <= EC->epsilon0) return (0);
|
---|
| 265 | if (rs >= 1) {
|
---|
| 266 | deps = 1.+EC->beta_1[up]*sqrt(rs)+EC->beta_2[up]*rs;
|
---|
| 267 | eps=EC->gamma[up]/deps;
|
---|
| 268 | beta11 = EC->beta_1[up];
|
---|
| 269 | beta12 = beta11*EC->beta_1[up];
|
---|
| 270 | beta13 = beta12*EC->beta_1[up];
|
---|
| 271 | beta21 = EC->beta_2[up];
|
---|
| 272 | beta22 = beta21*EC->beta_2[up];
|
---|
| 273 | beta23 = beta22*EC->beta_2[up];
|
---|
| 274 | // return (-2.*eps*PI*PI/(243.*deps*deps*deps)*(
|
---|
| 275 | // pow(rs,13./2.) * (35.*beta11)
|
---|
| 276 | // + pow(rs,7.) * (64.*beta21 + 76.*beta12)
|
---|
| 277 | // + pow(rs,15./2.) * (35.*beta13 + 234.*beta11*beta21)
|
---|
| 278 | // + pow(rs,8.) * (176.*beta22 + 140.*beta12*beta21)
|
---|
| 279 | // + pow(rs,17./2.) * (175.*beta11*beta22)
|
---|
| 280 | // + pow(rs,9.) * (64.*beta23) ));
|
---|
| 281 | return (-2.*eps*PI*PI/(243.*deps*deps*deps)*(
|
---|
| 282 | rs*rs*rs*rs*rs*rs*(
|
---|
| 283 | sqrt(rs) *((35.*beta11) // fractial ones
|
---|
| 284 | + rs *((35.*beta13 + 234.*beta11*beta21)
|
---|
| 285 | + rs * (175.*beta11*beta22))
|
---|
| 286 | )
|
---|
| 287 | + rs * ((64.*beta21 + 76.*beta12) // integer ones
|
---|
| 288 | + rs * ((176.*beta22 + 140.*beta12*beta21)
|
---|
| 289 | + rs * ((64.*beta23) ))))));
|
---|
| 290 | }
|
---|
| 291 | //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])));
|
---|
| 292 | 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])));
|
---|
| 293 | }
|
---|
| 294 |
|
---|
| 295 | /** Calculates spin polarisation \f$\zeta\f$
|
---|
| 296 | * \param *EC ExCor exchange correlation structure
|
---|
| 297 | * \param pUp SpinUp electron density
|
---|
| 298 | * \param pDown SpinDown electron density
|
---|
| 299 | * \return \f$\frac{pUp - pDown}{pUp + pDown}\f$
|
---|
| 300 | * \note zeta is never less than ExCor::epsilon0
|
---|
| 301 | */
|
---|
| 302 | #ifdef HAVE_INLINE
|
---|
| 303 | inline double CalcZeta(struct ExCor *EC, double pUp, double pDown) {
|
---|
| 304 | #else
|
---|
| 305 | double CalcZeta(struct ExCor *EC, double pUp, double pDown) {
|
---|
| 306 | #endif
|
---|
| 307 | if (fabs(pUp+pDown) < EC->epsilon0) return(0);
|
---|
| 308 | return ((pUp - pDown)/(pUp + pDown));
|
---|
| 309 | }
|
---|
| 310 |
|
---|
| 311 | /** Calculates derivative \f$\frac{\delta \zeta}{\delta p}\f$
|
---|
| 312 | * \param *EC ExCor exchange correlation structure
|
---|
| 313 | * \param ST SpinType
|
---|
| 314 | * \param pUp SpinUp density
|
---|
| 315 | * \param pDown SpinDown density
|
---|
| 316 | * \return p is pUp or pDown: \f$\pm 2 \frac{p}{ (pUp+pDown)^2 }\f$
|
---|
| 317 | * \sa CalcZeta()
|
---|
| 318 | */
|
---|
| 319 | #ifdef HAVE_INLINE
|
---|
| 320 | inline static double CalcDzeta(struct ExCor *EC, enum SpinType ST, double pUp, double pDown) {
|
---|
| 321 | #else
|
---|
| 322 | static double CalcDzeta(struct ExCor *EC, enum SpinType ST, double pUp, double pDown) {
|
---|
| 323 | #endif
|
---|
| 324 | double res = 0;
|
---|
| 325 | if (fabs(pUp+pDown) < EC->epsilon0) return(0);
|
---|
| 326 | //EC = EC;
|
---|
| 327 | switch(ST) {
|
---|
| 328 | case SpinUp:
|
---|
| 329 | res = 2.*pDown/(pUp+pDown)/(pUp+pDown); // formula checked
|
---|
| 330 | break;
|
---|
| 331 | case SpinDown:
|
---|
| 332 | res = -2.*pUp/(pUp+pDown)/(pUp+pDown); // formula checked
|
---|
| 333 | break;
|
---|
| 334 | case SpinDouble:
|
---|
| 335 | res = 0;
|
---|
| 336 | break;
|
---|
| 337 | }
|
---|
| 338 | return (res);
|
---|
| 339 | }
|
---|
| 340 |
|
---|
| 341 | /** Calculates second derivative \f$\frac{\delta^2 \zeta}{\delta p^2}\f$
|
---|
| 342 | * \param *EC ExCor exchange correlation structure
|
---|
| 343 | * \param ST SpinType
|
---|
| 344 | * \param pUp SpinUp density
|
---|
| 345 | * \param pDown SpinDown density
|
---|
| 346 | * \return p is pUp or pDown: \f$\pm 4 \frac{p}{ (pUp+pDown)^3 }\f$
|
---|
| 347 | * \sa CalcZeta(), CalcDZeta()
|
---|
| 348 | */
|
---|
| 349 | #ifdef HAVE_INLINE
|
---|
| 350 | inline static double CalcD2zeta(struct ExCor *EC, enum SpinType ST, double pUp, double pDown) {
|
---|
| 351 | #else
|
---|
| 352 | static double CalcD2zeta(struct ExCor *EC, enum SpinType ST, double pUp, double pDown) {
|
---|
| 353 | #endif
|
---|
| 354 | double res = 0;
|
---|
| 355 | if (fabs(pUp+pDown) < EC->epsilon0) return(0);
|
---|
| 356 | //EC = EC;
|
---|
| 357 | switch(ST) {
|
---|
| 358 | case SpinUp:
|
---|
| 359 | res = -4.*pDown/(pUp+pDown)/(pUp+pDown)/(pUp+pDown);
|
---|
| 360 | break;
|
---|
| 361 | case SpinDown:
|
---|
| 362 | res = 4.*pUp/(pUp+pDown)/(pUp+pDown)/(pUp+pDown);
|
---|
| 363 | break;
|
---|
| 364 | case SpinDouble:
|
---|
| 365 | res = 0;
|
---|
| 366 | break;
|
---|
| 367 | }
|
---|
| 368 | return (res);
|
---|
| 369 | }
|
---|
| 370 |
|
---|
| 371 | /** Calculates third derivative \f$\frac{\delta^3 \zeta}{\delta p^3}\f$
|
---|
| 372 | * \param *EC ExCor exchange correlation structure
|
---|
| 373 | * \param ST SpinType
|
---|
| 374 | * \param pUp SpinUp density
|
---|
| 375 | * \param pDown SpinDown density
|
---|
| 376 | * \return p is pUp or pDown: \f$\pm -12 \frac{p}{ (pUp+pDown)^4 }\f$
|
---|
| 377 | * \sa CalcZeta(), CalcDZeta(), CalcD2Zeta()
|
---|
| 378 | */
|
---|
| 379 | #ifdef HAVE_INLINE
|
---|
| 380 | inline static double CalcD3zeta(struct ExCor *EC, enum SpinType ST, double pUp, double pDown) {
|
---|
| 381 | #else
|
---|
| 382 | static double CalcD3zeta(struct ExCor *EC, enum SpinType ST, double pUp, double pDown) {
|
---|
| 383 | #endif
|
---|
| 384 | double res = 0;
|
---|
| 385 | if (fabs(pUp+pDown) < EC->epsilon0) return(0);
|
---|
| 386 | //EC = EC;
|
---|
| 387 | switch(ST) {
|
---|
| 388 | case SpinUp:
|
---|
| 389 | res = +12.*pDown/(pUp+pDown)/(pUp+pDown)/(pUp+pDown)/(pUp+pDown);
|
---|
| 390 | break;
|
---|
| 391 | case SpinDown:
|
---|
| 392 | res = -12.*pUp/(pUp+pDown)/(pUp+pDown)/(pUp+pDown)/(pUp+pDown);
|
---|
| 393 | break;
|
---|
| 394 | case SpinDouble:
|
---|
| 395 | res = 0;
|
---|
| 396 | break;
|
---|
| 397 | }
|
---|
| 398 | return (res);
|
---|
| 399 | }
|
---|
| 400 |
|
---|
| 401 | /** Calculates auxiliary factor \f$f(\zeta)\f$.
|
---|
| 402 | * This auxiliary factor is needed in the parametrized evaluation of the
|
---|
| 403 | * full spin-polarised correlation energy \f${\cal E}_c\f$ (see section 2.3.6)
|
---|
| 404 | * \param *EC ExCor exchange correlation structure
|
---|
| 405 | * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
|
---|
| 406 | * \return \f$\frac{(1+\zeta)^{4/3} + (1-\zeta)^{4/3} - 2} {2^{4/3}-2}\f$
|
---|
| 407 | */
|
---|
| 408 | #ifdef HAVE_INLINE
|
---|
| 409 | inline static double Calcf(struct ExCor *EC, double zeta) {
|
---|
| 410 | #else
|
---|
| 411 | static double Calcf(struct ExCor *EC, double zeta) {
|
---|
| 412 | #endif
|
---|
| 413 | double res =0;
|
---|
| 414 | EC = EC;
|
---|
| 415 | if (fabs(zeta) < EC->epsilon0) return(0.); // unpolarised case
|
---|
| 416 | //res = ((pow(1.+zeta,4./3.)+pow(1.-zeta,4./3.)-2.)/(EC->fac243-2.));
|
---|
| 417 | res = ((cbrt(1.+zeta)*(1.+zeta)+cbrt(1.-zeta)*(1.-zeta)-2.)/(EC->fac243-2.));
|
---|
| 418 | return (res);
|
---|
| 419 | }
|
---|
| 420 |
|
---|
| 421 | /** Calculates the derivative \f$\frac{\delta f}{\delta \zeta}\f$.
|
---|
| 422 | * \param *EC ExCor exchange correlation structure
|
---|
| 423 | * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
|
---|
| 424 | * \return \f$\frac{4}{3} \frac{(\zeta+1)^{1/3} - (1-\zeta)^{1/3}} {2^{4/3}-2}\f$
|
---|
| 425 | * \sa Calcf()
|
---|
| 426 | */
|
---|
| 427 | #ifdef HAVE_INLINE
|
---|
| 428 | inline static double CalcDf(struct ExCor *EC, double zeta) {
|
---|
| 429 | #else
|
---|
| 430 | static double CalcDf(struct ExCor *EC, double zeta) {
|
---|
| 431 | #endif
|
---|
| 432 | if (fabs(zeta) < EC->epsilon0) return(0.); // unpolarised case
|
---|
| 433 | //return((4./3.)*(pow(zeta+1.,1./3.)-pow(1.-zeta,1./3.))/(EC->fac243-2.)); // formula checked (16.5.06)
|
---|
| 434 | return((4./3.)*(cbrt(zeta+1.)-cbrt(1.-zeta))/(EC->fac243-2.)); // formula checked (16.5.06)
|
---|
| 435 | }
|
---|
| 436 |
|
---|
| 437 | /** Calculates second derivative \f$\frac{\delta^2 f}{\delta \zeta^2}\f$.
|
---|
| 438 | * \param *EC ExCor exchange correlation structure
|
---|
| 439 | * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
|
---|
| 440 | * \return \f$\frac{4}{9} \frac{(\zeta+1)^{-2/3}-(1-\zeta)^{-2/3} }{2^{4/3}-2}\f$
|
---|
| 441 | * \sa Calcf(), CalcDf()
|
---|
| 442 | */
|
---|
| 443 | #ifdef HAVE_INLINE
|
---|
| 444 | inline static double CalcD2f(struct ExCor *EC, double zeta) {
|
---|
| 445 | #else
|
---|
| 446 | static double CalcD2f(struct ExCor *EC, double zeta) {
|
---|
| 447 | #endif
|
---|
| 448 | if (fabs(zeta) < EC->epsilon0) return(0.); // unpolarised case
|
---|
| 449 | if (1.-fabs(zeta) < EC->epsilon0) return (0.); // second derivative not defined for these cases (which are needed!)
|
---|
| 450 | //if (1-fabs(zeta) < EC->epsilon0) { fprintf(stderr,"CalcD2f: zeta = %lg\n",zeta); return(0); }
|
---|
| 451 | //return((4./9.)*(pow(zeta+1.,-2./3.)+pow(1.-zeta,-2./3.))/(EC->fac243-2.)); // formula checked (16.5.06)
|
---|
[f2ce71c] | 452 | return((4./9.)*(1./cbrt((zeta+1.)*(zeta+1.))+1./cbrt((1.-zeta)*(1.-zeta)))/(EC->fac243-2.)); // formula checked (16.5.06)
|
---|
[a0bcf1] | 453 | }
|
---|
| 454 |
|
---|
| 455 | /** Calculates third derivative \f$\frac{\delta^3 f}{\delta \zeta^3}\f$.
|
---|
| 456 | * \param *EC ExCor exchange correlation structure
|
---|
| 457 | * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
|
---|
| 458 | * \return \f$\frac{8}{27} \frac{-(\zeta+1)^{-5/3}+(1-\zeta)^{-5/3} }{2^{4/3}-2}\f$
|
---|
| 459 | * \sa Calcf(), CalcDf()
|
---|
| 460 | */
|
---|
| 461 | #ifdef HAVE_INLINE
|
---|
| 462 | inline static double CalcD3f(struct ExCor *EC, double zeta) {
|
---|
| 463 | #else
|
---|
| 464 | static double CalcD3f(struct ExCor *EC, double zeta) {
|
---|
| 465 | #endif
|
---|
| 466 | if (fabs(zeta) < EC->epsilon0) return(0.); // unpolarised case
|
---|
| 467 | //if (1-fabs(zeta) < EC->epsilon0) return (0.); // second derivative not defined for these cases (which are needed!)
|
---|
| 468 | if (1-fabs(zeta) < EC->epsilon0) { fprintf(stderr,"CalcD2f: zeta = %lg\n",zeta); return(0.); }
|
---|
| 469 | //return((-8./27.)*(pow(zeta+1.,-5./3.)-pow(1.-zeta,-5./3.))/(EC->fac243-2.)); // formula checked (16.5.06)
|
---|
| 470 | 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)
|
---|
| 471 | }
|
---|
| 472 |
|
---|
| 473 | /** Calculates correlation energy with free spin-polarisation \f$E_c (n,\zeta)\f$.
|
---|
| 474 | * \param *EC exchange correlation structure
|
---|
| 475 | * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
|
---|
| 476 | * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
|
---|
| 477 | * \param p electron density \f$n\f$
|
---|
| 478 | * \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$
|
---|
| 479 | * \sa CalcECr()
|
---|
| 480 | */
|
---|
| 481 | #ifdef HAVE_INLINE
|
---|
| 482 | inline double CalcSECr(struct ExCor *EC, double rs, double zeta, double p) {
|
---|
| 483 | #else
|
---|
| 484 | double CalcSECr(struct ExCor *EC, double rs, double zeta, double p) {
|
---|
| 485 | #endif
|
---|
| 486 | double res =0;
|
---|
| 487 | double unpol, pol;
|
---|
| 488 | unpol = CalcECr(EC,rs,unpolarised);
|
---|
| 489 | pol = CalcECr(EC,rs,polarised);
|
---|
| 490 | res = (unpol+Calcf(EC,zeta)*(pol-unpol))*p;
|
---|
| 491 | return (res);
|
---|
| 492 | }
|
---|
| 493 |
|
---|
| 494 | /** Calculates SpinType-dependently correlation potential with free spin-polarisation \f$V_c (n,\zeta)\f$.
|
---|
| 495 | * \f[
|
---|
| 496 | * 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)
|
---|
| 497 | * + [{\cal E}_c(n,1)-{\cal E}_c(N,0)]\frac{\delta f(\zeta)}{\delta \zeta}\frac{\delta\zeta}{\delta n}
|
---|
| 498 | * \f]
|
---|
| 499 | * \param *EC ExCor exchange correlation structure
|
---|
| 500 | * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
|
---|
| 501 | * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
|
---|
| 502 | * \param ST SpinType
|
---|
| 503 | * \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$
|
---|
| 504 | */
|
---|
| 505 | #ifdef HAVE_INLINE
|
---|
| 506 | inline static double CalcSVCr(struct ExCor *EC, double rs, double zeta, enum SpinType ST) {
|
---|
| 507 | #else
|
---|
| 508 | static double CalcSVCr(struct ExCor *EC, double rs, double zeta, enum SpinType ST) {
|
---|
| 509 | #endif
|
---|
| 510 | double res = 0;
|
---|
| 511 | double VCR_unpol = CalcVCr(EC,rs,unpolarised);
|
---|
| 512 | switch (ST) {
|
---|
| 513 | case SpinUp:
|
---|
| 514 | 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)
|
---|
| 515 | break;
|
---|
| 516 | case SpinDown:
|
---|
| 517 | 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)
|
---|
| 518 | break;
|
---|
| 519 | case SpinDouble:
|
---|
[6c8ee7] | 520 | res = VCR_unpol; /*EC->fac1213**/
|
---|
[a0bcf1] | 521 | break;
|
---|
| 522 | }
|
---|
| 523 | return (res);
|
---|
| 524 | }
|
---|
| 525 |
|
---|
| 526 | /** Calculates complete exchange energy \f$E_x (n)\f$.
|
---|
| 527 | * \param *EC exchange correlation structure
|
---|
| 528 | * \param rsUp Wigner-Seitz-Radius \f$r_s\f$ for pUp
|
---|
| 529 | * \param rsDown Wigner-Seitz-Radius \f$r_s\f$ for pDown
|
---|
| 530 | * \param pUp SpinUp electron density
|
---|
| 531 | * \param pDown SpinDown electron density
|
---|
| 532 | * \return \f$E_x = {\cal E}_x \cdot n\f$
|
---|
| 533 | * \sa CalcEXrUP(), CalculateXCEnergyNoRT()
|
---|
| 534 | */
|
---|
| 535 | #ifdef HAVE_INLINE
|
---|
| 536 | inline double CalcSEXr(struct ExCor *EC, double rsUp, double rsDown, double pUp, double pDown) {
|
---|
| 537 | #else
|
---|
| 538 | double CalcSEXr(struct ExCor *EC, double rsUp, double rsDown, double pUp, double pDown) {
|
---|
| 539 | #endif
|
---|
| 540 | return(CalcEXrUP(EC,rsUp)*pUp+CalcEXrUP(EC,rsDown)*pDown);
|
---|
| 541 | }
|
---|
| 542 |
|
---|
| 543 | /** Calculates SpinType-dependently exchange potential \f$V_x (n)\f$.
|
---|
| 544 | * Essentially, just CalcVXrUP() is called, even in SpinDouble case
|
---|
| 545 | * \param *EC ExCor exchange correlation structure
|
---|
| 546 | * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
|
---|
| 547 | * \param ST SpinType
|
---|
| 548 | * \return CalcVXrUP()
|
---|
| 549 | */
|
---|
| 550 | #ifdef HAVE_INLINE
|
---|
| 551 | inline static double CalcSVXr(struct ExCor *EC, double rs, enum SpinType ST) {
|
---|
| 552 | #else
|
---|
| 553 | static double CalcSVXr(struct ExCor *EC, double rs, enum SpinType ST) {
|
---|
| 554 | #endif
|
---|
| 555 | double res = 0;
|
---|
| 556 | switch (ST) {
|
---|
| 557 | case SpinUp:
|
---|
| 558 | case SpinDown:
|
---|
| 559 | res = CalcVXrUP(EC,rs);
|
---|
| 560 | break;
|
---|
| 561 | case SpinDouble:
|
---|
| 562 | res = EC->fac1213*CalcVXrUP(EC,rs);
|
---|
| 563 | break;
|
---|
| 564 | }
|
---|
| 565 | return (res);
|
---|
| 566 | }
|
---|
| 567 |
|
---|
| 568 | /** Calculates SpinType-dependently derivative of exchange potential \f$\frac{\delta V_x}{\delta n}\f$.
|
---|
| 569 | * Essentially, just CalcVVXrUP() is called, even in SpinDouble case
|
---|
| 570 | * \param *EC ExCor exchange correlation structure
|
---|
| 571 | * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
|
---|
| 572 | * \param ST SpinType
|
---|
| 573 | * \return CalcVVXrUP()
|
---|
| 574 | */
|
---|
| 575 | #ifdef HAVE_INLINE
|
---|
| 576 | inline double CalcSVVXr(struct ExCor *EC, double rs, enum SpinType ST) {
|
---|
| 577 | #else
|
---|
| 578 | double CalcSVVXr(struct ExCor *EC, double rs, enum SpinType ST) {
|
---|
| 579 | #endif
|
---|
| 580 | double res = 0;
|
---|
| 581 | switch (ST) {
|
---|
| 582 | case SpinUp:
|
---|
| 583 | case SpinDown:
|
---|
| 584 | res = CalcVVXrUP(EC,rs);
|
---|
| 585 | break;
|
---|
| 586 | case SpinDouble:
|
---|
| 587 | res = EC->fac1213*CalcVVXrUP(EC,rs);
|
---|
| 588 | break;
|
---|
| 589 | }
|
---|
| 590 | return (res);
|
---|
| 591 | }
|
---|
| 592 |
|
---|
| 593 | /** Calculates SpinType-dependently second derivative of exchange potential \f$\frac{\delta^2 V_x}{\delta n^2}\f$.
|
---|
| 594 | * Essentially, just CalcVVVXrUP() is called, even in SpinDouble case
|
---|
| 595 | * \param *EC ExCor exchange correlation structure
|
---|
| 596 | * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
|
---|
| 597 | * \param ST SpinType
|
---|
| 598 | * \return CalcVVVXrUP()
|
---|
| 599 | */
|
---|
| 600 | #ifdef HAVE_INLINE
|
---|
| 601 | inline double CalcSVVVXr(struct ExCor *EC, double rs, enum SpinType ST) {
|
---|
| 602 | #else
|
---|
| 603 | double CalcSVVVXr(struct ExCor *EC, double rs, enum SpinType ST) {
|
---|
| 604 | #endif
|
---|
| 605 | double res = 0;
|
---|
| 606 | switch (ST) {
|
---|
| 607 | case SpinUp:
|
---|
| 608 | case SpinDown:
|
---|
| 609 | res = CalcVVVXrUP(EC,rs);
|
---|
| 610 | break;
|
---|
| 611 | case SpinDouble:
|
---|
| 612 | res = EC->fac1213*CalcVVVXrUP(EC,rs);
|
---|
| 613 | break;
|
---|
| 614 | }
|
---|
| 615 | return (res);
|
---|
| 616 | }
|
---|
| 617 |
|
---|
| 618 | /** Calculates SpinType-dependently first derivative of correlation potential with free spin-polarisation \f$\frac{\delta V_c}{\delta n} (n,\zeta)\f$.
|
---|
| 619 | * \param *EC ExCor exchange correlation structure
|
---|
| 620 | * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
|
---|
| 621 | * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
|
---|
| 622 | * \param ST SpinType
|
---|
| 623 | * \param pUp SpinUp density
|
---|
| 624 | * \param pDown SpinDown density
|
---|
| 625 | * \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 )
|
---|
| 626 | * + 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>
|
---|
| 627 | * \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
|
---|
| 628 | * \frac{\delta f}{\delta \zeta})\f$
|
---|
| 629 | * \sa CalcSVCr()
|
---|
| 630 | */
|
---|
| 631 | #ifdef HAVE_INLINE
|
---|
| 632 | inline double CalcSVVCr(struct ExCor *EC, double rs, double zeta, enum SpinType ST, double pUp, double pDown) {
|
---|
| 633 | #else
|
---|
| 634 | double CalcSVVCr(struct ExCor *EC, double rs, double zeta, enum SpinType ST, double pUp, double pDown) {
|
---|
| 635 | #endif
|
---|
| 636 | double res = 0;
|
---|
| 637 | double VVCR_pol, VVCR_unpol, VCR_pol, VCR_unpol, ECr_pol, ECr_unpol, f, Df, D2f, DZeta, D2Zeta;
|
---|
| 638 | switch (ST) {
|
---|
| 639 | case SpinUp:
|
---|
| 640 | case SpinDown:
|
---|
| 641 | // store some function calls for faster and easierly understandable operation
|
---|
| 642 | VVCR_unpol = CalcVVCr(EC,rs,unpolarised);
|
---|
| 643 | VVCR_pol = CalcVVCr(EC,rs,polarised);
|
---|
| 644 | ECr_pol = CalcECr(EC,rs,polarised);
|
---|
| 645 | ECr_unpol = CalcECr(EC,rs,unpolarised);
|
---|
| 646 | f = Calcf(EC,zeta);
|
---|
| 647 | Df = CalcDf(EC,zeta);
|
---|
| 648 | D2Zeta = CalcD2zeta(EC,ST,pUp,pDown);
|
---|
| 649 | if (CalcDzeta(EC,ST,pUp,pDown) != 0.0) {
|
---|
| 650 | VCR_unpol = CalcVCr(EC,rs,unpolarised);
|
---|
| 651 | VCR_pol = CalcVCr(EC,rs,polarised);
|
---|
| 652 | D2f = CalcD2f(EC,zeta);
|
---|
| 653 | DZeta = CalcDzeta(EC,ST,pUp,pDown);
|
---|
| 654 | 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)
|
---|
| 655 | } else {
|
---|
| 656 | res = VVCR_unpol + f*(VVCR_pol-VVCR_unpol) + (ECr_pol-ECr_unpol)*(pUp+pDown)*(D2Zeta*Df); //formula checked (15.5.06)
|
---|
| 657 | }
|
---|
| 658 | break;
|
---|
| 659 | case SpinDouble:
|
---|
| 660 | res = CalcVVCr(EC,rs,unpolarised);
|
---|
| 661 | break;
|
---|
| 662 | }
|
---|
| 663 | return (res);
|
---|
| 664 | }
|
---|
| 665 |
|
---|
| 666 | /** Calculates SpinType-dependently second derivative of correlation potential with free spin-polarisation \f$\frac{\delta V_c}{\delta n} (n,\zeta)\f$.
|
---|
| 667 | * \param *EC ExCor exchange correlation structure
|
---|
| 668 | * \param rs Wigner-Seitz-Radius \f$r_s\f$, see Calcrs()
|
---|
| 669 | * \param zeta spin polarisation \f$\zeta\f$, see CalcZeta()
|
---|
| 670 | * \param ST SpinType
|
---|
| 671 | * \param pUp SpinUp density
|
---|
| 672 | * \param pDown SpinDown density
|
---|
| 673 | * \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 )
|
---|
| 674 | * + 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>
|
---|
| 675 | * \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 ) +
|
---|
| 676 | * + ({\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$
|
---|
| 677 | * \sa CalcSVCr(), CalcSVVCr()
|
---|
| 678 | */
|
---|
| 679 | #ifdef HAVE_INLINE
|
---|
| 680 | inline double CalcSVVVCr(struct ExCor *EC, double rs, double zeta, enum SpinType ST, double pUp, double pDown) {
|
---|
| 681 | #else
|
---|
| 682 | double CalcSVVVCr(struct ExCor *EC, double rs, double zeta, enum SpinType ST, double pUp, double pDown) {
|
---|
| 683 | #endif
|
---|
| 684 | double res = 0;
|
---|
| 685 | double VVVCR_pol, VVVCR_unpol, VVCR_pol, VVCR_unpol, VCR_pol, VCR_unpol, ECr_pol, ECr_unpol, f, Df, D2f, D3f, DZeta, D2Zeta, D3Zeta;
|
---|
| 686 | switch (ST) {
|
---|
| 687 | case SpinUp:
|
---|
| 688 | case SpinDown:
|
---|
| 689 | // store some function calls for faster and easierly understandable operation
|
---|
| 690 | VVVCR_unpol = CalcVVVCr(EC,rs,unpolarised);
|
---|
| 691 | VVVCR_pol = CalcVVVCr(EC,rs,polarised);
|
---|
| 692 | VCR_unpol = CalcVCr(EC,rs,unpolarised);
|
---|
| 693 | VCR_pol = CalcVCr(EC,rs,polarised);
|
---|
| 694 | ECr_pol = CalcECr(EC,rs,polarised);
|
---|
| 695 | ECr_unpol = CalcECr(EC,rs,unpolarised);
|
---|
| 696 | f = Calcf(EC,zeta);
|
---|
| 697 | Df = CalcDf(EC,zeta);
|
---|
| 698 | D2Zeta = CalcD2zeta(EC,ST,pUp,pDown);
|
---|
| 699 | D3Zeta = CalcD3zeta(EC,ST,pUp,pDown);
|
---|
| 700 | if (CalcDzeta(EC,ST,pUp,pDown) != 0.0) {
|
---|
| 701 | VVCR_unpol = CalcVVCr(EC,rs,unpolarised);
|
---|
| 702 | VVCR_pol = CalcVVCr(EC,rs,polarised);
|
---|
| 703 | D2f = CalcD2f(EC,zeta);
|
---|
| 704 | D3f = CalcD3f(EC,zeta);
|
---|
| 705 | DZeta = CalcDzeta(EC,ST,pUp,pDown);
|
---|
| 706 | 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)
|
---|
| 707 | } else {
|
---|
| 708 | 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)
|
---|
| 709 | }
|
---|
| 710 | break;
|
---|
| 711 | case SpinDouble:
|
---|
| 712 | res = CalcVVVCr(EC,rs,unpolarised);
|
---|
| 713 | break;
|
---|
| 714 | }
|
---|
| 715 | return (res);
|
---|
| 716 | }
|
---|
| 717 |
|
---|
| 718 | /** SpinType-dependent calculation of exchange and correlation energy with free spin-polarisation without riemann tensor.
|
---|
| 719 | * Discretely does the integration of \f${\cal E}_{xc} (n,\zeta) \cdot n\f$ on
|
---|
| 720 | * the radial mesh of the respective real densities.<br>
|
---|
| 721 | * Takes either whole density or separated for each SpinType and calls
|
---|
| 722 | * Calcrs(), then summing each CalcSEXr() for each Density::LocalSizeR
|
---|
| 723 | * and times factor RunStruct::XCEnergyFactor / LatticeLevel::MaxN
|
---|
| 724 | * \param *P Problem at hand
|
---|
| 725 | */
|
---|
| 726 | void CalculateXCEnergyNoRT(struct Problem *P)
|
---|
| 727 | {
|
---|
| 728 | struct Lattice *Lat = &P->Lat;
|
---|
| 729 | struct Energy *E = Lat->E;
|
---|
| 730 | struct PseudoPot *PP = &P->PP;
|
---|
| 731 | struct RunStruct *R = &P->R;
|
---|
| 732 | struct LatticeLevel *Lev0 = R->Lev0;
|
---|
| 733 | struct Psis *Psi = &Lat->Psi;
|
---|
| 734 | struct Density *Dens = Lev0->Dens;
|
---|
| 735 | struct ExCor *EC = &P->ExCo;
|
---|
| 736 | double SumEc = 0;
|
---|
[6c8ee7] | 737 | //double SumEx_GC = 0.; // Gradient correction part according to Becke '92
|
---|
[a0bcf1] | 738 | double rs, p = 0.0, pUp = 0.0, pDown = 0.0, zeta, rsUp, rsDown, SumEx=0.0;
|
---|
| 739 | double Factor = R->XCEnergyFactor/Lev0->MaxN;
|
---|
[6c8ee7] | 740 | //double Dp, DpUp, DpDown;
|
---|
[a0bcf1] | 741 | int i;
|
---|
| 742 |
|
---|
| 743 | for (i = 0; i < Dens->LocalSizeR; i++) { // for each node in radial mesh
|
---|
| 744 | // put (corecorrected) densities in p, pUp, pDown
|
---|
| 745 | p = Dens->DensityArray[TotalDensity][i];
|
---|
[6c8ee7] | 746 | //Dp = DensityGradient(Dens->DensityArray[TotalDensity], i, Lev0, Lat);
|
---|
[a0bcf1] | 747 | if (R->CurrentMin > UnOccupied)
|
---|
| 748 | if (PP->corecorr == CoreCorrected)
|
---|
| 749 | p += Dens->DensityArray[CoreWaveDensity][i];
|
---|
| 750 | switch (Psi->PsiST) {
|
---|
[6c8ee7] | 751 | default:
|
---|
[a0bcf1] | 752 | case SpinDouble:
|
---|
| 753 | pUp = 0.5*p;
|
---|
| 754 | pDown = 0.5*p;
|
---|
[6c8ee7] | 755 | //DpUp = 0.5*Dp;
|
---|
| 756 | //DpDown = 0.5*Dp;
|
---|
[a0bcf1] | 757 | break;
|
---|
| 758 | case SpinUp:
|
---|
| 759 | case SpinDown:
|
---|
| 760 | pUp = Dens->DensityArray[TotalUpDensity][i];
|
---|
| 761 | pDown = Dens->DensityArray[TotalDownDensity][i];
|
---|
[6c8ee7] | 762 | //DpUp = DensityGradient(Dens->DensityArray[TotalUpDensity], i, Lev0, Lat);
|
---|
| 763 | //DpDown = DensityGradient(Dens->DensityArray[TotalDownDensity], i, Lev0, Lat);
|
---|
[a0bcf1] | 764 | if (PP->corecorr == CoreCorrected) {
|
---|
| 765 | pUp += 0.5*Dens->DensityArray[CoreWaveDensity][i];
|
---|
| 766 | pDown += 0.5*Dens->DensityArray[CoreWaveDensity][i];
|
---|
| 767 | }
|
---|
| 768 | break;
|
---|
| 769 | }
|
---|
| 770 | // set all to zero if one of them is negative
|
---|
| 771 | if ((p < 0) || (pUp < 0) || (pDown < 0)) {
|
---|
| 772 | /*fprintf(stderr,"index %i pc %g p %g\n",i,Dens->DensityArray[CoreWaveDensity][i],p);*/
|
---|
| 773 | p = 0.0;
|
---|
| 774 | pUp = 0.0;
|
---|
| 775 | pDown = 0.0;
|
---|
| 776 | }
|
---|
| 777 | // Calculation with the densities and summation
|
---|
| 778 | rs = Calcrs(EC,p);
|
---|
| 779 | rsUp = Calcrs(EC,pUp);
|
---|
| 780 | rsDown = Calcrs(EC,pDown);
|
---|
| 781 | SumEx += CalcSEXr(EC, rsUp, rsDown, pUp, pDown);
|
---|
| 782 | zeta = CalcZeta(EC,pUp,pDown);
|
---|
| 783 | SumEc += CalcSECr(EC, rs, zeta, p);
|
---|
[6c8ee7] | 784 | //SumEx_GC += CalcSE_GC(EC, pUp, DpUp);
|
---|
| 785 | //SumEx_GC += CalcSE_GC(EC, pDown, DpDown);
|
---|
[a0bcf1] | 786 | }
|
---|
| 787 | E->AllLocalDensityEnergy[CorrelationEnergy] = Factor*SumEc;
|
---|
[6c8ee7] | 788 | E->AllLocalDensityEnergy[ExchangeEnergy] = Factor*(SumEx); // - SumEx_GC);
|
---|
[a0bcf1] | 789 | }
|
---|
| 790 |
|
---|
[b70503] | 791 | /** Returns \f$y = \sinh^{-1}(x)\f$.
|
---|
| 792 | * Recall that \f$ x = sinh (y) = \frac{\exp(y)-\exp(-y)}{2}\f$. Then, solve the in \f$\exp(y)\f$ quadratic
|
---|
| 793 | * equation: \f$\exp(2y) - 1 - 2x\exp(y) = 0\f$ by pq-formula.
|
---|
| 794 | * \param arg argument \f$x\f$
|
---|
| 795 | * \return \f$\sinh^{-1}(x) = \ln(x+\sqrt{x^2+1})\f$
|
---|
| 796 | */
|
---|
| 797 | inline static double sinh_inverse(double arg)
|
---|
| 798 | {
|
---|
| 799 | return log(arg+sqrt(arg*arg+1));
|
---|
| 800 | }
|
---|
| 801 |
|
---|
| 802 | /** Gradient correction according to Becke '92.
|
---|
| 803 | * \param EC ExCor structure
|
---|
| 804 | * \param p \f$\rho\f$ local density
|
---|
| 805 | * \param Dp \f$|\nabla \rho|\f$ local magnitude of density gradient
|
---|
| 806 | * \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$
|
---|
| 807 | */
|
---|
| 808 | double CalcSE_GC(struct ExCor *EC, double p, double Dp)
|
---|
| 809 | {
|
---|
| 810 | double res = 0.;
|
---|
| 811 | double b = 0.0042; // in atomic units, empirical constant fitted to noble gas atoms
|
---|
| 812 | double x = Dp/pow(p, 4./3.);
|
---|
| 813 | res = b * pow(p, 4./3.) * x/(1+6.*b*x*sinh_inverse(x));
|
---|
| 814 |
|
---|
| 815 | return res;
|
---|
| 816 | }
|
---|
| 817 |
|
---|
| 818 | /** Evaluates magnitude of gradient by simple 3-star.
|
---|
| 819 | * \param *density density array
|
---|
| 820 | * \param i current node
|
---|
| 821 | * \param *Lev LatticeLevel structure for number of nodes per axis
|
---|
| 822 | * \param *Lat Lattice structure for axis lengths
|
---|
| 823 | * \return \f$|\nabla\rho|\f$
|
---|
| 824 | */
|
---|
| 825 | double DensityGradient(fftw_real *density, int i, struct LatticeLevel *Lev, struct Lattice *Lat)
|
---|
| 826 | {
|
---|
| 827 | double res=0., right_diff, left_diff;
|
---|
| 828 | int neighbour[NDIM], nodes[NDIM];
|
---|
| 829 | nodes[0] = Lev->Plan0.plan->local_nx;
|
---|
| 830 | nodes[1] = Lev->Plan0.plan->N[1];
|
---|
| 831 | nodes[2] = Lev->Plan0.plan->N[2];
|
---|
| 832 | neighbour[0] = Lev->Plan0.plan->N[1] * Lev->Plan0.plan->N[2];
|
---|
| 833 | neighbour[1] = Lev->Plan0.plan->N[2];
|
---|
| 834 | neighbour[2] = 1.;
|
---|
| 835 | int k, l, nr, i_check = i;
|
---|
| 836 | double h;
|
---|
| 837 |
|
---|
| 838 | //iS = n[2] + N[2]*(n[1] + N[1]*n0); // howto access the array ...
|
---|
| 839 |
|
---|
| 840 | for (k=NDIM-1;k>=0;k--) { // for each axis
|
---|
| 841 | h = 0.;
|
---|
| 842 | for (l=0;l<NDIM;l++)
|
---|
| 843 | h += Lat->RealBasis[k*NDIM+l]/Lev->Plan0.plan->N[k]; // finite distance
|
---|
| 844 | h = sqrt(h);
|
---|
| 845 | // check which limit exists: right, left, both?
|
---|
| 846 | right_diff = 0.;
|
---|
| 847 | left_diff = 0.;
|
---|
| 848 | nr = 0; // neighbour counter
|
---|
| 849 | if (i_check % nodes[k] != nodes[k]-1) {// right neighbour?
|
---|
| 850 | right_diff = density[i] - density[i+neighbour[k]];
|
---|
| 851 | nr++;
|
---|
| 852 | }
|
---|
| 853 | if (i_check % nodes[k] != 0) { // left neighbour?
|
---|
| 854 | left_diff = density[i] - density[i-neighbour[k]];
|
---|
| 855 | nr++;
|
---|
| 856 | }
|
---|
| 857 | res += (left_diff + right_diff)/(h*nr) * (left_diff + right_diff)/(h*nr);
|
---|
| 858 | i_check /= nodes[k]; // remove axis from i_check
|
---|
| 859 | }
|
---|
| 860 |
|
---|
| 861 | return sqrt(res);
|
---|
| 862 | }
|
---|
| 863 |
|
---|
[a0bcf1] | 864 | /** SpinType-dependent calculation of exchange and correlation with free spin-polarisation energy with riemann tensor.
|
---|
| 865 | * Discretely does the integration of \f${\cal E}_{xc} (n,\zeta) \cdot n\f$ on
|
---|
| 866 | * the radial mesh of the respective real densities.<br>
|
---|
| 867 | *
|
---|
| 868 | * Like CalculateXCEnergyNoRT(), only CalcSEXr() and CalcSECr() are
|
---|
| 869 | * divided by RTFactor Lattice->RT.DensityR
|
---|
| 870 | * \param *P Problem at hand
|
---|
| 871 | */
|
---|
| 872 | void CalculateXCEnergyUseRT(struct Problem *P)
|
---|
| 873 | {
|
---|
| 874 | struct Lattice *Lat = &P->Lat;
|
---|
| 875 | struct Energy *E = Lat->E;
|
---|
| 876 | struct PseudoPot *PP = &P->PP;
|
---|
| 877 | struct RunStruct *R = &P->R;
|
---|
| 878 | struct LatticeLevel *Lev0 = R->Lev0;
|
---|
| 879 | struct Psis *Psi = &Lat->Psi;
|
---|
| 880 | struct Density *Dens = Lev0->Dens;
|
---|
| 881 | struct ExCor *EC = &P->ExCo;
|
---|
| 882 | double SumEc = 0;
|
---|
| 883 | double rs, p = 0.0, pUp = 0.0, pDown = 0.0, zeta, rsUp, rsDown, SumEx = 0.0;
|
---|
| 884 | double Factor = R->XCEnergyFactor/Lev0->MaxN;
|
---|
| 885 | int i;
|
---|
| 886 | fftw_real *RTFactor = Lat->RT.DensityR[RTADetPreRT];
|
---|
| 887 |
|
---|
| 888 | for (i = 0; i < Dens->LocalSizeR; i++) { // for each point of radial mesh take density p[i]
|
---|
| 889 | p = Dens->DensityArray[TotalDensity][i];
|
---|
| 890 | if (PP->corecorr == CoreCorrected)
|
---|
| 891 | p += Dens->DensityArray[CoreWaveDensity][i];
|
---|
| 892 | switch (Psi->PsiST) {
|
---|
| 893 | case SpinDouble:
|
---|
| 894 | pUp = 0.5*p;
|
---|
| 895 | pDown = 0.5*p;
|
---|
| 896 | break;
|
---|
| 897 | case SpinUp:
|
---|
| 898 | case SpinDown:
|
---|
| 899 | pUp = Dens->DensityArray[TotalUpDensity][i];
|
---|
| 900 | pDown = Dens->DensityArray[TotalDownDensity][i];
|
---|
| 901 | if (PP->corecorr == CoreCorrected) {
|
---|
| 902 | pUp += 0.5*Dens->DensityArray[CoreWaveDensity][i];
|
---|
| 903 | pDown += 0.5*Dens->DensityArray[CoreWaveDensity][i];
|
---|
| 904 | }
|
---|
| 905 | break;
|
---|
| 906 | default:
|
---|
| 907 | ;
|
---|
| 908 | }
|
---|
| 909 | if ((p < 0) || (pUp < 0) || (pDown < 0)) {
|
---|
| 910 | /*fprintf(stderr,"index %i pc %g p %g\n",i,Dens->DensityArray[CoreWaveDensity][i],p);*/
|
---|
| 911 | p = 0.0;
|
---|
| 912 | pUp = 0.0;
|
---|
| 913 | pDown = 0.0;
|
---|
| 914 | }
|
---|
| 915 | // .. calculate Ec and Ex and sum them up ...
|
---|
| 916 | rs = Calcrs(EC,p);
|
---|
| 917 | rsUp = Calcrs(EC,pUp);
|
---|
| 918 | rsDown = Calcrs(EC,pDown);
|
---|
| 919 | SumEx += CalcSEXr(EC, rsUp, rsDown, pUp, pDown)/fabs(RTFactor[i]);
|
---|
| 920 | zeta = CalcZeta(EC,pUp,pDown);
|
---|
| 921 | SumEc += CalcSECr(EC, rs, zeta, p)/fabs(RTFactor[i]);
|
---|
| 922 | }
|
---|
| 923 | // ... and factorise with discrete integration width
|
---|
| 924 | E->AllLocalDensityEnergy[CorrelationEnergy] = Factor*SumEc;
|
---|
| 925 | E->AllLocalDensityEnergy[ExchangeEnergy] = Factor*SumEx;
|
---|
| 926 | }
|
---|
| 927 |
|
---|
| 928 | /** Calculates SpinType-dependently exchange correlation potential with free spin-polarisation without Riemann tensor.
|
---|
| 929 | * Goes through all possible nodes, calculates the potential and stores it in \a *HGR
|
---|
| 930 | * \param *P Problem at hand
|
---|
| 931 | * \param *HGR pointer storage array for (added!) result: \f$V_c (n,\zeta)(R) + V_x (n,\zeta)(R)\f$
|
---|
| 932 | * \sa CalcSVXr(), CalcSVCr()
|
---|
| 933 | */
|
---|
| 934 | void CalculateXCPotentialNoRT(struct Problem *P, fftw_real *HGR)
|
---|
| 935 | {
|
---|
| 936 | struct Lattice *Lat = &P->Lat;
|
---|
| 937 | struct PseudoPot *PP = &P->PP;
|
---|
| 938 | struct RunStruct *R = &P->R;
|
---|
| 939 | struct LatticeLevel *Lev0 = R->Lev0;
|
---|
| 940 | struct LatticeLevel *LevS = R->LevS;
|
---|
| 941 | struct Psis *Psi = &Lat->Psi;
|
---|
| 942 | struct Density *Dens = Lev0->Dens;
|
---|
| 943 | struct ExCor *EC = &P->ExCo;
|
---|
| 944 | double rsX, rsC, zeta, p = 0.0, pUp = 0.0, pDown = 0.0;
|
---|
| 945 | int nx,ny,nz,i;
|
---|
| 946 | const int Nx = LevS->Plan0.plan->local_nx;
|
---|
| 947 | const int Ny = LevS->Plan0.plan->N[1];
|
---|
| 948 | const int Nz = LevS->Plan0.plan->N[2];
|
---|
| 949 | const int NUpx = LevS->NUp[0]; // factors due to density being calculated on a finer grid
|
---|
| 950 | const int NUpy = LevS->NUp[1];
|
---|
| 951 | const int NUpz = LevS->NUp[2];
|
---|
| 952 | double tmp;
|
---|
| 953 | for (nx=0;nx<Nx;nx++)
|
---|
| 954 | for (ny=0;ny<Ny;ny++)
|
---|
| 955 | for (nz=0;nz<Nz;nz++) {
|
---|
| 956 | i = nz*NUpz+Nz*NUpz*(ny*NUpy+Ny*NUpy*nx*NUpx);
|
---|
| 957 | p = Dens->DensityArray[TotalDensity][i];
|
---|
[f7754e] | 958 | //if (isnan(p)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): p_%i = NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
|
---|
[a0bcf1] | 959 | if (PP->corecorr == CoreCorrected)
|
---|
| 960 | p += Dens->DensityArray[CoreWaveDensity][i];
|
---|
| 961 | switch (Psi->PsiST) {
|
---|
| 962 | case SpinDouble:
|
---|
| 963 | pUp = 0.5*p;
|
---|
| 964 | pDown = 0.5*p;
|
---|
| 965 | break;
|
---|
| 966 | case SpinUp:
|
---|
| 967 | case SpinDown:
|
---|
| 968 | pUp = Dens->DensityArray[TotalUpDensity][i];
|
---|
| 969 | pDown = Dens->DensityArray[TotalDownDensity][i];
|
---|
| 970 | if (PP->corecorr == CoreCorrected) { // additional factor due to PseudoPot'entials
|
---|
| 971 | pUp += 0.5*Dens->DensityArray[CoreWaveDensity][i];
|
---|
| 972 | pDown += 0.5*Dens->DensityArray[CoreWaveDensity][i];
|
---|
| 973 | }
|
---|
| 974 | break;
|
---|
| 975 | default:
|
---|
| 976 | ;
|
---|
| 977 | }
|
---|
[f7754e] | 978 | //if (isnan(pUp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): pUp_%i = NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
|
---|
| 979 | //if (isnan(pDown)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): pDown_%i = NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
|
---|
[a0bcf1] | 980 | if ((p < 0) || (pUp < 0) || (pDown < 0)) {
|
---|
| 981 | p = 0;
|
---|
| 982 | pUp = 0;
|
---|
| 983 | pDown = 0;
|
---|
| 984 | }
|
---|
| 985 | switch (Psi->PsiST) {
|
---|
| 986 | case SpinUp:
|
---|
| 987 | rsX = Calcrs(EC, pUp);
|
---|
| 988 | break;
|
---|
| 989 | case SpinDown:
|
---|
| 990 | rsX = Calcrs(EC, pDown);
|
---|
| 991 | break;
|
---|
| 992 | case SpinDouble:
|
---|
| 993 | //rsX = Calcrs(EC,p/2.); // why half of it???
|
---|
| 994 | rsX = Calcrs(EC, p);
|
---|
| 995 | break;
|
---|
| 996 | default:
|
---|
| 997 | rsX = 0.0;
|
---|
| 998 | }
|
---|
| 999 | rsC = Calcrs(EC,p);
|
---|
| 1000 | zeta = CalcZeta(EC,pUp,pDown);
|
---|
[f7754e] | 1001 | //if (isnan(rsC)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): rsC%i = NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
|
---|
| 1002 | //if (isnan(zeta)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): zeta%i = NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
|
---|
[a0bcf1] | 1003 | switch (R->CurrentMin) {
|
---|
| 1004 | case UnOccupied: // here epsilon appears instead of the potential in the integrand due to different variation
|
---|
[f7754e] | 1005 | tmp = CalcSEXr(EC,Calcrs(EC, pUp),Calcrs(EC, pDown),pUp,pDown); // \todo ../p was here before test!
|
---|
| 1006 | //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i = NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
|
---|
[a0bcf1] | 1007 | tmp += CalcSECr(EC, rsX,zeta, 1);
|
---|
[f7754e] | 1008 | //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i += NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
|
---|
[a0bcf1] | 1009 | break;
|
---|
| 1010 | default:
|
---|
| 1011 | tmp = CalcSVXr(EC,rsX,Psi->PsiST);
|
---|
| 1012 | //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i def= NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
|
---|
| 1013 | tmp += CalcSVCr(EC, rsC,zeta,Psi->PsiST);
|
---|
| 1014 | //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i def+= NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
|
---|
| 1015 | break;
|
---|
| 1016 | }
|
---|
| 1017 | //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i := NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
|
---|
| 1018 | HGR[i] += tmp;
|
---|
| 1019 | //if (isnan(HGR[i])) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): HGR[%i] = NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
|
---|
| 1020 | }
|
---|
| 1021 | }
|
---|
| 1022 |
|
---|
| 1023 | /** Calculates SpinType-dependently derivative of exchange correlation potential with free spin-polarisation without Riemann tensor.
|
---|
| 1024 | * \f[
|
---|
| 1025 | * A^{XC} = \sum_R \frac{\delta V^{XC}}{\delta n} \rho^2 \qquad (\textnormal{section 5.1, line search})
|
---|
| 1026 | * \f]
|
---|
| 1027 | * With the density from Dens::DensityArray calculates discretely the integral over the summed
|
---|
| 1028 | * derivative, SpinType is taken from Psi::PsiST. Due to the principle of the theory the wave
|
---|
| 1029 | * functions themselves are not needed explicitely, see also CalculateXCEnergyNoRT()
|
---|
| 1030 | * \param *P Problem at hand
|
---|
| 1031 | * \param *PsiCD \f$\rho (r)\f$
|
---|
| 1032 | * \return \f$A^{XC}\f$
|
---|
| 1033 | * \sa CalcSVVXr(), CalcSVVCr() - derivatives of exchange and correlation potential
|
---|
| 1034 | */
|
---|
| 1035 | double CalculateXCddEddt0NoRT(struct Problem *P, fftw_real *PsiCD)
|
---|
| 1036 | {
|
---|
| 1037 | struct Lattice *Lat = &P->Lat;
|
---|
| 1038 | struct PseudoPot *PP = &P->PP;
|
---|
| 1039 | struct RunStruct *R = &P->R;
|
---|
| 1040 | struct LatticeLevel *Lev0 = R->Lev0;
|
---|
| 1041 | struct Psis *Psi = &Lat->Psi;
|
---|
| 1042 | struct Density *Dens = Lev0->Dens;
|
---|
| 1043 | struct ExCor *EC = &P->ExCo;
|
---|
[f7754e] | 1044 | double SumExc = 0.;
|
---|
[a0bcf1] | 1045 | double rsX=0.0, rsC, p = 0.0, pUp = 0.0, pDown = 0.0, zeta;
|
---|
| 1046 | double Factor = R->XCEnergyFactor/Lev0->MaxN;
|
---|
| 1047 | int i;
|
---|
| 1048 |
|
---|
| 1049 | for (i = 0; i < Dens->LocalSizeR; i++) {
|
---|
| 1050 | p = Dens->DensityArray[TotalDensity][i];
|
---|
| 1051 | if (PP->corecorr == CoreCorrected)
|
---|
| 1052 | p += Dens->DensityArray[CoreWaveDensity][i];
|
---|
| 1053 | switch (Psi->PsiST) {
|
---|
| 1054 | case SpinDouble:
|
---|
| 1055 | pUp = 0.5*p;
|
---|
| 1056 | pDown = 0.5*p;
|
---|
| 1057 | break;
|
---|
| 1058 | case SpinUp:
|
---|
| 1059 | case SpinDown:
|
---|
| 1060 | pUp = Dens->DensityArray[TotalUpDensity][i];
|
---|
| 1061 | pDown = Dens->DensityArray[TotalDownDensity][i];
|
---|
| 1062 | if (PP->corecorr == CoreCorrected) {
|
---|
| 1063 | pUp += 0.5*Dens->DensityArray[CoreWaveDensity][i];
|
---|
| 1064 | pDown += 0.5*Dens->DensityArray[CoreWaveDensity][i];
|
---|
| 1065 | }
|
---|
| 1066 | break;
|
---|
| 1067 | default:
|
---|
| 1068 | ;
|
---|
| 1069 | }
|
---|
| 1070 | if ((p < 0) || (pUp < 0) || (pDown < 0)) {
|
---|
| 1071 | /*fprintf(stderr,"index %i pc %g p %g\n",i,Dens->DensityArray[CoreWaveDensity][i],p);*/
|
---|
| 1072 | p = 0.0;
|
---|
| 1073 | pUp = 0.0;
|
---|
| 1074 | pDown = 0.0;
|
---|
| 1075 | }
|
---|
| 1076 | switch (Psi->PsiST) {
|
---|
| 1077 | case SpinUp:
|
---|
| 1078 | rsX = Calcrs(EC,pUp);
|
---|
| 1079 | break;
|
---|
| 1080 | case SpinDown:
|
---|
| 1081 | rsX = Calcrs(EC,pDown);
|
---|
| 1082 | break;
|
---|
| 1083 | case SpinDouble:
|
---|
| 1084 | //rsX = Calcrs(EC,p/2.); // why half of it???
|
---|
| 1085 | rsX = Calcrs(EC,p);
|
---|
| 1086 | break;
|
---|
| 1087 | }
|
---|
| 1088 | rsC = Calcrs(EC,p);
|
---|
| 1089 | zeta = CalcZeta(EC,pUp,pDown);
|
---|
| 1090 | SumExc += (CalcSVVXr(EC,rsX,Psi->PsiST) + CalcSVVCr(EC, rsC,zeta,Psi->PsiST,pUp,pDown))*PsiCD[i]*PsiCD[i];
|
---|
| 1091 | }
|
---|
| 1092 | return (SumExc*Factor);
|
---|
| 1093 | }
|
---|
| 1094 |
|
---|
| 1095 | /** Initialises ExCor structure with parametrization values.
|
---|
| 1096 | * All of the entries in the structure are set to parametrization values, see [CA80].
|
---|
| 1097 | * \param *P Problem at hand
|
---|
| 1098 | * \param *EC Exchange correlation ExCor structure
|
---|
| 1099 | */
|
---|
| 1100 | void InitExchangeCorrelationEnergy(struct Problem *P, struct ExCor *EC)
|
---|
| 1101 | {
|
---|
| 1102 | EC->gamma[unpolarised] = -0.1423;
|
---|
| 1103 | EC->beta_1[unpolarised] = 1.0529;
|
---|
| 1104 | EC->beta_2[unpolarised] = 0.3334;
|
---|
| 1105 | EC->C[unpolarised] = 0.0020;
|
---|
| 1106 | EC->D[unpolarised] = -0.0116;
|
---|
| 1107 | EC->A[unpolarised] = 0.0311;
|
---|
| 1108 | EC->B[unpolarised] = -0.048;
|
---|
| 1109 | EC->gamma[polarised] = -0.0843;
|
---|
| 1110 | EC->beta_1[polarised] = 1.3981;
|
---|
| 1111 | EC->beta_2[polarised] = 0.2611;
|
---|
| 1112 | EC->C[polarised] = 0.0007;
|
---|
| 1113 | EC->D[polarised] = -0.0048;
|
---|
| 1114 | EC->A[polarised] = 0.01555;
|
---|
| 1115 | EC->B[polarised] = -0.0269;
|
---|
| 1116 | EC->fac34pi = -3./(4.*PI);
|
---|
| 1117 | EC->facexrs = EC->fac34pi*pow(9.*PI/4.,1./3.);
|
---|
| 1118 | EC->epsilon0 = MYEPSILON;
|
---|
| 1119 | EC->fac6PI23 = pow(6./PI,2./3.);
|
---|
| 1120 | EC->facPI213 = pow(PI/2.,1./3.);
|
---|
| 1121 | EC->fac3PI23 = pow(3.*PI,2./3.);
|
---|
| 1122 | EC->fac6PIPI23 = pow(6.*PI*PI,2./3.);
|
---|
| 1123 | EC->fac243 = pow(2.,4./3.);
|
---|
| 1124 | EC->fac1213 = pow(1./2.,1./3.);
|
---|
| 1125 | }
|
---|