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