1 /**********************************************************************
2 Dr_VH_AtomF.c:
3
4 Dr_VH_AtomF.c is a subroutine to calculate the radial derivative
5 of the Hartree potential of a free atom specified by "spe".
6
7 Log of Dr_VH_AtomF.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
17 #ifdef MAX
18 #undef MAX
19 #endif
20 #define MAX(a,b) ((a)>(b))? (a):(b)
21
22 #ifdef MIN
23 #undef MIN
24 #endif
25 #define MIN(a,b) ((a)<(b))? (a):(b)
26
Dr_VH_AtomF(int spe,int N,double x,double r,double * xv,double * rv,double * yv)27 double Dr_VH_AtomF(int spe, int N, double x, double r, double *xv, double *rv, double *yv)
28 {
29 int i;
30 double t,dt,y;
31 double xmin,xmax,tmp;
32
33 xmin = xv[0];
34 xmax = xv[N-1];
35
36 if (xmax<=x){
37 return -Spe_Core_Charge[spe]/r/r;
38 }
39 else if (r<Spe_VPS_RV[spe][0]){
40
41 int m;
42 double rm,h1,h2,h3,f1,f2,f3,f4,a,b;
43 double g1,g2,x1,x2,y1,y2,y12,y22,f,df;
44
45 m = 4;
46 rm = rv[m];
47
48 h1 = rv[m-1] - rv[m-2];
49 h2 = rv[m] - rv[m-1];
50 h3 = rv[m+1] - rv[m];
51
52 f1 = yv[m-2];
53 f2 = yv[m-1];
54 f3 = yv[m];
55 f4 = yv[m+1];
56
57 g1 = ((f3-f2)*h1/h2 + (f2-f1)*h2/h1)/(h1+h2);
58 g2 = ((f4-f3)*h2/h3 + (f3-f2)*h3/h2)/(h2+h3);
59
60 x1 = rm - rv[m-1];
61 x2 = rm - rv[m];
62 y1 = x1/h2;
63 y2 = x2/h2;
64 y12 = y1*y1;
65 y22 = y2*y2;
66
67 f = y22*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2)
68 + y12*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1);
69
70 df = 2.0*y2/h2*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2)
71 + y22*(2.0*f2 + h2*g1)/h2
72 + 2.0*y1/h2*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1)
73 - y12*(2.0*f3 - h2*g2)/h2;
74
75 a = 0.5*df/rm;
76 b = f - a*rm*rm;
77 return 2.0*a*r;
78
79 }
80 else{
81
82 x = MAX(x,xmin);
83 tmp = ((double)N-1.0)/(xmax-xmin);
84 t = (x-xmin)*tmp;
85 i = floor(t);
86 dt = t - (double)i;
87
88 return 0.5*(( 3.0*(yv[i+3]-yv[i]-3.0*(yv[i+2]-yv[i+1]))*dt
89 +2.0*(-yv[i+3]+4.0*yv[i+2]-5.0*yv[i+1]+2.0*yv[i]))*dt
90 +(yv[i+2]-yv[i]))*tmp/r;
91 }
92 }
93
94
95
96
97
98
99