(git:374b731)
Loading...
Searching...
No Matches
t_c_g0.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8!--------------------------------------------------------------------------------------------------!
9! Copyright (c) 2008, 2009, Joost VandeVondele and Manuel Guidon !
10! All rights reserved. !
11! !
12! Redistribution and use in source and binary forms, with or without !
13! modification, are permitted provided that the following conditions are met: !
14! * Redistributions of source code must retain the above copyright !
15! notice, this list of conditions and the following disclaimer. !
16! * Redistributions in binary form must reproduce the above copyright !
17! notice, this list of conditions and the following disclaimer in the !
18! documentation and/or other materials provided with the distribution. !
19! !
20! THIS SOFTWARE IS PROVIDED BY Joost VandeVondele and Manuel Guidon AS IS AND ANY !
21! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED !
22! WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE !
23! DISCLAIMED. IN NO EVENT SHALL Joost VandeVondele or Manuel Guidon BE LIABLE FOR ANY !
24! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES !
25! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; !
26! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND !
27! ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT !
28! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS !
29! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. !
30!--------------------------------------------------------------------------------------------------!
31
32! **************************************************************************************************
33!> \brief This module computes the basic integrals for the truncated coulomb operator
34!>
35!> res(1) =G_0(R,T)= ((2*erf(sqrt(t))+erf(R-sqrt(t))-erf(R+sqrt(t)))/sqrt(t))
36!>
37!> and up to 21 derivatives with respect to T
38!>
39!> res(n+1)=(-1)**n d^n/dT^n G_0(R,T)
40!>
41!> The function is only computed for values of R,T which fulfil
42!>
43!> R**2 - 11.0_dp*R + 0.0_dp < T < R**2 + 11.0_dp*R + 50.0_dp where R>=0 T>=0
44!>
45!> for T larger than the upper bound, 0 is returned
46!> (which is accurate at least up to 1.0E-16)
47!> while for T smaller than the lower bound, the caller is instructed
48!> to use the conventional gamma function instead
49!> (i.e. the limit of above expression for R to Infinity)
50!>
51!> \author Joost VandeVondele and Manuel Guidon
52!> \par History
53!> Nov 2008, 2009 Joost VandeVondele and Manuel Guidon
54!> May 2019 A. Bussy: Added a get_maxl_init function to get current status of nderiv_init and
55!> moved the file to common (made it accessible from aobasis, same place as gamma.F).
56! **************************************************************************************************
57MODULE t_c_g0
58 USE kinds, ONLY: dp
59 USE message_passing, ONLY: mp_comm_type
60#include "../base/base_uses.f90"
61
62 IMPLICIT NONE
63 PRIVATE
64
65 PUBLIC :: t_c_g0_n, init, free_c0, get_lmax_init
66 REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE, SAVE :: c0
67 INTEGER, PARAMETER :: degree = 13
68 REAL(KIND=dp), PARAMETER :: target_error = 0.100000e-08
69 INTEGER, PARAMETER :: nderiv_max = 21
70 INTEGER, SAVE :: nderiv_init = -1
71 INTEGER, SAVE :: patches = -1
72
73CONTAINS
74
75! **************************************************************************************************
76!> \brief ...
77!> \param RES ...
78!> \param use_gamma ...
79!> \param R ...
80!> \param T ...
81!> \param NDERIV ...
82! **************************************************************************************************
83 SUBROUTINE t_c_g0_n(RES, use_gamma, R, T, NDERIV)
84 REAL(kind=dp), INTENT(OUT) :: res(*)
85 LOGICAL, INTENT(OUT) :: use_gamma
86 REAL(kind=dp), INTENT(IN) :: r, t
87 INTEGER, INTENT(IN) :: nderiv
88
89 REAL(kind=dp) :: lower, tg1, tg2, upper, x1, x2
90
91 use_gamma = .false.
92 upper = r**2 + 11.0_dp*r + 50.0_dp
93 lower = r**2 - 11.0_dp*r + 0.0_dp
94 IF (t > upper) THEN
95 res(1:nderiv + 1) = 0.0_dp
96 RETURN
97 END IF
98 IF (r <= 11.0_dp) THEN
99 x2 = r/11.0_dp
100 upper = r**2 + 11.0_dp*r + 50.0_dp
101 lower = 0.0_dp
102 x1 = (t - lower)/(upper - lower)
103 IF (x1 <= 0.500000000000000000e+00_dp) THEN
104 IF (x2 <= 0.500000000000000000e+00_dp) THEN
105 IF (x2 <= 0.250000000000000000e+00_dp) THEN
106 IF (x2 <= 0.125000000000000000e+00_dp) THEN
107 IF (x1 <= 0.250000000000000000e+00_dp) THEN
108 IF (x2 <= 0.625000000000000000e-01_dp) THEN
109 IF (x1 <= 0.125000000000000000e+00_dp) THEN
110 IF (x2 <= 0.312500000000000000e-01_dp) THEN
111 IF (x1 <= 0.625000000000000000e-01_dp) THEN
112 IF (x2 <= 0.156250000000000000e-01_dp) THEN
113 IF (x1 <= 0.312500000000000000e-01_dp) THEN
114 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
115 tg2 = (2*x2 - 0.156250000000000000e-01_dp)*0.640000000000000000e+02_dp
116 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 1))
117 ELSE
118 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
119 tg2 = (2*x2 - 0.156250000000000000e-01_dp)*0.640000000000000000e+02_dp
120 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 2))
121 END IF
122 ELSE
123 IF (x1 <= 0.312500000000000000e-01_dp) THEN
124 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
125 tg2 = (2*x2 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
126 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 3))
127 ELSE
128 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
129 tg2 = (2*x2 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
130 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 4))
131 END IF
132 END IF
133 ELSE
134 IF (x2 <= 0.156250000000000000e-01_dp) THEN
135 tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
136 tg2 = (2*x2 - 0.156250000000000000e-01_dp)*0.640000000000000000e+02_dp
137 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 5))
138 ELSE
139 tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
140 tg2 = (2*x2 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
141 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 6))
142 END IF
143 END IF
144 ELSE
145 IF (x1 <= 0.625000000000000000e-01_dp) THEN
146 IF (x2 <= 0.468750000000000000e-01_dp) THEN
147 IF (x1 <= 0.312500000000000000e-01_dp) THEN
148 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
149 tg2 = (2*x2 - 0.781250000000000000e-01_dp)*0.640000000000000000e+02_dp
150 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 7))
151 ELSE
152 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
153 tg2 = (2*x2 - 0.781250000000000000e-01_dp)*0.640000000000000000e+02_dp
154 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 8))
155 END IF
156 ELSE
157 IF (x1 <= 0.312500000000000000e-01_dp) THEN
158 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
159 tg2 = (2*x2 - 0.109375000000000000e+00_dp)*0.640000000000000000e+02_dp
160 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 9))
161 ELSE
162 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
163 tg2 = (2*x2 - 0.109375000000000000e+00_dp)*0.640000000000000000e+02_dp
164 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 10))
165 END IF
166 END IF
167 ELSE
168 IF (x2 <= 0.468750000000000000e-01_dp) THEN
169 tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
170 tg2 = (2*x2 - 0.781250000000000000e-01_dp)*0.640000000000000000e+02_dp
171 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 11))
172 ELSE
173 tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
174 tg2 = (2*x2 - 0.109375000000000000e+00_dp)*0.640000000000000000e+02_dp
175 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 12))
176 END IF
177 END IF
178 END IF
179 ELSE
180 IF (x2 <= 0.312500000000000000e-01_dp) THEN
181 IF (x1 <= 0.187500000000000000e+00_dp) THEN
182 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
183 tg2 = (2*x2 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
184 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 13))
185 ELSE
186 tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
187 tg2 = (2*x2 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
188 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 14))
189 END IF
190 ELSE
191 IF (x1 <= 0.187500000000000000e+00_dp) THEN
192 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
193 tg2 = (2*x2 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
194 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 15))
195 ELSE
196 tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
197 tg2 = (2*x2 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
198 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 16))
199 END IF
200 END IF
201 END IF
202 ELSE
203 IF (x1 <= 0.125000000000000000e+00_dp) THEN
204 IF (x2 <= 0.937500000000000000e-01_dp) THEN
205 IF (x1 <= 0.625000000000000000e-01_dp) THEN
206 IF (x2 <= 0.781250000000000000e-01_dp) THEN
207 IF (x1 <= 0.312500000000000000e-01_dp) THEN
208 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
209 tg2 = (2*x2 - 0.140625000000000000e+00_dp)*0.640000000000000000e+02_dp
210 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 17))
211 ELSE
212 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
213 tg2 = (2*x2 - 0.140625000000000000e+00_dp)*0.640000000000000000e+02_dp
214 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 18))
215 END IF
216 ELSE
217 IF (x1 <= 0.312500000000000000e-01_dp) THEN
218 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
219 tg2 = (2*x2 - 0.171875000000000000e+00_dp)*0.640000000000000000e+02_dp
220 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 19))
221 ELSE
222 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
223 tg2 = (2*x2 - 0.171875000000000000e+00_dp)*0.640000000000000000e+02_dp
224 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 20))
225 END IF
226 END IF
227 ELSE
228 IF (x2 <= 0.781250000000000000e-01_dp) THEN
229 tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
230 tg2 = (2*x2 - 0.140625000000000000e+00_dp)*0.640000000000000000e+02_dp
231 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 21))
232 ELSE
233 tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
234 tg2 = (2*x2 - 0.171875000000000000e+00_dp)*0.640000000000000000e+02_dp
235 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 22))
236 END IF
237 END IF
238 ELSE
239 IF (x1 <= 0.625000000000000000e-01_dp) THEN
240 IF (x2 <= 0.109375000000000000e+00_dp) THEN
241 IF (x1 <= 0.312500000000000000e-01_dp) THEN
242 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
243 tg2 = (2*x2 - 0.203125000000000000e+00_dp)*0.640000000000000000e+02_dp
244 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 23))
245 ELSE
246 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
247 tg2 = (2*x2 - 0.203125000000000000e+00_dp)*0.640000000000000000e+02_dp
248 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 24))
249 END IF
250 ELSE
251 IF (x1 <= 0.312500000000000000e-01_dp) THEN
252 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
253 tg2 = (2*x2 - 0.234375000000000000e+00_dp)*0.640000000000000000e+02_dp
254 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 25))
255 ELSE
256 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
257 tg2 = (2*x2 - 0.234375000000000000e+00_dp)*0.640000000000000000e+02_dp
258 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 26))
259 END IF
260 END IF
261 ELSE
262 IF (x2 <= 0.109375000000000000e+00_dp) THEN
263 tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
264 tg2 = (2*x2 - 0.203125000000000000e+00_dp)*0.640000000000000000e+02_dp
265 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 27))
266 ELSE
267 tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
268 tg2 = (2*x2 - 0.234375000000000000e+00_dp)*0.640000000000000000e+02_dp
269 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 28))
270 END IF
271 END IF
272 END IF
273 ELSE
274 IF (x1 <= 0.187500000000000000e+00_dp) THEN
275 IF (x2 <= 0.937500000000000000e-01_dp) THEN
276 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
277 tg2 = (2*x2 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
278 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 29))
279 ELSE
280 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
281 tg2 = (2*x2 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
282 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 30))
283 END IF
284 ELSE
285 tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
286 tg2 = (2*x2 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
287 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 31))
288 END IF
289 END IF
290 END IF
291 ELSE
292 IF (x1 <= 0.375000000000000000e+00_dp) THEN
293 tg1 = (2*x1 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
294 tg2 = (2*x2 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
295 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 32))
296 ELSE
297 tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
298 tg2 = (2*x2 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
299 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 33))
300 END IF
301 END IF
302 ELSE
303 IF (x1 <= 0.250000000000000000e+00_dp) THEN
304 IF (x2 <= 0.187500000000000000e+00_dp) THEN
305 IF (x1 <= 0.125000000000000000e+00_dp) THEN
306 IF (x2 <= 0.156250000000000000e+00_dp) THEN
307 IF (x1 <= 0.625000000000000000e-01_dp) THEN
308 IF (x1 <= 0.312500000000000000e-01_dp) THEN
309 IF (x2 <= 0.140625000000000000e+00_dp) THEN
310 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
311 tg2 = (2*x2 - 0.265625000000000000e+00_dp)*0.640000000000000000e+02_dp
312 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 34))
313 ELSE
314 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
315 tg2 = (2*x2 - 0.296875000000000000e+00_dp)*0.640000000000000000e+02_dp
316 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 35))
317 END IF
318 ELSE
319 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
320 tg2 = (2*x2 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
321 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 36))
322 END IF
323 ELSE
324 IF (x1 <= 0.937500000000000000e-01_dp) THEN
325 tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
326 tg2 = (2*x2 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
327 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 37))
328 ELSE
329 tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
330 tg2 = (2*x2 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
331 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 38))
332 END IF
333 END IF
334 ELSE
335 IF (x1 <= 0.625000000000000000e-01_dp) THEN
336 IF (x1 <= 0.312500000000000000e-01_dp) THEN
337 IF (x2 <= 0.171875000000000000e+00_dp) THEN
338 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
339 tg2 = (2*x2 - 0.328125000000000000e+00_dp)*0.640000000000000000e+02_dp
340 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 39))
341 ELSE
342 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
343 tg2 = (2*x2 - 0.359375000000000000e+00_dp)*0.640000000000000000e+02_dp
344 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 40))
345 END IF
346 ELSE
347 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
348 tg2 = (2*x2 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
349 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 41))
350 END IF
351 ELSE
352 tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
353 tg2 = (2*x2 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
354 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 42))
355 END IF
356 END IF
357 ELSE
358 IF (x1 <= 0.187500000000000000e+00_dp) THEN
359 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
360 tg2 = (2*x2 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
361 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 43))
362 ELSE
363 tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
364 tg2 = (2*x2 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
365 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 44))
366 END IF
367 END IF
368 ELSE
369 IF (x1 <= 0.125000000000000000e+00_dp) THEN
370 IF (x1 <= 0.625000000000000000e-01_dp) THEN
371 IF (x2 <= 0.218750000000000000e+00_dp) THEN
372 IF (x1 <= 0.312500000000000000e-01_dp) THEN
373 IF (x2 <= 0.203125000000000000e+00_dp) THEN
374 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
375 tg2 = (2*x2 - 0.390625000000000000e+00_dp)*0.640000000000000000e+02_dp
376 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 45))
377 ELSE
378 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
379 tg2 = (2*x2 - 0.421875000000000000e+00_dp)*0.640000000000000000e+02_dp
380 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 46))
381 END IF
382 ELSE
383 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
384 tg2 = (2*x2 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
385 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 47))
386 END IF
387 ELSE
388 IF (x1 <= 0.312500000000000000e-01_dp) THEN
389 IF (x2 <= 0.234375000000000000e+00_dp) THEN
390 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
391 tg2 = (2*x2 - 0.453125000000000000e+00_dp)*0.640000000000000000e+02_dp
392 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 48))
393 ELSE
394 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
395 tg2 = (2*x2 - 0.484375000000000000e+00_dp)*0.640000000000000000e+02_dp
396 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 49))
397 END IF
398 ELSE
399 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
400 tg2 = (2*x2 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
401 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 50))
402 END IF
403 END IF
404 ELSE
405 IF (x2 <= 0.218750000000000000e+00_dp) THEN
406 tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
407 tg2 = (2*x2 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
408 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 51))
409 ELSE
410 tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
411 tg2 = (2*x2 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
412 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 52))
413 END IF
414 END IF
415 ELSE
416 IF (x1 <= 0.187500000000000000e+00_dp) THEN
417 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
418 tg2 = (2*x2 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
419 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 53))
420 ELSE
421 tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
422 tg2 = (2*x2 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
423 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 54))
424 END IF
425 END IF
426 END IF
427 ELSE
428 IF (x1 <= 0.375000000000000000e+00_dp) THEN
429 IF (x1 <= 0.312500000000000000e+00_dp) THEN
430 tg1 = (2*x1 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
431 tg2 = (2*x2 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
432 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 55))
433 ELSE
434 tg1 = (2*x1 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
435 tg2 = (2*x2 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
436 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 56))
437 END IF
438 ELSE
439 tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
440 tg2 = (2*x2 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
441 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 57))
442 END IF
443 END IF
444 END IF
445 ELSE
446 IF (x1 <= 0.250000000000000000e+00_dp) THEN
447 IF (x1 <= 0.125000000000000000e+00_dp) THEN
448 IF (x1 <= 0.625000000000000000e-01_dp) THEN
449 IF (x2 <= 0.375000000000000000e+00_dp) THEN
450 IF (x2 <= 0.312500000000000000e+00_dp) THEN
451 IF (x1 <= 0.312500000000000000e-01_dp) THEN
452 IF (x2 <= 0.281250000000000000e+00_dp) THEN
453 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
454 tg2 = (2*x2 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
455 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 58))
456 ELSE
457 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
458 tg2 = (2*x2 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
459 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 59))
460 END IF
461 ELSE
462 IF (x2 <= 0.281250000000000000e+00_dp) THEN
463 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
464 tg2 = (2*x2 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
465 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 60))
466 ELSE
467 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
468 tg2 = (2*x2 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
469 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 61))
470 END IF
471 END IF
472 ELSE
473 IF (x1 <= 0.312500000000000000e-01_dp) THEN
474 IF (x2 <= 0.343750000000000000e+00_dp) THEN
475 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
476 tg2 = (2*x2 - 0.656250000000000000e+00_dp)*0.320000000000000000e+02_dp
477 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 62))
478 ELSE
479 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
480 tg2 = (2*x2 - 0.718750000000000000e+00_dp)*0.320000000000000000e+02_dp
481 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 63))
482 END IF
483 ELSE
484 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
485 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
486 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 64))
487 END IF
488 END IF
489 ELSE
490 IF (x1 <= 0.312500000000000000e-01_dp) THEN
491 IF (x2 <= 0.437500000000000000e+00_dp) THEN
492 IF (x1 <= 0.156250000000000000e-01_dp) THEN
493 tg1 = (2*x1 - 0.156250000000000000e-01_dp)*0.640000000000000000e+02_dp
494 tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
495 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 65))
496 ELSE
497 tg1 = (2*x1 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
498 tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
499 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 66))
500 END IF
501 ELSE
502 IF (x1 <= 0.156250000000000000e-01_dp) THEN
503 tg1 = (2*x1 - 0.156250000000000000e-01_dp)*0.640000000000000000e+02_dp
504 tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
505 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 67))
506 ELSE
507 tg1 = (2*x1 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
508 tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
509 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 68))
510 END IF
511 END IF
512 ELSE
513 IF (x2 <= 0.437500000000000000e+00_dp) THEN
514 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
515 tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
516 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 69))
517 ELSE
518 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
519 tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
520 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 70))
521 END IF
522 END IF
523 END IF
524 ELSE
525 IF (x2 <= 0.375000000000000000e+00_dp) THEN
526 IF (x2 <= 0.312500000000000000e+00_dp) THEN
527 IF (x1 <= 0.937500000000000000e-01_dp) THEN
528 tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
529 tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
530 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 71))
531 ELSE
532 tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
533 tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
534 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 72))
535 END IF
536 ELSE
537 IF (x1 <= 0.937500000000000000e-01_dp) THEN
538 tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
539 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
540 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 73))
541 ELSE
542 tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
543 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
544 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 74))
545 END IF
546 END IF
547 ELSE
548 IF (x1 <= 0.937500000000000000e-01_dp) THEN
549 IF (x2 <= 0.437500000000000000e+00_dp) THEN
550 tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
551 tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
552 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 75))
553 ELSE
554 tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
555 tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
556 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 76))
557 END IF
558 ELSE
559 IF (x2 <= 0.437500000000000000e+00_dp) THEN
560 tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
561 tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
562 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 77))
563 ELSE
564 tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
565 tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
566 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 78))
567 END IF
568 END IF
569 END IF
570 END IF
571 ELSE
572 IF (x2 <= 0.375000000000000000e+00_dp) THEN
573 IF (x1 <= 0.187500000000000000e+00_dp) THEN
574 IF (x2 <= 0.312500000000000000e+00_dp) THEN
575 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
576 tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
577 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 79))
578 ELSE
579 IF (x1 <= 0.156250000000000000e+00_dp) THEN
580 tg1 = (2*x1 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
581 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
582 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 80))
583 ELSE
584 tg1 = (2*x1 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
585 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
586 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 81))
587 END IF
588 END IF
589 ELSE
590 IF (x2 <= 0.312500000000000000e+00_dp) THEN
591 tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
592 tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
593 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 82))
594 ELSE
595 tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
596 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
597 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 83))
598 END IF
599 END IF
600 ELSE
601 IF (x1 <= 0.187500000000000000e+00_dp) THEN
602 IF (x2 <= 0.437500000000000000e+00_dp) THEN
603 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
604 tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
605 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 84))
606 ELSE
607 IF (x1 <= 0.156250000000000000e+00_dp) THEN
608 tg1 = (2*x1 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
609 tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
610 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 85))
611 ELSE
612 tg1 = (2*x1 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
613 tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
614 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 86))
615 END IF
616 END IF
617 ELSE
618 IF (x1 <= 0.218750000000000000e+00_dp) THEN
619 tg1 = (2*x1 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
620 tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
621 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 87))
622 ELSE
623 tg1 = (2*x1 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
624 tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
625 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 88))
626 END IF
627 END IF
628 END IF
629 END IF
630 ELSE
631 IF (x1 <= 0.375000000000000000e+00_dp) THEN
632 IF (x2 <= 0.375000000000000000e+00_dp) THEN
633 IF (x1 <= 0.312500000000000000e+00_dp) THEN
634 tg1 = (2*x1 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
635 tg2 = (2*x2 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
636 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 89))
637 ELSE
638 tg1 = (2*x1 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
639 tg2 = (2*x2 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
640 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 90))
641 END IF
642 ELSE
643 IF (x1 <= 0.312500000000000000e+00_dp) THEN
644 IF (x1 <= 0.281250000000000000e+00_dp) THEN
645 tg1 = (2*x1 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
646 tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
647 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 91))
648 ELSE
649 tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
650 tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
651 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 92))
652 END IF
653 ELSE
654 tg1 = (2*x1 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
655 tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
656 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 93))
657 END IF
658 END IF
659 ELSE
660 IF (x1 <= 0.437500000000000000e+00_dp) THEN
661 tg1 = (2*x1 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
662 tg2 = (2*x2 - 0.750000000000000000e+00_dp)*0.400000000000000000e+01_dp
663 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 94))
664 ELSE
665 tg1 = (2*x1 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
666 tg2 = (2*x2 - 0.750000000000000000e+00_dp)*0.400000000000000000e+01_dp
667 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 95))
668 END IF
669 END IF
670 END IF
671 END IF
672 ELSE
673 IF (x1 <= 0.250000000000000000e+00_dp) THEN
674 IF (x1 <= 0.125000000000000000e+00_dp) THEN
675 IF (x1 <= 0.625000000000000000e-01_dp) THEN
676 IF (x1 <= 0.312500000000000000e-01_dp) THEN
677 IF (x1 <= 0.156250000000000000e-01_dp) THEN
678 IF (x1 <= 0.781250000000000000e-02_dp) THEN
679 IF (x2 <= 0.750000000000000000e+00_dp) THEN
680 tg1 = (2*x1 - 0.781250000000000000e-02_dp)*0.128000000000000000e+03_dp
681 tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
682 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 96))
683 ELSE
684 tg1 = (2*x1 - 0.781250000000000000e-02_dp)*0.128000000000000000e+03_dp
685 tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
686 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 97))
687 END IF
688 ELSE
689 IF (x2 <= 0.750000000000000000e+00_dp) THEN
690 tg1 = (2*x1 - 0.234375000000000000e-01_dp)*0.128000000000000000e+03_dp
691 tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
692 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 98))
693 ELSE
694 tg1 = (2*x1 - 0.234375000000000000e-01_dp)*0.128000000000000000e+03_dp
695 tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
696 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 99))
697 END IF
698 END IF
699 ELSE
700 IF (x2 <= 0.750000000000000000e+00_dp) THEN
701 tg1 = (2*x1 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
702 tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
703 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 100))
704 ELSE
705 tg1 = (2*x1 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
706 tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
707 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 101))
708 END IF
709 END IF
710 ELSE
711 IF (x2 <= 0.750000000000000000e+00_dp) THEN
712 IF (x2 <= 0.625000000000000000e+00_dp) THEN
713 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
714 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
715 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 102))
716 ELSE
717 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
718 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
719 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 103))
720 END IF
721 ELSE
722 IF (x1 <= 0.468750000000000000e-01_dp) THEN
723 tg1 = (2*x1 - 0.781250000000000000e-01_dp)*0.640000000000000000e+02_dp
724 tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
725 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 104))
726 ELSE
727 tg1 = (2*x1 - 0.109375000000000000e+00_dp)*0.640000000000000000e+02_dp
728 tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
729 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 105))
730 END IF
731 END IF
732 END IF
733 ELSE
734 IF (x2 <= 0.750000000000000000e+00_dp) THEN
735 IF (x2 <= 0.625000000000000000e+00_dp) THEN
736 IF (x1 <= 0.937500000000000000e-01_dp) THEN
737 tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
738 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
739 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 106))
740 ELSE
741 tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
742 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
743 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 107))
744 END IF
745 ELSE
746 IF (x1 <= 0.937500000000000000e-01_dp) THEN
747 tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
748 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
749 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 108))
750 ELSE
751 tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
752 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
753 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 109))
754 END IF
755 END IF
756 ELSE
757 IF (x1 <= 0.937500000000000000e-01_dp) THEN
758 tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
759 tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
760 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 110))
761 ELSE
762 tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
763 tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
764 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 111))
765 END IF
766 END IF
767 END IF
768 ELSE
769 IF (x2 <= 0.750000000000000000e+00_dp) THEN
770 IF (x2 <= 0.625000000000000000e+00_dp) THEN
771 IF (x1 <= 0.187500000000000000e+00_dp) THEN
772 IF (x1 <= 0.156250000000000000e+00_dp) THEN
773 tg1 = (2*x1 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
774 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
775 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 112))
776 ELSE
777 tg1 = (2*x1 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
778 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
779 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 113))
780 END IF
781 ELSE
782 IF (x1 <= 0.218750000000000000e+00_dp) THEN
783 tg1 = (2*x1 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
784 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
785 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 114))
786 ELSE
787 tg1 = (2*x1 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
788 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
789 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 115))
790 END IF
791 END IF
792 ELSE
793 IF (x1 <= 0.187500000000000000e+00_dp) THEN
794 IF (x1 <= 0.156250000000000000e+00_dp) THEN
795 tg1 = (2*x1 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
796 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
797 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 116))
798 ELSE
799 tg1 = (2*x1 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
800 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
801 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 117))
802 END IF
803 ELSE
804 IF (x1 <= 0.218750000000000000e+00_dp) THEN
805 tg1 = (2*x1 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
806 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
807 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 118))
808 ELSE
809 tg1 = (2*x1 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
810 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
811 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 119))
812 END IF
813 END IF
814 END IF
815 ELSE
816 IF (x1 <= 0.187500000000000000e+00_dp) THEN
817 IF (x2 <= 0.875000000000000000e+00_dp) THEN
818 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
819 tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
820 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 120))
821 ELSE
822 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
823 tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
824 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 121))
825 END IF
826 ELSE
827 IF (x2 <= 0.875000000000000000e+00_dp) THEN
828 IF (x1 <= 0.218750000000000000e+00_dp) THEN
829 tg1 = (2*x1 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
830 tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
831 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 122))
832 ELSE
833 tg1 = (2*x1 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
834 tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
835 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 123))
836 END IF
837 ELSE
838 tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
839 tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
840 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 124))
841 END IF
842 END IF
843 END IF
844 END IF
845 ELSE
846 IF (x1 <= 0.375000000000000000e+00_dp) THEN
847 IF (x2 <= 0.750000000000000000e+00_dp) THEN
848 IF (x1 <= 0.312500000000000000e+00_dp) THEN
849 IF (x2 <= 0.625000000000000000e+00_dp) THEN
850 IF (x1 <= 0.281250000000000000e+00_dp) THEN
851 tg1 = (2*x1 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
852 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
853 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 125))
854 ELSE
855 tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
856 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
857 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 126))
858 END IF
859 ELSE
860 IF (x1 <= 0.281250000000000000e+00_dp) THEN
861 tg1 = (2*x1 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
862 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
863 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 127))
864 ELSE
865 tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
866 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
867 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 128))
868 END IF
869 END IF
870 ELSE
871 IF (x2 <= 0.625000000000000000e+00_dp) THEN
872 tg1 = (2*x1 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
873 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
874 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 129))
875 ELSE
876 IF (x1 <= 0.343750000000000000e+00_dp) THEN
877 tg1 = (2*x1 - 0.656250000000000000e+00_dp)*0.320000000000000000e+02_dp
878 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
879 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 130))
880 ELSE
881 tg1 = (2*x1 - 0.718750000000000000e+00_dp)*0.320000000000000000e+02_dp
882 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
883 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 131))
884 END IF
885 END IF
886 END IF
887 ELSE
888 IF (x1 <= 0.312500000000000000e+00_dp) THEN
889 IF (x2 <= 0.875000000000000000e+00_dp) THEN
890 IF (x1 <= 0.281250000000000000e+00_dp) THEN
891 tg1 = (2*x1 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
892 tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
893 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 132))
894 ELSE
895 tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
896 tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
897 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 133))
898 END IF
899 ELSE
900 IF (x1 <= 0.281250000000000000e+00_dp) THEN
901 tg1 = (2*x1 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
902 tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
903 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 134))
904 ELSE
905 tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
906 tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
907 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 135))
908 END IF
909 END IF
910 ELSE
911 IF (x2 <= 0.875000000000000000e+00_dp) THEN
912 IF (x1 <= 0.343750000000000000e+00_dp) THEN
913 tg1 = (2*x1 - 0.656250000000000000e+00_dp)*0.320000000000000000e+02_dp
914 tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
915 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 136))
916 ELSE
917 tg1 = (2*x1 - 0.718750000000000000e+00_dp)*0.320000000000000000e+02_dp
918 tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
919 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 137))
920 END IF
921 ELSE
922 tg1 = (2*x1 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
923 tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
924 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 138))
925 END IF
926 END IF
927 END IF
928 ELSE
929 IF (x2 <= 0.750000000000000000e+00_dp) THEN
930 IF (x1 <= 0.437500000000000000e+00_dp) THEN
931 IF (x2 <= 0.625000000000000000e+00_dp) THEN
932 tg1 = (2*x1 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
933 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
934 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 139))
935 ELSE
936 tg1 = (2*x1 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
937 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
938 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 140))
939 END IF
940 ELSE
941 tg1 = (2*x1 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
942 tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
943 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 141))
944 END IF
945 ELSE
946 IF (x1 <= 0.437500000000000000e+00_dp) THEN
947 IF (x2 <= 0.875000000000000000e+00_dp) THEN
948 tg1 = (2*x1 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
949 tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
950 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 142))
951 ELSE
952 tg1 = (2*x1 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
953 tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
954 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 143))
955 END IF
956 ELSE
957 IF (x2 <= 0.875000000000000000e+00_dp) THEN
958 tg1 = (2*x1 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
959 tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
960 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 144))
961 ELSE
962 tg1 = (2*x1 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
963 tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
964 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 145))
965 END IF
966 END IF
967 END IF
968 END IF
969 END IF
970 END IF
971 ELSE
972 IF (x1 <= 0.750000000000000000e+00_dp) THEN
973 IF (x2 <= 0.500000000000000000e+00_dp) THEN
974 IF (x1 <= 0.625000000000000000e+00_dp) THEN
975 IF (x2 <= 0.250000000000000000e+00_dp) THEN
976 tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
977 tg2 = (2*x2 - 0.250000000000000000e+00_dp)*0.400000000000000000e+01_dp
978 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 146))
979 ELSE
980 tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
981 tg2 = (2*x2 - 0.750000000000000000e+00_dp)*0.400000000000000000e+01_dp
982 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 147))
983 END IF
984 ELSE
985 tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
986 tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
987 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 148))
988 END IF
989 ELSE
990 IF (x1 <= 0.625000000000000000e+00_dp) THEN
991 IF (x2 <= 0.750000000000000000e+00_dp) THEN
992 IF (x1 <= 0.562500000000000000e+00_dp) THEN
993 tg1 = (2*x1 - 0.106250000000000000e+01_dp)*0.160000000000000000e+02_dp
994 tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
995 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 149))
996 ELSE
997 tg1 = (2*x1 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
998 tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
999 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 150))
1000 END IF
1001 ELSE
1002 IF (x1 <= 0.562500000000000000e+00_dp) THEN
1003 tg1 = (2*x1 - 0.106250000000000000e+01_dp)*0.160000000000000000e+02_dp
1004 tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
1005 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 151))
1006 ELSE
1007 tg1 = (2*x1 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
1008 tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
1009 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 152))
1010 END IF
1011 END IF
1012 ELSE
1013 IF (x1 <= 0.687500000000000000e+00_dp) THEN
1014 tg1 = (2*x1 - 0.131250000000000000e+01_dp)*0.160000000000000000e+02_dp
1015 tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1016 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 153))
1017 ELSE
1018 tg1 = (2*x1 - 0.143750000000000000e+01_dp)*0.160000000000000000e+02_dp
1019 tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1020 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 154))
1021 END IF
1022 END IF
1023 END IF
1024 ELSE
1025 tg1 = (2*x1 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
1026 tg2 = (2*x2 - 0.100000000000000000e+01_dp)*0.100000000000000000e+01_dp
1027 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 155))
1028 END IF
1029 END IF
1030 ELSE
1031 IF (t < lower) THEN
1032 use_gamma = .true.
1033 RETURN
1034 END IF
1035 x2 = 11.0_dp/r
1036 x1 = (t - lower)/(upper - lower)
1037 IF (x1 <= 0.500000000000000000e+00_dp) THEN
1038 IF (x1 <= 0.250000000000000000e+00_dp) THEN
1039 IF (x2 <= 0.500000000000000000e+00_dp) THEN
1040 IF (x1 <= 0.125000000000000000e+00_dp) THEN
1041 IF (x2 <= 0.250000000000000000e+00_dp) THEN
1042 tg1 = (2*x1 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
1043 tg2 = (2*x2 - 0.250000000000000000e+00_dp)*0.400000000000000000e+01_dp
1044 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 156))
1045 ELSE
1046 tg1 = (2*x1 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
1047 tg2 = (2*x2 - 0.750000000000000000e+00_dp)*0.400000000000000000e+01_dp
1048 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 157))
1049 END IF
1050 ELSE
1051 tg1 = (2*x1 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
1052 tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
1053 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 158))
1054 END IF
1055 ELSE
1056 IF (x1 <= 0.125000000000000000e+00_dp) THEN
1057 IF (x2 <= 0.750000000000000000e+00_dp) THEN
1058 IF (x2 <= 0.625000000000000000e+00_dp) THEN
1059 tg1 = (2*x1 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
1060 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
1061 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 159))
1062 ELSE
1063 IF (x1 <= 0.625000000000000000e-01_dp) THEN
1064 tg1 = (2*x1 - 0.625000000000000000e-01_dp)*0.160000000000000000e+02_dp
1065 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
1066 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 160))
1067 ELSE
1068 tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
1069 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
1070 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 161))
1071 END IF
1072 END IF
1073 ELSE
1074 IF (x1 <= 0.625000000000000000e-01_dp) THEN
1075 IF (x2 <= 0.875000000000000000e+00_dp) THEN
1076 IF (x1 <= 0.312500000000000000e-01_dp) THEN
1077 IF (x2 <= 0.812500000000000000e+00_dp) THEN
1078 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
1079 tg2 = (2*x2 - 0.156250000000000000e+01_dp)*0.160000000000000000e+02_dp
1080 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 162))
1081 ELSE
1082 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
1083 tg2 = (2*x2 - 0.168750000000000000e+01_dp)*0.160000000000000000e+02_dp
1084 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 163))
1085 END IF
1086 ELSE
1087 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
1088 tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
1089 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 164))
1090 END IF
1091 ELSE
1092 IF (x1 <= 0.312500000000000000e-01_dp) THEN
1093 IF (x2 <= 0.937500000000000000e+00_dp) THEN
1094 IF (x1 <= 0.156250000000000000e-01_dp) THEN
1095 IF (x2 <= 0.906250000000000000e+00_dp) THEN
1096 tg1 = (2*x1 - 0.156250000000000000e-01_dp)*0.640000000000000000e+02_dp
1097 tg2 = (2*x2 - 0.178125000000000000e+01_dp)*0.320000000000000000e+02_dp
1098 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 165))
1099 ELSE
1100 tg1 = (2*x1 - 0.156250000000000000e-01_dp)*0.640000000000000000e+02_dp
1101 tg2 = (2*x2 - 0.184375000000000000e+01_dp)*0.320000000000000000e+02_dp
1102 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 166))
1103 END IF
1104 ELSE
1105 tg1 = (2*x1 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
1106 tg2 = (2*x2 - 0.181250000000000000e+01_dp)*0.160000000000000000e+02_dp
1107 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 167))
1108 END IF
1109 ELSE
1110 IF (x1 <= 0.156250000000000000e-01_dp) THEN
1111 IF (x2 <= 0.968750000000000000e+00_dp) THEN
1112 IF (x1 <= 0.781250000000000000e-02_dp) THEN
1113 IF (x2 <= 0.953125000000000000e+00_dp) THEN
1114 tg1 = (2*x1 - 0.781250000000000000e-02_dp)*0.128000000000000000e+03_dp
1115 tg2 = (2*x2 - 0.189062500000000000e+01_dp)*0.640000000000000000e+02_dp
1116 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 168))
1117 ELSE
1118 tg1 = (2*x1 - 0.781250000000000000e-02_dp)*0.128000000000000000e+03_dp
1119 tg2 = (2*x2 - 0.192187500000000000e+01_dp)*0.640000000000000000e+02_dp
1120 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 169))
1121 END IF
1122 ELSE
1123 tg1 = (2*x1 - 0.234375000000000000e-01_dp)*0.128000000000000000e+03_dp
1124 tg2 = (2*x2 - 0.190625000000000000e+01_dp)*0.320000000000000000e+02_dp
1125 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 170))
1126 END IF
1127 ELSE
1128 IF (x1 <= 0.781250000000000000e-02_dp) THEN
1129 IF (x2 <= 0.984375000000000000e+00_dp) THEN
1130 tg1 = (2*x1 - 0.781250000000000000e-02_dp)*0.128000000000000000e+03_dp
1131 tg2 = (2*x2 - 0.195312500000000000e+01_dp)*0.640000000000000000e+02_dp
1132 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 171))
1133 ELSE
1134 tg1 = (2*x1 - 0.781250000000000000e-02_dp)*0.128000000000000000e+03_dp
1135 tg2 = (2*x2 - 0.198437500000000000e+01_dp)*0.640000000000000000e+02_dp
1136 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 172))
1137 END IF
1138 ELSE
1139 IF (x2 <= 0.984375000000000000e+00_dp) THEN
1140 tg1 = (2*x1 - 0.234375000000000000e-01_dp)*0.128000000000000000e+03_dp
1141 tg2 = (2*x2 - 0.195312500000000000e+01_dp)*0.640000000000000000e+02_dp
1142 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 173))
1143 ELSE
1144 tg1 = (2*x1 - 0.234375000000000000e-01_dp)*0.128000000000000000e+03_dp
1145 tg2 = (2*x2 - 0.198437500000000000e+01_dp)*0.640000000000000000e+02_dp
1146 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 174))
1147 END IF
1148 END IF
1149 END IF
1150 ELSE
1151 IF (x2 <= 0.968750000000000000e+00_dp) THEN
1152 tg1 = (2*x1 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
1153 tg2 = (2*x2 - 0.190625000000000000e+01_dp)*0.320000000000000000e+02_dp
1154 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 175))
1155 ELSE
1156 IF (x1 <= 0.234375000000000000e-01_dp) THEN
1157 tg1 = (2*x1 - 0.390625000000000000e-01_dp)*0.128000000000000000e+03_dp
1158 tg2 = (2*x2 - 0.196875000000000000e+01_dp)*0.320000000000000000e+02_dp
1159 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 176))
1160 ELSE
1161 tg1 = (2*x1 - 0.546875000000000000e-01_dp)*0.128000000000000000e+03_dp
1162 tg2 = (2*x2 - 0.196875000000000000e+01_dp)*0.320000000000000000e+02_dp
1163 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 177))
1164 END IF
1165 END IF
1166 END IF
1167 END IF
1168 ELSE
1169 IF (x2 <= 0.937500000000000000e+00_dp) THEN
1170 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
1171 tg2 = (2*x2 - 0.181250000000000000e+01_dp)*0.160000000000000000e+02_dp
1172 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 178))
1173 ELSE
1174 IF (x1 <= 0.468750000000000000e-01_dp) THEN
1175 IF (x2 <= 0.968750000000000000e+00_dp) THEN
1176 tg1 = (2*x1 - 0.781250000000000000e-01_dp)*0.640000000000000000e+02_dp
1177 tg2 = (2*x2 - 0.190625000000000000e+01_dp)*0.320000000000000000e+02_dp
1178 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 179))
1179 ELSE
1180 tg1 = (2*x1 - 0.781250000000000000e-01_dp)*0.640000000000000000e+02_dp
1181 tg2 = (2*x2 - 0.196875000000000000e+01_dp)*0.320000000000000000e+02_dp
1182 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 180))
1183 END IF
1184 ELSE
1185 tg1 = (2*x1 - 0.109375000000000000e+00_dp)*0.640000000000000000e+02_dp
1186 tg2 = (2*x2 - 0.193750000000000000e+01_dp)*0.160000000000000000e+02_dp
1187 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 181))
1188 END IF
1189 END IF
1190 END IF
1191 END IF
1192 ELSE
1193 IF (x2 <= 0.875000000000000000e+00_dp) THEN
1194 tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
1195 tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
1196 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 182))
1197 ELSE
1198 IF (x1 <= 0.937500000000000000e-01_dp) THEN
1199 IF (x2 <= 0.937500000000000000e+00_dp) THEN
1200 tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
1201 tg2 = (2*x2 - 0.181250000000000000e+01_dp)*0.160000000000000000e+02_dp
1202 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 183))
1203 ELSE
1204 tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
1205 tg2 = (2*x2 - 0.193750000000000000e+01_dp)*0.160000000000000000e+02_dp
1206 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 184))
1207 END IF
1208 ELSE
1209 tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
1210 tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
1211 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 185))
1212 END IF
1213 END IF
1214 END IF
1215 END IF
1216 ELSE
1217 IF (x2 <= 0.750000000000000000e+00_dp) THEN
1218 tg1 = (2*x1 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
1219 tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
1220 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 186))
1221 ELSE
1222 IF (x1 <= 0.187500000000000000e+00_dp) THEN
1223 IF (x2 <= 0.875000000000000000e+00_dp) THEN
1224 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
1225 tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
1226 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 187))
1227 ELSE
1228 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
1229 tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
1230 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 188))
1231 END IF
1232 ELSE
1233 tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
1234 tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
1235 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 189))
1236 END IF
1237 END IF
1238 END IF
1239 END IF
1240 ELSE
1241 IF (x1 <= 0.375000000000000000e+00_dp) THEN
1242 IF (x1 <= 0.312500000000000000e+00_dp) THEN
1243 IF (x1 <= 0.281250000000000000e+00_dp) THEN
1244 tg1 = (2*x1 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
1245 tg2 = (2*x2 - 0.100000000000000000e+01_dp)*0.100000000000000000e+01_dp
1246 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 190))
1247 ELSE
1248 tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
1249 tg2 = (2*x2 - 0.100000000000000000e+01_dp)*0.100000000000000000e+01_dp
1250 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 191))
1251 END IF
1252 ELSE
1253 IF (x2 <= 0.500000000000000000e+00_dp) THEN
1254 tg1 = (2*x1 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
1255 tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
1256 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 192))
1257 ELSE
1258 tg1 = (2*x1 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
1259 tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1260 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 193))
1261 END IF
1262 END IF
1263 ELSE
1264 IF (x1 <= 0.437500000000000000e+00_dp) THEN
1265 IF (x2 <= 0.500000000000000000e+00_dp) THEN
1266 tg1 = (2*x1 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
1267 tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
1268 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 194))
1269 ELSE
1270 IF (x1 <= 0.406250000000000000e+00_dp) THEN
1271 tg1 = (2*x1 - 0.781250000000000000e+00_dp)*0.320000000000000000e+02_dp
1272 tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1273 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 195))
1274 ELSE
1275 tg1 = (2*x1 - 0.843750000000000000e+00_dp)*0.320000000000000000e+02_dp
1276 tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1277 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 196))
1278 END IF
1279 END IF
1280 ELSE
1281 IF (x2 <= 0.500000000000000000e+00_dp) THEN
1282 tg1 = (2*x1 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
1283 tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
1284 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 197))
1285 ELSE
1286 tg1 = (2*x1 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
1287 tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1288 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 198))
1289 END IF
1290 END IF
1291 END IF
1292 END IF
1293 ELSE
1294 IF (x1 <= 0.750000000000000000e+00_dp) THEN
1295 IF (x1 <= 0.625000000000000000e+00_dp) THEN
1296 IF (x2 <= 0.500000000000000000e+00_dp) THEN
1297 IF (x1 <= 0.562500000000000000e+00_dp) THEN
1298 tg1 = (2*x1 - 0.106250000000000000e+01_dp)*0.160000000000000000e+02_dp
1299 tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
1300 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 199))
1301 ELSE
1302 tg1 = (2*x1 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
1303 tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
1304 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 200))
1305 END IF
1306 ELSE
1307 IF (x1 <= 0.562500000000000000e+00_dp) THEN
1308 tg1 = (2*x1 - 0.106250000000000000e+01_dp)*0.160000000000000000e+02_dp
1309 tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1310 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 201))
1311 ELSE
1312 tg1 = (2*x1 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
1313 tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1314 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 202))
1315 END IF
1316 END IF
1317 ELSE
1318 IF (x2 <= 0.500000000000000000e+00_dp) THEN
1319 IF (x1 <= 0.687500000000000000e+00_dp) THEN
1320 tg1 = (2*x1 - 0.131250000000000000e+01_dp)*0.160000000000000000e+02_dp
1321 tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
1322 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 203))
1323 ELSE
1324 tg1 = (2*x1 - 0.143750000000000000e+01_dp)*0.160000000000000000e+02_dp
1325 tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
1326 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 204))
1327 END IF
1328 ELSE
1329 tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
1330 tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1331 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 205))
1332 END IF
1333 END IF
1334 ELSE
1335 IF (x1 <= 0.875000000000000000e+00_dp) THEN
1336 tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
1337 tg2 = (2*x2 - 0.100000000000000000e+01_dp)*0.100000000000000000e+01_dp
1338 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 206))
1339 ELSE
1340 tg1 = (2*x1 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
1341 tg2 = (2*x2 - 0.100000000000000000e+01_dp)*0.100000000000000000e+01_dp
1342 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 207))
1343 END IF
1344 END IF
1345 END IF
1346 END IF
1347 END SUBROUTINE t_c_g0_n
1348
1349! **************************************************************************************************
1350!> \brief ...
1351!> \param Nder the number of derivatives that will actually be used
1352!> \param iunit contains the data file to initialize the table
1353!> \param mepos ...
1354!> \param group ...
1355! **************************************************************************************************
1356 SUBROUTINE init(Nder, iunit, mepos, group)
1357 INTEGER, INTENT(IN) :: nder, iunit, mepos
1358
1359 CLASS(mp_comm_type), INTENT(IN) :: group
1360
1361 INTEGER :: i
1362 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: chunk
1363
1364 patches = 207
1365 IF (nder > nderiv_max) cpabort("T_C_G0 init failed")
1366 nderiv_init = nder
1367 IF (ALLOCATED(c0)) DEALLOCATE (c0)
1368 ! round up to a multiple of 32 to give some generous alignment for each C0
1369 ALLOCATE (c0(32*((31 + (nder + 1)*(degree + 1)*(degree + 2)/2)/32), patches))
1370 ! init mpi'ed buffers to silence warnings under valgrind
1371 c0 = 1.0e99_dp
1372 IF (mepos == 0) THEN
1373 ALLOCATE (chunk((nderiv_max + 1)*(degree + 1)*(degree + 2)/2))
1374 DO i = 1, patches
1375 READ (iunit, *) chunk
1376 c0(1:(nder + 1)*(degree + 1)*(degree + 2)/2, i) = chunk(1:(nder + 1)*(degree + 1)*(degree + 2)/2)
1377 END DO
1378 DEALLOCATE (chunk)
1379 END IF
1380 CALL group%bcast(c0, 0)
1381
1382 END SUBROUTINE init
1383
1384! **************************************************************************************************
1385!> \brief ...
1386! **************************************************************************************************
1387 SUBROUTINE free_c0()
1388 IF (ALLOCATED(c0)) DEALLOCATE (c0)
1389 nderiv_init = -1
1390 END SUBROUTINE free_c0
1391
1392! **************************************************************************************************
1393!> \brief ...
1394!> \param RES ...
1395!> \param NDERIV ...
1396!> \param TG1 ...
1397!> \param TG2 ...
1398!> \param C0 ...
1399! **************************************************************************************************
1400 SUBROUTINE pd2val(RES, NDERIV, TG1, TG2, C0)
1401 REAL(kind=dp), INTENT(OUT) :: res(*)
1402 INTEGER, INTENT(IN) :: nderiv
1403 REAL(kind=dp), INTENT(IN) :: tg1, tg2, c0(105, *)
1404
1405 REAL(kind=dp), PARAMETER :: sqrt2 = 1.4142135623730950488016887242096980785696718753_dp
1406
1407 INTEGER :: k
1408 REAL(kind=dp) :: t1(0:13), t2(0:13)
1409
1410 t1(0) = 1.0_dp
1411 t2(0) = 1.0_dp
1412 t1(1) = sqrt2*tg1
1413 t2(1) = sqrt2*tg2
1414 t1(2) = 2*tg1*t1(1) - sqrt2
1415 t2(2) = 2*tg2*t2(1) - sqrt2
1416 t1(3) = 2*tg1*t1(2) - t1(1)
1417 t2(3) = 2*tg2*t2(2) - t2(1)
1418 t1(4) = 2*tg1*t1(3) - t1(2)
1419 t2(4) = 2*tg2*t2(3) - t2(2)
1420 t1(5) = 2*tg1*t1(4) - t1(3)
1421 t2(5) = 2*tg2*t2(4) - t2(3)
1422 t1(6) = 2*tg1*t1(5) - t1(4)
1423 t2(6) = 2*tg2*t2(5) - t2(4)
1424 t1(7) = 2*tg1*t1(6) - t1(5)
1425 t2(7) = 2*tg2*t2(6) - t2(5)
1426 t1(8) = 2*tg1*t1(7) - t1(6)
1427 t2(8) = 2*tg2*t2(7) - t2(6)
1428 t1(9) = 2*tg1*t1(8) - t1(7)
1429 t2(9) = 2*tg2*t2(8) - t2(7)
1430 t1(10) = 2*tg1*t1(9) - t1(8)
1431 t2(10) = 2*tg2*t2(9) - t2(8)
1432 t1(11) = 2*tg1*t1(10) - t1(9)
1433 t2(11) = 2*tg2*t2(10) - t2(9)
1434 t1(12) = 2*tg1*t1(11) - t1(10)
1435 t2(12) = 2*tg2*t2(11) - t2(10)
1436 t1(13) = 2*tg1*t1(12) - t1(11)
1437 t2(13) = 2*tg2*t2(12) - t2(11)
1438 DO k = 1, nderiv + 1
1439 res(k) = 0.0_dp
1440 res(k) = res(k) + dot_product(t1(0:13), c0(1:14, k))*t2(0)
1441 res(k) = res(k) + dot_product(t1(0:12), c0(15:27, k))*t2(1)
1442 res(k) = res(k) + dot_product(t1(0:11), c0(28:39, k))*t2(2)
1443 res(k) = res(k) + dot_product(t1(0:10), c0(40:50, k))*t2(3)
1444 res(k) = res(k) + dot_product(t1(0:9), c0(51:60, k))*t2(4)
1445 res(k) = res(k) + dot_product(t1(0:8), c0(61:69, k))*t2(5)
1446 res(k) = res(k) + dot_product(t1(0:7), c0(70:77, k))*t2(6)
1447 res(k) = res(k) + dot_product(t1(0:6), c0(78:84, k))*t2(7)
1448 res(k) = res(k) + dot_product(t1(0:5), c0(85:90, k))*t2(8)
1449 res(k) = res(k) + dot_product(t1(0:4), c0(91:95, k))*t2(9)
1450 res(k) = res(k) + dot_product(t1(0:3), c0(96:99, k))*t2(10)
1451 res(k) = res(k) + dot_product(t1(0:2), c0(100:102, k))*t2(11)
1452 res(k) = res(k) + dot_product(t1(0:1), c0(103:104, k))*t2(12)
1453 res(k) = res(k) + dot_product(t1(0:0), c0(105:105, k))*t2(13)
1454 END DO
1455 END SUBROUTINE pd2val
1456
1457! **************************************************************************************************
1458!> \brief Returns the value of nderiv_init so that one can check if opening the potential file is
1459!> worhtwhile
1460!> \return ...
1461!> \author A. Bussy, 05.2019
1462! **************************************************************************************************
1463 FUNCTION get_lmax_init() RESULT(lmax_init)
1464
1465 INTEGER :: lmax_init
1466
1467 lmax_init = nderiv_init
1468
1469 END FUNCTION get_lmax_init
1470
1471END MODULE t_c_g0
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
t_c_g0
This module computes the basic integrals for the truncated coulomb operator.
Definition t_c_g0.F:57
t_c_g0::t_c_g0_n
subroutine, public t_c_g0_n(res, use_gamma, r, t, nderiv)
...
Definition t_c_g0.F:84
t_c_g0::init
subroutine, public init(nder, iunit, mepos, group)
...
Definition t_c_g0.F:1357
t_c_g0::get_lmax_init
integer function, public get_lmax_init()
Returns the value of nderiv_init so that one can check if opening the potential file is worhtwhile.
Definition t_c_g0.F:1464
t_c_g0::free_c0
subroutine, public free_c0()
...
Definition t_c_g0.F:1388
message_passing::mp_comm_type
Definition message_passing.F:189