1 /**********************************************************************
2 FT_PAO.c:
3
4 FT_PAO.c is a subroutine to Fourier transform pseudo atomic
5 orbitals.
6
7 Log of FT_PAO.c:
8
9 15/Sep/2002 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
28
FT_PAO()29 void FT_PAO()
30 {
31 int numprocs,myid,ID,tag=999;
32 int count,NumSpe;
33 int i,kj,num_k;
34 int Lspe,spe,GL,Mul;
35 int RestartRead_Succeed;
36 double dk,norm_k,h;
37 double rmin,rmax,r,r2,sum;
38 double sy,sjp,syp;
39 double Sr,Dr,dum0;
40 double **SphB;
41 double *tmp_SphB,*tmp_SphBp;
42 double TStime, TEtime;
43 /* for MPI */
44 MPI_Status stat;
45 MPI_Request request;
46 /* for OpenMP */
47 int OMPID,Nthrds,Nprocs;
48
49 char fileFT[YOUSO10];
50 char operate[300];
51 FILE *fp;
52 size_t size;
53
54 dtime(&TStime);
55
56 /* MPI */
57 MPI_Comm_size(mpi_comm_level1,&numprocs);
58 MPI_Comm_rank(mpi_comm_level1,&myid);
59
60 if (myid==Host_ID && 0<level_stdout) printf("<FT_PAO> Fourier transform of pseudo atomic orbitals\n");
61
62 RestartRead_Succeed = 0;
63
64 /***********************************************************
65 In case of Scf_RestartFromFile==1, read Spe_RF_Bessel
66 ***********************************************************/
67
68 if (Scf_RestartFromFile==1){
69
70 /****************************************************
71 generate radial grids in the k-space
72 ****************************************************/
73
74 dk = PAO_Nkmax/(double)Ngrid_NormK;
75 for (i=0; i<Ngrid_NormK; i++){
76 NormK[i] = (double)i*dk;
77 }
78
79 /***********************************************************
80 read Spe_RF_Bessel
81 ***********************************************************/
82
83 sprintf(fileFT,"%s%s_rst/%s.ftpao",filepath,filename,filename);
84
85 if ((fp = fopen(fileFT,"rb")) != NULL){
86
87 RestartRead_Succeed = 1;
88
89 for (spe=0; spe<SpeciesNum; spe++){
90 for (GL=0; GL<=Spe_MaxL_Basis[spe]; GL++){
91 for (Mul=0; Mul<Spe_Num_Basis[spe][GL]; Mul++){
92
93 size = fread(&Spe_RF_Bessel[spe][GL][Mul][0],sizeof(double),List_YOUSO[15],fp);
94 if (size!=List_YOUSO[15]) RestartRead_Succeed = 0;
95 }
96 }
97 }
98
99 fclose(fp);
100 }
101 else{
102 printf("Could not open a file %s in FT_PAO\n",fileFT);
103 }
104 }
105
106 /***********************************************************
107 if (RestartRead_Succeed==0), calculate Spe_RF_Bessel
108 ***********************************************************/
109
110 if (RestartRead_Succeed==0){
111
112 for (Lspe=0; Lspe<MSpeciesNum; Lspe++){
113
114 spe = Species_Top[myid] + Lspe;
115
116 num_k = Ngrid_NormK;
117 dk = PAO_Nkmax/(double)num_k;
118 rmin = Spe_PAO_RV[spe][0];
119 rmax = Spe_Atom_Cut1[spe] + 0.5;
120 h = (rmax - rmin)/(double)OneD_Grid;
121
122 /* kj loop */
123
124 #pragma omp parallel shared(num_k,dk,spe,rmin,rmax,h,Spe_PAO_RV,Spe_Atom_Cut1,OneD_Grid,Spe_MaxL_Basis,Spe_Num_Basis,Spe_RF_Bessel) private(OMPID,Nthrds,Nprocs,kj,norm_k,i,r,r2,tmp_SphB,tmp_SphBp,GL,SphB,Mul,sum)
125 {
126
127 /* allocate SphB */
128
129 SphB = (double**)malloc(sizeof(double*)*(Spe_MaxL_Basis[spe]+3));
130 for(GL=0; GL<(Spe_MaxL_Basis[spe]+3); GL++){
131 SphB[GL] = (double*)malloc(sizeof(double)*(OneD_Grid+1));
132 }
133
134 tmp_SphB = (double*)malloc(sizeof(double)*(Spe_MaxL_Basis[spe]+3));
135 tmp_SphBp = (double*)malloc(sizeof(double)*(Spe_MaxL_Basis[spe]+3));
136
137 /* get info. on OpenMP */
138
139 OMPID = omp_get_thread_num();
140 Nthrds = omp_get_num_threads();
141 Nprocs = omp_get_num_procs();
142
143 for ( kj=OMPID; kj<num_k; kj+=Nthrds ){
144
145 norm_k = (double)kj*dk;
146
147 /* calculate SphB */
148
149 for (i=0; i<=OneD_Grid; i++){
150
151 r = rmin + (double)i*h;
152
153 Spherical_Bessel(norm_k*r,Spe_MaxL_Basis[spe],tmp_SphB,tmp_SphBp);
154
155 r2 = r*r;
156 for(GL=0; GL<=Spe_MaxL_Basis[spe]; GL++){
157 SphB[GL][i] = tmp_SphB[GL]*r2;
158 }
159 }
160
161 for(GL=0; GL<=Spe_MaxL_Basis[spe]; GL++){
162 SphB[GL][0] *= 0.5;
163 SphB[GL][OneD_Grid] *= 0.5;
164 }
165
166 /* loof for GL and Mul */
167
168 for (GL=0; GL<=Spe_MaxL_Basis[spe]; GL++){
169 for (Mul=0; Mul<Spe_Num_Basis[spe][GL]; Mul++){
170
171 /****************************************************
172 \int jL(k*r)RL r^2 dr
173 ****************************************************/
174
175 /* trapezoidal rule */
176
177 sum = 0.0;
178
179 for (i=0; i<=OneD_Grid; i++){
180 r = rmin + (double)i*h;
181 sum += RadialF(spe,GL,Mul,r)*SphB[GL][i];
182 }
183 sum = sum*h;
184
185 Spe_RF_Bessel[spe][GL][Mul][kj] = sum;
186
187 } /* Mul */
188 } /* GL */
189 } /* kj */
190
191 /* free SphB */
192
193 for(GL=0; GL<(Spe_MaxL_Basis[spe]+3); GL++){
194 free(SphB[GL]);
195 }
196 free(SphB);
197
198 free(tmp_SphB);
199 free(tmp_SphBp);
200
201 } /* #pragma omp parallel */
202
203 } /* Lspe */
204
205 /****************************************************
206 generate radial grids in the k-space
207 ****************************************************/
208
209 dk = PAO_Nkmax/(double)Ngrid_NormK;
210 for (i=0; i<Ngrid_NormK; i++){
211 NormK[i] = (double)i*dk;
212 }
213
214 /***********************************************************
215 sending and receiving of Spe_RF_Bessel by MPI
216 ***********************************************************/
217
218 for (ID=0; ID<Num_Procs2; ID++){
219 NumSpe = Species_End[ID] - Species_Top[ID] + 1;
220 for (Lspe=0; Lspe<NumSpe; Lspe++){
221 spe = Species_Top[ID] + Lspe;
222 for (GL=0; GL<=Spe_MaxL_Basis[spe]; GL++){
223 for (Mul=0; Mul<Spe_Num_Basis[spe][GL]; Mul++){
224
225 MPI_Bcast(&Spe_RF_Bessel[spe][GL][Mul][0],
226 List_YOUSO[15],MPI_DOUBLE,ID, mpi_comm_level1);
227 MPI_Barrier(mpi_comm_level1);
228 }
229 }
230 }
231 }
232
233 /***********************************************************
234 save Spe_RF_Bessel
235 ***********************************************************/
236
237 if (myid==Host_ID){
238
239 sprintf(operate,"%s%s_rst",filepath,filename);
240 mkdir(operate,0775);
241
242 sprintf(fileFT,"%s%s_rst/%s.ftpao",filepath,filename,filename);
243
244 if ((fp = fopen(fileFT,"wb")) != NULL){
245
246 for (spe=0; spe<SpeciesNum; spe++){
247 for (GL=0; GL<=Spe_MaxL_Basis[spe]; GL++){
248 for (Mul=0; Mul<Spe_Num_Basis[spe][GL]; Mul++){
249 fwrite(&Spe_RF_Bessel[spe][GL][Mul][0],sizeof(double),List_YOUSO[15],fp);
250 }
251 }
252 }
253
254 fclose(fp);
255 }
256 else{
257 printf("Could not open a file %s in FT_PAO\n",fileFT);
258 }
259 }
260
261 } /* if (RestartRead_Succeed==0) */
262
263 /***********************************************************
264 elapsed time
265 ***********************************************************/
266
267 dtime(&TEtime);
268
269 /*
270 printf("myid=%2d Elapsed Time (s) = %15.12f\n",myid,TEtime-TStime);
271 MPI_Finalize();
272 exit(0);
273 */
274
275 }
276
277
278
279