1 /**********************************************************************
2 openmx_common.c:
3
4 openmx_common.c is a collective routine of subroutines
5 which are often used.
6
7 Log of openmx_common.c:
8
9 22/Nov/2001 Released by T.Ozaki
10
11 ***********************************************************************/
12
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <string.h>
16 #include <math.h>
17 #include <ctype.h>
18
19 #include "openmx_common.h"
20
21
22
Generation_ATV(int N)23 void Generation_ATV(int N)
24 {
25 int Rn,i,j,k;
26 double di,dj,dk;
27
28 Rn = 1;
29 di = -(N+1);
30 for (i=-N; i<=N; i++){
31 di = di + 1.0;
32 dj = -(N+1);
33 for (j=-N; j<=N; j++){
34 dj = dj + 1.0;
35 dk = -(N+1);
36 for (k=-N; k<=N; k++){
37
38 dk = dk + 1.0;
39 if (i==0 && j==0 && k==0){
40 atv[0][1] = 0.0;
41 atv[0][2] = 0.0;
42 atv[0][3] = 0.0;
43 atv_ijk[0][1] = 0;
44 atv_ijk[0][2] = 0;
45 atv_ijk[0][3] = 0;
46 ratv[i+N][j+N][k+N] = 0;
47 }
48 else{
49 atv[Rn][1] = di*tv[1][1] + dj*tv[2][1] + dk*tv[3][1];
50 atv[Rn][2] = di*tv[1][2] + dj*tv[2][2] + dk*tv[3][2];
51 atv[Rn][3] = di*tv[1][3] + dj*tv[2][3] + dk*tv[3][3];
52 atv_ijk[Rn][1] = i;
53 atv_ijk[Rn][2] = j;
54 atv_ijk[Rn][3] = k;
55 ratv[i+N][j+N][k+N] = Rn;
56 Rn = Rn + 1;
57 }
58 }
59 }
60 }
61
62 }
63
Cross_Product(double a[4],double b[4],double c[4])64 void Cross_Product(double a[4], double b[4], double c[4])
65 {
66 c[1] = a[2]*b[3] - a[3]*b[2];
67 c[2] = a[3]*b[1] - a[1]*b[3];
68 c[3] = a[1]*b[2] - a[2]*b[1];
69 }
70
Dot_Product(double a[4],double b[4])71 double Dot_Product(double a[4], double b[4])
72 {
73 double sum;
74 sum = a[1]*b[1] + a[2]*b[2] + a[3]*b[3];
75 return sum;
76 }
77
78
R_atv(int N,int i,int j,int k)79 int R_atv(int N, int i, int j, int k)
80 {
81 int Rn;
82 Rn = ratv[i+N][j+N][k+N];
83 return Rn;
84 }
85
86
Complex(double re,double im)87 dcomplex Complex(double re, double im)
88 {
89 dcomplex c;
90 c.r = re;
91 c.i = im;
92 return c;
93 }
94
Cadd(dcomplex a,dcomplex b)95 dcomplex Cadd(dcomplex a, dcomplex b)
96 {
97 dcomplex c;
98 c.r = a.r + b.r;
99 c.i = a.i + b.i;
100 return c;
101 }
102
Csub(dcomplex a,dcomplex b)103 dcomplex Csub(dcomplex a, dcomplex b)
104 {
105 dcomplex c;
106 c.r = a.r - b.r;
107 c.i = a.i - b.i;
108 return c;
109 }
110
Cmul(dcomplex a,dcomplex b)111 dcomplex Cmul(dcomplex a, dcomplex b)
112 {
113 dcomplex c;
114 c.r = a.r*b.r - a.i*b.i;
115 c.i = a.i*b.r + a.r*b.i;
116 return c;
117 }
118
Conjg(dcomplex z)119 dcomplex Conjg(dcomplex z)
120 {
121 dcomplex c;
122 c.r = z.r;
123 c.i = -z.i;
124 return c;
125 }
126
Cdiv(dcomplex a,dcomplex b)127 dcomplex Cdiv(dcomplex a, dcomplex b)
128 {
129 dcomplex c;
130 double r,den;
131 if (fabs(b.r) >= fabs(b.i)){
132 r = b.i/b.r;
133 den = b.r + r*b.i;
134 c.r = (a.r + r*a.i)/den;
135 c.i = (a.i - r*a.r)/den;
136 }
137 else{
138 r = b.r/b.i;
139 den = b.i + r*b.r;
140 c.r = (a.r*r + a.i)/den;
141 c.i = (a.i*r - a.r)/den;
142 }
143 return c;
144 }
145
Cabs(dcomplex z)146 double Cabs(dcomplex z)
147 {
148 double x,y,ans,temp;
149 x = fabs(z.r);
150 y = fabs(z.i);
151 if (x<1.0e-30)
152 ans = y;
153 else if (y<1.0e-30)
154 ans = x;
155 else if (x>y){
156 temp = y/x;
157 ans = x*sqrt(1.0+temp*temp);
158 } else{
159 temp = x/y;
160 ans = y*sqrt(1.0+temp*temp);
161 }
162 return ans;
163 }
164
Csqrt(dcomplex z)165 dcomplex Csqrt(dcomplex z)
166 {
167 dcomplex c;
168 double x,y,w,r;
169 if ( fabs(z.r)<1.0e-30 && z.i<1.0e-30 ){
170 c.r = 0.0;
171 c.i = 0.0;
172 return c;
173 }
174 else{
175 x = fabs(z.r);
176 y = fabs(z.i);
177 if (x>=y){
178 r = y/x;
179 w = sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
180 } else {
181 r = x/y;
182 w = sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
183 }
184 if (z.r>=0.0){
185 c.r = w;
186 c.i = z.i/(2.0*w);
187 } else {
188 c.i = (z.i>=0) ? w : -w;
189 c.r = z.i/(2.0*c.i);
190 }
191 return c;
192 }
193 }
194
195
196
Csin(dcomplex a)197 dcomplex Csin(dcomplex a)
198 {
199 double ar;
200
201 if(fabs(a.r) > 2. * PI)
202 ar = (int)(a.r / 2. / PI) * 2. * PI;
203 else
204 ar = a.r;
205
206 return Complex(sin(ar) * cosh(a.i), cos(ar) * sinh(a.i));
207 }
208
209
Ccos(dcomplex a)210 dcomplex Ccos(dcomplex a)
211 {
212 double ar;
213
214 if(fabs(a.r) > 2. * PI)
215 ar = (int)(a.r / 2. / M_PI) * 2. * M_PI;
216 else
217 ar = a.r;
218
219 return Complex(cos(ar) * cosh(a.i), - sin(ar) * sinh(a.i));
220 }
221
222
Cexp(dcomplex a)223 dcomplex Cexp(dcomplex a)
224 {
225 double x;
226
227 x = exp(a.r);
228 return Complex(x * cos(a.i), x * sin(a.i));
229 }
230
231
232
RCadd(double x,dcomplex a)233 dcomplex RCadd(double x, dcomplex a)
234 {
235 dcomplex c;
236 c.r = x + a.r;
237 c.i = a.i;
238 return c;
239 }
240
RCsub(double x,dcomplex a)241 dcomplex RCsub(double x, dcomplex a)
242 {
243 dcomplex c;
244 c.r = x - a.r;
245 c.i = -a.i;
246 return c;
247 }
248
RCmul(double x,dcomplex a)249 dcomplex RCmul(double x, dcomplex a)
250 {
251 dcomplex c;
252 c.r = x*a.r;
253 c.i = x*a.i;
254 return c;
255 }
256
CRmul(dcomplex a,double x)257 dcomplex CRmul(dcomplex a, double x)
258 {
259 dcomplex c;
260 c.r = x*a.r;
261 c.i = x*a.i;
262 return c;
263 }
264
RCdiv(double x,dcomplex a)265 dcomplex RCdiv(double x, dcomplex a)
266 {
267 dcomplex c;
268 double xx,yy,w;
269 xx = a.r;
270 yy = a.i;
271 w = xx*xx+yy*yy;
272 c.r = x*a.r/w;
273 c.i = -x*a.i/w;
274 return c;
275 }
276
CRC(dcomplex a,double x,dcomplex b)277 dcomplex CRC(dcomplex a, double x, dcomplex b)
278 {
279 dcomplex c;
280 c.r = a.r - x - b.r;
281 c.i = a.i - b.i;
282 return c;
283 }
284
285
Cswap(dcomplex * a,dcomplex * b)286 void Cswap(dcomplex *a, dcomplex *b)
287 {
288 dcomplex temp;
289 temp.r = a->r;
290 temp.i = a->i;
291 a->r = b->r;
292 a->i = b->i;
293 b->r = temp.r;
294 b->i = temp.i;
295 }
296
297
298
rnd(double width)299 double rnd(double width)
300 {
301 /****************************************************
302 This rnd() function generates random number
303 -width/2 to width/2
304 ****************************************************/
305
306 double result;
307
308 result = rand();
309 result = result*width/(double)RAND_MAX - 0.5*width;
310 return result;
311 }
312
rnd0to1()313 double rnd0to1()
314 {
315 /****************************************************
316 This rnd() function generates random number 0 to 1
317 ****************************************************/
318
319 double result;
320
321 result = rand();
322 result /= (double)RAND_MAX;
323 return result;
324 }
325
326
sgn(double nu)327 double sgn(double nu)
328 {
329 double result;
330 if (nu<0.0)
331 result = -1.0;
332 else
333 result = 1.0;
334 return result;
335 }
336
isgn(int nu)337 double isgn(int nu)
338 {
339 double result;
340 if (nu<0)
341 result = -1.0;
342 else
343 result = 1.0;
344 return result;
345 }
346
fnjoint(char name1[YOUSO10],char name2[YOUSO10],char name3[YOUSO10])347 void fnjoint(char name1[YOUSO10],char name2[YOUSO10],char name3[YOUSO10])
348 {
349 char name4[YOUSO10];
350 char *f1 = name1,
351 *f2 = name2,
352 *f3 = name3,
353 *f4 = name4;
354
355 while(*f1)
356 {
357 *f4 = *f1;
358 *f1++;
359 *f4++;
360 }
361 while(*f2)
362 {
363 *f4 = *f2;
364 *f2++;
365 *f4++;
366 }
367 while(*f3)
368 {
369 *f4 = *f3;
370 *f3++;
371 *f4++;
372 }
373 *f4 = *f3;
374 chcp(name3,name4);
375 }
376
fnjoint2(char name1[YOUSO10],char name2[YOUSO10],char name3[YOUSO10],char name4[YOUSO10])377 void fnjoint2(char name1[YOUSO10], char name2[YOUSO10],
378 char name3[YOUSO10], char name4[YOUSO10])
379 {
380 char *f1 = name1,
381 *f2 = name2,
382 *f3 = name3,
383 *f4 = name4;
384
385 while(*f1)
386 {
387 *f4 = *f1;
388 *f1++;
389 *f4++;
390 }
391 while(*f2)
392 {
393 *f4 = *f2;
394 *f2++;
395 *f4++;
396 }
397 while(*f3)
398 {
399 *f4 = *f3;
400 *f3++;
401 *f4++;
402 }
403 *f4 = *f3;
404 }
405
406
chcp(char name1[YOUSO10],char name2[YOUSO10])407 void chcp(char name1[YOUSO10],char name2[YOUSO10])
408 {
409
410 /****************************************************
411 name2 -> name1
412 ****************************************************/
413
414 char *f1 = name1,
415 *f2 = name2;
416 while(*f2){
417 *f1 = *f2;
418 *f1++;
419 *f2++;
420 }
421 *f1 = *f2;
422 }
423
SEQ(char str1[YOUSO10],char str2[YOUSO10])424 int SEQ(char str1[YOUSO10], char str2[YOUSO10])
425 {
426
427 int i,result,l1,l2;
428
429 l1 = strlen(str1);
430 l2 = strlen(str2);
431
432 result = 1;
433 if (l1 == l2){
434 for (i=0; i<=l1-1; i++){
435 if (str1[i]!=str2[i]) result = 0;
436 }
437 }
438 else
439 result = 0;
440
441 return result;
442 }
443
444
spline3(double r,double r1,double rcut,double g,double dg,double value[2])445 void spline3(double r, double r1, double rcut,
446 double g, double dg, double value[2])
447 {
448
449 /****************************************************
450 r -> a given distance
451 r1 -> a shortest distatnce in a spline function
452 rcut -> a cut-off distance in a spline function
453 g -> a function value at r1
454 dg -> a derivative at r1
455
456 a function value at r -> value[0]
457 a derivative at r -> value[1]
458 ****************************************************/
459
460 double a0,a1,a2,a3;
461 double rcut2,rcut3,r12,r13,dr;
462
463 rcut2 = rcut*rcut;
464 rcut3 = rcut2*rcut;
465 r12 = r1*r1;
466 r13 = r12*r1;
467 dr = r1 - rcut;
468 a3 = (2.0*g-dg*dr)/(rcut3-r13+3.0*r12*rcut-3.0*r1*rcut2);
469 a2 = 0.5*dg/dr - 1.5*(r1+rcut)*a3;
470 a1 = -rcut*dg/dr + 3.0*r1*rcut*a3;
471 a0 = -a1*rcut-a2*rcut2-a3*rcut3;
472 value[0] = a0+a1*r+a2*r*r+a3*r*r*r;
473 value[1] = a1+2.0*a2*r+3.0*a3*r*r;
474 }
475
largest(double a,double b)476 double largest(double a, double b)
477 {
478 double result;
479
480 if (b<=a) result = a;
481 else result = b;
482 return result;
483 }
484
smallest(double a,double b)485 double smallest(double a, double b)
486 {
487 double result;
488
489 if (b<=a) result = b;
490 else result = a;
491 return result;
492 }
493
494
495
496
497
498
499
500
501
502
503
504
505
asbessel(int n,double x,double sbe[2])506 void asbessel(int n, double x, double sbe[2])
507 {
508
509 /* This rourine suffers from numerical instabilities for a small x */
510
511 double x2,x3,x4,x5,x6,x7,x8;
512
513 if (6<n){
514 printf("n=%2d is not supported in asbessel.",n);
515 exit(0);
516 }
517
518 switch(n){
519
520 case 0:
521 x2 = x*x;
522 sbe[0] = sin(x)/x;
523 sbe[1] = (cos(x) - sin(x)/x)/x;
524 break;
525
526 case 1:
527 x2 = x*x;
528 x3 = x2*x;
529 sbe[0] = -(cos(x)/x) + sin(x)/(x*x);
530 sbe[1] = (2.0*cos(x))/x2 - (2.0*sin(x))/x3 + sin(x)/x;
531 break;
532
533 case 2:
534 x2 = x*x;
535 x3 = x2*x;
536 x4 = x2*x2;
537 sbe[0] = (-3.0*cos(x))/x2 + (3.0*sin(x))/x3 - sin(x)/x;
538 sbe[1] = (9.0*cos(x))/x3 - cos(x)/x - (9.0*sin(x))/x4
539 + (4.0*sin(x))/x2;
540 break;
541
542 case 3:
543 x2 = x*x;
544 x3 = x2*x;
545 x4 = x2*x2;
546 x5 = x4*x;
547 sbe[0] = (-15.0*cos(x))/x3 + cos(x)/x + (15.0*sin(x))/x4
548 - (6.0*sin(x))/x2;
549 sbe[1] = (60.0*cos(x))/x4 - (7.0*cos(x))/x2
550 - (60.0*sin(x))/x5 + (27.0*sin(x))/x3 - sin(x)/x;
551 break;
552
553 case 4:
554 x2 = x*x;
555 x3 = x2*x;
556 x4 = x2*x2;
557 x5 = x4*x;
558 x6 = x5*x;
559 sbe[0] = (-105.0*cos(x))/x4 + (10.0*cos(x))/x2 + (105.0*sin(x))/x5
560 - (45.0*sin(x))/x3 + sin(x)/x;
561 sbe[1] = (525.0*cos(x))/x5 - (65.0*cos(x))/x3 + cos(x)/x
562 - (525.0*sin(x))/x6 + (240.0*sin(x))/x4 - (11.0*sin(x))/x2;
563 break;
564
565 case 5:
566 x2 = x*x;
567 x3 = x2*x;
568 x4 = x2*x2;
569 x5 = x4*x;
570 x6 = x5*x;
571 x7 = x3*x4;
572 sbe[0] = (-945.0*cos(x))/x5 + (105.0*cos(x))/x3 - cos(x)/x
573 + (945.0*sin(x))/x6 - (420.0*sin(x))/x4 + (15.0*sin(x))/x2;
574 sbe[1] = (5670.0*cos(x))/x6 - (735.0*cos(x))/x4 + (16.0*cos(x))/x2
575 - (5670.0*sin(x))/x7 + (2625.0*sin(x))/x5 - (135.0*sin(x))/x3
576 + sin(x)/x;
577 break;
578
579 case 6:
580 x2 = x*x;
581 x3 = x2*x;
582 x4 = x2*x2;
583 x5 = x4*x;
584 x6 = x5*x;
585 x7 = x3*x4;
586 x8 = x4*x4;
587 sbe[0] = (-10395.0*cos(x))/x6 + (1260.0*cos(x))/x4
588 - (21.0*cos(x))/x2 + (10395.0*sin(x))/x7
589 - (4725.0*sin(x))/x5 + (210.0*sin(x))/x3 - sin(x)/x;
590
591 sbe[1] = (72765.0*cos(x))/x7 - (9765.0*cos(x))/x5
592 + (252.0*cos(x))/x3 - cos(x)/x - (72765.0*sin(x))/x8
593 + (34020.0*sin(x))/x6 - (1890.0*sin(x))/x4
594 + (22.0*sin(x))/x2;
595 break;
596
597 }
598 }
599
600
601
602
ComplexSH(int l,int m,double theta,double phi,double SH[2],double dSHt[2],double dSHp[2])603 void ComplexSH(int l, int m, double theta, double phi,
604 double SH[2], double dSHt[2], double dSHp[2])
605 {
606 int i;
607 long double fact0,fact1;
608 double co,si,tmp0,ALeg[2];
609
610 /* Compute (l-|m|)! */
611
612 fact0 = 1.0;
613 for (i=1; i<=(l-abs(m)); i++){
614 fact0 *= i;
615 }
616
617 /* Compute (l+|m|)! */
618 fact1 = 1.0;
619 for (i=1; i<=(l+abs(m)); i++){
620 fact1 *= i;
621 }
622
623 /* sqrt((2*l+1)/(4*PI)*(l-|m|)!/(l+|m|)!) */
624
625 tmp0 = sqrt((2.0*(double)l+1.0)/(4.0*PI)*fact0/fact1);
626
627 /* P_l^|m| */
628
629 Associated_Legendre(l,abs(m),cos(theta),ALeg);
630
631 /* Ylm */
632
633 co = cos((double)m*phi);
634 si = sin((double)m*phi);
635
636 if (0<=m){
637 SH[0] = tmp0*ALeg[0]*co;
638 SH[1] = tmp0*ALeg[0]*si;
639 dSHt[0] = tmp0*ALeg[1]*co;
640 dSHt[1] = tmp0*ALeg[1]*si;
641 dSHp[0] = -(double)m*tmp0*ALeg[0]*si;
642 dSHp[1] = (double)m*tmp0*ALeg[0]*co;
643 }
644 else{
645 if (abs(m)%2==0){
646 SH[0] = tmp0*ALeg[0]*co;
647 SH[1] = tmp0*ALeg[0]*si;
648 dSHt[0] = tmp0*ALeg[1]*co;
649 dSHt[1] = tmp0*ALeg[1]*si;
650 dSHp[0] = -(double)m*tmp0*ALeg[0]*si;
651 dSHp[1] = (double)m*tmp0*ALeg[0]*co;
652 }
653 else{
654 SH[0] = -tmp0*ALeg[0]*co;
655 SH[1] = -tmp0*ALeg[0]*si;
656 dSHt[0] = -tmp0*ALeg[1]*co;
657 dSHt[1] = -tmp0*ALeg[1]*si;
658 dSHp[0] = (double)m*tmp0*ALeg[0]*si;
659 dSHp[1] = -(double)m*tmp0*ALeg[0]*co;
660 }
661 }
662
663 }
664
665
666
667
Associated_Legendre(int l,int m,double x,double ALeg[2])668 void Associated_Legendre(int l, int m, double x, double ALeg[2])
669 {
670 /*****************************************************
671 associated Legendre polynomial Plm(x) with integers
672 m (0<=m<=l) and l. The range of x is -1<=x<=1.
673 Its derivative is given by
674 dP_l^m(x)/dtheta =
675 1/sqrt{1-x*x}*(l*x*Plm(x)-(l+m)*P{l-1}m(x))
676 where x=cos(theta)
677 ******************************************************/
678 double cut0=1.0e-24,cut1=1.0e-12;
679 double Pm,Pm1,f,p0,p1,dP,tmp0;
680 int i,ll;
681
682 if (m<0 || m>l || fabs(x)>1.0){
683 printf("Invalid arguments in routine Associated_Legendre\n");
684 exit(0);
685 }
686 else if ((1.0-cut0)<fabs(x)){
687 x = sgn(x)*(1.0-cut0);
688 }
689
690 /* calculate Pm */
691
692 Pm = 1.0;
693
694 if (m>0){
695
696 f = 1.0;
697 tmp0 = sqrt((1.0 - x)*(1.0 + x));
698 for (i=1; i<=m; i++){
699 Pm = -Pm*f*tmp0;
700 f += 2.0;
701 }
702 }
703
704 if (l==m){
705 p0 = Pm;
706 p1 = 0.0;
707
708 tmp0 = sqrt(1.0-x*x);
709 if (cut1<tmp0) dP = ((double)l*x*p0 - (double)(l+m)*p1)/tmp0;
710 else dP = 0.0;
711
712 ALeg[0] = p0;
713 ALeg[1] = dP;
714 }
715
716 else{
717
718 /* calculate Pm1 */
719
720 Pm1 = x*(2.0*(double)m + 1.0)*Pm;
721
722 if (l==(m+1)){
723 p0 = Pm1;
724 p1 = Pm;
725 tmp0 = sqrt(1.0-x*x);
726
727 if (cut1<tmp0) dP = ((double)l*x*p0 - (double)(l+m)*p1)/tmp0;
728 else dP = 0.0;
729
730 ALeg[0] = p0;
731 ALeg[1] = dP;
732 }
733
734 /* calculate Plm, l>m+1 */
735
736 else{
737
738 for (ll=m+2; ll<=l; ll++){
739 tmp0 = (x*(2.0*(double)ll-1.0)*Pm1 - ((double)ll+(double)m-1.0)*Pm)/(double)(ll-m);
740 Pm = Pm1;
741 Pm1 = tmp0;
742 }
743 p0 = Pm1;
744 p1 = Pm;
745
746 tmp0 = sqrt(1.0-x*x);
747
748 if (cut1<tmp0) dP = ((double)l*x*p0 - (double)(l+m)*p1)/tmp0;
749 else dP = 0.0;
750
751 ALeg[0] = p0;
752 ALeg[1] = dP;
753 }
754 }
755 }
756
757
758
Im_pow(int fu,int Ls)759 dcomplex Im_pow(int fu, int Ls)
760 {
761
762 dcomplex Cres;
763
764 if (fu!=1 && fu!=-1){
765 printf("Invalid arguments in Im_pow\n");
766 exit(0);
767 }
768 else{
769
770 if (fu==1){
771 if (Ls%2==0){
772 if (Ls%4==0) Cres = Complex( 1.0, 0.0);
773 else Cres = Complex(-1.0, 0.0);
774 }
775 else{
776 if ((Ls+1)%4==0) Cres = Complex( 0.0,-1.0);
777 else Cres = Complex( 0.0, 1.0);
778 }
779 }
780
781 else{
782 if (Ls%2==0){
783 if (Ls%4==0) Cres = Complex( 1.0, 0.0);
784 else Cres = Complex(-1.0, 0.0);
785 }
786 else{
787 if ((Ls+1)%4==0) Cres = Complex( 0.0, 1.0);
788 else Cres = Complex( 0.0,-1.0);
789 }
790 }
791
792 }
793
794 return Cres;
795 }
796
GN2N(int GN,int N3[4])797 void GN2N(int GN, int N3[4])
798 {
799 int n1,n2,n3;
800
801 n1 = GN/(Ngrid2*Ngrid3);
802 n2 = (GN - n1*(Ngrid2*Ngrid3))/Ngrid3;
803 n3 = GN - n1*(Ngrid2*Ngrid3) - n2*Ngrid3;
804 N3[1] = n1;
805 N3[2] = n2;
806 N3[3] = n3;
807 }
808
AproxFactN(int N0)809 int AproxFactN(int N0)
810 {
811 int N1,N18,po,N[5],i;
812 int result;
813
814 printf("AproxFactN: N0=%d\n",N0);
815
816 if (N0<=4){
817 result = 4;
818 }
819 else{
820
821 po = 0;
822 N1 = 1;
823 do{
824 N1 = 2*N1;
825 if (N0<N1){
826 N18 = N1/16;
827 po = 1;
828 }
829 } while (po==0);
830
831 printf("AproxFactN: N18=%d\n",N18);
832
833 N[0] = N18*4;
834 N[1] = N18*5;
835 N[2] = N18*6;
836 N[3] = N18*7;
837 N[4] = N18*8;
838
839 po = 0;
840 i = -1;
841 do{
842 i++;
843 printf("AproxFactN: i,N[i],N0=%d %d %d\n",i,N[i],N0);
844 if (0<=(N[i]-N0)) po = 1;
845 } while(po==0);
846
847 result = N[i];
848 }
849
850 return result;
851 }
852
853
854
Get_Grid_XYZ(int GN,double xyz[4])855 void Get_Grid_XYZ(int GN, double xyz[4])
856 {
857 int n1,n2,n3;
858
859 n1 = GN/(Ngrid2*Ngrid3);
860 n2 = (GN - n1*(Ngrid2*Ngrid3))/Ngrid3;
861 n3 = GN - n1*(Ngrid2*Ngrid3) - n2*Ngrid3;
862
863 xyz[1] = (double)n1*gtv[1][1] + (double)n2*gtv[2][1]
864 + (double)n3*gtv[3][1] + Grid_Origin[1];
865 xyz[2] = (double)n1*gtv[1][2] + (double)n2*gtv[2][2]
866 + (double)n3*gtv[3][2] + Grid_Origin[2];
867 xyz[3] = (double)n1*gtv[1][3] + (double)n2*gtv[2][3]
868 + (double)n3*gtv[3][3] + Grid_Origin[3];
869 }
870
871
k_inversion(int i,int j,int k,int mi,int mj,int mk,int * ii,int * ij,int * ik)872 void k_inversion(int i, int j, int k,
873 int mi, int mj, int mk,
874 int *ii, int *ij, int *ik )
875 {
876 *ii= mi-i-1;
877 *ij= mj-j-1;
878 *ik= mk-k-1;
879 }
880
881
string_tolower(char * buf,char * buf1)882 char *string_tolower(char *buf, char *buf1)
883 {
884 char *c=buf;
885 char *c1=buf1;
886
887 while (*c){
888 *c1=tolower(*c);
889 c++;
890 c1++;
891 }
892 return buf;
893 }
894