Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 ai_overlap_debug, ONLY: init_os_overlap2,&
31 : os_overlap2
32 : USE gamma, ONLY: fgamma => fgamma_0
33 : USE gfun, ONLY: gfun_values
34 : USE kinds, ONLY: dp
35 : USE mathconstants, ONLY: pi
36 : USE mathlib, ONLY: binomial
37 : USE orbital_pointers, ONLY: indco,&
38 : ncoset
39 : #include "../base/base_uses.f90"
40 :
41 : IMPLICIT NONE
42 :
43 : PRIVATE
44 :
45 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap_ppl'
46 :
47 : ! *** Public subroutines ***
48 :
49 : PUBLIC :: ecploc_integral, ppl_integral, ppl_integral_ri
50 :
51 : CONTAINS
52 :
53 : ! **************************************************************************************************
54 : !> \brief Calculation of three-center overlap integrals <a|c|b> over
55 : !> Cartesian Gaussian functions for the local part of the Goedecker
56 : !> pseudopotential (GTH). c is a primitive Gaussian-type function
57 : !> with a set of even angular momentum indices.
58 : !>
59 : !> <a|V(ppl)|b> = <a| (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
60 : !> C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
61 : !> zetc = alpha**2/2
62 : !>
63 : !> \param la_max_set ...
64 : !> \param la_min_set ...
65 : !> \param npgfa ...
66 : !> \param rpgfa ...
67 : !> \param zeta ...
68 : !> \param lb_max_set ...
69 : !> \param lb_min_set ...
70 : !> \param npgfb ...
71 : !> \param rpgfb ...
72 : !> \param zetb ...
73 : !> \param nexp_ppl ...
74 : !> \param alpha_ppl ...
75 : !> \param nct_ppl ...
76 : !> \param cexp_ppl ...
77 : !> \param rpgfc ...
78 : !> \param rab ...
79 : !> \param dab ...
80 : !> \param rac ...
81 : !> \param dac ...
82 : !> \param rbc ...
83 : !> \param dbc ...
84 : !> \param vab ...
85 : !> \param s ...
86 : !> \param pab ...
87 : !> \param force_a ...
88 : !> \param force_b ...
89 : !> \param fs ...
90 : !> \param hab2 The derivative of the ppl integrals according to the weighting factors deltaR
91 : !> \param hab2_work ...
92 : !> \param deltaR Weighting factors for the derivatives wrt. nuclear positions
93 : !> \param iatom ...
94 : !> \param jatom ...
95 : !> \param katom ...
96 : !> \date May 2011
97 : !> \author Juerg Hutter
98 : !> \version 1.0
99 : !> \note Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
100 : ! **************************************************************************************************
101 34881833 : SUBROUTINE ppl_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
102 34881833 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
103 69763666 : rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
104 69763666 : hab2, hab2_work, deltaR, iatom, jatom, katom)
105 : INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
106 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
107 : INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
108 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
109 : INTEGER, INTENT(IN) :: nexp_ppl
110 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: alpha_ppl
111 : INTEGER, DIMENSION(:), INTENT(IN) :: nct_ppl
112 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: cexp_ppl
113 : REAL(KIND=dp), INTENT(IN) :: rpgfc
114 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
115 : REAL(KIND=dp), INTENT(IN) :: dab
116 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
117 : REAL(KIND=dp), INTENT(IN) :: dac
118 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rbc
119 : REAL(KIND=dp), INTENT(IN) :: dbc
120 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
121 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: s
122 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
123 : OPTIONAL :: pab
124 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a, force_b
125 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
126 : OPTIONAL :: fs, hab2, hab2_work
127 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
128 : OPTIONAL :: deltaR
129 : INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom, katom
130 :
131 : INTEGER :: iexp, ij, ipgf, jpgf, mmax, nexp
132 : REAL(KIND=dp) :: rho, sab, t, zetc
133 34881833 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: auxint
134 : REAL(KIND=dp), DIMENSION(3) :: pci
135 :
136 34881833 : IF (PRESENT(pab)) THEN
137 8983339 : CPASSERT(PRESENT(force_a))
138 8983339 : CPASSERT(PRESENT(force_b))
139 8983339 : CPASSERT(PRESENT(fs))
140 8983339 : mmax = la_max_set + lb_max_set + 2
141 8983339 : force_a(:) = 0.0_dp
142 8983339 : force_b(:) = 0.0_dp
143 25898494 : ELSE IF (PRESENT(hab2)) THEN
144 870 : mmax = la_max_set + lb_max_set + 2
145 : ELSE
146 25897624 : mmax = la_max_set + lb_max_set
147 : END IF
148 :
149 139527332 : ALLOCATE (auxint(0:mmax, npgfa*npgfb))
150 2327914262 : auxint = 0._dp
151 :
152 : ! *** Calculate auxiliary integrals ***
153 :
154 154497928 : DO ipgf = 1, npgfa
155 : ! *** Screening ***
156 119616095 : IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
157 238513995 : DO jpgf = 1, npgfb
158 : ! *** Screening ***
159 159011237 : IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
160 : (rpgfa(ipgf) + rpgfb(jpgf) < dab)) CYCLE
161 58183250 : ij = (ipgf - 1)*npgfb + jpgf
162 58183250 : rho = zeta(ipgf) + zetb(jpgf)
163 232733000 : pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
164 58183250 : sab = EXP(-(zeta(ipgf)*zetb(jpgf)/rho*dab*dab))
165 232733000 : t = rho*SUM(pci(:)*pci(:))
166 :
167 116630430 : DO iexp = 1, nexp_ppl
168 58447180 : nexp = nct_ppl(iexp)
169 58447180 : zetc = alpha_ppl(iexp)
170 116630430 : CALL ppl_aux(auxint(0:mmax, ij), mmax, t, rho, nexp, cexp_ppl(:, iexp), zetc)
171 : END DO
172 :
173 391110681 : auxint(0:mmax, ij) = sab*auxint(0:mmax, ij)
174 :
175 : END DO
176 : END DO
177 :
178 : CALL os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
179 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
180 : rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
181 : vab2=hab2, vab2_work=hab2_work, &
182 191320840 : deltaR=deltaR, iatom=iatom, jatom=jatom, katom=katom)
183 :
184 34881833 : DEALLOCATE (auxint)
185 :
186 34881833 : END SUBROUTINE ppl_integral
187 :
188 : ! **************************************************************************************************
189 : !> \brief Calculation of three-center potential integrals <a|V(r)|b> over
190 : !> Cartesian Gaussian functions for the local part of ECP
191 : !> pseudopotential. Multiple terms C1-4 are possible.
192 : !>
193 : !> <a|V(ecploc)|b> = <a| C1/r*exp(-a1*r**2) + C2*exp(-a2*r**2) + C3*r*exp(-a3*r**2) +
194 : !> C4*r**2*exp(-a4*r**2)|b>
195 : !>
196 : !> \param la_max_set ...
197 : !> \param la_min_set ...
198 : !> \param npgfa ...
199 : !> \param rpgfa ...
200 : !> \param zeta ...
201 : !> \param lb_max_set ...
202 : !> \param lb_min_set ...
203 : !> \param npgfb ...
204 : !> \param rpgfb ...
205 : !> \param zetb ...
206 : !> \param nexp_ppl ...
207 : !> \param alpha_ppl ...
208 : !> \param nct_ppl ...
209 : !> \param cexp_ppl ...
210 : !> \param rpgfc ...
211 : !> \param rab ...
212 : !> \param dab ...
213 : !> \param rac ...
214 : !> \param dac ...
215 : !> \param rbc ...
216 : !> \param dbc ...
217 : !> \param vab ...
218 : !> \param s ...
219 : !> \param pab ...
220 : !> \param force_a ...
221 : !> \param force_b ...
222 : !> \param fs ...
223 : !> \param hab2 The derivative of the ppl integrals according to the weighting factors deltaR
224 : !> \param hab2_work ...
225 : !> \param deltaR Weighting factors for the derivatives wrt. nuclear positions
226 : !> \param iatom ...
227 : !> \param jatom ...
228 : !> \param katom ...
229 : !> \date 2025
230 : !> \author Juerg Hutter
231 : !> \version 1.0
232 : ! **************************************************************************************************
233 12600 : SUBROUTINE ecploc_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
234 12600 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
235 12600 : nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
236 25200 : rab, dab, rac, dac, rbc, dbc, vab, s, pab, &
237 25200 : force_a, force_b, fs, hab2, hab2_work, &
238 12600 : deltaR, iatom, jatom, katom)
239 : INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
240 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
241 : INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
242 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
243 : INTEGER, INTENT(IN) :: nexp_ppl
244 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: alpha_ppl
245 : INTEGER, DIMENSION(:), INTENT(IN) :: nct_ppl
246 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: cexp_ppl
247 : REAL(KIND=dp), INTENT(IN) :: rpgfc
248 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
249 : REAL(KIND=dp), INTENT(IN) :: dab
250 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
251 : REAL(KIND=dp), INTENT(IN) :: dac
252 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rbc
253 : REAL(KIND=dp), INTENT(IN) :: dbc
254 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
255 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: s
256 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
257 : OPTIONAL :: pab
258 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a, force_b
259 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
260 : OPTIONAL :: fs, hab2, hab2_work
261 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
262 : OPTIONAL :: deltaR
263 : INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom, katom
264 :
265 : INTEGER :: iexp, ij, ipgf, jpgf, mmax, nexp
266 : REAL(KIND=dp) :: rho, sab, t, zetc
267 12600 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: auxint
268 : REAL(KIND=dp), DIMENSION(3) :: pci
269 :
270 12600 : IF (PRESENT(pab)) THEN
271 2888 : CPASSERT(PRESENT(force_a))
272 2888 : CPASSERT(PRESENT(force_b))
273 2888 : CPASSERT(PRESENT(fs))
274 2888 : mmax = la_max_set + lb_max_set + 2
275 2888 : force_a(:) = 0.0_dp
276 2888 : force_b(:) = 0.0_dp
277 9712 : ELSE IF (PRESENT(hab2)) THEN
278 0 : mmax = la_max_set + lb_max_set + 2
279 : ELSE
280 9712 : mmax = la_max_set + lb_max_set
281 : END IF
282 :
283 50400 : ALLOCATE (auxint(0:mmax, npgfa*npgfb))
284 126232 : auxint = 0._dp
285 :
286 : ! *** Calculate auxiliary integrals ***
287 :
288 31259 : DO ipgf = 1, npgfa
289 : ! *** Screening ***
290 18659 : IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
291 52371 : DO jpgf = 1, npgfb
292 : ! *** Screening ***
293 24316 : IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
294 : (rpgfa(ipgf) + rpgfb(jpgf) < dab)) CYCLE
295 19271 : ij = (ipgf - 1)*npgfb + jpgf
296 19271 : rho = zeta(ipgf) + zetb(jpgf)
297 77084 : pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
298 19271 : sab = EXP(-(zeta(ipgf)*zetb(jpgf)/rho*dab*dab))
299 77084 : t = rho*SUM(pci(:)*pci(:))
300 :
301 48993 : DO iexp = 1, nexp_ppl
302 29722 : nexp = nct_ppl(iexp)
303 29722 : zetc = alpha_ppl(iexp)
304 48993 : CALL ecploc_aux(auxint(0:mmax, ij), mmax, t, rho, nexp, cexp_ppl(1, iexp), zetc)
305 : END DO
306 :
307 97574 : auxint(0:mmax, ij) = sab*auxint(0:mmax, ij)
308 :
309 : END DO
310 : END DO
311 :
312 : CALL os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
313 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
314 : rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
315 : vab2=hab2, vab2_work=hab2_work, &
316 69824 : deltaR=deltaR, iatom=iatom, jatom=jatom, katom=katom)
317 :
318 12600 : DEALLOCATE (auxint)
319 :
320 12600 : END SUBROUTINE ecploc_integral
321 : ! **************************************************************************************************
322 : !> \brief Calculation of two-center overlap integrals <a|c> over
323 : !> Cartesian Gaussian functions for the local part of the Goedecker
324 : !> pseudopotential (GTH). c is a primitive Gaussian-type function
325 : !> with a set of even angular momentum indices.
326 : !>
327 : !> <a|V(ppl)|b> = <a| (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
328 : !> C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
329 : !> zetc = alpha**2/2
330 : !>
331 : !> \param la_max_set ...
332 : !> \param la_min_set ...
333 : !> \param npgfa ...
334 : !> \param rpgfa ...
335 : !> \param zeta ...
336 : !> \param nexp_ppl ...
337 : !> \param alpha_ppl ...
338 : !> \param nct_ppl ...
339 : !> \param cexp_ppl ...
340 : !> \param rpgfc ...
341 : !> \param rac ...
342 : !> \param dac ...
343 : !> \param va ...
344 : !> \param dva ...
345 : !> \date December 2017
346 : !> \author Juerg Hutter
347 : !> \version 1.0
348 : ! **************************************************************************************************
349 328 : SUBROUTINE ppl_integral_ri(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
350 328 : nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
351 656 : rac, dac, va, dva)
352 : INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
353 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
354 : INTEGER, INTENT(IN) :: nexp_ppl
355 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: alpha_ppl
356 : INTEGER, DIMENSION(:), INTENT(IN) :: nct_ppl
357 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: cexp_ppl
358 : REAL(KIND=dp), INTENT(IN) :: rpgfc
359 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
360 : REAL(KIND=dp), INTENT(IN) :: dac
361 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: va
362 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
363 : OPTIONAL :: dva
364 :
365 : INTEGER :: i, iexp, ipgf, iw, mmax, na, nexp
366 : INTEGER, DIMENSION(3) :: ani, anm, anp
367 : LOGICAL :: debug
368 : REAL(KIND=dp) :: oint, oref, rho, t, zetc
369 328 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: auxint
370 : REAL(KIND=dp), DIMENSION(3) :: doint, doref
371 :
372 328 : debug = .FALSE.
373 :
374 328 : IF (PRESENT(dva)) THEN
375 164 : mmax = la_max_set + 1
376 : ELSE
377 164 : mmax = la_max_set
378 : END IF
379 :
380 1312 : ALLOCATE (auxint(0:mmax, npgfa))
381 1588 : auxint = 0._dp
382 :
383 : ! *** Calculate auxiliary integrals ***
384 656 : DO ipgf = 1, npgfa
385 328 : IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
386 328 : rho = zeta(ipgf)
387 328 : t = rho*dac*dac
388 :
389 984 : DO iexp = 1, nexp_ppl
390 328 : nexp = nct_ppl(iexp)
391 328 : zetc = alpha_ppl(iexp)
392 656 : CALL ppl_aux(auxint(0:mmax, ipgf), mmax, t, rho, nexp, cexp_ppl(:, iexp), zetc)
393 : END DO
394 :
395 : END DO
396 :
397 328 : IF (PRESENT(dva)) THEN
398 : CALL os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
399 164 : auxint, rpgfc, rac, dac, va, dva)
400 : ELSE
401 : CALL os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
402 164 : auxint, rpgfc, rac, dac, va)
403 : END IF
404 :
405 328 : DEALLOCATE (auxint)
406 :
407 : IF (debug) THEN
408 : iw = 6
409 : na = 0
410 : DO ipgf = 1, npgfa
411 : IF (rpgfa(ipgf) + rpgfc < dac) THEN
412 : na = na + ncoset(la_max_set)
413 : CYCLE
414 : END IF
415 : rho = zeta(ipgf)
416 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
417 : oref = va(na + i)
418 : ani(1:3) = indco(1:3, i)
419 : oint = ppl_ri_test(rho, ani, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
420 : ! test
421 : IF (ABS(oint - oref) > 1.0e-12_dp) THEN
422 : WRITE (iw, '(A,3i2,i5,F10.4,2G24.12)') "PPL int error ", ani, la_max_set, dac, oint, oref
423 : END IF
424 : IF (PRESENT(dva)) THEN
425 : anp = ani + [1, 0, 0]
426 : anm = ani - [1, 0, 0]
427 : doint(1) = 2._dp*rho*ppl_ri_test(rho, anp, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) &
428 : - ani(1)*ppl_ri_test(rho, anm, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
429 : anp = ani + [0, 1, 0]
430 : anm = ani - [0, 1, 0]
431 : doint(2) = 2._dp*rho*ppl_ri_test(rho, anp, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) &
432 : - ani(2)*ppl_ri_test(rho, anm, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
433 : anp = ani + [0, 0, 1]
434 : anm = ani - [0, 0, 1]
435 : doint(3) = 2._dp*rho*ppl_ri_test(rho, anp, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) &
436 : - ani(3)*ppl_ri_test(rho, anm, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
437 : doref(1:3) = dva(na + i, 1:3)
438 : IF (ANY(ABS(doint - doref) > 1.0e-6_dp)) THEN
439 : WRITE (iw, '(A,3i2,i5,F10.4,2G24.12)') " PPL dint error ", &
440 : ani, la_max_set, dac, SUM(ABS(doint)), SUM(ABS(doref))
441 : END IF
442 : END IF
443 : END DO
444 : na = na + ncoset(la_max_set)
445 : END DO
446 : END IF
447 :
448 328 : END SUBROUTINE ppl_integral_ri
449 :
450 : ! **************************************************************************************************
451 : !> \brief ...
452 : !> \param rho ...
453 : !> \param ani ...
454 : !> \param rac ...
455 : !> \param nexp_ppl ...
456 : !> \param nct_ppl ...
457 : !> \param alpha_ppl ...
458 : !> \param cexp_ppl ...
459 : !> \return ...
460 : ! **************************************************************************************************
461 0 : FUNCTION ppl_ri_test(rho, ani, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) RESULT(oint)
462 : REAL(KIND=dp), INTENT(IN) :: rho
463 : INTEGER, DIMENSION(3), INTENT(IN) :: ani
464 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
465 : INTEGER, INTENT(IN) :: nexp_ppl
466 : INTEGER, DIMENSION(:), INTENT(IN) :: nct_ppl
467 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: alpha_ppl
468 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: cexp_ppl
469 : REAL(KIND=dp) :: oint
470 :
471 : INTEGER :: iexp, nexp, ni
472 : REAL(KIND=dp) :: cn, zetc
473 : REAL(KIND=dp), DIMENSION(3) :: ra
474 :
475 0 : oint = 0.0_dp
476 0 : ra = 0.0_dp
477 0 : DO iexp = 1, nexp_ppl
478 0 : nexp = nct_ppl(iexp)
479 0 : zetc = alpha_ppl(iexp)
480 0 : CALL init_os_overlap2(rho, zetc, ra, -rac)
481 0 : DO ni = 1, nexp
482 0 : cn = cexp_ppl(ni, iexp)
483 0 : SELECT CASE (ni)
484 : CASE (1)
485 0 : oint = oint + cn*os_overlap2(ani, [0, 0, 0])
486 : CASE (2)
487 0 : oint = oint + cn*os_overlap2(ani, [2, 0, 0])
488 0 : oint = oint + cn*os_overlap2(ani, [0, 2, 0])
489 0 : oint = oint + cn*os_overlap2(ani, [0, 0, 2])
490 : CASE (3)
491 0 : oint = oint + cn*os_overlap2(ani, [4, 0, 0])
492 0 : oint = oint + cn*os_overlap2(ani, [0, 4, 0])
493 0 : oint = oint + cn*os_overlap2(ani, [0, 0, 4])
494 0 : oint = oint + 2.0_dp*cn*os_overlap2(ani, [2, 2, 0])
495 0 : oint = oint + 2.0_dp*cn*os_overlap2(ani, [0, 2, 2])
496 0 : oint = oint + 2.0_dp*cn*os_overlap2(ani, [2, 0, 2])
497 : CASE (4)
498 0 : oint = oint + cn*os_overlap2(ani, [6, 0, 0])
499 0 : oint = oint + cn*os_overlap2(ani, [0, 6, 0])
500 0 : oint = oint + cn*os_overlap2(ani, [0, 0, 6])
501 0 : oint = oint + 3.0_dp*cn*os_overlap2(ani, [4, 2, 0])
502 0 : oint = oint + 3.0_dp*cn*os_overlap2(ani, [4, 0, 2])
503 0 : oint = oint + 3.0_dp*cn*os_overlap2(ani, [2, 4, 0])
504 0 : oint = oint + 3.0_dp*cn*os_overlap2(ani, [0, 4, 2])
505 0 : oint = oint + 3.0_dp*cn*os_overlap2(ani, [2, 0, 4])
506 0 : oint = oint + 3.0_dp*cn*os_overlap2(ani, [0, 2, 4])
507 0 : oint = oint + 6.0_dp*cn*os_overlap2(ani, [2, 2, 2])
508 : CASE DEFAULT
509 0 : CPABORT("OVERLAP_PPL")
510 : END SELECT
511 : END DO
512 : END DO
513 :
514 0 : END FUNCTION ppl_ri_test
515 :
516 : ! **************************************************************************************************
517 : !> \brief ...
518 : !> \param auxint ...
519 : !> \param mmax ...
520 : !> \param t ...
521 : !> \param rho ...
522 : !> \param nexp_ppl ...
523 : !> \param cexp_ppl ...
524 : !> \param zetc ...
525 : ! **************************************************************************************************
526 58447508 : SUBROUTINE ppl_aux(auxint, mmax, t, rho, nexp_ppl, cexp_ppl, zetc)
527 : INTEGER, INTENT(IN) :: mmax
528 : REAL(KIND=dp), DIMENSION(0:mmax) :: auxint
529 : REAL(KIND=dp), INTENT(IN) :: t, rho
530 : INTEGER, INTENT(IN) :: nexp_ppl
531 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cexp_ppl
532 : REAL(KIND=dp), INTENT(IN) :: zetc
533 :
534 : INTEGER :: i, j, ke, kp, pmax
535 : REAL(KIND=dp) :: a2, a3, a4, cc, f, q, q2, q4, q6, rho2, &
536 : rho3, t2, t3
537 : REAL(KIND=dp), DIMENSION(0:6) :: polder
538 116895016 : REAL(KIND=dp), DIMENSION(0:mmax) :: expder
539 :
540 58447508 : CPASSERT(nexp_ppl > 0)
541 58447508 : q = rho + zetc
542 58447508 : polder = 0._dp
543 58447508 : pmax = 0
544 58447508 : IF (nexp_ppl > 0) THEN
545 58447508 : polder(0) = polder(0) + cexp_ppl(1)
546 58447508 : pmax = 0
547 : END IF
548 58447508 : IF (nexp_ppl > 1) THEN
549 44154892 : q2 = q*q
550 44154892 : a2 = 0.5_dp/q2*cexp_ppl(2)
551 44154892 : polder(0) = polder(0) + a2*(2._dp*rho*t + 3._dp*q)
552 44154892 : polder(1) = polder(1) - a2*2._dp*rho
553 44154892 : pmax = 1
554 : END IF
555 58447508 : IF (nexp_ppl > 2) THEN
556 944146 : q4 = q2*q2
557 944146 : rho2 = rho*rho
558 944146 : t2 = t*t
559 944146 : a3 = 0.25_dp/q4*cexp_ppl(3)
560 944146 : polder(0) = polder(0) + a3*(4._dp*rho2*t2 + 20._dp*rho*t*q + 15._dp*q2)
561 944146 : polder(1) = polder(1) - a3*(8._dp*rho2*t + 20._dp*rho*q)
562 944146 : polder(2) = polder(2) + a3*8._dp*rho2
563 944146 : pmax = 2
564 : END IF
565 58447508 : IF (nexp_ppl > 3) THEN
566 944146 : q6 = q4*q2
567 944146 : rho3 = rho2*rho
568 944146 : t3 = t2*t
569 944146 : a4 = 0.125_dp/q6*cexp_ppl(4)
570 944146 : polder(0) = polder(0) + a4*(8._dp*rho3*t3 + 84._dp*rho2*t2*q + 210._dp*rho*t*q2 + 105._dp*q*q2)
571 944146 : polder(1) = polder(1) - a4*(24._dp*rho3*t2 + 168._dp*rho2*t*q + 210._dp*rho*q2)
572 944146 : polder(2) = polder(2) + a4*(48._dp*rho3*t + 168._dp*rho2*q)
573 944146 : polder(3) = polder(3) - a4*48_dp*rho3
574 944146 : pmax = 3
575 : END IF
576 58447508 : IF (nexp_ppl > 4) THEN
577 0 : CPABORT("nexp_ppl > 4")
578 : END IF
579 :
580 58447508 : f = zetc/q
581 58447508 : cc = (pi/q)**1.5_dp*EXP(-t*f)
582 :
583 58447508 : IF (mmax >= 0) expder(0) = cc
584 213051552 : DO i = 1, mmax
585 213051552 : expder(i) = f*expder(i - 1)
586 : END DO
587 :
588 271499060 : DO i = 0, mmax
589 601263583 : DO j = 0, MIN(i, pmax)
590 329764523 : kp = j
591 329764523 : ke = i - j
592 542816075 : auxint(i) = auxint(i) + expder(ke)*polder(kp)*binomial(i, j)
593 : END DO
594 : END DO
595 :
596 58447508 : END SUBROUTINE ppl_aux
597 : ! **************************************************************************************************
598 : !> \brief ...
599 : !> \param auxint ...
600 : !> \param mmax ...
601 : !> \param t ...
602 : !> \param rho ...
603 : !> \param nexp ...
604 : !> \param cexp ...
605 : !> \param zetc ...
606 : ! **************************************************************************************************
607 29722 : SUBROUTINE ecploc_aux(auxint, mmax, t, rho, nexp, cexp, zetc)
608 : INTEGER, INTENT(IN) :: mmax
609 : REAL(KIND=dp), DIMENSION(0:mmax) :: auxint
610 : REAL(KIND=dp), INTENT(IN) :: t, rho
611 : INTEGER, INTENT(IN) :: nexp
612 : REAL(KIND=dp), INTENT(IN) :: cexp, zetc
613 :
614 : INTEGER :: i, j, ke, kf
615 : REAL(KIND=dp) :: c0, c1, cc, cval, fa, fr, q, ts
616 29722 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: expder, fdiff, funder, gfund
617 :
618 29722 : q = rho + zetc
619 29722 : fa = zetc/q
620 29722 : fr = rho/q
621 : !
622 148610 : ALLOCATE (expder(0:mmax), funder(0:mmax + 1))
623 : !
624 30262 : SELECT CASE (nexp)
625 : CASE (0)
626 540 : cval = 2.0_dp*cexp/SQRT(q)*pi**1.5_dp*EXP(-t*fa)
627 540 : expder(0) = cval
628 1296 : DO i = 1, mmax
629 1296 : expder(i) = fa*expder(i - 1)
630 : END DO
631 540 : ts = fr*t
632 1080 : ALLOCATE (gfund(0:mmax))
633 540 : CALL gfun_values(mmax, ts, gfund)
634 :
635 540 : funder(0) = gfund(0)
636 1296 : DO i = 1, mmax
637 756 : funder(i) = 0.0_dp
638 3267 : DO j = 0, i
639 2727 : funder(i) = funder(i) + (-1)**j*binomial(i, j)*gfund(j)
640 : END DO
641 : END DO
642 :
643 540 : DEALLOCATE (gfund)
644 1296 : DO i = 1, mmax
645 1296 : funder(i) = fr**i*funder(i)
646 : END DO
647 1836 : DO i = 0, mmax
648 4347 : DO j = 0, i
649 2511 : kf = j
650 2511 : ke = i - j
651 3807 : auxint(i) = auxint(i) + expder(ke)*funder(kf)*binomial(i, j)
652 : END DO
653 : END DO
654 : CASE (1)
655 1690 : cval = cexp*2._dp*pi/q*EXP(-t*fa)
656 1690 : expder(0) = cval
657 4286 : DO i = 1, mmax
658 4286 : expder(i) = fa*expder(i - 1)
659 : END DO
660 1690 : ts = fr*t
661 1690 : CALL fgamma(mmax, ts, funder)
662 4286 : DO i = 1, mmax
663 4286 : funder(i) = fr**i*funder(i)
664 : END DO
665 5976 : DO i = 0, mmax
666 14734 : DO j = 0, i
667 8758 : kf = j
668 8758 : ke = i - j
669 13044 : auxint(i) = auxint(i) + expder(ke)*funder(kf)*binomial(i, j)
670 : END DO
671 : END DO
672 : CASE (2)
673 26882 : cval = cexp*(pi/q)**1.5_dp*EXP(-t*fa)
674 26882 : expder(0) = cval
675 79886 : DO i = 1, mmax
676 79886 : expder(i) = fa*expder(i - 1)
677 : END DO
678 106768 : auxint(0:mmax) = auxint(0:mmax) + expder(0:mmax)
679 : CASE (3)
680 610 : cval = 2.*pi*cexp/q**2*EXP(-t*fa)
681 610 : expder(0) = cval
682 1694 : DO i = 1, mmax
683 1694 : expder(i) = fa*expder(i - 1)
684 : END DO
685 610 : ts = fr*t
686 610 : CALL fgamma(mmax + 1, ts, funder)
687 1220 : ALLOCATE (fdiff(0:mmax))
688 610 : fdiff(0) = (1.0_dp + ts)*funder(0) - ts*funder(1)
689 1694 : DO i = 1, mmax
690 : fdiff(i) = fr**i*(-i*funder(i - 1) + (1.0_dp + ts)*funder(i) &
691 1694 : + i*funder(i) - ts*funder(i + 1))
692 : END DO
693 2304 : DO i = 0, mmax
694 6040 : DO j = 0, i
695 3736 : kf = j
696 3736 : ke = i - j
697 5430 : auxint(i) = auxint(i) + expder(ke)*fdiff(kf)*binomial(i, j)
698 : END DO
699 : END DO
700 610 : DEALLOCATE (fdiff)
701 : CASE (4)
702 0 : cval = cexp/(4._dp*q**2)*(pi/q)**1.5_dp*EXP(-t*fa)
703 0 : expder(0) = cval
704 0 : DO i = 1, mmax
705 0 : expder(i) = fa*expder(i - 1)
706 : END DO
707 0 : c0 = 4._dp*rho/fa
708 0 : c1 = 6._dp*q + 4._dp*rho*t
709 0 : DO i = 0, mmax
710 0 : cc = -i*c0 + c1
711 0 : expder(i) = cc*expder(i)
712 : END DO
713 0 : auxint(0:mmax) = auxint(0:mmax) + expder(0:mmax)
714 : CASE DEFAULT
715 29722 : CPABORT("nexp out of range [1..4]")
716 : END SELECT
717 : !
718 29722 : DEALLOCATE (expder, funder)
719 :
720 29722 : END SUBROUTINE ecploc_aux
721 : ! **************************************************************************************************
722 :
723 : END MODULE ai_overlap_ppl
|
