(git:f2099e5)
Loading...
Searching...
No Matches
ai_overlap_ppl.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!> \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! **************************************************************************************************
27MODULE 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
48CONTAINS
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 SUBROUTINE ppl_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
99 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
100 rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
101 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: auxint
131 REAL(kind=dp), DIMENSION(3) :: pci
132
133 IF (PRESENT(pab)) THEN
134 cpassert(PRESENT(force_a))
135 cpassert(PRESENT(force_b))
136 cpassert(PRESENT(fs))
137 mmax = la_max_set + lb_max_set + 2
138 force_a(:) = 0.0_dp
139 force_b(:) = 0.0_dp
140 ELSE IF (PRESENT(hab2)) THEN
141 mmax = la_max_set + lb_max_set + 2
142 ELSE
143 mmax = la_max_set + lb_max_set
144 END IF
145
146 ALLOCATE (auxint(0:mmax, npgfa*npgfb))
147 auxint = 0._dp
148
149 ! *** Calculate auxiliary integrals ***
150
151 DO ipgf = 1, npgfa
152 ! *** Screening ***
153 IF (rpgfa(ipgf) + rpgfc < dac) cycle
154 DO jpgf = 1, npgfb
155 ! *** Screening ***
156 IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
157 (rpgfa(ipgf) + rpgfb(jpgf) < dab)) cycle
158 ij = (ipgf - 1)*npgfb + jpgf
159 rho = zeta(ipgf) + zetb(jpgf)
160 pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
161 sab = exp(-(zeta(ipgf)*zetb(jpgf)/rho*dab*dab))
162 t = rho*sum(pci(:)*pci(:))
163
164 DO iexp = 1, nexp_ppl
165 nexp = nct_ppl(iexp)
166 zetc = alpha_ppl(iexp)
167 CALL ppl_aux(auxint(0:mmax, ij), mmax, t, rho, nexp, cexp_ppl(:, iexp), zetc)
168 END DO
169
170 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 deltar=deltar, iatom=iatom, jatom=jatom, katom=katom)
180
181 DEALLOCATE (auxint)
182
183 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 SUBROUTINE ecploc_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
231 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
232 nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
233 rab, dab, rac, dac, rbc, dbc, vab, s, pab, &
234 force_a, force_b, fs, hab2, hab2_work, &
235 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: auxint
265 REAL(kind=dp), DIMENSION(3) :: pci
266
267 IF (PRESENT(pab)) THEN
268 cpassert(PRESENT(force_a))
269 cpassert(PRESENT(force_b))
270 cpassert(PRESENT(fs))
271 mmax = la_max_set + lb_max_set + 2
272 force_a(:) = 0.0_dp
273 force_b(:) = 0.0_dp
274 ELSE IF (PRESENT(hab2)) THEN
275 mmax = la_max_set + lb_max_set + 2
276 ELSE
277 mmax = la_max_set + lb_max_set
278 END IF
279
280 ALLOCATE (auxint(0:mmax, npgfa*npgfb))
281 auxint = 0._dp
282
283 ! *** Calculate auxiliary integrals ***
284
285 DO ipgf = 1, npgfa
286 ! *** Screening ***
287 IF (rpgfa(ipgf) + rpgfc < dac) cycle
288 DO jpgf = 1, npgfb
289 ! *** Screening ***
290 IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
291 (rpgfa(ipgf) + rpgfb(jpgf) < dab)) cycle
292 ij = (ipgf - 1)*npgfb + jpgf
293 rho = zeta(ipgf) + zetb(jpgf)
294 pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
295 sab = exp(-(zeta(ipgf)*zetb(jpgf)/rho*dab*dab))
296 t = rho*sum(pci(:)*pci(:))
297
298 DO iexp = 1, nexp_ppl
299 nexp = nct_ppl(iexp)
300 zetc = alpha_ppl(iexp)
301 CALL ecploc_aux(auxint(0:mmax, ij), mmax, t, rho, nexp, cexp_ppl(1, iexp), zetc)
302 END DO
303
304 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 deltar=deltar, iatom=iatom, jatom=jatom, katom=katom)
314
315 DEALLOCATE (auxint)
316
317 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 SUBROUTINE ppl_integral_ri(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
347 nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
348 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 IF (PRESENT(dva)) THEN
367 mmax = la_max_set + 1
368 ELSE
369 mmax = la_max_set
370 END IF
371
372 ALLOCATE (auxint(0:mmax, npgfa))
373 auxint = 0._dp
374
375 ! *** Calculate auxiliary integrals ***
376 DO ipgf = 1, npgfa
377 IF (rpgfa(ipgf) + rpgfc < dac) cycle
378 rho = zeta(ipgf)
379 t = rho*dac*dac
380
381 DO iexp = 1, nexp_ppl
382 nexp = nct_ppl(iexp)
383 zetc = alpha_ppl(iexp)
384 CALL ppl_aux(auxint(0:mmax, ipgf), mmax, t, rho, nexp, cexp_ppl(:, iexp), zetc)
385 END DO
386
387 END DO
388
389 IF (PRESENT(dva)) THEN
390 CALL os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
391 auxint, rpgfc, rac, dac, va, dva)
392 ELSE
393 CALL os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
394 auxint, rpgfc, rac, dac, va)
395 END IF
396
397 DEALLOCATE (auxint)
398
399 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 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 REAL(kind=dp), DIMENSION(0:mmax) :: expder
424
425 cpassert(nexp_ppl > 0)
426 q = rho + zetc
427 polder = 0._dp
428 pmax = 0
429 IF (nexp_ppl > 0) THEN
430 polder(0) = polder(0) + cexp_ppl(1)
431 pmax = 0
432 END IF
433 IF (nexp_ppl > 1) THEN
434 q2 = q*q
435 a2 = 0.5_dp/q2*cexp_ppl(2)
436 polder(0) = polder(0) + a2*(2._dp*rho*t + 3._dp*q)
437 polder(1) = polder(1) - a2*2._dp*rho
438 pmax = 1
439 END IF
440 IF (nexp_ppl > 2) THEN
441 q4 = q2*q2
442 rho2 = rho*rho
443 t2 = t*t
444 a3 = 0.25_dp/q4*cexp_ppl(3)
445 polder(0) = polder(0) + a3*(4._dp*rho2*t2 + 20._dp*rho*t*q + 15._dp*q2)
446 polder(1) = polder(1) - a3*(8._dp*rho2*t + 20._dp*rho*q)
447 polder(2) = polder(2) + a3*8._dp*rho2
448 pmax = 2
449 END IF
450 IF (nexp_ppl > 3) THEN
451 q6 = q4*q2
452 rho3 = rho2*rho
453 t3 = t2*t
454 a4 = 0.125_dp/q6*cexp_ppl(4)
455 polder(0) = polder(0) + a4*(8._dp*rho3*t3 + 84._dp*rho2*t2*q + 210._dp*rho*t*q2 + 105._dp*q*q2)
456 polder(1) = polder(1) - a4*(24._dp*rho3*t2 + 168._dp*rho2*t*q + 210._dp*rho*q2)
457 polder(2) = polder(2) + a4*(48._dp*rho3*t + 168._dp*rho2*q)
458 polder(3) = polder(3) - a4*48_dp*rho3
459 pmax = 3
460 END IF
461 IF (nexp_ppl > 4) THEN
462 cpabort("nexp_ppl > 4")
463 END IF
464
465 f = zetc/q
466 cc = (pi/q)**1.5_dp*exp(-t*f)
467
468 IF (mmax >= 0) expder(0) = cc
469 DO i = 1, mmax
470 expder(i) = f*expder(i - 1)
471 END DO
472
473 DO i = 0, mmax
474 DO j = 0, min(i, pmax)
475 kp = j
476 ke = i - j
477 auxint(i) = auxint(i) + expder(ke)*polder(kp)*binomial(i, j)
478 END DO
479 END DO
480
481 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: expder, fdiff, funder, gfund
502
503 q = rho + zetc
504 fa = zetc/q
505 fr = rho/q
506 !
507 ALLOCATE (expder(0:mmax), funder(0:mmax + 1))
508 !
509 SELECT CASE (nexp)
510 CASE (0)
511 cval = 2.0_dp*cexp/sqrt(q)*pi**1.5_dp*exp(-t*fa)
512 expder(0) = cval
513 DO i = 1, mmax
514 expder(i) = fa*expder(i - 1)
515 END DO
516 ts = fr*t
517 ALLOCATE (gfund(0:mmax))
518 CALL gfun_values(mmax, ts, gfund)
519
520 funder(0) = gfund(0)
521 DO i = 1, mmax
522 funder(i) = 0.0_dp
523 DO j = 0, i
524 funder(i) = funder(i) + (-1)**j*binomial(i, j)*gfund(j)
525 END DO
526 END DO
527
528 DEALLOCATE (gfund)
529 DO i = 1, mmax
530 funder(i) = fr**i*funder(i)
531 END DO
532 DO i = 0, mmax
533 DO j = 0, i
534 kf = j
535 ke = i - j
536 auxint(i) = auxint(i) + expder(ke)*funder(kf)*binomial(i, j)
537 END DO
538 END DO
539 CASE (1)
540 cval = cexp*2._dp*pi/q*exp(-t*fa)
541 expder(0) = cval
542 DO i = 1, mmax
543 expder(i) = fa*expder(i - 1)
544 END DO
545 ts = fr*t
546 CALL fgamma(mmax, ts, funder)
547 DO i = 1, mmax
548 funder(i) = fr**i*funder(i)
549 END DO
550 DO i = 0, mmax
551 DO j = 0, i
552 kf = j
553 ke = i - j
554 auxint(i) = auxint(i) + expder(ke)*funder(kf)*binomial(i, j)
555 END DO
556 END DO
557 CASE (2)
558 cval = cexp*(pi/q)**1.5_dp*exp(-t*fa)
559 expder(0) = cval
560 DO i = 1, mmax
561 expder(i) = fa*expder(i - 1)
562 END DO
563 auxint(0:mmax) = auxint(0:mmax) + expder(0:mmax)
564 CASE (3)
565 cval = twopi*cexp/q**2*exp(-t*fa)
566 expder(0) = cval
567 DO i = 1, mmax
568 expder(i) = fa*expder(i - 1)
569 END DO
570 ts = fr*t
571 CALL fgamma(mmax + 1, ts, funder)
572 ALLOCATE (fdiff(0:mmax))
573 fdiff(0) = (1.0_dp + ts)*funder(0) - ts*funder(1)
574 DO i = 1, mmax
575 fdiff(i) = fr**i*(-i*funder(i - 1) + (1.0_dp + ts)*funder(i) &
576 + i*funder(i) - ts*funder(i + 1))
577 END DO
578 DO i = 0, mmax
579 DO j = 0, i
580 kf = j
581 ke = i - j
582 auxint(i) = auxint(i) + expder(ke)*fdiff(kf)*binomial(i, j)
583 END DO
584 END DO
585 DEALLOCATE (fdiff)
586 CASE (4)
587 cval = cexp/(4._dp*q**2)*(pi/q)**1.5_dp*exp(-t*fa)
588 expder(0) = cval
589 DO i = 1, mmax
590 expder(i) = fa*expder(i - 1)
591 END DO
592 c0 = 4._dp*rho/fa
593 c1 = 6._dp*q + 4._dp*rho*t
594 DO i = 0, mmax
595 cc = -i*c0 + c1
596 expder(i) = cc*expder(i)
597 END DO
598 auxint(0:mmax) = auxint(0:mmax) + expder(0:mmax)
599 CASE DEFAULT
600 cpabort("nexp out of range [1..4]")
601 END SELECT
602 !
603 DEALLOCATE (expder, funder)
604
605 END SUBROUTINE ecploc_aux
606! **************************************************************************************************
607
608END MODULE ai_overlap_ppl
ai_oneelectron
Calculation of general three-center integrals over Cartesian Gaussian-type functions and a spherical ...
Definition ai_oneelectron.F:34
ai_oneelectron::os_3center
subroutine, public os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, vab2, vab2_work, deltar, iatom, jatom, katom)
Calculation of three-center integrals <a|c|b> over Cartesian Gaussian functions and a spherical poten...
Definition ai_oneelectron.F:96
ai_oneelectron::os_2center
subroutine, public os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, auxint, rpgfc, rac, dac, va, dva)
Calculation of two-center integrals <a|c> over Cartesian Gaussian functions and a spherical potential...
Definition ai_oneelectron.F:609
ai_overlap_ppl
Calculation of three-center overlap integrals over Cartesian Gaussian-type functions for the second t...
Definition ai_overlap_ppl.F:27
ai_overlap_ppl::ppl_integral
subroutine, public ppl_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, hab2, hab2_work, deltar, iatom, jatom, katom)
Calculation of three-center overlap integrals <a|c|b> over Cartesian Gaussian functions for the local...
Definition ai_overlap_ppl.F:102
ai_overlap_ppl::ecploc_integral
subroutine, public ecploc_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, hab2, hab2_work, deltar, iatom, jatom, katom)
Calculation of three-center potential integrals <a|V(r)|b> over Cartesian Gaussian functions for the ...
Definition ai_overlap_ppl.F:236
ai_overlap_ppl::ppl_integral_ri
subroutine, public ppl_integral_ri(la_max_set, la_min_set, npgfa, rpgfa, zeta, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, rac, dac, va, dva)
Calculation of two-center overlap integrals <a|c> over Cartesian Gaussian functions for the local par...
Definition ai_overlap_ppl.F:349
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
gamma::fgamma_0
subroutine, public fgamma_0(nmax, t, f)
Calculation of the incomplete Gamma function F(t) for multicenter integrals over Gaussian functions....
Definition gamma.F:154
gfun
Calculation of the G function G_n(t) for 1/R^2 operators.
Definition gfun.F:29
gfun::gfun_values
subroutine, public gfun_values(nmax, t, g)
...
Definition gfun.F:53
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::binomial
elemental real(kind=dp) function, public binomial(n, k)
The binomial coefficient n over k for 0 <= k <= n is calculated, otherwise zero is returned.
Definition mathlib.F:214