(git:d18deda)
Loading...
Searching...
No Matches
qs_collocate_density.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculate the plane wave density by collocating the primitive Gaussian
10!> functions (pgf).
11!> \par History
12!> - rewrote collocate for increased accuracy and speed
13!> - introduced the PGI hack for increased speed with that compiler
14!> (22.02.02)
15!> - Added Multiple Grid feature
16!> - new way to get over the grid (01.03.02)
17!> - removed timing calls since they were getting expensive
18!> - Updated with the new QS data structures (09.04.02,MK)
19!> - introduction of the real space grid type ( prelim. version JVdV 05.02)
20!> - parallel FFT (JGH 22.05.02)
21!> - multigrid arrays independent from density (JGH 30.08.02)
22!> - old density stored in g space (JGH 30.08.02)
23!> - distributed real space code (JGH 17.07.03)
24!> - refactoring and new loop ordering (JGH 23.11.03)
25!> - OpenMP parallelization (JGH 03.12.03)
26!> - Modified to compute tau (Joost 12.03)
27!> - removed the incremental density rebuild (Joost 01.04)
28!> - introduced realspace multigridding (Joost 02.04)
29!> - introduced map_consistent (Joost 02.04)
30!> - Addition of the subroutine calculate_atomic_charge_density (TdK, 08.05)
31!> - rewrite of the collocate/integrate kernels (Joost VandeVondele, 03.07)
32!> - Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
33!> \author Matthias Krack (03.04.2001)
34!> 1) Joost VandeVondele (01.2002)
35!> Thomas D. Kuehne (04.08.2005)
36!> Ole Schuett (2020)
37! **************************************************************************************************
38MODULE qs_collocate_density
39 USE admm_types, ONLY: get_admm_env
40 USE ao_util, ONLY: exp_radius_very_extended
41 USE atomic_kind_types, ONLY: atomic_kind_type, &
42 get_atomic_kind, &
43 get_atomic_kind_set
44 USE basis_set_types, ONLY: get_gto_basis_set, &
45 gto_basis_set_type
46 USE cell_types, ONLY: cell_type, &
47 pbc
48 USE cp_control_types, ONLY: dft_control_type
49 USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
50 USE cp_fm_types, ONLY: cp_fm_get_element, &
51 cp_fm_get_info, &
52 cp_fm_type
53 USE cp_dbcsr_api, ONLY: dbcsr_copy, &
54 dbcsr_get_block_p, &
55 dbcsr_p_type, &
56 dbcsr_type
57 USE external_potential_types, ONLY: get_potential, &
58 gth_potential_type
59 USE gaussian_gridlevels, ONLY: gaussian_gridlevel, &
60 gridlevel_info_type
61 USE grid_api, ONLY: &
62 grid_func_ab, grid_func_core_x, grid_func_core_y, grid_func_core_z, grid_func_dab_x, &
63 grid_func_dab_y, grid_func_dab_z, grid_func_dabpadb_x, grid_func_dabpadb_y, &
64 grid_func_dabpadb_z, grid_func_dadb, grid_func_dx, grid_func_dxdx, grid_func_dxdy, &
65 grid_func_dy, grid_func_dydy, grid_func_dydz, grid_func_dz, grid_func_dzdx, &
66 grid_func_dzdz, collocate_pgf_product, grid_collocate_task_list
67 USE input_constants, ONLY: &
68 orb_dx2, orb_dxy, orb_dy2, orb_dyz, orb_dz2, orb_dzx, orb_px, orb_py, orb_pz, orb_s
69 USE kinds, ONLY: default_string_length, &
70 dp
71 USE lri_environment_types, ONLY: lri_kind_type
72 USE memory_utilities, ONLY: reallocate
73 USE message_passing, ONLY: mp_comm_type
74 USE orbital_pointers, ONLY: coset, &
75 ncoset
76 USE particle_types, ONLY: particle_type
77 USE pw_env_types, ONLY: pw_env_get, &
78 pw_env_type
79 USE pw_methods, ONLY: pw_axpy, &
80 pw_integrate_function, &
81 pw_transfer, &
82 pw_zero
83 USE pw_pool_types, ONLY: pw_pool_p_type, &
84 pw_pool_type, &
85 pw_pools_create_pws, &
86 pw_pools_give_back_pws
87 USE pw_types, ONLY: pw_r3d_rs_type, &
88 pw_c1d_gs_type, &
89 pw_r3d_rs_type
90 USE qs_environment_types, ONLY: get_qs_env, &
91 qs_environment_type
92 USE qs_kind_types, ONLY: get_qs_kind, &
93 get_qs_kind_set, &
94 qs_kind_type
95 USE qs_ks_types, ONLY: get_ks_env, &
96 qs_ks_env_type
97 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
98 USE realspace_grid_types, ONLY: map_gaussian_here, &
99 realspace_grid_desc_p_type, &
100 realspace_grid_type, &
101 rs_grid_zero, &
102 transfer_rs2pw
103 USE rs_pw_interface, ONLY: density_rs2pw
104 USE task_list_methods, ONLY: rs_copy_to_buffer, &
105 rs_distribute_matrix, &
106 rs_scatter_matrices
107 USE task_list_types, ONLY: atom_pair_type, &
108 task_list_type, &
109 task_type
110
111!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
112
113#include "./base/base_uses.f90"
114
115 IMPLICIT NONE
116
117 PRIVATE
118
119 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_collocate_density'
120! *** Public subroutines ***
121
122 PUBLIC :: calculate_ppl_grid, &
123 calculate_rho_core, &
124 calculate_lri_rho_elec, &
125 calculate_rho_single_gaussian, &
126 calculate_rho_metal, &
127 calculate_rho_resp_single, &
128 calculate_rho_resp_all, &
129 calculate_rho_elec, &
130 calculate_drho_elec, &
131 calculate_wavefunction, &
132 collocate_function, &
133 calculate_rho_nlcc, &
134 calculate_drho_elec_dr, &
135 calculate_drho_core, &
136 collocate_single_gaussian
137
138 INTERFACE calculate_rho_core
139 MODULE PROCEDURE calculate_rho_core_r3d_rs
140 MODULE PROCEDURE calculate_rho_core_c1d_gs
141 END INTERFACE
142
143 INTERFACE calculate_rho_resp_all
144 MODULE PROCEDURE calculate_rho_resp_all_r3d_rs, calculate_rho_resp_all_c1d_gs
145 END INTERFACE
146
147CONTAINS
148
149! **************************************************************************************************
150!> \brief computes the density of the non-linear core correction on the grid
151!> \param rho_nlcc ...
152!> \param qs_env ...
153! **************************************************************************************************
154 SUBROUTINE calculate_rho_nlcc(rho_nlcc, qs_env)
155
156 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho_nlcc
157 TYPE(qs_environment_type), POINTER :: qs_env
158
159 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_rho_nlcc'
160
161 INTEGER :: atom_a, handle, iatom, iexp_nlcc, ikind, &
162 ithread, j, n, natom, nc, nexp_nlcc, &
163 ni, npme, nthread, subpatch_pattern
164 INTEGER, DIMENSION(:), POINTER :: atom_list, cores, nct_nlcc
165 LOGICAL :: nlcc
166 REAL(kind=dp) :: alpha, eps_rho_rspace, radius
167 REAL(kind=dp), DIMENSION(3) :: ra
168 REAL(kind=dp), DIMENSION(:), POINTER :: alpha_nlcc
169 REAL(kind=dp), DIMENSION(:, :), POINTER :: cval_nlcc, pab
170 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
171 TYPE(cell_type), POINTER :: cell
172 TYPE(dft_control_type), POINTER :: dft_control
173 TYPE(gth_potential_type), POINTER :: gth_potential
174 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
175 TYPE(pw_env_type), POINTER :: pw_env
176 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
177 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
178 TYPE(realspace_grid_type), POINTER :: rs_rho
179
180 CALL timeset(routinen, handle)
181
182 NULLIFY (cell, dft_control, pab, particle_set, atomic_kind_set, &
183 qs_kind_set, atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
184
185 CALL get_qs_env(qs_env=qs_env, &
186 atomic_kind_set=atomic_kind_set, &
187 qs_kind_set=qs_kind_set, &
188 cell=cell, &
189 dft_control=dft_control, &
190 particle_set=particle_set, &
191 pw_env=pw_env)
192 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
193 auxbas_pw_pool=auxbas_pw_pool)
194 ! be careful in parallel nsmax is chosen with multigrid in mind!
195 CALL rs_grid_zero(rs_rho)
196
197 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
198
199 DO ikind = 1, SIZE(atomic_kind_set)
200 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
201 CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
202
203 IF (.NOT. ASSOCIATED(gth_potential)) cycle
204 CALL get_potential(potential=gth_potential, nlcc_present=nlcc, nexp_nlcc=nexp_nlcc, &
205 alpha_nlcc=alpha_nlcc, nct_nlcc=nct_nlcc, cval_nlcc=cval_nlcc)
206
207 IF (.NOT. nlcc) cycle
208
209 DO iexp_nlcc = 1, nexp_nlcc
210
211 alpha = alpha_nlcc(iexp_nlcc)
212 nc = nct_nlcc(iexp_nlcc)
213
214 ni = ncoset(2*nc - 2)
215 ALLOCATE (pab(ni, 1))
216 pab = 0._dp
217
218 nthread = 1
219 ithread = 0
220
221 CALL reallocate(cores, 1, natom)
222 npme = 0
223 cores = 0
224
225 ! prepare core function
226 DO j = 1, nc
227 SELECT CASE (j)
228 CASE (1)
229 pab(1, 1) = cval_nlcc(1, iexp_nlcc)
230 CASE (2)
231 n = coset(2, 0, 0)
232 pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
233 n = coset(0, 2, 0)
234 pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
235 n = coset(0, 0, 2)
236 pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
237 CASE (3)
238 n = coset(4, 0, 0)
239 pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
240 n = coset(0, 4, 0)
241 pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
242 n = coset(0, 0, 4)
243 pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
244 n = coset(2, 2, 0)
245 pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
246 n = coset(2, 0, 2)
247 pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
248 n = coset(0, 2, 2)
249 pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
250 CASE (4)
251 n = coset(6, 0, 0)
252 pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
253 n = coset(0, 6, 0)
254 pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
255 n = coset(0, 0, 6)
256 pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
257 n = coset(4, 2, 0)
258 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
259 n = coset(4, 0, 2)
260 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
261 n = coset(2, 4, 0)
262 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
263 n = coset(2, 0, 4)
264 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
265 n = coset(0, 4, 2)
266 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
267 n = coset(0, 2, 4)
268 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
269 n = coset(2, 2, 2)
270 pab(n, 1) = 6._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
271 CASE DEFAULT
272 cpabort("")
273 END SELECT
274 END DO
275 IF (dft_control%nspins == 2) pab = pab*0.5_dp
276
277 DO iatom = 1, natom
278 atom_a = atom_list(iatom)
279 ra(:) = pbc(particle_set(atom_a)%r, cell)
280 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
281 ! replicated realspace grid, split the atoms up between procs
282 IF (modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
283 npme = npme + 1
284 cores(npme) = iatom
285 END IF
286 ELSE
287 npme = npme + 1
288 cores(npme) = iatom
289 END IF
290 END DO
291
292 DO j = 1, npme
293
294 iatom = cores(j)
295 atom_a = atom_list(iatom)
296 ra(:) = pbc(particle_set(atom_a)%r, cell)
297 subpatch_pattern = 0
298 ni = 2*nc - 2
299 radius = exp_radius_very_extended(la_min=0, la_max=ni, lb_min=0, lb_max=0, &
300 ra=ra, rb=ra, rp=ra, &
301 zetp=1/(2*alpha**2), eps=eps_rho_rspace, &
302 pab=pab, o1=0, o2=0, & ! without map_consistent
303 prefactor=1.0_dp, cutoff=0.0_dp)
304
305 CALL collocate_pgf_product(ni, 1/(2*alpha**2), 0, 0, 0.0_dp, 0, ra, &
306 (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, rs_rho, &
307 ga_gb_function=grid_func_ab, radius=radius, &
308 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
309
310 END DO
311
312 DEALLOCATE (pab)
313
314 END DO
315
316 END DO
317
318 IF (ASSOCIATED(cores)) THEN
319 DEALLOCATE (cores)
320 END IF
321
322 CALL transfer_rs2pw(rs_rho, rho_nlcc)
323
324 CALL timestop(handle)
325
326 END SUBROUTINE calculate_rho_nlcc
327
328! **************************************************************************************************
329!> \brief computes the local pseudopotential (without erf term) on the grid
330!> \param vppl ...
331!> \param qs_env ...
332! **************************************************************************************************
333 SUBROUTINE calculate_ppl_grid(vppl, qs_env)
334
335 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vppl
336 TYPE(qs_environment_type), POINTER :: qs_env
337
338 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_ppl_grid'
339
340 INTEGER :: atom_a, handle, iatom, ikind, ithread, &
341 j, lppl, n, natom, ni, npme, nthread, &
342 subpatch_pattern
343 INTEGER, DIMENSION(:), POINTER :: atom_list, cores
344 REAL(kind=dp) :: alpha, eps_rho_rspace, radius
345 REAL(kind=dp), DIMENSION(3) :: ra
346 REAL(kind=dp), DIMENSION(:), POINTER :: cexp_ppl
347 REAL(kind=dp), DIMENSION(:, :), POINTER :: pab
348 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
349 TYPE(cell_type), POINTER :: cell
350 TYPE(dft_control_type), POINTER :: dft_control
351 TYPE(gth_potential_type), POINTER :: gth_potential
352 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
353 TYPE(pw_env_type), POINTER :: pw_env
354 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
355 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
356 TYPE(realspace_grid_type), POINTER :: rs_rho
357
358 CALL timeset(routinen, handle)
359
360 NULLIFY (cell, dft_control, pab, atomic_kind_set, qs_kind_set, particle_set, &
361 atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
362
363 CALL get_qs_env(qs_env=qs_env, &
364 atomic_kind_set=atomic_kind_set, &
365 qs_kind_set=qs_kind_set, &
366 cell=cell, &
367 dft_control=dft_control, &
368 particle_set=particle_set, &
369 pw_env=pw_env)
370 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
371 auxbas_pw_pool=auxbas_pw_pool)
372 ! be careful in parallel nsmax is chosen with multigrid in mind!
373 CALL rs_grid_zero(rs_rho)
374
375 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
376
377 DO ikind = 1, SIZE(atomic_kind_set)
378 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
379 CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
380
381 IF (.NOT. ASSOCIATED(gth_potential)) cycle
382 CALL get_potential(potential=gth_potential, alpha_ppl=alpha, nexp_ppl=lppl, cexp_ppl=cexp_ppl)
383
384 IF (lppl <= 0) cycle
385
386 ni = ncoset(2*lppl - 2)
387 ALLOCATE (pab(ni, 1))
388 pab = 0._dp
389
390 nthread = 1
391 ithread = 0
392
393 CALL reallocate(cores, 1, natom)
394 npme = 0
395 cores = 0
396
397 ! prepare core function
398 DO j = 1, lppl
399 SELECT CASE (j)
400 CASE (1)
401 pab(1, 1) = cexp_ppl(1)
402 CASE (2)
403 n = coset(2, 0, 0)
404 pab(n, 1) = cexp_ppl(2)
405 n = coset(0, 2, 0)
406 pab(n, 1) = cexp_ppl(2)
407 n = coset(0, 0, 2)
408 pab(n, 1) = cexp_ppl(2)
409 CASE (3)
410 n = coset(4, 0, 0)
411 pab(n, 1) = cexp_ppl(3)
412 n = coset(0, 4, 0)
413 pab(n, 1) = cexp_ppl(3)
414 n = coset(0, 0, 4)
415 pab(n, 1) = cexp_ppl(3)
416 n = coset(2, 2, 0)
417 pab(n, 1) = 2._dp*cexp_ppl(3)
418 n = coset(2, 0, 2)
419 pab(n, 1) = 2._dp*cexp_ppl(3)
420 n = coset(0, 2, 2)
421 pab(n, 1) = 2._dp*cexp_ppl(3)
422 CASE (4)
423 n = coset(6, 0, 0)
424 pab(n, 1) = cexp_ppl(4)
425 n = coset(0, 6, 0)
426 pab(n, 1) = cexp_ppl(4)
427 n = coset(0, 0, 6)
428 pab(n, 1) = cexp_ppl(4)
429 n = coset(4, 2, 0)
430 pab(n, 1) = 3._dp*cexp_ppl(4)
431 n = coset(4, 0, 2)
432 pab(n, 1) = 3._dp*cexp_ppl(4)
433 n = coset(2, 4, 0)
434 pab(n, 1) = 3._dp*cexp_ppl(4)
435 n = coset(2, 0, 4)
436 pab(n, 1) = 3._dp*cexp_ppl(4)
437 n = coset(0, 4, 2)
438 pab(n, 1) = 3._dp*cexp_ppl(4)
439 n = coset(0, 2, 4)
440 pab(n, 1) = 3._dp*cexp_ppl(4)
441 n = coset(2, 2, 2)
442 pab(n, 1) = 6._dp*cexp_ppl(4)
443 CASE DEFAULT
444 cpabort("")
445 END SELECT
446 END DO
447
448 DO iatom = 1, natom
449 atom_a = atom_list(iatom)
450 ra(:) = pbc(particle_set(atom_a)%r, cell)
451 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
452 ! replicated realspace grid, split the atoms up between procs
453 IF (modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
454 npme = npme + 1
455 cores(npme) = iatom
456 END IF
457 ELSE
458 npme = npme + 1
459 cores(npme) = iatom
460 END IF
461 END DO
462
463 IF (npme .GT. 0) THEN
464 DO j = 1, npme
465
466 iatom = cores(j)
467 atom_a = atom_list(iatom)
468 ra(:) = pbc(particle_set(atom_a)%r, cell)
469 subpatch_pattern = 0
470 ni = 2*lppl - 2
471
472 radius = exp_radius_very_extended(la_min=0, la_max=ni, &
473 lb_min=0, lb_max=0, &
474 ra=ra, rb=ra, rp=ra, &
475 zetp=alpha, eps=eps_rho_rspace, &
476 pab=pab, o1=0, o2=0, & ! without map_consistent
477 prefactor=1.0_dp, cutoff=0.0_dp)
478
479 CALL collocate_pgf_product(ni, alpha, 0, 0, 0.0_dp, 0, ra, &
480 (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, rs_rho, &
481 radius=radius, ga_gb_function=grid_func_ab, &
482 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
483
484 END DO
485 END IF
486
487 DEALLOCATE (pab)
488
489 END DO
490
491 IF (ASSOCIATED(cores)) THEN
492 DEALLOCATE (cores)
493 END IF
494
495 CALL transfer_rs2pw(rs_rho, vppl)
496
497 CALL timestop(handle)
498
499 END SUBROUTINE calculate_ppl_grid
500
501! **************************************************************************************************
502!> \brief Collocates the fitted lri density on a grid.
503!> \param lri_rho_g ...
504!> \param lri_rho_r ...
505!> \param qs_env ...
506!> \param lri_coef ...
507!> \param total_rho ...
508!> \param basis_type ...
509!> \param exact_1c_terms ...
510!> \param pmat replicated block diagonal density matrix (optional)
511!> \param atomlist list of atoms to be included (optional)
512!> \par History
513!> 04.2013
514!> \author Dorothea Golze
515! **************************************************************************************************
516 SUBROUTINE calculate_lri_rho_elec(lri_rho_g, lri_rho_r, qs_env, &
517 lri_coef, total_rho, basis_type, exact_1c_terms, pmat, atomlist)
518
519 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: lri_rho_g
520 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: lri_rho_r
521 TYPE(qs_environment_type), POINTER :: qs_env
522 TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
523 REAL(kind=dp), INTENT(OUT) :: total_rho
524 CHARACTER(len=*), INTENT(IN) :: basis_type
525 LOGICAL, INTENT(IN) :: exact_1c_terms
526 TYPE(dbcsr_type), OPTIONAL :: pmat
527 INTEGER, DIMENSION(:), OPTIONAL :: atomlist
528
529 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_lri_rho_elec'
530
531 INTEGER :: atom_a, group_size, handle, iatom, igrid_level, ikind, ipgf, iset, jpgf, jset, &
532 m1, maxco, maxsgf_set, my_pos, na1, natom, nb1, ncoa, ncob, nseta, offset, sgfa, sgfb
533 INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, nsgfa
534 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
535 LOGICAL :: found
536 LOGICAL, ALLOCATABLE, DIMENSION(:) :: map_it
537 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: map_it2
538 REAL(kind=dp) :: eps_rho_rspace, radius, zetp
539 REAL(kind=dp), DIMENSION(3) :: ra
540 REAL(kind=dp), DIMENSION(:), POINTER :: aci
541 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block, pab, sphi_a, work, zeta
542 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
543 TYPE(cell_type), POINTER :: cell
544 TYPE(dft_control_type), POINTER :: dft_control
545 TYPE(gridlevel_info_type), POINTER :: gridlevel_info
546 TYPE(gto_basis_set_type), POINTER :: lri_basis_set, orb_basis_set
547 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
548 TYPE(pw_env_type), POINTER :: pw_env
549 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
550 TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_gspace
551 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_rspace
552 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
553 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
554 TYPE(realspace_grid_type), POINTER :: rs_grid
555
556 NULLIFY (aci, atomic_kind_set, qs_kind_set, atom_list, cell, &
557 dft_control, first_sgfa, gridlevel_info, la_max, &
558 la_min, lri_basis_set, npgfa, nsgfa, &
559 pab, particle_set, pw_env, pw_pools, rs_grid, rs_rho, sphi_a, &
560 work, zeta)
561
562 CALL timeset(routinen, handle)
563
564 IF (exact_1c_terms) THEN
565 cpassert(PRESENT(pmat))
566 END IF
567
568 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
569 atomic_kind_set=atomic_kind_set, &
570 cell=cell, particle_set=particle_set, &
571 pw_env=pw_env, &
572 dft_control=dft_control)
573
574 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
575 gridlevel_info => pw_env%gridlevel_info
576
577 ! *** set up the pw multi-grids *** !
578 cpassert(ASSOCIATED(pw_env))
579 CALL pw_env_get(pw_env=pw_env, rs_grids=rs_rho, pw_pools=pw_pools)
580
581 CALL pw_pools_create_pws(pw_pools, mgrid_rspace)
582
583 CALL pw_pools_create_pws(pw_pools, mgrid_gspace)
584
585 ! *** set up the rs multi-grids *** !
586 DO igrid_level = 1, gridlevel_info%ngrid_levels
587 CALL rs_grid_zero(rs_rho(igrid_level))
588 END DO
589
590 !take maxco from the LRI basis set!
591 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
592 maxco=maxco, basis_type=basis_type)
593
594 ALLOCATE (pab(maxco, 1))
595 offset = 0
596 my_pos = mgrid_rspace(1)%pw_grid%para%group%mepos
597 group_size = mgrid_rspace(1)%pw_grid%para%group%num_pe
598
599 DO ikind = 1, SIZE(atomic_kind_set)
600
601 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
602 CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type=basis_type)
603
604 !Take the lri basis set here!
605 CALL get_gto_basis_set(gto_basis_set=lri_basis_set, lmax=la_max, &
606 lmin=la_min, zet=zeta, nset=nseta, npgf=npgfa, &
607 sphi=sphi_a, first_sgf=first_sgfa, nsgf_set=nsgfa)
608
609 DO iatom = 1, natom
610 atom_a = atom_list(iatom)
611 IF (PRESENT(atomlist)) THEN
612 IF (atomlist(atom_a) == 0) cycle
613 END IF
614 ra(:) = pbc(particle_set(atom_a)%r, cell)
615 aci => lri_coef(ikind)%acoef(iatom, :)
616
617 m1 = maxval(npgfa(1:nseta))
618 ALLOCATE (map_it(m1))
619 DO iset = 1, nseta
620 ! collocate this set locally?
621 map_it = .false.
622 DO ipgf = 1, npgfa(iset)
623 igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
624 rs_grid => rs_rho(igrid_level)
625 map_it(ipgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
626 END DO
627 offset = offset + 1
628
629 IF (any(map_it(1:npgfa(iset)))) THEN
630 sgfa = first_sgfa(1, iset)
631 ncoa = npgfa(iset)*ncoset(la_max(iset))
632 m1 = sgfa + nsgfa(iset) - 1
633 ALLOCATE (work(nsgfa(iset), 1))
634 work(1:nsgfa(iset), 1) = aci(sgfa:m1)
635 pab = 0._dp
636
637 CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), 1.0_dp, lri_basis_set%sphi(1, sgfa), &
638 SIZE(lri_basis_set%sphi, 1), work(1, 1), SIZE(work, 1), 0.0_dp, pab(1, 1), &
639 SIZE(pab, 1))
640
641 DO ipgf = 1, npgfa(iset)
642 na1 = (ipgf - 1)*ncoset(la_max(iset))
643 igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
644 rs_grid => rs_rho(igrid_level)
645 IF (map_it(ipgf)) THEN
646 radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
647 lb_min=0, lb_max=0, &
648 ra=ra, rb=ra, rp=ra, &
649 zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
650 prefactor=1.0_dp, cutoff=1.0_dp)
651
652 CALL collocate_pgf_product(la_max=la_max(iset), &
653 zeta=zeta(ipgf, iset), &
654 la_min=la_min(iset), &
655 lb_max=0, zetb=0.0_dp, lb_min=0, &
656 ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
657 scale=1._dp, &
658 pab=pab, o1=na1, o2=0, &
659 rsgrid=rs_grid, &
660 radius=radius, &
661 ga_gb_function=grid_func_ab)
662 END IF
663 END DO
664 DEALLOCATE (work)
665 END IF
666 END DO
667 DEALLOCATE (map_it)
668 END DO
669 END DO
670
671 DEALLOCATE (pab)
672
673 ! process the one-center terms
674 IF (exact_1c_terms) THEN
675 ! find maximum numbers
676 offset = 0
677 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
678 maxco=maxco, &
679 maxsgf_set=maxsgf_set, &
680 basis_type="ORB")
681 ALLOCATE (pab(maxco, maxco), work(maxco, maxsgf_set))
682
683 DO ikind = 1, SIZE(atomic_kind_set)
684 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
685 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
686 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, lmax=la_max, &
687 lmin=la_min, zet=zeta, nset=nseta, npgf=npgfa, &
688 sphi=sphi_a, first_sgf=first_sgfa, nsgf_set=nsgfa)
689 DO iatom = 1, natom
690 atom_a = atom_list(iatom)
691 ra(:) = pbc(particle_set(atom_a)%r, cell)
692 CALL dbcsr_get_block_p(matrix=pmat, row=atom_a, col=atom_a, block=p_block, found=found)
693 m1 = maxval(npgfa(1:nseta))
694 ALLOCATE (map_it2(m1, m1))
695 DO iset = 1, nseta
696 DO jset = 1, nseta
697 ! processor mappint
698 map_it2 = .false.
699 DO ipgf = 1, npgfa(iset)
700 DO jpgf = 1, npgfa(jset)
701 zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
702 igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
703 rs_grid => rs_rho(igrid_level)
704 map_it2(ipgf, jpgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
705 END DO
706 END DO
707 offset = offset + 1
708 !
709 IF (any(map_it2(1:npgfa(iset), 1:npgfa(jset)))) THEN
710 ncoa = npgfa(iset)*ncoset(la_max(iset))
711 sgfa = first_sgfa(1, iset)
712 ncob = npgfa(jset)*ncoset(la_max(jset))
713 sgfb = first_sgfa(1, jset)
714 ! decontract density block
715 CALL dgemm("N", "N", ncoa, nsgfa(jset), nsgfa(iset), &
716 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
717 p_block(sgfa, sgfb), SIZE(p_block, 1), &
718 0.0_dp, work(1, 1), maxco)
719 CALL dgemm("N", "T", ncoa, ncob, nsgfa(jset), &
720 1.0_dp, work(1, 1), maxco, &
721 sphi_a(1, sgfb), SIZE(sphi_a, 1), &
722 0.0_dp, pab(1, 1), maxco)
723 DO ipgf = 1, npgfa(iset)
724 DO jpgf = 1, npgfa(jset)
725 zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
726 igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
727 rs_grid => rs_rho(igrid_level)
728
729 na1 = (ipgf - 1)*ncoset(la_max(iset))
730 nb1 = (jpgf - 1)*ncoset(la_max(jset))
731
732 IF (map_it2(ipgf, jpgf)) THEN
733 radius = exp_radius_very_extended(la_min=la_min(iset), &
734 la_max=la_max(iset), &
735 lb_min=la_min(jset), &
736 lb_max=la_max(jset), &
737 ra=ra, rb=ra, rp=ra, &
738 zetp=zetp, eps=eps_rho_rspace, &
739 prefactor=1.0_dp, cutoff=1.0_dp)
740
741 CALL collocate_pgf_product( &
742 la_max(iset), zeta(ipgf, iset), la_min(iset), &
743 la_max(jset), zeta(jpgf, jset), la_min(jset), &
744 ra, (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, na1, nb1, &
745 rs_grid, &
746 radius=radius, ga_gb_function=grid_func_ab)
747 END IF
748 END DO
749 END DO
750 END IF
751 END DO
752 END DO
753 DEALLOCATE (map_it2)
754 !
755 END DO
756 END DO
757 DEALLOCATE (pab, work)
758 END IF
759
760 CALL pw_zero(lri_rho_g)
761 CALL pw_zero(lri_rho_r)
762
763 DO igrid_level = 1, gridlevel_info%ngrid_levels
764 CALL pw_zero(mgrid_rspace(igrid_level))
765 CALL transfer_rs2pw(rs=rs_rho(igrid_level), &
766 pw=mgrid_rspace(igrid_level))
767 END DO
768
769 DO igrid_level = 1, gridlevel_info%ngrid_levels
770 CALL pw_zero(mgrid_gspace(igrid_level))
771 CALL pw_transfer(mgrid_rspace(igrid_level), &
772 mgrid_gspace(igrid_level))
773 CALL pw_axpy(mgrid_gspace(igrid_level), lri_rho_g)
774 END DO
775 CALL pw_transfer(lri_rho_g, lri_rho_r)
776 total_rho = pw_integrate_function(lri_rho_r, isign=-1)
777
778 ! *** give back the multi-grids *** !
779 CALL pw_pools_give_back_pws(pw_pools, mgrid_gspace)
780 CALL pw_pools_give_back_pws(pw_pools, mgrid_rspace)
781
782 CALL timestop(handle)
783
784 END SUBROUTINE calculate_lri_rho_elec
785
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 SUBROUTINE calculate_rho_core_r3d_rs (rho_core, total_rho, qs_env, calpha, ccore, only_nopaw)
796
797 TYPE(pw_r3d_rs_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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: pab
813 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
814 TYPE(cell_type), POINTER :: cell
815 TYPE(dft_control_type), POINTER :: dft_control
816 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 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
821 TYPE(realspace_grid_type), POINTER :: rs_rho
822
823 CALL timeset(routinen, handle)
824 NULLIFY (cell, dft_control, pab, atomic_kind_set, qs_kind_set, particle_set, &
825 atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
826 ALLOCATE (pab(1, 1))
827
828 my_only_nopaw = .false.
829 IF (PRESENT(only_nopaw)) my_only_nopaw = only_nopaw
830 IF (PRESENT(calpha)) THEN
831 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 pw_env=pw_env)
841 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
842 auxbas_pw_pool=auxbas_pw_pool)
843 ! be careful in parallel nsmax is chosen with multigrid in mind!
844 CALL rs_grid_zero(rs_rho)
845
846 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
847
848 DO ikind = 1, SIZE(atomic_kind_set)
849 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
850 IF (PRESENT(calpha)) THEN
851 alpha = calpha(ikind)
852 pab(1, 1) = ccore(ikind)
853 ELSE
854 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom, &
855 alpha_core_charge=alpha, ccore_charge=pab(1, 1))
856 END IF
857
858 IF (my_only_nopaw .AND. paw_atom) cycle
859 IF (alpha == 0.0_dp .OR. pab(1, 1) == 0.0_dp) cycle
860
861 nthread = 1
862 ithread = 0
863
864 CALL reallocate(cores, 1, natom)
865 npme = 0
866 cores = 0
867
868 DO iatom = 1, natom
869 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
870 ! replicated realspace grid, split the atoms up between procs
871 IF (modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
872 npme = npme + 1
873 cores(npme) = iatom
874 END IF
875 ELSE
876 npme = npme + 1
877 cores(npme) = iatom
878 END IF
879 END DO
880
881 IF (npme .GT. 0) THEN
882 DO j = 1, npme
883
884 iatom = cores(j)
885 atom_a = atom_list(iatom)
886 ra(:) = pbc(particle_set(atom_a)%r, cell)
887 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 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 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
899
900 END DO
901 END IF
902
903 END DO
904
905 IF (ASSOCIATED(cores)) THEN
906 DEALLOCATE (cores)
907 END IF
908 DEALLOCATE (pab)
909
910 CALL auxbas_pw_pool%create_pw(rhoc_r)
911
912 CALL transfer_rs2pw(rs_rho, rhoc_r)
913
914 total_rho = pw_integrate_function(rhoc_r, isign=-1)
915
916 CALL pw_transfer(rhoc_r, rho_core)
917
918 CALL auxbas_pw_pool%give_back_pw(rhoc_r)
919
920 CALL timestop(handle)
921
922 END SUBROUTINE calculate_rho_core_r3d_rs
923! **************************************************************************************************
924!> \brief computes the density of the core charges on the grid
925!> \param rho_core ...
926!> \param total_rho ...
927!> \param qs_env ...
928!> \param calpha ...
929!> \param ccore ...
930!> \param only_nopaw ...
931! **************************************************************************************************
932 SUBROUTINE calculate_rho_core_c1d_gs (rho_core, total_rho, qs_env, calpha, ccore, only_nopaw)
933
934 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_core
935 REAL(KIND=dp), INTENT(OUT) :: total_rho
936 TYPE(qs_environment_type), POINTER :: qs_env
937 REAL(kind=dp), DIMENSION(:), OPTIONAL :: calpha, ccore
938 LOGICAL, INTENT(IN), OPTIONAL :: only_nopaw
939
940 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_rho_core'
941
942 INTEGER :: atom_a, handle, iatom, ikind, ithread, &
943 j, natom, npme, nthread, &
944 subpatch_pattern
945 INTEGER, DIMENSION(:), POINTER :: atom_list, cores
946 LOGICAL :: my_only_nopaw, paw_atom
947 REAL(kind=dp) :: alpha, eps_rho_rspace, radius
948 REAL(kind=dp), DIMENSION(3) :: ra
949 REAL(kind=dp), DIMENSION(:, :), POINTER :: pab
950 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
951 TYPE(cell_type), POINTER :: cell
952 TYPE(dft_control_type), POINTER :: dft_control
953 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
954 TYPE(pw_env_type), POINTER :: pw_env
955 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
956 TYPE(pw_r3d_rs_type) :: rhoc_r
957 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
958 TYPE(realspace_grid_type), POINTER :: rs_rho
959
960 CALL timeset(routinen, handle)
961 NULLIFY (cell, dft_control, pab, atomic_kind_set, qs_kind_set, particle_set, &
962 atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
963 ALLOCATE (pab(1, 1))
964
965 my_only_nopaw = .false.
966 IF (PRESENT(only_nopaw)) my_only_nopaw = only_nopaw
967 IF (PRESENT(calpha)) THEN
968 cpassert(PRESENT(ccore))
969 END IF
970
971 CALL get_qs_env(qs_env=qs_env, &
972 atomic_kind_set=atomic_kind_set, &
973 qs_kind_set=qs_kind_set, &
974 cell=cell, &
975 dft_control=dft_control, &
976 particle_set=particle_set, &
977 pw_env=pw_env)
978 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
979 auxbas_pw_pool=auxbas_pw_pool)
980 ! be careful in parallel nsmax is chosen with multigrid in mind!
981 CALL rs_grid_zero(rs_rho)
982
983 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
984
985 DO ikind = 1, SIZE(atomic_kind_set)
986 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
987 IF (PRESENT(calpha)) THEN
988 alpha = calpha(ikind)
989 pab(1, 1) = ccore(ikind)
990 ELSE
991 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom, &
992 alpha_core_charge=alpha, ccore_charge=pab(1, 1))
993 END IF
994
995 IF (my_only_nopaw .AND. paw_atom) cycle
996 IF (alpha == 0.0_dp .OR. pab(1, 1) == 0.0_dp) cycle
997
998 nthread = 1
999 ithread = 0
1000
1001 CALL reallocate(cores, 1, natom)
1002 npme = 0
1003 cores = 0
1004
1005 DO iatom = 1, natom
1006 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
1007 ! replicated realspace grid, split the atoms up between procs
1008 IF (modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
1009 npme = npme + 1
1010 cores(npme) = iatom
1011 END IF
1012 ELSE
1013 npme = npme + 1
1014 cores(npme) = iatom
1015 END IF
1016 END DO
1017
1018 IF (npme .GT. 0) THEN
1019 DO j = 1, npme
1020
1021 iatom = cores(j)
1022 atom_a = atom_list(iatom)
1023 ra(:) = pbc(particle_set(atom_a)%r, cell)
1024 subpatch_pattern = 0
1025 radius = exp_radius_very_extended(la_min=0, la_max=0, &
1026 lb_min=0, lb_max=0, &
1027 ra=ra, rb=ra, rp=ra, &
1028 zetp=alpha, eps=eps_rho_rspace, &
1029 pab=pab, o1=0, o2=0, & ! without map_consistent
1030 prefactor=-1.0_dp, cutoff=0.0_dp)
1031
1032 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
1033 (/0.0_dp, 0.0_dp, 0.0_dp/), -1.0_dp, pab, 0, 0, rs_rho, &
1034 radius=radius, ga_gb_function=grid_func_ab, &
1035 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
1036
1037 END DO
1038 END IF
1039
1040 END DO
1041
1042 IF (ASSOCIATED(cores)) THEN
1043 DEALLOCATE (cores)
1044 END IF
1045 DEALLOCATE (pab)
1046
1047 CALL auxbas_pw_pool%create_pw(rhoc_r)
1048
1049 CALL transfer_rs2pw(rs_rho, rhoc_r)
1050
1051 total_rho = pw_integrate_function(rhoc_r, isign=-1)
1052
1053 CALL pw_transfer(rhoc_r, rho_core)
1054
1055 CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1056
1057 CALL timestop(handle)
1058
1059 END SUBROUTINE calculate_rho_core_c1d_gs
1060
1061! *****************************************************************************
1062!> \brief Computes the derivative of the density of the core charges with
1063!> respect to the nuclear coordinates on the grid.
1064!> \param drho_core The resulting density derivative
1065!> \param qs_env ...
1066!> \param beta Derivative direction
1067!> \param lambda Atom index
1068!> \note SL November 2014, ED 2021
1069! **************************************************************************************************
1070 SUBROUTINE calculate_drho_core(drho_core, qs_env, beta, lambda)
1071
1072 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: drho_core
1073 TYPE(qs_environment_type), POINTER :: qs_env
1074 INTEGER, INTENT(IN) :: beta, lambda
1075
1076 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_drho_core'
1077
1078 INTEGER :: atom_a, dabqadb_func, handle, iatom, &
1079 ikind, ithread, j, natom, npme, &
1080 nthread, subpatch_pattern
1081 INTEGER, DIMENSION(:), POINTER :: atom_list, cores
1082 REAL(kind=dp) :: alpha, eps_rho_rspace, radius
1083 REAL(kind=dp), DIMENSION(3) :: ra
1084 REAL(kind=dp), DIMENSION(:, :), POINTER :: pab
1085 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1086 TYPE(cell_type), POINTER :: cell
1087 TYPE(dft_control_type), POINTER :: dft_control
1088 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1089 TYPE(pw_env_type), POINTER :: pw_env
1090 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1091 TYPE(pw_r3d_rs_type) :: rhoc_r
1092 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1093 TYPE(realspace_grid_type), POINTER :: rs_rho
1094
1095 CALL timeset(routinen, handle)
1096 NULLIFY (cell, dft_control, pab, atomic_kind_set, qs_kind_set, particle_set, &
1097 atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
1098 ALLOCATE (pab(1, 1))
1099
1100 CALL get_qs_env(qs_env=qs_env, &
1101 atomic_kind_set=atomic_kind_set, &
1102 qs_kind_set=qs_kind_set, &
1103 cell=cell, &
1104 dft_control=dft_control, &
1105 particle_set=particle_set, &
1106 pw_env=pw_env)
1107 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
1108 auxbas_pw_pool=auxbas_pw_pool)
1109 ! be careful in parallel nsmax is chosen with multigrid in mind!
1110 CALL rs_grid_zero(rs_rho)
1111
1112 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1113
1114 SELECT CASE (beta)
1115 CASE (1)
1116 dabqadb_func = grid_func_core_x
1117 CASE (2)
1118 dabqadb_func = grid_func_core_y
1119 CASE (3)
1120 dabqadb_func = grid_func_core_z
1121 CASE DEFAULT
1122 cpabort("invalid beta")
1123 END SELECT
1124 DO ikind = 1, SIZE(atomic_kind_set)
1125 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
1126 CALL get_qs_kind(qs_kind_set(ikind), &
1127 alpha_core_charge=alpha, ccore_charge=pab(1, 1))
1128
1129 IF (alpha == 0.0_dp .OR. pab(1, 1) == 0.0_dp) cycle
1130
1131 nthread = 1
1132 ithread = 0
1133
1134 CALL reallocate(cores, 1, natom)
1135 npme = 0
1136 cores = 0
1137
1138 DO iatom = 1, natom
1139 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
1140 ! replicated realspace grid, split the atoms up between procs
1141 IF (modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
1142 npme = npme + 1
1143 cores(npme) = iatom
1144 END IF
1145 ELSE
1146 npme = npme + 1
1147 cores(npme) = iatom
1148 END IF
1149 END DO
1150
1151 IF (npme .GT. 0) THEN
1152 DO j = 1, npme
1153
1154 iatom = cores(j)
1155 atom_a = atom_list(iatom)
1156 IF (atom_a /= lambda) cycle
1157 ra(:) = pbc(particle_set(atom_a)%r, cell)
1158 subpatch_pattern = 0
1159 radius = exp_radius_very_extended(la_min=0, la_max=0, &
1160 lb_min=0, lb_max=0, &
1161 ra=ra, rb=ra, rp=ra, &
1162 zetp=alpha, eps=eps_rho_rspace, &
1163 pab=pab, o1=0, o2=0, & ! without map_consistent
1164 prefactor=-1.0_dp, cutoff=0.0_dp)
1165
1166 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
1167 (/0.0_dp, 0.0_dp, 0.0_dp/), -1.0_dp, pab, 0, 0, rs_rho, &
1168 radius=radius, ga_gb_function=dabqadb_func, &
1169 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
1170
1171 END DO
1172 END IF
1173
1174 END DO
1175
1176 IF (ASSOCIATED(cores)) THEN
1177 DEALLOCATE (cores)
1178 END IF
1179 DEALLOCATE (pab)
1180
1181 CALL auxbas_pw_pool%create_pw(rhoc_r)
1182
1183 CALL transfer_rs2pw(rs_rho, rhoc_r)
1184
1185 CALL pw_transfer(rhoc_r, drho_core)
1186
1187 CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1188
1189 CALL timestop(handle)
1190
1191 END SUBROUTINE calculate_drho_core
1192
1193! **************************************************************************************************
1194!> \brief collocate a single Gaussian on the grid
1195!> \param rho_gb charge density generated by a single gaussian
1196!> \param qs_env qs environment
1197!> \param iatom_in atom index
1198!> \par History
1199!> 12.2011 created
1200!> \author Dorothea Golze
1201! **************************************************************************************************
1202 SUBROUTINE calculate_rho_single_gaussian(rho_gb, qs_env, iatom_in)
1203
1204 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_gb
1205 TYPE(qs_environment_type), POINTER :: qs_env
1206 INTEGER, INTENT(IN) :: iatom_in
1207
1208 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_rho_single_gaussian'
1209
1210 INTEGER :: atom_a, handle, iatom, npme, &
1211 subpatch_pattern
1212 REAL(kind=dp) :: eps_rho_rspace, radius
1213 REAL(kind=dp), DIMENSION(3) :: ra
1214 REAL(kind=dp), DIMENSION(:, :), POINTER :: pab
1215 TYPE(cell_type), POINTER :: cell
1216 TYPE(dft_control_type), POINTER :: dft_control
1217 TYPE(pw_env_type), POINTER :: pw_env
1218 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1219 TYPE(pw_r3d_rs_type) :: rhoc_r
1220 TYPE(realspace_grid_type), POINTER :: rs_rho
1221
1222 CALL timeset(routinen, handle)
1223 NULLIFY (cell, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool)
1224
1225 ALLOCATE (pab(1, 1))
1226
1227 CALL get_qs_env(qs_env=qs_env, &
1228 cell=cell, &
1229 dft_control=dft_control, &
1230 pw_env=pw_env)
1231 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
1232 auxbas_pw_pool=auxbas_pw_pool)
1233 CALL rs_grid_zero(rs_rho)
1234
1235 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1236 pab(1, 1) = 1.0_dp
1237 iatom = iatom_in
1238
1239 npme = 0
1240
1241 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
1242 IF (modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
1243 npme = npme + 1
1244 END IF
1245 ELSE
1246 npme = npme + 1
1247 END IF
1248
1249 IF (npme .GT. 0) THEN
1250 atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom)
1251 ra(:) = pbc(qs_env%qmmm_env_qm%image_charge_pot%particles_all(atom_a)%r, cell)
1252 subpatch_pattern = 0
1253 radius = exp_radius_very_extended(la_min=0, la_max=0, &
1254 lb_min=0, lb_max=0, &
1255 ra=ra, rb=ra, rp=ra, &
1256 zetp=qs_env%qmmm_env_qm%image_charge_pot%eta, &
1257 eps=eps_rho_rspace, &
1258 pab=pab, o1=0, o2=0, & ! without map_consistent
1259 prefactor=1.0_dp, cutoff=0.0_dp)
1260
1261 CALL collocate_pgf_product(0, qs_env%qmmm_env_qm%image_charge_pot%eta, &
1262 0, 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, rs_rho, &
1263 radius=radius, ga_gb_function=grid_func_ab, &
1264 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
1265 END IF
1266
1267 DEALLOCATE (pab)
1268
1269 CALL auxbas_pw_pool%create_pw(rhoc_r)
1270
1271 CALL transfer_rs2pw(rs_rho, rhoc_r)
1272
1273 CALL pw_transfer(rhoc_r, rho_gb)
1274
1275 CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1276
1277 CALL timestop(handle)
1278
1279 END SUBROUTINE calculate_rho_single_gaussian
1280
1281! **************************************************************************************************
1282!> \brief computes the image charge density on the grid (including coeffcients)
1283!> \param rho_metal image charge density
1284!> \param coeff expansion coefficients of the image charge density, i.e.
1285!> rho_metal=sum_a c_a*g_a
1286!> \param total_rho_metal total induced image charge density
1287!> \param qs_env qs environment
1288!> \par History
1289!> 01.2012 created
1290!> \author Dorothea Golze
1291! **************************************************************************************************
1292 SUBROUTINE calculate_rho_metal(rho_metal, coeff, total_rho_metal, qs_env)
1293
1294 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_metal
1295 REAL(kind=dp), DIMENSION(:), POINTER :: coeff
1296 REAL(kind=dp), INTENT(OUT), OPTIONAL :: total_rho_metal
1297 TYPE(qs_environment_type), POINTER :: qs_env
1298
1299 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_rho_metal'
1300
1301 INTEGER :: atom_a, handle, iatom, j, natom, npme, &
1302 subpatch_pattern
1303 INTEGER, DIMENSION(:), POINTER :: cores
1304 REAL(kind=dp) :: eps_rho_rspace, radius
1305 REAL(kind=dp), DIMENSION(3) :: ra
1306 REAL(kind=dp), DIMENSION(:, :), POINTER :: pab
1307 TYPE(cell_type), POINTER :: cell
1308 TYPE(dft_control_type), POINTER :: dft_control
1309 TYPE(pw_env_type), POINTER :: pw_env
1310 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1311 TYPE(pw_r3d_rs_type) :: rhoc_r
1312 TYPE(realspace_grid_type), POINTER :: rs_rho
1313
1314 CALL timeset(routinen, handle)
1315
1316 NULLIFY (cell, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool, cores)
1317
1318 ALLOCATE (pab(1, 1))
1319
1320 CALL get_qs_env(qs_env=qs_env, &
1321 cell=cell, &
1322 dft_control=dft_control, &
1323 pw_env=pw_env)
1324 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
1325 auxbas_pw_pool=auxbas_pw_pool)
1326 CALL rs_grid_zero(rs_rho)
1327
1328 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1329 pab(1, 1) = 1.0_dp
1330
1331 natom = SIZE(qs_env%qmmm_env_qm%image_charge_pot%image_mm_list)
1332
1333 CALL reallocate(cores, 1, natom)
1334 npme = 0
1335 cores = 0
1336
1337 DO iatom = 1, natom
1338 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
1339 IF (modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
1340 npme = npme + 1
1341 cores(npme) = iatom
1342 END IF
1343 ELSE
1344 npme = npme + 1
1345 cores(npme) = iatom
1346 END IF
1347 END DO
1348
1349 IF (npme .GT. 0) THEN
1350 DO j = 1, npme
1351 iatom = cores(j)
1352 atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom)
1353 ra(:) = pbc(qs_env%qmmm_env_qm%image_charge_pot%particles_all(atom_a)%r, cell)
1354 subpatch_pattern = 0
1355 radius = exp_radius_very_extended(la_min=0, la_max=0, &
1356 lb_min=0, lb_max=0, &
1357 ra=ra, rb=ra, rp=ra, &
1358 zetp=qs_env%qmmm_env_qm%image_charge_pot%eta, &
1359 eps=eps_rho_rspace, &
1360 pab=pab, o1=0, o2=0, & ! without map_consistent
1361 prefactor=coeff(iatom), cutoff=0.0_dp)
1362
1363 CALL collocate_pgf_product( &
1364 0, qs_env%qmmm_env_qm%image_charge_pot%eta, &
1365 0, 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), coeff(iatom), pab, 0, 0, rs_rho, &
1366 radius=radius, ga_gb_function=grid_func_ab, &
1367 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
1368 END DO
1369 END IF
1370
1371 DEALLOCATE (pab, cores)
1372
1373 CALL auxbas_pw_pool%create_pw(rhoc_r)
1374
1375 CALL transfer_rs2pw(rs_rho, rhoc_r)
1376
1377 IF (PRESENT(total_rho_metal)) &
1378 !minus sign: account for the fact that rho_metal has opposite sign
1379 total_rho_metal = pw_integrate_function(rhoc_r, isign=-1)
1380
1381 CALL pw_transfer(rhoc_r, rho_metal)
1382 CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1383
1384 CALL timestop(handle)
1385
1386 END SUBROUTINE calculate_rho_metal
1387
1388! **************************************************************************************************
1389!> \brief collocate a single Gaussian on the grid for periodic RESP fitting
1390!> \param rho_gb charge density generated by a single gaussian
1391!> \param qs_env qs environment
1392!> \param eta width of single Gaussian
1393!> \param iatom_in atom index
1394!> \par History
1395!> 06.2012 created
1396!> \author Dorothea Golze
1397! **************************************************************************************************
1398 SUBROUTINE calculate_rho_resp_single(rho_gb, qs_env, eta, iatom_in)
1399
1400 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_gb
1401 TYPE(qs_environment_type), POINTER :: qs_env
1402 REAL(kind=dp), INTENT(IN) :: eta
1403 INTEGER, INTENT(IN) :: iatom_in
1404
1405 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_rho_resp_single'
1406
1407 INTEGER :: handle, iatom, npme, subpatch_pattern
1408 REAL(kind=dp) :: eps_rho_rspace, radius
1409 REAL(kind=dp), DIMENSION(3) :: ra
1410 REAL(kind=dp), DIMENSION(:, :), POINTER :: pab
1411 TYPE(cell_type), POINTER :: cell
1412 TYPE(dft_control_type), POINTER :: dft_control
1413 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1414 TYPE(pw_env_type), POINTER :: pw_env
1415 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1416 TYPE(pw_r3d_rs_type) :: rhoc_r
1417 TYPE(realspace_grid_type), POINTER :: rs_rho
1418
1419 CALL timeset(routinen, handle)
1420 NULLIFY (cell, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool, &
1421 particle_set)
1422
1423 ALLOCATE (pab(1, 1))
1424
1425 CALL get_qs_env(qs_env=qs_env, &
1426 cell=cell, &
1427 dft_control=dft_control, &
1428 particle_set=particle_set, &
1429 pw_env=pw_env)
1430 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
1431 auxbas_pw_pool=auxbas_pw_pool)
1432 CALL rs_grid_zero(rs_rho)
1433
1434 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1435 pab(1, 1) = 1.0_dp
1436 iatom = iatom_in
1437
1438 npme = 0
1439
1440 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
1441 IF (modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
1442 npme = npme + 1
1443 END IF
1444 ELSE
1445 npme = npme + 1
1446 END IF
1447
1448 IF (npme .GT. 0) THEN
1449 ra(:) = pbc(particle_set(iatom)%r, cell)
1450 subpatch_pattern = 0
1451 radius = exp_radius_very_extended(la_min=0, la_max=0, &
1452 lb_min=0, lb_max=0, &
1453 ra=ra, rb=ra, rp=ra, &
1454 zetp=eta, eps=eps_rho_rspace, &
1455 pab=pab, o1=0, o2=0, & ! without map_consistent
1456 prefactor=1.0_dp, cutoff=0.0_dp)
1457
1458 CALL collocate_pgf_product(0, eta, 0, 0, 0.0_dp, 0, ra, &
1459 (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, rs_rho, &
1460 radius=radius, ga_gb_function=grid_func_ab, &
1461 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
1462 END IF
1463
1464 DEALLOCATE (pab)
1465
1466 CALL auxbas_pw_pool%create_pw(rhoc_r)
1467
1468 CALL transfer_rs2pw(rs_rho, rhoc_r)
1469
1470 CALL pw_transfer(rhoc_r, rho_gb)
1471
1472 CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1473
1474 CALL timestop(handle)
1475
1476 END SUBROUTINE calculate_rho_resp_single
1477
1478! **************************************************************************************************
1479!> \brief computes the RESP charge density on a grid based on the RESP charges
1480!> \param rho_resp RESP charge density
1481!> \param coeff RESP charges, take care of normalization factor
1482!> (eta/pi)**1.5 later
1483!> \param natom number of atoms
1484!> \param eta width of single Gaussian
1485!> \param qs_env qs environment
1486!> \par History
1487!> 01.2012 created
1488!> \author Dorothea Golze
1489! **************************************************************************************************
1490 SUBROUTINE calculate_rho_resp_all_r3d_rs (rho_resp, coeff, natom, eta, qs_env)
1491
1492 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho_resp
1493 REAL(KIND=dp), DIMENSION(:), POINTER :: coeff
1494 INTEGER, INTENT(IN) :: natom
1495 REAL(KIND=dp), INTENT(IN) :: eta
1496 TYPE(qs_environment_type), POINTER :: qs_env
1497
1498 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_rho_resp_all'
1499
1500 INTEGER :: handle, iatom, j, npme, subpatch_pattern
1501 INTEGER, DIMENSION(:), POINTER :: cores
1502 REAL(kind=dp) :: eps_rho_rspace, radius
1503 REAL(kind=dp), DIMENSION(3) :: ra
1504 REAL(kind=dp), DIMENSION(:, :), POINTER :: pab
1505 TYPE(cell_type), POINTER :: cell
1506 TYPE(dft_control_type), POINTER :: dft_control
1507 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1508 TYPE(pw_env_type), POINTER :: pw_env
1509 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1510 TYPE(pw_r3d_rs_type) :: rhoc_r
1511 TYPE(realspace_grid_type), POINTER :: rs_rho
1512
1513 CALL timeset(routinen, handle)
1514
1515 NULLIFY (cell, cores, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool, &
1516 particle_set)
1517
1518 ALLOCATE (pab(1, 1))
1519
1520 CALL get_qs_env(qs_env=qs_env, &
1521 cell=cell, &
1522 dft_control=dft_control, &
1523 particle_set=particle_set, &
1524 pw_env=pw_env)
1525 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
1526 auxbas_pw_pool=auxbas_pw_pool)
1527 CALL rs_grid_zero(rs_rho)
1528
1529 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1530 pab(1, 1) = 1.0_dp
1531
1532 CALL reallocate(cores, 1, natom)
1533 npme = 0
1534 cores = 0
1535
1536 DO iatom = 1, natom
1537 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
1538 IF (modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
1539 npme = npme + 1
1540 cores(npme) = iatom
1541 END IF
1542 ELSE
1543 npme = npme + 1
1544 cores(npme) = iatom
1545 END IF
1546 END DO
1547
1548 IF (npme .GT. 0) THEN
1549 DO j = 1, npme
1550 iatom = cores(j)
1551 ra(:) = pbc(particle_set(iatom)%r, cell)
1552 subpatch_pattern = 0
1553 radius = exp_radius_very_extended(la_min=0, la_max=0, &
1554 lb_min=0, lb_max=0, &
1555 ra=ra, rb=ra, rp=ra, &
1556 zetp=eta, eps=eps_rho_rspace, &
1557 pab=pab, o1=0, o2=0, & ! without map_consistent
1558 prefactor=coeff(iatom), cutoff=0.0_dp)
1559
1560 CALL collocate_pgf_product( &
1561 0, eta, &
1562 0, 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), coeff(iatom), pab, 0, 0, rs_rho, &
1563 radius=radius, ga_gb_function=grid_func_ab, &
1564 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
1565 END DO
1566 END IF
1567
1568 DEALLOCATE (pab, cores)
1569
1570 CALL auxbas_pw_pool%create_pw(rhoc_r)
1571
1572 CALL transfer_rs2pw(rs_rho, rhoc_r)
1573
1574 CALL pw_transfer(rhoc_r, rho_resp)
1575 CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1576
1577 CALL timestop(handle)
1578
1579 END SUBROUTINE calculate_rho_resp_all_r3d_rs
1580! **************************************************************************************************
1581!> \brief computes the RESP charge density on a grid based on the RESP charges
1582!> \param rho_resp RESP charge density
1583!> \param coeff RESP charges, take care of normalization factor
1584!> (eta/pi)**1.5 later
1585!> \param natom number of atoms
1586!> \param eta width of single Gaussian
1587!> \param qs_env qs environment
1588!> \par History
1589!> 01.2012 created
1590!> \author Dorothea Golze
1591! **************************************************************************************************
1592 SUBROUTINE calculate_rho_resp_all_c1d_gs (rho_resp, coeff, natom, eta, qs_env)
1593
1594 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_resp
1595 REAL(KIND=dp), DIMENSION(:), POINTER :: coeff
1596 INTEGER, INTENT(IN) :: natom
1597 REAL(KIND=dp), INTENT(IN) :: eta
1598 TYPE(qs_environment_type), POINTER :: qs_env
1599
1600 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_rho_resp_all'
1601
1602 INTEGER :: handle, iatom, j, npme, subpatch_pattern
1603 INTEGER, DIMENSION(:), POINTER :: cores
1604 REAL(kind=dp) :: eps_rho_rspace, radius
1605 REAL(kind=dp), DIMENSION(3) :: ra
1606 REAL(kind=dp), DIMENSION(:, :), POINTER :: pab
1607 TYPE(cell_type), POINTER :: cell
1608 TYPE(dft_control_type), POINTER :: dft_control
1609 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1610 TYPE(pw_env_type), POINTER :: pw_env
1611 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1612 TYPE(pw_r3d_rs_type) :: rhoc_r
1613 TYPE(realspace_grid_type), POINTER :: rs_rho
1614
1615 CALL timeset(routinen, handle)
1616
1617 NULLIFY (cell, cores, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool, &
1618 particle_set)
1619
1620 ALLOCATE (pab(1, 1))
1621
1622 CALL get_qs_env(qs_env=qs_env, &
1623 cell=cell, &
1624 dft_control=dft_control, &
1625 particle_set=particle_set, &
1626 pw_env=pw_env)
1627 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
1628 auxbas_pw_pool=auxbas_pw_pool)
1629 CALL rs_grid_zero(rs_rho)
1630
1631 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1632 pab(1, 1) = 1.0_dp
1633
1634 CALL reallocate(cores, 1, natom)
1635 npme = 0
1636 cores = 0
1637
1638 DO iatom = 1, natom
1639 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
1640 IF (modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
1641 npme = npme + 1
1642 cores(npme) = iatom
1643 END IF
1644 ELSE
1645 npme = npme + 1
1646 cores(npme) = iatom
1647 END IF
1648 END DO
1649
1650 IF (npme .GT. 0) THEN
1651 DO j = 1, npme
1652 iatom = cores(j)
1653 ra(:) = pbc(particle_set(iatom)%r, cell)
1654 subpatch_pattern = 0
1655 radius = exp_radius_very_extended(la_min=0, la_max=0, &
1656 lb_min=0, lb_max=0, &
1657 ra=ra, rb=ra, rp=ra, &
1658 zetp=eta, eps=eps_rho_rspace, &
1659 pab=pab, o1=0, o2=0, & ! without map_consistent
1660 prefactor=coeff(iatom), cutoff=0.0_dp)
1661
1662 CALL collocate_pgf_product( &
1663 0, eta, &
1664 0, 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), coeff(iatom), pab, 0, 0, rs_rho, &
1665 radius=radius, ga_gb_function=grid_func_ab, &
1666 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
1667 END DO
1668 END IF
1669
1670 DEALLOCATE (pab, cores)
1671
1672 CALL auxbas_pw_pool%create_pw(rhoc_r)
1673
1674 CALL transfer_rs2pw(rs_rho, rhoc_r)
1675
1676 CALL pw_transfer(rhoc_r, rho_resp)
1677 CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1678
1679 CALL timestop(handle)
1680
1681 END SUBROUTINE calculate_rho_resp_all_c1d_gs
1682
1683! **************************************************************************************************
1684!> \brief computes the density corresponding to a given density matrix on the grid
1685!> \param matrix_p ...
1686!> \param matrix_p_kp ...
1687!> \param rho ...
1688!> \param rho_gspace ...
1689!> \param total_rho ...
1690!> \param ks_env ...
1691!> \param soft_valid ...
1692!> \param compute_tau ...
1693!> \param compute_grad ...
1694!> \param basis_type ...
1695!> \param der_type ...
1696!> \param idir ...
1697!> \param task_list_external ...
1698!> \param pw_env_external ...
1699!> \par History
1700!> IAB (15-Feb-2010): Added OpenMP parallelisation to task loop
1701!> (c) The Numerical Algorithms Group (NAG) Ltd, 2010 on behalf of the HECToR project
1702!> Anything that is not the default ORB basis_type requires an external_task_list 12.2019, (A.Bussy)
1703!> Ole Schuett (2020): Migrated to C, see grid_api.F
1704!> \note
1705!> both rho and rho_gspace contain the new rho
1706!> (in real and g-space respectively)
1707! **************************************************************************************************
1708 SUBROUTINE calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, &
1709 ks_env, soft_valid, compute_tau, compute_grad, &
1710 basis_type, der_type, idir, task_list_external, pw_env_external)
1711
1712 TYPE(dbcsr_type), OPTIONAL, TARGET :: matrix_p
1713 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
1714 POINTER :: matrix_p_kp
1715 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho
1716 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_gspace
1717 REAL(kind=dp), INTENT(OUT), OPTIONAL :: total_rho
1718 TYPE(qs_ks_env_type), POINTER :: ks_env
1719 LOGICAL, INTENT(IN), OPTIONAL :: soft_valid, compute_tau, compute_grad
1720 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: basis_type
1721 INTEGER, INTENT(IN), OPTIONAL :: der_type, idir
1722 TYPE(task_list_type), OPTIONAL, POINTER :: task_list_external
1723 TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
1724
1725 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_rho_elec'
1726
1727 CHARACTER(LEN=default_string_length) :: my_basis_type
1728 INTEGER :: ga_gb_function, handle, ilevel, img, &
1729 nimages, nlevels
1730 LOGICAL :: any_distributed, my_compute_grad, &
1731 my_compute_tau, my_soft_valid
1732 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_images
1733 TYPE(dft_control_type), POINTER :: dft_control
1734 TYPE(mp_comm_type) :: group
1735 TYPE(pw_env_type), POINTER :: pw_env
1736 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
1737 TYPE(task_list_type), POINTER :: task_list
1738
1739 CALL timeset(routinen, handle)
1740
1741 NULLIFY (matrix_images, dft_control, pw_env, rs_rho, task_list)
1742
1743 ! Figure out which function to collocate.
1744 my_compute_tau = .false.
1745 IF (PRESENT(compute_tau)) my_compute_tau = compute_tau
1746 my_compute_grad = .false.
1747 IF (PRESENT(compute_grad)) my_compute_grad = compute_grad
1748 IF (PRESENT(der_type)) THEN
1749 SELECT CASE (der_type)
1750 CASE (orb_s)
1751 ga_gb_function = grid_func_ab
1752 CASE (orb_px)
1753 ga_gb_function = grid_func_dx
1754 CASE (orb_py)
1755 ga_gb_function = grid_func_dy
1756 CASE (orb_pz)
1757 ga_gb_function = grid_func_dz
1758 CASE (orb_dxy)
1759 ga_gb_function = grid_func_dxdy
1760 CASE (orb_dyz)
1761 ga_gb_function = grid_func_dydz
1762 CASE (orb_dzx)
1763 ga_gb_function = grid_func_dzdx
1764 CASE (orb_dx2)
1765 ga_gb_function = grid_func_dxdx
1766 CASE (orb_dy2)
1767 ga_gb_function = grid_func_dydy
1768 CASE (orb_dz2)
1769 ga_gb_function = grid_func_dzdz
1770 CASE DEFAULT
1771 cpabort("Unknown der_type")
1772 END SELECT
1773 ELSE IF (my_compute_tau) THEN
1774 ga_gb_function = grid_func_dadb
1775 ELSE IF (my_compute_grad) THEN
1776 cpassert(PRESENT(idir))
1777 SELECT CASE (idir)
1778 CASE (1)
1779 ga_gb_function = grid_func_dabpadb_x
1780 CASE (2)
1781 ga_gb_function = grid_func_dabpadb_y
1782 CASE (3)
1783 ga_gb_function = grid_func_dabpadb_z
1784 CASE DEFAULT
1785 cpabort("invalid idir")
1786 END SELECT
1787 ELSE
1788 ga_gb_function = grid_func_ab
1789 END IF
1790
1791 ! Figure out which basis_type to use.
1792 my_basis_type = "ORB" ! by default, the full density is calculated
1793 IF (PRESENT(basis_type)) my_basis_type = basis_type
1794 cpassert(my_basis_type == "ORB" .OR. PRESENT(task_list_external))
1795
1796 ! Figure out which task_list to use.
1797 my_soft_valid = .false.
1798 IF (PRESENT(soft_valid)) my_soft_valid = soft_valid
1799 IF (PRESENT(task_list_external)) THEN
1800 task_list => task_list_external
1801 ELSEIF (my_soft_valid) THEN
1802 CALL get_ks_env(ks_env, task_list_soft=task_list)
1803 ELSE
1804 CALL get_ks_env(ks_env, task_list=task_list)
1805 END IF
1806 cpassert(ASSOCIATED(task_list))
1807
1808 ! Figure out which pw_env to use.
1809 IF (PRESENT(pw_env_external)) THEN
1810 pw_env => pw_env_external
1811 ELSE
1812 CALL get_ks_env(ks_env, pw_env=pw_env)
1813 END IF
1814 cpassert(ASSOCIATED(pw_env))
1815
1816 ! Get grids.
1817 CALL pw_env_get(pw_env, rs_grids=rs_rho)
1818 nlevels = SIZE(rs_rho)
1819 group = rs_rho(1)%desc%group
1820
1821 ! Check if any of the grids is distributed.
1822 any_distributed = .false.
1823 DO ilevel = 1, nlevels
1824 any_distributed = any_distributed .OR. rs_rho(ilevel)%desc%distributed
1825 END DO
1826
1827 ! Gather all matrix images in a single array.
1828 CALL get_ks_env(ks_env, dft_control=dft_control)
1829 nimages = dft_control%nimages
1830 ALLOCATE (matrix_images(nimages))
1831 IF (PRESENT(matrix_p_kp)) THEN
1832 cpassert(.NOT. PRESENT(matrix_p))
1833 DO img = 1, nimages
1834 matrix_images(img)%matrix => matrix_p_kp(img)%matrix
1835 END DO
1836 ELSE
1837 cpassert(PRESENT(matrix_p) .AND. nimages == 1)
1838 matrix_images(1)%matrix => matrix_p
1839 END IF
1840
1841 ! Distribute matrix blocks.
1842 IF (any_distributed) THEN
1843 CALL rs_scatter_matrices(matrix_images, task_list%pab_buffer, task_list, group)
1844 ELSE
1845 CALL rs_copy_to_buffer(matrix_images, task_list%pab_buffer, task_list)
1846 END IF
1847 DEALLOCATE (matrix_images)
1848
1849 ! Map all tasks onto the grids
1850 CALL grid_collocate_task_list(task_list=task_list%grid_task_list, &
1851 ga_gb_function=ga_gb_function, &
1852 pab_blocks=task_list%pab_buffer, &
1853 rs_grids=rs_rho)
1854
1855 ! Merge realspace multi-grids into single planewave grid.
1856 CALL density_rs2pw(pw_env, rs_rho, rho, rho_gspace)
1857 IF (PRESENT(total_rho)) total_rho = pw_integrate_function(rho, isign=-1)
1858
1859 CALL timestop(handle)
1860
1861 END SUBROUTINE calculate_rho_elec
1862
1863! **************************************************************************************************
1864!> \brief computes the gradient of the density corresponding to a given
1865!> density matrix on the grid
1866!> \param matrix_p ...
1867!> \param matrix_p_kp ...
1868!> \param drho ...
1869!> \param drho_gspace ...
1870!> \param qs_env ...
1871!> \param soft_valid ...
1872!> \param basis_type ...
1873!> \note this is an alternative to calculate the gradient through FFTs
1874! **************************************************************************************************
1875 SUBROUTINE calculate_drho_elec(matrix_p, matrix_p_kp, drho, drho_gspace, qs_env, &
1876 soft_valid, basis_type)
1877
1878 TYPE(dbcsr_type), OPTIONAL, TARGET :: matrix_p
1879 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
1880 POINTER :: matrix_p_kp
1881 TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT) :: drho
1882 TYPE(pw_c1d_gs_type), DIMENSION(3), INTENT(INOUT) :: drho_gspace
1883 TYPE(qs_environment_type), POINTER :: qs_env
1884 LOGICAL, INTENT(IN), OPTIONAL :: soft_valid
1885 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: basis_type
1886
1887 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_drho_elec'
1888
1889 CHARACTER(LEN=default_string_length) :: my_basis_type
1890 INTEGER :: bcol, brow, dabqadb_func, handle, iatom, iatom_old, idir, igrid_level, ikind, &
1891 ikind_old, img, img_old, ipgf, iset, iset_old, itask, ithread, jatom, jatom_old, jkind, &
1892 jkind_old, jpgf, jset, jset_old, maxco, maxsgf_set, na1, na2, natoms, nb1, nb2, ncoa, &
1893 ncob, nimages, nseta, nsetb, ntasks, nthread, sgfa, sgfb
1894 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
1895 npgfb, nsgfa, nsgfb
1896 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1897 LOGICAL :: atom_pair_changed, distributed_rs_grids, &
1898 do_kp, found, my_soft, use_subpatch
1899 REAL(kind=dp) :: eps_rho_rspace, f, prefactor, radius, &
1900 scale, zetp
1901 REAL(kind=dp), DIMENSION(3) :: ra, rab, rab_inv, rb, rp
1902 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block, pab, sphi_a, sphi_b, work, &
1903 zeta, zetb
1904 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: pabt, workt
1905 TYPE(atom_pair_type), DIMENSION(:), POINTER :: atom_pair_recv, atom_pair_send
1906 TYPE(cell_type), POINTER :: cell
1907 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: deltap
1908 TYPE(dft_control_type), POINTER :: dft_control
1909 TYPE(gridlevel_info_type), POINTER :: gridlevel_info
1910 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1911 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1912 POINTER :: sab_orb
1913 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1914 TYPE(pw_env_type), POINTER :: pw_env
1915 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1916 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
1917 POINTER :: rs_descs
1918 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
1919 TYPE(task_list_type), POINTER :: task_list, task_list_soft
1920 TYPE(task_type), DIMENSION(:), POINTER :: tasks
1921
1922 CALL timeset(routinen, handle)
1923
1924 cpassert(PRESENT(matrix_p) .OR. PRESENT(matrix_p_kp))
1925 do_kp = PRESENT(matrix_p_kp)
1926
1927 NULLIFY (cell, dft_control, orb_basis_set, deltap, qs_kind_set, &
1928 sab_orb, particle_set, rs_rho, pw_env, rs_descs, la_max, la_min, &
1929 lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, p_block, sphi_a, &
1930 sphi_b, zeta, zetb, first_sgfa, first_sgfb, tasks, pabt, workt)
1931
1932 ! by default, the full density is calculated
1933 my_soft = .false.
1934 IF (PRESENT(soft_valid)) my_soft = soft_valid
1935
1936 IF (PRESENT(basis_type)) THEN
1937 my_basis_type = basis_type
1938 ELSE
1939 my_basis_type = "ORB"
1940 END IF
1941
1942 CALL get_qs_env(qs_env=qs_env, &
1943 qs_kind_set=qs_kind_set, &
1944 cell=cell, &
1945 dft_control=dft_control, &
1946 particle_set=particle_set, &
1947 sab_orb=sab_orb, &
1948 pw_env=pw_env)
1949
1950 SELECT CASE (my_basis_type)
1951 CASE ("ORB")
1952 CALL get_qs_env(qs_env=qs_env, &
1953 task_list=task_list, &
1954 task_list_soft=task_list_soft)
1955 CASE ("AUX_FIT")
1956 CALL get_qs_env(qs_env=qs_env, &
1957 task_list_soft=task_list_soft)
1958 CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list)
1959 END SELECT
1960
1961 ! *** assign from pw_env
1962 gridlevel_info => pw_env%gridlevel_info
1963
1964 ! *** Allocate work storage ***
1965 nthread = 1
1966 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
1967 maxco=maxco, &
1968 maxsgf_set=maxsgf_set, &
1969 basis_type=my_basis_type)
1970 CALL reallocate(pabt, 1, maxco, 1, maxco, 0, nthread - 1)
1971 CALL reallocate(workt, 1, maxco, 1, maxsgf_set, 0, nthread - 1)
1972
1973 ! find maximum numbers
1974 nimages = dft_control%nimages
1975 cpassert(nimages == 1 .OR. do_kp)
1976
1977 natoms = SIZE(particle_set)
1978
1979 ! get the task lists
1980 IF (my_soft) task_list => task_list_soft
1981 cpassert(ASSOCIATED(task_list))
1982 tasks => task_list%tasks
1983 atom_pair_send => task_list%atom_pair_send
1984 atom_pair_recv => task_list%atom_pair_recv
1985 ntasks = task_list%ntasks
1986
1987 ! *** set up the rs multi-grids
1988 cpassert(ASSOCIATED(pw_env))
1989 CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho)
1990 DO igrid_level = 1, gridlevel_info%ngrid_levels
1991 distributed_rs_grids = rs_rho(igrid_level)%desc%distributed
1992 END DO
1993
1994 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1995
1996 ! *** Initialize working density matrix ***
1997 ! distributed rs grids require a matrix that will be changed
1998 ! whereas this is not the case for replicated grids
1999 ALLOCATE (deltap(nimages))
2000 IF (distributed_rs_grids) THEN
2001 DO img = 1, nimages
2002 END DO
2003 ! this matrix has no strict sparsity pattern in parallel
2004 ! deltap%sparsity_id=-1
2005 IF (do_kp) THEN
2006 DO img = 1, nimages
2007 CALL dbcsr_copy(deltap(img)%matrix, matrix_p_kp(img)%matrix, &
2008 name="DeltaP")
2009 END DO
2010 ELSE
2011 CALL dbcsr_copy(deltap(1)%matrix, matrix_p, name="DeltaP")
2012 END IF
2013 ELSE
2014 IF (do_kp) THEN
2015 DO img = 1, nimages
2016 deltap(img)%matrix => matrix_p_kp(img)%matrix
2017 END DO
2018 ELSE
2019 deltap(1)%matrix => matrix_p
2020 END IF
2021 END IF
2022
2023 ! distribute the matrix
2024 IF (distributed_rs_grids) THEN
2025 CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltap, &
2026 atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
2027 nimages=nimages, scatter=.true.)
2028 END IF
2029
2030 ! map all tasks on the grids
2031
2032 ithread = 0
2033 pab => pabt(:, :, ithread)
2034 work => workt(:, :, ithread)
2035
2036 loop_xyz: DO idir = 1, 3
2037
2038 DO igrid_level = 1, gridlevel_info%ngrid_levels
2039 CALL rs_grid_zero(rs_rho(igrid_level))
2040 END DO
2041
2042 iatom_old = -1; jatom_old = -1; iset_old = -1; jset_old = -1
2043 ikind_old = -1; jkind_old = -1; img_old = -1
2044 loop_tasks: DO itask = 1, ntasks
2045
2046 !decode the atom pair and basis info
2047 igrid_level = tasks(itask)%grid_level
2048 img = tasks(itask)%image
2049 iatom = tasks(itask)%iatom
2050 jatom = tasks(itask)%jatom
2051 iset = tasks(itask)%iset
2052 jset = tasks(itask)%jset
2053 ipgf = tasks(itask)%ipgf
2054 jpgf = tasks(itask)%jpgf
2055
2056 ikind = particle_set(iatom)%atomic_kind%kind_number
2057 jkind = particle_set(jatom)%atomic_kind%kind_number
2058
2059 IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old .OR. img .NE. img_old) THEN
2060
2061 IF (iatom .NE. iatom_old) ra(:) = pbc(particle_set(iatom)%r, cell)
2062
2063 IF (iatom <= jatom) THEN
2064 brow = iatom
2065 bcol = jatom
2066 ELSE
2067 brow = jatom
2068 bcol = iatom
2069 END IF
2070
2071 IF (ikind .NE. ikind_old) THEN
2072 IF (my_soft) THEN
2073 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
2074 basis_type="ORB_SOFT")
2075 ELSE
2076 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
2077 basis_type=my_basis_type)
2078 END IF
2079 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
2080 first_sgf=first_sgfa, &
2081 lmax=la_max, &
2082 lmin=la_min, &
2083 npgf=npgfa, &
2084 nset=nseta, &
2085 nsgf_set=nsgfa, &
2086 sphi=sphi_a, &
2087 zet=zeta)
2088 END IF
2089
2090 IF (jkind .NE. jkind_old) THEN
2091 IF (my_soft) THEN
2092 CALL get_qs_kind(qs_kind_set(jkind), basis_set=orb_basis_set, &
2093 basis_type="ORB_SOFT")
2094 ELSE
2095 CALL get_qs_kind(qs_kind_set(jkind), basis_set=orb_basis_set, &
2096 basis_type=my_basis_type)
2097 END IF
2098 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
2099 first_sgf=first_sgfb, &
2100 lmax=lb_max, &
2101 lmin=lb_min, &
2102 npgf=npgfb, &
2103 nset=nsetb, &
2104 nsgf_set=nsgfb, &
2105 sphi=sphi_b, &
2106 zet=zetb)
2107 END IF
2108
2109 CALL dbcsr_get_block_p(matrix=deltap(img)%matrix, &
2110 row=brow, col=bcol, block=p_block, found=found)
2111 cpassert(found)
2112
2113 iatom_old = iatom
2114 jatom_old = jatom
2115 ikind_old = ikind
2116 jkind_old = jkind
2117 img_old = img
2118 atom_pair_changed = .true.
2119
2120 ELSE
2121
2122 atom_pair_changed = .false.
2123
2124 END IF
2125
2126 IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset) THEN
2127
2128 ncoa = npgfa(iset)*ncoset(la_max(iset))
2129 sgfa = first_sgfa(1, iset)
2130 ncob = npgfb(jset)*ncoset(lb_max(jset))
2131 sgfb = first_sgfb(1, jset)
2132
2133 IF (iatom <= jatom) THEN
2134 CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
2135 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
2136 p_block(sgfa, sgfb), SIZE(p_block, 1), &
2137 0.0_dp, work(1, 1), maxco)
2138 CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
2139 1.0_dp, work(1, 1), maxco, &
2140 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
2141 0.0_dp, pab(1, 1), maxco)
2142 ELSE
2143 CALL dgemm("N", "N", ncob, nsgfa(iset), nsgfb(jset), &
2144 1.0_dp, sphi_b(1, sgfb), SIZE(sphi_b, 1), &
2145 p_block(sgfb, sgfa), SIZE(p_block, 1), &
2146 0.0_dp, work(1, 1), maxco)
2147 CALL dgemm("N", "T", ncob, ncoa, nsgfa(iset), &
2148 1.0_dp, work(1, 1), maxco, &
2149 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
2150 0.0_dp, pab(1, 1), maxco)
2151 END IF
2152
2153 iset_old = iset
2154 jset_old = jset
2155
2156 END IF
2157
2158 rab(:) = tasks(itask)%rab
2159 rb(:) = ra(:) + rab(:)
2160 zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
2161
2162 f = zetb(jpgf, jset)/zetp
2163 rp(:) = ra(:) + f*rab(:)
2164 prefactor = exp(-zeta(ipgf, iset)*f*dot_product(rab, rab))
2165 radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
2166 lb_min=lb_min(jset), lb_max=lb_max(jset), &
2167 ra=ra, rb=rb, rp=rp, &
2168 zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
2169 prefactor=prefactor, cutoff=1.0_dp)
2170
2171 na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
2172 na2 = ipgf*ncoset(la_max(iset))
2173 nb1 = (jpgf - 1)*ncoset(lb_max(jset)) + 1
2174 nb2 = jpgf*ncoset(lb_max(jset))
2175
2176 ! takes the density matrix symmetry in account, i.e. off-diagonal blocks need to be mapped 'twice'
2177 IF (iatom == jatom .AND. img == 1) THEN
2178 scale = 1.0_dp
2179 ELSE
2180 scale = 2.0_dp
2181 END IF
2182
2183 ! check whether we need to use fawzi's generalised collocation scheme
2184 IF (rs_rho(igrid_level)%desc%distributed) THEN
2185 !tasks(4,:) is 0 for replicated, 1 for distributed 2 for exceptional distributed tasks
2186 IF (tasks(itask)%dist_type .EQ. 2) THEN
2187 use_subpatch = .true.
2188 ELSE
2189 use_subpatch = .false.
2190 END IF
2191 ELSE
2192 use_subpatch = .false.
2193 END IF
2194
2195 SELECT CASE (idir)
2196 CASE (1)
2197 dabqadb_func = grid_func_dabpadb_x
2198 CASE (2)
2199 dabqadb_func = grid_func_dabpadb_y
2200 CASE (3)
2201 dabqadb_func = grid_func_dabpadb_z
2202 CASE DEFAULT
2203 cpabort("invalid idir")
2204 END SELECT
2205
2206 IF (iatom <= jatom) THEN
2207 CALL collocate_pgf_product( &
2208 la_max(iset), zeta(ipgf, iset), la_min(iset), &
2209 lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
2210 ra, rab, scale, pab, na1 - 1, nb1 - 1, &
2211 rs_rho(igrid_level), &
2212 radius=radius, ga_gb_function=dabqadb_func, &
2213 use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern)
2214 ELSE
2215 rab_inv = -rab
2216 CALL collocate_pgf_product( &
2217 lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
2218 la_max(iset), zeta(ipgf, iset), la_min(iset), &
2219 rb, rab_inv, scale, pab, nb1 - 1, na1 - 1, &
2220 rs_rho(igrid_level), &
2221 radius=radius, ga_gb_function=dabqadb_func, &
2222 use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern)
2223 END IF
2224
2225 END DO loop_tasks
2226
2227 CALL density_rs2pw(pw_env, rs_rho, drho(idir), drho_gspace(idir))
2228
2229 END DO loop_xyz
2230
2231 ! *** Release work storage ***
2232 IF (distributed_rs_grids) THEN
2233 CALL dbcsr_deallocate_matrix_set(deltap)
2234 ELSE
2235 DO img = 1, nimages
2236 NULLIFY (deltap(img)%matrix)
2237 END DO
2238 DEALLOCATE (deltap)
2239 END IF
2240
2241 DEALLOCATE (pabt, workt)
2242
2243 CALL timestop(handle)
2244
2245 END SUBROUTINE calculate_drho_elec
2246
2247! **************************************************************************************************
2248!> \brief Computes the gradient wrt. nuclear coordinates of a density on the grid
2249!> The density is given in terms of the density matrix_p
2250!> \param matrix_p Density matrix
2251!> \param matrix_p_kp ...
2252!> \param drho Density gradient on the grid
2253!> \param drho_gspace Density gradient on the reciprocal grid
2254!> \param qs_env ...
2255!> \param soft_valid ...
2256!> \param basis_type ...
2257!> \param beta Derivative direction
2258!> \param lambda Atom index
2259!> \note SL, ED 2021
2260!> Adapted from calculate_drho_elec
2261! **************************************************************************************************
2262 SUBROUTINE calculate_drho_elec_dr(matrix_p, matrix_p_kp, drho, drho_gspace, qs_env, &
2263 soft_valid, basis_type, beta, lambda)
2264
2265 TYPE(dbcsr_type), OPTIONAL, TARGET :: matrix_p
2266 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
2267 POINTER :: matrix_p_kp
2268 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: drho
2269 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: drho_gspace
2270 TYPE(qs_environment_type), POINTER :: qs_env
2271 LOGICAL, INTENT(IN), OPTIONAL :: soft_valid
2272 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: basis_type
2273 INTEGER, INTENT(IN) :: beta, lambda
2274
2275 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_drho_elec_dR'
2276
2277 CHARACTER(LEN=default_string_length) :: my_basis_type
2278 INTEGER :: bcol, brow, dabqadb_func, handle, iatom, iatom_old, igrid_level, ikind, &
2279 ikind_old, img, img_old, ipgf, iset, iset_old, itask, ithread, jatom, jatom_old, jkind, &
2280 jkind_old, jpgf, jset, jset_old, maxco, maxsgf_set, na1, na2, natoms, nb1, nb2, ncoa, &
2281 ncob, nimages, nseta, nsetb, ntasks, nthread, sgfa, sgfb
2282 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
2283 npgfb, nsgfa, nsgfb
2284 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
2285 LOGICAL :: atom_pair_changed, distributed_rs_grids, &
2286 do_kp, found, my_soft, use_subpatch
2287 REAL(kind=dp) :: eps_rho_rspace, f, prefactor, radius, &
2288 scale, zetp
2289 REAL(kind=dp), DIMENSION(3) :: ra, rab, rab_inv, rb, rp
2290 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block, pab, sphi_a, sphi_b, work, &
2291 zeta, zetb
2292 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: pabt, workt
2293 TYPE(atom_pair_type), DIMENSION(:), POINTER :: atom_pair_recv, atom_pair_send
2294 TYPE(cell_type), POINTER :: cell
2295 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: deltap
2296 TYPE(dft_control_type), POINTER :: dft_control
2297 TYPE(gridlevel_info_type), POINTER :: gridlevel_info
2298 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
2299 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2300 TYPE(pw_env_type), POINTER :: pw_env
2301 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2302 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
2303 POINTER :: rs_descs
2304 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
2305 TYPE(task_list_type), POINTER :: task_list, task_list_soft
2306 TYPE(task_type), DIMENSION(:), POINTER :: tasks
2307
2308 CALL timeset(routinen, handle)
2309
2310 cpassert(PRESENT(matrix_p) .OR. PRESENT(matrix_p_kp))
2311 do_kp = PRESENT(matrix_p_kp)
2312
2313 NULLIFY (cell, dft_control, orb_basis_set, deltap, qs_kind_set, &
2314 particle_set, rs_rho, pw_env, rs_descs, la_max, la_min, lb_max, &
2315 lb_min, npgfa, npgfb, nsgfa, nsgfb, p_block, sphi_a, sphi_b, &
2316 zeta, zetb, first_sgfa, first_sgfb, tasks, pabt, workt)
2317
2318 ! by default, the full density is calculated
2319 my_soft = .false.
2320 IF (PRESENT(soft_valid)) my_soft = soft_valid
2321
2322 IF (PRESENT(basis_type)) THEN
2323 my_basis_type = basis_type
2324 ELSE
2325 my_basis_type = "ORB"
2326 END IF
2327
2328 CALL get_qs_env(qs_env=qs_env, &
2329 qs_kind_set=qs_kind_set, &
2330 cell=cell, &
2331 dft_control=dft_control, &
2332 particle_set=particle_set, &
2333 pw_env=pw_env)
2334
2335 SELECT CASE (my_basis_type)
2336 CASE ("ORB")
2337 CALL get_qs_env(qs_env=qs_env, &
2338 task_list=task_list, &
2339 task_list_soft=task_list_soft)
2340 CASE ("AUX_FIT")
2341 CALL get_qs_env(qs_env=qs_env, &
2342 task_list_soft=task_list_soft)
2343 CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list)
2344 END SELECT
2345
2346 ! *** assign from pw_env
2347 gridlevel_info => pw_env%gridlevel_info
2348
2349 ! *** Allocate work storage ***
2350 nthread = 1
2351 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
2352 maxco=maxco, &
2353 maxsgf_set=maxsgf_set, &
2354 basis_type=my_basis_type)
2355 CALL reallocate(pabt, 1, maxco, 1, maxco, 0, nthread - 1)
2356 CALL reallocate(workt, 1, maxco, 1, maxsgf_set, 0, nthread - 1)
2357
2358 ! find maximum numbers
2359 nimages = dft_control%nimages
2360 cpassert(nimages == 1 .OR. do_kp)
2361
2362 natoms = SIZE(particle_set)
2363
2364 ! get the task lists
2365 IF (my_soft) task_list => task_list_soft
2366 cpassert(ASSOCIATED(task_list))
2367 tasks => task_list%tasks
2368 atom_pair_send => task_list%atom_pair_send
2369 atom_pair_recv => task_list%atom_pair_recv
2370 ntasks = task_list%ntasks
2371
2372 ! *** set up the rs multi-grids
2373 cpassert(ASSOCIATED(pw_env))
2374 CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho)
2375 DO igrid_level = 1, gridlevel_info%ngrid_levels
2376 distributed_rs_grids = rs_rho(igrid_level)%desc%distributed
2377 END DO
2378
2379 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
2380
2381 ! *** Initialize working density matrix ***
2382 ! distributed rs grids require a matrix that will be changed
2383 ! whereas this is not the case for replicated grids
2384 ALLOCATE (deltap(nimages))
2385 IF (distributed_rs_grids) THEN
2386 DO img = 1, nimages
2387 END DO
2388 ! this matrix has no strict sparsity pattern in parallel
2389 ! deltap%sparsity_id=-1
2390 IF (do_kp) THEN
2391 DO img = 1, nimages
2392 CALL dbcsr_copy(deltap(img)%matrix, matrix_p_kp(img)%matrix, &
2393 name="DeltaP")
2394 END DO
2395 ELSE
2396 CALL dbcsr_copy(deltap(1)%matrix, matrix_p, name="DeltaP")
2397 END IF
2398 ELSE
2399 IF (do_kp) THEN
2400 DO img = 1, nimages
2401 deltap(img)%matrix => matrix_p_kp(img)%matrix
2402 END DO
2403 ELSE
2404 deltap(1)%matrix => matrix_p
2405 END IF
2406 END IF
2407
2408 ! distribute the matrix
2409 IF (distributed_rs_grids) THEN
2410 CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltap, &
2411 atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
2412 nimages=nimages, scatter=.true.)
2413 END IF
2414
2415 ! map all tasks on the grids
2416
2417 ithread = 0
2418 pab => pabt(:, :, ithread)
2419 work => workt(:, :, ithread)
2420
2421 DO igrid_level = 1, gridlevel_info%ngrid_levels
2422 CALL rs_grid_zero(rs_rho(igrid_level))
2423 END DO
2424
2425 iatom_old = -1; jatom_old = -1; iset_old = -1; jset_old = -1
2426 ikind_old = -1; jkind_old = -1; img_old = -1
2427 loop_tasks: DO itask = 1, ntasks
2428
2429 !decode the atom pair and basis info
2430 igrid_level = tasks(itask)%grid_level
2431 img = tasks(itask)%image
2432 iatom = tasks(itask)%iatom
2433 jatom = tasks(itask)%jatom
2434 iset = tasks(itask)%iset
2435 jset = tasks(itask)%jset
2436 ipgf = tasks(itask)%ipgf
2437 jpgf = tasks(itask)%jpgf
2438
2439 ikind = particle_set(iatom)%atomic_kind%kind_number
2440 jkind = particle_set(jatom)%atomic_kind%kind_number
2441
2442 IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old .OR. img .NE. img_old) THEN
2443
2444 IF (iatom .NE. iatom_old) ra(:) = pbc(particle_set(iatom)%r, cell)
2445
2446 IF (iatom <= jatom) THEN
2447 brow = iatom
2448 bcol = jatom
2449 ELSE
2450 brow = jatom
2451 bcol = iatom
2452 END IF
2453
2454 IF (ikind .NE. ikind_old) THEN
2455 IF (my_soft) THEN
2456 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
2457 basis_type="ORB_SOFT")
2458 ELSE
2459 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
2460 basis_type=my_basis_type)
2461 END IF
2462 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
2463 first_sgf=first_sgfa, &
2464 lmax=la_max, &
2465 lmin=la_min, &
2466 npgf=npgfa, &
2467 nset=nseta, &
2468 nsgf_set=nsgfa, &
2469 sphi=sphi_a, &
2470 zet=zeta)
2471 END IF
2472
2473 IF (jkind .NE. jkind_old) THEN
2474 IF (my_soft) THEN
2475 CALL get_qs_kind(qs_kind_set(jkind), basis_set=orb_basis_set, &
2476 basis_type="ORB_SOFT")
2477 ELSE
2478 CALL get_qs_kind(qs_kind_set(jkind), basis_set=orb_basis_set, &
2479 basis_type=my_basis_type)
2480 END IF
2481 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
2482 first_sgf=first_sgfb, &
2483 lmax=lb_max, &
2484 lmin=lb_min, &
2485 npgf=npgfb, &
2486 nset=nsetb, &
2487 nsgf_set=nsgfb, &
2488 sphi=sphi_b, &
2489 zet=zetb)
2490 END IF
2491
2492 CALL dbcsr_get_block_p(matrix=deltap(img)%matrix, &
2493 row=brow, col=bcol, block=p_block, found=found)
2494 cpassert(found)
2495
2496 iatom_old = iatom
2497 jatom_old = jatom
2498 ikind_old = ikind
2499 jkind_old = jkind
2500 img_old = img
2501 atom_pair_changed = .true.
2502
2503 ELSE
2504
2505 atom_pair_changed = .false.
2506
2507 END IF
2508
2509 IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset) THEN
2510
2511 ncoa = npgfa(iset)*ncoset(la_max(iset))
2512 sgfa = first_sgfa(1, iset)
2513 ncob = npgfb(jset)*ncoset(lb_max(jset))
2514 sgfb = first_sgfb(1, jset)
2515
2516 IF (iatom <= jatom) THEN
2517 CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
2518 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
2519 p_block(sgfa, sgfb), SIZE(p_block, 1), &
2520 0.0_dp, work(1, 1), maxco)
2521 CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
2522 1.0_dp, work(1, 1), maxco, &
2523 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
2524 0.0_dp, pab(1, 1), maxco)
2525 ELSE
2526 CALL dgemm("N", "N", ncob, nsgfa(iset), nsgfb(jset), &
2527 1.0_dp, sphi_b(1, sgfb), SIZE(sphi_b, 1), &
2528 p_block(sgfb, sgfa), SIZE(p_block, 1), &
2529 0.0_dp, work(1, 1), maxco)
2530 CALL dgemm("N", "T", ncob, ncoa, nsgfa(iset), &
2531 1.0_dp, work(1, 1), maxco, &
2532 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
2533 0.0_dp, pab(1, 1), maxco)
2534 END IF
2535
2536 iset_old = iset
2537 jset_old = jset
2538
2539 END IF
2540
2541 rab(:) = tasks(itask)%rab
2542 rb(:) = ra(:) + rab(:)
2543 zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
2544
2545 f = zetb(jpgf, jset)/zetp
2546 rp(:) = ra(:) + f*rab(:)
2547 prefactor = exp(-zeta(ipgf, iset)*f*dot_product(rab, rab))
2548 radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
2549 lb_min=lb_min(jset), lb_max=lb_max(jset), &
2550 ra=ra, rb=rb, rp=rp, &
2551 zetp=zetp, eps=eps_rho_rspace, &
2552 prefactor=prefactor, cutoff=1.0_dp)
2553
2554 na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
2555 na2 = ipgf*ncoset(la_max(iset))
2556 nb1 = (jpgf - 1)*ncoset(lb_max(jset)) + 1
2557 nb2 = jpgf*ncoset(lb_max(jset))
2558
2559 ! takes the density matrix symmetry in account, i.e. off-diagonal blocks need to be mapped 'twice'
2560 IF (iatom == jatom .AND. img == 1) THEN
2561 scale = 1.0_dp
2562 ELSE
2563 scale = 2.0_dp
2564 END IF
2565
2566 ! check whether we need to use fawzi's generalised collocation scheme
2567 IF (rs_rho(igrid_level)%desc%distributed) THEN
2568 !tasks(4,:) is 0 for replicated, 1 for distributed 2 for exceptional distributed tasks
2569 IF (tasks(itask)%dist_type .EQ. 2) THEN
2570 use_subpatch = .true.
2571 ELSE
2572 use_subpatch = .false.
2573 END IF
2574 ELSE
2575 use_subpatch = .false.
2576 END IF
2577
2578 SELECT CASE (beta)
2579 CASE (1)
2580 dabqadb_func = grid_func_dab_x
2581 CASE (2)
2582 dabqadb_func = grid_func_dab_y
2583 CASE (3)
2584 dabqadb_func = grid_func_dab_z
2585 CASE DEFAULT
2586 cpabort("invalid beta")
2587 END SELECT
2588
2589 IF (iatom <= jatom) THEN
2590 IF (iatom == lambda) &
2591 CALL collocate_pgf_product( &
2592 la_max(iset), zeta(ipgf, iset), la_min(iset), &
2593 lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
2594 ra, rab, scale, pab, na1 - 1, nb1 - 1, &
2595 rsgrid=rs_rho(igrid_level), &
2596 ga_gb_function=dabqadb_func, radius=radius, &
2597 use_subpatch=use_subpatch, &
2598 subpatch_pattern=tasks(itask)%subpatch_pattern)
2599 IF (jatom == lambda) &
2600 CALL collocate_pgf_product( &
2601 la_max(iset), zeta(ipgf, iset), la_min(iset), &
2602 lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
2603 ra, rab, scale, pab, na1 - 1, nb1 - 1, &
2604 rsgrid=rs_rho(igrid_level), &
2605 ga_gb_function=dabqadb_func + 3, radius=radius, &
2606 use_subpatch=use_subpatch, &
2607 subpatch_pattern=tasks(itask)%subpatch_pattern)
2608 ELSE
2609 rab_inv = -rab
2610 IF (jatom == lambda) &
2611 CALL collocate_pgf_product( &
2612 lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
2613 la_max(iset), zeta(ipgf, iset), la_min(iset), &
2614 rb, rab_inv, scale, pab, nb1 - 1, na1 - 1, &
2615 rs_rho(igrid_level), &
2616 ga_gb_function=dabqadb_func, radius=radius, &
2617 use_subpatch=use_subpatch, &
2618 subpatch_pattern=tasks(itask)%subpatch_pattern)
2619 IF (iatom == lambda) &
2620 CALL collocate_pgf_product( &
2621 lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
2622 la_max(iset), zeta(ipgf, iset), la_min(iset), &
2623 rb, rab_inv, scale, pab, nb1 - 1, na1 - 1, &
2624 rs_rho(igrid_level), &
2625 ga_gb_function=dabqadb_func + 3, radius=radius, &
2626 use_subpatch=use_subpatch, &
2627 subpatch_pattern=tasks(itask)%subpatch_pattern)
2628 END IF
2629
2630 END DO loop_tasks
2631
2632 CALL density_rs2pw(pw_env, rs_rho, drho, drho_gspace)
2633
2634 ! *** Release work storage ***
2635 IF (distributed_rs_grids) THEN
2636 CALL dbcsr_deallocate_matrix_set(deltap)
2637 ELSE
2638 DO img = 1, nimages
2639 NULLIFY (deltap(img)%matrix)
2640 END DO
2641 DEALLOCATE (deltap)
2642 END IF
2643
2644 DEALLOCATE (pabt, workt)
2645
2646 CALL timestop(handle)
2647
2648 END SUBROUTINE calculate_drho_elec_dr
2649
2650! **************************************************************************************************
2651!> \brief maps a single gaussian on the grid
2652!> \param rho ...
2653!> \param rho_gspace ...
2654!> \param atomic_kind_set ...
2655!> \param qs_kind_set ...
2656!> \param cell ...
2657!> \param dft_control ...
2658!> \param particle_set ...
2659!> \param pw_env ...
2660!> \param required_function ...
2661!> \param basis_type ...
2662!> \par History
2663!> 08.2022 created from calculate_wavefunction
2664!> \note
2665!> modified calculate_wave function assuming that the collocation of only a single Gaussian is required.
2666!> chooses a basis function (in contrast to calculate_rho_core or calculate_rho_single_gaussian)
2667! **************************************************************************************************
2668 SUBROUTINE collocate_single_gaussian(rho, rho_gspace, &
2669 atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
2670 pw_env, required_function, basis_type)
2671
2672 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho
2673 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_gspace
2674 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2675 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2676 TYPE(cell_type), POINTER :: cell
2677 TYPE(dft_control_type), POINTER :: dft_control
2678 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2679 TYPE(pw_env_type), POINTER :: pw_env
2680 INTEGER, INTENT(IN) :: required_function
2681 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: basis_type
2682
2683 CHARACTER(LEN=*), PARAMETER :: routinen = 'collocate_single_gaussian'
2684
2685 CHARACTER(LEN=default_string_length) :: my_basis_type
2686 INTEGER :: group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, maxco, maxsgf_set, &
2687 my_index, my_pos, na1, na2, natom, ncoa, nseta, offset, sgfa
2688 INTEGER, ALLOCATABLE, DIMENSION(:) :: where_is_the_point
2689 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
2690 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
2691 LOGICAL :: found
2692 REAL(kind=dp) :: dab, eps_rho_rspace, radius, scale
2693 REAL(kind=dp), DIMENSION(3) :: ra
2694 REAL(kind=dp), DIMENSION(:, :), POINTER :: pab, sphi_a, zeta
2695 TYPE(gridlevel_info_type), POINTER :: gridlevel_info
2696 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
2697 TYPE(mp_comm_type) :: group
2698 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
2699 TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_gspace
2700 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_rspace
2701 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
2702
2703 IF (PRESENT(basis_type)) THEN
2704 my_basis_type = basis_type
2705 ELSE
2706 my_basis_type = "ORB"
2707 END IF
2708
2709 CALL timeset(routinen, handle)
2710
2711 NULLIFY (orb_basis_set, pab, la_max, la_min, npgfa, nsgfa, sphi_a, &
2712 zeta, first_sgfa, rs_rho, pw_pools)
2713
2714 ! *** set up the pw multi-grids
2715 cpassert(ASSOCIATED(pw_env))
2716 CALL pw_env_get(pw_env, rs_grids=rs_rho, pw_pools=pw_pools, &
2717 gridlevel_info=gridlevel_info)
2718
2719 CALL pw_pools_create_pws(pw_pools, mgrid_gspace)
2720 CALL pw_pools_create_pws(pw_pools, mgrid_rspace)
2721
2722 ! *** set up rs multi-grids
2723 DO igrid_level = 1, gridlevel_info%ngrid_levels
2724 CALL rs_grid_zero(rs_rho(igrid_level))
2725 END DO
2726
2727 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
2728! *** Allocate work storage ***
2729 CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
2730 CALL get_qs_kind_set(qs_kind_set, &
2731 maxco=maxco, &
2732 maxsgf_set=maxsgf_set, &
2733 basis_type=my_basis_type)
2734
2735 ALLOCATE (pab(maxco, 1))
2736
2737 offset = 0
2738 group = mgrid_rspace(1)%pw_grid%para%group
2739 my_pos = mgrid_rspace(1)%pw_grid%para%group%mepos
2740 group_size = mgrid_rspace(1)%pw_grid%para%group%num_pe
2741 ALLOCATE (where_is_the_point(0:group_size - 1))
2742
2743 DO iatom = 1, natom
2744 ikind = particle_set(iatom)%atomic_kind%kind_number
2745 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=my_basis_type)
2746 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
2747 first_sgf=first_sgfa, &
2748 lmax=la_max, &
2749 lmin=la_min, &
2750 npgf=npgfa, &
2751 nset=nseta, &
2752 nsgf_set=nsgfa, &
2753 sphi=sphi_a, &
2754 zet=zeta)
2755 ra(:) = pbc(particle_set(iatom)%r, cell)
2756 dab = 0.0_dp
2757
2758 DO iset = 1, nseta
2759
2760 ncoa = npgfa(iset)*ncoset(la_max(iset))
2761 sgfa = first_sgfa(1, iset)
2762
2763 found = .false.
2764 my_index = 0
2765 DO i = 1, nsgfa(iset)
2766 IF (offset + i == required_function) THEN
2767 my_index = i
2768 found = .true.
2769 EXIT
2770 END IF
2771 END DO
2772
2773 IF (found) THEN
2774
2775 pab(1:ncoa, 1) = sphi_a(1:ncoa, sgfa + my_index - 1)
2776
2777 DO ipgf = 1, npgfa(iset)
2778
2779 na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
2780 na2 = ipgf*ncoset(la_max(iset))
2781
2782 scale = 1.0_dp
2783 igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
2784
2785 IF (map_gaussian_here(rs_rho(igrid_level), cell%h_inv, ra, offset, group_size, my_pos)) THEN
2786 radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
2787 lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
2788 zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
2789 prefactor=1.0_dp, cutoff=1.0_dp)
2790
2791 CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), la_min(iset), &
2792 0, 0.0_dp, 0, &
2793 ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
2794 scale, pab, na1 - 1, 0, rs_rho(igrid_level), &
2795 radius=radius, ga_gb_function=grid_func_ab)
2796 END IF
2797
2798 END DO
2799
2800 END IF
2801
2802 offset = offset + nsgfa(iset)
2803
2804 END DO
2805
2806 END DO
2807
2808 DO igrid_level = 1, gridlevel_info%ngrid_levels
2809 CALL transfer_rs2pw(rs_rho(igrid_level), &
2810 mgrid_rspace(igrid_level))
2811 END DO
2812
2813 CALL pw_zero(rho_gspace)
2814 DO igrid_level = 1, gridlevel_info%ngrid_levels
2815 CALL pw_transfer(mgrid_rspace(igrid_level), &
2816 mgrid_gspace(igrid_level))
2817 CALL pw_axpy(mgrid_gspace(igrid_level), rho_gspace)
2818 END DO
2819
2820 CALL pw_transfer(rho_gspace, rho)
2821
2822 ! Release work storage
2823 DEALLOCATE (pab)
2824
2825 ! give back the pw multi-grids
2826 CALL pw_pools_give_back_pws(pw_pools, mgrid_gspace)
2827 CALL pw_pools_give_back_pws(pw_pools, mgrid_rspace)
2828
2829 CALL timestop(handle)
2830
2831 END SUBROUTINE collocate_single_gaussian
2832
2833! **************************************************************************************************
2834!> \brief maps a given wavefunction on the grid
2835!> \param mo_vectors ...
2836!> \param ivector ...
2837!> \param rho ...
2838!> \param rho_gspace ...
2839!> \param atomic_kind_set ...
2840!> \param qs_kind_set ...
2841!> \param cell ...
2842!> \param dft_control ...
2843!> \param particle_set ...
2844!> \param pw_env ...
2845!> \param basis_type ...
2846!> \par History
2847!> 08.2002 created [Joost VandeVondele]
2848!> 03.2006 made independent of qs_env [Joost VandeVondele]
2849!> 08.2024 call collocate_function [JGH]
2850! **************************************************************************************************
2851 SUBROUTINE calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, &
2852 atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
2853 pw_env, basis_type)
2854 TYPE(cp_fm_type), INTENT(IN) :: mo_vectors
2855 INTEGER, INTENT(IN) :: ivector
2856 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho
2857 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_gspace
2858 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2859 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2860 TYPE(cell_type), POINTER :: cell
2861 TYPE(dft_control_type), POINTER :: dft_control
2862 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2863 TYPE(pw_env_type), POINTER :: pw_env
2864 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: basis_type
2865
2866 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_wavefunction'
2867
2868 INTEGER :: handle, i, nao
2869 LOGICAL :: local
2870 REAL(kind=dp) :: eps_rho_rspace
2871 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvector
2872
2873 CALL timeset(routinen, handle)
2874
2875 CALL cp_fm_get_info(matrix=mo_vectors, nrow_global=nao)
2876 ALLOCATE (eigenvector(nao))
2877 DO i = 1, nao
2878 CALL cp_fm_get_element(mo_vectors, i, ivector, eigenvector(i), local)
2879 END DO
2880
2881 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
2882
2883 CALL collocate_function(eigenvector, rho, rho_gspace, &
2884 atomic_kind_set, qs_kind_set, cell, particle_set, pw_env, &
2885 eps_rho_rspace, basis_type)
2886
2887 DEALLOCATE (eigenvector)
2888
2889 CALL timestop(handle)
2890
2891 END SUBROUTINE calculate_wavefunction
2892
2893! **************************************************************************************************
2894!> \brief maps a given function on the grid
2895!> \param vector ...
2896!> \param rho ...
2897!> \param rho_gspace ...
2898!> \param atomic_kind_set ...
2899!> \param qs_kind_set ...
2900!> \param cell ...
2901!> \param particle_set ...
2902!> \param pw_env ...
2903!> \param eps_rho_rspace ...
2904!> \param basis_type ...
2905!> \par History
2906!> 08.2002 created [Joost VandeVondele]
2907!> 03.2006 made independent of qs_env [Joost VandeVondele]
2908!> 08.2024 specialized version from calculate_wavefunction [JGH]
2909!> \notes
2910!> modified calculate_rho_elec, should write the wavefunction represented by vector
2911!> it's presumably dominated by the FFT and the rs->pw and back routines
2912! **************************************************************************************************
2913 SUBROUTINE collocate_function(vector, rho, rho_gspace, &
2914 atomic_kind_set, qs_kind_set, cell, particle_set, pw_env, &
2915 eps_rho_rspace, basis_type)
2916 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: vector
2917 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho
2918 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_gspace
2919 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2920 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2921 TYPE(cell_type), POINTER :: cell
2922 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2923 TYPE(pw_env_type), POINTER :: pw_env
2924 REAL(kind=dp), INTENT(IN) :: eps_rho_rspace
2925 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: basis_type
2926
2927 CHARACTER(LEN=*), PARAMETER :: routinen = 'collocate_function'
2928
2929 CHARACTER(LEN=default_string_length) :: my_basis_type
2930 INTEGER :: group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, maxco, maxsgf_set, &
2931 my_pos, na1, na2, natom, ncoa, nseta, offset, sgfa
2932 INTEGER, ALLOCATABLE, DIMENSION(:) :: where_is_the_point
2933 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
2934 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
2935 REAL(kind=dp) :: dab, radius, scale
2936 REAL(kind=dp), DIMENSION(3) :: ra
2937 REAL(kind=dp), DIMENSION(:, :), POINTER :: pab, sphi_a, work, zeta
2938 TYPE(gridlevel_info_type), POINTER :: gridlevel_info
2939 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
2940 TYPE(mp_comm_type) :: group
2941 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
2942 TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_gspace
2943 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_rspace
2944 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
2945
2946 CALL timeset(routinen, handle)
2947
2948 IF (PRESENT(basis_type)) THEN
2949 my_basis_type = basis_type
2950 ELSE
2951 my_basis_type = "ORB"
2952 END IF
2953
2954 NULLIFY (orb_basis_set, pab, work, la_max, la_min, &
2955 npgfa, nsgfa, sphi_a, zeta, first_sgfa, rs_rho, pw_pools)
2956
2957 ! *** set up the pw multi-grids
2958 cpassert(ASSOCIATED(pw_env))
2959 CALL pw_env_get(pw_env, rs_grids=rs_rho, pw_pools=pw_pools, &
2960 gridlevel_info=gridlevel_info)
2961
2962 CALL pw_pools_create_pws(pw_pools, mgrid_gspace)
2963 CALL pw_pools_create_pws(pw_pools, mgrid_rspace)
2964
2965 ! *** set up rs multi-grids
2966 DO igrid_level = 1, gridlevel_info%ngrid_levels
2967 CALL rs_grid_zero(rs_rho(igrid_level))
2968 END DO
2969
2970 ! *** Allocate work storage ***
2971 CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
2972 CALL get_qs_kind_set(qs_kind_set, &
2973 maxco=maxco, &
2974 maxsgf_set=maxsgf_set, &
2975 basis_type=my_basis_type)
2976
2977 ALLOCATE (pab(maxco, 1))
2978 ALLOCATE (work(maxco, 1))
2979
2980 offset = 0
2981 group = mgrid_rspace(1)%pw_grid%para%group
2982 my_pos = mgrid_rspace(1)%pw_grid%para%group%mepos
2983 group_size = mgrid_rspace(1)%pw_grid%para%group%num_pe
2984 ALLOCATE (where_is_the_point(0:group_size - 1))
2985
2986 DO iatom = 1, natom
2987 ikind = particle_set(iatom)%atomic_kind%kind_number
2988 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=my_basis_type)
2989 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
2990 first_sgf=first_sgfa, &
2991 lmax=la_max, &
2992 lmin=la_min, &
2993 npgf=npgfa, &
2994 nset=nseta, &
2995 nsgf_set=nsgfa, &
2996 sphi=sphi_a, &
2997 zet=zeta)
2998 ra(:) = pbc(particle_set(iatom)%r, cell)
2999 dab = 0.0_dp
3000
3001 DO iset = 1, nseta
3002
3003 ncoa = npgfa(iset)*ncoset(la_max(iset))
3004 sgfa = first_sgfa(1, iset)
3005
3006 DO i = 1, nsgfa(iset)
3007 work(i, 1) = vector(offset + i)
3008 END DO
3009
3010 CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), &
3011 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
3012 work(1, 1), SIZE(work, 1), &
3013 0.0_dp, pab(1, 1), SIZE(pab, 1))
3014
3015 DO ipgf = 1, npgfa(iset)
3016
3017 na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
3018 na2 = ipgf*ncoset(la_max(iset))
3019
3020 scale = 1.0_dp
3021 igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
3022
3023 IF (map_gaussian_here(rs_rho(igrid_level), cell%h_inv, ra, offset, group_size, my_pos)) THEN
3024 radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
3025 lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
3026 zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
3027 prefactor=1.0_dp, cutoff=1.0_dp)
3028
3029 CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), la_min(iset), &
3030 0, 0.0_dp, 0, &
3031 ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
3032 scale, pab, na1 - 1, 0, rs_rho(igrid_level), &
3033 radius=radius, ga_gb_function=grid_func_ab)
3034 END IF
3035
3036 END DO
3037
3038 offset = offset + nsgfa(iset)
3039
3040 END DO
3041
3042 END DO
3043
3044 DO igrid_level = 1, gridlevel_info%ngrid_levels
3045 CALL transfer_rs2pw(rs_rho(igrid_level), &
3046 mgrid_rspace(igrid_level))
3047 END DO
3048
3049 CALL pw_zero(rho_gspace)
3050 DO igrid_level = 1, gridlevel_info%ngrid_levels
3051 CALL pw_transfer(mgrid_rspace(igrid_level), &
3052 mgrid_gspace(igrid_level))
3053 CALL pw_axpy(mgrid_gspace(igrid_level), rho_gspace)
3054 END DO
3055
3056 CALL pw_transfer(rho_gspace, rho)
3057
3058 ! Release work storage
3059 DEALLOCATE (pab)
3060 DEALLOCATE (work)
3061
3062 ! give back the pw multi-grids
3063 CALL pw_pools_give_back_pws(pw_pools, mgrid_gspace)
3064 CALL pw_pools_give_back_pws(pw_pools, mgrid_rspace)
3065
3066 CALL timestop(handle)
3067
3068 END SUBROUTINE collocate_function
3069
3070END MODULE qs_collocate_density
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:120
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:206
cell_types::pbc
Definition cell_types.F:92
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
external_potential_types::get_potential
Definition external_potential_types.F:276
memory_utilities::reallocate
Definition memory_utilities.F:27
pw_methods::pw_axpy
Definition pw_methods.F:164
pw_methods::pw_integrate_function
Definition pw_methods.F:115
pw_methods::pw_transfer
Definition pw_methods.F:303
pw_methods::pw_zero
Definition pw_methods.F:82
pw_pool_types::pw_pools_create_pws
Definition pw_pool_types.F:139
pw_pool_types::pw_pools_give_back_pws
Definition pw_pool_types.F:150
qs_collocate_density::calculate_rho_core
Definition qs_collocate_density.F:138
qs_collocate_density::calculate_rho_resp_all
Definition qs_collocate_density.F:143
admm_types
Types and set/get functions for auxiliary density matrix methods.
Definition admm_types.F:15
admm_types::get_admm_env
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition admm_types.F:593
ao_util
All kind of helpful little routines.
Definition ao_util.F:14
ao_util::exp_radius_very_extended
real(kind=dp) function, public exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, zetp, eps, prefactor, cutoff, epsabs)
computes the radius of the Gaussian outside of which it is smaller than eps
Definition ao_util.F:209
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
Definition basis_set_types.F:624
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:679
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1009
cp_fm_types::cp_fm_get_element
subroutine, public cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
returns an element of a fm this value is valid on every cpu using this call is expensive
Definition cp_fm_types.F:636
external_potential_types
Definition of the atomic potential types.
Definition external_potential_types.F:14
gaussian_gridlevels
Definition gaussian_gridlevels.F:13
gaussian_gridlevels::gaussian_gridlevel
integer function, public gaussian_gridlevel(gridlevel_info, exponent)
...
Definition gaussian_gridlevels.F:148
grid_api
Fortran API for the grid package, which is written in C.
Definition grid_api.F:12
grid_api::grid_func_core_x
integer, parameter, public grid_func_core_x
Definition grid_api.F:62
grid_api::grid_func_dab_z
integer, parameter, public grid_func_dab_z
Definition grid_api.F:57
grid_api::grid_collocate_task_list
subroutine, public grid_collocate_task_list(task_list, ga_gb_function, pab_blocks, rs_grids)
Collocate all tasks of in given list onto given grids.
Definition grid_api.F:968
grid_api::grid_func_dzdx
integer, parameter, public grid_func_dzdx
Definition grid_api.F:51
grid_api::grid_func_dzdz
integer, parameter, public grid_func_dzdz
Definition grid_api.F:54
grid_api::grid_func_dydz
integer, parameter, public grid_func_dydz
Definition grid_api.F:50
grid_api::grid_func_dxdy
integer, parameter, public grid_func_dxdy
Definition grid_api.F:49
grid_api::grid_func_dabpadb_y
integer, parameter, public grid_func_dabpadb_y
Definition grid_api.F:44
grid_api::grid_func_dab_y
integer, parameter, public grid_func_dab_y
Definition grid_api.F:56
grid_api::grid_func_dxdx
integer, parameter, public grid_func_dxdx
Definition grid_api.F:52
grid_api::grid_func_dadb
integer, parameter, public grid_func_dadb
Definition grid_api.F:30
grid_api::grid_func_dydy
integer, parameter, public grid_func_dydy
Definition grid_api.F:53
grid_api::grid_func_dabpadb_z
integer, parameter, public grid_func_dabpadb_z
Definition grid_api.F:45
grid_api::grid_func_dabpadb_x
integer, parameter, public grid_func_dabpadb_x
Definition grid_api.F:43
grid_api::grid_func_dx
integer, parameter, public grid_func_dx
Definition grid_api.F:46
grid_api::grid_func_dz
integer, parameter, public grid_func_dz
Definition grid_api.F:48
grid_api::grid_func_ab
integer, parameter, public grid_func_ab
Definition grid_api.F:29
grid_api::grid_func_core_y
integer, parameter, public grid_func_core_y
Definition grid_api.F:63
grid_api::grid_func_dab_x
integer, parameter, public grid_func_dab_x
Definition grid_api.F:55
grid_api::collocate_pgf_product
subroutine, public collocate_pgf_product(la_max, zeta, la_min, lb_max, zetb, lb_min, ra, rab, scale, pab, o1, o2, rsgrid, ga_gb_function, radius, use_subpatch, subpatch_pattern)
low level collocation of primitive gaussian functions
Definition grid_api.F:119
grid_api::grid_func_core_z
integer, parameter, public grid_func_core_z
Definition grid_api.F:64
grid_api::grid_func_dy
integer, parameter, public grid_func_dy
Definition grid_api.F:47
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::orb_dxy
integer, parameter, public orb_dxy
Definition input_constants.F:486
input_constants::orb_pz
integer, parameter, public orb_pz
Definition input_constants.F:486
input_constants::orb_dz2
integer, parameter, public orb_dz2
Definition input_constants.F:486
input_constants::orb_s
integer, parameter, public orb_s
Definition input_constants.F:486
input_constants::orb_py
integer, parameter, public orb_py
Definition input_constants.F:486
input_constants::orb_dyz
integer, parameter, public orb_dyz
Definition input_constants.F:486
input_constants::orb_px
integer, parameter, public orb_px
Definition input_constants.F:486
input_constants::orb_dzx
integer, parameter, public orb_dzx
Definition input_constants.F:486
input_constants::orb_dy2
integer, parameter, public orb_dy2
Definition input_constants.F:486
input_constants::orb_dx2
integer, parameter, public orb_dx2
Definition input_constants.F:486
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
lri_environment_types
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Definition lri_environment_types.F:18
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition orbital_pointers.F:45
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition qs_collocate_density.F:38
qs_collocate_density::calculate_drho_elec_dr
subroutine, public calculate_drho_elec_dr(matrix_p, matrix_p_kp, drho, drho_gspace, qs_env, soft_valid, basis_type, beta, lambda)
Computes the gradient wrt. nuclear coordinates of a density on the grid The density is given in terms...
Definition qs_collocate_density.F:2264
qs_collocate_density::calculate_drho_elec
subroutine, public calculate_drho_elec(matrix_p, matrix_p_kp, drho, drho_gspace, qs_env, soft_valid, basis_type)
computes the gradient of the density corresponding to a given density matrix on the grid
Definition qs_collocate_density.F:1877
qs_collocate_density::calculate_rho_elec
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
Definition qs_collocate_density.F:1711
qs_collocate_density::collocate_function
subroutine, public collocate_function(vector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, particle_set, pw_env, eps_rho_rspace, basis_type)
maps a given function on the grid
Definition qs_collocate_density.F:2916
qs_collocate_density::calculate_rho_metal
subroutine, public calculate_rho_metal(rho_metal, coeff, total_rho_metal, qs_env)
computes the image charge density on the grid (including coeffcients)
Definition qs_collocate_density.F:1293
qs_collocate_density::calculate_rho_resp_single
subroutine, public calculate_rho_resp_single(rho_gb, qs_env, eta, iatom_in)
collocate a single Gaussian on the grid for periodic RESP fitting
Definition qs_collocate_density.F:1399
qs_collocate_density::calculate_rho_nlcc
subroutine, public calculate_rho_nlcc(rho_nlcc, qs_env)
computes the density of the non-linear core correction on the grid
Definition qs_collocate_density.F:155
qs_collocate_density::calculate_lri_rho_elec
subroutine, public calculate_lri_rho_elec(lri_rho_g, lri_rho_r, qs_env, lri_coef, total_rho, basis_type, exact_1c_terms, pmat, atomlist)
Collocates the fitted lri density on a grid.
Definition qs_collocate_density.F:518
qs_collocate_density::collocate_single_gaussian
subroutine, public collocate_single_gaussian(rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, required_function, basis_type)
maps a single gaussian on the grid
Definition qs_collocate_density.F:2671
qs_collocate_density::calculate_rho_single_gaussian
subroutine, public calculate_rho_single_gaussian(rho_gb, qs_env, iatom_in)
collocate a single Gaussian on the grid
Definition qs_collocate_density.F:1203
qs_collocate_density::calculate_ppl_grid
subroutine, public calculate_ppl_grid(vppl, qs_env)
computes the local pseudopotential (without erf term) on the grid
Definition qs_collocate_density.F:334
qs_collocate_density::calculate_wavefunction
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type)
maps a given wavefunction on the grid
Definition qs_collocate_density.F:2854
qs_collocate_density::calculate_drho_core
subroutine, public calculate_drho_core(drho_core, qs_env, beta, lambda)
Computes the derivative of the density of the core charges with respect to the nuclear coordinates on...
Definition qs_collocate_density.F:1071
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:454
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:880
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition qs_ks_types.F:333
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
realspace_grid_types
Definition realspace_grid_types.F:18
realspace_grid_types::map_gaussian_here
pure logical function, public map_gaussian_here(rs_grid, h_inv, ra, offset, group_size, my_pos)
...
Definition realspace_grid_types.F:2144
realspace_grid_types::transfer_rs2pw
subroutine, public transfer_rs2pw(rs, pw)
...
Definition realspace_grid_types.F:652
realspace_grid_types::rs_grid_zero
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
Definition realspace_grid_types.F:1905
rs_pw_interface
Transfers densities from PW to RS grids and potentials from PW to RS.
Definition rs_pw_interface.F:14
rs_pw_interface::density_rs2pw
subroutine, public density_rs2pw(pw_env, rs_rho, rho, rho_gspace)
given partial densities on the realspace multigrids, computes the full density on the plane wave grid...
Definition rs_pw_interface.F:65
task_list_methods
generate the tasks lists used by collocate and integrate routines
Definition task_list_methods.F:14
task_list_methods::rs_scatter_matrices
subroutine, public rs_scatter_matrices(src_matrices, dest_buffer, task_list, group)
Scatters dbcsr matrix blocks and receives them into a buffer as needed before collocation.
Definition task_list_methods.F:2525
task_list_methods::rs_distribute_matrix
subroutine, public rs_distribute_matrix(rs_descs, pmats, atom_pair_send, atom_pair_recv, nimages, scatter, hmats)
redistributes the matrix so that it can be used in realspace operations i.e. according to the task li...
Definition task_list_methods.F:2187
task_list_methods::rs_copy_to_buffer
subroutine, public rs_copy_to_buffer(src_matrices, dest_buffer, task_list)
Copies the DBCSR blocks into buffer, replaces rs_scatter_matrix for non-distributed grids.
Definition task_list_methods.F:2596
task_list_types
types for task lists
Definition task_list_types.F:14
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
external_potential_types::gth_potential_type
Definition external_potential_types.F:116
gaussian_gridlevels::gridlevel_info_type
Definition gaussian_gridlevels.F:34
lri_environment_types::lri_kind_type
Definition lri_environment_types.F:332
message_passing::mp_comm_type
Definition message_passing.F:189
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_p_type
to create arrays of pools
Definition pw_pool_types.F:135
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:136
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
realspace_grid_types::realspace_grid_desc_p_type
Definition realspace_grid_types.F:158
realspace_grid_types::realspace_grid_type
Definition realspace_grid_types.F:136
task_list_types::atom_pair_type
Definition task_list_types.F:51
task_list_types::task_list_type
Definition task_list_types.F:58
task_list_types::task_type
Definition task_list_types.F:31