Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : ! **************************************************************************************************
8 : MODULE qs_rho0_ggrid
9 : USE atomic_kind_types, ONLY: atomic_kind_type,&
10 : get_atomic_kind
11 : USE basis_set_types, ONLY: get_gto_basis_set,&
12 : gto_basis_set_type
13 : USE cell_types, ONLY: cell_type,&
14 : pbc
15 : USE cp_control_types, ONLY: dft_control_type
16 : USE gaussian_gridlevels, ONLY: gaussian_gridlevel
17 : USE grid_api, ONLY: GRID_FUNC_AB,&
18 : collocate_pgf_product
19 : USE kinds, ONLY: dp
20 : USE memory_utilities, ONLY: reallocate
21 : USE message_passing, ONLY: mp_para_env_type
22 : USE orbital_pointers, ONLY: indco,&
23 : nco,&
24 : ncoset,&
25 : nso,&
26 : nsoset
27 : USE orbital_transformation_matrices, ONLY: orbtramat
28 : USE particle_types, ONLY: particle_type
29 : USE pw_env_types, ONLY: pw_env_get,&
30 : pw_env_type
31 : USE pw_methods, ONLY: pw_axpy,&
32 : pw_copy,&
33 : pw_integrate_function,&
34 : pw_scale,&
35 : pw_transfer,&
36 : pw_zero
37 : USE pw_pool_types, ONLY: pw_pool_p_type,&
38 : pw_pool_type
39 : USE pw_types, ONLY: pw_c1d_gs_type,&
40 : pw_r3d_rs_type
41 : USE qs_environment_types, ONLY: get_qs_env,&
42 : qs_environment_type
43 : USE qs_force_types, ONLY: qs_force_type
44 : USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
45 : harmonics_atom_type
46 : USE qs_integrate_potential, ONLY: integrate_pgf_product
47 : USE qs_kind_types, ONLY: get_qs_kind,&
48 : qs_kind_type
49 : USE qs_local_rho_types, ONLY: get_local_rho,&
50 : local_rho_type
51 : USE qs_rho0_types, ONLY: get_rho0_mpole,&
52 : rho0_mpole_type
53 : USE qs_rho_atom_types, ONLY: get_rho_atom,&
54 : rho_atom_coeff,&
55 : rho_atom_type
56 : USE realspace_grid_types, ONLY: realspace_grid_desc_p_type,&
57 : realspace_grid_desc_type,&
58 : realspace_grid_type,&
59 : rs_grid_create,&
60 : rs_grid_release,&
61 : rs_grid_zero,&
62 : transfer_pw2rs,&
63 : transfer_rs2pw
64 : USE util, ONLY: get_limit
65 : USE virial_types, ONLY: virial_type
66 : #include "./base/base_uses.f90"
67 :
68 : IMPLICIT NONE
69 :
70 : PRIVATE
71 :
72 : ! Global parameters (only in this module)
73 :
74 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho0_ggrid'
75 :
76 : ! Public subroutines
77 :
78 : PUBLIC :: put_rho0_on_grid, rho0_s_grid_create, integrate_vhg0_rspace
79 :
80 : CONTAINS
81 :
82 : ! **************************************************************************************************
83 : !> \brief ...
84 : !> \param qs_env ...
85 : !> \param rho0 ...
86 : !> \param tot_rs_int ...
87 : !> \param my_pools ...
88 : !> \param my_rs_grids ...
89 : !> \param my_rs_descs ...
90 : ! **************************************************************************************************
91 16010 : SUBROUTINE put_rho0_on_grid(qs_env, rho0, tot_rs_int, my_pools, my_rs_grids, my_rs_descs)
92 :
93 : TYPE(qs_environment_type), POINTER :: qs_env
94 : TYPE(rho0_mpole_type), POINTER :: rho0
95 : REAL(KIND=dp), INTENT(OUT) :: tot_rs_int
96 : TYPE(pw_pool_p_type), DIMENSION(:), OPTIONAL, &
97 : POINTER :: my_pools
98 : TYPE(realspace_grid_type), DIMENSION(:), &
99 : OPTIONAL, POINTER :: my_rs_grids
100 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
101 : OPTIONAL, POINTER :: my_rs_descs
102 :
103 : CHARACTER(LEN=*), PARAMETER :: routineN = 'put_rho0_on_grid'
104 :
105 : INTEGER :: auxbas_grid, handle, iat, iatom, igrid, &
106 : ikind, ithread, j, l0_ikind, lmax0, &
107 : nat, nch_ik, nch_max, npme
108 16010 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
109 : LOGICAL :: paw_atom
110 : REAL(KIND=dp) :: eps_rho_rspace, rpgf0, zet0
111 : REAL(KIND=dp), DIMENSION(3) :: ra
112 16010 : REAL(KIND=dp), DIMENSION(:), POINTER :: Qlm_c
113 16010 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
114 16010 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
115 : TYPE(cell_type), POINTER :: cell
116 : TYPE(dft_control_type), POINTER :: dft_control
117 : TYPE(mp_para_env_type), POINTER :: para_env
118 16010 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
119 : TYPE(pw_c1d_gs_type) :: coeff_gspace
120 : TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs
121 : TYPE(pw_env_type), POINTER :: pw_env
122 16010 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
123 : TYPE(pw_pool_type), POINTER :: pw_pool
124 : TYPE(pw_r3d_rs_type) :: coeff_rspace, rho0_r_tmp
125 : TYPE(pw_r3d_rs_type), POINTER :: rho0_s_rs
126 16010 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
127 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
128 16010 : POINTER :: descs
129 : TYPE(realspace_grid_desc_type), POINTER :: desc
130 16010 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: grids
131 : TYPE(realspace_grid_type), POINTER :: rs_grid
132 :
133 16010 : CALL timeset(routineN, handle)
134 :
135 16010 : NULLIFY (atomic_kind_set, qs_kind_set, cores, pab, Qlm_c)
136 :
137 16010 : NULLIFY (dft_control, pw_env, particle_set, para_env, cell, rho0_s_gs, rho0_s_rs)
138 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
139 : particle_set=particle_set, &
140 : atomic_kind_set=atomic_kind_set, &
141 : qs_kind_set=qs_kind_set, &
142 : para_env=para_env, &
143 16010 : pw_env=pw_env, cell=cell)
144 16010 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
145 :
146 16010 : NULLIFY (descs, pw_pools)
147 16010 : CALL pw_env_get(pw_env=pw_env, rs_descs=descs, rs_grids=grids, pw_pools=pw_pools)
148 16010 : auxbas_grid = pw_env%auxbas_grid
149 :
150 16010 : NULLIFY (rho0_s_gs, rho0_s_rs)
151 : CALL get_rho0_mpole(rho0_mpole=rho0, lmax_0=lmax0, &
152 : zet0_h=zet0, igrid_zet0_s=igrid, &
153 : rho0_s_gs=rho0_s_gs, &
154 16010 : rho0_s_rs=rho0_s_rs)
155 :
156 : ! *** set up the rs grid at level igrid
157 16010 : NULLIFY (rs_grid, desc, pw_pool)
158 : ! IF present, overwrite qs grid for new pool
159 16010 : IF (PRESENT(my_pools)) THEN
160 972 : desc => my_rs_descs(igrid)%rs_desc
161 972 : rs_grid => my_rs_grids(igrid)
162 972 : pw_pool => my_pools(igrid)%pool
163 : ELSE
164 15038 : desc => descs(igrid)%rs_desc
165 15038 : rs_grid => grids(igrid)
166 15038 : pw_pool => pw_pools(igrid)%pool
167 : END IF
168 :
169 16010 : CPASSERT(ASSOCIATED(desc))
170 16010 : CPASSERT(ASSOCIATED(pw_pool))
171 :
172 16010 : IF (igrid /= auxbas_grid) THEN
173 12 : CALL pw_pool%create_pw(coeff_rspace)
174 12 : CALL pw_pool%create_pw(coeff_gspace)
175 : END IF
176 16010 : CALL rs_grid_zero(rs_grid)
177 :
178 16010 : nch_max = ncoset(lmax0)
179 :
180 48030 : ALLOCATE (pab(nch_max, 1))
181 :
182 48580 : DO ikind = 1, SIZE(atomic_kind_set)
183 32570 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
184 32570 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
185 :
186 32570 : IF (.NOT. paw_atom .AND. dft_control%qs_control%gapw_control%nopaw_as_gpw) CYCLE
187 :
188 : CALL get_rho0_mpole(rho0_mpole=rho0, ikind=ikind, l0_ikind=l0_ikind, &
189 30982 : rpgf0_s=rpgf0)
190 :
191 30982 : nch_ik = ncoset(l0_ikind)
192 387820 : pab = 0.0_dp
193 :
194 30982 : CALL reallocate(cores, 1, nat)
195 30982 : npme = 0
196 79972 : cores = 0
197 :
198 79972 : DO iat = 1, nat
199 48990 : iatom = atom_list(iat)
200 48990 : ra(:) = pbc(particle_set(iatom)%r, cell)
201 79972 : IF (rs_grid%desc%parallel .AND. .NOT. rs_grid%desc%distributed) THEN
202 : ! replicated realspace grid, split the atoms up between procs
203 48990 : IF (MODULO(nat, rs_grid%desc%group_size) == rs_grid%desc%my_pos) THEN
204 24495 : npme = npme + 1
205 24495 : cores(npme) = iat
206 : END IF
207 : ELSE
208 0 : npme = npme + 1
209 0 : cores(npme) = iat
210 : END IF
211 :
212 : END DO
213 :
214 : ithread = 0
215 135039 : DO j = 1, npme
216 :
217 24495 : iat = cores(j)
218 24495 : iatom = atom_list(iat)
219 :
220 24495 : CALL get_rho0_mpole(rho0_mpole=rho0, iat=iatom, Qlm_car=Qlm_c)
221 :
222 465586 : pab(1:nch_ik, 1) = Qlm_c(1:nch_ik)
223 :
224 24495 : ra(:) = pbc(particle_set(iatom)%r, cell)
225 :
226 : CALL collocate_pgf_product( &
227 : l0_ikind, zet0, 0, 0, 0.0_dp, 0, &
228 : ra, (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, &
229 : rs_grid, ga_gb_function=GRID_FUNC_AB, radius=rpgf0, &
230 57065 : use_subpatch=.TRUE., subpatch_pattern=0)
231 :
232 : END DO ! j
233 :
234 : END DO ! ikind
235 :
236 16010 : IF (ASSOCIATED(cores)) THEN
237 15866 : DEALLOCATE (cores)
238 : END IF
239 :
240 16010 : DEALLOCATE (pab)
241 :
242 16010 : IF (igrid /= auxbas_grid) THEN
243 12 : CALL transfer_rs2pw(rs_grid, coeff_rspace)
244 12 : CALL pw_zero(rho0_s_gs)
245 12 : CALL pw_transfer(coeff_rspace, coeff_gspace)
246 12 : CALL pw_axpy(coeff_gspace, rho0_s_gs)
247 :
248 12 : tot_rs_int = pw_integrate_function(coeff_rspace, isign=-1)
249 :
250 12 : CALL pw_pool%give_back_pw(coeff_rspace)
251 12 : CALL pw_pool%give_back_pw(coeff_gspace)
252 : ELSE
253 :
254 15998 : CALL pw_pool%create_pw(rho0_r_tmp)
255 :
256 15998 : CALL transfer_rs2pw(rs_grid, rho0_r_tmp)
257 :
258 15998 : tot_rs_int = pw_integrate_function(rho0_r_tmp, isign=-1)
259 :
260 15998 : CALL pw_transfer(rho0_r_tmp, rho0_s_rs)
261 15998 : CALL pw_pool%give_back_pw(rho0_r_tmp)
262 :
263 15998 : CALL pw_zero(rho0_s_gs)
264 15998 : CALL pw_transfer(rho0_s_rs, rho0_s_gs)
265 : END IF
266 16010 : CALL timestop(handle)
267 :
268 32020 : END SUBROUTINE put_rho0_on_grid
269 :
270 : ! **************************************************************************************************
271 : !> \brief ...
272 : !> \param pw_env ...
273 : !> \param rho0_mpole ...
274 : ! **************************************************************************************************
275 2616 : SUBROUTINE rho0_s_grid_create(pw_env, rho0_mpole)
276 :
277 : TYPE(pw_env_type), POINTER :: pw_env
278 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
279 :
280 : CHARACTER(len=*), PARAMETER :: routineN = 'rho0_s_grid_create'
281 :
282 : INTEGER :: handle
283 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
284 :
285 2616 : CALL timeset(routineN, handle)
286 :
287 2616 : CPASSERT(ASSOCIATED(pw_env))
288 :
289 2616 : NULLIFY (auxbas_pw_pool)
290 2616 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
291 2616 : CPASSERT(ASSOCIATED(auxbas_pw_pool))
292 :
293 : ! reallocate rho0 on the global grid in real and reciprocal space
294 2616 : CPASSERT(ASSOCIATED(rho0_mpole))
295 :
296 : ! rho0 density in real space
297 2616 : IF (ASSOCIATED(rho0_mpole%rho0_s_rs)) THEN
298 814 : CALL rho0_mpole%rho0_s_rs%release()
299 : ELSE
300 1802 : ALLOCATE (rho0_mpole%rho0_s_rs)
301 : END IF
302 2616 : CALL auxbas_pw_pool%create_pw(rho0_mpole%rho0_s_rs)
303 :
304 : ! rho0 density in reciprocal space
305 2616 : IF (ASSOCIATED(rho0_mpole%rho0_s_gs)) THEN
306 814 : CALL rho0_mpole%rho0_s_gs%release()
307 : ELSE
308 1802 : ALLOCATE (rho0_mpole%rho0_s_gs)
309 : END IF
310 2616 : CALL auxbas_pw_pool%create_pw(rho0_mpole%rho0_s_gs)
311 :
312 : ! Find the grid level suitable for rho0_soft
313 2616 : rho0_mpole%igrid_zet0_s = gaussian_gridlevel(pw_env%gridlevel_info, 2.0_dp*rho0_mpole%zet0_h)
314 :
315 2616 : CALL timestop(handle)
316 :
317 2616 : END SUBROUTINE rho0_s_grid_create
318 :
319 : ! **************************************************************************************************
320 : !> \brief ...
321 : !> \param qs_env ...
322 : !> \param v_rspace ...
323 : !> \param para_env ...
324 : !> \param calculate_forces ...
325 : !> \param local_rho_set ...
326 : !> \param local_rho_set_2nd ...
327 : !> \param atener ...
328 : !> \param kforce ...
329 : !> \param my_pools ...
330 : !> \param my_rs_descs ...
331 : ! **************************************************************************************************
332 14922 : SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set, &
333 14922 : local_rho_set_2nd, atener, kforce, my_pools, my_rs_descs)
334 :
335 : TYPE(qs_environment_type), POINTER :: qs_env
336 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
337 : TYPE(mp_para_env_type), POINTER :: para_env
338 : LOGICAL, INTENT(IN) :: calculate_forces
339 : TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set, local_rho_set_2nd
340 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atener
341 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: kforce
342 : TYPE(pw_pool_p_type), DIMENSION(:), OPTIONAL, &
343 : POINTER :: my_pools
344 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
345 : OPTIONAL, POINTER :: my_rs_descs
346 :
347 : CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_vhg0_rspace'
348 :
349 : INTEGER :: auxbas_grid, bo(2), handle, i, iat, iatom, ic, icg, ico, ig1, ig2, igrid, ii, &
350 : ikind, ipgf1, ipgf2, is, iset1, iset2, iso, iso1, iso2, ispin, j, l0_ikind, llmax, lmax0, &
351 : lshell, lx, ly, lz, m1, m2, max_iso_not0_local, max_s_harm, maxl, maxso, mepos, n1, n2, &
352 : nat, nch_ik, nch_max, ncurr, nset, nsotot, nspins, num_pe
353 14922 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
354 14922 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
355 14922 : INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf
356 : LOGICAL :: grid_distributed, paw_atom, use_virial
357 : REAL(KIND=dp) :: eps_rho_rspace, force_tmp(3), fscale, &
358 : ra(3), rpgf0, zet0
359 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
360 14922 : REAL(KIND=dp), DIMENSION(:), POINTER :: hab_sph, norm_l, Qlm
361 14922 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, hdab_sph, intloc, pab
362 14922 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: a_hdab_sph, hdab, Qlm_gg
363 14922 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: a_hdab
364 14922 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
365 : TYPE(cell_type), POINTER :: cell
366 : TYPE(dft_control_type), POINTER :: dft_control
367 : TYPE(gto_basis_set_type), POINTER :: basis_1c_set
368 : TYPE(harmonics_atom_type), POINTER :: harmonics
369 14922 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
370 : TYPE(pw_c1d_gs_type) :: coeff_gaux, coeff_gspace
371 : TYPE(pw_env_type), POINTER :: pw_env
372 14922 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
373 : TYPE(pw_pool_type), POINTER :: pw_aux, pw_pool
374 : TYPE(pw_r3d_rs_type) :: coeff_raux, coeff_rspace
375 14922 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
376 14922 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
377 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
378 14922 : POINTER :: rs_descs
379 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
380 298440 : TYPE(realspace_grid_type) :: rs_v
381 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
382 14922 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
383 14922 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
384 : TYPE(rho_atom_type), POINTER :: rho_atom
385 : TYPE(virial_type), POINTER :: virial
386 :
387 14922 : CALL timeset(routineN, handle)
388 :
389 : ! In case of the external density computed forces probably also
390 : ! need to be stored outside qs_env. We can then remove the
391 : ! attribute 'OPTIONAL' from the argument 'local_rho_set'.
392 :
393 : ! CPASSERT(.NOT. (calculate_forces .AND. PRESENT(local_rho_set)))
394 : ! IF (calculate_forces .AND. PRESENT(local_rho_set)) THEN
395 : ! CPWARN("Forces and External Density!")
396 : ! END IF
397 :
398 14922 : NULLIFY (atomic_kind_set, qs_kind_set, dft_control, particle_set)
399 14922 : NULLIFY (cell, force, pw_env, rho0_mpole, rho_atom_set)
400 :
401 : CALL get_qs_env(qs_env=qs_env, &
402 : atomic_kind_set=atomic_kind_set, &
403 : qs_kind_set=qs_kind_set, &
404 : cell=cell, &
405 : dft_control=dft_control, &
406 : force=force, pw_env=pw_env, &
407 : rho0_mpole=rho0_mpole, &
408 : rho_atom_set=rho_atom_set, &
409 : particle_set=particle_set, &
410 14922 : virial=virial)
411 :
412 14922 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
413 :
414 14922 : nspins = dft_control%nspins
415 :
416 : ! The aim of the following code was to return immediately if the subroutine
417 : ! was called for triplet excited states in spin-restricted case. This check
418 : ! is also performed before invocation of this subroutine. It should be save
419 : ! to remove the optional argument 'do_triplet' from the subroutine interface.
420 : !my_tddft = PRESENT(local_rho_set)
421 : !IF (my_tddft) THEN
422 : ! IF (PRESENT(do_triplet)) THEN
423 : ! IF (nspins == 1 .AND. do_triplet) RETURN
424 : ! ELSE
425 : ! IF (nspins == 1 .AND. dft_control%tddfpt_control%res_etype /= tddfpt_singlet) RETURN
426 : ! END IF
427 : !END IF
428 :
429 14922 : IF (PRESENT(local_rho_set)) &
430 1910 : CALL get_local_rho(local_rho_set, rho0_mpole=rho0_mpole, rho_atom_set=rho_atom_set)
431 : ! Q from rho0_mpole of local_rho_set
432 : ! for TDDFT forces we need mixed potential / integral space
433 : ! potential stored on local_rho_set_2nd
434 14922 : IF (PRESENT(local_rho_set_2nd)) THEN
435 128 : CALL get_local_rho(local_rho_set_2nd, rho_atom_set=rho_atom_set)
436 : END IF
437 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, lmax_0=lmax0, &
438 : zet0_h=zet0, igrid_zet0_s=igrid, &
439 14922 : norm_g0l_h=norm_l)
440 :
441 : ! Setup of the potential on the multigrids
442 14922 : NULLIFY (rs_descs, pw_pools)
443 14922 : CPASSERT(ASSOCIATED(pw_env))
444 14922 : CALL pw_env_get(pw_env, rs_descs=rs_descs, pw_pools=pw_pools)
445 :
446 : ! Assign from pw_env
447 14922 : auxbas_grid = pw_env%auxbas_grid
448 :
449 : ! IF present, overwrite qs grid for new pool
450 : ! Get the potential on the right grid
451 14922 : IF (PRESENT(my_pools)) THEN
452 50 : rs_desc => my_rs_descs(igrid)%rs_desc
453 50 : pw_pool => my_pools(igrid)%pool
454 : ELSE
455 14872 : rs_desc => rs_descs(igrid)%rs_desc
456 14872 : pw_pool => pw_pools(igrid)%pool
457 : END IF
458 :
459 14922 : CALL pw_pool%create_pw(coeff_gspace)
460 14922 : CALL pw_pool%create_pw(coeff_rspace)
461 :
462 14922 : IF (igrid /= auxbas_grid) THEN
463 12 : pw_aux => pw_pools(auxbas_grid)%pool
464 12 : CALL pw_aux%create_pw(coeff_gaux)
465 12 : CALL pw_transfer(v_rspace, coeff_gaux)
466 12 : CALL pw_copy(coeff_gaux, coeff_gspace)
467 12 : CALL pw_transfer(coeff_gspace, coeff_rspace)
468 12 : CALL pw_aux%give_back_pw(coeff_gaux)
469 12 : CALL pw_aux%create_pw(coeff_raux)
470 12 : fscale = coeff_rspace%pw_grid%dvol/coeff_raux%pw_grid%dvol
471 12 : CALL pw_scale(coeff_rspace, fscale)
472 12 : CALL pw_aux%give_back_pw(coeff_raux)
473 : ELSE
474 :
475 14910 : IF (coeff_gspace%pw_grid%spherical) THEN
476 0 : CALL pw_transfer(v_rspace, coeff_gspace)
477 0 : CALL pw_transfer(coeff_gspace, coeff_rspace)
478 : ELSE
479 14910 : CALL pw_copy(v_rspace, coeff_rspace)
480 : END IF
481 : END IF
482 14922 : CALL pw_pool%give_back_pw(coeff_gspace)
483 :
484 : ! Setup the rs grid at level igrid
485 14922 : CALL rs_grid_create(rs_v, rs_desc)
486 14922 : CALL rs_grid_zero(rs_v)
487 14922 : CALL transfer_pw2rs(rs_v, coeff_rspace)
488 :
489 14922 : CALL pw_pool%give_back_pw(coeff_rspace)
490 :
491 : ! Now the potential is on the right grid => integration
492 :
493 14922 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
494 :
495 : ! Allocate work storage
496 :
497 14922 : NULLIFY (hab, hab_sph, hdab, hdab_sph, pab, a_hdab, a_hdab_sph)
498 14922 : nch_max = ncoset(lmax0)
499 14922 : CALL reallocate(hab, 1, nch_max, 1, 1)
500 14922 : CALL reallocate(hab_sph, 1, nch_max)
501 14922 : CALL reallocate(hdab, 1, 3, 1, nch_max, 1, 1)
502 14922 : CALL reallocate(hdab_sph, 1, 3, 1, nch_max)
503 14922 : CALL reallocate(a_hdab, 1, 3, 1, 3, 1, nch_max, 1, 1)
504 14922 : CALL reallocate(a_hdab_sph, 1, 3, 1, 3, 1, nch_max)
505 14922 : CALL reallocate(pab, 1, nch_max, 1, 1)
506 :
507 14922 : ncurr = -1
508 :
509 14922 : grid_distributed = rs_v%desc%distributed
510 :
511 14922 : fscale = 1.0_dp
512 14922 : IF (PRESENT(kforce)) THEN
513 40 : fscale = kforce
514 : END IF
515 :
516 44954 : DO ikind = 1, SIZE(atomic_kind_set, 1)
517 30032 : NULLIFY (basis_1c_set, atom_list, harmonics)
518 30032 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
519 : CALL get_qs_kind(qs_kind_set(ikind), &
520 : basis_set=basis_1c_set, basis_type="GAPW_1C", &
521 : paw_atom=paw_atom, &
522 30032 : harmonics=harmonics)
523 :
524 30032 : IF (.NOT. paw_atom) CYCLE
525 :
526 28492 : NULLIFY (Qlm_gg, lmax, npgf)
527 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, &
528 : l0_ikind=l0_ikind, Qlm_gg=Qlm_gg, & ! Qs different
529 28492 : rpgf0_s=rpgf0)
530 :
531 : CALL get_gto_basis_set(gto_basis_set=basis_1c_set, &
532 : lmax=lmax, lmin=lmin, &
533 : maxso=maxso, maxl=maxl, &
534 28492 : nset=nset, npgf=npgf)
535 :
536 28492 : nsotot = maxso*nset
537 113968 : ALLOCATE (intloc(nsotot, nsotot))
538 :
539 : ! Initialize the local KS integrals
540 :
541 28492 : nch_ik = ncoset(l0_ikind)
542 353740 : pab = 1.0_dp
543 28492 : max_s_harm = harmonics%max_s_harm
544 28492 : llmax = harmonics%llmax
545 :
546 170952 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
547 :
548 28492 : num_pe = para_env%num_pe
549 28492 : mepos = para_env%mepos
550 85476 : DO j = 0, num_pe - 1
551 56984 : bo = get_limit(nat, num_pe, j)
552 56984 : IF (.NOT. grid_distributed .AND. j /= mepos) CYCLE
553 :
554 79702 : DO iat = bo(1), bo(2)
555 22718 : iatom = atom_list(iat)
556 22718 : ra(:) = pbc(particle_set(iatom)%r, cell)
557 :
558 22718 : NULLIFY (Qlm)
559 22718 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, iat=iatom, Qlm_tot=Qlm)
560 :
561 278888 : hab = 0.0_dp
562 979244 : hdab = 0.0_dp
563 65185276 : intloc = 0._dp
564 22718 : IF (use_virial) THEN
565 274 : my_virial_a = 0.0_dp
566 274 : my_virial_b = 0.0_dp
567 36168 : a_hdab = 0.0_dp
568 : END IF
569 :
570 : CALL integrate_pgf_product( &
571 : l0_ikind, zet0, 0, 0, 0.0_dp, 0, &
572 : ra, (/0.0_dp, 0.0_dp, 0.0_dp/), rs_v, &
573 : hab, pab, o1=0, o2=0, &
574 : radius=rpgf0, &
575 : calculate_forces=calculate_forces, &
576 : use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b, &
577 22718 : hdab=hdab, a_hdab=a_hdab, use_subpatch=.TRUE., subpatch_pattern=0)
578 :
579 : ! Convert from cartesian to spherical
580 82148 : DO lshell = 0, l0_ikind
581 253506 : DO is = 1, nso(lshell)
582 171358 : iso = is + nsoset(lshell - 1)
583 171358 : hab_sph(iso) = 0.0_dp
584 685432 : hdab_sph(1:3, iso) = 0.0_dp
585 2227654 : a_hdab_sph(1:3, 1:3, iso) = 0.0_dp
586 1006574 : DO ic = 1, nco(lshell)
587 775786 : ico = ic + ncoset(lshell - 1)
588 775786 : lx = indco(1, ico)
589 775786 : ly = indco(2, ico)
590 775786 : lz = indco(3, ico)
591 : hab_sph(iso) = hab_sph(iso) + &
592 : norm_l(lshell)* &
593 : orbtramat(lshell)%slm(is, ic)* &
594 775786 : hab(ico, 1)
595 775786 : IF (calculate_forces) THEN
596 : hdab_sph(1:3, iso) = hdab_sph(1:3, iso) + &
597 : norm_l(lshell)* &
598 : orbtramat(lshell)%slm(is, ic)* &
599 293184 : hdab(1:3, ico, 1)
600 : END IF
601 947144 : IF (use_virial) THEN
602 43840 : DO ii = 1, 3
603 142480 : DO i = 1, 3
604 : a_hdab_sph(i, ii, iso) = a_hdab_sph(i, ii, iso) + &
605 : norm_l(lshell)* &
606 : orbtramat(lshell)%slm(is, ic)* &
607 131520 : a_hdab(i, ii, ico, 1)
608 : END DO
609 : END DO
610 : END IF
611 :
612 : END DO ! ic
613 : END DO ! is
614 : END DO ! lshell
615 :
616 : m1 = 0
617 78771 : DO iset1 = 1, nset
618 :
619 : m2 = 0
620 237938 : DO iset2 = 1, nset
621 : CALL get_none0_cg_list(harmonics%my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
622 181885 : max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
623 181885 : n1 = nsoset(lmax(iset1))
624 628269 : DO ipgf1 = 1, npgf(iset1)
625 446384 : n2 = nsoset(lmax(iset2))
626 1925312 : DO ipgf2 = 1, npgf(iset2)
627 :
628 8212123 : DO iso = 1, MIN(nsoset(l0_ikind), max_iso_not0_local)
629 18262291 : DO icg = 1, cg_n_list(iso)
630 10496552 : iso1 = cg_list(1, icg, iso)
631 10496552 : iso2 = cg_list(2, icg, iso)
632 :
633 10496552 : ig1 = iso1 + n1*(ipgf1 - 1) + m1
634 10496552 : ig2 = iso2 + n2*(ipgf2 - 1) + m2
635 :
636 16965248 : intloc(ig1, ig2) = intloc(ig1, ig2) + Qlm_gg(ig1, ig2, iso)*hab_sph(iso) ! potential times Q
637 :
638 : END DO ! icg
639 : END DO ! iso
640 :
641 : END DO ! ipgf2
642 : END DO ! ipgf1
643 419823 : m2 = m2 + maxso
644 : END DO ! iset2
645 78771 : m1 = m1 + maxso
646 : END DO ! iset1
647 :
648 22718 : IF (grid_distributed) THEN
649 : ! Sum result over all processors
650 0 : CALL para_env%sum(intloc)
651 : END IF
652 :
653 22718 : IF (j == mepos) THEN
654 22718 : rho_atom => rho_atom_set(iatom)
655 22718 : CALL get_rho_atom(rho_atom=rho_atom, ga_Vlocal_gb_h=int_local_h, ga_Vlocal_gb_s=int_local_s)
656 49381 : DO ispin = 1, nspins
657 147977518 : int_local_h(ispin)%r_coef = int_local_h(ispin)%r_coef + intloc
658 148000236 : int_local_s(ispin)%r_coef = int_local_s(ispin)%r_coef + intloc
659 : END DO
660 : END IF
661 :
662 22718 : IF (PRESENT(atener)) THEN
663 178 : DO iso = 1, nsoset(l0_ikind)
664 178 : atener(iatom) = atener(iatom) + 0.5_dp*Qlm(iso)*hab_sph(iso)
665 : END DO
666 : END IF
667 :
668 22718 : IF (calculate_forces) THEN
669 920 : force_tmp(1:3) = 0.0_dp
670 8752 : DO iso = 1, nsoset(l0_ikind)
671 7832 : force_tmp(1) = force_tmp(1) + Qlm(iso)*hdab_sph(1, iso) ! Q here is from local_rho_set
672 7832 : force_tmp(2) = force_tmp(2) + Qlm(iso)*hdab_sph(2, iso)
673 8752 : force_tmp(3) = force_tmp(3) + Qlm(iso)*hdab_sph(3, iso)
674 : END DO
675 3680 : force(ikind)%g0s_Vh_elec(1:3, iat) = force(ikind)%g0s_Vh_elec(1:3, iat) + fscale*force_tmp(1:3)
676 : END IF
677 79702 : IF (use_virial) THEN
678 274 : my_virial_a = 0.0_dp
679 2740 : DO iso = 1, nsoset(l0_ikind)
680 10138 : DO ii = 1, 3
681 32058 : DO i = 1, 3
682 : ! Q from local_rho_set
683 22194 : virial%pv_gapw(i, ii) = virial%pv_gapw(i, ii) + fscale*Qlm(iso)*a_hdab_sph(i, ii, iso)
684 29592 : virial%pv_virial(i, ii) = virial%pv_virial(i, ii) + fscale*Qlm(iso)*a_hdab_sph(i, ii, iso)
685 : END DO
686 : END DO
687 : END DO
688 : END IF
689 :
690 : END DO
691 : END DO
692 :
693 28492 : DEALLOCATE (intloc)
694 103478 : DEALLOCATE (cg_list, cg_n_list)
695 :
696 : END DO ! ikind
697 :
698 14922 : CALL rs_grid_release(rs_v)
699 :
700 14922 : DEALLOCATE (hab, hdab, hab_sph, hdab_sph, pab, a_hdab, a_hdab_sph)
701 :
702 14922 : CALL timestop(handle)
703 :
704 29844 : END SUBROUTINE integrate_vhg0_rspace
705 :
706 : END MODULE qs_rho0_ggrid
|
