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 : ! **************************************************************************************************
9 : !> \brief Calculate the plane wave density by collocating the primitive Gaussian
10 : !> functions (pgf).
11 : !> \par History
12 : !> - rewrote collocate for increased accuracy and speed
13 : !> - introduced the PGI hack for increased speed with that compiler
14 : !> (22.02.02)
15 : !> - Added Multiple Grid feature
16 : !> - new way to get over the grid (01.03.02)
17 : !> - removed timing calls since they were getting expensive
18 : !> - Updated with the new QS data structures (09.04.02,MK)
19 : !> - introduction of the real space grid type ( prelim. version JVdV 05.02)
20 : !> - parallel FFT (JGH 22.05.02)
21 : !> - multigrid arrays independent from density (JGH 30.08.02)
22 : !> - old density stored in g space (JGH 30.08.02)
23 : !> - distributed real space code (JGH 17.07.03)
24 : !> - refactoring and new loop ordering (JGH 23.11.03)
25 : !> - OpenMP parallelization (JGH 03.12.03)
26 : !> - Modified to compute tau (Joost 12.03)
27 : !> - removed the incremental density rebuild (Joost 01.04)
28 : !> - introduced realspace multigridding (Joost 02.04)
29 : !> - introduced map_consistent (Joost 02.04)
30 : !> - Addition of the subroutine calculate_atomic_charge_density (TdK, 08.05)
31 : !> - rewrite of the collocate/integrate kernels (Joost VandeVondele, 03.07)
32 : !> - Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
33 : !> \author Matthias Krack (03.04.2001)
34 : !> 1) Joost VandeVondele (01.2002)
35 : !> Thomas D. Kuehne (04.08.2005)
36 : !> Ole Schuett (2020)
37 : ! **************************************************************************************************
38 : MODULE qs_collocate_density
39 : USE admm_types, ONLY: get_admm_env
40 : USE ao_util, ONLY: exp_radius_very_extended
41 : USE atomic_kind_types, ONLY: atomic_kind_type, &
42 : get_atomic_kind, &
43 : get_atomic_kind_set
44 : USE basis_set_types, ONLY: get_gto_basis_set, &
45 : gto_basis_set_type
46 : USE cell_types, ONLY: cell_type, &
47 : pbc
48 : USE cp_control_types, ONLY: dft_control_type
49 : USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
50 : USE cp_fm_types, ONLY: cp_fm_get_element, &
51 : cp_fm_get_info, &
52 : cp_fm_type
53 : USE cp_dbcsr_api, ONLY: dbcsr_copy, &
54 : dbcsr_get_block_p, &
55 : dbcsr_p_type, &
56 : dbcsr_type
57 : USE external_potential_types, ONLY: get_potential, &
58 : gth_potential_type
59 : USE gaussian_gridlevels, ONLY: gaussian_gridlevel, &
60 : gridlevel_info_type
61 : USE grid_api, ONLY: &
62 : GRID_FUNC_AB, GRID_FUNC_CORE_X, GRID_FUNC_CORE_Y, GRID_FUNC_CORE_Z, GRID_FUNC_DAB_X, &
63 : GRID_FUNC_DAB_Y, GRID_FUNC_DAB_Z, GRID_FUNC_DABpADB_X, GRID_FUNC_DABpADB_Y, &
64 : GRID_FUNC_DABpADB_Z, GRID_FUNC_DADB, GRID_FUNC_DX, GRID_FUNC_DXDX, GRID_FUNC_DXDY, &
65 : GRID_FUNC_DY, GRID_FUNC_DYDY, GRID_FUNC_DYDZ, GRID_FUNC_DZ, GRID_FUNC_DZDX, &
66 : GRID_FUNC_DZDZ, collocate_pgf_product, grid_collocate_task_list
67 : USE input_constants, ONLY: &
68 : orb_dx2, orb_dxy, orb_dy2, orb_dyz, orb_dz2, orb_dzx, orb_px, orb_py, orb_pz, orb_s
69 : USE kinds, ONLY: default_string_length, &
70 : dp
71 : USE lri_environment_types, ONLY: lri_kind_type
72 : USE memory_utilities, ONLY: reallocate
73 : USE message_passing, ONLY: mp_comm_type
74 : USE orbital_pointers, ONLY: coset, &
75 : ncoset
76 : USE particle_types, ONLY: particle_type
77 : USE pw_env_types, ONLY: pw_env_get, &
78 : pw_env_type
79 : USE pw_methods, ONLY: pw_axpy, &
80 : pw_integrate_function, &
81 : pw_transfer, &
82 : pw_zero
83 : USE pw_pool_types, ONLY: pw_pool_p_type, &
84 : pw_pool_type, &
85 : pw_pools_create_pws, &
86 : pw_pools_give_back_pws
87 : USE pw_types, ONLY: pw_r3d_rs_type, &
88 : pw_c1d_gs_type, &
89 : pw_r3d_rs_type
90 : USE qs_environment_types, ONLY: get_qs_env, &
91 : qs_environment_type
92 : USE qs_kind_types, ONLY: get_qs_kind, &
93 : get_qs_kind_set, &
94 : qs_kind_type
95 : USE qs_ks_types, ONLY: get_ks_env, &
96 : qs_ks_env_type
97 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
98 : USE realspace_grid_types, ONLY: map_gaussian_here, &
99 : realspace_grid_desc_p_type, &
100 : realspace_grid_type, &
101 : rs_grid_zero, &
102 : transfer_rs2pw
103 : USE rs_pw_interface, ONLY: density_rs2pw
104 : USE task_list_methods, ONLY: rs_copy_to_buffer, &
105 : rs_distribute_matrix, &
106 : rs_scatter_matrices
107 : USE task_list_types, ONLY: atom_pair_type, &
108 : task_list_type, &
109 : task_type
110 :
111 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
112 :
113 : #include "./base/base_uses.f90"
114 :
115 : IMPLICIT NONE
116 :
117 : PRIVATE
118 :
119 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_collocate_density'
120 : ! *** Public subroutines ***
121 :
122 : PUBLIC :: calculate_ppl_grid, &
123 : calculate_rho_core, &
124 : calculate_lri_rho_elec, &
125 : calculate_rho_single_gaussian, &
126 : calculate_rho_metal, &
127 : calculate_rho_resp_single, &
128 : calculate_rho_resp_all, &
129 : calculate_rho_elec, &
130 : calculate_drho_elec, &
131 : calculate_wavefunction, &
132 : collocate_function, &
133 : calculate_rho_nlcc, &
134 : calculate_drho_elec_dR, &
135 : calculate_drho_core, &
136 : collocate_single_gaussian
137 :
138 : INTERFACE calculate_rho_core
139 : MODULE PROCEDURE calculate_rho_core_r3d_rs
140 : MODULE PROCEDURE calculate_rho_core_c1d_gs
141 : END INTERFACE
142 :
143 : INTERFACE calculate_rho_resp_all
144 : MODULE PROCEDURE calculate_rho_resp_all_r3d_rs, calculate_rho_resp_all_c1d_gs
145 : END INTERFACE
146 :
147 : CONTAINS
148 :
149 : ! **************************************************************************************************
150 : !> \brief computes the density of the non-linear core correction on the grid
151 : !> \param rho_nlcc ...
152 : !> \param qs_env ...
153 : ! **************************************************************************************************
154 36 : SUBROUTINE calculate_rho_nlcc(rho_nlcc, qs_env)
155 :
156 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho_nlcc
157 : TYPE(qs_environment_type), POINTER :: qs_env
158 :
159 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_nlcc'
160 :
161 : INTEGER :: atom_a, handle, iatom, iexp_nlcc, ikind, &
162 : ithread, j, n, natom, nc, nexp_nlcc, &
163 : ni, npme, nthread, subpatch_pattern
164 36 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores, nct_nlcc
165 : LOGICAL :: nlcc
166 : REAL(KIND=dp) :: alpha, eps_rho_rspace, radius
167 : REAL(KIND=dp), DIMENSION(3) :: ra
168 36 : REAL(KIND=dp), DIMENSION(:), POINTER :: alpha_nlcc
169 36 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_nlcc, pab
170 36 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
171 : TYPE(cell_type), POINTER :: cell
172 : TYPE(dft_control_type), POINTER :: dft_control
173 : TYPE(gth_potential_type), POINTER :: gth_potential
174 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
175 : TYPE(pw_env_type), POINTER :: pw_env
176 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
177 36 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
178 : TYPE(realspace_grid_type), POINTER :: rs_rho
179 :
180 36 : CALL timeset(routineN, handle)
181 :
182 36 : NULLIFY (cell, dft_control, pab, particle_set, atomic_kind_set, &
183 36 : qs_kind_set, atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
184 :
185 : CALL get_qs_env(qs_env=qs_env, &
186 : atomic_kind_set=atomic_kind_set, &
187 : qs_kind_set=qs_kind_set, &
188 : cell=cell, &
189 : dft_control=dft_control, &
190 : particle_set=particle_set, &
191 36 : pw_env=pw_env)
192 : CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
193 36 : auxbas_pw_pool=auxbas_pw_pool)
194 : ! be careful in parallel nsmax is chosen with multigrid in mind!
195 36 : CALL rs_grid_zero(rs_rho)
196 :
197 36 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
198 :
199 92 : DO ikind = 1, SIZE(atomic_kind_set)
200 56 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
201 56 : CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
202 :
203 56 : IF (.NOT. ASSOCIATED(gth_potential)) CYCLE
204 : CALL get_potential(potential=gth_potential, nlcc_present=nlcc, nexp_nlcc=nexp_nlcc, &
205 56 : alpha_nlcc=alpha_nlcc, nct_nlcc=nct_nlcc, cval_nlcc=cval_nlcc)
206 :
207 56 : IF (.NOT. nlcc) CYCLE
208 :
209 256 : DO iexp_nlcc = 1, nexp_nlcc
210 :
211 54 : alpha = alpha_nlcc(iexp_nlcc)
212 54 : nc = nct_nlcc(iexp_nlcc)
213 :
214 54 : ni = ncoset(2*nc - 2)
215 162 : ALLOCATE (pab(ni, 1))
216 306 : pab = 0._dp
217 :
218 54 : nthread = 1
219 54 : ithread = 0
220 :
221 54 : CALL reallocate(cores, 1, natom)
222 54 : npme = 0
223 232 : cores = 0
224 :
225 : ! prepare core function
226 124 : DO j = 1, nc
227 54 : SELECT CASE (j)
228 : CASE (1)
229 54 : pab(1, 1) = cval_nlcc(1, iexp_nlcc)
230 : CASE (2)
231 16 : n = coset(2, 0, 0)
232 16 : pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
233 16 : n = coset(0, 2, 0)
234 16 : pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
235 16 : n = coset(0, 0, 2)
236 16 : pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
237 : CASE (3)
238 0 : n = coset(4, 0, 0)
239 0 : pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
240 0 : n = coset(0, 4, 0)
241 0 : pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
242 0 : n = coset(0, 0, 4)
243 0 : pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
244 0 : n = coset(2, 2, 0)
245 0 : pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
246 0 : n = coset(2, 0, 2)
247 0 : pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
248 0 : n = coset(0, 2, 2)
249 0 : pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
250 : CASE (4)
251 0 : n = coset(6, 0, 0)
252 0 : pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
253 0 : n = coset(0, 6, 0)
254 0 : pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
255 0 : n = coset(0, 0, 6)
256 0 : pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
257 0 : n = coset(4, 2, 0)
258 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
259 0 : n = coset(4, 0, 2)
260 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
261 0 : n = coset(2, 4, 0)
262 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
263 0 : n = coset(2, 0, 4)
264 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
265 0 : n = coset(0, 4, 2)
266 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
267 0 : n = coset(0, 2, 4)
268 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
269 0 : n = coset(2, 2, 2)
270 0 : pab(n, 1) = 6._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
271 : CASE DEFAULT
272 70 : CPABORT("")
273 : END SELECT
274 : END DO
275 54 : IF (dft_control%nspins == 2) pab = pab*0.5_dp
276 :
277 232 : DO iatom = 1, natom
278 178 : atom_a = atom_list(iatom)
279 178 : ra(:) = pbc(particle_set(atom_a)%r, cell)
280 232 : IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
281 : ! replicated realspace grid, split the atoms up between procs
282 178 : IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
283 89 : npme = npme + 1
284 89 : cores(npme) = iatom
285 : END IF
286 : ELSE
287 0 : npme = npme + 1
288 0 : cores(npme) = iatom
289 : END IF
290 : END DO
291 :
292 143 : DO j = 1, npme
293 :
294 89 : iatom = cores(j)
295 89 : atom_a = atom_list(iatom)
296 89 : ra(:) = pbc(particle_set(atom_a)%r, cell)
297 89 : subpatch_pattern = 0
298 89 : ni = 2*nc - 2
299 : radius = exp_radius_very_extended(la_min=0, la_max=ni, lb_min=0, lb_max=0, &
300 : ra=ra, rb=ra, rp=ra, &
301 : zetp=1/(2*alpha**2), eps=eps_rho_rspace, &
302 : pab=pab, o1=0, o2=0, & ! without map_consistent
303 89 : prefactor=1.0_dp, cutoff=0.0_dp)
304 :
305 : CALL collocate_pgf_product(ni, 1/(2*alpha**2), 0, 0, 0.0_dp, 0, ra, &
306 : [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, pab, 0, 0, rs_rho, &
307 : ga_gb_function=GRID_FUNC_AB, radius=radius, &
308 143 : use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
309 :
310 : END DO
311 :
312 110 : DEALLOCATE (pab)
313 :
314 : END DO
315 :
316 : END DO
317 :
318 36 : IF (ASSOCIATED(cores)) THEN
319 36 : DEALLOCATE (cores)
320 : END IF
321 :
322 36 : CALL transfer_rs2pw(rs_rho, rho_nlcc)
323 :
324 36 : CALL timestop(handle)
325 :
326 36 : END SUBROUTINE calculate_rho_nlcc
327 :
328 : ! **************************************************************************************************
329 : !> \brief computes the local pseudopotential (without erf term) on the grid
330 : !> \param vppl ...
331 : !> \param qs_env ...
332 : ! **************************************************************************************************
333 12 : SUBROUTINE calculate_ppl_grid(vppl, qs_env)
334 :
335 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vppl
336 : TYPE(qs_environment_type), POINTER :: qs_env
337 :
338 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ppl_grid'
339 :
340 : INTEGER :: atom_a, handle, iatom, ikind, ithread, &
341 : j, lppl, n, natom, ni, npme, nthread, &
342 : subpatch_pattern
343 12 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
344 : REAL(KIND=dp) :: alpha, eps_rho_rspace, radius
345 : REAL(KIND=dp), DIMENSION(3) :: ra
346 12 : REAL(KIND=dp), DIMENSION(:), POINTER :: cexp_ppl
347 12 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
348 12 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
349 : TYPE(cell_type), POINTER :: cell
350 : TYPE(dft_control_type), POINTER :: dft_control
351 : TYPE(gth_potential_type), POINTER :: gth_potential
352 12 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
353 : TYPE(pw_env_type), POINTER :: pw_env
354 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
355 12 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
356 : TYPE(realspace_grid_type), POINTER :: rs_rho
357 :
358 12 : CALL timeset(routineN, handle)
359 :
360 12 : NULLIFY (cell, dft_control, pab, atomic_kind_set, qs_kind_set, particle_set, &
361 12 : atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
362 :
363 : CALL get_qs_env(qs_env=qs_env, &
364 : atomic_kind_set=atomic_kind_set, &
365 : qs_kind_set=qs_kind_set, &
366 : cell=cell, &
367 : dft_control=dft_control, &
368 : particle_set=particle_set, &
369 12 : pw_env=pw_env)
370 : CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
371 12 : auxbas_pw_pool=auxbas_pw_pool)
372 : ! be careful in parallel nsmax is chosen with multigrid in mind!
373 12 : CALL rs_grid_zero(rs_rho)
374 :
375 12 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
376 :
377 28 : DO ikind = 1, SIZE(atomic_kind_set)
378 16 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
379 16 : CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
380 :
381 16 : IF (.NOT. ASSOCIATED(gth_potential)) CYCLE
382 16 : CALL get_potential(potential=gth_potential, alpha_ppl=alpha, nexp_ppl=lppl, cexp_ppl=cexp_ppl)
383 :
384 16 : IF (lppl <= 0) CYCLE
385 :
386 16 : ni = ncoset(2*lppl - 2)
387 48 : ALLOCATE (pab(ni, 1))
388 192 : pab = 0._dp
389 :
390 16 : nthread = 1
391 16 : ithread = 0
392 :
393 16 : CALL reallocate(cores, 1, natom)
394 16 : npme = 0
395 60 : cores = 0
396 :
397 : ! prepare core function
398 48 : DO j = 1, lppl
399 16 : SELECT CASE (j)
400 : CASE (1)
401 16 : pab(1, 1) = cexp_ppl(1)
402 : CASE (2)
403 16 : n = coset(2, 0, 0)
404 16 : pab(n, 1) = cexp_ppl(2)
405 16 : n = coset(0, 2, 0)
406 16 : pab(n, 1) = cexp_ppl(2)
407 16 : n = coset(0, 0, 2)
408 16 : pab(n, 1) = cexp_ppl(2)
409 : CASE (3)
410 0 : n = coset(4, 0, 0)
411 0 : pab(n, 1) = cexp_ppl(3)
412 0 : n = coset(0, 4, 0)
413 0 : pab(n, 1) = cexp_ppl(3)
414 0 : n = coset(0, 0, 4)
415 0 : pab(n, 1) = cexp_ppl(3)
416 0 : n = coset(2, 2, 0)
417 0 : pab(n, 1) = 2._dp*cexp_ppl(3)
418 0 : n = coset(2, 0, 2)
419 0 : pab(n, 1) = 2._dp*cexp_ppl(3)
420 0 : n = coset(0, 2, 2)
421 0 : pab(n, 1) = 2._dp*cexp_ppl(3)
422 : CASE (4)
423 0 : n = coset(6, 0, 0)
424 0 : pab(n, 1) = cexp_ppl(4)
425 0 : n = coset(0, 6, 0)
426 0 : pab(n, 1) = cexp_ppl(4)
427 0 : n = coset(0, 0, 6)
428 0 : pab(n, 1) = cexp_ppl(4)
429 0 : n = coset(4, 2, 0)
430 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
431 0 : n = coset(4, 0, 2)
432 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
433 0 : n = coset(2, 4, 0)
434 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
435 0 : n = coset(2, 0, 4)
436 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
437 0 : n = coset(0, 4, 2)
438 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
439 0 : n = coset(0, 2, 4)
440 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
441 0 : n = coset(2, 2, 2)
442 0 : pab(n, 1) = 6._dp*cexp_ppl(4)
443 : CASE DEFAULT
444 32 : CPABORT("")
445 : END SELECT
446 : END DO
447 :
448 60 : DO iatom = 1, natom
449 44 : atom_a = atom_list(iatom)
450 44 : ra(:) = pbc(particle_set(atom_a)%r, cell)
451 60 : IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
452 : ! replicated realspace grid, split the atoms up between procs
453 44 : IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
454 22 : npme = npme + 1
455 22 : cores(npme) = iatom
456 : END IF
457 : ELSE
458 0 : npme = npme + 1
459 0 : cores(npme) = iatom
460 : END IF
461 : END DO
462 :
463 16 : IF (npme > 0) THEN
464 36 : DO j = 1, npme
465 :
466 22 : iatom = cores(j)
467 22 : atom_a = atom_list(iatom)
468 22 : ra(:) = pbc(particle_set(atom_a)%r, cell)
469 22 : subpatch_pattern = 0
470 22 : ni = 2*lppl - 2
471 :
472 : radius = exp_radius_very_extended(la_min=0, la_max=ni, &
473 : lb_min=0, lb_max=0, &
474 : ra=ra, rb=ra, rp=ra, &
475 : zetp=alpha, eps=eps_rho_rspace, &
476 : pab=pab, o1=0, o2=0, & ! without map_consistent
477 22 : prefactor=1.0_dp, cutoff=0.0_dp)
478 :
479 : CALL collocate_pgf_product(ni, alpha, 0, 0, 0.0_dp, 0, ra, &
480 : [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, pab, 0, 0, rs_rho, &
481 : radius=radius, ga_gb_function=GRID_FUNC_AB, &
482 36 : use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
483 :
484 : END DO
485 : END IF
486 :
487 60 : DEALLOCATE (pab)
488 :
489 : END DO
490 :
491 12 : IF (ASSOCIATED(cores)) THEN
492 12 : DEALLOCATE (cores)
493 : END IF
494 :
495 12 : CALL transfer_rs2pw(rs_rho, vppl)
496 :
497 12 : CALL timestop(handle)
498 :
499 12 : END SUBROUTINE calculate_ppl_grid
500 :
501 : ! **************************************************************************************************
502 : !> \brief Collocates the fitted lri density on a grid.
503 : !> \param lri_rho_g ...
504 : !> \param lri_rho_r ...
505 : !> \param qs_env ...
506 : !> \param lri_coef ...
507 : !> \param total_rho ...
508 : !> \param basis_type ...
509 : !> \param exact_1c_terms ...
510 : !> \param pmat replicated block diagonal density matrix (optional)
511 : !> \param atomlist list of atoms to be included (optional)
512 : !> \par History
513 : !> 04.2013
514 : !> \author Dorothea Golze
515 : ! **************************************************************************************************
516 908 : SUBROUTINE calculate_lri_rho_elec(lri_rho_g, lri_rho_r, qs_env, &
517 908 : lri_coef, total_rho, basis_type, exact_1c_terms, pmat, atomlist)
518 :
519 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: lri_rho_g
520 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: lri_rho_r
521 : TYPE(qs_environment_type), POINTER :: qs_env
522 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
523 : REAL(KIND=dp), INTENT(OUT) :: total_rho
524 : CHARACTER(len=*), INTENT(IN) :: basis_type
525 : LOGICAL, INTENT(IN) :: exact_1c_terms
526 : TYPE(dbcsr_type), OPTIONAL :: pmat
527 : INTEGER, DIMENSION(:), OPTIONAL :: atomlist
528 :
529 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_lri_rho_elec'
530 :
531 : INTEGER :: atom_a, group_size, handle, iatom, igrid_level, ikind, ipgf, iset, jpgf, jset, &
532 : m1, maxco, maxsgf_set, my_pos, na1, natom, nb1, ncoa, ncob, nseta, offset, sgfa, sgfb
533 908 : INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, nsgfa
534 908 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
535 : LOGICAL :: found
536 908 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: map_it
537 908 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: map_it2
538 : REAL(KIND=dp) :: eps_rho_rspace, radius, zetp
539 : REAL(KIND=dp), DIMENSION(3) :: ra
540 908 : REAL(KIND=dp), DIMENSION(:), POINTER :: aci
541 908 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, pab, sphi_a, work, zeta
542 908 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
543 : TYPE(cell_type), POINTER :: cell
544 : TYPE(dft_control_type), POINTER :: dft_control
545 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
546 : TYPE(gto_basis_set_type), POINTER :: lri_basis_set, orb_basis_set
547 908 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
548 : TYPE(pw_env_type), POINTER :: pw_env
549 908 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
550 908 : TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_gspace
551 908 : TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_rspace
552 908 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
553 908 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
554 : TYPE(realspace_grid_type), POINTER :: rs_grid
555 :
556 908 : NULLIFY (aci, atomic_kind_set, qs_kind_set, atom_list, cell, &
557 908 : dft_control, first_sgfa, gridlevel_info, la_max, &
558 908 : la_min, lri_basis_set, npgfa, nsgfa, &
559 908 : pab, particle_set, pw_env, pw_pools, rs_grid, rs_rho, sphi_a, &
560 908 : work, zeta)
561 :
562 908 : CALL timeset(routineN, handle)
563 :
564 908 : IF (exact_1c_terms) THEN
565 48 : CPASSERT(PRESENT(pmat))
566 : END IF
567 :
568 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
569 : atomic_kind_set=atomic_kind_set, &
570 : cell=cell, particle_set=particle_set, &
571 : pw_env=pw_env, &
572 908 : dft_control=dft_control)
573 :
574 908 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
575 908 : gridlevel_info => pw_env%gridlevel_info
576 :
577 : ! *** set up the pw multi-grids *** !
578 908 : CPASSERT(ASSOCIATED(pw_env))
579 908 : CALL pw_env_get(pw_env=pw_env, rs_grids=rs_rho, pw_pools=pw_pools)
580 :
581 908 : CALL pw_pools_create_pws(pw_pools, mgrid_rspace)
582 :
583 908 : CALL pw_pools_create_pws(pw_pools, mgrid_gspace)
584 :
585 : ! *** set up the rs multi-grids *** !
586 4480 : DO igrid_level = 1, gridlevel_info%ngrid_levels
587 4480 : CALL rs_grid_zero(rs_rho(igrid_level))
588 : END DO
589 :
590 : !take maxco from the LRI basis set!
591 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
592 908 : maxco=maxco, basis_type=basis_type)
593 :
594 2724 : ALLOCATE (pab(maxco, 1))
595 908 : offset = 0
596 908 : my_pos = mgrid_rspace(1)%pw_grid%para%group%mepos
597 908 : group_size = mgrid_rspace(1)%pw_grid%para%group%num_pe
598 :
599 2718 : DO ikind = 1, SIZE(atomic_kind_set)
600 :
601 1810 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
602 1810 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type=basis_type)
603 :
604 : !Take the lri basis set here!
605 : CALL get_gto_basis_set(gto_basis_set=lri_basis_set, lmax=la_max, &
606 : lmin=la_min, zet=zeta, nset=nseta, npgf=npgfa, &
607 1810 : sphi=sphi_a, first_sgf=first_sgfa, nsgf_set=nsgfa)
608 :
609 7950 : DO iatom = 1, natom
610 3422 : atom_a = atom_list(iatom)
611 3422 : IF (PRESENT(ATOMLIST)) THEN
612 500 : IF (atomlist(atom_a) == 0) CYCLE
613 : END IF
614 3222 : ra(:) = pbc(particle_set(atom_a)%r, cell)
615 3222 : aci => lri_coef(ikind)%acoef(iatom, :)
616 :
617 48280 : m1 = MAXVAL(npgfa(1:nseta))
618 9666 : ALLOCATE (map_it(m1))
619 48280 : DO iset = 1, nseta
620 : ! collocate this set locally?
621 95420 : map_it = .FALSE.
622 95228 : DO ipgf = 1, npgfa(iset)
623 50170 : igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
624 50170 : rs_grid => rs_rho(igrid_level)
625 95228 : map_it(ipgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
626 : END DO
627 45058 : offset = offset + 1
628 :
629 73365 : IF (ANY(map_it(1:npgfa(iset)))) THEN
630 22529 : sgfa = first_sgfa(1, iset)
631 22529 : ncoa = npgfa(iset)*ncoset(la_max(iset))
632 22529 : m1 = sgfa + nsgfa(iset) - 1
633 67587 : ALLOCATE (work(nsgfa(iset), 1))
634 373257 : work(1:nsgfa(iset), 1) = aci(sgfa:m1)
635 668770 : pab = 0._dp
636 :
637 : CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), 1.0_dp, lri_basis_set%sphi(1, sgfa), &
638 : SIZE(lri_basis_set%sphi, 1), work(1, 1), SIZE(work, 1), 0.0_dp, pab(1, 1), &
639 22529 : SIZE(pab, 1))
640 :
641 47614 : DO ipgf = 1, npgfa(iset)
642 25085 : na1 = (ipgf - 1)*ncoset(la_max(iset))
643 25085 : igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
644 25085 : rs_grid => rs_rho(igrid_level)
645 47614 : IF (map_it(ipgf)) THEN
646 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
647 : lb_min=0, lb_max=0, &
648 : ra=ra, rb=ra, rp=ra, &
649 : zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
650 25085 : prefactor=1.0_dp, cutoff=1.0_dp)
651 :
652 : CALL collocate_pgf_product(la_max=la_max(iset), &
653 : zeta=zeta(ipgf, iset), &
654 : la_min=la_min(iset), &
655 : lb_max=0, zetb=0.0_dp, lb_min=0, &
656 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], &
657 : scale=1._dp, &
658 : pab=pab, o1=na1, o2=0, &
659 : rsgrid=rs_grid, &
660 : radius=radius, &
661 25085 : ga_gb_function=GRID_FUNC_AB)
662 : END IF
663 : END DO
664 22529 : DEALLOCATE (work)
665 : END IF
666 : END DO
667 5232 : DEALLOCATE (map_it)
668 : END DO
669 : END DO
670 :
671 908 : DEALLOCATE (pab)
672 :
673 : ! process the one-center terms
674 908 : IF (exact_1c_terms) THEN
675 : ! find maximum numbers
676 48 : offset = 0
677 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
678 : maxco=maxco, &
679 : maxsgf_set=maxsgf_set, &
680 48 : basis_type="ORB")
681 336 : ALLOCATE (pab(maxco, maxco), work(maxco, maxsgf_set))
682 :
683 144 : DO ikind = 1, SIZE(atomic_kind_set)
684 96 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
685 96 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
686 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, lmax=la_max, &
687 : lmin=la_min, zet=zeta, nset=nseta, npgf=npgfa, &
688 96 : sphi=sphi_a, first_sgf=first_sgfa, nsgf_set=nsgfa)
689 528 : DO iatom = 1, natom
690 288 : atom_a = atom_list(iatom)
691 288 : ra(:) = pbc(particle_set(atom_a)%r, cell)
692 288 : CALL dbcsr_get_block_p(matrix=pmat, row=atom_a, col=atom_a, BLOCK=p_block, found=found)
693 576 : m1 = MAXVAL(npgfa(1:nseta))
694 1152 : ALLOCATE (map_it2(m1, m1))
695 576 : DO iset = 1, nseta
696 864 : DO jset = 1, nseta
697 : ! processor mappint
698 16416 : map_it2 = .FALSE.
699 2304 : DO ipgf = 1, npgfa(iset)
700 16416 : DO jpgf = 1, npgfa(jset)
701 14112 : zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
702 14112 : igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
703 14112 : rs_grid => rs_rho(igrid_level)
704 16128 : map_it2(ipgf, jpgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
705 : END DO
706 : END DO
707 288 : offset = offset + 1
708 : !
709 8640 : IF (ANY(map_it2(1:npgfa(iset), 1:npgfa(jset)))) THEN
710 144 : ncoa = npgfa(iset)*ncoset(la_max(iset))
711 144 : sgfa = first_sgfa(1, iset)
712 144 : ncob = npgfa(jset)*ncoset(la_max(jset))
713 144 : sgfb = first_sgfa(1, jset)
714 : ! decontract density block
715 : CALL dgemm("N", "N", ncoa, nsgfa(jset), nsgfa(iset), &
716 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
717 : p_block(sgfa, sgfb), SIZE(p_block, 1), &
718 144 : 0.0_dp, work(1, 1), maxco)
719 : CALL dgemm("N", "T", ncoa, ncob, nsgfa(jset), &
720 : 1.0_dp, work(1, 1), maxco, &
721 : sphi_a(1, sgfb), SIZE(sphi_a, 1), &
722 144 : 0.0_dp, pab(1, 1), maxco)
723 1152 : DO ipgf = 1, npgfa(iset)
724 8208 : DO jpgf = 1, npgfa(jset)
725 7056 : zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
726 7056 : igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
727 7056 : rs_grid => rs_rho(igrid_level)
728 :
729 7056 : na1 = (ipgf - 1)*ncoset(la_max(iset))
730 7056 : nb1 = (jpgf - 1)*ncoset(la_max(jset))
731 :
732 8064 : IF (map_it2(ipgf, jpgf)) THEN
733 : radius = exp_radius_very_extended(la_min=la_min(iset), &
734 : la_max=la_max(iset), &
735 : lb_min=la_min(jset), &
736 : lb_max=la_max(jset), &
737 : ra=ra, rb=ra, rp=ra, &
738 : zetp=zetp, eps=eps_rho_rspace, &
739 7056 : prefactor=1.0_dp, cutoff=1.0_dp)
740 :
741 : CALL collocate_pgf_product( &
742 : la_max(iset), zeta(ipgf, iset), la_min(iset), &
743 : la_max(jset), zeta(jpgf, jset), la_min(jset), &
744 : ra, [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, pab, na1, nb1, &
745 : rs_grid, &
746 7056 : radius=radius, ga_gb_function=GRID_FUNC_AB)
747 : END IF
748 : END DO
749 : END DO
750 : END IF
751 : END DO
752 : END DO
753 672 : DEALLOCATE (map_it2)
754 : !
755 : END DO
756 : END DO
757 96 : DEALLOCATE (pab, work)
758 : END IF
759 :
760 908 : CALL pw_zero(lri_rho_g)
761 908 : CALL pw_zero(lri_rho_r)
762 :
763 4480 : DO igrid_level = 1, gridlevel_info%ngrid_levels
764 3572 : CALL pw_zero(mgrid_rspace(igrid_level))
765 : CALL transfer_rs2pw(rs=rs_rho(igrid_level), &
766 4480 : pw=mgrid_rspace(igrid_level))
767 : END DO
768 :
769 4480 : DO igrid_level = 1, gridlevel_info%ngrid_levels
770 3572 : CALL pw_zero(mgrid_gspace(igrid_level))
771 : CALL pw_transfer(mgrid_rspace(igrid_level), &
772 3572 : mgrid_gspace(igrid_level))
773 4480 : CALL pw_axpy(mgrid_gspace(igrid_level), lri_rho_g)
774 : END DO
775 908 : CALL pw_transfer(lri_rho_g, lri_rho_r)
776 908 : total_rho = pw_integrate_function(lri_rho_r, isign=-1)
777 :
778 : ! *** give back the multi-grids *** !
779 908 : CALL pw_pools_give_back_pws(pw_pools, mgrid_gspace)
780 908 : CALL pw_pools_give_back_pws(pw_pools, mgrid_rspace)
781 :
782 908 : CALL timestop(handle)
783 :
784 4540 : END SUBROUTINE calculate_lri_rho_elec
785 :
786 : #:for kind in ["r3d_rs", "c1d_gs"]
787 : ! **************************************************************************************************
788 : !> \brief computes the density of the core charges on the grid
789 : !> \param rho_core ...
790 : !> \param total_rho ...
791 : !> \param qs_env ...
792 : !> \param calpha ...
793 : !> \param ccore ...
794 : !> \param only_nopaw ...
795 : ! **************************************************************************************************
796 9466 : SUBROUTINE calculate_rho_core_${kind}$ (rho_core, total_rho, qs_env, calpha, ccore, only_nopaw)
797 :
798 : TYPE(pw_${kind}$_type), INTENT(INOUT) :: rho_core
799 : REAL(KIND=dp), INTENT(OUT) :: total_rho
800 : TYPE(qs_environment_type), POINTER :: qs_env
801 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: calpha, ccore
802 : LOGICAL, INTENT(IN), OPTIONAL :: only_nopaw
803 :
804 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_core'
805 :
806 : INTEGER :: atom_a, handle, iatom, ikind, ithread, &
807 : j, natom, npme, nthread, &
808 : subpatch_pattern
809 9466 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
810 : LOGICAL :: my_only_nopaw, paw_atom
811 : REAL(KIND=dp) :: alpha, eps_rho_rspace, radius
812 : REAL(KIND=dp), DIMENSION(3) :: ra
813 9466 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
814 9466 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
815 : TYPE(cell_type), POINTER :: cell
816 : TYPE(dft_control_type), POINTER :: dft_control
817 9466 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
818 : TYPE(pw_env_type), POINTER :: pw_env
819 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
820 : TYPE(pw_r3d_rs_type) :: rhoc_r
821 9466 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
822 : TYPE(realspace_grid_type), POINTER :: rs_rho
823 :
824 9466 : CALL timeset(routineN, handle)
825 9466 : NULLIFY (cell, dft_control, pab, atomic_kind_set, qs_kind_set, particle_set, &
826 9466 : atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
827 9466 : ALLOCATE (pab(1, 1))
828 :
829 9466 : my_only_nopaw = .FALSE.
830 9466 : IF (PRESENT(only_nopaw)) my_only_nopaw = only_nopaw
831 9466 : IF (PRESENT(calpha)) THEN
832 2 : CPASSERT(PRESENT(ccore))
833 : END IF
834 :
835 : CALL get_qs_env(qs_env=qs_env, &
836 : atomic_kind_set=atomic_kind_set, &
837 : qs_kind_set=qs_kind_set, &
838 : cell=cell, &
839 : dft_control=dft_control, &
840 : particle_set=particle_set, &
841 9466 : pw_env=pw_env)
842 : CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
843 9466 : auxbas_pw_pool=auxbas_pw_pool)
844 : ! be careful in parallel nsmax is chosen with multigrid in mind!
845 9466 : CALL rs_grid_zero(rs_rho)
846 :
847 9466 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
848 :
849 26347 : DO ikind = 1, SIZE(atomic_kind_set)
850 16881 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
851 16881 : IF (PRESENT(calpha)) THEN
852 4 : alpha = calpha(ikind)
853 4 : pab(1, 1) = ccore(ikind)
854 : ELSE
855 16877 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
856 16877 : IF (my_only_nopaw .AND. paw_atom) CYCLE
857 : CALL get_qs_kind(qs_kind_set(ikind), alpha_core_charge=alpha, &
858 16717 : ccore_charge=pab(1, 1))
859 : END IF
860 :
861 16721 : IF (my_only_nopaw .AND. paw_atom) CYCLE
862 16721 : IF (alpha == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE
863 :
864 16545 : nthread = 1
865 16545 : ithread = 0
866 :
867 16545 : CALL reallocate(cores, 1, natom)
868 16545 : npme = 0
869 53532 : cores = 0
870 :
871 53532 : DO iatom = 1, natom
872 53532 : IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
873 : ! replicated realspace grid, split the atoms up between procs
874 36160 : IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
875 18080 : npme = npme + 1
876 18080 : cores(npme) = iatom
877 : END IF
878 : ELSE
879 827 : npme = npme + 1
880 827 : cores(npme) = iatom
881 : END IF
882 : END DO
883 :
884 42892 : IF (npme > 0) THEN
885 32052 : DO j = 1, npme
886 :
887 18907 : iatom = cores(j)
888 18907 : atom_a = atom_list(iatom)
889 18907 : ra(:) = pbc(particle_set(atom_a)%r, cell)
890 18907 : subpatch_pattern = 0
891 : radius = exp_radius_very_extended(la_min=0, la_max=0, &
892 : lb_min=0, lb_max=0, &
893 : ra=ra, rb=ra, rp=ra, &
894 : zetp=alpha, eps=eps_rho_rspace, &
895 : pab=pab, o1=0, o2=0, & ! without map_consistent
896 18907 : prefactor=-1.0_dp, cutoff=0.0_dp)
897 :
898 : CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
899 : [0.0_dp, 0.0_dp, 0.0_dp], -1.0_dp, pab, 0, 0, rs_rho, &
900 : radius=radius, ga_gb_function=GRID_FUNC_AB, &
901 32052 : use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
902 :
903 : END DO
904 : END IF
905 :
906 : END DO
907 :
908 9466 : IF (ASSOCIATED(cores)) THEN
909 9458 : DEALLOCATE (cores)
910 : END IF
911 9466 : DEALLOCATE (pab)
912 :
913 9466 : CALL auxbas_pw_pool%create_pw(rhoc_r)
914 :
915 9466 : CALL transfer_rs2pw(rs_rho, rhoc_r)
916 :
917 9466 : total_rho = pw_integrate_function(rhoc_r, isign=-1)
918 :
919 9466 : CALL pw_transfer(rhoc_r, rho_core)
920 :
921 9466 : CALL auxbas_pw_pool%give_back_pw(rhoc_r)
922 :
923 9466 : CALL timestop(handle)
924 :
925 9466 : END SUBROUTINE calculate_rho_core_${kind}$
926 : #:endfor
927 :
928 : ! *****************************************************************************
929 : !> \brief Computes the derivative of the density of the core charges with
930 : !> respect to the nuclear coordinates on the grid.
931 : !> \param drho_core The resulting density derivative
932 : !> \param qs_env ...
933 : !> \param beta Derivative direction
934 : !> \param lambda Atom index
935 : !> \note SL November 2014, ED 2021
936 : ! **************************************************************************************************
937 216 : SUBROUTINE calculate_drho_core(drho_core, qs_env, beta, lambda)
938 :
939 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: drho_core
940 : TYPE(qs_environment_type), POINTER :: qs_env
941 : INTEGER, INTENT(IN) :: beta, lambda
942 :
943 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_drho_core'
944 :
945 : INTEGER :: atom_a, dabqadb_func, handle, iatom, &
946 : ikind, ithread, j, natom, npme, &
947 : nthread, subpatch_pattern
948 216 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
949 : REAL(KIND=dp) :: alpha, eps_rho_rspace, radius
950 : REAL(KIND=dp), DIMENSION(3) :: ra
951 216 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
952 216 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
953 : TYPE(cell_type), POINTER :: cell
954 : TYPE(dft_control_type), POINTER :: dft_control
955 216 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
956 : TYPE(pw_env_type), POINTER :: pw_env
957 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
958 : TYPE(pw_r3d_rs_type) :: rhoc_r
959 216 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
960 : TYPE(realspace_grid_type), POINTER :: rs_rho
961 :
962 216 : CALL timeset(routineN, handle)
963 216 : NULLIFY (cell, dft_control, pab, atomic_kind_set, qs_kind_set, particle_set, &
964 216 : atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
965 216 : ALLOCATE (pab(1, 1))
966 :
967 : CALL get_qs_env(qs_env=qs_env, &
968 : atomic_kind_set=atomic_kind_set, &
969 : qs_kind_set=qs_kind_set, &
970 : cell=cell, &
971 : dft_control=dft_control, &
972 : particle_set=particle_set, &
973 216 : pw_env=pw_env)
974 : CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
975 216 : auxbas_pw_pool=auxbas_pw_pool)
976 : ! be careful in parallel nsmax is chosen with multigrid in mind!
977 216 : CALL rs_grid_zero(rs_rho)
978 :
979 216 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
980 :
981 288 : SELECT CASE (beta)
982 : CASE (1)
983 72 : dabqadb_func = GRID_FUNC_CORE_X
984 : CASE (2)
985 72 : dabqadb_func = GRID_FUNC_CORE_Y
986 : CASE (3)
987 72 : dabqadb_func = GRID_FUNC_CORE_Z
988 : CASE DEFAULT
989 216 : CPABORT("invalid beta")
990 : END SELECT
991 648 : DO ikind = 1, SIZE(atomic_kind_set)
992 432 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
993 : CALL get_qs_kind(qs_kind_set(ikind), &
994 432 : alpha_core_charge=alpha, ccore_charge=pab(1, 1))
995 :
996 432 : IF (alpha == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE
997 :
998 432 : nthread = 1
999 432 : ithread = 0
1000 :
1001 432 : CALL reallocate(cores, 1, natom)
1002 432 : npme = 0
1003 1080 : cores = 0
1004 :
1005 1080 : DO iatom = 1, natom
1006 1080 : IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
1007 : ! replicated realspace grid, split the atoms up between procs
1008 648 : IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
1009 324 : npme = npme + 1
1010 324 : cores(npme) = iatom
1011 : END IF
1012 : ELSE
1013 0 : npme = npme + 1
1014 0 : cores(npme) = iatom
1015 : END IF
1016 : END DO
1017 :
1018 1080 : IF (npme > 0) THEN
1019 648 : DO j = 1, npme
1020 :
1021 324 : iatom = cores(j)
1022 324 : atom_a = atom_list(iatom)
1023 324 : IF (atom_a /= lambda) CYCLE
1024 108 : ra(:) = pbc(particle_set(atom_a)%r, cell)
1025 108 : subpatch_pattern = 0
1026 : radius = exp_radius_very_extended(la_min=0, la_max=0, &
1027 : lb_min=0, lb_max=0, &
1028 : ra=ra, rb=ra, rp=ra, &
1029 : zetp=alpha, eps=eps_rho_rspace, &
1030 : pab=pab, o1=0, o2=0, & ! without map_consistent
1031 108 : prefactor=-1.0_dp, cutoff=0.0_dp)
1032 :
1033 : CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
1034 : [0.0_dp, 0.0_dp, 0.0_dp], -1.0_dp, pab, 0, 0, rs_rho, &
1035 : radius=radius, ga_gb_function=dabqadb_func, &
1036 648 : use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
1037 :
1038 : END DO
1039 : END IF
1040 :
1041 : END DO
1042 :
1043 216 : IF (ASSOCIATED(cores)) THEN
1044 216 : DEALLOCATE (cores)
1045 : END IF
1046 216 : DEALLOCATE (pab)
1047 :
1048 216 : CALL auxbas_pw_pool%create_pw(rhoc_r)
1049 :
1050 216 : CALL transfer_rs2pw(rs_rho, rhoc_r)
1051 :
1052 216 : CALL pw_transfer(rhoc_r, drho_core)
1053 :
1054 216 : CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1055 :
1056 216 : CALL timestop(handle)
1057 :
1058 216 : END SUBROUTINE calculate_drho_core
1059 :
1060 : ! **************************************************************************************************
1061 : !> \brief collocate a single Gaussian on the grid
1062 : !> \param rho_gb charge density generated by a single gaussian
1063 : !> \param qs_env qs environment
1064 : !> \param iatom_in atom index
1065 : !> \par History
1066 : !> 12.2011 created
1067 : !> \author Dorothea Golze
1068 : ! **************************************************************************************************
1069 4 : SUBROUTINE calculate_rho_single_gaussian(rho_gb, qs_env, iatom_in)
1070 :
1071 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_gb
1072 : TYPE(qs_environment_type), POINTER :: qs_env
1073 : INTEGER, INTENT(IN) :: iatom_in
1074 :
1075 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_single_gaussian'
1076 :
1077 : INTEGER :: atom_a, handle, iatom, npme, &
1078 : subpatch_pattern
1079 : REAL(KIND=dp) :: eps_rho_rspace, radius
1080 : REAL(KIND=dp), DIMENSION(3) :: ra
1081 4 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
1082 : TYPE(cell_type), POINTER :: cell
1083 : TYPE(dft_control_type), POINTER :: dft_control
1084 : TYPE(pw_env_type), POINTER :: pw_env
1085 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1086 : TYPE(pw_r3d_rs_type) :: rhoc_r
1087 : TYPE(realspace_grid_type), POINTER :: rs_rho
1088 :
1089 4 : CALL timeset(routineN, handle)
1090 4 : NULLIFY (cell, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool)
1091 :
1092 4 : ALLOCATE (pab(1, 1))
1093 :
1094 : CALL get_qs_env(qs_env=qs_env, &
1095 : cell=cell, &
1096 : dft_control=dft_control, &
1097 4 : pw_env=pw_env)
1098 : CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
1099 4 : auxbas_pw_pool=auxbas_pw_pool)
1100 4 : CALL rs_grid_zero(rs_rho)
1101 :
1102 4 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1103 4 : pab(1, 1) = 1.0_dp
1104 4 : iatom = iatom_in
1105 :
1106 4 : npme = 0
1107 :
1108 4 : IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
1109 4 : IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
1110 : npme = npme + 1
1111 : END IF
1112 : ELSE
1113 : npme = npme + 1
1114 : END IF
1115 :
1116 : IF (npme > 0) THEN
1117 2 : atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom)
1118 2 : ra(:) = pbc(qs_env%qmmm_env_qm%image_charge_pot%particles_all(atom_a)%r, cell)
1119 2 : subpatch_pattern = 0
1120 : radius = exp_radius_very_extended(la_min=0, la_max=0, &
1121 : lb_min=0, lb_max=0, &
1122 : ra=ra, rb=ra, rp=ra, &
1123 : zetp=qs_env%qmmm_env_qm%image_charge_pot%eta, &
1124 : eps=eps_rho_rspace, &
1125 : pab=pab, o1=0, o2=0, & ! without map_consistent
1126 2 : prefactor=1.0_dp, cutoff=0.0_dp)
1127 :
1128 : CALL collocate_pgf_product(0, qs_env%qmmm_env_qm%image_charge_pot%eta, &
1129 : 0, 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, pab, 0, 0, rs_rho, &
1130 : radius=radius, ga_gb_function=GRID_FUNC_AB, &
1131 2 : use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
1132 : END IF
1133 :
1134 4 : DEALLOCATE (pab)
1135 :
1136 4 : CALL auxbas_pw_pool%create_pw(rhoc_r)
1137 :
1138 4 : CALL transfer_rs2pw(rs_rho, rhoc_r)
1139 :
1140 4 : CALL pw_transfer(rhoc_r, rho_gb)
1141 :
1142 4 : CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1143 :
1144 4 : CALL timestop(handle)
1145 :
1146 4 : END SUBROUTINE calculate_rho_single_gaussian
1147 :
1148 : ! **************************************************************************************************
1149 : !> \brief computes the image charge density on the grid (including coeffcients)
1150 : !> \param rho_metal image charge density
1151 : !> \param coeff expansion coefficients of the image charge density, i.e.
1152 : !> rho_metal=sum_a c_a*g_a
1153 : !> \param total_rho_metal total induced image charge density
1154 : !> \param qs_env qs environment
1155 : !> \par History
1156 : !> 01.2012 created
1157 : !> \author Dorothea Golze
1158 : ! **************************************************************************************************
1159 90 : SUBROUTINE calculate_rho_metal(rho_metal, coeff, total_rho_metal, qs_env)
1160 :
1161 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_metal
1162 : REAL(KIND=dp), DIMENSION(:), POINTER :: coeff
1163 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: total_rho_metal
1164 : TYPE(qs_environment_type), POINTER :: qs_env
1165 :
1166 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_metal'
1167 :
1168 : INTEGER :: atom_a, handle, iatom, j, natom, npme, &
1169 : subpatch_pattern
1170 90 : INTEGER, DIMENSION(:), POINTER :: cores
1171 : REAL(KIND=dp) :: eps_rho_rspace, radius
1172 : REAL(KIND=dp), DIMENSION(3) :: ra
1173 90 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
1174 : TYPE(cell_type), POINTER :: cell
1175 : TYPE(dft_control_type), POINTER :: dft_control
1176 : TYPE(pw_env_type), POINTER :: pw_env
1177 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1178 : TYPE(pw_r3d_rs_type) :: rhoc_r
1179 : TYPE(realspace_grid_type), POINTER :: rs_rho
1180 :
1181 90 : CALL timeset(routineN, handle)
1182 :
1183 90 : NULLIFY (cell, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool, cores)
1184 :
1185 90 : ALLOCATE (pab(1, 1))
1186 :
1187 : CALL get_qs_env(qs_env=qs_env, &
1188 : cell=cell, &
1189 : dft_control=dft_control, &
1190 90 : pw_env=pw_env)
1191 : CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
1192 90 : auxbas_pw_pool=auxbas_pw_pool)
1193 90 : CALL rs_grid_zero(rs_rho)
1194 :
1195 90 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1196 90 : pab(1, 1) = 1.0_dp
1197 :
1198 90 : natom = SIZE(qs_env%qmmm_env_qm%image_charge_pot%image_mm_list)
1199 :
1200 90 : CALL reallocate(cores, 1, natom)
1201 90 : npme = 0
1202 270 : cores = 0
1203 :
1204 270 : DO iatom = 1, natom
1205 270 : IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
1206 180 : IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
1207 90 : npme = npme + 1
1208 90 : cores(npme) = iatom
1209 : END IF
1210 : ELSE
1211 0 : npme = npme + 1
1212 0 : cores(npme) = iatom
1213 : END IF
1214 : END DO
1215 :
1216 90 : IF (npme > 0) THEN
1217 180 : DO j = 1, npme
1218 90 : iatom = cores(j)
1219 90 : atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom)
1220 90 : ra(:) = pbc(qs_env%qmmm_env_qm%image_charge_pot%particles_all(atom_a)%r, cell)
1221 90 : subpatch_pattern = 0
1222 : radius = exp_radius_very_extended(la_min=0, la_max=0, &
1223 : lb_min=0, lb_max=0, &
1224 : ra=ra, rb=ra, rp=ra, &
1225 : zetp=qs_env%qmmm_env_qm%image_charge_pot%eta, &
1226 : eps=eps_rho_rspace, &
1227 : pab=pab, o1=0, o2=0, & ! without map_consistent
1228 90 : prefactor=coeff(iatom), cutoff=0.0_dp)
1229 :
1230 : CALL collocate_pgf_product( &
1231 : 0, qs_env%qmmm_env_qm%image_charge_pot%eta, &
1232 : 0, 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], coeff(iatom), pab, 0, 0, rs_rho, &
1233 : radius=radius, ga_gb_function=GRID_FUNC_AB, &
1234 180 : use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
1235 : END DO
1236 : END IF
1237 :
1238 90 : DEALLOCATE (pab, cores)
1239 :
1240 90 : CALL auxbas_pw_pool%create_pw(rhoc_r)
1241 :
1242 90 : CALL transfer_rs2pw(rs_rho, rhoc_r)
1243 :
1244 90 : IF (PRESENT(total_rho_metal)) &
1245 : !minus sign: account for the fact that rho_metal has opposite sign
1246 90 : total_rho_metal = pw_integrate_function(rhoc_r, isign=-1)
1247 :
1248 90 : CALL pw_transfer(rhoc_r, rho_metal)
1249 90 : CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1250 :
1251 90 : CALL timestop(handle)
1252 :
1253 90 : END SUBROUTINE calculate_rho_metal
1254 :
1255 : ! **************************************************************************************************
1256 : !> \brief collocate a single Gaussian on the grid for periodic RESP fitting
1257 : !> \param rho_gb charge density generated by a single gaussian
1258 : !> \param qs_env qs environment
1259 : !> \param eta width of single Gaussian
1260 : !> \param iatom_in atom index
1261 : !> \par History
1262 : !> 06.2012 created
1263 : !> \author Dorothea Golze
1264 : ! **************************************************************************************************
1265 66 : SUBROUTINE calculate_rho_resp_single(rho_gb, qs_env, eta, iatom_in)
1266 :
1267 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_gb
1268 : TYPE(qs_environment_type), POINTER :: qs_env
1269 : REAL(KIND=dp), INTENT(IN) :: eta
1270 : INTEGER, INTENT(IN) :: iatom_in
1271 :
1272 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_resp_single'
1273 :
1274 : INTEGER :: handle, iatom, npme, subpatch_pattern
1275 : REAL(KIND=dp) :: eps_rho_rspace, radius
1276 : REAL(KIND=dp), DIMENSION(3) :: ra
1277 66 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
1278 : TYPE(cell_type), POINTER :: cell
1279 : TYPE(dft_control_type), POINTER :: dft_control
1280 66 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1281 : TYPE(pw_env_type), POINTER :: pw_env
1282 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1283 : TYPE(pw_r3d_rs_type) :: rhoc_r
1284 : TYPE(realspace_grid_type), POINTER :: rs_rho
1285 :
1286 66 : CALL timeset(routineN, handle)
1287 66 : NULLIFY (cell, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool, &
1288 66 : particle_set)
1289 :
1290 66 : ALLOCATE (pab(1, 1))
1291 :
1292 : CALL get_qs_env(qs_env=qs_env, &
1293 : cell=cell, &
1294 : dft_control=dft_control, &
1295 : particle_set=particle_set, &
1296 66 : pw_env=pw_env)
1297 : CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
1298 66 : auxbas_pw_pool=auxbas_pw_pool)
1299 66 : CALL rs_grid_zero(rs_rho)
1300 :
1301 66 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1302 66 : pab(1, 1) = 1.0_dp
1303 66 : iatom = iatom_in
1304 :
1305 66 : npme = 0
1306 :
1307 66 : IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
1308 66 : IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
1309 : npme = npme + 1
1310 : END IF
1311 : ELSE
1312 : npme = npme + 1
1313 : END IF
1314 :
1315 : IF (npme > 0) THEN
1316 33 : ra(:) = pbc(particle_set(iatom)%r, cell)
1317 33 : subpatch_pattern = 0
1318 : radius = exp_radius_very_extended(la_min=0, la_max=0, &
1319 : lb_min=0, lb_max=0, &
1320 : ra=ra, rb=ra, rp=ra, &
1321 : zetp=eta, eps=eps_rho_rspace, &
1322 : pab=pab, o1=0, o2=0, & ! without map_consistent
1323 33 : prefactor=1.0_dp, cutoff=0.0_dp)
1324 :
1325 : CALL collocate_pgf_product(0, eta, 0, 0, 0.0_dp, 0, ra, &
1326 : [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, pab, 0, 0, rs_rho, &
1327 : radius=radius, ga_gb_function=GRID_FUNC_AB, &
1328 33 : use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
1329 : END IF
1330 :
1331 66 : DEALLOCATE (pab)
1332 :
1333 66 : CALL auxbas_pw_pool%create_pw(rhoc_r)
1334 :
1335 66 : CALL transfer_rs2pw(rs_rho, rhoc_r)
1336 :
1337 66 : CALL pw_transfer(rhoc_r, rho_gb)
1338 :
1339 66 : CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1340 :
1341 66 : CALL timestop(handle)
1342 :
1343 66 : END SUBROUTINE calculate_rho_resp_single
1344 :
1345 : #:for kind in ["r3d_rs", "c1d_gs"]
1346 : ! **************************************************************************************************
1347 : !> \brief computes the RESP charge density on a grid based on the RESP charges
1348 : !> \param rho_resp RESP charge density
1349 : !> \param coeff RESP charges, take care of normalization factor
1350 : !> (eta/pi)**1.5 later
1351 : !> \param natom number of atoms
1352 : !> \param eta width of single Gaussian
1353 : !> \param qs_env qs environment
1354 : !> \par History
1355 : !> 01.2012 created
1356 : !> \author Dorothea Golze
1357 : ! **************************************************************************************************
1358 24 : SUBROUTINE calculate_rho_resp_all_${kind}$ (rho_resp, coeff, natom, eta, qs_env)
1359 :
1360 : TYPE(pw_${kind}$_type), INTENT(INOUT) :: rho_resp
1361 : REAL(KIND=dp), DIMENSION(:), POINTER :: coeff
1362 : INTEGER, INTENT(IN) :: natom
1363 : REAL(KIND=dp), INTENT(IN) :: eta
1364 : TYPE(qs_environment_type), POINTER :: qs_env
1365 :
1366 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_resp_all'
1367 :
1368 : INTEGER :: handle, iatom, j, npme, subpatch_pattern
1369 24 : INTEGER, DIMENSION(:), POINTER :: cores
1370 : REAL(KIND=dp) :: eps_rho_rspace, radius
1371 : REAL(KIND=dp), DIMENSION(3) :: ra
1372 24 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
1373 : TYPE(cell_type), POINTER :: cell
1374 : TYPE(dft_control_type), POINTER :: dft_control
1375 24 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1376 : TYPE(pw_env_type), POINTER :: pw_env
1377 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1378 : TYPE(pw_r3d_rs_type) :: rhoc_r
1379 : TYPE(realspace_grid_type), POINTER :: rs_rho
1380 :
1381 24 : CALL timeset(routineN, handle)
1382 :
1383 24 : NULLIFY (cell, cores, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool, &
1384 24 : particle_set)
1385 :
1386 24 : ALLOCATE (pab(1, 1))
1387 :
1388 : CALL get_qs_env(qs_env=qs_env, &
1389 : cell=cell, &
1390 : dft_control=dft_control, &
1391 : particle_set=particle_set, &
1392 24 : pw_env=pw_env)
1393 : CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
1394 24 : auxbas_pw_pool=auxbas_pw_pool)
1395 24 : CALL rs_grid_zero(rs_rho)
1396 :
1397 24 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1398 24 : pab(1, 1) = 1.0_dp
1399 :
1400 24 : CALL reallocate(cores, 1, natom)
1401 24 : npme = 0
1402 142 : cores = 0
1403 :
1404 142 : DO iatom = 1, natom
1405 142 : IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
1406 118 : IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
1407 59 : npme = npme + 1
1408 59 : cores(npme) = iatom
1409 : END IF
1410 : ELSE
1411 0 : npme = npme + 1
1412 0 : cores(npme) = iatom
1413 : END IF
1414 : END DO
1415 :
1416 24 : IF (npme > 0) THEN
1417 83 : DO j = 1, npme
1418 59 : iatom = cores(j)
1419 59 : ra(:) = pbc(particle_set(iatom)%r, cell)
1420 59 : subpatch_pattern = 0
1421 : radius = exp_radius_very_extended(la_min=0, la_max=0, &
1422 : lb_min=0, lb_max=0, &
1423 : ra=ra, rb=ra, rp=ra, &
1424 : zetp=eta, eps=eps_rho_rspace, &
1425 : pab=pab, o1=0, o2=0, & ! without map_consistent
1426 59 : prefactor=coeff(iatom), cutoff=0.0_dp)
1427 :
1428 : CALL collocate_pgf_product( &
1429 : 0, eta, &
1430 : 0, 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], coeff(iatom), pab, 0, 0, rs_rho, &
1431 : radius=radius, ga_gb_function=GRID_FUNC_AB, &
1432 83 : use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
1433 : END DO
1434 : END IF
1435 :
1436 24 : DEALLOCATE (pab, cores)
1437 :
1438 24 : CALL auxbas_pw_pool%create_pw(rhoc_r)
1439 :
1440 24 : CALL transfer_rs2pw(rs_rho, rhoc_r)
1441 :
1442 24 : CALL pw_transfer(rhoc_r, rho_resp)
1443 24 : CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1444 :
1445 24 : CALL timestop(handle)
1446 :
1447 24 : END SUBROUTINE calculate_rho_resp_all_${kind}$
1448 : #:endfor
1449 :
1450 : ! **************************************************************************************************
1451 : !> \brief computes the density corresponding to a given density matrix on the grid
1452 : !> \param matrix_p ...
1453 : !> \param matrix_p_kp ...
1454 : !> \param rho ...
1455 : !> \param rho_gspace ...
1456 : !> \param total_rho ...
1457 : !> \param ks_env ...
1458 : !> \param soft_valid ...
1459 : !> \param compute_tau ...
1460 : !> \param compute_grad ...
1461 : !> \param basis_type ...
1462 : !> \param der_type ...
1463 : !> \param idir ...
1464 : !> \param task_list_external ...
1465 : !> \param pw_env_external ...
1466 : !> \par History
1467 : !> IAB (15-Feb-2010): Added OpenMP parallelisation to task loop
1468 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2010 on behalf of the HECToR project
1469 : !> Anything that is not the default ORB basis_type requires an external_task_list 12.2019, (A.Bussy)
1470 : !> Ole Schuett (2020): Migrated to C, see grid_api.F
1471 : !> \note
1472 : !> both rho and rho_gspace contain the new rho
1473 : !> (in real and g-space respectively)
1474 : ! **************************************************************************************************
1475 208186 : SUBROUTINE calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, &
1476 : ks_env, soft_valid, compute_tau, compute_grad, &
1477 : basis_type, der_type, idir, task_list_external, pw_env_external)
1478 :
1479 : TYPE(dbcsr_type), OPTIONAL, TARGET :: matrix_p
1480 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
1481 : POINTER :: matrix_p_kp
1482 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho
1483 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_gspace
1484 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: total_rho
1485 : TYPE(qs_ks_env_type), POINTER :: ks_env
1486 : LOGICAL, INTENT(IN), OPTIONAL :: soft_valid, compute_tau, compute_grad
1487 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: basis_type
1488 : INTEGER, INTENT(IN), OPTIONAL :: der_type, idir
1489 : TYPE(task_list_type), OPTIONAL, POINTER :: task_list_external
1490 : TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
1491 :
1492 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_elec'
1493 :
1494 : CHARACTER(LEN=default_string_length) :: my_basis_type
1495 : INTEGER :: ga_gb_function, handle, ilevel, img, &
1496 : nimages, nlevels
1497 : LOGICAL :: any_distributed, my_compute_grad, &
1498 : my_compute_tau, my_soft_valid
1499 208186 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_images
1500 : TYPE(dft_control_type), POINTER :: dft_control
1501 : TYPE(mp_comm_type) :: group
1502 : TYPE(pw_env_type), POINTER :: pw_env
1503 208186 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
1504 : TYPE(task_list_type), POINTER :: task_list
1505 :
1506 208186 : CALL timeset(routineN, handle)
1507 :
1508 208186 : NULLIFY (matrix_images, dft_control, pw_env, rs_rho, task_list)
1509 :
1510 : ! Figure out which function to collocate.
1511 208186 : my_compute_tau = .FALSE.
1512 208186 : IF (PRESENT(compute_tau)) my_compute_tau = compute_tau
1513 208186 : my_compute_grad = .FALSE.
1514 208186 : IF (PRESENT(compute_grad)) my_compute_grad = compute_grad
1515 208186 : IF (PRESENT(der_type)) THEN
1516 84 : SELECT CASE (der_type)
1517 : CASE (orb_s)
1518 36 : ga_gb_function = GRID_FUNC_AB
1519 : CASE (orb_px)
1520 0 : ga_gb_function = GRID_FUNC_DX
1521 : CASE (orb_py)
1522 0 : ga_gb_function = GRID_FUNC_DY
1523 : CASE (orb_pz)
1524 12 : ga_gb_function = GRID_FUNC_DZ
1525 : CASE (orb_dxy)
1526 0 : ga_gb_function = GRID_FUNC_DXDY
1527 : CASE (orb_dyz)
1528 0 : ga_gb_function = GRID_FUNC_DYDZ
1529 : CASE (orb_dzx)
1530 0 : ga_gb_function = GRID_FUNC_DZDX
1531 : CASE (orb_dx2)
1532 0 : ga_gb_function = GRID_FUNC_DXDX
1533 : CASE (orb_dy2)
1534 0 : ga_gb_function = GRID_FUNC_DYDY
1535 : CASE (orb_dz2)
1536 0 : ga_gb_function = GRID_FUNC_DZDZ
1537 : CASE DEFAULT
1538 48 : CPABORT("Unknown der_type")
1539 : END SELECT
1540 208138 : ELSE IF (my_compute_tau) THEN
1541 5060 : ga_gb_function = GRID_FUNC_DADB
1542 203078 : ELSE IF (my_compute_grad) THEN
1543 258 : CPASSERT(PRESENT(idir))
1544 344 : SELECT CASE (idir)
1545 : CASE (1)
1546 86 : ga_gb_function = GRID_FUNC_DABpADB_X
1547 : CASE (2)
1548 86 : ga_gb_function = GRID_FUNC_DABpADB_Y
1549 : CASE (3)
1550 86 : ga_gb_function = GRID_FUNC_DABpADB_Z
1551 : CASE DEFAULT
1552 258 : CPABORT("invalid idir")
1553 : END SELECT
1554 : ELSE
1555 202820 : ga_gb_function = GRID_FUNC_AB
1556 : END IF
1557 :
1558 : ! Figure out which basis_type to use.
1559 208186 : my_basis_type = "ORB" ! by default, the full density is calculated
1560 208186 : IF (PRESENT(basis_type)) my_basis_type = basis_type
1561 208186 : CPASSERT(my_basis_type == "ORB" .OR. PRESENT(task_list_external))
1562 :
1563 : ! Figure out which task_list to use.
1564 208186 : my_soft_valid = .FALSE.
1565 208186 : IF (PRESENT(soft_valid)) my_soft_valid = soft_valid
1566 208186 : IF (PRESENT(task_list_external)) THEN
1567 40272 : task_list => task_list_external
1568 167914 : ELSEIF (my_soft_valid) THEN
1569 26540 : CALL get_ks_env(ks_env, task_list_soft=task_list)
1570 : ELSE
1571 141374 : CALL get_ks_env(ks_env, task_list=task_list)
1572 : END IF
1573 208186 : CPASSERT(ASSOCIATED(task_list))
1574 :
1575 : ! Figure out which pw_env to use.
1576 208186 : IF (PRESENT(pw_env_external)) THEN
1577 22144 : pw_env => pw_env_external
1578 : ELSE
1579 186042 : CALL get_ks_env(ks_env, pw_env=pw_env)
1580 : END IF
1581 208186 : CPASSERT(ASSOCIATED(pw_env))
1582 :
1583 : ! Get grids.
1584 208186 : CALL pw_env_get(pw_env, rs_grids=rs_rho)
1585 208186 : nlevels = SIZE(rs_rho)
1586 208186 : group = rs_rho(1)%desc%group
1587 :
1588 : ! Check if any of the grids is distributed.
1589 208186 : any_distributed = .FALSE.
1590 1032580 : DO ilevel = 1, nlevels
1591 1856066 : any_distributed = any_distributed .OR. rs_rho(ilevel)%desc%distributed
1592 : END DO
1593 :
1594 : ! Gather all matrix images in a single array.
1595 208186 : CALL get_ks_env(ks_env, dft_control=dft_control)
1596 208186 : nimages = dft_control%nimages
1597 959004 : ALLOCATE (matrix_images(nimages))
1598 208186 : IF (PRESENT(matrix_p_kp)) THEN
1599 175932 : CPASSERT(.NOT. PRESENT(matrix_p))
1600 478124 : DO img = 1, nimages
1601 478124 : matrix_images(img)%matrix => matrix_p_kp(img)%matrix
1602 : END DO
1603 : ELSE
1604 32254 : CPASSERT(PRESENT(matrix_p) .AND. nimages == 1)
1605 32254 : matrix_images(1)%matrix => matrix_p
1606 : END IF
1607 :
1608 : ! Distribute matrix blocks.
1609 208186 : IF (any_distributed) THEN
1610 230 : CALL rs_scatter_matrices(matrix_images, task_list%pab_buffer, task_list, group)
1611 : ELSE
1612 207956 : CALL rs_copy_to_buffer(matrix_images, task_list%pab_buffer, task_list)
1613 : END IF
1614 208186 : DEALLOCATE (matrix_images)
1615 :
1616 : ! Map all tasks onto the grids
1617 : CALL grid_collocate_task_list(task_list=task_list%grid_task_list, &
1618 : ga_gb_function=ga_gb_function, &
1619 : pab_blocks=task_list%pab_buffer, &
1620 208186 : rs_grids=rs_rho)
1621 :
1622 : ! Merge realspace multi-grids into single planewave grid.
1623 208186 : CALL density_rs2pw(pw_env, rs_rho, rho, rho_gspace)
1624 208186 : IF (PRESENT(total_rho)) total_rho = pw_integrate_function(rho, isign=-1)
1625 :
1626 208186 : CALL timestop(handle)
1627 :
1628 208186 : END SUBROUTINE calculate_rho_elec
1629 :
1630 : ! **************************************************************************************************
1631 : !> \brief computes the gradient of the density corresponding to a given
1632 : !> density matrix on the grid
1633 : !> \param matrix_p ...
1634 : !> \param matrix_p_kp ...
1635 : !> \param drho ...
1636 : !> \param drho_gspace ...
1637 : !> \param qs_env ...
1638 : !> \param soft_valid ...
1639 : !> \param basis_type ...
1640 : !> \note this is an alternative to calculate the gradient through FFTs
1641 : ! **************************************************************************************************
1642 0 : SUBROUTINE calculate_drho_elec(matrix_p, matrix_p_kp, drho, drho_gspace, qs_env, &
1643 : soft_valid, basis_type)
1644 :
1645 : TYPE(dbcsr_type), OPTIONAL, TARGET :: matrix_p
1646 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
1647 : POINTER :: matrix_p_kp
1648 : TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT) :: drho
1649 : TYPE(pw_c1d_gs_type), DIMENSION(3), INTENT(INOUT) :: drho_gspace
1650 : TYPE(qs_environment_type), POINTER :: qs_env
1651 : LOGICAL, INTENT(IN), OPTIONAL :: soft_valid
1652 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: basis_type
1653 :
1654 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_drho_elec'
1655 :
1656 : CHARACTER(LEN=default_string_length) :: my_basis_type
1657 : INTEGER :: bcol, brow, dabqadb_func, handle, iatom, iatom_old, idir, igrid_level, ikind, &
1658 : ikind_old, img, img_old, ipgf, iset, iset_old, itask, ithread, jatom, jatom_old, jkind, &
1659 : jkind_old, jpgf, jset, jset_old, maxco, maxsgf_set, na1, na2, natoms, nb1, nb2, ncoa, &
1660 : ncob, nimages, nseta, nsetb, ntasks, nthread, sgfa, sgfb
1661 0 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
1662 0 : npgfb, nsgfa, nsgfb
1663 0 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1664 : LOGICAL :: atom_pair_changed, distributed_rs_grids, &
1665 : do_kp, found, my_soft, use_subpatch
1666 : REAL(KIND=dp) :: eps_rho_rspace, f, prefactor, radius, &
1667 : scale, zetp
1668 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rab_inv, rb, rp
1669 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, pab, sphi_a, sphi_b, work, &
1670 0 : zeta, zetb
1671 0 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pabt, workt
1672 0 : TYPE(atom_pair_type), DIMENSION(:), POINTER :: atom_pair_recv, atom_pair_send
1673 : TYPE(cell_type), POINTER :: cell
1674 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: deltap
1675 : TYPE(dft_control_type), POINTER :: dft_control
1676 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
1677 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1678 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1679 0 : POINTER :: sab_orb
1680 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1681 : TYPE(pw_env_type), POINTER :: pw_env
1682 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1683 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
1684 0 : POINTER :: rs_descs
1685 0 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
1686 : TYPE(task_list_type), POINTER :: task_list, task_list_soft
1687 0 : TYPE(task_type), DIMENSION(:), POINTER :: tasks
1688 :
1689 0 : CALL timeset(routineN, handle)
1690 :
1691 0 : CPASSERT(PRESENT(matrix_p) .OR. PRESENT(matrix_p_kp))
1692 0 : do_kp = PRESENT(matrix_p_kp)
1693 :
1694 0 : NULLIFY (cell, dft_control, orb_basis_set, deltap, qs_kind_set, &
1695 0 : sab_orb, particle_set, rs_rho, pw_env, rs_descs, la_max, la_min, &
1696 0 : lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, p_block, sphi_a, &
1697 0 : sphi_b, zeta, zetb, first_sgfa, first_sgfb, tasks, pabt, workt)
1698 :
1699 : ! by default, the full density is calculated
1700 0 : my_soft = .FALSE.
1701 0 : IF (PRESENT(soft_valid)) my_soft = soft_valid
1702 :
1703 0 : IF (PRESENT(basis_type)) THEN
1704 0 : my_basis_type = basis_type
1705 : ELSE
1706 0 : my_basis_type = "ORB"
1707 : END IF
1708 :
1709 : CALL get_qs_env(qs_env=qs_env, &
1710 : qs_kind_set=qs_kind_set, &
1711 : cell=cell, &
1712 : dft_control=dft_control, &
1713 : particle_set=particle_set, &
1714 : sab_orb=sab_orb, &
1715 0 : pw_env=pw_env)
1716 :
1717 0 : SELECT CASE (my_basis_type)
1718 : CASE ("ORB")
1719 : CALL get_qs_env(qs_env=qs_env, &
1720 : task_list=task_list, &
1721 0 : task_list_soft=task_list_soft)
1722 : CASE ("AUX_FIT")
1723 : CALL get_qs_env(qs_env=qs_env, &
1724 0 : task_list_soft=task_list_soft)
1725 0 : CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list)
1726 : END SELECT
1727 :
1728 : ! *** assign from pw_env
1729 0 : gridlevel_info => pw_env%gridlevel_info
1730 :
1731 : ! *** Allocate work storage ***
1732 0 : nthread = 1
1733 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
1734 : maxco=maxco, &
1735 : maxsgf_set=maxsgf_set, &
1736 0 : basis_type=my_basis_type)
1737 0 : CALL reallocate(pabt, 1, maxco, 1, maxco, 0, nthread - 1)
1738 0 : CALL reallocate(workt, 1, maxco, 1, maxsgf_set, 0, nthread - 1)
1739 :
1740 : ! find maximum numbers
1741 0 : nimages = dft_control%nimages
1742 0 : CPASSERT(nimages == 1 .OR. do_kp)
1743 :
1744 0 : natoms = SIZE(particle_set)
1745 :
1746 : ! get the task lists
1747 0 : IF (my_soft) task_list => task_list_soft
1748 0 : CPASSERT(ASSOCIATED(task_list))
1749 0 : tasks => task_list%tasks
1750 0 : atom_pair_send => task_list%atom_pair_send
1751 0 : atom_pair_recv => task_list%atom_pair_recv
1752 0 : ntasks = task_list%ntasks
1753 :
1754 : ! *** set up the rs multi-grids
1755 0 : CPASSERT(ASSOCIATED(pw_env))
1756 0 : CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho)
1757 0 : DO igrid_level = 1, gridlevel_info%ngrid_levels
1758 0 : distributed_rs_grids = rs_rho(igrid_level)%desc%distributed
1759 : END DO
1760 :
1761 0 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1762 :
1763 : ! *** Initialize working density matrix ***
1764 : ! distributed rs grids require a matrix that will be changed
1765 : ! whereas this is not the case for replicated grids
1766 0 : ALLOCATE (deltap(nimages))
1767 0 : IF (distributed_rs_grids) THEN
1768 0 : DO img = 1, nimages
1769 : END DO
1770 : ! this matrix has no strict sparsity pattern in parallel
1771 : ! deltap%sparsity_id=-1
1772 0 : IF (do_kp) THEN
1773 0 : DO img = 1, nimages
1774 : CALL dbcsr_copy(deltap(img)%matrix, matrix_p_kp(img)%matrix, &
1775 0 : name="DeltaP")
1776 : END DO
1777 : ELSE
1778 0 : CALL dbcsr_copy(deltap(1)%matrix, matrix_p, name="DeltaP")
1779 : END IF
1780 : ELSE
1781 0 : IF (do_kp) THEN
1782 0 : DO img = 1, nimages
1783 0 : deltap(img)%matrix => matrix_p_kp(img)%matrix
1784 : END DO
1785 : ELSE
1786 0 : deltap(1)%matrix => matrix_p
1787 : END IF
1788 : END IF
1789 :
1790 : ! distribute the matrix
1791 0 : IF (distributed_rs_grids) THEN
1792 : CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltap, &
1793 : atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
1794 0 : nimages=nimages, scatter=.TRUE.)
1795 : END IF
1796 :
1797 : ! map all tasks on the grids
1798 :
1799 0 : ithread = 0
1800 0 : pab => pabt(:, :, ithread)
1801 0 : work => workt(:, :, ithread)
1802 :
1803 0 : loop_xyz: DO idir = 1, 3
1804 :
1805 0 : DO igrid_level = 1, gridlevel_info%ngrid_levels
1806 0 : CALL rs_grid_zero(rs_rho(igrid_level))
1807 : END DO
1808 :
1809 : iatom_old = -1; jatom_old = -1; iset_old = -1; jset_old = -1
1810 : ikind_old = -1; jkind_old = -1; img_old = -1
1811 0 : loop_tasks: DO itask = 1, ntasks
1812 :
1813 : !decode the atom pair and basis info
1814 0 : igrid_level = tasks(itask)%grid_level
1815 0 : img = tasks(itask)%image
1816 0 : iatom = tasks(itask)%iatom
1817 0 : jatom = tasks(itask)%jatom
1818 0 : iset = tasks(itask)%iset
1819 0 : jset = tasks(itask)%jset
1820 0 : ipgf = tasks(itask)%ipgf
1821 0 : jpgf = tasks(itask)%jpgf
1822 :
1823 0 : ikind = particle_set(iatom)%atomic_kind%kind_number
1824 0 : jkind = particle_set(jatom)%atomic_kind%kind_number
1825 :
1826 0 : IF (iatom /= iatom_old .OR. jatom /= jatom_old .OR. img /= img_old) THEN
1827 :
1828 0 : IF (iatom /= iatom_old) ra(:) = pbc(particle_set(iatom)%r, cell)
1829 :
1830 0 : IF (iatom <= jatom) THEN
1831 0 : brow = iatom
1832 0 : bcol = jatom
1833 : ELSE
1834 0 : brow = jatom
1835 0 : bcol = iatom
1836 : END IF
1837 :
1838 0 : IF (ikind /= ikind_old) THEN
1839 0 : IF (my_soft) THEN
1840 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
1841 0 : basis_type="ORB_SOFT")
1842 : ELSE
1843 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
1844 0 : basis_type=my_basis_type)
1845 : END IF
1846 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1847 : first_sgf=first_sgfa, &
1848 : lmax=la_max, &
1849 : lmin=la_min, &
1850 : npgf=npgfa, &
1851 : nset=nseta, &
1852 : nsgf_set=nsgfa, &
1853 : sphi=sphi_a, &
1854 0 : zet=zeta)
1855 : END IF
1856 :
1857 0 : IF (jkind /= jkind_old) THEN
1858 0 : IF (my_soft) THEN
1859 : CALL get_qs_kind(qs_kind_set(jkind), basis_set=orb_basis_set, &
1860 0 : basis_type="ORB_SOFT")
1861 : ELSE
1862 : CALL get_qs_kind(qs_kind_set(jkind), basis_set=orb_basis_set, &
1863 0 : basis_type=my_basis_type)
1864 : END IF
1865 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1866 : first_sgf=first_sgfb, &
1867 : lmax=lb_max, &
1868 : lmin=lb_min, &
1869 : npgf=npgfb, &
1870 : nset=nsetb, &
1871 : nsgf_set=nsgfb, &
1872 : sphi=sphi_b, &
1873 0 : zet=zetb)
1874 : END IF
1875 :
1876 : CALL dbcsr_get_block_p(matrix=deltap(img)%matrix, &
1877 0 : row=brow, col=bcol, BLOCK=p_block, found=found)
1878 0 : CPASSERT(found)
1879 :
1880 : iatom_old = iatom
1881 : jatom_old = jatom
1882 : ikind_old = ikind
1883 : jkind_old = jkind
1884 : img_old = img
1885 : atom_pair_changed = .TRUE.
1886 :
1887 : ELSE
1888 :
1889 : atom_pair_changed = .FALSE.
1890 :
1891 : END IF
1892 :
1893 0 : IF (atom_pair_changed .OR. iset_old /= iset .OR. jset_old /= jset) THEN
1894 :
1895 0 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1896 0 : sgfa = first_sgfa(1, iset)
1897 0 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1898 0 : sgfb = first_sgfb(1, jset)
1899 :
1900 0 : IF (iatom <= jatom) THEN
1901 : CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
1902 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1903 : p_block(sgfa, sgfb), SIZE(p_block, 1), &
1904 0 : 0.0_dp, work(1, 1), maxco)
1905 : CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
1906 : 1.0_dp, work(1, 1), maxco, &
1907 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1908 0 : 0.0_dp, pab(1, 1), maxco)
1909 : ELSE
1910 : CALL dgemm("N", "N", ncob, nsgfa(iset), nsgfb(jset), &
1911 : 1.0_dp, sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1912 : p_block(sgfb, sgfa), SIZE(p_block, 1), &
1913 0 : 0.0_dp, work(1, 1), maxco)
1914 : CALL dgemm("N", "T", ncob, ncoa, nsgfa(iset), &
1915 : 1.0_dp, work(1, 1), maxco, &
1916 : sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1917 0 : 0.0_dp, pab(1, 1), maxco)
1918 : END IF
1919 :
1920 : iset_old = iset
1921 : jset_old = jset
1922 :
1923 : END IF
1924 :
1925 0 : rab(:) = tasks(itask)%rab
1926 0 : rb(:) = ra(:) + rab(:)
1927 0 : zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
1928 :
1929 0 : f = zetb(jpgf, jset)/zetp
1930 0 : rp(:) = ra(:) + f*rab(:)
1931 0 : prefactor = EXP(-zeta(ipgf, iset)*f*DOT_PRODUCT(rab, rab))
1932 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
1933 : lb_min=lb_min(jset), lb_max=lb_max(jset), &
1934 : ra=ra, rb=rb, rp=rp, &
1935 : zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
1936 0 : prefactor=prefactor, cutoff=1.0_dp)
1937 :
1938 0 : na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
1939 0 : na2 = ipgf*ncoset(la_max(iset))
1940 0 : nb1 = (jpgf - 1)*ncoset(lb_max(jset)) + 1
1941 0 : nb2 = jpgf*ncoset(lb_max(jset))
1942 :
1943 : ! takes the density matrix symmetry in account, i.e. off-diagonal blocks need to be mapped 'twice'
1944 0 : IF (iatom == jatom .AND. img == 1) THEN
1945 0 : scale = 1.0_dp
1946 : ELSE
1947 0 : scale = 2.0_dp
1948 : END IF
1949 :
1950 : ! check whether we need to use fawzi's generalised collocation scheme
1951 0 : IF (rs_rho(igrid_level)%desc%distributed) THEN
1952 : !tasks(4,:) is 0 for replicated, 1 for distributed 2 for exceptional distributed tasks
1953 0 : IF (tasks(itask)%dist_type == 2) THEN
1954 0 : use_subpatch = .TRUE.
1955 : ELSE
1956 0 : use_subpatch = .FALSE.
1957 : END IF
1958 : ELSE
1959 0 : use_subpatch = .FALSE.
1960 : END IF
1961 :
1962 0 : SELECT CASE (idir)
1963 : CASE (1)
1964 0 : dabqadb_func = GRID_FUNC_DABpADB_X
1965 : CASE (2)
1966 0 : dabqadb_func = GRID_FUNC_DABpADB_Y
1967 : CASE (3)
1968 0 : dabqadb_func = GRID_FUNC_DABpADB_Z
1969 : CASE DEFAULT
1970 0 : CPABORT("invalid idir")
1971 : END SELECT
1972 :
1973 0 : IF (iatom <= jatom) THEN
1974 : CALL collocate_pgf_product( &
1975 : la_max(iset), zeta(ipgf, iset), la_min(iset), &
1976 : lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1977 : ra, rab, scale, pab, na1 - 1, nb1 - 1, &
1978 : rs_rho(igrid_level), &
1979 : radius=radius, ga_gb_function=dabqadb_func, &
1980 0 : use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern)
1981 : ELSE
1982 0 : rab_inv = -rab
1983 : CALL collocate_pgf_product( &
1984 : lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1985 : la_max(iset), zeta(ipgf, iset), la_min(iset), &
1986 : rb, rab_inv, scale, pab, nb1 - 1, na1 - 1, &
1987 : rs_rho(igrid_level), &
1988 : radius=radius, ga_gb_function=dabqadb_func, &
1989 0 : use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern)
1990 : END IF
1991 :
1992 : END DO loop_tasks
1993 :
1994 0 : CALL density_rs2pw(pw_env, rs_rho, drho(idir), drho_gspace(idir))
1995 :
1996 : END DO loop_xyz
1997 :
1998 : ! *** Release work storage ***
1999 0 : IF (distributed_rs_grids) THEN
2000 0 : CALL dbcsr_deallocate_matrix_set(deltap)
2001 : ELSE
2002 0 : DO img = 1, nimages
2003 0 : NULLIFY (deltap(img)%matrix)
2004 : END DO
2005 0 : DEALLOCATE (deltap)
2006 : END IF
2007 :
2008 0 : DEALLOCATE (pabt, workt)
2009 :
2010 0 : CALL timestop(handle)
2011 :
2012 0 : END SUBROUTINE calculate_drho_elec
2013 :
2014 : ! **************************************************************************************************
2015 : !> \brief Computes the gradient wrt. nuclear coordinates of a density on the grid
2016 : !> The density is given in terms of the density matrix_p
2017 : !> \param matrix_p Density matrix
2018 : !> \param matrix_p_kp ...
2019 : !> \param drho Density gradient on the grid
2020 : !> \param drho_gspace Density gradient on the reciprocal grid
2021 : !> \param qs_env ...
2022 : !> \param soft_valid ...
2023 : !> \param basis_type ...
2024 : !> \param beta Derivative direction
2025 : !> \param lambda Atom index
2026 : !> \note SL, ED 2021
2027 : !> Adapted from calculate_drho_elec
2028 : ! **************************************************************************************************
2029 252 : SUBROUTINE calculate_drho_elec_dR(matrix_p, matrix_p_kp, drho, drho_gspace, qs_env, &
2030 : soft_valid, basis_type, beta, lambda)
2031 :
2032 : TYPE(dbcsr_type), OPTIONAL, TARGET :: matrix_p
2033 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
2034 : POINTER :: matrix_p_kp
2035 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: drho
2036 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: drho_gspace
2037 : TYPE(qs_environment_type), POINTER :: qs_env
2038 : LOGICAL, INTENT(IN), OPTIONAL :: soft_valid
2039 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: basis_type
2040 : INTEGER, INTENT(IN) :: beta, lambda
2041 :
2042 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_drho_elec_dR'
2043 :
2044 : CHARACTER(LEN=default_string_length) :: my_basis_type
2045 : INTEGER :: bcol, brow, dabqadb_func, handle, iatom, iatom_old, igrid_level, ikind, &
2046 : ikind_old, img, img_old, ipgf, iset, iset_old, itask, ithread, jatom, jatom_old, jkind, &
2047 : jkind_old, jpgf, jset, jset_old, maxco, maxsgf_set, na1, na2, natoms, nb1, nb2, ncoa, &
2048 : ncob, nimages, nseta, nsetb, ntasks, nthread, sgfa, sgfb
2049 252 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
2050 252 : npgfb, nsgfa, nsgfb
2051 252 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
2052 : LOGICAL :: atom_pair_changed, distributed_rs_grids, &
2053 : do_kp, found, my_soft, use_subpatch
2054 : REAL(KIND=dp) :: eps_rho_rspace, f, prefactor, radius, &
2055 : scale, zetp
2056 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rab_inv, rb, rp
2057 252 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, pab, sphi_a, sphi_b, work, &
2058 252 : zeta, zetb
2059 252 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pabt, workt
2060 252 : TYPE(atom_pair_type), DIMENSION(:), POINTER :: atom_pair_recv, atom_pair_send
2061 : TYPE(cell_type), POINTER :: cell
2062 252 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: deltap
2063 : TYPE(dft_control_type), POINTER :: dft_control
2064 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
2065 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
2066 252 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2067 : TYPE(pw_env_type), POINTER :: pw_env
2068 252 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2069 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
2070 252 : POINTER :: rs_descs
2071 252 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
2072 : TYPE(task_list_type), POINTER :: task_list, task_list_soft
2073 252 : TYPE(task_type), DIMENSION(:), POINTER :: tasks
2074 :
2075 252 : CALL timeset(routineN, handle)
2076 :
2077 252 : CPASSERT(PRESENT(matrix_p) .OR. PRESENT(matrix_p_kp))
2078 252 : do_kp = PRESENT(matrix_p_kp)
2079 :
2080 252 : NULLIFY (cell, dft_control, orb_basis_set, deltap, qs_kind_set, &
2081 252 : particle_set, rs_rho, pw_env, rs_descs, la_max, la_min, lb_max, &
2082 252 : lb_min, npgfa, npgfb, nsgfa, nsgfb, p_block, sphi_a, sphi_b, &
2083 252 : zeta, zetb, first_sgfa, first_sgfb, tasks, pabt, workt)
2084 :
2085 : ! by default, the full density is calculated
2086 252 : my_soft = .FALSE.
2087 252 : IF (PRESENT(soft_valid)) my_soft = soft_valid
2088 :
2089 252 : IF (PRESENT(basis_type)) THEN
2090 0 : my_basis_type = basis_type
2091 : ELSE
2092 252 : my_basis_type = "ORB"
2093 : END IF
2094 :
2095 : CALL get_qs_env(qs_env=qs_env, &
2096 : qs_kind_set=qs_kind_set, &
2097 : cell=cell, &
2098 : dft_control=dft_control, &
2099 : particle_set=particle_set, &
2100 252 : pw_env=pw_env)
2101 :
2102 252 : SELECT CASE (my_basis_type)
2103 : CASE ("ORB")
2104 : CALL get_qs_env(qs_env=qs_env, &
2105 : task_list=task_list, &
2106 252 : task_list_soft=task_list_soft)
2107 : CASE ("AUX_FIT")
2108 : CALL get_qs_env(qs_env=qs_env, &
2109 0 : task_list_soft=task_list_soft)
2110 252 : CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list)
2111 : END SELECT
2112 :
2113 : ! *** assign from pw_env
2114 252 : gridlevel_info => pw_env%gridlevel_info
2115 :
2116 : ! *** Allocate work storage ***
2117 252 : nthread = 1
2118 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
2119 : maxco=maxco, &
2120 : maxsgf_set=maxsgf_set, &
2121 252 : basis_type=my_basis_type)
2122 252 : CALL reallocate(pabt, 1, maxco, 1, maxco, 0, nthread - 1)
2123 252 : CALL reallocate(workt, 1, maxco, 1, maxsgf_set, 0, nthread - 1)
2124 :
2125 : ! find maximum numbers
2126 252 : nimages = dft_control%nimages
2127 252 : CPASSERT(nimages == 1 .OR. do_kp)
2128 :
2129 252 : natoms = SIZE(particle_set)
2130 :
2131 : ! get the task lists
2132 252 : IF (my_soft) task_list => task_list_soft
2133 252 : CPASSERT(ASSOCIATED(task_list))
2134 252 : tasks => task_list%tasks
2135 252 : atom_pair_send => task_list%atom_pair_send
2136 252 : atom_pair_recv => task_list%atom_pair_recv
2137 252 : ntasks = task_list%ntasks
2138 :
2139 : ! *** set up the rs multi-grids
2140 252 : CPASSERT(ASSOCIATED(pw_env))
2141 252 : CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho)
2142 774 : DO igrid_level = 1, gridlevel_info%ngrid_levels
2143 774 : distributed_rs_grids = rs_rho(igrid_level)%desc%distributed
2144 : END DO
2145 :
2146 252 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
2147 :
2148 : ! *** Initialize working density matrix ***
2149 : ! distributed rs grids require a matrix that will be changed
2150 : ! whereas this is not the case for replicated grids
2151 1008 : ALLOCATE (deltap(nimages))
2152 252 : IF (distributed_rs_grids) THEN
2153 0 : DO img = 1, nimages
2154 : END DO
2155 : ! this matrix has no strict sparsity pattern in parallel
2156 : ! deltap%sparsity_id=-1
2157 0 : IF (do_kp) THEN
2158 0 : DO img = 1, nimages
2159 : CALL dbcsr_copy(deltap(img)%matrix, matrix_p_kp(img)%matrix, &
2160 0 : name="DeltaP")
2161 : END DO
2162 : ELSE
2163 0 : CALL dbcsr_copy(deltap(1)%matrix, matrix_p, name="DeltaP")
2164 : END IF
2165 : ELSE
2166 252 : IF (do_kp) THEN
2167 0 : DO img = 1, nimages
2168 0 : deltap(img)%matrix => matrix_p_kp(img)%matrix
2169 : END DO
2170 : ELSE
2171 252 : deltap(1)%matrix => matrix_p
2172 : END IF
2173 : END IF
2174 :
2175 : ! distribute the matrix
2176 252 : IF (distributed_rs_grids) THEN
2177 : CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltap, &
2178 : atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
2179 0 : nimages=nimages, scatter=.TRUE.)
2180 : END IF
2181 :
2182 : ! map all tasks on the grids
2183 :
2184 252 : ithread = 0
2185 252 : pab => pabt(:, :, ithread)
2186 252 : work => workt(:, :, ithread)
2187 :
2188 774 : DO igrid_level = 1, gridlevel_info%ngrid_levels
2189 774 : CALL rs_grid_zero(rs_rho(igrid_level))
2190 : END DO
2191 :
2192 : iatom_old = -1; jatom_old = -1; iset_old = -1; jset_old = -1
2193 : ikind_old = -1; jkind_old = -1; img_old = -1
2194 16506 : loop_tasks: DO itask = 1, ntasks
2195 :
2196 : !decode the atom pair and basis info
2197 16254 : igrid_level = tasks(itask)%grid_level
2198 16254 : img = tasks(itask)%image
2199 16254 : iatom = tasks(itask)%iatom
2200 16254 : jatom = tasks(itask)%jatom
2201 16254 : iset = tasks(itask)%iset
2202 16254 : jset = tasks(itask)%jset
2203 16254 : ipgf = tasks(itask)%ipgf
2204 16254 : jpgf = tasks(itask)%jpgf
2205 :
2206 16254 : ikind = particle_set(iatom)%atomic_kind%kind_number
2207 16254 : jkind = particle_set(jatom)%atomic_kind%kind_number
2208 :
2209 16254 : IF (iatom /= iatom_old .OR. jatom /= jatom_old .OR. img /= img_old) THEN
2210 :
2211 1296 : IF (iatom /= iatom_old) ra(:) = pbc(particle_set(iatom)%r, cell)
2212 :
2213 1296 : IF (iatom <= jatom) THEN
2214 864 : brow = iatom
2215 864 : bcol = jatom
2216 : ELSE
2217 432 : brow = jatom
2218 432 : bcol = iatom
2219 : END IF
2220 :
2221 1296 : IF (ikind /= ikind_old) THEN
2222 252 : IF (my_soft) THEN
2223 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
2224 0 : basis_type="ORB_SOFT")
2225 : ELSE
2226 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
2227 252 : basis_type=my_basis_type)
2228 : END IF
2229 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
2230 : first_sgf=first_sgfa, &
2231 : lmax=la_max, &
2232 : lmin=la_min, &
2233 : npgf=npgfa, &
2234 : nset=nseta, &
2235 : nsgf_set=nsgfa, &
2236 : sphi=sphi_a, &
2237 252 : zet=zeta)
2238 : END IF
2239 :
2240 1296 : IF (jkind /= jkind_old) THEN
2241 864 : IF (my_soft) THEN
2242 : CALL get_qs_kind(qs_kind_set(jkind), basis_set=orb_basis_set, &
2243 0 : basis_type="ORB_SOFT")
2244 : ELSE
2245 : CALL get_qs_kind(qs_kind_set(jkind), basis_set=orb_basis_set, &
2246 864 : basis_type=my_basis_type)
2247 : END IF
2248 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
2249 : first_sgf=first_sgfb, &
2250 : lmax=lb_max, &
2251 : lmin=lb_min, &
2252 : npgf=npgfb, &
2253 : nset=nsetb, &
2254 : nsgf_set=nsgfb, &
2255 : sphi=sphi_b, &
2256 864 : zet=zetb)
2257 : END IF
2258 :
2259 : CALL dbcsr_get_block_p(matrix=deltap(img)%matrix, &
2260 1296 : row=brow, col=bcol, BLOCK=p_block, found=found)
2261 1296 : CPASSERT(found)
2262 :
2263 : iatom_old = iatom
2264 : jatom_old = jatom
2265 : ikind_old = ikind
2266 : jkind_old = jkind
2267 : img_old = img
2268 : atom_pair_changed = .TRUE.
2269 :
2270 : ELSE
2271 :
2272 : atom_pair_changed = .FALSE.
2273 :
2274 : END IF
2275 :
2276 16254 : IF (atom_pair_changed .OR. iset_old /= iset .OR. jset_old /= jset) THEN
2277 :
2278 1296 : ncoa = npgfa(iset)*ncoset(la_max(iset))
2279 1296 : sgfa = first_sgfa(1, iset)
2280 1296 : ncob = npgfb(jset)*ncoset(lb_max(jset))
2281 1296 : sgfb = first_sgfb(1, jset)
2282 :
2283 1296 : IF (iatom <= jatom) THEN
2284 : CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
2285 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
2286 : p_block(sgfa, sgfb), SIZE(p_block, 1), &
2287 864 : 0.0_dp, work(1, 1), maxco)
2288 : CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
2289 : 1.0_dp, work(1, 1), maxco, &
2290 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
2291 864 : 0.0_dp, pab(1, 1), maxco)
2292 : ELSE
2293 : CALL dgemm("N", "N", ncob, nsgfa(iset), nsgfb(jset), &
2294 : 1.0_dp, sphi_b(1, sgfb), SIZE(sphi_b, 1), &
2295 : p_block(sgfb, sgfa), SIZE(p_block, 1), &
2296 432 : 0.0_dp, work(1, 1), maxco)
2297 : CALL dgemm("N", "T", ncob, ncoa, nsgfa(iset), &
2298 : 1.0_dp, work(1, 1), maxco, &
2299 : sphi_a(1, sgfa), SIZE(sphi_a, 1), &
2300 432 : 0.0_dp, pab(1, 1), maxco)
2301 : END IF
2302 :
2303 : iset_old = iset
2304 : jset_old = jset
2305 :
2306 : END IF
2307 :
2308 65016 : rab(:) = tasks(itask)%rab
2309 65016 : rb(:) = ra(:) + rab(:)
2310 16254 : zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
2311 :
2312 16254 : f = zetb(jpgf, jset)/zetp
2313 65016 : rp(:) = ra(:) + f*rab(:)
2314 65016 : prefactor = EXP(-zeta(ipgf, iset)*f*DOT_PRODUCT(rab, rab))
2315 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
2316 : lb_min=lb_min(jset), lb_max=lb_max(jset), &
2317 : ra=ra, rb=rb, rp=rp, &
2318 : zetp=zetp, eps=eps_rho_rspace, &
2319 16254 : prefactor=prefactor, cutoff=1.0_dp)
2320 :
2321 16254 : na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
2322 16254 : na2 = ipgf*ncoset(la_max(iset))
2323 16254 : nb1 = (jpgf - 1)*ncoset(lb_max(jset)) + 1
2324 16254 : nb2 = jpgf*ncoset(lb_max(jset))
2325 :
2326 : ! takes the density matrix symmetry in account, i.e. off-diagonal blocks need to be mapped 'twice'
2327 16254 : IF (iatom == jatom .AND. img == 1) THEN
2328 8100 : scale = 1.0_dp
2329 : ELSE
2330 8154 : scale = 2.0_dp
2331 : END IF
2332 :
2333 : ! check whether we need to use fawzi's generalised collocation scheme
2334 16254 : IF (rs_rho(igrid_level)%desc%distributed) THEN
2335 : !tasks(4,:) is 0 for replicated, 1 for distributed 2 for exceptional distributed tasks
2336 0 : IF (tasks(itask)%dist_type == 2) THEN
2337 0 : use_subpatch = .TRUE.
2338 : ELSE
2339 0 : use_subpatch = .FALSE.
2340 : END IF
2341 : ELSE
2342 16254 : use_subpatch = .FALSE.
2343 : END IF
2344 :
2345 21672 : SELECT CASE (beta)
2346 : CASE (1)
2347 5418 : dabqadb_func = GRID_FUNC_DAB_X
2348 : CASE (2)
2349 5418 : dabqadb_func = GRID_FUNC_DAB_Y
2350 : CASE (3)
2351 5418 : dabqadb_func = GRID_FUNC_DAB_Z
2352 : CASE DEFAULT
2353 16254 : CPABORT("invalid beta")
2354 : END SELECT
2355 :
2356 16506 : IF (iatom <= jatom) THEN
2357 10854 : IF (iatom == lambda) &
2358 : CALL collocate_pgf_product( &
2359 : la_max(iset), zeta(ipgf, iset), la_min(iset), &
2360 : lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
2361 : ra, rab, scale, pab, na1 - 1, nb1 - 1, &
2362 : rsgrid=rs_rho(igrid_level), &
2363 : ga_gb_function=dabqadb_func, radius=radius, &
2364 : use_subpatch=use_subpatch, &
2365 3618 : subpatch_pattern=tasks(itask)%subpatch_pattern)
2366 10854 : IF (jatom == lambda) &
2367 : CALL collocate_pgf_product( &
2368 : la_max(iset), zeta(ipgf, iset), la_min(iset), &
2369 : lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
2370 : ra, rab, scale, pab, na1 - 1, nb1 - 1, &
2371 : rsgrid=rs_rho(igrid_level), &
2372 : ga_gb_function=dabqadb_func + 3, radius=radius, &
2373 : use_subpatch=use_subpatch, &
2374 3618 : subpatch_pattern=tasks(itask)%subpatch_pattern)
2375 : ELSE
2376 21600 : rab_inv = -rab
2377 5400 : IF (jatom == lambda) &
2378 : CALL collocate_pgf_product( &
2379 : lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
2380 : la_max(iset), zeta(ipgf, iset), la_min(iset), &
2381 : rb, rab_inv, scale, pab, nb1 - 1, na1 - 1, &
2382 : rs_rho(igrid_level), &
2383 : ga_gb_function=dabqadb_func, radius=radius, &
2384 : use_subpatch=use_subpatch, &
2385 1800 : subpatch_pattern=tasks(itask)%subpatch_pattern)
2386 5400 : IF (iatom == lambda) &
2387 : CALL collocate_pgf_product( &
2388 : lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
2389 : la_max(iset), zeta(ipgf, iset), la_min(iset), &
2390 : rb, rab_inv, scale, pab, nb1 - 1, na1 - 1, &
2391 : rs_rho(igrid_level), &
2392 : ga_gb_function=dabqadb_func + 3, radius=radius, &
2393 : use_subpatch=use_subpatch, &
2394 1800 : subpatch_pattern=tasks(itask)%subpatch_pattern)
2395 : END IF
2396 :
2397 : END DO loop_tasks
2398 :
2399 252 : CALL density_rs2pw(pw_env, rs_rho, drho, drho_gspace)
2400 :
2401 : ! *** Release work storage ***
2402 252 : IF (distributed_rs_grids) THEN
2403 0 : CALL dbcsr_deallocate_matrix_set(deltap)
2404 : ELSE
2405 504 : DO img = 1, nimages
2406 504 : NULLIFY (deltap(img)%matrix)
2407 : END DO
2408 252 : DEALLOCATE (deltap)
2409 : END IF
2410 :
2411 252 : DEALLOCATE (pabt, workt)
2412 :
2413 252 : CALL timestop(handle)
2414 :
2415 504 : END SUBROUTINE calculate_drho_elec_dR
2416 :
2417 : ! **************************************************************************************************
2418 : !> \brief maps a single gaussian on the grid
2419 : !> \param rho ...
2420 : !> \param rho_gspace ...
2421 : !> \param atomic_kind_set ...
2422 : !> \param qs_kind_set ...
2423 : !> \param cell ...
2424 : !> \param dft_control ...
2425 : !> \param particle_set ...
2426 : !> \param pw_env ...
2427 : !> \param required_function ...
2428 : !> \param basis_type ...
2429 : !> \par History
2430 : !> 08.2022 created from calculate_wavefunction
2431 : !> \note
2432 : !> modified calculate_wave function assuming that the collocation of only a single Gaussian is required.
2433 : !> chooses a basis function (in contrast to calculate_rho_core or calculate_rho_single_gaussian)
2434 : ! **************************************************************************************************
2435 28573 : SUBROUTINE collocate_single_gaussian(rho, rho_gspace, &
2436 : atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
2437 : pw_env, required_function, basis_type)
2438 :
2439 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho
2440 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_gspace
2441 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2442 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2443 : TYPE(cell_type), POINTER :: cell
2444 : TYPE(dft_control_type), POINTER :: dft_control
2445 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2446 : TYPE(pw_env_type), POINTER :: pw_env
2447 : INTEGER, INTENT(IN) :: required_function
2448 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: basis_type
2449 :
2450 : CHARACTER(LEN=*), PARAMETER :: routineN = 'collocate_single_gaussian'
2451 :
2452 : CHARACTER(LEN=default_string_length) :: my_basis_type
2453 : INTEGER :: group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, maxco, maxsgf_set, &
2454 : my_index, my_pos, na1, na2, natom, ncoa, nseta, offset, sgfa
2455 28573 : INTEGER, ALLOCATABLE, DIMENSION(:) :: where_is_the_point
2456 28573 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
2457 28573 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
2458 : LOGICAL :: found
2459 : REAL(KIND=dp) :: dab, eps_rho_rspace, radius, scale
2460 : REAL(KIND=dp), DIMENSION(3) :: ra
2461 28573 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab, sphi_a, zeta
2462 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
2463 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
2464 : TYPE(mp_comm_type) :: group
2465 28573 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
2466 28573 : TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_gspace
2467 28573 : TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_rspace
2468 28573 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
2469 :
2470 28573 : IF (PRESENT(basis_type)) THEN
2471 28573 : my_basis_type = basis_type
2472 : ELSE
2473 0 : my_basis_type = "ORB"
2474 : END IF
2475 :
2476 28573 : CALL timeset(routineN, handle)
2477 :
2478 28573 : NULLIFY (orb_basis_set, pab, la_max, la_min, npgfa, nsgfa, sphi_a, &
2479 28573 : zeta, first_sgfa, rs_rho, pw_pools)
2480 :
2481 : ! *** set up the pw multi-grids
2482 28573 : CPASSERT(ASSOCIATED(pw_env))
2483 : CALL pw_env_get(pw_env, rs_grids=rs_rho, pw_pools=pw_pools, &
2484 28573 : gridlevel_info=gridlevel_info)
2485 :
2486 28573 : CALL pw_pools_create_pws(pw_pools, mgrid_gspace)
2487 28573 : CALL pw_pools_create_pws(pw_pools, mgrid_rspace)
2488 :
2489 : ! *** set up rs multi-grids
2490 142865 : DO igrid_level = 1, gridlevel_info%ngrid_levels
2491 142865 : CALL rs_grid_zero(rs_rho(igrid_level))
2492 : END DO
2493 :
2494 28573 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
2495 : ! *** Allocate work storage ***
2496 28573 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
2497 : CALL get_qs_kind_set(qs_kind_set, &
2498 : maxco=maxco, &
2499 : maxsgf_set=maxsgf_set, &
2500 28573 : basis_type=my_basis_type)
2501 :
2502 85719 : ALLOCATE (pab(maxco, 1))
2503 :
2504 28573 : offset = 0
2505 28573 : group = mgrid_rspace(1)%pw_grid%para%group
2506 28573 : my_pos = mgrid_rspace(1)%pw_grid%para%group%mepos
2507 28573 : group_size = mgrid_rspace(1)%pw_grid%para%group%num_pe
2508 85719 : ALLOCATE (where_is_the_point(0:group_size - 1))
2509 :
2510 117551 : DO iatom = 1, natom
2511 88978 : ikind = particle_set(iatom)%atomic_kind%kind_number
2512 88978 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=my_basis_type)
2513 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
2514 : first_sgf=first_sgfa, &
2515 : lmax=la_max, &
2516 : lmin=la_min, &
2517 : npgf=npgfa, &
2518 : nset=nseta, &
2519 : nsgf_set=nsgfa, &
2520 : sphi=sphi_a, &
2521 88978 : zet=zeta)
2522 88978 : ra(:) = pbc(particle_set(iatom)%r, cell)
2523 88978 : dab = 0.0_dp
2524 :
2525 1047571 : DO iset = 1, nseta
2526 :
2527 841042 : ncoa = npgfa(iset)*ncoset(la_max(iset))
2528 841042 : sgfa = first_sgfa(1, iset)
2529 :
2530 841042 : found = .FALSE.
2531 841042 : my_index = 0
2532 3176013 : DO i = 1, nsgfa(iset)
2533 3176013 : IF (offset + i == required_function) THEN
2534 : my_index = i
2535 : found = .TRUE.
2536 : EXIT
2537 : END IF
2538 : END DO
2539 :
2540 841042 : IF (found) THEN
2541 :
2542 523769 : pab(1:ncoa, 1) = sphi_a(1:ncoa, sgfa + my_index - 1)
2543 :
2544 58202 : DO ipgf = 1, npgfa(iset)
2545 :
2546 29629 : na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
2547 29629 : na2 = ipgf*ncoset(la_max(iset))
2548 :
2549 29629 : scale = 1.0_dp
2550 29629 : igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
2551 :
2552 58202 : IF (map_gaussian_here(rs_rho(igrid_level), cell%h_inv, ra, offset, group_size, my_pos)) THEN
2553 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
2554 : lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
2555 : zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
2556 27872 : prefactor=1.0_dp, cutoff=1.0_dp)
2557 :
2558 : CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), la_min(iset), &
2559 : 0, 0.0_dp, 0, &
2560 : ra, [0.0_dp, 0.0_dp, 0.0_dp], &
2561 : scale, pab, na1 - 1, 0, rs_rho(igrid_level), &
2562 27872 : radius=radius, ga_gb_function=GRID_FUNC_AB)
2563 : END IF
2564 :
2565 : END DO
2566 :
2567 : END IF
2568 :
2569 930020 : offset = offset + nsgfa(iset)
2570 :
2571 : END DO
2572 :
2573 : END DO
2574 :
2575 142865 : DO igrid_level = 1, gridlevel_info%ngrid_levels
2576 : CALL transfer_rs2pw(rs_rho(igrid_level), &
2577 142865 : mgrid_rspace(igrid_level))
2578 : END DO
2579 :
2580 28573 : CALL pw_zero(rho_gspace)
2581 142865 : DO igrid_level = 1, gridlevel_info%ngrid_levels
2582 : CALL pw_transfer(mgrid_rspace(igrid_level), &
2583 114292 : mgrid_gspace(igrid_level))
2584 142865 : CALL pw_axpy(mgrid_gspace(igrid_level), rho_gspace)
2585 : END DO
2586 :
2587 28573 : CALL pw_transfer(rho_gspace, rho)
2588 :
2589 : ! Release work storage
2590 28573 : DEALLOCATE (pab)
2591 :
2592 : ! give back the pw multi-grids
2593 28573 : CALL pw_pools_give_back_pws(pw_pools, mgrid_gspace)
2594 28573 : CALL pw_pools_give_back_pws(pw_pools, mgrid_rspace)
2595 :
2596 28573 : CALL timestop(handle)
2597 :
2598 171438 : END SUBROUTINE collocate_single_gaussian
2599 :
2600 : ! **************************************************************************************************
2601 : !> \brief maps a given wavefunction on the grid
2602 : !> \param mo_vectors ...
2603 : !> \param ivector ...
2604 : !> \param rho ...
2605 : !> \param rho_gspace ...
2606 : !> \param atomic_kind_set ...
2607 : !> \param qs_kind_set ...
2608 : !> \param cell ...
2609 : !> \param dft_control ...
2610 : !> \param particle_set ...
2611 : !> \param pw_env ...
2612 : !> \param basis_type ...
2613 : !> \par History
2614 : !> 08.2002 created [Joost VandeVondele]
2615 : !> 03.2006 made independent of qs_env [Joost VandeVondele]
2616 : !> 08.2024 call collocate_function [JGH]
2617 : ! **************************************************************************************************
2618 1350 : SUBROUTINE calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, &
2619 : atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
2620 : pw_env, basis_type)
2621 : TYPE(cp_fm_type), INTENT(IN) :: mo_vectors
2622 : INTEGER, INTENT(IN) :: ivector
2623 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho
2624 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_gspace
2625 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2626 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2627 : TYPE(cell_type), POINTER :: cell
2628 : TYPE(dft_control_type), POINTER :: dft_control
2629 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2630 : TYPE(pw_env_type), POINTER :: pw_env
2631 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: basis_type
2632 :
2633 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_wavefunction'
2634 :
2635 : INTEGER :: handle, i, nao
2636 : LOGICAL :: local
2637 : REAL(KIND=dp) :: eps_rho_rspace
2638 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvector
2639 :
2640 1350 : CALL timeset(routineN, handle)
2641 :
2642 1350 : CALL cp_fm_get_info(matrix=mo_vectors, nrow_global=nao)
2643 4050 : ALLOCATE (eigenvector(nao))
2644 24846 : DO i = 1, nao
2645 24846 : CALL cp_fm_get_element(mo_vectors, i, ivector, eigenvector(i), local)
2646 : END DO
2647 :
2648 1350 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
2649 :
2650 : CALL collocate_function(eigenvector, rho, rho_gspace, &
2651 : atomic_kind_set, qs_kind_set, cell, particle_set, pw_env, &
2652 2700 : eps_rho_rspace, basis_type)
2653 :
2654 1350 : DEALLOCATE (eigenvector)
2655 :
2656 1350 : CALL timestop(handle)
2657 :
2658 1350 : END SUBROUTINE calculate_wavefunction
2659 :
2660 : ! **************************************************************************************************
2661 : !> \brief maps a given function on the grid
2662 : !> \param vector ...
2663 : !> \param rho ...
2664 : !> \param rho_gspace ...
2665 : !> \param atomic_kind_set ...
2666 : !> \param qs_kind_set ...
2667 : !> \param cell ...
2668 : !> \param particle_set ...
2669 : !> \param pw_env ...
2670 : !> \param eps_rho_rspace ...
2671 : !> \param basis_type ...
2672 : !> \par History
2673 : !> 08.2002 created [Joost VandeVondele]
2674 : !> 03.2006 made independent of qs_env [Joost VandeVondele]
2675 : !> 08.2024 specialized version from calculate_wavefunction [JGH]
2676 : !> \notes
2677 : !> modified calculate_rho_elec, should write the wavefunction represented by vector
2678 : !> it's presumably dominated by the FFT and the rs->pw and back routines
2679 : ! **************************************************************************************************
2680 39616 : SUBROUTINE collocate_function(vector, rho, rho_gspace, &
2681 : atomic_kind_set, qs_kind_set, cell, particle_set, pw_env, &
2682 : eps_rho_rspace, basis_type)
2683 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vector
2684 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho
2685 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_gspace
2686 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2687 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2688 : TYPE(cell_type), POINTER :: cell
2689 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2690 : TYPE(pw_env_type), POINTER :: pw_env
2691 : REAL(KIND=dp), INTENT(IN) :: eps_rho_rspace
2692 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: basis_type
2693 :
2694 : CHARACTER(LEN=*), PARAMETER :: routineN = 'collocate_function'
2695 :
2696 : CHARACTER(LEN=default_string_length) :: my_basis_type
2697 : INTEGER :: group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, maxco, maxsgf_set, &
2698 : my_pos, na1, na2, natom, ncoa, nseta, offset, sgfa
2699 19808 : INTEGER, ALLOCATABLE, DIMENSION(:) :: where_is_the_point
2700 19808 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
2701 19808 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
2702 : REAL(KIND=dp) :: dab, radius, scale
2703 : REAL(KIND=dp), DIMENSION(3) :: ra
2704 19808 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab, sphi_a, work, zeta
2705 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
2706 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
2707 : TYPE(mp_comm_type) :: group
2708 19808 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
2709 19808 : TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_gspace
2710 19808 : TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_rspace
2711 19808 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
2712 :
2713 19808 : CALL timeset(routineN, handle)
2714 :
2715 19808 : IF (PRESENT(basis_type)) THEN
2716 18154 : my_basis_type = basis_type
2717 : ELSE
2718 1654 : my_basis_type = "ORB"
2719 : END IF
2720 :
2721 19808 : NULLIFY (orb_basis_set, pab, work, la_max, la_min, &
2722 19808 : npgfa, nsgfa, sphi_a, zeta, first_sgfa, rs_rho, pw_pools)
2723 :
2724 : ! *** set up the pw multi-grids
2725 19808 : CPASSERT(ASSOCIATED(pw_env))
2726 : CALL pw_env_get(pw_env, rs_grids=rs_rho, pw_pools=pw_pools, &
2727 19808 : gridlevel_info=gridlevel_info)
2728 :
2729 19808 : CALL pw_pools_create_pws(pw_pools, mgrid_gspace)
2730 19808 : CALL pw_pools_create_pws(pw_pools, mgrid_rspace)
2731 :
2732 : ! *** set up rs multi-grids
2733 98824 : DO igrid_level = 1, gridlevel_info%ngrid_levels
2734 98824 : CALL rs_grid_zero(rs_rho(igrid_level))
2735 : END DO
2736 :
2737 : ! *** Allocate work storage ***
2738 19808 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
2739 : CALL get_qs_kind_set(qs_kind_set, &
2740 : maxco=maxco, &
2741 : maxsgf_set=maxsgf_set, &
2742 19808 : basis_type=my_basis_type)
2743 :
2744 59424 : ALLOCATE (pab(maxco, 1))
2745 39616 : ALLOCATE (work(maxco, 1))
2746 :
2747 19808 : offset = 0
2748 19808 : group = mgrid_rspace(1)%pw_grid%para%group
2749 19808 : my_pos = mgrid_rspace(1)%pw_grid%para%group%mepos
2750 19808 : group_size = mgrid_rspace(1)%pw_grid%para%group%num_pe
2751 59424 : ALLOCATE (where_is_the_point(0:group_size - 1))
2752 :
2753 82054 : DO iatom = 1, natom
2754 62246 : ikind = particle_set(iatom)%atomic_kind%kind_number
2755 62246 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=my_basis_type)
2756 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
2757 : first_sgf=first_sgfa, &
2758 : lmax=la_max, &
2759 : lmin=la_min, &
2760 : npgf=npgfa, &
2761 : nset=nseta, &
2762 : nsgf_set=nsgfa, &
2763 : sphi=sphi_a, &
2764 62246 : zet=zeta)
2765 62246 : ra(:) = pbc(particle_set(iatom)%r, cell)
2766 62246 : dab = 0.0_dp
2767 :
2768 689557 : DO iset = 1, nseta
2769 :
2770 545257 : ncoa = npgfa(iset)*ncoset(la_max(iset))
2771 545257 : sgfa = first_sgfa(1, iset)
2772 :
2773 2107976 : DO i = 1, nsgfa(iset)
2774 2107976 : work(i, 1) = vector(offset + i)
2775 : END DO
2776 :
2777 : CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), &
2778 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
2779 : work(1, 1), SIZE(work, 1), &
2780 545257 : 0.0_dp, pab(1, 1), SIZE(pab, 1))
2781 :
2782 1115751 : DO ipgf = 1, npgfa(iset)
2783 :
2784 570494 : na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
2785 570494 : na2 = ipgf*ncoset(la_max(iset))
2786 :
2787 570494 : scale = 1.0_dp
2788 570494 : igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
2789 :
2790 1115751 : IF (map_gaussian_here(rs_rho(igrid_level), cell%h_inv, ra, offset, group_size, my_pos)) THEN
2791 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
2792 : lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
2793 : zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
2794 521622 : prefactor=1.0_dp, cutoff=1.0_dp)
2795 :
2796 : CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), la_min(iset), &
2797 : 0, 0.0_dp, 0, &
2798 : ra, [0.0_dp, 0.0_dp, 0.0_dp], &
2799 : scale, pab, na1 - 1, 0, rs_rho(igrid_level), &
2800 521622 : radius=radius, ga_gb_function=GRID_FUNC_AB)
2801 : END IF
2802 :
2803 : END DO
2804 :
2805 607503 : offset = offset + nsgfa(iset)
2806 :
2807 : END DO
2808 :
2809 : END DO
2810 :
2811 98824 : DO igrid_level = 1, gridlevel_info%ngrid_levels
2812 : CALL transfer_rs2pw(rs_rho(igrid_level), &
2813 98824 : mgrid_rspace(igrid_level))
2814 : END DO
2815 :
2816 19808 : CALL pw_zero(rho_gspace)
2817 98824 : DO igrid_level = 1, gridlevel_info%ngrid_levels
2818 : CALL pw_transfer(mgrid_rspace(igrid_level), &
2819 79016 : mgrid_gspace(igrid_level))
2820 98824 : CALL pw_axpy(mgrid_gspace(igrid_level), rho_gspace)
2821 : END DO
2822 :
2823 19808 : CALL pw_transfer(rho_gspace, rho)
2824 :
2825 : ! Release work storage
2826 19808 : DEALLOCATE (pab)
2827 19808 : DEALLOCATE (work)
2828 :
2829 : ! give back the pw multi-grids
2830 19808 : CALL pw_pools_give_back_pws(pw_pools, mgrid_gspace)
2831 19808 : CALL pw_pools_give_back_pws(pw_pools, mgrid_rspace)
2832 :
2833 19808 : CALL timestop(handle)
2834 :
2835 99040 : END SUBROUTINE collocate_function
2836 :
2837 : END MODULE qs_collocate_density
|
