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