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