(git:ccc2433)
qs_rho0_ggrid.F
Go to the documentation of this file.
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  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  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  REAL(kind=dp), DIMENSION(:), POINTER :: qlm_c
104  REAL(kind=dp), DIMENSION(:, :), POINTER :: pab
105  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  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  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  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
118  TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
119  POINTER :: descs
120  TYPE(realspace_grid_desc_type), POINTER :: desc
121  TYPE(realspace_grid_type), DIMENSION(:), POINTER :: grids
122  TYPE(realspace_grid_type), POINTER :: rs_grid
123 
124  CALL timeset(routinen, handle)
125 
126  NULLIFY (atomic_kind_set, qs_kind_set, cores, pab, qlm_c)
127 
128  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  pw_env=pw_env, cell=cell)
135  eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
136 
137  NULLIFY (descs, pw_pools)
138  CALL pw_env_get(pw_env=pw_env, rs_descs=descs, rs_grids=grids, pw_pools=pw_pools)
139  auxbas_grid = pw_env%auxbas_grid
140 
141  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  rho0_s_rs=rho0_s_rs)
146 
147  ! *** set up the rs grid at level igrid
148  NULLIFY (rs_grid, desc, pw_pool)
149  desc => descs(igrid)%rs_desc
150  rs_grid => grids(igrid)
151  pw_pool => pw_pools(igrid)%pool
152 
153  cpassert(ASSOCIATED(desc))
154  cpassert(ASSOCIATED(pw_pool))
155 
156  IF (igrid /= auxbas_grid) THEN
157  CALL pw_pool%create_pw(coeff_rspace)
158  CALL pw_pool%create_pw(coeff_gspace)
159  END IF
160  CALL rs_grid_zero(rs_grid)
161 
162  nch_max = ncoset(lmax0)
163 
164  ALLOCATE (pab(nch_max, 1))
165 
166  DO ikind = 1, SIZE(atomic_kind_set)
167  CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
168  CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
169 
170  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  rpgf0_s=rpgf0)
174 
175  nch_ik = ncoset(l0_ikind)
176  pab = 0.0_dp
177 
178  CALL reallocate(cores, 1, nat)
179  npme = 0
180  cores = 0
181 
182  DO iat = 1, nat
183  iatom = atom_list(iat)
184  ra(:) = pbc(particle_set(iatom)%r, cell)
185  IF (rs_grid%desc%parallel .AND. .NOT. rs_grid%desc%distributed) THEN
186  ! replicated realspace grid, split the atoms up between procs
187  IF (modulo(nat, rs_grid%desc%group_size) == rs_grid%desc%my_pos) THEN
188  npme = npme + 1
189  cores(npme) = iat
190  END IF
191  ELSE
192  npme = npme + 1
193  cores(npme) = iat
194  END IF
195 
196  END DO
197 
198  ithread = 0
199  DO j = 1, npme
200 
201  iat = cores(j)
202  iatom = atom_list(iat)
203 
204  CALL get_rho0_mpole(rho0_mpole=rho0, iat=iatom, qlm_car=qlm_c)
205 
206  pab(1:nch_ik, 1) = qlm_c(1:nch_ik)
207 
208  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  use_subpatch=.true., subpatch_pattern=0)
215 
216  END DO ! j
217 
218  END DO ! ikind
219 
220  IF (ASSOCIATED(cores)) THEN
221  DEALLOCATE (cores)
222  END IF
223 
224  DEALLOCATE (pab)
225 
226  IF (igrid /= auxbas_grid) THEN
227  CALL transfer_rs2pw(rs_grid, coeff_rspace)
228  CALL pw_zero(rho0_s_gs)
229  CALL pw_transfer(coeff_rspace, coeff_gspace)
230  CALL pw_axpy(coeff_gspace, rho0_s_gs)
231 
232  tot_rs_int = pw_integrate_function(coeff_rspace, isign=-1)
233 
234  CALL pw_pool%give_back_pw(coeff_rspace)
235  CALL pw_pool%give_back_pw(coeff_gspace)
236  ELSE
237 
238  CALL pw_pool%create_pw(rho0_r_tmp)
239 
240  CALL transfer_rs2pw(rs_grid, rho0_r_tmp)
241 
242  tot_rs_int = pw_integrate_function(rho0_r_tmp, isign=-1)
243 
244  CALL pw_transfer(rho0_r_tmp, rho0_s_rs)
245  CALL pw_pool%give_back_pw(rho0_r_tmp)
246 
247  CALL pw_zero(rho0_s_gs)
248  CALL pw_transfer(rho0_s_rs, rho0_s_gs)
249  END IF
250  CALL timestop(handle)
251 
252  END SUBROUTINE put_rho0_on_grid
253 
254 ! **************************************************************************************************
255 !> \brief ...
256 !> \param pw_env ...
257 !> \param rho0_mpole ...
258 ! **************************************************************************************************
259  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  CALL timeset(routinen, handle)
270 
271  cpassert(ASSOCIATED(pw_env))
272 
273  NULLIFY (auxbas_pw_pool)
274  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
275  cpassert(ASSOCIATED(auxbas_pw_pool))
276 
277  ! reallocate rho0 on the global grid in real and reciprocal space
278  cpassert(ASSOCIATED(rho0_mpole))
279 
280  ! rho0 density in real space
281  IF (ASSOCIATED(rho0_mpole%rho0_s_rs)) THEN
282  CALL rho0_mpole%rho0_s_rs%release()
283  ELSE
284  ALLOCATE (rho0_mpole%rho0_s_rs)
285  END IF
286  CALL auxbas_pw_pool%create_pw(rho0_mpole%rho0_s_rs)
287 
288  ! rho0 density in reciprocal space
289  IF (ASSOCIATED(rho0_mpole%rho0_s_gs)) THEN
290  CALL rho0_mpole%rho0_s_gs%release()
291  ELSE
292  ALLOCATE (rho0_mpole%rho0_s_gs)
293  END IF
294  CALL auxbas_pw_pool%create_pw(rho0_mpole%rho0_s_gs)
295 
296  ! Find the grid level suitable for rho0_soft
297  rho0_mpole%igrid_zet0_s = gaussian_gridlevel(pw_env%gridlevel_info, 2.0_dp*rho0_mpole%zet0_h)
298 
299  CALL timestop(handle)
300 
301  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  SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set, &
315  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  INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
332  INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
333  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  REAL(kind=dp), DIMENSION(:), POINTER :: hab_sph, norm_l, qlm
339  REAL(kind=dp), DIMENSION(:, :), POINTER :: hab, hdab_sph, intloc, pab
340  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: a_hdab_sph, hdab, qlm_gg
341  REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: a_hdab
342  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  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  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  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
354  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
355  TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
356  POINTER :: rs_descs
357  TYPE(realspace_grid_desc_type), POINTER :: rs_desc
358  TYPE(realspace_grid_type) :: rs_v
359  TYPE(rho0_mpole_type), POINTER :: rho0_mpole
360  TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
361  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  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  NULLIFY (atomic_kind_set, qs_kind_set, dft_control, particle_set)
377  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  virial=virial)
389 
390  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
391 
392  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  IF (PRESENT(local_rho_set)) &
408  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  IF (PRESENT(local_rho_set_2nd)) THEN
413  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  norm_g0l_h=norm_l)
418 
419  ! Setup of the potential on the multigrids
420  NULLIFY (rs_descs, pw_pools)
421  cpassert(ASSOCIATED(pw_env))
422  CALL pw_env_get(pw_env, rs_descs=rs_descs, pw_pools=pw_pools)
423 
424  ! Assign from pw_env
425  auxbas_grid = pw_env%auxbas_grid
426 
427  ! Get the potential on the right grid
428  NULLIFY (rs_desc, pw_pool, pw_aux)
429  rs_desc => rs_descs(igrid)%rs_desc
430  pw_pool => pw_pools(igrid)%pool
431 
432  CALL pw_pool%create_pw(coeff_gspace)
433  CALL pw_pool%create_pw(coeff_rspace)
434 
435  IF (igrid /= auxbas_grid) THEN
436  pw_aux => pw_pools(auxbas_grid)%pool
437  CALL pw_aux%create_pw(coeff_gaux)
438  CALL pw_transfer(v_rspace, coeff_gaux)
439  CALL pw_copy(coeff_gaux, coeff_gspace)
440  CALL pw_transfer(coeff_gspace, coeff_rspace)
441  CALL pw_aux%give_back_pw(coeff_gaux)
442  CALL pw_aux%create_pw(coeff_raux)
443  fscale = coeff_rspace%pw_grid%dvol/coeff_raux%pw_grid%dvol
444  CALL pw_scale(coeff_rspace, fscale)
445  CALL pw_aux%give_back_pw(coeff_raux)
446  ELSE
447 
448  IF (coeff_gspace%pw_grid%spherical) THEN
449  CALL pw_transfer(v_rspace, coeff_gspace)
450  CALL pw_transfer(coeff_gspace, coeff_rspace)
451  ELSE
452  CALL pw_copy(v_rspace, coeff_rspace)
453  END IF
454  END IF
455  CALL pw_pool%give_back_pw(coeff_gspace)
456 
457  ! Setup the rs grid at level igrid
458  CALL rs_grid_create(rs_v, rs_desc)
459  CALL rs_grid_zero(rs_v)
460  CALL transfer_pw2rs(rs_v, coeff_rspace)
461 
462  CALL pw_pool%give_back_pw(coeff_rspace)
463 
464  ! Now the potential is on the right grid => integration
465 
466  eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
467 
468  ! Allocate work storage
469 
470  NULLIFY (hab, hab_sph, hdab, hdab_sph, pab, a_hdab, a_hdab_sph)
471  nch_max = ncoset(lmax0)
472  CALL reallocate(hab, 1, nch_max, 1, 1)
473  CALL reallocate(hab_sph, 1, nch_max)
474  CALL reallocate(hdab, 1, 3, 1, nch_max, 1, 1)
475  CALL reallocate(hdab_sph, 1, 3, 1, nch_max)
476  CALL reallocate(a_hdab, 1, 3, 1, 3, 1, nch_max, 1, 1)
477  CALL reallocate(a_hdab_sph, 1, 3, 1, 3, 1, nch_max)
478  CALL reallocate(pab, 1, nch_max, 1, 1)
479 
480  ncurr = -1
481 
482  grid_distributed = rs_v%desc%distributed
483 
484  fscale = 1.0_dp
485  IF (PRESENT(kforce)) THEN
486  fscale = kforce
487  END IF
488 
489  DO ikind = 1, SIZE(atomic_kind_set, 1)
490  NULLIFY (basis_1c_set, atom_list, harmonics)
491  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  harmonics=harmonics)
496 
497  IF (.NOT. paw_atom) cycle
498 
499  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  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  nset=nset, npgf=npgf)
508 
509  nsotot = maxso*nset
510  ALLOCATE (intloc(nsotot, nsotot))
511 
512  ! Initialize the local KS integrals
513 
514  nch_ik = ncoset(l0_ikind)
515  pab = 1.0_dp
516  max_s_harm = harmonics%max_s_harm
517  llmax = harmonics%llmax
518 
519  ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
520 
521  num_pe = para_env%num_pe
522  mepos = para_env%mepos
523  DO j = 0, num_pe - 1
524  bo = get_limit(nat, num_pe, j)
525  IF (.NOT. grid_distributed .AND. j /= mepos) cycle
526 
527  DO iat = bo(1), bo(2)
528  iatom = atom_list(iat)
529  ra(:) = pbc(particle_set(iatom)%r, cell)
530 
531  NULLIFY (qlm)
532  CALL get_rho0_mpole(rho0_mpole=rho0_mpole, iat=iatom, qlm_tot=qlm)
533 
534  hab = 0.0_dp
535  hdab = 0.0_dp
536  intloc = 0._dp
537  IF (use_virial) THEN
538  my_virial_a = 0.0_dp
539  my_virial_b = 0.0_dp
540  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  hdab=hdab, a_hdab=a_hdab, use_subpatch=.true., subpatch_pattern=0)
551 
552  ! Convert from cartesian to spherical
553  DO lshell = 0, l0_ikind
554  DO is = 1, nso(lshell)
555  iso = is + nsoset(lshell - 1)
556  hab_sph(iso) = 0.0_dp
557  hdab_sph(1:3, iso) = 0.0_dp
558  a_hdab_sph(1:3, 1:3, iso) = 0.0_dp
559  DO ic = 1, nco(lshell)
560  ico = ic + ncoset(lshell - 1)
561  lx = indco(1, ico)
562  ly = indco(2, ico)
563  lz = indco(3, ico)
564  hab_sph(iso) = hab_sph(iso) + &
565  norm_l(lshell)* &
566  orbtramat(lshell)%slm(is, ic)* &
567  hab(ico, 1)
568  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  hdab(1:3, ico, 1)
573  END IF
574  IF (use_virial) THEN
575  DO ii = 1, 3
576  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  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  DO iset1 = 1, nset
591 
592  m2 = 0
593  DO iset2 = 1, nset
594  CALL get_none0_cg_list(harmonics%my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
595  max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
596  n1 = nsoset(lmax(iset1))
597  DO ipgf1 = 1, npgf(iset1)
598  n2 = nsoset(lmax(iset2))
599  DO ipgf2 = 1, npgf(iset2)
600 
601  DO iso = 1, min(nsoset(l0_ikind), max_iso_not0_local)
602  DO icg = 1, cg_n_list(iso)
603  iso1 = cg_list(1, icg, iso)
604  iso2 = cg_list(2, icg, iso)
605 
606  ig1 = iso1 + n1*(ipgf1 - 1) + m1
607  ig2 = iso2 + n2*(ipgf2 - 1) + m2
608 
609  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  m2 = m2 + maxso
617  END DO ! iset2
618  m1 = m1 + maxso
619  END DO ! iset1
620 
621  IF (grid_distributed) THEN
622  ! Sum result over all processors
623  CALL para_env%sum(intloc)
624  END IF
625 
626  IF (j == mepos) THEN
627  rho_atom => rho_atom_set(iatom)
628  CALL get_rho_atom(rho_atom=rho_atom, ga_vlocal_gb_h=int_local_h, ga_vlocal_gb_s=int_local_s)
629  DO ispin = 1, nspins
630  int_local_h(ispin)%r_coef = int_local_h(ispin)%r_coef + intloc
631  int_local_s(ispin)%r_coef = int_local_s(ispin)%r_coef + intloc
632  END DO
633  END IF
634 
635  IF (PRESENT(atener)) THEN
636  DO iso = 1, nsoset(l0_ikind)
637  atener(iatom) = atener(iatom) + 0.5_dp*qlm(iso)*hab_sph(iso)
638  END DO
639  END IF
640 
641  IF (calculate_forces) THEN
642  force_tmp(1:3) = 0.0_dp
643  DO iso = 1, nsoset(l0_ikind)
644  force_tmp(1) = force_tmp(1) + qlm(iso)*hdab_sph(1, iso) ! Q here is from local_rho_set
645  force_tmp(2) = force_tmp(2) + qlm(iso)*hdab_sph(2, iso)
646  force_tmp(3) = force_tmp(3) + qlm(iso)*hdab_sph(3, iso)
647  END DO
648  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  IF (use_virial) THEN
651  my_virial_a = 0.0_dp
652  DO iso = 1, nsoset(l0_ikind)
653  DO ii = 1, 3
654  DO i = 1, 3
655  ! Q from local_rho_set
656  virial%pv_gapw(i, ii) = virial%pv_gapw(i, ii) + fscale*qlm(iso)*a_hdab_sph(i, ii, iso)
657  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  DEALLOCATE (intloc)
667  DEALLOCATE (cg_list, cg_n_list)
668 
669  END DO ! ikind
670 
671  CALL rs_grid_release(rs_v)
672 
673  DEALLOCATE (hab, hdab, hab_sph, hdab_sph, pab, a_hdab, a_hdab_sph)
674 
675  CALL timestop(handle)
676 
677  END SUBROUTINE integrate_vhg0_rspace
678 
679 END MODULE qs_rho0_ggrid
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
basis_set_types
Definition: basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius)
...
Definition: basis_set_types.F:615
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
gaussian_gridlevels
Definition: gaussian_gridlevels.F:13
gaussian_gridlevels::gaussian_gridlevel
integer function, public gaussian_gridlevel(gridlevel_info, exponent)
...
Definition: gaussian_gridlevels.F:148
grid_api
Fortran API for the grid package, which is written in C.
Definition: grid_api.F:12
grid_api::grid_func_ab
integer, parameter, public grid_func_ab
Definition: grid_api.F:29
grid_api::collocate_pgf_product
subroutine, public collocate_pgf_product(la_max, zeta, la_min, lb_max, zetb, lb_min, ra, rab, scale, pab, o1, o2, rsgrid, ga_gb_function, radius, use_subpatch, subpatch_pattern)
low level collocation of primitive gaussian functions
Definition: grid_api.F:119
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
memory_utilities
Utility routines for the memory handling.
Definition: memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition: orbital_pointers.F:42
orbital_pointers::nsoset
integer, dimension(:), allocatable, public nsoset
Definition: orbital_pointers.F:42
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition: orbital_pointers.F:42
orbital_pointers::indco
integer, dimension(:, :), allocatable, public indco
Definition: orbital_pointers.F:43
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition: orbital_pointers.F:42
orbital_transformation_matrices
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
Definition: orbital_transformation_matrices.F:17
orbital_transformation_matrices::orbtramat
type(orbtramat_type), dimension(:), pointer, public orbtramat
Definition: orbital_transformation_matrices.F:47
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition: pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
pw_methods
Definition: pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
pw_types
Definition: pw_types.F:24
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_force_types
Definition: qs_force_types.F:13
qs_harmonics_atom
Definition: qs_harmonics_atom.F:8
qs_integrate_potential
Integrate single or product functions over a potential on a RS grid.
Definition: qs_integrate_potential.F:22
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_local_rho_types
Definition: qs_local_rho_types.F:8
qs_local_rho_types::get_local_rho
subroutine, public get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set)
...
Definition: qs_local_rho_types.F:173
qs_rho0_ggrid
Definition: qs_rho0_ggrid.F:8
qs_rho0_ggrid::integrate_vhg0_rspace
subroutine, public integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set, local_rho_set_2nd, atener, kforce)
...
Definition: qs_rho0_ggrid.F:316
qs_rho0_ggrid::put_rho0_on_grid
subroutine, public put_rho0_on_grid(qs_env, rho0, tot_rs_int)
...
Definition: qs_rho0_ggrid.F:89
qs_rho0_ggrid::rho0_s_grid_create
subroutine, public rho0_s_grid_create(pw_env, rho0_mpole)
...
Definition: qs_rho0_ggrid.F:260
qs_rho0_types
Definition: qs_rho0_types.F:8
qs_rho0_types::get_rho0_mpole
subroutine, public get_rho0_mpole(rho0_mpole, g0_h, vg0_h, iat, ikind, lmax_0, l0_ikind, mp_gau_ikind, mp_rho, norm_g0l_h, Qlm_gg, Qlm_car, Qlm_tot, zet0_h, igrid_zet0_s, rpgf0_h, rpgf0_s, max_rpgf0_s, rho0_s_rs, rho0_s_gs)
...
Definition: qs_rho0_types.F:425
qs_rho_atom_types
Definition: qs_rho_atom_types.F:8
qs_rho_atom_types::get_rho_atom
subroutine, public get_rho_atom(rho_atom, cpc_h, cpc_s, rho_rad_h, rho_rad_s, drho_rad_h, drho_rad_s, vrho_rad_h, vrho_rad_s, rho_rad_h_d, rho_rad_s_d, ga_Vlocal_gb_h, ga_Vlocal_gb_s, int_scr_h, int_scr_s)
...
Definition: qs_rho_atom_types.F:195
realspace_grid_types
Definition: realspace_grid_types.F:18
realspace_grid_types::rs_grid_create
subroutine, public rs_grid_create(rs, desc)
...
Definition: realspace_grid_types.F:482
realspace_grid_types::transfer_pw2rs
subroutine, public transfer_pw2rs(rs, pw)
...
Definition: realspace_grid_types.F:698
realspace_grid_types::transfer_rs2pw
subroutine, public transfer_rs2pw(rs, pw)
...
Definition: realspace_grid_types.F:656
realspace_grid_types::rs_grid_release
subroutine, public rs_grid_release(rs_grid)
releases the given rs grid (see doc/ReferenceCounting.html)
Definition: realspace_grid_types.F:2016
realspace_grid_types::rs_grid_zero
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
Definition: realspace_grid_types.F:1909
util
All kind of helpful little routines.
Definition: util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition: util.F:333
virial_types
Definition: virial_types.F:13