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 7221832 : 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 7221832 : Rc_00 = 1.0_dp
62 7221832 : Rs_00 = 0.0_dp
63 :
64 7221832 : Rlm_c(0, 0) = Rc_00
65 7221832 : Rlm_s(0, 0) = Rs_00
66 :
67 : ! generate elements Rmm
68 : ! start
69 7221832 : IF (l > 0) THEN
70 4140372 : Rc = -0.5_dp*r(1)*Rc_00
71 4140372 : Rs = -0.5_dp*r(2)*Rc_00
72 4140372 : Rlm_c(1, 1) = Rc
73 4140372 : Rlm_s(1, 1) = Rs
74 4140372 : Rlm_c(1, -1) = -Rc
75 4140372 : Rlm_s(1, -1) = Rs
76 : END IF
77 7784524 : DO li = 2, l
78 562692 : temp_c = (-r(1)*Rc + r(2)*Rs)/(REAL(2*(li - 1) + 2, dp))
79 562692 : Rs = (-r(2)*Rc - r(1)*Rs)/(REAL(2*(li - 1) + 2, dp))
80 562692 : Rc = temp_c
81 562692 : Rlm_c(li, li) = Rc
82 562692 : Rlm_s(li, li) = Rs
83 7784524 : IF (MODULO(li, 2) /= 0) THEN
84 148626 : Rlm_c(li, -li) = -Rc
85 148626 : Rlm_s(li, -li) = Rs
86 : ELSE
87 414066 : Rlm_c(li, -li) = Rc
88 414066 : Rlm_s(li, -li) = -Rs
89 : END IF
90 : END DO
91 :
92 11924896 : DO mi = 0, l - 1
93 4703064 : Rmlm = Rlm_c(mi, mi)
94 4703064 : Rlm = r(3)*Rlm_c(mi, mi)
95 4703064 : Rlm_c(mi + 1, mi) = Rlm
96 4703064 : IF (MODULO(mi, 2) /= 0) THEN
97 414066 : Rlm_c(mi + 1, -mi) = -Rlm
98 : ELSE
99 4288998 : Rlm_c(mi + 1, -mi) = Rlm
100 : END IF
101 12760814 : DO li = mi + 2, l
102 835918 : prefac = (li + mi)*(li - mi)
103 835918 : Rplm = (REAL(2*li - 1, dp)*r(3)*Rlm - r2*Rmlm)/REAL(prefac, dp)
104 835918 : Rmlm = Rlm
105 835918 : Rlm = Rplm
106 835918 : Rlm_c(li, mi) = Rlm
107 5538982 : IF (MODULO(mi, 2) /= 0) THEN
108 210926 : Rlm_c(li, -mi) = -Rlm
109 : ELSE
110 624992 : Rlm_c(li, -mi) = Rlm
111 : END IF
112 : END DO
113 : END DO
114 7784524 : DO mi = 1, l - 1
115 562692 : Rmlm = Rlm_s(mi, mi)
116 562692 : Rlm = r(3)*Rlm_s(mi, mi)
117 562692 : Rlm_s(mi + 1, mi) = Rlm
118 562692 : IF (MODULO(mi, 2) /= 0) THEN
119 414066 : Rlm_s(mi + 1, -mi) = Rlm
120 : ELSE
121 148626 : Rlm_s(mi + 1, -mi) = -Rlm
122 : END IF
123 8057750 : DO li = mi + 2, l
124 273226 : prefac = (li + mi)*(li - mi)
125 273226 : Rplm = (REAL(2*li - 1, dp)*r(3)*Rlm - r2*Rmlm)/REAL(prefac, dp)
126 273226 : Rmlm = Rlm
127 273226 : Rlm = Rplm
128 273226 : Rlm_s(li, mi) = Rlm
129 835918 : IF (MODULO(mi, 2) /= 0) THEN
130 210926 : Rlm_s(li, -mi) = Rlm
131 : ELSE
132 62300 : Rlm_s(li, -mi) = -Rlm
133 : END IF
134 : END DO
135 : END DO
136 :
137 7221832 : 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 7221832 : 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 19146728 : DO l = 0, lmax
154 36610606 : DO m = 0, l
155 17463878 : temp = SQRT(fac(l + m)*fac(l - m))
156 17463878 : IF (MODULO(m, 2) /= 0) temp = -temp
157 17463878 : IF (m /= 0) temp = temp*SQRT(2.0_dp)
158 29388774 : A(l, m) = temp
159 : END DO
160 : END DO
161 :
162 7221832 : 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 12196 : 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 85948 : DO l = 0, lmax
184 348036 : DO m = 0, l
185 262088 : bm = 1.0_dp
186 262088 : bm_m = 1.0_dp
187 262088 : bm_p = 1.0_dp
188 262088 : IF (m /= 0) bm = SQRT(2.0_dp)
189 188336 : IF (m - 1 /= 0) bm_m = SQRT(2.0_dp)
190 200532 : IF (m + 1 /= 0) bm_p = SQRT(2.0_dp)
191 262088 : dA_p(l, m) = -bm/bm_p*SQRT(REAL((l - m)*(l - m - 1), dp))
192 262088 : dA_m(l, m) = -bm/bm_m*SQRT(REAL((l + m)*(l + m - 1), dp))
193 262088 : dA(l, m) = 2.0_dp*SQRT(REAL((l + m)*(l - m), dp))
194 335840 : IF (m == 0) dA_p(l, m) = 2.0_dp*dA_p(l, m)
195 : END DO
196 : END DO
197 12196 : 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 7221832 : 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 7221832 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: A
223 :
224 7221832 : Wa(:) = 0.0_dp
225 7221832 : Wb(:) = 0.0_dp
226 7221832 : Wmat(:) = 0.0_dp
227 :
228 28887328 : ALLOCATE (A(0:lmax, 0:lmax))
229 7221832 : CALL get_Alm(lmax, A)
230 :
231 17419487 : DO lb = 0, lbmax
232 10197655 : nlb = nsoset(lb - 1)
233 32898176 : DO la = 0, lamax(lb)
234 15478689 : nla = nsoset(la - 1)
235 15478689 : labmin = MIN(la, lb)
236 47958283 : DO mb = 0, lb
237 22281939 : A_lbmb = A(lb, mb)
238 22281939 : IF (MODULO(lb, 2) /= 0) A_lbmb = -A_lbmb
239 73346937 : DO ma = 0, la
240 35586309 : A_lama = A(la, ma)
241 35586309 : Alm_fac = A_lama*A_lbmb
242 116068385 : DO j = 0, labmin
243 58200137 : laj = la - j
244 58200137 : lbj = lb - j
245 58200137 : prefac = Alm_fac*REAL(2**(la + lb - j), dp)*dfac(2*j - 1)
246 58200137 : delta_k = 0.5_dp
247 58200137 : Wmat = 0.0_dp
248 147932553 : DO k = 0, j
249 89732416 : ma_m = ma - k
250 89732416 : ma_p = ma + k
251 89732416 : IF (laj < ABS(ma_m) .AND. laj < ABS(ma_p)) CYCLE
252 68996612 : mb_m = mb - k
253 68996612 : mb_p = mb + k
254 68996612 : IF (lbj < ABS(mb_m) .AND. lbj < ABS(mb_p)) CYCLE
255 56309780 : IF (k /= 0) delta_k = 1.0_dp
256 56309780 : A_jk = fac(j + k)*fac(j - k)
257 56309780 : IF (k /= 0) A_jk = 2.0_dp*A_jk
258 13629711 : IF (MODULO(k, 2) /= 0) THEN
259 : sign_fac = -1.0_dp
260 : ELSE
261 : sign_fac = 1.0_dp
262 : END IF
263 56309780 : Rca_m = Rc(laj, ma_m)
264 56309780 : Rsa_m = Rs(laj, ma_m)
265 56309780 : Rca_p = Rc(laj, ma_p)
266 56309780 : Rsa_p = Rs(laj, ma_p)
267 56309780 : Rcb_m = Rc(lbj, mb_m)
268 56309780 : Rsb_m = Rs(lbj, mb_m)
269 56309780 : Rcb_p = Rc(lbj, mb_p)
270 56309780 : Rsb_p = Rs(lbj, mb_p)
271 56309780 : Wa(1) = delta_k*(Rca_m + sign_fac*Rca_p)
272 56309780 : Wb(1) = delta_k*(Rcb_m + sign_fac*Rcb_p)
273 56309780 : Wa(2) = -Rsa_m + sign_fac*Rsa_p
274 56309780 : Wb(2) = -Rsb_m + sign_fac*Rsb_p
275 56309780 : Wmat(1) = Wmat(1) + prefac/A_jk*(Wa(1)*Wb(1) + Wa(2)*Wb(2))
276 56309780 : IF (mb > 0) THEN
277 26234851 : Wb(3) = delta_k*(Rsb_m + sign_fac*Rsb_p)
278 26234851 : Wb(4) = Rcb_m - sign_fac*Rcb_p
279 26234851 : Wmat(2) = Wmat(2) + prefac/A_jk*(Wa(1)*Wb(3) + Wa(2)*Wb(4))
280 : END IF
281 56309780 : IF (ma > 0) THEN
282 27374508 : Wa(3) = delta_k*(Rsa_m + sign_fac*Rsa_p)
283 27374508 : Wa(4) = Rca_m - sign_fac*Rca_p
284 27374508 : Wmat(3) = Wmat(3) + prefac/A_jk*(Wa(3)*Wb(1) + Wa(4)*Wb(2))
285 : END IF
286 114509917 : IF (ma > 0 .AND. mb > 0) THEN
287 16198894 : Wmat(4) = Wmat(4) + prefac/A_jk*(Wa(3)*Wb(3) + Wa(4)*Wb(4))
288 : END IF
289 : END DO
290 58200137 : Waux_mat(j + 1, nla + la + 1 + ma, nlb + lb + 1 + mb) = Wmat(1)
291 58200137 : IF (mb > 0) Waux_mat(j + 1, nla + la + 1 + ma, nlb + lb + 1 - mb) = Wmat(2)
292 58200137 : IF (ma > 0) Waux_mat(j + 1, nla + la + 1 - ma, nlb + lb + 1 + mb) = Wmat(3)
293 93786446 : 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 7221832 : DEALLOCATE (A)
301 :
302 7221832 : 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 12196 : 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 12196 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dA, dA_m, dA_p, Wam, Wamm, Wamp, Wbm, &
326 12196 : Wbmm, Wbmp
327 :
328 61674 : jmax = MIN(MAXVAL(lamax), lbmax)
329 73176 : ALLOCATE (Wam(0:jmax, 4), Wamm(0:jmax, 4), Wamp(0:jmax, 4))
330 48784 : 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 61674 : lmax = MAX(MAXVAL(lamax), lbmax)
336 97568 : ALLOCATE (dA_p(0:lmax, 0:lmax), dA_m(0:lmax, 0:lmax), dA(0:lmax, 0:lmax))
337 12196 : CALL get_dA_prefactors(lmax, dA_p, dA_m, dA)
338 :
339 61674 : DO lb = 0, lbmax
340 49478 : nlb = nsoset(lb - 1)
341 49478 : nlbm = 0
342 49478 : IF (lb > 0) nlbm = nsoset(lb - 2)
343 336518 : DO la = 0, lamax(lb)
344 274844 : nla = nsoset(la - 1)
345 274844 : nlam = 0
346 274844 : IF (la > 0) nlam = nsoset(la - 2)
347 274844 : labmin = MIN(la, lb)
348 274844 : lamb = MIN(la - 1, lb)
349 274844 : labm = MIN(la, lb - 1)
350 989466 : DO mb = 0, lb
351 665144 : dAb = dA(lb, mb)
352 665144 : dAb_p = dA_p(lb, mb)
353 665144 : dAb_m = dA_m(lb, mb)
354 665144 : ipb = nlb + lb + mb + 1
355 665144 : imb = nlb + lb - mb + 1
356 665144 : ipbm = nlbm + lb + mb
357 665144 : imbm = nlbm + lb - mb
358 3121896 : DO ma = 0, la
359 2181908 : dAa = dA(la, ma)
360 2181908 : dAa_p = dA_p(la, ma)
361 2181908 : dAa_m = dA_m(la, ma)
362 2181908 : ipa = nla + la + ma + 1
363 2181908 : ima = nla + la - ma + 1
364 2181908 : ipam = nlam + la + ma
365 2181908 : imam = nlam + la - ma
366 2181908 : Wam(:, :) = 0.0_dp
367 2181908 : Wamm(:, :) = 0.0_dp
368 2181908 : Wamp(:, :) = 0.0_dp
369 : !*** Wam: la-1, ma
370 2181908 : IF (ma <= la - 1) THEN
371 5045198 : Wam(0:lamb, 1) = Waux_mat(1:lamb + 1, ipam, ipb)
372 3730347 : IF (mb > 0) Wam(0:lamb, 2) = Waux_mat(1:lamb + 1, ipam, imb)
373 3945703 : IF (ma > 0) Wam(0:lamb, 3) = Waux_mat(1:lamb + 1, imam, ipb)
374 3039402 : 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 2181908 : IF (ma - 1 >= 0) THEN
378 5045198 : Wamm(0:lamb, 1) = Waux_mat(1:lamb + 1, ipam - 1, ipb)
379 3730347 : 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 3945703 : IF (ma - 1 > 0) Wamm(0:lamb, 3) = Waux_mat(1:lamb + 1, imam + 1, ipb)
382 3039402 : 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 2181908 : IF (ma + 1 <= la - 1) THEN
386 3406421 : Wamp(0:lamb, 1) = Waux_mat(1:lamb + 1, ipam + 1, ipb)
387 2500120 : IF (mb > 0) Wamp(0:lamb, 2) = Waux_mat(1:lamb + 1, ipam + 1, imb)
388 3406421 : IF (ma + 1 > 0) Wamp(0:lamb, 3) = Waux_mat(1:lamb + 1, imam - 1, ipb)
389 2500120 : IF (ma + 1 > 0 .AND. mb > 0) Wamp(0:lamb, 4) = Waux_mat(1:lamb + 1, imam - 1, imb)
390 : END IF
391 2181908 : Wbm(:, :) = 0.0_dp
392 2181908 : Wbmm(:, :) = 0.0_dp
393 2181908 : Wbmp(:, :) = 0.0_dp
394 : !*** Wbm: lb-1, mb
395 2181908 : IF (mb <= lb - 1) THEN
396 3854568 : Wbm(0:labm, 1) = Waux_mat(1:labm + 1, ipa, ipbm)
397 2675339 : IF (mb > 0) Wbm(0:labm, 2) = Waux_mat(1:labm + 1, ipa, imbm)
398 3107485 : IF (ma > 0) Wbm(0:labm, 3) = Waux_mat(1:labm + 1, ima, ipbm)
399 2263921 : 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 2181908 : IF (mb - 1 >= 0) THEN
403 3854568 : Wbmm(0:labm, 1) = Waux_mat(1:labm + 1, ipa, ipbm - 1)
404 2675339 : IF (mb - 1 > 0) Wbmm(0:labm, 2) = Waux_mat(1:labm + 1, ipa, imbm + 1)
405 3107485 : IF (ma > 0) Wbmm(0:labm, 3) = Waux_mat(1:labm + 1, ima, ipbm - 1)
406 2263921 : 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 2181908 : IF (mb + 1 <= lb - 1) THEN
410 2005776 : Wbmp(0:labm, 1) = Waux_mat(1:labm + 1, ipa, ipbm + 1)
411 2005776 : IF (mb + 1 > 0) Wbmp(0:labm, 2) = Waux_mat(1:labm + 1, ipa, imbm - 1)
412 1594358 : IF (ma > 0) Wbmp(0:labm, 3) = Waux_mat(1:labm + 1, ima, ipbm + 1)
413 1594358 : IF (ma > 0 .AND. mb + 1 > 0) Wbmp(0:labm, 4) = Waux_mat(1:labm + 1, ima, imbm - 1)
414 : END IF
415 8332898 : 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 5485846 : - dAb_p*Wbmp(j, 1) + dAb_m*Wbmm(j, 1)
419 5485846 : IF (mb > 0) THEN
420 : dWaux_mat(1, j + 1, ipa, imb) = dAa_p*Wamp(j, 2) - dAa_m*Wamm(j, 2) &
421 3507105 : - dAb_p*Wbmp(j, 2) + dAb_m*Wbmm(j, 2)
422 : END IF
423 5485846 : IF (ma > 0) THEN
424 : dWaux_mat(1, j + 1, ima, ipb) = dAa_p*Wamp(j, 3) - dAa_m*Wamm(j, 3) &
425 3999992 : - dAb_p*Wbmp(j, 3) + dAb_m*Wbmm(j, 3)
426 : END IF
427 5485846 : 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 2555376 : - 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 5485846 : - dAb_p*Wbmp(j, 2) - dAb_m*Wbmm(j, 2)
435 5485846 : IF (mb > 0) THEN
436 : dWaux_mat(2, j + 1, ipa, imb) = dAa_p*Wamp(j, 4) + dAa_m*Wamm(j, 4) &
437 3507105 : + dAb_p*Wbmp(j, 1) + dAb_m*Wbmm(j, 1)
438 : END IF
439 5485846 : IF (ma > 0) THEN
440 : dWaux_mat(2, j + 1, ima, ipb) = -dAa_p*Wamp(j, 1) - dAa_m*Wamm(j, 1) &
441 3999992 : - dAb_p*Wbmp(j, 4) - dAb_m*Wbmm(j, 4)
442 : END IF
443 5485846 : 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 2555376 : + dAb_p*Wbmp(j, 3) + dAb_m*Wbmm(j, 3)
446 : END IF
447 : !**** z compnent
448 5485846 : dWaux_mat(3, j + 1, ipa, ipb) = dAa*Wam(j, 1) - dAb*Wbm(j, 1)
449 5485846 : IF (mb > 0) THEN
450 3507105 : dWaux_mat(3, j + 1, ipa, imb) = dAa*Wam(j, 2) - dAb*Wbm(j, 2)
451 : END IF
452 5485846 : IF (ma > 0) THEN
453 3999992 : dWaux_mat(3, j + 1, ima, ipb) = dAa*Wam(j, 3) - dAb*Wbm(j, 3)
454 : END IF
455 7667754 : IF (ma > 0 .AND. mb > 0) THEN
456 2555376 : 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 12196 : DEALLOCATE (Wam, Wamm, Wamp)
466 12196 : DEALLOCATE (Wbm, Wbmm, Wbmp)
467 12196 : DEALLOCATE (dA, dA_p, dA_m)
468 :
469 12196 : 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 28775265 : SUBROUTINE construct_int_shg_ab(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
485 28775265 : 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 57859013 : DO jshellb = 1, nshellb
500 29083748 : lbj = lb(jshellb)
501 29083748 : fnlb = nsoset(lbj - 1) + 1
502 29083748 : lnlb = nsoset(lbj)
503 29083748 : fsgfb = first_sgfb(jshellb)
504 29083748 : lsgfb = fsgfb + 2*lbj
505 87938243 : DO ishella = 1, nshella
506 30079230 : lai = la(ishella)
507 30079230 : fnla = nsoset(lai - 1) + 1
508 30079230 : lnla = nsoset(lai)
509 30079230 : fsgfa = first_sgfa(ishella)
510 30079230 : lsgfa = fsgfa + 2*lai
511 30079230 : labmin = MIN(lai, lbj)
512 129280828 : DO mbj = 0, 2*lbj
513 277722290 : DO mai = 0, 2*lai
514 593758890 : DO j = 0, labmin
515 346115830 : prefac = swork_cont(lai + lbj - j + 1, ishella, jshellb)
516 : sab(fsgfa + mai, fsgfb + mbj) = sab(fsgfa + mai, fsgfb + mbj) &
517 523641040 : + 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 28775265 : 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 113782 : SUBROUTINE construct_dev_shg_ab(la, first_sgfa, nshella, lb, first_sgfb, nshellb, rab, &
542 113782 : 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 455128 : rabx2(:) = 2.0_dp*rab
559 463662 : DO jshellb = 1, nshellb
560 349880 : lbj = lb(jshellb)
561 349880 : fnlb = nsoset(lbj - 1) + 1
562 349880 : lnlb = nsoset(lbj)
563 349880 : fsgfb = first_sgfb(jshellb)
564 349880 : lsgfb = fsgfb + 2*lbj
565 1572454 : DO ishella = 1, nshella
566 1108792 : lai = la(ishella)
567 1108792 : fnla = nsoset(lai - 1) + 1
568 1108792 : lnla = nsoset(lai)
569 1108792 : fsgfa = first_sgfa(ishella)
570 1108792 : lsgfa = fsgfa + 2*lai
571 1108792 : labmin = MIN(lai, lbj)
572 3272555 : DO j = 0, labmin
573 1813883 : prefac = swork_cont(lai + lbj - j + 1, ishella, jshellb)
574 1813883 : dprefac = swork_cont(lai + lbj - j + 2, ishella, jshellb) !j+1
575 8364324 : 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 124504580 : + prefac*dWaux_mat(i, j + 1, fnla:lnla, fnlb:lnlb)
579 : END DO
580 : END DO
581 : END DO
582 : END DO
583 :
584 113782 : 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 69930 : SUBROUTINE construct_overlap_shg_aba(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
606 69930 : lca, first_sgfca, nshellca, cg_coeff, cg_none0_list, &
607 69930 : 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 246667 : DO kshella = 1, nshellca
628 176737 : lak = lca(kshella)
629 176737 : sgfca = first_sgfca(kshella)
630 176737 : ka = sgfca + lak
631 1070985 : DO jshellb = 1, nshellb
632 824318 : lbj = lb(jshellb)
633 824318 : nlb = nsoset(lbj - 1) + lbj + 1
634 824318 : sgfb = first_sgfb(jshellb)
635 824318 : jb = sgfb + lbj
636 4989601 : DO ishella = 1, nshella
637 3988546 : lai = la(ishella)
638 3988546 : sgfa = first_sgfa(ishella)
639 3988546 : ia = sgfa + lai
640 15043850 : DO mai = -lai, lai, 1
641 49875758 : DO mak = -lak, lak, 1
642 35656226 : isoa1 = indso_inv(lai, mai)
643 35656226 : isoa2 = indso_inv(lak, mak)
644 136770380 : DO mbj = -lbj, lbj, 1
645 312211669 : DO ilist = 1, ncg_none0(isoa1, isoa2)
646 185672275 : isoaa = cg_none0_list(isoa1, isoa2, ilist)
647 185672275 : laa = indso(1, isoaa)
648 185672275 : maa = indso(2, isoaa)
649 185672275 : nla = nsoset(laa - 1) + laa + 1
650 185672275 : labmin = MIN(laa, lbj)
651 185672275 : il = INT((lai + lak - laa)/2)
652 185672275 : stemp = 0.0_dp
653 574909559 : DO j = 0, labmin
654 389237284 : prefac = swork_cont(laa + lbj - j + 1, il, ishella, jshellb, kshella)
655 574909559 : stemp = stemp + prefac*Waux_mat(j + 1, nla + maa, nlb + mbj)
656 : END DO
657 276555443 : 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 69930 : 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 53452 : SUBROUTINE dev_overlap_shg_aba(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
691 106904 : lca, first_sgfca, nshellca, cg_coeff, cg_none0_list, &
692 53452 : 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 213808 : rabx2(:) = 2.0_dp*rab
716 :
717 186499 : DO kshella = 1, nshellca
718 133047 : lak = lca(kshella)
719 133047 : sgfca = first_sgfca(kshella)
720 133047 : ka = sgfca + lak
721 841744 : DO jshellb = 1, nshellb
722 655245 : lbj = lb(jshellb)
723 655245 : nlb = nsoset(lbj - 1) + lbj + 1
724 655245 : sgfb = first_sgfb(jshellb)
725 655245 : jb = sgfb + lbj
726 4037919 : DO ishella = 1, nshella
727 3249627 : lai = la(ishella)
728 3249627 : sgfa = first_sgfa(ishella)
729 3249627 : ia = sgfa + lai
730 12336561 : DO mai = -lai, lai, 1
731 40536187 : DO mak = -lak, lak, 1
732 28854871 : isoa1 = indso_inv(lai, mai)
733 28854871 : isoa2 = indso_inv(lak, mak)
734 112147403 : DO mbj = -lbj, lbj, 1
735 255801794 : DO ilist = 1, ncg_none0(isoa1, isoa2)
736 152086080 : isoaa = cg_none0_list(isoa1, isoa2, ilist)
737 152086080 : laa = indso(1, isoaa)
738 152086080 : maa = indso(2, isoaa)
739 152086080 : nla = nsoset(laa - 1) + laa + 1
740 152086080 : labmin = MIN(laa, lbj)
741 152086080 : il = (lai + lak - laa)/2 ! lai+lak-laa always even
742 152086080 : dtemp = 0.0_dp
743 473029574 : DO j = 0, labmin
744 320943494 : prefac = swork_cont(laa + lbj - j + 1, il, ishella, jshellb, kshella)
745 320943494 : dprefac = swork_cont(laa + lbj - j + 2, il, ishella, jshellb, kshella)
746 1435860056 : DO i = 1, 3
747 : dtemp(i) = dtemp(i) + rabx2(i)*dprefac*Waux_mat(j + 1, nla + maa, nlb + mbj) &
748 1283773976 : + prefac*dWaux_mat(i, j + 1, nla + maa, nlb + mbj)
749 : END DO
750 : END DO
751 683205163 : DO i = 1, 3
752 : dsaba(ia + mai, jb + mbj, ka + mak, i) = dsaba(ia + mai, jb + mbj, ka + mak, i) &
753 608344320 : + 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 53452 : 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 3349 : SUBROUTINE construct_overlap_shg_abb(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
785 3349 : lcb, first_sgfcb, nshellcb, cg_coeff, cg_none0_list, &
786 3349 : 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 15158 : DO kshellb = 1, nshellcb
807 11809 : lbk = lcb(kshellb)
808 11809 : sgfcb = first_sgfcb(kshellb)
809 11809 : kb = sgfcb + lbk
810 63289 : DO jshellb = 1, nshellb
811 48131 : lbj = lb(jshellb)
812 48131 : sgfb = first_sgfb(jshellb)
813 48131 : jb = sgfb + lbj
814 214794 : DO ishella = 1, nshella
815 154854 : lai = la(ishella)
816 154854 : nla = nsoset(lai - 1) + lai + 1
817 154854 : sgfa = first_sgfa(ishella)
818 154854 : ia = sgfa + lai
819 572645 : DO mbj = -lbj, lbj, 1
820 2185584 : DO mbk = -lbk, lbk, 1
821 1661070 : isob1 = indso_inv(lbj, mbj)
822 1661070 : isob2 = indso_inv(lbk, mbk)
823 5009744 : DO mai = -lai, lai, 1
824 11411921 : DO ilist = 1, ncg_none0(isob1, isob2)
825 6771837 : isobb = cg_none0_list(isob1, isob2, ilist)
826 6771837 : lbb = indso(1, isobb)
827 6771837 : mbb = indso(2, isobb)
828 6771837 : nlb = nsoset(lbb - 1) + lbb + 1
829 : ! tsgin: because we take the transpose of auxmat (calculated for (la,lb))
830 6771837 : tsign = 1.0_dp
831 6771837 : IF (MODULO(lbb - lai, 2) /= 0) tsign = -1.0_dp
832 6771837 : labmin = MIN(lai, lbb)
833 6771837 : il = INT((lbj + lbk - lbb)/2)
834 6771837 : stemp = 0.0_dp
835 18284074 : DO j = 0, labmin
836 11512237 : prefac = swork_cont(lai + lbb - j + 1, il, ishella, jshellb, kshellb)
837 18284074 : stemp = stemp + prefac*Waux_mat(j + 1, nlb + mbb, nla + mai)
838 : END DO
839 9750851 : 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 3349 : 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 1039 : SUBROUTINE dev_overlap_shg_abb(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
873 2078 : lcb, first_sgfcb, nshellcb, cg_coeff, cg_none0_list, &
874 1039 : 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 4156 : rabx2(:) = 2.0_dp*rab
899 :
900 4000 : DO kshellb = 1, nshellcb
901 2961 : lbk = lcb(kshellb)
902 2961 : sgfcb = first_sgfcb(kshellb)
903 2961 : kb = sgfcb + lbk
904 13207 : DO jshellb = 1, nshellb
905 9207 : lbj = lb(jshellb)
906 9207 : sgfb = first_sgfb(jshellb)
907 9207 : jb = sgfb + lbj
908 37581 : DO ishella = 1, nshella
909 25413 : lai = la(ishella)
910 25413 : nla = nsoset(lai - 1) + lai + 1
911 25413 : sgfa = first_sgfa(ishella)
912 25413 : ia = sgfa + lai
913 99495 : DO mbj = -lbj, lbj, 1
914 379793 : DO mbk = -lbk, lbk, 1
915 289505 : isob1 = indso_inv(lbj, mbj)
916 289505 : isob2 = indso_inv(lbk, mbk)
917 852145 : DO mai = -lai, lai, 1
918 2075895 : DO ilist = 1, ncg_none0(isob1, isob2)
919 1288625 : isobb = cg_none0_list(isob1, isob2, ilist)
920 1288625 : lbb = indso(1, isobb)
921 1288625 : mbb = indso(2, isobb)
922 1288625 : nlb = nsoset(lbb - 1) + lbb + 1
923 : ! tsgin: because we take the transpose of auxmat (calculated for (la,lb))
924 1288625 : tsign = 1.0_dp
925 1288625 : IF (MODULO(lbb - lai, 2) /= 0) tsign = -1.0_dp
926 1288625 : labmin = MIN(lai, lbb)
927 1288625 : il = (lbj + lbk - lbb)/2
928 1288625 : dtemp = 0.0_dp
929 3629345 : DO j = 0, labmin
930 2340720 : prefac = swork_cont(lai + lbb - j + 1, il, ishella, jshellb, kshellb)
931 2340720 : dprefac = swork_cont(lai + lbb - j + 2, il, ishella, jshellb, kshellb)
932 10651505 : DO i = 1, 3
933 : dtemp(i) = dtemp(i) + rabx2(i)*dprefac*Waux_mat(j + 1, nlb + mbb, nla + mai) &
934 9362880 : + prefac*dWaux_mat(i, j + 1, nlb + mbb, nla + mai)
935 : END DO
936 : END DO
937 5652265 : DO i = 1, 3
938 : dsabb(ia + mai, jb + mbj, kb + mbk, i) = dsabb(ia + mai, jb + mbj, kb + mbk, i) &
939 5154500 : + 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 1039 : END SUBROUTINE dev_overlap_shg_abb
950 :
951 : END MODULE construct_shg
952 :
|
