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