1 /*
2 $Id$
3 dirac_exchange.c - 6/9/95
4 author - Eric Bylaska
5
6 This file contains a routine for finding the Dirac
7 exchange potential and energy of a spin density rho, defined on
8 a log grid.
9
10 */
11
12 #include <stdio.h>
13 #include "loggrid.h"
14 #include "dirac_exchange.h"
15
16 static double alpha=0.6666666666666666667;
17
18 /********************************
19 * *
20 * set_Dirac_alpha *
21 * *
22 ********************************/
23
set_Dirac_alpha(aa)24 void set_Dirac_alpha(aa)
25 double aa;
26 {
27 alpha = aa;
28 }
29
30 /********************************
31 * *
32 * Dirac_alpha *
33 * *
34 ********************************/
35
Dirac_alpha()36 double Dirac_alpha()
37 {
38
39 return alpha;
40 }
41
42
43
44 /********************************
45 * *
46 * R_Dirac_Exchange *
47 * *
48 ********************************/
49
50 /* this routine calculates the spin
51 polarized Dirac exchange functional.
52
53 Entry - rho[]: the density
54 Exit - Vx[]: the exchange functional
55 Ex: the exchange energy
56 Px: The variational exchange corrections for the eigenvalues.
57 */
58
R_Dirac_Exchange(rho,Vx,Ex,Px)59 void R_Dirac_Exchange(rho,Vx,Ex,Px)
60
61 double rho[],
62 Vx[],
63 *Ex,
64 *Px;
65 {
66 int i;
67 double onethird;
68 double pi;
69 double n;
70 double n_onethird;
71 double ex_p;
72 double ux_p;
73
74 /* loggrid variables */
75 int Ngrid;
76
77 /* tempory local grids */
78 double *ex_functional;
79 double *tmp;
80
81 /* define constants */
82 onethird = 1.0/3.0;
83 pi = 4.0*atan(1.0);
84
85
86 /* access the loggrid variables */
87 Ngrid = N_LogGrid();
88
89 /* allocate temporary memory */
90 ex_functional = alloc_LogGrid();
91 tmp = alloc_LogGrid();
92
93 for (i=0; i<Ngrid; ++i)
94 {
95 n = rho[i]/(4.0*pi);
96 n_onethird = pow((3.0*n/pi),onethird);
97
98 ex_p = -(9.0/8.0)*alpha*n_onethird;
99 ux_p = -(3.0/2.0)*alpha*n_onethird;
100
101 ex_functional[i] = ex_p;
102 Vx[i] = ux_p;
103 } /*for i*/
104
105
106 /* cacluate Ex, and Px */
107 /* note that the integration is weird, because */
108 /* we are integrating from 0 to infinity, and */
109 /* our log grid goes from r0 to 45.0 */
110
111 /* integrate Ex = integrate(rho*ex_functional) */
112 for (i=0; i<Ngrid; ++i)
113 tmp[i] = rho[i]*ex_functional[i];
114 *Ex = Integrate_LogGrid(tmp);
115
116 /* integrate px = integrate(rho*Vx) */
117 for (i=0; i<Ngrid; ++i)
118 tmp[i] = (rho[i])*Vx[i];
119 *Px = Integrate_LogGrid(tmp);
120
121
122 /* deallocate temporary memory */
123 dealloc_LogGrid(ex_functional);
124 dealloc_LogGrid(tmp);
125
126 } /* Dirac_Exchange */
127