source: pcp/src/excor.c@ b70503

Last change on this file since b70503 was b70503, checked in by Frederik Heber <heber@…>, 17 years ago

sinh_inverse(), CalcSE_GC() and DensityGradient() added for Gradient Correction XC
  • Property mode set to 100644
File size: 44.4 KB
Line 
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
49inline double Calcrs(struct ExCor *EC, double p) {
50#else
51double 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
68inline static double CalcVXrUP(struct ExCor *EC, double rs) {
69#else
70static 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
83inline static double CalcVVXrUP(struct ExCor *EC, double rs) {
84#else
85static 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
97inline static double CalcVVVXrUP(struct ExCor *EC, double rs) {
98#else
99static 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
118inline double CalcEXr(struct ExCor *EC, double rs, double zeta) {
119#else
120double 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
140inline static double CalcEXrUP(struct ExCor *EC, double rs) {
141#else
142static 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
163inline double CalcECr(struct ExCor *EC, double rs, enum UnPolarised up) {
164#else
165double 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
195inline static double CalcVCr(struct ExCor *EC, double rs, enum UnPolarised up) {
196#else
197static 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
223inline static double CalcVVCr(struct ExCor *EC, double rs, enum UnPolarised up) {
224#else
225static 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
258inline static double CalcVVVCr(struct ExCor *EC, double rs, enum UnPolarised up) {
259#else
260static 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
303inline double CalcZeta(struct ExCor *EC, double pUp, double pDown) {
304#else
305double 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
320inline static double CalcDzeta(struct ExCor *EC, enum SpinType ST, double pUp, double pDown) {
321#else
322static 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
350inline static double CalcD2zeta(struct ExCor *EC, enum SpinType ST, double pUp, double pDown) {
351#else
352static 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
380inline static double CalcD3zeta(struct ExCor *EC, enum SpinType ST, double pUp, double pDown) {
381#else
382static 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
409inline static double Calcf(struct ExCor *EC, double zeta) {
410#else
411static 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
428inline static double CalcDf(struct ExCor *EC, double zeta) {
429#else
430static 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
444inline static double CalcD2f(struct ExCor *EC, double zeta) {
445#else
446static 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)
452 return((4./9.)*(1./cbrt((zeta+1.)*(zeta+1.))+1./cbrt((1.-zeta)*(1.-zeta)))/(EC->fac243-2.)); // formula checked (16.5.06)
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
462inline static double CalcD3f(struct ExCor *EC, double zeta) {
463#else
464static 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
482inline double CalcSECr(struct ExCor *EC, double rs, double zeta, double p) {
483#else
484double 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
506inline static double CalcSVCr(struct ExCor *EC, double rs, double zeta, enum SpinType ST) {
507#else
508static 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:
520 res = VCR_unpol; /*EC->fac1213**/
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
536inline double CalcSEXr(struct ExCor *EC, double rsUp, double rsDown, double pUp, double pDown) {
537#else
538double 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
551inline static double CalcSVXr(struct ExCor *EC, double rs, enum SpinType ST) {
552#else
553static 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
576inline double CalcSVVXr(struct ExCor *EC, double rs, enum SpinType ST) {
577#else
578double 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
601inline double CalcSVVVXr(struct ExCor *EC, double rs, enum SpinType ST) {
602#else
603double 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
632inline double CalcSVVCr(struct ExCor *EC, double rs, double zeta, enum SpinType ST, double pUp, double pDown) {
633#else
634double 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
680inline double CalcSVVVCr(struct ExCor *EC, double rs, double zeta, enum SpinType ST, double pUp, double pDown) {
681#else
682double 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 */
726void 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;
737 //double SumEx_GC = 0.; // Gradient correction part according to Becke '92
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;
740 //double Dp, DpUp, DpDown;
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];
746 //Dp = DensityGradient(Dens->DensityArray[TotalDensity], i, Lev0, Lat);
747 if (R->CurrentMin > UnOccupied)
748 if (PP->corecorr == CoreCorrected)
749 p += Dens->DensityArray[CoreWaveDensity][i];
750 switch (Psi->PsiST) {
751 default:
752 case SpinDouble:
753 pUp = 0.5*p;
754 pDown = 0.5*p;
755 //DpUp = 0.5*Dp;
756 //DpDown = 0.5*Dp;
757 break;
758 case SpinUp:
759 case SpinDown:
760 pUp = Dens->DensityArray[TotalUpDensity][i];
761 pDown = Dens->DensityArray[TotalDownDensity][i];
762 //DpUp = DensityGradient(Dens->DensityArray[TotalUpDensity], i, Lev0, Lat);
763 //DpDown = DensityGradient(Dens->DensityArray[TotalDownDensity], i, Lev0, Lat);
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);
784 //SumEx_GC += CalcSE_GC(EC, pUp, DpUp);
785 //SumEx_GC += CalcSE_GC(EC, pDown, DpDown);
786 }
787 E->AllLocalDensityEnergy[CorrelationEnergy] = Factor*SumEc;
788 E->AllLocalDensityEnergy[ExchangeEnergy] = Factor*(SumEx); // - SumEx_GC);
789}
790
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 */
797inline 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 */
808double 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 */
825double 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
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 */
872void 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 */
934void 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];
958 if (PP->corecorr == CoreCorrected)
959 p += Dens->DensityArray[CoreWaveDensity][i];
960 switch (Psi->PsiST) {
961 case SpinDouble:
962 pUp = 0.5*p;
963 pDown = 0.5*p;
964 break;
965 case SpinUp:
966 case SpinDown:
967 pUp = Dens->DensityArray[TotalUpDensity][i];
968 pDown = Dens->DensityArray[TotalDownDensity][i];
969 if (PP->corecorr == CoreCorrected) { // additional factor due to PseudoPot'entials
970 pUp += 0.5*Dens->DensityArray[CoreWaveDensity][i];
971 pDown += 0.5*Dens->DensityArray[CoreWaveDensity][i];
972 }
973 break;
974 default:
975 ;
976 }
977 if ((p < 0) || (pUp < 0) || (pDown < 0)) {
978 p = 0;
979 pUp = 0;
980 pDown = 0;
981 }
982 switch (Psi->PsiST) {
983 case SpinUp:
984 rsX = Calcrs(EC, pUp);
985 break;
986 case SpinDown:
987 rsX = Calcrs(EC, pDown);
988 break;
989 case SpinDouble:
990 //rsX = Calcrs(EC,p/2.); // why half of it???
991 rsX = Calcrs(EC, p);
992 break;
993 default:
994 rsX = 0.0;
995 }
996 rsC = Calcrs(EC,p);
997 zeta = CalcZeta(EC,pUp,pDown);
998 switch (R->CurrentMin) {
999 case UnOccupied: // here epsilon appears instead of the potential in the integrand due to different variation
1000 tmp = CalcSEXr(EC,Calcrs(EC, pUp),Calcrs(EC, pDown),pUp,pDown)/p;
1001 //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i un= NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
1002 tmp += CalcSECr(EC, rsX,zeta, 1);
1003 //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i un+= NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
1004 break;
1005 default:
1006 tmp = CalcSVXr(EC,rsX,Psi->PsiST);
1007 //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i def= NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
1008 tmp += CalcSVCr(EC, rsC,zeta,Psi->PsiST);
1009 //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i def+= NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
1010 break;
1011 }
1012 //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i := NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
1013 HGR[i] += tmp;
1014 //if (isnan(HGR[i])) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): HGR[%i] = NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
1015 }
1016}
1017
1018/** Calculates SpinType-dependently derivative of exchange correlation potential with free spin-polarisation without Riemann tensor.
1019 * \f[
1020 * A^{XC} = \sum_R \frac{\delta V^{XC}}{\delta n} \rho^2 \qquad (\textnormal{section 5.1, line search})
1021 * \f]
1022 * With the density from Dens::DensityArray calculates discretely the integral over the summed
1023 * derivative, SpinType is taken from Psi::PsiST. Due to the principle of the theory the wave
1024 * functions themselves are not needed explicitely, see also CalculateXCEnergyNoRT()
1025 * \param *P Problem at hand
1026 * \param *PsiCD \f$\rho (r)\f$
1027 * \return \f$A^{XC}\f$
1028 * \sa CalcSVVXr(), CalcSVVCr() - derivatives of exchange and correlation potential
1029 */
1030double CalculateXCddEddt0NoRT(struct Problem *P, fftw_real *PsiCD)
1031{
1032 struct Lattice *Lat = &P->Lat;
1033 struct PseudoPot *PP = &P->PP;
1034 struct RunStruct *R = &P->R;
1035 struct LatticeLevel *Lev0 = R->Lev0;
1036 struct Psis *Psi = &Lat->Psi;
1037 struct Density *Dens = Lev0->Dens;
1038 struct ExCor *EC = &P->ExCo;
1039 double SumExc = 0;
1040 double rsX=0.0, rsC, p = 0.0, pUp = 0.0, pDown = 0.0, zeta;
1041 double Factor = R->XCEnergyFactor/Lev0->MaxN;
1042 int i;
1043
1044 for (i = 0; i < Dens->LocalSizeR; i++) {
1045 p = Dens->DensityArray[TotalDensity][i];
1046 if (PP->corecorr == CoreCorrected)
1047 p += Dens->DensityArray[CoreWaveDensity][i];
1048 switch (Psi->PsiST) {
1049 case SpinDouble:
1050 pUp = 0.5*p;
1051 pDown = 0.5*p;
1052 break;
1053 case SpinUp:
1054 case SpinDown:
1055 pUp = Dens->DensityArray[TotalUpDensity][i];
1056 pDown = Dens->DensityArray[TotalDownDensity][i];
1057 if (PP->corecorr == CoreCorrected) {
1058 pUp += 0.5*Dens->DensityArray[CoreWaveDensity][i];
1059 pDown += 0.5*Dens->DensityArray[CoreWaveDensity][i];
1060 }
1061 break;
1062 default:
1063 ;
1064 }
1065 if ((p < 0) || (pUp < 0) || (pDown < 0)) {
1066 /*fprintf(stderr,"index %i pc %g p %g\n",i,Dens->DensityArray[CoreWaveDensity][i],p);*/
1067 p = 0.0;
1068 pUp = 0.0;
1069 pDown = 0.0;
1070 }
1071 switch (Psi->PsiST) {
1072 case SpinUp:
1073 rsX = Calcrs(EC,pUp);
1074 break;
1075 case SpinDown:
1076 rsX = Calcrs(EC,pDown);
1077 break;
1078 case SpinDouble:
1079 //rsX = Calcrs(EC,p/2.); // why half of it???
1080 rsX = Calcrs(EC,p);
1081 break;
1082 }
1083 rsC = Calcrs(EC,p);
1084 zeta = CalcZeta(EC,pUp,pDown);
1085 SumExc += (CalcSVVXr(EC,rsX,Psi->PsiST) + CalcSVVCr(EC, rsC,zeta,Psi->PsiST,pUp,pDown))*PsiCD[i]*PsiCD[i];
1086 }
1087 return (SumExc*Factor);
1088}
1089
1090/** Initialises ExCor structure with parametrization values.
1091 * All of the entries in the structure are set to parametrization values, see [CA80].
1092 * \param *P Problem at hand
1093 * \param *EC Exchange correlation ExCor structure
1094 */
1095void InitExchangeCorrelationEnergy(struct Problem *P, struct ExCor *EC)
1096{
1097 EC->gamma[unpolarised] = -0.1423;
1098 EC->beta_1[unpolarised] = 1.0529;
1099 EC->beta_2[unpolarised] = 0.3334;
1100 EC->C[unpolarised] = 0.0020;
1101 EC->D[unpolarised] = -0.0116;
1102 EC->A[unpolarised] = 0.0311;
1103 EC->B[unpolarised] = -0.048;
1104 EC->gamma[polarised] = -0.0843;
1105 EC->beta_1[polarised] = 1.3981;
1106 EC->beta_2[polarised] = 0.2611;
1107 EC->C[polarised] = 0.0007;
1108 EC->D[polarised] = -0.0048;
1109 EC->A[polarised] = 0.01555;
1110 EC->B[polarised] = -0.0269;
1111 EC->fac34pi = -3./(4.*PI);
1112 EC->facexrs = EC->fac34pi*pow(9.*PI/4.,1./3.);
1113 EC->epsilon0 = MYEPSILON;
1114 EC->fac6PI23 = pow(6./PI,2./3.);
1115 EC->facPI213 = pow(PI/2.,1./3.);
1116 EC->fac3PI23 = pow(3.*PI,2./3.);
1117 EC->fac6PIPI23 = pow(6.*PI*PI,2./3.);
1118 EC->fac243 = pow(2.,4./3.);
1119 EC->fac1213 = pow(1./2.,1./3.);
1120}
Note: See TracBrowser for help on using the repository browser.