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 : ! **************************************************************************************************
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_cneo_ggrid, ONLY: integrate_vhgg_rspace,&
42 : rhoz_cneo_s_grid_create
43 : USE qs_cneo_types, ONLY: cneo_potential_type,&
44 : rhoz_cneo_type
45 : USE qs_environment_types, ONLY: get_qs_env,&
46 : qs_environment_type
47 : USE qs_force_types, ONLY: qs_force_type
48 : USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
49 : harmonics_atom_type
50 : USE qs_integrate_potential, ONLY: integrate_pgf_product
51 : USE qs_kind_types, ONLY: get_qs_kind,&
52 : qs_kind_type
53 : USE qs_local_rho_types, ONLY: get_local_rho,&
54 : local_rho_type
55 : USE qs_rho0_types, ONLY: get_rho0_mpole,&
56 : rho0_mpole_type
57 : USE qs_rho_atom_types, ONLY: get_rho_atom,&
58 : rho_atom_coeff,&
59 : rho_atom_type
60 : USE realspace_grid_types, ONLY: realspace_grid_desc_p_type,&
61 : realspace_grid_desc_type,&
62 : realspace_grid_type,&
63 : rs_grid_create,&
64 : rs_grid_release,&
65 : rs_grid_zero,&
66 : transfer_pw2rs,&
67 : transfer_rs2pw
68 : USE util, ONLY: get_limit
69 : USE virial_types, ONLY: virial_type
70 : #include "./base/base_uses.f90"
71 :
72 : IMPLICIT NONE
73 :
74 : PRIVATE
75 :
76 : ! Global parameters (only in this module)
77 :
78 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho0_ggrid'
79 :
80 : ! Public subroutines
81 :
82 : PUBLIC :: put_rho0_on_grid, rho0_s_grid_create, integrate_vhg0_rspace
83 :
84 : CONTAINS
85 :
86 : ! **************************************************************************************************
87 : !> \brief ...
88 : !> \param qs_env ...
89 : !> \param rho0 ...
90 : !> \param tot_rs_int ...
91 : !> \param my_pools ...
92 : !> \param my_rs_grids ...
93 : !> \param my_rs_descs ...
94 : ! **************************************************************************************************
95 31924 : SUBROUTINE put_rho0_on_grid(qs_env, rho0, tot_rs_int, my_pools, my_rs_grids, my_rs_descs)
96 :
97 : TYPE(qs_environment_type), POINTER :: qs_env
98 : TYPE(rho0_mpole_type), POINTER :: rho0
99 : REAL(KIND=dp), INTENT(OUT) :: tot_rs_int
100 : TYPE(pw_pool_p_type), DIMENSION(:), OPTIONAL, &
101 : POINTER :: my_pools
102 : TYPE(realspace_grid_type), DIMENSION(:), &
103 : OPTIONAL, POINTER :: my_rs_grids
104 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
105 : OPTIONAL, POINTER :: my_rs_descs
106 :
107 : CHARACTER(LEN=*), PARAMETER :: routineN = 'put_rho0_on_grid'
108 :
109 : INTEGER :: auxbas_grid, handle, iat, iatom, igrid, &
110 : ikind, ithread, j, l0_ikind, lmax0, &
111 : nat, nch_ik, nch_max, npme
112 31924 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
113 : LOGICAL :: paw_atom
114 : REAL(KIND=dp) :: eps_rho_rspace, rpgf0, zet0
115 : REAL(KIND=dp), DIMENSION(3) :: ra
116 31924 : REAL(KIND=dp), DIMENSION(:), POINTER :: Qlm_c
117 31924 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
118 31924 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
119 : TYPE(cell_type), POINTER :: cell
120 : TYPE(dft_control_type), POINTER :: dft_control
121 : TYPE(mp_para_env_type), POINTER :: para_env
122 31924 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
123 : TYPE(pw_c1d_gs_type) :: coeff_gspace
124 : TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs
125 : TYPE(pw_env_type), POINTER :: pw_env
126 31924 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
127 : TYPE(pw_pool_type), POINTER :: pw_pool
128 : TYPE(pw_r3d_rs_type) :: coeff_rspace
129 : TYPE(pw_r3d_rs_type), POINTER :: rho0_s_rs
130 31924 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
131 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
132 31924 : POINTER :: descs
133 : TYPE(realspace_grid_desc_type), POINTER :: desc
134 31924 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: grids
135 : TYPE(realspace_grid_type), POINTER :: rs_grid
136 :
137 31924 : CALL timeset(routineN, handle)
138 :
139 31924 : NULLIFY (atomic_kind_set, qs_kind_set, cores, pab, Qlm_c)
140 :
141 31924 : NULLIFY (dft_control, pw_env, particle_set, para_env, cell, rho0_s_gs, rho0_s_rs)
142 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
143 : particle_set=particle_set, &
144 : atomic_kind_set=atomic_kind_set, &
145 : qs_kind_set=qs_kind_set, &
146 : para_env=para_env, &
147 31924 : pw_env=pw_env, cell=cell)
148 31924 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
149 :
150 31924 : NULLIFY (descs, pw_pools)
151 31924 : CALL pw_env_get(pw_env=pw_env, rs_descs=descs, rs_grids=grids, pw_pools=pw_pools)
152 31924 : auxbas_grid = pw_env%auxbas_grid
153 :
154 31924 : NULLIFY (rho0_s_gs, rho0_s_rs)
155 : CALL get_rho0_mpole(rho0_mpole=rho0, lmax_0=lmax0, &
156 : zet0_h=zet0, igrid_zet0_s=igrid, &
157 : rho0_s_gs=rho0_s_gs, &
158 31924 : rho0_s_rs=rho0_s_rs)
159 :
160 : ! *** set up the rs grid at level igrid
161 31924 : NULLIFY (rs_grid, desc, pw_pool)
162 : ! IF present, overwrite qs grid for new pool
163 31924 : IF (PRESENT(my_pools)) THEN
164 2288 : desc => my_rs_descs(igrid)%rs_desc
165 2288 : rs_grid => my_rs_grids(igrid)
166 2288 : pw_pool => my_pools(igrid)%pool
167 : ELSE
168 29636 : desc => descs(igrid)%rs_desc
169 29636 : rs_grid => grids(igrid)
170 29636 : pw_pool => pw_pools(igrid)%pool
171 : END IF
172 :
173 31924 : CPASSERT(ASSOCIATED(desc))
174 31924 : CPASSERT(ASSOCIATED(pw_pool))
175 :
176 31924 : IF (igrid /= auxbas_grid) THEN
177 12 : CALL pw_pool%create_pw(coeff_rspace)
178 12 : CALL pw_pool%create_pw(coeff_gspace)
179 : END IF
180 31924 : CALL rs_grid_zero(rs_grid)
181 :
182 31924 : nch_max = ncoset(lmax0)
183 :
184 95772 : ALLOCATE (pab(nch_max, 1))
185 :
186 93080 : DO ikind = 1, SIZE(atomic_kind_set)
187 61156 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
188 61156 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
189 :
190 61156 : IF (.NOT. paw_atom .AND. dft_control%qs_control%gapw_control%nopaw_as_gpw) CYCLE
191 :
192 : CALL get_rho0_mpole(rho0_mpole=rho0, ikind=ikind, l0_ikind=l0_ikind, &
193 57704 : rpgf0_s=rpgf0)
194 :
195 57704 : nch_ik = ncoset(l0_ikind)
196 730954 : pab = 0.0_dp
197 :
198 57704 : CALL reallocate(cores, 1, nat)
199 57704 : npme = 0
200 148280 : cores = 0
201 :
202 148280 : DO iat = 1, nat
203 90576 : iatom = atom_list(iat)
204 90576 : ra(:) = pbc(particle_set(iatom)%r, cell)
205 148280 : IF (rs_grid%desc%parallel .AND. .NOT. rs_grid%desc%distributed) THEN
206 : ! replicated realspace grid, split the atoms up between procs
207 90576 : IF (MODULO(iatom, rs_grid%desc%group_size) == rs_grid%desc%my_pos) THEN
208 45288 : npme = npme + 1
209 45288 : cores(npme) = iat
210 : END IF
211 : ELSE
212 0 : npme = npme + 1
213 0 : cores(npme) = iat
214 : END IF
215 :
216 : END DO
217 :
218 : ithread = 0
219 253776 : DO j = 1, npme
220 :
221 45288 : iat = cores(j)
222 45288 : iatom = atom_list(iat)
223 :
224 45288 : CALL get_rho0_mpole(rho0_mpole=rho0, iat=iatom, Qlm_car=Qlm_c)
225 :
226 882074 : pab(1:nch_ik, 1) = Qlm_c(1:nch_ik)
227 :
228 45288 : ra(:) = pbc(particle_set(iatom)%r, cell)
229 :
230 : CALL collocate_pgf_product( &
231 : l0_ikind, zet0, 0, 0, 0.0_dp, 0, &
232 : ra, [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, pab, 0, 0, &
233 : rs_grid, ga_gb_function=GRID_FUNC_AB, radius=rpgf0, &
234 106444 : use_subpatch=.TRUE., subpatch_pattern=0)
235 :
236 : END DO ! j
237 :
238 : END DO ! ikind
239 :
240 31924 : IF (ASSOCIATED(cores)) THEN
241 30680 : DEALLOCATE (cores)
242 : END IF
243 :
244 31924 : DEALLOCATE (pab)
245 :
246 31924 : IF (igrid /= auxbas_grid) THEN
247 12 : CALL transfer_rs2pw(rs_grid, coeff_rspace)
248 12 : CALL pw_zero(rho0_s_gs)
249 12 : CALL pw_transfer(coeff_rspace, coeff_gspace)
250 12 : CALL pw_axpy(coeff_gspace, rho0_s_gs)
251 :
252 12 : tot_rs_int = pw_integrate_function(coeff_rspace, isign=-1)
253 :
254 12 : CALL pw_pool%give_back_pw(coeff_rspace)
255 12 : CALL pw_pool%give_back_pw(coeff_gspace)
256 : ELSE
257 31912 : IF (ASSOCIATED(rs_grid%desc%pw, rho0_s_rs%pw_grid)) THEN
258 31912 : CALL transfer_rs2pw(rs_grid, rho0_s_rs)
259 31912 : tot_rs_int = pw_integrate_function(rho0_s_rs, isign=-1)
260 : ELSE
261 : ! Custom pools may use a distinct grid descriptor at the same level.
262 0 : CALL pw_pool%create_pw(coeff_rspace)
263 0 : CALL transfer_rs2pw(rs_grid, coeff_rspace)
264 0 : tot_rs_int = pw_integrate_function(coeff_rspace, isign=-1)
265 0 : CALL pw_transfer(coeff_rspace, rho0_s_rs)
266 0 : CALL pw_pool%give_back_pw(coeff_rspace)
267 : END IF
268 :
269 31912 : CALL pw_zero(rho0_s_gs)
270 31912 : CALL pw_transfer(rho0_s_rs, rho0_s_gs)
271 : END IF
272 31924 : CALL timestop(handle)
273 :
274 63848 : END SUBROUTINE put_rho0_on_grid
275 :
276 : ! **************************************************************************************************
277 : !> \brief ...
278 : !> \param pw_env ...
279 : !> \param rho0_mpole ...
280 : ! **************************************************************************************************
281 5294 : SUBROUTINE rho0_s_grid_create(pw_env, rho0_mpole)
282 :
283 : TYPE(pw_env_type), POINTER :: pw_env
284 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
285 :
286 : CHARACTER(len=*), PARAMETER :: routineN = 'rho0_s_grid_create'
287 :
288 : INTEGER :: handle
289 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
290 :
291 5294 : CALL timeset(routineN, handle)
292 :
293 5294 : CPASSERT(ASSOCIATED(pw_env))
294 :
295 5294 : NULLIFY (auxbas_pw_pool)
296 5294 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
297 5294 : CPASSERT(ASSOCIATED(auxbas_pw_pool))
298 :
299 : ! reallocate rho0 on the global grid in real and reciprocal space
300 5294 : CPASSERT(ASSOCIATED(rho0_mpole))
301 :
302 : ! rho0 density in real space
303 5294 : IF (ASSOCIATED(rho0_mpole%rho0_s_rs)) THEN
304 1480 : CALL rho0_mpole%rho0_s_rs%release()
305 : ELSE
306 3814 : ALLOCATE (rho0_mpole%rho0_s_rs)
307 : END IF
308 5294 : CALL auxbas_pw_pool%create_pw(rho0_mpole%rho0_s_rs)
309 :
310 : ! rho0 density in reciprocal space
311 5294 : IF (ASSOCIATED(rho0_mpole%rho0_s_gs)) THEN
312 1480 : CALL rho0_mpole%rho0_s_gs%release()
313 : ELSE
314 3814 : ALLOCATE (rho0_mpole%rho0_s_gs)
315 : END IF
316 5294 : CALL auxbas_pw_pool%create_pw(rho0_mpole%rho0_s_gs)
317 :
318 : ! Find the grid level suitable for rho0_soft
319 5294 : rho0_mpole%igrid_zet0_s = gaussian_gridlevel(pw_env%gridlevel_info, 2.0_dp*rho0_mpole%zet0_h)
320 :
321 5294 : CALL timestop(handle)
322 :
323 5294 : IF (rho0_mpole%do_cneo) THEN
324 16 : CALL rhoz_cneo_s_grid_create(pw_env, rho0_mpole)
325 : END IF
326 :
327 5294 : END SUBROUTINE rho0_s_grid_create
328 :
329 : ! **************************************************************************************************
330 : !> \brief ...
331 : !> \param qs_env ...
332 : !> \param v_rspace ...
333 : !> \param para_env ...
334 : !> \param calculate_forces ...
335 : !> \param local_rho_set ...
336 : !> \param local_rho_set_2nd ...
337 : !> \param atener ...
338 : !> \param kforce ...
339 : !> \param my_pools ...
340 : !> \param my_rs_descs ...
341 : ! **************************************************************************************************
342 29254 : SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set, &
343 29254 : local_rho_set_2nd, atener, kforce, my_pools, my_rs_descs)
344 :
345 : TYPE(qs_environment_type), POINTER :: qs_env
346 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
347 : TYPE(mp_para_env_type), POINTER :: para_env
348 : LOGICAL, INTENT(IN) :: calculate_forces
349 : TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set, local_rho_set_2nd
350 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atener
351 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: kforce
352 : TYPE(pw_pool_p_type), DIMENSION(:), OPTIONAL, &
353 : POINTER :: my_pools
354 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
355 : OPTIONAL, POINTER :: my_rs_descs
356 :
357 : CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_vhg0_rspace'
358 :
359 : INTEGER :: auxbas_grid, bo(2), handle, i, iat, iatom, ic, icg, ico, ig1, ig2, igrid, ii, &
360 : ikind, ipgf1, ipgf2, is, iset1, iset2, iso, iso1, iso2, ispin, j, l0_ikind, llmax, &
361 : llmax_nuc, lmax0, lshell, lx, ly, lz, m1, m2, max_iso_not0_local, max_s_harm, &
362 : max_s_harm_nuc, maxl, maxl_nuc, maxso, maxso_nuc, mepos, n1, n2, nat, nch_ik, nch_max, &
363 : ncurr, nset, nset_nuc, nsotot, nsotot_nuc, nspins, num_pe
364 29254 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
365 29254 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
366 29254 : INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmax_nuc, lmin, &
367 29254 : lmin_nuc, npgf, npgf_nuc
368 : LOGICAL :: direct_transfer, grid_distributed, &
369 : paw_atom, use_virial
370 : REAL(KIND=dp) :: eps_rho_rspace, force_tmp(3), fscale, &
371 : ra(3), rpgf0, zet0
372 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
373 29254 : REAL(KIND=dp), DIMENSION(:), POINTER :: hab_sph, norm_l, Qlm
374 29254 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, hdab_sph, intloc, intloc_nuc, pab
375 29254 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: a_hdab_sph, hdab, Qlm_gg
376 29254 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: a_hdab
377 29254 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
378 : TYPE(cell_type), POINTER :: cell
379 : TYPE(cneo_potential_type), POINTER :: cneo_potential
380 : TYPE(dft_control_type), POINTER :: dft_control
381 : TYPE(gto_basis_set_type), POINTER :: basis_1c_set, nuc_basis_set
382 : TYPE(harmonics_atom_type), POINTER :: harmonics
383 29254 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
384 : TYPE(pw_c1d_gs_type) :: coeff_gaux, coeff_gspace
385 : TYPE(pw_env_type), POINTER :: pw_env
386 29254 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
387 : TYPE(pw_pool_type), POINTER :: pw_aux, pw_pool
388 : TYPE(pw_r3d_rs_type) :: coeff_rspace
389 29254 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
390 29254 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
391 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
392 29254 : POINTER :: rs_descs
393 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
394 585080 : TYPE(realspace_grid_type) :: rs_v
395 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
396 29254 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
397 29254 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
398 : TYPE(rho_atom_type), POINTER :: rho_atom
399 29254 : TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
400 : TYPE(rhoz_cneo_type), POINTER :: rhoz_cneo
401 : TYPE(virial_type), POINTER :: virial
402 :
403 29254 : CALL timeset(routineN, handle)
404 :
405 : ! In case of the external density computed forces probably also
406 : ! need to be stored outside qs_env. We can then remove the
407 : ! attribute 'OPTIONAL' from the argument 'local_rho_set'.
408 :
409 : ! CPASSERT(.NOT. (calculate_forces .AND. PRESENT(local_rho_set)))
410 : ! IF (calculate_forces .AND. PRESENT(local_rho_set)) THEN
411 : ! CPWARN("Forces and External Density!")
412 : ! END IF
413 :
414 29254 : NULLIFY (atomic_kind_set, qs_kind_set, dft_control, particle_set)
415 29254 : NULLIFY (cell, force, pw_env, rho0_mpole, rho_atom_set, rhoz_cneo_set)
416 :
417 : CALL get_qs_env(qs_env=qs_env, &
418 : atomic_kind_set=atomic_kind_set, &
419 : qs_kind_set=qs_kind_set, &
420 : cell=cell, &
421 : dft_control=dft_control, &
422 : force=force, pw_env=pw_env, &
423 : rho0_mpole=rho0_mpole, &
424 : rho_atom_set=rho_atom_set, &
425 : rhoz_cneo_set=rhoz_cneo_set, &
426 : particle_set=particle_set, &
427 29254 : virial=virial)
428 :
429 29254 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
430 :
431 29254 : nspins = dft_control%nspins
432 :
433 : ! The aim of the following code was to return immediately if the subroutine
434 : ! was called for triplet excited states in spin-restricted case. This check
435 : ! is also performed before invocation of this subroutine. It should be save
436 : ! to remove the optional argument 'do_triplet' from the subroutine interface.
437 : !my_tddft = PRESENT(local_rho_set)
438 : !IF (my_tddft) THEN
439 : ! IF (PRESENT(do_triplet)) THEN
440 : ! IF (nspins == 1 .AND. do_triplet) RETURN
441 : ! ELSE
442 : ! IF (nspins == 1 .AND. dft_control%tddfpt_control%res_etype /= tddfpt_singlet) RETURN
443 : ! END IF
444 : !END IF
445 :
446 29254 : IF (PRESENT(local_rho_set)) THEN
447 : CALL get_local_rho(local_rho_set, rho0_mpole=rho0_mpole, rho_atom_set=rho_atom_set, &
448 4792 : rhoz_cneo_set=rhoz_cneo_set)
449 : END IF
450 : ! Q from rho0_mpole of local_rho_set
451 : ! for TDDFT forces we need mixed potential / integral space
452 : ! potential stored on local_rho_set_2nd
453 29254 : IF (PRESENT(local_rho_set_2nd)) THEN
454 348 : CALL get_local_rho(local_rho_set_2nd, rho_atom_set=rho_atom_set)
455 : END IF
456 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, lmax_0=lmax0, &
457 : zet0_h=zet0, igrid_zet0_s=igrid, &
458 29254 : norm_g0l_h=norm_l)
459 :
460 : ! Setup of the potential on the multigrids
461 29254 : NULLIFY (rs_descs, pw_pools)
462 29254 : CPASSERT(ASSOCIATED(pw_env))
463 29254 : CALL pw_env_get(pw_env, rs_descs=rs_descs, pw_pools=pw_pools)
464 :
465 : ! Assign from pw_env
466 29254 : auxbas_grid = pw_env%auxbas_grid
467 :
468 : ! IF present, overwrite qs grid for new pool
469 : ! Get the potential on the right grid
470 29254 : IF (PRESENT(my_pools)) THEN
471 50 : rs_desc => my_rs_descs(igrid)%rs_desc
472 50 : pw_pool => my_pools(igrid)%pool
473 : ELSE
474 29204 : rs_desc => rs_descs(igrid)%rs_desc
475 29204 : pw_pool => pw_pools(igrid)%pool
476 : END IF
477 :
478 : direct_transfer = igrid == auxbas_grid .AND. .NOT. pw_pool%pw_grid%spherical .AND. &
479 29254 : ASSOCIATED(rs_desc%pw, v_rspace%pw_grid)
480 : IF (.NOT. direct_transfer) THEN
481 12 : CALL pw_pool%create_pw(coeff_rspace)
482 12 : IF (igrid /= auxbas_grid .OR. pw_pool%pw_grid%spherical) THEN
483 12 : CALL pw_pool%create_pw(coeff_gspace)
484 12 : IF (igrid /= auxbas_grid) THEN
485 12 : pw_aux => pw_pools(auxbas_grid)%pool
486 12 : CALL pw_aux%create_pw(coeff_gaux)
487 12 : CALL pw_transfer(v_rspace, coeff_gaux)
488 12 : CALL pw_copy(coeff_gaux, coeff_gspace)
489 12 : CALL pw_aux%give_back_pw(coeff_gaux)
490 : ELSE
491 : ! Preserve the reciprocal-space projection on spherical grids.
492 0 : CALL pw_transfer(v_rspace, coeff_gspace)
493 : END IF
494 12 : CALL pw_transfer(coeff_gspace, coeff_rspace)
495 12 : CALL pw_pool%give_back_pw(coeff_gspace)
496 12 : IF (igrid /= auxbas_grid) THEN
497 12 : fscale = coeff_rspace%pw_grid%dvol/pw_aux%pw_grid%dvol
498 12 : CALL pw_scale(coeff_rspace, fscale)
499 : END IF
500 : ELSE
501 0 : CALL pw_copy(v_rspace, coeff_rspace)
502 : END IF
503 : END IF
504 :
505 : ! Setup the rs grid at level igrid
506 29254 : CALL rs_grid_create(rs_v, rs_desc)
507 29254 : CALL rs_grid_zero(rs_v)
508 29254 : IF (direct_transfer) THEN
509 : ! The potential is already on the target grid; no PW copy is needed.
510 29242 : CALL transfer_pw2rs(rs_v, v_rspace)
511 : ELSE
512 12 : CALL transfer_pw2rs(rs_v, coeff_rspace)
513 12 : CALL pw_pool%give_back_pw(coeff_rspace)
514 : END IF
515 :
516 : ! Now the potential is on the right grid => integration
517 :
518 29254 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
519 :
520 : ! Allocate work storage
521 :
522 29254 : NULLIFY (hab, hab_sph, hdab, hdab_sph, pab, a_hdab, a_hdab_sph)
523 29254 : nch_max = ncoset(lmax0)
524 29254 : CALL reallocate(hab, 1, nch_max, 1, 1)
525 29254 : CALL reallocate(hab_sph, 1, nch_max)
526 29254 : CALL reallocate(hdab, 1, 3, 1, nch_max, 1, 1)
527 29254 : CALL reallocate(hdab_sph, 1, 3, 1, nch_max)
528 29254 : CALL reallocate(a_hdab, 1, 3, 1, 3, 1, nch_max, 1, 1)
529 29254 : CALL reallocate(a_hdab_sph, 1, 3, 1, 3, 1, nch_max)
530 29254 : CALL reallocate(pab, 1, nch_max, 1, 1)
531 :
532 29254 : ncurr = -1
533 :
534 29254 : grid_distributed = rs_v%desc%distributed
535 :
536 29254 : fscale = 1.0_dp
537 29254 : IF (PRESENT(kforce)) THEN
538 88 : fscale = kforce
539 : END IF
540 :
541 84612 : DO ikind = 1, SIZE(atomic_kind_set, 1)
542 55358 : NULLIFY (basis_1c_set, atom_list, harmonics, cneo_potential)
543 55358 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
544 : CALL get_qs_kind(qs_kind_set(ikind), &
545 : basis_set=basis_1c_set, basis_type="GAPW_1C", &
546 : paw_atom=paw_atom, &
547 : harmonics=harmonics, &
548 55358 : cneo_potential=cneo_potential)
549 :
550 55358 : IF (.NOT. paw_atom) CYCLE
551 :
552 52004 : NULLIFY (Qlm_gg, lmax, npgf)
553 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, &
554 : l0_ikind=l0_ikind, Qlm_gg=Qlm_gg, & ! Qs different
555 52004 : rpgf0_s=rpgf0)
556 :
557 : CALL get_gto_basis_set(gto_basis_set=basis_1c_set, &
558 : lmax=lmax, lmin=lmin, &
559 : maxso=maxso, maxl=maxl, &
560 52004 : nset=nset, npgf=npgf)
561 :
562 52004 : nsotot = maxso*nset
563 208016 : ALLOCATE (intloc(nsotot, nsotot))
564 :
565 : ! Initialize the local KS integrals
566 :
567 52004 : nch_ik = ncoset(l0_ikind)
568 654154 : pab = 1.0_dp
569 52004 : max_s_harm = harmonics%max_s_harm
570 52004 : llmax = harmonics%llmax
571 :
572 52004 : NULLIFY (intloc_nuc)
573 52004 : maxl_nuc = -1
574 52004 : max_s_harm_nuc = 0
575 52004 : llmax_nuc = -1
576 52004 : IF (ASSOCIATED(cneo_potential)) THEN
577 48 : NULLIFY (nuc_basis_set)
578 : CALL get_qs_kind(qs_kind_set(ikind), &
579 : basis_set=nuc_basis_set, &
580 48 : basis_type="NUC")
581 :
582 : CALL get_gto_basis_set(gto_basis_set=nuc_basis_set, &
583 : lmax=lmax_nuc, lmin=lmin_nuc, &
584 : maxso=maxso_nuc, maxl=maxl_nuc, &
585 48 : nset=nset_nuc, npgf=npgf_nuc)
586 48 : nsotot_nuc = maxso_nuc*nset_nuc
587 192 : ALLOCATE (intloc_nuc(nsotot_nuc, nsotot_nuc))
588 48 : max_s_harm_nuc = cneo_potential%harmonics%max_s_harm
589 96 : llmax_nuc = cneo_potential%harmonics%llmax
590 : END IF
591 :
592 0 : ALLOCATE (cg_list(2, nsoset(MAX(maxl, maxl_nuc))**2, &
593 : MAX(max_s_harm, max_s_harm_nuc)), &
594 312024 : cg_n_list(MAX(max_s_harm, max_s_harm_nuc)))
595 :
596 52004 : num_pe = para_env%num_pe
597 52004 : mepos = para_env%mepos
598 156012 : DO j = 0, num_pe - 1
599 104008 : bo = get_limit(nat, num_pe, j)
600 104008 : IF (.NOT. grid_distributed .AND. j /= mepos) CYCLE
601 :
602 145153 : DO iat = bo(1), bo(2)
603 41145 : iatom = atom_list(iat)
604 41145 : ra(:) = pbc(particle_set(iatom)%r, cell)
605 :
606 41145 : NULLIFY (Qlm)
607 41145 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, iat=iatom, Qlm_tot=Qlm)
608 :
609 507254 : hab = 0.0_dp
610 1782146 : hdab = 0.0_dp
611 137891159 : intloc = 0._dp
612 41145 : IF (use_virial) THEN
613 949 : my_virial_a = 0.0_dp
614 949 : my_virial_b = 0.0_dp
615 119652 : a_hdab = 0.0_dp
616 : END IF
617 :
618 : CALL integrate_pgf_product( &
619 : l0_ikind, zet0, 0, 0, 0.0_dp, 0, &
620 : ra, [0.0_dp, 0.0_dp, 0.0_dp], rs_v, &
621 : hab, pab, o1=0, o2=0, &
622 : radius=rpgf0, &
623 : calculate_forces=calculate_forces, &
624 : use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b, &
625 41145 : hdab=hdab, a_hdab=a_hdab, use_subpatch=.TRUE., subpatch_pattern=0)
626 :
627 : ! Convert from cartesian to spherical
628 150350 : DO lshell = 0, l0_ikind
629 468399 : DO is = 1, nso(lshell)
630 318049 : iso = is + nsoset(lshell - 1)
631 318049 : hab_sph(iso) = 0.0_dp
632 1272196 : hdab_sph(1:3, iso) = 0.0_dp
633 4134637 : a_hdab_sph(1:3, 1:3, iso) = 0.0_dp
634 1892347 : DO ic = 1, nco(lshell)
635 1465093 : ico = ic + ncoset(lshell - 1)
636 1465093 : lx = indco(1, ico)
637 1465093 : ly = indco(2, ico)
638 1465093 : lz = indco(3, ico)
639 : hab_sph(iso) = hab_sph(iso) + &
640 : norm_l(lshell)* &
641 : orbtramat(lshell)%slm(is, ic)* &
642 1465093 : hab(ico, 1)
643 1465093 : IF (calculate_forces) THEN
644 : hdab_sph(1:3, iso) = hdab_sph(1:3, iso) + &
645 : norm_l(lshell)* &
646 : orbtramat(lshell)%slm(is, ic)* &
647 541824 : hdab(1:3, ico, 1)
648 : END IF
649 1783142 : IF (use_virial) THEN
650 142792 : DO ii = 1, 3
651 464074 : DO i = 1, 3
652 : a_hdab_sph(i, ii, iso) = a_hdab_sph(i, ii, iso) + &
653 : norm_l(lshell)* &
654 : orbtramat(lshell)%slm(is, ic)* &
655 428376 : a_hdab(i, ii, ico, 1)
656 : END DO
657 : END DO
658 : END IF
659 :
660 : END DO ! ic
661 : END DO ! is
662 : END DO ! lshell
663 :
664 : m1 = 0
665 143391 : DO iset1 = 1, nset
666 :
667 : m2 = 0
668 439816 : DO iset2 = 1, nset
669 : CALL get_none0_cg_list(harmonics%my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
670 337570 : max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
671 337570 : n1 = nsoset(lmax(iset1))
672 1187948 : DO ipgf1 = 1, npgf(iset1)
673 850378 : n2 = nsoset(lmax(iset2))
674 3666475 : DO ipgf2 = 1, npgf(iset2)
675 :
676 16538809 : DO iso = 1, MIN(nsoset(l0_ikind), max_iso_not0_local)
677 39980413 : DO icg = 1, cg_n_list(iso)
678 24291982 : iso1 = cg_list(1, icg, iso)
679 24291982 : iso2 = cg_list(2, icg, iso)
680 :
681 24291982 : ig1 = iso1 + n1*(ipgf1 - 1) + m1
682 24291982 : ig2 = iso2 + n2*(ipgf2 - 1) + m2
683 :
684 37501886 : intloc(ig1, ig2) = intloc(ig1, ig2) + Qlm_gg(ig1, ig2, iso)*hab_sph(iso) ! potential times Q
685 :
686 : END DO ! icg
687 : END DO ! iso
688 :
689 : END DO ! ipgf2
690 : END DO ! ipgf1
691 777386 : m2 = m2 + maxso
692 : END DO ! iset2
693 143391 : m1 = m1 + maxso
694 : END DO ! iset1
695 :
696 41145 : IF (ASSOCIATED(cneo_potential)) THEN
697 305578 : intloc_nuc = 0.0_dp
698 46 : m1 = 0
699 460 : DO iset1 = 1, nset_nuc
700 414 : n1 = nsoset(lmax_nuc(iset1))
701 414 : m2 = 0
702 4140 : DO iset2 = 1, nset_nuc
703 3726 : n2 = nsoset(lmax_nuc(iset2))
704 : CALL get_none0_cg_list(cneo_potential%harmonics%my_CG, lmin_nuc(iset1), &
705 : lmax_nuc(iset1), lmin_nuc(iset2), lmax_nuc(iset2), &
706 : max_s_harm_nuc, llmax_nuc, cg_list, cg_n_list, &
707 3726 : max_iso_not0_local)
708 7452 : DO ipgf1 = 1, npgf_nuc(iset1)
709 11178 : DO ipgf2 = 1, npgf_nuc(iset2)
710 :
711 29578 : DO iso = 1, MIN(nsoset(l0_ikind), max_iso_not0_local)
712 50968 : DO icg = 1, cg_n_list(iso)
713 25116 : iso1 = cg_list(1, icg, iso)
714 25116 : iso2 = cg_list(2, icg, iso)
715 :
716 25116 : ig1 = iso1 + n1*(ipgf1 - 1) + m1
717 25116 : ig2 = iso2 + n2*(ipgf2 - 1) + m2
718 :
719 : intloc_nuc(ig1, ig2) = intloc_nuc(ig1, ig2) - cneo_potential%zeff* &
720 47242 : cneo_potential%Qlm_gg(ig1, ig2, iso)*hab_sph(iso)
721 :
722 : END DO ! icg
723 : END DO ! iso
724 :
725 : END DO ! ipgf2
726 : END DO ! ipgf1
727 7866 : m2 = m2 + maxso_nuc
728 : END DO ! iset2
729 460 : m1 = m1 + maxso_nuc
730 : END DO ! iset1
731 : END IF
732 :
733 41145 : IF (grid_distributed) THEN
734 : ! Sum result over all processors
735 0 : CALL para_env%sum(intloc)
736 0 : IF (ASSOCIATED(cneo_potential)) THEN
737 0 : CALL para_env%sum(intloc_nuc)
738 : END IF
739 : END IF
740 :
741 41145 : IF (j == mepos) THEN
742 41145 : rho_atom => rho_atom_set(iatom)
743 41145 : CALL get_rho_atom(rho_atom=rho_atom, ga_Vlocal_gb_h=int_local_h, ga_Vlocal_gb_s=int_local_s)
744 88301 : DO ispin = 1, nspins
745 305007172 : int_local_h(ispin)%r_coef = int_local_h(ispin)%r_coef + intloc
746 305048317 : int_local_s(ispin)%r_coef = int_local_s(ispin)%r_coef + intloc
747 : END DO
748 41145 : IF (ASSOCIATED(cneo_potential)) THEN
749 46 : rhoz_cneo => rhoz_cneo_set(iatom)
750 611156 : rhoz_cneo%ga_Vlocal_gb_h = rhoz_cneo%ga_Vlocal_gb_h + intloc_nuc
751 611156 : rhoz_cneo%ga_Vlocal_gb_s = rhoz_cneo%ga_Vlocal_gb_s + intloc_nuc
752 : END IF
753 : END IF
754 :
755 41145 : IF (PRESENT(atener)) THEN
756 178 : DO iso = 1, nsoset(l0_ikind)
757 178 : atener(iatom) = atener(iatom) + 0.5_dp*Qlm(iso)*hab_sph(iso)
758 : END DO
759 : END IF
760 :
761 41145 : IF (calculate_forces) THEN
762 1662 : force_tmp(1:3) = 0.0_dp
763 15860 : DO iso = 1, nsoset(l0_ikind)
764 14198 : force_tmp(1) = force_tmp(1) + Qlm(iso)*hdab_sph(1, iso) ! Q here is from local_rho_set
765 14198 : force_tmp(2) = force_tmp(2) + Qlm(iso)*hdab_sph(2, iso)
766 15860 : force_tmp(3) = force_tmp(3) + Qlm(iso)*hdab_sph(3, iso)
767 : END DO
768 6648 : force(ikind)%g0s_Vh_elec(1:3, iat) = force(ikind)%g0s_Vh_elec(1:3, iat) + fscale*force_tmp(1:3)
769 : END IF
770 145153 : IF (use_virial) THEN
771 949 : my_virial_a = 0.0_dp
772 9026 : DO iso = 1, nsoset(l0_ikind)
773 33257 : DO ii = 1, 3
774 105001 : DO i = 1, 3
775 : ! Q from local_rho_set
776 72693 : virial%pv_gapw(i, ii) = virial%pv_gapw(i, ii) + fscale*Qlm(iso)*a_hdab_sph(i, ii, iso)
777 96924 : virial%pv_virial(i, ii) = virial%pv_virial(i, ii) + fscale*Qlm(iso)*a_hdab_sph(i, ii, iso)
778 : END DO
779 : END DO
780 : END DO
781 : END IF
782 :
783 : END DO
784 : END DO
785 :
786 52004 : DEALLOCATE (intloc)
787 52004 : IF (ASSOCIATED(intloc_nuc)) DEALLOCATE (intloc_nuc)
788 191974 : DEALLOCATE (cg_list, cg_n_list)
789 :
790 : END DO ! ikind
791 :
792 29254 : CALL rs_grid_release(rs_v)
793 :
794 29254 : DEALLOCATE (hab, hdab, hab_sph, hdab_sph, pab, a_hdab, a_hdab_sph)
795 :
796 29254 : CALL timestop(handle)
797 :
798 29254 : IF (rho0_mpole%do_cneo) THEN
799 48 : IF (PRESENT(kforce)) THEN
800 : CALL integrate_vhgg_rspace(qs_env, v_rspace, para_env, calculate_forces, &
801 0 : rhoz_cneo_set, kforce)
802 : ELSE
803 : CALL integrate_vhgg_rspace(qs_env, v_rspace, para_env, calculate_forces, &
804 48 : rhoz_cneo_set)
805 : END IF
806 : END IF
807 :
808 58508 : END SUBROUTINE integrate_vhg0_rspace
809 :
810 : END MODULE qs_rho0_ggrid
|
