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 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 14388254 : 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 14388254 : INTEGER, DIMENSION(:), POINTER :: npgf, nshell
56 : REAL(KIND=dp) :: aif, gcc
57 14388254 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: norm
58 14388254 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
59 :
60 14388254 : nset = basis%nset
61 14388254 : npgf => basis%npgf
62 14388254 : nshell => basis%nshell
63 14388254 : zet => basis%zet
64 :
65 14388254 : maxpgf = SIZE(basis%gcc, 1)
66 14388254 : maxshell = SIZE(basis%gcc, 2)
67 57553016 : ALLOCATE (norm(basis%nset, maxshell))
68 71941270 : ALLOCATE (scon_shg(maxpgf, maxshell, nset))
69 80896824 : scon_shg = 0.0_dp
70 :
71 14388254 : CALL basis_norm_shg(basis, norm)
72 :
73 36553374 : DO iset = 1, nset
74 58721484 : DO ishell = 1, nshell(iset)
75 22168110 : l = basis%l(ishell, iset)
76 66504964 : DO ipgf = 1, npgf(iset)
77 22171734 : aif = 1.0_dp/((2._dp*zet(ipgf, iset))**l)
78 22171734 : gcc = basis%gcc(ipgf, ishell, iset)
79 44339844 : scon_shg(ipgf, ishell, iset) = norm(iset, ishell)*gcc*aif
80 : END DO
81 : END DO
82 : END DO
83 :
84 14388254 : DEALLOCATE (norm)
85 :
86 14388254 : END SUBROUTINE contraction_matrix_shg
87 :
88 : !***************************************************************************************************
89 : !> \brief normalization solid harmonic Gaussians (SHG)
90 : !> \param basis ...
91 : !> \param norm ...
92 : ! **************************************************************************************************
93 14388554 : 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 50954126 : norm = 0._dp
102 :
103 36555610 : DO iset = 1, basis%nset
104 58728684 : DO ishell = 1, basis%nshell(iset)
105 22173074 : l = basis%l(ishell, iset)
106 22173074 : expa = 0.5_dp*REAL(2*l + 3, dp)
107 22173074 : ppl = fac(2*l + 2)*pi**(1.5_dp)/fac(l + 1)
108 22173074 : ppl = ppl/(2._dp**REAL(2*l + 1, dp))
109 22173074 : ppl = ppl/REAL(2*l + 1, dp)
110 44353396 : DO ipgf = 1, basis%npgf(iset)
111 22180322 : cci = basis%gcc(ipgf, ishell, iset)
112 22180322 : aai = basis%zet(ipgf, iset)
113 66581974 : DO jpgf = 1, basis%npgf(iset)
114 22228578 : ccj = basis%gcc(jpgf, ishell, iset)
115 22228578 : aaj = basis%zet(jpgf, iset)
116 44408900 : norm(iset, ishell) = norm(iset, ishell) + cci*ccj*ppl/(aai + aaj)**expa
117 : END DO
118 : END DO
119 44340130 : norm(iset, ishell) = 1.0_dp/SQRT(norm(iset, ishell))
120 : END DO
121 : END DO
122 :
123 14388554 : 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 150 : 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 150 : 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 150 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: shg_fac
147 150 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: norm_orb, norm_ri, zet_orb, zet_ri
148 :
149 150 : nset_orb = orb_basis%nset
150 150 : npgf_orb => orb_basis%npgf
151 150 : nshell_orb => orb_basis%nshell
152 150 : zet_orb => orb_basis%zet
153 150 : nset_ri = ri_basis%nset
154 150 : npgf_ri => ri_basis%npgf
155 150 : nshell_ri => ri_basis%nshell
156 150 : zet_ri => ri_basis%zet
157 :
158 150 : maxpgf_orb = SIZE(orb_basis%gcc, 1)
159 150 : maxshell_orb = SIZE(orb_basis%gcc, 2)
160 600 : ALLOCATE (norm_orb(nset_orb, maxshell_orb))
161 150 : maxpgf_ri = SIZE(ri_basis%gcc, 1)
162 150 : maxshell_ri = SIZE(ri_basis%gcc, 2)
163 600 : ALLOCATE (norm_ri(nset_ri, maxshell_ri))
164 :
165 150 : CALL basis_norm_shg(orb_basis, norm_orb)
166 150 : CALL basis_norm_shg(ri_basis, norm_ri)
167 :
168 750 : ALLOCATE (orb_index(maxpgf_orb, maxshell_orb, nset_orb))
169 750 : ALLOCATE (ri_index(maxpgf_ri, maxshell_ri, nset_ri))
170 :
171 : !** index orbital basis set
172 150 : nf_orb = 0
173 344 : DO iset = 1, nset_orb
174 864 : DO ishell = 1, nshell_orb(iset)
175 4094 : DO ipgf = 1, npgf_orb(iset)
176 3380 : nf_orb = nf_orb + 1
177 3900 : 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 1892 : DO iset = 1, nset_ri
185 6336 : DO ishell = 1, nshell_ri(iset)
186 11394 : DO ipgf = 1, npgf_ri(iset)
187 5208 : nf_ri = nf_ri + 1
188 9652 : ri_index(ipgf, ishell, iset) = nf_ri
189 : END DO
190 : END DO
191 : END DO
192 :
193 344 : lmax_orb = MAXVAL(orb_basis%lmax)
194 1892 : lmax_ri = MAXVAL(ri_basis%lmax)
195 150 : nl_max = INT((lmax_orb + lmax_ri)/2) + 1
196 900 : ALLOCATE (scon_mix(nl_max, nf_ri, nf_orb, nl_max))
197 2043224 : scon_mix = 0.0_dp
198 :
199 450 : ALLOCATE (shg_fac(0:nl_max - 1))
200 150 : shg_fac(0) = 1.0_dp
201 :
202 344 : DO iset = 1, nset_orb
203 864 : DO ishell = 1, nshell_orb(iset)
204 520 : l1 = orb_basis%l(ishell, iset)
205 520 : const1 = SQRT(1.0_dp/REAL(2*l1 + 1, dp))
206 7054 : DO jset = 1, nset_ri
207 23662 : DO jshell = 1, nshell_ri(jset)
208 16802 : l2 = ri_basis%l(jshell, jset)
209 16802 : const2 = SQRT(1.0_dp/REAL(2*l2 + 1, dp))
210 16802 : nl = INT((l1 + l2)/2)
211 16802 : IF (l1 == 0 .OR. l2 == 0) nl = 0
212 48748 : DO il = 0, nl
213 25606 : l = l1 + l2 - 2*il
214 25606 : const = const1*const2*2.0_dp*SQRT(pi*REAL(2*l + 1, dp))
215 37330 : DO iil = 1, il
216 : shg_fac(iil) = fac(l + iil - 1)*ifac(l)*REAL(l, dp) &
217 37330 : *fac(il)/fac(il - iil)/fac(iil)
218 : END DO
219 216144 : DO ipgf = 1, npgf_orb(iset)
220 173736 : forb = orb_index(ipgf, ishell, iset)
221 173736 : gcc_orb = orb_basis%gcc(ipgf, ishell, iset)
222 173736 : scon1 = norm_orb(iset, ishell)*gcc_orb
223 397072 : DO jpgf = 1, npgf_ri(jset)
224 197730 : fri = ri_index(jpgf, jshell, jset)
225 197730 : gcc_ri = ri_basis%gcc(jpgf, jshell, jset)
226 197730 : scon2 = norm_ri(jset, jshell)*gcc_ri
227 197730 : zet = zet_orb(ipgf, iset) + zet_ri(jpgf, jset)
228 197730 : cjf = 1.0_dp/((2._dp*zet)**l)
229 197730 : prefac = const*cjf*scon1*scon2
230 655656 : DO iil = 0, il
231 481920 : 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 150 : DEALLOCATE (norm_orb, norm_ri, shg_fac)
242 :
243 150 : 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 2 : 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 72 : 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 72 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rga
319 :
320 72 : nlist_max = 6
321 72 : nsfunc1 = nsoset(maxl1)
322 72 : nsfunc2 = nsoset(maxl2)
323 72 : maxl = maxl1 + maxl2
324 72 : nsfunc = nsoset(maxl)
325 :
326 72 : CALL clebsch_gordon_init(maxl)
327 :
328 360 : ALLOCATE (my_cg(nsfunc1, nsfunc2, nsfunc))
329 849138 : my_cg = 0.0_dp
330 288 : ALLOCATE (ncg_none0(nsfunc1, nsfunc2))
331 15110 : ncg_none0 = 0
332 360 : ALLOCATE (cg_none0_list(nsfunc1, nsfunc2, nlist_max))
333 90732 : cg_none0_list = 0
334 :
335 216 : ALLOCATE (rga(maxl, 2))
336 72 : rga = 0.0_dp
337 278 : DO lc1 = 0, maxl1
338 888 : DO iso1 = nsoset(lc1 - 1) + 1, nsoset(lc1)
339 610 : l1 = indso(1, iso1)
340 610 : m1 = indso(2, iso1)
341 3630 : DO lc2 = 0, maxl2
342 16994 : DO iso2 = nsoset(lc2 - 1) + 1, nsoset(lc2)
343 13570 : l2 = indso(1, iso2)
344 13570 : m2 = indso(2, iso2)
345 13570 : CALL clebsch_gordon(l1, m1, l2, m2, rga)
346 13570 : l1l2 = l1 + l2
347 13570 : mp = m1 + m2
348 13570 : mm = m1 - m2
349 13570 : IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN
350 6320 : mp = -ABS(mp)
351 6320 : mm = -ABS(mm)
352 : ELSE
353 7250 : mp = ABS(mp)
354 7250 : mm = ABS(mm)
355 : END IF
356 54148 : DO lp = MOD(l1 + l2, 2), l1l2, 2
357 40578 : il = lp/2 + 1
358 40578 : IF (ABS(mp) <= lp) THEN
359 27742 : IF (mp >= 0) THEN
360 16504 : iso = nsoset(lp - 1) + lp + 1 + mp
361 : ELSE
362 11238 : iso = nsoset(lp - 1) + lp + 1 - ABS(mp)
363 : END IF
364 27742 : my_cg(iso1, iso2, iso) = rga(il, 1)
365 : END IF
366 54148 : IF (mp /= mm .AND. ABS(mm) <= lp) THEN
367 15928 : IF (mm >= 0) THEN
368 10296 : iso = nsoset(lp - 1) + lp + 1 + mm
369 : ELSE
370 5632 : iso = nsoset(lp - 1) + lp + 1 - ABS(mm)
371 : END IF
372 15928 : my_cg(iso1, iso2, iso) = rga(il, 2)
373 : END IF
374 : END DO
375 13570 : nlist = 0
376 790740 : DO ilist = 1, nsfunc
377 790740 : IF (ABS(my_cg(iso1, iso2, ilist)) > 1.E-8_dp) THEN
378 37734 : nlist = nlist + 1
379 37734 : 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 37734 : cg_none0_list(iso1, iso2, nlist) = ilist
384 : END IF
385 : END DO
386 16384 : ncg_none0(iso1, iso2) = nlist
387 : END DO ! iso2
388 : END DO ! lc2
389 : END DO ! iso1
390 : END DO ! lc1
391 :
392 72 : DEALLOCATE (rga)
393 72 : CALL clebsch_gordon_deallocate()
394 :
395 72 : END SUBROUTINE get_clebsch_gordon_coefficients
396 :
397 : END MODULE generic_shg_integrals_init
|
