1 /**********************************************************************
2 XC_PW92C.c:
3
4 XC_PW92C.c is a subroutine to calculate the correlation
5 energy density and potential developed by Perdew and Wang
6 (Ref: J.P.Perdew and Y.Wang, PRB, 45, 13244 (1992))
7 for given up (dens[0]) and down (dens[1]) densities.
8
9 This routine was written by T.Ozaki, based on the original fortran
10 code provided by the SIESTA group through their website.
11 Thanks to them.
12
13 Log of PW92C.c:
14
15 22/Nov/2001 Released by T.Ozaki
16
17 ***********************************************************************/
18
19 #include <stdio.h>
20 #include <math.h>
21 #include <stdlib.h>
22 #include "openmx_common.h"
23
24 #define den_min 1.0e-14
25 #define den_min_half 0.5*1.0e-14
26
XC_PW92C(double dens[2],double Ec[1],double Vc[2])27 void XC_PW92C(double dens[2], double Ec[1], double Vc[2])
28 {
29 int i;
30 double dtot,rs,srs,zeta,coe;
31 double dum,dum1,dum2,b,c,dbdrs,dcdrs;
32 double tmp0,tmp1,tmp2,tmp12,tmp14,tmp22,tmp24;
33 double fpp0,f,dfdz;
34 double G[3],dGdrs[3];
35 double dEcdd[2];
36 double dEcdrs,dEcdz;
37 double drsdd,dzdd[2];
38
39 /****************************************************
40 parameters from Table I of
41 Perdew and Wang, PRB, 45, 13244 (92)
42 ****************************************************/
43
44 double p[3] = {1.0000000, 1.0000000, 1.0000000};
45 double A[3] = {0.0310910, 0.0155450, 0.0168870};
46 double alpha1[3] = {0.2137000, 0.2054800, 0.1112500};
47 double beta1[3] = {7.5957000, 14.1189000, 10.3570000};
48 double beta2[3] = {3.5876000, 6.1977000, 3.6231000};
49 double beta3[3] = {1.6382000, 3.3662000, 0.8802600};
50 double beta4[3] = {0.4929400, 0.6251700, 0.4967100};
51
52 /****************************************************
53 zeta and rs
54 ****************************************************/
55
56 coe = 0.6203504908994; /* pow(3.0/4.0/PI,1.0/3.0); */
57 dtot = dens[0] + dens[1];
58
59 if (dtot<=den_min){
60 rs = 6203.504908994; /* coe*pow(1.0e-12,-1.0/3.0); */
61 dens[0] = den_min_half;
62 dens[1] = den_min_half;
63 dtot = den_min;
64 }
65 else
66 rs = coe*pow(dtot,-0.33333333333333333333333);
67
68 tmp0 = 1.0/dtot;
69 zeta = tmp0*(dens[0] - dens[1]);
70
71 if (0.99<zeta) zeta = 0.99;
72 if (zeta<-1.0) zeta = -0.99;
73
74 drsdd = -0.3333333333333333333333*rs*tmp0;
75 dzdd[0] = tmp0*( 1.0 - zeta);
76 dzdd[1] = tmp0*(-1.0 - zeta);
77
78 /****************************************************
79 eps_c(rs,0)=G(0)
80 eps_c(rs,1)=G(1)
81 -alpha_c(rs)=G(2)
82 using eq.(10) in
83 Perdew and Wang, PRB, 45, 13244 (1992))
84 ****************************************************/
85
86 srs = sqrt(rs);
87
88 for (i=0; i<=2; i++){
89 b = beta1[i]*srs + rs*(beta2[i] + beta3[i]*srs + beta4[i]*rs);
90
91 dbdrs = beta1[i]*0.50/srs
92 + beta2[i]
93 + beta3[i]*1.50*srs
94 + beta4[i]*2.0*rs;
95
96 c = 1.0 + 1.0/(2.0*A[i]*b);
97 dcdrs = -(c - 1.0)*dbdrs/b;
98 dum = log(c);
99 dum1 = 1.0 + alpha1[i]*rs;
100 G[i] = -2.0*A[i]*dum1*dum;
101 dGdrs[i] = -2.0*A[i]*(alpha1[i]*dum + dum1*dcdrs/c);
102 }
103
104 /****************************************************
105 f''(0) and f(zeta) from eq.(9)
106 ****************************************************/
107
108 c = 1.92366105093154; /* 1/(2*(2^{1/3}-1)) */
109 fpp0 = 1.70992093416137;
110 dum1 = 1.0 + zeta;
111 dum2 = 1.0 - zeta;
112
113 tmp1 = pow(dum1,0.333333333333333333);
114 tmp2 = pow(dum2,0.333333333333333333);
115
116 tmp12 = tmp1*tmp1;
117 tmp22 = tmp2*tmp2;
118 tmp14 = tmp12*tmp12;
119 tmp24 = tmp22*tmp22;
120
121 f = (tmp14 + tmp24 - 2.0)*c;
122 dfdz = 1.333333333333333333*(tmp1 - tmp2)*c;
123
124 /****************************************************
125 eps_c(rs,zeta) from eq.(8)
126 ****************************************************/
127
128 dum1 = zeta*zeta*zeta;
129 dum = dum1*zeta;
130
131 Ec[0] = G[0] - G[2]*f/fpp0*(1.0 - dum) + (G[1] - G[0])*f*dum;
132
133 dEcdrs = dGdrs[0] - dGdrs[2]*f/fpp0*(1.0 - dum)
134 + (dGdrs[1] - dGdrs[0])*f*dum;
135
136 dEcdz = - G[2]/fpp0*(dfdz*(1.0 - dum) - f*4.0*dum1)
137 + (G[1] - G[0])*(dfdz*dum + f*4.0*dum1);
138
139 /*
140 printf("rs = %18.15f\n",rs);
141 printf("zeta = %18.15f\n",zeta);
142 printf("Ec = %18.15f\n",Ec[0]);
143 printf("Ec2 = %18.15f %18.15f\n",
144 -G[2]*f/fpp0*(1.0 - dum),-dGdrs[2]*f/fpp0*(1.0 - dum));
145 printf("Ec3 = %18.15f %18.15f\n",
146 (G[1] - G[0])*f*dum,(dGdrs[1] - dGdrs[0])*f*dum);
147
148 printf("G0 dGdrs0 = %18.15f %18.15f\n",G[0],dGdrs[0]);
149 printf("G1 dGdrs1 = %18.15f %18.15f\n",G[1],dGdrs[1]);
150 printf("G2 dGdrs2 = %18.15f %18.15f\n",G[2],dGdrs[2]);
151
152 printf("dEcdrs = %18.15f\n",dEcdrs);
153 printf("dEcdz = %18.15f\n",dEcdz);
154 */
155
156 /****************************************************
157 Find correlation potential
158 ****************************************************/
159
160 dum = dEcdrs*drsdd;
161 dEcdd[0] = dum + dEcdz*dzdd[0];
162 dEcdd[1] = dum + dEcdz*dzdd[1];
163 Vc[0] = Ec[0] + dtot*dEcdd[0];
164 Vc[1] = Ec[0] + dtot*dEcdd[1];
165
166 }
167