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