Line data Source code
1 : MODULE m_grdchlh
2 : use m_juDFT
3 : ! -----------------------------------------------------------------
4 : ! input: rv present: exponential mesh. otherwise dx interval mesh.
5 : ! ro: charge or quantity to be derivated.
6 : ! evaluates d(ro)/dr,d{d(ro)/dr}/dr.
7 : ! drr=d(ro)/dr, ddrr=d(drr)/dr.
8 : ! coded by t.asada. 1996.
9 : ! ------------------------------------------------------------------
10 : CONTAINS
11 :
12 43813 : SUBROUTINE grdchlh(dx,ro, drr,ddrr,rv,order)
13 :
14 : USE m_constants
15 :
16 : IMPLICIT NONE
17 :
18 : REAL, INTENT (IN) :: dx
19 : REAL, INTENT (IN) :: ro(:)
20 : REAL, INTENT (OUT) :: drr(:),ddrr(:)
21 :
22 : real,intent(in),optional :: rv(:)
23 : INTEGER, INTENT (IN),OPTIONAL :: order
24 :
25 : REAL drx,drx0,drx1,drx2,drx3,drxx,drxx0,drxx1,drxx2,drxx3
26 : INTEGER i,i1,i2,i3,i4,i5,i6,j,nred,ndvgrd
27 :
28 : !call timestart("grdchlh")
29 :
30 43813 : ndvgrd=6
31 43813 : if (present(order)) ndvgrd=order
32 :
33 15299983 : drr = 0.0
34 15299983 : ddrr = 0.0
35 :
36 43813 : IF (size(ro) < 4) RETURN
37 :
38 43813 : IF (ndvgrd < 3 .OR. ndvgrd>6) THEN
39 0 : CALL juDFT_error("ndvgrd<3 .or. ndvgrd>6",calledby="grdchlh")
40 : ENDIF
41 : 126 FORMAT (/,' ndvgrd should be ge.4 .or. le.6. ndvgrd=',i3)
42 :
43 43813 : i1 = 1
44 43813 : i2 = 2
45 43813 : i3 = 3
46 43813 : i4 = 4
47 43813 : i5 = 5
48 43813 : i6 = 6
49 :
50 43813 : IF (ndvgrd==3) THEN
51 :
52 0 : drx1 = f131(ro(i1),ro(i2),ro(i3),dx)
53 0 : drxx1 = f231(ro(i1),ro(i2),ro(i3),dx)
54 :
55 43813 : ELSEIF (ndvgrd==4) THEN
56 :
57 0 : drx1 = f141(ro(i1),ro(i2),ro(i3),ro(i4),dx)
58 0 : drxx1 = f241(ro(i1),ro(i2),ro(i3),ro(i4),dx)
59 0 : drx2 = f142(ro(i1),ro(i2),ro(i3),ro(i4),dx)
60 0 : drxx2 = f242(ro(i1),ro(i2),ro(i3),ro(i4),dx)
61 :
62 43813 : ELSEIF (ndvgrd==5) THEN
63 :
64 0 : drx1 = f151(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),dx)
65 0 : drxx1 = f251(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),dx)
66 0 : drx2 = f152(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),dx)
67 0 : drxx2 = f252(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),dx)
68 :
69 43813 : ELSEIF (ndvgrd==6) THEN
70 :
71 43813 : drx1 = f161(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),ro(i6),dx)
72 43813 : drxx1 = f261(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),ro(i6),dx)
73 43813 : drx2 = f162(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),ro(i6),dx)
74 43813 : drxx2 = f262(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),ro(i6),dx)
75 43813 : drx3 = f163(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),ro(i6),dx)
76 43813 : drxx3 = f263(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),ro(i6),dx)
77 :
78 : ENDIF
79 :
80 43813 : IF (present(rv)) THEN
81 17555 : drr(i1) = drx1/rv(i1)
82 17555 : ddrr(i1) = (drxx1-drx1)/rv(i1)**2
83 : ELSE
84 26258 : drr(i1) = drx1
85 26258 : ddrr(i1) = drxx1
86 : ENDIF
87 :
88 43813 : IF (ndvgrd>3) THEN
89 :
90 43813 : IF (present(rv)) THEN
91 17555 : drr(i2) = drx2/rv(i2)
92 17555 : ddrr(i2) = (drxx2-drx2)/rv(i2)**2
93 : ELSE
94 26258 : drr(i2) = drx2
95 26258 : ddrr(i2) = drxx2
96 : ENDIF
97 :
98 43813 : IF (ndvgrd==6) THEN
99 43813 : IF (present(rv)) THEN
100 17555 : drr(i3) = drx3/rv(i3)
101 17555 : ddrr(i3) = (drxx3-drx3)/rv(i3)**2
102 : ELSE
103 26258 : drr(i3) = drx3
104 26258 : ddrr(i3) = drxx3
105 : ENDIF
106 : ENDIF
107 :
108 : ENDIF
109 :
110 43813 : nred = REAL(ndvgrd)/2 + 0.1
111 :
112 43813 : IF (size(ro)-nred.LE.1) THEN
113 0 : WRITE(oUnit,fmt='(/'' size(ro)-nred < 1. size(ro),nred='',3i4)') size(ro),nred
114 0 : CALL juDFT_error("size(ro)-nred.le.1",calledby="grdchlh")
115 : ENDIF
116 :
117 14828193 : DO j = nred + 1,size(ro) - nred
118 :
119 14784380 : IF (ndvgrd==3) THEN
120 :
121 0 : drx = f132(ro(j-1),ro(j),ro(j+1),dx)
122 0 : drxx = f232(ro(j-1),ro(j),ro(j+1),dx)
123 :
124 14784380 : ELSEIF (ndvgrd==4) THEN
125 :
126 0 : drx = f142(ro(j-1),ro(j),ro(j+1),ro(j+2),dx)
127 0 : drxx = f242(ro(j-1),ro(j),ro(j+1),ro(j+2),dx)
128 :
129 14784380 : ELSEIF (ndvgrd==5) THEN
130 :
131 0 : drx = f153(ro(j-2),ro(j-1),ro(j),ro(j+1),ro(j+2),dx)
132 0 : drxx = f253(ro(j-2),ro(j-1),ro(j),ro(j+1),ro(j+2),dx)
133 :
134 14784380 : ELSEIF (ndvgrd==6) THEN
135 :
136 14784380 : drx = f164(ro(j-3),ro(j-2),ro(j-1),ro(j),ro(j+1),ro(j+2), dx)
137 14784380 : drxx = f264(ro(j-3),ro(j-2),ro(j-1),ro(j),ro(j+1),ro(j+2), dx)
138 :
139 : ENDIF
140 :
141 14828193 : IF (present(rv)) THEN
142 12311328 : drr(j) = drx/rv(j)
143 12311328 : ddrr(j) = (drxx-drx)/rv(j)**2
144 : ELSE
145 2473052 : drr(j) = drx
146 2473052 : ddrr(j) = drxx
147 : ENDIF
148 :
149 : ENDDO ! j
150 :
151 43813 : IF (ndvgrd==3) THEN
152 :
153 0 : drx0 = f133(ro(size(ro)-2),ro(size(ro)-1),ro(size(ro)),dx)
154 0 : drxx0 = f233(ro(size(ro)-2),ro(size(ro)-1),ro(size(ro)),dx)
155 :
156 43813 : ELSEIF (ndvgrd==4) THEN
157 :
158 0 : drx1 = f143(ro(size(ro)-3),ro(size(ro)-2),ro(size(ro)-1),ro(size(ro)),dx)
159 0 : drxx1 = f243(ro(size(ro)-3),ro(size(ro)-2),ro(size(ro)-1),ro(size(ro)),dx)
160 0 : drx0 = f144(ro(size(ro)-3),ro(size(ro)-2),ro(size(ro)-1),ro(size(ro)),dx)
161 0 : drxx0 = f244(ro(size(ro)-3),ro(size(ro)-2),ro(size(ro)-1),ro(size(ro)),dx)
162 :
163 43813 : ELSEIF (ndvgrd==5) THEN
164 :
165 0 : drx1 = f154(ro(size(ro)-4),ro(size(ro)-3),ro(size(ro)-2),ro(size(ro)-1),ro(size(ro)), dx)
166 0 : drxx1 = f254(ro(size(ro)-4),ro(size(ro)-3),ro(size(ro)-2),ro(size(ro)-1),ro(size(ro)), dx)
167 0 : drx0 = f155(ro(size(ro)-4),ro(size(ro)-3),ro(size(ro)-2),ro(size(ro)-1),ro(size(ro)), dx)
168 0 : drxx0 = f255(ro(size(ro)-4),ro(size(ro)-3),ro(size(ro)-2),ro(size(ro)-1),ro(size(ro)), dx)
169 :
170 43813 : ELSEIF (ndvgrd==6) THEN
171 :
172 43813 : drx2 = f164(ro(size(ro)-5),ro(size(ro)-4),ro(size(ro)-3),ro(size(ro)-2), ro(size(ro)-1),ro(size(ro)),dx)
173 43813 : drxx2 = f264(ro(size(ro)-5),ro(size(ro)-4),ro(size(ro)-3),ro(size(ro)-2), ro(size(ro)-1),ro(size(ro)),dx)
174 :
175 43813 : drx1 = f165(ro(size(ro)-5),ro(size(ro)-4),ro(size(ro)-3),ro(size(ro)-2), ro(size(ro)-1),ro(size(ro)),dx)
176 43813 : drxx1 = f265(ro(size(ro)-5),ro(size(ro)-4),ro(size(ro)-3),ro(size(ro)-2), ro(size(ro)-1),ro(size(ro)),dx)
177 :
178 43813 : drx0 = f166(ro(size(ro)-5),ro(size(ro)-4),ro(size(ro)-3),ro(size(ro)-2), ro(size(ro)-1),ro(size(ro)),dx)
179 43813 : drxx0 = f266(ro(size(ro)-5),ro(size(ro)-4),ro(size(ro)-3),ro(size(ro)-2), ro(size(ro)-1),ro(size(ro)),dx)
180 :
181 : ENDIF
182 :
183 43813 : IF (ndvgrd>3) THEN
184 :
185 43813 : IF (ndvgrd==6) THEN
186 43813 : IF (present(rv)) THEN
187 17555 : drr(size(ro)-2) = drx2/rv(size(ro)-2)
188 17555 : ddrr(size(ro)-2) = (drxx2-drx2)/rv(size(ro)-2)**2
189 : ELSE
190 26258 : drr(size(ro)-2) = drx2
191 26258 : ddrr(size(ro)-2) = drxx2
192 : ENDIF
193 : ENDIF
194 :
195 43813 : IF (present(rv)) THEN
196 17555 : drr(size(ro)-1) = drx1/rv(size(ro)-1)
197 17555 : ddrr(size(ro)-1) = (drxx1-drx1)/rv(size(ro)-1)**2
198 : ELSE
199 26258 : drr(size(ro)-1) = drx1
200 26258 : ddrr(size(ro)-1) = drxx1
201 : ENDIF
202 :
203 : ENDIF
204 :
205 43813 : IF (present(rv)) THEN
206 17555 : drr(size(ro)) = drx0/rv(size(ro))
207 17555 : ddrr(size(ro)) = (drxx0-drx0)/rv(size(ro))**2
208 : ELSE
209 26258 : drr(size(ro)) = drx0
210 26258 : ddrr(size(ro)) = drxx0
211 : ENDIF
212 :
213 : !call timestop("grdchlh")
214 : END SUBROUTINE grdchlh
215 : !--------------------------------------------------------------------
216 : ! Functions: 3 point formula :
217 : !
218 0 : REAL FUNCTION f131(f0,f1,f2,d)
219 : REAL f0,f1,f2,d
220 0 : f131 = (-3*f0+4*f1-f2)/ (2*d)
221 0 : END FUNCTION f131
222 0 : REAL FUNCTION f132(g1,f0,f1,d)
223 : REAL g1,f0,f1,d
224 0 : f132 = (-1*g1-0*f0+f1)/ (2*d)
225 0 : END FUNCTION f132
226 0 : REAL FUNCTION f133(g2,g1,f0,d)
227 : REAL g2,g1,f0,d
228 0 : f133 = (g2-4*g1+3*f0)/ (2*d)
229 0 : END FUNCTION f133
230 : !
231 : !.....four point formula for the 1st deriv.
232 0 : REAL FUNCTION f141(f0,f1,f2,f3,d)
233 : REAL f0,f1,f2,f3,d
234 0 : f141 = (-11*f0+18*f1-9*f2+2*f3)/ (6*d)
235 0 : END FUNCTION f141
236 0 : REAL FUNCTION f142(g1,f0,f1,f2,d)
237 : REAL g1,f0,f1,f2,d
238 0 : f142 = (-2*g1-3*f0+6*f1-f2)/ (6*d)
239 0 : END FUNCTION f142
240 0 : REAL FUNCTION f143(g2,g1,f0,f1,d)
241 : REAL g2,g1,f0,f1,d
242 0 : f143 = (g2-6*g1+3*f0+2*f1)/ (6*d)
243 0 : END FUNCTION f143
244 0 : REAL FUNCTION f144(g3,g2,g1,f0,d)
245 : REAL g3,g2,g1,f0,d
246 0 : f144 = (-2*g3+9*g2-18*g1+11*f0)/ (6*d)
247 0 : END FUNCTION f144
248 : !
249 : !.....five point formula for the 1st deriv.
250 0 : REAL FUNCTION f151(f0,f1,f2,f3,f4,d)
251 : REAL f0,f1,f2,f3,f4,d
252 0 : f151 = (-50*f0+96*f1-72*f2+32*f3-6*f4)/ (24*d)
253 0 : END FUNCTION f151
254 0 : REAL FUNCTION f152(g1,f0,f1,f2,f3,d)
255 : REAL g1,f0,f1,f2,f3,d
256 0 : f152 = (-6*g1-20*f0+36*f1-12*f2+2*f3)/ (24*d)
257 0 : END FUNCTION f152
258 0 : REAL FUNCTION f153(g2,g1,f0,f1,f2,d)
259 : REAL g2,g1,f0,f1,f2,d
260 0 : f153 = (2*g2-16*g1-0*f0+16*f1-2*f2)/ (24*d)
261 0 : END FUNCTION f153
262 0 : REAL FUNCTION f154(g3,g2,g1,f0,f1,d)
263 : REAL g3,g2,g1,f0,f1,d
264 0 : f154 = (-2*g3+12*g2-36*g1+20*f0+6*f1)/ (24*d)
265 0 : END FUNCTION f154
266 0 : REAL FUNCTION f155(g4,g3,g2,g1,f0,d)
267 : REAL g4,g3,g2,g1,f0,d
268 0 : f155 = (6*g4-32*g3+72*g2-96*g1+50*f0)/ (24*d)
269 0 : END FUNCTION f155
270 : !
271 : !.....six point formula for the 1st deriv.
272 43813 : REAL FUNCTION f161(f0,f1,f2,f3,f4,f5,d)
273 : REAL f0,f1,f2,f3,f4,f5,d
274 43813 : f161 = (-274*f0+600*f1-600*f2+400*f3-150*f4+24*f5)/ (120*d)
275 0 : END FUNCTION f161
276 43813 : REAL FUNCTION f162(g1,f0,f1,f2,f3,f4,d)
277 : REAL g1,f0,f1,f2,f3,f4,d
278 43813 : f162 = (-24*g1-130*f0+240*f1-120*f2+40*f3-6*f4)/ (120*d)
279 0 : END FUNCTION f162
280 43813 : REAL FUNCTION f163(g2,g1,f0,f1,f2,f3,d)
281 : REAL g2,g1,f0,f1,f2,f3,d
282 43813 : f163 = (6*g2-60*g1-40*f0+120*f1-30*f2+4*f3)/ (120*d)
283 0 : END FUNCTION f163
284 14828193 : REAL FUNCTION f164(g3,g2,g1,f0,f1,f2,d)
285 : REAL g3,g2,g1,f0,f1,f2,d
286 14828193 : f164 = (-4*g3+30*g2-120*g1+40*f0+60*f1-6*f2)/ (120*d)
287 0 : END FUNCTION f164
288 43813 : REAL FUNCTION f165(g4,g3,g2,g1,f0,f1,d)
289 : REAL g4,g3,g2,g1,f0,f1,d
290 43813 : f165 = (6*g4-40*g3+120*g2-240*g1+130*f0+24*f1)/ (120*d)
291 0 : END FUNCTION f165
292 43813 : REAL FUNCTION f166(g5,g4,g3,g2,g1,f0,d)
293 : REAL g5,g4,g3,g2,g1,f0,d
294 43813 : f166 = (-24*g5+150*g4-400*g3+600*g2-600*g1+274*f0)/ (120*d)
295 0 : END FUNCTION f166
296 : !
297 : !.....three point formula for the 2nd deriv.
298 0 : REAL FUNCTION f231(f0,f1,f2,d)
299 : REAL f0,f1,f2,d
300 0 : f231 = (f0-2*f1+f2)/ (d*d)
301 0 : END FUNCTION f231
302 0 : REAL FUNCTION f232(g1,f0,f1,d)
303 : REAL g1,f0,f1,d
304 0 : f232 = (g1-2*f0+f1)/ (d*d)
305 0 : END FUNCTION f232
306 0 : REAL FUNCTION f233(g2,g1,f0,d)
307 : REAL g2,g1,f0,d
308 0 : f233 = (g2-2*g1+f0)/ (d*d)
309 0 : END FUNCTION f233
310 : !
311 : !.....four point formula for the 2nd deriv.
312 0 : REAL FUNCTION f241(f0,f1,f2,f3,d)
313 : REAL f0,f1,f2,f3,d
314 0 : f241 = (6*f0-15*f1+12*f2-3*f3)/ (3*d*d)
315 0 : END FUNCTION f241
316 0 : REAL FUNCTION f242(g1,f0,f1,f2,d)
317 : REAL g1,f0,f1,f2,d
318 0 : f242 = (3*g1-6*f0+3*f1+0*f2)/ (3*d*d)
319 0 : END FUNCTION f242
320 0 : REAL FUNCTION f243(g2,g1,f0,f1,d)
321 : REAL g2,g1,f0,f1,d
322 0 : f243 = (0*g2+3*g1-6*f0+3*f1)/ (3*d*d)
323 0 : END FUNCTION f243
324 0 : REAL FUNCTION f244(g3,g2,g1,f0,d)
325 : REAL g3,g2,g1,f0,d
326 0 : f244 = (-3*g3+2*g2+15*g1+6*f0)/ (3*d*d)
327 0 : END FUNCTION f244
328 : !
329 : !.....five point formula for the 2nd deriv.
330 0 : REAL FUNCTION f251(f0,f1,f2,f3,f4,d)
331 : REAL f0,f1,f2,f3,f4,d
332 0 : f251 = (35*f0-104*f1+114*f2-56*f3+11*f4)/ (12*d*d)
333 0 : END FUNCTION f251
334 0 : REAL FUNCTION f252(g1,f0,f1,f2,f3,d)
335 : REAL g1,f0,f1,f2,f3,d
336 0 : f252 = (11*g1-20*f0+6*f1+4*f2-f3)/ (12*d*d)
337 0 : END FUNCTION f252
338 0 : REAL FUNCTION f253(g2,g1,f0,f1,f2,d)
339 : REAL g2,g1,f0,f1,f2,d
340 0 : f253 = (-g2+16*g1-30*f0+16*f1-f2)/ (12*d*d)
341 0 : END FUNCTION f253
342 0 : REAL FUNCTION f254(g3,g2,g1,f0,f1,d)
343 : REAL g3,g2,g1,f0,f1,d
344 0 : f254 = (-g3+4*g2+6*g1-20*f0+11*f1)/ (12*d*d)
345 0 : END FUNCTION f254
346 0 : REAL FUNCTION f255(g4,g3,g2,g1,f0,d)
347 : REAL g4,g3,g2,g1,f0,d
348 0 : f255 = (11*g4-56*g3+114*g2-104*g1+35*f0)/ (12*d*d)
349 0 : END FUNCTION f255
350 : !
351 : !.....six point formula for the 2nd deriv.
352 43813 : REAL FUNCTION f261(f0,f1,f2,f3,f4,f5,d)
353 : REAL f0,f1,f2,f3,f4,f5,d
354 43813 : f261 = (225*f0-770*f1+1070*f2-780*f3+305*f4-50*f5)/ (60*d*d)
355 0 : END FUNCTION f261
356 43813 : REAL FUNCTION f262(g1,f0,f1,f2,f3,f4,d)
357 : REAL g1,f0,f1,f2,f3,f4,d
358 43813 : f262 = (50*g1-75*f0-20*f1+70*f2-30*f3+5*f4)/ (60*d*d)
359 0 : END FUNCTION f262
360 43813 : REAL FUNCTION f263(g2,g1,f0,f1,f2,f3,d)
361 : REAL g2,g1,f0,f1,f2,f3,d
362 43813 : f263 = (-5*g2+80*g1-150*f0+80*f1-5*f2+0*f3)/ (60*d*d)
363 0 : END FUNCTION f263
364 14828193 : REAL FUNCTION f264(g3,g2,g1,f0,f1,f2,d)
365 : REAL g3,g2,g1,f0,f1,f2,d
366 14784380 : f264 = (0*g3-5*g2+80*g1-150*f0+80*f1-5*f2)/ (60*d*d)
367 0 : END FUNCTION f264
368 43813 : REAL FUNCTION f265(g4,g3,g2,g1,f0,f1,d)
369 : REAL g4,g3,g2,g1,f0,f1,d
370 43813 : f265 = (5*g4-30*g3+70*g2-20*g1-75*f0+50*f1)/ (60*d*d)
371 0 : END FUNCTION f265
372 43813 : REAL FUNCTION f266(g5,g4,g3,g2,g1,f0,d)
373 : REAL g5,g4,g3,g2,g1,f0,d
374 43813 : f266 = (-50*g5+305*g4-780*g3+1070*g2-770*g1+225*f0)/ (60*d*d)
375 0 : END FUNCTION f266
376 : !--------------------------------------------------------------------
377 : END MODULE m_grdchlh
|
