1 #include <assert.h>
2 #include <math.h>
3 #include <fftw3.h>
4
5 #include "vdw_radial.h"
6 #include "vdwxc.h"
7 #include "vdw_kernel.h"
8
9 #if defined HAVE_CONFIG_H
10 #include "config.h"
11 #endif
12
13 #include "vdw_core.h"
14
15 #if defined VDW_TIMING
16 #include "vdw_timing.h"
17 #endif
18
19 #define true 1
20 #define false 0
21
22 static const int SPLINE_ALPHAS = 0;
23
vdwxc_write_data(int N,double dr,double * data,char * file)24 void vdwxc_write_data(int N, double dr, double* data, char* file)
25 {
26 FILE *fp;
27 fp = fopen(file, "w");
28 int g;
29 for (g=0; g<N; g++) {
30 fprintf(fp, "%.10f %.10f\n", g*dr, data[g]);
31 }
32 fclose(fp);
33 }
34
vdwxc_fourier_bessel_transform(int Ninput,int Noutput,int lda_input,int lda_output,double dr,int premul_input,double * input,double * output,int direction)35 void vdwxc_fourier_bessel_transform(int Ninput, int Noutput,
36 int lda_input, int lda_output,
37 double dr,
38 int premul_input,
39 double* input, double* output,
40 int direction)
41 {
42 assert (direction == VDW_FBT_FORWARD || direction == VDW_FBT_BACKWARD);
43
44 int k;
45 int g;
46
47 // For now, just a reference implementation without FFTW.
48 if (direction == VDW_FBT_FORWARD) {
49 vdwxc_write_data(Ninput, dr, input, "fbtin.txt");
50 for (k=0; k<Noutput; k++) {
51 output[k*lda_output] = 0.0;
52 for (g=0; g<Ninput; g++) {
53 if (k == 0) {
54 // What to do to avoid division by zero?
55 //if (premul_input) {
56 // output[k*lda_output] += 2.0*VDW_PI*g*1.0/Noutput*input[g*lda_input];
57 //} else {
58 output[k*lda_output] += 2.0*VDW_PI*g*g*1.0/Noutput*input[g*lda_input]*dr*dr;
59 //}
60 } else {
61 //if (premul_input) {
62 // output[k*lda_output] += 2*sin(VDW_PI*g*k*1.0/(Noutput))/k*input[g*lda_input];
63 //} else {
64 output[k*lda_output] += 2*g*sin(VDW_PI*g*k*1.0/(Noutput))/k*input[g*lda_input]*dr;
65 //}
66 }
67 }
68 output[k*lda_output] *= 4*VDW_PI;
69 }
70 vdwxc_write_data(Noutput, dr, output, "fbtout.txt");
71 }
72 if (direction == VDW_FBT_BACKWARD)
73 {
74 vdwxc_write_data(Noutput, dr, output, "ifbtin.txt");
75 double dk = 2*VDW_PI/(Noutput*dr);
76 for (g=0; g<Ninput; g++) {
77 input[g*lda_input] = 0.0;
78 for (k=0; k<Noutput; k++) {
79 if (g == 0) {
80 // What to do to avoid division by zero?
81 input[g*lda_input] += (4*VDW_PI)*VDW_PI*k*k*dk*1.0/Noutput*output[k*lda_output];
82 //double c = k*VDW_PI*k*1.0/Noutput*output[k];
83 //printf("contribution %.10f \n", c);
84 } else {
85 input[g*lda_input] += (4*VDW_PI)*k*dk*sin(VDW_PI*g*k*1.0/(Noutput))/g*output[k*lda_output];
86 }
87
88 }
89 input[g*lda_input] /= (2*VDW_PI)*(2*VDW_PI)*(2*VDW_PI); //((2*VDW_PI)*(2*VDW_PI)*(2*VDW_PI));
90 }
91 vdwxc_write_data(Ninput, dr, input, "ifbtout.txt");
92 }
93 }
94
95
96 // Debug function, to test the spherical Bessel transform
97 // OBSOLETE!
vdwxc_calculate_radial_hartree(int N,double dr,double * rho_i,double * dedn_i)98 double vdwxc_calculate_radial_hartree(int N, double dr, double* rho_i, double* dedn_i)
99 {
100 double* rho_k = fftw_malloc(2 * sizeof(double) * N);
101 double* kernel_i = fftw_malloc(2 * sizeof(double) * N);
102 double* kernel_k = fftw_malloc(2 * sizeof(double) * N);
103 int i;
104 for (i=0; i<2*N; i++) {
105 kernel_i[i] = 1.0;
106 }
107 vdwxc_fourier_bessel_transform(N, 2*N, 1, 1, dr, false, rho_i, rho_k, VDW_FBT_FORWARD);
108 vdwxc_fourier_bessel_transform(2*N, 2*N, 1, 1, dr, true, kernel_i, kernel_k, VDW_FBT_FORWARD);
109 for (i=0; i<2*N; i++) {
110 rho_k[i] *= kernel_k[i];
111 }
112
113 vdwxc_fourier_bessel_transform(N, 2*N, 1, 1, dr, false, rho_i, rho_k, VDW_FBT_BACKWARD);
114 fftw_free(rho_k);
115 fftw_free(kernel_i);
116 fftw_free(kernel_k);
117 return 0.0;
118 }
119
120
121
vdwxc_calculate_radial(vdwxc_data data,int N,double dr,double * rho_i,double * sigma_i,double * dedn_i,double * dedsigma_i)122 double vdwxc_calculate_radial(vdwxc_data data, int N, double dr, double* rho_i,
123 double* sigma_i, double* dedn_i,
124 double* dedsigma_i)
125 {
126 START_TIMER(radial, "vdw_df_calculate_radial");
127 int pad = 2;
128 assert(0);
129
130 printf("N: %d\n", N);
131 double* radial_q0_g = fftw_malloc(sizeof(double) * N);
132 double* radial_dq0_drho_g = fftw_malloc(sizeof(double) * N);
133 double* radial_dq0_dsigma_g = fftw_malloc(sizeof(double) * N);
134
135 double* radial_theta_ga = fftw_malloc(sizeof(double) * SPLINE_ALPHAS * N);
136 double* radial_F_ga = fftw_malloc(sizeof(double) * SPLINE_ALPHAS * N);
137
138 double* radial_theta_ka = fftw_malloc(pad * sizeof(double) * SPLINE_ALPHAS * N);
139 double* radial_F_ka = fftw_malloc(pad * sizeof(double) * SPLINE_ALPHAS * N);
140
141 printf("Allocated buffers\n");
142
143 vdwxc_write_data(N, dr, rho_i, "rho_input.dat");
144 printf("Wrote data\n");
145
146 START_TIMER(q0, "q0");
147 printf("Started q0 timer\n");
148 //printf("%d\n",data);
149 //vdwxc_calculate_q0(N, data->Z_ab, data->q_cut,
150 // rho_i, sigma_i, radial_q0_g,
151 // radial_dq0_drho_g, radial_dq0_dsigma_g);
152 END_TIMER(q0, "q0");
153
154 START_TIMER(thetas, "thetas");
155 //vdwxc_calculate_thetas(1, 1, N, 1,
156 // rho_i, radial_q0_g,
157 // radial_theta_ga);
158 END_TIMER(thetas, "thetas");
159
160 int a1;
161 int a2;
162 for (a1=0; a1<SPLINE_ALPHAS; a1++) {
163 vdwxc_fourier_bessel_transform(N, pad*N,
164 SPLINE_ALPHAS, SPLINE_ALPHAS,
165 dr, false,
166 radial_theta_ga+a1, radial_theta_ka+a1,
167 VDW_FBT_FORWARD);
168
169 }
170
171 // Convolution
172 double kernel_aa[SPLINE_ALPHAS*SPLINE_ALPHAS];
173 double energy = 0.0;
174 complex double F_a[SPLINE_ALPHAS];
175
176 int k;
177 int g;
178 double L = dr * N;
179 double dk = VDW_PI / (L * pad);
180
181 for (k=0; k<pad*N; k++) {
182 double kk = k * dk;
183 double weight = 4*VDW_PI*kk*kk*dk;
184 //vdwxc_interpolate_kernels(kk, kernel_aa);
185 assert(0);
186 double* ikernel_aa = kernel_aa;
187 for (a1=0; a1 < SPLINE_ALPHAS; a1++) {
188 double F = 0.0;
189 for (a2=0; a2 < SPLINE_ALPHAS; a2++) {
190 // This array should be completely contiguous and needs no fancy indexing.
191 double a2val = radial_theta_ka[k * SPLINE_ALPHAS + a2];
192 double kernel = *ikernel_aa++;
193 double dF = a2val * kernel;
194 F += dF;
195 }
196 //assert(kindex < data->icell.Nlocal[0]*data->icell.Nlocal[1]*data->icell.Nlocal[2]);
197 double a1val = radial_theta_ka[k * SPLINE_ALPHAS + a1];
198 energy += a1val * F * weight;
199 F_a[a1] = F;
200 }
201 for(a1=0; a1 < SPLINE_ALPHAS; a1++) {
202 radial_F_ka[k * SPLINE_ALPHAS + a1] = F_a[a1];
203 }
204 }
205
206 // Todo: Figure out what is going on in the prefactor.
207 energy *= 0.5 * (L)*(L)*dr*dr/(4*VDW_PI)/(4*VDW_PI)/pow(2*VDW_PI,3.0/2);
208 /*for(a1=0; a1 < SPLINE_ALPHAS; a1++) {
209 for (k=0; k<data->rgd.Nr*2; k++) {
210 assert (radial_theta_ka[k * SPLINE_ALPHAS + a1] == 0);
211 }
212 }*/
213
214 for (a1=0; a1<SPLINE_ALPHAS; a1++) {
215 vdwxc_fourier_bessel_transform(N, pad*N,
216 SPLINE_ALPHAS, SPLINE_ALPHAS,
217 dr, false,
218 radial_F_ga+a1, radial_F_ka+a1,
219 VDW_FBT_BACKWARD);
220 }
221
222 /*
223 for(a1=0; a1 < SPLINE_ALPHAS; a1++) {
224 for (k=0; k<data->rgd.Nr; k++) {
225 //printf("%d %d %.15f\n",a1,k, radial_theta_ga[k * SPLINE_ALPHAS + a1]);
226 assert (radial_theta_ga[k * SPLINE_ALPHAS + a1] == 0);
227 }
228 }*/
229
230
231 int alpha;
232
233 double p_a[SPLINE_ALPHAS];
234 double dpdq_a[SPLINE_ALPHAS];
235 double prefactor = 0.5 / pad;
236 for (g=0; g<N; g++) {
237 double rho_dq0_drho = rho_i[g] * radial_dq0_drho_g[g];
238 double rho_dq0_dsigma = rho_i[g] * radial_dq0_dsigma_g[g];
239 dedn_i[g] = 0.0;
240 dedsigma_i[g] = 0.0; // zero or not?
241 //vdw_evaluate_palpha_splines(1.0, radial_q0_g[g], p_a);
242 //vdw_evaluate_palpha_splines_derivative(radial_q0_g[g], dpdq_a);
243 assert(0);
244 for(alpha=0; alpha<SPLINE_ALPHAS; alpha++) {
245 dedn_i[g] += prefactor * radial_F_ga[alpha+g*SPLINE_ALPHAS] *
246 (p_a[alpha] + dpdq_a[alpha] * rho_dq0_drho);
247 dedsigma_i[g] += prefactor * radial_F_ga[alpha+g*SPLINE_ALPHAS] *
248 dpdq_a[alpha] * rho_dq0_dsigma;
249 }
250 }
251
252
253 fftw_free(radial_q0_g);
254 fftw_free(radial_dq0_drho_g);
255 fftw_free(radial_dq0_dsigma_g);
256 fftw_free(radial_theta_ga);
257 fftw_free(radial_theta_ka);
258 fftw_free(radial_F_ka);
259 fftw_free(radial_F_ga);
260 END_TIMER(radial, "vdw_df_calculate_radial");
261
262 printf("vdw-energy %.10f \n", energy);
263 return energy;
264
265 // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
266 // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
267 // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
268 // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
269 // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
270 /* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXxx
271 data->radial_plan_r2r =
272 fftw_plan_many_dft_r2r(1, // 1D FFT
273 &data->rgd.Nr, // Table of 1 number, number of grid points
274 SPLINE_ALPHAS, // how many FFTs
275 radial_theta_ga, // Input array
276 NULL, // full array (no subarray)
277 SPLINE_ALPHAS, // Concecutive grid points are within 20 doubles
278 1, // Individual transforms are next to eachother in memory
279 radial_theta_ka, // Output array
280 NULL, // full array (no subarray)
281 SPLINE_ALPHAS,
282 1,
283 FFTW_RODFT10,
284 FFTW_ESTIMATE);
285
286 */
287
288
289 //START_TIMER(potential, "potential");
290 //END_TIMER(potential, "potential");
291
292
293 // store integrals somehow for checking/debugging?
294 return energy;
295 }
296
297
298