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