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