1 #include <stdio.h>
2
3 /*****************************************************************
4 Improved analytical tetrahedron method
5 Idea: PRB 49, 16223 (1994), by Bloechl, Jepsen and Andersen
6 *******************************************************************/
7
8 /*
9 input e[0:3],n
10 output: e[0]<e[1]<e[2]<e[4]
11 */
OrderE0(double * e,int n)12 void OrderE0(double *e, int n)
13 {
14
15 int i,j;
16 double t;
17
18 for (i=0;i<n-1;i++) {
19 for (j=i;j<n;j++) {
20 if (e[j]<e[i]) {
21 t = e[j];
22 e[j]=e[i];
23 e[i]=t;
24
25 }
26 }
27 }
28 }
29
30
31
32 /*
33 input e[0:3],a[0:3],n
34 output: e[0]<e[1]<e[2]<e[4]
35 a[] is also interchanged corresponding to e[]:
36 */
OrderE(double * e,double * a,int n)37 void OrderE(double *e,double *a, int n)
38 {
39
40 int i,j;
41 double t;
42
43 for (i=0;i<n-1;i++) {
44 for (j=i;j<n;j++) {
45 if (e[j]<e[i]) {
46 t = e[j];
47 e[j]=e[i];
48 e[i]=t;
49
50 t = a[j];
51 a[j]=a[i];
52 a[i]=t;
53
54 }
55 }
56 }
57 }
58
59
60 /*
61 density of states
62 Appendix C of the paper
63
64 assume that et[] is orderd
65
66 input et[0:3], *e
67 output dos
68 */
69
ATM_Dos(double * et,double * e,double * dos)70 void ATM_Dos(double *et, double *e, double *dos)
71 {
72 double e21, e31, e32, e41, e42, e43;
73 double e1,e2,e4;
74
75 /* Parameter adjustments */
76 --et;
77 #if 0
78 printf("ATM: %lf %lf %lf %lf %lf\n",*e,et[1],et[2],et[3],et[4]);
79 #endif
80
81 e21 = et[2] - et[1];
82 e31 = et[3] - et[1];
83 e32 = et[3] - et[2];
84 e41 = et[4] - et[1];
85 e42 = et[4] - et[2];
86 e43 = et[4] - et[3];
87 if (*e < et[1]) {
88 *dos = 0.;
89 } else if (*e > et[1] && *e < et[2]) {
90 e1=(*e - et[1]);
91 *dos = 3.0 * e1*e1 / (e21 * e31 * e41);
92 } else if (*e > et[2] && *e < et[3]) {
93 e2=(*e - et[2]);
94 *dos = (e21 * 3. + e2 * 6. - (e31 + e42) * 3. * e2 * e2
95 / (e32 * e42)) / (e31 * e41);
96 } else if (*e > et[3] && *e < et[4]) {
97 e4= (et[4] - *e);
98 *dos = 3.0* e4*e4 / (e41 * e42 * e43);
99 } else if (*e > et[4]) {
100 *dos = 0.;
101 }
102 #if 0
103 printf("ATM_Dos: %lf %lf %lf %lf %lf->%lf\n",*e, et[1],et[2],et[3],et[4],*dos);
104 #endif
105
106 }
107
108
109 /*
110 Appendix B of the paper
111 In the Appendix B, integrated weight is written, differenciate it.
112 assume that et[] is ordered.
113
114 input et[0:3], at[0:3],*e
115 output spectrum
116
117 <at> = sum_i at[i] w[i]
118 w[i] is a function of et[i]
119 An integrated w[i] is written in the paper.
120
121 */
ATM_Spectrum(double * et,double * at,double * e,double * spectrum)122 void ATM_Spectrum(double *et,double *at, double *e, double *spectrum)
123 {
124 double a21, a31, a32, a41;
125 double e21, e31, e32, e41, e42, e43;
126 double a42, a43;
127 double dos;
128
129 /* Parameter adjustments , --- f2c technique */
130 --at;
131
132 ATM_Dos(et, e, &dos);
133 --et;
134 #if 0
135 printf("DOS-> %lf %lf\n",*e,dos);
136 printf("%lf %lf %lf %lf\n",et[1],et[2],et[3],et[4]);
137 #endif
138 e21 = et[2] - et[1];
139 e31 = et[3] - et[1];
140 e32 = et[3] - et[2];
141 e41 = et[4] - et[1];
142 e42 = et[4] - et[2];
143 e43 = et[4] - et[3];
144 a21 = at[2] - at[1];
145 a31 = at[3] - at[1];
146 a32 = at[3] - at[2];
147 a41 = at[4] - at[1];
148 a42 = at[4] - at[2];
149 a43 = at[4] - at[3];
150 if (*e < et[1]) {
151 *spectrum = 0.;
152 } else if (*e > et[1] && *e < et[2]) {
153 *spectrum = dos * (at[1] + (*e - et[1]) * .33333333333333331 * (a21 /
154 e21 + a31 / e31 + a41 / e41));
155 } else if (*e > et[2] && *e < et[3]) {
156 *spectrum = dos * (at[1] + (a21 + e21 * a31 / e31 + e21 * a41 / e41) *
157 .33333333333333331 * (et[3] - *e) / e32 + (at[4] - at[1] - (
158 e43 * a41 / e41 + e43 * a42 / e42 + a43) * .33333333333333331)
159 * (*e - et[2]) / e32);
160 } else if (*e > et[3] && *e < et[4]) {
161 *spectrum = dos * (at[4] + (*e - et[4]) * .33333333333333331 * (a41 /
162 e41 + a42 / e42 + a43 / e43));
163 } else if (*e > et[4]) {
164 *spectrum = 0.;
165 }
166
167 }
168
169