Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Initialization for solid harmonic Gaussian (SHG) integral scheme. Scheme for calculation
10 : !> of contracted, spherical Gaussian integrals using the solid harmonics. Initialization of
11 : !> the contraction matrices
12 : !> \par History
13 : !> created [08.2016]
14 : !> \author Dorothea Golze
15 : ! **************************************************************************************************
16 : MODULE generic_shg_integrals_init
17 : USE basis_set_types, ONLY: gto_basis_set_type
18 : USE kinds, ONLY: dp
19 : USE mathconstants, ONLY: fac,&
20 : ifac,&
21 : pi
22 : USE memory_utilities, ONLY: reallocate
23 : USE orbital_pointers, ONLY: indso,&
24 : nsoset
25 : USE spherical_harmonics, ONLY: clebsch_gordon,&
26 : clebsch_gordon_deallocate,&
27 : clebsch_gordon_init
28 : #include "../base/base_uses.f90"
29 :
30 : IMPLICIT NONE
31 :
32 : PRIVATE
33 :
34 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_shg_integrals_init'
35 :
36 : PUBLIC :: contraction_matrix_shg, contraction_matrix_shg_mix, contraction_matrix_shg_rx2m, &
37 : get_clebsch_gordon_coefficients
38 :
39 : ! **************************************************************************************************
40 :
41 : CONTAINS
42 :
43 : ! **************************************************************************************************
44 : !> \brief contraction matrix for SHG integrals
45 : !> \param basis ...
46 : !> \param scon_shg contraction matrix
47 : ! **************************************************************************************************
48 5055952 : SUBROUTINE contraction_matrix_shg(basis, scon_shg)
49 :
50 : TYPE(gto_basis_set_type), POINTER :: basis
51 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_shg
52 :
53 : INTEGER :: ipgf, iset, ishell, l, maxpgf, maxshell, &
54 : nset
55 5055952 : INTEGER, DIMENSION(:), POINTER :: npgf, nshell
56 : REAL(KIND=dp) :: aif, gcc
57 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: norm
58 5055952 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
59 :
60 5055952 : nset = basis%nset
61 5055952 : npgf => basis%npgf
62 5055952 : nshell => basis%nshell
63 5055952 : zet => basis%zet
64 :
65 5055952 : maxpgf = SIZE(basis%gcc, 1)
66 5055952 : maxshell = SIZE(basis%gcc, 2)
67 20223808 : ALLOCATE (norm(basis%nset, maxshell))
68 25279760 : ALLOCATE (scon_shg(maxpgf, maxshell, nset))
69 34210144 : scon_shg = 0.0_dp
70 :
71 5055952 : CALL basis_norm_shg(basis, norm)
72 :
73 14770082 : DO iset = 1, nset
74 24486870 : DO ishell = 1, nshell(iset)
75 9716788 : l = basis%l(ishell, iset)
76 29150906 : DO ipgf = 1, npgf(iset)
77 9719988 : aif = 1.0_dp/((2._dp*zet(ipgf, iset))**l)
78 9719988 : gcc = basis%gcc(ipgf, ishell, iset)
79 19436776 : scon_shg(ipgf, ishell, iset) = norm(iset, ishell)*gcc*aif
80 : END DO
81 : END DO
82 : END DO
83 :
84 5055952 : DEALLOCATE (norm)
85 :
86 5055952 : END SUBROUTINE contraction_matrix_shg
87 :
88 : !***************************************************************************************************
89 : !> \brief normalization solid harmonic Gaussians (SHG)
90 : !> \param basis ...
91 : !> \param norm ...
92 : ! **************************************************************************************************
93 5056216 : SUBROUTINE basis_norm_shg(basis, norm)
94 :
95 : TYPE(gto_basis_set_type), POINTER :: basis
96 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: norm
97 :
98 : INTEGER :: ipgf, iset, ishell, jpgf, l
99 : REAL(KIND=dp) :: aai, aaj, cci, ccj, expa, ppl
100 :
101 19836984 : norm = 0._dp
102 :
103 14771958 : DO iset = 1, basis%nset
104 24493054 : DO ishell = 1, basis%nshell(iset)
105 9721096 : l = basis%l(ishell, iset)
106 9721096 : expa = 0.5_dp*REAL(2*l + 3, dp)
107 9721096 : ppl = fac(2*l + 2)*pi**(1.5_dp)/fac(l + 1)
108 9721096 : ppl = ppl/(2._dp**REAL(2*l + 1, dp))
109 9721096 : ppl = ppl/REAL(2*l + 1, dp)
110 19448592 : DO ipgf = 1, basis%npgf(iset)
111 9727496 : cci = basis%gcc(ipgf, ishell, iset)
112 9727496 : aai = basis%zet(ipgf, iset)
113 29218568 : DO jpgf = 1, basis%npgf(iset)
114 9769976 : ccj = basis%gcc(jpgf, ishell, iset)
115 9769976 : aaj = basis%zet(jpgf, iset)
116 19497472 : norm(iset, ishell) = norm(iset, ishell) + cci*ccj*ppl/(aai + aaj)**expa
117 : END DO
118 : END DO
119 19436838 : norm(iset, ishell) = 1.0_dp/SQRT(norm(iset, ishell))
120 : END DO
121 : END DO
122 :
123 5056216 : END SUBROUTINE basis_norm_shg
124 :
125 : ! **************************************************************************************************
126 : !> \brief mixed contraction matrix for SHG integrals [aba] and [abb] for orbital and ri basis
127 : !> at the same atom
128 : !> \param orb_basis orbital basis
129 : !> \param ri_basis ...
130 : !> \param orb_index index for orbital basis
131 : !> \param ri_index index for ri basis
132 : !> \param scon_mix mixed contraction matrix
133 : ! **************************************************************************************************
134 132 : SUBROUTINE contraction_matrix_shg_mix(orb_basis, ri_basis, orb_index, ri_index, scon_mix)
135 :
136 : TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
137 : INTEGER, DIMENSION(:, :, :), POINTER :: orb_index, ri_index
138 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: scon_mix
139 :
140 : INTEGER :: forb, fri, iil, il, ipgf, iset, ishell, jpgf, jset, jshell, l, l1, l2, lmax_orb, &
141 : lmax_ri, maxpgf_orb, maxpgf_ri, maxshell_orb, maxshell_ri, nf_orb, nf_ri, nl, nl_max, &
142 : nset_orb, nset_ri
143 132 : INTEGER, DIMENSION(:), POINTER :: npgf_orb, npgf_ri, nshell_orb, nshell_ri
144 : REAL(KIND=dp) :: cjf, const, const1, const2, gcc_orb, &
145 : gcc_ri, prefac, scon1, scon2, zet
146 132 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: shg_fac
147 132 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: norm_orb, norm_ri, zet_orb, zet_ri
148 :
149 132 : nset_orb = orb_basis%nset
150 132 : npgf_orb => orb_basis%npgf
151 132 : nshell_orb => orb_basis%nshell
152 132 : zet_orb => orb_basis%zet
153 132 : nset_ri = ri_basis%nset
154 132 : npgf_ri => ri_basis%npgf
155 132 : nshell_ri => ri_basis%nshell
156 132 : zet_ri => ri_basis%zet
157 :
158 132 : maxpgf_orb = SIZE(orb_basis%gcc, 1)
159 132 : maxshell_orb = SIZE(orb_basis%gcc, 2)
160 528 : ALLOCATE (norm_orb(nset_orb, maxshell_orb))
161 132 : maxpgf_ri = SIZE(ri_basis%gcc, 1)
162 132 : maxshell_ri = SIZE(ri_basis%gcc, 2)
163 528 : ALLOCATE (norm_ri(nset_ri, maxshell_ri))
164 :
165 132 : CALL basis_norm_shg(orb_basis, norm_orb)
166 132 : CALL basis_norm_shg(ri_basis, norm_ri)
167 :
168 660 : ALLOCATE (orb_index(maxpgf_orb, maxshell_orb, nset_orb))
169 660 : ALLOCATE (ri_index(maxpgf_ri, maxshell_ri, nset_ri))
170 :
171 : !** index orbital basis set
172 132 : nf_orb = 0
173 308 : DO iset = 1, nset_orb
174 754 : DO ishell = 1, nshell_orb(iset)
175 3504 : DO ipgf = 1, npgf_orb(iset)
176 2882 : nf_orb = nf_orb + 1
177 3328 : orb_index(ipgf, ishell, iset) = nf_orb
178 : END DO
179 : END DO
180 : END DO
181 :
182 : !** index ri basis set
183 : nf_ri = 0
184 1568 : DO iset = 1, nset_ri
185 5430 : DO ishell = 1, nshell_ri(iset)
186 9924 : DO ipgf = 1, npgf_ri(iset)
187 4626 : nf_ri = nf_ri + 1
188 8488 : ri_index(ipgf, ishell, iset) = nf_ri
189 : END DO
190 : END DO
191 : END DO
192 :
193 308 : lmax_orb = MAXVAL(orb_basis%lmax)
194 1568 : lmax_ri = MAXVAL(ri_basis%lmax)
195 132 : nl_max = INT((lmax_orb + lmax_ri)/2) + 1
196 792 : ALLOCATE (scon_mix(nl_max, nf_ri, nf_orb, nl_max))
197 1880530 : scon_mix = 0.0_dp
198 :
199 396 : ALLOCATE (shg_fac(0:nl_max - 1))
200 132 : shg_fac(0) = 1.0_dp
201 :
202 308 : DO iset = 1, nset_orb
203 754 : DO ishell = 1, nshell_orb(iset)
204 446 : l1 = orb_basis%l(ishell, iset)
205 446 : const1 = SQRT(1.0_dp/REAL(2*l1 + 1, dp))
206 5784 : DO jset = 1, nset_ri
207 20028 : DO jshell = 1, nshell_ri(jset)
208 14420 : l2 = ri_basis%l(jshell, jset)
209 14420 : const2 = SQRT(1.0_dp/REAL(2*l2 + 1, dp))
210 14420 : nl = INT((l1 + l2)/2)
211 14420 : IF (l1 == 0 .OR. l2 == 0) nl = 0
212 41756 : DO il = 0, nl
213 22174 : l = l1 + l2 - 2*il
214 22174 : const = const1*const2*2.0_dp*SQRT(pi*REAL(2*l + 1, dp))
215 32646 : DO iil = 1, il
216 : shg_fac(iil) = fac(l + iil - 1)*ifac(l)*REAL(l, dp) &
217 32646 : *fac(il)/fac(il - iil)/fac(iil)
218 : END DO
219 187010 : DO ipgf = 1, npgf_orb(iset)
220 150416 : forb = orb_index(ipgf, ishell, iset)
221 150416 : gcc_orb = orb_basis%gcc(ipgf, ishell, iset)
222 150416 : scon1 = norm_orb(iset, ishell)*gcc_orb
223 347000 : DO jpgf = 1, npgf_ri(jset)
224 174410 : fri = ri_index(jpgf, jshell, jset)
225 174410 : gcc_ri = ri_basis%gcc(jpgf, jshell, jset)
226 174410 : scon2 = norm_ri(jset, jshell)*gcc_ri
227 174410 : zet = zet_orb(ipgf, iset) + zet_ri(jpgf, jset)
228 174410 : cjf = 1.0_dp/((2._dp*zet)**l)
229 174410 : prefac = const*cjf*scon1*scon2
230 577228 : DO iil = 0, il
231 426812 : scon_mix(iil + 1, fri, forb, il + 1) = prefac*shg_fac(iil)/zet**iil
232 : END DO
233 : END DO
234 : END DO
235 : END DO
236 : END DO
237 : END DO
238 : END DO
239 : END DO
240 :
241 132 : DEALLOCATE (norm_orb, norm_ri, shg_fac)
242 :
243 132 : END SUBROUTINE contraction_matrix_shg_mix
244 :
245 : ! **************************************************************************************************
246 : !> \brief ...
247 : !> \param basis ...
248 : !> \param m ...
249 : !> \param scon_shg ...
250 : !> \param scon_rx2m ...
251 : ! **************************************************************************************************
252 2 : SUBROUTINE contraction_matrix_shg_rx2m(basis, m, scon_shg, scon_rx2m)
253 :
254 : TYPE(gto_basis_set_type), POINTER :: basis
255 : INTEGER, INTENT(IN) :: m
256 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scon_shg
257 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: scon_rx2m
258 :
259 : INTEGER :: ipgf, iset, ishell, j, l, maxpgf, &
260 : maxshell, nset
261 2 : INTEGER, DIMENSION(:), POINTER :: npgf, nshell
262 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: shg_fac
263 2 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
264 :
265 2 : npgf => basis%npgf
266 2 : nshell => basis%nshell
267 2 : zet => basis%zet
268 2 : nset = basis%nset
269 :
270 2 : maxpgf = SIZE(basis%gcc, 1)
271 2 : maxshell = SIZE(basis%gcc, 2)
272 12 : ALLOCATE (scon_rx2m(maxpgf, m + 1, maxshell, nset))
273 742 : scon_rx2m = 0.0_dp
274 6 : ALLOCATE (shg_fac(0:m))
275 2 : shg_fac(0) = 1.0_dp
276 :
277 22 : DO iset = 1, nset
278 88 : DO ishell = 1, nshell(iset)
279 66 : l = basis%l(ishell, iset)
280 264 : DO j = 1, m
281 : shg_fac(j) = fac(l + j - 1)*ifac(l)*REAL(l, dp) &
282 264 : *fac(m)/fac(m - j)/fac(j)
283 : END DO
284 152 : DO ipgf = 1, npgf(iset)
285 396 : DO j = 0, m
286 : scon_rx2m(ipgf, j + 1, ishell, iset) = scon_shg(ipgf, ishell, iset) &
287 330 : *shg_fac(j)/zet(ipgf, iset)**j
288 : END DO
289 : END DO
290 : END DO
291 : END DO
292 :
293 2 : DEALLOCATE (shg_fac)
294 :
295 2 : END SUBROUTINE contraction_matrix_shg_rx2m
296 :
297 : ! **************************************************************************************************
298 : !> \brief calculate the Clebsch-Gordon (CG) coefficients for expansion of the
299 : !> product of two spherical harmonic Gaussians
300 : !> \param my_cg matrix storing CG coefficients
301 : !> \param cg_none0_list list of none-zero CG coefficients
302 : !> \param ncg_none0 number of none-zero CG coefficients
303 : !> \param maxl1 maximal l quantum number of 1st spherical function
304 : !> \param maxl2 maximal l quantum number of 2nd spherical function
305 : ! **************************************************************************************************
306 62 : SUBROUTINE get_clebsch_gordon_coefficients(my_cg, cg_none0_list, ncg_none0, &
307 : maxl1, maxl2)
308 :
309 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: my_cg
310 : INTEGER, DIMENSION(:, :, :), POINTER :: cg_none0_list
311 : INTEGER, DIMENSION(:, :), POINTER :: ncg_none0
312 : INTEGER, INTENT(IN) :: maxl1, maxl2
313 :
314 : INTEGER :: il, ilist, iso, iso1, iso2, l1, l1l2, &
315 : l2, lc1, lc2, lp, m1, m2, maxl, mm, &
316 : mp, nlist, nlist_max, nsfunc, nsfunc1, &
317 : nsfunc2
318 62 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rga
319 :
320 62 : nlist_max = 6
321 62 : nsfunc1 = nsoset(maxl1)
322 62 : nsfunc2 = nsoset(maxl2)
323 62 : maxl = maxl1 + maxl2
324 62 : nsfunc = nsoset(maxl)
325 :
326 62 : CALL clebsch_gordon_init(maxl)
327 :
328 310 : ALLOCATE (my_cg(nsfunc1, nsfunc2, nsfunc))
329 791168 : my_cg = 0.0_dp
330 248 : ALLOCATE (ncg_none0(nsfunc1, nsfunc2))
331 13500 : ncg_none0 = 0
332 310 : ALLOCATE (cg_none0_list(nsfunc1, nsfunc2, nlist_max))
333 81062 : cg_none0_list = 0
334 :
335 186 : ALLOCATE (rga(maxl, 2))
336 850 : rga = 0.0_dp
337 238 : DO lc1 = 0, maxl1
338 758 : DO iso1 = nsoset(lc1 - 1) + 1, nsoset(lc1)
339 520 : l1 = indso(1, iso1)
340 520 : m1 = indso(2, iso1)
341 3150 : DO lc2 = 0, maxl2
342 15104 : DO iso2 = nsoset(lc2 - 1) + 1, nsoset(lc2)
343 12130 : l2 = indso(1, iso2)
344 12130 : m2 = indso(2, iso2)
345 12130 : CALL clebsch_gordon(l1, m1, l2, m2, rga)
346 12130 : l1l2 = l1 + l2
347 12130 : mp = m1 + m2
348 12130 : mm = m1 - m2
349 12130 : IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN
350 5660 : mp = -ABS(mp)
351 5660 : mm = -ABS(mm)
352 : ELSE
353 6470 : mp = ABS(mp)
354 6470 : mm = ABS(mm)
355 : END IF
356 49088 : DO lp = MOD(l1 + l2, 2), l1l2, 2
357 36958 : il = lp/2 + 1
358 36958 : IF (ABS(mp) <= lp) THEN
359 25182 : IF (mp >= 0) THEN
360 14894 : iso = nsoset(lp - 1) + lp + 1 + mp
361 : ELSE
362 10288 : iso = nsoset(lp - 1) + lp + 1 - ABS(mp)
363 : END IF
364 25182 : my_cg(iso1, iso2, iso) = rga(il, 1)
365 : END IF
366 49088 : IF (mp /= mm .AND. ABS(mm) <= lp) THEN
367 14588 : IF (mm >= 0) THEN
368 9416 : iso = nsoset(lp - 1) + lp + 1 + mm
369 : ELSE
370 5172 : iso = nsoset(lp - 1) + lp + 1 - ABS(mm)
371 : END IF
372 14588 : my_cg(iso1, iso2, iso) = rga(il, 2)
373 : END IF
374 : END DO
375 12130 : nlist = 0
376 737460 : DO ilist = 1, nsfunc
377 737460 : IF (ABS(my_cg(iso1, iso2, ilist)) > 1.E-8_dp) THEN
378 34354 : nlist = nlist + 1
379 34354 : IF (nlist > nlist_max) THEN
380 8 : CALL reallocate(cg_none0_list, 1, nsfunc1, 1, nsfunc2, 1, nlist)
381 8 : nlist_max = nlist
382 : END IF
383 34354 : cg_none0_list(iso1, iso2, nlist) = ilist
384 : END IF
385 : END DO
386 14584 : ncg_none0(iso1, iso2) = nlist
387 : END DO ! iso2
388 : END DO ! lc2
389 : END DO ! iso1
390 : END DO ! lc1
391 :
392 62 : DEALLOCATE (rga)
393 62 : CALL clebsch_gordon_deallocate()
394 :
395 62 : END SUBROUTINE get_clebsch_gordon_coefficients
396 :
397 : END MODULE generic_shg_integrals_init
|
