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