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 : ! **************************************************************************************************
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 : ! **************************************************************************************************
88 16018 : SUBROUTINE put_rho0_on_grid(qs_env, rho0, tot_rs_int)
89 :
90 : TYPE(qs_environment_type), POINTER :: qs_env
91 : TYPE(rho0_mpole_type), POINTER :: rho0
92 : REAL(KIND=dp), INTENT(OUT) :: tot_rs_int
93 :
94 : CHARACTER(LEN=*), PARAMETER :: routineN = 'put_rho0_on_grid'
95 :
96 : INTEGER :: auxbas_grid, handle, iat, iatom, igrid, &
97 : ikind, ithread, j, l0_ikind, lmax0, &
98 : nat, nch_ik, nch_max, npme
99 16018 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
100 : LOGICAL :: paw_atom
101 : REAL(KIND=dp) :: eps_rho_rspace, rpgf0, zet0
102 : REAL(KIND=dp), DIMENSION(3) :: ra
103 16018 : REAL(KIND=dp), DIMENSION(:), POINTER :: Qlm_c
104 16018 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
105 16018 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
106 : TYPE(cell_type), POINTER :: cell
107 : TYPE(dft_control_type), POINTER :: dft_control
108 : TYPE(mp_para_env_type), POINTER :: para_env
109 16018 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
110 : TYPE(pw_c1d_gs_type) :: coeff_gspace
111 : TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs
112 : TYPE(pw_env_type), POINTER :: pw_env
113 16018 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
114 : TYPE(pw_pool_type), POINTER :: pw_pool
115 : TYPE(pw_r3d_rs_type) :: coeff_rspace, rho0_r_tmp
116 : TYPE(pw_r3d_rs_type), POINTER :: rho0_s_rs
117 16018 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
118 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
119 16018 : POINTER :: descs
120 : TYPE(realspace_grid_desc_type), POINTER :: desc
121 16018 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: grids
122 : TYPE(realspace_grid_type), POINTER :: rs_grid
123 :
124 16018 : CALL timeset(routineN, handle)
125 :
126 16018 : NULLIFY (atomic_kind_set, qs_kind_set, cores, pab, Qlm_c)
127 :
128 16018 : NULLIFY (dft_control, pw_env, particle_set, para_env, cell, rho0_s_gs, rho0_s_rs)
129 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
130 : particle_set=particle_set, &
131 : atomic_kind_set=atomic_kind_set, &
132 : qs_kind_set=qs_kind_set, &
133 : para_env=para_env, &
134 16018 : pw_env=pw_env, cell=cell)
135 16018 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
136 :
137 16018 : NULLIFY (descs, pw_pools)
138 16018 : CALL pw_env_get(pw_env=pw_env, rs_descs=descs, rs_grids=grids, pw_pools=pw_pools)
139 16018 : auxbas_grid = pw_env%auxbas_grid
140 :
141 16018 : NULLIFY (rho0_s_gs, rho0_s_rs)
142 : CALL get_rho0_mpole(rho0_mpole=rho0, lmax_0=lmax0, &
143 : zet0_h=zet0, igrid_zet0_s=igrid, &
144 : rho0_s_gs=rho0_s_gs, &
145 16018 : rho0_s_rs=rho0_s_rs)
146 :
147 : ! *** set up the rs grid at level igrid
148 16018 : NULLIFY (rs_grid, desc, pw_pool)
149 16018 : desc => descs(igrid)%rs_desc
150 16018 : rs_grid => grids(igrid)
151 16018 : pw_pool => pw_pools(igrid)%pool
152 :
153 16018 : CPASSERT(ASSOCIATED(desc))
154 16018 : CPASSERT(ASSOCIATED(pw_pool))
155 :
156 16018 : IF (igrid /= auxbas_grid) THEN
157 12 : CALL pw_pool%create_pw(coeff_rspace)
158 12 : CALL pw_pool%create_pw(coeff_gspace)
159 : END IF
160 16018 : CALL rs_grid_zero(rs_grid)
161 :
162 16018 : nch_max = ncoset(lmax0)
163 :
164 48054 : ALLOCATE (pab(nch_max, 1))
165 :
166 48504 : DO ikind = 1, SIZE(atomic_kind_set)
167 32486 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
168 32486 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
169 :
170 32486 : IF (.NOT. paw_atom .AND. dft_control%qs_control%gapw_control%nopaw_as_gpw) CYCLE
171 :
172 : CALL get_rho0_mpole(rho0_mpole=rho0, ikind=ikind, l0_ikind=l0_ikind, &
173 30940 : rpgf0_s=rpgf0)
174 :
175 30940 : nch_ik = ncoset(l0_ikind)
176 386110 : pab = 0.0_dp
177 :
178 30940 : CALL reallocate(cores, 1, nat)
179 30940 : npme = 0
180 79932 : cores = 0
181 :
182 79932 : DO iat = 1, nat
183 48992 : iatom = atom_list(iat)
184 48992 : ra(:) = pbc(particle_set(iatom)%r, cell)
185 79932 : IF (rs_grid%desc%parallel .AND. .NOT. rs_grid%desc%distributed) THEN
186 : ! replicated realspace grid, split the atoms up between procs
187 48992 : IF (MODULO(nat, rs_grid%desc%group_size) == rs_grid%desc%my_pos) THEN
188 24496 : npme = npme + 1
189 24496 : cores(npme) = iat
190 : END IF
191 : ELSE
192 0 : npme = npme + 1
193 0 : cores(npme) = iat
194 : END IF
195 :
196 : END DO
197 :
198 : ithread = 0
199 134880 : DO j = 1, npme
200 :
201 24496 : iat = cores(j)
202 24496 : iatom = atom_list(iat)
203 :
204 24496 : CALL get_rho0_mpole(rho0_mpole=rho0, iat=iatom, Qlm_car=Qlm_c)
205 :
206 461828 : pab(1:nch_ik, 1) = Qlm_c(1:nch_ik)
207 :
208 24496 : ra(:) = pbc(particle_set(iatom)%r, cell)
209 :
210 : CALL collocate_pgf_product( &
211 : l0_ikind, zet0, 0, 0, 0.0_dp, 0, &
212 : ra, (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, &
213 : rs_grid, ga_gb_function=GRID_FUNC_AB, radius=rpgf0, &
214 56982 : use_subpatch=.TRUE., subpatch_pattern=0)
215 :
216 : END DO ! j
217 :
218 : END DO ! ikind
219 :
220 16018 : IF (ASSOCIATED(cores)) THEN
221 15916 : DEALLOCATE (cores)
222 : END IF
223 :
224 16018 : DEALLOCATE (pab)
225 :
226 16018 : IF (igrid /= auxbas_grid) THEN
227 12 : CALL transfer_rs2pw(rs_grid, coeff_rspace)
228 12 : CALL pw_zero(rho0_s_gs)
229 12 : CALL pw_transfer(coeff_rspace, coeff_gspace)
230 12 : CALL pw_axpy(coeff_gspace, rho0_s_gs)
231 :
232 12 : tot_rs_int = pw_integrate_function(coeff_rspace, isign=-1)
233 :
234 12 : CALL pw_pool%give_back_pw(coeff_rspace)
235 12 : CALL pw_pool%give_back_pw(coeff_gspace)
236 : ELSE
237 :
238 16006 : CALL pw_pool%create_pw(rho0_r_tmp)
239 :
240 16006 : CALL transfer_rs2pw(rs_grid, rho0_r_tmp)
241 :
242 16006 : tot_rs_int = pw_integrate_function(rho0_r_tmp, isign=-1)
243 :
244 16006 : CALL pw_transfer(rho0_r_tmp, rho0_s_rs)
245 16006 : CALL pw_pool%give_back_pw(rho0_r_tmp)
246 :
247 16006 : CALL pw_zero(rho0_s_gs)
248 16006 : CALL pw_transfer(rho0_s_rs, rho0_s_gs)
249 : END IF
250 16018 : CALL timestop(handle)
251 :
252 32036 : END SUBROUTINE put_rho0_on_grid
253 :
254 : ! **************************************************************************************************
255 : !> \brief ...
256 : !> \param pw_env ...
257 : !> \param rho0_mpole ...
258 : ! **************************************************************************************************
259 2620 : SUBROUTINE rho0_s_grid_create(pw_env, rho0_mpole)
260 :
261 : TYPE(pw_env_type), POINTER :: pw_env
262 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
263 :
264 : CHARACTER(len=*), PARAMETER :: routineN = 'rho0_s_grid_create'
265 :
266 : INTEGER :: handle
267 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
268 :
269 2620 : CALL timeset(routineN, handle)
270 :
271 2620 : CPASSERT(ASSOCIATED(pw_env))
272 :
273 2620 : NULLIFY (auxbas_pw_pool)
274 2620 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
275 2620 : CPASSERT(ASSOCIATED(auxbas_pw_pool))
276 :
277 : ! reallocate rho0 on the global grid in real and reciprocal space
278 2620 : CPASSERT(ASSOCIATED(rho0_mpole))
279 :
280 : ! rho0 density in real space
281 2620 : IF (ASSOCIATED(rho0_mpole%rho0_s_rs)) THEN
282 828 : CALL rho0_mpole%rho0_s_rs%release()
283 : ELSE
284 1792 : ALLOCATE (rho0_mpole%rho0_s_rs)
285 : END IF
286 2620 : CALL auxbas_pw_pool%create_pw(rho0_mpole%rho0_s_rs)
287 :
288 : ! rho0 density in reciprocal space
289 2620 : IF (ASSOCIATED(rho0_mpole%rho0_s_gs)) THEN
290 828 : CALL rho0_mpole%rho0_s_gs%release()
291 : ELSE
292 1792 : ALLOCATE (rho0_mpole%rho0_s_gs)
293 : END IF
294 2620 : CALL auxbas_pw_pool%create_pw(rho0_mpole%rho0_s_gs)
295 :
296 : ! Find the grid level suitable for rho0_soft
297 2620 : rho0_mpole%igrid_zet0_s = gaussian_gridlevel(pw_env%gridlevel_info, 2.0_dp*rho0_mpole%zet0_h)
298 :
299 2620 : CALL timestop(handle)
300 :
301 2620 : END SUBROUTINE rho0_s_grid_create
302 :
303 : ! **************************************************************************************************
304 : !> \brief ...
305 : !> \param qs_env ...
306 : !> \param v_rspace ...
307 : !> \param para_env ...
308 : !> \param calculate_forces ...
309 : !> \param local_rho_set ...
310 : !> \param local_rho_set_2nd ...
311 : !> \param atener ...
312 : !> \param kforce ...
313 : ! **************************************************************************************************
314 14954 : SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set, &
315 14954 : local_rho_set_2nd, atener, kforce)
316 :
317 : TYPE(qs_environment_type), POINTER :: qs_env
318 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
319 : TYPE(mp_para_env_type), POINTER :: para_env
320 : LOGICAL, INTENT(IN) :: calculate_forces
321 : TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set, local_rho_set_2nd
322 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atener
323 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: kforce
324 :
325 : CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_vhg0_rspace'
326 :
327 : INTEGER :: auxbas_grid, bo(2), handle, i, iat, iatom, ic, icg, ico, ig1, ig2, igrid, ii, &
328 : ikind, ipgf1, ipgf2, is, iset1, iset2, iso, iso1, iso2, ispin, j, l0_ikind, llmax, lmax0, &
329 : lshell, lx, ly, lz, m1, m2, max_iso_not0_local, max_s_harm, maxl, maxso, mepos, n1, n2, &
330 : nat, nch_ik, nch_max, ncurr, nset, nsotot, nspins, num_pe
331 14954 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
332 14954 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
333 14954 : INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf
334 : LOGICAL :: grid_distributed, paw_atom, use_virial
335 : REAL(KIND=dp) :: eps_rho_rspace, force_tmp(3), fscale, &
336 : ra(3), rpgf0, zet0
337 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
338 14954 : REAL(KIND=dp), DIMENSION(:), POINTER :: hab_sph, norm_l, Qlm
339 14954 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, hdab_sph, intloc, pab
340 14954 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: a_hdab_sph, hdab, Qlm_gg
341 14954 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: a_hdab
342 14954 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
343 : TYPE(cell_type), POINTER :: cell
344 : TYPE(dft_control_type), POINTER :: dft_control
345 : TYPE(gto_basis_set_type), POINTER :: basis_1c_set
346 : TYPE(harmonics_atom_type), POINTER :: harmonics
347 14954 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
348 : TYPE(pw_c1d_gs_type) :: coeff_gaux, coeff_gspace
349 : TYPE(pw_env_type), POINTER :: pw_env
350 14954 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
351 : TYPE(pw_pool_type), POINTER :: pw_aux, pw_pool
352 : TYPE(pw_r3d_rs_type) :: coeff_raux, coeff_rspace
353 14954 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
354 14954 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
355 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
356 14954 : POINTER :: rs_descs
357 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
358 299080 : TYPE(realspace_grid_type) :: rs_v
359 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
360 14954 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
361 14954 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
362 : TYPE(rho_atom_type), POINTER :: rho_atom
363 : TYPE(virial_type), POINTER :: virial
364 :
365 14954 : CALL timeset(routineN, handle)
366 :
367 : ! In case of the external density computed forces probably also
368 : ! need to be stored outside qs_env. We can then remove the
369 : ! attribute 'OPTIONAL' from the argument 'local_rho_set'.
370 :
371 : ! CPASSERT(.NOT. (calculate_forces .AND. PRESENT(local_rho_set)))
372 : ! IF (calculate_forces .AND. PRESENT(local_rho_set)) THEN
373 : ! CPWARN("Forces and External Density!")
374 : ! END IF
375 :
376 14954 : NULLIFY (atomic_kind_set, qs_kind_set, dft_control, particle_set)
377 14954 : NULLIFY (cell, force, pw_env, rho0_mpole, rho_atom_set)
378 :
379 : CALL get_qs_env(qs_env=qs_env, &
380 : atomic_kind_set=atomic_kind_set, &
381 : qs_kind_set=qs_kind_set, &
382 : cell=cell, &
383 : dft_control=dft_control, &
384 : force=force, pw_env=pw_env, &
385 : rho0_mpole=rho0_mpole, &
386 : rho_atom_set=rho_atom_set, &
387 : particle_set=particle_set, &
388 14954 : virial=virial)
389 :
390 14954 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
391 :
392 14954 : nspins = dft_control%nspins
393 :
394 : ! The aim of the following code was to return immediately if the subroutine
395 : ! was called for triplet excited states in spin-restricted case. This check
396 : ! is also performed before invocation of this subroutine. It should be save
397 : ! to remove the optional argument 'do_triplet' from the subroutine interface.
398 : !my_tddft = PRESENT(local_rho_set)
399 : !IF (my_tddft) THEN
400 : ! IF (PRESENT(do_triplet)) THEN
401 : ! IF (nspins == 1 .AND. do_triplet) RETURN
402 : ! ELSE
403 : ! IF (nspins == 1 .AND. dft_control%tddfpt_control%res_etype /= tddfpt_singlet) RETURN
404 : ! END IF
405 : !END IF
406 :
407 14954 : IF (PRESENT(local_rho_set)) &
408 1860 : CALL get_local_rho(local_rho_set, rho0_mpole=rho0_mpole, rho_atom_set=rho_atom_set)
409 : ! Q from rho0_mpole of local_rho_set
410 : ! for TDDFT forces we need mixed potential / integral space
411 : ! potential stored on local_rho_set_2nd
412 14954 : IF (PRESENT(local_rho_set_2nd)) THEN
413 124 : CALL get_local_rho(local_rho_set_2nd, rho_atom_set=rho_atom_set)
414 : END IF
415 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, lmax_0=lmax0, &
416 : zet0_h=zet0, igrid_zet0_s=igrid, &
417 14954 : norm_g0l_h=norm_l)
418 :
419 : ! Setup of the potential on the multigrids
420 14954 : NULLIFY (rs_descs, pw_pools)
421 14954 : CPASSERT(ASSOCIATED(pw_env))
422 14954 : CALL pw_env_get(pw_env, rs_descs=rs_descs, pw_pools=pw_pools)
423 :
424 : ! Assign from pw_env
425 14954 : auxbas_grid = pw_env%auxbas_grid
426 :
427 : ! Get the potential on the right grid
428 14954 : NULLIFY (rs_desc, pw_pool, pw_aux)
429 14954 : rs_desc => rs_descs(igrid)%rs_desc
430 14954 : pw_pool => pw_pools(igrid)%pool
431 :
432 14954 : CALL pw_pool%create_pw(coeff_gspace)
433 14954 : CALL pw_pool%create_pw(coeff_rspace)
434 :
435 14954 : IF (igrid /= auxbas_grid) THEN
436 12 : pw_aux => pw_pools(auxbas_grid)%pool
437 12 : CALL pw_aux%create_pw(coeff_gaux)
438 12 : CALL pw_transfer(v_rspace, coeff_gaux)
439 12 : CALL pw_copy(coeff_gaux, coeff_gspace)
440 12 : CALL pw_transfer(coeff_gspace, coeff_rspace)
441 12 : CALL pw_aux%give_back_pw(coeff_gaux)
442 12 : CALL pw_aux%create_pw(coeff_raux)
443 12 : fscale = coeff_rspace%pw_grid%dvol/coeff_raux%pw_grid%dvol
444 12 : CALL pw_scale(coeff_rspace, fscale)
445 12 : CALL pw_aux%give_back_pw(coeff_raux)
446 : ELSE
447 :
448 14942 : IF (coeff_gspace%pw_grid%spherical) THEN
449 0 : CALL pw_transfer(v_rspace, coeff_gspace)
450 0 : CALL pw_transfer(coeff_gspace, coeff_rspace)
451 : ELSE
452 14942 : CALL pw_copy(v_rspace, coeff_rspace)
453 : END IF
454 : END IF
455 14954 : CALL pw_pool%give_back_pw(coeff_gspace)
456 :
457 : ! Setup the rs grid at level igrid
458 14954 : CALL rs_grid_create(rs_v, rs_desc)
459 14954 : CALL rs_grid_zero(rs_v)
460 14954 : CALL transfer_pw2rs(rs_v, coeff_rspace)
461 :
462 14954 : CALL pw_pool%give_back_pw(coeff_rspace)
463 :
464 : ! Now the potential is on the right grid => integration
465 :
466 14954 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
467 :
468 : ! Allocate work storage
469 :
470 14954 : NULLIFY (hab, hab_sph, hdab, hdab_sph, pab, a_hdab, a_hdab_sph)
471 14954 : nch_max = ncoset(lmax0)
472 14954 : CALL reallocate(hab, 1, nch_max, 1, 1)
473 14954 : CALL reallocate(hab_sph, 1, nch_max)
474 14954 : CALL reallocate(hdab, 1, 3, 1, nch_max, 1, 1)
475 14954 : CALL reallocate(hdab_sph, 1, 3, 1, nch_max)
476 14954 : CALL reallocate(a_hdab, 1, 3, 1, 3, 1, nch_max, 1, 1)
477 14954 : CALL reallocate(a_hdab_sph, 1, 3, 1, 3, 1, nch_max)
478 14954 : CALL reallocate(pab, 1, nch_max, 1, 1)
479 :
480 14954 : ncurr = -1
481 :
482 14954 : grid_distributed = rs_v%desc%distributed
483 :
484 14954 : fscale = 1.0_dp
485 14954 : IF (PRESENT(kforce)) THEN
486 38 : fscale = kforce
487 : END IF
488 :
489 44950 : DO ikind = 1, SIZE(atomic_kind_set, 1)
490 29996 : NULLIFY (basis_1c_set, atom_list, harmonics)
491 29996 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
492 : CALL get_qs_kind(qs_kind_set(ikind), &
493 : basis_set=basis_1c_set, basis_type="GAPW_1C", &
494 : paw_atom=paw_atom, &
495 29996 : harmonics=harmonics)
496 :
497 29996 : IF (.NOT. paw_atom) CYCLE
498 :
499 28498 : NULLIFY (Qlm_gg, lmax, npgf)
500 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, &
501 : l0_ikind=l0_ikind, Qlm_gg=Qlm_gg, & ! Qs different
502 28498 : rpgf0_s=rpgf0)
503 :
504 : CALL get_gto_basis_set(gto_basis_set=basis_1c_set, &
505 : lmax=lmax, lmin=lmin, &
506 : maxso=maxso, maxl=maxl, &
507 28498 : nset=nset, npgf=npgf)
508 :
509 28498 : nsotot = maxso*nset
510 113992 : ALLOCATE (intloc(nsotot, nsotot))
511 :
512 : ! Initialize the local KS integrals
513 :
514 28498 : nch_ik = ncoset(l0_ikind)
515 352606 : pab = 1.0_dp
516 28498 : max_s_harm = harmonics%max_s_harm
517 28498 : llmax = harmonics%llmax
518 :
519 170988 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
520 :
521 28498 : num_pe = para_env%num_pe
522 28498 : mepos = para_env%mepos
523 85494 : DO j = 0, num_pe - 1
524 56996 : bo = get_limit(nat, num_pe, j)
525 56996 : IF (.NOT. grid_distributed .AND. j /= mepos) CYCLE
526 :
527 79751 : DO iat = bo(1), bo(2)
528 22755 : iatom = atom_list(iat)
529 22755 : ra(:) = pbc(particle_set(iatom)%r, cell)
530 :
531 22755 : NULLIFY (Qlm)
532 22755 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, iat=iatom, Qlm_tot=Qlm)
533 :
534 278126 : hab = 0.0_dp
535 975974 : hdab = 0.0_dp
536 65060861 : intloc = 0._dp
537 22755 : IF (use_virial) THEN
538 274 : my_virial_a = 0.0_dp
539 274 : my_virial_b = 0.0_dp
540 36168 : a_hdab = 0.0_dp
541 : END IF
542 :
543 : CALL integrate_pgf_product( &
544 : l0_ikind, zet0, 0, 0, 0.0_dp, 0, &
545 : ra, (/0.0_dp, 0.0_dp, 0.0_dp/), rs_v, &
546 : hab, pab, o1=0, o2=0, &
547 : radius=rpgf0, &
548 : calculate_forces=calculate_forces, &
549 : use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b, &
550 22755 : hdab=hdab, a_hdab=a_hdab, use_subpatch=.TRUE., subpatch_pattern=0)
551 :
552 : ! Convert from cartesian to spherical
553 81876 : DO lshell = 0, l0_ikind
554 251887 : DO is = 1, nso(lshell)
555 170011 : iso = is + nsoset(lshell - 1)
556 170011 : hab_sph(iso) = 0.0_dp
557 680044 : hdab_sph(1:3, iso) = 0.0_dp
558 2210143 : a_hdab_sph(1:3, 1:3, iso) = 0.0_dp
559 998208 : DO ic = 1, nco(lshell)
560 769076 : ico = ic + ncoset(lshell - 1)
561 769076 : lx = indco(1, ico)
562 769076 : ly = indco(2, ico)
563 769076 : lz = indco(3, ico)
564 : hab_sph(iso) = hab_sph(iso) + &
565 : norm_l(lshell)* &
566 : orbtramat(lshell)%slm(is, ic)* &
567 769076 : hab(ico, 1)
568 769076 : IF (calculate_forces) THEN
569 : hdab_sph(1:3, iso) = hdab_sph(1:3, iso) + &
570 : norm_l(lshell)* &
571 : orbtramat(lshell)%slm(is, ic)* &
572 289344 : hdab(1:3, ico, 1)
573 : END IF
574 939087 : IF (use_virial) THEN
575 43840 : DO ii = 1, 3
576 142480 : DO i = 1, 3
577 : a_hdab_sph(i, ii, iso) = a_hdab_sph(i, ii, iso) + &
578 : norm_l(lshell)* &
579 : orbtramat(lshell)%slm(is, ic)* &
580 131520 : a_hdab(i, ii, ico, 1)
581 : END DO
582 : END DO
583 : END IF
584 :
585 : END DO ! ic
586 : END DO ! is
587 : END DO ! lshell
588 :
589 : m1 = 0
590 78597 : DO iset1 = 1, nset
591 :
592 : m2 = 0
593 236134 : DO iset2 = 1, nset
594 : CALL get_none0_cg_list(harmonics%my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
595 180292 : max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
596 180292 : n1 = nsoset(lmax(iset1))
597 624200 : DO ipgf1 = 1, npgf(iset1)
598 443908 : n2 = nsoset(lmax(iset2))
599 1918049 : DO ipgf2 = 1, npgf(iset2)
600 :
601 8183787 : DO iso = 1, MIN(nsoset(l0_ikind), max_iso_not0_local)
602 18201354 : DO icg = 1, cg_n_list(iso)
603 10461475 : iso1 = cg_list(1, icg, iso)
604 10461475 : iso2 = cg_list(2, icg, iso)
605 :
606 10461475 : ig1 = iso1 + n1*(ipgf1 - 1) + m1
607 10461475 : ig2 = iso2 + n2*(ipgf2 - 1) + m2
608 :
609 16907505 : intloc(ig1, ig2) = intloc(ig1, ig2) + Qlm_gg(ig1, ig2, iso)*hab_sph(iso) ! potential times Q
610 :
611 : END DO ! icg
612 : END DO ! iso
613 :
614 : END DO ! ipgf2
615 : END DO ! ipgf1
616 416426 : m2 = m2 + maxso
617 : END DO ! iset2
618 78597 : m1 = m1 + maxso
619 : END DO ! iset1
620 :
621 22755 : IF (grid_distributed) THEN
622 : ! Sum result over all processors
623 0 : CALL para_env%sum(intloc)
624 : END IF
625 :
626 22755 : IF (j == mepos) THEN
627 22755 : rho_atom => rho_atom_set(iatom)
628 22755 : CALL get_rho_atom(rho_atom=rho_atom, ga_Vlocal_gb_h=int_local_h, ga_Vlocal_gb_s=int_local_s)
629 49536 : DO ispin = 1, nspins
630 147860394 : int_local_h(ispin)%r_coef = int_local_h(ispin)%r_coef + intloc
631 147883149 : int_local_s(ispin)%r_coef = int_local_s(ispin)%r_coef + intloc
632 : END DO
633 : END IF
634 :
635 22755 : IF (PRESENT(atener)) THEN
636 178 : DO iso = 1, nsoset(l0_ikind)
637 178 : atener(iatom) = atener(iatom) + 0.5_dp*Qlm(iso)*hab_sph(iso)
638 : END DO
639 : END IF
640 :
641 22755 : IF (calculate_forces) THEN
642 908 : force_tmp(1:3) = 0.0_dp
643 8632 : DO iso = 1, nsoset(l0_ikind)
644 7724 : force_tmp(1) = force_tmp(1) + Qlm(iso)*hdab_sph(1, iso) ! Q here is from local_rho_set
645 7724 : force_tmp(2) = force_tmp(2) + Qlm(iso)*hdab_sph(2, iso)
646 8632 : force_tmp(3) = force_tmp(3) + Qlm(iso)*hdab_sph(3, iso)
647 : END DO
648 3632 : force(ikind)%g0s_Vh_elec(1:3, iat) = force(ikind)%g0s_Vh_elec(1:3, iat) + fscale*force_tmp(1:3)
649 : END IF
650 79751 : IF (use_virial) THEN
651 274 : my_virial_a = 0.0_dp
652 2740 : DO iso = 1, nsoset(l0_ikind)
653 10138 : DO ii = 1, 3
654 32058 : DO i = 1, 3
655 : ! Q from local_rho_set
656 22194 : virial%pv_gapw(i, ii) = virial%pv_gapw(i, ii) + fscale*Qlm(iso)*a_hdab_sph(i, ii, iso)
657 29592 : virial%pv_virial(i, ii) = virial%pv_virial(i, ii) + fscale*Qlm(iso)*a_hdab_sph(i, ii, iso)
658 : END DO
659 : END DO
660 : END DO
661 : END IF
662 :
663 : END DO
664 : END DO
665 :
666 28498 : DEALLOCATE (intloc)
667 103444 : DEALLOCATE (cg_list, cg_n_list)
668 :
669 : END DO ! ikind
670 :
671 14954 : CALL rs_grid_release(rs_v)
672 :
673 14954 : DEALLOCATE (hab, hdab, hab_sph, hdab_sph, pab, a_hdab, a_hdab_sph)
674 :
675 14954 : CALL timestop(handle)
676 :
677 29908 : END SUBROUTINE integrate_vhg0_rspace
678 :
679 : END MODULE qs_rho0_ggrid
|
