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