Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
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 920177 : 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 920177 : Rc_00 = 1.0_dp
62 920177 : Rs_00 = 0.0_dp
63 :
64 920177 : Rlm_c(0, 0) = Rc_00
65 920177 : Rlm_s(0, 0) = Rs_00
66 :
67 : ! generate elements Rmm
68 : ! start
69 920177 : IF (l > 0) THEN
70 711993 : Rc = -0.5_dp*r(1)*Rc_00
71 711993 : Rs = -0.5_dp*r(2)*Rc_00
72 711993 : Rlm_c(1, 1) = Rc
73 711993 : Rlm_s(1, 1) = Rs
74 711993 : Rlm_c(1, -1) = -Rc
75 711993 : Rlm_s(1, -1) = Rs
76 : END IF
77 1349628 : DO li = 2, l
78 429451 : temp_c = (-r(1)*Rc + r(2)*Rs)/(REAL(2*(li - 1) + 2, dp))
79 429451 : Rs = (-r(2)*Rc - r(1)*Rs)/(REAL(2*(li - 1) + 2, dp))
80 429451 : Rc = temp_c
81 429451 : Rlm_c(li, li) = Rc
82 429451 : Rlm_s(li, li) = Rs
83 1349628 : IF (MODULO(li, 2) /= 0) THEN
84 122986 : Rlm_c(li, -li) = -Rc
85 122986 : Rlm_s(li, -li) = Rs
86 : ELSE
87 306465 : Rlm_c(li, -li) = Rc
88 306465 : Rlm_s(li, -li) = -Rs
89 : END IF
90 : END DO
91 :
92 2061621 : DO mi = 0, l - 1
93 1141444 : Rmlm = Rlm_c(mi, mi)
94 1141444 : Rlm = r(3)*Rlm_c(mi, mi)
95 1141444 : Rlm_c(mi + 1, mi) = Rlm
96 1141444 : IF (MODULO(mi, 2) /= 0) THEN
97 306465 : Rlm_c(mi + 1, -mi) = -Rlm
98 : ELSE
99 834979 : Rlm_c(mi + 1, -mi) = Rlm
100 : END IF
101 2687770 : DO li = mi + 2, l
102 626149 : prefac = (li + mi)*(li - mi)
103 626149 : Rplm = (REAL(2*li - 1, dp)*r(3)*Rlm - r2*Rmlm)/REAL(prefac, dp)
104 626149 : Rmlm = Rlm
105 626149 : Rlm = Rplm
106 626149 : Rlm_c(li, mi) = Rlm
107 1767593 : IF (MODULO(mi, 2) /= 0) THEN
108 159842 : Rlm_c(li, -mi) = -Rlm
109 : ELSE
110 466307 : Rlm_c(li, -mi) = Rlm
111 : END IF
112 : END DO
113 : END DO
114 1349628 : DO mi = 1, l - 1
115 429451 : Rmlm = Rlm_s(mi, mi)
116 429451 : Rlm = r(3)*Rlm_s(mi, mi)
117 429451 : Rlm_s(mi + 1, mi) = Rlm
118 429451 : IF (MODULO(mi, 2) /= 0) THEN
119 306465 : Rlm_s(mi + 1, -mi) = Rlm
120 : ELSE
121 122986 : Rlm_s(mi + 1, -mi) = -Rlm
122 : END IF
123 1546326 : DO li = mi + 2, l
124 196698 : prefac = (li + mi)*(li - mi)
125 196698 : Rplm = (REAL(2*li - 1, dp)*r(3)*Rlm - r2*Rmlm)/REAL(prefac, dp)
126 196698 : Rmlm = Rlm
127 196698 : Rlm = Rplm
128 196698 : Rlm_s(li, mi) = Rlm
129 626149 : IF (MODULO(mi, 2) /= 0) THEN
130 159842 : Rlm_s(li, -mi) = Rlm
131 : ELSE
132 36856 : Rlm_s(li, -mi) = -Rlm
133 : END IF
134 : END DO
135 : END DO
136 :
137 920177 : 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 920177 : 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 2981798 : DO l = 0, lmax
154 6811012 : DO m = 0, l
155 3829214 : temp = SQRT(fac(l + m)*fac(l - m))
156 3829214 : IF (MODULO(m, 2) /= 0) temp = -temp
157 3829214 : IF (m /= 0) temp = temp*SQRT(2.0_dp)
158 5890835 : A(l, m) = temp
159 : END DO
160 : END DO
161 :
162 920177 : 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 6396 : 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 45356 : DO l = 0, lmax
184 185688 : DO m = 0, l
185 140332 : bm = 1.0_dp
186 140332 : bm_m = 1.0_dp
187 140332 : bm_p = 1.0_dp
188 140332 : IF (m /= 0) bm = SQRT(2.0_dp)
189 140332 : IF (m - 1 /= 0) bm_m = SQRT(2.0_dp)
190 140332 : IF (m + 1 /= 0) bm_p = SQRT(2.0_dp)
191 140332 : dA_p(l, m) = -bm/bm_p*SQRT(REAL((l - m)*(l - m - 1), dp))
192 140332 : dA_m(l, m) = -bm/bm_m*SQRT(REAL((l + m)*(l + m - 1), dp))
193 140332 : dA(l, m) = 2.0_dp*SQRT(REAL((l + m)*(l - m), dp))
194 179292 : IF (m == 0) dA_p(l, m) = 2.0_dp*dA_p(l, m)
195 : END DO
196 : END DO
197 6396 : 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 920177 : 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 920177 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: A
223 :
224 920177 : Wa(:) = 0.0_dp
225 920177 : Wb(:) = 0.0_dp
226 920177 : Wmat(:) = 0.0_dp
227 :
228 3680708 : ALLOCATE (A(0:lmax, 0:lmax))
229 920177 : CALL get_Alm(lmax, A)
230 :
231 2701885 : DO lb = 0, lbmax
232 1781708 : nlb = nsoset(lb - 1)
233 6777893 : DO la = 0, lamax(lb)
234 4076008 : nla = nsoset(la - 1)
235 4076008 : labmin = MIN(la, lb)
236 13522267 : DO mb = 0, lb
237 7664551 : A_lbmb = A(lb, mb)
238 7664551 : IF (MODULO(lb, 2) /= 0) A_lbmb = -A_lbmb
239 27675772 : DO ma = 0, la
240 15935213 : A_lama = A(la, ma)
241 15935213 : Alm_fac = A_lama*A_lbmb
242 56198583 : DO j = 0, labmin
243 32598819 : laj = la - j
244 32598819 : lbj = lb - j
245 32598819 : prefac = Alm_fac*REAL(2**(la + lb - j), dp)*dfac(2*j - 1)
246 32598819 : delta_k = 0.5_dp
247 32598819 : Wmat = 0.0_dp
248 89205538 : DO k = 0, j
249 56606719 : ma_m = ma - k
250 56606719 : ma_p = ma + k
251 56606719 : IF (laj < ABS(ma_m) .AND. laj < ABS(ma_p)) CYCLE
252 41145420 : mb_m = mb - k
253 41145420 : mb_p = mb + k
254 41145420 : IF (lbj < ABS(mb_m) .AND. lbj < ABS(mb_p)) CYCLE
255 31929230 : IF (k /= 0) delta_k = 1.0_dp
256 31929230 : A_jk = fac(j + k)*fac(j - k)
257 31929230 : IF (k /= 0) A_jk = 2.0_dp*A_jk
258 31929230 : IF (MODULO(k, 2) /= 0) THEN
259 : sign_fac = -1.0_dp
260 : ELSE
261 22848101 : sign_fac = 1.0_dp
262 : END IF
263 31929230 : Rca_m = Rc(laj, ma_m)
264 31929230 : Rsa_m = Rs(laj, ma_m)
265 31929230 : Rca_p = Rc(laj, ma_p)
266 31929230 : Rsa_p = Rs(laj, ma_p)
267 31929230 : Rcb_m = Rc(lbj, mb_m)
268 31929230 : Rsb_m = Rs(lbj, mb_m)
269 31929230 : Rcb_p = Rc(lbj, mb_p)
270 31929230 : Rsb_p = Rs(lbj, mb_p)
271 31929230 : Wa(1) = delta_k*(Rca_m + sign_fac*Rca_p)
272 31929230 : Wb(1) = delta_k*(Rcb_m + sign_fac*Rcb_p)
273 31929230 : Wa(2) = -Rsa_m + sign_fac*Rsa_p
274 31929230 : Wb(2) = -Rsb_m + sign_fac*Rsb_p
275 31929230 : Wmat(1) = Wmat(1) + prefac/A_jk*(Wa(1)*Wb(1) + Wa(2)*Wb(2))
276 31929230 : IF (mb > 0) THEN
277 18546860 : Wb(3) = delta_k*(Rsb_m + sign_fac*Rsb_p)
278 18546860 : Wb(4) = Rcb_m - sign_fac*Rcb_p
279 18546860 : Wmat(2) = Wmat(2) + prefac/A_jk*(Wa(1)*Wb(3) + Wa(2)*Wb(4))
280 : END IF
281 31929230 : IF (ma > 0) THEN
282 19215685 : Wa(3) = delta_k*(Rsa_m + sign_fac*Rsa_p)
283 19215685 : Wa(4) = Rca_m - sign_fac*Rca_p
284 19215685 : Wmat(3) = Wmat(3) + prefac/A_jk*(Wa(3)*Wb(1) + Wa(4)*Wb(2))
285 : END IF
286 64528049 : IF (ma > 0 .AND. mb > 0) THEN
287 12277587 : Wmat(4) = Wmat(4) + prefac/A_jk*(Wa(3)*Wb(3) + Wa(4)*Wb(4))
288 : END IF
289 : END DO
290 32598819 : Waux_mat(j + 1, nla + la + 1 + ma, nlb + lb + 1 + mb) = Wmat(1)
291 32598819 : IF (mb > 0) Waux_mat(j + 1, nla + la + 1 + ma, nlb + lb + 1 - mb) = Wmat(2)
292 32598819 : IF (ma > 0) Waux_mat(j + 1, nla + la + 1 - ma, nlb + lb + 1 + mb) = Wmat(3)
293 48534032 : 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 920177 : DEALLOCATE (A)
301 :
302 920177 : 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 6396 : 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 6396 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dA, dA_m, dA_p, Wam, Wamm, Wamp, Wbm, &
326 6396 : Wbmm, Wbmp
327 :
328 32680 : jmax = MIN(MAXVAL(lamax), lbmax)
329 63960 : ALLOCATE (Wam(0:jmax, 4), Wamm(0:jmax, 4), Wamp(0:jmax, 4))
330 63960 : 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 32680 : lmax = MAX(MAXVAL(lamax), lbmax)
336 63960 : ALLOCATE (dA_p(0:lmax, 0:lmax), dA_m(0:lmax, 0:lmax), dA(0:lmax, 0:lmax))
337 6396 : CALL get_dA_prefactors(lmax, dA_p, dA_m, dA)
338 :
339 32680 : DO lb = 0, lbmax
340 26284 : nlb = nsoset(lb - 1)
341 26284 : nlbm = 0
342 26284 : IF (lb > 0) nlbm = nsoset(lb - 2)
343 179986 : DO la = 0, lamax(lb)
344 147306 : nla = nsoset(la - 1)
345 147306 : nlam = 0
346 147306 : IF (la > 0) nlam = nsoset(la - 2)
347 147306 : labmin = MIN(la, lb)
348 147306 : lamb = MIN(la - 1, lb)
349 147306 : labm = MIN(la, lb - 1)
350 537322 : DO mb = 0, lb
351 363732 : dAb = dA(lb, mb)
352 363732 : dAb_p = dA_p(lb, mb)
353 363732 : dAb_m = dA_m(lb, mb)
354 363732 : ipb = nlb + lb + mb + 1
355 363732 : imb = nlb + lb - mb + 1
356 363732 : ipbm = nlbm + lb + mb
357 363732 : imbm = nlbm + lb - mb
358 1730862 : DO ma = 0, la
359 1219824 : dAa = dA(la, ma)
360 1219824 : dAa_p = dA_p(la, ma)
361 1219824 : dAa_m = dA_m(la, ma)
362 1219824 : ipa = nla + la + ma + 1
363 1219824 : ima = nla + la - ma + 1
364 1219824 : ipam = nlam + la + ma
365 1219824 : imam = nlam + la - ma
366 27537348 : Wam(:, :) = 0.0_dp
367 27537348 : Wamm(:, :) = 0.0_dp
368 27537348 : Wamp(:, :) = 0.0_dp
369 : !*** Wam: la-1, ma
370 1219824 : IF (ma <= la - 1) THEN
371 2924270 : Wam(0:lamb, 1) = Waux_mat(1:lamb + 1, ipam, ipb)
372 2183157 : IF (mb > 0) Wam(0:lamb, 2) = Waux_mat(1:lamb + 1, ipam, imb)
373 2305787 : IF (ma > 0) Wam(0:lamb, 3) = Waux_mat(1:lamb + 1, imam, ipb)
374 1787750 : 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 1219824 : IF (ma - 1 >= 0) THEN
378 2924270 : Wamm(0:lamb, 1) = Waux_mat(1:lamb + 1, ipam - 1, ipb)
379 2183157 : IF (mb > 0) Wamm(0:lamb, 2) = Waux_mat(1:lamb + 1, ipam - 1, imb)
380 2305787 : 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 1787750 : 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 1219824 : IF (ma + 1 <= la - 1) THEN
385 2009941 : Wamp(0:lamb, 1) = Waux_mat(1:lamb + 1, ipam + 1, ipb)
386 1491904 : IF (mb > 0) Wamp(0:lamb, 2) = Waux_mat(1:lamb + 1, ipam + 1, imb)
387 2009941 : IF (ma + 1 > 0) Wamp(0:lamb, 3) = Waux_mat(1:lamb + 1, imam - 1, ipb)
388 1491904 : IF (ma + 1 > 0 .AND. mb > 0) Wamp(0:lamb, 4) = Waux_mat(1:lamb + 1, imam - 1, imb)
389 : END IF
390 27537348 : Wbm(:, :) = 0.0_dp
391 27537348 : Wbmm(:, :) = 0.0_dp
392 27537348 : Wbmp(:, :) = 0.0_dp
393 : !*** Wbm: lb-1, mb
394 1219824 : IF (mb <= lb - 1) THEN
395 2266824 : Wbm(0:labm, 1) = Waux_mat(1:labm + 1, ipa, ipbm)
396 1597563 : IF (mb > 0) Wbm(0:labm, 2) = Waux_mat(1:labm + 1, ipa, imbm)
397 1838473 : IF (ma > 0) Wbm(0:labm, 3) = Waux_mat(1:labm + 1, ima, ipbm)
398 1354185 : 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 1219824 : IF (mb - 1 >= 0) THEN
402 2266824 : Wbmm(0:labm, 1) = Waux_mat(1:labm + 1, ipa, ipbm - 1)
403 1597563 : IF (mb - 1 > 0) Wbmm(0:labm, 2) = Waux_mat(1:labm + 1, ipa, imbm + 1)
404 1838473 : IF (ma > 0) Wbmm(0:labm, 3) = Waux_mat(1:labm + 1, ima, ipbm - 1)
405 1354185 : 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 1219824 : IF (mb + 1 <= lb - 1) THEN
409 1229384 : Wbmp(0:labm, 1) = Waux_mat(1:labm + 1, ipa, ipbm + 1)
410 1229384 : IF (mb + 1 > 0) Wbmp(0:labm, 2) = Waux_mat(1:labm + 1, ipa, imbm - 1)
411 986006 : IF (ma > 0) Wbmp(0:labm, 3) = Waux_mat(1:labm + 1, ima, ipbm + 1)
412 986006 : IF (ma > 0 .AND. mb + 1 > 0) Wbmp(0:labm, 4) = Waux_mat(1:labm + 1, ima, imbm - 1)
413 : END IF
414 4751406 : 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 3167850 : - dAb_p*Wbmp(j, 1) + dAb_m*Wbmm(j, 1)
418 3167850 : IF (mb > 0) THEN
419 : dWaux_mat(1, j + 1, ipa, imb) = dAa_p*Wamp(j, 2) - dAa_m*Wamm(j, 2) &
420 2064253 : - dAb_p*Wbmp(j, 2) + dAb_m*Wbmm(j, 2)
421 : END IF
422 3167850 : IF (ma > 0) THEN
423 : dWaux_mat(1, j + 1, ima, ipb) = dAa_p*Wamp(j, 3) - dAa_m*Wamm(j, 3) &
424 2336904 : - dAb_p*Wbmp(j, 3) + dAb_m*Wbmm(j, 3)
425 : END IF
426 3167850 : 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 1523984 : - 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 3167850 : - dAb_p*Wbmp(j, 2) - dAb_m*Wbmm(j, 2)
434 3167850 : IF (mb > 0) THEN
435 : dWaux_mat(2, j + 1, ipa, imb) = dAa_p*Wamp(j, 4) + dAa_m*Wamm(j, 4) &
436 2064253 : + dAb_p*Wbmp(j, 1) + dAb_m*Wbmm(j, 1)
437 : END IF
438 3167850 : IF (ma > 0) THEN
439 : dWaux_mat(2, j + 1, ima, ipb) = -dAa_p*Wamp(j, 1) - dAa_m*Wamm(j, 1) &
440 2336904 : - dAb_p*Wbmp(j, 4) - dAb_m*Wbmm(j, 4)
441 : END IF
442 3167850 : 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 1523984 : + dAb_p*Wbmp(j, 3) + dAb_m*Wbmm(j, 3)
445 : END IF
446 : !**** z compnent
447 3167850 : dWaux_mat(3, j + 1, ipa, ipb) = dAa*Wam(j, 1) - dAb*Wbm(j, 1)
448 3167850 : IF (mb > 0) THEN
449 2064253 : dWaux_mat(3, j + 1, ipa, imb) = dAa*Wam(j, 2) - dAb*Wbm(j, 2)
450 : END IF
451 3167850 : IF (ma > 0) THEN
452 2336904 : dWaux_mat(3, j + 1, ima, ipb) = dAa*Wam(j, 3) - dAb*Wbm(j, 3)
453 : END IF
454 4387674 : IF (ma > 0 .AND. mb > 0) THEN
455 1523984 : 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 6396 : DEALLOCATE (Wam, Wamm, Wamp)
465 6396 : DEALLOCATE (Wbm, Wbmm, Wbmp)
466 6396 : DEALLOCATE (dA, dA_p, dA_m)
467 :
468 6396 : 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 17837956 : SUBROUTINE construct_int_shg_ab(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
484 17837956 : 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 35900335 : DO jshellb = 1, nshellb
499 18062379 : lbj = lb(jshellb)
500 18062379 : fnlb = nsoset(lbj - 1) + 1
501 18062379 : lnlb = nsoset(lbj)
502 18062379 : fsgfb = first_sgfb(jshellb)
503 18062379 : lsgfb = fsgfb + 2*lbj
504 54682454 : DO ishella = 1, nshella
505 18782119 : lai = la(ishella)
506 18782119 : fnla = nsoset(lai - 1) + 1
507 18782119 : lnla = nsoset(lai)
508 18782119 : fsgfa = first_sgfa(ishella)
509 18782119 : lsgfa = fsgfa + 2*lai
510 18782119 : labmin = MIN(lai, lbj)
511 89464391 : DO mbj = 0, 2*lbj
512 220786247 : DO mai = 0, 2*lai
513 509253584 : DO j = 0, labmin
514 307249456 : prefac = swork_cont(lai + lbj - j + 1, ishella, jshellb)
515 : sab(fsgfa + mai, fsgfb + mbj) = sab(fsgfa + mai, fsgfb + mbj) &
516 456633691 : + 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 17837956 : 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 77242 : SUBROUTINE construct_dev_shg_ab(la, first_sgfa, nshella, lb, first_sgfb, nshellb, rab, &
541 77242 : 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 308968 : rabx2(:) = 2.0_dp*rab
558 324154 : DO jshellb = 1, nshellb
559 246912 : lbj = lb(jshellb)
560 246912 : fnlb = nsoset(lbj - 1) + 1
561 246912 : lnlb = nsoset(lbj)
562 246912 : fsgfb = first_sgfb(jshellb)
563 246912 : lsgfb = fsgfb + 2*lbj
564 1108796 : DO ishella = 1, nshella
565 784642 : lai = la(ishella)
566 784642 : fnla = nsoset(lai - 1) + 1
567 784642 : lnla = nsoset(lai)
568 784642 : fsgfa = first_sgfa(ishella)
569 784642 : lsgfa = fsgfa + 2*lai
570 784642 : labmin = MIN(lai, lbj)
571 2345125 : DO j = 0, labmin
572 1313571 : prefac = swork_cont(lai + lbj - j + 1, ishella, jshellb)
573 1313571 : dprefac = swork_cont(lai + lbj - j + 2, ishella, jshellb) !j+1
574 6038926 : 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 97498992 : + prefac*dWaux_mat(i, j + 1, fnla:lnla, fnlb:lnlb)
578 : END DO
579 : END DO
580 : END DO
581 : END DO
582 :
583 77242 : 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 40040 : SUBROUTINE construct_overlap_shg_aba(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
605 40040 : lca, first_sgfca, nshellca, cg_coeff, cg_none0_list, &
606 40040 : 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 143017 : DO kshella = 1, nshellca
627 102977 : lak = lca(kshella)
628 102977 : sgfca = first_sgfca(kshella)
629 102977 : ka = sgfca + lak
630 599463 : DO jshellb = 1, nshellb
631 456446 : lbj = lb(jshellb)
632 456446 : nlb = nsoset(lbj - 1) + lbj + 1
633 456446 : sgfb = first_sgfb(jshellb)
634 456446 : jb = sgfb + lbj
635 2711645 : DO ishella = 1, nshella
636 2152222 : lai = la(ishella)
637 2152222 : sgfa = first_sgfa(ishella)
638 2152222 : ia = sgfa + lai
639 8069462 : DO mai = -lai, lai, 1
640 26951258 : DO mak = -lak, lak, 1
641 19338242 : isoa1 = indso_inv(lai, mai)
642 19338242 : isoa2 = indso_inv(lak, mak)
643 73280668 : DO mbj = -lbj, lbj, 1
644 167368307 : DO ilist = 1, ncg_none0(isoa1, isoa2)
645 99548433 : isoaa = cg_none0_list(isoa1, isoa2, ilist)
646 99548433 : laa = indso(1, isoaa)
647 99548433 : maa = indso(2, isoaa)
648 99548433 : nla = nsoset(laa - 1) + laa + 1
649 99548433 : labmin = MIN(laa, lbj)
650 99548433 : il = INT((lai + lak - laa)/2)
651 99548433 : stemp = 0.0_dp
652 306983753 : DO j = 0, labmin
653 207435320 : prefac = swork_cont(laa + lbj - j + 1, il, ishella, jshellb, kshella)
654 306983753 : stemp = stemp + prefac*Waux_mat(j + 1, nla + maa, nlb + mbj)
655 : END DO
656 148030065 : 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 40040 : 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 27388 : SUBROUTINE dev_overlap_shg_aba(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
690 54776 : lca, first_sgfca, nshellca, cg_coeff, cg_none0_list, &
691 27388 : 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 109552 : rabx2(:) = 2.0_dp*rab
715 :
716 95801 : DO kshella = 1, nshellca
717 68413 : lak = lca(kshella)
718 68413 : sgfca = first_sgfca(kshella)
719 68413 : ka = sgfca + lak
720 428276 : DO jshellb = 1, nshellb
721 332475 : lbj = lb(jshellb)
722 332475 : nlb = nsoset(lbj - 1) + lbj + 1
723 332475 : sgfb = first_sgfb(jshellb)
724 332475 : jb = sgfb + lbj
725 2038117 : DO ishella = 1, nshella
726 1637229 : lai = la(ishella)
727 1637229 : sgfa = first_sgfa(ishella)
728 1637229 : ia = sgfa + lai
729 6211191 : DO mai = -lai, lai, 1
730 20427013 : DO mak = -lak, lak, 1
731 14548297 : isoa1 = indso_inv(lai, mai)
732 14548297 : isoa2 = indso_inv(lak, mak)
733 56467669 : DO mbj = -lbj, lbj, 1
734 128860480 : DO ilist = 1, ncg_none0(isoa1, isoa2)
735 76634298 : isoaa = cg_none0_list(isoa1, isoa2, ilist)
736 76634298 : laa = indso(1, isoaa)
737 76634298 : maa = indso(2, isoaa)
738 76634298 : nla = nsoset(laa - 1) + laa + 1
739 76634298 : labmin = MIN(laa, lbj)
740 76634298 : il = (lai + lak - laa)/2 ! lai+lak-laa always even
741 76634298 : dtemp = 0.0_dp
742 238324798 : DO j = 0, labmin
743 161690500 : prefac = swork_cont(laa + lbj - j + 1, il, ishella, jshellb, kshella)
744 161690500 : dprefac = swork_cont(laa + lbj - j + 2, il, ishella, jshellb, kshella)
745 723396298 : DO i = 1, 3
746 : dtemp(i) = dtemp(i) + rabx2(i)*dprefac*Waux_mat(j + 1, nla + maa, nlb + mbj) &
747 646762000 : + prefac*dWaux_mat(i, j + 1, nla + maa, nlb + mbj)
748 : END DO
749 : END DO
750 344215077 : DO i = 1, 3
751 : dsaba(ia + mai, jb + mbj, ka + mak, i) = dsaba(ia + mai, jb + mbj, ka + mak, i) &
752 306537192 : + 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 27388 : 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 3277 : SUBROUTINE construct_overlap_shg_abb(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
784 3277 : lcb, first_sgfcb, nshellcb, cg_coeff, cg_none0_list, &
785 3277 : 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 14952 : DO kshellb = 1, nshellcb
806 11675 : lbk = lcb(kshellb)
807 11675 : sgfcb = first_sgfcb(kshellb)
808 11675 : kb = sgfcb + lbk
809 62545 : DO jshellb = 1, nshellb
810 47593 : lbj = lb(jshellb)
811 47593 : sgfb = first_sgfb(jshellb)
812 47593 : jb = sgfb + lbj
813 212112 : DO ishella = 1, nshella
814 152844 : lai = la(ishella)
815 152844 : nla = nsoset(lai - 1) + lai + 1
816 152844 : sgfa = first_sgfa(ishella)
817 152844 : ia = sgfa + lai
818 565795 : DO mbj = -lbj, lbj, 1
819 2165118 : DO mbk = -lbk, lbk, 1
820 1646916 : isob1 = indso_inv(lbj, mbj)
821 1646916 : isob2 = indso_inv(lbk, mbk)
822 4963078 : DO mai = -lai, lai, 1
823 11318995 : DO ilist = 1, ncg_none0(isob1, isob2)
824 6721275 : isobb = cg_none0_list(isob1, isob2, ilist)
825 6721275 : lbb = indso(1, isobb)
826 6721275 : mbb = indso(2, isobb)
827 6721275 : nlb = nsoset(lbb - 1) + lbb + 1
828 : ! tsgin: because we take the transpose of auxmat (calculated for (la,lb))
829 6721275 : tsign = 1.0_dp
830 6721275 : IF (MODULO(lbb - lai, 2) /= 0) tsign = -1.0_dp
831 6721275 : labmin = MIN(lai, lbb)
832 6721275 : il = INT((lbj + lbk - lbb)/2)
833 6721275 : stemp = 0.0_dp
834 18145266 : DO j = 0, labmin
835 11423991 : prefac = swork_cont(lai + lbb - j + 1, il, ishella, jshellb, kshellb)
836 18145266 : stemp = stemp + prefac*Waux_mat(j + 1, nlb + mbb, nla + mai)
837 : END DO
838 9672079 : 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 3277 : 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 967 : SUBROUTINE dev_overlap_shg_abb(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
872 1934 : lcb, first_sgfcb, nshellcb, cg_coeff, cg_none0_list, &
873 967 : 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 3868 : rabx2(:) = 2.0_dp*rab
898 :
899 3794 : DO kshellb = 1, nshellcb
900 2827 : lbk = lcb(kshellb)
901 2827 : sgfcb = first_sgfcb(kshellb)
902 2827 : kb = sgfcb + lbk
903 12463 : DO jshellb = 1, nshellb
904 8669 : lbj = lb(jshellb)
905 8669 : sgfb = first_sgfb(jshellb)
906 8669 : jb = sgfb + lbj
907 34899 : DO ishella = 1, nshella
908 23403 : lai = la(ishella)
909 23403 : nla = nsoset(lai - 1) + lai + 1
910 23403 : sgfa = first_sgfa(ishella)
911 23403 : ia = sgfa + lai
912 92645 : DO mbj = -lbj, lbj, 1
913 359327 : DO mbk = -lbk, lbk, 1
914 275351 : isob1 = indso_inv(lbj, mbj)
915 275351 : isob2 = indso_inv(lbk, mbk)
916 805479 : DO mai = -lai, lai, 1
917 1982969 : DO ilist = 1, ncg_none0(isob1, isob2)
918 1238063 : isobb = cg_none0_list(isob1, isob2, ilist)
919 1238063 : lbb = indso(1, isobb)
920 1238063 : mbb = indso(2, isobb)
921 1238063 : nlb = nsoset(lbb - 1) + lbb + 1
922 : ! tsgin: because we take the transpose of auxmat (calculated for (la,lb))
923 1238063 : tsign = 1.0_dp
924 1238063 : IF (MODULO(lbb - lai, 2) /= 0) tsign = -1.0_dp
925 1238063 : labmin = MIN(lai, lbb)
926 1238063 : il = (lbj + lbk - lbb)/2
927 1238063 : dtemp = 0.0_dp
928 3490537 : DO j = 0, labmin
929 2252474 : prefac = swork_cont(lai + lbb - j + 1, il, ishella, jshellb, kshellb)
930 2252474 : dprefac = swork_cont(lai + lbb - j + 2, il, ishella, jshellb, kshellb)
931 10247959 : DO i = 1, 3
932 : dtemp(i) = dtemp(i) + rabx2(i)*dprefac*Waux_mat(j + 1, nlb + mbb, nla + mai) &
933 9009896 : + prefac*dWaux_mat(i, j + 1, nlb + mbb, nla + mai)
934 : END DO
935 : END DO
936 5421807 : DO i = 1, 3
937 : dsabb(ia + mai, jb + mbj, kb + mbk, i) = dsabb(ia + mai, jb + mbj, kb + mbk, i) &
938 4952252 : + 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 967 : END SUBROUTINE dev_overlap_shg_abb
949 :
950 : END MODULE construct_shg
951 :
|
