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