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