1 /* revpbe_exchange.c - 6/9/95
2 author - Eric Bylaska
3
4 This file contains a routine for finding the PBE GGA
5 exchange potential and energy of a spin density rho, defined on
6 a log grid.
7
8 THis GGA was developed by Norskov
9
10 */
11 #include <stdio.h>
12 #include "loggrid.h"
13 #include "revpbe_exchange.h"
14
15
16
17 /********************************
18 * *
19 * R_revPBE_Exchange *
20 * *
21 ********************************/
22
23 /* this routine calculates the PBE
24 GGA exchange functional.
25
26 Entry - rho[]: the density
27 Exit - Vx[]: the exchange functional
28 Ex: the exchange energy
29 Px[]: The variational exchange corrections for the eigenvalues.
30 */
31
R_revPBE_Exchange(rho,Vx,Ex,Px)32 void R_revPBE_Exchange(rho,Vx,Ex,Px)
33
34 double rho[],
35 Vx[],
36 *Ex,
37 *Px;
38 {
39 int i;
40 double mu,kappa;
41 double onethird,fourthird,twothird;
42 double pi;
43 double n;
44 double n_onethird;
45 double ex_p;
46 double ux_p;
47 double s,F,Fs,Fss,kf,P0;
48 double u,v;
49 double agr,lap,delgr;
50
51 /* loggrid variables */
52 int Ngrid;
53
54 /* tempory local grids */
55 double *ex_functional;
56 double *tmp;
57 double *drho;
58 double *ddrho;
59 double *dadrho;
60 double *rgrid;
61
62 /* define constants */
63 mu = 0.2195149727645171;
64 /*kappa = 0.8040000000000000;*/
65 kappa = 1.2450000000000000;
66 pi = 4.0*atan(1.0);
67 onethird = 1.0/3.0;
68 fourthird = 4.0/3.0;
69 twothird = 2.0/3.0;
70
71
72 /* access the loggrid variables */
73 Ngrid = N_LogGrid();
74 rgrid = r_LogGrid();
75
76 /* allocate temporary memory */
77 ex_functional = alloc_LogGrid();
78 tmp = alloc_LogGrid();
79 drho = alloc_LogGrid();
80 ddrho = alloc_LogGrid();
81 dadrho = alloc_LogGrid();
82
83 /* calculate drho,ddrho, */
84 for (i=0; i<Ngrid; ++i)
85 tmp[i] = rho[i]/(4.0*pi);
86 Derivative_LogGrid(tmp,drho);
87 Derivative_LogGrid(drho,ddrho);
88
89 /* calculate dadrho */
90 for (i=0; i<Ngrid; ++i)
91 tmp[i] = fabs(drho[i]);
92 Derivative_LogGrid(tmp,dadrho);
93
94
95
96 for (i=0; i<Ngrid; ++i)
97 {
98 /* regular inputs to GGA */
99 n = rho[i]/(4.0*pi);
100 if (n > 1.0e-18)
101 {
102 agr = fabs(drho[i]);
103 delgr = drho[i]*dadrho[i];
104 lap = ddrho[i] + (2.0/rgrid[i])*drho[i];
105
106 n_onethird = pow((3.0*n/pi),onethird);
107
108 kf = pow( (3.0*pi*pi*n), onethird);
109 s = agr/(2.0*kf*n);
110 u = delgr/(n*n*(8.0*kf*kf*kf));
111 v = lap/(n*(4.0*kf*kf));
112 P0 = 1.0 + (mu/kappa)*s*s;
113
114
115 F = (1.0 + kappa - kappa/P0);
116 Fs = 2.0*mu/(P0*P0);
117 Fss = -4.0*(mu/kappa)*s*Fs/P0;
118
119 ex_p = -(3.0/4.0)*n_onethird*F;
120 ux_p = -(3.0/4.0)*n_onethird*(
121 fourthird*F
122 - v*Fs
123 - (u - fourthird*s*s*s)*Fss
124 );
125
126 ex_functional[i] = ex_p;
127 Vx[i] = ux_p;
128 }
129 else
130 {
131 ex_functional[i] = 0.0;
132 Vx[i] = 0.0;
133 }
134 } /*for i*/
135
136
137 /* cacluate Ex, and Px */
138 /* note that the integration is weird, because */
139 /* we are integrating from 0 to infinity, and */
140 /* our log grid goes from r0 to 45.0 */
141
142 /* integrate Ex = integrate((rho_down+rho_down)*ex_functional) */
143 for (i=0; i<Ngrid; ++i)
144 tmp[i] = (rho[i])*ex_functional[i];
145 *Ex = Integrate_LogGrid(tmp);
146
147 /* integrate px_up = integrate(rho*Vx) */
148 for (i=0; i<Ngrid; ++i)
149 tmp[i] = (rho[i])*Vx[i];
150 *Px = Integrate_LogGrid(tmp);
151
152
153
154 /* deallocate temporary memory */
155 dealloc_LogGrid(ex_functional);
156 dealloc_LogGrid(tmp);
157 dealloc_LogGrid(drho);
158 dealloc_LogGrid(ddrho);
159 dealloc_LogGrid(dadrho);
160
161 } /* R_revPBE_Exchange */
162 /* $Id$ */
163