source: pcp/src/excor.c@ 774ae8

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

-initial commit
-Minimum set of files needed from ESPACK SVN repository
-Switch to three tantamount package parts instead of all relating to pcp (as at some time Ralf's might find inclusion as well)
  • Property mode set to 100644
File size: 41.5 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 = /*EC->fac1213**/VCR_unpol;
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 rs, p = 0.0, pUp = 0.0, pDown = 0.0, zeta, rsUp, rsDown, SumEx=0.0;
738 double Factor = R->XCEnergyFactor/Lev0->MaxN;
739 int i;
740
741 for (i = 0; i < Dens->LocalSizeR; i++) { // for each node in radial mesh
742 // put (corecorrected) densities in p, pUp, pDown
743 p = Dens->DensityArray[TotalDensity][i];
744 if (R->CurrentMin > UnOccupied)
745 if (PP->corecorr == CoreCorrected)
746 p += Dens->DensityArray[CoreWaveDensity][i];
747 switch (Psi->PsiST) {
748 case SpinDouble:
749 pUp = 0.5*p;
750 pDown = 0.5*p;
751 break;
752 case SpinUp:
753 case SpinDown:
754 pUp = Dens->DensityArray[TotalUpDensity][i];
755 pDown = Dens->DensityArray[TotalDownDensity][i];
756 if (PP->corecorr == CoreCorrected) {
757 pUp += 0.5*Dens->DensityArray[CoreWaveDensity][i];
758 pDown += 0.5*Dens->DensityArray[CoreWaveDensity][i];
759 }
760 break;
761 default:
762 ;
763 }
764 // set all to zero if one of them is negative
765 if ((p < 0) || (pUp < 0) || (pDown < 0)) {
766 /*fprintf(stderr,"index %i pc %g p %g\n",i,Dens->DensityArray[CoreWaveDensity][i],p);*/
767 p = 0.0;
768 pUp = 0.0;
769 pDown = 0.0;
770 }
771 // Calculation with the densities and summation
772 rs = Calcrs(EC,p);
773 rsUp = Calcrs(EC,pUp);
774 rsDown = Calcrs(EC,pDown);
775 SumEx += CalcSEXr(EC, rsUp, rsDown, pUp, pDown);
776 zeta = CalcZeta(EC,pUp,pDown);
777 SumEc += CalcSECr(EC, rs, zeta, p);
778 }
779 E->AllLocalDensityEnergy[CorrelationEnergy] = Factor*SumEc;
780 E->AllLocalDensityEnergy[ExchangeEnergy] = Factor*SumEx;
781}
782
783/** SpinType-dependent calculation of exchange and correlation with free spin-polarisation energy with riemann tensor.
784 * Discretely does the integration of \f${\cal E}_{xc} (n,\zeta) \cdot n\f$ on
785 * the radial mesh of the respective real densities.<br>
786 *
787 * Like CalculateXCEnergyNoRT(), only CalcSEXr() and CalcSECr() are
788 * divided by RTFactor Lattice->RT.DensityR
789 * \param *P Problem at hand
790 */
791void CalculateXCEnergyUseRT(struct Problem *P)
792{
793 struct Lattice *Lat = &P->Lat;
794 struct Energy *E = Lat->E;
795 struct PseudoPot *PP = &P->PP;
796 struct RunStruct *R = &P->R;
797 struct LatticeLevel *Lev0 = R->Lev0;
798 struct Psis *Psi = &Lat->Psi;
799 struct Density *Dens = Lev0->Dens;
800 struct ExCor *EC = &P->ExCo;
801 double SumEc = 0;
802 double rs, p = 0.0, pUp = 0.0, pDown = 0.0, zeta, rsUp, rsDown, SumEx = 0.0;
803 double Factor = R->XCEnergyFactor/Lev0->MaxN;
804 int i;
805 fftw_real *RTFactor = Lat->RT.DensityR[RTADetPreRT];
806
807 for (i = 0; i < Dens->LocalSizeR; i++) { // for each point of radial mesh take density p[i]
808 p = Dens->DensityArray[TotalDensity][i];
809 if (PP->corecorr == CoreCorrected)
810 p += Dens->DensityArray[CoreWaveDensity][i];
811 switch (Psi->PsiST) {
812 case SpinDouble:
813 pUp = 0.5*p;
814 pDown = 0.5*p;
815 break;
816 case SpinUp:
817 case SpinDown:
818 pUp = Dens->DensityArray[TotalUpDensity][i];
819 pDown = Dens->DensityArray[TotalDownDensity][i];
820 if (PP->corecorr == CoreCorrected) {
821 pUp += 0.5*Dens->DensityArray[CoreWaveDensity][i];
822 pDown += 0.5*Dens->DensityArray[CoreWaveDensity][i];
823 }
824 break;
825 default:
826 ;
827 }
828 if ((p < 0) || (pUp < 0) || (pDown < 0)) {
829 /*fprintf(stderr,"index %i pc %g p %g\n",i,Dens->DensityArray[CoreWaveDensity][i],p);*/
830 p = 0.0;
831 pUp = 0.0;
832 pDown = 0.0;
833 }
834 // .. calculate Ec and Ex and sum them up ...
835 rs = Calcrs(EC,p);
836 rsUp = Calcrs(EC,pUp);
837 rsDown = Calcrs(EC,pDown);
838 SumEx += CalcSEXr(EC, rsUp, rsDown, pUp, pDown)/fabs(RTFactor[i]);
839 zeta = CalcZeta(EC,pUp,pDown);
840 SumEc += CalcSECr(EC, rs, zeta, p)/fabs(RTFactor[i]);
841 }
842 // ... and factorise with discrete integration width
843 E->AllLocalDensityEnergy[CorrelationEnergy] = Factor*SumEc;
844 E->AllLocalDensityEnergy[ExchangeEnergy] = Factor*SumEx;
845}
846
847/** Calculates SpinType-dependently exchange correlation potential with free spin-polarisation without Riemann tensor.
848 * Goes through all possible nodes, calculates the potential and stores it in \a *HGR
849 * \param *P Problem at hand
850 * \param *HGR pointer storage array for (added!) result: \f$V_c (n,\zeta)(R) + V_x (n,\zeta)(R)\f$
851 * \sa CalcSVXr(), CalcSVCr()
852 */
853void CalculateXCPotentialNoRT(struct Problem *P, fftw_real *HGR)
854{
855 struct Lattice *Lat = &P->Lat;
856 struct PseudoPot *PP = &P->PP;
857 struct RunStruct *R = &P->R;
858 struct LatticeLevel *Lev0 = R->Lev0;
859 struct LatticeLevel *LevS = R->LevS;
860 struct Psis *Psi = &Lat->Psi;
861 struct Density *Dens = Lev0->Dens;
862 struct ExCor *EC = &P->ExCo;
863 double rsX, rsC, zeta, p = 0.0, pUp = 0.0, pDown = 0.0;
864 int nx,ny,nz,i;
865 const int Nx = LevS->Plan0.plan->local_nx;
866 const int Ny = LevS->Plan0.plan->N[1];
867 const int Nz = LevS->Plan0.plan->N[2];
868 const int NUpx = LevS->NUp[0]; // factors due to density being calculated on a finer grid
869 const int NUpy = LevS->NUp[1];
870 const int NUpz = LevS->NUp[2];
871 double tmp;
872 for (nx=0;nx<Nx;nx++)
873 for (ny=0;ny<Ny;ny++)
874 for (nz=0;nz<Nz;nz++) {
875 i = nz*NUpz+Nz*NUpz*(ny*NUpy+Ny*NUpy*nx*NUpx);
876 p = Dens->DensityArray[TotalDensity][i];
877 if (PP->corecorr == CoreCorrected)
878 p += Dens->DensityArray[CoreWaveDensity][i];
879 switch (Psi->PsiST) {
880 case SpinDouble:
881 pUp = 0.5*p;
882 pDown = 0.5*p;
883 break;
884 case SpinUp:
885 case SpinDown:
886 pUp = Dens->DensityArray[TotalUpDensity][i];
887 pDown = Dens->DensityArray[TotalDownDensity][i];
888 if (PP->corecorr == CoreCorrected) { // additional factor due to PseudoPot'entials
889 pUp += 0.5*Dens->DensityArray[CoreWaveDensity][i];
890 pDown += 0.5*Dens->DensityArray[CoreWaveDensity][i];
891 }
892 break;
893 default:
894 ;
895 }
896 if ((p < 0) || (pUp < 0) || (pDown < 0)) {
897 p = 0;
898 pUp = 0;
899 pDown = 0;
900 }
901 switch (Psi->PsiST) {
902 case SpinUp:
903 rsX = Calcrs(EC, pUp);
904 break;
905 case SpinDown:
906 rsX = Calcrs(EC, pDown);
907 break;
908 case SpinDouble:
909 //rsX = Calcrs(EC,p/2.); // why half of it???
910 rsX = Calcrs(EC, p);
911 break;
912 default:
913 rsX = 0.0;
914 }
915 rsC = Calcrs(EC,p);
916 zeta = CalcZeta(EC,pUp,pDown);
917 switch (R->CurrentMin) {
918 case UnOccupied: // here epsilon appears instead of the potential in the integrand due to different variation
919 tmp = CalcSEXr(EC,Calcrs(EC, pUp),Calcrs(EC, pDown),pUp,pDown)/p;
920 //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i un= NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
921 tmp += CalcSECr(EC, rsX,zeta, 1);
922 //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i un+= NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
923 break;
924 default:
925 tmp = CalcSVXr(EC,rsX,Psi->PsiST);
926 //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i def= NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
927 tmp += CalcSVCr(EC, rsC,zeta,Psi->PsiST);
928 //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i def+= NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
929 break;
930 }
931 //if (isnan(tmp)) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): tmp_%i := NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
932 HGR[i] += tmp;
933 //if (isnan(HGR[i])) { fprintf(stderr,"WARNGING: CalculateXCPotentialNoRT(): HGR[%i] = NaN!\n", i); Error(SomeError, "NaN-Fehler!"); }
934 }
935}
936
937/** Calculates SpinType-dependently derivative of exchange correlation potential with free spin-polarisation without Riemann tensor.
938 * \f[
939 * A^{XC} = \sum_R \frac{\delta V^{XC}}{\delta n} \rho^2 \qquad (\textnormal{section 5.1, line search})
940 * \f]
941 * With the density from Dens::DensityArray calculates discretely the integral over the summed
942 * derivative, SpinType is taken from Psi::PsiST. Due to the principle of the theory the wave
943 * functions themselves are not needed explicitely, see also CalculateXCEnergyNoRT()
944 * \param *P Problem at hand
945 * \param *PsiCD \f$\rho (r)\f$
946 * \return \f$A^{XC}\f$
947 * \sa CalcSVVXr(), CalcSVVCr() - derivatives of exchange and correlation potential
948 */
949double CalculateXCddEddt0NoRT(struct Problem *P, fftw_real *PsiCD)
950{
951 struct Lattice *Lat = &P->Lat;
952 struct PseudoPot *PP = &P->PP;
953 struct RunStruct *R = &P->R;
954 struct LatticeLevel *Lev0 = R->Lev0;
955 struct Psis *Psi = &Lat->Psi;
956 struct Density *Dens = Lev0->Dens;
957 struct ExCor *EC = &P->ExCo;
958 double SumExc = 0;
959 double rsX=0.0, rsC, p = 0.0, pUp = 0.0, pDown = 0.0, zeta;
960 double Factor = R->XCEnergyFactor/Lev0->MaxN;
961 int i;
962
963 for (i = 0; i < Dens->LocalSizeR; i++) {
964 p = Dens->DensityArray[TotalDensity][i];
965 if (PP->corecorr == CoreCorrected)
966 p += Dens->DensityArray[CoreWaveDensity][i];
967 switch (Psi->PsiST) {
968 case SpinDouble:
969 pUp = 0.5*p;
970 pDown = 0.5*p;
971 break;
972 case SpinUp:
973 case SpinDown:
974 pUp = Dens->DensityArray[TotalUpDensity][i];
975 pDown = Dens->DensityArray[TotalDownDensity][i];
976 if (PP->corecorr == CoreCorrected) {
977 pUp += 0.5*Dens->DensityArray[CoreWaveDensity][i];
978 pDown += 0.5*Dens->DensityArray[CoreWaveDensity][i];
979 }
980 break;
981 default:
982 ;
983 }
984 if ((p < 0) || (pUp < 0) || (pDown < 0)) {
985 /*fprintf(stderr,"index %i pc %g p %g\n",i,Dens->DensityArray[CoreWaveDensity][i],p);*/
986 p = 0.0;
987 pUp = 0.0;
988 pDown = 0.0;
989 }
990 switch (Psi->PsiST) {
991 case SpinUp:
992 rsX = Calcrs(EC,pUp);
993 break;
994 case SpinDown:
995 rsX = Calcrs(EC,pDown);
996 break;
997 case SpinDouble:
998 //rsX = Calcrs(EC,p/2.); // why half of it???
999 rsX = Calcrs(EC,p);
1000 break;
1001 }
1002 rsC = Calcrs(EC,p);
1003 zeta = CalcZeta(EC,pUp,pDown);
1004 SumExc += (CalcSVVXr(EC,rsX,Psi->PsiST) + CalcSVVCr(EC, rsC,zeta,Psi->PsiST,pUp,pDown))*PsiCD[i]*PsiCD[i];
1005 }
1006 return (SumExc*Factor);
1007}
1008
1009/** Initialises ExCor structure with parametrization values.
1010 * All of the entries in the structure are set to parametrization values, see [CA80].
1011 * \param *P Problem at hand
1012 * \param *EC Exchange correlation ExCor structure
1013 */
1014void InitExchangeCorrelationEnergy(struct Problem *P, struct ExCor *EC)
1015{
1016 EC->gamma[unpolarised] = -0.1423;
1017 EC->beta_1[unpolarised] = 1.0529;
1018 EC->beta_2[unpolarised] = 0.3334;
1019 EC->C[unpolarised] = 0.0020;
1020 EC->D[unpolarised] = -0.0116;
1021 EC->A[unpolarised] = 0.0311;
1022 EC->B[unpolarised] = -0.048;
1023 EC->gamma[polarised] = -0.0843;
1024 EC->beta_1[polarised] = 1.3981;
1025 EC->beta_2[polarised] = 0.2611;
1026 EC->C[polarised] = 0.0007;
1027 EC->D[polarised] = -0.0048;
1028 EC->A[polarised] = 0.01555;
1029 EC->B[polarised] = -0.0269;
1030 EC->fac34pi = -3./(4.*PI);
1031 EC->facexrs = EC->fac34pi*pow(9.*PI/4.,1./3.);
1032 EC->epsilon0 = MYEPSILON;
1033 EC->fac6PI23 = pow(6./PI,2./3.);
1034 EC->facPI213 = pow(PI/2.,1./3.);
1035 EC->fac3PI23 = pow(3.*PI,2./3.);
1036 EC->fac6PIPI23 = pow(6.*PI*PI,2./3.);
1037 EC->fac243 = pow(2.,4./3.);
1038 EC->fac1213 = pow(1./2.,1./3.);
1039}
Note: See TracBrowser for help on using the repository browser.