(git:98357aa)
Loading...
Searching...
No Matches
qs_rho_atom_methods.F
Go to the documentation of this file.
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!--------------------------------------------------------------------------------------------------!
8
19 USE kinds, ONLY: dp
20 USE kpoint_types, ONLY: get_kpoint_info,&
26 USE mathconstants, ONLY: fourpi,&
27 pi
30 USE orbital_pointers, ONLY: indso,&
31 nsoset
41 USE qs_kind_types, ONLY: get_qs_kind,&
50 USE qs_oce_methods, ONLY: proj_blk
61 USE util, ONLY: get_limit
62 USE whittaker, ONLY: whittaker_c0a,&
64
65!$ USE OMP_LIB, ONLY: omp_get_max_threads, &
66!$ omp_get_thread_num, &
67!$ omp_lock_kind, &
68!$ omp_init_lock, omp_set_lock, &
69!$ omp_unset_lock, omp_destroy_lock
70
71#include "./base/base_uses.f90"
72
73 IMPLICIT NONE
74
75 PRIVATE
76
77! *** Global parameters (only in this module)
78
79 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho_atom_methods'
80
81! *** Public subroutines ***
82
88
89CONTAINS
90
91! **************************************************************************************************
92!> \brief ...
93!> \param para_env ...
94!> \param rho_atom_set ...
95!> \param qs_kind ...
96!> \param atom_list ...
97!> \param natom ...
98!> \param nspins ...
99!> \param tot_rho1_h ...
100!> \param tot_rho1_s ...
101!> \param rho1_h_spin ...
102!> \param rho1_s_spin ...
103!> \param rho1_h_aspin ...
104!> \param rho1_s_aspin ...
105! **************************************************************************************************
106 SUBROUTINE calculate_rho_atom(para_env, rho_atom_set, qs_kind, atom_list, &
107 natom, nspins, tot_rho1_h, tot_rho1_s, &
108 rho1_h_spin, rho1_s_spin, rho1_h_aspin, rho1_s_aspin)
109
110 TYPE(mp_para_env_type), POINTER :: para_env
111 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
112 TYPE(qs_kind_type), INTENT(IN) :: qs_kind
113 INTEGER, DIMENSION(:), INTENT(IN) :: atom_list
114 INTEGER, INTENT(IN) :: natom, nspins
115 REAL(dp), DIMENSION(:), INTENT(INOUT) :: tot_rho1_h, tot_rho1_s
116 REAL(dp), INTENT(INOUT) :: rho1_h_spin, rho1_s_spin, rho1_h_aspin, &
117 rho1_s_aspin
118
119 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_rho_atom'
120
121 INTEGER :: damax_iso_not0_local, handle, i, i1, i2, iat, iatom, icg, ipgf1, ipgf2, ir, &
122 iset1, iset2, iso, iso1, iso1_coeff, iso1_first, iso1_last, iso2, iso2_coeff, iso2_first, &
123 iso2_last, j, l, l_iso, l_sub, l_sum, lmax12, lmax_expansion, lmin12, m1s, m2s, &
124 max_iso_not0, max_iso_not0_local, max_npgf, max_s_harm, maxl, maxso, mepos, n1s, n2s, na, &
125 nr, nset, num_pe, size1, size2
126 INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list, dacg_n_list
127 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list, dacg_list
128 INTEGER, DIMENSION(2) :: bo
129 INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, o2nindex
130 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: done_vgg
131 REAL(dp) :: c1, c2, cpc_h, cpc_s, rfun, rho_h, &
132 rho_s, root_zet12, zet12
133 REAL(dp), ALLOCATABLE, DIMENSION(:) :: erf_zet12, g1, g2, gg0, int1, int2
134 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dgg, gg, gg_lm1, sfun_h, sfun_s
135 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: g_rad, vgg
136 REAL(dp), DIMENSION(:, :), POINTER :: coeff_h, coeff_s, zet
137 REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg
138 REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_cg_dxyz
139 TYPE(grid_atom_type), POINTER :: grid_atom
140 TYPE(gto_basis_set_type), POINTER :: basis_1c
141 TYPE(harmonics_atom_type), POINTER :: harmonics
142
143 CALL timeset(routinen, handle)
144
145 !Note: tau is taken care of separately in qs_vxc_atom.F
146
147 NULLIFY (basis_1c)
148 NULLIFY (harmonics, grid_atom)
149 NULLIFY (lmin, lmax, npgf, zet, my_cg, my_cg_dxyz, coeff_h, coeff_s)
150
151 CALL get_qs_kind(qs_kind, grid_atom=grid_atom, harmonics=harmonics)
152 CALL get_qs_kind(qs_kind, basis_set=basis_1c, basis_type="GAPW_1C")
153
154 CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
155 maxl=maxl, npgf=npgf, nset=nset, zet=zet, &
156 maxso=maxso)
157
158 CALL get_paw_basis_info(basis_1c, o2nindex=o2nindex)
159
160 max_iso_not0 = harmonics%max_iso_not0
161 max_s_harm = harmonics%max_s_harm
162
163 nr = grid_atom%nr
164 max_npgf = maxval(npgf(1:nset))
165 lmax_expansion = indso(1, max_iso_not0)
166 ! Distribute the atoms of this kind
167 num_pe = para_env%num_pe
168 mepos = para_env%mepos
169 bo = get_limit(natom, num_pe, mepos)
170
171 my_cg => harmonics%my_CG
172 my_cg_dxyz => harmonics%my_CG_dxyz
173
174 ALLOCATE (g1(nr), g2(nr), gg0(nr), gg(nr, 0:2*maxl), dgg(nr, 0:2*maxl), gg_lm1(nr, 0:2*maxl))
175 ALLOCATE (erf_zet12(nr), vgg(nr, 0:2*maxl, 0:indso(1, max_iso_not0)))
176 ALLOCATE (done_vgg(0:2*maxl, 0:indso(1, max_iso_not0)))
177 ALLOCATE (int1(nr), int2(nr))
178 ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm), &
179 dacg_list(2, nsoset(maxl)**2, max_s_harm), dacg_n_list(max_s_harm))
180 ALLOCATE (g_rad(nr, max_npgf, nset))
181
182 DO iset1 = 1, nset
183 DO ipgf1 = 1, npgf(iset1)
184 g_rad(1:nr, ipgf1, iset1) = exp(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
185 END DO
186 END DO
187
188 DO iat = bo(1), bo(2)
189 iatom = atom_list(iat)
190 DO i = 1, nspins
191 IF (.NOT. ASSOCIATED(rho_atom_set(iatom)%rho_rad_h(i)%r_coef)) THEN
192 CALL allocate_rho_atom_rad(rho_atom_set, iatom, i, nr, max_iso_not0)
193 ELSE
194 CALL set2zero_rho_atom_rad(rho_atom_set, iatom, i)
195 END IF
196 END DO
197 END DO
198
199 m1s = 0
200 DO iset1 = 1, nset
201 m2s = 0
202 DO iset2 = 1, nset
203
204 CALL get_none0_cg_list(my_cg, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
205 max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
206 cpassert(max_iso_not0_local <= max_iso_not0)
207 CALL get_none0_cg_list(my_cg_dxyz, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
208 max_s_harm, lmax_expansion, dacg_list, dacg_n_list, damax_iso_not0_local)
209 n1s = nsoset(lmax(iset1))
210
211 DO ipgf1 = 1, npgf(iset1)
212 iso1_first = nsoset(lmin(iset1) - 1) + 1 + n1s*(ipgf1 - 1) + m1s
213 iso1_last = nsoset(lmax(iset1)) + n1s*(ipgf1 - 1) + m1s
214 size1 = iso1_last - iso1_first + 1
215 iso1_first = o2nindex(iso1_first)
216 iso1_last = o2nindex(iso1_last)
217 i1 = iso1_last - iso1_first + 1
218 cpassert(size1 == i1)
219 i1 = nsoset(lmin(iset1) - 1) + 1
220
221 g1(1:nr) = g_rad(1:nr, ipgf1, iset1)
222
223 n2s = nsoset(lmax(iset2))
224 DO ipgf2 = 1, npgf(iset2)
225 iso2_first = nsoset(lmin(iset2) - 1) + 1 + n2s*(ipgf2 - 1) + m2s
226 iso2_last = nsoset(lmax(iset2)) + n2s*(ipgf2 - 1) + m2s
227 size2 = iso2_last - iso2_first + 1
228 iso2_first = o2nindex(iso2_first)
229 iso2_last = o2nindex(iso2_last)
230 i2 = iso2_last - iso2_first + 1
231 cpassert(size2 == i2)
232 i2 = nsoset(lmin(iset2) - 1) + 1
233
234 g2(1:nr) = g_rad(1:nr, ipgf2, iset2)
235 lmin12 = lmin(iset1) + lmin(iset2)
236 lmax12 = lmax(iset1) + lmax(iset2)
237
238 zet12 = zet(ipgf1, iset1) + zet(ipgf2, iset2)
239 root_zet12 = sqrt(zet(ipgf1, iset1) + zet(ipgf2, iset2))
240 DO ir = 1, nr
241 erf_zet12(ir) = erf(root_zet12*grid_atom%rad(ir))
242 END DO
243
244 gg = 0.0_dp
245 dgg = 0.0_dp
246 gg_lm1 = 0.0_dp
247 vgg = 0.0_dp
248 done_vgg = .false.
249 ! reduce the number of terms in the expansion local densities
250 IF (lmin12 <= lmax_expansion) THEN
251 IF (lmin12 == 0) THEN
252 gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
253 gg_lm1(1:nr, lmin12) = 0.0_dp
254 gg0(1:nr) = gg(1:nr, lmin12)
255 ELSE
256 gg0(1:nr) = g1(1:nr)*g2(1:nr)
257 gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
258 gg_lm1(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12 - 1)*g1(1:nr)*g2(1:nr)
259 END IF
260
261 ! reduce the number of terms in the expansion local densities
262 IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
263
264 DO l = lmin12 + 1, lmax12
265 gg(1:nr, l) = grid_atom%rad(1:nr)*gg(1:nr, l - 1)
266 gg_lm1(1:nr, l) = gg(1:nr, l - 1)
267 dgg(1:nr, l - 1) = -2.0_dp*(zet(ipgf1, iset1) + zet(ipgf2, iset2))*gg(1:nr, l)
268
269 END DO
270 dgg(1:nr, lmax12) = -2.0_dp*(zet(ipgf1, iset1) + &
271 zet(ipgf2, iset2))*grid_atom%rad(1:nr)*gg(1:nr, lmax12)
272
273 c2 = sqrt(pi*pi*pi/(zet12*zet12*zet12))
274
275 DO iso = 1, max_iso_not0_local
276 l_iso = indso(1, iso)
277 c1 = fourpi/(2._dp*real(l_iso, dp) + 1._dp)
278 DO icg = 1, cg_n_list(iso)
279 iso1 = cg_list(1, icg, iso)
280 iso2 = cg_list(2, icg, iso)
281
282 l = indso(1, iso1) + indso(1, iso2)
283 cpassert(l <= lmax_expansion)
284 IF (done_vgg(l, l_iso)) cycle
285 l_sum = l + l_iso
286 l_sub = l - l_iso
287
288 IF (l_sum == 0) THEN
289 vgg(1:nr, l, l_iso) = erf_zet12(1:nr)*grid_atom%oorad2l(1:nr, 1)*c2
290 ELSE
291 CALL whittaker_c0a(int1, grid_atom%rad, gg0, erf_zet12, zet12, l, l_iso, nr)
292 CALL whittaker_ci(int2, grid_atom%rad, gg0, zet12, l_sub, nr)
293
294 DO ir = 1, nr
295 int2(ir) = grid_atom%rad2l(ir, l_iso)*int2(ir)
296 vgg(ir, l, l_iso) = c1*(int1(ir) + int2(ir))
297 END DO
298 END IF
299 done_vgg(l, l_iso) = .true.
300 END DO
301 END DO
302 END IF ! lmax_expansion
303
304 DO iat = bo(1), bo(2)
305 iatom = atom_list(iat)
306
307 DO i = 1, nspins
308 coeff_h => rho_atom_set(iatom)%cpc_h(i)%r_coef
309 coeff_s => rho_atom_set(iatom)%cpc_s(i)%r_coef
310
311 DO iso = 1, max_iso_not0_local
312 l_iso = indso(1, iso)
313 DO icg = 1, cg_n_list(iso)
314 iso1 = cg_list(1, icg, iso)
315 iso2 = cg_list(2, icg, iso)
316
317 l = indso(1, iso1) + indso(1, iso2)
318 cpassert(l <= lmax_expansion)
319 iso1_coeff = iso1_first + iso1 - i1
320 iso2_coeff = iso2_first + iso2 - i2
321 cpc_h = coeff_h(iso1_coeff, iso2_coeff)*my_cg(iso1, iso2, iso)
322 cpc_s = coeff_s(iso1_coeff, iso2_coeff)*my_cg(iso1, iso2, iso)
323
324 rho_atom_set(iatom)%rho_rad_h(i)%r_coef(1:nr, iso) = &
325 rho_atom_set(iatom)%rho_rad_h(i)%r_coef(1:nr, iso) + &
326 gg(1:nr, l)*cpc_h
327
328 rho_atom_set(iatom)%rho_rad_s(i)%r_coef(1:nr, iso) = &
329 rho_atom_set(iatom)%rho_rad_s(i)%r_coef(1:nr, iso) + &
330 gg(1:nr, l)*cpc_s
331
332 rho_atom_set(iatom)%drho_rad_h(i)%r_coef(1:nr, iso) = &
333 rho_atom_set(iatom)%drho_rad_h(i)%r_coef(1:nr, iso) + &
334 dgg(1:nr, l)*cpc_h
335
336 rho_atom_set(iatom)%drho_rad_s(i)%r_coef(1:nr, iso) = &
337 rho_atom_set(iatom)%drho_rad_s(i)%r_coef(1:nr, iso) + &
338 dgg(1:nr, l)*cpc_s
339
340 rho_atom_set(iatom)%vrho_rad_h(i)%r_coef(1:nr, iso) = &
341 rho_atom_set(iatom)%vrho_rad_h(i)%r_coef(1:nr, iso) + &
342 vgg(1:nr, l, l_iso)*cpc_h
343
344 rho_atom_set(iatom)%vrho_rad_s(i)%r_coef(1:nr, iso) = &
345 rho_atom_set(iatom)%vrho_rad_s(i)%r_coef(1:nr, iso) + &
346 vgg(1:nr, l, l_iso)*cpc_s
347
348 END DO ! icg
349
350 END DO ! iso
351
352 DO iso = 1, max_iso_not0 !damax_iso_not0_local
353 l_iso = indso(1, iso)
354 DO icg = 1, dacg_n_list(iso)
355 iso1 = dacg_list(1, icg, iso)
356 iso2 = dacg_list(2, icg, iso)
357 l = indso(1, iso1) + indso(1, iso2)
358 cpassert(l <= lmax_expansion)
359 iso1_coeff = iso1_first + iso1 - i1
360 iso2_coeff = iso2_first + iso2 - i2
361 cpc_h = coeff_h(iso1_coeff, iso2_coeff)
362 cpc_s = coeff_s(iso1_coeff, iso2_coeff)
363 DO j = 1, 3
364 rho_atom_set(iatom)%rho_rad_h_d(j, i)%r_coef(1:nr, iso) = &
365 rho_atom_set(iatom)%rho_rad_h_d(j, i)%r_coef(1:nr, iso) + &
366 gg_lm1(1:nr, l)*cpc_h*my_cg_dxyz(j, iso1, iso2, iso)
367
368 rho_atom_set(iatom)%rho_rad_s_d(j, i)%r_coef(1:nr, iso) = &
369 rho_atom_set(iatom)%rho_rad_s_d(j, i)%r_coef(1:nr, iso) + &
370 gg_lm1(1:nr, l)*cpc_s*my_cg_dxyz(j, iso1, iso2, iso)
371 END DO
372 END DO ! icg
373
374 END DO ! iso
375
376 END DO ! i
377 END DO ! iat
378
379 END DO ! ipgf2
380 END DO ! ipgf1
381 m2s = m2s + maxso
382 END DO ! iset2
383 m1s = m1s + maxso
384 END DO ! iset1
385
386 DO iat = bo(1), bo(2)
387 iatom = atom_list(iat)
388
389 DO i = 1, nspins
390
391 DO iso = 1, max_iso_not0
392 rho_s = 0.0_dp
393 rho_h = 0.0_dp
394 DO ir = 1, nr
395 rho_h = rho_h + rho_atom_set(iatom)%rho_rad_h(i)%r_coef(ir, iso)*grid_atom%wr(ir)
396 rho_s = rho_s + rho_atom_set(iatom)%rho_rad_s(i)%r_coef(ir, iso)*grid_atom%wr(ir)
397 END DO ! ir
398 tot_rho1_h(i) = tot_rho1_h(i) + rho_h*harmonics%slm_int(iso)
399 tot_rho1_s(i) = tot_rho1_s(i) + rho_s*harmonics%slm_int(iso)
400 END DO ! iso
401
402 END DO ! ispin
403
404 IF (nspins == 2) THEN
405 na = SIZE(harmonics%slm, 1)
406 ALLOCATE (sfun_h(nr, na), sfun_s(nr, na))
407 sfun_h = 0.0_dp
408 sfun_s = 0.0_dp
409 DO iso = 1, max_iso_not0
410 DO ir = 1, nr
411 rfun = grid_atom%wr(ir)*(rho_atom_set(iatom)%rho_rad_h(1)%r_coef(ir, iso) - &
412 rho_atom_set(iatom)%rho_rad_h(2)%r_coef(ir, iso))
413 sfun_h(ir, 1:na) = sfun_h(ir, 1:na) + rfun*harmonics%slm(1:na, iso)*grid_atom%wa(1:na)
414 rfun = grid_atom%wr(ir)*(rho_atom_set(iatom)%rho_rad_s(1)%r_coef(ir, iso) - &
415 rho_atom_set(iatom)%rho_rad_s(2)%r_coef(ir, iso))
416 sfun_s(ir, 1:na) = sfun_s(ir, 1:na) + rfun*harmonics%slm(1:na, iso)*grid_atom%wa(1:na)
417 END DO
418 END DO
419 rho1_h_spin = rho1_h_spin + sum(sfun_h(1:nr, 1:na))
420 rho1_s_spin = rho1_s_spin + sum(sfun_s(1:nr, 1:na))
421 rho1_h_aspin = rho1_h_aspin + sum(abs(sfun_h(1:nr, 1:na)))
422 rho1_s_aspin = rho1_s_aspin + sum(abs(sfun_s(1:nr, 1:na)))
423 DEALLOCATE (sfun_h, sfun_s)
424 END IF
425
426 END DO ! iat
427
428 DEALLOCATE (g1, g2, gg0, gg, gg_lm1, dgg, vgg, done_vgg, erf_zet12, int1, int2, g_rad)
429 DEALLOCATE (cg_list, cg_n_list, dacg_list, dacg_n_list)
430 DEALLOCATE (o2nindex)
431
432 CALL timestop(handle)
433
434 END SUBROUTINE calculate_rho_atom
435
436! **************************************************************************************************
437!> \brief Replicate the radial hard/soft density data needed to evaluate one-center tails on
438!> rank-local target grids. The compact one-center density matrices are already global;
439!> this routine performs one packed reduction for the derived radial fields of a kind.
440!> \param para_env ...
441!> \param rho_atom_set ...
442!> \param qs_kind ...
443!> \param atom_list ...
444!> \param natom ...
445!> \param nspins ...
446! **************************************************************************************************
447 SUBROUTINE replicate_rho_atom_radial(para_env, rho_atom_set, qs_kind, atom_list, natom, nspins)
448 TYPE(mp_para_env_type), POINTER :: para_env
449 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
450 TYPE(qs_kind_type), INTENT(IN) :: qs_kind
451 INTEGER, DIMENSION(:), INTENT(IN) :: atom_list
452 INTEGER, INTENT(IN) :: natom, nspins
453
454 CHARACTER(len=*), PARAMETER :: routinen = 'replicate_rho_atom_radial'
455
456 INTEGER :: block_size, bo(2), cursor, handle, iat, &
457 iatom, ispin, j, max_iso_not0, ncoeff, &
458 nr
459 REAL(dp), ALLOCATABLE, DIMENSION(:) :: buffer
460 TYPE(grid_atom_type), POINTER :: grid_atom
461 TYPE(harmonics_atom_type), POINTER :: harmonics
462
463 CALL timeset(routinen, handle)
464
465 NULLIFY (grid_atom, harmonics)
466 CALL get_qs_kind(qs_kind, grid_atom=grid_atom, harmonics=harmonics)
467 cpassert(ASSOCIATED(grid_atom))
468 cpassert(ASSOCIATED(harmonics))
469 nr = grid_atom%nr
470 max_iso_not0 = harmonics%max_iso_not0
471 ncoeff = nr*max_iso_not0
472 block_size = 10*nspins*ncoeff
473 ALLOCATE (buffer(block_size*natom))
474 buffer = 0.0_dp
475
476 bo = get_limit(natom, para_env%num_pe, para_env%mepos)
477 DO iat = bo(1), bo(2)
478 iatom = atom_list(iat)
479 cursor = (iat - 1)*block_size
480 DO ispin = 1, nspins
481 buffer(cursor + 1:cursor + ncoeff) = &
482 reshape(rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef, [ncoeff])
483 cursor = cursor + ncoeff
484 buffer(cursor + 1:cursor + ncoeff) = &
485 reshape(rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef, [ncoeff])
486 cursor = cursor + ncoeff
487 buffer(cursor + 1:cursor + ncoeff) = &
488 reshape(rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef, [ncoeff])
489 cursor = cursor + ncoeff
490 buffer(cursor + 1:cursor + ncoeff) = &
491 reshape(rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef, [ncoeff])
492 cursor = cursor + ncoeff
493 DO j = 1, 3
494 buffer(cursor + 1:cursor + ncoeff) = &
495 reshape(rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef, [ncoeff])
496 cursor = cursor + ncoeff
497 buffer(cursor + 1:cursor + ncoeff) = &
498 reshape(rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef, [ncoeff])
499 cursor = cursor + ncoeff
500 END DO
501 END DO
502 cpassert(cursor == iat*block_size)
503 END DO
504 CALL para_env%sum(buffer)
505
506 DO iat = 1, natom
507 iatom = atom_list(iat)
508 cursor = (iat - 1)*block_size
509 DO ispin = 1, nspins
510 IF (.NOT. ASSOCIATED(rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef)) THEN
511 CALL allocate_rho_atom_rad(rho_atom_set, iatom, ispin, nr, max_iso_not0)
512 END IF
513 rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef = &
514 reshape(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
515 cursor = cursor + ncoeff
516 rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef = &
517 reshape(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
518 cursor = cursor + ncoeff
519 rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef = &
520 reshape(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
521 cursor = cursor + ncoeff
522 rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef = &
523 reshape(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
524 cursor = cursor + ncoeff
525 DO j = 1, 3
526 rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef = &
527 reshape(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
528 cursor = cursor + ncoeff
529 rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef = &
530 reshape(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
531 cursor = cursor + ncoeff
532 END DO
533 END DO
534 cpassert(cursor == iat*block_size)
535 END DO
536 DEALLOCATE (buffer)
537
538 CALL timestop(handle)
539
540 END SUBROUTINE replicate_rho_atom_radial
541
542! **************************************************************************************************
543!> \brief ...
544!> \param qs_env QuickStep environment
545!> (accessed components: atomic_kind_set, dft_control%nimages,
546!> dft_control%nspins, kpoints%cell_to_index)
547!> \param rho_ao density matrix in atomic basis set
548!> \param rho_atom_set ...
549!> \param qs_kind_set list of QuickStep kinds
550!> \param oce one-centre expansion coefficients
551!> \param sab neighbour pair list
552!> \param para_env parallel environment
553!> \par History
554!> Add OpenMP [Apr 2016, EPCC]
555!> Use automatic arrays [Sep 2016, M Tucker]
556!> Allow for external non-default kind_set, oce and sab [Dec 2019, A Bussy]
557!> \note Consider to declare 'rho_ao' dummy argument as a pointer to the two-dimensional
558!> (1:nspins, 1:nimages) set of matrices.
559! **************************************************************************************************
560 SUBROUTINE calculate_rho_atom_coeff(qs_env, rho_ao, rho_atom_set, qs_kind_set, oce, sab, para_env)
561
562 TYPE(qs_environment_type), POINTER :: qs_env
563 TYPE(dbcsr_p_type), DIMENSION(*) :: rho_ao
564 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
565 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
566 TYPE(oce_matrix_type), POINTER :: oce
567 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
568 POINTER :: sab
569 TYPE(mp_para_env_type), POINTER :: para_env
570
571 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_rho_atom_coeff'
572
573 INTEGER :: bo(2), handle, i, iac, iatom, ibc, icol, ikind, img, irow, ispin, jatom, jkind, &
574 kac, katom, kbc, kkind, len_cpc, len_pc1, max_gau, max_nsgf, mepos, n_cont_a, n_cont_b, &
575 nat_kind, natom, nimages, nkind, nsoctot, nspins, num_pe
576 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, nsatbas_kind
577 INTEGER, DIMENSION(3) :: cell_b
578 INTEGER, DIMENSION(:), POINTER :: a_list, list_a, list_b
579 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
580 LOGICAL :: dista, distab, distb, found, paw_atom
581 LOGICAL, ALLOCATABLE, DIMENSION(:) :: has_intac, paw_kind
582 REAL(dp), ALLOCATABLE, DIMENSION(:) :: proj_work1, proj_work2
583 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: p_matrix
584 REAL(kind=dp) :: eps_cpc, factor, pmax
585 REAL(kind=dp), DIMENSION(3) :: rab
586 REAL(kind=dp), DIMENSION(:, :), POINTER :: c_coeff_hh_a, c_coeff_hh_b, &
587 c_coeff_ss_a, c_coeff_ss_b, r_coef_h, &
588 r_coef_s
589 TYPE(alist_type), POINTER :: alist_ac, alist_bc
590 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
591 TYPE(dft_control_type), POINTER :: dft_control
592 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
593 TYPE(gto_basis_set_type), POINTER :: basis_1c, basis_set_a, basis_set_b
594 TYPE(kpoint_type), POINTER :: kpoints
596 DIMENSION(:), POINTER :: nl_iterator
597 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: p_block_spin
598
599!$ INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: locks
600!$ INTEGER :: lock, number_of_locks
601
602 CALL timeset(routinen, handle)
603
604 CALL get_qs_env(qs_env=qs_env, &
605 dft_control=dft_control, &
606 atomic_kind_set=atomic_kind_set)
607
608 eps_cpc = dft_control%qs_control%gapw_control%eps_cpc
609
610 cpassert(ASSOCIATED(qs_kind_set))
611 cpassert(ASSOCIATED(rho_atom_set))
612 cpassert(ASSOCIATED(oce))
613 cpassert(ASSOCIATED(sab))
614
615 nspins = dft_control%nspins
616 nimages = dft_control%nimages
617
618 NULLIFY (cell_to_index)
619 IF (nimages > 1) THEN
620 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
621 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
622 END IF
623
624 CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
625 CALL get_qs_kind_set(qs_kind_set, maxsgf=max_nsgf, maxgtops=max_gau, basis_type='GAPW_1C')
626
627 nkind = SIZE(atomic_kind_set)
628 ! Initialize to 0 the CPC coefficients and the local density arrays
629 DO ikind = 1, nkind
630 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=a_list, natom=nat_kind)
631 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
632
633 IF (.NOT. paw_atom) cycle
634 DO i = 1, nat_kind
635 iatom = a_list(i)
636 DO ispin = 1, nspins
637 rho_atom_set(iatom)%cpc_h(ispin)%r_coef = 0.0_dp
638 rho_atom_set(iatom)%cpc_s(ispin)%r_coef = 0.0_dp
639 END DO ! ispin
640 END DO ! i
641
642 num_pe = para_env%num_pe
643 mepos = para_env%mepos
644 bo = get_limit(nat_kind, num_pe, mepos)
645 DO i = bo(1), bo(2)
646 iatom = a_list(i)
647 DO ispin = 1, nspins
648 rho_atom_set(iatom)%ga_Vlocal_gb_h(ispin)%r_coef = 0.0_dp
649 rho_atom_set(iatom)%ga_Vlocal_gb_s(ispin)%r_coef = 0.0_dp
650 END DO ! ispin
651 END DO ! i
652 END DO ! ikind
653
654 ALLOCATE (basis_set_list(nkind))
655 ALLOCATE (paw_kind(nkind), nsatbas_kind(nkind), has_intac(nkind*nkind))
656 paw_kind(:) = .false.
657 nsatbas_kind(:) = 0
658 has_intac(:) = .false.
659 DO ikind = 1, nkind
660 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a)
661 IF (ASSOCIATED(basis_set_a)) THEN
662 basis_set_list(ikind)%gto_basis_set => basis_set_a
663 ELSE
664 NULLIFY (basis_set_list(ikind)%gto_basis_set)
665 END IF
666 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C", &
667 paw_atom=paw_kind(ikind))
668 IF (paw_kind(ikind)) CALL get_paw_basis_info(basis_1c, nsatbas=nsatbas_kind(ikind))
669 END DO
670 DO ikind = 1, nkind*nkind
671 has_intac(ikind) = ASSOCIATED(oce%intac(ikind)%alist)
672 END DO
673
674 len_pc1 = max_nsgf*max_gau
675 len_cpc = max_gau*max_gau
676
677 num_pe = 1
678!$ num_pe = omp_get_max_threads()
679 CALL neighbor_list_iterator_create(nl_iterator, sab, nthread=num_pe)
680
681!$OMP PARALLEL DEFAULT( NONE ) &
682!$OMP SHARED( max_nsgf, max_gau &
683!$OMP , len_PC1, len_CPC &
684!$OMP , nl_iterator, basis_set_list &
685!$OMP , nimages, cell_to_index &
686!$OMP , nspins, rho_ao &
687!$OMP , nkind, qs_kind_set &
688!$OMP , oce, eps_cpc &
689!$OMP , rho_atom_set &
690!$OMP , natom, locks, number_of_locks &
691!$OMP , paw_kind, nsatbas_kind, has_intac &
692!$OMP ) &
693!$OMP PRIVATE( p_block_spin, ispin &
694!$OMP , p_matrix, proj_work1, proj_work2 &
695!$OMP , mepos &
696!$OMP , ikind, jkind, iatom, jatom &
697!$OMP , cell_b, rab &
698!$OMP , basis_set_a, basis_set_b &
699!$OMP , pmax, irow, icol, img &
700!$OMP , found &
701!$OMP , kkind &
702!$OMP , nsoctot, katom &
703!$OMP , iac , alist_ac, kac, n_cont_a, list_a &
704!$OMP , ibc , alist_bc, kbc, n_cont_b, list_b &
705!$OMP , C_coeff_hh_a, C_coeff_ss_a, dista &
706!$OMP , C_coeff_hh_b, C_coeff_ss_b, distb &
707!$OMP , distab &
708!$OMP , factor, r_coef_h, r_coef_s &
709!$OMP )
710
711 ALLOCATE (p_block_spin(nspins))
712 ALLOCATE (p_matrix(max_nsgf, max_nsgf))
713 ALLOCATE (proj_work1(len_pc1), proj_work2(len_cpc))
714
715!$OMP SINGLE
716!$ number_of_locks = nspins*natom
717!$ ALLOCATE (locks(number_of_locks))
718!$OMP END SINGLE
719
720!$OMP DO
721!$ DO lock = 1, number_of_locks
722!$ call omp_init_lock(locks(lock))
723!$ END DO
724!$OMP END DO
725
726 mepos = 0
727!$ mepos = omp_get_thread_num()
728 DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
729
730 CALL get_iterator_info(nl_iterator, mepos=mepos, &
731 ikind=ikind, jkind=jkind, &
732 iatom=iatom, jatom=jatom, &
733 cell=cell_b, r=rab)
734
735 basis_set_a => basis_set_list(ikind)%gto_basis_set
736 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
737 basis_set_b => basis_set_list(jkind)%gto_basis_set
738 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
739
740 pmax = 0._dp
741 IF (iatom <= jatom) THEN
742 irow = iatom
743 icol = jatom
744 ELSE
745 irow = jatom
746 icol = iatom
747 END IF
748
749 IF (nimages > 1) THEN
750 img = cell_to_index(cell_b(1), cell_b(2), cell_b(3))
751 cpassert(img > 0)
752 ELSE
753 img = 1
754 END IF
755
756 DO ispin = 1, nspins
757 CALL dbcsr_get_block_p(matrix=rho_ao(nspins*(img - 1) + ispin)%matrix, &
758 row=irow, col=icol, block=p_block_spin(ispin)%r_coef, &
759 found=found)
760 pmax = pmax + maxval(abs(p_block_spin(ispin)%r_coef))
761 END DO
762
763 DO kkind = 1, nkind
764 IF (.NOT. paw_kind(kkind)) cycle
765
766 nsoctot = nsatbas_kind(kkind)
767
768 iac = ikind + nkind*(kkind - 1)
769 ibc = jkind + nkind*(kkind - 1)
770 IF (.NOT. has_intac(iac)) cycle
771 IF (.NOT. has_intac(ibc)) cycle
772
773 CALL get_alist(oce%intac(iac), alist_ac, iatom)
774 CALL get_alist(oce%intac(ibc), alist_bc, jatom)
775 IF (.NOT. ASSOCIATED(alist_ac)) cycle
776 IF (.NOT. ASSOCIATED(alist_bc)) cycle
777
778 DO kac = 1, alist_ac%nclist
779 DO kbc = 1, alist_bc%nclist
780 IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
781 IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
782 IF (pmax*alist_bc%clist(kbc)%maxac*alist_ac%clist(kac)%maxac < eps_cpc) cycle
783
784 n_cont_a = alist_ac%clist(kac)%nsgf_cnt
785 n_cont_b = alist_bc%clist(kbc)%nsgf_cnt
786 IF (n_cont_a == 0 .OR. n_cont_b == 0) cycle
787
788 list_a => alist_ac%clist(kac)%sgf_list
789 list_b => alist_bc%clist(kbc)%sgf_list
790
791 katom = alist_ac%clist(kac)%catom
792
793 IF (iatom == katom .AND. all(alist_ac%clist(kac)%cell == 0)) THEN
794 c_coeff_hh_a => alist_ac%clist(kac)%achint(:, :, 1)
795 c_coeff_ss_a => alist_ac%clist(kac)%acint(:, :, 1)
796 dista = .false.
797 ELSE
798 c_coeff_hh_a => alist_ac%clist(kac)%acint(:, :, 1)
799 c_coeff_ss_a => alist_ac%clist(kac)%acint(:, :, 1)
800 dista = .true.
801 END IF
802 IF (jatom == katom .AND. all(alist_bc%clist(kbc)%cell == 0)) THEN
803 c_coeff_hh_b => alist_bc%clist(kbc)%achint(:, :, 1)
804 c_coeff_ss_b => alist_bc%clist(kbc)%acint(:, :, 1)
805 distb = .false.
806 ELSE
807 c_coeff_hh_b => alist_bc%clist(kbc)%acint(:, :, 1)
808 c_coeff_ss_b => alist_bc%clist(kbc)%acint(:, :, 1)
809 distb = .true.
810 END IF
811
812 distab = dista .AND. distb
813
814 DO ispin = 1, nspins
815
816 IF (iatom <= jatom) THEN
817 CALL alist_pre_align_blk(p_block_spin(ispin)%r_coef, &
818 SIZE(p_block_spin(ispin)%r_coef, 1), p_matrix, SIZE(p_matrix, 1), &
819 list_a, n_cont_a, list_b, n_cont_b)
820 ELSE
821 CALL alist_pre_align_blk(p_block_spin(ispin)%r_coef, &
822 SIZE(p_block_spin(ispin)%r_coef, 1), p_matrix, SIZE(p_matrix, 1), &
823 list_b, n_cont_b, list_a, n_cont_a)
824 END IF
825
826 factor = 1.0_dp
827 IF (iatom == jatom) factor = 0.5_dp
828
829 r_coef_h => rho_atom_set(katom)%cpc_h(ispin)%r_coef
830 r_coef_s => rho_atom_set(katom)%cpc_s(ispin)%r_coef
831
832!$ CALL omp_set_lock(locks((katom - 1)*nspins + ispin))
833 IF (iatom <= jatom) THEN
834 CALL proj_blk(c_coeff_hh_a, c_coeff_ss_a, n_cont_a, &
835 c_coeff_hh_b, c_coeff_ss_b, n_cont_b, &
836 p_matrix, max_nsgf, r_coef_h, r_coef_s, nsoctot, &
837 len_pc1, len_cpc, factor, distab, proj_work1, proj_work2)
838 ELSE
839 CALL proj_blk(c_coeff_hh_b, c_coeff_ss_b, n_cont_b, &
840 c_coeff_hh_a, c_coeff_ss_a, n_cont_a, &
841 p_matrix, max_nsgf, r_coef_h, r_coef_s, nsoctot, &
842 len_pc1, len_cpc, factor, distab, proj_work1, proj_work2)
843 END IF
844!$ CALL omp_unset_lock(locks((katom - 1)*nspins + ispin))
845
846 END DO
847 EXIT !search loop over jatom-katom list
848 END IF
849 END DO
850 END DO
851 END DO
852 END DO
853 ! Wait for all threads to finish the loop before locks can be freed
854!$OMP BARRIER
855
856!$OMP DO
857!$ DO lock = 1, number_of_locks
858!$ call omp_destroy_lock(locks(lock))
859!$ END DO
860!$OMP END DO
861!$OMP SINGLE
862!$ DEALLOCATE (locks)
863!$OMP END SINGLE NOWAIT
864
865 DEALLOCATE (p_block_spin, p_matrix, proj_work1, proj_work2)
866!$OMP END PARALLEL
867
868 CALL neighbor_list_iterator_release(nl_iterator)
869
870 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
871
872 DO iatom = 1, natom
873 ikind = kind_of(iatom)
874
875 DO ispin = 1, nspins
876 IF (ASSOCIATED(rho_atom_set(iatom)%cpc_h(ispin)%r_coef)) THEN
877 CALL para_env%sum(rho_atom_set(iatom)%cpc_h(ispin)%r_coef)
878 CALL para_env%sum(rho_atom_set(iatom)%cpc_s(ispin)%r_coef)
879 r_coef_h => rho_atom_set(iatom)%cpc_h(ispin)%r_coef
880 r_coef_s => rho_atom_set(iatom)%cpc_s(ispin)%r_coef
881 r_coef_h(:, :) = r_coef_h(:, :) + transpose(r_coef_h(:, :))
882 r_coef_s(:, :) = r_coef_s(:, :) + transpose(r_coef_s(:, :))
883 END IF
884 END DO
885
886 END DO
887
888 DEALLOCATE (kind_of, basis_set_list, paw_kind, nsatbas_kind, has_intac)
889
890 CALL timestop(handle)
891
892 END SUBROUTINE calculate_rho_atom_coeff
893
894! **************************************************************************************************
895!> \brief ...
896!> \param rho_atom_set the type to initialize
897!> \param atomic_kind_set list of atomic kinds
898!> \param qs_kind_set the kind set from which to take quantum numbers and basis info
899!> \param dft_control DFT control type
900!> \param para_env parallel environment
901!> \par History:
902!> - Generalised by providing the rho_atom_set and the qs_kind_set 12.2019 (A.Bussy)
903! **************************************************************************************************
904 SUBROUTINE init_rho_atom(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
905
906 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
907 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
908 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
909 TYPE(dft_control_type), POINTER :: dft_control
910 TYPE(mp_para_env_type), POINTER :: para_env
911
912 CHARACTER(len=*), PARAMETER :: routinen = 'init_rho_atom'
913
914 INTEGER :: handle, ikind, il, iso, iso1, iso2, l1, l1l2, l2, la, lc1, lc2, lcleb, ll, llmax, &
915 lmax_sphere, lp, m1, m2, max_s_harm, max_s_set, maxl, maxlgto, maxs, mm, mp, na, nat, &
916 natom, nr, nspins, quadrature
917 INTEGER, DIMENSION(:), POINTER :: atom_list
918 LOGICAL :: paw_atom
919 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rga
920 REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg
921 TYPE(gapw_control_type), POINTER :: gapw_control
922 TYPE(grid_atom_type), POINTER :: grid_atom
923 TYPE(gto_basis_set_type), POINTER :: basis_1c_set
924 TYPE(harmonics_atom_type), POINTER :: harmonics
925
926 CALL timeset(routinen, handle)
927
928 NULLIFY (basis_1c_set)
929 NULLIFY (my_cg, grid_atom, harmonics, atom_list)
930
931 cpassert(ASSOCIATED(atomic_kind_set))
932 cpassert(ASSOCIATED(dft_control))
933 cpassert(ASSOCIATED(para_env))
934 cpassert(ASSOCIATED(qs_kind_set))
935
936 CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
937
938 CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="GAPW_1C")
939
940 nspins = dft_control%nspins
941 gapw_control => dft_control%qs_control%gapw_control
942
943 lmax_sphere = gapw_control%lmax_sphere
944
945 llmax = min(lmax_sphere, 2*maxlgto)
946 max_s_harm = nsoset(llmax)
947 max_s_set = nsoset(maxlgto)
948
949 lcleb = max(llmax, 2*maxlgto, 1)
950
951! *** allocate calculate the CG coefficients up to the maxl ***
952 CALL clebsch_gordon_init(lcleb)
953 CALL reallocate(my_cg, 1, max_s_set, 1, max_s_set, 1, max_s_harm)
954
955 ALLOCATE (rga(lcleb, 2))
956 DO lc1 = 0, maxlgto
957 DO iso1 = nsoset(lc1 - 1) + 1, nsoset(lc1)
958 l1 = indso(1, iso1)
959 m1 = indso(2, iso1)
960 DO lc2 = 0, maxlgto
961 DO iso2 = nsoset(lc2 - 1) + 1, nsoset(lc2)
962 l2 = indso(1, iso2)
963 m2 = indso(2, iso2)
964 CALL clebsch_gordon(l1, m1, l2, m2, rga)
965 IF (l1 + l2 > llmax) THEN
966 l1l2 = llmax
967 ELSE
968 l1l2 = l1 + l2
969 END IF
970 mp = m1 + m2
971 mm = m1 - m2
972 IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN
973 mp = -abs(mp)
974 mm = -abs(mm)
975 ELSE
976 mp = abs(mp)
977 mm = abs(mm)
978 END IF
979 DO lp = mod(l1 + l2, 2), l1l2, 2
980 il = lp/2 + 1
981 IF (abs(mp) <= lp) THEN
982 IF (mp >= 0) THEN
983 iso = nsoset(lp - 1) + lp + 1 + mp
984 ELSE
985 iso = nsoset(lp - 1) + lp + 1 - abs(mp)
986 END IF
987 my_cg(iso1, iso2, iso) = rga(il, 1)
988 END IF
989 IF (mp /= mm .AND. abs(mm) <= lp) THEN
990 IF (mm >= 0) THEN
991 iso = nsoset(lp - 1) + lp + 1 + mm
992 ELSE
993 iso = nsoset(lp - 1) + lp + 1 - abs(mm)
994 END IF
995 my_cg(iso1, iso2, iso) = rga(il, 2)
996 END IF
997 END DO
998 END DO ! iso2
999 END DO ! lc2
1000 END DO ! iso1
1001 END DO ! lc1
1002 DEALLOCATE (rga)
1004
1005! *** initialize the Lebedev grids ***
1006 CALL init_lebedev_grids()
1007 quadrature = gapw_control%quadrature
1008
1009 DO ikind = 1, SIZE(atomic_kind_set)
1010 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
1011 CALL get_qs_kind(qs_kind_set(ikind), &
1012 paw_atom=paw_atom, &
1013 grid_atom=grid_atom, &
1014 harmonics=harmonics, &
1015 ngrid_rad=nr, ngrid_ang=na)
1016
1017! *** determine the Lebedev grid for this kind ***
1018
1020 na = lebedev_grid(ll)%n
1021 la = lebedev_grid(ll)%l
1022 grid_atom%ng_sphere = na
1023 grid_atom%nr = nr
1024
1025 IF (llmax > la) THEN
1026 WRITE (*, '(/,72("*"))')
1027 WRITE (*, '(T2,A,T66,I4)') &
1028 "WARNING: the lebedev grid is built for angular momentum l up to ", la, &
1029 " the max l of spherical harmonics is larger, l_max = ", llmax, &
1030 " good integration is guaranteed only for l <= ", la
1031 WRITE (*, '(72("*"),/)')
1032 END IF
1033
1034! *** calculate the radial grid ***
1035 CALL create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
1036
1037! *** calculate the spherical harmonics on the grid ***
1038
1039 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_1c_set, basis_type="GAPW_1C")
1040 CALL get_gto_basis_set(gto_basis_set=basis_1c_set, maxl=maxl)
1041 maxs = nsoset(maxl)
1042 CALL create_harmonics_atom(harmonics, &
1043 my_cg, na, llmax, maxs, max_s_harm, ll, grid_atom%wa, &
1044 grid_atom%azi, grid_atom%pol)
1045 CALL get_maxl_cg(harmonics, basis_1c_set, llmax, max_s_harm)
1046
1047 END DO
1048
1050 DEALLOCATE (my_cg)
1051
1052 CALL allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
1053
1054 CALL timestop(handle)
1055
1056 END SUBROUTINE init_rho_atom
1057
1058! **************************************************************************************************
1059!> \brief ...
1060!> \param rho_atom_set ...
1061!> \param atomic_kind_set list of atomic kinds
1062!> \param qs_kind_set the kind set from which to take quantum numbers and basis info
1063!> \param dft_control DFT control type
1064!> \param para_env parallel environment
1065! **************************************************************************************************
1066 SUBROUTINE allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
1067
1068 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
1069 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1070 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1071 TYPE(dft_control_type), POINTER :: dft_control
1072 TYPE(mp_para_env_type), POINTER :: para_env
1073
1074 CHARACTER(len=*), PARAMETER :: routinen = 'allocate_rho_atom_internals'
1075
1076 INTEGER :: bo(2), handle, iat, iatom, ikind, ispin, &
1077 max_iso_not0, maxso, mepos, nat, &
1078 natom, nsatbas, nset, nsotot, nspins, &
1079 num_pe
1080 INTEGER, DIMENSION(:), POINTER :: atom_list
1081 LOGICAL :: paw_atom
1082 TYPE(gto_basis_set_type), POINTER :: basis_1c
1083 TYPE(harmonics_atom_type), POINTER :: harmonics
1084
1085 CALL timeset(routinen, handle)
1086
1087 cpassert(ASSOCIATED(atomic_kind_set))
1088 cpassert(ASSOCIATED(dft_control))
1089 cpassert(ASSOCIATED(para_env))
1090 cpassert(ASSOCIATED(qs_kind_set))
1091
1092 CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
1093
1094 nspins = dft_control%nspins
1095
1096 IF (ASSOCIATED(rho_atom_set)) THEN
1097 CALL deallocate_rho_atom_set(rho_atom_set)
1098 END IF
1099 ALLOCATE (rho_atom_set(natom))
1100
1101 DO ikind = 1, SIZE(atomic_kind_set)
1102
1103 NULLIFY (atom_list, harmonics)
1104 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
1105 CALL get_qs_kind(qs_kind_set(ikind), &
1106 paw_atom=paw_atom, &
1107 harmonics=harmonics)
1108
1109 IF (paw_atom) THEN
1110 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
1111 CALL get_gto_basis_set(gto_basis_set=basis_1c, nset=nset, maxso=maxso)
1112 nsotot = nset*maxso
1113 CALL get_paw_basis_info(basis_1c, nsatbas=nsatbas)
1114 END IF
1115
1116 max_iso_not0 = harmonics%max_iso_not0
1117 DO iat = 1, nat
1118 iatom = atom_list(iat)
1119 ! *** allocate the radial density for each LM,for each atom ***
1120
1121 ALLOCATE (rho_atom_set(iatom)%rho_rad_h(nspins))
1122 ALLOCATE (rho_atom_set(iatom)%rho_rad_s(nspins))
1123 ALLOCATE (rho_atom_set(iatom)%vrho_rad_h(nspins))
1124 ALLOCATE (rho_atom_set(iatom)%vrho_rad_s(nspins))
1125
1126 ALLOCATE (rho_atom_set(iatom)%cpc_h(nspins), &
1127 rho_atom_set(iatom)%cpc_s(nspins), &
1128 rho_atom_set(iatom)%drho_rad_h(nspins), &
1129 rho_atom_set(iatom)%drho_rad_s(nspins), &
1130 rho_atom_set(iatom)%rho_rad_h_d(3, nspins), &
1131 rho_atom_set(iatom)%rho_rad_s_d(3, nspins))
1132 ALLOCATE (rho_atom_set(iatom)%int_scr_h(nspins), &
1133 rho_atom_set(iatom)%int_scr_s(nspins))
1134
1135 IF (paw_atom) THEN
1136 DO ispin = 1, nspins
1137 ALLOCATE (rho_atom_set(iatom)%cpc_h(ispin)%r_coef(1:nsatbas, 1:nsatbas), &
1138 rho_atom_set(iatom)%cpc_s(ispin)%r_coef(1:nsatbas, 1:nsatbas))
1139 ALLOCATE (rho_atom_set(iatom)%int_scr_h(ispin)%r_coef(1:nsatbas, 1:nsatbas), &
1140 rho_atom_set(iatom)%int_scr_s(ispin)%r_coef(1:nsatbas, 1:nsatbas))
1141
1142 rho_atom_set(iatom)%cpc_h(ispin)%r_coef = 0.0_dp
1143 rho_atom_set(iatom)%cpc_s(ispin)%r_coef = 0.0_dp
1144 END DO
1145 END IF
1146
1147 END DO ! iat
1148
1149 num_pe = para_env%num_pe
1150 mepos = para_env%mepos
1151 bo = get_limit(nat, num_pe, mepos)
1152 DO iat = bo(1), bo(2)
1153 iatom = atom_list(iat)
1154 ALLOCATE (rho_atom_set(iatom)%ga_Vlocal_gb_h(nspins), &
1155 rho_atom_set(iatom)%ga_Vlocal_gb_s(nspins))
1156 IF (paw_atom) THEN
1157 DO ispin = 1, nspins
1158 CALL reallocate(rho_atom_set(iatom)%ga_Vlocal_gb_h(ispin)%r_coef, &
1159 1, nsotot, 1, nsotot)
1160 CALL reallocate(rho_atom_set(iatom)%ga_Vlocal_gb_s(ispin)%r_coef, &
1161 1, nsotot, 1, nsotot)
1162
1163 rho_atom_set(iatom)%ga_Vlocal_gb_h(ispin)%r_coef = 0.0_dp
1164 rho_atom_set(iatom)%ga_Vlocal_gb_s(ispin)%r_coef = 0.0_dp
1165 END DO
1166 END IF
1167
1168 END DO ! iat
1169
1170 END DO
1171
1172 CALL timestop(handle)
1173
1174 END SUBROUTINE allocate_rho_atom_internals
1175
1176! **************************************************************************************************
1177!> \brief ...
1178!> \param rho_atom_set ...
1179!> \param iatom ...
1180!> \param ispin ...
1181!> \param nr ...
1182!> \param max_iso_not0 ...
1183! **************************************************************************************************
1184 SUBROUTINE allocate_rho_atom_rad(rho_atom_set, iatom, ispin, nr, max_iso_not0)
1185
1186 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
1187 INTEGER, INTENT(IN) :: iatom, ispin, nr, max_iso_not0
1188
1189 CHARACTER(len=*), PARAMETER :: routinen = 'allocate_rho_atom_rad'
1190
1191 INTEGER :: handle, j
1192
1193 CALL timeset(routinen, handle)
1194
1195 ALLOCATE (rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef(1:nr, 1:max_iso_not0), &
1196 rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef(1:nr, 1:max_iso_not0), &
1197 rho_atom_set(iatom)%vrho_rad_h(ispin)%r_coef(1:nr, 1:max_iso_not0), &
1198 rho_atom_set(iatom)%vrho_rad_s(ispin)%r_coef(1:nr, 1:max_iso_not0))
1199
1200 rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef = 0.0_dp
1201 rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef = 0.0_dp
1202 rho_atom_set(iatom)%vrho_rad_h(ispin)%r_coef = 0.0_dp
1203 rho_atom_set(iatom)%vrho_rad_s(ispin)%r_coef = 0.0_dp
1204
1205 ALLOCATE (rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef(nr, max_iso_not0), &
1206 rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef(nr, max_iso_not0))
1207 rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef = 0.0_dp
1208 rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef = 0.0_dp
1209
1210 DO j = 1, 3
1211 ALLOCATE (rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef(nr, max_iso_not0), &
1212 rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef(nr, max_iso_not0))
1213 rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef = 0.0_dp
1214 rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef = 0.0_dp
1215 END DO
1216
1217 CALL timestop(handle)
1218
1219 END SUBROUTINE allocate_rho_atom_rad
1220
1221! **************************************************************************************************
1222!> \brief ...
1223!> \param rho_atom_set ...
1224!> \param iatom ...
1225!> \param ispin ...
1226! **************************************************************************************************
1227 SUBROUTINE set2zero_rho_atom_rad(rho_atom_set, iatom, ispin)
1228
1229 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
1230 INTEGER, INTENT(IN) :: iatom, ispin
1231
1232 INTEGER :: j
1233
1234 rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef = 0.0_dp
1235 rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef = 0.0_dp
1236
1237 rho_atom_set(iatom)%vrho_rad_h(ispin)%r_coef = 0.0_dp
1238 rho_atom_set(iatom)%vrho_rad_s(ispin)%r_coef = 0.0_dp
1239
1240 rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef = 0.0_dp
1241 rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef = 0.0_dp
1242
1243 DO j = 1, 3
1244 rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef = 0.0_dp
1245 rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef = 0.0_dp
1246 END DO
1247
1248 END SUBROUTINE set2zero_rho_atom_rad
1249
1250END MODULE qs_rho_atom_methods
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered)
Retrieve information from a kpoint environment.
Generation of the spherical Lebedev grids. All Lebedev grids were generated with a precision of at le...
Definition lebedev.F:57
subroutine, public deallocate_lebedev_grids()
...
Definition lebedev.F:324
type(oh_grid), dimension(nlg), target, public lebedev_grid
Definition lebedev.F:85
integer function, public get_number_of_lebedev_grid(l, n)
Get the number of the Lebedev grid, which has the requested angular momentum quantnum number l or siz...
Definition lebedev.F:114
subroutine, public init_lebedev_grids()
Load the coordinates and weights of the nonredundant Lebedev grid points.
Definition lebedev.F:344
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
real(kind=dp), parameter, public fourpi
Utility routines for the memory handling.
Interface to the message passing library MPI.
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public nsoset
integer, dimension(:, :), allocatable, public indso
subroutine, public get_paw_basis_info(basis_1c, o2nindex, n2oindex, nsatbas)
Return some info on the PAW basis derived from a GTO basis set.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
subroutine, public create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
...
subroutine, public get_maxl_cg(harmonics, orb_basis, llmax, max_s_harm)
...
subroutine, public create_harmonics_atom(harmonics, my_cg, na, llmax, maxs, max_s_harm, ll, wa, azi, pol)
...
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Define the neighbor list data types and the corresponding functionality.
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Routines for the construction of the coefficients for the expansion of the atomic densities rho1_hard...
subroutine, public proj_blk(h_a, s_a, na, h_b, s_b, nb, blk, ldb, proj_h, proj_s, nso, len1, len2, fac, distab, work1, work2)
Project a matrix block onto the local atomic functions.
subroutine, public calculate_rho_atom(para_env, rho_atom_set, qs_kind, atom_list, natom, nspins, tot_rho1_h, tot_rho1_s, rho1_h_spin, rho1_s_spin, rho1_h_aspin, rho1_s_aspin)
...
subroutine, public allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
...
subroutine, public replicate_rho_atom_radial(para_env, rho_atom_set, qs_kind, atom_list, natom, nspins)
Replicate the radial hard/soft density data needed to evaluate one-center tails on rank-local target ...
subroutine, public init_rho_atom(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
...
subroutine, public calculate_rho_atom_coeff(qs_env, rho_ao, rho_atom_set, qs_kind_set, oce, sab, para_env)
...
subroutine, public deallocate_rho_atom_set(rho_atom_set)
...
General overlap type integrals containers.
subroutine, public alist_pre_align_blk(blk_in, ldin, blk_out, ldout, ilist, in, jlist, jn)
...
subroutine, public get_alist(sap_int, alist, atom)
...
Calculate spherical harmonics.
subroutine, public clebsch_gordon_init(l)
...
subroutine, public clebsch_gordon_deallocate()
...
All kind of helpful little routines.
Definition util.F:14
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition util.F:333
Calculates special integrals.
Definition whittaker.F:12
subroutine, public whittaker_c0a(wc, r, expa, erfa, alpha, l1, l2, n)
int(y^(2+l1+l2) * exp(-alpha*y*y),y=0..x) / x^(l2+1); wc(:) :: output r(:) :: coordinate expa(:) :: e...
Definition whittaker.F:52
subroutine, public whittaker_ci(wc, r, expa, alpha, l, n)
int(y^(l+1) * exp(-alpha*y*y),y=x..infinity);
Definition whittaker.F:340
Provides all information about an atomic kind.
Contains information about kpoints.
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.