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