(git:b195825)
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-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 MODULE qs_rho_atom_methods
8 
9  USE atomic_kind_types, ONLY: atomic_kind_type,&
10  get_atomic_kind,&
11  get_atomic_kind_set
12  USE basis_set_types, ONLY: get_gto_basis_set,&
13  gto_basis_set_p_type,&
14  gto_basis_set_type
15  USE cp_control_types, ONLY: dft_control_type,&
16  gapw_control_type
17  USE dbcsr_api, ONLY: dbcsr_get_block_p,&
18  dbcsr_p_type
19  USE kinds, ONLY: dp
20  USE kpoint_types, ONLY: get_kpoint_info,&
21  kpoint_type
22  USE lebedev, ONLY: deallocate_lebedev_grids,&
23  get_number_of_lebedev_grid,&
24  init_lebedev_grids,&
25  lebedev_grid
26  USE mathconstants, ONLY: fourpi,&
27  pi
28  USE memory_utilities, ONLY: reallocate
29  USE message_passing, ONLY: mp_para_env_type
30  USE orbital_pointers, ONLY: indso,&
31  nsoset
32  USE paw_basis_types, ONLY: get_paw_basis_info
33  USE qs_environment_types, ONLY: get_qs_env,&
34  qs_environment_type
35  USE qs_grid_atom, ONLY: create_grid_atom,&
36  grid_atom_type
37  USE qs_harmonics_atom, ONLY: create_harmonics_atom,&
38  get_maxl_cg,&
39  get_none0_cg_list,&
40  harmonics_atom_type
41  USE qs_kind_types, ONLY: get_qs_kind,&
42  get_qs_kind_set,&
43  qs_kind_type
44  USE qs_neighbor_list_types, ONLY: get_iterator_info,&
45  neighbor_list_iterate,&
46  neighbor_list_iterator_create,&
47  neighbor_list_iterator_p_type,&
48  neighbor_list_iterator_release,&
49  neighbor_list_set_p_type
50  USE qs_oce_methods, ONLY: proj_blk
51  USE qs_oce_types, ONLY: oce_matrix_type
52  USE qs_rho_atom_types, ONLY: deallocate_rho_atom_set,&
53  rho_atom_coeff,&
54  rho_atom_type
55  USE sap_kind_types, ONLY: alist_pre_align_blk,&
56  alist_type,&
57  get_alist
58  USE spherical_harmonics, ONLY: clebsch_gordon,&
59  clebsch_gordon_deallocate,&
60  clebsch_gordon_init
61  USE util, ONLY: get_limit
62  USE whittaker, ONLY: whittaker_c0a,&
63  whittaker_ci
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 
83  PUBLIC :: allocate_rho_atom_internals, &
84  calculate_rho_atom, &
85  calculate_rho_atom_coeff, &
86  init_rho_atom
87 
88 CONTAINS
89 
90 ! **************************************************************************************************
91 !> \brief ...
92 !> \param para_env ...
93 !> \param rho_atom_set ...
94 !> \param qs_kind ...
95 !> \param atom_list ...
96 !> \param natom ...
97 !> \param nspins ...
98 !> \param tot_rho1_h ...
99 !> \param tot_rho1_s ...
100 ! **************************************************************************************************
101  SUBROUTINE calculate_rho_atom(para_env, rho_atom_set, qs_kind, atom_list, &
102  natom, nspins, tot_rho1_h, tot_rho1_s)
103 
104  TYPE(mp_para_env_type), POINTER :: para_env
105  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
106  TYPE(qs_kind_type), INTENT(IN) :: qs_kind
107  INTEGER, DIMENSION(:), INTENT(IN) :: atom_list
108  INTEGER, INTENT(IN) :: natom, nspins
109  REAL(dp), DIMENSION(:), INTENT(INOUT) :: tot_rho1_h, tot_rho1_s
110 
111  CHARACTER(len=*), PARAMETER :: routinen = 'calculate_rho_atom'
112 
113  INTEGER :: damax_iso_not0_local, handle, i, i1, i2, iat, iatom, icg, ipgf1, ipgf2, ir, &
114  iset1, iset2, iso, iso1, iso1_first, iso1_last, iso2, iso2_first, iso2_last, j, l, l1, &
115  l2, l_iso, l_sub, l_sum, lmax12, lmax_expansion, lmin12, m1s, m2s, max_iso_not0, &
116  max_iso_not0_local, max_s_harm, maxl, maxso, mepos, n1s, n2s, nr, nset, num_pe, size1, &
117  size2
118  INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list, dacg_n_list
119  INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list, dacg_list
120  INTEGER, DIMENSION(2) :: bo
121  INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, o2nindex
122  LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: done_vgg
123  REAL(dp) :: c1, c2, rho_h, rho_s, root_zet12, zet12
124  REAL(dp), ALLOCATABLE, DIMENSION(:) :: erf_zet12, g1, g2, gg0, int1, int2
125  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: cpch_sphere, cpcs_sphere, dgg, gg, gg_lm1
126  REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: vgg
127  REAL(dp), DIMENSION(:, :), POINTER :: coeff, zet
128  REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg
129  REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_cg_dxyz
130  TYPE(grid_atom_type), POINTER :: grid_atom
131  TYPE(gto_basis_set_type), POINTER :: basis_1c
132  TYPE(harmonics_atom_type), POINTER :: harmonics
133 
134  CALL timeset(routinen, handle)
135 
136  !Note: tau is taken care of separately in qs_vxc_atom.F
137 
138  NULLIFY (basis_1c)
139  NULLIFY (harmonics, grid_atom)
140  NULLIFY (lmin, lmax, npgf, zet, my_cg, my_cg_dxyz, coeff)
141 
142  CALL get_qs_kind(qs_kind, grid_atom=grid_atom, harmonics=harmonics)
143  CALL get_qs_kind(qs_kind, basis_set=basis_1c, basis_type="GAPW_1C")
144 
145  CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
146  maxl=maxl, npgf=npgf, nset=nset, zet=zet, &
147  maxso=maxso)
148 
149  CALL get_paw_basis_info(basis_1c, o2nindex=o2nindex)
150 
151  max_iso_not0 = harmonics%max_iso_not0
152  max_s_harm = harmonics%max_s_harm
153 
154  nr = grid_atom%nr
155  lmax_expansion = indso(1, max_iso_not0)
156  ! Distribute the atoms of this kind
157  num_pe = para_env%num_pe
158  mepos = para_env%mepos
159  bo = get_limit(natom, num_pe, mepos)
160 
161  my_cg => harmonics%my_CG
162  my_cg_dxyz => harmonics%my_CG_dxyz
163 
164  ALLOCATE (cpch_sphere(nsoset(maxl), nsoset(maxl)))
165  ALLOCATE (cpcs_sphere(nsoset(maxl), nsoset(maxl)))
166  ALLOCATE (g1(nr), g2(nr), gg0(nr), gg(nr, 0:2*maxl), dgg(nr, 0:2*maxl), gg_lm1(nr, 0:2*maxl))
167  ALLOCATE (erf_zet12(nr), vgg(nr, 0:2*maxl, 0:indso(1, max_iso_not0)))
168  ALLOCATE (done_vgg(0:2*maxl, 0:indso(1, max_iso_not0)))
169  ALLOCATE (int1(nr), int2(nr))
170  ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm), &
171  dacg_list(2, nsoset(maxl)**2, max_s_harm), dacg_n_list(max_s_harm))
172 
173  DO iat = bo(1), bo(2)
174  iatom = atom_list(iat)
175  DO i = 1, nspins
176  IF (.NOT. ASSOCIATED(rho_atom_set(iatom)%rho_rad_h(i)%r_coef)) THEN
177  CALL allocate_rho_atom_rad(rho_atom_set, iatom, i, nr, max_iso_not0)
178  ELSE
179  CALL set2zero_rho_atom_rad(rho_atom_set, iatom, i)
180  END IF
181  END DO
182  END DO
183 
184  m1s = 0
185  DO iset1 = 1, nset
186  m2s = 0
187  DO iset2 = 1, nset
188 
189  CALL get_none0_cg_list(my_cg, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
190  max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
191  cpassert(max_iso_not0_local .LE. max_iso_not0)
192  CALL get_none0_cg_list(my_cg_dxyz, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
193  max_s_harm, lmax_expansion, dacg_list, dacg_n_list, damax_iso_not0_local)
194  n1s = nsoset(lmax(iset1))
195 
196  DO ipgf1 = 1, npgf(iset1)
197  iso1_first = nsoset(lmin(iset1) - 1) + 1 + n1s*(ipgf1 - 1) + m1s
198  iso1_last = nsoset(lmax(iset1)) + n1s*(ipgf1 - 1) + m1s
199  size1 = iso1_last - iso1_first + 1
200  iso1_first = o2nindex(iso1_first)
201  iso1_last = o2nindex(iso1_last)
202  i1 = iso1_last - iso1_first + 1
203  cpassert(size1 == i1)
204  i1 = nsoset(lmin(iset1) - 1) + 1
205 
206  g1(1:nr) = exp(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
207 
208  n2s = nsoset(lmax(iset2))
209  DO ipgf2 = 1, npgf(iset2)
210  iso2_first = nsoset(lmin(iset2) - 1) + 1 + n2s*(ipgf2 - 1) + m2s
211  iso2_last = nsoset(lmax(iset2)) + n2s*(ipgf2 - 1) + m2s
212  size2 = iso2_last - iso2_first + 1
213  iso2_first = o2nindex(iso2_first)
214  iso2_last = o2nindex(iso2_last)
215  i2 = iso2_last - iso2_first + 1
216  cpassert(size2 == i2)
217  i2 = nsoset(lmin(iset2) - 1) + 1
218 
219  g2(1:nr) = exp(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
220  lmin12 = lmin(iset1) + lmin(iset2)
221  lmax12 = lmax(iset1) + lmax(iset2)
222 
223  zet12 = zet(ipgf1, iset1) + zet(ipgf2, iset2)
224  root_zet12 = sqrt(zet(ipgf1, iset1) + zet(ipgf2, iset2))
225  DO ir = 1, nr
226  erf_zet12(ir) = erf(root_zet12*grid_atom%rad(ir))
227  END DO
228 
229  gg = 0.0_dp
230  dgg = 0.0_dp
231  gg_lm1 = 0.0_dp
232  vgg = 0.0_dp
233  done_vgg = .false.
234  ! reduce the number of terms in the expansion local densities
235  IF (lmin12 .LE. lmax_expansion) THEN
236  IF (lmin12 == 0) THEN
237  gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
238  gg_lm1(1:nr, lmin12) = 0.0_dp
239  gg0(1:nr) = gg(1:nr, lmin12)
240  ELSE
241  gg0(1:nr) = g1(1:nr)*g2(1:nr)
242  gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
243  gg_lm1(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12 - 1)*g1(1:nr)*g2(1:nr)
244  END IF
245 
246  ! reduce the number of terms in the expansion local densities
247  IF (lmax12 .GT. lmax_expansion) lmax12 = lmax_expansion
248 
249  DO l = lmin12 + 1, lmax12
250  gg(1:nr, l) = grid_atom%rad(1:nr)*gg(1:nr, l - 1)
251  gg_lm1(1:nr, l) = gg(1:nr, l - 1)
252  dgg(1:nr, l - 1) = -2.0_dp*(zet(ipgf1, iset1) + zet(ipgf2, iset2))*gg(1:nr, l)
253 
254  END DO
255  dgg(1:nr, lmax12) = -2.0_dp*(zet(ipgf1, iset1) + &
256  zet(ipgf2, iset2))*grid_atom%rad(1:nr)*gg(1:nr, lmax12)
257 
258  c2 = sqrt(pi*pi*pi/(zet12*zet12*zet12))
259 
260  DO iso = 1, max_iso_not0_local
261  l_iso = indso(1, iso)
262  c1 = fourpi/(2._dp*real(l_iso, dp) + 1._dp)
263  DO icg = 1, cg_n_list(iso)
264  iso1 = cg_list(1, icg, iso)
265  iso2 = cg_list(2, icg, iso)
266 
267  l = indso(1, iso1) + indso(1, iso2)
268  cpassert(l <= lmax_expansion)
269  IF (done_vgg(l, l_iso)) cycle
270  l_sum = l + l_iso
271  l_sub = l - l_iso
272 
273  IF (l_sum == 0) THEN
274  vgg(1:nr, l, l_iso) = erf_zet12(1:nr)*grid_atom%oorad2l(1:nr, 1)*c2
275  ELSE
276  CALL whittaker_c0a(int1, grid_atom%rad, gg0, erf_zet12, zet12, l, l_iso, nr)
277  CALL whittaker_ci(int2, grid_atom%rad, gg0, zet12, l_sub, nr)
278 
279  DO ir = 1, nr
280  int2(ir) = grid_atom%rad2l(ir, l_iso)*int2(ir)
281  vgg(ir, l, l_iso) = c1*(int1(ir) + int2(ir))
282  END DO
283  END IF
284  done_vgg(l, l_iso) = .true.
285  END DO
286  END DO
287  END IF ! lmax_expansion
288 
289  DO iat = bo(1), bo(2)
290  iatom = atom_list(iat)
291 
292  DO i = 1, nspins
293  cpch_sphere = 0.0_dp
294  cpcs_sphere = 0.0_dp
295 
296  coeff => rho_atom_set(iatom)%cpc_h(i)%r_coef
297  cpch_sphere(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = coeff(iso1_first:iso1_last, iso2_first:iso2_last)
298 
299  coeff => rho_atom_set(iatom)%cpc_s(i)%r_coef
300  cpcs_sphere(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = coeff(iso1_first:iso1_last, iso2_first:iso2_last)
301 
302  DO iso = 1, max_iso_not0_local
303  l_iso = indso(1, iso)
304  DO icg = 1, cg_n_list(iso)
305  iso1 = cg_list(1, icg, iso)
306  iso2 = cg_list(2, icg, iso)
307 
308  l1 = indso(1, iso1)
309  l2 = indso(1, iso2)
310 
311  l = indso(1, iso1) + indso(1, iso2)
312  cpassert(l <= lmax_expansion)
313 
314  rho_atom_set(iatom)%rho_rad_h(i)%r_coef(1:nr, iso) = &
315  rho_atom_set(iatom)%rho_rad_h(i)%r_coef(1:nr, iso) + &
316  gg(1:nr, l)*cpch_sphere(iso1, iso2)*my_cg(iso1, iso2, iso)
317 
318  rho_atom_set(iatom)%rho_rad_s(i)%r_coef(1:nr, iso) = &
319  rho_atom_set(iatom)%rho_rad_s(i)%r_coef(1:nr, iso) + &
320  gg(1:nr, l)*cpcs_sphere(iso1, iso2)*my_cg(iso1, iso2, iso)
321 
322  rho_atom_set(iatom)%drho_rad_h(i)%r_coef(1:nr, iso) = &
323  rho_atom_set(iatom)%drho_rad_h(i)%r_coef(1:nr, iso) + &
324  dgg(1:nr, l)*cpch_sphere(iso1, iso2)*my_cg(iso1, iso2, iso)
325 
326  rho_atom_set(iatom)%drho_rad_s(i)%r_coef(1:nr, iso) = &
327  rho_atom_set(iatom)%drho_rad_s(i)%r_coef(1:nr, iso) + &
328  dgg(1:nr, l)*cpcs_sphere(iso1, iso2)*my_cg(iso1, iso2, iso)
329 
330  rho_atom_set(iatom)%vrho_rad_h(i)%r_coef(1:nr, iso) = &
331  rho_atom_set(iatom)%vrho_rad_h(i)%r_coef(1:nr, iso) + &
332  vgg(1:nr, l, l_iso)*cpch_sphere(iso1, iso2)*my_cg(iso1, iso2, iso)
333 
334  rho_atom_set(iatom)%vrho_rad_s(i)%r_coef(1:nr, iso) = &
335  rho_atom_set(iatom)%vrho_rad_s(i)%r_coef(1:nr, iso) + &
336  vgg(1:nr, l, l_iso)*cpcs_sphere(iso1, iso2)*my_cg(iso1, iso2, iso)
337 
338  END DO ! icg
339 
340  END DO ! iso
341 
342  DO iso = 1, max_iso_not0 !damax_iso_not0_local
343  l_iso = indso(1, iso)
344  DO icg = 1, dacg_n_list(iso)
345  iso1 = dacg_list(1, icg, iso)
346  iso2 = dacg_list(2, icg, iso)
347  l = indso(1, iso1) + indso(1, iso2)
348  cpassert(l <= lmax_expansion)
349  DO j = 1, 3
350  rho_atom_set(iatom)%rho_rad_h_d(j, i)%r_coef(1:nr, iso) = &
351  rho_atom_set(iatom)%rho_rad_h_d(j, i)%r_coef(1:nr, iso) + &
352  gg_lm1(1:nr, l)*cpch_sphere(iso1, iso2)*my_cg_dxyz(j, iso1, iso2, iso)
353 
354  rho_atom_set(iatom)%rho_rad_s_d(j, i)%r_coef(1:nr, iso) = &
355  rho_atom_set(iatom)%rho_rad_s_d(j, i)%r_coef(1:nr, iso) + &
356  gg_lm1(1:nr, l)*cpcs_sphere(iso1, iso2)*my_cg_dxyz(j, iso1, iso2, iso)
357  END DO
358  END DO ! icg
359 
360  END DO ! iso
361 
362  END DO ! i
363  END DO ! iat
364 
365  END DO ! ipgf2
366  END DO ! ipgf1
367  m2s = m2s + maxso
368  END DO ! iset2
369  m1s = m1s + maxso
370  END DO ! iset1
371 
372  DO iat = bo(1), bo(2)
373  iatom = atom_list(iat)
374 
375  DO i = 1, nspins
376 
377  DO iso = 1, max_iso_not0
378  rho_s = 0.0_dp
379  rho_h = 0.0_dp
380  DO ir = 1, nr
381  rho_h = rho_h + rho_atom_set(iatom)%rho_rad_h(i)%r_coef(ir, iso)*grid_atom%wr(ir)
382  rho_s = rho_s + rho_atom_set(iatom)%rho_rad_s(i)%r_coef(ir, iso)*grid_atom%wr(ir)
383  END DO ! ir
384  tot_rho1_h(i) = tot_rho1_h(i) + rho_h*harmonics%slm_int(iso)
385  tot_rho1_s(i) = tot_rho1_s(i) + rho_s*harmonics%slm_int(iso)
386  END DO ! iso
387 
388  END DO ! ispin
389 
390  END DO ! iat
391 
392  DEALLOCATE (cpch_sphere, cpcs_sphere)
393  DEALLOCATE (g1, g2, gg0, gg, gg_lm1, dgg, vgg, done_vgg, erf_zet12, int1, int2)
394  DEALLOCATE (cg_list, cg_n_list, dacg_list, dacg_n_list)
395  DEALLOCATE (o2nindex)
396 
397  CALL timestop(handle)
398 
399  END SUBROUTINE calculate_rho_atom
400 
401 ! **************************************************************************************************
402 !> \brief ...
403 !> \param qs_env QuickStep environment
404 !> (accessed components: atomic_kind_set, dft_control%nimages,
405 !> dft_control%nspins, kpoints%cell_to_index)
406 !> \param rho_ao density matrix in atomic basis set
407 !> \param rho_atom_set ...
408 !> \param qs_kind_set list of QuickStep kinds
409 !> \param oce one-centre expansion coefficients
410 !> \param sab neighbour pair list
411 !> \param para_env parallel environment
412 !> \par History
413 !> Add OpenMP [Apr 2016, EPCC]
414 !> Use automatic arrays [Sep 2016, M Tucker]
415 !> Allow for external non-default kind_set, oce and sab [Dec 2019, A Bussy]
416 !> \note Consider to declare 'rho_ao' dummy argument as a pointer to the two-dimensional
417 !> (1:nspins, 1:nimages) set of matrices.
418 ! **************************************************************************************************
419  SUBROUTINE calculate_rho_atom_coeff(qs_env, rho_ao, rho_atom_set, qs_kind_set, oce, sab, para_env)
420 
421  TYPE(qs_environment_type), POINTER :: qs_env
422  TYPE(dbcsr_p_type), DIMENSION(*) :: rho_ao
423  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
424  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
425  TYPE(oce_matrix_type), POINTER :: oce
426  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
427  POINTER :: sab
428  TYPE(mp_para_env_type), POINTER :: para_env
429 
430  CHARACTER(len=*), PARAMETER :: routinen = 'calculate_rho_atom_coeff'
431 
432  INTEGER :: bo(2), handle, i, iac, iatom, ibc, icol, ikind, img, irow, ispin, jatom, jkind, &
433  kac, katom, kbc, kkind, len_cpc, len_pc1, max_gau, max_nsgf, mepos, n_cont_a, n_cont_b, &
434  nat_kind, natom, nimages, nkind, nsatbas, nsoctot, nspins, num_pe
435  INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
436  INTEGER, DIMENSION(3) :: cell_b
437  INTEGER, DIMENSION(:), POINTER :: a_list, list_a, list_b
438  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
439  LOGICAL :: dista, distab, distb, found, paw_atom
440  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: p_matrix
441  REAL(kind=dp) :: eps_cpc, factor, pmax
442  REAL(kind=dp), DIMENSION(3) :: rab
443  REAL(kind=dp), DIMENSION(:, :), POINTER :: c_coeff_hh_a, c_coeff_hh_b, &
444  c_coeff_ss_a, c_coeff_ss_b, r_coef_h, &
445  r_coef_s
446  TYPE(alist_type), POINTER :: alist_ac, alist_bc
447  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
448  TYPE(dft_control_type), POINTER :: dft_control
449  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
450  TYPE(gto_basis_set_type), POINTER :: basis_1c, basis_set_a, basis_set_b
451  TYPE(kpoint_type), POINTER :: kpoints
452  TYPE(neighbor_list_iterator_p_type), &
453  DIMENSION(:), POINTER :: nl_iterator
454  TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: p_block_spin
455 
456 !$ INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: locks
457 !$ INTEGER :: lock, number_of_locks
458 
459  CALL timeset(routinen, handle)
460 
461  CALL get_qs_env(qs_env=qs_env, &
462  dft_control=dft_control, &
463  atomic_kind_set=atomic_kind_set)
464 
465  eps_cpc = dft_control%qs_control%gapw_control%eps_cpc
466 
467  cpassert(ASSOCIATED(qs_kind_set))
468  cpassert(ASSOCIATED(rho_atom_set))
469  cpassert(ASSOCIATED(oce))
470  cpassert(ASSOCIATED(sab))
471 
472  nspins = dft_control%nspins
473  nimages = dft_control%nimages
474 
475  NULLIFY (cell_to_index)
476  IF (nimages > 1) THEN
477  CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
478  CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
479  END IF
480 
481  CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
482  CALL get_qs_kind_set(qs_kind_set, maxsgf=max_nsgf, maxgtops=max_gau, basis_type='GAPW_1C')
483 
484  nkind = SIZE(atomic_kind_set)
485  ! Initialize to 0 the CPC coefficients and the local density arrays
486  DO ikind = 1, nkind
487  CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=a_list, natom=nat_kind)
488  CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
489 
490  IF (.NOT. paw_atom) cycle
491  DO i = 1, nat_kind
492  iatom = a_list(i)
493  DO ispin = 1, nspins
494  rho_atom_set(iatom)%cpc_h(ispin)%r_coef = 0.0_dp
495  rho_atom_set(iatom)%cpc_s(ispin)%r_coef = 0.0_dp
496  END DO ! ispin
497  END DO ! i
498 
499  num_pe = para_env%num_pe
500  mepos = para_env%mepos
501  bo = get_limit(nat_kind, num_pe, mepos)
502  DO i = bo(1), bo(2)
503  iatom = a_list(i)
504  DO ispin = 1, nspins
505  rho_atom_set(iatom)%ga_Vlocal_gb_h(ispin)%r_coef = 0.0_dp
506  rho_atom_set(iatom)%ga_Vlocal_gb_s(ispin)%r_coef = 0.0_dp
507  END DO ! ispin
508  END DO ! i
509  END DO ! ikind
510 
511  ALLOCATE (basis_set_list(nkind))
512  DO ikind = 1, nkind
513  CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a)
514  IF (ASSOCIATED(basis_set_a)) THEN
515  basis_set_list(ikind)%gto_basis_set => basis_set_a
516  ELSE
517  NULLIFY (basis_set_list(ikind)%gto_basis_set)
518  END IF
519  END DO
520 
521  len_pc1 = max_nsgf*max_gau
522  len_cpc = max_gau*max_gau
523 
524  num_pe = 1
525 !$ num_pe = omp_get_max_threads()
526  CALL neighbor_list_iterator_create(nl_iterator, sab, nthread=num_pe)
527 
528 !$OMP PARALLEL DEFAULT( NONE ) &
529 !$OMP SHARED( max_nsgf, max_gau &
530 !$OMP , len_PC1, len_CPC &
531 !$OMP , nl_iterator, basis_set_list &
532 !$OMP , nimages, cell_to_index &
533 !$OMP , nspins, rho_ao &
534 !$OMP , nkind, qs_kind_set &
535 !$OMP , oce, eps_cpc &
536 !$OMP , rho_atom_set &
537 !$OMP , natom, locks, number_of_locks &
538 !$OMP ) &
539 !$OMP PRIVATE( p_block_spin, ispin &
540 !$OMP , p_matrix, mepos &
541 !$OMP , ikind, jkind, iatom, jatom &
542 !$OMP , cell_b, rab, basis_1c &
543 !$OMP , basis_set_a, basis_set_b &
544 !$OMP , pmax, irow, icol, img &
545 !$OMP , found &
546 !$OMP , kkind &
547 !$OMP , paw_atom, nsatbas &
548 !$OMP , nsoctot, katom &
549 !$OMP , iac , alist_ac, kac, n_cont_a, list_a &
550 !$OMP , ibc , alist_bc, kbc, n_cont_b, list_b &
551 !$OMP , C_coeff_hh_a, C_coeff_ss_a, dista &
552 !$OMP , C_coeff_hh_b, C_coeff_ss_b, distb &
553 !$OMP , distab &
554 !$OMP , factor, r_coef_h, r_coef_s &
555 !$OMP )
556 
557  ALLOCATE (p_block_spin(nspins))
558  ALLOCATE (p_matrix(max_nsgf, max_nsgf))
559 
560 !$OMP SINGLE
561 !$ number_of_locks = nspins*natom
562 !$ ALLOCATE (locks(number_of_locks))
563 !$OMP END SINGLE
564 
565 !$OMP DO
566 !$ DO lock = 1, number_of_locks
567 !$ call omp_init_lock(locks(lock))
568 !$ END DO
569 !$OMP END DO
570 
571  mepos = 0
572 !$ mepos = omp_get_thread_num()
573  DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
574 
575  CALL get_iterator_info(nl_iterator, mepos=mepos, &
576  ikind=ikind, jkind=jkind, &
577  iatom=iatom, jatom=jatom, &
578  cell=cell_b, r=rab)
579 
580  basis_set_a => basis_set_list(ikind)%gto_basis_set
581  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
582  basis_set_b => basis_set_list(jkind)%gto_basis_set
583  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
584 
585  pmax = 0._dp
586  IF (iatom <= jatom) THEN
587  irow = iatom
588  icol = jatom
589  ELSE
590  irow = jatom
591  icol = iatom
592  END IF
593 
594  IF (nimages > 1) THEN
595  img = cell_to_index(cell_b(1), cell_b(2), cell_b(3))
596  cpassert(img > 0)
597  ELSE
598  img = 1
599  END IF
600 
601  DO ispin = 1, nspins
602  CALL dbcsr_get_block_p(matrix=rho_ao(nspins*(img - 1) + ispin)%matrix, &
603  row=irow, col=icol, block=p_block_spin(ispin)%r_coef, &
604  found=found)
605  pmax = pmax + maxval(abs(p_block_spin(ispin)%r_coef))
606  END DO
607 
608  DO kkind = 1, nkind
609  CALL get_qs_kind(qs_kind_set(kkind), basis_set=basis_1c, basis_type="GAPW_1C", &
610  paw_atom=paw_atom)
611 
612  IF (.NOT. paw_atom) cycle
613 
614  CALL get_paw_basis_info(basis_1c, nsatbas=nsatbas)
615  nsoctot = nsatbas
616 
617  iac = ikind + nkind*(kkind - 1)
618  ibc = jkind + nkind*(kkind - 1)
619  IF (.NOT. ASSOCIATED(oce%intac(iac)%alist)) cycle
620  IF (.NOT. ASSOCIATED(oce%intac(ibc)%alist)) cycle
621 
622  CALL get_alist(oce%intac(iac), alist_ac, iatom)
623  CALL get_alist(oce%intac(ibc), alist_bc, jatom)
624  IF (.NOT. ASSOCIATED(alist_ac)) cycle
625  IF (.NOT. ASSOCIATED(alist_bc)) cycle
626 
627  DO kac = 1, alist_ac%nclist
628  DO kbc = 1, alist_bc%nclist
629  IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
630  IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
631  IF (pmax*alist_bc%clist(kbc)%maxac*alist_ac%clist(kac)%maxac < eps_cpc) cycle
632 
633  n_cont_a = alist_ac%clist(kac)%nsgf_cnt
634  n_cont_b = alist_bc%clist(kbc)%nsgf_cnt
635  IF (n_cont_a .EQ. 0 .OR. n_cont_b .EQ. 0) cycle
636 
637  list_a => alist_ac%clist(kac)%sgf_list
638  list_b => alist_bc%clist(kbc)%sgf_list
639 
640  katom = alist_ac%clist(kac)%catom
641 
642  IF (iatom == katom .AND. all(alist_ac%clist(kac)%cell == 0)) THEN
643  c_coeff_hh_a => alist_ac%clist(kac)%achint(:, :, 1)
644  c_coeff_ss_a => alist_ac%clist(kac)%acint(:, :, 1)
645  dista = .false.
646  ELSE
647  c_coeff_hh_a => alist_ac%clist(kac)%acint(:, :, 1)
648  c_coeff_ss_a => alist_ac%clist(kac)%acint(:, :, 1)
649  dista = .true.
650  END IF
651  IF (jatom == katom .AND. all(alist_bc%clist(kbc)%cell == 0)) THEN
652  c_coeff_hh_b => alist_bc%clist(kbc)%achint(:, :, 1)
653  c_coeff_ss_b => alist_bc%clist(kbc)%acint(:, :, 1)
654  distb = .false.
655  ELSE
656  c_coeff_hh_b => alist_bc%clist(kbc)%acint(:, :, 1)
657  c_coeff_ss_b => alist_bc%clist(kbc)%acint(:, :, 1)
658  distb = .true.
659  END IF
660 
661  distab = dista .AND. distb
662 
663  DO ispin = 1, nspins
664 
665  IF (iatom <= jatom) THEN
666  CALL alist_pre_align_blk(p_block_spin(ispin)%r_coef, &
667  SIZE(p_block_spin(ispin)%r_coef, 1), p_matrix, SIZE(p_matrix, 1), &
668  list_a, n_cont_a, list_b, n_cont_b)
669  ELSE
670  CALL alist_pre_align_blk(p_block_spin(ispin)%r_coef, &
671  SIZE(p_block_spin(ispin)%r_coef, 1), p_matrix, SIZE(p_matrix, 1), &
672  list_b, n_cont_b, list_a, n_cont_a)
673  END IF
674 
675  factor = 1.0_dp
676  IF (iatom == jatom) factor = 0.5_dp
677 
678  r_coef_h => rho_atom_set(katom)%cpc_h(ispin)%r_coef
679  r_coef_s => rho_atom_set(katom)%cpc_s(ispin)%r_coef
680 
681 !$ CALL omp_set_lock(locks((katom - 1)*nspins + ispin))
682  IF (iatom <= jatom) THEN
683  CALL proj_blk(c_coeff_hh_a, c_coeff_ss_a, n_cont_a, &
684  c_coeff_hh_b, c_coeff_ss_b, n_cont_b, &
685  p_matrix, max_nsgf, r_coef_h, r_coef_s, nsoctot, &
686  len_pc1, len_cpc, factor, distab)
687  ELSE
688  CALL proj_blk(c_coeff_hh_b, c_coeff_ss_b, n_cont_b, &
689  c_coeff_hh_a, c_coeff_ss_a, n_cont_a, &
690  p_matrix, max_nsgf, r_coef_h, r_coef_s, nsoctot, &
691  len_pc1, len_cpc, factor, distab)
692  END IF
693 !$ CALL omp_unset_lock(locks((katom - 1)*nspins + ispin))
694 
695  END DO
696  EXIT !search loop over jatom-katom list
697  END IF
698  END DO
699  END DO
700  END DO
701  END DO
702  ! Wait for all threads to finish the loop before locks can be freed
703 !$OMP BARRIER
704 
705 !$OMP DO
706 !$ DO lock = 1, number_of_locks
707 !$ call omp_destroy_lock(locks(lock))
708 !$ END DO
709 !$OMP END DO
710 !$OMP SINGLE
711 !$ DEALLOCATE (locks)
712 !$OMP END SINGLE NOWAIT
713 
714  DEALLOCATE (p_block_spin, p_matrix)
715 !$OMP END PARALLEL
716 
717  CALL neighbor_list_iterator_release(nl_iterator)
718 
719  CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
720 
721  DO iatom = 1, natom
722  ikind = kind_of(iatom)
723 
724  DO ispin = 1, nspins
725  IF (ASSOCIATED(rho_atom_set(iatom)%cpc_h(ispin)%r_coef)) THEN
726  CALL para_env%sum(rho_atom_set(iatom)%cpc_h(ispin)%r_coef)
727  CALL para_env%sum(rho_atom_set(iatom)%cpc_s(ispin)%r_coef)
728  r_coef_h => rho_atom_set(iatom)%cpc_h(ispin)%r_coef
729  r_coef_s => rho_atom_set(iatom)%cpc_s(ispin)%r_coef
730  r_coef_h(:, :) = r_coef_h(:, :) + transpose(r_coef_h(:, :))
731  r_coef_s(:, :) = r_coef_s(:, :) + transpose(r_coef_s(:, :))
732  END IF
733  END DO
734 
735  END DO
736 
737  DEALLOCATE (kind_of, basis_set_list)
738 
739  CALL timestop(handle)
740 
741  END SUBROUTINE calculate_rho_atom_coeff
742 
743 ! **************************************************************************************************
744 !> \brief ...
745 !> \param rho_atom_set the type to initialize
746 !> \param atomic_kind_set list of atomic kinds
747 !> \param qs_kind_set the kind set from which to take quantum numbers and basis info
748 !> \param dft_control DFT control type
749 !> \param para_env parallel environment
750 !> \par History:
751 !> - Generalised by providing the rho_atom_set and the qs_kind_set 12.2019 (A.Bussy)
752 ! **************************************************************************************************
753  SUBROUTINE init_rho_atom(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
754 
755  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
756  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
757  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
758  TYPE(dft_control_type), POINTER :: dft_control
759  TYPE(mp_para_env_type), POINTER :: para_env
760 
761  CHARACTER(len=*), PARAMETER :: routinen = 'init_rho_atom'
762 
763  INTEGER :: handle, ikind, il, iso, iso1, iso2, l1, l1l2, l2, la, lc1, lc2, lcleb, ll, llmax, &
764  lmax_sphere, lp, m1, m2, max_s_harm, max_s_set, maxl, maxlgto, maxs, mm, mp, na, nat, &
765  natom, nr, nspins, quadrature
766  INTEGER, DIMENSION(:), POINTER :: atom_list
767  LOGICAL :: paw_atom
768  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rga
769  REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg
770  TYPE(gapw_control_type), POINTER :: gapw_control
771  TYPE(grid_atom_type), POINTER :: grid_atom
772  TYPE(gto_basis_set_type), POINTER :: basis_1c_set
773  TYPE(harmonics_atom_type), POINTER :: harmonics
774 
775  CALL timeset(routinen, handle)
776 
777  NULLIFY (basis_1c_set)
778  NULLIFY (my_cg, grid_atom, harmonics, atom_list)
779 
780  cpassert(ASSOCIATED(atomic_kind_set))
781  cpassert(ASSOCIATED(dft_control))
782  cpassert(ASSOCIATED(para_env))
783  cpassert(ASSOCIATED(qs_kind_set))
784 
785  CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
786 
787  CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="GAPW_1C")
788 
789  nspins = dft_control%nspins
790  gapw_control => dft_control%qs_control%gapw_control
791 
792  lmax_sphere = gapw_control%lmax_sphere
793 
794  llmax = min(lmax_sphere, 2*maxlgto)
795  max_s_harm = nsoset(llmax)
796  max_s_set = nsoset(maxlgto)
797 
798  lcleb = max(llmax, 2*maxlgto, 1)
799 
800 ! *** allocate calculate the CG coefficients up to the maxl ***
801  CALL clebsch_gordon_init(lcleb)
802  CALL reallocate(my_cg, 1, max_s_set, 1, max_s_set, 1, max_s_harm)
803 
804  ALLOCATE (rga(lcleb, 2))
805  DO lc1 = 0, maxlgto
806  DO iso1 = nsoset(lc1 - 1) + 1, nsoset(lc1)
807  l1 = indso(1, iso1)
808  m1 = indso(2, iso1)
809  DO lc2 = 0, maxlgto
810  DO iso2 = nsoset(lc2 - 1) + 1, nsoset(lc2)
811  l2 = indso(1, iso2)
812  m2 = indso(2, iso2)
813  CALL clebsch_gordon(l1, m1, l2, m2, rga)
814  IF (l1 + l2 > llmax) THEN
815  l1l2 = llmax
816  ELSE
817  l1l2 = l1 + l2
818  END IF
819  mp = m1 + m2
820  mm = m1 - m2
821  IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN
822  mp = -abs(mp)
823  mm = -abs(mm)
824  ELSE
825  mp = abs(mp)
826  mm = abs(mm)
827  END IF
828  DO lp = mod(l1 + l2, 2), l1l2, 2
829  il = lp/2 + 1
830  IF (abs(mp) <= lp) THEN
831  IF (mp >= 0) THEN
832  iso = nsoset(lp - 1) + lp + 1 + mp
833  ELSE
834  iso = nsoset(lp - 1) + lp + 1 - abs(mp)
835  END IF
836  my_cg(iso1, iso2, iso) = rga(il, 1)
837  END IF
838  IF (mp /= mm .AND. abs(mm) <= lp) THEN
839  IF (mm >= 0) THEN
840  iso = nsoset(lp - 1) + lp + 1 + mm
841  ELSE
842  iso = nsoset(lp - 1) + lp + 1 - abs(mm)
843  END IF
844  my_cg(iso1, iso2, iso) = rga(il, 2)
845  END IF
846  END DO
847  END DO ! iso2
848  END DO ! lc2
849  END DO ! iso1
850  END DO ! lc1
851  DEALLOCATE (rga)
852  CALL clebsch_gordon_deallocate()
853 
854 ! *** initialize the Lebedev grids ***
855  CALL init_lebedev_grids()
856  quadrature = gapw_control%quadrature
857 
858  DO ikind = 1, SIZE(atomic_kind_set)
859  CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
860  CALL get_qs_kind(qs_kind_set(ikind), &
861  paw_atom=paw_atom, &
862  grid_atom=grid_atom, &
863  harmonics=harmonics, &
864  ngrid_rad=nr, ngrid_ang=na)
865 
866 ! *** determine the Lebedev grid for this kind ***
867 
868  ll = get_number_of_lebedev_grid(n=na)
869  na = lebedev_grid(ll)%n
870  la = lebedev_grid(ll)%l
871  grid_atom%ng_sphere = na
872  grid_atom%nr = nr
873 
874  IF (llmax > la) THEN
875  WRITE (*, '(/,72("*"))')
876  WRITE (*, '(T2,A,T66,I4)') &
877  "WARNING: the lebedev grid is built for angular momentum l up to ", la, &
878  " the max l of spherical harmonics is larger, l_max = ", llmax, &
879  " good integration is guaranteed only for l <= ", la
880  WRITE (*, '(72("*"),/)')
881  END IF
882 
883 ! *** calculate the radial grid ***
884  CALL create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
885 
886 ! *** calculate the spherical harmonics on the grid ***
887 
888  CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_1c_set, basis_type="GAPW_1C")
889  CALL get_gto_basis_set(gto_basis_set=basis_1c_set, maxl=maxl)
890  maxs = nsoset(maxl)
891  CALL create_harmonics_atom(harmonics, &
892  my_cg, na, llmax, maxs, max_s_harm, ll, grid_atom%wa, &
893  grid_atom%azi, grid_atom%pol)
894  CALL get_maxl_cg(harmonics, basis_1c_set, llmax, max_s_harm)
895 
896  END DO
897 
898  CALL deallocate_lebedev_grids()
899  DEALLOCATE (my_cg)
900 
901  CALL allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
902 
903  CALL timestop(handle)
904 
905  END SUBROUTINE init_rho_atom
906 
907 ! **************************************************************************************************
908 !> \brief ...
909 !> \param rho_atom_set ...
910 !> \param atomic_kind_set list of atomic kinds
911 !> \param qs_kind_set the kind set from which to take quantum numbers and basis info
912 !> \param dft_control DFT control type
913 !> \param para_env parallel environment
914 ! **************************************************************************************************
915  SUBROUTINE allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
916 
917  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
918  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
919  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
920  TYPE(dft_control_type), POINTER :: dft_control
921  TYPE(mp_para_env_type), POINTER :: para_env
922 
923  CHARACTER(len=*), PARAMETER :: routinen = 'allocate_rho_atom_internals'
924 
925  INTEGER :: bo(2), handle, iat, iatom, ikind, ispin, &
926  max_iso_not0, maxso, mepos, nat, &
927  natom, nsatbas, nset, nsotot, nspins, &
928  num_pe
929  INTEGER, DIMENSION(:), POINTER :: atom_list
930  LOGICAL :: paw_atom
931  TYPE(gto_basis_set_type), POINTER :: basis_1c
932  TYPE(harmonics_atom_type), POINTER :: harmonics
933 
934  CALL timeset(routinen, handle)
935 
936  cpassert(ASSOCIATED(atomic_kind_set))
937  cpassert(ASSOCIATED(dft_control))
938  cpassert(ASSOCIATED(para_env))
939  cpassert(ASSOCIATED(qs_kind_set))
940 
941  CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
942 
943  nspins = dft_control%nspins
944 
945  IF (ASSOCIATED(rho_atom_set)) THEN
946  CALL deallocate_rho_atom_set(rho_atom_set)
947  END IF
948  ALLOCATE (rho_atom_set(natom))
949 
950  DO ikind = 1, SIZE(atomic_kind_set)
951 
952  NULLIFY (atom_list, harmonics)
953  CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
954  CALL get_qs_kind(qs_kind_set(ikind), &
955  paw_atom=paw_atom, &
956  harmonics=harmonics)
957 
958  IF (paw_atom) THEN
959  CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
960  CALL get_gto_basis_set(gto_basis_set=basis_1c, nset=nset, maxso=maxso)
961  nsotot = nset*maxso
962  CALL get_paw_basis_info(basis_1c, nsatbas=nsatbas)
963  END IF
964 
965  max_iso_not0 = harmonics%max_iso_not0
966  DO iat = 1, nat
967  iatom = atom_list(iat)
968  ! *** allocate the radial density for each LM,for each atom ***
969 
970  ALLOCATE (rho_atom_set(iatom)%rho_rad_h(nspins))
971  ALLOCATE (rho_atom_set(iatom)%rho_rad_s(nspins))
972  ALLOCATE (rho_atom_set(iatom)%vrho_rad_h(nspins))
973  ALLOCATE (rho_atom_set(iatom)%vrho_rad_s(nspins))
974 
975  ALLOCATE (rho_atom_set(iatom)%cpc_h(nspins), &
976  rho_atom_set(iatom)%cpc_s(nspins), &
977  rho_atom_set(iatom)%drho_rad_h(nspins), &
978  rho_atom_set(iatom)%drho_rad_s(nspins), &
979  rho_atom_set(iatom)%rho_rad_h_d(3, nspins), &
980  rho_atom_set(iatom)%rho_rad_s_d(3, nspins))
981  ALLOCATE (rho_atom_set(iatom)%int_scr_h(nspins), &
982  rho_atom_set(iatom)%int_scr_s(nspins))
983 
984  IF (paw_atom) THEN
985  DO ispin = 1, nspins
986  ALLOCATE (rho_atom_set(iatom)%cpc_h(ispin)%r_coef(1:nsatbas, 1:nsatbas), &
987  rho_atom_set(iatom)%cpc_s(ispin)%r_coef(1:nsatbas, 1:nsatbas))
988  ALLOCATE (rho_atom_set(iatom)%int_scr_h(ispin)%r_coef(1:nsatbas, 1:nsatbas), &
989  rho_atom_set(iatom)%int_scr_s(ispin)%r_coef(1:nsatbas, 1:nsatbas))
990 
991  rho_atom_set(iatom)%cpc_h(ispin)%r_coef = 0.0_dp
992  rho_atom_set(iatom)%cpc_s(ispin)%r_coef = 0.0_dp
993  END DO
994  END IF
995 
996  END DO ! iat
997 
998  num_pe = para_env%num_pe
999  mepos = para_env%mepos
1000  bo = get_limit(nat, num_pe, mepos)
1001  DO iat = bo(1), bo(2)
1002  iatom = atom_list(iat)
1003  ALLOCATE (rho_atom_set(iatom)%ga_Vlocal_gb_h(nspins), &
1004  rho_atom_set(iatom)%ga_Vlocal_gb_s(nspins))
1005  IF (paw_atom) THEN
1006  DO ispin = 1, nspins
1007  CALL reallocate(rho_atom_set(iatom)%ga_Vlocal_gb_h(ispin)%r_coef, &
1008  1, nsotot, 1, nsotot)
1009  CALL reallocate(rho_atom_set(iatom)%ga_Vlocal_gb_s(ispin)%r_coef, &
1010  1, nsotot, 1, nsotot)
1011 
1012  rho_atom_set(iatom)%ga_Vlocal_gb_h(ispin)%r_coef = 0.0_dp
1013  rho_atom_set(iatom)%ga_Vlocal_gb_s(ispin)%r_coef = 0.0_dp
1014  END DO
1015  END IF
1016 
1017  END DO ! iat
1018 
1019  END DO
1020 
1021  CALL timestop(handle)
1022 
1023  END SUBROUTINE allocate_rho_atom_internals
1024 
1025 ! **************************************************************************************************
1026 !> \brief ...
1027 !> \param rho_atom_set ...
1028 !> \param iatom ...
1029 !> \param ispin ...
1030 !> \param nr ...
1031 !> \param max_iso_not0 ...
1032 ! **************************************************************************************************
1033  SUBROUTINE allocate_rho_atom_rad(rho_atom_set, iatom, ispin, nr, max_iso_not0)
1034 
1035  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
1036  INTEGER, INTENT(IN) :: iatom, ispin, nr, max_iso_not0
1037 
1038  CHARACTER(len=*), PARAMETER :: routinen = 'allocate_rho_atom_rad'
1039 
1040  INTEGER :: handle, j
1041 
1042  CALL timeset(routinen, handle)
1043 
1044  ALLOCATE (rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef(1:nr, 1:max_iso_not0), &
1045  rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef(1:nr, 1:max_iso_not0), &
1046  rho_atom_set(iatom)%vrho_rad_h(ispin)%r_coef(1:nr, 1:max_iso_not0), &
1047  rho_atom_set(iatom)%vrho_rad_s(ispin)%r_coef(1:nr, 1:max_iso_not0))
1048 
1049  rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef = 0.0_dp
1050  rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef = 0.0_dp
1051  rho_atom_set(iatom)%vrho_rad_h(ispin)%r_coef = 0.0_dp
1052  rho_atom_set(iatom)%vrho_rad_s(ispin)%r_coef = 0.0_dp
1053 
1054  ALLOCATE (rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef(nr, max_iso_not0), &
1055  rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef(nr, max_iso_not0))
1056  rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef = 0.0_dp
1057  rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef = 0.0_dp
1058 
1059  DO j = 1, 3
1060  ALLOCATE (rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef(nr, max_iso_not0), &
1061  rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef(nr, max_iso_not0))
1062  rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef = 0.0_dp
1063  rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef = 0.0_dp
1064  END DO
1065 
1066  CALL timestop(handle)
1067 
1068  END SUBROUTINE allocate_rho_atom_rad
1069 
1070 ! **************************************************************************************************
1071 !> \brief ...
1072 !> \param rho_atom_set ...
1073 !> \param iatom ...
1074 !> \param ispin ...
1075 ! **************************************************************************************************
1076  SUBROUTINE set2zero_rho_atom_rad(rho_atom_set, iatom, ispin)
1077 
1078  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
1079  INTEGER, INTENT(IN) :: iatom, ispin
1080 
1081  INTEGER :: j
1082 
1083  rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef = 0.0_dp
1084  rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef = 0.0_dp
1085 
1086  rho_atom_set(iatom)%vrho_rad_h(ispin)%r_coef = 0.0_dp
1087  rho_atom_set(iatom)%vrho_rad_s(ispin)%r_coef = 0.0_dp
1088 
1089  rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef = 0.0_dp
1090  rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef = 0.0_dp
1091 
1092  DO j = 1, 3
1093  rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef = 0.0_dp
1094  rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef = 0.0_dp
1095  END DO
1096 
1097  END SUBROUTINE set2zero_rho_atom_rad
1098 
1099 END MODULE qs_rho_atom_methods
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
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.
Definition: atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
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.
Definition: atomic_kind_types.F:138
basis_set_types
Definition: basis_set_types.F:15
basis_set_types::get_gto_basis_set
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)
...
Definition: basis_set_types.F:615
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
kpoint_types::get_kpoint_info
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)
Retrieve information from a kpoint environment.
Definition: kpoint_types.F:333
lebedev
Generation of the spherical Lebedev grids. All Lebedev grids were generated with a precision of at le...
Definition: lebedev.F:57
lebedev::deallocate_lebedev_grids
subroutine, public deallocate_lebedev_grids()
...
Definition: lebedev.F:324
lebedev::lebedev_grid
type(oh_grid), dimension(nlg), target, public lebedev_grid
Definition: lebedev.F:85
lebedev::get_number_of_lebedev_grid
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
lebedev::init_lebedev_grids
subroutine, public init_lebedev_grids()
Load the coordinates and weights of the nonredundant Lebedev grid points.
Definition: lebedev.F:344
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition: mathconstants.F:113
memory_utilities
Utility routines for the memory handling.
Definition: memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::nsoset
integer, dimension(:), allocatable, public nsoset
Definition: orbital_pointers.F:42
orbital_pointers::indso
integer, dimension(:, :), allocatable, public indso
Definition: orbital_pointers.F:43
paw_basis_types
Definition: paw_basis_types.F:13
paw_basis_types::get_paw_basis_info
subroutine, public get_paw_basis_info(basis_1c, o2nindex, n2oindex, nsatbas)
Return some info on the PAW basis derived from a GTO basis set.
Definition: paw_basis_types.F:40
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
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, 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_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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, 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, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_grid_atom
Definition: qs_grid_atom.F:8
qs_grid_atom::create_grid_atom
subroutine, public create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
...
Definition: qs_grid_atom.F:183
qs_harmonics_atom
Definition: qs_harmonics_atom.F:8
qs_harmonics_atom::get_maxl_cg
subroutine, public get_maxl_cg(harmonics, orb_basis, llmax, max_s_harm)
...
Definition: qs_harmonics_atom.F:403
qs_harmonics_atom::create_harmonics_atom
subroutine, public create_harmonics_atom(harmonics, my_CG, na, llmax, maxs, max_s_harm, ll, wa, azi, pol)
...
Definition: qs_harmonics_atom.F:158
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, 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_r3d_rs_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, 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, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_kind_types::get_qs_kind_set
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, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Definition: qs_kind_types.F:864
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition: qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition: qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition: qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition: qs_neighbor_list_types.F:549
qs_oce_methods
Routines for the construction of the coefficients for the expansion of the atomic densities rho1_hard...
Definition: qs_oce_methods.F:16
qs_oce_methods::proj_blk
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)
Project a matrix block onto the local atomic functions.
Definition: qs_oce_methods.F:784
qs_oce_types
Definition: qs_oce_types.F:8
qs_rho_atom_methods
Definition: qs_rho_atom_methods.F:7
qs_rho_atom_methods::calculate_rho_atom
subroutine, public calculate_rho_atom(para_env, rho_atom_set, qs_kind, atom_list, natom, nspins, tot_rho1_h, tot_rho1_s)
...
Definition: qs_rho_atom_methods.F:103
qs_rho_atom_methods::allocate_rho_atom_internals
subroutine, public allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
...
Definition: qs_rho_atom_methods.F:916
qs_rho_atom_methods::init_rho_atom
subroutine, public init_rho_atom(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
...
Definition: qs_rho_atom_methods.F:754
qs_rho_atom_methods::calculate_rho_atom_coeff
subroutine, public calculate_rho_atom_coeff(qs_env, rho_ao, rho_atom_set, qs_kind_set, oce, sab, para_env)
...
Definition: qs_rho_atom_methods.F:420
qs_rho_atom_types
Definition: qs_rho_atom_types.F:8
qs_rho_atom_types::deallocate_rho_atom_set
subroutine, public deallocate_rho_atom_set(rho_atom_set)
...
Definition: qs_rho_atom_types.F:57
sap_kind_types
General overlap type integrals containers.
Definition: sap_kind_types.F:12
sap_kind_types::alist_pre_align_blk
subroutine, public alist_pre_align_blk(blk_in, ldin, blk_out, ldout, ilist, in, jlist, jn)
...
Definition: sap_kind_types.F:172
sap_kind_types::get_alist
subroutine, public get_alist(sap_int, alist, atom)
...
Definition: sap_kind_types.F:140
spherical_harmonics
Calculate spherical harmonics.
Definition: spherical_harmonics.F:27
spherical_harmonics::clebsch_gordon_init
subroutine, public clebsch_gordon_init(l)
...
Definition: spherical_harmonics.F:238
spherical_harmonics::clebsch_gordon_deallocate
subroutine, public clebsch_gordon_deallocate()
...
Definition: spherical_harmonics.F:220
util
All kind of helpful little routines.
Definition: util.F:14
util::get_limit
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
whittaker
Calculates special integrals.
Definition: whittaker.F:12
whittaker::whittaker_c0a
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
whittaker::whittaker_ci
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