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 ...
76 : !> \param s ...
77 : !> \param lds ...
78 : !> \date 06.2014
79 : !> \author Dorothea Golze
80 : ! **************************************************************************************************
81 9 : SUBROUTINE overlap_aabb(la_max_set1, la_min_set1, npgfa1, rpgfa1, zeta1, &
82 9 : la_max_set2, la_min_set2, npgfa2, rpgfa2, zeta2, &
83 9 : lb_max_set1, lb_min_set1, npgfb1, rpgfb1, zetb1, &
84 18 : lb_max_set2, lb_min_set2, npgfb2, rpgfb2, zetb2, &
85 9 : asets_equal, bsets_equal, rab, dab, saabb, s, lds)
86 :
87 : INTEGER, INTENT(IN) :: la_max_set1, la_min_set1, npgfa1
88 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa1, zeta1
89 : INTEGER, INTENT(IN) :: la_max_set2, la_min_set2, npgfa2
90 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa2, zeta2
91 : INTEGER, INTENT(IN) :: lb_max_set1, lb_min_set1, npgfb1
92 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb1, zetb1
93 : INTEGER, INTENT(IN) :: lb_max_set2, lb_min_set2, npgfb2
94 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb2, zetb2
95 : LOGICAL, INTENT(IN) :: asets_equal, bsets_equal
96 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
97 : REAL(KIND=dp), INTENT(IN) :: dab
98 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
99 : INTENT(INOUT) :: saabb
100 : INTEGER, INTENT(IN) :: lds
101 : REAL(KIND=dp), DIMENSION(lds, lds), INTENT(INOUT) :: s
102 :
103 : CHARACTER(len=*), PARAMETER :: routineN = 'overlap_aabb'
104 :
105 : INTEGER :: ax, ay, az, bx, by, bz, coa, cob, handle, i, ia, ib, ipgf, j, ja, jb, jpgf, &
106 : jpgf_start, kpgf, la, la_max, la_min, lb, lb_max, lb_min, ldrr, lpgf, lpgf_start, ncoa1, &
107 : ncoa2, ncob1, ncob2
108 : INTEGER, DIMENSION(3) :: na, naa, nb, nbb, nia, nib, nja, njb
109 : REAL(KIND=dp) :: f0, zeta, zetb, zetp
110 9 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
111 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp
112 :
113 9 : CALL timeset(routineN, handle)
114 9 : ldrr = MAX(la_max_set1 + la_max_set2, lb_max_set1 + lb_max_set2) + 1
115 45 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
116 :
117 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
118 :
119 9 : ncoa1 = 0
120 9 : ncoa2 = 0
121 9 : ncob1 = 0
122 9 : ncob2 = 0
123 :
124 72 : DO ipgf = 1, npgfa1
125 :
126 63 : ncoa2 = 0
127 :
128 63 : IF (asets_equal) THEN
129 252 : jpgf_start = ipgf
130 252 : DO i = 1, jpgf_start - 1
131 252 : ncoa2 = ncoa2 + ncoset(la_max_set2)
132 : END DO
133 : ELSE
134 : jpgf_start = 1
135 : END IF
136 :
137 315 : DO jpgf = jpgf_start, npgfa2
138 :
139 252 : ncob1 = 0
140 252 : zeta = zeta1(ipgf) + zeta2(jpgf)
141 252 : la_max = la_max_set1 + la_max_set2
142 252 : la_min = la_min_set1 + la_min_set2
143 :
144 2016 : DO kpgf = 1, npgfb1
145 :
146 1764 : ncob2 = 0
147 :
148 1764 : IF (bsets_equal) THEN
149 7056 : lpgf_start = kpgf
150 7056 : DO i = 1, lpgf_start - 1
151 7056 : ncob2 = ncob2 + ncoset(lb_max_set2)
152 : END DO
153 : ELSE
154 : lpgf_start = 1
155 : END IF
156 :
157 8820 : DO lpgf = lpgf_start, npgfb2
158 :
159 : ! *** Screening ***
160 : IF ((rpgfa1(ipgf) + rpgfb1(kpgf) < dab) .OR. &
161 : (rpgfa2(jpgf) + rpgfb1(kpgf) < dab) .OR. &
162 7056 : (rpgfa1(ipgf) + rpgfb2(lpgf) < dab) .OR. &
163 : (rpgfa2(jpgf) + rpgfb2(lpgf) < dab)) THEN
164 735 : DO jb = ncoset(lb_min_set2 - 1) + 1, ncoset(lb_max_set2)
165 3087 : DO ib = ncoset(lb_min_set1 - 1) + 1, ncoset(lb_max_set1)
166 12348 : DO ja = ncoset(la_min_set2 - 1) + 1, ncoset(la_max_set2)
167 49392 : DO ia = ncoset(la_min_set1 - 1) + 1, ncoset(la_max_set1)
168 37632 : saabb(ncoa1 + ia, ncoa2 + ja, ncob1 + ib, ncob2 + jb) = 0._dp
169 37632 : IF (asets_equal) saabb(ncoa2 + ja, ncoa1 + ia, ncob1 + ib, ncob2 + jb) = 0._dp
170 37632 : IF (bsets_equal) saabb(ncoa1 + ia, ncoa2 + ja, ncob2 + jb, ncob1 + ib) = 0._dp
171 47040 : IF (asets_equal .AND. bsets_equal) THEN
172 37632 : saabb(ncoa2 + ja, ncoa1 + ia, ncob2 + jb, ncob1 + ib) = 0._dp
173 : END IF
174 : END DO
175 : END DO
176 : END DO
177 : END DO
178 147 : ncob2 = ncob2 + ncoset(lb_max_set2)
179 147 : CYCLE
180 : END IF
181 :
182 6909 : zetb = zetb1(kpgf) + zetb2(lpgf)
183 6909 : lb_max = lb_max_set1 + lb_max_set2
184 6909 : lb_min = lb_min_set1 + lb_min_set2
185 :
186 : ! *** Calculate some prefactors ***
187 :
188 6909 : zetp = 1.0_dp/(zeta + zetb)
189 :
190 6909 : f0 = SQRT((pi*zetp)**3)*EXP(-zeta*zetb*zetp*dab*dab)
191 27636 : rap(:) = zetb*zetp*rab(:)
192 27636 : rbp(:) = -zeta*zetp*rab(:)
193 :
194 6909 : CALL os_rr_ovlp(rap, la_max, rbp, lb_max, 1.0_dp/zetp, ldrr, rr)
195 :
196 27636 : DO lb = 0, lb_max
197 69090 : DO bx = 0, lb
198 131271 : DO by = 0, lb - bx
199 69090 : bz = lb - bx - by
200 69090 : cob = coset(bx, by, bz)
201 317814 : DO la = 0, la_max
202 690900 : DO ax = 0, la
203 1312710 : DO ay = 0, la - ax
204 690900 : az = la - ax - ay
205 690900 : coa = coset(ax, ay, az)
206 1105440 : s(coa, cob) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az, bz, 3)
207 : END DO
208 : END DO
209 : END DO
210 : END DO
211 : END DO
212 : END DO
213 :
214 : ! *** Store the primitive overlap integrals ***
215 34545 : DO jb = ncoset(lb_min_set2 - 1) + 1, ncoset(lb_max_set2)
216 110544 : njb(1:3) = indco(1:3, jb)
217 145089 : DO ib = ncoset(lb_min_set1 - 1) + 1, ncoset(lb_max_set1)
218 442176 : nib(1:3) = indco(1:3, ib)
219 442176 : nbb(1:3) = nib + njb
220 580356 : DO ja = ncoset(la_min_set2 - 1) + 1, ncoset(la_max_set2)
221 1768704 : nja(1:3) = indco(1:3, ja)
222 2321424 : DO ia = ncoset(la_min_set1 - 1) + 1, ncoset(la_max_set1)
223 7074816 : nia(1:3) = indco(1:3, ia)
224 7074816 : naa(1:3) = nia + nja
225 : ! now loop over all elements of s
226 19897920 : DO j = ncoset(lb_min - 1) + 1, ncoset(lb_max)
227 70748160 : nb(1:3) = indco(1:3, j)
228 196326144 : DO i = ncoset(la_min - 1) + 1, ncoset(la_max)
229 707481600 : na(1:3) = indco(1:3, i)
230 638944320 : IF (ALL(na == naa) .AND. ALL(nb == nbb)) THEN
231 1768704 : saabb(ncoa1 + ia, ncoa2 + ja, ncob1 + ib, ncob2 + jb) = s(i, j)
232 1768704 : IF (asets_equal) saabb(ncoa2 + ja, ncoa1 + ia, ncob1 + ib, ncob2 + jb) = s(i, j)
233 1768704 : IF (bsets_equal) saabb(ncoa1 + ia, ncoa2 + ja, ncob2 + jb, ncob1 + ib) = s(i, j)
234 1768704 : IF (asets_equal .AND. bsets_equal) THEN
235 1768704 : saabb(ncoa2 + ja, ncoa1 + ia, ncob2 + jb, ncob1 + ib) = s(i, j)
236 : END IF
237 : END IF
238 : END DO
239 : END DO
240 : END DO
241 : END DO
242 : END DO
243 : END DO
244 :
245 8673 : ncob2 = ncob2 + ncoset(lb_max_set2)
246 :
247 : END DO
248 :
249 2016 : ncob1 = ncob1 + ncoset(lb_max_set1)
250 :
251 : END DO
252 :
253 315 : ncoa2 = ncoa2 + ncoset(la_max_set2)
254 :
255 : END DO
256 :
257 72 : ncoa1 = ncoa1 + ncoset(la_max_set1)
258 :
259 : END DO
260 :
261 9 : DEALLOCATE (rr)
262 9 : CALL timestop(handle)
263 :
264 9 : END SUBROUTINE overlap_aabb
265 :
266 : END MODULE ai_overlap_aabb
|
