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 7221708 : 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 7221708 : Rc_00 = 1.0_dp
62 7221708 : Rs_00 = 0.0_dp
63 :
64 7221708 : Rlm_c(0, 0) = Rc_00
65 7221708 : Rlm_s(0, 0) = Rs_00
66 :
67 : ! generate elements Rmm
68 : ! start
69 7221708 : IF (l > 0) THEN
70 4140248 : Rc = -0.5_dp*r(1)*Rc_00
71 4140248 : Rs = -0.5_dp*r(2)*Rc_00
72 4140248 : Rlm_c(1, 1) = Rc
73 4140248 : Rlm_s(1, 1) = Rs
74 4140248 : Rlm_c(1, -1) = -Rc
75 4140248 : Rlm_s(1, -1) = Rs
76 : END IF
77 7784356 : DO li = 2, l
78 562648 : temp_c = (-r(1)*Rc + r(2)*Rs)/(REAL(2*(li - 1) + 2, dp))
79 562648 : Rs = (-r(2)*Rc - r(1)*Rs)/(REAL(2*(li - 1) + 2, dp))
80 562648 : Rc = temp_c
81 562648 : Rlm_c(li, li) = Rc
82 562648 : Rlm_s(li, li) = Rs
83 7784356 : IF (MODULO(li, 2) /= 0) THEN
84 148674 : Rlm_c(li, -li) = -Rc
85 148674 : Rlm_s(li, -li) = Rs
86 : ELSE
87 413974 : Rlm_c(li, -li) = Rc
88 413974 : Rlm_s(li, -li) = -Rs
89 : END IF
90 : END DO
91 :
92 11924604 : DO mi = 0, l - 1
93 4702896 : Rmlm = Rlm_c(mi, mi)
94 4702896 : Rlm = r(3)*Rlm_c(mi, mi)
95 4702896 : Rlm_c(mi + 1, mi) = Rlm
96 4702896 : IF (MODULO(mi, 2) /= 0) THEN
97 413974 : Rlm_c(mi + 1, -mi) = -Rlm
98 : ELSE
99 4288922 : Rlm_c(mi + 1, -mi) = Rlm
100 : END IF
101 12760590 : DO li = mi + 2, l
102 835986 : prefac = (li + mi)*(li - mi)
103 835986 : Rplm = (REAL(2*li - 1, dp)*r(3)*Rlm - r2*Rmlm)/REAL(prefac, dp)
104 835986 : Rmlm = Rlm
105 835986 : Rlm = Rplm
106 835986 : Rlm_c(li, mi) = Rlm
107 5538882 : IF (MODULO(mi, 2) /= 0) THEN
108 211006 : Rlm_c(li, -mi) = -Rlm
109 : ELSE
110 624980 : Rlm_c(li, -mi) = Rlm
111 : END IF
112 : END DO
113 : END DO
114 7784356 : DO mi = 1, l - 1
115 562648 : Rmlm = Rlm_s(mi, mi)
116 562648 : Rlm = r(3)*Rlm_s(mi, mi)
117 562648 : Rlm_s(mi + 1, mi) = Rlm
118 562648 : IF (MODULO(mi, 2) /= 0) THEN
119 413974 : Rlm_s(mi + 1, -mi) = Rlm
120 : ELSE
121 148674 : Rlm_s(mi + 1, -mi) = -Rlm
122 : END IF
123 8057694 : 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 835986 : 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 7221708 : 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 7221708 : 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 19146312 : DO l = 0, lmax
154 36609798 : DO m = 0, l
155 17463486 : temp = SQRT(fac(l + m)*fac(l - m))
156 17463486 : IF (MODULO(m, 2) /= 0) temp = -temp
157 17463486 : IF (m /= 0) temp = temp*SQRT(2.0_dp)
158 29388090 : A(l, m) = temp
159 : END DO
160 : END DO
161 :
162 7221708 : 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 7221708 : 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 7221708 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: A
223 :
224 7221708 : Wa(:) = 0.0_dp
225 7221708 : Wb(:) = 0.0_dp
226 7221708 : Wmat(:) = 0.0_dp
227 :
228 28886832 : ALLOCATE (A(0:lmax, 0:lmax))
229 7221708 : CALL get_Alm(lmax, A)
230 :
231 17419011 : DO lb = 0, lbmax
232 10197303 : nlb = nsoset(lb - 1)
233 32896860 : DO la = 0, lamax(lb)
234 15477849 : nla = nsoset(la - 1)
235 15477849 : labmin = MIN(la, lb)
236 47955435 : DO mb = 0, lb
237 22280283 : A_lbmb = A(lb, mb)
238 22280283 : IF (MODULO(lb, 2) /= 0) A_lbmb = -A_lbmb
239 73342017 : DO ma = 0, la
240 35583885 : A_lama = A(la, ma)
241 35583885 : Alm_fac = A_lama*A_lbmb
242 116059985 : DO j = 0, labmin
243 58195817 : laj = la - j
244 58195817 : lbj = lb - j
245 58195817 : prefac = Alm_fac*REAL(2**(la + lb - j), dp)*dfac(2*j - 1)
246 58195817 : delta_k = 0.5_dp
247 58195817 : Wmat = 0.0_dp
248 147921621 : DO k = 0, j
249 89725804 : ma_m = ma - k
250 89725804 : ma_p = ma + k
251 89725804 : IF (laj < ABS(ma_m) .AND. laj < ABS(ma_p)) CYCLE
252 68992740 : mb_m = mb - k
253 68992740 : mb_p = mb + k
254 68992740 : IF (lbj < ABS(mb_m) .AND. lbj < ABS(mb_p)) CYCLE
255 56306184 : IF (k /= 0) delta_k = 1.0_dp
256 56306184 : A_jk = fac(j + k)*fac(j - k)
257 56306184 : IF (k /= 0) A_jk = 2.0_dp*A_jk
258 13628967 : IF (MODULO(k, 2) /= 0) THEN
259 : sign_fac = -1.0_dp
260 : ELSE
261 : sign_fac = 1.0_dp
262 : END IF
263 56306184 : Rca_m = Rc(laj, ma_m)
264 56306184 : Rsa_m = Rs(laj, ma_m)
265 56306184 : Rca_p = Rc(laj, ma_p)
266 56306184 : Rsa_p = Rs(laj, ma_p)
267 56306184 : Rcb_m = Rc(lbj, mb_m)
268 56306184 : Rsb_m = Rs(lbj, mb_m)
269 56306184 : Rcb_p = Rc(lbj, mb_p)
270 56306184 : Rsb_p = Rs(lbj, mb_p)
271 56306184 : Wa(1) = delta_k*(Rca_m + sign_fac*Rca_p)
272 56306184 : Wb(1) = delta_k*(Rcb_m + sign_fac*Rcb_p)
273 56306184 : Wa(2) = -Rsa_m + sign_fac*Rsa_p
274 56306184 : Wb(2) = -Rsb_m + sign_fac*Rsb_p
275 56306184 : Wmat(1) = Wmat(1) + prefac/A_jk*(Wa(1)*Wb(1) + Wa(2)*Wb(2))
276 56306184 : IF (mb > 0) THEN
277 26232943 : Wb(3) = delta_k*(Rsb_m + sign_fac*Rsb_p)
278 26232943 : Wb(4) = Rcb_m - sign_fac*Rcb_p
279 26232943 : Wmat(2) = Wmat(2) + prefac/A_jk*(Wa(1)*Wb(3) + Wa(2)*Wb(4))
280 : END IF
281 56306184 : IF (ma > 0) THEN
282 27373292 : Wa(3) = delta_k*(Rsa_m + sign_fac*Rsa_p)
283 27373292 : Wa(4) = Rca_m - sign_fac*Rca_p
284 27373292 : Wmat(3) = Wmat(3) + prefac/A_jk*(Wa(3)*Wb(1) + Wa(4)*Wb(2))
285 : END IF
286 114502001 : IF (ma > 0 .AND. mb > 0) THEN
287 16198070 : Wmat(4) = Wmat(4) + prefac/A_jk*(Wa(3)*Wb(3) + Wa(4)*Wb(4))
288 : END IF
289 : END DO
290 58195817 : Waux_mat(j + 1, nla + la + 1 + ma, nlb + lb + 1 + mb) = Wmat(1)
291 58195817 : IF (mb > 0) Waux_mat(j + 1, nla + la + 1 + ma, nlb + lb + 1 - mb) = Wmat(2)
292 58195817 : IF (ma > 0) Waux_mat(j + 1, nla + la + 1 - ma, nlb + lb + 1 + mb) = Wmat(3)
293 93779702 : 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 7221708 : DEALLOCATE (A)
301 :
302 7221708 : 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 2182520 : Wam(:, :) = 0.0_dp
367 2182520 : Wamm(:, :) = 0.0_dp
368 2182520 : 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 : ! order: e.g. -1 0 1, if < 0 |m|, -1 means -m+1
381 3946483 : IF (ma - 1 > 0) Wamm(0:lamb, 3) = Waux_mat(1:lamb + 1, imam + 1, ipb)
382 3039982 : IF (ma - 1 > 0 .AND. mb > 0) Wamm(0:lamb, 4) = Waux_mat(1:lamb + 1, imam + 1, imb)
383 : END IF
384 : !*** Wamp: la-1, ma+1
385 2182520 : IF (ma + 1 <= la - 1) THEN
386 3407029 : Wamp(0:lamb, 1) = Waux_mat(1:lamb + 1, ipam + 1, ipb)
387 2500528 : IF (mb > 0) Wamp(0:lamb, 2) = Waux_mat(1:lamb + 1, ipam + 1, imb)
388 3407029 : IF (ma + 1 > 0) Wamp(0:lamb, 3) = Waux_mat(1:lamb + 1, imam - 1, ipb)
389 2500528 : IF (ma + 1 > 0 .AND. mb > 0) Wamp(0:lamb, 4) = Waux_mat(1:lamb + 1, imam - 1, imb)
390 : END IF
391 2182520 : Wbm(:, :) = 0.0_dp
392 2182520 : Wbmm(:, :) = 0.0_dp
393 2182520 : Wbmp(:, :) = 0.0_dp
394 : !*** Wbm: lb-1, mb
395 2182520 : IF (mb <= lb - 1) THEN
396 3855304 : Wbm(0:labm, 1) = Waux_mat(1:labm + 1, ipa, ipbm)
397 2675803 : IF (mb > 0) Wbm(0:labm, 2) = Waux_mat(1:labm + 1, ipa, imbm)
398 3108049 : IF (ma > 0) Wbm(0:labm, 3) = Waux_mat(1:labm + 1, ima, ipbm)
399 2264313 : IF (ma > 0 .AND. mb > 0) Wbm(0:labm, 4) = Waux_mat(1:labm + 1, ima, imbm)
400 : END IF
401 : !*** Wbmm: lb-1, mb-1
402 2182520 : IF (mb - 1 >= 0) THEN
403 3855304 : Wbmm(0:labm, 1) = Waux_mat(1:labm + 1, ipa, ipbm - 1)
404 2675803 : IF (mb - 1 > 0) Wbmm(0:labm, 2) = Waux_mat(1:labm + 1, ipa, imbm + 1)
405 3108049 : IF (ma > 0) Wbmm(0:labm, 3) = Waux_mat(1:labm + 1, ima, ipbm - 1)
406 2264313 : IF (ma > 0 .AND. mb - 1 > 0) Wbmm(0:labm, 4) = Waux_mat(1:labm + 1, ima, imbm + 1)
407 : END IF
408 : !*** Wbmp: lb-1, mb+1
409 2182520 : IF (mb + 1 <= lb - 1) THEN
410 2006040 : Wbmp(0:labm, 1) = Waux_mat(1:labm + 1, ipa, ipbm + 1)
411 2006040 : IF (mb + 1 > 0) Wbmp(0:labm, 2) = Waux_mat(1:labm + 1, ipa, imbm - 1)
412 1594550 : IF (ma > 0) Wbmp(0:labm, 3) = Waux_mat(1:labm + 1, ima, ipbm + 1)
413 1594550 : IF (ma > 0 .AND. mb + 1 > 0) Wbmp(0:labm, 4) = Waux_mat(1:labm + 1, ima, imbm - 1)
414 : END IF
415 8334934 : DO j = 0, labmin
416 : !*** x component
417 : dWaux_mat(1, j + 1, ipa, ipb) = dAa_p*Wamp(j, 1) - dAa_m*Wamm(j, 1) &
418 5487042 : - dAb_p*Wbmp(j, 1) + dAb_m*Wbmm(j, 1)
419 5487042 : IF (mb > 0) THEN
420 : dWaux_mat(1, j + 1, ipa, imb) = dAa_p*Wamp(j, 2) - dAa_m*Wamm(j, 2) &
421 3507749 : - dAb_p*Wbmp(j, 2) + dAb_m*Wbmm(j, 2)
422 : END IF
423 5487042 : IF (ma > 0) THEN
424 : dWaux_mat(1, j + 1, ima, ipb) = dAa_p*Wamp(j, 3) - dAa_m*Wamm(j, 3) &
425 4000776 : - dAb_p*Wbmp(j, 3) + dAb_m*Wbmm(j, 3)
426 : END IF
427 5487042 : IF (ma > 0 .AND. mb > 0) THEN
428 : dWaux_mat(1, j + 1, ima, imb) = dAa_p*Wamp(j, 4) - dAa_m*Wamm(j, 4) &
429 2555792 : - dAb_p*Wbmp(j, 4) + dAb_m*Wbmm(j, 4)
430 : END IF
431 :
432 : !**** y component
433 : dWaux_mat(2, j + 1, ipa, ipb) = dAa_p*Wamp(j, 3) + dAa_m*Wamm(j, 3) &
434 5487042 : - dAb_p*Wbmp(j, 2) - dAb_m*Wbmm(j, 2)
435 5487042 : IF (mb > 0) THEN
436 : dWaux_mat(2, j + 1, ipa, imb) = dAa_p*Wamp(j, 4) + dAa_m*Wamm(j, 4) &
437 3507749 : + dAb_p*Wbmp(j, 1) + dAb_m*Wbmm(j, 1)
438 : END IF
439 5487042 : IF (ma > 0) THEN
440 : dWaux_mat(2, j + 1, ima, ipb) = -dAa_p*Wamp(j, 1) - dAa_m*Wamm(j, 1) &
441 4000776 : - dAb_p*Wbmp(j, 4) - dAb_m*Wbmm(j, 4)
442 : END IF
443 5487042 : IF (ma > 0 .AND. mb > 0) THEN
444 : dWaux_mat(2, j + 1, ima, imb) = -dAa_p*Wamp(j, 2) - dAa_m*Wamm(j, 2) &
445 2555792 : + dAb_p*Wbmp(j, 3) + dAb_m*Wbmm(j, 3)
446 : END IF
447 : !**** z compnent
448 5487042 : dWaux_mat(3, j + 1, ipa, ipb) = dAa*Wam(j, 1) - dAb*Wbm(j, 1)
449 5487042 : IF (mb > 0) THEN
450 3507749 : dWaux_mat(3, j + 1, ipa, imb) = dAa*Wam(j, 2) - dAb*Wbm(j, 2)
451 : END IF
452 5487042 : IF (ma > 0) THEN
453 4000776 : dWaux_mat(3, j + 1, ima, ipb) = dAa*Wam(j, 3) - dAb*Wbm(j, 3)
454 : END IF
455 7669562 : IF (ma > 0 .AND. mb > 0) THEN
456 2555792 : dWaux_mat(3, j + 1, ima, imb) = dAa*Wam(j, 4) - dAb*Wbm(j, 4)
457 : END IF
458 :
459 : END DO
460 : END DO
461 : END DO
462 : END DO
463 : END DO
464 :
465 12204 : DEALLOCATE (Wam, Wamm, Wamp)
466 12204 : DEALLOCATE (Wbm, Wbmm, Wbmp)
467 12204 : DEALLOCATE (dA, dA_p, dA_m)
468 :
469 12204 : END SUBROUTINE get_dW_matrix
470 :
471 : ! **************************************************************************************************
472 : !> \brief calculates [ab] SHG overlap integrals using precomputed angular-
473 : !> dependent part
474 : !> \param la set of l quantum number on a
475 : !> \param first_sgfa indexing
476 : !> \param nshella number of shells for a
477 : !> \param lb set of l quantum number on b
478 : !> \param first_sgfb indexing
479 : !> \param nshellb number of shells for b
480 : !> \param swork_cont contracted and normalized [s|s] integrals
481 : !> \param Waux_mat precomputed angular-dependent part
482 : !> \param sab contracted integral of spherical harmonic Gaussianslm
483 : ! **************************************************************************************************
484 28768449 : SUBROUTINE construct_int_shg_ab(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
485 28768449 : swork_cont, Waux_mat, sab)
486 :
487 : INTEGER, DIMENSION(:), INTENT(IN) :: la, first_sgfa
488 : INTEGER, INTENT(IN) :: nshella
489 : INTEGER, DIMENSION(:), INTENT(IN) :: lb, first_sgfb
490 : INTEGER, INTENT(IN) :: nshellb
491 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: swork_cont, Waux_mat
492 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
493 :
494 : INTEGER :: fnla, fnlb, fsgfa, fsgfb, ishella, j, &
495 : jshellb, labmin, lai, lbj, lnla, lnlb, &
496 : lsgfa, lsgfb, mai, mbj
497 : REAL(KIND=dp) :: prefac
498 :
499 57847553 : DO jshellb = 1, nshellb
500 29079104 : lbj = lb(jshellb)
501 29079104 : fnlb = nsoset(lbj - 1) + 1
502 29079104 : lnlb = nsoset(lbj)
503 29079104 : fsgfb = first_sgfb(jshellb)
504 29079104 : lsgfb = fsgfb + 2*lbj
505 87925721 : DO ishella = 1, nshella
506 30078168 : lai = la(ishella)
507 30078168 : fnla = nsoset(lai - 1) + 1
508 30078168 : lnla = nsoset(lai)
509 30078168 : fsgfa = first_sgfa(ishella)
510 30078168 : lsgfa = fsgfa + 2*lai
511 30078168 : labmin = MIN(lai, lbj)
512 129277684 : DO mbj = 0, 2*lbj
513 277746304 : DO mai = 0, 2*lai
514 593850006 : DO j = 0, labmin
515 346181870 : prefac = swork_cont(lai + lbj - j + 1, ishella, jshellb)
516 : sab(fsgfa + mai, fsgfb + mbj) = sab(fsgfa + mai, fsgfb + mbj) &
517 523729594 : + prefac*Waux_mat(j + 1, fnla + mai, fnlb + mbj)
518 : END DO
519 : END DO
520 : END DO
521 : END DO
522 : END DO
523 :
524 28768449 : END SUBROUTINE construct_int_shg_ab
525 :
526 : ! **************************************************************************************************
527 : !> \brief calculates derivatives of [ab] SHG overlap integrals using precomputed
528 : !> angular-dependent part
529 : !> \param la set of l quantum number on a
530 : !> \param first_sgfa indexing
531 : !> \param nshella number of shells for a
532 : !> \param lb set of l quantum number on b
533 : !> \param first_sgfb indexing
534 : !> \param nshellb number of shells for b
535 : !> \param rab distance vector Ra-Rb
536 : !> \param swork_cont contracted and normalized [s|s] integrals
537 : !> \param Waux_mat precomputed angular-dependent part
538 : !> \param dWaux_mat ...
539 : !> \param dsab derivative of contracted integral of spherical harmonic Gaussians
540 : ! **************************************************************************************************
541 116946 : SUBROUTINE construct_dev_shg_ab(la, first_sgfa, nshella, lb, first_sgfb, nshellb, rab, &
542 116946 : swork_cont, Waux_mat, dWaux_mat, dsab)
543 :
544 : INTEGER, DIMENSION(:), INTENT(IN) :: la, first_sgfa
545 : INTEGER, INTENT(IN) :: nshella
546 : INTEGER, DIMENSION(:), INTENT(IN) :: lb, first_sgfb
547 : INTEGER, INTENT(IN) :: nshellb
548 : REAL(KIND=dp), INTENT(IN) :: rab(3)
549 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: swork_cont, Waux_mat
550 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dWaux_mat
551 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: dsab
552 :
553 : INTEGER :: fnla, fnlb, fsgfa, fsgfb, i, ishella, j, &
554 : jshellb, labmin, lai, lbj, lnla, lnlb, &
555 : lsgfa, lsgfb
556 : REAL(KIND=dp) :: dprefac, prefac, rabx2(3)
557 :
558 467784 : rabx2(:) = 2.0_dp*rab
559 472162 : DO jshellb = 1, nshellb
560 355216 : lbj = lb(jshellb)
561 355216 : fnlb = nsoset(lbj - 1) + 1
562 355216 : lnlb = nsoset(lbj)
563 355216 : fsgfb = first_sgfb(jshellb)
564 355216 : lsgfb = fsgfb + 2*lbj
565 1589872 : DO ishella = 1, nshella
566 1117710 : lai = la(ishella)
567 1117710 : fnla = nsoset(lai - 1) + 1
568 1117710 : lnla = nsoset(lai)
569 1117710 : fsgfa = first_sgfa(ishella)
570 1117710 : lsgfa = fsgfa + 2*lai
571 1117710 : labmin = MIN(lai, lbj)
572 3300881 : DO j = 0, labmin
573 1827955 : prefac = swork_cont(lai + lbj - j + 1, ishella, jshellb)
574 1827955 : dprefac = swork_cont(lai + lbj - j + 2, ishella, jshellb) !j+1
575 8429530 : DO i = 1, 3
576 : dsab(fsgfa:lsgfa, fsgfb:lsgfb, i) = dsab(fsgfa:lsgfa, fsgfb:lsgfb, i) &
577 : + rabx2(i)*dprefac*Waux_mat(j + 1, fnla:lnla, fnlb:lnlb) &
578 125264344 : + prefac*dWaux_mat(i, j + 1, fnla:lnla, fnlb:lnlb)
579 : END DO
580 : END DO
581 : END DO
582 : END DO
583 :
584 116946 : END SUBROUTINE construct_dev_shg_ab
585 :
586 : ! **************************************************************************************************
587 : !> \brief calculates [aba] SHG overlap integrals using precomputed angular-
588 : !> dependent part
589 : !> \param la set of l quantum number on a, orbital basis
590 : !> \param first_sgfa indexing
591 : !> \param nshella number of shells for a, orbital basis
592 : !> \param lb set of l quantum number on b. orbital basis
593 : !> \param first_sgfb indexing
594 : !> \param nshellb number of shells for b, orbital basis
595 : !> \param lca of l quantum number on a, aux basis
596 : !> \param first_sgfca indexing
597 : !> \param nshellca number of shells for a, aux basis
598 : !> \param cg_coeff Clebsch-Gordon coefficients
599 : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
600 : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
601 : !> \param swork_cont contracted and normalized [s|ra^n|s] integrals
602 : !> \param Waux_mat precomputed angular-dependent part
603 : !> \param saba contracted overlap [aba] of spherical harmonic Gaussians
604 : ! **************************************************************************************************
605 70378 : SUBROUTINE construct_overlap_shg_aba(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
606 70378 : lca, first_sgfca, nshellca, cg_coeff, cg_none0_list, &
607 70378 : ncg_none0, swork_cont, Waux_mat, saba)
608 :
609 : INTEGER, DIMENSION(:), INTENT(IN) :: la, first_sgfa
610 : INTEGER, INTENT(IN) :: nshella
611 : INTEGER, DIMENSION(:), INTENT(IN) :: lb, first_sgfb
612 : INTEGER, INTENT(IN) :: nshellb
613 : INTEGER, DIMENSION(:), INTENT(IN) :: lca, first_sgfca
614 : INTEGER, INTENT(IN) :: nshellca
615 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff
616 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list
617 : INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0
618 : REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), &
619 : INTENT(IN) :: swork_cont
620 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Waux_mat
621 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: saba
622 :
623 : INTEGER :: ia, il, ilist, ishella, isoa1, isoa2, isoaa, j, jb, jshellb, ka, kshella, laa, &
624 : labmin, lai, lak, lbj, maa, mai, mak, mbj, nla, nlb, sgfa, sgfb, sgfca
625 : REAL(KIND=dp) :: prefac, stemp
626 :
627 247917 : DO kshella = 1, nshellca
628 177539 : lak = lca(kshella)
629 177539 : sgfca = first_sgfca(kshella)
630 177539 : ka = sgfca + lak
631 1075317 : DO jshellb = 1, nshellb
632 827400 : lbj = lb(jshellb)
633 827400 : nlb = nsoset(lbj - 1) + lbj + 1
634 827400 : sgfb = first_sgfb(jshellb)
635 827400 : jb = sgfb + lbj
636 5005859 : DO ishella = 1, nshella
637 4000920 : lai = la(ishella)
638 4000920 : sgfa = first_sgfa(ishella)
639 4000920 : ia = sgfa + lai
640 15087228 : DO mai = -lai, lai, 1
641 50009388 : DO mak = -lak, lak, 1
642 35749560 : isoa1 = indso_inv(lai, mai)
643 35749560 : isoa2 = indso_inv(lak, mak)
644 137109082 : DO mbj = -lbj, lbj, 1
645 312943791 : DO ilist = 1, ncg_none0(isoa1, isoa2)
646 186093617 : isoaa = cg_none0_list(isoa1, isoa2, ilist)
647 186093617 : laa = indso(1, isoaa)
648 186093617 : maa = indso(2, isoaa)
649 186093617 : nla = nsoset(laa - 1) + laa + 1
650 186093617 : labmin = MIN(laa, lbj)
651 186093617 : il = INT((lai + lak - laa)/2)
652 186093617 : stemp = 0.0_dp
653 576178399 : DO j = 0, labmin
654 390084782 : prefac = swork_cont(laa + lbj - j + 1, il, ishella, jshellb, kshella)
655 576178399 : stemp = stemp + prefac*Waux_mat(j + 1, nla + maa, nlb + mbj)
656 : END DO
657 277194231 : saba(ia + mai, jb + mbj, ka + mak) = saba(ia + mai, jb + mbj, ka + mak) + cg_coeff(isoa1, isoa2, isoaa)*stemp
658 : END DO
659 : END DO
660 : END DO
661 : END DO
662 : END DO
663 : END DO
664 : END DO
665 :
666 70378 : END SUBROUTINE construct_overlap_shg_aba
667 :
668 : ! **************************************************************************************************
669 : !> \brief calculates derivatives of [aba] SHG overlap integrals using
670 : !> precomputed angular-dependent part
671 : !> \param la set of l quantum number on a, orbital basis
672 : !> \param first_sgfa indexing
673 : !> \param nshella number of shells for a, orbital basis
674 : !> \param lb set of l quantum number on b. orbital basis
675 : !> \param first_sgfb indexing
676 : !> \param nshellb number of shells for b, orbital basis
677 : !> \param lca of l quantum number on a, aux basis
678 : !> \param first_sgfca indexing
679 : !> \param nshellca number of shells for a, aux basis
680 : !> \param cg_coeff Clebsch-Gordon coefficients
681 : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
682 : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
683 : !> \param rab distance vector Ra-Rb
684 : !> \param swork_cont contracted and normalized [s|ra^n|s] integrals
685 : !> \param Waux_mat precomputed angular-dependent part
686 : !> \param dWaux_mat derivatives of precomputed angular-dependent part
687 : !> \param dsaba derivative of contracted overlap [aba] of spherical harmonic
688 : !> Gaussians
689 : ! **************************************************************************************************
690 53612 : SUBROUTINE dev_overlap_shg_aba(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
691 107224 : lca, first_sgfca, nshellca, cg_coeff, cg_none0_list, &
692 53612 : ncg_none0, rab, swork_cont, Waux_mat, dWaux_mat, dsaba)
693 :
694 : INTEGER, DIMENSION(:), INTENT(IN) :: la, first_sgfa
695 : INTEGER, INTENT(IN) :: nshella
696 : INTEGER, DIMENSION(:), INTENT(IN) :: lb, first_sgfb
697 : INTEGER, INTENT(IN) :: nshellb
698 : INTEGER, DIMENSION(:), INTENT(IN) :: lca, first_sgfca
699 : INTEGER, INTENT(IN) :: nshellca
700 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff
701 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list
702 : INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0
703 : REAL(KIND=dp), INTENT(IN) :: rab(3)
704 : REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), &
705 : INTENT(IN) :: swork_cont
706 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Waux_mat
707 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dWaux_mat
708 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
709 : INTENT(INOUT) :: dsaba
710 :
711 : INTEGER :: i, ia, il, ilist, ishella, isoa1, isoa2, isoaa, j, jb, jshellb, ka, kshella, laa, &
712 : labmin, lai, lak, lbj, maa, mai, mak, mbj, nla, nlb, sgfa, sgfb, sgfca
713 : REAL(KIND=dp) :: dprefac, dtemp(3), prefac, rabx2(3)
714 :
715 214448 : rabx2(:) = 2.0_dp*rab
716 :
717 186925 : DO kshella = 1, nshellca
718 133313 : lak = lca(kshella)
719 133313 : sgfca = first_sgfca(kshella)
720 133313 : ka = sgfca + lak
721 843100 : DO jshellb = 1, nshellb
722 656175 : lbj = lb(jshellb)
723 656175 : nlb = nsoset(lbj - 1) + lbj + 1
724 656175 : sgfb = first_sgfb(jshellb)
725 656175 : jb = sgfb + lbj
726 4042313 : DO ishella = 1, nshella
727 3252825 : lai = la(ishella)
728 3252825 : sgfa = first_sgfa(ishella)
729 3252825 : ia = sgfa + lai
730 12346971 : DO mai = -lai, lai, 1
731 40565761 : DO mak = -lak, lak, 1
732 28874965 : isoa1 = indso_inv(lai, mai)
733 28874965 : isoa2 = indso_inv(lak, mak)
734 112211889 : DO mbj = -lbj, lbj, 1
735 255925300 : DO ilist = 1, ncg_none0(isoa1, isoa2)
736 152151382 : isoaa = cg_none0_list(isoa1, isoa2, ilist)
737 152151382 : laa = indso(1, isoaa)
738 152151382 : maa = indso(2, isoaa)
739 152151382 : nla = nsoset(laa - 1) + laa + 1
740 152151382 : labmin = MIN(laa, lbj)
741 152151382 : il = (lai + lak - laa)/2 ! lai+lak-laa always even
742 152151382 : dtemp = 0.0_dp
743 473206046 : DO j = 0, labmin
744 321054664 : prefac = swork_cont(laa + lbj - j + 1, il, ishella, jshellb, kshella)
745 321054664 : dprefac = swork_cont(laa + lbj - j + 2, il, ishella, jshellb, kshella)
746 1436370038 : DO i = 1, 3
747 : dtemp(i) = dtemp(i) + rabx2(i)*dprefac*Waux_mat(j + 1, nla + maa, nlb + mbj) &
748 1284218656 : + prefac*dWaux_mat(i, j + 1, nla + maa, nlb + mbj)
749 : END DO
750 : END DO
751 683504481 : DO i = 1, 3
752 : dsaba(ia + mai, jb + mbj, ka + mak, i) = dsaba(ia + mai, jb + mbj, ka + mak, i) &
753 608605528 : + cg_coeff(isoa1, isoa2, isoaa)*dtemp(i)
754 : END DO
755 : END DO
756 : END DO
757 : END DO
758 : END DO
759 : END DO
760 : END DO
761 : END DO
762 :
763 53612 : END SUBROUTINE dev_overlap_shg_aba
764 :
765 : ! **************************************************************************************************
766 : !> \brief calculates [abb] SHG overlap integrals using precomputed angular-
767 : !> dependent part
768 : !> \param la set of l quantum number on a, orbital basis
769 : !> \param first_sgfa indexing
770 : !> \param nshella number of shells for a, orbital basis
771 : !> \param lb set of l quantum number on b. orbital basis
772 : !> \param first_sgfb indexing
773 : !> \param nshellb number of shells for b, orbital basis
774 : !> \param lcb l quantum number on b, aux basis
775 : !> \param first_sgfcb indexing
776 : !> \param nshellcb number of shells for b, aux basis
777 : !> \param cg_coeff Clebsch-Gordon coefficients
778 : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
779 : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
780 : !> \param swork_cont contracted and normalized [s|rb^n|s] integrals
781 : !> \param Waux_mat precomputed angular-dependent part
782 : !> \param sabb contracted overlap [abb] of spherical harmonic Gaussians
783 : ! **************************************************************************************************
784 3421 : SUBROUTINE construct_overlap_shg_abb(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
785 3421 : lcb, first_sgfcb, nshellcb, cg_coeff, cg_none0_list, &
786 3421 : ncg_none0, swork_cont, Waux_mat, sabb)
787 :
788 : INTEGER, DIMENSION(:), INTENT(IN) :: la, first_sgfa
789 : INTEGER, INTENT(IN) :: nshella
790 : INTEGER, DIMENSION(:), INTENT(IN) :: lb, first_sgfb
791 : INTEGER, INTENT(IN) :: nshellb
792 : INTEGER, DIMENSION(:), INTENT(IN) :: lcb, first_sgfcb
793 : INTEGER, INTENT(IN) :: nshellcb
794 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff
795 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list
796 : INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0
797 : REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), &
798 : INTENT(IN) :: swork_cont
799 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Waux_mat
800 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: sabb
801 :
802 : INTEGER :: ia, il, ilist, ishella, isob1, isob2, isobb, j, jb, jshellb, kb, kshellb, labmin, &
803 : lai, lbb, lbj, lbk, mai, mbb, mbj, mbk, nla, nlb, sgfa, sgfb, sgfcb
804 : REAL(KIND=dp) :: prefac, stemp, tsign
805 :
806 15364 : DO kshellb = 1, nshellcb
807 11943 : lbk = lcb(kshellb)
808 11943 : sgfcb = first_sgfcb(kshellb)
809 11943 : kb = sgfcb + lbk
810 64033 : DO jshellb = 1, nshellb
811 48669 : lbj = lb(jshellb)
812 48669 : sgfb = first_sgfb(jshellb)
813 48669 : jb = sgfb + lbj
814 217476 : DO ishella = 1, nshella
815 156864 : lai = la(ishella)
816 156864 : nla = nsoset(lai - 1) + lai + 1
817 156864 : sgfa = first_sgfa(ishella)
818 156864 : ia = sgfa + lai
819 579495 : DO mbj = -lbj, lbj, 1
820 2206050 : DO mbk = -lbk, lbk, 1
821 1675224 : isob1 = indso_inv(lbj, mbj)
822 1675224 : isob2 = indso_inv(lbk, mbk)
823 5056410 : DO mai = -lai, lai, 1
824 11504847 : DO ilist = 1, ncg_none0(isob1, isob2)
825 6822399 : isobb = cg_none0_list(isob1, isob2, ilist)
826 6822399 : lbb = indso(1, isobb)
827 6822399 : mbb = indso(2, isobb)
828 6822399 : nlb = nsoset(lbb - 1) + lbb + 1
829 : ! tsgin: because we take the transpose of auxmat (calculated for (la,lb))
830 6822399 : tsign = 1.0_dp
831 6822399 : IF (MODULO(lbb - lai, 2) /= 0) tsign = -1.0_dp
832 6822399 : labmin = MIN(lai, lbb)
833 6822399 : il = INT((lbj + lbk - lbb)/2)
834 6822399 : stemp = 0.0_dp
835 18422882 : DO j = 0, labmin
836 11600483 : prefac = swork_cont(lai + lbb - j + 1, il, ishella, jshellb, kshellb)
837 18422882 : stemp = stemp + prefac*Waux_mat(j + 1, nlb + mbb, nla + mai)
838 : END DO
839 9829623 : sabb(ia + mai, jb + mbj, kb + mbk) = sabb(ia + mai, jb + mbj, kb + mbk) + tsign*cg_coeff(isob1, isob2, isobb)*stemp
840 : END DO
841 : END DO
842 : END DO
843 : END DO
844 : END DO
845 : END DO
846 : END DO
847 :
848 3421 : END SUBROUTINE construct_overlap_shg_abb
849 :
850 : ! **************************************************************************************************
851 : !> \brief calculates derivatives of [abb] SHG overlap integrals using
852 : !> precomputed angular-dependent part
853 : !> \param la set of l quantum number on a, orbital basis
854 : !> \param first_sgfa indexing
855 : !> \param nshella number of shells for a, orbital basis
856 : !> \param lb set of l quantum number on b. orbital basis
857 : !> \param first_sgfb indexing
858 : !> \param nshellb number of shells for b, orbital basis
859 : !> \param lcb l quantum number on b, aux basis
860 : !> \param first_sgfcb indexing
861 : !> \param nshellcb number of shells for b, aux basis
862 : !> \param cg_coeff Clebsch-Gordon coefficients
863 : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
864 : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
865 : !> \param rab distance vector Ra-Rb
866 : !> \param swork_cont contracted and normalized [s|rb^n|s] integrals
867 : !> \param Waux_mat precomputed angular-dependent part
868 : !> \param dWaux_mat derivatives of precomputed angular-dependent part
869 : !> \param dsabb derivative of contracted overlap [abb] of spherical harmonic
870 : !> Gaussians
871 : ! **************************************************************************************************
872 1111 : SUBROUTINE dev_overlap_shg_abb(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
873 2222 : lcb, first_sgfcb, nshellcb, cg_coeff, cg_none0_list, &
874 1111 : ncg_none0, rab, swork_cont, Waux_mat, dWaux_mat, dsabb)
875 :
876 : INTEGER, DIMENSION(:), INTENT(IN) :: la, first_sgfa
877 : INTEGER, INTENT(IN) :: nshella
878 : INTEGER, DIMENSION(:), INTENT(IN) :: lb, first_sgfb
879 : INTEGER, INTENT(IN) :: nshellb
880 : INTEGER, DIMENSION(:), INTENT(IN) :: lcb, first_sgfcb
881 : INTEGER, INTENT(IN) :: nshellcb
882 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff
883 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list
884 : INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0
885 : REAL(KIND=dp), INTENT(IN) :: rab(3)
886 : REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), &
887 : INTENT(IN) :: swork_cont
888 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Waux_mat
889 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dWaux_mat
890 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
891 : INTENT(INOUT) :: dsabb
892 :
893 : INTEGER :: i, ia, il, ilist, ishella, isob1, isob2, isobb, j, jb, jshellb, kb, kshellb, &
894 : labmin, lai, lbb, lbj, lbk, mai, mbb, mbj, mbk, nla, nlb, sgfa, sgfb, sgfcb
895 : REAL(KIND=dp) :: dprefac, dtemp(3), prefac, rabx2(3), &
896 : tsign
897 :
898 4444 : rabx2(:) = 2.0_dp*rab
899 :
900 4206 : DO kshellb = 1, nshellcb
901 3095 : lbk = lcb(kshellb)
902 3095 : sgfcb = first_sgfcb(kshellb)
903 3095 : kb = sgfcb + lbk
904 13951 : DO jshellb = 1, nshellb
905 9745 : lbj = lb(jshellb)
906 9745 : sgfb = first_sgfb(jshellb)
907 9745 : jb = sgfb + lbj
908 40263 : DO ishella = 1, nshella
909 27423 : lai = la(ishella)
910 27423 : nla = nsoset(lai - 1) + lai + 1
911 27423 : sgfa = first_sgfa(ishella)
912 27423 : ia = sgfa + lai
913 106345 : DO mbj = -lbj, lbj, 1
914 400259 : DO mbk = -lbk, lbk, 1
915 303659 : isob1 = indso_inv(lbj, mbj)
916 303659 : isob2 = indso_inv(lbk, mbk)
917 898811 : DO mai = -lai, lai, 1
918 2168821 : DO ilist = 1, ncg_none0(isob1, isob2)
919 1339187 : isobb = cg_none0_list(isob1, isob2, ilist)
920 1339187 : lbb = indso(1, isobb)
921 1339187 : mbb = indso(2, isobb)
922 1339187 : nlb = nsoset(lbb - 1) + lbb + 1
923 : ! tsgin: because we take the transpose of auxmat (calculated for (la,lb))
924 1339187 : tsign = 1.0_dp
925 1339187 : IF (MODULO(lbb - lai, 2) /= 0) tsign = -1.0_dp
926 1339187 : labmin = MIN(lai, lbb)
927 1339187 : il = (lbj + lbk - lbb)/2
928 1339187 : dtemp = 0.0_dp
929 3768153 : DO j = 0, labmin
930 2428966 : prefac = swork_cont(lai + lbb - j + 1, il, ishella, jshellb, kshellb)
931 2428966 : dprefac = swork_cont(lai + lbb - j + 2, il, ishella, jshellb, kshellb)
932 11055051 : DO i = 1, 3
933 : dtemp(i) = dtemp(i) + rabx2(i)*dprefac*Waux_mat(j + 1, nlb + mbb, nla + mai) &
934 9715864 : + prefac*dWaux_mat(i, j + 1, nlb + mbb, nla + mai)
935 : END DO
936 : END DO
937 5882723 : DO i = 1, 3
938 : dsabb(ia + mai, jb + mbj, kb + mbk, i) = dsabb(ia + mai, jb + mbj, kb + mbk, i) &
939 5356748 : + tsign*cg_coeff(isob1, isob2, isobb)*dtemp(i)
940 : END DO
941 : END DO
942 : END DO
943 : END DO
944 : END DO
945 : END DO
946 : END DO
947 : END DO
948 :
949 1111 : END SUBROUTINE dev_overlap_shg_abb
950 :
951 : END MODULE construct_shg
952 :
|
