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 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 : replicate_rho_atom_radial
88 :
89 : CONTAINS
90 :
91 : ! **************************************************************************************************
92 : !> \brief ...
93 : !> \param para_env ...
94 : !> \param rho_atom_set ...
95 : !> \param qs_kind ...
96 : !> \param atom_list ...
97 : !> \param natom ...
98 : !> \param nspins ...
99 : !> \param tot_rho1_h ...
100 : !> \param tot_rho1_s ...
101 : !> \param rho1_h_spin ...
102 : !> \param rho1_s_spin ...
103 : !> \param rho1_h_aspin ...
104 : !> \param rho1_s_aspin ...
105 : ! **************************************************************************************************
106 76292 : SUBROUTINE calculate_rho_atom(para_env, rho_atom_set, qs_kind, atom_list, &
107 76292 : natom, nspins, tot_rho1_h, tot_rho1_s, &
108 : rho1_h_spin, rho1_s_spin, rho1_h_aspin, rho1_s_aspin)
109 :
110 : TYPE(mp_para_env_type), POINTER :: para_env
111 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
112 : TYPE(qs_kind_type), INTENT(IN) :: qs_kind
113 : INTEGER, DIMENSION(:), INTENT(IN) :: atom_list
114 : INTEGER, INTENT(IN) :: natom, nspins
115 : REAL(dp), DIMENSION(:), INTENT(INOUT) :: tot_rho1_h, tot_rho1_s
116 : REAL(dp), INTENT(INOUT) :: rho1_h_spin, rho1_s_spin, rho1_h_aspin, &
117 : rho1_s_aspin
118 :
119 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_atom'
120 :
121 : INTEGER :: damax_iso_not0_local, handle, i, i1, i2, iat, iatom, icg, ipgf1, ipgf2, ir, &
122 : iset1, iset2, iso, iso1, iso1_coeff, iso1_first, iso1_last, iso2, iso2_coeff, iso2_first, &
123 : iso2_last, j, l, l_iso, l_sub, l_sum, lmax12, lmax_expansion, lmin12, m1s, m2s, &
124 : max_iso_not0, max_iso_not0_local, max_npgf, max_s_harm, maxl, maxso, mepos, n1s, n2s, na, &
125 : nr, nset, num_pe, size1, size2
126 76292 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list, dacg_n_list
127 76292 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list, dacg_list
128 : INTEGER, DIMENSION(2) :: bo
129 76292 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, o2nindex
130 76292 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: done_vgg
131 : REAL(dp) :: c1, c2, cpc_h, cpc_s, rfun, rho_h, &
132 : rho_s, root_zet12, zet12
133 76292 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: erf_zet12, g1, g2, gg0, int1, int2
134 76292 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dgg, gg, gg_lm1, sfun_h, sfun_s
135 76292 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: g_rad, vgg
136 76292 : REAL(dp), DIMENSION(:, :), POINTER :: coeff_h, coeff_s, zet
137 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
138 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_CG_dxyz
139 : TYPE(grid_atom_type), POINTER :: grid_atom
140 : TYPE(gto_basis_set_type), POINTER :: basis_1c
141 : TYPE(harmonics_atom_type), POINTER :: harmonics
142 :
143 76292 : CALL timeset(routineN, handle)
144 :
145 : !Note: tau is taken care of separately in qs_vxc_atom.F
146 :
147 76292 : NULLIFY (basis_1c)
148 76292 : NULLIFY (harmonics, grid_atom)
149 76292 : NULLIFY (lmin, lmax, npgf, zet, my_CG, my_CG_dxyz, coeff_h, coeff_s)
150 :
151 76292 : CALL get_qs_kind(qs_kind, grid_atom=grid_atom, harmonics=harmonics)
152 76292 : CALL get_qs_kind(qs_kind, basis_set=basis_1c, basis_type="GAPW_1C")
153 :
154 : CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
155 : maxl=maxl, npgf=npgf, nset=nset, zet=zet, &
156 76292 : maxso=maxso)
157 :
158 76292 : CALL get_paw_basis_info(basis_1c, o2nindex=o2nindex)
159 :
160 76292 : max_iso_not0 = harmonics%max_iso_not0
161 76292 : max_s_harm = harmonics%max_s_harm
162 :
163 76292 : nr = grid_atom%nr
164 277798 : max_npgf = MAXVAL(npgf(1:nset))
165 76292 : lmax_expansion = indso(1, max_iso_not0)
166 : ! Distribute the atoms of this kind
167 76292 : num_pe = para_env%num_pe
168 76292 : mepos = para_env%mepos
169 76292 : bo = get_limit(natom, num_pe, mepos)
170 :
171 76292 : my_CG => harmonics%my_CG
172 76292 : my_CG_dxyz => harmonics%my_CG_dxyz
173 :
174 915504 : ALLOCATE (g1(nr), g2(nr), gg0(nr), gg(nr, 0:2*maxl), dgg(nr, 0:2*maxl), gg_lm1(nr, 0:2*maxl))
175 457752 : ALLOCATE (erf_zet12(nr), vgg(nr, 0:2*maxl, 0:indso(1, max_iso_not0)))
176 305168 : ALLOCATE (done_vgg(0:2*maxl, 0:indso(1, max_iso_not0)))
177 228876 : ALLOCATE (int1(nr), int2(nr))
178 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm), &
179 686628 : dacg_list(2, nsoset(maxl)**2, max_s_harm), dacg_n_list(max_s_harm))
180 381460 : ALLOCATE (g_rad(nr, max_npgf, nset))
181 :
182 277798 : DO iset1 = 1, nset
183 780618 : DO ipgf1 = 1, npgf(iset1)
184 27436922 : g_rad(1:nr, ipgf1, iset1) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
185 : END DO
186 : END DO
187 :
188 134474 : DO iat = bo(1), bo(2)
189 58182 : iatom = atom_list(iat)
190 199487 : DO i = 1, nspins
191 123195 : IF (.NOT. ASSOCIATED(rho_atom_set(iatom)%rho_rad_h(i)%r_coef)) THEN
192 7971 : CALL allocate_rho_atom_rad(rho_atom_set, iatom, i, nr, max_iso_not0)
193 : ELSE
194 57042 : CALL set2zero_rho_atom_rad(rho_atom_set, iatom, i)
195 : END IF
196 : END DO
197 : END DO
198 :
199 : m1s = 0
200 277798 : DO iset1 = 1, nset
201 : m2s = 0
202 914664 : DO iset2 = 1, nset
203 :
204 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
205 713158 : max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
206 713158 : CPASSERT(max_iso_not0_local <= max_iso_not0)
207 : CALL get_none0_cg_list(my_CG_dxyz, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
208 713158 : max_s_harm, lmax_expansion, dacg_list, dacg_n_list, damax_iso_not0_local)
209 713158 : n1s = nsoset(lmax(iset1))
210 :
211 2337254 : DO ipgf1 = 1, npgf(iset1)
212 1624096 : iso1_first = nsoset(lmin(iset1) - 1) + 1 + n1s*(ipgf1 - 1) + m1s
213 1624096 : iso1_last = nsoset(lmax(iset1)) + n1s*(ipgf1 - 1) + m1s
214 1624096 : size1 = iso1_last - iso1_first + 1
215 1624096 : iso1_first = o2nindex(iso1_first)
216 1624096 : iso1_last = o2nindex(iso1_last)
217 1624096 : i1 = iso1_last - iso1_first + 1
218 1624096 : CPASSERT(size1 == i1)
219 1624096 : i1 = nsoset(lmin(iset1) - 1) + 1
220 :
221 87385952 : g1(1:nr) = g_rad(1:nr, ipgf1, iset1)
222 :
223 1624096 : n2s = nsoset(lmax(iset2))
224 6758782 : DO ipgf2 = 1, npgf(iset2)
225 4421528 : iso2_first = nsoset(lmin(iset2) - 1) + 1 + n2s*(ipgf2 - 1) + m2s
226 4421528 : iso2_last = nsoset(lmax(iset2)) + n2s*(ipgf2 - 1) + m2s
227 4421528 : size2 = iso2_last - iso2_first + 1
228 4421528 : iso2_first = o2nindex(iso2_first)
229 4421528 : iso2_last = o2nindex(iso2_last)
230 4421528 : i2 = iso2_last - iso2_first + 1
231 4421528 : CPASSERT(size2 == i2)
232 4421528 : i2 = nsoset(lmin(iset2) - 1) + 1
233 :
234 239103332 : g2(1:nr) = g_rad(1:nr, ipgf2, iset2)
235 4421528 : lmin12 = lmin(iset1) + lmin(iset2)
236 4421528 : lmax12 = lmax(iset1) + lmax(iset2)
237 :
238 4421528 : zet12 = zet(ipgf1, iset1) + zet(ipgf2, iset2)
239 4421528 : root_zet12 = SQRT(zet(ipgf1, iset1) + zet(ipgf2, iset2))
240 239103332 : DO ir = 1, nr
241 239103332 : erf_zet12(ir) = erf(root_zet12*grid_atom%rad(ir))
242 : END DO
243 :
244 4421528 : gg = 0.0_dp
245 4421528 : dgg = 0.0_dp
246 4421528 : gg_lm1 = 0.0_dp
247 4421528 : vgg = 0.0_dp
248 4421528 : done_vgg = .FALSE.
249 : ! reduce the number of terms in the expansion local densities
250 4421528 : IF (lmin12 <= lmax_expansion) THEN
251 4419038 : IF (lmin12 == 0) THEN
252 131076198 : gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
253 131076198 : gg_lm1(1:nr, lmin12) = 0.0_dp
254 131076198 : gg0(1:nr) = gg(1:nr, lmin12)
255 : ELSE
256 107900144 : gg0(1:nr) = g1(1:nr)*g2(1:nr)
257 107900144 : gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
258 107900144 : gg_lm1(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12 - 1)*g1(1:nr)*g2(1:nr)
259 : END IF
260 :
261 : ! reduce the number of terms in the expansion local densities
262 4419038 : IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
263 :
264 6967130 : DO l = lmin12 + 1, lmax12
265 140492212 : gg(1:nr, l) = grid_atom%rad(1:nr)*gg(1:nr, l - 1)
266 140492212 : gg_lm1(1:nr, l) = gg(1:nr, l - 1)
267 144911250 : dgg(1:nr, l - 1) = -2.0_dp*(zet(ipgf1, iset1) + zet(ipgf2, iset2))*gg(1:nr, l)
268 :
269 : END DO
270 : dgg(1:nr, lmax12) = -2.0_dp*(zet(ipgf1, iset1) + &
271 238976342 : zet(ipgf2, iset2))*grid_atom%rad(1:nr)*gg(1:nr, lmax12)
272 :
273 4419038 : c2 = SQRT(pi*pi*pi/(zet12*zet12*zet12))
274 :
275 37170372 : DO iso = 1, max_iso_not0_local
276 32751334 : l_iso = indso(1, iso)
277 32751334 : c1 = fourpi/(2._dp*REAL(l_iso, dp) + 1._dp)
278 109493144 : DO icg = 1, cg_n_list(iso)
279 72322772 : iso1 = cg_list(1, icg, iso)
280 72322772 : iso2 = cg_list(2, icg, iso)
281 :
282 72322772 : l = indso(1, iso1) + indso(1, iso2)
283 72322772 : CPASSERT(l <= lmax_expansion)
284 72322772 : IF (done_vgg(l, l_iso)) CYCLE
285 9129672 : L_sum = l + l_iso
286 9129672 : L_sub = l - l_iso
287 :
288 9129672 : IF (l_sum == 0) THEN
289 131076198 : vgg(1:nr, l, l_iso) = erf_zet12(1:nr)*grid_atom%oorad2l(1:nr, 1)*c2
290 : ELSE
291 6787278 : CALL whittaker_c0a(int1, grid_atom%rad, gg0, erf_zet12, zet12, l, l_iso, nr)
292 6787278 : CALL whittaker_ci(int2, grid_atom%rad, gg0, zet12, L_sub, nr)
293 :
294 362402538 : DO ir = 1, nr
295 355615260 : int2(ir) = grid_atom%rad2l(ir, l_iso)*int2(ir)
296 362402538 : vgg(ir, l, l_iso) = c1*(int1(ir) + int2(ir))
297 : END DO
298 : END IF
299 105074106 : done_vgg(l, l_iso) = .TRUE.
300 : END DO
301 : END DO
302 : END IF ! lmax_expansion
303 :
304 9069892 : DO iat = bo(1), bo(2)
305 3024268 : iatom = atom_list(iat)
306 :
307 10951831 : DO i = 1, nspins
308 3506035 : coeff_h => rho_atom_set(iatom)%cpc_h(i)%r_coef
309 3506035 : coeff_s => rho_atom_set(iatom)%cpc_s(i)%r_coef
310 :
311 25884617 : DO iso = 1, max_iso_not0_local
312 22378582 : l_iso = indso(1, iso)
313 73127258 : DO icg = 1, cg_n_list(iso)
314 47242641 : iso1 = cg_list(1, icg, iso)
315 47242641 : iso2 = cg_list(2, icg, iso)
316 :
317 47242641 : l = indso(1, iso1) + indso(1, iso2)
318 47242641 : CPASSERT(l <= lmax_expansion)
319 47242641 : iso1_coeff = iso1_first + iso1 - i1
320 47242641 : iso2_coeff = iso2_first + iso2 - i2
321 47242641 : cpc_h = coeff_h(iso1_coeff, iso2_coeff)*my_CG(iso1, iso2, iso)
322 47242641 : cpc_s = coeff_s(iso1_coeff, iso2_coeff)*my_CG(iso1, iso2, iso)
323 :
324 : rho_atom_set(iatom)%rho_rad_h(i)%r_coef(1:nr, iso) = &
325 : rho_atom_set(iatom)%rho_rad_h(i)%r_coef(1:nr, iso) + &
326 2517555005 : gg(1:nr, l)*cpc_h
327 :
328 : rho_atom_set(iatom)%rho_rad_s(i)%r_coef(1:nr, iso) = &
329 : rho_atom_set(iatom)%rho_rad_s(i)%r_coef(1:nr, iso) + &
330 2517555005 : gg(1:nr, l)*cpc_s
331 :
332 : rho_atom_set(iatom)%drho_rad_h(i)%r_coef(1:nr, iso) = &
333 : rho_atom_set(iatom)%drho_rad_h(i)%r_coef(1:nr, iso) + &
334 2517555005 : dgg(1:nr, l)*cpc_h
335 :
336 : rho_atom_set(iatom)%drho_rad_s(i)%r_coef(1:nr, iso) = &
337 : rho_atom_set(iatom)%drho_rad_s(i)%r_coef(1:nr, iso) + &
338 2517555005 : dgg(1:nr, l)*cpc_s
339 :
340 : rho_atom_set(iatom)%vrho_rad_h(i)%r_coef(1:nr, iso) = &
341 : rho_atom_set(iatom)%vrho_rad_h(i)%r_coef(1:nr, iso) + &
342 2517555005 : vgg(1:nr, l, l_iso)*cpc_h
343 :
344 : rho_atom_set(iatom)%vrho_rad_s(i)%r_coef(1:nr, iso) = &
345 : rho_atom_set(iatom)%vrho_rad_s(i)%r_coef(1:nr, iso) + &
346 2539933587 : vgg(1:nr, l, l_iso)*cpc_s
347 :
348 : END DO ! icg
349 :
350 : END DO ! iso
351 :
352 76530858 : DO iso = 1, max_iso_not0 !damax_iso_not0_local
353 70000555 : l_iso = indso(1, iso)
354 150885200 : DO icg = 1, dacg_n_list(iso)
355 77378610 : iso1 = dacg_list(1, icg, iso)
356 77378610 : iso2 = dacg_list(2, icg, iso)
357 77378610 : l = indso(1, iso1) + indso(1, iso2)
358 77378610 : CPASSERT(l <= lmax_expansion)
359 77378610 : iso1_coeff = iso1_first + iso1 - i1
360 77378610 : iso2_coeff = iso2_first + iso2 - i2
361 77378610 : cpc_h = coeff_h(iso1_coeff, iso2_coeff)
362 77378610 : cpc_s = coeff_s(iso1_coeff, iso2_coeff)
363 379514995 : DO j = 1, 3
364 : rho_atom_set(iatom)%rho_rad_h_d(j, i)%r_coef(1:nr, iso) = &
365 : rho_atom_set(iatom)%rho_rad_h_d(j, i)%r_coef(1:nr, iso) + &
366 12233568900 : gg_lm1(1:nr, l)*cpc_h*my_CG_dxyz(j, iso1, iso2, iso)
367 :
368 : rho_atom_set(iatom)%rho_rad_s_d(j, i)%r_coef(1:nr, iso) = &
369 : rho_atom_set(iatom)%rho_rad_s_d(j, i)%r_coef(1:nr, iso) + &
370 12310947510 : gg_lm1(1:nr, l)*cpc_s*my_CG_dxyz(j, iso1, iso2, iso)
371 : END DO
372 : END DO ! icg
373 :
374 : END DO ! iso
375 :
376 : END DO ! i
377 : END DO ! iat
378 :
379 : END DO ! ipgf2
380 : END DO ! ipgf1
381 2340980 : m2s = m2s + maxso
382 : END DO ! iset2
383 277798 : m1s = m1s + maxso
384 : END DO ! iset1
385 :
386 134474 : DO iat = bo(1), bo(2)
387 58182 : iatom = atom_list(iat)
388 :
389 123195 : DO i = 1, nspins
390 :
391 994720 : DO iso = 1, max_iso_not0
392 : rho_s = 0.0_dp
393 : rho_h = 0.0_dp
394 47193319 : DO ir = 1, nr
395 46321794 : rho_h = rho_h + rho_atom_set(iatom)%rho_rad_h(i)%r_coef(ir, iso)*grid_atom%wr(ir)
396 47193319 : rho_s = rho_s + rho_atom_set(iatom)%rho_rad_s(i)%r_coef(ir, iso)*grid_atom%wr(ir)
397 : END DO ! ir
398 871525 : tot_rho1_h(i) = tot_rho1_h(i) + rho_h*harmonics%slm_int(iso)
399 936538 : tot_rho1_s(i) = tot_rho1_s(i) + rho_s*harmonics%slm_int(iso)
400 : END DO ! iso
401 :
402 : END DO ! ispin
403 :
404 134474 : IF (nspins == 2) THEN
405 6831 : na = SIZE(harmonics%slm, 1)
406 40986 : ALLOCATE (sfun_h(nr, na), sfun_s(nr, na))
407 6831 : sfun_h = 0.0_dp
408 6831 : sfun_s = 0.0_dp
409 97150 : DO iso = 1, max_iso_not0
410 5793820 : DO ir = 1, nr
411 : rfun = grid_atom%wr(ir)*(rho_atom_set(iatom)%rho_rad_h(1)%r_coef(ir, iso) - &
412 5696670 : rho_atom_set(iatom)%rho_rad_h(2)%r_coef(ir, iso))
413 290386170 : sfun_h(ir, 1:na) = sfun_h(ir, 1:na) + rfun*harmonics%slm(1:na, iso)*grid_atom%wa(1:na)
414 : rfun = grid_atom%wr(ir)*(rho_atom_set(iatom)%rho_rad_s(1)%r_coef(ir, iso) - &
415 5696670 : rho_atom_set(iatom)%rho_rad_s(2)%r_coef(ir, iso))
416 290476489 : sfun_s(ir, 1:na) = sfun_s(ir, 1:na) + rfun*harmonics%slm(1:na, iso)*grid_atom%wa(1:na)
417 : END DO
418 : END DO
419 20831857 : rho1_h_spin = rho1_h_spin + SUM(sfun_h(1:nr, 1:na))
420 20831857 : rho1_s_spin = rho1_s_spin + SUM(sfun_s(1:nr, 1:na))
421 20831857 : rho1_h_aspin = rho1_h_aspin + SUM(ABS(sfun_h(1:nr, 1:na)))
422 20831857 : rho1_s_aspin = rho1_s_aspin + SUM(ABS(sfun_s(1:nr, 1:na)))
423 6831 : DEALLOCATE (sfun_h, sfun_s)
424 : END IF
425 :
426 : END DO ! iat
427 :
428 76292 : DEALLOCATE (g1, g2, gg0, gg, gg_lm1, dgg, vgg, done_vgg, erf_zet12, int1, int2, g_rad)
429 76292 : DEALLOCATE (cg_list, cg_n_list, dacg_list, dacg_n_list)
430 76292 : DEALLOCATE (o2nindex)
431 :
432 76292 : CALL timestop(handle)
433 :
434 228876 : END SUBROUTINE calculate_rho_atom
435 :
436 : ! **************************************************************************************************
437 : !> \brief Replicate the radial hard/soft density data needed to evaluate one-center tails on
438 : !> rank-local target grids. The compact one-center density matrices are already global;
439 : !> this routine performs one packed reduction for the derived radial fields of a kind.
440 : !> \param para_env ...
441 : !> \param rho_atom_set ...
442 : !> \param qs_kind ...
443 : !> \param atom_list ...
444 : !> \param natom ...
445 : !> \param nspins ...
446 : ! **************************************************************************************************
447 32 : SUBROUTINE replicate_rho_atom_radial(para_env, rho_atom_set, qs_kind, atom_list, natom, nspins)
448 : TYPE(mp_para_env_type), POINTER :: para_env
449 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
450 : TYPE(qs_kind_type), INTENT(IN) :: qs_kind
451 : INTEGER, DIMENSION(:), INTENT(IN) :: atom_list
452 : INTEGER, INTENT(IN) :: natom, nspins
453 :
454 : CHARACTER(len=*), PARAMETER :: routineN = 'replicate_rho_atom_radial'
455 :
456 : INTEGER :: block_size, bo(2), cursor, handle, iat, &
457 : iatom, ispin, j, max_iso_not0, ncoeff, &
458 : nr
459 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: buffer
460 : TYPE(grid_atom_type), POINTER :: grid_atom
461 : TYPE(harmonics_atom_type), POINTER :: harmonics
462 :
463 32 : CALL timeset(routineN, handle)
464 :
465 32 : NULLIFY (grid_atom, harmonics)
466 32 : CALL get_qs_kind(qs_kind, grid_atom=grid_atom, harmonics=harmonics)
467 32 : CPASSERT(ASSOCIATED(grid_atom))
468 32 : CPASSERT(ASSOCIATED(harmonics))
469 32 : nr = grid_atom%nr
470 32 : max_iso_not0 = harmonics%max_iso_not0
471 32 : ncoeff = nr*max_iso_not0
472 32 : block_size = 10*nspins*ncoeff
473 96 : ALLOCATE (buffer(block_size*natom))
474 32 : buffer = 0.0_dp
475 :
476 32 : bo = get_limit(natom, para_env%num_pe, para_env%mepos)
477 60 : DO iat = bo(1), bo(2)
478 28 : iatom = atom_list(iat)
479 28 : cursor = (iat - 1)*block_size
480 56 : DO ispin = 1, nspins
481 : buffer(cursor + 1:cursor + ncoeff) = &
482 56 : RESHAPE(rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef, [ncoeff])
483 28 : cursor = cursor + ncoeff
484 : buffer(cursor + 1:cursor + ncoeff) = &
485 56 : RESHAPE(rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef, [ncoeff])
486 28 : cursor = cursor + ncoeff
487 : buffer(cursor + 1:cursor + ncoeff) = &
488 56 : RESHAPE(rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef, [ncoeff])
489 28 : cursor = cursor + ncoeff
490 : buffer(cursor + 1:cursor + ncoeff) = &
491 56 : RESHAPE(rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef, [ncoeff])
492 28 : cursor = cursor + ncoeff
493 140 : DO j = 1, 3
494 : buffer(cursor + 1:cursor + ncoeff) = &
495 168 : RESHAPE(rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef, [ncoeff])
496 84 : cursor = cursor + ncoeff
497 : buffer(cursor + 1:cursor + ncoeff) = &
498 168 : RESHAPE(rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef, [ncoeff])
499 112 : cursor = cursor + ncoeff
500 : END DO
501 : END DO
502 60 : CPASSERT(cursor == iat*block_size)
503 : END DO
504 32 : CALL para_env%sum(buffer)
505 :
506 88 : DO iat = 1, natom
507 56 : iatom = atom_list(iat)
508 56 : cursor = (iat - 1)*block_size
509 112 : DO ispin = 1, nspins
510 56 : IF (.NOT. ASSOCIATED(rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef)) THEN
511 4 : CALL allocate_rho_atom_rad(rho_atom_set, iatom, ispin, nr, max_iso_not0)
512 : END IF
513 : rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef = &
514 24240 : RESHAPE(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
515 56 : cursor = cursor + ncoeff
516 : rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef = &
517 24240 : RESHAPE(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
518 56 : cursor = cursor + ncoeff
519 : rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef = &
520 24240 : RESHAPE(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
521 56 : cursor = cursor + ncoeff
522 : rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef = &
523 24240 : RESHAPE(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
524 56 : cursor = cursor + ncoeff
525 280 : DO j = 1, 3
526 : rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef = &
527 72720 : RESHAPE(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
528 168 : cursor = cursor + ncoeff
529 : rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef = &
530 72720 : RESHAPE(buffer(cursor + 1:cursor + ncoeff), [nr, max_iso_not0])
531 224 : cursor = cursor + ncoeff
532 : END DO
533 : END DO
534 88 : CPASSERT(cursor == iat*block_size)
535 : END DO
536 32 : DEALLOCATE (buffer)
537 :
538 32 : CALL timestop(handle)
539 :
540 32 : END SUBROUTINE replicate_rho_atom_radial
541 :
542 : ! **************************************************************************************************
543 : !> \brief ...
544 : !> \param qs_env QuickStep environment
545 : !> (accessed components: atomic_kind_set, dft_control%nimages,
546 : !> dft_control%nspins, kpoints%cell_to_index)
547 : !> \param rho_ao density matrix in atomic basis set
548 : !> \param rho_atom_set ...
549 : !> \param qs_kind_set list of QuickStep kinds
550 : !> \param oce one-centre expansion coefficients
551 : !> \param sab neighbour pair list
552 : !> \param para_env parallel environment
553 : !> \par History
554 : !> Add OpenMP [Apr 2016, EPCC]
555 : !> Use automatic arrays [Sep 2016, M Tucker]
556 : !> Allow for external non-default kind_set, oce and sab [Dec 2019, A Bussy]
557 : !> \note Consider to declare 'rho_ao' dummy argument as a pointer to the two-dimensional
558 : !> (1:nspins, 1:nimages) set of matrices.
559 : ! **************************************************************************************************
560 43790 : SUBROUTINE calculate_rho_atom_coeff(qs_env, rho_ao, rho_atom_set, qs_kind_set, oce, sab, para_env)
561 :
562 : TYPE(qs_environment_type), POINTER :: qs_env
563 : TYPE(dbcsr_p_type), DIMENSION(*) :: rho_ao
564 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
565 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
566 : TYPE(oce_matrix_type), POINTER :: oce
567 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
568 : POINTER :: sab
569 : TYPE(mp_para_env_type), POINTER :: para_env
570 :
571 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_atom_coeff'
572 :
573 : INTEGER :: bo(2), handle, i, iac, iatom, ibc, icol, ikind, img, irow, ispin, jatom, jkind, &
574 : kac, katom, kbc, kkind, len_CPC, len_PC1, max_gau, max_nsgf, mepos, n_cont_a, n_cont_b, &
575 : nat_kind, natom, nimages, nkind, nsoctot, nspins, num_pe
576 43790 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, nsatbas_kind
577 : INTEGER, DIMENSION(3) :: cell_b
578 43790 : INTEGER, DIMENSION(:), POINTER :: a_list, list_a, list_b
579 43790 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
580 : LOGICAL :: dista, distab, distb, found, paw_atom
581 43790 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: has_intac, paw_kind
582 43790 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: proj_work1, proj_work2
583 43790 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: p_matrix
584 : REAL(KIND=dp) :: eps_cpc, factor, pmax
585 : REAL(KIND=dp), DIMENSION(3) :: rab
586 43790 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: C_coeff_hh_a, C_coeff_hh_b, &
587 43790 : C_coeff_ss_a, C_coeff_ss_b, r_coef_h, &
588 43790 : r_coef_s
589 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
590 43790 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
591 : TYPE(dft_control_type), POINTER :: dft_control
592 43790 : 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 43790 : DIMENSION(:), POINTER :: nl_iterator
597 43790 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: p_block_spin
598 :
599 43790 : !$ INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: locks
600 : !$ INTEGER :: lock, number_of_locks
601 :
602 43790 : CALL timeset(routineN, handle)
603 :
604 : CALL get_qs_env(qs_env=qs_env, &
605 : dft_control=dft_control, &
606 43790 : atomic_kind_set=atomic_kind_set)
607 :
608 43790 : eps_cpc = dft_control%qs_control%gapw_control%eps_cpc
609 :
610 43790 : CPASSERT(ASSOCIATED(qs_kind_set))
611 43790 : CPASSERT(ASSOCIATED(rho_atom_set))
612 43790 : CPASSERT(ASSOCIATED(oce))
613 43790 : CPASSERT(ASSOCIATED(sab))
614 :
615 43790 : nspins = dft_control%nspins
616 43790 : nimages = dft_control%nimages
617 :
618 43790 : NULLIFY (cell_to_index)
619 43790 : IF (nimages > 1) THEN
620 506 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
621 506 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
622 : END IF
623 :
624 43790 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
625 43790 : CALL get_qs_kind_set(qs_kind_set, maxsgf=max_nsgf, maxgtops=max_gau, basis_type='GAPW_1C')
626 :
627 43790 : nkind = SIZE(atomic_kind_set)
628 : ! Initialize to 0 the CPC coefficients and the local density arrays
629 130532 : DO ikind = 1, nkind
630 86742 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=a_list, natom=nat_kind)
631 86742 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
632 :
633 86742 : IF (.NOT. paw_atom) CYCLE
634 200124 : DO i = 1, nat_kind
635 121576 : iatom = a_list(i)
636 335990 : DO ispin = 1, nspins
637 61745750 : rho_atom_set(iatom)%cpc_h(ispin)%r_coef = 0.0_dp
638 61867326 : rho_atom_set(iatom)%cpc_s(ispin)%r_coef = 0.0_dp
639 : END DO ! ispin
640 : END DO ! i
641 :
642 78548 : num_pe = para_env%num_pe
643 78548 : mepos = para_env%mepos
644 78548 : bo = get_limit(nat_kind, num_pe, mepos)
645 269868 : DO i = bo(1), bo(2)
646 60788 : iatom = a_list(i)
647 215463 : DO ispin = 1, nspins
648 188968375 : rho_atom_set(iatom)%ga_Vlocal_gb_h(ispin)%r_coef = 0.0_dp
649 189029163 : 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 218112 : ALLOCATE (basis_set_list(nkind))
655 262740 : ALLOCATE (paw_kind(nkind), nsatbas_kind(nkind), has_intac(nkind*nkind))
656 43790 : paw_kind(:) = .FALSE.
657 43790 : nsatbas_kind(:) = 0
658 43790 : has_intac(:) = .FALSE.
659 130532 : DO ikind = 1, nkind
660 86742 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a)
661 86742 : IF (ASSOCIATED(basis_set_a)) THEN
662 86742 : 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 86742 : paw_atom=paw_kind(ikind))
668 130532 : IF (paw_kind(ikind)) CALL get_paw_basis_info(basis_1c, nsatbas=nsatbas_kind(ikind))
669 : END DO
670 232628 : DO ikind = 1, nkind*nkind
671 232628 : has_intac(ikind) = ASSOCIATED(oce%intac(ikind)%alist)
672 : END DO
673 :
674 43790 : len_PC1 = max_nsgf*max_gau
675 43790 : len_CPC = max_gau*max_gau
676 :
677 : num_pe = 1
678 43790 : !$ num_pe = omp_get_max_threads()
679 43790 : CALL neighbor_list_iterator_create(nl_iterator, sab, nthread=num_pe)
680 :
681 : !$OMP PARALLEL DEFAULT( NONE ) &
682 : !$OMP SHARED( max_nsgf, max_gau &
683 : !$OMP , len_PC1, len_CPC &
684 : !$OMP , nl_iterator, basis_set_list &
685 : !$OMP , nimages, cell_to_index &
686 : !$OMP , nspins, rho_ao &
687 : !$OMP , nkind, qs_kind_set &
688 : !$OMP , oce, eps_cpc &
689 : !$OMP , rho_atom_set &
690 : !$OMP , natom, locks, number_of_locks &
691 : !$OMP , paw_kind, nsatbas_kind, has_intac &
692 : !$OMP ) &
693 : !$OMP PRIVATE( p_block_spin, ispin &
694 : !$OMP , p_matrix, proj_work1, proj_work2 &
695 : !$OMP , mepos &
696 : !$OMP , ikind, jkind, iatom, jatom &
697 : !$OMP , cell_b, rab &
698 : !$OMP , basis_set_a, basis_set_b &
699 : !$OMP , pmax, irow, icol, img &
700 : !$OMP , found &
701 : !$OMP , kkind &
702 : !$OMP , nsoctot, katom &
703 : !$OMP , iac , alist_ac, kac, n_cont_a, list_a &
704 : !$OMP , ibc , alist_bc, kbc, n_cont_b, list_b &
705 : !$OMP , C_coeff_hh_a, C_coeff_ss_a, dista &
706 : !$OMP , C_coeff_hh_b, C_coeff_ss_b, distb &
707 : !$OMP , distab &
708 : !$OMP , factor, r_coef_h, r_coef_s &
709 43790 : !$OMP )
710 :
711 : ALLOCATE (p_block_spin(nspins))
712 : ALLOCATE (p_matrix(max_nsgf, max_nsgf))
713 : ALLOCATE (proj_work1(len_PC1), proj_work2(len_CPC))
714 :
715 : !$OMP SINGLE
716 : !$ number_of_locks = nspins*natom
717 : !$ ALLOCATE (locks(number_of_locks))
718 : !$OMP END SINGLE
719 :
720 : !$OMP DO
721 : !$ DO lock = 1, number_of_locks
722 : !$ call omp_init_lock(locks(lock))
723 : !$ END DO
724 : !$OMP END DO
725 :
726 : mepos = 0
727 : !$ mepos = omp_get_thread_num()
728 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
729 :
730 : CALL get_iterator_info(nl_iterator, mepos=mepos, &
731 : ikind=ikind, jkind=jkind, &
732 : iatom=iatom, jatom=jatom, &
733 : cell=cell_b, r=rab)
734 :
735 : basis_set_a => basis_set_list(ikind)%gto_basis_set
736 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
737 : basis_set_b => basis_set_list(jkind)%gto_basis_set
738 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
739 :
740 : pmax = 0._dp
741 : IF (iatom <= jatom) THEN
742 : irow = iatom
743 : icol = jatom
744 : ELSE
745 : irow = jatom
746 : icol = iatom
747 : END IF
748 :
749 : IF (nimages > 1) THEN
750 : img = cell_to_index(cell_b(1), cell_b(2), cell_b(3))
751 : CPASSERT(img > 0)
752 : ELSE
753 : img = 1
754 : END IF
755 :
756 : DO ispin = 1, nspins
757 : CALL dbcsr_get_block_p(matrix=rho_ao(nspins*(img - 1) + ispin)%matrix, &
758 : row=irow, col=icol, BLOCK=p_block_spin(ispin)%r_coef, &
759 : found=found)
760 : pmax = pmax + MAXVAL(ABS(p_block_spin(ispin)%r_coef))
761 : END DO
762 :
763 : DO kkind = 1, nkind
764 : IF (.NOT. paw_kind(kkind)) CYCLE
765 :
766 : nsoctot = nsatbas_kind(kkind)
767 :
768 : iac = ikind + nkind*(kkind - 1)
769 : ibc = jkind + nkind*(kkind - 1)
770 : IF (.NOT. has_intac(iac)) CYCLE
771 : IF (.NOT. has_intac(ibc)) CYCLE
772 :
773 : CALL get_alist(oce%intac(iac), alist_ac, iatom)
774 : CALL get_alist(oce%intac(ibc), alist_bc, jatom)
775 : IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
776 : IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
777 :
778 : DO kac = 1, alist_ac%nclist
779 : DO kbc = 1, alist_bc%nclist
780 : IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
781 : IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
782 : IF (pmax*alist_bc%clist(kbc)%maxac*alist_ac%clist(kac)%maxac < eps_cpc) CYCLE
783 :
784 : n_cont_a = alist_ac%clist(kac)%nsgf_cnt
785 : n_cont_b = alist_bc%clist(kbc)%nsgf_cnt
786 : IF (n_cont_a == 0 .OR. n_cont_b == 0) CYCLE
787 :
788 : list_a => alist_ac%clist(kac)%sgf_list
789 : list_b => alist_bc%clist(kbc)%sgf_list
790 :
791 : katom = alist_ac%clist(kac)%catom
792 :
793 : IF (iatom == katom .AND. ALL(alist_ac%clist(kac)%cell == 0)) THEN
794 : C_coeff_hh_a => alist_ac%clist(kac)%achint(:, :, 1)
795 : C_coeff_ss_a => alist_ac%clist(kac)%acint(:, :, 1)
796 : dista = .FALSE.
797 : ELSE
798 : C_coeff_hh_a => alist_ac%clist(kac)%acint(:, :, 1)
799 : C_coeff_ss_a => alist_ac%clist(kac)%acint(:, :, 1)
800 : dista = .TRUE.
801 : END IF
802 : IF (jatom == katom .AND. ALL(alist_bc%clist(kbc)%cell == 0)) THEN
803 : C_coeff_hh_b => alist_bc%clist(kbc)%achint(:, :, 1)
804 : C_coeff_ss_b => alist_bc%clist(kbc)%acint(:, :, 1)
805 : distb = .FALSE.
806 : ELSE
807 : C_coeff_hh_b => alist_bc%clist(kbc)%acint(:, :, 1)
808 : C_coeff_ss_b => alist_bc%clist(kbc)%acint(:, :, 1)
809 : distb = .TRUE.
810 : END IF
811 :
812 : distab = dista .AND. distb
813 :
814 : DO ispin = 1, nspins
815 :
816 : IF (iatom <= jatom) THEN
817 : CALL alist_pre_align_blk(p_block_spin(ispin)%r_coef, &
818 : SIZE(p_block_spin(ispin)%r_coef, 1), p_matrix, SIZE(p_matrix, 1), &
819 : list_a, n_cont_a, list_b, n_cont_b)
820 : ELSE
821 : CALL alist_pre_align_blk(p_block_spin(ispin)%r_coef, &
822 : SIZE(p_block_spin(ispin)%r_coef, 1), p_matrix, SIZE(p_matrix, 1), &
823 : list_b, n_cont_b, list_a, n_cont_a)
824 : END IF
825 :
826 : factor = 1.0_dp
827 : IF (iatom == jatom) factor = 0.5_dp
828 :
829 : r_coef_h => rho_atom_set(katom)%cpc_h(ispin)%r_coef
830 : r_coef_s => rho_atom_set(katom)%cpc_s(ispin)%r_coef
831 :
832 : !$ CALL omp_set_lock(locks((katom - 1)*nspins + ispin))
833 : IF (iatom <= jatom) THEN
834 : CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, &
835 : C_coeff_hh_b, C_coeff_ss_b, n_cont_b, &
836 : p_matrix, max_nsgf, r_coef_h, r_coef_s, nsoctot, &
837 : len_PC1, len_CPC, factor, distab, proj_work1, proj_work2)
838 : ELSE
839 : CALL proj_blk(C_coeff_hh_b, C_coeff_ss_b, n_cont_b, &
840 : C_coeff_hh_a, C_coeff_ss_a, n_cont_a, &
841 : p_matrix, max_nsgf, r_coef_h, r_coef_s, nsoctot, &
842 : len_PC1, len_CPC, factor, distab, proj_work1, proj_work2)
843 : END IF
844 : !$ CALL omp_unset_lock(locks((katom - 1)*nspins + ispin))
845 :
846 : END DO
847 : EXIT !search loop over jatom-katom list
848 : END IF
849 : END DO
850 : END DO
851 : END DO
852 : END DO
853 : ! Wait for all threads to finish the loop before locks can be freed
854 : !$OMP BARRIER
855 :
856 : !$OMP DO
857 : !$ DO lock = 1, number_of_locks
858 : !$ call omp_destroy_lock(locks(lock))
859 : !$ END DO
860 : !$OMP END DO
861 : !$OMP SINGLE
862 : !$ DEALLOCATE (locks)
863 : !$OMP END SINGLE NOWAIT
864 :
865 : DEALLOCATE (p_block_spin, p_matrix, proj_work1, proj_work2)
866 : !$OMP END PARALLEL
867 :
868 43790 : CALL neighbor_list_iterator_release(nl_iterator)
869 :
870 43790 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
871 :
872 182854 : DO iatom = 1, natom
873 294214 : ikind = kind_of(iatom)
874 :
875 338004 : DO ispin = 1, nspins
876 294214 : IF (ASSOCIATED(rho_atom_set(iatom)%cpc_h(ispin)%r_coef)) THEN
877 123355634 : CALL para_env%sum(rho_atom_set(iatom)%cpc_h(ispin)%r_coef)
878 123355634 : CALL para_env%sum(rho_atom_set(iatom)%cpc_s(ispin)%r_coef)
879 135866 : r_coef_h => rho_atom_set(iatom)%cpc_h(ispin)%r_coef
880 135866 : r_coef_s => rho_atom_set(iatom)%cpc_s(ispin)%r_coef
881 123491500 : r_coef_h(:, :) = r_coef_h(:, :) + TRANSPOSE(r_coef_h(:, :))
882 123491500 : r_coef_s(:, :) = r_coef_s(:, :) + TRANSPOSE(r_coef_s(:, :))
883 : END IF
884 : END DO
885 :
886 : END DO
887 :
888 43790 : DEALLOCATE (kind_of, basis_set_list, paw_kind, nsatbas_kind, has_intac)
889 :
890 43790 : CALL timestop(handle)
891 :
892 131370 : END SUBROUTINE calculate_rho_atom_coeff
893 :
894 : ! **************************************************************************************************
895 : !> \brief ...
896 : !> \param rho_atom_set the type to initialize
897 : !> \param atomic_kind_set list of atomic kinds
898 : !> \param qs_kind_set the kind set from which to take quantum numbers and basis info
899 : !> \param dft_control DFT control type
900 : !> \param para_env parallel environment
901 : !> \par History:
902 : !> - Generalised by providing the rho_atom_set and the qs_kind_set 12.2019 (A.Bussy)
903 : ! **************************************************************************************************
904 1592 : SUBROUTINE init_rho_atom(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
905 :
906 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
907 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
908 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
909 : TYPE(dft_control_type), POINTER :: dft_control
910 : TYPE(mp_para_env_type), POINTER :: para_env
911 :
912 : CHARACTER(len=*), PARAMETER :: routineN = 'init_rho_atom'
913 :
914 : INTEGER :: handle, ikind, il, iso, iso1, iso2, l1, l1l2, l2, la, lc1, lc2, lcleb, ll, llmax, &
915 : lmax_sphere, lp, m1, m2, max_s_harm, max_s_set, maxl, maxlgto, maxs, mm, mp, na, nat, &
916 : natom, nr, nspins, quadrature
917 1592 : INTEGER, DIMENSION(:), POINTER :: atom_list
918 : LOGICAL :: paw_atom
919 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rga
920 1592 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
921 : TYPE(gapw_control_type), POINTER :: gapw_control
922 : TYPE(grid_atom_type), POINTER :: grid_atom
923 : TYPE(gto_basis_set_type), POINTER :: basis_1c_set
924 : TYPE(harmonics_atom_type), POINTER :: harmonics
925 :
926 1592 : CALL timeset(routineN, handle)
927 :
928 1592 : NULLIFY (basis_1c_set)
929 1592 : NULLIFY (my_CG, grid_atom, harmonics, atom_list)
930 :
931 1592 : CPASSERT(ASSOCIATED(atomic_kind_set))
932 1592 : CPASSERT(ASSOCIATED(dft_control))
933 1592 : CPASSERT(ASSOCIATED(para_env))
934 1592 : CPASSERT(ASSOCIATED(qs_kind_set))
935 :
936 1592 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
937 :
938 1592 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="GAPW_1C")
939 :
940 1592 : nspins = dft_control%nspins
941 1592 : gapw_control => dft_control%qs_control%gapw_control
942 :
943 1592 : lmax_sphere = gapw_control%lmax_sphere
944 :
945 1592 : llmax = MIN(lmax_sphere, 2*maxlgto)
946 1592 : max_s_harm = nsoset(llmax)
947 1592 : max_s_set = nsoset(maxlgto)
948 :
949 1592 : lcleb = MAX(llmax, 2*maxlgto, 1)
950 :
951 : ! *** allocate calculate the CG coefficients up to the maxl ***
952 1592 : CALL clebsch_gordon_init(lcleb)
953 1592 : CALL reallocate(my_CG, 1, max_s_set, 1, max_s_set, 1, max_s_harm)
954 :
955 4776 : ALLOCATE (rga(lcleb, 2))
956 5710 : DO lc1 = 0, maxlgto
957 17152 : DO iso1 = nsoset(lc1 - 1) + 1, nsoset(lc1)
958 11442 : l1 = indso(1, iso1)
959 11442 : m1 = indso(2, iso1)
960 48886 : DO lc2 = 0, maxlgto
961 145586 : DO iso2 = nsoset(lc2 - 1) + 1, nsoset(lc2)
962 100818 : l2 = indso(1, iso2)
963 100818 : m2 = indso(2, iso2)
964 100818 : CALL clebsch_gordon(l1, m1, l2, m2, rga)
965 100818 : IF (l1 + l2 > llmax) THEN
966 : l1l2 = llmax
967 : ELSE
968 : l1l2 = l1 + l2
969 : END IF
970 100818 : mp = m1 + m2
971 100818 : mm = m1 - m2
972 100818 : IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN
973 44688 : mp = -ABS(mp)
974 44688 : mm = -ABS(mm)
975 : ELSE
976 56130 : mp = ABS(mp)
977 56130 : mm = ABS(mm)
978 : END IF
979 366130 : DO lp = MOD(l1 + l2, 2), l1l2, 2
980 231986 : il = lp/2 + 1
981 231986 : IF (ABS(mp) <= lp) THEN
982 163894 : IF (mp >= 0) THEN
983 108534 : iso = nsoset(lp - 1) + lp + 1 + mp
984 : ELSE
985 55360 : iso = nsoset(lp - 1) + lp + 1 - ABS(mp)
986 : END IF
987 163894 : my_CG(iso1, iso2, iso) = rga(il, 1)
988 : END IF
989 332804 : IF (mp /= mm .AND. ABS(mm) <= lp) THEN
990 80044 : IF (mm >= 0) THEN
991 53212 : iso = nsoset(lp - 1) + lp + 1 + mm
992 : ELSE
993 26832 : iso = nsoset(lp - 1) + lp + 1 - ABS(mm)
994 : END IF
995 80044 : my_CG(iso1, iso2, iso) = rga(il, 2)
996 : END IF
997 : END DO
998 : END DO ! iso2
999 : END DO ! lc2
1000 : END DO ! iso1
1001 : END DO ! lc1
1002 1592 : DEALLOCATE (rga)
1003 1592 : CALL clebsch_gordon_deallocate()
1004 :
1005 : ! *** initialize the Lebedev grids ***
1006 1592 : CALL init_lebedev_grids()
1007 1592 : quadrature = gapw_control%quadrature
1008 :
1009 4480 : DO ikind = 1, SIZE(atomic_kind_set)
1010 2888 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
1011 : CALL get_qs_kind(qs_kind_set(ikind), &
1012 : paw_atom=paw_atom, &
1013 : grid_atom=grid_atom, &
1014 : harmonics=harmonics, &
1015 2888 : ngrid_rad=nr, ngrid_ang=na)
1016 :
1017 : ! *** determine the Lebedev grid for this kind ***
1018 :
1019 2888 : ll = get_number_of_lebedev_grid(n=na)
1020 2888 : na = lebedev_grid(ll)%n
1021 2888 : la = lebedev_grid(ll)%l
1022 2888 : grid_atom%ng_sphere = na
1023 2888 : grid_atom%nr = nr
1024 :
1025 2888 : IF (llmax > la) THEN
1026 0 : WRITE (*, '(/,72("*"))')
1027 : WRITE (*, '(T2,A,T66,I4)') &
1028 0 : "WARNING: the lebedev grid is built for angular momentum l up to ", la, &
1029 0 : " the max l of spherical harmonics is larger, l_max = ", llmax, &
1030 0 : " good integration is guaranteed only for l <= ", la
1031 0 : WRITE (*, '(72("*"),/)')
1032 : END IF
1033 :
1034 : ! *** calculate the radial grid ***
1035 2888 : CALL create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
1036 :
1037 : ! *** calculate the spherical harmonics on the grid ***
1038 :
1039 2888 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_1c_set, basis_type="GAPW_1C")
1040 2888 : CALL get_gto_basis_set(gto_basis_set=basis_1c_set, maxl=maxl)
1041 2888 : maxs = nsoset(maxl)
1042 : CALL create_harmonics_atom(harmonics, &
1043 : my_CG, na, llmax, maxs, max_s_harm, ll, grid_atom%wa, &
1044 2888 : grid_atom%azi, grid_atom%pol)
1045 10256 : CALL get_maxl_CG(harmonics, basis_1c_set, llmax, max_s_harm)
1046 :
1047 : END DO
1048 :
1049 1592 : CALL deallocate_lebedev_grids()
1050 1592 : DEALLOCATE (my_CG)
1051 :
1052 1592 : CALL allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
1053 :
1054 1592 : CALL timestop(handle)
1055 :
1056 3184 : END SUBROUTINE init_rho_atom
1057 :
1058 : ! **************************************************************************************************
1059 : !> \brief ...
1060 : !> \param rho_atom_set ...
1061 : !> \param atomic_kind_set list of atomic kinds
1062 : !> \param qs_kind_set the kind set from which to take quantum numbers and basis info
1063 : !> \param dft_control DFT control type
1064 : !> \param para_env parallel environment
1065 : ! **************************************************************************************************
1066 5504 : SUBROUTINE allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
1067 :
1068 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
1069 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1070 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1071 : TYPE(dft_control_type), POINTER :: dft_control
1072 : TYPE(mp_para_env_type), POINTER :: para_env
1073 :
1074 : CHARACTER(len=*), PARAMETER :: routineN = 'allocate_rho_atom_internals'
1075 :
1076 : INTEGER :: bo(2), handle, iat, iatom, ikind, ispin, &
1077 : max_iso_not0, maxso, mepos, nat, &
1078 : natom, nsatbas, nset, nsotot, nspins, &
1079 : num_pe
1080 5504 : INTEGER, DIMENSION(:), POINTER :: atom_list
1081 : LOGICAL :: paw_atom
1082 : TYPE(gto_basis_set_type), POINTER :: basis_1c
1083 : TYPE(harmonics_atom_type), POINTER :: harmonics
1084 :
1085 5504 : CALL timeset(routineN, handle)
1086 :
1087 5504 : CPASSERT(ASSOCIATED(atomic_kind_set))
1088 5504 : CPASSERT(ASSOCIATED(dft_control))
1089 5504 : CPASSERT(ASSOCIATED(para_env))
1090 5504 : CPASSERT(ASSOCIATED(qs_kind_set))
1091 :
1092 5504 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
1093 :
1094 5504 : nspins = dft_control%nspins
1095 :
1096 5504 : IF (ASSOCIATED(rho_atom_set)) THEN
1097 0 : CALL deallocate_rho_atom_set(rho_atom_set)
1098 : END IF
1099 33612 : ALLOCATE (rho_atom_set(natom))
1100 :
1101 16482 : DO ikind = 1, SIZE(atomic_kind_set)
1102 :
1103 10978 : NULLIFY (atom_list, harmonics)
1104 10978 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
1105 : CALL get_qs_kind(qs_kind_set(ikind), &
1106 : paw_atom=paw_atom, &
1107 10978 : harmonics=harmonics)
1108 :
1109 10978 : IF (paw_atom) THEN
1110 10020 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
1111 10020 : CALL get_gto_basis_set(gto_basis_set=basis_1c, nset=nset, maxso=maxso)
1112 10020 : nsotot = nset*maxso
1113 10020 : CALL get_paw_basis_info(basis_1c, nsatbas=nsatbas)
1114 : END IF
1115 :
1116 10978 : max_iso_not0 = harmonics%max_iso_not0
1117 28078 : DO iat = 1, nat
1118 17100 : iatom = atom_list(iat)
1119 : ! *** allocate the radial density for each LM,for each atom ***
1120 :
1121 69736 : ALLOCATE (rho_atom_set(iatom)%rho_rad_h(nspins))
1122 52636 : ALLOCATE (rho_atom_set(iatom)%rho_rad_s(nspins))
1123 52636 : ALLOCATE (rho_atom_set(iatom)%vrho_rad_h(nspins))
1124 52636 : ALLOCATE (rho_atom_set(iatom)%vrho_rad_s(nspins))
1125 :
1126 : ALLOCATE (rho_atom_set(iatom)%cpc_h(nspins), &
1127 : rho_atom_set(iatom)%cpc_s(nspins), &
1128 : rho_atom_set(iatom)%drho_rad_h(nspins), &
1129 : rho_atom_set(iatom)%drho_rad_s(nspins), &
1130 : rho_atom_set(iatom)%rho_rad_h_d(3, nspins), &
1131 358032 : rho_atom_set(iatom)%rho_rad_s_d(3, nspins))
1132 : ALLOCATE (rho_atom_set(iatom)%int_scr_h(nspins), &
1133 88172 : rho_atom_set(iatom)%int_scr_s(nspins))
1134 :
1135 28078 : IF (paw_atom) THEN
1136 31582 : DO ispin = 1, nspins
1137 : ALLOCATE (rho_atom_set(iatom)%cpc_h(ispin)%r_coef(1:nsatbas, 1:nsatbas), &
1138 98508 : rho_atom_set(iatom)%cpc_s(ispin)%r_coef(1:nsatbas, 1:nsatbas))
1139 : ALLOCATE (rho_atom_set(iatom)%int_scr_h(ispin)%r_coef(1:nsatbas, 1:nsatbas), &
1140 82090 : rho_atom_set(iatom)%int_scr_s(ispin)%r_coef(1:nsatbas, 1:nsatbas))
1141 :
1142 7579602 : rho_atom_set(iatom)%cpc_h(ispin)%r_coef = 0.0_dp
1143 7594766 : rho_atom_set(iatom)%cpc_s(ispin)%r_coef = 0.0_dp
1144 : END DO
1145 : END IF
1146 :
1147 : END DO ! iat
1148 :
1149 10978 : num_pe = para_env%num_pe
1150 10978 : mepos = para_env%mepos
1151 10978 : bo = get_limit(nat, num_pe, mepos)
1152 36010 : DO iat = bo(1), bo(2)
1153 8550 : iatom = atom_list(iat)
1154 : ALLOCATE (rho_atom_set(iatom)%ga_Vlocal_gb_h(nspins), &
1155 52636 : rho_atom_set(iatom)%ga_Vlocal_gb_s(nspins))
1156 19528 : IF (paw_atom) THEN
1157 15791 : DO ispin = 1, nspins
1158 : CALL reallocate(rho_atom_set(iatom)%ga_Vlocal_gb_h(ispin)%r_coef, &
1159 8209 : 1, nsotot, 1, nsotot)
1160 : CALL reallocate(rho_atom_set(iatom)%ga_Vlocal_gb_s(ispin)%r_coef, &
1161 8209 : 1, nsotot, 1, nsotot)
1162 :
1163 22556625 : rho_atom_set(iatom)%ga_Vlocal_gb_h(ispin)%r_coef = 0.0_dp
1164 22564207 : rho_atom_set(iatom)%ga_Vlocal_gb_s(ispin)%r_coef = 0.0_dp
1165 : END DO
1166 : END IF
1167 :
1168 : END DO ! iat
1169 :
1170 : END DO
1171 :
1172 5504 : CALL timestop(handle)
1173 :
1174 11008 : END SUBROUTINE allocate_rho_atom_internals
1175 :
1176 : ! **************************************************************************************************
1177 : !> \brief ...
1178 : !> \param rho_atom_set ...
1179 : !> \param iatom ...
1180 : !> \param ispin ...
1181 : !> \param nr ...
1182 : !> \param max_iso_not0 ...
1183 : ! **************************************************************************************************
1184 7975 : SUBROUTINE allocate_rho_atom_rad(rho_atom_set, iatom, ispin, nr, max_iso_not0)
1185 :
1186 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
1187 : INTEGER, INTENT(IN) :: iatom, ispin, nr, max_iso_not0
1188 :
1189 : CHARACTER(len=*), PARAMETER :: routineN = 'allocate_rho_atom_rad'
1190 :
1191 : INTEGER :: handle, j
1192 :
1193 7975 : CALL timeset(routineN, handle)
1194 :
1195 : ALLOCATE (rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef(1:nr, 1:max_iso_not0), &
1196 : rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef(1:nr, 1:max_iso_not0), &
1197 : rho_atom_set(iatom)%vrho_rad_h(ispin)%r_coef(1:nr, 1:max_iso_not0), &
1198 79750 : rho_atom_set(iatom)%vrho_rad_s(ispin)%r_coef(1:nr, 1:max_iso_not0))
1199 :
1200 6124838 : rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef = 0.0_dp
1201 6124838 : rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef = 0.0_dp
1202 6124838 : rho_atom_set(iatom)%vrho_rad_h(ispin)%r_coef = 0.0_dp
1203 6124838 : rho_atom_set(iatom)%vrho_rad_s(ispin)%r_coef = 0.0_dp
1204 :
1205 : ALLOCATE (rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef(nr, max_iso_not0), &
1206 47850 : rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef(nr, max_iso_not0))
1207 6124838 : rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef = 0.0_dp
1208 6124838 : rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef = 0.0_dp
1209 :
1210 31900 : DO j = 1, 3
1211 : ALLOCATE (rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef(nr, max_iso_not0), &
1212 119625 : rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef(nr, max_iso_not0))
1213 18374514 : rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef = 0.0_dp
1214 18382489 : rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef = 0.0_dp
1215 : END DO
1216 :
1217 7975 : CALL timestop(handle)
1218 :
1219 7975 : END SUBROUTINE allocate_rho_atom_rad
1220 :
1221 : ! **************************************************************************************************
1222 : !> \brief ...
1223 : !> \param rho_atom_set ...
1224 : !> \param iatom ...
1225 : !> \param ispin ...
1226 : ! **************************************************************************************************
1227 57042 : SUBROUTINE set2zero_rho_atom_rad(rho_atom_set, iatom, ispin)
1228 :
1229 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
1230 : INTEGER, INTENT(IN) :: iatom, ispin
1231 :
1232 : INTEGER :: j
1233 :
1234 41134926 : rho_atom_set(iatom)%rho_rad_h(ispin)%r_coef = 0.0_dp
1235 41134926 : rho_atom_set(iatom)%rho_rad_s(ispin)%r_coef = 0.0_dp
1236 :
1237 41134926 : rho_atom_set(iatom)%vrho_rad_h(ispin)%r_coef = 0.0_dp
1238 41134926 : rho_atom_set(iatom)%vrho_rad_s(ispin)%r_coef = 0.0_dp
1239 :
1240 41134926 : rho_atom_set(iatom)%drho_rad_h(ispin)%r_coef = 0.0_dp
1241 41134926 : rho_atom_set(iatom)%drho_rad_s(ispin)%r_coef = 0.0_dp
1242 :
1243 228168 : DO j = 1, 3
1244 123404778 : rho_atom_set(iatom)%rho_rad_h_d(j, ispin)%r_coef = 0.0_dp
1245 123461820 : rho_atom_set(iatom)%rho_rad_s_d(j, ispin)%r_coef = 0.0_dp
1246 : END DO
1247 :
1248 57042 : END SUBROUTINE set2zero_rho_atom_rad
1249 :
1250 : END MODULE qs_rho_atom_methods
|
