1 /**********************************************************************
2 TRAN_Set_Electrode_Grid.c:
3
4 TRAN_Set_Electrode_Grid.c is a subroutine to calculate charge density
5 near the boundary region of the extended system, contributed from
6 the electrodes.
7
8 The contribution is added to the regions [0:TRAN_grid_bound[0]] and
9 [TRAN_grid_bound[1]:Ngrid1-1].
10
11 Log of TRAN_Set_Electrode_Grid.c:
12
13 24/July/2008 Released by T.Ozaki
14
15 ***********************************************************************/
16 /* revised by Y. Xiao for Noncollinear NEGF calculations */
17
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <math.h>
21 #include "openmx_common.h"
22 #include "mpi.h"
23 #include <fftw3.h>
24 #include "tran_variables.h"
25 #include "tran_prototypes.h"
26
27 #define print_stdout 0
28
29 double Interpolated_Func(double R, double *xg, double *v, int N);
30 double Dot_Product(double a[4], double b[4]);
31
32
33
34 /*
35 * input
36 * dDen_Grid_e
37 * dVHart_Grid_e
38 *
39 * output
40 * double **ElectrodeDensity_Grid
41 * double **VHart_Boundary
42 */
43
44
45 /*
46 * This is also used for Vtot_Grid
47 */
48
49
50 /* implicit input from trans_variables.h
51 * double **dDen_Grid_e
52 * double **dVHart_Grid_e
53 */
54
55
56
57
58
59
TRAN_Set_Electrode_Grid(MPI_Comm comm1,int * TRAN_Poisson_flag2,double * Grid_Origin,double tv[4][4],double Left_tv[4][4],double Right_tv[4][4],double gtv[4][4],int Ngrid1,int Ngrid2,int Ngrid3)60 void TRAN_Set_Electrode_Grid(MPI_Comm comm1,
61 int *TRAN_Poisson_flag2,
62 double *Grid_Origin, /* origin of the grid */
63 double tv[4][4], /* unit vector of the cell*/
64 double Left_tv[4][4], /* unit vector left */
65 double Right_tv[4][4], /* unit vector right */
66 double gtv[4][4], /* unit vector of the grid point, which is gtv*integer */
67 int Ngrid1,
68 int Ngrid2,
69 int Ngrid3 /* # of c grid points */
70 )
71 {
72 int l1[2];
73 int i,j,k,k2,k3;
74 int tnoA,tnoB,wanB,wanA;
75 int GA_AN,MA_AN,Nc,GNc,LB_AN_e;
76 int GB_AN;
77 int Rn_e,GA_AN_e,GB_AN_e,GRc;
78 int side,direction;
79 int spin,p2,p3;
80 int n1,n2,n3;
81 int id,gidx;
82 double rcutA,rcutB,r,r1,r2;
83 double dx,dy,dz,xx;
84 double dx1,dy1,dz1;
85 double dx2,dy2,dz2;
86 double x1,y1,z1;
87 double x2,y2,z2;
88 double sum,tmp;
89 double offset[4];
90 double R[4];
91 double xyz[4];
92 double *Chi0,*Chi1,*Chi2;
93 int idim=1;
94 int myid;
95 fftw_complex *in, *out;
96 fftw_plan p;
97
98 MPI_Comm_rank(comm1, &myid);
99
100 /* for passing TRAN_Poisson_flag to "DFT" */
101 *TRAN_Poisson_flag2 = TRAN_Poisson_flag;
102
103 if (myid==Host_ID){
104 printf("<TRAN_Set_Electrode_Grid>\n");
105 }
106
107 /* allocation of array */
108
109 Chi0 = (double*)malloc(sizeof(double)*List_YOUSO[7]);
110 Chi1 = (double*)malloc(sizeof(double)*List_YOUSO[7]);
111 Chi2 = (double*)malloc(sizeof(double)*List_YOUSO[7]);
112
113 in = fftw_malloc(sizeof(fftw_complex)*List_YOUSO[17]);
114 out = fftw_malloc(sizeof(fftw_complex)*List_YOUSO[17]);
115
116 /* allocation of arrays */
117 if (print_stdout){
118 printf("%d %d %d %d %d\n",Ngrid1,Ngrid2,Ngrid3,TRAN_grid_bound[0], TRAN_grid_bound[1]);
119 }
120
121 for (side=0; side<2; side++) {
122 ElectrodeDensity_Grid[side] = (double**)malloc(sizeof(double*)*(SpinP_switch_e[side]+1));
123 }
124
125 /* left lead */
126
127 side = 0;
128 l1[0] = 0;
129 l1[1] = TRAN_grid_bound[0];
130 for (spin=0; spin<=SpinP_switch_e[side]; spin++) {
131 ElectrodeDensity_Grid[side][spin] = (double*)malloc(sizeof(double)*Ngrid3*Ngrid2*(l1[1]-l1[0]+1));
132 }
133
134 ElectrodeADensity_Grid[side] = (double*)malloc(sizeof(double)*Ngrid3*Ngrid2*(l1[1]-l1[0]+1));
135 ElectrodedVHart_Grid[side] = (double*)malloc(sizeof(double)*Ngrid3*Ngrid2*(l1[1]-l1[0]+1));
136
137 VHart_Boundary[side] = (dcomplex***)malloc(sizeof(dcomplex**)*Ngrid1_e[side]);
138 for (n1=0; n1<Ngrid1_e[side]; n1++){
139 VHart_Boundary[side][n1] = (dcomplex**)malloc(sizeof(dcomplex*)*Ngrid2);
140 for (n2=0; n2<Ngrid2; n2++){
141 VHart_Boundary[side][n1][n2] = (dcomplex*)malloc(sizeof(dcomplex)*Ngrid3);
142 }
143 }
144
145 /* right lead */
146
147 side = 1;
148 l1[0] = TRAN_grid_bound[1];
149 l1[1] = Ngrid1 - 1;
150 for (spin=0; spin<=SpinP_switch_e[side]; spin++) {
151 ElectrodeDensity_Grid[side][spin] = (double*)malloc(sizeof(double)*Ngrid3*Ngrid2*(l1[1]-l1[0]+1));
152 }
153
154 ElectrodeADensity_Grid[side] = (double*)malloc(sizeof(double)*Ngrid3*Ngrid2*(l1[1]-l1[0]+1));
155 ElectrodedVHart_Grid[side] = (double*)malloc(sizeof(double)*Ngrid3*Ngrid2*(l1[1]-l1[0]+1));
156
157 VHart_Boundary[side] = (dcomplex***)malloc(sizeof(dcomplex**)*Ngrid1_e[side]);
158 for (n1=0; n1<Ngrid1_e[side]; n1++){
159 VHart_Boundary[side][n1] = (dcomplex**)malloc(sizeof(dcomplex*)*Ngrid2);
160 for (n2=0; n2<Ngrid2; n2++){
161 VHart_Boundary[side][n1][n2] = (dcomplex*)malloc(sizeof(dcomplex)*Ngrid3);
162 }
163 }
164
165 /*******************************************************
166 charge density contributed by the left and right sides
167 *******************************************************/
168
169 for (side=0; side<=1; side++){
170
171 if (side==0){
172 direction = -1;
173 l1[0] = 0;
174 l1[1] = TRAN_grid_bound[0];
175 }
176 else{
177 direction = 1;
178 l1[0] = TRAN_grid_bound[1];
179 l1[1] = Ngrid1-1;
180 }
181
182 /* initialize ElectrodeDensity_Grid and ElectrodeADensity_Grid */
183
184 for (spin=0; spin<=SpinP_switch_e[side]; spin++) {
185 for (i=0; i<Ngrid3*Ngrid2*(l1[1]-l1[0]+1); i++){
186 ElectrodeDensity_Grid[side][spin][i] = 0.0;
187 }
188 }
189 for (i=0; i<Ngrid3*Ngrid2*(l1[1]-l1[0]+1); i++){
190 ElectrodeADensity_Grid[side][i] = 0.0;
191 }
192
193 /* calculate charge density */
194
195 for (n1=l1[0]; n1<=l1[1]; n1++){
196 for (n2=0; n2<Ngrid2; n2++){
197 for (n3=0; n3<Ngrid3; n3++){
198
199 GNc = n1*Ngrid2*Ngrid3 + n2*Ngrid3 + n3;
200 Get_Grid_XYZ(GNc, xyz);
201
202 for (p2=-1; p2<=1; p2++){
203 for (p3=-1; p3<=1; p3++){
204 for (GA_AN_e=1; GA_AN_e<=atomnum_e[side]; GA_AN_e++){
205
206 if (side==0) GA_AN = GA_AN_e;
207 else GA_AN = GA_AN_e + Latomnum + Catomnum;
208
209 wanA = WhatSpecies[GA_AN];
210 tnoA = Spe_Total_CNO[wanA];
211 rcutA = Spe_Atom_Cut1[wanA];
212
213 x1 = Gxyz[GA_AN][1]
214 + (double)direction*tv_e[side][1][1]
215 + (double)p2*tv_e[side][2][1]
216 + (double)p3*tv_e[side][3][1];
217
218 y1 = Gxyz[GA_AN][2]
219 + (double)direction*tv_e[side][1][2]
220 + (double)p2*tv_e[side][2][2]
221 + (double)p3*tv_e[side][3][2];
222
223 z1 = Gxyz[GA_AN][3]
224 + (double)direction*tv_e[side][1][3]
225 + (double)p2*tv_e[side][2][3]
226 + (double)p3*tv_e[side][3][3];
227
228 dx1 = xyz[1] - x1;
229 dy1 = xyz[2] - y1;
230 dz1 = xyz[3] - z1;
231 r1 = sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
232
233 if (r1<=rcutA){
234
235 /******************************
236 ElectrodeDensity_Grid
237 ******************************/
238
239 Get_Orbitals(wanA,dx1,dy1,dz1,Chi1);
240
241 for (LB_AN_e=0; LB_AN_e<=FNAN_e[side][GA_AN_e]; LB_AN_e++){
242
243 GB_AN_e = natn_e[side][GA_AN_e][LB_AN_e];
244 Rn_e = ncn_e[side][GA_AN_e][LB_AN_e];
245
246 if (side==0) GB_AN = GB_AN_e;
247 else GB_AN = GB_AN_e + Latomnum + Catomnum;
248
249 wanB = WhatSpecies[GB_AN];
250 tnoB = Spe_Total_CNO[wanB];
251 rcutB = Spe_Atom_Cut1[wanB];
252
253 x2 = Gxyz[GB_AN][1]
254 + (double)(direction+atv_ijk_e[side][Rn_e][1])*tv_e[side][1][1]
255 + (double)(p2 +atv_ijk_e[side][Rn_e][2])*tv_e[side][2][1]
256 + (double)(p3 +atv_ijk_e[side][Rn_e][3])*tv_e[side][3][1];
257
258 y2 = Gxyz[GB_AN][2]
259 + (double)(direction+atv_ijk_e[side][Rn_e][1])*tv_e[side][1][2]
260 + (double)(p2 +atv_ijk_e[side][Rn_e][2])*tv_e[side][2][2]
261 + (double)(p3 +atv_ijk_e[side][Rn_e][3])*tv_e[side][3][2];
262
263 z2 = Gxyz[GB_AN][3]
264 + (double)(direction+atv_ijk_e[side][Rn_e][1])*tv_e[side][1][3]
265 + (double)(p2 +atv_ijk_e[side][Rn_e][2])*tv_e[side][2][3]
266 + (double)(p3 +atv_ijk_e[side][Rn_e][3])*tv_e[side][3][3];
267
268 dx2 = xyz[1] - x2;
269 dy2 = xyz[2] - y2;
270 dz2 = xyz[3] - z2;
271 r2 = sqrt(dx2*dx2 + dy2*dy2 + dz2*dz2);
272
273 if (r2<=rcutB){
274
275 Get_Orbitals(wanB,dx2,dy2,dz2,Chi2);
276
277 for (spin=0; spin<=SpinP_switch; spin++) {
278
279 if (atv_ijk_e[side][Rn_e][1]==(-direction)){
280
281 sum = 0.0;
282 for (i=0; i<tnoA; i++){
283 for (j=0; j<tnoB; j++){
284 /* revised by Y. Xiao for Noncollinear NEGF calculations */
285 /* sum += 2.0*DM_e[side][0][spin][GA_AN_e][LB_AN_e][i][j]*Chi1[i]*Chi2[j]; */
286 if(spin==3){
287 sum -= 2.0*DM_e[side][0][spin][GA_AN_e][LB_AN_e][i][j]*Chi1[i]*Chi2[j];
288 } else {
289 sum += 2.0*DM_e[side][0][spin][GA_AN_e][LB_AN_e][i][j]*Chi1[i]*Chi2[j];
290 }
291 /* until here by Y. Xiao for Noncollinear NEGF calculations */
292 }
293 }
294 }
295
296 else if (atv_ijk_e[side][Rn_e][1]==0){
297
298 sum = 0.0;
299 for (i=0; i<tnoA; i++){
300 for (j=0; j<tnoB; j++){
301 /* revised by Y. Xiao for Noncollinear NEGF calculations */
302 /* sum += DM_e[side][0][spin][GA_AN_e][LB_AN_e][i][j]*Chi1[i]*Chi2[j]; */
303 if(spin==3){
304 sum -= DM_e[side][0][spin][GA_AN_e][LB_AN_e][i][j]*Chi1[i]*Chi2[j];
305 } else {
306 sum += DM_e[side][0][spin][GA_AN_e][LB_AN_e][i][j]*Chi1[i]*Chi2[j];
307 }
308 /* until here by Y. Xiao for Noncollinear NEGF calculations */
309 }
310 }
311 }
312
313 gidx = (n1-l1[0])*Ngrid2*Ngrid3 + n2*Ngrid3 + n3;
314 ElectrodeDensity_Grid[side][spin][gidx] += sum;
315
316 }
317 }
318 }
319
320 /******************************
321 ElectrodeADensity_Grid
322 ******************************/
323
324 xx = 0.5*log(dx1*dx1 + dy1*dy1 + dz1*dz1);
325 gidx = (n1-l1[0])*Ngrid2*Ngrid3 + n2*Ngrid3 + n3;
326 ElectrodeADensity_Grid[side][gidx] += 0.5*KumoF( Spe_Num_Mesh_PAO[wanA], xx,
327 Spe_PAO_XV[wanA], Spe_PAO_RV[wanA],
328 Spe_Atomic_Den[wanA]);
329
330 } /* if (r1<=rcutA) */
331 }
332 }
333 }
334 }
335 }
336 }
337
338 } /* side */
339
340 /*******************************************************
341 2D FFT of dDen_Grid_e, which is difference between
342 charge density calculated by KS wave functions and
343 the superposition of atomic charge density, of the
344 left and right leads on the bc plane for TRAN_Poisson
345 *******************************************************/
346
347 for (side=0; side<=1; side++){
348
349 /* set VHart_Boundary in real space */
350
351 for (n1=0; n1<Ngrid1_e[side]; n1++){
352 for (n2=0; n2<Ngrid2; n2++){
353 for (n3=0; n3<Ngrid3; n3++){
354
355 /* borrow the array "VHart_Boundary" */
356 VHart_Boundary[side][n1][n2][n3].r = dDen_Grid_e[side][n1*Ngrid2*Ngrid3+n2*Ngrid3+n3];
357 VHart_Boundary[side][n1][n2][n3].i = 0.0;
358 }
359 }
360 }
361
362 /* FFT of VHart_Boundary for c-axis */
363
364 p = fftw_plan_dft_1d(Ngrid3,in,out,-1,FFTW_ESTIMATE);
365
366 for (n1=0; n1<Ngrid1_e[side]; n1++){
367 for (n2=0; n2<Ngrid2; n2++){
368
369 for (n3=0; n3<Ngrid3; n3++){
370
371 in[n3][0] = VHart_Boundary[side][n1][n2][n3].r;
372 in[n3][1] = VHart_Boundary[side][n1][n2][n3].i;
373 }
374
375 fftw_execute(p);
376
377 for (k3=0; k3<Ngrid3; k3++){
378
379 VHart_Boundary[side][n1][n2][k3].r = out[k3][0];
380 VHart_Boundary[side][n1][n2][k3].i = out[k3][1];
381 }
382 }
383 }
384
385 fftw_destroy_plan(p);
386
387 /* FFT of VHart_Boundary for b-axis */
388
389 p = fftw_plan_dft_1d(Ngrid2,in,out,-1,FFTW_ESTIMATE);
390
391 for (n1=0; n1<Ngrid1_e[side]; n1++){
392 for (k3=0; k3<Ngrid3; k3++){
393
394 for (n2=0; n2<Ngrid2; n2++){
395
396 in[n2][0] = VHart_Boundary[side][n1][n2][k3].r;
397 in[n2][1] = VHart_Boundary[side][n1][n2][k3].i;
398 }
399
400 fftw_execute(p);
401
402 for (k2=0; k2<Ngrid2; k2++){
403
404 VHart_Boundary[side][n1][k2][k3].r = out[k2][0];
405 VHart_Boundary[side][n1][k2][k3].i = out[k2][1];
406 }
407 }
408 }
409
410 fftw_destroy_plan(p);
411
412 tmp = 1.0/(double)(Ngrid2*Ngrid3);
413
414 for (n1=0; n1<Ngrid1_e[side]; n1++){
415 for (k2=0; k2<Ngrid2; k2++){
416 for (k3=0; k3<Ngrid3; k3++){
417 VHart_Boundary[side][n1][k2][k3].r *= tmp;
418 VHart_Boundary[side][n1][k2][k3].i *= tmp;
419 }
420 }
421 }
422
423 } /* side */
424
425 /*******************************************************
426 interpolation of 2D Fourier transformed difference
427 charge of the left and right leads on the bc plane for
428 TRAN_Poisson_FFT_Extended
429 *******************************************************/
430
431 {
432 int ip,Num;
433 double x;
434 double *vr,*vi,*xg;
435
436 if (1.0<fabs(tv[1][1])) ip = 1;
437 else if (1.0<fabs(tv[1][2])) ip = 2;
438 else if (1.0<fabs(tv[1][3])) ip = 3;
439
440 side = 0;
441 IntNgrid1_e[side] = (int)(fabs((double)TRAN_FFTE_CpyNum*tv_e[side][1][ip]+1.0e-8)/length_gtv[1]);
442
443 dDen_IntBoundary[side] = (dcomplex***)malloc(sizeof(dcomplex**)*IntNgrid1_e[side]);
444 for (n1=0; n1<IntNgrid1_e[side]; n1++){
445 dDen_IntBoundary[side][n1] = (dcomplex**)malloc(sizeof(dcomplex*)*Ngrid2);
446 for (n2=0; n2<Ngrid2; n2++){
447 dDen_IntBoundary[side][n1][n2] = (dcomplex*)malloc(sizeof(dcomplex)*Ngrid3);
448 }
449 }
450
451 side = 1;
452 IntNgrid1_e[side] = (int)(fabs((double)TRAN_FFTE_CpyNum*tv_e[side][1][ip]+1.0e-8)/length_gtv[1]);
453
454 dDen_IntBoundary[side] = (dcomplex***)malloc(sizeof(dcomplex**)*IntNgrid1_e[side]);
455 for (n1=0; n1<IntNgrid1_e[side]; n1++){
456 dDen_IntBoundary[side][n1] = (dcomplex**)malloc(sizeof(dcomplex*)*Ngrid2);
457 for (n2=0; n2<Ngrid2; n2++){
458 dDen_IntBoundary[side][n1][n2] = (dcomplex*)malloc(sizeof(dcomplex)*Ngrid3);
459 }
460 }
461
462 for (side=0; side<=1; side++){
463
464 vr = (double*)malloc(sizeof(double)*(TRAN_FFTE_CpyNum+2)*Ngrid1_e[side]);
465 vi = (double*)malloc(sizeof(double)*(TRAN_FFTE_CpyNum+2)*Ngrid1_e[side]);
466 xg = (double*)malloc(sizeof(double)*(TRAN_FFTE_CpyNum+2)*Ngrid1_e[side]);
467
468 for (k2=0; k2<Ngrid2; k2++){
469 for (k3=0; k3<Ngrid3; k3++){
470
471 /* set up the 1D data */
472
473 for (i=0; i<(TRAN_FFTE_CpyNum+2); i++){
474
475 for (n1=0; n1<Ngrid1_e[side]; n1++){
476 vr[n1+i*Ngrid1_e[side]] = VHart_Boundary[side][n1][k2][k3].r;
477 vi[n1+i*Ngrid1_e[side]] = VHart_Boundary[side][n1][k2][k3].i;
478 }
479
480 for (n1=0; n1<Ngrid1_e[side]; n1++){
481 xg[n1+i*Ngrid1_e[side]] =
482 (double)(n1+i*Ngrid1_e[side])*fabs(tv_e[side][1][ip]/(double)Ngrid1_e[side])
483 - (double)(Ngrid1_e[side]*(TRAN_FFTE_CpyNum+1))*fabs(tv_e[side][1][ip]/(double)Ngrid1_e[side])
484 + Grid_Origin[ip];
485 }
486 }
487
488 /*
489 if (myid==0 && side==0 && k2==1 && k3==0){
490 for (n1=0; n1<Ngrid1_e[side]*(TRAN_FFTE_CpyNum+2); n1++){
491 printf("A %15.12f %15.12f %15.12f\n",xg[n1],vr[n1],vi[n1]);
492 }
493 }
494 */
495
496 /* interpolation */
497
498 for (n1=0; n1<IntNgrid1_e[side]; n1++){
499
500 x = (double)n1*length_gtv[1] - (double)IntNgrid1_e[side]*length_gtv[1] + Grid_Origin[ip];
501
502 dDen_IntBoundary[side][n1][k2][k3].r = Interpolated_Func(x, xg, vr, (TRAN_FFTE_CpyNum+2)*Ngrid1_e[side]);
503 dDen_IntBoundary[side][n1][k2][k3].i = Interpolated_Func(x, xg, vi, (TRAN_FFTE_CpyNum+2)*Ngrid1_e[side]);
504
505
506 /*
507 if (myid==0 && side==0 && k2==1 && k3==0){
508 printf("B %15.12f %15.12f %15.12f\n",
509 x,dDen_IntBoundary[side][n1][k2][k3].r,
510 dDen_IntBoundary[side][n1][k2][k3].i);
511 }
512 */
513
514 }
515 }
516 }
517
518 free(xg);
519 free(vi);
520 free(vr);
521
522 }
523 }
524
525 /****************************************************
526 2D FFT of dVHartree of the left and right leads
527 on the bc plane for TRAN_Poisson
528 ****************************************************/
529
530 for (side=0; side<=1; side++){
531
532 /* set VHart_Boundary in real space */
533
534 for (n1=0; n1<Ngrid1_e[side]; n1++){
535 for (n2=0; n2<Ngrid2; n2++){
536 for (n3=0; n3<Ngrid3; n3++){
537
538 VHart_Boundary[side][n1][n2][n3].r = dVHart_Grid_e[side][n1*Ngrid2*Ngrid3+n2*Ngrid3+n3];
539 VHart_Boundary[side][n1][n2][n3].i = 0.0;
540 }
541 }
542 }
543
544 /* FFT of VHart_Boundary for c-axis */
545
546 p = fftw_plan_dft_1d(Ngrid3,in,out,-1,FFTW_ESTIMATE);
547
548 for (n1=0; n1<Ngrid1_e[side]; n1++){
549 for (n2=0; n2<Ngrid2; n2++){
550
551 for (n3=0; n3<Ngrid3; n3++){
552
553 in[n3][0] = VHart_Boundary[side][n1][n2][n3].r;
554 in[n3][1] = VHart_Boundary[side][n1][n2][n3].i;
555 }
556
557 fftw_execute(p);
558
559 for (k3=0; k3<Ngrid3; k3++){
560
561 VHart_Boundary[side][n1][n2][k3].r = out[k3][0];
562 VHart_Boundary[side][n1][n2][k3].i = out[k3][1];
563 }
564 }
565 }
566
567 fftw_destroy_plan(p);
568
569 /* FFT of VHart_Boundary for b-axis */
570
571 p = fftw_plan_dft_1d(Ngrid2,in,out,-1,FFTW_ESTIMATE);
572
573 for (n1=0; n1<Ngrid1_e[side]; n1++){
574 for (k3=0; k3<Ngrid3; k3++){
575
576 for (n2=0; n2<Ngrid2; n2++){
577
578 in[n2][0] = VHart_Boundary[side][n1][n2][k3].r;
579 in[n2][1] = VHart_Boundary[side][n1][n2][k3].i;
580 }
581
582 fftw_execute(p);
583
584 for (k2=0; k2<Ngrid2; k2++){
585
586 VHart_Boundary[side][n1][k2][k3].r = out[k2][0];
587 VHart_Boundary[side][n1][k2][k3].i = out[k2][1];
588 }
589 }
590 }
591
592 fftw_destroy_plan(p);
593
594 tmp = 1.0/(double)(Ngrid2*Ngrid3);
595
596 for (n1=0; n1<Ngrid1_e[side]; n1++){
597 for (k2=0; k2<Ngrid2; k2++){
598 for (k3=0; k3<Ngrid3; k3++){
599 VHart_Boundary[side][n1][k2][k3].r *= tmp;
600 VHart_Boundary[side][n1][k2][k3].i *= tmp;
601 }
602 }
603 }
604
605 } /* side */
606
607 /* freeing of arrays */
608 free(Chi0);
609 free(Chi1);
610 free(Chi2);
611
612 fftw_free(in);
613 fftw_free(out);
614 }
615
616
617
618
Interpolated_Func(double R,double * xg,double * v,int N)619 double Interpolated_Func(double R, double *xg, double *v, int N)
620 {
621 int mp_min,mp_max,m;
622 double h1,h2,h3,f1,f2,f3,f4;
623 double g1,g2,x1,x2,y1,y2,f;
624 double result;
625
626 mp_min = 0;
627 mp_max = N - 1;
628
629 if (R<xg[0]){
630 printf("Error #1 in Interpolated_Func\n");
631 exit(0);
632 }
633 else if (xg[N-1]<R){
634 printf("Error #2 in Interpolated_Func\n");
635 exit(0);
636 }
637 else{
638
639 do{
640 m = (mp_min + mp_max)/2;
641 if (xg[m]<R)
642 mp_min = m;
643 else
644 mp_max = m;
645 }
646 while((mp_max-mp_min)!=1);
647 m = mp_max;
648
649 if (m<2)
650 m = 2;
651 else if (N<=m)
652 m = N - 2;
653
654 /****************************************************
655 spline like interpolation
656 ****************************************************/
657
658 if (m==1){
659
660 h2 = xg[m] - xg[m-1];
661 h3 = xg[m+1] - xg[m];
662
663 f2 = v[m-1];
664 f3 = v[m];
665 f4 = v[m+1];
666
667 h1 = -(h2+h3);
668 f1 = f4;
669 }
670 else if (m==(N-1)){
671
672 h1 = xg[m-1] - xg[m-2];
673 h2 = xg[m] - xg[m-1];
674
675 f1 = v[m-2];
676 f2 = v[m-1];
677 f3 = v[m];
678
679 h3 = -(h1+h2);
680 f4 = f1;
681 }
682 else{
683 h1 = xg[m-1] - xg[m-2];
684 h2 = xg[m] - xg[m-1];
685 h3 = xg[m+1] - xg[m];
686
687 f1 = v[m-2];
688 f2 = v[m-1];
689 f3 = v[m];
690 f4 = v[m+1];
691 }
692
693 /****************************************************
694 calculate the value at R
695 ****************************************************/
696
697 g1 = ((f3-f2)*h1/h2 + (f2-f1)*h2/h1)/(h1+h2);
698 g2 = ((f4-f3)*h2/h3 + (f3-f2)*h3/h2)/(h2+h3);
699
700 x1 = R - xg[m-1];
701 x2 = R - xg[m];
702 y1 = x1/h2;
703 y2 = x2/h2;
704
705 f = y2*y2*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2)
706 + y1*y1*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1);
707 }
708 result = f;
709 return result;
710 }
711