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 : MODULE qs_harmonics_atom
9 :
10 : USE basis_set_types, ONLY: get_gto_basis_set,&
11 : gto_basis_set_type
12 : USE kinds, ONLY: dp
13 : USE lebedev, ONLY: lebedev_grid
14 : USE memory_utilities, ONLY: reallocate
15 : USE orbital_pointers, ONLY: indco,&
16 : indso,&
17 : nco,&
18 : ncoset,&
19 : nso,&
20 : nsoset
21 : USE orbital_transformation_matrices, ONLY: orbtramat
22 : USE spherical_harmonics, ONLY: dy_lm,&
23 : y_lm
24 : #include "./base/base_uses.f90"
25 :
26 : IMPLICIT NONE
27 :
28 : PRIVATE
29 :
30 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_harmonics_atom'
31 :
32 : TYPE harmonics_atom_type
33 : INTEGER :: max_s_harm = -1, llmax = -1, &
34 : max_iso_not0 = -1, &
35 : dmax_iso_not0 = -1, &
36 : damax_iso_not0 = -1, &
37 : ngrid = -1
38 : REAL(dp), DIMENSION(:, :), POINTER :: a => NULL(), slm => NULL()
39 : REAL(dp), DIMENSION(:, :, :), POINTER :: dslm => NULL(), dslm_dxyz => NULL()
40 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG => NULL()
41 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_CG_dxyz => NULL()
42 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_CG_dxyz_asym => NULL()
43 : REAL(dp), DIMENSION(:), POINTER :: slm_int => NULL()
44 :
45 : END TYPE harmonics_atom_type
46 :
47 : PUBLIC :: allocate_harmonics_atom, &
48 : create_harmonics_atom, &
49 : deallocate_harmonics_atom, &
50 : get_none0_cg_list
51 :
52 : PUBLIC :: harmonics_atom_type, get_maxl_CG
53 :
54 : INTERFACE get_none0_cg_list
55 : MODULE PROCEDURE get_none0_cg_list3, get_none0_cg_list4
56 : END INTERFACE
57 :
58 : CONTAINS
59 :
60 : ! **************************************************************************************************
61 : !> \brief Allocate a spherical harmonics set for the atom grid.
62 : !> \param harmonics ...
63 : !> \version 1.0
64 : ! **************************************************************************************************
65 3038 : SUBROUTINE allocate_harmonics_atom(harmonics)
66 :
67 : TYPE(harmonics_atom_type), POINTER :: harmonics
68 :
69 3038 : IF (ASSOCIATED(harmonics)) CALL deallocate_harmonics_atom(harmonics)
70 :
71 3038 : ALLOCATE (harmonics)
72 :
73 3038 : harmonics%max_s_harm = 0
74 3038 : harmonics%llmax = 0
75 3038 : harmonics%max_iso_not0 = 0
76 3038 : harmonics%dmax_iso_not0 = 0
77 3038 : harmonics%damax_iso_not0 = 0
78 3038 : harmonics%ngrid = 0
79 :
80 : NULLIFY (harmonics%slm)
81 : NULLIFY (harmonics%dslm)
82 : NULLIFY (harmonics%dslm_dxyz)
83 : NULLIFY (harmonics%slm_int)
84 : NULLIFY (harmonics%my_CG)
85 : NULLIFY (harmonics%my_CG_dxyz)
86 : NULLIFY (harmonics%my_CG_dxyz_asym)
87 : NULLIFY (harmonics%a)
88 :
89 3038 : END SUBROUTINE allocate_harmonics_atom
90 :
91 : ! **************************************************************************************************
92 : !> \brief Deallocate the spherical harmonics set for the atom grid.
93 : !> \param harmonics ...
94 : !> \version 1.0
95 : ! **************************************************************************************************
96 3038 : SUBROUTINE deallocate_harmonics_atom(harmonics)
97 :
98 : TYPE(harmonics_atom_type), POINTER :: harmonics
99 :
100 3038 : IF (ASSOCIATED(harmonics)) THEN
101 :
102 3038 : IF (ASSOCIATED(harmonics%slm)) THEN
103 3030 : DEALLOCATE (harmonics%slm)
104 : END IF
105 :
106 3038 : IF (ASSOCIATED(harmonics%dslm)) THEN
107 3030 : DEALLOCATE (harmonics%dslm)
108 : END IF
109 :
110 3038 : IF (ASSOCIATED(harmonics%dslm_dxyz)) THEN
111 3030 : DEALLOCATE (harmonics%dslm_dxyz)
112 : END IF
113 :
114 3038 : IF (ASSOCIATED(harmonics%slm_int)) THEN
115 3030 : DEALLOCATE (harmonics%slm_int)
116 : END IF
117 :
118 3038 : IF (ASSOCIATED(harmonics%my_CG)) THEN
119 3038 : DEALLOCATE (harmonics%my_CG)
120 : END IF
121 :
122 3038 : IF (ASSOCIATED(harmonics%my_CG_dxyz)) THEN
123 3030 : DEALLOCATE (harmonics%my_CG_dxyz)
124 : END IF
125 :
126 3038 : IF (ASSOCIATED(harmonics%my_CG_dxyz_asym)) THEN
127 3030 : DEALLOCATE (harmonics%my_CG_dxyz_asym)
128 : END IF
129 :
130 3038 : IF (ASSOCIATED(harmonics%a)) THEN
131 3030 : DEALLOCATE (harmonics%a)
132 : END IF
133 :
134 3038 : DEALLOCATE (harmonics)
135 : ELSE
136 : CALL cp_abort(__LOCATION__, &
137 : "The pointer harmonics is not associated and "// &
138 0 : "cannot be deallocated")
139 : END IF
140 :
141 3038 : END SUBROUTINE deallocate_harmonics_atom
142 :
143 : ! **************************************************************************************************
144 : !> \brief ...
145 : !> \param harmonics ...
146 : !> \param my_CG ...
147 : !> \param na ...
148 : !> \param llmax ...
149 : !> \param maxs ...
150 : !> \param max_s_harm ...
151 : !> \param ll ...
152 : !> \param wa ...
153 : !> \param azi ...
154 : !> \param pol ...
155 : !> \note Slight refactoring + OMP parallelized (03.2020 A. Bussy)
156 : ! **************************************************************************************************
157 3030 : SUBROUTINE create_harmonics_atom(harmonics, my_CG, na, llmax, maxs, max_s_harm, ll, wa, azi, pol)
158 :
159 : TYPE(harmonics_atom_type), POINTER :: harmonics
160 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
161 : INTEGER, INTENT(IN) :: na, llmax, maxs, max_s_harm, ll
162 : REAL(dp), DIMENSION(:), INTENT(IN) :: wa, azi, pol
163 :
164 : CHARACTER(len=*), PARAMETER :: routineN = 'create_harmonics_atom'
165 :
166 : INTEGER :: handle, i, ia, ic, is, is1, is2, iso, &
167 : iso1, iso2, l, l1, l2, lmax_grid, lx, &
168 : ly, lz, m, m1, m2, max_s_grid, n
169 : REAL(dp) :: drx, dry, drz, rx, ry, rz
170 : REAL(dp), DIMENSION(2) :: cin, dylm
171 3030 : REAL(dp), DIMENSION(:), POINTER :: slm_int, y
172 3030 : REAL(dp), DIMENSION(:, :), POINTER :: dc, slm
173 3030 : REAL(dp), DIMENSION(:, :, :), POINTER :: dslm_dxyz
174 :
175 3030 : CALL timeset(routineN, handle)
176 :
177 3030 : NULLIFY (y, slm, dslm_dxyz, dc)
178 :
179 3030 : CPASSERT(ASSOCIATED(harmonics))
180 :
181 3030 : max_s_grid = MAX(maxs, max_s_harm)
182 3030 : lmax_grid = indso(1, max_s_grid)
183 :
184 3030 : harmonics%max_s_harm = max_s_harm
185 3030 : harmonics%llmax = llmax
186 3030 : harmonics%ngrid = na
187 :
188 3030 : NULLIFY (harmonics%my_CG, harmonics%my_CG_dxyz, harmonics%my_CG_dxyz_asym)
189 3030 : CALL reallocate(harmonics%my_CG, 1, maxs, 1, maxs, 1, max_s_harm)
190 3030 : CALL reallocate(harmonics%my_CG_dxyz, 1, 3, 1, maxs, 1, maxs, 1, max_s_harm)
191 3030 : CALL reallocate(harmonics%my_CG_dxyz_asym, 1, 3, 1, maxs, 1, maxs, 1, max_s_harm)
192 :
193 67772 : DO i = 1, max_s_harm
194 553024 : DO is1 = 1, maxs
195 12876166 : harmonics%my_CG(1:maxs, is1, i) = my_CG(1:maxs, is1, i)
196 : END DO
197 : END DO
198 :
199 : ! allocate and calculate the spherical harmonics LM for this grid
200 : ! and their derivatives
201 3030 : NULLIFY (harmonics%slm, harmonics%dslm, harmonics%dslm_dxyz, harmonics%a, harmonics%slm_int)
202 3030 : CALL reallocate(harmonics%slm, 1, na, 1, max_s_grid)
203 3030 : CALL reallocate(harmonics%dslm, 1, 2, 1, na, 1, maxs)
204 3030 : CALL reallocate(harmonics%dslm_dxyz, 1, 3, 1, na, 1, max_s_grid)
205 3030 : CALL reallocate(harmonics%a, 1, 3, 1, na)
206 3030 : CALL reallocate(harmonics%slm_int, 1, max_s_grid)
207 :
208 : NULLIFY (slm, dslm_dxyz, slm_int)
209 3030 : slm => harmonics%slm
210 3030 : dslm_dxyz => harmonics%dslm_dxyz
211 15058124 : dslm_dxyz = 0.0_dp
212 3030 : slm_int => harmonics%slm_int
213 67772 : slm_int = 0.0_dp
214 :
215 : !$OMP PARALLEL DEFAULT(NONE), &
216 : !$OMP SHARED (slm,dslm_dxyz,slm_int,max_s_harm,max_s_grid,ll,lebedev_grid,na,harmonics,wa,indco,orbtramat) &
217 : !$OMP SHARED (nso,nsoset,nco,maxs,indso,ncoset,pol,azi,llmax,lmax_grid) &
218 : !$OMP PRIVATE(ia,iso,l,m,i,lx,ly,lz,rx,ry,rz,drx,dry,drz,ic,dc,iso1,iso2,cin,dylm) &
219 3030 : !$OMP PRIVATE(is1,l1,m1,is2,l2,m2,is,n,y)
220 :
221 : ALLOCATE (y(na))
222 : !$OMP DO
223 : DO iso = 1, max_s_grid
224 : l = indso(1, iso)
225 : m = indso(2, iso)
226 : CALL y_lm(lebedev_grid(ll)%r, y, l, m)
227 :
228 : DO ia = 1, na
229 : slm(ia, iso) = y(ia)
230 : slm_int(iso) = slm_int(iso) + slm(ia, iso)*wa(ia)
231 : END DO ! ia
232 : END DO ! iso
233 : !$OMP END DO
234 : DEALLOCATE (y)
235 :
236 : !$OMP DO
237 : DO ia = 1, na
238 : harmonics%a(:, ia) = lebedev_grid(ll)%r(:, ia)
239 : END DO
240 : !$OMP END DO
241 :
242 : !
243 : ! The derivatives dslm_dxyz and its expansions my_CG_dxyz and my_CG_dxyz_asymm
244 : ! are NOT the dSlm/dx but the scaled by r**(l-1) derivatives of the monomial
245 : ! terms x^n1 y^n2 z^n3 transformed by spherical harmonics expansion coefficients
246 : !
247 :
248 : ALLOCATE (dc(nco(lmax_grid), 3))
249 : !$OMP DO
250 : DO ia = 1, na
251 : DO l = 0, lmax_grid
252 : DO ic = 1, nco(l)
253 : lx = indco(1, ic + ncoset(l - 1))
254 : ly = indco(2, ic + ncoset(l - 1))
255 : lz = indco(3, ic + ncoset(l - 1))
256 :
257 : IF (lx == 0) THEN
258 : rx = 1.0_dp
259 : drx = 0.0_dp
260 : ELSE IF (lx == 1) THEN
261 : rx = lebedev_grid(ll)%r(1, ia)
262 : drx = 1.0_dp
263 : ELSE
264 : rx = lebedev_grid(ll)%r(1, ia)**lx
265 : drx = REAL(lx, dp)*lebedev_grid(ll)%r(1, ia)**(lx - 1)
266 : END IF
267 : IF (ly == 0) THEN
268 : ry = 1.0_dp
269 : dry = 0.0_dp
270 : ELSE IF (ly == 1) THEN
271 : ry = lebedev_grid(ll)%r(2, ia)
272 : dry = 1.0_dp
273 : ELSE
274 : ry = lebedev_grid(ll)%r(2, ia)**ly
275 : dry = REAL(ly, dp)*lebedev_grid(ll)%r(2, ia)**(ly - 1)
276 : END IF
277 : IF (lz == 0) THEN
278 : rz = 1.0_dp
279 : drz = 0.0_dp
280 : ELSE IF (lz == 1) THEN
281 : rz = lebedev_grid(ll)%r(3, ia)
282 : drz = 1.0_dp
283 : ELSE
284 : rz = lebedev_grid(ll)%r(3, ia)**lz
285 : drz = REAL(lz, dp)*lebedev_grid(ll)%r(3, ia)**(lz - 1)
286 : END IF
287 : dc(ic, 1) = drx*ry*rz
288 : dc(ic, 2) = rx*dry*rz
289 : dc(ic, 3) = rx*ry*drz
290 : END DO
291 : n = nsoset(l - 1)
292 : DO is = 1, nso(l)
293 : iso = is + n
294 : DO ic = 1, nco(l)
295 : dslm_dxyz(:, ia, iso) = dslm_dxyz(:, ia, iso) + &
296 : orbtramat(l)%slm(is, ic)*dc(ic, :)
297 : END DO
298 : END DO
299 : END DO ! l
300 : END DO !ia
301 : !$OMP END DO
302 : DEALLOCATE (dc)
303 :
304 : ! Expansion coefficients of the cartesian derivatives
305 : ! of the product of two harmonics :
306 : ! d(Y(l1m1) * Y(l2m2))/dx ; d(Y(l1m1) * Y(l2m2))/dy ; d(Y(l1m1) * Y(l2m2))/dz
307 :
308 : !$OMP DO COLLAPSE(3)
309 : DO iso1 = 1, maxs
310 : DO iso2 = 1, maxs
311 : DO iso = 1, max_s_harm
312 : rx = 0.0_dp
313 : ry = 0.0_dp
314 : rz = 0.0_dp
315 :
316 : DO ia = 1, na
317 : rx = rx + wa(ia)*slm(ia, iso)* &
318 : (dslm_dxyz(1, ia, iso1)*slm(ia, iso2) + slm(ia, iso1)*dslm_dxyz(1, ia, iso2))
319 : ry = ry + wa(ia)*slm(ia, iso)* &
320 : (dslm_dxyz(2, ia, iso1)*slm(ia, iso2) + slm(ia, iso1)*dslm_dxyz(2, ia, iso2))
321 : rz = rz + wa(ia)*slm(ia, iso)* &
322 : (dslm_dxyz(3, ia, iso1)*slm(ia, iso2) + slm(ia, iso1)*dslm_dxyz(3, ia, iso2))
323 : END DO
324 :
325 : harmonics%my_CG_dxyz(1, iso1, iso2, iso) = rx
326 : harmonics%my_CG_dxyz(2, iso1, iso2, iso) = ry
327 : harmonics%my_CG_dxyz(3, iso1, iso2, iso) = rz
328 :
329 : END DO
330 : END DO
331 : END DO
332 : !$OMP END DO
333 :
334 : ! Expansion coefficients of the cartesian of the combinations
335 : ! Y(l1m1) * d(Y(l2m2))/dx - d(Y(l1m1))/dx * Y(l2m2)
336 : ! Y(l1m1) * d(Y(l2m2))/dy - d(Y(l1m1))/dy * Y(l2m2)
337 : ! Y(l1m1) * d(Y(l2m2))/dz - d(Y(l1m1))/dz * Y(l2m2)
338 :
339 : !$OMP DO COLLAPSE(3)
340 : DO iso1 = 1, maxs
341 : DO iso2 = 1, maxs
342 : DO iso = 1, max_s_harm
343 : drx = 0.0_dp
344 : dry = 0.0_dp
345 : drz = 0.0_dp
346 :
347 : DO ia = 1, na
348 : drx = drx + wa(ia)*slm(ia, iso)* &
349 : (-dslm_dxyz(1, ia, iso1)*slm(ia, iso2) + &
350 : slm(ia, iso1)*dslm_dxyz(1, ia, iso2))
351 : dry = dry + wa(ia)*slm(ia, iso)* &
352 : (-dslm_dxyz(2, ia, iso1)*slm(ia, iso2) + &
353 : slm(ia, iso1)*dslm_dxyz(2, ia, iso2))
354 : drz = drz + wa(ia)*slm(ia, iso)* &
355 : (-dslm_dxyz(3, ia, iso1)*slm(ia, iso2) + &
356 : slm(ia, iso1)*dslm_dxyz(3, ia, iso2))
357 : END DO
358 :
359 : harmonics%my_CG_dxyz_asym(1, iso1, iso2, iso) = drx
360 : harmonics%my_CG_dxyz_asym(2, iso1, iso2, iso) = dry
361 : harmonics%my_CG_dxyz_asym(3, iso1, iso2, iso) = drz
362 :
363 : END DO ! iso
364 : END DO ! iso2
365 : END DO ! iso1
366 : !$OMP END DO
367 :
368 : ! Calculate the derivatives of the harmonics with respect of the 2 angles
369 : ! the first angle (polar) is acos(lebedev_grid(ll)%r(3))
370 : ! the second angle (azimutal) is atan(lebedev_grid(ll)%r(2)/lebedev_grid(ll)%r(1))
371 : !$OMP DO
372 : DO iso = 1, maxs
373 : l = indso(1, iso)
374 : m = indso(2, iso)
375 : DO ia = 1, na
376 : cin(1) = pol(ia)
377 : cin(2) = azi(ia)
378 : CALL dy_lm(cin, dylm, l, m)
379 : harmonics%dslm(:, ia, iso) = dylm(:)
380 : END DO
381 : END DO
382 : !$OMP END DO
383 :
384 : ! expansion coefficients of product of polar angle derivatives (dslm(1...)) in
385 : ! spherical harmonics (used for tau functionals)
386 : !$OMP END PARALLEL
387 :
388 3030 : CALL timestop(handle)
389 :
390 3030 : END SUBROUTINE create_harmonics_atom
391 :
392 : ! **************************************************************************************************
393 : !> \brief ...
394 : !> \param harmonics ...
395 : !> \param orb_basis ...
396 : !> \param llmax ...
397 : !> \param max_s_harm ...
398 : ! **************************************************************************************************
399 6060 : SUBROUTINE get_maxl_CG(harmonics, orb_basis, llmax, max_s_harm)
400 :
401 : TYPE(harmonics_atom_type), POINTER :: harmonics
402 : TYPE(gto_basis_set_type), POINTER :: orb_basis
403 : INTEGER, INTENT(IN) :: llmax, max_s_harm
404 :
405 : CHARACTER(len=*), PARAMETER :: routineN = 'get_maxl_CG'
406 :
407 : INTEGER :: damax_iso_not0, dmax_iso_not0, handle, &
408 : is1, is2, itmp, max_iso_not0, nset
409 3030 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin
410 :
411 3030 : CALL timeset(routineN, handle)
412 :
413 3030 : CPASSERT(ASSOCIATED(harmonics))
414 :
415 3030 : CALL get_gto_basis_set(gto_basis_set=orb_basis, lmax=lmax, lmin=lmin, nset=nset)
416 :
417 : ! *** Assign indexes for the non null CG coefficients ***
418 3030 : max_iso_not0 = 0
419 3030 : dmax_iso_not0 = 0
420 3030 : damax_iso_not0 = 0
421 11936 : DO is1 = 1, nset
422 49302 : DO is2 = 1, nset
423 : CALL get_none0_cg_list(harmonics%my_CG, &
424 : lmin(is1), lmax(is1), lmin(is2), lmax(is2), &
425 37366 : max_s_harm, llmax, max_iso_not0=itmp)
426 37366 : max_iso_not0 = MAX(max_iso_not0, itmp)
427 : CALL get_none0_cg_list(harmonics%my_CG_dxyz, &
428 : lmin(is1), lmax(is1), lmin(is2), lmax(is2), &
429 37366 : max_s_harm, llmax, max_iso_not0=itmp)
430 37366 : dmax_iso_not0 = MAX(dmax_iso_not0, itmp)
431 : CALL get_none0_cg_list(harmonics%my_CG_dxyz_asym, &
432 : lmin(is1), lmax(is1), lmin(is2), lmax(is2), &
433 37366 : max_s_harm, llmax, max_iso_not0=itmp)
434 46272 : damax_iso_not0 = MAX(damax_iso_not0, itmp)
435 : END DO ! is2
436 : END DO ! is1
437 3030 : harmonics%max_iso_not0 = max_iso_not0
438 3030 : harmonics%dmax_iso_not0 = dmax_iso_not0
439 3030 : harmonics%damax_iso_not0 = damax_iso_not0
440 :
441 3030 : CALL timestop(handle)
442 :
443 3030 : END SUBROUTINE get_maxl_CG
444 :
445 : ! **************************************************************************************************
446 : !> \brief ...
447 : !> \param cgc ...
448 : !> \param lmin1 ...
449 : !> \param lmax1 ...
450 : !> \param lmin2 ...
451 : !> \param lmax2 ...
452 : !> \param max_s_harm ...
453 : !> \param llmax ...
454 : !> \param list ...
455 : !> \param n_list ...
456 : !> \param max_iso_not0 ...
457 : ! **************************************************************************************************
458 1103280 : SUBROUTINE get_none0_cg_list4(cgc, lmin1, lmax1, lmin2, lmax2, max_s_harm, llmax, &
459 1103280 : list, n_list, max_iso_not0)
460 :
461 : REAL(dp), DIMENSION(:, :, :, :), INTENT(IN) :: cgc
462 : INTEGER, INTENT(IN) :: lmin1, lmax1, lmin2, lmax2, max_s_harm, &
463 : llmax
464 : INTEGER, DIMENSION(:, :, :), INTENT(OUT), OPTIONAL :: list
465 : INTEGER, DIMENSION(:), INTENT(OUT), OPTIONAL :: n_list
466 : INTEGER, INTENT(OUT) :: max_iso_not0
467 :
468 : INTEGER :: iso, iso1, iso2, l1, l2, nlist
469 :
470 1103280 : CPASSERT(nsoset(lmax1) <= SIZE(cgc, 2))
471 1103280 : CPASSERT(nsoset(lmax2) <= SIZE(cgc, 3))
472 1103280 : CPASSERT(max_s_harm <= SIZE(cgc, 4))
473 1103280 : IF (PRESENT(n_list) .AND. PRESENT(list)) THEN
474 1028548 : CPASSERT(max_s_harm <= SIZE(list, 3))
475 : END IF
476 1103280 : max_iso_not0 = 0
477 25792408 : IF (PRESENT(n_list) .AND. PRESENT(list)) n_list = 0
478 27140704 : DO iso = 1, max_s_harm
479 26037424 : nlist = 0
480 56888382 : DO l1 = lmin1, lmax1
481 134667362 : DO iso1 = nsoset(l1 - 1) + 1, nsoset(l1)
482 206121051 : DO l2 = lmin2, lmax2
483 97491113 : IF (l1 + l2 > llmax) CYCLE
484 438265168 : DO iso2 = nsoset(l2 - 1) + 1, nsoset(l2)
485 263035525 : IF (ABS(cgc(1, iso1, iso2, iso)) + &
486 : ABS(cgc(2, iso1, iso2, iso)) + &
487 97491113 : ABS(cgc(3, iso1, iso2, iso)) > 1.E-8_dp) THEN
488 33257751 : nlist = nlist + 1
489 33257751 : IF (PRESENT(n_list) .AND. PRESENT(list)) THEN
490 28804337 : list(1, nlist, iso) = iso1
491 28804337 : list(2, nlist, iso) = iso2
492 : END IF
493 33257751 : max_iso_not0 = MAX(max_iso_not0, iso)
494 : END IF
495 : END DO
496 : END DO
497 : END DO
498 : END DO
499 27140704 : IF (PRESENT(n_list) .AND. PRESENT(list)) n_list(iso) = nlist
500 : END DO
501 1103280 : END SUBROUTINE get_none0_cg_list4
502 :
503 : ! **************************************************************************************************
504 : !> \brief ...
505 : !> \param cgc ...
506 : !> \param lmin1 ...
507 : !> \param lmax1 ...
508 : !> \param lmin2 ...
509 : !> \param lmax2 ...
510 : !> \param max_s_harm ...
511 : !> \param llmax ...
512 : !> \param list ...
513 : !> \param n_list ...
514 : !> \param max_iso_not0 ...
515 : ! **************************************************************************************************
516 1912388 : SUBROUTINE get_none0_cg_list3(cgc, lmin1, lmax1, lmin2, lmax2, max_s_harm, llmax, &
517 1912388 : list, n_list, max_iso_not0)
518 :
519 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: cgc
520 : INTEGER, INTENT(IN) :: lmin1, lmax1, lmin2, lmax2, max_s_harm, &
521 : llmax
522 : INTEGER, DIMENSION(:, :, :), INTENT(OUT), OPTIONAL :: list
523 : INTEGER, DIMENSION(:), INTENT(OUT), OPTIONAL :: n_list
524 : INTEGER, INTENT(OUT) :: max_iso_not0
525 :
526 : INTEGER :: iso, iso1, iso2, l1, l2, nlist
527 :
528 1912388 : CPASSERT(nsoset(lmax1) <= SIZE(cgc, 1))
529 1912388 : CPASSERT(nsoset(lmax2) <= SIZE(cgc, 2))
530 1912388 : CPASSERT(max_s_harm <= SIZE(cgc, 3))
531 1912388 : IF (PRESENT(n_list) .AND. PRESENT(list)) THEN
532 1874374 : CPASSERT(max_s_harm <= SIZE(list, 3))
533 : END IF
534 1912388 : max_iso_not0 = 0
535 49040712 : IF (PRESENT(n_list) .AND. PRESENT(list)) n_list = 0
536 48370384 : DO iso = 1, max_s_harm
537 46457996 : nlist = 0
538 101482459 : DO l1 = lmin1, lmax1
539 238244242 : DO iso1 = nsoset(l1 - 1) + 1, nsoset(l1)
540 363383709 : DO l2 = lmin2, lmax2
541 171597463 : IF (l1 + l2 > llmax) CYCLE
542 764710203 : DO iso2 = nsoset(l2 - 1) + 1, nsoset(l2)
543 628010445 : IF (ABS(cgc(iso1, iso2, iso)) > 1.E-8_dp) THEN
544 30510191 : nlist = nlist + 1
545 30510191 : IF (PRESENT(n_list) .AND. PRESENT(list)) THEN
546 29692873 : list(1, nlist, iso) = iso1
547 29692873 : list(2, nlist, iso) = iso2
548 : END IF
549 30510191 : max_iso_not0 = MAX(max_iso_not0, iso)
550 : END IF
551 : END DO
552 : END DO
553 : END DO
554 : END DO
555 48370384 : IF (PRESENT(n_list) .AND. PRESENT(list)) n_list(iso) = nlist
556 : END DO
557 1912388 : END SUBROUTINE get_none0_cg_list3
558 :
559 0 : END MODULE qs_harmonics_atom
|
