Line data Source code
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 : !> \brief Calculation of three-center overlap integrals over Cartesian
10 : !> Gaussian-type functions for the second term V(ppl) of the local
11 : !> part of the Goedecker pseudopotential (GTH):
12 : !>
13 : !> <a|V(local)|b> = <a|V(erf) + V(ppl)|b>
14 : !> = <a|V(erf)|b> + <a|V(ppl)|b>
15 : !> = <a|-Z(eff)*erf(SQRT(2)*alpha*r)/r +
16 : !> (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
17 : !> C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
18 : !> \par Literature
19 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
20 : !> S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B 54, 1703 (1996)
21 : !> C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B 58, 3641 (1998)
22 : !> \par History
23 : !> - Derivatives added (17.05.2002,MK)
24 : !> - Complete refactoring (05.2011,jhu)
25 : !> \author Matthias Krack (04.10.2000)
26 : ! **************************************************************************************************
27 : MODULE ai_overlap_ppl
28 : USE ai_oneelectron, ONLY: os_2center,&
29 : os_3center
30 : USE gamma, ONLY: fgamma => fgamma_0
31 : USE gfun, ONLY: gfun_values
32 : USE kinds, ONLY: dp
33 : USE mathconstants, ONLY: pi
34 : USE mathlib, ONLY: binomial
35 : #include "../base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 :
39 : PRIVATE
40 :
41 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap_ppl'
42 :
43 : ! *** Public subroutines ***
44 :
45 : PUBLIC :: ecploc_integral, ppl_integral, ppl_integral_ri
46 :
47 : CONTAINS
48 :
49 : ! **************************************************************************************************
50 : !> \brief Calculation of three-center overlap integrals <a|c|b> over
51 : !> Cartesian Gaussian functions for the local part of the Goedecker
52 : !> pseudopotential (GTH). c is a primitive Gaussian-type function
53 : !> with a set of even angular momentum indices.
54 : !>
55 : !> <a|V(ppl)|b> = <a| (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
56 : !> C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
57 : !> zetc = alpha**2/2
58 : !>
59 : !> \param la_max_set ...
60 : !> \param la_min_set ...
61 : !> \param npgfa ...
62 : !> \param rpgfa ...
63 : !> \param zeta ...
64 : !> \param lb_max_set ...
65 : !> \param lb_min_set ...
66 : !> \param npgfb ...
67 : !> \param rpgfb ...
68 : !> \param zetb ...
69 : !> \param nexp_ppl ...
70 : !> \param alpha_ppl ...
71 : !> \param nct_ppl ...
72 : !> \param cexp_ppl ...
73 : !> \param rpgfc ...
74 : !> \param rab ...
75 : !> \param dab ...
76 : !> \param rac ...
77 : !> \param dac ...
78 : !> \param rbc ...
79 : !> \param dbc ...
80 : !> \param vab ...
81 : !> \param s ...
82 : !> \param pab ...
83 : !> \param force_a ...
84 : !> \param force_b ...
85 : !> \param fs ...
86 : !> \param hab2 The derivative of the ppl integrals according to the weighting factors deltaR
87 : !> \param hab2_work ...
88 : !> \param deltaR Weighting factors for the derivatives wrt. nuclear positions
89 : !> \param iatom ...
90 : !> \param jatom ...
91 : !> \param katom ...
92 : !> \date May 2011
93 : !> \author Juerg Hutter
94 : !> \version 1.0
95 : !> \note Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
96 : ! **************************************************************************************************
97 51005201 : SUBROUTINE ppl_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
98 51005201 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
99 102010402 : rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
100 102010402 : hab2, hab2_work, deltaR, iatom, jatom, katom)
101 : INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
102 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
103 : INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
104 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
105 : INTEGER, INTENT(IN) :: nexp_ppl
106 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: alpha_ppl
107 : INTEGER, DIMENSION(:), INTENT(IN) :: nct_ppl
108 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: cexp_ppl
109 : REAL(KIND=dp), INTENT(IN) :: rpgfc
110 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
111 : REAL(KIND=dp), INTENT(IN) :: dab
112 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
113 : REAL(KIND=dp), INTENT(IN) :: dac
114 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rbc
115 : REAL(KIND=dp), INTENT(IN) :: dbc
116 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
117 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: s
118 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
119 : OPTIONAL :: pab
120 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a, force_b
121 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
122 : OPTIONAL :: fs, hab2, hab2_work
123 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
124 : OPTIONAL :: deltaR
125 : INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom, katom
126 :
127 : INTEGER :: iexp, ij, ipgf, jpgf, mmax, nexp
128 : REAL(KIND=dp) :: rho, sab, t, zetc
129 51005201 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: auxint
130 : REAL(KIND=dp), DIMENSION(3) :: pci
131 :
132 51005201 : IF (PRESENT(pab)) THEN
133 14687953 : CPASSERT(PRESENT(force_a))
134 14687953 : CPASSERT(PRESENT(force_b))
135 14687953 : CPASSERT(PRESENT(fs))
136 14687953 : mmax = la_max_set + lb_max_set + 2
137 14687953 : force_a(:) = 0.0_dp
138 14687953 : force_b(:) = 0.0_dp
139 36317248 : ELSE IF (PRESENT(hab2)) THEN
140 870 : mmax = la_max_set + lb_max_set + 2
141 : ELSE
142 36316378 : mmax = la_max_set + lb_max_set
143 : END IF
144 :
145 204020804 : ALLOCATE (auxint(0:mmax, npgfa*npgfb))
146 51005201 : auxint = 0._dp
147 :
148 : ! *** Calculate auxiliary integrals ***
149 :
150 236259897 : DO ipgf = 1, npgfa
151 : ! *** Screening ***
152 185254696 : IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
153 365949469 : DO jpgf = 1, npgfb
154 : ! *** Screening ***
155 248573414 : IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
156 : (rpgfa(ipgf) + rpgfb(jpgf) < dab)) CYCLE
157 87811445 : ij = (ipgf - 1)*npgfb + jpgf
158 87811445 : rho = zeta(ipgf) + zetb(jpgf)
159 351245780 : pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
160 87811445 : sab = EXP(-(zeta(ipgf)*zetb(jpgf)/rho*dab*dab))
161 351245780 : t = rho*SUM(pci(:)*pci(:))
162 :
163 175890328 : DO iexp = 1, nexp_ppl
164 88078883 : nexp = nct_ppl(iexp)
165 88078883 : zetc = alpha_ppl(iexp)
166 175890328 : CALL ppl_aux(auxint(0:mmax, ij), mmax, t, rho, nexp, cexp_ppl(:, iexp), zetc)
167 : END DO
168 :
169 596586755 : auxint(0:mmax, ij) = sab*auxint(0:mmax, ij)
170 :
171 : END DO
172 : END DO
173 :
174 : CALL os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
175 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
176 : rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
177 : vab2=hab2, vab2_work=hab2_work, &
178 276651820 : deltaR=deltaR, iatom=iatom, jatom=jatom, katom=katom)
179 :
180 51005201 : DEALLOCATE (auxint)
181 :
182 51005201 : END SUBROUTINE ppl_integral
183 :
184 : ! **************************************************************************************************
185 : !> \brief Calculation of three-center potential integrals <a|V(r)|b> over
186 : !> Cartesian Gaussian functions for the local part of ECP
187 : !> pseudopotential. Multiple terms C1-4 are possible.
188 : !>
189 : !> <a|V(ecploc)|b> = <a| C1/r*exp(-a1*r**2) + C2*exp(-a2*r**2) + C3*r*exp(-a3*r**2) +
190 : !> C4*r**2*exp(-a4*r**2)|b>
191 : !>
192 : !> \param la_max_set ...
193 : !> \param la_min_set ...
194 : !> \param npgfa ...
195 : !> \param rpgfa ...
196 : !> \param zeta ...
197 : !> \param lb_max_set ...
198 : !> \param lb_min_set ...
199 : !> \param npgfb ...
200 : !> \param rpgfb ...
201 : !> \param zetb ...
202 : !> \param nexp_ppl ...
203 : !> \param alpha_ppl ...
204 : !> \param nct_ppl ...
205 : !> \param cexp_ppl ...
206 : !> \param rpgfc ...
207 : !> \param rab ...
208 : !> \param dab ...
209 : !> \param rac ...
210 : !> \param dac ...
211 : !> \param rbc ...
212 : !> \param dbc ...
213 : !> \param vab ...
214 : !> \param s ...
215 : !> \param pab ...
216 : !> \param force_a ...
217 : !> \param force_b ...
218 : !> \param fs ...
219 : !> \param hab2 The derivative of the ppl integrals according to the weighting factors deltaR
220 : !> \param hab2_work ...
221 : !> \param deltaR Weighting factors for the derivatives wrt. nuclear positions
222 : !> \param iatom ...
223 : !> \param jatom ...
224 : !> \param katom ...
225 : !> \date 2025
226 : !> \author Juerg Hutter
227 : !> \version 1.0
228 : ! **************************************************************************************************
229 37648 : SUBROUTINE ecploc_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
230 37648 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
231 37648 : nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
232 75296 : rab, dab, rac, dac, rbc, dbc, vab, s, pab, &
233 75296 : force_a, force_b, fs, hab2, hab2_work, &
234 37648 : deltaR, iatom, jatom, katom)
235 : INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
236 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
237 : INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
238 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
239 : INTEGER, INTENT(IN) :: nexp_ppl
240 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: alpha_ppl
241 : INTEGER, DIMENSION(:), INTENT(IN) :: nct_ppl
242 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: cexp_ppl
243 : REAL(KIND=dp), INTENT(IN) :: rpgfc
244 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
245 : REAL(KIND=dp), INTENT(IN) :: dab
246 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
247 : REAL(KIND=dp), INTENT(IN) :: dac
248 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rbc
249 : REAL(KIND=dp), INTENT(IN) :: dbc
250 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
251 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: s
252 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
253 : OPTIONAL :: pab
254 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a, force_b
255 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
256 : OPTIONAL :: fs, hab2, hab2_work
257 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
258 : OPTIONAL :: deltaR
259 : INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom, katom
260 :
261 : INTEGER :: iexp, ij, ipgf, jpgf, mmax, nexp
262 : REAL(KIND=dp) :: rho, sab, t, zetc
263 37648 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: auxint
264 : REAL(KIND=dp), DIMENSION(3) :: pci
265 :
266 37648 : IF (PRESENT(pab)) THEN
267 12394 : CPASSERT(PRESENT(force_a))
268 12394 : CPASSERT(PRESENT(force_b))
269 12394 : CPASSERT(PRESENT(fs))
270 12394 : mmax = la_max_set + lb_max_set + 2
271 12394 : force_a(:) = 0.0_dp
272 12394 : force_b(:) = 0.0_dp
273 25254 : ELSE IF (PRESENT(hab2)) THEN
274 0 : mmax = la_max_set + lb_max_set + 2
275 : ELSE
276 25254 : mmax = la_max_set + lb_max_set
277 : END IF
278 :
279 150592 : ALLOCATE (auxint(0:mmax, npgfa*npgfb))
280 37648 : auxint = 0._dp
281 :
282 : ! *** Calculate auxiliary integrals ***
283 :
284 141595 : DO ipgf = 1, npgfa
285 : ! *** Screening ***
286 103947 : IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
287 244770 : DO jpgf = 1, npgfb
288 : ! *** Screening ***
289 153392 : IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
290 : (rpgfa(ipgf) + rpgfb(jpgf) < dab)) CYCLE
291 67479 : ij = (ipgf - 1)*npgfb + jpgf
292 67479 : rho = zeta(ipgf) + zetb(jpgf)
293 269916 : pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
294 67479 : sab = EXP(-(zeta(ipgf)*zetb(jpgf)/rho*dab*dab))
295 269916 : t = rho*SUM(pci(:)*pci(:))
296 :
297 291853 : DO iexp = 1, nexp_ppl
298 224374 : nexp = nct_ppl(iexp)
299 224374 : zetc = alpha_ppl(iexp)
300 291853 : CALL ecploc_aux(auxint(0:mmax, ij), mmax, t, rho, nexp, cexp_ppl(1, iexp), zetc)
301 : END DO
302 :
303 380480 : auxint(0:mmax, ij) = sab*auxint(0:mmax, ij)
304 :
305 : END DO
306 : END DO
307 :
308 : CALL os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
309 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
310 : rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
311 : vab2=hab2, vab2_work=hab2_work, &
312 201100 : deltaR=deltaR, iatom=iatom, jatom=jatom, katom=katom)
313 :
314 37648 : DEALLOCATE (auxint)
315 :
316 37648 : END SUBROUTINE ecploc_integral
317 : ! **************************************************************************************************
318 : !> \brief Calculation of two-center overlap integrals <a|c> over
319 : !> Cartesian Gaussian functions for the local part of the Goedecker
320 : !> pseudopotential (GTH). c is a primitive Gaussian-type function
321 : !> with a set of even angular momentum indices.
322 : !>
323 : !> <a|V(ppl)|b> = <a| (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
324 : !> C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
325 : !> zetc = alpha**2/2
326 : !>
327 : !> \param la_max_set ...
328 : !> \param la_min_set ...
329 : !> \param npgfa ...
330 : !> \param rpgfa ...
331 : !> \param zeta ...
332 : !> \param nexp_ppl ...
333 : !> \param alpha_ppl ...
334 : !> \param nct_ppl ...
335 : !> \param cexp_ppl ...
336 : !> \param rpgfc ...
337 : !> \param rac ...
338 : !> \param dac ...
339 : !> \param va ...
340 : !> \param dva ...
341 : !> \date December 2017
342 : !> \author Juerg Hutter
343 : !> \version 1.0
344 : ! **************************************************************************************************
345 328 : SUBROUTINE ppl_integral_ri(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
346 328 : nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
347 656 : rac, dac, va, dva)
348 : INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
349 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
350 : INTEGER, INTENT(IN) :: nexp_ppl
351 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: alpha_ppl
352 : INTEGER, DIMENSION(:), INTENT(IN) :: nct_ppl
353 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: cexp_ppl
354 : REAL(KIND=dp), INTENT(IN) :: rpgfc
355 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
356 : REAL(KIND=dp), INTENT(IN) :: dac
357 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: va
358 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
359 : OPTIONAL :: dva
360 :
361 : INTEGER :: iexp, ipgf, mmax, nexp
362 : REAL(KIND=dp) :: rho, t, zetc
363 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: auxint
364 :
365 328 : IF (PRESENT(dva)) THEN
366 164 : mmax = la_max_set + 1
367 : ELSE
368 164 : mmax = la_max_set
369 : END IF
370 :
371 1312 : ALLOCATE (auxint(0:mmax, npgfa))
372 328 : auxint = 0._dp
373 :
374 : ! *** Calculate auxiliary integrals ***
375 656 : DO ipgf = 1, npgfa
376 328 : IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
377 328 : rho = zeta(ipgf)
378 328 : t = rho*dac*dac
379 :
380 984 : DO iexp = 1, nexp_ppl
381 328 : nexp = nct_ppl(iexp)
382 328 : zetc = alpha_ppl(iexp)
383 656 : CALL ppl_aux(auxint(0:mmax, ipgf), mmax, t, rho, nexp, cexp_ppl(:, iexp), zetc)
384 : END DO
385 :
386 : END DO
387 :
388 328 : IF (PRESENT(dva)) THEN
389 : CALL os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
390 164 : auxint, rpgfc, rac, dac, va, dva)
391 : ELSE
392 : CALL os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
393 164 : auxint, rpgfc, rac, dac, va)
394 : END IF
395 :
396 328 : DEALLOCATE (auxint)
397 :
398 328 : END SUBROUTINE ppl_integral_ri
399 :
400 : ! **************************************************************************************************
401 : !> \brief ...
402 : !> \param auxint ...
403 : !> \param mmax ...
404 : !> \param t ...
405 : !> \param rho ...
406 : !> \param nexp_ppl ...
407 : !> \param cexp_ppl ...
408 : !> \param zetc ...
409 : ! **************************************************************************************************
410 88079211 : SUBROUTINE ppl_aux(auxint, mmax, t, rho, nexp_ppl, cexp_ppl, zetc)
411 : INTEGER, INTENT(IN) :: mmax
412 : REAL(KIND=dp), DIMENSION(0:mmax) :: auxint
413 : REAL(KIND=dp), INTENT(IN) :: t, rho
414 : INTEGER, INTENT(IN) :: nexp_ppl
415 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cexp_ppl
416 : REAL(KIND=dp), INTENT(IN) :: zetc
417 :
418 : INTEGER :: i, j, ke, kp, pmax
419 : REAL(KIND=dp) :: a2, a3, a4, cc, f, q, q2, q4, q6, rho2, &
420 : rho3, t2, t3
421 : REAL(KIND=dp), DIMENSION(0:6) :: polder
422 176158422 : REAL(KIND=dp), DIMENSION(0:mmax) :: expder
423 :
424 88079211 : CPASSERT(nexp_ppl > 0)
425 88079211 : q = rho + zetc
426 88079211 : polder = 0._dp
427 88079211 : pmax = 0
428 88079211 : IF (nexp_ppl > 0) THEN
429 88079211 : polder(0) = polder(0) + cexp_ppl(1)
430 88079211 : pmax = 0
431 : END IF
432 88079211 : IF (nexp_ppl > 1) THEN
433 73550033 : q2 = q*q
434 73550033 : a2 = 0.5_dp/q2*cexp_ppl(2)
435 73550033 : polder(0) = polder(0) + a2*(2._dp*rho*t + 3._dp*q)
436 73550033 : polder(1) = polder(1) - a2*2._dp*rho
437 73550033 : pmax = 1
438 : END IF
439 88079211 : IF (nexp_ppl > 2) THEN
440 1007012 : q4 = q2*q2
441 1007012 : rho2 = rho*rho
442 1007012 : t2 = t*t
443 1007012 : a3 = 0.25_dp/q4*cexp_ppl(3)
444 1007012 : polder(0) = polder(0) + a3*(4._dp*rho2*t2 + 20._dp*rho*t*q + 15._dp*q2)
445 1007012 : polder(1) = polder(1) - a3*(8._dp*rho2*t + 20._dp*rho*q)
446 1007012 : polder(2) = polder(2) + a3*8._dp*rho2
447 1007012 : pmax = 2
448 : END IF
449 88079211 : IF (nexp_ppl > 3) THEN
450 1007012 : q6 = q4*q2
451 1007012 : rho3 = rho2*rho
452 1007012 : t3 = t2*t
453 1007012 : a4 = 0.125_dp/q6*cexp_ppl(4)
454 1007012 : polder(0) = polder(0) + a4*(8._dp*rho3*t3 + 84._dp*rho2*t2*q + 210._dp*rho*t*q2 + 105._dp*q*q2)
455 1007012 : polder(1) = polder(1) - a4*(24._dp*rho3*t2 + 168._dp*rho2*t*q + 210._dp*rho*q2)
456 1007012 : polder(2) = polder(2) + a4*(48._dp*rho3*t + 168._dp*rho2*q)
457 1007012 : polder(3) = polder(3) - a4*48_dp*rho3
458 1007012 : pmax = 3
459 : END IF
460 88079211 : IF (nexp_ppl > 4) THEN
461 0 : CPABORT("nexp_ppl > 4")
462 : END IF
463 :
464 88079211 : f = zetc/q
465 88079211 : cc = (pi/q)**1.5_dp*EXP(-t*f)
466 :
467 88079211 : IF (mmax >= 0) expder(0) = cc
468 322796280 : DO i = 1, mmax
469 322796280 : expder(i) = f*expder(i - 1)
470 : END DO
471 :
472 410875491 : DO i = 0, mmax
473 930402174 : DO j = 0, MIN(i, pmax)
474 519526683 : kp = j
475 519526683 : ke = i - j
476 842322963 : auxint(i) = auxint(i) + expder(ke)*polder(kp)*binomial(i, j)
477 : END DO
478 : END DO
479 :
480 88079211 : END SUBROUTINE ppl_aux
481 : ! **************************************************************************************************
482 : !> \brief ...
483 : !> \param auxint ...
484 : !> \param mmax ...
485 : !> \param t ...
486 : !> \param rho ...
487 : !> \param nexp ...
488 : !> \param cexp ...
489 : !> \param zetc ...
490 : ! **************************************************************************************************
491 224374 : SUBROUTINE ecploc_aux(auxint, mmax, t, rho, nexp, cexp, zetc)
492 : INTEGER, INTENT(IN) :: mmax
493 : REAL(KIND=dp), DIMENSION(0:mmax) :: auxint
494 : REAL(KIND=dp), INTENT(IN) :: t, rho
495 : INTEGER, INTENT(IN) :: nexp
496 : REAL(KIND=dp), INTENT(IN) :: cexp, zetc
497 :
498 : INTEGER :: i, j, ke, kf
499 : REAL(KIND=dp) :: c0, c1, cc, cval, fa, fr, q, ts
500 224374 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: expder, fdiff, funder, gfund
501 :
502 224374 : q = rho + zetc
503 224374 : fa = zetc/q
504 224374 : fr = rho/q
505 : !
506 1121870 : ALLOCATE (expder(0:mmax), funder(0:mmax + 1))
507 : !
508 240334 : SELECT CASE (nexp)
509 : CASE (0)
510 15960 : cval = 2.0_dp*cexp/SQRT(q)*pi**1.5_dp*EXP(-t*fa)
511 15960 : expder(0) = cval
512 45392 : DO i = 1, mmax
513 45392 : expder(i) = fa*expder(i - 1)
514 : END DO
515 15960 : ts = fr*t
516 31920 : ALLOCATE (gfund(0:mmax))
517 15960 : CALL gfun_values(mmax, ts, gfund)
518 :
519 15960 : funder(0) = gfund(0)
520 45392 : DO i = 1, mmax
521 29432 : funder(i) = 0.0_dp
522 127561 : DO j = 0, i
523 111601 : funder(i) = funder(i) + (-1)**j*binomial(i, j)*gfund(j)
524 : END DO
525 : END DO
526 :
527 15960 : DEALLOCATE (gfund)
528 45392 : DO i = 1, mmax
529 45392 : funder(i) = fr**i*funder(i)
530 : END DO
531 61352 : DO i = 0, mmax
532 159481 : DO j = 0, i
533 98129 : kf = j
534 98129 : ke = i - j
535 143521 : auxint(i) = auxint(i) + expder(ke)*funder(kf)*binomial(i, j)
536 : END DO
537 : END DO
538 : CASE (1)
539 49898 : cval = cexp*2._dp*pi/q*EXP(-t*fa)
540 49898 : expder(0) = cval
541 143982 : DO i = 1, mmax
542 143982 : expder(i) = fa*expder(i - 1)
543 : END DO
544 49898 : ts = fr*t
545 49898 : CALL fgamma(mmax, ts, funder)
546 143982 : DO i = 1, mmax
547 143982 : funder(i) = fr**i*funder(i)
548 : END DO
549 193880 : DO i = 0, mmax
550 518574 : DO j = 0, i
551 324694 : kf = j
552 324694 : ke = i - j
553 468676 : auxint(i) = auxint(i) + expder(ke)*funder(kf)*binomial(i, j)
554 : END DO
555 : END DO
556 : CASE (2)
557 134712 : cval = cexp*(pi/q)**1.5_dp*EXP(-t*fa)
558 134712 : expder(0) = cval
559 386632 : DO i = 1, mmax
560 386632 : expder(i) = fa*expder(i - 1)
561 : END DO
562 521344 : auxint(0:mmax) = auxint(0:mmax) + expder(0:mmax)
563 : CASE (3)
564 23804 : cval = 2.*pi*cexp/q**2*EXP(-t*fa)
565 23804 : expder(0) = cval
566 71008 : DO i = 1, mmax
567 71008 : expder(i) = fa*expder(i - 1)
568 : END DO
569 23804 : ts = fr*t
570 23804 : CALL fgamma(mmax + 1, ts, funder)
571 47608 : ALLOCATE (fdiff(0:mmax))
572 23804 : fdiff(0) = (1.0_dp + ts)*funder(0) - ts*funder(1)
573 71008 : DO i = 1, mmax
574 : fdiff(i) = fr**i*(-i*funder(i - 1) + (1.0_dp + ts)*funder(i) &
575 71008 : + i*funder(i) - ts*funder(i + 1))
576 : END DO
577 94812 : DO i = 0, mmax
578 263102 : DO j = 0, i
579 168290 : kf = j
580 168290 : ke = i - j
581 239298 : auxint(i) = auxint(i) + expder(ke)*fdiff(kf)*binomial(i, j)
582 : END DO
583 : END DO
584 23804 : DEALLOCATE (fdiff)
585 : CASE (4)
586 0 : cval = cexp/(4._dp*q**2)*(pi/q)**1.5_dp*EXP(-t*fa)
587 0 : expder(0) = cval
588 0 : DO i = 1, mmax
589 0 : expder(i) = fa*expder(i - 1)
590 : END DO
591 0 : c0 = 4._dp*rho/fa
592 0 : c1 = 6._dp*q + 4._dp*rho*t
593 0 : DO i = 0, mmax
594 0 : cc = -i*c0 + c1
595 0 : expder(i) = cc*expder(i)
596 : END DO
597 0 : auxint(0:mmax) = auxint(0:mmax) + expder(0:mmax)
598 : CASE DEFAULT
599 224374 : CPABORT("nexp out of range [1..4]")
600 : END SELECT
601 : !
602 224374 : DEALLOCATE (expder, funder)
603 :
604 224374 : END SUBROUTINE ecploc_aux
605 : ! **************************************************************************************************
606 :
607 : END MODULE ai_overlap_ppl
|
