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