1 /**********************************************************************
2 XC_Ceperly_Alder.c:
3
4 XC_Ceperly_Alder.c is a subroutine to calculate an exchange-
5 correlation potential for a given density "den" by the local
6 density approximation, which is based on the original works
7 by Ceperly and Alder and parametrized by Perdew and Zunger.
8
9 Log of XC_Ceperly_Alder.c:
10
11 22/Nov/2001 Released by T.Ozaki
12
13 ***********************************************************************/
14
15 #include <stdio.h>
16 #include <math.h>
17 #include "openmx_common.h"
18
XC_Ceperly_Alder(double den,int P_switch)19 double XC_Ceperly_Alder(double den, int P_switch)
20 {
21
22 /****************************************************
23 P_switch:
24 0 \epsilon_XC (XC energy density)
25 1 \mu_XC (XC potential)
26 2 \epsilon_XC - \mu_XC
27 3 derivative of XC energy density w.r.t. den
28 ****************************************************/
29
30 double dum,rs,coe;
31 double Ex,Ec,dEx,dEc;
32 double tmp0,tmp1;
33 double result;
34
35 /****************************************************
36 Non-relativisic
37 ****************************************************/
38
39 if (den<=1.0e-15){
40 result = 0.0;
41 }
42 else{
43
44 coe = 0.6203504908994; /* pow(3.0/4.0/PI,1.0/3.0); */
45 rs = coe*pow(den,-0.3333333333333333333);
46
47 tmp0 = 0.458165293632163/rs;
48 Ex = -tmp0;
49 dEx = tmp0/rs;
50
51 if (1.0<=rs){
52 tmp0 = sqrt(rs);
53 dum = (1.0 + 1.0529*tmp0 + 0.3334*rs);
54 tmp1 = 0.1423/dum;
55 Ec = -tmp1;
56 dEc = tmp1/dum*(0.52645/tmp0 + 0.3334);
57 }
58 else{
59 tmp0 = log(rs);
60 Ec = -0.0480 + 0.0311*tmp0 + rs*(0.0020*tmp0 - 0.0116);
61 dEc = 0.0311/rs + 0.0020*tmp0 - 0.0096;
62 }
63
64 /*
65 printf("Ex=%15.12f %15.12f\n",Ex,Ex-0.33333333333333333333*rs*dEx);
66 */
67
68 if (P_switch==0)
69 result = Ex + Ec;
70 else if (P_switch==1)
71 result = Ex + Ec - 0.33333333333333333333*rs*(dEx + dEc);
72 else if (P_switch==2)
73 result = 0.3333333333333333333*rs*(dEx + dEc);
74 else if (P_switch==3)
75 result = -0.3333333333333333333/(coe*coe*coe)*rs*rs*rs*rs*(dEx + dEc);
76
77 }
78
79 return result;
80 }
81
82
83
84
85
86
87
88
89
90
91
92