Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of 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 Parameters
14 : !> - ax,ay,az : Angular momentum index numbers of orbital a.
15 : !> - bx,by,bz : Angular momentum index numbers of orbital b.
16 : !> - coset : Cartesian orbital set pointer.
17 : !> - dab : Distance between the atomic centers a and b.
18 : !> - l{a,b} : Angular momentum quantum number of shell a or b.
19 : !> - l{a,b}_max: Maximum angular momentum quantum number of shell a or b.
20 : !> - l{a,b}_min: Minimum angular momentum quantum number of shell a or b.
21 : !> - rab : Distance vector between the atomic centers a and b.
22 : !> - rpgf{a,b} : Radius of the primitive Gaussian-type function a or b.
23 : !> - sab : Shell set of overlap integrals.
24 : !> - zet{a,b} : Exponents of the Gaussian-type functions a or b.
25 : !> - zetp : Reciprocal of the sum of the exponents of orbital a and b.
26 : ! **************************************************************************************************
27 : MODULE ai_overlap_aabb
28 :
29 : USE ai_os_rr, ONLY: os_rr_ovlp
30 : USE kinds, ONLY: dp
31 : USE mathconstants, ONLY: pi
32 : USE orbital_pointers, ONLY: coset,&
33 : indco,&
34 : ncoset
35 : #include "../base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 :
39 : PRIVATE
40 :
41 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap_aabb'
42 :
43 : ! *** Public subroutines ***
44 : PUBLIC :: overlap_aabb
45 :
46 : CONTAINS
47 :
48 : ! **************************************************************************************************
49 : !> \brief Purpose: Calculation of the two-center overlap integrals [aa|bb]
50 : !> over Cartesian Gaussian-type functions.
51 : !> \param la_max_set1 ...
52 : !> \param la_min_set1 ...
53 : !> \param npgfa1 ...
54 : !> \param rpgfa1 ...
55 : !> \param zeta1 ...
56 : !> \param la_max_set2 ...
57 : !> \param la_min_set2 ...
58 : !> \param npgfa2 ...
59 : !> \param rpgfa2 ...
60 : !> \param zeta2 ...
61 : !> \param lb_max_set1 ...
62 : !> \param lb_min_set1 ...
63 : !> \param npgfb1 ...
64 : !> \param rpgfb1 ...
65 : !> \param zetb1 ...
66 : !> \param lb_max_set2 ...
67 : !> \param lb_min_set2 ...
68 : !> \param npgfb2 ...
69 : !> \param rpgfb2 ...
70 : !> \param zetb2 ...
71 : !> \param asets_equal ...
72 : !> \param bsets_equal ...
73 : !> \param rab ...
74 : !> \param dab ...
75 : !> \param saabb primitive integrals, including zero screened blocks and lower-shell padding
76 : !> \date 06.2014
77 : !> \author Dorothea Golze
78 : ! **************************************************************************************************
79 45 : SUBROUTINE overlap_aabb(la_max_set1, la_min_set1, npgfa1, rpgfa1, zeta1, &
80 9 : la_max_set2, la_min_set2, npgfa2, rpgfa2, zeta2, &
81 9 : lb_max_set1, lb_min_set1, npgfb1, rpgfb1, zetb1, &
82 9 : lb_max_set2, lb_min_set2, npgfb2, rpgfb2, zetb2, &
83 9 : asets_equal, bsets_equal, rab, dab, saabb)
84 :
85 : INTEGER, INTENT(IN) :: la_max_set1, la_min_set1, npgfa1
86 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa1, zeta1
87 : INTEGER, INTENT(IN) :: la_max_set2, la_min_set2, npgfa2
88 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa2, zeta2
89 : INTEGER, INTENT(IN) :: lb_max_set1, lb_min_set1, npgfb1
90 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb1, zetb1
91 : INTEGER, INTENT(IN) :: lb_max_set2, lb_min_set2, npgfb2
92 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb2, zetb2
93 : LOGICAL, INTENT(IN) :: asets_equal, bsets_equal
94 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
95 : REAL(KIND=dp), INTENT(IN) :: dab
96 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(OUT) :: saabb
97 :
98 : CHARACTER(len=*), PARAMETER :: routineN = 'overlap_aabb'
99 :
100 : INTEGER :: coa, cob, handle, ia, ib, ipgf, ja, jb, jpgf, jpgf_start, kpgf, la_max, la_min, &
101 : lb_max, lb_min, ldrr, lpgf, lpgf_start, ncoa1, ncoa2, ncob1, ncob2
102 9 : INTEGER, ALLOCATABLE :: map_a(:, :), map_b(:, :)
103 : INTEGER, DIMENSION(3) :: na, nb
104 : REAL(KIND=dp) :: f0, zeta, zetb, zetp
105 9 : REAL(KIND=dp), ALLOCATABLE :: s(:, :)
106 9 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
107 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp
108 :
109 9 : CALL timeset(routineN, handle)
110 5736789 : saabb = 0.0_dp
111 9 : la_max = la_max_set1 + la_max_set2
112 9 : la_min = la_min_set1 + la_min_set2
113 9 : lb_max = lb_max_set1 + lb_max_set2
114 9 : lb_min = lb_min_set1 + lb_min_set2
115 9 : ldrr = MAX(la_max, lb_max) + 1
116 45 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
117 36 : ALLOCATE (s(ncoset(la_min - 1) + 1:ncoset(la_max), ncoset(lb_min - 1) + 1:ncoset(lb_max)))
118 : ALLOCATE (map_a(ncoset(la_min_set1 - 1) + 1:ncoset(la_max_set1), &
119 : ncoset(la_min_set2 - 1) + 1:ncoset(la_max_set2)), &
120 : map_b(ncoset(lb_min_set1 - 1) + 1:ncoset(lb_max_set1), &
121 63 : ncoset(lb_min_set2 - 1) + 1:ncoset(lb_max_set2)))
122 :
123 : ! Same-center products add Cartesian degrees; the maps are shared by all primitives.
124 63 : DO ja = LBOUND(map_a, 2), UBOUND(map_a, 2)
125 261 : DO ia = LBOUND(map_a, 1), UBOUND(map_a, 1)
126 576 : na = indco(:, ia) + indco(:, ja)
127 180 : map_a(ia, ja) = coset(na(1), na(2), na(3))
128 : END DO
129 : END DO
130 63 : DO jb = LBOUND(map_b, 2), UBOUND(map_b, 2)
131 261 : DO ib = LBOUND(map_b, 1), UBOUND(map_b, 1)
132 576 : nb = indco(:, ib) + indco(:, jb)
133 180 : map_b(ib, jb) = coset(nb(1), nb(2), nb(3))
134 : END DO
135 : END DO
136 :
137 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
138 :
139 72 : DO ipgf = 1, npgfa1
140 63 : ncoa1 = (ipgf - 1)*ncoset(la_max_set1)
141 63 : jpgf_start = MERGE(ipgf, 1, asets_equal)
142 :
143 324 : DO jpgf = jpgf_start, npgfa2
144 :
145 252 : ncoa2 = (jpgf - 1)*ncoset(la_max_set2)
146 252 : zeta = zeta1(ipgf) + zeta2(jpgf)
147 :
148 2079 : DO kpgf = 1, npgfb1
149 :
150 1764 : ncob1 = (kpgf - 1)*ncoset(lb_max_set1)
151 1764 : lpgf_start = MERGE(kpgf, 1, bsets_equal)
152 :
153 9072 : DO lpgf = lpgf_start, npgfb2
154 7056 : ncob2 = (lpgf - 1)*ncoset(lb_max_set2)
155 :
156 : ! *** Screening ***
157 7056 : IF (MIN(rpgfa1(ipgf), rpgfa2(jpgf)) + &
158 : MIN(rpgfb1(kpgf), rpgfb2(lpgf)) < dab) CYCLE
159 :
160 6909 : zetb = zetb1(kpgf) + zetb2(lpgf)
161 :
162 : ! *** Calculate some prefactors ***
163 :
164 6909 : zetp = 1.0_dp/(zeta + zetb)
165 :
166 6909 : f0 = SQRT((pi*zetp)**3)*EXP(-zeta*zetb*zetp*dab*dab)
167 27636 : rap(:) = zetb*zetp*rab(:)
168 27636 : rbp(:) = -zeta*zetp*rab(:)
169 :
170 6909 : CALL os_rr_ovlp(rap, la_max, rbp, lb_max, 1.0_dp/zetp, ldrr, rr)
171 :
172 75999 : DO cob = ncoset(lb_min - 1) + 1, ncoset(lb_max)
173 276360 : nb = indco(:, cob)
174 766899 : DO coa = ncoset(la_min - 1) + 1, ncoset(la_max)
175 2763600 : na = indco(:, coa)
176 759990 : s(coa, cob) = f0*rr(na(1), nb(1), 1)*rr(na(2), nb(2), 2)*rr(na(3), nb(3), 3)
177 : END DO
178 : END DO
179 :
180 : ! *** Store the primitive overlap integrals ***
181 36309 : DO jb = ncoset(lb_min_set2 - 1) + 1, ncoset(lb_max_set2)
182 145236 : DO ib = ncoset(lb_min_set1 - 1) + 1, ncoset(lb_max_set1)
183 110544 : cob = map_b(ib, jb)
184 580356 : DO ja = ncoset(la_min_set2 - 1) + 1, ncoset(la_max_set2)
185 2321424 : DO ia = ncoset(la_min_set1 - 1) + 1, ncoset(la_max_set1)
186 1768704 : coa = map_a(ia, ja)
187 1768704 : saabb(ncoa1 + ia, ncoa2 + ja, ncob1 + ib, ncob2 + jb) = s(coa, cob)
188 1768704 : IF (asets_equal) saabb(ncoa2 + ja, ncoa1 + ia, ncob1 + ib, ncob2 + jb) = s(coa, cob)
189 1768704 : IF (bsets_equal) saabb(ncoa1 + ia, ncoa2 + ja, ncob2 + jb, ncob1 + ib) = s(coa, cob)
190 2210880 : IF (asets_equal .AND. bsets_equal) THEN
191 1768704 : saabb(ncoa2 + ja, ncoa1 + ia, ncob2 + jb, ncob1 + ib) = s(coa, cob)
192 : END IF
193 : END DO
194 : END DO
195 : END DO
196 : END DO
197 :
198 : END DO
199 : END DO
200 : END DO
201 : END DO
202 :
203 9 : DEALLOCATE (rr, s, map_a, map_b)
204 9 : CALL timestop(handle)
205 :
206 9 : END SUBROUTINE overlap_aabb
207 :
208 : END MODULE ai_overlap_aabb
|
