1 /**********************************************************************
2 RF_BesselF.c:
3
4 RF_BesselF.c is a subroutine to calculate radial part of PAO of
5 atom specified by "Gensi" in k-space.
6
7 Log of RF_BesselF.c:
8
9 22/Nov/2001 Released by T.Ozaki
10
11 ***********************************************************************/
12
13 #include <stdio.h>
14 #include <math.h>
15 #include "openmx_common.h"
16
RF_BesselF(int Gensi,int GL,int Mul,double R)17 double RF_BesselF(int Gensi, int GL, int Mul, double R)
18 {
19 int mp_min,mp_max,m,po;
20 double h1,h2,h3,f1,f2,f3,f4;
21 double g1,g2,x1,x2,y1,y2,f;
22 double result;
23
24 mp_min = 0;
25 mp_max = Ngrid_NormK - 1;
26 po = 0;
27
28 if (R<NormK[0]){
29 m = 1;
30 }
31 else if (NormK[mp_max]<R){
32 result = 0.0;
33 po = 1;
34 }
35 else{
36 do{
37 m = (mp_min + mp_max)/2;
38 if (NormK[m]<R)
39 mp_min = m;
40 else
41 mp_max = m;
42 }
43 while((mp_max-mp_min)!=1);
44 m = mp_max;
45
46 if (m<2)
47 m = 2;
48 else if (Ngrid_NormK<=m)
49 m = Ngrid_NormK - 2;
50 }
51
52 /****************************************************
53 Spline like interpolation
54 ****************************************************/
55
56 if (po==0){
57
58 if (m==1){
59 h2 = NormK[m] - NormK[m-1];
60 h3 = NormK[m+1] - NormK[m];
61
62 f2 = Spe_RF_Bessel[Gensi][GL][Mul][m-1];
63 f3 = Spe_RF_Bessel[Gensi][GL][Mul][m];
64 f4 = Spe_RF_Bessel[Gensi][GL][Mul][m+1];
65
66 h1 = -(h2+h3);
67 f1 = f4;
68 }
69 else if (m==(Ngrid_NormK-1)){
70 h1 = NormK[m-1] - NormK[m-2];
71 h2 = NormK[m] - NormK[m-1];
72
73 f1 = Spe_RF_Bessel[Gensi][GL][Mul][m-2];
74 f2 = Spe_RF_Bessel[Gensi][GL][Mul][m-1];
75 f3 = Spe_RF_Bessel[Gensi][GL][Mul][m];
76
77 h3 = -(h1+h2);
78 f4 = f1;
79 }
80 else{
81 h1 = NormK[m-1] - NormK[m-2];
82 h2 = NormK[m] - NormK[m-1];
83 h3 = NormK[m+1] - NormK[m];
84
85 f1 = Spe_RF_Bessel[Gensi][GL][Mul][m-2];
86 f2 = Spe_RF_Bessel[Gensi][GL][Mul][m-1];
87 f3 = Spe_RF_Bessel[Gensi][GL][Mul][m];
88 f4 = Spe_RF_Bessel[Gensi][GL][Mul][m+1];
89 }
90
91 /****************************************************
92 Calculate the value at R
93 ****************************************************/
94
95 g1 = ((f3-f2)*h1/h2 + (f2-f1)*h2/h1)/(h1+h2);
96 g2 = ((f4-f3)*h2/h3 + (f3-f2)*h3/h2)/(h2+h3);
97
98 x1 = R - NormK[m-1];
99 x2 = R - NormK[m];
100 y1 = x1/h2;
101 y2 = x2/h2;
102
103 f = y2*y2*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2)
104 + y1*y1*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1);
105
106 result = f;
107 }
108
109 return result;
110 }
111