(git:ec11232)
Loading...
Searching...
No Matches
ai_overlap.F
Go to the documentation of this file.
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 the overlap integrals over Cartesian Gaussian-type
10!> functions.
11!> \par Literature
12!> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
13!> \par History
14!> - Derivatives added (02.05.2002,MK)
15!> - New OS routine with simpler logic (11.07.2014, JGH)
16!> \author Matthias Krack (08.10.1999)
17! **************************************************************************************************
18MODULE ai_overlap
19 USE ai_os_rr, ONLY: os_rr_ovlp
20 USE kinds, ONLY: dp
21 USE mathconstants, ONLY: pi,&
22 twopi,&
23 z_one
24 USE orbital_pointers, ONLY: coset,&
25 nco,&
26 ncoset,&
27 nso
28 USE orbital_transformation_matrices, ONLY: orbtramat
29#include "../base/base_uses.f90"
30
31 IMPLICIT NONE
32
33 PRIVATE
34
35 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap'
36
37! *** Public subroutines ***
38 PUBLIC :: overlap, overlap_ab, overlap_aab, overlap_ab_s, overlap_ab_sp, &
39 overlap_abb
40
41CONTAINS
42
43! **************************************************************************************************
44!> \brief Purpose: Calculation of the two-center overlap integrals [a|b] over
45!> Cartesian Gaussian-type functions.
46!> \param la_max_set Max L on center A
47!> \param la_min_set Min L on center A
48!> \param npgfa Number of primitives on center A
49!> \param rpgfa Range of functions on A, used for screening
50!> \param zeta Exponents on center A
51!> \param lb_max_set Max L on center B
52!> \param lb_min_set Min L on center B
53!> \param npgfb Number of primitives on center B
54!> \param rpgfb Range of functions on B, used for screening
55!> \param zetb Exponents on center B
56!> \param rab Distance vector A-B
57!> \param dab Distance A-B
58!> \param sab Final Integrals, basic and derivatives
59!> \param da_max_set Some additional derivative information
60!> \param return_derivatives Return integral derivatives
61!> \param s Work space
62!> \param lds Leading dimension of s
63!> \param sdab Return additional derivative integrals
64!> \param pab Density matrix block, used to calculate forces
65!> \param force_a Force vector [da/dR|b]
66!> \date 19.09.2000
67!> \author MK
68!> \version 1.0
69! **************************************************************************************************
70 SUBROUTINE overlap(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
71 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
72 rab, dab, sab, da_max_set, return_derivatives, s, lds, &
73 sdab, pab, force_a)
74 INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
75 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
76 INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
77 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
78 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
79 REAL(kind=dp), INTENT(IN) :: dab
80 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
81 INTEGER, INTENT(IN) :: da_max_set
82 LOGICAL, INTENT(IN) :: return_derivatives
83 INTEGER, INTENT(IN) :: lds
84 REAL(kind=dp), DIMENSION(lds, lds, *), &
85 INTENT(INOUT) :: s
86 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
87 OPTIONAL :: sdab
88 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
89 OPTIONAL :: pab
90 REAL(kind=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a
91
92 INTEGER :: ax, ay, az, bx, by, bz, cda, cdax, cday, cdaz, coa, coamx, coamy, coamz, coapx, &
93 coapy, coapz, cob, cobm2x, cobm2y, cobm2z, cobmx, cobmy, cobmz, da, da_max, dax, day, &
94 daz, i, ipgf, j, jk, jpgf, jstart, k, la, la_max, la_min, la_start, lb, lb_max, lb_min, &
95 lb_start, na, nb, nda
96 LOGICAL :: calculate_force_a
97 REAL(kind=dp) :: f0, f1, f2, f3, f4, fax, fay, faz, ftz, &
98 zetp
99 REAL(kind=dp), DIMENSION(3) :: rap, rbp
100
101 IF (PRESENT(pab) .AND. PRESENT(force_a)) THEN
102 calculate_force_a = .true.
103 force_a(:) = 0.0_dp
104 ELSE
105 calculate_force_a = .false.
106 END IF
107
108 IF (PRESENT(sdab) .OR. calculate_force_a) THEN
109 IF (da_max_set == 0) THEN
110 da_max = 1
111 la_max = la_max_set + 1
112 la_min = max(0, la_min_set - 1)
113 ELSE
114 da_max = da_max_set
115 la_max = la_max_set + da_max_set + 1
116 la_min = max(0, la_min_set - da_max_set - 1)
117 END IF
118 ELSE
119 da_max = da_max_set
120 la_max = la_max_set + da_max_set
121 la_min = max(0, la_min_set - da_max_set)
122 END IF
123
124 lb_max = lb_max_set
125 lb_min = lb_min_set
126
127! *** Loop over all pairs of primitive Gaussian-type functions ***
128
129 na = 0
130 nda = 0
131
132 DO ipgf = 1, npgfa
133
134 nb = 0
135
136 DO jpgf = 1, npgfb
137
138! *** Screening ***
139
140 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
141 DO j = nb + 1, nb + ncoset(lb_max_set)
142 DO i = na + 1, na + ncoset(la_max_set)
143 sab(i, j) = 0.0_dp
144 END DO
145 END DO
146 IF (return_derivatives) THEN
147 DO k = 2, ncoset(da_max_set)
148 jstart = (k - 1)*SIZE(sab, 1)
149 DO j = jstart + nb + 1, jstart + nb + ncoset(lb_max_set)
150 DO i = na + 1, na + ncoset(la_max_set)
151 sab(i, j) = 0.0_dp
152 END DO
153 END DO
154 END DO
155 END IF
156 nb = nb + ncoset(lb_max_set)
157 cycle
158 END IF
159
160! *** Calculate some prefactors ***
161
162 zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
163
164 f0 = sqrt((pi*zetp)**3)
165 f1 = zetb(jpgf)*zetp
166 f2 = 0.5_dp*zetp
167
168! *** Calculate the basic two-center overlap integral [s|s] ***
169
170 s(1, 1, 1) = f0*exp(-zeta(ipgf)*f1*dab*dab)
171
172! *** Recurrence steps: [s|s] -> [a|b] ***
173
174 IF (la_max > 0) THEN
175
176! *** Vertical recurrence steps: [s|s] -> [a|s] ***
177
178 rap(:) = f1*rab(:)
179
180! *** [p|s] = (Pi - Ai)*[s|s] (i = x,y,z) ***
181
182 s(2, 1, 1) = rap(1)*s(1, 1, 1) ! [px|s]
183 s(3, 1, 1) = rap(2)*s(1, 1, 1) ! [py|s]
184 s(4, 1, 1) = rap(3)*s(1, 1, 1) ! [pz|s]
185
186 IF (la_max > 1) THEN
187
188! *** [d|s] ***
189
190 f3 = f2*s(1, 1, 1)
191
192 s(5, 1, 1) = rap(1)*s(2, 1, 1) + f3 ! [dx2|s]
193 s(6, 1, 1) = rap(1)*s(3, 1, 1) ! [dxy|s]
194 s(7, 1, 1) = rap(1)*s(4, 1, 1) ! [dxz|s]
195 s(8, 1, 1) = rap(2)*s(3, 1, 1) + f3 ! [dy2|s]
196 s(9, 1, 1) = rap(2)*s(4, 1, 1) ! [dyz|s]
197 s(10, 1, 1) = rap(3)*s(4, 1, 1) + f3 ! [dz2|s]
198
199 IF (la_max > 2) THEN
200
201! *** [f|s] ***
202
203 f3 = 2.0_dp*f2
204
205 s(11, 1, 1) = rap(1)*s(5, 1, 1) + f3*s(2, 1, 1) ! [fx3 |s]
206 s(12, 1, 1) = rap(1)*s(6, 1, 1) + f2*s(3, 1, 1) ! [fx2y|s]
207 s(13, 1, 1) = rap(1)*s(7, 1, 1) + f2*s(4, 1, 1) ! [fx2z|s]
208 s(14, 1, 1) = rap(1)*s(8, 1, 1) ! [fxy2|s]
209 s(15, 1, 1) = rap(1)*s(9, 1, 1) ! [fxyz|s]
210 s(16, 1, 1) = rap(1)*s(10, 1, 1) ! [fxz2|s]
211 s(17, 1, 1) = rap(2)*s(8, 1, 1) + f3*s(3, 1, 1) ! [fy3 |s]
212 s(18, 1, 1) = rap(2)*s(9, 1, 1) + f2*s(4, 1, 1) ! [fy2z|s]
213 s(19, 1, 1) = rap(2)*s(10, 1, 1) ! [fyz2|s]
214 s(20, 1, 1) = rap(3)*s(10, 1, 1) + f3*s(4, 1, 1) ! [fz3 |s]
215
216 IF (la_max > 3) THEN
217
218! *** [g|s] ***
219
220 f4 = 3.0_dp*f2
221
222 s(21, 1, 1) = rap(1)*s(11, 1, 1) + f4*s(5, 1, 1) ! [gx4 |s]
223 s(22, 1, 1) = rap(1)*s(12, 1, 1) + f3*s(6, 1, 1) ! [gx3y |s]
224 s(23, 1, 1) = rap(1)*s(13, 1, 1) + f3*s(7, 1, 1) ! [gx3z |s]
225 s(24, 1, 1) = rap(1)*s(14, 1, 1) + f2*s(8, 1, 1) ! [gx2y2|s]
226 s(25, 1, 1) = rap(1)*s(15, 1, 1) + f2*s(9, 1, 1) ! [gx2yz|s]
227 s(26, 1, 1) = rap(1)*s(16, 1, 1) + f2*s(10, 1, 1) ! [gx2z2|s]
228 s(27, 1, 1) = rap(1)*s(17, 1, 1) ! [gxy3 |s]
229 s(28, 1, 1) = rap(1)*s(18, 1, 1) ! [gxy2z|s]
230 s(29, 1, 1) = rap(1)*s(19, 1, 1) ! [gxyz2|s]
231 s(30, 1, 1) = rap(1)*s(20, 1, 1) ! [gxz3 |s]
232 s(31, 1, 1) = rap(2)*s(17, 1, 1) + f4*s(8, 1, 1) ! [gy4 |s]
233 s(32, 1, 1) = rap(2)*s(18, 1, 1) + f3*s(9, 1, 1) ! [gy3z |s]
234 s(33, 1, 1) = rap(2)*s(19, 1, 1) + f2*s(10, 1, 1) ! [gy2z2|s]
235 s(34, 1, 1) = rap(2)*s(20, 1, 1) ! [gyz3 |s]
236 s(35, 1, 1) = rap(3)*s(20, 1, 1) + f4*s(10, 1, 1) ! [gz4 |s]
237
238! *** [a|s] = (Pi - Ai)*[a-1i|s] + f2*Ni(a-1i)*[a-2i|s] ***
239
240 DO la = 5, la_max
241
242! *** Increase the angular momentum component z of a ***
243
244 s(coset(0, 0, la), 1, 1) = &
245 rap(3)*s(coset(0, 0, la - 1), 1, 1) + &
246 f2*real(la - 1, dp)*s(coset(0, 0, la - 2), 1, 1)
247
248! *** Increase the angular momentum component y of a ***
249
250 az = la - 1
251 s(coset(0, 1, az), 1, 1) = rap(2)*s(coset(0, 0, az), 1, 1)
252 DO ay = 2, la
253 az = la - ay
254 s(coset(0, ay, az), 1, 1) = &
255 rap(2)*s(coset(0, ay - 1, az), 1, 1) + &
256 f2*real(ay - 1, dp)*s(coset(0, ay - 2, az), 1, 1)
257 END DO
258
259! *** Increase the angular momentum component x of a ***
260
261 DO ay = 0, la - 1
262 az = la - 1 - ay
263 s(coset(1, ay, az), 1, 1) = rap(1)*s(coset(0, ay, az), 1, 1)
264 END DO
265 DO ax = 2, la
266 f3 = f2*real(ax - 1, dp)
267 DO ay = 0, la - ax
268 az = la - ax - ay
269 s(coset(ax, ay, az), 1, 1) = &
270 rap(1)*s(coset(ax - 1, ay, az), 1, 1) + &
271 f3*s(coset(ax - 2, ay, az), 1, 1)
272 END DO
273 END DO
274
275 END DO
276
277 END IF
278
279 END IF
280
281 END IF
282
283! *** Recurrence steps: [a|s] -> [a|b] ***
284
285 IF (lb_max > 0) THEN
286
287 DO j = 2, ncoset(lb_max)
288 DO i = 1, ncoset(la_min)
289 s(i, j, 1) = 0.0_dp
290 END DO
291 END DO
292
293! *** Horizontal recurrence steps ***
294
295 rbp(:) = rap(:) - rab(:)
296
297! *** [a|p] = [a+1i|s] - (Bi - Ai)*[a|s] ***
298
299 IF (lb_max == 1) THEN
300 la_start = la_min
301 ELSE
302 la_start = max(0, la_min - 1)
303 END IF
304
305 DO la = la_start, la_max - 1
306 DO ax = 0, la
307 DO ay = 0, la - ax
308 az = la - ax - ay
309 coa = coset(ax, ay, az)
310 coapx = coset(ax + 1, ay, az)
311 coapy = coset(ax, ay + 1, az)
312 coapz = coset(ax, ay, az + 1)
313 s(coa, 2, 1) = s(coapx, 1, 1) - rab(1)*s(coa, 1, 1)
314 s(coa, 3, 1) = s(coapy, 1, 1) - rab(2)*s(coa, 1, 1)
315 s(coa, 4, 1) = s(coapz, 1, 1) - rab(3)*s(coa, 1, 1)
316 END DO
317 END DO
318 END DO
319
320! *** Vertical recurrence step ***
321
322! *** [a|p] = (Pi - Bi)*[a|s] + f2*Ni(a)*[a-1i|s] ***
323
324 DO ax = 0, la_max
325 fax = f2*real(ax, dp)
326 DO ay = 0, la_max - ax
327 fay = f2*real(ay, dp)
328 az = la_max - ax - ay
329 faz = f2*real(az, dp)
330 coa = coset(ax, ay, az)
331 coamx = coset(ax - 1, ay, az)
332 coamy = coset(ax, ay - 1, az)
333 coamz = coset(ax, ay, az - 1)
334 s(coa, 2, 1) = rbp(1)*s(coa, 1, 1) + fax*s(coamx, 1, 1)
335 s(coa, 3, 1) = rbp(2)*s(coa, 1, 1) + fay*s(coamy, 1, 1)
336 s(coa, 4, 1) = rbp(3)*s(coa, 1, 1) + faz*s(coamz, 1, 1)
337 END DO
338 END DO
339
340! *** Recurrence steps: [a|p] -> [a|b] ***
341
342 DO lb = 2, lb_max
343
344! *** Horizontal recurrence steps ***
345
346! *** [a|b] = [a+1i|b-1i] - (Bi - Ai)*[a|b-1i] ***
347
348 IF (lb == lb_max) THEN
349 la_start = la_min
350 ELSE
351 la_start = max(0, la_min - 1)
352 END IF
353
354 DO la = la_start, la_max - 1
355 DO ax = 0, la
356 DO ay = 0, la - ax
357 az = la - ax - ay
358 coa = coset(ax, ay, az)
359 coapx = coset(ax + 1, ay, az)
360 coapy = coset(ax, ay + 1, az)
361 coapz = coset(ax, ay, az + 1)
362
363! *** Shift of angular momentum component z from a to b ***
364
365 cob = coset(0, 0, lb)
366 cobmz = coset(0, 0, lb - 1)
367 s(coa, cob, 1) = s(coapz, cobmz, 1) - rab(3)*s(coa, cobmz, 1)
368
369! *** Shift of angular momentum component y from a to b ***
370
371 DO by = 1, lb
372 bz = lb - by
373 cob = coset(0, by, bz)
374 cobmy = coset(0, by - 1, bz)
375 s(coa, cob, 1) = s(coapy, cobmy, 1) - rab(2)*s(coa, cobmy, 1)
376 END DO
377
378! *** Shift of angular momentum component x from a to b ***
379
380 DO bx = 1, lb
381 DO by = 0, lb - bx
382 bz = lb - bx - by
383 cob = coset(bx, by, bz)
384 cobmx = coset(bx - 1, by, bz)
385 s(coa, cob, 1) = s(coapx, cobmx, 1) - rab(1)*s(coa, cobmx, 1)
386 END DO
387 END DO
388
389 END DO
390 END DO
391 END DO
392
393! *** Vertical recurrence step ***
394
395! *** [a|b] = (Pi - Bi)*[a|b-1i] + f2*Ni(a)*[a-1i|b-1i] + ***
396! *** f2*Ni(b-1i)*[a|b-2i] ***
397
398 DO ax = 0, la_max
399 fax = f2*real(ax, dp)
400 DO ay = 0, la_max - ax
401 fay = f2*real(ay, dp)
402 az = la_max - ax - ay
403 faz = f2*real(az, dp)
404 coa = coset(ax, ay, az)
405 coamx = coset(ax - 1, ay, az)
406 coamy = coset(ax, ay - 1, az)
407 coamz = coset(ax, ay, az - 1)
408
409! *** Increase the angular momentum component z of b ***
410
411 f3 = f2*real(lb - 1, dp)
412 cob = coset(0, 0, lb)
413 cobmz = coset(0, 0, lb - 1)
414 cobm2z = coset(0, 0, lb - 2)
415 s(coa, cob, 1) = rbp(3)*s(coa, cobmz, 1) + &
416 faz*s(coamz, cobmz, 1) + &
417 f3*s(coa, cobm2z, 1)
418
419! *** Increase the angular momentum component y of b ***
420
421 bz = lb - 1
422 cob = coset(0, 1, bz)
423 cobmy = coset(0, 0, bz)
424 s(coa, cob, 1) = rbp(2)*s(coa, cobmy, 1) + &
425 fay*s(coamy, cobmy, 1)
426 DO by = 2, lb
427 bz = lb - by
428 f3 = f2*real(by - 1, dp)
429 cob = coset(0, by, bz)
430 cobmy = coset(0, by - 1, bz)
431 cobm2y = coset(0, by - 2, bz)
432 s(coa, cob, 1) = rbp(2)*s(coa, cobmy, 1) + &
433 fay*s(coamy, cobmy, 1) + &
434 f3*s(coa, cobm2y, 1)
435 END DO
436
437! *** Increase the angular momentum component x of b ***
438
439 DO by = 0, lb - 1
440 bz = lb - 1 - by
441 cob = coset(1, by, bz)
442 cobmx = coset(0, by, bz)
443 s(coa, cob, 1) = rbp(1)*s(coa, cobmx, 1) + &
444 fax*s(coamx, cobmx, 1)
445 END DO
446 DO bx = 2, lb
447 f3 = f2*real(bx - 1, dp)
448 DO by = 0, lb - bx
449 bz = lb - bx - by
450 cob = coset(bx, by, bz)
451 cobmx = coset(bx - 1, by, bz)
452 cobm2x = coset(bx - 2, by, bz)
453 s(coa, cob, 1) = rbp(1)*s(coa, cobmx, 1) + &
454 fax*s(coamx, cobmx, 1) + &
455 f3*s(coa, cobm2x, 1)
456 END DO
457 END DO
458
459 END DO
460 END DO
461
462 END DO
463
464 END IF
465
466 ELSE
467
468 IF (lb_max > 0) THEN
469
470! *** Vertical recurrence steps: [s|s] -> [s|b] ***
471
472 rbp(:) = (f1 - 1.0_dp)*rab(:)
473
474! *** [s|p] = (Pi - Bi)*[s|s] ***
475
476 s(1, 2, 1) = rbp(1)*s(1, 1, 1) ! [s|px]
477 s(1, 3, 1) = rbp(2)*s(1, 1, 1) ! [s|py]
478 s(1, 4, 1) = rbp(3)*s(1, 1, 1) ! [s|pz]
479
480 IF (lb_max > 1) THEN
481
482! *** [s|d] ***
483
484 f3 = f2*s(1, 1, 1)
485
486 s(1, 5, 1) = rbp(1)*s(1, 2, 1) + f3 ! [s|dx2]
487 s(1, 6, 1) = rbp(1)*s(1, 3, 1) ! [s|dxy]
488 s(1, 7, 1) = rbp(1)*s(1, 4, 1) ! [s|dxz]
489 s(1, 8, 1) = rbp(2)*s(1, 3, 1) + f3 ! [s|dy2]
490 s(1, 9, 1) = rbp(2)*s(1, 4, 1) ! [s|dyz]
491 s(1, 10, 1) = rbp(3)*s(1, 4, 1) + f3 ! [s|dz2]
492
493! *** [s|b] = (Pi - Bi)*[s|b-1i] + f2*Ni(b-1i)*[s|b-2i] ***
494
495 DO lb = 3, lb_max
496
497! *** Increase the angular momentum component z of b ***
498
499 s(1, coset(0, 0, lb), 1) = &
500 rbp(3)*s(1, coset(0, 0, lb - 1), 1) + &
501 f2*real(lb - 1, dp)*s(1, coset(0, 0, lb - 2), 1)
502
503! *** Increase the angular momentum component y of b ***
504
505 bz = lb - 1
506 s(1, coset(0, 1, bz), 1) = rbp(2)*s(1, coset(0, 0, bz), 1)
507 DO by = 2, lb
508 bz = lb - by
509 s(1, coset(0, by, bz), 1) = &
510 rbp(2)*s(1, coset(0, by - 1, bz), 1) + &
511 f2*real(by - 1, dp)*s(1, coset(0, by - 2, bz), 1)
512 END DO
513
514! *** Increase the angular momentum component x of b ***
515
516 DO by = 0, lb - 1
517 bz = lb - 1 - by
518 s(1, coset(1, by, bz), 1) = rbp(1)*s(1, coset(0, by, bz), 1)
519 END DO
520 DO bx = 2, lb
521 f3 = f2*real(bx - 1, dp)
522 DO by = 0, lb - bx
523 bz = lb - bx - by
524 s(1, coset(bx, by, bz), 1) = &
525 rbp(1)*s(1, coset(bx - 1, by, bz), 1) + &
526 f3*s(1, coset(bx - 2, by, bz), 1)
527 END DO
528 END DO
529
530 END DO
531
532 END IF
533
534 END IF
535
536 END IF
537
538! *** Store the primitive overlap integrals ***
539
540 DO j = 1, ncoset(lb_max_set)
541 DO i = 1, ncoset(la_max_set)
542 sab(na + i, nb + j) = s(i, j, 1)
543 END DO
544 END DO
545
546! *** Calculate the requested derivatives with respect ***
547! *** to the nuclear coordinates of the atomic center a ***
548
549 IF (PRESENT(sdab) .OR. return_derivatives) THEN
550 la_start = 0
551 lb_start = 0
552 ELSE
553 la_start = la_min_set
554 lb_start = lb_min_set
555 END IF
556
557 DO da = 0, da_max - 1
558 ftz = 2.0_dp*zeta(ipgf)
559 DO dax = 0, da
560 DO day = 0, da - dax
561 daz = da - dax - day
562 cda = coset(dax, day, daz)
563 cdax = coset(dax + 1, day, daz)
564 cday = coset(dax, day + 1, daz)
565 cdaz = coset(dax, day, daz + 1)
566
567! *** [da/dAi|b] = 2*zeta*[a+1i|b] - Ni(a)[a-1i|b] ***
568
569 DO la = la_start, la_max - da - 1
570 DO ax = 0, la
571 fax = real(ax, dp)
572 DO ay = 0, la - ax
573 fay = real(ay, dp)
574 az = la - ax - ay
575 faz = real(az, dp)
576 coa = coset(ax, ay, az)
577 coamx = coset(ax - 1, ay, az)
578 coamy = coset(ax, ay - 1, az)
579 coamz = coset(ax, ay, az - 1)
580 coapx = coset(ax + 1, ay, az)
581 coapy = coset(ax, ay + 1, az)
582 coapz = coset(ax, ay, az + 1)
583 DO lb = lb_start, lb_max_set
584 DO bx = 0, lb
585 DO by = 0, lb - bx
586 bz = lb - bx - by
587 cob = coset(bx, by, bz)
588 s(coa, cob, cdax) = ftz*s(coapx, cob, cda) - &
589 fax*s(coamx, cob, cda)
590 s(coa, cob, cday) = ftz*s(coapy, cob, cda) - &
591 fay*s(coamy, cob, cda)
592 s(coa, cob, cdaz) = ftz*s(coapz, cob, cda) - &
593 faz*s(coamz, cob, cda)
594 END DO
595 END DO
596 END DO
597 END DO
598 END DO
599 END DO
600
601 END DO
602 END DO
603 END DO
604
605! *** Return all the calculated derivatives of the ***
606! *** primitive overlap integrals, if requested ***
607
608 IF (return_derivatives) THEN
609 DO k = 2, ncoset(da_max_set)
610 jstart = (k - 1)*SIZE(sab, 1)
611 DO j = 1, ncoset(lb_max_set)
612 jk = jstart + j
613 DO i = 1, ncoset(la_max_set)
614 sab(na + i, nb + jk) = s(i, j, k)
615 END DO
616 END DO
617 END DO
618 END IF
619
620! *** Calculate the force contribution for the atomic center a ***
621
622 IF (calculate_force_a) THEN
623 DO k = 1, 3
624 DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
625 DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
626 force_a(k) = force_a(k) + pab(na + i, nb + j)*s(i, j, k + 1)
627 END DO
628 END DO
629 END DO
630 END IF
631
632! *** Store the first derivatives of the primitive overlap integrals ***
633! *** which are used as auxiliary integrals for the calculation of ***
634! *** the kinetic energy integrals if requested ***
635
636 IF (PRESENT(sdab)) THEN
637 sdab(nda + 1, nb + 1, 1) = s(1, 1, 1)
638 DO k = 2, 4
639 DO j = 1, ncoset(lb_max_set)
640 DO i = 1, ncoset(la_max - 1)
641 sdab(nda + i, nb + j, k) = s(i, j, k)
642 END DO
643 END DO
644 END DO
645 END IF
646
647 nb = nb + ncoset(lb_max_set)
648
649 END DO
650
651 na = na + ncoset(la_max_set)
652 nda = nda + ncoset(la_max - 1)
653
654 END DO
655
656 END SUBROUTINE overlap
657
658! **************************************************************************************************
659!> \brief Calculation of the two-center overlap integrals [a|b] over
660!> Cartesian Gaussian-type functions. First and second derivatives
661!> \param la_max Max L on center A
662!> \param la_min Min L on center A
663!> \param npgfa Number of primitives on center A
664!> \param rpgfa Range of functions on A, used for screening
665!> \param zeta Exponents on center A
666!> \param lb_max Max L on center B
667!> \param lb_min Min L on center B
668!> \param npgfb Number of primitives on center B
669!> \param rpgfb Range of functions on B, used for screening
670!> \param zetb Exponents on center B
671!> \param rab Distance vector A-B
672!> \param sab Final overlap integrals
673!> \param dab First derivative overlap integrals
674!> \param ddab Second derivative overlap integrals
675!> \date 01.07.2014
676!> \author JGH
677! **************************************************************************************************
678 SUBROUTINE overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, &
679 lb_max, lb_min, npgfb, rpgfb, zetb, &
680 rab, sab, dab, ddab)
681 INTEGER, INTENT(IN) :: la_max, la_min, npgfa
682 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
683 INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
684 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
685 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
686 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT), &
687 OPTIONAL :: sab
688 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
689 OPTIONAL :: dab, ddab
690
691 INTEGER :: ax, ay, az, bx, by, bz, coa, cob, ia, &
692 ib, ipgf, jpgf, la, lb, ldrr, lma, &
693 lmb, ma, mb, na, nb, ofa, ofb
694 REAL(kind=dp) :: a, ambm, ambp, apbm, apbp, b, dumx, &
695 dumy, dumz, f0, rab2, tab, xhi, zet
696 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
697 REAL(kind=dp), DIMENSION(3) :: rap, rbp
698
699 ! Distance of the centers a and b
700
701 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
702 tab = sqrt(rab2)
703
704 ! Maximum l for auxiliary integrals
705 IF (PRESENT(sab)) THEN
706 lma = la_max
707 lmb = lb_max
708 END IF
709 IF (PRESENT(dab)) THEN
710 lma = la_max + 1
711 lmb = lb_max
712 END IF
713 IF (PRESENT(ddab)) THEN
714 lma = la_max + 1
715 lmb = lb_max + 1
716 END IF
717 ldrr = max(lma, lmb) + 1
718
719 ! Allocate space for auxiliary integrals
720 ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
721
722 ! Number of integrals, check size of arrays
723 ofa = ncoset(la_min - 1)
724 ofb = ncoset(lb_min - 1)
725 na = ncoset(la_max) - ofa
726 nb = ncoset(lb_max) - ofb
727 IF (PRESENT(sab)) THEN
728 cpassert((SIZE(sab, 1) >= na*npgfa))
729 cpassert((SIZE(sab, 2) >= nb*npgfb))
730 END IF
731 IF (PRESENT(dab)) THEN
732 cpassert((SIZE(dab, 1) >= na*npgfa))
733 cpassert((SIZE(dab, 2) >= nb*npgfb))
734 cpassert((SIZE(dab, 3) >= 3))
735 END IF
736 IF (PRESENT(ddab)) THEN
737 cpassert((SIZE(ddab, 1) >= na*npgfa))
738 cpassert((SIZE(ddab, 2) >= nb*npgfb))
739 cpassert((SIZE(ddab, 3) >= 6))
740 END IF
741
742 ! Loops over all pairs of primitive Gaussian-type functions
743 ma = 0
744 DO ipgf = 1, npgfa
745 mb = 0
746 DO jpgf = 1, npgfb
747 ! Distance Screening
748 IF (rpgfa(ipgf) + rpgfb(jpgf) < tab) THEN
749 IF (PRESENT(sab)) sab(ma + 1:ma + na, mb + 1:mb + nb) = 0.0_dp
750 IF (PRESENT(dab)) dab(ma + 1:ma + na, mb + 1:mb + nb, 1:3) = 0.0_dp
751 IF (PRESENT(ddab)) ddab(ma + 1:ma + na, mb + 1:mb + nb, 1:6) = 0.0_dp
752 mb = mb + nb
753 cycle
754 END IF
755
756 ! Calculate some prefactors
757 a = zeta(ipgf)
758 b = zetb(jpgf)
759 zet = a + b
760 xhi = a*b/zet
761 rap = b*rab/zet
762 rbp = -a*rab/zet
763
764 ! [s|s] integral
765 f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
766
767 ! Calculate the recurrence relation
768 CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
769
770 DO lb = lb_min, lb_max
771 DO bx = 0, lb
772 DO by = 0, lb - bx
773 bz = lb - bx - by
774 cob = coset(bx, by, bz) - ofb
775 ib = mb + cob
776 DO la = la_min, la_max
777 DO ax = 0, la
778 DO ay = 0, la - ax
779 az = la - ax - ay
780 coa = coset(ax, ay, az) - ofa
781 ia = ma + coa
782 ! integrals
783 IF (PRESENT(sab)) THEN
784 sab(ia, ib) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az, bz, 3)
785 END IF
786 ! first derivatives
787 IF (PRESENT(dab)) THEN
788 ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
789 ! dx
790 dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
791 IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
792 dab(ia, ib, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
793 ! dy
794 dumy = 2.0_dp*a*rr(ay + 1, by, 2)
795 IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
796 dab(ia, ib, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
797 ! dz
798 dumz = 2.0_dp*a*rr(az + 1, bz, 3)
799 IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
800 dab(ia, ib, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
801 END IF
802 ! 2nd derivatives
803 IF (PRESENT(ddab)) THEN
804 ! (dda|b) = -4*a*b*(a+1|b+1) + 2*a*N(b)*(a+1|b-1)
805 ! + 2*b*N(a)*(a-1|b+1) - N(a)*N(b)*(a-1|b-1)
806 ! dx dx
807 apbp = f0*rr(ax + 1, bx + 1, 1)*rr(ay, by, 2)*rr(az, bz, 3)
808 IF (bx > 0) THEN
809 apbm = f0*rr(ax + 1, bx - 1, 1)*rr(ay, by, 2)*rr(az, bz, 3)
810 ELSE
811 apbm = 0.0_dp
812 END IF
813 IF (ax > 0) THEN
814 ambp = f0*rr(ax - 1, bx + 1, 1)*rr(ay, by, 2)*rr(az, bz, 3)
815 ELSE
816 ambp = 0.0_dp
817 END IF
818 IF (ax > 0 .AND. bx > 0) THEN
819 ambm = f0*rr(ax - 1, bx - 1, 1)*rr(ay, by, 2)*rr(az, bz, 3)
820 ELSE
821 ambm = 0.0_dp
822 END IF
823 ddab(ia, ib, 1) = -4.0_dp*a*b*apbm + 2.0_dp*a*real(bx, dp)*apbm &
824 + 2.0_dp*b*real(ax, dp)*ambp - real(ax, dp)*real(bx, dp)*ambm
825 ! dx dy
826 apbp = f0*rr(ax + 1, bx, 1)*rr(ay, by + 1, 2)*rr(az, bz, 3)
827 IF (by > 0) THEN
828 apbm = f0*rr(ax + 1, bx, 1)*rr(ay, by - 1, 2)*rr(az, bz, 3)
829 ELSE
830 apbm = 0.0_dp
831 END IF
832 IF (ax > 0) THEN
833 ambp = f0*rr(ax - 1, bx, 1)*rr(ay, by + 1, 2)*rr(az, bz, 3)
834 ELSE
835 ambp = 0.0_dp
836 END IF
837 IF (ax > 0 .AND. by > 0) THEN
838 ambm = f0*rr(ax - 1, bx, 1)*rr(ay, by - 1, 2)*rr(az, bz, 3)
839 ELSE
840 ambm = 0.0_dp
841 END IF
842 ddab(ia, ib, 2) = -4.0_dp*a*b*apbm + 2.0_dp*a*real(by, dp)*apbm &
843 + 2.0_dp*b*real(ax, dp)*ambp - real(ax, dp)*real(by, dp)*ambm
844 ! dx dz
845 apbp = f0*rr(ax + 1, bx, 1)*rr(ay, by, 2)*rr(az, bz + 1, 3)
846 IF (bz > 0) THEN
847 apbm = f0*rr(ax + 1, bx, 1)*rr(ay, by, 2)*rr(az, bz - 1, 3)
848 ELSE
849 apbm = 0.0_dp
850 END IF
851 IF (ax > 0) THEN
852 ambp = f0*rr(ax - 1, bx, 1)*rr(ay, by, 2)*rr(az, bz + 1, 3)
853 ELSE
854 ambp = 0.0_dp
855 END IF
856 IF (ax > 0 .AND. bz > 0) THEN
857 ambm = f0*rr(ax - 1, bx, 1)*rr(ay, by, 2)*rr(az, bz - 1, 3)
858 ELSE
859 ambm = 0.0_dp
860 END IF
861 ddab(ia, ib, 3) = -4.0_dp*a*b*apbm + 2.0_dp*a*real(bz, dp)*apbm &
862 + 2.0_dp*b*real(ax, dp)*ambp - real(ax, dp)*real(bz, dp)*ambm
863 ! dy dy
864 apbp = f0*rr(ax, bx, 1)*rr(ay + 1, by + 1, 2)*rr(az, bz, 3)
865 IF (by > 0) THEN
866 apbm = f0*rr(ax, bx, 1)*rr(ay + 1, by - 1, 2)*rr(az, bz, 3)
867 ELSE
868 apbm = 0.0_dp
869 END IF
870 IF (ay > 0) THEN
871 ambp = f0*rr(ax, bx, 1)*rr(ay - 1, by + 1, 2)*rr(az, bz, 3)
872 ELSE
873 ambp = 0.0_dp
874 END IF
875 IF (ay > 0 .AND. by > 0) THEN
876 ambm = f0*rr(ax, bx, 1)*rr(ay - 1, by - 1, 2)*rr(az, bz, 3)
877 ELSE
878 ambm = 0.0_dp
879 END IF
880 ddab(ia, ib, 4) = -4.0_dp*a*b*apbm + 2.0_dp*a*real(by, dp)*apbm &
881 + 2.0_dp*b*real(ay, dp)*ambp - real(ay, dp)*real(by, dp)*ambm
882 ! dy dz
883 apbp = f0*rr(ax, bx, 1)*rr(ay + 1, by, 2)*rr(az, bz + 1, 3)
884 IF (bz > 0) THEN
885 apbm = f0*rr(ax, bx, 1)*rr(ay + 1, by, 2)*rr(az, bz - 1, 3)
886 ELSE
887 apbm = 0.0_dp
888 END IF
889 IF (ay > 0) THEN
890 ambp = f0*rr(ax, bx, 1)*rr(ay - 1, by, 2)*rr(az, bz + 1, 3)
891 ELSE
892 ambp = 0.0_dp
893 END IF
894 IF (ay > 0 .AND. bz > 0) THEN
895 ambm = f0*rr(ax, bx, 1)*rr(ay - 1, by, 2)*rr(az, bz - 1, 3)
896 ELSE
897 ambm = 0.0_dp
898 END IF
899 ddab(ia, ib, 5) = -4.0_dp*a*b*apbm + 2.0_dp*a*real(bz, dp)*apbm &
900 + 2.0_dp*b*real(ay, dp)*ambp - real(ay, dp)*real(bz, dp)*ambm
901 ! dz dz
902 apbp = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az + 1, bz + 1, 3)
903 IF (bz > 0) THEN
904 apbm = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az + 1, bz - 1, 3)
905 ELSE
906 apbm = 0.0_dp
907 END IF
908 IF (az > 0) THEN
909 ambp = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az - 1, bz + 1, 3)
910 ELSE
911 ambp = 0.0_dp
912 END IF
913 IF (az > 0 .AND. bz > 0) THEN
914 ambm = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az - 1, bz - 1, 3)
915 ELSE
916 ambm = 0.0_dp
917 END IF
918 ddab(ia, ib, 6) = -4.0_dp*a*b*apbm + 2.0_dp*a*real(bz, dp)*apbm &
919 + 2.0_dp*b*real(az, dp)*ambp - real(az, dp)*real(bz, dp)*ambm
920 END IF
921 !
922 END DO
923 END DO
924 END DO !la
925 END DO
926 END DO
927 END DO !lb
928
929 mb = mb + nb
930 END DO
931 ma = ma + na
932 END DO
933
934 DEALLOCATE (rr)
935
936 END SUBROUTINE overlap_ab
937
938! **************************************************************************************************
939!> \brief Calculation of the two-center overlap integrals [aa|b] over
940!> Cartesian Gaussian-type functions.
941!> \param la1_max Max L on center A (basis 1)
942!> \param la1_min Min L on center A (basis 1)
943!> \param npgfa1 Number of primitives on center A (basis 1)
944!> \param rpgfa1 Range of functions on A, used for screening (basis 1)
945!> \param zeta1 Exponents on center A (basis 1)
946!> \param la2_max Max L on center A (basis 2)
947!> \param la2_min Min L on center A (basis 2)
948!> \param npgfa2 Number of primitives on center A (basis 2)
949!> \param rpgfa2 Range of functions on A, used for screening (basis 2)
950!> \param zeta2 Exponents on center A (basis 2)
951!> \param lb_max Max L on center B
952!> \param lb_min Min L on center B
953!> \param npgfb Number of primitives on center B
954!> \param rpgfb Range of functions on B, used for screening
955!> \param zetb Exponents on center B
956!> \param rab Distance vector A-B
957!> \param saab Final overlap integrals
958!> \param daab First derivative overlap integrals
959!> \param saba Final overlap integrals; different order
960!> \param daba First derivative overlap integrals; different order
961!> \date 01.07.2014
962!> \author JGH
963! **************************************************************************************************
964 SUBROUTINE overlap_aab(la1_max, la1_min, npgfa1, rpgfa1, zeta1, &
965 la2_max, la2_min, npgfa2, rpgfa2, zeta2, &
966 lb_max, lb_min, npgfb, rpgfb, zetb, &
967 rab, saab, daab, saba, daba)
968 INTEGER, INTENT(IN) :: la1_max, la1_min, npgfa1
969 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa1, zeta1
970 INTEGER, INTENT(IN) :: la2_max, la2_min, npgfa2
971 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa2, zeta2
972 INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
973 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
974 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
975 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
976 OPTIONAL :: saab
977 REAL(kind=dp), DIMENSION(:, :, :, :), &
978 INTENT(INOUT), OPTIONAL :: daab
979 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
980 OPTIONAL :: saba
981 REAL(kind=dp), DIMENSION(:, :, :, :), &
982 INTENT(INOUT), OPTIONAL :: daba
983
984 INTEGER :: ax, ax1, ax2, ay, ay1, ay2, az, az1, az2, bx, by, bz, coa1, coa2, cob, i1pgf, &
985 i2pgf, ia1, ia2, ib, jpgf, la1, la2, lb, ldrr, lma, lmb, ma1, ma2, mb, na1, na2, nb, &
986 ofa1, ofa2, ofb
987 REAL(kind=dp) :: a, b, dumx, dumy, dumz, f0, rab2, rpgfa, &
988 tab, xhi, zet
989 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
990 REAL(kind=dp), DIMENSION(3) :: rap, rbp
991
992 ! Distance of the centers a and b
993
994 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
995 tab = sqrt(rab2)
996
997 ! Maximum l for auxiliary integrals
998 IF (PRESENT(saab) .OR. PRESENT(saba)) THEN
999 lma = la1_max + la2_max
1000 lmb = lb_max
1001 END IF
1002 IF (PRESENT(daab) .OR. PRESENT(daba)) THEN
1003 lma = la1_max + la2_max + 1
1004 lmb = lb_max
1005 END IF
1006 ldrr = max(lma, lmb) + 1
1007
1008 ! Allocate space for auxiliary integrals
1009 ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
1010
1011 ! Number of integrals, check size of arrays
1012 ofa1 = ncoset(la1_min - 1)
1013 ofa2 = ncoset(la2_min - 1)
1014 ofb = ncoset(lb_min - 1)
1015 na1 = ncoset(la1_max) - ofa1
1016 na2 = ncoset(la2_max) - ofa2
1017 nb = ncoset(lb_max) - ofb
1018 IF (PRESENT(saab)) THEN
1019 cpassert((SIZE(saab, 1) >= na1*npgfa1))
1020 cpassert((SIZE(saab, 2) >= na2*npgfa2))
1021 cpassert((SIZE(saab, 3) >= nb*npgfb))
1022 END IF
1023 IF (PRESENT(daab)) THEN
1024 cpassert((SIZE(daab, 1) >= na1*npgfa1))
1025 cpassert((SIZE(daab, 2) >= na2*npgfa2))
1026 cpassert((SIZE(daab, 3) >= nb*npgfb))
1027 cpassert((SIZE(daab, 4) >= 3))
1028 END IF
1029 IF (PRESENT(saba)) THEN
1030 cpassert((SIZE(saba, 1) >= na1*npgfa1))
1031 cpassert((SIZE(saba, 2) >= nb*npgfb))
1032 cpassert((SIZE(saba, 3) >= na2*npgfa2))
1033 END IF
1034 IF (PRESENT(daba)) THEN
1035 cpassert((SIZE(daba, 1) >= na1*npgfa1))
1036 cpassert((SIZE(daba, 2) >= nb*npgfb))
1037 cpassert((SIZE(daba, 3) >= na2*npgfa2))
1038 cpassert((SIZE(daba, 4) >= 3))
1039 END IF
1040
1041 ! Loops over all primitive Gaussian-type functions
1042 ma1 = 0
1043 DO i1pgf = 1, npgfa1
1044 ma2 = 0
1045 DO i2pgf = 1, npgfa2
1046 rpgfa = min(rpgfa1(i1pgf), rpgfa2(i2pgf))
1047 mb = 0
1048 DO jpgf = 1, npgfb
1049 ! Distance Screening
1050 IF (rpgfa + rpgfb(jpgf) < tab) THEN
1051 IF (PRESENT(saab)) saab(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, mb + 1:mb + nb) = 0.0_dp
1052 IF (PRESENT(daab)) daab(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, mb + 1:mb + nb, 1:3) = 0.0_dp
1053 IF (PRESENT(saba)) saba(ma1 + 1:ma1 + na1, mb + 1:mb + nb, ma2 + 1:ma2 + na2) = 0.0_dp
1054 IF (PRESENT(daba)) daba(ma1 + 1:ma1 + na1, mb + 1:mb + nb, ma2 + 1:ma2 + na2, 1:3) = 0.0_dp
1055 mb = mb + nb
1056 cycle
1057 END IF
1058
1059 ! Calculate some prefactors
1060 a = zeta1(i1pgf) + zeta2(i2pgf)
1061 b = zetb(jpgf)
1062 zet = a + b
1063 xhi = a*b/zet
1064 rap = b*rab/zet
1065 rbp = -a*rab/zet
1066
1067 ! [ss|s] integral
1068 f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
1069
1070 ! Calculate the recurrence relation
1071 CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
1072
1073 DO lb = lb_min, lb_max
1074 DO bx = 0, lb
1075 DO by = 0, lb - bx
1076 bz = lb - bx - by
1077 cob = coset(bx, by, bz) - ofb
1078 ib = mb + cob
1079 DO la2 = la2_min, la2_max
1080 DO ax2 = 0, la2
1081 DO ay2 = 0, la2 - ax2
1082 az2 = la2 - ax2 - ay2
1083 coa2 = coset(ax2, ay2, az2) - ofa2
1084 ia2 = ma2 + coa2
1085 DO la1 = la1_min, la1_max
1086 DO ax1 = 0, la1
1087 DO ay1 = 0, la1 - ax1
1088 az1 = la1 - ax1 - ay1
1089 coa1 = coset(ax1, ay1, az1) - ofa1
1090 ia1 = ma1 + coa1
1091 ! integrals
1092 IF (PRESENT(saab)) THEN
1093 saab(ia1, ia2, ib) = f0*rr(ax1 + ax2, bx, 1)*rr(ay1 + ay2, by, 2)*rr(az1 + az2, bz, 3)
1094 END IF
1095 IF (PRESENT(saba)) THEN
1096 saba(ia1, ib, ia2) = f0*rr(ax1 + ax2, bx, 1)*rr(ay1 + ay2, by, 2)*rr(az1 + az2, bz, 3)
1097 END IF
1098 ! first derivatives
1099 IF (PRESENT(daab)) THEN
1100 ax = ax1 + ax2
1101 ay = ay1 + ay2
1102 az = az1 + az2
1103 ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
1104 ! dx
1105 dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
1106 IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
1107 daab(ia1, ia2, ib, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
1108 ! dy
1109 dumy = 2.0_dp*a*rr(ay + 1, by, 2)
1110 IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
1111 daab(ia1, ia2, ib, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
1112 ! dz
1113 dumz = 2.0_dp*a*rr(az + 1, bz, 3)
1114 IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
1115 daab(ia1, ia2, ib, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
1116 END IF
1117 IF (PRESENT(daba)) THEN
1118 ax = ax1 + ax2
1119 ay = ay1 + ay2
1120 az = az1 + az2
1121 ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
1122 ! dx
1123 dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
1124 IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
1125 daba(ia1, ib, ia2, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
1126 ! dy
1127 dumy = 2.0_dp*a*rr(ay + 1, by, 2)
1128 IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
1129 daba(ia1, ib, ia2, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
1130 ! dz
1131 dumz = 2.0_dp*a*rr(az + 1, bz, 3)
1132 IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
1133 daba(ia1, ib, ia2, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
1134 END IF
1135 !
1136 END DO
1137 END DO
1138 END DO !la1
1139 END DO
1140 END DO
1141 END DO !la2
1142 END DO
1143 END DO
1144 END DO !lb
1145
1146 mb = mb + nb
1147 END DO
1148 ma2 = ma2 + na2
1149 END DO
1150 ma1 = ma1 + na1
1151 END DO
1152
1153 DEALLOCATE (rr)
1154
1155 END SUBROUTINE overlap_aab
1156
1157! **************************************************************************************************
1158!> \brief Calculation of the two-center overlap integrals [a|bb] over
1159!> Cartesian Gaussian-type functions.
1160!> \param la_max Max L on center A
1161!> \param la_min Min L on center A
1162!> \param npgfa Number of primitives on center A
1163!> \param rpgfa Range of functions on A, used for screening
1164!> \param zeta Exponents on center A
1165!> \param lb1_max Max L on center B (basis 1)
1166!> \param lb1_min Min L on center B (basis 1)
1167!> \param npgfb1 Number of primitives on center B (basis 1)
1168!> \param rpgfb1 Range of functions on B, used for screening (basis 1)
1169!> \param zetb1 Exponents on center B (basis 1)
1170!> \param lb2_max Max L on center B (basis 2)
1171!> \param lb2_min Min L on center B (basis 2)
1172!> \param npgfb2 Number of primitives on center B (basis 2)
1173!> \param rpgfb2 Range of functions on B, used for screening (basis 2)
1174!> \param zetb2 Exponents on center B (basis 2)
1175!> \param rab Distance vector A-B
1176!> \param sabb Final overlap integrals
1177!> \param dabb First derivative overlap integrals
1178!> \date 01.07.2014
1179!> \author JGH
1180! **************************************************************************************************
1181 SUBROUTINE overlap_abb(la_max, la_min, npgfa, rpgfa, zeta, &
1182 lb1_max, lb1_min, npgfb1, rpgfb1, zetb1, &
1183 lb2_max, lb2_min, npgfb2, rpgfb2, zetb2, &
1184 rab, sabb, dabb)
1185 INTEGER, INTENT(IN) :: la_max, la_min, npgfa
1186 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
1187 INTEGER, INTENT(IN) :: lb1_max, lb1_min, npgfb1
1188 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb1, zetb1
1189 INTEGER, INTENT(IN) :: lb2_max, lb2_min, npgfb2
1190 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb2, zetb2
1191 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
1192 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
1193 OPTIONAL :: sabb
1194 REAL(kind=dp), DIMENSION(:, :, :, :), &
1195 INTENT(INOUT), OPTIONAL :: dabb
1196
1197 INTEGER :: ax, ay, az, bx, bx1, bx2, by, by1, by2, bz, bz1, bz2, coa, cob1, cob2, ia, ib1, &
1198 ib2, ipgf, j1pgf, j2pgf, la, lb1, lb2, ldrr, lma, lmb, ma, mb1, mb2, na, nb1, nb2, ofa, &
1199 ofb1, ofb2
1200 REAL(kind=dp) :: a, b, dumx, dumy, dumz, f0, rab2, rpgfb, &
1201 tab, xhi, zet
1202 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
1203 REAL(kind=dp), DIMENSION(3) :: rap, rbp
1204
1205 ! Distance of the centers a and b
1206
1207 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1208 tab = sqrt(rab2)
1209
1210 ! Maximum l for auxiliary integrals
1211 IF (PRESENT(sabb)) THEN
1212 lma = la_max
1213 lmb = lb1_max + lb2_max
1214 END IF
1215 IF (PRESENT(dabb)) THEN
1216 lma = la_max + 1
1217 lmb = lb1_max + lb2_max
1218 END IF
1219 ldrr = max(lma, lmb) + 1
1220
1221 ! Allocate space for auxiliary integrals
1222 ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
1223
1224 ! Number of integrals, check size of arrays
1225 ofa = ncoset(la_min - 1)
1226 ofb1 = ncoset(lb1_min - 1)
1227 ofb2 = ncoset(lb2_min - 1)
1228 na = ncoset(la_max) - ofa
1229 nb1 = ncoset(lb1_max) - ofb1
1230 nb2 = ncoset(lb2_max) - ofb2
1231 IF (PRESENT(sabb)) THEN
1232 cpassert((SIZE(sabb, 1) >= na*npgfa))
1233 cpassert((SIZE(sabb, 2) >= nb1*npgfb1))
1234 cpassert((SIZE(sabb, 3) >= nb2*npgfb2))
1235 END IF
1236 IF (PRESENT(dabb)) THEN
1237 cpassert((SIZE(dabb, 1) >= na*npgfa))
1238 cpassert((SIZE(dabb, 2) >= nb1*npgfb1))
1239 cpassert((SIZE(dabb, 3) >= nb2*npgfb2))
1240 cpassert((SIZE(dabb, 4) >= 3))
1241 END IF
1242
1243 ! Loops over all pairs of primitive Gaussian-type functions
1244 ma = 0
1245 DO ipgf = 1, npgfa
1246 mb1 = 0
1247 DO j1pgf = 1, npgfb1
1248 mb2 = 0
1249 DO j2pgf = 1, npgfb2
1250 ! Distance Screening
1251 rpgfb = min(rpgfb1(j1pgf), rpgfb2(j2pgf))
1252 IF (rpgfa(ipgf) + rpgfb < tab) THEN
1253 IF (PRESENT(sabb)) sabb(ma + 1:ma + na, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2) = 0.0_dp
1254 IF (PRESENT(dabb)) dabb(ma + 1:ma + na, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2, 1:3) = 0.0_dp
1255 mb2 = mb2 + nb2
1256 cycle
1257 END IF
1258
1259 ! Calculate some prefactors
1260 a = zeta(ipgf)
1261 b = zetb1(j1pgf) + zetb2(j2pgf)
1262 zet = a + b
1263 xhi = a*b/zet
1264 rap = b*rab/zet
1265 rbp = -a*rab/zet
1266
1267 ! [s|s] integral
1268 f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
1269
1270 ! Calculate the recurrence relation
1271 CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
1272
1273 DO lb2 = lb2_min, lb2_max
1274 DO bx2 = 0, lb2
1275 DO by2 = 0, lb2 - bx2
1276 bz2 = lb2 - bx2 - by2
1277 cob2 = coset(bx2, by2, bz2) - ofb2
1278 ib2 = mb2 + cob2
1279 DO lb1 = lb1_min, lb1_max
1280 DO bx1 = 0, lb1
1281 DO by1 = 0, lb1 - bx1
1282 bz1 = lb1 - bx1 - by1
1283 cob1 = coset(bx1, by1, bz1) - ofb1
1284 ib1 = mb1 + cob1
1285 DO la = la_min, la_max
1286 DO ax = 0, la
1287 DO ay = 0, la - ax
1288 az = la - ax - ay
1289 coa = coset(ax, ay, az) - ofa
1290 ia = ma + coa
1291 ! integrals
1292 IF (PRESENT(sabb)) THEN
1293 sabb(ia, ib1, ib2) = f0*rr(ax, bx1 + bx2, 1)*rr(ay, by1 + by2, 2)*rr(az, bz1 + bz2, 3)
1294 END IF
1295 ! first derivatives
1296 IF (PRESENT(dabb)) THEN
1297 bx = bx1 + bx2
1298 by = by1 + by2
1299 bz = bz1 + bz2
1300 ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
1301 ! dx
1302 dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
1303 IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
1304 dabb(ia, ib1, ib2, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
1305 ! dy
1306 dumy = 2.0_dp*a*rr(ay + 1, by, 2)
1307 IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
1308 dabb(ia, ib1, ib2, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
1309 ! dz
1310 dumz = 2.0_dp*a*rr(az + 1, bz, 3)
1311 IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
1312 dabb(ia, ib1, ib2, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
1313 END IF
1314 !
1315 END DO
1316 END DO
1317 END DO !la
1318 END DO
1319 END DO
1320 END DO !lb1
1321 END DO
1322 END DO
1323 END DO !lb2
1324
1325 mb2 = mb2 + nb2
1326 END DO
1327 mb1 = mb1 + nb1
1328 END DO
1329 ma = ma + na
1330 END DO
1331
1332 DEALLOCATE (rr)
1333
1334 END SUBROUTINE overlap_abb
1335
1336! **************************************************************************************************
1337!> \brief Calculation of the two-center overlap integrals [aaa|b] over
1338!> Cartesian Gaussian-type functions.
1339!> \param la1_max Max L on center A (basis 1)
1340!> \param la1_min Min L on center A (basis 1)
1341!> \param npgfa1 Number of primitives on center A (basis 1)
1342!> \param rpgfa1 Range of functions on A, used for screening (basis 1)
1343!> \param zeta1 Exponents on center A (basis 1)
1344!> \param la2_max Max L on center A (basis 2)
1345!> \param la2_min Min L on center A (basis 2)
1346!> \param npgfa2 Number of primitives on center A (basis 2)
1347!> \param rpgfa2 Range of functions on A, used for screening (basis 2)
1348!> \param zeta2 Exponents on center A (basis 2)
1349!> \param la3_max Max L on center A (basis 3)
1350!> \param la3_min Min L on center A (basis 3)
1351!> \param npgfa3 Number of primitives on center A (basis 3)
1352!> \param rpgfa3 Range of functions on A, used for screening (basis 3)
1353!> \param zeta3 Exponents on center A (basis 3)
1354!> \param lb_max Max L on center B
1355!> \param lb_min Min L on center B
1356!> \param npgfb Number of primitives on center B
1357!> \param rpgfb Range of functions on B, used for screening
1358!> \param zetb Exponents on center B
1359!> \param rab Distance vector A-B
1360!> \param saaab Final overlap integrals
1361!> \param daaab First derivative overlap integrals
1362!> \date 01.07.2014
1363!> \author JGH
1364! **************************************************************************************************
1365 SUBROUTINE overlap_aaab(la1_max, la1_min, npgfa1, rpgfa1, zeta1, &
1366 la2_max, la2_min, npgfa2, rpgfa2, zeta2, &
1367 la3_max, la3_min, npgfa3, rpgfa3, zeta3, &
1368 lb_max, lb_min, npgfb, rpgfb, zetb, &
1369 rab, saaab, daaab)
1370 INTEGER, INTENT(IN) :: la1_max, la1_min, npgfa1
1371 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa1, zeta1
1372 INTEGER, INTENT(IN) :: la2_max, la2_min, npgfa2
1373 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa2, zeta2
1374 INTEGER, INTENT(IN) :: la3_max, la3_min, npgfa3
1375 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa3, zeta3
1376 INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
1377 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
1378 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
1379 REAL(kind=dp), DIMENSION(:, :, :, :), &
1380 INTENT(INOUT), OPTIONAL :: saaab
1381 REAL(kind=dp), DIMENSION(:, :, :, :, :), &
1382 INTENT(INOUT), OPTIONAL :: daaab
1383
1384 INTEGER :: ax, ax1, ax2, ax3, ay, ay1, ay2, ay3, az, az1, az2, az3, bx, by, bz, coa1, coa2, &
1385 coa3, cob, i1pgf, i2pgf, i3pgf, ia1, ia2, ia3, ib, jpgf, la1, la2, la3, lb, ldrr, lma, &
1386 lmb, ma1, ma2, ma3, mb, na1, na2, na3, nb, ofa1, ofa2, ofa3, ofb
1387 REAL(kind=dp) :: a, b, dumx, dumy, dumz, f0, rab2, rpgfa, &
1388 tab, xhi, zet
1389 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
1390 REAL(kind=dp), DIMENSION(3) :: rap, rbp
1391
1392! Distance of the centers a and b
1393
1394 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1395 tab = sqrt(rab2)
1396
1397 ! Maximum l for auxiliary integrals
1398 IF (PRESENT(saaab)) THEN
1399 lma = la1_max + la2_max + la3_max
1400 lmb = lb_max
1401 END IF
1402 IF (PRESENT(daaab)) THEN
1403 lma = la1_max + la2_max + la3_max + 1
1404 lmb = lb_max
1405 END IF
1406 ldrr = max(lma, lmb) + 1
1407
1408 ! Allocate space for auxiliary integrals
1409 ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
1410
1411 ! Number of integrals, check size of arrays
1412 ofa1 = ncoset(la1_min - 1)
1413 ofa2 = ncoset(la2_min - 1)
1414 ofa3 = ncoset(la3_min - 1)
1415 ofb = ncoset(lb_min - 1)
1416 na1 = ncoset(la1_max) - ofa1
1417 na2 = ncoset(la2_max) - ofa2
1418 na3 = ncoset(la3_max) - ofa3
1419 nb = ncoset(lb_max) - ofb
1420 IF (PRESENT(saaab)) THEN
1421 cpassert((SIZE(saaab, 1) >= na1*npgfa1))
1422 cpassert((SIZE(saaab, 2) >= na2*npgfa2))
1423 cpassert((SIZE(saaab, 3) >= na3*npgfa3))
1424 cpassert((SIZE(saaab, 4) >= nb*npgfb))
1425 END IF
1426 IF (PRESENT(daaab)) THEN
1427 cpassert((SIZE(daaab, 1) >= na1*npgfa1))
1428 cpassert((SIZE(daaab, 2) >= na2*npgfa2))
1429 cpassert((SIZE(daaab, 3) >= na3*npgfa3))
1430 cpassert((SIZE(daaab, 4) >= nb*npgfb))
1431 cpassert((SIZE(daaab, 5) >= 3))
1432 END IF
1433
1434 ! Loops over all primitive Gaussian-type functions
1435 ma1 = 0
1436 DO i1pgf = 1, npgfa1
1437 ma2 = 0
1438 DO i2pgf = 1, npgfa2
1439 ma3 = 0
1440 DO i3pgf = 1, npgfa3
1441 rpgfa = min(rpgfa1(i1pgf), rpgfa2(i2pgf), rpgfa3(i3pgf))
1442 mb = 0
1443 DO jpgf = 1, npgfb
1444 ! Distance Screening
1445 IF (rpgfa + rpgfb(jpgf) < tab) THEN
1446 IF (PRESENT(saaab)) saaab(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, ma3 + 1:ma3 + na3, mb + 1:mb + nb) = 0.0_dp
1447 IF (PRESENT(daaab)) daaab(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, ma3 + 1:ma3 + na3, mb + 1:mb + nb, 1:3) = 0.0_dp
1448 mb = mb + nb
1449 cycle
1450 END IF
1451
1452 ! Calculate some prefactors
1453 a = zeta1(i1pgf) + zeta2(i2pgf) + zeta3(i3pgf)
1454 b = zetb(jpgf)
1455 zet = a + b
1456 xhi = a*b/zet
1457 rap = b*rab/zet
1458 rbp = -a*rab/zet
1459
1460 ! [sss|s] integral
1461 f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
1462
1463 ! Calculate the recurrence relation
1464 CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
1465
1466 DO lb = lb_min, lb_max
1467 DO bx = 0, lb
1468 DO by = 0, lb - bx
1469 bz = lb - bx - by
1470 cob = coset(bx, by, bz) - ofb
1471 ib = mb + cob
1472 DO la3 = la3_min, la3_max
1473 DO ax3 = 0, la3
1474 DO ay3 = 0, la3 - ax3
1475 az3 = la3 - ax3 - ay3
1476 coa3 = coset(ax3, ay3, az3) - ofa3
1477 ia3 = ma3 + coa3
1478 DO la2 = la2_min, la2_max
1479 DO ax2 = 0, la2
1480 DO ay2 = 0, la2 - ax2
1481 az2 = la2 - ax2 - ay2
1482 coa2 = coset(ax2, ay2, az2) - ofa2
1483 ia2 = ma2 + coa2
1484 DO la1 = la1_min, la1_max
1485 DO ax1 = 0, la1
1486 DO ay1 = 0, la1 - ax1
1487 az1 = la1 - ax1 - ay1
1488 coa1 = coset(ax1, ay1, az1) - ofa1
1489 ia1 = ma1 + coa1
1490 ! integrals
1491 IF (PRESENT(saaab)) THEN
1492 saaab(ia1, ia2, ia3, ib) = f0*rr(ax1 + ax2 + ax3, bx, 1)* &
1493 rr(ay1 + ay2 + ay3, by, 2)*rr(az1 + az2 + az3, bz, 3)
1494 END IF
1495 ! first derivatives
1496 IF (PRESENT(daaab)) THEN
1497 ax = ax1 + ax2 + ax3
1498 ay = ay1 + ay2 + ay3
1499 az = az1 + az2 + az3
1500 ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
1501 ! dx
1502 dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
1503 IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
1504 daaab(ia1, ia2, ia3, ib, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
1505 ! dy
1506 dumy = 2.0_dp*a*rr(ay + 1, by, 2)
1507 IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
1508 daaab(ia1, ia2, ia3, ib, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
1509 ! dz
1510 dumz = 2.0_dp*a*rr(az + 1, bz, 3)
1511 IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
1512 daaab(ia1, ia2, ia3, ib, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
1513 END IF
1514 !
1515 END DO
1516 END DO
1517 END DO !la1
1518 END DO
1519 END DO
1520 END DO !la2
1521 END DO
1522 END DO
1523 END DO !la3
1524 END DO
1525 END DO
1526 END DO !lb
1527
1528 mb = mb + nb
1529 END DO
1530 ma3 = ma3 + na3
1531 END DO
1532 ma2 = ma2 + na2
1533 END DO
1534 ma1 = ma1 + na1
1535 END DO
1536
1537 DEALLOCATE (rr)
1538
1539 END SUBROUTINE overlap_aaab
1540! **************************************************************************************************
1541!> \brief Calculation of the two-center overlap integrals [aa|bb] over
1542!> Cartesian Gaussian-type functions.
1543!> \param la1_max Max L on center A (basis 1)
1544!> \param la1_min Min L on center A (basis 1)
1545!> \param npgfa1 Number of primitives on center A (basis 1)
1546!> \param rpgfa1 Range of functions on A, used for screening (basis 1)
1547!> \param zeta1 Exponents on center A (basis 1)
1548!> \param la2_max Max L on center A (basis 2)
1549!> \param la2_min Min L on center A (basis 2)
1550!> \param npgfa2 Number of primitives on center A (basis 2)
1551!> \param rpgfa2 Range of functions on A, used for screening (basis 2)
1552!> \param zeta2 Exponents on center A (basis 2)
1553!> \param lb1_max Max L on center B (basis 3)
1554!> \param lb1_min Min L on center B (basis 3)
1555!> \param npgfb1 Number of primitives on center B (basis 3)
1556!> \param rpgfb1 Range of functions on B, used for screening (basis 3)
1557!> \param zetb1 Exponents on center B (basis 3)
1558!> \param lb2_max Max L on center B (basis 4)
1559!> \param lb2_min Min L on center B (basis 4)
1560!> \param npgfb2 Number of primitives on center B (basis 4)
1561!> \param rpgfb2 Range of functions on B, used for screening (basis 4)
1562!> \param zetb2 Exponents on center B (basis 4)
1563!> \param rab Distance vector A-B
1564!> \param saabb Final overlap integrals
1565!> \param daabb First derivative overlap integrals
1566!> \date 01.07.2014
1567!> \author JGH
1568! **************************************************************************************************
1569 SUBROUTINE overlap_aabb(la1_max, la1_min, npgfa1, rpgfa1, zeta1, &
1570 la2_max, la2_min, npgfa2, rpgfa2, zeta2, &
1571 lb1_max, lb1_min, npgfb1, rpgfb1, zetb1, &
1572 lb2_max, lb2_min, npgfb2, rpgfb2, zetb2, &
1573 rab, saabb, daabb)
1574 INTEGER, INTENT(IN) :: la1_max, la1_min, npgfa1
1575 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa1, zeta1
1576 INTEGER, INTENT(IN) :: la2_max, la2_min, npgfa2
1577 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa2, zeta2
1578 INTEGER, INTENT(IN) :: lb1_max, lb1_min, npgfb1
1579 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb1, zetb1
1580 INTEGER, INTENT(IN) :: lb2_max, lb2_min, npgfb2
1581 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb2, zetb2
1582 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
1583 REAL(kind=dp), DIMENSION(:, :, :, :), &
1584 INTENT(INOUT), OPTIONAL :: saabb
1585 REAL(kind=dp), DIMENSION(:, :, :, :, :), &
1586 INTENT(INOUT), OPTIONAL :: daabb
1587
1588 INTEGER :: ax, ax1, ax2, ay, ay1, ay2, az, az1, az2, bx, bx1, bx2, by, by1, by2, bz, bz1, &
1589 bz2, coa1, coa2, cob1, cob2, i1pgf, i2pgf, ia1, ia2, ib1, ib2, j1pgf, j2pgf, la1, la2, &
1590 lb1, lb2, ldrr, lma, lmb, ma1, ma2, mb1, mb2, na1, na2, nb1, nb2, ofa1, ofa2, ofb1, ofb2
1591 REAL(kind=dp) :: a, b, dumx, dumy, dumz, f0, rab2, rpgfa, &
1592 rpgfb, tab, xhi, zet
1593 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
1594 REAL(kind=dp), DIMENSION(3) :: rap, rbp
1595
1596! Distance of the centers a and b
1597
1598 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1599 tab = sqrt(rab2)
1600
1601 ! Maximum l for auxiliary integrals
1602 IF (PRESENT(saabb)) THEN
1603 lma = la1_max + la2_max
1604 lmb = lb1_max + lb2_max
1605 END IF
1606 IF (PRESENT(daabb)) THEN
1607 lma = la1_max + la2_max + 1
1608 lmb = lb1_max + lb2_max
1609 END IF
1610 ldrr = max(lma, lmb) + 1
1611
1612 ! Allocate space for auxiliary integrals
1613 ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
1614
1615 ! Number of integrals, check size of arrays
1616 ofa1 = ncoset(la1_min - 1)
1617 ofa2 = ncoset(la2_min - 1)
1618 ofb1 = ncoset(lb1_min - 1)
1619 ofb2 = ncoset(lb2_min - 1)
1620 na1 = ncoset(la1_max) - ofa1
1621 na2 = ncoset(la2_max) - ofa2
1622 nb1 = ncoset(lb1_max) - ofb1
1623 nb2 = ncoset(lb2_max) - ofb2
1624 IF (PRESENT(saabb)) THEN
1625 cpassert((SIZE(saabb, 1) >= na1*npgfa1))
1626 cpassert((SIZE(saabb, 2) >= na2*npgfa2))
1627 cpassert((SIZE(saabb, 3) >= nb1*npgfb1))
1628 cpassert((SIZE(saabb, 4) >= nb2*npgfb2))
1629 END IF
1630 IF (PRESENT(daabb)) THEN
1631 cpassert((SIZE(daabb, 1) >= na1*npgfa1))
1632 cpassert((SIZE(daabb, 2) >= na2*npgfa2))
1633 cpassert((SIZE(daabb, 3) >= nb1*npgfb1))
1634 cpassert((SIZE(daabb, 4) >= nb2*npgfb2))
1635 cpassert((SIZE(daabb, 5) >= 3))
1636 END IF
1637
1638 ! Loops over all primitive Gaussian-type functions
1639 ma1 = 0
1640 DO i1pgf = 1, npgfa1
1641 ma2 = 0
1642 DO i2pgf = 1, npgfa2
1643 mb1 = 0
1644 rpgfa = min(rpgfa1(i1pgf), rpgfa2(i2pgf))
1645 DO j1pgf = 1, npgfb1
1646 mb2 = 0
1647 DO j2pgf = 1, npgfb2
1648 rpgfb = min(rpgfb1(j1pgf), rpgfb2(j2pgf))
1649 ! Distance Screening
1650 IF (rpgfa + rpgfb < tab) THEN
1651 IF (PRESENT(saabb)) saabb(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2) = 0.0_dp
1652 IF (PRESENT(daabb)) daabb(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2, 1:3) = 0.0_dp
1653 mb2 = mb2 + nb2
1654 cycle
1655 END IF
1656
1657 ! Calculate some prefactors
1658 a = zeta1(i1pgf) + zeta2(i2pgf)
1659 b = zetb1(j1pgf) + zetb2(j2pgf)
1660 zet = a + b
1661 xhi = a*b/zet
1662 rap = b*rab/zet
1663 rbp = -a*rab/zet
1664
1665 ! [ss|ss] integral
1666 f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
1667
1668 ! Calculate the recurrence relation
1669 CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
1670
1671 DO lb2 = lb2_min, lb2_max
1672 DO bx2 = 0, lb2
1673 DO by2 = 0, lb2 - bx2
1674 bz2 = lb2 - bx2 - by2
1675 cob2 = coset(bx2, by2, bz2) - ofb2
1676 ib2 = mb2 + cob2
1677 DO lb1 = lb1_min, lb1_max
1678 DO bx1 = 0, lb1
1679 DO by1 = 0, lb1 - bx1
1680 bz1 = lb1 - bx1 - by1
1681 cob1 = coset(bx1, by1, bz1) - ofb1
1682 ib1 = mb1 + cob1
1683 DO la2 = la2_min, la2_max
1684 DO ax2 = 0, la2
1685 DO ay2 = 0, la2 - ax2
1686 az2 = la2 - ax2 - ay2
1687 coa2 = coset(ax2, ay2, az2) - ofa2
1688 ia2 = ma2 + coa2
1689 DO la1 = la1_min, la1_max
1690 DO ax1 = 0, la1
1691 DO ay1 = 0, la1 - ax1
1692 az1 = la1 - ax1 - ay1
1693 coa1 = coset(ax1, ay1, az1) - ofa1
1694 ia1 = ma1 + coa1
1695 ! integrals
1696 IF (PRESENT(saabb)) THEN
1697 saabb(ia1, ia2, ib1, ib2) = f0*rr(ax1 + ax2, bx1 + bx2, 1)* &
1698 rr(ay1 + ay2, by1 + by2, 2)*rr(az1 + az2, bz1 + bz2, 3)
1699 END IF
1700 ! first derivatives
1701 IF (PRESENT(daabb)) THEN
1702 ax = ax1 + ax2
1703 ay = ay1 + ay2
1704 az = az1 + az2
1705 bx = bx1 + bx2
1706 by = by1 + by2
1707 bz = bz1 + bz2
1708 ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
1709 ! dx
1710 dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
1711 IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
1712 daabb(ia1, ia2, ib1, ib2, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
1713 ! dy
1714 dumy = 2.0_dp*a*rr(ay + 1, by, 2)
1715 IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
1716 daabb(ia1, ia2, ib1, ib2, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
1717 ! dz
1718 dumz = 2.0_dp*a*rr(az + 1, bz, 3)
1719 IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
1720 daabb(ia1, ia2, ib1, ib2, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
1721 END IF
1722 !
1723 END DO
1724 END DO
1725 END DO !la1
1726 END DO
1727 END DO
1728 END DO !la2
1729 END DO
1730 END DO
1731 END DO !lb1
1732 END DO
1733 END DO
1734 END DO !lb2
1735
1736 mb2 = mb2 + nb2
1737 END DO
1738 mb1 = mb1 + nb1
1739 END DO
1740 ma2 = ma2 + na2
1741 END DO
1742 ma1 = ma1 + na1
1743 END DO
1744
1745 DEALLOCATE (rr)
1746
1747 END SUBROUTINE overlap_aabb
1748! **************************************************************************************************
1749!> \brief Calculation of the two-center overlap integrals [a|bbb] over
1750!> Cartesian Gaussian-type functions.
1751!> \param la_max Max L on center A
1752!> \param la_min Min L on center A
1753!> \param npgfa Number of primitives on center A
1754!> \param rpgfa Range of functions on A, used for screening
1755!> \param zeta Exponents on center A
1756!> \param lb1_max Max L on center B (basis 1)
1757!> \param lb1_min Min L on center B (basis 1)
1758!> \param npgfb1 Number of primitives on center B (basis 1)
1759!> \param rpgfb1 Range of functions on B, used for screening (basis 1)
1760!> \param zetb1 Exponents on center B (basis 1)
1761!> \param lb2_max Max L on center B (basis 2)
1762!> \param lb2_min Min L on center B (basis 2)
1763!> \param npgfb2 Number of primitives on center B (basis 2)
1764!> \param rpgfb2 Range of functions on B, used for screening (basis 2)
1765!> \param zetb2 Exponents on center B (basis 2)
1766!> \param lb3_max Max L on center B (basis 3)
1767!> \param lb3_min Min L on center B (basis 3)
1768!> \param npgfb3 Number of primitives on center B (basis 3)
1769!> \param rpgfb3 Range of functions on B, used for screening (basis 3)
1770!> \param zetb3 Exponents on center B (basis 3)
1771!> \param rab Distance vector A-B
1772!> \param sabbb Final overlap integrals
1773!> \param dabbb First derivative overlap integrals
1774!> \date 01.07.2014
1775!> \author JGH
1776! **************************************************************************************************
1777 SUBROUTINE overlap_abbb(la_max, la_min, npgfa, rpgfa, zeta, &
1778 lb1_max, lb1_min, npgfb1, rpgfb1, zetb1, &
1779 lb2_max, lb2_min, npgfb2, rpgfb2, zetb2, &
1780 lb3_max, lb3_min, npgfb3, rpgfb3, zetb3, &
1781 rab, sabbb, dabbb)
1782 INTEGER, INTENT(IN) :: la_max, la_min, npgfa
1783 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
1784 INTEGER, INTENT(IN) :: lb1_max, lb1_min, npgfb1
1785 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb1, zetb1
1786 INTEGER, INTENT(IN) :: lb2_max, lb2_min, npgfb2
1787 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb2, zetb2
1788 INTEGER, INTENT(IN) :: lb3_max, lb3_min, npgfb3
1789 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb3, zetb3
1790 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
1791 REAL(kind=dp), DIMENSION(:, :, :, :), &
1792 INTENT(INOUT), OPTIONAL :: sabbb
1793 REAL(kind=dp), DIMENSION(:, :, :, :, :), &
1794 INTENT(INOUT), OPTIONAL :: dabbb
1795
1796 INTEGER :: ax, ay, az, bx, bx1, bx2, bx3, by, by1, by2, by3, bz, bz1, bz2, bz3, coa, cob1, &
1797 cob2, cob3, ia, ib1, ib2, ib3, ipgf, j1pgf, j2pgf, j3pgf, la, lb1, lb2, lb3, ldrr, lma, &
1798 lmb, ma, mb1, mb2, mb3, na, nb1, nb2, nb3, ofa, ofb1, ofb2, ofb3
1799 REAL(kind=dp) :: a, b, dumx, dumy, dumz, f0, rab2, rpgfb, &
1800 tab, xhi, zet
1801 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
1802 REAL(kind=dp), DIMENSION(3) :: rap, rbp
1803
1804! Distance of the centers a and b
1805
1806 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1807 tab = sqrt(rab2)
1808
1809 ! Maximum l for auxiliary integrals
1810 IF (PRESENT(sabbb)) THEN
1811 lma = la_max
1812 lmb = lb1_max + lb2_max + lb3_max
1813 END IF
1814 IF (PRESENT(dabbb)) THEN
1815 lma = la_max + 1
1816 lmb = lb1_max + lb2_max + lb3_max
1817 END IF
1818 ldrr = max(lma, lmb) + 1
1819
1820 ! Allocate space for auxiliary integrals
1821 ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
1822
1823 ! Number of integrals, check size of arrays
1824 ofa = ncoset(la_min - 1)
1825 ofb1 = ncoset(lb1_min - 1)
1826 ofb2 = ncoset(lb2_min - 1)
1827 ofb3 = ncoset(lb3_min - 1)
1828 na = ncoset(la_max) - ofa
1829 nb1 = ncoset(lb1_max) - ofb1
1830 nb2 = ncoset(lb2_max) - ofb2
1831 nb3 = ncoset(lb3_max) - ofb3
1832 IF (PRESENT(sabbb)) THEN
1833 cpassert((SIZE(sabbb, 1) >= na*npgfa))
1834 cpassert((SIZE(sabbb, 2) >= nb1*npgfb1))
1835 cpassert((SIZE(sabbb, 3) >= nb2*npgfb2))
1836 cpassert((SIZE(sabbb, 4) >= nb3*npgfb3))
1837 END IF
1838 IF (PRESENT(dabbb)) THEN
1839 cpassert((SIZE(dabbb, 1) >= na*npgfa))
1840 cpassert((SIZE(dabbb, 2) >= nb1*npgfb1))
1841 cpassert((SIZE(dabbb, 3) >= nb2*npgfb2))
1842 cpassert((SIZE(dabbb, 4) >= nb3*npgfb3))
1843 cpassert((SIZE(dabbb, 5) >= 3))
1844 END IF
1845
1846 ! Loops over all pairs of primitive Gaussian-type functions
1847 ma = 0
1848 DO ipgf = 1, npgfa
1849 mb1 = 0
1850 DO j1pgf = 1, npgfb1
1851 mb2 = 0
1852 DO j2pgf = 1, npgfb2
1853 mb3 = 0
1854 DO j3pgf = 1, npgfb3
1855 ! Distance Screening
1856 rpgfb = min(rpgfb1(j1pgf), rpgfb2(j2pgf), rpgfb3(j3pgf))
1857 IF (rpgfa(ipgf) + rpgfb < tab) THEN
1858 IF (PRESENT(sabbb)) sabbb(ma + 1:ma + na, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2, mb3 + 1:mb3 + nb3) = 0.0_dp
1859 IF (PRESENT(dabbb)) dabbb(ma + 1:ma + na, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2, mb3 + 1:mb3 + nb3, 1:3) = 0.0_dp
1860 mb3 = mb3 + nb3
1861 cycle
1862 END IF
1863
1864 ! Calculate some prefactors
1865 a = zeta(ipgf)
1866 b = zetb1(j1pgf) + zetb2(j2pgf) + zetb3(j3pgf)
1867 zet = a + b
1868 xhi = a*b/zet
1869 rap = b*rab/zet
1870 rbp = -a*rab/zet
1871
1872 ! [s|s] integral
1873 f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
1874
1875 ! Calculate the recurrence relation
1876 CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
1877
1878 DO lb3 = lb3_min, lb3_max
1879 DO bx3 = 0, lb3
1880 DO by3 = 0, lb3 - bx3
1881 bz3 = lb3 - bx3 - by3
1882 cob3 = coset(bx3, by3, bz3) - ofb3
1883 ib3 = mb3 + cob3
1884 DO lb2 = lb2_min, lb2_max
1885 DO bx2 = 0, lb2
1886 DO by2 = 0, lb2 - bx2
1887 bz2 = lb2 - bx2 - by2
1888 cob2 = coset(bx2, by2, bz2) - ofb2
1889 ib2 = mb2 + cob2
1890 DO lb1 = lb1_min, lb1_max
1891 DO bx1 = 0, lb1
1892 DO by1 = 0, lb1 - bx1
1893 bz1 = lb1 - bx1 - by1
1894 cob1 = coset(bx1, by1, bz1) - ofb1
1895 ib1 = mb1 + cob1
1896 DO la = la_min, la_max
1897 DO ax = 0, la
1898 DO ay = 0, la - ax
1899 az = la - ax - ay
1900 coa = coset(ax, ay, az) - ofa
1901 ia = ma + coa
1902 ! integrals
1903 IF (PRESENT(sabbb)) THEN
1904 sabbb(ia, ib1, ib2, ib3) = f0*rr(ax, bx1 + bx2 + bx3, 1)* &
1905 rr(ay, by1 + by2 + by3, 2)*rr(az, bz1 + bz2 + bz3, 3)
1906 END IF
1907 ! first derivatives
1908 IF (PRESENT(dabbb)) THEN
1909 bx = bx1 + bx2 + bx3
1910 by = by1 + by2 + by3
1911 bz = bz1 + bz2 + bz3
1912 ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
1913 ! dx
1914 dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
1915 IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
1916 dabbb(ia, ib1, ib2, ib3, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
1917 ! dy
1918 dumy = 2.0_dp*a*rr(ay + 1, by, 2)
1919 IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
1920 dabbb(ia, ib1, ib2, ib3, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
1921 ! dz
1922 dumz = 2.0_dp*a*rr(az + 1, bz, 3)
1923 IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
1924 dabbb(ia, ib1, ib2, ib3, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
1925 END IF
1926 !
1927 END DO
1928 END DO
1929 END DO !la
1930 END DO
1931 END DO
1932 END DO !lb1
1933 END DO
1934 END DO
1935 END DO !lb2
1936 END DO
1937 END DO
1938 END DO !lb3
1939
1940 mb3 = mb3 + nb3
1941 END DO
1942 mb2 = mb2 + nb2
1943 END DO
1944 mb1 = mb1 + nb1
1945 END DO
1946 ma = ma + na
1947 END DO
1948
1949 DEALLOCATE (rr)
1950
1951 END SUBROUTINE overlap_abbb
1952
1953! **************************************************************************************************
1954!> \brief Calculation of the two-center overlap integrals [a|b] over
1955!> Spherical Gaussian-type functions.
1956!> \param la Max L on center A
1957!> \param zeta Exponents on center A
1958!> \param lb Max L on center B
1959!> \param zetb Exponents on center B
1960!> \param rab Distance vector A-B
1961!> \param sab Final overlap integrals
1962!> \date 01.03.2016
1963!> \author JGH
1964! **************************************************************************************************
1965 SUBROUTINE overlap_ab_s(la, zeta, lb, zetb, rab, sab)
1966 INTEGER, INTENT(IN) :: la
1967 REAL(kind=dp), INTENT(IN) :: zeta
1968 INTEGER, INTENT(IN) :: lb
1969 REAL(kind=dp), INTENT(IN) :: zetb
1970 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
1971 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
1972
1973 REAL(kind=dp), PARAMETER :: huge4 = huge(1._dp)/4._dp
1974
1975 INTEGER :: nca, ncb, nsa, nsb
1976 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cab
1977 REAL(kind=dp), DIMENSION(1) :: rpgf, za, zb
1978 REAL(kind=dp), DIMENSION(:, :), POINTER :: c2sa, c2sb
1979
1980 rpgf(1) = huge4
1981 za(1) = zeta
1982 zb(1) = zetb
1983
1984 nca = nco(la)
1985 ncb = nco(lb)
1986 ALLOCATE (cab(nca, ncb))
1987 nsa = nso(la)
1988 nsb = nso(lb)
1989
1990 CALL overlap_ab(la, la, 1, rpgf, za, lb, lb, 1, rpgf, zb, rab, cab)
1991
1992 c2sa => orbtramat(la)%c2s
1993 c2sb => orbtramat(lb)%c2s
1994 sab(1:nsa, 1:nsb) = matmul(c2sa(1:nsa, 1:nca), &
1995 matmul(cab(1:nca, 1:ncb), transpose(c2sb(1:nsb, 1:ncb))))
1996
1997 DEALLOCATE (cab)
1998
1999 END SUBROUTINE overlap_ab_s
2000
2001! **************************************************************************************************
2002!> \brief Calculation of the overlap integrals [a|b] over
2003!> cubic periodic Spherical Gaussian-type functions.
2004!> \param la Max L on center A
2005!> \param zeta Exponents on center A
2006!> \param lb Max L on center B
2007!> \param zetb Exponents on center B
2008!> \param alat Lattice constant
2009!> \param sab Final overlap integrals
2010!> \date 01.03.2016
2011!> \author JGH
2012! **************************************************************************************************
2013 SUBROUTINE overlap_ab_sp(la, zeta, lb, zetb, alat, sab)
2014 INTEGER, INTENT(IN) :: la
2015 REAL(kind=dp), INTENT(IN) :: zeta
2016 INTEGER, INTENT(IN) :: lb
2017 REAL(kind=dp), INTENT(IN) :: zetb, alat
2018 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
2019
2020 COMPLEX(KIND=dp) :: zfg
2021 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: fun, gun
2022 INTEGER :: ax, ay, az, bx, by, bz, i, ia, ib, l, &
2023 l1, l2, na, nb, nca, ncb, nmax, nsa, &
2024 nsb
2025 REAL(kind=dp) :: oa, ob, ovol, zm
2026 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fexp, gexp, gval
2027 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cab
2028 REAL(kind=dp), DIMENSION(0:3, 0:3) :: fgsum
2029 REAL(kind=dp), DIMENSION(:, :), POINTER :: c2sa, c2sb
2030
2031 nca = nco(la)
2032 ncb = nco(lb)
2033 ALLOCATE (cab(nca, ncb))
2034 cab = 0.0_dp
2035 nsa = nso(la)
2036 nsb = nso(lb)
2037
2038 zm = min(zeta, zetb)
2039 nmax = nint(1.81_dp*alat*sqrt(zm) + 1.0_dp)
2040 ALLOCATE (fun(-nmax:nmax, 0:la), gun(-nmax:nmax, 0:lb), &
2041 fexp(-nmax:nmax), gexp(-nmax:nmax), gval(-nmax:nmax))
2042
2043 oa = 1._dp/zeta
2044 ob = 1._dp/zetb
2045 DO i = -nmax, nmax
2046 gval(i) = twopi/alat*real(i, kind=dp)
2047 fexp(i) = sqrt(oa*pi)*exp(-0.25_dp*oa*gval(i)**2)
2048 gexp(i) = sqrt(ob*pi)*exp(-0.25_dp*ob*gval(i)**2)
2049 END DO
2050 DO l = 0, la
2051 IF (l == 0) THEN
2052 fun(:, l) = z_one
2053 ELSEIF (l == 1) THEN
2054 fun(:, l) = cmplx(0.0_dp, 0.5_dp*oa*gval(:), kind=dp)
2055 ELSEIF (l == 2) THEN
2056 fun(:, l) = cmplx(-(0.5_dp*oa*gval(:))**2, 0.0_dp, kind=dp)
2057 fun(:, l) = fun(:, l) + cmplx(0.5_dp*oa, 0.0_dp, kind=dp)
2058 ELSEIF (l == 3) THEN
2059 fun(:, l) = cmplx(0.0_dp, -(0.5_dp*oa*gval(:))**3, kind=dp)
2060 fun(:, l) = fun(:, l) + cmplx(0.0_dp, 0.75_dp*oa*oa*gval(:), kind=dp)
2061 ELSE
2062 cpabort("l value too high")
2063 END IF
2064 END DO
2065 DO l = 0, lb
2066 IF (l == 0) THEN
2067 gun(:, l) = z_one
2068 ELSEIF (l == 1) THEN
2069 gun(:, l) = cmplx(0.0_dp, 0.5_dp*ob*gval(:), kind=dp)
2070 ELSEIF (l == 2) THEN
2071 gun(:, l) = cmplx(-(0.5_dp*ob*gval(:))**2, 0.0_dp, kind=dp)
2072 gun(:, l) = gun(:, l) + cmplx(0.5_dp*ob, 0.0_dp, kind=dp)
2073 ELSEIF (l == 3) THEN
2074 gun(:, l) = cmplx(0.0_dp, -(0.5_dp*ob*gval(:))**3, kind=dp)
2075 gun(:, l) = gun(:, l) + cmplx(0.0_dp, 0.75_dp*ob*ob*gval(:), kind=dp)
2076 ELSE
2077 cpabort("l value too high")
2078 END IF
2079 END DO
2080
2081 fgsum = 0.0_dp
2082 DO l1 = 0, la
2083 DO l2 = 0, lb
2084 zfg = sum(conjg(fun(:, l1))*fexp(:)*gun(:, l2)*gexp(:))
2085 fgsum(l1, l2) = real(zfg, kind=dp)
2086 END DO
2087 END DO
2088
2089 na = ncoset(la - 1)
2090 nb = ncoset(lb - 1)
2091 DO ax = 0, la
2092 DO ay = 0, la - ax
2093 az = la - ax - ay
2094 ia = coset(ax, ay, az) - na
2095 DO bx = 0, lb
2096 DO by = 0, lb - bx
2097 bz = lb - bx - by
2098 ib = coset(bx, by, bz) - nb
2099 cab(ia, ib) = fgsum(ax, bx)*fgsum(ay, by)*fgsum(az, bz)
2100 END DO
2101 END DO
2102 END DO
2103 END DO
2104
2105 c2sa => orbtramat(la)%c2s
2106 c2sb => orbtramat(lb)%c2s
2107 sab(1:nsa, 1:nsb) = matmul(c2sa(1:nsa, 1:nca), &
2108 matmul(cab(1:nca, 1:ncb), transpose(c2sb(1:nsb, 1:ncb))))
2109 ovol = 1._dp/(alat**3)
2110 sab(1:nsa, 1:nsb) = ovol*sab(1:nsa, 1:nsb)
2111
2112 DEALLOCATE (cab, fun, gun, fexp, gexp, gval)
2113
2114 END SUBROUTINE overlap_ab_sp
2115
2116END MODULE ai_overlap
ai_os_rr
Definition ai_os_rr.F:8
ai_os_rr::os_rr_ovlp
subroutine, public os_rr_ovlp(rap, la_max, rbp, lb_max, zet, ldrr, rr)
Calculation of the basic Obara-Saika recurrence relation.
Definition ai_os_rr.F:39
ai_overlap
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition ai_overlap.F:18
ai_overlap::overlap_ab
subroutine, public overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, sab, dab, ddab)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
Definition ai_overlap.F:681
ai_overlap::overlap
subroutine, public overlap(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, rab, dab, sab, da_max_set, return_derivatives, s, lds, sdab, pab, force_a)
Purpose: Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions...
Definition ai_overlap.F:74
ai_overlap::overlap_ab_s
subroutine, public overlap_ab_s(la, zeta, lb, zetb, rab, sab)
Calculation of the two-center overlap integrals [a|b] over Spherical Gaussian-type functions.
Definition ai_overlap.F:1966
ai_overlap::overlap_ab_sp
subroutine, public overlap_ab_sp(la, zeta, lb, zetb, alat, sab)
Calculation of the overlap integrals [a|b] over cubic periodic Spherical Gaussian-type functions.
Definition ai_overlap.F:2014
ai_overlap::overlap_aab
subroutine, public overlap_aab(la1_max, la1_min, npgfa1, rpgfa1, zeta1, la2_max, la2_min, npgfa2, rpgfa2, zeta2, lb_max, lb_min, npgfb, rpgfb, zetb, rab, saab, daab, saba, daba)
Calculation of the two-center overlap integrals [aa|b] over Cartesian Gaussian-type functions.
Definition ai_overlap.F:968
ai_overlap::overlap_abb
subroutine, public overlap_abb(la_max, la_min, npgfa, rpgfa, zeta, lb1_max, lb1_min, npgfb1, rpgfb1, zetb1, lb2_max, lb2_min, npgfb2, rpgfb2, zetb2, rab, sabb, dabb)
Calculation of the two-center overlap integrals [a|bb] over Cartesian Gaussian-type functions.
Definition ai_overlap.F:1185
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::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:143
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition orbital_pointers.F:42
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition orbital_pointers.F:45
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition orbital_pointers.F:42
orbital_transformation_matrices
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
Definition orbital_transformation_matrices.F:17
orbital_transformation_matrices::orbtramat
type(orbtramat_type), dimension(:), pointer, public orbtramat
Definition orbital_transformation_matrices.F:47