1 /**********************************************************************
2 FT_ProExpn_VNA.c:
3
4 FT_ProExpn_VNA.c is a subroutine to Fourier transform
5 VNA separable projectors.
6
7 Log of FT_ProExpn_VNA.c:
8
9 7/Apr/2004 Released by T.Ozaki
10
11 ***********************************************************************/
12
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <math.h>
16 #include <time.h>
17 /* stat section */
18 #include <sys/types.h>
19 #include <sys/stat.h>
20 #include <unistd.h>
21 /* end stat section */
22 #include "openmx_common.h"
23 #include "mpi.h"
24 #include <omp.h>
25
26
27
FT_ProExpn_VNA()28 void FT_ProExpn_VNA()
29 {
30 int numprocs,myid,ID,tag=999;
31 int count,NumSpe;
32 int L,i,kj;
33 int Lspe,spe,GL,Mul;
34 int RestartRead_Succeed;
35 double Sr,Dr;
36 double norm_k,h,dum0;
37 double rmin,rmax,r,sum;
38 double kmin,kmax,Sk,Dk;
39 double RGL[GL_Mesh + 2];
40 double *SumTmp;
41 double tmp0,tmp1;
42 double **SphB;
43 double *tmp_SphB,*tmp_SphBp;
44 double TStime, TEtime;
45 /* for MPI */
46 MPI_Status stat;
47 MPI_Request request;
48 /* for OpenMP */
49 int OMPID,Nthrds,Nprocs;
50
51 char fileFT[YOUSO10];
52 char operate[300];
53 FILE *fp;
54 size_t size;
55
56 dtime(&TStime);
57
58 /* MPI */
59 MPI_Comm_size(mpi_comm_level1,&numprocs);
60 MPI_Comm_rank(mpi_comm_level1,&myid);
61
62 if (myid==Host_ID && 0<level_stdout) printf("<FT_ProExpn_VNA> Fourier transform of VNA separable projectors\n");
63
64 RestartRead_Succeed = 0;
65
66 /***********************************************************
67 In case of Scf_RestartFromFile==1, read Spe_VNA_Bessel
68 ***********************************************************/
69
70 if (Scf_RestartFromFile==1){
71
72 /****************************************************
73 regenerate radial grids in the k-space
74 for the MPI calculation
75 ****************************************************/
76
77 for (kj=0; kj<GL_Mesh; kj++){
78 kmin = Radial_kmin;
79 kmax = PAO_Nkmax;
80 Sk = kmax + kmin;
81 Dk = kmax - kmin;
82 norm_k = 0.50*(Dk*GL_Abscissae[kj] + Sk);
83 GL_NormK[kj] = norm_k;
84 }
85
86 /***********************************************************
87 read Spe_VNA_Bessel
88 ***********************************************************/
89
90 sprintf(fileFT,"%s%s_rst/%s.ftPEvna",filepath,filename,filename);
91
92 if ((fp = fopen(fileFT,"rb")) != NULL){
93
94 RestartRead_Succeed = 1;
95
96 for (spe=0; spe<SpeciesNum; spe++){
97 for (L=0; L<=List_YOUSO[35]; L++){
98 for (Mul=0; Mul<List_YOUSO[34]; Mul++){
99
100 size = fread(&Spe_VNA_Bessel[spe][L][Mul][0],sizeof(double),GL_Mesh,fp);
101 if (size!=GL_Mesh) RestartRead_Succeed = 0;
102 }
103 }
104 }
105
106 fclose(fp);
107 }
108 else{
109 printf("Could not open a file %s in FT_ProExpn_VNA\n",fileFT);
110 }
111 }
112
113 /***********************************************************
114 if (RestartRead_Succeed==0), calculate Spe_VNA_Bessel
115 ***********************************************************/
116
117 if (RestartRead_Succeed==0){
118
119 for (Lspe=0; Lspe<MSpeciesNum; Lspe++){
120
121 spe = Species_Top[myid] + Lspe;
122
123 /* initalize */
124 /* tabulation on Gauss-Legendre radial grid */
125
126 rmin = Spe_VPS_RV[spe][0];
127 rmax = Spe_Atom_Cut1[spe] + 0.5;
128 Sr = rmax + rmin;
129 Dr = rmax - rmin;
130 for (i=0; i<GL_Mesh; i++){
131 RGL[i] = 0.50*(Dr*GL_Abscissae[i] + Sr);
132 }
133
134 kmin = Radial_kmin;
135 kmax = PAO_Nkmax;
136 Sk = kmax + kmin;
137 Dk = kmax - kmin;
138
139 /* loop for kj */
140
141 #pragma omp parallel shared(spe,List_YOUSO,GL_Weight,GL_Abscissae,Dr,Dk,Sk,RGL,Projector_VNA,Spe_VPS_RV,Spe_Num_Mesh_VPS,Spe_VNA_Bessel) private(SumTmp,SphB,tmp_SphB,tmp_SphBp,OMPID,Nthrds,Nprocs,kj,norm_k,i,r,L,Mul,tmp0,dum0)
142
143 {
144
145 /* allocate arrays */
146
147 SumTmp = (double*)malloc(sizeof(double)*List_YOUSO[34]);
148
149 SphB = (double**)malloc(sizeof(double*)*(List_YOUSO[35]+3));
150 for(L=0; L<(List_YOUSO[35]+3); L++){
151 SphB[L] = (double*)malloc(sizeof(double)*GL_Mesh);
152 }
153
154 tmp_SphB = (double*)malloc(sizeof(double)*(List_YOUSO[35]+3));
155 tmp_SphBp = (double*)malloc(sizeof(double)*(List_YOUSO[35]+3));
156
157 /* get info. on OpenMP */
158
159 OMPID = omp_get_thread_num();
160 Nthrds = omp_get_num_threads();
161 Nprocs = omp_get_num_procs();
162
163 for ( kj=OMPID; kj<GL_Mesh; kj+=Nthrds ){
164
165 norm_k = 0.50*(Dk*GL_Abscissae[kj] + Sk);
166
167 /* calculate SphB */
168
169 for (i=0; i<GL_Mesh; i++){
170
171 r = RGL[i];
172 Spherical_Bessel(norm_k*r,List_YOUSO[35],tmp_SphB,tmp_SphBp);
173
174 for(L=0; L<=List_YOUSO[35]; L++){
175 SphB[L][i] = tmp_SphB[L];
176 }
177 }
178
179 /* loop for L */
180
181 for (L=0; L<=List_YOUSO[35]; L++){
182
183 /****************************************************
184 \int jL(k*r)RL r^2 dr
185 ****************************************************/
186
187 for (Mul=0; Mul<List_YOUSO[34]; Mul++) SumTmp[Mul] = 0.0;
188
189 /* Gauss-Legendre quadrature */
190
191 for (i=0; i<GL_Mesh; i++){
192
193 r = RGL[i];
194
195 tmp0 = r*r*GL_Weight[i]*SphB[L][i];
196 for (Mul=0; Mul<List_YOUSO[34]; Mul++){
197 dum0 = PhiF(r, Projector_VNA[spe][L][Mul], Spe_VPS_RV[spe], Spe_Num_Mesh_VPS[spe]);
198 SumTmp[Mul] += dum0*tmp0;
199 }
200 }
201
202 for (Mul=0; Mul<List_YOUSO[34]; Mul++){
203 Spe_VNA_Bessel[spe][L][Mul][kj] = 0.5*Dr*SumTmp[Mul];
204 }
205
206 } /* L */
207
208 } /* kj */
209
210 /* free arrays */
211
212 free(SumTmp);
213
214 for(L=0; L<(List_YOUSO[35]+3); L++){
215 free(SphB[L]);
216 }
217 free(SphB);
218
219 free(tmp_SphB);
220 free(tmp_SphBp);
221
222 #pragma omp flush(Spe_VNA_Bessel)
223
224 } /* #pragma omp parallel */
225 } /* Lspe */
226
227 /****************************************************
228 regenerate radial grids in the k-space
229 for the MPI calculation
230 ****************************************************/
231
232 for (kj=0; kj<GL_Mesh; kj++){
233 kmin = Radial_kmin;
234 kmax = PAO_Nkmax;
235 Sk = kmax + kmin;
236 Dk = kmax - kmin;
237 norm_k = 0.50*(Dk*GL_Abscissae[kj] + Sk);
238 GL_NormK[kj] = norm_k;
239 }
240
241 /***********************************************************
242 sending and receiving of Spe_VNA_Bessel by MPI
243 ***********************************************************/
244
245 for (ID=0; ID<Num_Procs2; ID++){
246 NumSpe = Species_End[ID] - Species_Top[ID] + 1;
247 for (Lspe=0; Lspe<NumSpe; Lspe++){
248 spe = Species_Top[ID] + Lspe;
249 for (L=0; L<=List_YOUSO[35]; L++){
250 for (Mul=0; Mul<List_YOUSO[34]; Mul++){
251 MPI_Bcast(&Spe_VNA_Bessel[spe][L][Mul][0],
252 GL_Mesh,MPI_DOUBLE,ID,mpi_comm_level1);
253 }
254 }
255 }
256 }
257
258 /***********************************************************
259 save Spe_VNA_Bessel
260 ***********************************************************/
261
262 if (myid==Host_ID){
263
264 sprintf(fileFT,"%s%s_rst/%s.ftPEvna",filepath,filename,filename);
265
266 if ((fp = fopen(fileFT,"wb")) != NULL){
267
268 for (spe=0; spe<SpeciesNum; spe++){
269 for (L=0; L<=List_YOUSO[35]; L++){
270 for (Mul=0; Mul<List_YOUSO[34]; Mul++){
271 fwrite(&Spe_VNA_Bessel[spe][L][Mul][0],sizeof(double),GL_Mesh,fp);
272 }
273 }
274 }
275
276 fclose(fp);
277 }
278 else{
279 printf("Could not open a file %s in FT_ProExpn_VNA\n",fileFT);
280 }
281 }
282
283 } /* if (RestartRead_Succeed==0) */
284
285 /***********************************************************
286 elapsed time
287 ***********************************************************/
288
289 dtime(&TEtime);
290
291 /*
292 printf("myid=%2d Elapsed Time (s) = %15.12f\n",myid,TEtime-TStime);
293 MPI_Finalize();
294 exit(0);
295 */
296
297 }
298
299
300
301
302