(git:374b731)
Loading...
Searching...
No Matches
ai_overlap3.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!!****** cp2k/ai_overlap3 [1.0] *
8!!
9!! NAME
10!! ai_overlap3
11!!
12!! FUNCTION
13!! Calculation of three-center overlap integrals over Cartesian
14!! Gaussian-type functions.
15!!
16!! AUTHOR
17!! Matthias Krack (26.06.2001)
18!!
19!! LITERATURE
20!! S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
21!!
22!******************************************************************************
23
25
26! **************************************************************************************************
27
28! ax,ay,az : Angular momentum index numbers of orbital a.
29! bx,by,bz : Angular momentum index numbers of orbital b.
30! coset : Cartesian orbital set pointer.
31! dab : Distance between the atomic centers a and b.
32! dac : Distance between the atomic centers a and c.
33! dbc : Distance between the atomic centers b and c.
34! l{a,b,c} : Angular momentum quantum number of shell a, b or c.
35! l{a,b}_max : Maximum angular momentum quantum number of shell a, b or c.
36! ncoset : Number of Cartesian orbitals up to l.
37! rab : Distance vector between the atomic centers a and b.
38! rac : Distance vector between the atomic centers a and c.
39! rbc : Distance vector between the atomic centers b and c.
40! rpgf{a,b,c}: Radius of the primitive Gaussian-type function a or b.
41! zet{a,b,c} : Exponents of the Gaussian-type functions a or b.
42! zetg : Reciprocal of the sum of the exponents of orbital a, b and c.
43! zetp : Reciprocal of the sum of the exponents of orbital a and b.
44
45! **************************************************************************************************
46
47 USE kinds, ONLY: dp
48 USE mathconstants, ONLY: pi
49 USE orbital_pointers, ONLY: coset,&
50 ncoset
51#include "../base/base_uses.f90"
52
53 IMPLICIT NONE
54
55 PRIVATE
56
57 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap3'
58
59! *** Public subroutines ***
60
61 PUBLIC :: overlap3
62
63!!***
64! **************************************************************************************************
65
66CONTAINS
67
68! ***************************************************************************************************
69!> \brief Calculation of three-center overlap integrals [a|b|c] over primitive
70!> Cartesian Gaussian functions
71!> \param la_max_set ...
72!> \param npgfa ...
73!> \param zeta ...
74!> \param rpgfa ...
75!> \param la_min_set ...
76!> \param lb_max_set ...
77!> \param npgfb ...
78!> \param zetb ...
79!> \param rpgfb ...
80!> \param lb_min_set ...
81!> \param lc_max_set ...
82!> \param npgfc ...
83!> \param zetc ...
84!> \param rpgfc ...
85!> \param lc_min_set ...
86!> \param rab ...
87!> \param dab ...
88!> \param rac ...
89!> \param dac ...
90!> \param rbc ...
91!> \param dbc ...
92!> \param sabc integrals [a|b|c]
93!> \param sdabc derivative [da/dAi|b|c]
94!> \param sabdc derivative [a|b|dc/dCi]
95!> \param int_abc_ext the extremal value of sabc, i.e., MAXVAL(ABS(sabc))
96!> \par History
97!> 05.2014 created (Dorothea Golze)
98!> \author Dorothea Golze
99!> \note overlap3 essentially uses the setup of overlap3_old
100! **************************************************************************************************
101
102 SUBROUTINE overlap3(la_max_set, npgfa, zeta, rpgfa, la_min_set, &
103 lb_max_set, npgfb, zetb, rpgfb, lb_min_set, &
104 lc_max_set, npgfc, zetc, rpgfc, lc_min_set, &
105 rab, dab, rac, dac, rbc, dbc, sabc, &
106 sdabc, sabdc, int_abc_ext)
107
108 INTEGER, INTENT(IN) :: la_max_set, npgfa
109 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
110 INTEGER, INTENT(IN) :: la_min_set, lb_max_set, npgfb
111 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
112 INTEGER, INTENT(IN) :: lb_min_set, lc_max_set, npgfc
113 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc
114 INTEGER, INTENT(IN) :: lc_min_set
115 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
116 REAL(kind=dp), INTENT(IN) :: dab
117 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
118 REAL(kind=dp), INTENT(IN) :: dac
119 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rbc
120 REAL(kind=dp), INTENT(IN) :: dbc
121 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: sabc
122 REAL(kind=dp), DIMENSION(:, :, :, :), &
123 INTENT(INOUT), OPTIONAL :: sdabc, sabdc
124 REAL(dp), INTENT(OUT), OPTIONAL :: int_abc_ext
125
126 CHARACTER(len=*), PARAMETER :: routinen = 'overlap3'
127
128 INTEGER :: ax, ay, az, bx, by, bz, coa, coax, coay, coaz, coc, cocx, cocy, cocz, cx, cy, cz, &
129 handle, i, ipgf, j, jpgf, k, kpgf, l, la, la_max, la_min, la_start, lai, lb, lb_max, &
130 lb_min, lbi, lc, lc_max, lc_min, lci, na, nb, nc, nda, ndc
131 REAL(kind=dp) :: f0, f1, f2, f3, fcx, fcy, fcz, fx, fy, &
132 fz, rcp2, zetg, zetp
133 REAL(kind=dp), DIMENSION(3) :: rag, rbg, rcg, rcp
134 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: s
135 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: sda, sdc
136
137! ---------------------------------------------------------------------------
138
139 CALL timeset(routinen, handle)
140
141 NULLIFY (s, sda, sdc)
142
143 lai = 0
144 lbi = 0
145 lci = 0
146
147 IF (PRESENT(sdabc)) lai = 1
148 IF (PRESENT(sabdc)) lci = 1
149
150 la_max = la_max_set + lai
151 la_min = max(0, la_min_set - lai)
152 lb_max = lb_max_set
153 lb_min = lb_min_set
154 lc_max = lc_max_set + lci
155 lc_min = max(0, lc_min_set - lci)
156
157 ALLOCATE (s(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
158 s = 0._dp
159 IF (PRESENT(sdabc)) THEN
160 ALLOCATE (sda(ncoset(la_max), ncoset(lb_max), ncoset(lc_max), 3))
161 sda = 0._dp
162 END IF
163 IF (PRESENT(sabdc)) THEN
164 ALLOCATE (sdc(ncoset(la_max), ncoset(lb_max), ncoset(lc_max), 3))
165 sdc = 0._dp
166 END IF
167 IF (PRESENT(int_abc_ext)) THEN
168 int_abc_ext = 0.0_dp
169 END IF
170
171! *** Loop over all pairs of primitive Gaussian-type functions ***
172
173 na = 0
174 nda = 0
175 DO ipgf = 1, npgfa
176
177 nb = 0
178 DO jpgf = 1, npgfb
179
180 ! *** Screening ***
181 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
182 nb = nb + ncoset(lb_max_set)
183 cycle
184 END IF
185
186 nc = 0
187 ndc = 0
188 DO kpgf = 1, npgfc
189
190 ! *** Screening ***
191 IF ((rpgfb(jpgf) + rpgfc(kpgf) < dbc) .OR. &
192 (rpgfa(ipgf) + rpgfc(kpgf) < dac)) THEN
193 nc = nc + ncoset(lc_max_set)
194 ndc = ndc + ncoset(lc_max_set)
195 cycle
196 END IF
197
198 ! *** Calculate some prefactors ***
199 zetg = 1.0_dp/(zeta(ipgf) + zetb(jpgf) + zetc(kpgf))
200 zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
201 f0 = (pi*zetg)**1.5_dp
202 f1 = zetb(jpgf)*zetp
203 f2 = 0.5_dp*zetg
204 rcp(:) = f1*rab(:) - rac(:)
205 rcp2 = rcp(1)*rcp(1) + rcp(2)*rcp(2) + rcp(3)*rcp(3)
206
207 ! *** Calculate the basic three-center overlap integral [s|s|s] ***
208 s(1, 1, 1) = f0*exp(-(zeta(ipgf)*f1*dab*dab + zetc(kpgf)*zetg*rcp2/zetp))
209
210! *** Recurrence steps: [s|s|s] -> [a|s|s] ***
211
212 IF (la_max > 0) THEN
213
214! *** Vertical recurrence steps: [s|s|s] -> [a|s|s] ***
215
216 rag(:) = zetg*(zetb(jpgf)*rab(:) + zetc(kpgf)*rac(:))
217
218! *** [p|s|s] = (Gi - Ai)*[s|s|s] (i = x,y,z) ***
219
220 s(2, 1, 1) = rag(1)*s(1, 1, 1)
221 s(3, 1, 1) = rag(2)*s(1, 1, 1)
222 s(4, 1, 1) = rag(3)*s(1, 1, 1)
223
224! *** [a|s|s] = (Gi - Ai)*[a-1i|s|s] + f2*Ni(a-1i)*[a-2i|s|s] ***
225
226 DO la = 2, la_max
227
228! *** Increase the angular momentum component z of function a ***
229
230 s(coset(0, 0, la), 1, 1) = rag(3)*s(coset(0, 0, la - 1), 1, 1) + &
231 f2*real(la - 1, dp)*s(coset(0, 0, la - 2), 1, 1)
232
233! *** Increase the angular momentum component y of function a ***
234
235 az = la - 1
236 s(coset(0, 1, az), 1, 1) = rag(2)*s(coset(0, 0, az), 1, 1)
237
238 DO ay = 2, la
239 az = la - ay
240 s(coset(0, ay, az), 1, 1) = rag(2)*s(coset(0, ay - 1, az), 1, 1) + &
241 f2*real(ay - 1, dp)*s(coset(0, ay - 2, az), 1, 1)
242 END DO
243
244! *** Increase the angular momentum component x of function a ***
245
246 DO ay = 0, la - 1
247 az = la - 1 - ay
248 s(coset(1, ay, az), 1, 1) = rag(1)*s(coset(0, ay, az), 1, 1)
249 END DO
250
251 DO ax = 2, la
252 f3 = f2*real(ax - 1, dp)
253 DO ay = 0, la - ax
254 az = la - ax - ay
255 s(coset(ax, ay, az), 1, 1) = rag(1)*s(coset(ax - 1, ay, az), 1, 1) + &
256 f3*s(coset(ax - 2, ay, az), 1, 1)
257 END DO
258 END DO
259
260 END DO
261
262! *** Recurrence steps: [a|s|s] -> [a|s|b] ***
263
264 IF (lb_max > 0) THEN
265
266! *** Horizontal recurrence steps ***
267
268 rbg(:) = rag(:) - rab(:)
269
270! *** [a|s|p] = [a+1i|s|s] - (Bi - Ai)*[a|s|s] ***
271
272 IF (lb_max == 1) THEN
273 la_start = la_min
274 ELSE
275 la_start = max(0, la_min - 1)
276 END IF
277
278 DO la = la_start, la_max - 1
279 DO ax = 0, la
280 DO ay = 0, la - ax
281 az = la - ax - ay
282 coa = coset(ax, ay, az)
283 coax = coset(ax + 1, ay, az)
284 coay = coset(ax, ay + 1, az)
285 coaz = coset(ax, ay, az + 1)
286 s(coset(ax, ay, az), 2, 1) = s(coax, 1, 1) - rab(1)*s(coa, 1, 1)
287 s(coset(ax, ay, az), 3, 1) = s(coay, 1, 1) - rab(2)*s(coa, 1, 1)
288 s(coset(ax, ay, az), 4, 1) = s(coaz, 1, 1) - rab(3)*s(coa, 1, 1)
289 END DO
290 END DO
291 END DO
292
293! *** Vertical recurrence step ***
294
295! *** [a|s|p] = (Gi - Bi)*[a|s|s] + f2*Ni(a)*[a-1i|s|s] ***
296
297 DO ax = 0, la_max
298 fx = f2*real(ax, dp)
299 DO ay = 0, la_max - ax
300 fy = f2*real(ay, dp)
301 az = la_max - ax - ay
302 fz = f2*real(az, dp)
303 coa = coset(ax, ay, az)
304 IF (ax == 0) THEN
305 s(coa, 2, 1) = rbg(1)*s(coa, 1, 1)
306 ELSE
307 s(coa, 2, 1) = rbg(1)*s(coa, 1, 1) + fx*s(coset(ax - 1, ay, az), 1, 1)
308 END IF
309 IF (ay == 0) THEN
310 s(coa, 3, 1) = rbg(2)*s(coa, 1, 1)
311 ELSE
312 s(coa, 3, 1) = rbg(2)*s(coa, 1, 1) + fy*s(coset(ax, ay - 1, az), 1, 1)
313 END IF
314 IF (az == 0) THEN
315 s(coa, 4, 1) = rbg(3)*s(coa, 1, 1)
316 ELSE
317 s(coa, 4, 1) = rbg(3)*s(coa, 1, 1) + fz*s(coset(ax, ay, az - 1), 1, 1)
318 END IF
319 END DO
320 END DO
321
322! *** Recurrence steps: [a|s|p] -> [a|s|b] ***
323
324 DO lb = 2, lb_max
325
326! *** Horizontal recurrence steps ***
327
328! *** [a|s|b] = [a+1i|s|b-1i] - (Bi - Ai)*[a|s|b-1i] ***
329
330 IF (lb == lb_max) THEN
331 la_start = la_min
332 ELSE
333 la_start = max(0, la_min - 1)
334 END IF
335
336 DO la = la_start, la_max - 1
337 DO ax = 0, la
338 DO ay = 0, la - ax
339 az = la - ax - ay
340
341 coa = coset(ax, ay, az)
342 coax = coset(ax + 1, ay, az)
343 coay = coset(ax, ay + 1, az)
344 coaz = coset(ax, ay, az + 1)
345
346! *** Shift of angular momentum component z from a to b ***
347
348 s(coa, coset(0, 0, lb), 1) = &
349 s(coaz, coset(0, 0, lb - 1), 1) - &
350 rab(3)*s(coa, coset(0, 0, lb - 1), 1)
351
352! *** Shift of angular momentum component y from a to b ***
353
354 DO by = 1, lb
355 bz = lb - by
356 s(coa, coset(0, by, bz), 1) = &
357 s(coay, coset(0, by - 1, bz), 1) - &
358 rab(2)*s(coa, coset(0, by - 1, bz), 1)
359 END DO
360
361! *** Shift of angular momentum component x from a to b ***
362
363 DO bx = 1, lb
364 DO by = 0, lb - bx
365 bz = lb - bx - by
366 s(coa, coset(bx, by, bz), 1) = &
367 s(coax, coset(bx - 1, by, bz), 1) - &
368 rab(1)*s(coa, coset(bx - 1, by, bz), 1)
369 END DO
370 END DO
371
372 END DO
373 END DO
374 END DO
375
376! *** Vertical recurrence step ***
377
378! *** [a|s|b] = (Gi - Bi)*[a|s|b-1i] + ***
379! *** f2*Ni(a)*[a-1i|s|b-1i] + ***
380! *** f2*Ni(b-1i)*[a|s|b-2i] ***
381
382 DO ax = 0, la_max
383 fx = f2*real(ax, dp)
384 DO ay = 0, la_max - ax
385 fy = f2*real(ay, dp)
386 az = la_max - ax - ay
387 fz = f2*real(az, dp)
388
389 coa = coset(ax, ay, az)
390
391 f3 = f2*real(lb - 1, dp)
392
393! *** Shift of angular momentum component z from a to b ***
394
395 IF (az == 0) THEN
396 s(coa, coset(0, 0, lb), 1) = &
397 rbg(3)*s(coa, coset(0, 0, lb - 1), 1) + &
398 f3*s(coa, coset(0, 0, lb - 2), 1)
399 ELSE
400 coaz = coset(ax, ay, az - 1)
401 s(coa, coset(0, 0, lb), 1) = &
402 rbg(3)*s(coa, coset(0, 0, lb - 1), 1) + &
403 fz*s(coaz, coset(0, 0, lb - 1), 1) + &
404 f3*s(coa, coset(0, 0, lb - 2), 1)
405 END IF
406
407! *** Shift of angular momentum component y from a to b ***
408
409 IF (ay == 0) THEN
410 bz = lb - 1
411 s(coa, coset(0, 1, bz), 1) = &
412 rbg(2)*s(coa, coset(0, 0, bz), 1)
413 DO by = 2, lb
414 bz = lb - by
415 f3 = f2*real(by - 1, dp)
416 s(coa, coset(0, by, bz), 1) = &
417 rbg(2)*s(coa, coset(0, by - 1, bz), 1) + &
418 f3*s(coa, coset(0, by - 2, bz), 1)
419 END DO
420 ELSE
421 coay = coset(ax, ay - 1, az)
422 bz = lb - 1
423 s(coa, coset(0, 1, bz), 1) = &
424 rbg(2)*s(coa, coset(0, 0, bz), 1) + &
425 fy*s(coay, coset(0, 0, bz), 1)
426 DO by = 2, lb
427 bz = lb - by
428 f3 = f2*real(by - 1, dp)
429 s(coa, coset(0, by, bz), 1) = &
430 rbg(2)*s(coa, coset(0, by - 1, bz), 1) + &
431 fy*s(coay, coset(0, by - 1, bz), 1) + &
432 f3*s(coa, coset(0, by - 2, bz), 1)
433 END DO
434 END IF
435
436! *** Shift of angular momentum component x from a to b ***
437
438 IF (ax == 0) THEN
439 DO by = 0, lb - 1
440 bz = lb - 1 - by
441 s(coa, coset(1, by, bz), 1) = &
442 rbg(1)*s(coa, coset(0, by, bz), 1)
443 END DO
444 DO bx = 2, lb
445 f3 = f2*real(bx - 1, dp)
446 DO by = 0, lb - bx
447 bz = lb - bx - by
448 s(coa, coset(bx, by, bz), 1) = &
449 rbg(1)*s(coa, coset(bx - 1, by, bz), 1) + &
450 f3*s(coa, coset(bx - 2, by, bz), 1)
451 END DO
452 END DO
453 ELSE
454 coax = coset(ax - 1, ay, az)
455 DO by = 0, lb - 1
456 bz = lb - 1 - by
457 s(coa, coset(1, by, bz), 1) = &
458 rbg(1)*s(coa, coset(0, by, bz), 1) + &
459 fx*s(coax, coset(0, by, bz), 1)
460 END DO
461 DO bx = 2, lb
462 f3 = f2*real(bx - 1, dp)
463 DO by = 0, lb - bx
464 bz = lb - bx - by
465 s(coa, coset(bx, by, bz), 1) = &
466 rbg(1)*s(coa, coset(bx - 1, by, bz), 1) + &
467 fx*s(coax, coset(bx - 1, by, bz), 1) + &
468 f3*s(coa, coset(bx - 2, by, bz), 1)
469 END DO
470 END DO
471 END IF
472
473 END DO
474 END DO
475
476 END DO
477
478 END IF
479
480 ELSE
481
482 IF (lb_max > 0) THEN
483
484! *** Vertical recurrence steps: [s|s|s] -> [s|s|b] ***
485
486 rbg(:) = -zetg*(zeta(ipgf)*rab(:) - zetc(kpgf)*rbc(:))
487
488! *** [s|s|p] = (Gi - Bi)*[s|s|s] ***
489
490 s(1, 2, 1) = rbg(1)*s(1, 1, 1)
491 s(1, 3, 1) = rbg(2)*s(1, 1, 1)
492 s(1, 4, 1) = rbg(3)*s(1, 1, 1)
493
494! *** [s|s|b] = (Gi - Bi)*[s|s|b-1i] + f2*Ni(b-1i)*[s|s|b-2i] ***
495
496 DO lb = 2, lb_max
497
498! *** Increase the angular momentum component z of function b ***
499
500 s(1, coset(0, 0, lb), 1) = rbg(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 function b ***
504
505 bz = lb - 1
506 s(1, coset(0, 1, bz), 1) = rbg(2)*s(1, coset(0, 0, bz), 1)
507
508 DO by = 2, lb
509 bz = lb - by
510 s(1, coset(0, by, bz), 1) = &
511 rbg(2)*s(1, coset(0, by - 1, bz), 1) + &
512 f2*real(by - 1, dp)*s(1, coset(0, by - 2, bz), 1)
513 END DO
514
515! *** Increase the angular momentum component x of function b ***
516
517 DO by = 0, lb - 1
518 bz = lb - 1 - by
519 s(1, coset(1, by, bz), 1) = rbg(1)*s(1, coset(0, by, bz), 1)
520 END DO
521
522 DO bx = 2, lb
523 f3 = f2*real(bx - 1, dp)
524 DO by = 0, lb - bx
525 bz = lb - bx - by
526 s(1, coset(bx, by, bz), 1) = rbg(1)*s(1, coset(bx - 1, by, bz), 1) + &
527 f3*s(1, coset(bx - 2, by, bz), 1)
528 END DO
529 END DO
530
531 END DO
532
533 END IF
534
535 END IF
536
537! *** Recurrence steps: [a|s|b] -> [a|c|b] ***
538
539 IF (lc_max > 0) THEN
540
541! *** Vertical recurrence steps: [s|s|s] -> [s|c|s] ***
542
543 rcg(:) = -zetg*(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))
544
545! *** [s|p|s] = (Gi - Ci)*[s|s|s] (i = x,y,z) ***
546
547 s(1, 1, 2) = rcg(1)*s(1, 1, 1)
548 s(1, 1, 3) = rcg(2)*s(1, 1, 1)
549 s(1, 1, 4) = rcg(3)*s(1, 1, 1)
550
551! *** [s|c|s] = (Gi - Ci)*[s|c-1i|s] + f2*Ni(c-1i)*[s|c-2i|s] ***
552
553 DO lc = 2, lc_max
554
555! *** Increase the angular momentum component z of function c ***
556
557 s(1, 1, coset(0, 0, lc)) = rcg(3)*s(1, 1, coset(0, 0, lc - 1)) + &
558 f2*real(lc - 1, dp)*s(1, 1, coset(0, 0, lc - 2))
559
560! *** Increase the angular momentum component y of function c ***
561
562 cz = lc - 1
563 s(1, 1, coset(0, 1, cz)) = rcg(2)*s(1, 1, coset(0, 0, cz))
564
565 DO cy = 2, lc
566 cz = lc - cy
567 s(1, 1, coset(0, cy, cz)) = rcg(2)*s(1, 1, coset(0, cy - 1, cz)) + &
568 f2*real(cy - 1, dp)*s(1, 1, coset(0, cy - 2, cz))
569 END DO
570
571! *** Increase the angular momentum component x of function c ***
572
573 DO cy = 0, lc - 1
574 cz = lc - 1 - cy
575 s(1, 1, coset(1, cy, cz)) = rcg(1)*s(1, 1, coset(0, cy, cz))
576 END DO
577
578 DO cx = 2, lc
579 f3 = f2*real(cx - 1, dp)
580 DO cy = 0, lc - cx
581 cz = lc - cx - cy
582 s(1, 1, coset(cx, cy, cz)) = rcg(1)*s(1, 1, coset(cx - 1, cy, cz)) + &
583 f3*s(1, 1, coset(cx - 2, cy, cz))
584 END DO
585 END DO
586
587 END DO
588
589! *** Recurrence steps: [s|c|s] -> [a|c|b] ***
590
591 DO lc = 1, lc_max
592
593 DO cx = 0, lc
594 DO cy = 0, lc - cx
595 cz = lc - cx - cy
596
597 coc = coset(cx, cy, cz)
598 cocx = coset(max(0, cx - 1), cy, cz)
599 cocy = coset(cx, max(0, cy - 1), cz)
600 cocz = coset(cx, cy, max(0, cz - 1))
601
602 fcx = f2*real(cx, dp)
603 fcy = f2*real(cy, dp)
604 fcz = f2*real(cz, dp)
605
606! *** Recurrence steps: [s|c|s] -> [a|c|s] ***
607
608 IF (la_max > 0) THEN
609
610! *** Vertical recurrence steps: [s|c|s] -> [a|c|s] ***
611
612 rag(:) = rcg(:) + rac(:)
613
614! *** [p|c|s] = (Gi - Ai)*[s|c|s] + f2*Ni(c)*[s|c-1i|s] ***
615
616 s(2, 1, coc) = rag(1)*s(1, 1, coc) + fcx*s(1, 1, cocx)
617 s(3, 1, coc) = rag(2)*s(1, 1, coc) + fcy*s(1, 1, cocy)
618 s(4, 1, coc) = rag(3)*s(1, 1, coc) + fcz*s(1, 1, cocz)
619
620! *** [a|c|s] = (Gi - Ai)*[a-1i|c|s] + ***
621! *** f2*Ni(a-1i)*[a-2i|c|s] + ***
622! *** f2*Ni(c)*[a-1i|c-1i|s] ***
623
624 DO la = 2, la_max
625
626! *** Increase the angular momentum component z of a ***
627
628 s(coset(0, 0, la), 1, coc) = &
629 rag(3)*s(coset(0, 0, la - 1), 1, coc) + &
630 f2*real(la - 1, dp)*s(coset(0, 0, la - 2), 1, coc) + &
631 fcz*s(coset(0, 0, la - 1), 1, cocz)
632
633! *** Increase the angular momentum component y of a ***
634
635 az = la - 1
636 s(coset(0, 1, az), 1, coc) = &
637 rag(2)*s(coset(0, 0, az), 1, coc) + &
638 fcy*s(coset(0, 0, az), 1, cocy)
639
640 DO ay = 2, la
641 az = la - ay
642 s(coset(0, ay, az), 1, coc) = &
643 rag(2)*s(coset(0, ay - 1, az), 1, coc) + &
644 f2*real(ay - 1, dp)*s(coset(0, ay - 2, az), 1, coc) + &
645 fcy*s(coset(0, ay - 1, az), 1, cocy)
646 END DO
647
648! *** Increase the angular momentum component x of a ***
649
650 DO ay = 0, la - 1
651 az = la - 1 - ay
652 s(coset(1, ay, az), 1, coc) = &
653 rag(1)*s(coset(0, ay, az), 1, coc) + &
654 fcx*s(coset(0, ay, az), 1, cocx)
655 END DO
656
657 DO ax = 2, la
658 f3 = f2*real(ax - 1, dp)
659 DO ay = 0, la - ax
660 az = la - ax - ay
661 s(coset(ax, ay, az), 1, coc) = &
662 rag(1)*s(coset(ax - 1, ay, az), 1, coc) + &
663 f3*s(coset(ax - 2, ay, az), 1, coc) + &
664 fcx*s(coset(ax - 1, ay, az), 1, cocx)
665 END DO
666 END DO
667
668 END DO
669
670! *** Recurrence steps: [a|c|s] -> [a|c|b] ***
671
672 IF (lb_max > 0) THEN
673
674! *** Horizontal recurrence steps ***
675
676 rbg(:) = rag(:) - rab(:)
677
678! *** [a|c|p] = [a+1i|c|s] - (Bi - Ai)*[a|c|s] ***
679
680 IF (lb_max == 1) THEN
681 la_start = la_min
682 ELSE
683 la_start = max(0, la_min - 1)
684 END IF
685
686 DO la = la_start, la_max - 1
687 DO ax = 0, la
688 DO ay = 0, la - ax
689 az = la - ax - ay
690 coa = coset(ax, ay, az)
691 coax = coset(ax + 1, ay, az)
692 coay = coset(ax, ay + 1, az)
693 coaz = coset(ax, ay, az + 1)
694 s(coa, 2, coc) = s(coax, 1, coc) - rab(1)*s(coa, 1, coc)
695 s(coa, 3, coc) = s(coay, 1, coc) - rab(2)*s(coa, 1, coc)
696 s(coa, 4, coc) = s(coaz, 1, coc) - rab(3)*s(coa, 1, coc)
697 END DO
698 END DO
699 END DO
700
701! *** Vertical recurrence step ***
702
703! *** [a|c|p] = (Gi - Bi)*[a|c|s] + ***
704! f2*Ni(a)*[a-1i|c|s] + ***
705! f2*Ni(c)*[a|c-1i|s] ***
706
707 DO ax = 0, la_max
708 fx = f2*real(ax, dp)
709 DO ay = 0, la_max - ax
710 fy = f2*real(ay, dp)
711 az = la_max - ax - ay
712 fz = f2*real(az, dp)
713 coa = coset(ax, ay, az)
714 IF (ax == 0) THEN
715 s(coa, 2, coc) = rbg(1)*s(coa, 1, coc) + &
716 fcx*s(coa, 1, cocx)
717 ELSE
718 s(coa, 2, coc) = rbg(1)*s(coa, 1, coc) + &
719 fx*s(coset(ax - 1, ay, az), 1, coc) + &
720 fcx*s(coa, 1, cocx)
721 END IF
722 IF (ay == 0) THEN
723 s(coa, 3, coc) = rbg(2)*s(coa, 1, coc) + &
724 fcy*s(coa, 1, cocy)
725 ELSE
726 s(coa, 3, coc) = rbg(2)*s(coa, 1, coc) + &
727 fy*s(coset(ax, ay - 1, az), 1, coc) + &
728 fcy*s(coa, 1, cocy)
729 END IF
730 IF (az == 0) THEN
731 s(coa, 4, coc) = rbg(3)*s(coa, 1, coc) + &
732 fcz*s(coa, 1, cocz)
733 ELSE
734 s(coa, 4, coc) = rbg(3)*s(coa, 1, coc) + &
735 fz*s(coset(ax, ay, az - 1), 1, coc) + &
736 fcz*s(coa, 1, cocz)
737 END IF
738 END DO
739 END DO
740
741! *** Recurrence steps: [a|c|p] -> [a|c|b] ***
742
743 DO lb = 2, lb_max
744
745! *** Horizontal recurrence steps ***
746
747! *** [a|c|b] = [a+1i|c|b-1i] - (Bi - Ai)*[a|c|b-1i] ***
748
749 IF (lb == lb_max) THEN
750 la_start = la_min
751 ELSE
752 la_start = max(0, la_min - 1)
753 END IF
754
755 DO la = la_start, la_max - 1
756 DO ax = 0, la
757 DO ay = 0, la - ax
758 az = la - ax - ay
759
760 coa = coset(ax, ay, az)
761 coax = coset(ax + 1, ay, az)
762 coay = coset(ax, ay + 1, az)
763 coaz = coset(ax, ay, az + 1)
764
765! *** Shift of angular momentum ***
766! *** component z from a to b ***
767
768 s(coa, coset(0, 0, lb), coc) = &
769 s(coaz, coset(0, 0, lb - 1), coc) - &
770 rab(3)*s(coa, coset(0, 0, lb - 1), coc)
771
772! *** Shift of angular momentum ***
773! *** component y from a to b ***
774
775 DO by = 1, lb
776 bz = lb - by
777 s(coa, coset(0, by, bz), coc) = &
778 s(coay, coset(0, by - 1, bz), coc) - &
779 rab(2)*s(coa, coset(0, by - 1, bz), coc)
780 END DO
781
782! *** Shift of angular momentum ***
783! *** component x from a to b ***
784
785 DO bx = 1, lb
786 DO by = 0, lb - bx
787 bz = lb - bx - by
788 s(coa, coset(bx, by, bz), coc) = &
789 s(coax, coset(bx - 1, by, bz), coc) - &
790 rab(1)*s(coa, coset(bx - 1, by, bz), coc)
791 END DO
792 END DO
793
794 END DO
795 END DO
796 END DO
797
798! *** Vertical recurrence step ***
799
800! *** [a|c|b] = (Gi - Bi)*[a|c|b-1i] + ***
801! *** f2*Ni(a)*[a-1i|c|b-1i] + ***
802! *** f2*Ni(b-1i)*[a|c|b-2i] + ***
803! *** f2*Ni(c)*[a|c-1i|b-1i] ***
804
805 DO ax = 0, la_max
806 fx = f2*real(ax, dp)
807 DO ay = 0, la_max - ax
808 fy = f2*real(ay, dp)
809 az = la_max - ax - ay
810 fz = f2*real(az, dp)
811
812 coa = coset(ax, ay, az)
813 coax = coset(max(0, ax - 1), ay, az)
814 coay = coset(ax, max(0, ay - 1), az)
815 coaz = coset(ax, ay, max(0, az - 1))
816
817 f3 = f2*real(lb - 1, dp)
818
819! *** Shift of angular momentum ***
820! *** component z from a to b ***
821
822 IF (az == 0) THEN
823 s(coa, coset(0, 0, lb), coc) = &
824 rbg(3)*s(coa, coset(0, 0, lb - 1), coc) + &
825 f3*s(coa, coset(0, 0, lb - 2), coc) + &
826 fcz*s(coa, coset(0, 0, lb - 1), cocz)
827 ELSE
828 s(coa, coset(0, 0, lb), coc) = &
829 rbg(3)*s(coa, coset(0, 0, lb - 1), coc) + &
830 fz*s(coaz, coset(0, 0, lb - 1), coc) + &
831 f3*s(coa, coset(0, 0, lb - 2), coc) + &
832 fcz*s(coa, coset(0, 0, lb - 1), cocz)
833 END IF
834
835! *** Shift of angular momentum ***
836! *** component y from a to b ***
837
838 IF (ay == 0) THEN
839 bz = lb - 1
840 s(coa, coset(0, 1, bz), coc) = &
841 rbg(2)*s(coa, coset(0, 0, bz), coc) + &
842 fcy*s(coa, coset(0, 0, bz), cocy)
843 DO by = 2, lb
844 bz = lb - by
845 f3 = f2*real(by - 1, dp)
846 s(coa, coset(0, by, bz), coc) = &
847 rbg(2)*s(coa, coset(0, by - 1, bz), coc) + &
848 f3*s(coa, coset(0, by - 2, bz), coc) + &
849 fcy*s(coa, coset(0, by - 1, bz), cocy)
850 END DO
851 ELSE
852 bz = lb - 1
853 s(coa, coset(0, 1, bz), coc) = &
854 rbg(2)*s(coa, coset(0, 0, bz), coc) + &
855 fy*s(coay, coset(0, 0, bz), coc) + &
856 fcy*s(coa, coset(0, 0, bz), cocy)
857 DO by = 2, lb
858 bz = lb - by
859 f3 = f2*real(by - 1, dp)
860 s(coa, coset(0, by, bz), coc) = &
861 rbg(2)*s(coa, coset(0, by - 1, bz), coc) + &
862 fy*s(coay, coset(0, by - 1, bz), coc) + &
863 f3*s(coa, coset(0, by - 2, bz), coc) + &
864 fcy*s(coa, coset(0, by - 1, bz), cocy)
865 END DO
866 END IF
867
868! *** Shift of angular momentum ***
869! *** component x from a to b ***
870
871 IF (ax == 0) THEN
872 DO by = 0, lb - 1
873 bz = lb - 1 - by
874 s(coa, coset(1, by, bz), coc) = &
875 rbg(1)*s(coa, coset(0, by, bz), coc) + &
876 fcx*s(coa, coset(0, by, bz), cocx)
877 END DO
878 DO bx = 2, lb
879 f3 = f2*real(bx - 1, dp)
880 DO by = 0, lb - bx
881 bz = lb - bx - by
882 s(coa, coset(bx, by, bz), coc) = &
883 rbg(1)*s(coa, coset(bx - 1, by, bz), coc) + &
884 f3*s(coa, coset(bx - 2, by, bz), coc) + &
885 fcx*s(coa, coset(bx - 1, by, bz), cocx)
886 END DO
887 END DO
888 ELSE
889 DO by = 0, lb - 1
890 bz = lb - 1 - by
891 s(coa, coset(1, by, bz), coc) = &
892 rbg(1)*s(coa, coset(0, by, bz), coc) + &
893 fx*s(coax, coset(0, by, bz), coc) + &
894 fcx*s(coa, coset(0, by, bz), cocx)
895 END DO
896 DO bx = 2, lb
897 f3 = f2*real(bx - 1, dp)
898 DO by = 0, lb - bx
899 bz = lb - bx - by
900 s(coa, coset(bx, by, bz), coc) = &
901 rbg(1)*s(coa, coset(bx - 1, by, bz), coc) + &
902 fx*s(coax, coset(bx - 1, by, bz), coc) + &
903 f3*s(coa, coset(bx - 2, by, bz), coc) + &
904 fcx*s(coa, coset(bx - 1, by, bz), cocx)
905 END DO
906 END DO
907 END IF
908
909 END DO
910 END DO
911
912 END DO
913
914 END IF
915
916 ELSE
917
918 IF (lb_max > 0) THEN
919
920! *** Vertical recurrence steps: [s|c|s] -> [s|c|b] ***
921
922 rbg(:) = rcg(:) + rbc(:)
923
924! *** [s|c|p] = (Gi - Bi)*[s|c|s] + f2*Ni(c)*[s|c-1i|s] ***
925
926 s(1, 2, coc) = rbg(1)*s(1, 1, coc) + fcx*s(1, 1, cocx)
927 s(1, 3, coc) = rbg(2)*s(1, 1, coc) + fcy*s(1, 1, cocy)
928 s(1, 4, coc) = rbg(3)*s(1, 1, coc) + fcz*s(1, 1, cocz)
929
930! *** [s|c|b] = (Gi - Bi)*[s|c|b-1i] + ***
931! *** f2*Ni(b-1i)*[s|c|b-2i] ***
932! *** f2*Ni(c)*[s|c-1i|b-1i] ***
933
934 DO lb = 2, lb_max
935
936! *** Increase the angular momentum component z of b ***
937
938 s(1, coset(0, 0, lb), coc) = &
939 rbg(3)*s(1, coset(0, 0, lb - 1), coc) + &
940 f2*real(lb - 1, dp)*s(1, coset(0, 0, lb - 2), coc) + &
941 fcz*s(1, coset(0, 0, lb - 1), cocz)
942
943! *** Increase the angular momentum component y of b ***
944
945 bz = lb - 1
946 s(1, coset(0, 1, bz), coc) = &
947 rbg(2)*s(1, coset(0, 0, bz), coc) + &
948 fcy*s(1, coset(0, 0, bz), cocy)
949
950 DO by = 2, lb
951 bz = lb - by
952 s(1, coset(0, by, bz), coc) = &
953 rbg(2)*s(1, coset(0, by - 1, bz), coc) + &
954 f2*real(by - 1, dp)*s(1, coset(0, by - 2, bz), coc) + &
955 fcy*s(1, coset(0, by - 1, bz), cocy)
956 END DO
957
958! *** Increase the angular momentum component x of b ***
959
960 DO by = 0, lb - 1
961 bz = lb - 1 - by
962 s(1, coset(1, by, bz), coc) = &
963 rbg(1)*s(1, coset(0, by, bz), coc) + &
964 fcx*s(1, coset(0, by, bz), cocx)
965 END DO
966
967 DO bx = 2, lb
968 f3 = f2*real(bx - 1, dp)
969 DO by = 0, lb - bx
970 bz = lb - bx - by
971 s(1, coset(bx, by, bz), coc) = &
972 rbg(1)*s(1, coset(bx - 1, by, bz), coc) + &
973 f3*s(1, coset(bx - 2, by, bz), coc) + &
974 fcx*s(1, coset(bx - 1, by, bz), cocx)
975 END DO
976 END DO
977
978 END DO
979
980 END IF
981
982 END IF
983
984 END DO
985 END DO
986
987 END DO
988
989 END IF
990
991! *** Store integrals
992
993 IF (PRESENT(int_abc_ext)) THEN
994 DO k = ncoset(lc_min_set - 1) + 1, ncoset(lc_max_set)
995 DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
996 DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
997 sabc(na + i, nb + j, nc + k) = s(i, j, k)
998 int_abc_ext = max(int_abc_ext, abs(s(i, j, k)))
999 END DO
1000 END DO
1001 END DO
1002 ELSE
1003 DO k = ncoset(lc_min_set - 1) + 1, ncoset(lc_max_set)
1004 DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
1005 DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
1006 sabc(na + i, nb + j, nc + k) = s(i, j, k)
1007 END DO
1008 END DO
1009 END DO
1010 END IF
1011
1012! *** Calculate the requested derivatives with respect to ***
1013! *** the nuclear coordinates of the atomic center a and c ***
1014
1015 IF (PRESENT(sdabc) .OR. PRESENT(sabdc)) THEN
1016 CALL derivatives_overlap3(la_max_set, la_min_set, lb_max_set, lb_min_set, &
1017 lc_max_set, lc_min_set, zeta(ipgf), zetc(kpgf), &
1018 s, sda, sdc)
1019 END IF
1020
1021! *** Store the first derivatives of the primitive overlap integrals ***
1022
1023 IF (PRESENT(sdabc)) THEN
1024 DO k = 1, 3
1025 DO l = 1, ncoset(lc_max_set)
1026 DO j = 1, ncoset(lb_max_set)
1027 DO i = 1, ncoset(la_max_set)
1028 sdabc(nda + i, nb + j, nc + l, k) = sda(i, j, l, k)
1029 END DO
1030 END DO
1031 END DO
1032 END DO
1033 END IF
1034
1035 IF (PRESENT(sabdc)) THEN
1036 DO k = 1, 3
1037 DO l = 1, ncoset(lc_max_set)
1038 DO j = 1, ncoset(lb_max_set)
1039 DO i = 1, ncoset(la_max_set)
1040 sabdc(na + i, nb + j, ndc + l, k) = sdc(i, j, l, k)
1041 END DO
1042 END DO
1043 END DO
1044 END DO
1045 END IF
1046
1047 nc = nc + ncoset(lc_max_set)
1048 ndc = ndc + ncoset(lc_max_set)
1049 END DO
1050
1051 nb = nb + ncoset(lb_max)
1052 END DO
1053
1054 na = na + ncoset(la_max_set)
1055 nda = nda + ncoset(la_max_set)
1056 END DO
1057
1058 DEALLOCATE (s)
1059 IF (PRESENT(sdabc)) THEN
1060 DEALLOCATE (sda)
1061 END IF
1062 IF (PRESENT(sabdc)) THEN
1063 DEALLOCATE (sdc)
1064 END IF
1065
1066 CALL timestop(handle)
1067
1068 END SUBROUTINE overlap3
1069
1070! **************************************************************************************************
1071!> \brief Calculates the derivatives of the three-center overlap integral [a|b|c]
1072!> with respect to the nuclear coordinates of the atomic center a and c
1073!> \param la_max_set ...
1074!> \param la_min_set ...
1075!> \param lb_max_set ...
1076!> \param lb_min_set ...
1077!> \param lc_max_set ...
1078!> \param lc_min_set ...
1079!> \param zeta ...
1080!> \param zetc ...
1081!> \param s integrals [a|b|c]
1082!> \param sda derivative [da/dAi|b|c]
1083!> \param sdc derivative [a|b|dc/dCi]
1084! **************************************************************************************************
1085 SUBROUTINE derivatives_overlap3(la_max_set, la_min_set, lb_max_set, lb_min_set, &
1086 lc_max_set, lc_min_set, zeta, zetc, s, sda, sdc)
1087
1088 INTEGER, INTENT(IN) :: la_max_set, la_min_set, lb_max_set, &
1089 lb_min_set, lc_max_set, lc_min_set
1090 REAL(kind=dp), INTENT(IN) :: zeta, zetc
1091 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: s
1092 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: sda, sdc
1093
1094 CHARACTER(len=*), PARAMETER :: routinen = 'derivatives_overlap3'
1095
1096 INTEGER :: ax, ay, az, bx, by, bz, coa, coamx, coamy, coamz, coapx, coapy, coapz, cob, coc, &
1097 cocmx, cocmy, cocmz, cocpx, cocpy, cocpz, cx, cy, cz, devx, devy, devz, handle, la, lb, lc
1098 REAL(kind=dp) :: fax, fay, faz, fcx, fcy, fcz, fexpa, &
1099 fexpc
1100
1101 CALL timeset(routinen, handle)
1102
1103 fexpa = 2.0_dp*zeta
1104 fexpc = 2.0_dp*zetc
1105
1106! derivative with respec to x,y,z
1107
1108 devx = 1
1109 devy = 2
1110 devz = 3
1111
1112! *** [da/dAi|b|c] = 2*zeta*[a+1i|b|c] - Ni(a)[a-1i|b|c] ***
1113! *** [a|b|dc/dCi] = 2*zetc*[a|b|c+1i] - Ni(c)[a|b|c-1i] ***
1114
1115 DO la = la_min_set, la_max_set
1116 DO ax = 0, la
1117 fax = real(ax, dp)
1118 DO ay = 0, la - ax
1119 fay = real(ay, dp)
1120 az = la - ax - ay
1121 faz = real(az, dp)
1122 coa = coset(ax, ay, az)
1123 coamx = coset(ax - 1, ay, az)
1124 coamy = coset(ax, ay - 1, az)
1125 coamz = coset(ax, ay, az - 1)
1126 coapx = coset(ax + 1, ay, az)
1127 coapy = coset(ax, ay + 1, az)
1128 coapz = coset(ax, ay, az + 1)
1129 DO lb = lb_min_set, lb_max_set
1130 DO bx = 0, lb
1131 DO by = 0, lb - bx
1132 bz = lb - bx - by
1133 cob = coset(bx, by, bz)
1134 DO lc = lc_min_set, lc_max_set
1135 DO cx = 0, lc
1136 fcx = real(cx, dp)
1137 DO cy = 0, lc - cx
1138 fcy = real(cy, dp)
1139 cz = lc - cx - cy
1140 fcz = real(cz, dp)
1141 coc = coset(cx, cy, cz)
1142 cocmx = coset(cx - 1, cy, cz)
1143 cocmy = coset(cx, cy - 1, cz)
1144 cocmz = coset(cx, cy, cz - 1)
1145 cocpx = coset(cx + 1, cy, cz)
1146 cocpy = coset(cx, cy + 1, cz)
1147 cocpz = coset(cx, cy, cz + 1)
1148 IF (ASSOCIATED(sda)) THEN
1149 sda(coa, cob, coc, devx) = fexpa*s(coapx, cob, coc) - &
1150 fax*s(coamx, cob, coc)
1151 sda(coa, cob, coc, devy) = fexpa*s(coapy, cob, coc) - &
1152 fay*s(coamy, cob, coc)
1153 sda(coa, cob, coc, devz) = fexpa*s(coapz, cob, coc) - &
1154 faz*s(coamz, cob, coc)
1155 END IF
1156 IF (ASSOCIATED(sdc)) THEN
1157 sdc(coa, cob, coc, devx) = fexpc*s(coa, cob, cocpx) - &
1158 fcx*s(coa, cob, cocmx)
1159 sdc(coa, cob, coc, devy) = fexpc*s(coa, cob, cocpy) - &
1160 fcy*s(coa, cob, cocmy)
1161 sdc(coa, cob, coc, devz) = fexpc*s(coa, cob, cocpz) - &
1162 fcz*s(coa, cob, cocmz)
1163 END IF
1164 END DO
1165 END DO
1166 END DO
1167 END DO
1168 END DO
1169 END DO
1170 END DO
1171 END DO
1172 END DO
1173
1174 CALL timestop(handle)
1175
1176 END SUBROUTINE derivatives_overlap3
1177
1178END MODULE ai_overlap3
subroutine, public overlap3(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max_set, npgfb, zetb, rpgfb, lb_min_set, lc_max_set, npgfc, zetc, rpgfc, lc_min_set, rab, dab, rac, dac, rbc, dbc, sabc, sdabc, sabdc, int_abc_ext)
Calculation of three-center overlap integrals [a|b|c] over primitive Cartesian Gaussian functions.
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
integer, dimension(:, :, :), allocatable, public coset