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