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