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 integrals over solid harmonic Gaussian(SHG) functions.
10 : !> Routines for (a|O(r12)|b) and overlap integrals (ab), (aba) and (abb).
11 : !> \par Literature (partly)
12 : !> T.J. Giese and D. M. York, J. Chem. Phys, 128, 064104 (2008)
13 : !> T. Helgaker, P Joergensen, J. Olsen, Molecular Electronic-Structure
14 : !> Theory, Wiley
15 : !> \par History
16 : !> created [04.2015]
17 : !> \author Dorothea Golze
18 : ! **************************************************************************************************
19 : MODULE construct_shg
20 : USE kinds, ONLY: dp
21 : USE mathconstants, ONLY: dfac,&
22 : fac
23 : USE orbital_pointers, ONLY: indso,&
24 : indso_inv,&
25 : nsoset
26 : #include "../base/base_uses.f90"
27 :
28 : IMPLICIT NONE
29 :
30 : PRIVATE
31 :
32 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'construct_shg'
33 :
34 : ! *** Public subroutines ***
35 : PUBLIC :: get_real_scaled_solid_harmonic, get_W_matrix, get_dW_matrix, &
36 : construct_int_shg_ab, construct_dev_shg_ab, construct_overlap_shg_aba, &
37 : dev_overlap_shg_aba, construct_overlap_shg_abb, dev_overlap_shg_abb
38 :
39 : CONTAINS
40 :
41 : ! **************************************************************************************************
42 : !> \brief computes the real scaled solid harmonics Rlm up to a given l
43 : !> \param Rlm_c cosine part of real scaled soldi harmonics
44 : !> \param Rlm_s sine part of real scaled soldi harmonics
45 : !> \param l maximal l quantum up to where Rlm is calculated
46 : !> \param r distance vector between a and b
47 : !> \param r2 square of distance vector
48 : ! **************************************************************************************************
49 3489179 : SUBROUTINE get_real_scaled_solid_harmonic(Rlm_c, Rlm_s, l, r, r2)
50 :
51 : INTEGER, INTENT(IN) :: l
52 : REAL(KIND=dp), DIMENSION(0:l, -2*l:2*l), &
53 : INTENT(OUT) :: Rlm_s, Rlm_c
54 : REAL(KIND=dp), DIMENSION(3) :: r
55 : REAL(KIND=dp) :: r2
56 :
57 : INTEGER :: li, mi, prefac
58 : REAL(KIND=dp) :: Rc, Rc_00, Rlm, Rmlm, Rplm, Rs, Rs_00, &
59 : temp_c
60 :
61 3489179 : Rc_00 = 1.0_dp
62 3489179 : Rs_00 = 0.0_dp
63 :
64 3489179 : Rlm_c(0, 0) = Rc_00
65 3489179 : Rlm_s(0, 0) = Rs_00
66 :
67 : ! generate elements Rmm
68 : ! start
69 3489179 : IF (l > 0) THEN
70 2066599 : Rc = -0.5_dp*r(1)*Rc_00
71 2066599 : Rs = -0.5_dp*r(2)*Rc_00
72 2066599 : Rlm_c(1, 1) = Rc
73 2066599 : Rlm_s(1, 1) = Rs
74 2066599 : Rlm_c(1, -1) = -Rc
75 2066599 : Rlm_s(1, -1) = Rs
76 : END IF
77 4051798 : DO li = 2, l
78 562619 : temp_c = (-r(1)*Rc + r(2)*Rs)/(REAL(2*(li - 1) + 2, dp))
79 562619 : Rs = (-r(2)*Rc - r(1)*Rs)/(REAL(2*(li - 1) + 2, dp))
80 562619 : Rc = temp_c
81 562619 : Rlm_c(li, li) = Rc
82 562619 : Rlm_s(li, li) = Rs
83 4051798 : IF (MODULO(li, 2) /= 0) THEN
84 148674 : Rlm_c(li, -li) = -Rc
85 148674 : Rlm_s(li, -li) = Rs
86 : ELSE
87 413945 : Rlm_c(li, -li) = Rc
88 413945 : Rlm_s(li, -li) = -Rs
89 : END IF
90 : END DO
91 :
92 6118397 : DO mi = 0, l - 1
93 2629218 : Rmlm = Rlm_c(mi, mi)
94 2629218 : Rlm = r(3)*Rlm_c(mi, mi)
95 2629218 : Rlm_c(mi + 1, mi) = Rlm
96 2629218 : IF (MODULO(mi, 2) /= 0) THEN
97 413945 : Rlm_c(mi + 1, -mi) = -Rlm
98 : ELSE
99 2215273 : Rlm_c(mi + 1, -mi) = Rlm
100 : END IF
101 6954354 : DO li = mi + 2, l
102 835957 : prefac = (li + mi)*(li - mi)
103 835957 : Rplm = (REAL(2*li - 1, dp)*r(3)*Rlm - r2*Rmlm)/REAL(prefac, dp)
104 835957 : Rmlm = Rlm
105 835957 : Rlm = Rplm
106 835957 : Rlm_c(li, mi) = Rlm
107 3465175 : IF (MODULO(mi, 2) /= 0) THEN
108 211006 : Rlm_c(li, -mi) = -Rlm
109 : ELSE
110 624951 : Rlm_c(li, -mi) = Rlm
111 : END IF
112 : END DO
113 : END DO
114 4051798 : DO mi = 1, l - 1
115 562619 : Rmlm = Rlm_s(mi, mi)
116 562619 : Rlm = r(3)*Rlm_s(mi, mi)
117 562619 : Rlm_s(mi + 1, mi) = Rlm
118 562619 : IF (MODULO(mi, 2) /= 0) THEN
119 413945 : Rlm_s(mi + 1, -mi) = Rlm
120 : ELSE
121 148674 : Rlm_s(mi + 1, -mi) = -Rlm
122 : END IF
123 4325136 : DO li = mi + 2, l
124 273338 : prefac = (li + mi)*(li - mi)
125 273338 : Rplm = (REAL(2*li - 1, dp)*r(3)*Rlm - r2*Rmlm)/REAL(prefac, dp)
126 273338 : Rmlm = Rlm
127 273338 : Rlm = Rplm
128 273338 : Rlm_s(li, mi) = Rlm
129 835957 : IF (MODULO(mi, 2) /= 0) THEN
130 211006 : Rlm_s(li, -mi) = Rlm
131 : ELSE
132 62332 : Rlm_s(li, -mi) = -Rlm
133 : END IF
134 : END DO
135 : END DO
136 :
137 3489179 : END SUBROUTINE get_real_scaled_solid_harmonic
138 :
139 : ! **************************************************************************************************
140 : !> \brief Calculate the prefactor A(l,m) = (-1)^m \sqrt[(2-delta(m,0))(l+m)!(l-m)!]
141 : !> \param lmax maximal l quantum number
142 : !> \param A matrix storing the prefactor for a given l and m
143 : ! **************************************************************************************************
144 3489179 : SUBROUTINE get_Alm(lmax, A)
145 :
146 : INTEGER, INTENT(IN) :: lmax
147 : REAL(KIND=dp), DIMENSION(0:lmax, 0:lmax), &
148 : INTENT(INOUT) :: A
149 :
150 : INTEGER :: l, m
151 : REAL(KIND=dp) :: temp
152 :
153 9607576 : DO l = 0, lmax
154 19191148 : DO m = 0, l
155 9583572 : temp = SQRT(fac(l + m)*fac(l - m))
156 9583572 : IF (MODULO(m, 2) /= 0) temp = -temp
157 9583572 : IF (m /= 0) temp = temp*SQRT(2.0_dp)
158 15701969 : A(l, m) = temp
159 : END DO
160 : END DO
161 :
162 3489179 : END SUBROUTINE get_Alm
163 :
164 : ! **************************************************************************************************
165 : !> \brief calculates the prefactors for the derivatives of the W matrix
166 : !> \param lmax maximal l quantum number
167 : !> \param dA_p = A(l,m)/A(l-1,m+1)
168 : !> \param dA_m = A(l,m)/A(l-1,m-1)
169 : !> \param dA = A(l,m)/A(l-1,m)
170 : !> \note for m=0, W_l-1,-1 can't be read from Waux_mat, but we use
171 : !> W_l-1,-1 = -W_l-1,1 [cc(1), cs(2)] or W_l-1,-1 = W_l-1,1 [[sc(3), ss(4)], i.e.
172 : !> effectively we multiply dA_p by 2
173 : ! **************************************************************************************************
174 12204 : SUBROUTINE get_dA_prefactors(lmax, dA_p, dA_m, dA)
175 :
176 : INTEGER, INTENT(IN) :: lmax
177 : REAL(KIND=dp), DIMENSION(0:lmax, 0:lmax), &
178 : INTENT(INOUT) :: dA_p, dA_m, dA
179 :
180 : INTEGER :: l, m
181 : REAL(KIND=dp) :: bm, bm_m, bm_p
182 :
183 85996 : DO l = 0, lmax
184 348208 : DO m = 0, l
185 262212 : bm = 1.0_dp
186 262212 : bm_m = 1.0_dp
187 262212 : bm_p = 1.0_dp
188 262212 : IF (m /= 0) bm = SQRT(2.0_dp)
189 188420 : IF (m - 1 /= 0) bm_m = SQRT(2.0_dp)
190 200624 : IF (m + 1 /= 0) bm_p = SQRT(2.0_dp)
191 262212 : dA_p(l, m) = -bm/bm_p*SQRT(REAL((l - m)*(l - m - 1), dp))
192 262212 : dA_m(l, m) = -bm/bm_m*SQRT(REAL((l + m)*(l + m - 1), dp))
193 262212 : dA(l, m) = 2.0_dp*SQRT(REAL((l + m)*(l - m), dp))
194 336004 : IF (m == 0) dA_p(l, m) = 2.0_dp*dA_p(l, m)
195 : END DO
196 : END DO
197 12204 : END SUBROUTINE get_dA_prefactors
198 :
199 : ! **************************************************************************************************
200 : !> \brief calculates the angular dependent-part of the SHG integrals,
201 : !> transformation matrix W, see literature above
202 : !> \param lamax array of maximal l quantum number on a;
203 : !> lamax(lb) with lb= 0..lbmax
204 : !> \param lbmax maximal l quantum number on b
205 : !> \param lmax maximal l quantum number
206 : !> \param Rc cosine part of real scaled solid harmonics
207 : !> \param Rs sine part of real scaled solid harmonics
208 : !> \param Waux_mat stores the angular-dependent part of the SHG integrals
209 : ! **************************************************************************************************
210 3489179 : SUBROUTINE get_W_matrix(lamax, lbmax, lmax, Rc, Rs, Waux_mat)
211 :
212 : INTEGER, DIMENSION(:), POINTER :: lamax
213 : INTEGER, INTENT(IN) :: lbmax, lmax
214 : REAL(KIND=dp), DIMENSION(0:lmax, -2*lmax:2*lmax), &
215 : INTENT(IN) :: Rc, Rs
216 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: Waux_mat
217 :
218 : INTEGER :: j, k, la, labmin, laj, lb, lbj, ma, &
219 : ma_m, ma_p, mb, mb_m, mb_p, nla, nlb
220 : REAL(KIND=dp) :: A_jk, A_lama, A_lbmb, Alm_fac, delta_k, prefac, Rca_m, Rca_p, Rcb_m, Rcb_p, &
221 : Rsa_m, Rsa_p, Rsb_m, Rsb_p, sign_fac, Wa(4), Wb(4), Wmat(4)
222 3489179 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: A
223 :
224 3489179 : Wa(:) = 0.0_dp
225 3489179 : Wb(:) = 0.0_dp
226 3489179 : Wmat(:) = 0.0_dp
227 :
228 13956716 : ALLOCATE (A(0:lmax, 0:lmax))
229 3489179 : CALL get_Alm(lmax, A)
230 :
231 8709725 : DO lb = 0, lbmax
232 5220546 : nlb = nsoset(lb - 1)
233 17551773 : DO la = 0, lamax(lb)
234 8842048 : nla = nsoset(la - 1)
235 8842048 : labmin = MIN(la, lb)
236 28047985 : DO mb = 0, lb
237 13985391 : A_lbmb = A(lb, mb)
238 13985391 : IF (MODULO(lb, 2) /= 0) A_lbmb = -A_lbmb
239 48042460 : DO ma = 0, la
240 25215021 : A_lama = A(la, ma)
241 25215021 : Alm_fac = A_lama*A_lbmb
242 85367899 : DO j = 0, labmin
243 46167487 : laj = la - j
244 46167487 : lbj = lb - j
245 46167487 : prefac = Alm_fac*REAL(2**(la + lb - j), dp)*dfac(2*j - 1)
246 46167487 : delta_k = 0.5_dp
247 46167487 : Wmat = 0.0_dp
248 122205414 : DO k = 0, j
249 76037927 : ma_m = ma - k
250 76037927 : ma_p = ma + k
251 76037927 : IF (laj < ABS(ma_m) .AND. laj < ABS(ma_p)) CYCLE
252 56964280 : mb_m = mb - k
253 56964280 : mb_p = mb + k
254 56964280 : IF (lbj < ABS(mb_m) .AND. lbj < ABS(mb_p)) CYCLE
255 45107448 : IF (k /= 0) delta_k = 1.0_dp
256 45107448 : A_jk = fac(j + k)*fac(j - k)
257 45107448 : IF (k /= 0) A_jk = 2.0_dp*A_jk
258 45107448 : IF (MODULO(k, 2) /= 0) THEN
259 : sign_fac = -1.0_dp
260 : ELSE
261 33864797 : sign_fac = 1.0_dp
262 : END IF
263 45107448 : Rca_m = Rc(laj, ma_m)
264 45107448 : Rsa_m = Rs(laj, ma_m)
265 45107448 : Rca_p = Rc(laj, ma_p)
266 45107448 : Rsa_p = Rs(laj, ma_p)
267 45107448 : Rcb_m = Rc(lbj, mb_m)
268 45107448 : Rsb_m = Rs(lbj, mb_m)
269 45107448 : Rcb_p = Rc(lbj, mb_p)
270 45107448 : Rsb_p = Rs(lbj, mb_p)
271 45107448 : Wa(1) = delta_k*(Rca_m + sign_fac*Rca_p)
272 45107448 : Wb(1) = delta_k*(Rcb_m + sign_fac*Rcb_p)
273 45107448 : Wa(2) = -Rsa_m + sign_fac*Rsa_p
274 45107448 : Wb(2) = -Rsb_m + sign_fac*Rsb_p
275 45107448 : Wmat(1) = Wmat(1) + prefac/A_jk*(Wa(1)*Wb(1) + Wa(2)*Wb(2))
276 45107448 : IF (mb > 0) THEN
277 23743998 : Wb(3) = delta_k*(Rsb_m + sign_fac*Rsb_p)
278 23743998 : Wb(4) = Rcb_m - sign_fac*Rcb_p
279 23743998 : Wmat(2) = Wmat(2) + prefac/A_jk*(Wa(1)*Wb(3) + Wa(2)*Wb(4))
280 : END IF
281 45107448 : IF (ma > 0) THEN
282 24884347 : Wa(3) = delta_k*(Rsa_m + sign_fac*Rsa_p)
283 24884347 : Wa(4) = Rca_m - sign_fac*Rca_p
284 24884347 : Wmat(3) = Wmat(3) + prefac/A_jk*(Wa(3)*Wb(1) + Wa(4)*Wb(2))
285 : END IF
286 91274935 : IF (ma > 0 .AND. mb > 0) THEN
287 15368281 : Wmat(4) = Wmat(4) + prefac/A_jk*(Wa(3)*Wb(3) + Wa(4)*Wb(4))
288 : END IF
289 : END DO
290 46167487 : Waux_mat(j + 1, nla + la + 1 + ma, nlb + lb + 1 + mb) = Wmat(1)
291 46167487 : IF (mb > 0) Waux_mat(j + 1, nla + la + 1 + ma, nlb + lb + 1 - mb) = Wmat(2)
292 46167487 : IF (ma > 0) Waux_mat(j + 1, nla + la + 1 - ma, nlb + lb + 1 + mb) = Wmat(3)
293 71382508 : IF (ma > 0 .AND. mb > 0) Waux_mat(j + 1, nla + la + 1 - ma, nlb + lb + 1 - mb) = Wmat(4)
294 : END DO
295 : END DO
296 : END DO
297 : END DO
298 : END DO
299 :
300 3489179 : DEALLOCATE (A)
301 :
302 3489179 : END SUBROUTINE get_W_matrix
303 :
304 : ! **************************************************************************************************
305 : !> \brief calculates derivatives of transformation matrix W,
306 : !> \param lamax array of maximal l quantum number on a;
307 : !> lamax(lb) with lb= 0..lbmax
308 : !> \param lbmax maximal l quantum number on b
309 : !> \param Waux_mat stores the angular-dependent part of the SHG integrals
310 : !> \param dWaux_mat stores the derivatives of the angular-dependent part of
311 : !> the SHG integrals
312 : ! **************************************************************************************************
313 12204 : SUBROUTINE get_dW_matrix(lamax, lbmax, Waux_mat, dWaux_mat)
314 :
315 : INTEGER, DIMENSION(:), POINTER :: lamax
316 : INTEGER, INTENT(IN) :: lbmax
317 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: Waux_mat
318 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
319 : INTENT(INOUT) :: dWaux_mat
320 :
321 : INTEGER :: ima, imam, imb, imbm, ipa, ipam, ipb, &
322 : ipbm, j, jmax, la, labm, labmin, lamb, &
323 : lb, lmax, ma, mb, nla, nlam, nlb, nlbm
324 : REAL(KIND=dp) :: dAa, dAa_m, dAa_p, dAb, dAb_m, dAb_p
325 12204 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dA, dA_m, dA_p, Wam, Wamm, Wamp, Wbm, &
326 12204 : Wbmm, Wbmp
327 :
328 61708 : jmax = MIN(MAXVAL(lamax), lbmax)
329 73224 : ALLOCATE (Wam(0:jmax, 4), Wamm(0:jmax, 4), Wamp(0:jmax, 4))
330 48816 : ALLOCATE (Wbm(0:jmax, 4), Wbmm(0:jmax, 4), Wbmp(0:jmax, 4))
331 :
332 : !*** get dA_p=A(l,m)/A(l-1,m+1)
333 : !*** get dA_m=A(l,m)/A(l-1,m-1)
334 : !*** get dA=2*A(l,m)/A(l-1,m)
335 61708 : lmax = MAX(MAXVAL(lamax), lbmax)
336 97632 : ALLOCATE (dA_p(0:lmax, 0:lmax), dA_m(0:lmax, 0:lmax), dA(0:lmax, 0:lmax))
337 12204 : CALL get_dA_prefactors(lmax, dA_p, dA_m, dA)
338 :
339 61708 : DO lb = 0, lbmax
340 49504 : nlb = nsoset(lb - 1)
341 49504 : nlbm = 0
342 49504 : IF (lb > 0) nlbm = nsoset(lb - 2)
343 336666 : DO la = 0, lamax(lb)
344 274958 : nla = nsoset(la - 1)
345 274958 : nlam = 0
346 274958 : IF (la > 0) nlam = nsoset(la - 2)
347 274958 : labmin = MIN(la, lb)
348 274958 : lamb = MIN(la - 1, lb)
349 274958 : labm = MIN(la, lb - 1)
350 989834 : DO mb = 0, lb
351 665372 : dAb = dA(lb, mb)
352 665372 : dAb_p = dA_p(lb, mb)
353 665372 : dAb_m = dA_m(lb, mb)
354 665372 : ipb = nlb + lb + mb + 1
355 665372 : imb = nlb + lb - mb + 1
356 665372 : ipbm = nlbm + lb + mb
357 665372 : imbm = nlbm + lb - mb
358 3122850 : DO ma = 0, la
359 2182520 : dAa = dA(la, ma)
360 2182520 : dAa_p = dA_p(la, ma)
361 2182520 : dAa_m = dA_m(la, ma)
362 2182520 : ipa = nla + la + ma + 1
363 2182520 : ima = nla + la - ma + 1
364 2182520 : ipam = nlam + la + ma
365 2182520 : imam = nlam + la - ma
366 47750940 : Wam(:, :) = 0.0_dp
367 47750940 : Wamm(:, :) = 0.0_dp
368 47750940 : Wamp(:, :) = 0.0_dp
369 : !*** Wam: la-1, ma
370 2182520 : IF (ma <= la - 1) THEN
371 5046254 : Wam(0:lamb, 1) = Waux_mat(1:lamb + 1, ipam, ipb)
372 3731073 : IF (mb > 0) Wam(0:lamb, 2) = Waux_mat(1:lamb + 1, ipam, imb)
373 3946483 : IF (ma > 0) Wam(0:lamb, 3) = Waux_mat(1:lamb + 1, imam, ipb)
374 3039982 : IF (ma > 0 .AND. mb > 0) Wam(0:lamb, 4) = Waux_mat(1:lamb + 1, imam, imb)
375 : END IF
376 : !*** Wamm: la-1, ma-1
377 2182520 : IF (ma - 1 >= 0) THEN
378 5046254 : Wamm(0:lamb, 1) = Waux_mat(1:lamb + 1, ipam - 1, ipb)
379 3731073 : IF (mb > 0) Wamm(0:lamb, 2) = Waux_mat(1:lamb + 1, ipam - 1, imb)
380 3946483 : IF (ma - 1 > 0) Wamm(0:lamb, 3) = Waux_mat(1:lamb + 1, imam + 1, ipb) !order: e.g. -1 0 1, if < 0 |m|, -1 means -m+1
381 3039982 : IF (ma - 1 > 0 .AND. mb > 0) Wamm(0:lamb, 4) = Waux_mat(1:lamb + 1, imam + 1, imb)
382 : END IF
383 : !*** Wamp: la-1, ma+1
384 2182520 : IF (ma + 1 <= la - 1) THEN
385 3407029 : Wamp(0:lamb, 1) = Waux_mat(1:lamb + 1, ipam + 1, ipb)
386 2500528 : IF (mb > 0) Wamp(0:lamb, 2) = Waux_mat(1:lamb + 1, ipam + 1, imb)
387 3407029 : IF (ma + 1 > 0) Wamp(0:lamb, 3) = Waux_mat(1:lamb + 1, imam - 1, ipb)
388 2500528 : IF (ma + 1 > 0 .AND. mb > 0) Wamp(0:lamb, 4) = Waux_mat(1:lamb + 1, imam - 1, imb)
389 : END IF
390 47750940 : Wbm(:, :) = 0.0_dp
391 47750940 : Wbmm(:, :) = 0.0_dp
392 47750940 : Wbmp(:, :) = 0.0_dp
393 : !*** Wbm: lb-1, mb
394 2182520 : IF (mb <= lb - 1) THEN
395 3855304 : Wbm(0:labm, 1) = Waux_mat(1:labm + 1, ipa, ipbm)
396 2675803 : IF (mb > 0) Wbm(0:labm, 2) = Waux_mat(1:labm + 1, ipa, imbm)
397 3108049 : IF (ma > 0) Wbm(0:labm, 3) = Waux_mat(1:labm + 1, ima, ipbm)
398 2264313 : IF (ma > 0 .AND. mb > 0) Wbm(0:labm, 4) = Waux_mat(1:labm + 1, ima, imbm)
399 : END IF
400 : !*** Wbmm: lb-1, mb-1
401 2182520 : IF (mb - 1 >= 0) THEN
402 3855304 : Wbmm(0:labm, 1) = Waux_mat(1:labm + 1, ipa, ipbm - 1)
403 2675803 : IF (mb - 1 > 0) Wbmm(0:labm, 2) = Waux_mat(1:labm + 1, ipa, imbm + 1)
404 3108049 : IF (ma > 0) Wbmm(0:labm, 3) = Waux_mat(1:labm + 1, ima, ipbm - 1)
405 2264313 : IF (ma > 0 .AND. mb - 1 > 0) Wbmm(0:labm, 4) = Waux_mat(1:labm + 1, ima, imbm + 1)
406 : END IF
407 : !*** Wbmp: lb-1, mb+1
408 2182520 : IF (mb + 1 <= lb - 1) THEN
409 2006040 : Wbmp(0:labm, 1) = Waux_mat(1:labm + 1, ipa, ipbm + 1)
410 2006040 : IF (mb + 1 > 0) Wbmp(0:labm, 2) = Waux_mat(1:labm + 1, ipa, imbm - 1)
411 1594550 : IF (ma > 0) Wbmp(0:labm, 3) = Waux_mat(1:labm + 1, ima, ipbm + 1)
412 1594550 : IF (ma > 0 .AND. mb + 1 > 0) Wbmp(0:labm, 4) = Waux_mat(1:labm + 1, ima, imbm - 1)
413 : END IF
414 8334934 : DO j = 0, labmin
415 : !*** x component
416 : dWaux_mat(1, j + 1, ipa, ipb) = dAa_p*Wamp(j, 1) - dAa_m*Wamm(j, 1) &
417 5487042 : - dAb_p*Wbmp(j, 1) + dAb_m*Wbmm(j, 1)
418 5487042 : IF (mb > 0) THEN
419 : dWaux_mat(1, j + 1, ipa, imb) = dAa_p*Wamp(j, 2) - dAa_m*Wamm(j, 2) &
420 3507749 : - dAb_p*Wbmp(j, 2) + dAb_m*Wbmm(j, 2)
421 : END IF
422 5487042 : IF (ma > 0) THEN
423 : dWaux_mat(1, j + 1, ima, ipb) = dAa_p*Wamp(j, 3) - dAa_m*Wamm(j, 3) &
424 4000776 : - dAb_p*Wbmp(j, 3) + dAb_m*Wbmm(j, 3)
425 : END IF
426 5487042 : IF (ma > 0 .AND. mb > 0) THEN
427 : dWaux_mat(1, j + 1, ima, imb) = dAa_p*Wamp(j, 4) - dAa_m*Wamm(j, 4) &
428 2555792 : - dAb_p*Wbmp(j, 4) + dAb_m*Wbmm(j, 4)
429 : END IF
430 :
431 : !**** y component
432 : dWaux_mat(2, j + 1, ipa, ipb) = dAa_p*Wamp(j, 3) + dAa_m*Wamm(j, 3) &
433 5487042 : - dAb_p*Wbmp(j, 2) - dAb_m*Wbmm(j, 2)
434 5487042 : IF (mb > 0) THEN
435 : dWaux_mat(2, j + 1, ipa, imb) = dAa_p*Wamp(j, 4) + dAa_m*Wamm(j, 4) &
436 3507749 : + dAb_p*Wbmp(j, 1) + dAb_m*Wbmm(j, 1)
437 : END IF
438 5487042 : IF (ma > 0) THEN
439 : dWaux_mat(2, j + 1, ima, ipb) = -dAa_p*Wamp(j, 1) - dAa_m*Wamm(j, 1) &
440 4000776 : - dAb_p*Wbmp(j, 4) - dAb_m*Wbmm(j, 4)
441 : END IF
442 5487042 : IF (ma > 0 .AND. mb > 0) THEN
443 : dWaux_mat(2, j + 1, ima, imb) = -dAa_p*Wamp(j, 2) - dAa_m*Wamm(j, 2) &
444 2555792 : + dAb_p*Wbmp(j, 3) + dAb_m*Wbmm(j, 3)
445 : END IF
446 : !**** z compnent
447 5487042 : dWaux_mat(3, j + 1, ipa, ipb) = dAa*Wam(j, 1) - dAb*Wbm(j, 1)
448 5487042 : IF (mb > 0) THEN
449 3507749 : dWaux_mat(3, j + 1, ipa, imb) = dAa*Wam(j, 2) - dAb*Wbm(j, 2)
450 : END IF
451 5487042 : IF (ma > 0) THEN
452 4000776 : dWaux_mat(3, j + 1, ima, ipb) = dAa*Wam(j, 3) - dAb*Wbm(j, 3)
453 : END IF
454 7669562 : IF (ma > 0 .AND. mb > 0) THEN
455 2555792 : dWaux_mat(3, j + 1, ima, imb) = dAa*Wam(j, 4) - dAb*Wbm(j, 4)
456 : END IF
457 :
458 : END DO
459 : END DO
460 : END DO
461 : END DO
462 : END DO
463 :
464 12204 : DEALLOCATE (Wam, Wamm, Wamp)
465 12204 : DEALLOCATE (Wbm, Wbmm, Wbmp)
466 12204 : DEALLOCATE (dA, dA_p, dA_m)
467 :
468 12204 : END SUBROUTINE get_dW_matrix
469 :
470 : ! **************************************************************************************************
471 : !> \brief calculates [ab] SHG overlap integrals using precomputed angular-
472 : !> dependent part
473 : !> \param la set of l quantum number on a
474 : !> \param first_sgfa indexing
475 : !> \param nshella number of shells for a
476 : !> \param lb set of l quantum number on b
477 : !> \param first_sgfb indexing
478 : !> \param nshellb number of shells for b
479 : !> \param swork_cont contracted and normalized [s|s] integrals
480 : !> \param Waux_mat precomputed angular-dependent part
481 : !> \param sab contracted integral of spherical harmonic Gaussianslm
482 : ! **************************************************************************************************
483 22132224 : SUBROUTINE construct_int_shg_ab(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
484 22132224 : swork_cont, Waux_mat, sab)
485 :
486 : INTEGER, DIMENSION(:), INTENT(IN) :: la, first_sgfa
487 : INTEGER, INTENT(IN) :: nshella
488 : INTEGER, DIMENSION(:), INTENT(IN) :: lb, first_sgfb
489 : INTEGER, INTENT(IN) :: nshellb
490 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: swork_cont, Waux_mat
491 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
492 :
493 : INTEGER :: fnla, fnlb, fsgfa, fsgfb, ishella, j, &
494 : jshellb, labmin, lai, lbj, lnla, lnlb, &
495 : lsgfa, lsgfb, mai, mbj
496 : REAL(KIND=dp) :: prefac
497 :
498 44575103 : DO jshellb = 1, nshellb
499 22442879 : lbj = lb(jshellb)
500 22442879 : fnlb = nsoset(lbj - 1) + 1
501 22442879 : lnlb = nsoset(lbj)
502 22442879 : fsgfb = first_sgfb(jshellb)
503 22442879 : lsgfb = fsgfb + 2*lbj
504 68017046 : DO ishella = 1, nshella
505 23441943 : lai = la(ishella)
506 23441943 : fnla = nsoset(lai - 1) + 1
507 23441943 : lnla = nsoset(lai)
508 23441943 : fsgfa = first_sgfa(ishella)
509 23441943 : lsgfa = fsgfa + 2*lai
510 23441943 : labmin = MIN(lai, lbj)
511 106050579 : DO mbj = 0, 2*lbj
512 246222815 : DO mai = 0, 2*lai
513 550296192 : DO j = 0, labmin
514 327515320 : prefac = swork_cont(lai + lbj - j + 1, ishella, jshellb)
515 : sab(fsgfa + mai, fsgfb + mbj) = sab(fsgfa + mai, fsgfb + mbj) &
516 490130435 : + prefac*Waux_mat(j + 1, fnla + mai, fnlb + mbj)
517 : END DO
518 : END DO
519 : END DO
520 : END DO
521 : END DO
522 :
523 22132224 : END SUBROUTINE construct_int_shg_ab
524 :
525 : ! **************************************************************************************************
526 : !> \brief calculates derivatives of [ab] SHG overlap integrals using precomputed
527 : !> angular-dependent part
528 : !> \param la set of l quantum number on a
529 : !> \param first_sgfa indexing
530 : !> \param nshella number of shells for a
531 : !> \param lb set of l quantum number on b
532 : !> \param first_sgfb indexing
533 : !> \param nshellb number of shells for b
534 : !> \param rab distance vector Ra-Rb
535 : !> \param swork_cont contracted and normalized [s|s] integrals
536 : !> \param Waux_mat precomputed angular-dependent part
537 : !> \param dWaux_mat ...
538 : !> \param dsab derivative of contracted integral of spherical harmonic Gaussians
539 : ! **************************************************************************************************
540 116946 : SUBROUTINE construct_dev_shg_ab(la, first_sgfa, nshella, lb, first_sgfb, nshellb, rab, &
541 116946 : swork_cont, Waux_mat, dWaux_mat, dsab)
542 :
543 : INTEGER, DIMENSION(:), INTENT(IN) :: la, first_sgfa
544 : INTEGER, INTENT(IN) :: nshella
545 : INTEGER, DIMENSION(:), INTENT(IN) :: lb, first_sgfb
546 : INTEGER, INTENT(IN) :: nshellb
547 : REAL(KIND=dp), INTENT(IN) :: rab(3)
548 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: swork_cont, Waux_mat
549 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dWaux_mat
550 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: dsab
551 :
552 : INTEGER :: fnla, fnlb, fsgfa, fsgfb, i, ishella, j, &
553 : jshellb, labmin, lai, lbj, lnla, lnlb, &
554 : lsgfa, lsgfb
555 : REAL(KIND=dp) :: dprefac, prefac, rabx2(3)
556 :
557 467784 : rabx2(:) = 2.0_dp*rab
558 472162 : DO jshellb = 1, nshellb
559 355216 : lbj = lb(jshellb)
560 355216 : fnlb = nsoset(lbj - 1) + 1
561 355216 : lnlb = nsoset(lbj)
562 355216 : fsgfb = first_sgfb(jshellb)
563 355216 : lsgfb = fsgfb + 2*lbj
564 1589872 : DO ishella = 1, nshella
565 1117710 : lai = la(ishella)
566 1117710 : fnla = nsoset(lai - 1) + 1
567 1117710 : lnla = nsoset(lai)
568 1117710 : fsgfa = first_sgfa(ishella)
569 1117710 : lsgfa = fsgfa + 2*lai
570 1117710 : labmin = MIN(lai, lbj)
571 3300881 : DO j = 0, labmin
572 1827955 : prefac = swork_cont(lai + lbj - j + 1, ishella, jshellb)
573 1827955 : dprefac = swork_cont(lai + lbj - j + 2, ishella, jshellb) !j+1
574 8429530 : DO i = 1, 3
575 : dsab(fsgfa:lsgfa, fsgfb:lsgfb, i) = dsab(fsgfa:lsgfa, fsgfb:lsgfb, i) &
576 : + rabx2(i)*dprefac*Waux_mat(j + 1, fnla:lnla, fnlb:lnlb) &
577 125264344 : + prefac*dWaux_mat(i, j + 1, fnla:lnla, fnlb:lnlb)
578 : END DO
579 : END DO
580 : END DO
581 : END DO
582 :
583 116946 : END SUBROUTINE construct_dev_shg_ab
584 :
585 : ! **************************************************************************************************
586 : !> \brief calculates [aba] SHG overlap integrals using precomputed angular-
587 : !> dependent part
588 : !> \param la set of l quantum number on a, orbital basis
589 : !> \param first_sgfa indexing
590 : !> \param nshella number of shells for a, orbital basis
591 : !> \param lb set of l quantum number on b. orbital basis
592 : !> \param first_sgfb indexing
593 : !> \param nshellb number of shells for b, orbital basis
594 : !> \param lca of l quantum number on a, aux basis
595 : !> \param first_sgfca indexing
596 : !> \param nshellca number of shells for a, aux basis
597 : !> \param cg_coeff Clebsch-Gordon coefficients
598 : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
599 : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
600 : !> \param swork_cont contracted and normalized [s|ra^n|s] integrals
601 : !> \param Waux_mat precomputed angular-dependent part
602 : !> \param saba contracted overlap [aba] of spherical harmonic Gaussians
603 : ! **************************************************************************************************
604 70378 : SUBROUTINE construct_overlap_shg_aba(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
605 70378 : lca, first_sgfca, nshellca, cg_coeff, cg_none0_list, &
606 70378 : ncg_none0, swork_cont, Waux_mat, saba)
607 :
608 : INTEGER, DIMENSION(:), INTENT(IN) :: la, first_sgfa
609 : INTEGER, INTENT(IN) :: nshella
610 : INTEGER, DIMENSION(:), INTENT(IN) :: lb, first_sgfb
611 : INTEGER, INTENT(IN) :: nshellb
612 : INTEGER, DIMENSION(:), INTENT(IN) :: lca, first_sgfca
613 : INTEGER, INTENT(IN) :: nshellca
614 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff
615 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list
616 : INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0
617 : REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), &
618 : INTENT(IN) :: swork_cont
619 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Waux_mat
620 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: saba
621 :
622 : INTEGER :: ia, il, ilist, ishella, isoa1, isoa2, isoaa, j, jb, jshellb, ka, kshella, laa, &
623 : labmin, lai, lak, lbj, maa, mai, mak, mbj, nla, nlb, sgfa, sgfb, sgfca
624 : REAL(KIND=dp) :: prefac, stemp
625 :
626 247917 : DO kshella = 1, nshellca
627 177539 : lak = lca(kshella)
628 177539 : sgfca = first_sgfca(kshella)
629 177539 : ka = sgfca + lak
630 1075317 : DO jshellb = 1, nshellb
631 827400 : lbj = lb(jshellb)
632 827400 : nlb = nsoset(lbj - 1) + lbj + 1
633 827400 : sgfb = first_sgfb(jshellb)
634 827400 : jb = sgfb + lbj
635 5005859 : DO ishella = 1, nshella
636 4000920 : lai = la(ishella)
637 4000920 : sgfa = first_sgfa(ishella)
638 4000920 : ia = sgfa + lai
639 15087228 : DO mai = -lai, lai, 1
640 50009388 : DO mak = -lak, lak, 1
641 35749560 : isoa1 = indso_inv(lai, mai)
642 35749560 : isoa2 = indso_inv(lak, mak)
643 137109082 : DO mbj = -lbj, lbj, 1
644 312943791 : DO ilist = 1, ncg_none0(isoa1, isoa2)
645 186093617 : isoaa = cg_none0_list(isoa1, isoa2, ilist)
646 186093617 : laa = indso(1, isoaa)
647 186093617 : maa = indso(2, isoaa)
648 186093617 : nla = nsoset(laa - 1) + laa + 1
649 186093617 : labmin = MIN(laa, lbj)
650 186093617 : il = INT((lai + lak - laa)/2)
651 186093617 : stemp = 0.0_dp
652 576178399 : DO j = 0, labmin
653 390084782 : prefac = swork_cont(laa + lbj - j + 1, il, ishella, jshellb, kshella)
654 576178399 : stemp = stemp + prefac*Waux_mat(j + 1, nla + maa, nlb + mbj)
655 : END DO
656 277194231 : saba(ia + mai, jb + mbj, ka + mak) = saba(ia + mai, jb + mbj, ka + mak) + cg_coeff(isoa1, isoa2, isoaa)*stemp
657 : END DO
658 : END DO
659 : END DO
660 : END DO
661 : END DO
662 : END DO
663 : END DO
664 :
665 70378 : END SUBROUTINE construct_overlap_shg_aba
666 :
667 : ! **************************************************************************************************
668 : !> \brief calculates derivatives of [aba] SHG overlap integrals using
669 : !> precomputed angular-dependent part
670 : !> \param la set of l quantum number on a, orbital basis
671 : !> \param first_sgfa indexing
672 : !> \param nshella number of shells for a, orbital basis
673 : !> \param lb set of l quantum number on b. orbital basis
674 : !> \param first_sgfb indexing
675 : !> \param nshellb number of shells for b, orbital basis
676 : !> \param lca of l quantum number on a, aux basis
677 : !> \param first_sgfca indexing
678 : !> \param nshellca number of shells for a, aux basis
679 : !> \param cg_coeff Clebsch-Gordon coefficients
680 : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
681 : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
682 : !> \param rab distance vector Ra-Rb
683 : !> \param swork_cont contracted and normalized [s|ra^n|s] integrals
684 : !> \param Waux_mat precomputed angular-dependent part
685 : !> \param dWaux_mat derivatives of precomputed angular-dependent part
686 : !> \param dsaba derivative of contracted overlap [aba] of spherical harmonic
687 : !> Gaussians
688 : ! **************************************************************************************************
689 53612 : SUBROUTINE dev_overlap_shg_aba(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
690 107224 : lca, first_sgfca, nshellca, cg_coeff, cg_none0_list, &
691 53612 : ncg_none0, rab, swork_cont, Waux_mat, dWaux_mat, dsaba)
692 :
693 : INTEGER, DIMENSION(:), INTENT(IN) :: la, first_sgfa
694 : INTEGER, INTENT(IN) :: nshella
695 : INTEGER, DIMENSION(:), INTENT(IN) :: lb, first_sgfb
696 : INTEGER, INTENT(IN) :: nshellb
697 : INTEGER, DIMENSION(:), INTENT(IN) :: lca, first_sgfca
698 : INTEGER, INTENT(IN) :: nshellca
699 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff
700 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list
701 : INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0
702 : REAL(KIND=dp), INTENT(IN) :: rab(3)
703 : REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), &
704 : INTENT(IN) :: swork_cont
705 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Waux_mat
706 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dWaux_mat
707 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
708 : INTENT(INOUT) :: dsaba
709 :
710 : INTEGER :: i, ia, il, ilist, ishella, isoa1, isoa2, isoaa, j, jb, jshellb, ka, kshella, laa, &
711 : labmin, lai, lak, lbj, maa, mai, mak, mbj, nla, nlb, sgfa, sgfb, sgfca
712 : REAL(KIND=dp) :: dprefac, dtemp(3), prefac, rabx2(3)
713 :
714 214448 : rabx2(:) = 2.0_dp*rab
715 :
716 186925 : DO kshella = 1, nshellca
717 133313 : lak = lca(kshella)
718 133313 : sgfca = first_sgfca(kshella)
719 133313 : ka = sgfca + lak
720 843100 : DO jshellb = 1, nshellb
721 656175 : lbj = lb(jshellb)
722 656175 : nlb = nsoset(lbj - 1) + lbj + 1
723 656175 : sgfb = first_sgfb(jshellb)
724 656175 : jb = sgfb + lbj
725 4042313 : DO ishella = 1, nshella
726 3252825 : lai = la(ishella)
727 3252825 : sgfa = first_sgfa(ishella)
728 3252825 : ia = sgfa + lai
729 12346971 : DO mai = -lai, lai, 1
730 40565761 : DO mak = -lak, lak, 1
731 28874965 : isoa1 = indso_inv(lai, mai)
732 28874965 : isoa2 = indso_inv(lak, mak)
733 112211889 : DO mbj = -lbj, lbj, 1
734 255925300 : DO ilist = 1, ncg_none0(isoa1, isoa2)
735 152151382 : isoaa = cg_none0_list(isoa1, isoa2, ilist)
736 152151382 : laa = indso(1, isoaa)
737 152151382 : maa = indso(2, isoaa)
738 152151382 : nla = nsoset(laa - 1) + laa + 1
739 152151382 : labmin = MIN(laa, lbj)
740 152151382 : il = (lai + lak - laa)/2 ! lai+lak-laa always even
741 152151382 : dtemp = 0.0_dp
742 473206046 : DO j = 0, labmin
743 321054664 : prefac = swork_cont(laa + lbj - j + 1, il, ishella, jshellb, kshella)
744 321054664 : dprefac = swork_cont(laa + lbj - j + 2, il, ishella, jshellb, kshella)
745 1436370038 : DO i = 1, 3
746 : dtemp(i) = dtemp(i) + rabx2(i)*dprefac*Waux_mat(j + 1, nla + maa, nlb + mbj) &
747 1284218656 : + prefac*dWaux_mat(i, j + 1, nla + maa, nlb + mbj)
748 : END DO
749 : END DO
750 683504481 : DO i = 1, 3
751 : dsaba(ia + mai, jb + mbj, ka + mak, i) = dsaba(ia + mai, jb + mbj, ka + mak, i) &
752 608605528 : + cg_coeff(isoa1, isoa2, isoaa)*dtemp(i)
753 : END DO
754 : END DO
755 : END DO
756 : END DO
757 : END DO
758 : END DO
759 : END DO
760 : END DO
761 :
762 53612 : END SUBROUTINE dev_overlap_shg_aba
763 :
764 : ! **************************************************************************************************
765 : !> \brief calculates [abb] SHG overlap integrals using precomputed angular-
766 : !> dependent part
767 : !> \param la set of l quantum number on a, orbital basis
768 : !> \param first_sgfa indexing
769 : !> \param nshella number of shells for a, orbital basis
770 : !> \param lb set of l quantum number on b. orbital basis
771 : !> \param first_sgfb indexing
772 : !> \param nshellb number of shells for b, orbital basis
773 : !> \param lcb l quantum number on b, aux basis
774 : !> \param first_sgfcb indexing
775 : !> \param nshellcb number of shells for b, aux basis
776 : !> \param cg_coeff Clebsch-Gordon coefficients
777 : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
778 : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
779 : !> \param swork_cont contracted and normalized [s|rb^n|s] integrals
780 : !> \param Waux_mat precomputed angular-dependent part
781 : !> \param sabb contracted overlap [abb] of spherical harmonic Gaussians
782 : ! **************************************************************************************************
783 3421 : SUBROUTINE construct_overlap_shg_abb(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
784 3421 : lcb, first_sgfcb, nshellcb, cg_coeff, cg_none0_list, &
785 3421 : ncg_none0, swork_cont, Waux_mat, sabb)
786 :
787 : INTEGER, DIMENSION(:), INTENT(IN) :: la, first_sgfa
788 : INTEGER, INTENT(IN) :: nshella
789 : INTEGER, DIMENSION(:), INTENT(IN) :: lb, first_sgfb
790 : INTEGER, INTENT(IN) :: nshellb
791 : INTEGER, DIMENSION(:), INTENT(IN) :: lcb, first_sgfcb
792 : INTEGER, INTENT(IN) :: nshellcb
793 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff
794 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list
795 : INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0
796 : REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), &
797 : INTENT(IN) :: swork_cont
798 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Waux_mat
799 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: sabb
800 :
801 : INTEGER :: ia, il, ilist, ishella, isob1, isob2, isobb, j, jb, jshellb, kb, kshellb, labmin, &
802 : lai, lbb, lbj, lbk, mai, mbb, mbj, mbk, nla, nlb, sgfa, sgfb, sgfcb
803 : REAL(KIND=dp) :: prefac, stemp, tsign
804 :
805 15364 : DO kshellb = 1, nshellcb
806 11943 : lbk = lcb(kshellb)
807 11943 : sgfcb = first_sgfcb(kshellb)
808 11943 : kb = sgfcb + lbk
809 64033 : DO jshellb = 1, nshellb
810 48669 : lbj = lb(jshellb)
811 48669 : sgfb = first_sgfb(jshellb)
812 48669 : jb = sgfb + lbj
813 217476 : DO ishella = 1, nshella
814 156864 : lai = la(ishella)
815 156864 : nla = nsoset(lai - 1) + lai + 1
816 156864 : sgfa = first_sgfa(ishella)
817 156864 : ia = sgfa + lai
818 579495 : DO mbj = -lbj, lbj, 1
819 2206050 : DO mbk = -lbk, lbk, 1
820 1675224 : isob1 = indso_inv(lbj, mbj)
821 1675224 : isob2 = indso_inv(lbk, mbk)
822 5056410 : DO mai = -lai, lai, 1
823 11504847 : DO ilist = 1, ncg_none0(isob1, isob2)
824 6822399 : isobb = cg_none0_list(isob1, isob2, ilist)
825 6822399 : lbb = indso(1, isobb)
826 6822399 : mbb = indso(2, isobb)
827 6822399 : nlb = nsoset(lbb - 1) + lbb + 1
828 : ! tsgin: because we take the transpose of auxmat (calculated for (la,lb))
829 6822399 : tsign = 1.0_dp
830 6822399 : IF (MODULO(lbb - lai, 2) /= 0) tsign = -1.0_dp
831 6822399 : labmin = MIN(lai, lbb)
832 6822399 : il = INT((lbj + lbk - lbb)/2)
833 6822399 : stemp = 0.0_dp
834 18422882 : DO j = 0, labmin
835 11600483 : prefac = swork_cont(lai + lbb - j + 1, il, ishella, jshellb, kshellb)
836 18422882 : stemp = stemp + prefac*Waux_mat(j + 1, nlb + mbb, nla + mai)
837 : END DO
838 9829623 : sabb(ia + mai, jb + mbj, kb + mbk) = sabb(ia + mai, jb + mbj, kb + mbk) + tsign*cg_coeff(isob1, isob2, isobb)*stemp
839 : END DO
840 : END DO
841 : END DO
842 : END DO
843 : END DO
844 : END DO
845 : END DO
846 :
847 3421 : END SUBROUTINE construct_overlap_shg_abb
848 :
849 : ! **************************************************************************************************
850 : !> \brief calculates derivatives of [abb] SHG overlap integrals using
851 : !> precomputed angular-dependent part
852 : !> \param la set of l quantum number on a, orbital basis
853 : !> \param first_sgfa indexing
854 : !> \param nshella number of shells for a, orbital basis
855 : !> \param lb set of l quantum number on b. orbital basis
856 : !> \param first_sgfb indexing
857 : !> \param nshellb number of shells for b, orbital basis
858 : !> \param lcb l quantum number on b, aux basis
859 : !> \param first_sgfcb indexing
860 : !> \param nshellcb number of shells for b, aux basis
861 : !> \param cg_coeff Clebsch-Gordon coefficients
862 : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
863 : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
864 : !> \param rab distance vector Ra-Rb
865 : !> \param swork_cont contracted and normalized [s|rb^n|s] integrals
866 : !> \param Waux_mat precomputed angular-dependent part
867 : !> \param dWaux_mat derivatives of precomputed angular-dependent part
868 : !> \param dsabb derivative of contracted overlap [abb] of spherical harmonic
869 : !> Gaussians
870 : ! **************************************************************************************************
871 1111 : SUBROUTINE dev_overlap_shg_abb(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
872 2222 : lcb, first_sgfcb, nshellcb, cg_coeff, cg_none0_list, &
873 1111 : ncg_none0, rab, swork_cont, Waux_mat, dWaux_mat, dsabb)
874 :
875 : INTEGER, DIMENSION(:), INTENT(IN) :: la, first_sgfa
876 : INTEGER, INTENT(IN) :: nshella
877 : INTEGER, DIMENSION(:), INTENT(IN) :: lb, first_sgfb
878 : INTEGER, INTENT(IN) :: nshellb
879 : INTEGER, DIMENSION(:), INTENT(IN) :: lcb, first_sgfcb
880 : INTEGER, INTENT(IN) :: nshellcb
881 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff
882 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list
883 : INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0
884 : REAL(KIND=dp), INTENT(IN) :: rab(3)
885 : REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), &
886 : INTENT(IN) :: swork_cont
887 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Waux_mat
888 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dWaux_mat
889 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
890 : INTENT(INOUT) :: dsabb
891 :
892 : INTEGER :: i, ia, il, ilist, ishella, isob1, isob2, isobb, j, jb, jshellb, kb, kshellb, &
893 : labmin, lai, lbb, lbj, lbk, mai, mbb, mbj, mbk, nla, nlb, sgfa, sgfb, sgfcb
894 : REAL(KIND=dp) :: dprefac, dtemp(3), prefac, rabx2(3), &
895 : tsign
896 :
897 4444 : rabx2(:) = 2.0_dp*rab
898 :
899 4206 : DO kshellb = 1, nshellcb
900 3095 : lbk = lcb(kshellb)
901 3095 : sgfcb = first_sgfcb(kshellb)
902 3095 : kb = sgfcb + lbk
903 13951 : DO jshellb = 1, nshellb
904 9745 : lbj = lb(jshellb)
905 9745 : sgfb = first_sgfb(jshellb)
906 9745 : jb = sgfb + lbj
907 40263 : DO ishella = 1, nshella
908 27423 : lai = la(ishella)
909 27423 : nla = nsoset(lai - 1) + lai + 1
910 27423 : sgfa = first_sgfa(ishella)
911 27423 : ia = sgfa + lai
912 106345 : DO mbj = -lbj, lbj, 1
913 400259 : DO mbk = -lbk, lbk, 1
914 303659 : isob1 = indso_inv(lbj, mbj)
915 303659 : isob2 = indso_inv(lbk, mbk)
916 898811 : DO mai = -lai, lai, 1
917 2168821 : DO ilist = 1, ncg_none0(isob1, isob2)
918 1339187 : isobb = cg_none0_list(isob1, isob2, ilist)
919 1339187 : lbb = indso(1, isobb)
920 1339187 : mbb = indso(2, isobb)
921 1339187 : nlb = nsoset(lbb - 1) + lbb + 1
922 : ! tsgin: because we take the transpose of auxmat (calculated for (la,lb))
923 1339187 : tsign = 1.0_dp
924 1339187 : IF (MODULO(lbb - lai, 2) /= 0) tsign = -1.0_dp
925 1339187 : labmin = MIN(lai, lbb)
926 1339187 : il = (lbj + lbk - lbb)/2
927 1339187 : dtemp = 0.0_dp
928 3768153 : DO j = 0, labmin
929 2428966 : prefac = swork_cont(lai + lbb - j + 1, il, ishella, jshellb, kshellb)
930 2428966 : dprefac = swork_cont(lai + lbb - j + 2, il, ishella, jshellb, kshellb)
931 11055051 : DO i = 1, 3
932 : dtemp(i) = dtemp(i) + rabx2(i)*dprefac*Waux_mat(j + 1, nlb + mbb, nla + mai) &
933 9715864 : + prefac*dWaux_mat(i, j + 1, nlb + mbb, nla + mai)
934 : END DO
935 : END DO
936 5882723 : DO i = 1, 3
937 : dsabb(ia + mai, jb + mbj, kb + mbk, i) = dsabb(ia + mai, jb + mbj, kb + mbk, i) &
938 5356748 : + tsign*cg_coeff(isob1, isob2, isobb)*dtemp(i)
939 : END DO
940 : END DO
941 : END DO
942 : END DO
943 : END DO
944 : END DO
945 : END DO
946 : END DO
947 :
948 1111 : END SUBROUTINE dev_overlap_shg_abb
949 :
950 : END MODULE construct_shg
951 :
|
