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