1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
|
/* Calculate XC pot */
/* Copyright (c) 2007-2019 MJ Rutter
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3
* of the Licence, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, see http://www.gnu.org/licenses/
*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include "c2xsf.h"
/* Calculate XC pot according to Perdew and Zunger
*
* PRB 23 5048 (1981) https://doi.org/10.1103/PhysRevB.23.5048
*/
double lda(double dens){
double r_ws,Ec,uc,ux;
double B1,B2,gamma,A,B,C,D;
double alpha,rscfcv;
if (dens<1e-30) return 0;
A= 0.0311;
B=-0.048;
C= 0.0020;
D=-0.0116;
B1=1.0529;
B2=0.3334; /* Castep has 0.630035 = 0.5291 * 0.3334 */
gamma=-0.1423; /* Castep has -3.872211 = 27.2116 * gamma */
alpha=pow(9*M_PI/4,1./3.);
rscfcv=3*alpha/(4*M_PI); /* = 0.4581653 */
/* Castep has 6.59747 = 27.211 * 0.5291 * 0.4581653 */
/* Find Wigner-Seitz radius */
r_ws=pow(3/(4*M_PI*dens),1./3.);
/* and convert to Bohr */
r_ws/=BOHR;
if (debug>2) fprintf(stderr,"Wigner-Seitz radius: %lf Bohr\n",r_ws);
ux=-4*rscfcv/(3*r_ws);
if (debug>1) fprintf(stderr,"Exchange potential: %lf V\n",ux*H_eV);
if (r_ws>=1.0){
Ec=gamma/(1+B1*sqrt(r_ws)+B2*r_ws);
uc=Ec*(1+(7*B1/6)*sqrt(r_ws)+(4*B2/3)*r_ws)/(1+B1*sqrt(r_ws)+B2*r_ws);
}
else
uc=A*log(r_ws)+(B-A/3)+(2*C*r_ws/3)*log(r_ws)+(2*D-C)*r_ws/3;
return (ux+uc)*H_eV;
}
#ifdef TEST
int debug;
int main(int argc, char **argv)
{
double dens,pot;
debug=3;
sscanf(argv[1],"%lf",&dens);
printf("Density: %lf eA^-3\n",dens);
pot=lda(dens);
printf("Potential: %lf V\n",pot);
return 0;
}
#endif
|
