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