1 /**********************************************************************
2 Set_Orbitals_Grid.c:
3
4 Set_Orbitals_Grid.c is a subroutine to calculate the value of basis
5 functions on each grid point.
6
7 Log of Set_Orbitals_Grid.c:
8
9 22/Nov/2001 Released by T.Ozaki
10
11 ***********************************************************************/
12
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <time.h>
16 #include <math.h>
17 #include "openmx_common.h"
18 #include "mpi.h"
19 #include <omp.h>
20
21
22
23
Set_Orbitals_Grid(int Cnt_kind)24 double Set_Orbitals_Grid(int Cnt_kind)
25 {
26 int i,j,n,Mc_AN,Gc_AN,Cwan,NO0,GNc,GRc;
27 int Gh_AN,Mh_AN,Rnh,Hwan,NO1,Nog,h_AN;
28 long int k,Nc;
29 double time0;
30 double x,y,z,dx,dy,dz;
31 double TStime,TEtime;
32 double Cxyz[4];
33 int numprocs,myid,tag=999,ID,IDS,IDR;
34 double Stime_atom,Etime_atom;
35
36 MPI_Status stat;
37 MPI_Request request;
38
39 /* for OpenMP */
40 int OMPID,Nthrds,Nprocs;
41
42 /* MPI */
43 MPI_Comm_size(mpi_comm_level1,&numprocs);
44 MPI_Comm_rank(mpi_comm_level1,&myid);
45
46 dtime(&TStime);
47
48 /*****************************************************
49 Calculate orbitals on grids
50 *****************************************************/
51
52 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
53
54 dtime(&Stime_atom);
55
56 Gc_AN = M2G[Mc_AN];
57 Cwan = WhatSpecies[Gc_AN];
58
59 if (Cnt_kind==0) NO0 = Spe_Total_NO[Cwan];
60 else NO0 = Spe_Total_CNO[Cwan];
61
62 #pragma omp parallel shared(Comp2Real,Spe_PAO_RWF,Spe_Num_Basis,Spe_MaxL_Basis,Spe_PAO_RV,Spe_Num_Mesh_PAO,List_YOUSO,Orbs_Grid,Cnt_kind,Gxyz,atv,CellListAtom,GridListAtom,GridN_Atom,Gc_AN,Cwan,Mc_AN,NO0) private(OMPID,Nthrds,Nprocs,Nc,GNc,GRc,Cxyz,x,y,z,dx,dy,dz,i,j)
63 {
64 double *Chi0;
65 double Cxyz0[4];
66 double **RF;
67 double **AF;
68 int i,L0,Mul0,M0,i1;
69 double S_coordinate[3];
70 double dum,dum1,dum2,dum3,dum4,a,b,c,d;
71 double siQ,coQ,siP,coP,Q,P,R;
72 double rm,df,sum0,sum1;
73 double SH[Supported_MaxL*2+1][2];
74 double dSHt[Supported_MaxL*2+1][2];
75 double dSHp[Supported_MaxL*2+1][2];
76
77 /* Radial */
78 int mp_min,mp_max,m,po,wan;
79 double h1,h2,h3,f1,f2,f3,f4;
80 double g1,g2,x1,x2,y1,y2,y12,y22,f;
81 double r,r1,theta,phi,Min_r;
82
83 /* allocation of array */
84
85 Chi0 = (double*)malloc(sizeof(double)*List_YOUSO[7]);
86
87 RF = (double**)malloc(sizeof(double*)*(List_YOUSO[25]+1));
88 for (i=0; i<(List_YOUSO[25]+1); i++){
89 RF[i] = (double*)malloc(sizeof(double)*List_YOUSO[24]);
90 }
91
92 AF = (double**)malloc(sizeof(double*)*(List_YOUSO[25]+1));
93 for (i=0; i<(List_YOUSO[25]+1); i++){
94 AF[i] = (double*)malloc(sizeof(double)*(2*(List_YOUSO[25]+1)+1));
95 }
96
97 /* get info. on OpenMP */
98
99 OMPID = omp_get_thread_num();
100 Nthrds = omp_get_num_threads();
101 Nprocs = omp_get_num_procs();
102
103 for (Nc=OMPID*GridN_Atom[Gc_AN]/Nthrds; Nc<(OMPID+1)*GridN_Atom[Gc_AN]/Nthrds; Nc++){
104
105 GNc = GridListAtom[Mc_AN][Nc];
106 GRc = CellListAtom[Mc_AN][Nc];
107
108 Get_Grid_XYZ(GNc,Cxyz);
109 x = Cxyz[1] + atv[GRc][1] - Gxyz[Gc_AN][1];
110 y = Cxyz[2] + atv[GRc][2] - Gxyz[Gc_AN][2];
111 z = Cxyz[3] + atv[GRc][3] - Gxyz[Gc_AN][3];
112
113 if (Cnt_kind==0){
114
115 /* Get_Orbitals(Cwan,x,y,z,Chi0); */
116 /* start of inlining of Get_Orbitals */
117
118 wan = Cwan;
119
120 /* xyz2spherical(x,y,z,0.0,0.0,0.0,S_coordinate); */
121 /* start of inlining of xyz2spherical */
122
123 Min_r = 10e-15;
124 dum = x*x + y*y;
125 r = sqrt(dum + z*z);
126 r1 = sqrt(dum);
127
128 if (Min_r<=r){
129
130 if (r<fabs(z))
131 dum1 = sgn(z)*1.0;
132 else
133 dum1 = z/r;
134
135 theta = acos(dum1);
136
137 if (Min_r<=r1){
138 if (0.0<=x){
139
140 if (r1<fabs(y))
141 dum1 = sgn(y)*1.0;
142 else
143 dum1 = y/r1;
144
145 phi = asin(dum1);
146 }
147 else{
148
149 if (r1<fabs(y))
150 dum1 = sgn(y)*1.0;
151 else
152 dum1 = y/r1;
153
154 phi = PI - asin(dum1);
155 }
156 }
157 else{
158 phi = 0.0;
159 }
160 }
161 else{
162 theta = 0.5*PI;
163 phi = 0.0;
164 }
165
166 R = r;
167 Q = theta;
168 P = phi;
169
170 /* end of inlining of xyz2spherical */
171
172 po = 0;
173 mp_min = 0;
174 mp_max = Spe_Num_Mesh_PAO[wan] - 1;
175
176 if (Spe_PAO_RV[wan][Spe_Num_Mesh_PAO[wan]-1]<R){
177
178 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
179 for (Mul0=0; Mul0<Spe_Num_Basis[wan][L0]; Mul0++){
180 RF[L0][Mul0] = 0.0;
181 }
182 }
183
184 po = 1;
185 }
186
187 else if (R<Spe_PAO_RV[wan][0]){
188
189 m = 4;
190 rm = Spe_PAO_RV[wan][m];
191
192 h1 = Spe_PAO_RV[wan][m-1] - Spe_PAO_RV[wan][m-2];
193 h2 = Spe_PAO_RV[wan][m] - Spe_PAO_RV[wan][m-1];
194 h3 = Spe_PAO_RV[wan][m+1] - Spe_PAO_RV[wan][m];
195
196 x1 = rm - Spe_PAO_RV[wan][m-1];
197 x2 = rm - Spe_PAO_RV[wan][m];
198 y1 = x1/h2;
199 y2 = x2/h2;
200 y12 = y1*y1;
201 y22 = y2*y2;
202
203 dum = h1 + h2;
204 dum1 = h1/h2/dum;
205 dum2 = h2/h1/dum;
206 dum = h2 + h3;
207 dum3 = h2/h3/dum;
208 dum4 = h3/h2/dum;
209
210 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
211 for (Mul0=0; Mul0<Spe_Num_Basis[wan][L0]; Mul0++){
212
213 f1 = Spe_PAO_RWF[wan][L0][Mul0][m-2];
214 f2 = Spe_PAO_RWF[wan][L0][Mul0][m-1];
215 f3 = Spe_PAO_RWF[wan][L0][Mul0][m];
216 f4 = Spe_PAO_RWF[wan][L0][Mul0][m+1];
217
218 if (m==1){
219 h1 = -(h2+h3);
220 f1 = f4;
221 }
222 else if (m==(Spe_Num_Mesh_PAO[wan]-1)){
223 h3 = -(h1+h2);
224 f4 = f1;
225 }
226
227 dum = f3 - f2;
228 g1 = dum*dum1 + (f2-f1)*dum2;
229 g2 = (f4-f3)*dum3 + dum*dum4;
230
231 f = y22*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2)
232 + y12*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1);
233
234 df = 2.0*y2/h2*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2)
235 + y22*(2.0*f2 + h2*g1)/h2
236 + 2.0*y1/h2*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1)
237 - y12*(2.0*f3 - h2*g2)/h2;
238
239 if (L0==0){
240 a = 0.0;
241 b = 0.5*df/rm;
242 c = 0.0;
243 d = f - b*rm*rm;
244 }
245
246 else if (L0==1){
247 a = (rm*df - f)/(2.0*rm*rm*rm);
248 b = 0.0;
249 c = df - 3.0*a*rm*rm;
250 d = 0.0;
251 }
252
253 else{
254 b = (3.0*f - rm*df)/(rm*rm);
255 a = (f - b*rm*rm)/(rm*rm*rm);
256 c = 0.0;
257 d = 0.0;
258 }
259
260 RF[L0][Mul0] = a*R*R*R + b*R*R + c*R + d;
261
262 }
263 }
264
265 }
266
267 else{
268
269 do{
270 m = (mp_min + mp_max)/2;
271 if (Spe_PAO_RV[wan][m]<R)
272 mp_min = m;
273 else
274 mp_max = m;
275 }
276 while((mp_max-mp_min)!=1);
277 m = mp_max;
278
279 h1 = Spe_PAO_RV[wan][m-1] - Spe_PAO_RV[wan][m-2];
280 h2 = Spe_PAO_RV[wan][m] - Spe_PAO_RV[wan][m-1];
281 h3 = Spe_PAO_RV[wan][m+1] - Spe_PAO_RV[wan][m];
282
283 x1 = R - Spe_PAO_RV[wan][m-1];
284 x2 = R - Spe_PAO_RV[wan][m];
285 y1 = x1/h2;
286 y2 = x2/h2;
287 y12 = y1*y1;
288 y22 = y2*y2;
289
290 dum = h1 + h2;
291 dum1 = h1/h2/dum;
292 dum2 = h2/h1/dum;
293 dum = h2 + h3;
294 dum3 = h2/h3/dum;
295 dum4 = h3/h2/dum;
296
297 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
298 for (Mul0=0; Mul0<Spe_Num_Basis[wan][L0]; Mul0++){
299
300 f1 = Spe_PAO_RWF[wan][L0][Mul0][m-2];
301 f2 = Spe_PAO_RWF[wan][L0][Mul0][m-1];
302 f3 = Spe_PAO_RWF[wan][L0][Mul0][m];
303 f4 = Spe_PAO_RWF[wan][L0][Mul0][m+1];
304
305 if (m==1){
306 h1 = -(h2+h3);
307 f1 = f4;
308 }
309 else if (m==(Spe_Num_Mesh_PAO[wan]-1)){
310 h3 = -(h1+h2);
311 f4 = f1;
312 }
313
314 dum = f3 - f2;
315 g1 = dum*dum1 + (f2-f1)*dum2;
316 g2 = (f4-f3)*dum3 + dum*dum4;
317
318 f = y22*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2)
319 + y12*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1);
320
321 RF[L0][Mul0] = f;
322
323 }
324 }
325
326 }
327
328 if (po==0){
329
330 /* Angular */
331 siQ = sin(Q);
332 coQ = cos(Q);
333 siP = sin(P);
334 coP = cos(P);
335
336 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
337
338 if (L0==0){
339 AF[0][0] = 0.282094791773878;
340 }
341 else if (L0==1){
342 dum = 0.48860251190292*siQ;
343 AF[1][0] = dum*coP;
344 AF[1][1] = dum*siP;
345 AF[1][2] = 0.48860251190292*coQ;
346 }
347 else if (L0==2){
348 dum1 = siQ*siQ;
349 dum2 = 1.09254843059208*siQ*coQ;
350 AF[2][0] = 0.94617469575756*coQ*coQ - 0.31539156525252;
351 AF[2][1] = 0.54627421529604*dum1*(1.0 - 2.0*siP*siP);
352 AF[2][2] = 1.09254843059208*dum1*siP*coP;
353 AF[2][3] = dum2*coP;
354 AF[2][4] = dum2*siP;
355 }
356
357 else if (L0==3){
358 AF[3][0] = 0.373176332590116*(5.0*coQ*coQ*coQ - 3.0*coQ);
359 AF[3][1] = 0.457045799464466*coP*siQ*(5.0*coQ*coQ - 1.0);
360 AF[3][2] = 0.457045799464466*siP*siQ*(5.0*coQ*coQ - 1.0);
361 AF[3][3] = 1.44530572132028*siQ*siQ*coQ*(coP*coP-siP*siP);
362 AF[3][4] = 2.89061144264055*siQ*siQ*coQ*siP*coP;
363 AF[3][5] = 0.590043589926644*siQ*siQ*siQ*(4.0*coP*coP*coP - 3.0*coP);
364 AF[3][6] = 0.590043589926644*siQ*siQ*siQ*(3.0*siP - 4.0*siP*siP*siP);
365 }
366
367 else if (4<=L0){
368
369 /* calculation of complex spherical harmonics functions */
370 for(m=-L0; m<=L0; m++){
371 ComplexSH(L0,m,Q,P,SH[L0+m],dSHt[L0+m],dSHp[L0+m]);
372 }
373
374 /* transformation of complex to real */
375 for (i=0; i<(L0*2+1); i++){
376
377 sum0 = 0.0;
378 sum1 = 0.0;
379 for (j=0; j<(L0*2+1); j++){
380 sum0 += Comp2Real[L0][i][j].r*SH[j][0] - Comp2Real[L0][i][j].i*SH[j][1];
381 sum1 += Comp2Real[L0][i][j].r*SH[j][1] + Comp2Real[L0][i][j].i*SH[j][0];
382 }
383 AF[L0][i] = sum0 + sum1;
384 }
385
386 }
387 }
388 }
389
390 /* Chi0 */
391 i1 = -1;
392 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
393 for (Mul0=0; Mul0<Spe_Num_Basis[wan][L0]; Mul0++){
394 for (M0=0; M0<=2*L0; M0++){
395 i1++;
396 Chi0[i1] = RF[L0][Mul0]*AF[L0][M0];
397 }
398 }
399 }
400 /* end of inlining of Get_Orbitals */
401
402 }
403 else{
404 Get_Cnt_Orbitals(Mc_AN,x,y,z,Chi0);
405 }
406
407 for (i=0; i<NO0; i++){
408 Orbs_Grid[Mc_AN][Nc][i] = (Type_Orbs_Grid)Chi0[i];/* AITUNE */
409 }
410
411 } /* Nc */
412
413 /* freeing of array */
414
415 free(Chi0);
416
417 for (i=0; i<(List_YOUSO[25]+1); i++){
418 free(RF[i]);
419 }
420 free(RF);
421
422 for (i=0; i<(List_YOUSO[25]+1); i++){
423 free(AF[i]);
424 }
425 free(AF);
426
427 } /* #pragma omp parallel */
428
429 dtime(&Etime_atom);
430 time_per_atom[Gc_AN] += Etime_atom - Stime_atom;
431 }
432
433 /****************************************************
434 Calculate Orbs_Grid_FNAN
435 ****************************************************/
436
437 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
438
439 Gc_AN = M2G[Mc_AN];
440
441 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
442
443 Gh_AN = natn[Gc_AN][h_AN];
444
445 if (G2ID[Gh_AN]!=myid){
446
447 Mh_AN = F_G2M[Gh_AN];
448 Rnh = ncn[Gc_AN][h_AN];
449 Hwan = WhatSpecies[Gh_AN];
450
451 if (Cnt_kind==0) NO1 = Spe_Total_NO[Hwan];
452 else NO1 = Spe_Total_CNO[Hwan];
453
454 #pragma omp parallel shared(List_YOUSO,Orbs_Grid_FNAN,NO1,Mh_AN,Hwan,Cnt_kind,Rnh,Gh_AN,Gxyz,atv,NumOLG,Mc_AN,h_AN,GListTAtoms1,GridListAtom,CellListAtom) private(OMPID,Nthrds,Nprocs,Nog,Nc,GNc,GRc,x,y,z,j)
455 {
456
457 double *Chi0;
458 double Cxyz0[4];
459 double **RF;
460 double **AF;
461 int i,L0,Mul0,M0,i1;
462 double S_coordinate[3];
463 double dum,dum1,dum2,dum3,dum4,a,b,c,d;
464 double siQ,coQ,siP,coP,Q,P,R;
465 double rm,df,sum0,sum1;
466 double SH[Supported_MaxL*2+1][2];
467 double dSHt[Supported_MaxL*2+1][2];
468 double dSHp[Supported_MaxL*2+1][2];
469
470 /* Radial */
471 int mp_min,mp_max,m,po,wan;
472 double h1,h2,h3,f1,f2,f3,f4;
473 double g1,g2,x1,x2,y1,y2,y12,y22,f;
474 double r,r1,theta,phi,Min_r;
475
476 /* allocation of arrays */
477
478 Chi0 = (double*)malloc(sizeof(double)*List_YOUSO[7]);
479
480 RF = (double**)malloc(sizeof(double*)*(List_YOUSO[25]+1));
481 for (i=0; i<(List_YOUSO[25]+1); i++){
482 RF[i] = (double*)malloc(sizeof(double)*List_YOUSO[24]);
483 }
484
485 AF = (double**)malloc(sizeof(double*)*(List_YOUSO[25]+1));
486 for (i=0; i<(List_YOUSO[25]+1); i++){
487 AF[i] = (double*)malloc(sizeof(double)*(2*(List_YOUSO[25]+1)+1));
488 }
489
490 /* get info. on OpenMP */
491
492 OMPID = omp_get_thread_num();
493 Nthrds = omp_get_num_threads();
494 Nprocs = omp_get_num_procs();
495
496 for (Nog=OMPID*NumOLG[Mc_AN][h_AN]/Nthrds; Nog<(OMPID+1)*NumOLG[Mc_AN][h_AN]/Nthrds; Nog++){
497
498 Nc = GListTAtoms1[Mc_AN][h_AN][Nog];
499 GNc = GridListAtom[Mc_AN][Nc];
500 GRc = CellListAtom[Mc_AN][Nc];
501
502 Get_Grid_XYZ(GNc,Cxyz0);
503
504 x = Cxyz0[1] + atv[GRc][1] - Gxyz[Gh_AN][1] - atv[Rnh][1];
505 y = Cxyz0[2] + atv[GRc][2] - Gxyz[Gh_AN][2] - atv[Rnh][2];
506 z = Cxyz0[3] + atv[GRc][3] - Gxyz[Gh_AN][3] - atv[Rnh][3];
507
508 if (Cnt_kind==0){
509
510 /* Get_Orbitals(Hwan,x,y,z,Chi0); */
511 /* start of inlining of Get_Orbitals */
512
513 wan = Hwan;
514
515 /* xyz2spherical(x,y,z,0.0,0.0,0.0,S_coordinate); */
516 /* start of inlining of xyz2spherical */
517
518 Min_r = 10e-15;
519 dum = x*x + y*y;
520 r = sqrt(dum + z*z);
521 r1 = sqrt(dum);
522
523 if (Min_r<=r){
524
525 if (r<fabs(z))
526 dum1 = sgn(z)*1.0;
527 else
528 dum1 = z/r;
529
530 theta = acos(dum1);
531
532 if (Min_r<=r1){
533 if (0.0<=x){
534
535 if (r1<fabs(y))
536 dum1 = sgn(y)*1.0;
537 else
538 dum1 = y/r1;
539
540 phi = asin(dum1);
541 }
542 else{
543
544 if (r1<fabs(y))
545 dum1 = sgn(y)*1.0;
546 else
547 dum1 = y/r1;
548
549 phi = PI - asin(dum1);
550 }
551 }
552 else{
553 phi = 0.0;
554 }
555 }
556 else{
557 theta = 0.5*PI;
558 phi = 0.0;
559 }
560
561 R = r;
562 Q = theta;
563 P = phi;
564
565 /* end of inlining of xyz2spherical */
566
567 po = 0;
568 mp_min = 0;
569 mp_max = Spe_Num_Mesh_PAO[wan] - 1;
570
571 if (Spe_PAO_RV[wan][Spe_Num_Mesh_PAO[wan]-1]<R){
572
573 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
574 for (Mul0=0; Mul0<Spe_Num_Basis[wan][L0]; Mul0++){
575 RF[L0][Mul0] = 0.0;
576 }
577 }
578
579 po = 1;
580 }
581
582 else if (R<Spe_PAO_RV[wan][0]){
583
584 m = 4;
585 rm = Spe_PAO_RV[wan][m];
586
587 h1 = Spe_PAO_RV[wan][m-1] - Spe_PAO_RV[wan][m-2];
588 h2 = Spe_PAO_RV[wan][m] - Spe_PAO_RV[wan][m-1];
589 h3 = Spe_PAO_RV[wan][m+1] - Spe_PAO_RV[wan][m];
590
591 x1 = rm - Spe_PAO_RV[wan][m-1];
592 x2 = rm - Spe_PAO_RV[wan][m];
593 y1 = x1/h2;
594 y2 = x2/h2;
595 y12 = y1*y1;
596 y22 = y2*y2;
597
598 dum = h1 + h2;
599 dum1 = h1/h2/dum;
600 dum2 = h2/h1/dum;
601 dum = h2 + h3;
602 dum3 = h2/h3/dum;
603 dum4 = h3/h2/dum;
604
605 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
606 for (Mul0=0; Mul0<Spe_Num_Basis[wan][L0]; Mul0++){
607
608 f1 = Spe_PAO_RWF[wan][L0][Mul0][m-2];
609 f2 = Spe_PAO_RWF[wan][L0][Mul0][m-1];
610 f3 = Spe_PAO_RWF[wan][L0][Mul0][m];
611 f4 = Spe_PAO_RWF[wan][L0][Mul0][m+1];
612
613 if (m==1){
614 h1 = -(h2+h3);
615 f1 = f4;
616 }
617 else if (m==(Spe_Num_Mesh_PAO[wan]-1)){
618 h3 = -(h1+h2);
619 f4 = f1;
620 }
621
622 dum = f3 - f2;
623 g1 = dum*dum1 + (f2-f1)*dum2;
624 g2 = (f4-f3)*dum3 + dum*dum4;
625
626 f = y22*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2)
627 + y12*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1);
628
629 df = 2.0*y2/h2*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2)
630 + y22*(2.0*f2 + h2*g1)/h2
631 + 2.0*y1/h2*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1)
632 - y12*(2.0*f3 - h2*g2)/h2;
633
634 if (L0==0){
635 a = 0.0;
636 b = 0.5*df/rm;
637 c = 0.0;
638 d = f - b*rm*rm;
639 }
640
641 else if (L0==1){
642 a = (rm*df - f)/(2.0*rm*rm*rm);
643 b = 0.0;
644 c = df - 3.0*a*rm*rm;
645 d = 0.0;
646 }
647
648 else{
649 b = (3.0*f - rm*df)/(rm*rm);
650 a = (f - b*rm*rm)/(rm*rm*rm);
651 c = 0.0;
652 d = 0.0;
653 }
654
655 RF[L0][Mul0] = a*R*R*R + b*R*R + c*R + d;
656
657 }
658 }
659
660 }
661
662 else{
663
664 do{
665 m = (mp_min + mp_max)/2;
666 if (Spe_PAO_RV[wan][m]<R)
667 mp_min = m;
668 else
669 mp_max = m;
670 }
671 while((mp_max-mp_min)!=1);
672 m = mp_max;
673
674 h1 = Spe_PAO_RV[wan][m-1] - Spe_PAO_RV[wan][m-2];
675 h2 = Spe_PAO_RV[wan][m] - Spe_PAO_RV[wan][m-1];
676 h3 = Spe_PAO_RV[wan][m+1] - Spe_PAO_RV[wan][m];
677
678 x1 = R - Spe_PAO_RV[wan][m-1];
679 x2 = R - Spe_PAO_RV[wan][m];
680 y1 = x1/h2;
681 y2 = x2/h2;
682 y12 = y1*y1;
683 y22 = y2*y2;
684
685 dum = h1 + h2;
686 dum1 = h1/h2/dum;
687 dum2 = h2/h1/dum;
688 dum = h2 + h3;
689 dum3 = h2/h3/dum;
690 dum4 = h3/h2/dum;
691
692 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
693 for (Mul0=0; Mul0<Spe_Num_Basis[wan][L0]; Mul0++){
694
695 f1 = Spe_PAO_RWF[wan][L0][Mul0][m-2];
696 f2 = Spe_PAO_RWF[wan][L0][Mul0][m-1];
697 f3 = Spe_PAO_RWF[wan][L0][Mul0][m];
698 f4 = Spe_PAO_RWF[wan][L0][Mul0][m+1];
699
700 if (m==1){
701 h1 = -(h2+h3);
702 f1 = f4;
703 }
704 else if (m==(Spe_Num_Mesh_PAO[wan]-1)){
705 h3 = -(h1+h2);
706 f4 = f1;
707 }
708
709 dum = f3 - f2;
710 g1 = dum*dum1 + (f2-f1)*dum2;
711 g2 = (f4-f3)*dum3 + dum*dum4;
712
713 f = y22*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2)
714 + y12*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1);
715
716 RF[L0][Mul0] = f;
717
718 }
719 }
720
721 }
722
723 if (po==0){
724
725 /* Angular */
726 siQ = sin(Q);
727 coQ = cos(Q);
728 siP = sin(P);
729 coP = cos(P);
730
731 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
732
733 if (L0==0){
734 AF[0][0] = 0.282094791773878;
735 }
736 else if (L0==1){
737 dum = 0.48860251190292*siQ;
738 AF[1][0] = dum*coP;
739 AF[1][1] = dum*siP;
740 AF[1][2] = 0.48860251190292*coQ;
741 }
742 else if (L0==2){
743 dum1 = siQ*siQ;
744 dum2 = 1.09254843059208*siQ*coQ;
745 AF[2][0] = 0.94617469575756*coQ*coQ - 0.31539156525252;
746 AF[2][1] = 0.54627421529604*dum1*(1.0 - 2.0*siP*siP);
747 AF[2][2] = 1.09254843059208*dum1*siP*coP;
748 AF[2][3] = dum2*coP;
749 AF[2][4] = dum2*siP;
750 }
751
752 else if (L0==3){
753 AF[3][0] = 0.373176332590116*(5.0*coQ*coQ*coQ - 3.0*coQ);
754 AF[3][1] = 0.457045799464466*coP*siQ*(5.0*coQ*coQ - 1.0);
755 AF[3][2] = 0.457045799464466*siP*siQ*(5.0*coQ*coQ - 1.0);
756 AF[3][3] = 1.44530572132028*siQ*siQ*coQ*(coP*coP-siP*siP);
757 AF[3][4] = 2.89061144264055*siQ*siQ*coQ*siP*coP;
758 AF[3][5] = 0.590043589926644*siQ*siQ*siQ*(4.0*coP*coP*coP - 3.0*coP);
759 AF[3][6] = 0.590043589926644*siQ*siQ*siQ*(3.0*siP - 4.0*siP*siP*siP);
760 }
761
762 else if (4<=L0){
763
764 /* calculation of complex spherical harmonics functions */
765 for(m=-L0; m<=L0; m++){
766 ComplexSH(L0,m,Q,P,SH[L0+m],dSHt[L0+m],dSHp[L0+m]);
767 }
768
769 /* transformation of complex to real */
770 for (i=0; i<(L0*2+1); i++){
771
772 sum0 = 0.0;
773 sum1 = 0.0;
774 for (j=0; j<(L0*2+1); j++){
775 sum0 += Comp2Real[L0][i][j].r*SH[j][0] - Comp2Real[L0][i][j].i*SH[j][1];
776 sum1 += Comp2Real[L0][i][j].r*SH[j][1] + Comp2Real[L0][i][j].i*SH[j][0];
777 }
778 AF[L0][i] = sum0 + sum1;
779 }
780
781 }
782 }
783 }
784
785 /* Chi0 */
786 i1 = -1;
787 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
788 for (Mul0=0; Mul0<Spe_Num_Basis[wan][L0]; Mul0++){
789 for (M0=0; M0<=2*L0; M0++){
790 i1++;
791 Chi0[i1] = RF[L0][Mul0]*AF[L0][M0];
792 }
793 }
794 }
795
796 /* end of inlining of Get_Orbitals */
797
798 } /* if (Cnt_kind==0) */
799
800 else{
801 Get_Cnt_Orbitals(Mh_AN,x,y,z,Chi0);
802 }
803
804 for (j=0; j<NO1; j++){
805 Orbs_Grid_FNAN[Mc_AN][h_AN][Nog][j] = (Type_Orbs_Grid)Chi0[j];/* AITUNE */
806 }
807
808 } /* Nog */
809
810 /* freeing of arrays */
811
812 free(Chi0);
813
814 for (i=0; i<(List_YOUSO[25]+1); i++){
815 free(RF[i]);
816 }
817 free(RF);
818
819 for (i=0; i<(List_YOUSO[25]+1); i++){
820 free(AF[i]);
821 }
822 free(AF);
823
824 }
825 }
826 }
827 }
828
829 /* time */
830 dtime(&TEtime);
831 time0 = TEtime - TStime;
832
833 return time0;
834 }
835