Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : 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 cp_dbcsr_api, ONLY: dbcsr_get_block_p,&
18 : dbcsr_p_type
19 : USE gaunt_coefficients, ONLY: real_gaunt_coefficients
20 : USE kinds, ONLY: dp
21 : USE kpoint_types, ONLY: get_kpoint_info,&
22 : kpoint_type
23 : USE lebedev, ONLY: deallocate_lebedev_grids,&
24 : get_number_of_lebedev_grid,&
25 : init_lebedev_grids,&
26 : lebedev_grid
27 : USE mathconstants, ONLY: fourpi,&
28 : pi
29 : USE memory_utilities, ONLY: reallocate
30 : USE message_passing, ONLY: mp_para_env_type
31 : USE orbital_pointers, ONLY: indso,&
32 : nsoset
33 : USE paw_basis_types, ONLY: get_paw_basis_info
34 : USE qs_environment_types, ONLY: get_qs_env,&
35 : qs_environment_type
36 : USE qs_grid_atom, ONLY: create_grid_atom,&
37 : grid_atom_type
38 : USE qs_harmonics_atom, ONLY: create_harmonics_atom,&
39 : get_maxl_CG,&
40 : get_none0_cg_list,&
41 : harmonics_atom_type
42 : USE qs_kind_types, ONLY: get_qs_kind,&
43 : get_qs_kind_set,&
44 : qs_kind_type
45 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
46 : neighbor_list_iterate,&
47 : neighbor_list_iterator_create,&
48 : neighbor_list_iterator_p_type,&
49 : neighbor_list_iterator_release,&
50 : neighbor_list_set_p_type
51 : USE qs_oce_methods, ONLY: proj_blk
52 : USE qs_oce_types, ONLY: oce_matrix_type
53 : USE qs_rho_atom_types, ONLY: deallocate_rho_atom_set,&
54 : rho_atom_coeff,&
55 : rho_atom_type
56 : USE sap_kind_types, ONLY: alist_pre_align_blk,&
57 : alist_type,&
58 : clist_type,&
59 : get_alist
60 : USE util, ONLY: get_limit
61 : USE whittaker, ONLY: whittaker_c0a,&
62 : whittaker_ci
63 :
64 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, &
65 : !$ omp_get_thread_num, &
66 : !$ omp_lock_kind, &
67 : !$ omp_init_lock, omp_set_lock, &
68 : !$ omp_unset_lock, omp_destroy_lock
69 :
70 : #include "./base/base_uses.f90"
71 :
72 : IMPLICIT NONE
73 :
74 : PRIVATE
75 :
76 : ! *** Global parameters (only in this module)
77 :
78 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho_atom_methods'
79 :
80 : ! *** Public subroutines ***
81 :
82 : PUBLIC :: allocate_rho_atom_internals, &
83 : calculate_rho_atom, &
84 : calculate_rho_atom_coeff, &
85 : init_rho_atom, &
86 : replicate_rho_atom_radial
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 : !> \param rho1_h_spin ...
101 : !> \param rho1_s_spin ...
102 : !> \param rho1_h_aspin ...
103 : !> \param rho1_s_aspin ...
104 : ! **************************************************************************************************
105 79648 : SUBROUTINE calculate_rho_atom(para_env, rho_atom_set, qs_kind, atom_list, &
106 79648 : natom, nspins, tot_rho1_h, tot_rho1_s, &
107 : rho1_h_spin, rho1_s_spin, rho1_h_aspin, rho1_s_aspin)
108 :
109 : TYPE(mp_para_env_type), POINTER :: para_env
110 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
111 : TYPE(qs_kind_type), INTENT(IN) :: qs_kind
112 : INTEGER, DIMENSION(:), INTENT(IN) :: atom_list
113 : INTEGER, INTENT(IN) :: natom, nspins
114 : REAL(dp), DIMENSION(:), INTENT(INOUT) :: tot_rho1_h, tot_rho1_s
115 : REAL(dp), INTENT(INOUT) :: rho1_h_spin, rho1_s_spin, rho1_h_aspin, &
116 : rho1_s_aspin
117 :
118 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_atom'
119 :
120 : INTEGER :: damax_iso_not0_local, handle, i, i1, i2, iat, iatom, icg, ipgf1, ipgf2, ir, &
121 : iset1, iset2, iso, iso1, iso1_coeff, iso1_first, iso1_last, iso2, iso2_coeff, iso2_first, &
122 : iso2_last, j, l, l_iso, l_sub, l_sum, lmax12, lmax_expansion, lmin12, m1s, m2s, &
123 : max_iso_not0, max_iso_not0_local, max_npgf, max_s_harm, maxl, maxso, mepos, n1s, n2s, na, &
124 : nr, nset, num_pe, size1, size2
125 79648 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list, dacg_n_list
126 79648 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list, dacg_list
127 : INTEGER, DIMENSION(2) :: bo
128 79648 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, o2nindex
129 79648 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: done_vgg
130 : REAL(dp) :: c1, c2, cpc_h, cpc_s, rfun, rho_h, &
131 : rho_s, root_zet12, zet12
132 79648 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: erf_zet12, g1, g2, gg0, int1, int2
133 79648 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dgg, gg, gg_lm1, sfun_h, sfun_s
134 79648 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: g_rad, vgg
135 79648 : REAL(dp), DIMENSION(:, :), POINTER :: coeff_h, coeff_s, zet
136 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
137 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_CG_dxyz
138 : TYPE(grid_atom_type), POINTER :: grid_atom
139 : TYPE(gto_basis_set_type), POINTER :: basis_1c
140 : TYPE(harmonics_atom_type), POINTER :: harmonics
141 :
142 79648 : CALL timeset(routineN, handle)
143 :
144 : !Note: tau is taken care of separately in qs_vxc_atom.F
145 :
146 79648 : NULLIFY (basis_1c)
147 79648 : NULLIFY (harmonics, grid_atom)
148 79648 : NULLIFY (lmin, lmax, npgf, zet, my_CG, my_CG_dxyz, coeff_h, coeff_s)
149 :
150 79648 : CALL get_qs_kind(qs_kind, grid_atom=grid_atom, harmonics=harmonics)
151 79648 : CALL get_qs_kind(qs_kind, basis_set=basis_1c, basis_type="GAPW_1C")
152 :
153 : CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
154 : maxl=maxl, npgf=npgf, nset=nset, zet=zet, &
155 79648 : maxso=maxso)
156 :
157 79648 : CALL get_paw_basis_info(basis_1c, o2nindex=o2nindex)
158 :
159 79648 : max_iso_not0 = harmonics%max_iso_not0
160 79648 : max_s_harm = harmonics%max_s_harm
161 :
162 79648 : nr = grid_atom%nr
163 289410 : max_npgf = MAXVAL(npgf(1:nset))
164 79648 : lmax_expansion = indso(1, max_iso_not0)
165 : ! Distribute the atoms of this kind
166 79648 : num_pe = para_env%num_pe
167 79648 : mepos = para_env%mepos
168 79648 : bo = get_limit(natom, num_pe, mepos)
169 :
170 79648 : my_CG => harmonics%my_CG
171 79648 : my_CG_dxyz => harmonics%my_CG_dxyz
172 :
173 955776 : ALLOCATE (g1(nr), g2(nr), gg0(nr), gg(nr, 0:2*maxl), dgg(nr, 0:2*maxl), gg_lm1(nr, 0:2*maxl))
174 477888 : ALLOCATE (erf_zet12(nr), vgg(nr, 0:2*maxl, 0:indso(1, max_iso_not0)))
175 318592 : ALLOCATE (done_vgg(0:2*maxl, 0:indso(1, max_iso_not0)))
176 238944 : ALLOCATE (int1(nr), int2(nr))
177 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm), &
178 716832 : dacg_list(2, nsoset(maxl)**2, max_s_harm), dacg_n_list(max_s_harm))
179 398240 : ALLOCATE (g_rad(nr, max_npgf, nset))
180 :
181 289410 : DO iset1 = 1, nset
182 823954 : DO ipgf1 = 1, npgf(iset1)
183 29119362 : g_rad(1:nr, ipgf1, iset1) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
184 : END DO
185 : END DO
186 :
187 140053 : DO iat = bo(1), bo(2)
188 60405 : iatom = atom_list(iat)
189 207784 : DO i = 1, nspins
190 128136 : IF (.NOT. ASSOCIATED(rho_atom_set(iatom)%rho_rad_h(i)%r_coef)) THEN
191 8283 : CALL allocate_rho_atom_rad(rho_atom_set, iatom, i, nr, max_iso_not0)
192 : ELSE
193 59448 : CALL set2zero_rho_atom_rad(rho_atom_set, iatom, i)
194 : END IF
195 : END DO
196 : END DO
197 :
198 : m1s = 0
199 289410 : DO iset1 = 1, nset
200 : m2s = 0
201 949916 : DO iset2 = 1, nset
202 :
203 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
204 740154 : max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
205 740154 : CPASSERT(max_iso_not0_local <= max_iso_not0)
206 : CALL get_none0_cg_list(my_CG_dxyz, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
207 740154 : max_s_harm, lmax_expansion, dacg_list, dacg_n_list, damax_iso_not0_local)
208 740154 : n1s = nsoset(lmax(iset1))
209 :
210 2465860 : DO ipgf1 = 1, npgf(iset1)
211 1725706 : iso1_first = nsoset(lmin(iset1) - 1) + 1 + n1s*(ipgf1 - 1) + m1s
212 1725706 : iso1_last = nsoset(lmax(iset1)) + n1s*(ipgf1 - 1) + m1s
213 1725706 : size1 = iso1_last - iso1_first + 1
214 1725706 : iso1_first = o2nindex(iso1_first)
215 1725706 : iso1_last = o2nindex(iso1_last)
216 1725706 : i1 = iso1_last - iso1_first + 1
217 1725706 : CPASSERT(size1 == i1)
218 1725706 : i1 = nsoset(lmin(iset1) - 1) + 1
219 :
220 92709062 : g1(1:nr) = g_rad(1:nr, ipgf1, iset1)
221 :
222 1725706 : n2s = nsoset(lmax(iset2))
223 7331196 : DO ipgf2 = 1, npgf(iset2)
224 4865336 : iso2_first = nsoset(lmin(iset2) - 1) + 1 + n2s*(ipgf2 - 1) + m2s
225 4865336 : iso2_last = nsoset(lmax(iset2)) + n2s*(ipgf2 - 1) + m2s
226 4865336 : size2 = iso2_last - iso2_first + 1
227 4865336 : iso2_first = o2nindex(iso2_first)
228 4865336 : iso2_last = o2nindex(iso2_last)
229 4865336 : i2 = iso2_last - iso2_first + 1
230 4865336 : CPASSERT(size2 == i2)
231 4865336 : i2 = nsoset(lmin(iset2) - 1) + 1
232 :
233 262175360 : g2(1:nr) = g_rad(1:nr, ipgf2, iset2)
234 4865336 : lmin12 = lmin(iset1) + lmin(iset2)
235 4865336 : lmax12 = lmax(iset1) + lmax(iset2)
236 :
237 4865336 : zet12 = zet(ipgf1, iset1) + zet(ipgf2, iset2)
238 4865336 : root_zet12 = SQRT(zet(ipgf1, iset1) + zet(ipgf2, iset2))
239 262175360 : DO ir = 1, nr
240 262175360 : erf_zet12(ir) = erf(root_zet12*grid_atom%rad(ir))
241 : END DO
242 :
243 4865336 : gg = 0.0_dp
244 4865336 : dgg = 0.0_dp
245 4865336 : gg_lm1 = 0.0_dp
246 4865336 : vgg = 0.0_dp
247 4865336 : done_vgg = .FALSE.
248 : ! reduce the number of terms in the expansion local densities
249 4865336 : IF (lmin12 <= lmax_expansion) THEN
250 4862846 : IF (lmin12 == 0) THEN
251 140484540 : gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
252 140484540 : gg_lm1(1:nr, lmin12) = 0.0_dp
253 140484540 : gg0(1:nr) = gg(1:nr, lmin12)
254 : ELSE
255 121563830 : gg0(1:nr) = g1(1:nr)*g2(1:nr)
256 121563830 : gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
257 121563830 : gg_lm1(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12 - 1)*g1(1:nr)*g2(1:nr)
258 : END IF
259 :
260 : ! reduce the number of terms in the expansion local densities
261 4862846 : IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
262 :
263 7627010 : DO l = lmin12 + 1, lmax12
264 152150324 : gg(1:nr, l) = grid_atom%rad(1:nr)*gg(1:nr, l - 1)
265 152150324 : gg_lm1(1:nr, l) = gg(1:nr, l - 1)
266 157013170 : dgg(1:nr, l - 1) = -2.0_dp*(zet(ipgf1, iset1) + zet(ipgf2, iset2))*gg(1:nr, l)
267 :
268 : END DO
269 : dgg(1:nr, lmax12) = -2.0_dp*(zet(ipgf1, iset1) + &
270 262048370 : zet(ipgf2, iset2))*grid_atom%rad(1:nr)*gg(1:nr, lmax12)
271 :
272 4862846 : c2 = SQRT(pi*pi*pi/(zet12*zet12*zet12))
273 :
274 41240604 : DO iso = 1, max_iso_not0_local
275 36377758 : l_iso = indso(1, iso)
276 36377758 : c1 = fourpi/(2._dp*REAL(l_iso, dp) + 1._dp)
277 120956394 : DO icg = 1, cg_n_list(iso)
278 79715790 : iso1 = cg_list(1, icg, iso)
279 79715790 : iso2 = cg_list(2, icg, iso)
280 :
281 79715790 : l = indso(1, iso1) + indso(1, iso2)
282 79715790 : CPASSERT(l <= lmax_expansion)
283 79715790 : IF (done_vgg(l, l_iso)) CYCLE
284 10023112 : L_sum = l + l_iso
285 10023112 : L_sub = l - l_iso
286 :
287 10023112 : IF (l_sum == 0) THEN
288 140484540 : vgg(1:nr, l, l_iso) = erf_zet12(1:nr)*grid_atom%oorad2l(1:nr, 1)*c2
289 : ELSE
290 7504396 : CALL whittaker_c0a(int1, grid_atom%rad, gg0, erf_zet12, zet12, l, l_iso, nr)
291 7504396 : CALL whittaker_ci(int2, grid_atom%rad, gg0, zet12, L_sub, nr)
292 :
293 399995496 : DO ir = 1, nr
294 392491100 : int2(ir) = grid_atom%rad2l(ir, l_iso)*int2(ir)
295 399995496 : vgg(ir, l, l_iso) = c1*(int1(ir) + int2(ir))
296 : END DO
297 : END IF
298 116093548 : done_vgg(l, l_iso) = .TRUE.
299 : END DO
300 : END DO
301 : END IF ! lmax_expansion
302 :
303 9868063 : DO iat = bo(1), bo(2)
304 3277021 : iatom = atom_list(iat)
305 :
306 11976105 : DO i = 1, nspins
307 3833748 : coeff_h => rho_atom_set(iatom)%cpc_h(i)%r_coef
308 3833748 : coeff_s => rho_atom_set(iatom)%cpc_s(i)%r_coef
309 :
310 28700385 : DO iso = 1, max_iso_not0_local
311 24866637 : l_iso = indso(1, iso)
312 80876796 : DO icg = 1, cg_n_list(iso)
313 52176411 : iso1 = cg_list(1, icg, iso)
314 52176411 : iso2 = cg_list(2, icg, iso)
315 :
316 52176411 : l = indso(1, iso1) + indso(1, iso2)
317 52176411 : CPASSERT(l <= lmax_expansion)
318 52176411 : iso1_coeff = iso1_first + iso1 - i1
319 52176411 : iso2_coeff = iso2_first + iso2 - i2
320 52176411 : cpc_h = coeff_h(iso1_coeff, iso2_coeff)*my_CG(iso1, iso2, iso)
321 52176411 : cpc_s = coeff_s(iso1_coeff, iso2_coeff)*my_CG(iso1, iso2, iso)
322 :
323 : rho_atom_set(iatom)%rho_rad_h(i)%r_coef(1:nr, iso) = &
324 : rho_atom_set(iatom)%rho_rad_h(i)%r_coef(1:nr, iso) + &
325 2780796595 : gg(1:nr, l)*cpc_h
326 :
327 : rho_atom_set(iatom)%rho_rad_s(i)%r_coef(1:nr, iso) = &
328 : rho_atom_set(iatom)%rho_rad_s(i)%r_coef(1:nr, iso) + &
329 2780796595 : gg(1:nr, l)*cpc_s
330 :
331 : rho_atom_set(iatom)%drho_rad_h(i)%r_coef(1:nr, iso) = &
332 : rho_atom_set(iatom)%drho_rad_h(i)%r_coef(1:nr, iso) + &
333 2780796595 : dgg(1:nr, l)*cpc_h
334 :
335 : rho_atom_set(iatom)%drho_rad_s(i)%r_coef(1:nr, iso) = &
336 : rho_atom_set(iatom)%drho_rad_s(i)%r_coef(1:nr, iso) + &
337 2780796595 : dgg(1:nr, l)*cpc_s
338 :
339 : rho_atom_set(iatom)%vrho_rad_h(i)%r_coef(1:nr, iso) = &
340 : rho_atom_set(iatom)%vrho_rad_h(i)%r_coef(1:nr, iso) + &
341 2780796595 : vgg(1:nr, l, l_iso)*cpc_h
342 :
343 : rho_atom_set(iatom)%vrho_rad_s(i)%r_coef(1:nr, iso) = &
344 : rho_atom_set(iatom)%vrho_rad_s(i)%r_coef(1:nr, iso) + &
345 2805663232 : vgg(1:nr, l, l_iso)*cpc_s
346 :
347 : END DO ! icg
348 :
349 : END DO ! iso
350 :
351 85193861 : DO iso = 1, max_iso_not0 !damax_iso_not0_local
352 78083092 : l_iso = indso(1, iso)
353 167455703 : DO icg = 1, dacg_n_list(iso)
354 85538863 : iso1 = dacg_list(1, icg, iso)
355 85538863 : iso2 = dacg_list(2, icg, iso)
356 85538863 : l = indso(1, iso1) + indso(1, iso2)
357 85538863 : CPASSERT(l <= lmax_expansion)
358 85538863 : iso1_coeff = iso1_first + iso1 - i1
359 85538863 : iso2_coeff = iso2_first + iso2 - i2
360 85538863 : cpc_h = coeff_h(iso1_coeff, iso2_coeff)
361 85538863 : cpc_s = coeff_s(iso1_coeff, iso2_coeff)
362 420238544 : DO j = 1, 3
363 : rho_atom_set(iatom)%rho_rad_h_d(j, i)%r_coef(1:nr, iso) = &
364 : rho_atom_set(iatom)%rho_rad_h_d(j, i)%r_coef(1:nr, iso) + &
365 13534290849 : gg_lm1(1:nr, l)*cpc_h*my_CG_dxyz(j, iso1, iso2, iso)
366 :
367 : rho_atom_set(iatom)%rho_rad_s_d(j, i)%r_coef(1:nr, iso) = &
368 : rho_atom_set(iatom)%rho_rad_s_d(j, i)%r_coef(1:nr, iso) + &
369 13619829712 : gg_lm1(1:nr, l)*cpc_s*my_CG_dxyz(j, iso1, iso2, iso)
370 : END DO
371 : END DO ! icg
372 :
373 : END DO ! iso
374 :
375 : END DO ! i
376 : END DO ! iat
377 :
378 : END DO ! ipgf2
379 : END DO ! ipgf1
380 2430224 : m2s = m2s + maxso
381 : END DO ! iset2
382 289410 : m1s = m1s + maxso
383 : END DO ! iset1
384 :
385 140053 : DO iat = bo(1), bo(2)
386 60405 : iatom = atom_list(iat)
387 :
388 128136 : DO i = 1, nspins
389 :
390 1049451 : DO iso = 1, max_iso_not0
391 : rho_s = 0.0_dp
392 : rho_h = 0.0_dp
393 49948029 : DO ir = 1, nr
394 49026714 : rho_h = rho_h + rho_atom_set(iatom)%rho_rad_h(i)%r_coef(ir, iso)*grid_atom%wr(ir)
395 49948029 : rho_s = rho_s + rho_atom_set(iatom)%rho_rad_s(i)%r_coef(ir, iso)*grid_atom%wr(ir)
396 : END DO ! ir
397 921315 : tot_rho1_h(i) = tot_rho1_h(i) + rho_h*harmonics%slm_int(iso)
398 989046 : tot_rho1_s(i) = tot_rho1_s(i) + rho_s*harmonics%slm_int(iso)
399 : END DO ! iso
400 :
401 : END DO ! ispin
402 :
403 140053 : IF (nspins == 2) THEN
404 7326 : na = SIZE(harmonics%slm, 1)
405 43956 : ALLOCATE (sfun_h(nr, na), sfun_s(nr, na))
406 7326 : sfun_h = 0.0_dp
407 7326 : sfun_s = 0.0_dp
408 109612 : DO iso = 1, max_iso_not0
409 6493682 : DO ir = 1, nr
410 : rfun = grid_atom%wr(ir)*(rho_atom_set(iatom)%rho_rad_h(1)%r_coef(ir, iso) - &
411 6384070 : rho_atom_set(iatom)%rho_rad_h(2)%r_coef(ir, iso))
412 356469570 : sfun_h(ir, 1:na) = sfun_h(ir, 1:na) + rfun*harmonics%slm(1:na, iso)*grid_atom%wa(1:na)
413 : rfun = grid_atom%wr(ir)*(rho_atom_set(iatom)%rho_rad_s(1)%r_coef(ir, iso) - &
414 6384070 : rho_atom_set(iatom)%rho_rad_s(2)%r_coef(ir, iso))
415 356571856 : sfun_s(ir, 1:na) = sfun_s(ir, 1:na) + rfun*harmonics%slm(1:na, iso)*grid_atom%wa(1:na)
416 : END DO
417 : END DO
418 23624802 : rho1_h_spin = rho1_h_spin + SUM(sfun_h(1:nr, 1:na))
419 23624802 : rho1_s_spin = rho1_s_spin + SUM(sfun_s(1:nr, 1:na))
420 23624802 : rho1_h_aspin = rho1_h_aspin + SUM(ABS(sfun_h(1:nr, 1:na)))
421 23624802 : rho1_s_aspin = rho1_s_aspin + SUM(ABS(sfun_s(1:nr, 1:na)))
422 7326 : DEALLOCATE (sfun_h, sfun_s)
423 : END IF
424 :
425 : END DO ! iat
426 :
427 79648 : DEALLOCATE (g1, g2, gg0, gg, gg_lm1, dgg, vgg, done_vgg, erf_zet12, int1, int2, g_rad)
428 79648 : DEALLOCATE (cg_list, cg_n_list, dacg_list, dacg_n_list)
429 79648 : DEALLOCATE (o2nindex)
430 :
431 79648 : CALL timestop(handle)
432 :
433 238944 : END SUBROUTINE calculate_rho_atom
434 :
435 : ! **************************************************************************************************
436 : !> \brief Replicate the radial hard/soft density data needed to evaluate one-center tails on
437 : !> rank-local target grids. The compact one-center density matrices are already global;
438 : !> this routine performs one packed reduction for the derived radial fields of a kind.
439 : !> \param para_env ...
440 : !> \param rho_atom_set ...
441 : !> \param qs_kind ...
442 : !> \param atom_list ...
443 : !> \param natom ...
444 : !> \param nspins ...
445 : ! **************************************************************************************************
446 160 : SUBROUTINE replicate_rho_atom_radial(para_env, rho_atom_set, qs_kind, atom_list, natom, nspins)
447 : TYPE(mp_para_env_type), POINTER :: para_env
448 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
449 : TYPE(qs_kind_type), INTENT(IN) :: qs_kind
450 : INTEGER, DIMENSION(:), INTENT(IN) :: atom_list
451 : INTEGER, INTENT(IN) :: natom, nspins
452 :
453 : CHARACTER(len=*), PARAMETER :: routineN = 'replicate_rho_atom_radial'
454 :
455 : INTEGER :: block_size, bo(2), cursor, handle, iat, &
456 : iatom, ispin, j, max_iso_not0, ncoeff, &
457 : nr
458 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: buffer
459 : TYPE(grid_atom_type), POINTER :: grid_atom
460 : TYPE(harmonics_atom_type), POINTER :: harmonics
461 :
462 160 : CALL timeset(routineN, handle)
463 :
464 160 : NULLIFY (grid_atom, harmonics)
465 160 : CALL get_qs_kind(qs_kind, grid_atom=grid_atom, harmonics=harmonics)
466 160 : CPASSERT(ASSOCIATED(grid_atom))
467 160 : CPASSERT(ASSOCIATED(harmonics))
468 160 : nr = grid_atom%nr
469 160 : max_iso_not0 = harmonics%max_iso_not0
470 160 : ncoeff = nr*max_iso_not0
471 160 : block_size = 10*nspins*ncoeff
472 480 : ALLOCATE (buffer(block_size*natom))
473 160 : buffer = 0.0_dp
474 :
475 160 : bo = get_limit(natom, para_env%num_pe, para_env%mepos)
476 286 : DO iat = bo(1), bo(2)
477 126 : iatom = atom_list(iat)
478 126 : cursor = (iat - 1)*block_size
479 252 : DO ispin = 1, nspins
480 : buffer(cursor + 1:cursor + ncoeff) = &
481 252 : RESHAPE(rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef, [ncoeff])
482 126 : cursor = cursor + ncoeff
483 : buffer(cursor + 1:cursor + ncoeff) = &
484 252 : RESHAPE(rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef, [ncoeff])
485 126 : cursor = cursor + ncoeff
486 : buffer(cursor + 1:cursor + ncoeff) = &
487 252 : RESHAPE(rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef, [ncoeff])
488 126 : cursor = cursor + ncoeff
489 : buffer(cursor + 1:cursor + ncoeff) = &
490 252 : RESHAPE(rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef, [ncoeff])
491 126 : cursor = cursor + ncoeff
492 630 : DO j = 1, 3
493 : buffer(cursor + 1:cursor + ncoeff) = &
494 756 : RESHAPE(rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef, [ncoeff])
495 378 : cursor = cursor + ncoeff
496 : buffer(cursor + 1:cursor + ncoeff) = &
497 756 : RESHAPE(rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef, [ncoeff])
498 504 : cursor = cursor + ncoeff
499 : END DO
500 : END DO
501 286 : CPASSERT(cursor == iat*block_size)
502 : END DO
503 160 : CALL para_env%sum(buffer)
504 :
505 412 : DO iat = 1, natom
506 252 : iatom = atom_list(iat)
507 252 : cursor = (iat - 1)*block_size
508 504 : DO ispin = 1, nspins
509 252 : IF (.NOT. ASSOCIATED(rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef)) THEN
510 35 : CALL allocate_rho_atom_rad(rho_atom_set, iatom, ispin, nr, max_iso_not0)
511 : END IF
512 : rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef = &
513 33360 : RESHAPE(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
514 252 : cursor = cursor + ncoeff
515 : rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef = &
516 33360 : RESHAPE(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
517 252 : cursor = cursor + ncoeff
518 : rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef = &
519 33360 : RESHAPE(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
520 252 : cursor = cursor + ncoeff
521 : rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef = &
522 33360 : RESHAPE(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
523 252 : cursor = cursor + ncoeff
524 1260 : DO j = 1, 3
525 : rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef = &
526 100080 : RESHAPE(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
527 756 : cursor = cursor + ncoeff
528 : rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef = &
529 100080 : RESHAPE(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
530 1008 : cursor = cursor + ncoeff
531 : END DO
532 : END DO
533 412 : CPASSERT(cursor == iat*block_size)
534 : END DO
535 160 : DEALLOCATE (buffer)
536 :
537 160 : CALL timestop(handle)
538 :
539 160 : END SUBROUTINE replicate_rho_atom_radial
540 :
541 : ! **************************************************************************************************
542 : !> \brief ...
543 : !> \param qs_env QuickStep environment
544 : !> (accessed components: atomic_kind_set, dft_control%nimages,
545 : !> dft_control%nspins, kpoints%cell_to_index)
546 : !> \param rho_ao density matrix in atomic basis set
547 : !> \param rho_atom_set ...
548 : !> \param qs_kind_set list of QuickStep kinds
549 : !> \param oce one-centre expansion coefficients
550 : !> \param sab neighbour pair list
551 : !> \param para_env parallel environment
552 : !> \par History
553 : !> Add OpenMP [Apr 2016, EPCC]
554 : !> Use automatic arrays [Sep 2016, M Tucker]
555 : !> Allow for external non-default kind_set, oce and sab [Dec 2019, A Bussy]
556 : !> \note Consider to declare 'rho_ao' dummy argument as a pointer to the two-dimensional
557 : !> (1:nspins, 1:nimages) set of matrices.
558 : ! **************************************************************************************************
559 46238 : SUBROUTINE calculate_rho_atom_coeff(qs_env, rho_ao, rho_atom_set, qs_kind_set, oce, sab, para_env)
560 :
561 : TYPE(qs_environment_type), POINTER :: qs_env
562 : TYPE(dbcsr_p_type), DIMENSION(*) :: rho_ao
563 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
564 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
565 : TYPE(oce_matrix_type), POINTER :: oce
566 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
567 : POINTER :: sab
568 : TYPE(mp_para_env_type), POINTER :: para_env
569 :
570 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_atom_coeff'
571 :
572 : INTEGER :: bo(2), handle, i, iac, iatom, ibc, icol, ikind, img, irow, ispin, jatom, jkind, &
573 : kac, katom, kbc, kkind, len_CPC, len_PC1, max_gau, max_nsgf, mepos, n_cont_a, n_cont_b, &
574 : nat_kind, natom, nimages, nkind, nsoctot, nspins, num_pe
575 46238 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, nsatbas_kind
576 : INTEGER, DIMENSION(3) :: cell_b
577 46238 : INTEGER, DIMENSION(:), POINTER :: a_list, list_a, list_b
578 46238 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
579 : LOGICAL :: dista, distab, distb, found, paw_atom
580 46238 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: has_intac, paw_kind
581 46238 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: proj_work1, proj_work2
582 46238 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: p_matrix
583 : REAL(KIND=dp) :: eps_cpc, factor, pmax
584 : REAL(KIND=dp), DIMENSION(3) :: rab
585 46238 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: C_coeff_hh_a, C_coeff_hh_b, &
586 46238 : C_coeff_ss_a, C_coeff_ss_b, r_coef_h, &
587 46238 : r_coef_s
588 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
589 46238 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
590 : TYPE(clist_type), POINTER :: clist_a, clist_b
591 : TYPE(dft_control_type), POINTER :: dft_control
592 46238 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
593 : TYPE(gto_basis_set_type), POINTER :: basis_1c, basis_set_a, basis_set_b
594 : TYPE(kpoint_type), POINTER :: kpoints
595 : TYPE(neighbor_list_iterator_p_type), &
596 46238 : DIMENSION(:), POINTER :: nl_iterator
597 46238 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: p_block_spin
598 :
599 46238 : !$ INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: locks
600 : !$ INTEGER :: lock, number_of_locks
601 :
602 46238 : CALL timeset(routineN, handle)
603 :
604 : CALL get_qs_env(qs_env=qs_env, &
605 : dft_control=dft_control, &
606 46238 : atomic_kind_set=atomic_kind_set)
607 :
608 46238 : eps_cpc = dft_control%qs_control%gapw_control%eps_cpc
609 :
610 46238 : CPASSERT(ASSOCIATED(qs_kind_set))
611 46238 : CPASSERT(ASSOCIATED(rho_atom_set))
612 46238 : CPASSERT(ASSOCIATED(oce))
613 46238 : CPASSERT(ASSOCIATED(sab))
614 :
615 46238 : nspins = dft_control%nspins
616 46238 : nimages = dft_control%nimages
617 :
618 46238 : NULLIFY (cell_to_index)
619 46238 : IF (nimages > 1) THEN
620 680 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
621 680 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
622 : END IF
623 :
624 46238 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
625 46238 : CALL get_qs_kind_set(qs_kind_set, maxsgf=max_nsgf, maxgtops=max_gau, basis_type='GAPW_1C')
626 :
627 46238 : nkind = SIZE(atomic_kind_set)
628 : ! Initialize to 0 the CPC coefficients and the local density arrays
629 136858 : DO ikind = 1, nkind
630 90620 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=a_list, natom=nat_kind)
631 90620 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
632 :
633 90620 : IF (.NOT. paw_atom) CYCLE
634 207372 : DO i = 1, nat_kind
635 125484 : iatom = a_list(i)
636 348200 : DO ispin = 1, nspins
637 67625556 : rho_atom_set(iatom)%cpc_h(ispin)%r_coef = 0.0_dp
638 67751040 : rho_atom_set(iatom)%cpc_s(ispin)%r_coef = 0.0_dp
639 : END DO ! ispin
640 : END DO ! i
641 :
642 81888 : num_pe = para_env%num_pe
643 81888 : mepos = para_env%mepos
644 81888 : bo = get_limit(nat_kind, num_pe, mepos)
645 281488 : DO i = bo(1), bo(2)
646 62742 : iatom = a_list(i)
647 223776 : DO ispin = 1, nspins
648 193861446 : rho_atom_set(iatom)%ga_Vlocal_gb_h(ispin)%r_coef = 0.0_dp
649 193924188 : rho_atom_set(iatom)%ga_Vlocal_gb_s(ispin)%r_coef = 0.0_dp
650 : END DO ! ispin
651 : END DO ! i
652 : END DO ! ikind
653 :
654 229334 : ALLOCATE (basis_set_list(nkind))
655 277428 : ALLOCATE (paw_kind(nkind), nsatbas_kind(nkind), has_intac(nkind*nkind))
656 46238 : paw_kind(:) = .FALSE.
657 46238 : nsatbas_kind(:) = 0
658 46238 : has_intac(:) = .FALSE.
659 136858 : DO ikind = 1, nkind
660 90620 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a)
661 90620 : IF (ASSOCIATED(basis_set_a)) THEN
662 90620 : basis_set_list(ikind)%gto_basis_set => basis_set_a
663 : ELSE
664 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
665 : END IF
666 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C", &
667 90620 : paw_atom=paw_kind(ikind))
668 136858 : IF (paw_kind(ikind)) CALL get_paw_basis_info(basis_1c, nsatbas=nsatbas_kind(ikind))
669 : END DO
670 241798 : DO ikind = 1, nkind*nkind
671 241798 : has_intac(ikind) = ASSOCIATED(oce%intac(ikind)%alist)
672 : END DO
673 :
674 46238 : len_PC1 = max_nsgf*max_gau
675 46238 : len_CPC = max_gau*max_gau
676 :
677 : num_pe = 1
678 46238 : !$ num_pe = omp_get_max_threads()
679 46238 : CALL neighbor_list_iterator_create(nl_iterator, sab, nthread=num_pe)
680 :
681 : !$OMP PARALLEL DEFAULT( NONE ) &
682 : !$OMP SHARED( max_nsgf, max_gau &
683 : !$OMP , len_PC1, len_CPC &
684 : !$OMP , nl_iterator, basis_set_list &
685 : !$OMP , nimages, cell_to_index &
686 : !$OMP , nspins, rho_ao &
687 : !$OMP , nkind, qs_kind_set &
688 : !$OMP , oce, eps_cpc &
689 : !$OMP , rho_atom_set &
690 : !$OMP , natom, locks, number_of_locks &
691 : !$OMP , paw_kind, nsatbas_kind, has_intac &
692 : !$OMP ) &
693 : !$OMP PRIVATE( p_block_spin, ispin &
694 : !$OMP , p_matrix, proj_work1, proj_work2 &
695 : !$OMP , mepos &
696 : !$OMP , ikind, jkind, iatom, jatom &
697 : !$OMP , cell_b, rab &
698 : !$OMP , basis_set_a, basis_set_b &
699 : !$OMP , pmax, irow, icol, img &
700 : !$OMP , found &
701 : !$OMP , kkind &
702 : !$OMP , nsoctot, katom &
703 : !$OMP , iac , alist_ac, kac, n_cont_a, list_a &
704 : !$OMP , ibc , alist_bc, kbc, n_cont_b, list_b &
705 : !$OMP , clist_a, clist_b &
706 : !$OMP , C_coeff_hh_a, C_coeff_ss_a, dista &
707 : !$OMP , C_coeff_hh_b, C_coeff_ss_b, distb &
708 : !$OMP , distab &
709 : !$OMP , factor, r_coef_h, r_coef_s &
710 46238 : !$OMP )
711 :
712 : ALLOCATE (p_block_spin(nspins))
713 : ALLOCATE (p_matrix(max_nsgf, max_nsgf))
714 : ALLOCATE (proj_work1(len_PC1), proj_work2(len_CPC))
715 :
716 : !$OMP SINGLE
717 : !$ number_of_locks = nspins*natom
718 : !$ ALLOCATE (locks(number_of_locks))
719 : !$OMP END SINGLE
720 :
721 : !$OMP DO
722 : !$ DO lock = 1, number_of_locks
723 : !$ call omp_init_lock(locks(lock))
724 : !$ END DO
725 : !$OMP END DO
726 :
727 : mepos = 0
728 : !$ mepos = omp_get_thread_num()
729 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
730 :
731 : CALL get_iterator_info(nl_iterator, mepos=mepos, &
732 : ikind=ikind, jkind=jkind, &
733 : iatom=iatom, jatom=jatom, &
734 : cell=cell_b, r=rab)
735 :
736 : basis_set_a => basis_set_list(ikind)%gto_basis_set
737 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
738 : basis_set_b => basis_set_list(jkind)%gto_basis_set
739 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
740 :
741 : pmax = 0._dp
742 : IF (iatom <= jatom) THEN
743 : irow = iatom
744 : icol = jatom
745 : ELSE
746 : irow = jatom
747 : icol = iatom
748 : END IF
749 :
750 : IF (nimages > 1) THEN
751 : img = cell_to_index(cell_b(1), cell_b(2), cell_b(3))
752 : CPASSERT(img > 0)
753 : ELSE
754 : img = 1
755 : END IF
756 :
757 : DO ispin = 1, nspins
758 : CALL dbcsr_get_block_p(matrix=rho_ao(nspins*(img - 1) + ispin)%matrix, &
759 : row=irow, col=icol, BLOCK=p_block_spin(ispin)%r_coef, &
760 : found=found)
761 : pmax = pmax + MAXVAL(ABS(p_block_spin(ispin)%r_coef))
762 : END DO
763 :
764 : DO kkind = 1, nkind
765 : IF (.NOT. paw_kind(kkind)) CYCLE
766 :
767 : nsoctot = nsatbas_kind(kkind)
768 :
769 : iac = ikind + nkind*(kkind - 1)
770 : ibc = jkind + nkind*(kkind - 1)
771 : IF (.NOT. has_intac(iac)) CYCLE
772 : IF (.NOT. has_intac(ibc)) CYCLE
773 :
774 : CALL get_alist(oce%intac(iac), alist_ac, iatom)
775 : CALL get_alist(oce%intac(ibc), alist_bc, jatom)
776 : IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
777 : IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
778 :
779 : DO kac = 1, alist_ac%nclist
780 : DO kbc = 1, alist_bc%nclist
781 : IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
782 : IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
783 : IF (pmax*alist_bc%clist(kbc)%maxac*alist_ac%clist(kac)%maxac < eps_cpc) CYCLE
784 :
785 : n_cont_a = alist_ac%clist(kac)%nsgf_cnt
786 : n_cont_b = alist_bc%clist(kbc)%nsgf_cnt
787 : IF (n_cont_a == 0 .OR. n_cont_b == 0) CYCLE
788 :
789 : list_a => alist_ac%clist(kac)%sgf_list
790 : list_b => alist_bc%clist(kbc)%sgf_list
791 :
792 : katom = alist_ac%clist(kac)%catom
793 :
794 : ! Scalar aliases avoid nested component-slice bounds-check temporaries.
795 : clist_a => alist_ac%clist(kac)
796 : clist_b => alist_bc%clist(kbc)
797 : IF (iatom == katom .AND. ALL(alist_ac%clist(kac)%cell == 0)) THEN
798 : C_coeff_hh_a => clist_a%achint(:, :, 1)
799 : C_coeff_ss_a => clist_a%acint(:, :, 1)
800 : dista = .FALSE.
801 : ELSE
802 : C_coeff_hh_a => clist_a%acint(:, :, 1)
803 : C_coeff_ss_a => clist_a%acint(:, :, 1)
804 : dista = .TRUE.
805 : END IF
806 : IF (jatom == katom .AND. ALL(alist_bc%clist(kbc)%cell == 0)) THEN
807 : C_coeff_hh_b => clist_b%achint(:, :, 1)
808 : C_coeff_ss_b => clist_b%acint(:, :, 1)
809 : distb = .FALSE.
810 : ELSE
811 : C_coeff_hh_b => clist_b%acint(:, :, 1)
812 : C_coeff_ss_b => clist_b%acint(:, :, 1)
813 : distb = .TRUE.
814 : END IF
815 :
816 : distab = dista .AND. distb
817 :
818 : DO ispin = 1, nspins
819 :
820 : IF (iatom <= jatom) THEN
821 : CALL alist_pre_align_blk(p_block_spin(ispin)%r_coef, &
822 : SIZE(p_block_spin(ispin)%r_coef, 1), p_matrix, SIZE(p_matrix, 1), &
823 : list_a, n_cont_a, list_b, n_cont_b)
824 : ELSE
825 : CALL alist_pre_align_blk(p_block_spin(ispin)%r_coef, &
826 : SIZE(p_block_spin(ispin)%r_coef, 1), p_matrix, SIZE(p_matrix, 1), &
827 : list_b, n_cont_b, list_a, n_cont_a)
828 : END IF
829 :
830 : factor = 1.0_dp
831 : IF (iatom == jatom) factor = 0.5_dp
832 :
833 : r_coef_h => rho_atom_set(katom)%cpc_h(ispin)%r_coef
834 : r_coef_s => rho_atom_set(katom)%cpc_s(ispin)%r_coef
835 :
836 : !$ CALL omp_set_lock(locks((katom - 1)*nspins + ispin))
837 : IF (iatom <= jatom) THEN
838 : CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, &
839 : C_coeff_hh_b, C_coeff_ss_b, n_cont_b, &
840 : p_matrix, max_nsgf, r_coef_h, r_coef_s, nsoctot, &
841 : len_PC1, len_CPC, factor, distab, proj_work1, proj_work2)
842 : ELSE
843 : CALL proj_blk(C_coeff_hh_b, C_coeff_ss_b, n_cont_b, &
844 : C_coeff_hh_a, C_coeff_ss_a, n_cont_a, &
845 : p_matrix, max_nsgf, r_coef_h, r_coef_s, nsoctot, &
846 : len_PC1, len_CPC, factor, distab, proj_work1, proj_work2)
847 : END IF
848 : !$ CALL omp_unset_lock(locks((katom - 1)*nspins + ispin))
849 :
850 : END DO
851 : EXIT !search loop over jatom-katom list
852 : END IF
853 : END DO
854 : END DO
855 : END DO
856 : END DO
857 : ! Wait for all threads to finish the loop before locks can be freed
858 : !$OMP BARRIER
859 :
860 : !$OMP DO
861 : !$ DO lock = 1, number_of_locks
862 : !$ call omp_destroy_lock(locks(lock))
863 : !$ END DO
864 : !$OMP END DO
865 : !$OMP SINGLE
866 : !$ DEALLOCATE (locks)
867 : !$OMP END SINGLE NOWAIT
868 :
869 : DEALLOCATE (p_block_spin, p_matrix, proj_work1, proj_work2)
870 : !$OMP END PARALLEL
871 :
872 46238 : CALL neighbor_list_iterator_release(nl_iterator)
873 :
874 46238 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
875 :
876 190406 : DO iatom = 1, natom
877 305592 : ikind = kind_of(iatom)
878 :
879 351830 : DO ispin = 1, nspins
880 305592 : IF (ASSOCIATED(rho_atom_set(iatom)%cpc_h(ispin)%r_coef)) THEN
881 135110284 : CALL para_env%sum(rho_atom_set(iatom)%cpc_h(ispin)%r_coef)
882 135110284 : CALL para_env%sum(rho_atom_set(iatom)%cpc_s(ispin)%r_coef)
883 140828 : r_coef_h => rho_atom_set(iatom)%cpc_h(ispin)%r_coef
884 140828 : r_coef_s => rho_atom_set(iatom)%cpc_s(ispin)%r_coef
885 135251112 : r_coef_h(:, :) = r_coef_h(:, :) + TRANSPOSE(r_coef_h(:, :))
886 135251112 : r_coef_s(:, :) = r_coef_s(:, :) + TRANSPOSE(r_coef_s(:, :))
887 : END IF
888 : END DO
889 :
890 : END DO
891 :
892 46238 : DEALLOCATE (kind_of, basis_set_list, paw_kind, nsatbas_kind, has_intac)
893 :
894 46238 : CALL timestop(handle)
895 :
896 138714 : END SUBROUTINE calculate_rho_atom_coeff
897 :
898 : ! **************************************************************************************************
899 : !> \brief ...
900 : !> \param rho_atom_set the type to initialize
901 : !> \param atomic_kind_set list of atomic kinds
902 : !> \param qs_kind_set the kind set from which to take quantum numbers and basis info
903 : !> \param dft_control DFT control type
904 : !> \param para_env parallel environment
905 : !> \par History:
906 : !> - Generalised by providing the rho_atom_set and the qs_kind_set 12.2019 (A.Bussy)
907 : ! **************************************************************************************************
908 5088 : SUBROUTINE init_rho_atom(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
909 :
910 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
911 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
912 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
913 : TYPE(dft_control_type), POINTER :: dft_control
914 : TYPE(mp_para_env_type), POINTER :: para_env
915 :
916 : CHARACTER(len=*), PARAMETER :: routineN = 'init_rho_atom'
917 :
918 : CHARACTER(LEN=4) :: my_ikind, my_la, my_llmax
919 : INTEGER :: handle, ikind, la, ll, llmax, &
920 : lmax_sphere, max_s_harm, max_s_set, &
921 : maxl, maxlgto, maxs, na, nat, natom, &
922 : nr, nspins, quadrature
923 1696 : INTEGER, DIMENSION(:), POINTER :: atom_list
924 : LOGICAL :: paw_atom
925 1696 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
926 : TYPE(gapw_control_type), POINTER :: gapw_control
927 : TYPE(grid_atom_type), POINTER :: grid_atom
928 : TYPE(gto_basis_set_type), POINTER :: basis_1c_set
929 : TYPE(harmonics_atom_type), POINTER :: harmonics
930 :
931 1696 : CALL timeset(routineN, handle)
932 :
933 1696 : NULLIFY (basis_1c_set)
934 1696 : NULLIFY (my_CG, grid_atom, harmonics, atom_list)
935 :
936 1696 : CPASSERT(ASSOCIATED(atomic_kind_set))
937 1696 : CPASSERT(ASSOCIATED(dft_control))
938 1696 : CPASSERT(ASSOCIATED(para_env))
939 1696 : CPASSERT(ASSOCIATED(qs_kind_set))
940 :
941 1696 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
942 :
943 1696 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="GAPW_1C")
944 :
945 1696 : nspins = dft_control%nspins
946 1696 : gapw_control => dft_control%qs_control%gapw_control
947 :
948 1696 : lmax_sphere = gapw_control%lmax_sphere
949 :
950 1696 : llmax = MIN(lmax_sphere, 2*maxlgto)
951 1696 : max_s_harm = nsoset(llmax)
952 1696 : max_s_set = nsoset(maxlgto)
953 :
954 : ! *** calculate the Gaunt coefficients up to the maxl ***
955 1696 : CALL real_gaunt_coefficients(my_CG, maxlgto, maxlgto, llmax)
956 :
957 : ! *** initialize the Lebedev grids ***
958 1696 : CALL init_lebedev_grids()
959 1696 : quadrature = gapw_control%quadrature
960 :
961 4722 : DO ikind = 1, SIZE(atomic_kind_set)
962 3026 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
963 : CALL get_qs_kind(qs_kind_set(ikind), &
964 : paw_atom=paw_atom, &
965 : grid_atom=grid_atom, &
966 : harmonics=harmonics, &
967 3026 : ngrid_rad=nr, ngrid_ang=na)
968 :
969 : ! *** determine the Lebedev grid for this kind ***
970 :
971 3026 : ll = get_number_of_lebedev_grid(n=na)
972 3026 : na = lebedev_grid(ll)%n
973 3026 : la = lebedev_grid(ll)%l
974 3026 : grid_atom%ng_sphere = na
975 3026 : grid_atom%nr = nr
976 :
977 3026 : IF (llmax > la) THEN
978 0 : IF (para_env%is_source()) THEN
979 0 : WRITE (my_la, "(I4)") la
980 0 : WRITE (my_ikind, "(I4)") ikind
981 0 : WRITE (my_llmax, "(I4)") llmax
982 : CALL cp_warn(__LOCATION__, &
983 : "The Lebedev grid has been built for angular momentum l up "// &
984 : "to "//TRIM(my_la)//" for atom kind "//TRIM(my_ikind)//", "// &
985 : "but the max l of spherical harmonics "//TRIM(my_llmax)//" is "// &
986 0 : "larger and good integration cannot be guaranteed")
987 : END IF
988 : END IF
989 :
990 : ! *** calculate the radial grid ***
991 3026 : CALL create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
992 :
993 : ! *** calculate the spherical harmonics on the grid ***
994 :
995 3026 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_1c_set, basis_type="GAPW_1C")
996 3026 : CALL get_gto_basis_set(gto_basis_set=basis_1c_set, maxl=maxl)
997 3026 : maxs = nsoset(maxl)
998 : CALL create_harmonics_atom(harmonics, &
999 : my_CG, na, llmax, maxs, max_s_harm, ll, grid_atom%wa, &
1000 3026 : grid_atom%azi, grid_atom%pol)
1001 10774 : CALL get_maxl_CG(harmonics, basis_1c_set, llmax, max_s_harm)
1002 :
1003 : END DO
1004 :
1005 1696 : CALL deallocate_lebedev_grids()
1006 1696 : DEALLOCATE (my_CG)
1007 :
1008 1696 : CALL allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
1009 :
1010 1696 : CALL timestop(handle)
1011 :
1012 1696 : END SUBROUTINE init_rho_atom
1013 :
1014 : ! **************************************************************************************************
1015 : !> \brief ...
1016 : !> \param rho_atom_set ...
1017 : !> \param atomic_kind_set list of atomic kinds
1018 : !> \param qs_kind_set the kind set from which to take quantum numbers and basis info
1019 : !> \param dft_control DFT control type
1020 : !> \param para_env parallel environment
1021 : ! **************************************************************************************************
1022 5704 : SUBROUTINE allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
1023 :
1024 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
1025 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1026 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1027 : TYPE(dft_control_type), POINTER :: dft_control
1028 : TYPE(mp_para_env_type), POINTER :: para_env
1029 :
1030 : CHARACTER(len=*), PARAMETER :: routineN = 'allocate_rho_atom_internals'
1031 :
1032 : INTEGER :: bo(2), handle, iat, iatom, ikind, ispin, &
1033 : max_iso_not0, maxso, mepos, nat, &
1034 : natom, nsatbas, nset, nsotot, nspins, &
1035 : num_pe
1036 5704 : INTEGER, DIMENSION(:), POINTER :: atom_list
1037 : LOGICAL :: paw_atom
1038 : TYPE(gto_basis_set_type), POINTER :: basis_1c
1039 : TYPE(harmonics_atom_type), POINTER :: harmonics
1040 :
1041 5704 : CALL timeset(routineN, handle)
1042 :
1043 5704 : CPASSERT(ASSOCIATED(atomic_kind_set))
1044 5704 : CPASSERT(ASSOCIATED(dft_control))
1045 5704 : CPASSERT(ASSOCIATED(para_env))
1046 5704 : CPASSERT(ASSOCIATED(qs_kind_set))
1047 :
1048 5704 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
1049 :
1050 5704 : nspins = dft_control%nspins
1051 :
1052 5704 : IF (ASSOCIATED(rho_atom_set)) THEN
1053 0 : CALL deallocate_rho_atom_set(rho_atom_set)
1054 : END IF
1055 34584 : ALLOCATE (rho_atom_set(natom))
1056 :
1057 16908 : DO ikind = 1, SIZE(atomic_kind_set)
1058 :
1059 11204 : NULLIFY (atom_list, harmonics)
1060 11204 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
1061 : CALL get_qs_kind(qs_kind_set(ikind), &
1062 : paw_atom=paw_atom, &
1063 11204 : harmonics=harmonics)
1064 :
1065 11204 : IF (paw_atom) THEN
1066 10226 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
1067 10226 : CALL get_gto_basis_set(gto_basis_set=basis_1c, nset=nset, maxso=maxso)
1068 10226 : nsotot = nset*maxso
1069 10226 : CALL get_paw_basis_info(basis_1c, nsatbas=nsatbas)
1070 : END IF
1071 :
1072 11204 : max_iso_not0 = harmonics%max_iso_not0
1073 28676 : DO iat = 1, nat
1074 17472 : iatom = atom_list(iat)
1075 : ! *** allocate the radial density for each LM,for each atom ***
1076 :
1077 71522 : ALLOCATE (rho_atom_set(iatom)%rho_rad_h(nspins))
1078 54050 : ALLOCATE (rho_atom_set(iatom)%rho_rad_s(nspins))
1079 54050 : ALLOCATE (rho_atom_set(iatom)%vrho_rad_h(nspins))
1080 54050 : ALLOCATE (rho_atom_set(iatom)%vrho_rad_s(nspins))
1081 :
1082 : ALLOCATE (rho_atom_set(iatom)%cpc_h(nspins), &
1083 : rho_atom_set(iatom)%cpc_s(nspins), &
1084 : rho_atom_set(iatom)%drho_rad_h(nspins), &
1085 : rho_atom_set(iatom)%drho_rad_s(nspins), &
1086 : rho_atom_set(iatom)%rho_rad_h_d(3, nspins), &
1087 369048 : rho_atom_set(iatom)%rho_rad_s_d(3, nspins))
1088 : ALLOCATE (rho_atom_set(iatom)%int_scr_h(nspins), &
1089 90628 : rho_atom_set(iatom)%int_scr_s(nspins))
1090 :
1091 28676 : IF (paw_atom) THEN
1092 32536 : DO ispin = 1, nspins
1093 : ALLOCATE (rho_atom_set(iatom)%cpc_h(ispin)%r_coef(1:nsatbas, 1:nsatbas), &
1094 102252 : rho_atom_set(iatom)%cpc_s(ispin)%r_coef(1:nsatbas, 1:nsatbas))
1095 : ALLOCATE (rho_atom_set(iatom)%int_scr_h(ispin)%r_coef(1:nsatbas, 1:nsatbas), &
1096 85210 : rho_atom_set(iatom)%int_scr_s(ispin)%r_coef(1:nsatbas, 1:nsatbas))
1097 :
1098 8522026 : rho_atom_set(iatom)%cpc_h(ispin)%r_coef = 0.0_dp
1099 8537520 : rho_atom_set(iatom)%cpc_s(ispin)%r_coef = 0.0_dp
1100 : END DO
1101 : END IF
1102 :
1103 : END DO ! iat
1104 :
1105 11204 : num_pe = para_env%num_pe
1106 11204 : mepos = para_env%mepos
1107 11204 : bo = get_limit(nat, num_pe, mepos)
1108 36848 : DO iat = bo(1), bo(2)
1109 8736 : iatom = atom_list(iat)
1110 : ALLOCATE (rho_atom_set(iatom)%ga_Vlocal_gb_h(nspins), &
1111 54050 : rho_atom_set(iatom)%ga_Vlocal_gb_s(nspins))
1112 19940 : IF (paw_atom) THEN
1113 16268 : DO ispin = 1, nspins
1114 : CALL reallocate(rho_atom_set(iatom)%ga_Vlocal_gb_h(ispin)%r_coef, &
1115 8521 : 1, nsotot, 1, nsotot)
1116 : CALL reallocate(rho_atom_set(iatom)%ga_Vlocal_gb_s(ispin)%r_coef, &
1117 8521 : 1, nsotot, 1, nsotot)
1118 :
1119 23805935 : rho_atom_set(iatom)%ga_Vlocal_gb_h(ispin)%r_coef = 0.0_dp
1120 23813682 : rho_atom_set(iatom)%ga_Vlocal_gb_s(ispin)%r_coef = 0.0_dp
1121 : END DO
1122 : END IF
1123 :
1124 : END DO ! iat
1125 :
1126 : END DO
1127 :
1128 5704 : CALL timestop(handle)
1129 :
1130 11408 : END SUBROUTINE allocate_rho_atom_internals
1131 :
1132 : ! **************************************************************************************************
1133 : !> \brief ...
1134 : !> \param rho_atom_set ...
1135 : !> \param iatom ...
1136 : !> \param ispin ...
1137 : !> \param nr ...
1138 : !> \param max_iso_not0 ...
1139 : ! **************************************************************************************************
1140 8318 : SUBROUTINE allocate_rho_atom_rad(rho_atom_set, iatom, ispin, nr, max_iso_not0)
1141 :
1142 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
1143 : INTEGER, INTENT(IN) :: iatom, ispin, nr, max_iso_not0
1144 :
1145 : CHARACTER(len=*), PARAMETER :: routineN = 'allocate_rho_atom_rad'
1146 :
1147 : INTEGER :: handle, j
1148 :
1149 8318 : CALL timeset(routineN, handle)
1150 :
1151 : ALLOCATE (rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef(1:nr, 1:max_iso_not0), &
1152 : rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef(1:nr, 1:max_iso_not0), &
1153 : rho_atom_set(iatom)%vrho_rad_h(ispin)%r_coef(1:nr, 1:max_iso_not0), &
1154 83180 : rho_atom_set(iatom)%vrho_rad_s(ispin)%r_coef(1:nr, 1:max_iso_not0))
1155 :
1156 6598742 : rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef = 0.0_dp
1157 6598742 : rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef = 0.0_dp
1158 6598742 : rho_atom_set(iatom)%vrho_rad_h(ispin)%r_coef = 0.0_dp
1159 6598742 : rho_atom_set(iatom)%vrho_rad_s(ispin)%r_coef = 0.0_dp
1160 :
1161 : ALLOCATE (rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef(nr, max_iso_not0), &
1162 49908 : rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef(nr, max_iso_not0))
1163 6598742 : rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef = 0.0_dp
1164 6598742 : rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef = 0.0_dp
1165 :
1166 33272 : DO j = 1, 3
1167 : ALLOCATE (rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef(nr, max_iso_not0), &
1168 124770 : rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef(nr, max_iso_not0))
1169 19796226 : rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef = 0.0_dp
1170 19804544 : rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef = 0.0_dp
1171 : END DO
1172 :
1173 8318 : CALL timestop(handle)
1174 :
1175 8318 : END SUBROUTINE allocate_rho_atom_rad
1176 :
1177 : ! **************************************************************************************************
1178 : !> \brief ...
1179 : !> \param rho_atom_set ...
1180 : !> \param iatom ...
1181 : !> \param ispin ...
1182 : ! **************************************************************************************************
1183 59448 : SUBROUTINE set2zero_rho_atom_rad(rho_atom_set, iatom, ispin)
1184 :
1185 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
1186 : INTEGER, INTENT(IN) :: iatom, ispin
1187 :
1188 : INTEGER :: j
1189 :
1190 43424616 : rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef = 0.0_dp
1191 43424616 : rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef = 0.0_dp
1192 :
1193 43424616 : rho_atom_set(iatom)%vrho_rad_h(ispin)%r_coef = 0.0_dp
1194 43424616 : rho_atom_set(iatom)%vrho_rad_s(ispin)%r_coef = 0.0_dp
1195 :
1196 43424616 : rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef = 0.0_dp
1197 43424616 : rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef = 0.0_dp
1198 :
1199 237792 : DO j = 1, 3
1200 130273848 : rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef = 0.0_dp
1201 130333296 : rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef = 0.0_dp
1202 : END DO
1203 :
1204 59448 : END SUBROUTINE set2zero_rho_atom_rad
1205 :
1206 : END MODULE qs_rho_atom_methods
|
