1 /**********************************************************************
2   XC_PBE.c:
3 
4      XC_PBE.c is a subroutine to calculate the exchange-correlation
5      potential developed by Perdew, Burke and Ernzerhof within
6      generalized gradient approximation.
7 
8      This routine was written by T.Ozaki, based on the original fortran
9      code provided by the SIESTA group through their website.
10      Thanks to them.
11 
12      Ref: J.P.Perdew, K.Burke & M.Ernzerhof, PRL 77, 3865 (1996)
13 
14   Log of XC_PBE.c:
15 
16      22/Nov/2001  Released by T.Ozaki
17 ***********************************************************************/
18 
19 #include <stdio.h>
20 #include <math.h>
21 #include <stdlib.h>
22 #include "openmx_common.h"
23 
24 #define FOUTHD  (4.0/3.0)
25 #define HALF    0.50
26 #define THD     (1.0/3.0)
27 #define THRHLF  1.50
28 #define TWO     2.0
29 #define TWOTHD  (2.0/3.0)
30 #define beta    0.06672455060314922
31 #define kappa   0.8040
32 
33 #pragma optimization_level 1
XC_PBE(double dens[2],double GDENS[3][2],double Exc[2],double DEXDD[2],double DECDD[2],double DEXDGD[3][2],double DECDGD[3][2])34 void XC_PBE(double dens[2], double GDENS[3][2], double Exc[2],
35             double DEXDD[2], double DECDD[2],
36             double DEXDGD[3][2], double DECDGD[3][2])
37 {
38   int IS,IX;
39   double Ec_unif[1],Vc_unif[2];
40   double dt,rs,zeta;
41   double den_min,gd_min,phi,t,ks,kF,f1,f2,f3,f4;
42   double A,H,Fc,Fx;
43   double GDMT,GDT[3];
44   double DRSDD,DKFDD,DKSDD,DZDD[2],DPDZ;
45   double DECUDD,DPDD,DTDD,DF1DD,DF2DD,DF3DD,DF4DD,DADD;
46   double DHDD,DFCDD[2],DTDGD,DF3DGD,DF4DGD,DHDGD,DFCDGD[3][2];
47   double DS[2],GDMS,KFS,s,f,DFDD,DFXDD[2],Vx_unif[2],Ex_unif[1];
48   double GDS,DSDGD,DSDD,DF1DGD,DFDGD,DFXDGD[3][2];
49   double D[2],GD[3][2],GDM[2];
50   double gamma,mu;
51 
52   gamma = (1.0 - log(TWO))/(PI*PI);
53   mu = beta*PI*PI/3.0;
54 
55   /****************************************************
56          Lower bounds of density and its gradient
57               to avoid divisions by zero
58   ****************************************************/
59 
60   den_min = 1.0e-14;
61   gd_min  = 1.0e-14;
62 
63   /****************************************************
64    Translate density and its gradient to new variables
65   ****************************************************/
66 
67   dens[0] = largest(0.5*den_min,dens[0]);
68   dens[1] = largest(0.5*den_min,dens[1]);
69 
70   D[0] = dens[0];
71   D[1] = dens[1];
72   dt = dens[0] + dens[1];
73 
74   for (IX=0; IX<=2; IX++){
75     GD[IX][0] = GDENS[IX][0];
76     GD[IX][1] = GDENS[IX][1];
77     GDT[IX] = GDENS[IX][0] + GDENS[IX][1];
78   }
79   GDM[0] = sqrt(GD[0][0]*GD[0][0] + GD[1][0]*GD[1][0] + GD[2][0]*GD[2][0]);
80   GDM[1] = sqrt(GD[0][1]*GD[0][1] + GD[1][1]*GD[1][1] + GD[2][1]*GD[2][1]);
81   GDMT   = sqrt(GDT[0]*GDT[0] + GDT[1]*GDT[1] + GDT[2]*GDT[2]);
82   GDMT = largest(gd_min, GDMT);
83 
84   /*
85   printf("GDM0=%15.12f\n",GDM[0]);
86   printf("GDM1=%15.12f\n",GDM[1]);
87   printf("GDMT=%15.12f\n",GDMT);
88   */
89 
90   /****************************************************
91           Local correlation energy and potential
92   ****************************************************/
93 
94   XC_PW92C(dens, Ec_unif, Vc_unif);
95 
96   /*
97   printf("PW92C Ec %15.12f\n",Ec_unif[0]);
98   printf("PW92C Vc %15.12f %15.12f\n",Vc_unif[0],Vc_unif[1]);
99   */
100 
101   /****************************************************
102                 Total correlation energy
103   ****************************************************/
104 
105   rs = pow(3.0/(4.0*PI*dt),THD);
106   kF = pow(3.0*PI*PI*dt,THD);
107   ks = sqrt(4.0*kF/PI);
108   zeta = (dens[0] - dens[1])/dt;
109 
110   if (0.99<zeta) zeta =  0.99;
111   if (zeta<-1.0) zeta = -0.99;
112 
113   /*
114   printf("zeta=%50.45f\n",zeta);
115   */
116 
117   phi = 0.50*(pow(1.0 + zeta,TWOTHD)
118             + pow(1.0 - zeta,TWOTHD));
119   t = GDMT/(2.0*phi*ks*dt);
120   f1 = Ec_unif[0]/(gamma*phi*phi*phi);
121   f2 = exp(-f1);
122 
123   /*
124   printf("ks=%15.12f\n",ks);
125   printf("t=%15.12f\n",t);
126   printf("f2=%15.12f\n",f2);
127   printf("phi^3=%15.12f\n",phi*phi*phi);
128   */
129 
130   A = beta/gamma/(f2 - 1.0);
131 
132   /*
133   printf("A=%15.12f\n",A);
134   */
135 
136   f3 = t*t + A*t*t*t*t;
137   f4 = beta/gamma * f3/(1.0 + A*f3);
138   H = gamma*phi*phi*phi*log(1.0 + f4);
139   Fc = Ec_unif[0] + H;
140 
141   /*
142   printf("At^2=%15.12f\n",A*t*t);
143   printf("A^2t^4=%15.12f\n",A*A*t*t*t*t);
144   printf("beta=%15.12f\n",beta);
145   printf("gamma=%15.12f\n",gamma);
146   printf("t=%15.12f\n",t);
147 
148   { double coe;
149 
150   coe = beta/gamma*t*t;
151   printf("coe=%15.12f\n",beta/gamma);
152   }
153   printf("H=%15.12f\n",H);
154   */
155 
156   /****************************************************
157               Correlation energy derivatives
158   ****************************************************/
159 
160   DRSDD = -(THD*rs/dt);
161   DKFDD =   THD*kF/dt;
162   DKSDD = HALF*ks*DKFDD/kF;
163   DZDD[0] = 1.0/dt - zeta/dt;
164   DZDD[1] = -(1.0/dt) - zeta/dt;
165 
166   DPDZ = HALF*TWOTHD*(1.0/pow(1.0+zeta,THD) - 1.0/pow(1.0-zeta,THD));
167 
168   for (IS=0; IS<=1; IS++){
169     DECUDD = (Vc_unif[IS] - Ec_unif[0])/dt;
170     DPDD = DPDZ*DZDD[IS];
171 
172     /*
173     printf("IS=%2d DPDZ=%15.12f DZDD=%15.12f\n",IS,DPDZ,DZDD[IS]);
174     */
175 
176     DTDD = (-t)*(DPDD/phi + DKSDD/ks + 1.0/dt);
177     DF1DD = f1*(DECUDD/Ec_unif[0] - 3.0*DPDD/phi);
178     DF2DD = (-f2)*DF1DD;
179     DADD = (-A)*DF2DD/(f2 - 1.0);
180     DF3DD = (2.0*t + 4.0*A*t*t*t) * DTDD + DADD*t*t*t*t;
181     DF4DD = f4*(DF3DD/f3 - (DADD*f3+A*DF3DD)/(1.0 + A*f3));
182 
183     /*
184     printf("DPDD=%15.12f phi=%15.12f\n",DPDD,phi);
185     */
186 
187     DHDD = 3.0*H*DPDD/phi;
188     DHDD = DHDD + gamma*phi*phi*phi*DF4DD/(1.0 + f4);
189     DFCDD[IS] = Vc_unif[IS] + H + dt * DHDD;
190 
191     /*
192     printf("IS=%2d Vc_unif=%15.12f H=%15.12f dt=%15.12f DHDD=%15.12f\n",IS,Vc_unif[IS],H,dt,DHDD);
193     */
194 
195     for (IX=0; IX<=2; IX++){
196       DTDGD = (t/GDMT)*GDT[IX]/GDMT;
197       DF3DGD = DTDGD*(2.0*t + 4.0*A*t*t*t);
198       DF4DGD = f4*DF3DGD*(1.0/f3 - A/(1.0 + A*f3));
199       DHDGD = gamma*phi*phi*phi*DF4DGD/(1.0 + f4);
200       DFCDGD[IX][IS] = dt*DHDGD;
201     }
202   }
203 
204   /****************************************************
205               Exchange energy and potential
206   ****************************************************/
207 
208   Fx = 0.0;
209   for (IS=0; IS<=1; IS++){
210 
211     DS[IS] = largest(den_min,2.0*D[IS]);
212     GDMS = largest(gd_min, 2.0*GDM[IS]);
213     KFS = pow(3.0*PI*PI*DS[IS],THD);
214     s = GDMS/(2.0*KFS*DS[IS]);
215     f1 = 1.0 + mu*s*s/kappa;
216     f = 1.0 + kappa - kappa/f1;
217 
218     /****************************************************
219                 Note nspin=1 in call to XC_EX
220     ****************************************************/
221 
222     XC_EX(1, DS[IS], DS, Ex_unif, Vx_unif);
223 
224     Fx = Fx + DS[IS]*Ex_unif[0]*f;
225     DKFDD = THD * KFS/DS[IS];
226     DSDD = s*(-(DKFDD/KFS) - 1.0/DS[IS]);
227     DF1DD = 2.0*(f1 - 1.0)*DSDD/s;
228     DFDD = kappa*DF1DD/(f1*f1);
229     DFXDD[IS] = Vx_unif[0]*f + DS[IS]*Ex_unif[0]*DFDD;
230     for (IX=0; IX<=2; IX++){
231       GDS = 2.0*GD[IX][IS];
232       DSDGD = (s/GDMS)*GDS/GDMS;
233       DF1DGD = 2.0*mu*s*DSDGD/kappa;
234       DFDGD = kappa*DF1DGD/(f1*f1);
235       DFXDGD[IX][IS] = DS[IS]*Ex_unif[0]*DFDGD;
236     }
237   }
238   Fx = HALF*Fx/dt;
239 
240   /****************************************************
241                    Set output arguments
242   ****************************************************/
243 
244   Exc[0] = Fx;
245   Exc[1] = Fc;
246   for (IS=0; IS<=1; IS++){
247     DEXDD[IS] = DFXDD[IS];
248     DECDD[IS] = DFCDD[IS];
249     for (IX=0; IX<=2; IX++){
250       DEXDGD[IX][IS] = DFXDGD[IX][IS];
251       DECDGD[IX][IS] = DFCDGD[IX][IS];
252     }
253   }
254 }
255