(git:f2099e5)
Loading...
Searching...
No Matches
kpoint_methods.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 Routines needed for kpoint calculation
10!> \par History
11!> 2014.07 created [JGH]
12!> 2014.11 unified k-point and gamma-point code [Ole Schuett]
13!> \author JGH
14! **************************************************************************************************
15MODULE kpoint_methods
16 USE atomic_kind_types, ONLY: get_atomic_kind
17 USE cell_types, ONLY: cell_type,&
18 pbc_stable,&
19 real_to_scaled
20 USE cp_blacs_env, ONLY: cp_blacs_env_create,&
21 cp_blacs_env_type
22 USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,&
23 cp_cfm_scale_and_add,&
24 cp_cfm_transpose
25 USE cp_cfm_types, ONLY: cp_cfm_create,&
26 cp_cfm_get_info,&
27 cp_cfm_get_submatrix,&
28 cp_cfm_release,&
29 cp_cfm_to_cfm,&
30 cp_cfm_to_fm,&
31 cp_cfm_type,&
32 cp_fm_to_cfm
33 USE cp_control_types, ONLY: hairy_probes_type
34 USE cp_dbcsr_api, ONLY: &
35 dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribute, &
36 dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_get_block_p, dbcsr_get_info, &
37 dbcsr_get_readonly_block_p, dbcsr_get_stored_coordinates, dbcsr_iterator_blocks_left, &
38 dbcsr_iterator_next_block, dbcsr_iterator_readonly_start, dbcsr_iterator_start, &
39 dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_replicate_all, dbcsr_set, &
40 dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
41 USE cp_dbcsr_contrib, ONLY: dbcsr_dot
42 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
43 USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr
44 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
45 USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
46 fm_pool_create_fm,&
47 fm_pool_give_back_fm
48 USE cp_fm_struct, ONLY: cp_fm_struct_equivalent,&
49 cp_fm_struct_type
50 USE cp_fm_types, ONLY: &
51 copy_info_type, cp_fm_cleanup_copy_general, cp_fm_create, cp_fm_finish_copy_general, &
52 cp_fm_get_diag, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_p_type, cp_fm_release, &
53 cp_fm_start_copy_general, cp_fm_to_fm, cp_fm_type
54 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
55 USE cryssym, ONLY: crys_sym_gen,&
56 csym_type,&
57 kpoint_gen,&
58 kpoint_gen_general,&
59 print_crys_symmetry,&
60 print_kp_symmetry,&
61 release_csym_type
62 USE hairy_probes, ONLY: probe_occupancy_kp
63 USE input_constants, ONLY: smear_fermi_dirac,&
64 smear_gaussian,&
65 smear_mp,&
66 smear_mv
67 USE input_cp2k_kpoints, ONLY: lattice_fft_auto,&
68 lattice_fft_off,&
69 use_spglib_kpoint_backend,&
70 use_spglib_kpoint_symmetry
71 USE kinds, ONLY: dp,&
72 int_8
73 USE kpoint_lattice_fft, ONLY: cell_to_k_grid_fft,&
74 k_grid_to_cell_execute,&
75 k_grid_to_cell_prepare,&
76 k_grid_to_cell_release,&
77 k_grid_to_cell_work_type,&
78 lattice_fft_shape,&
79 regular_kpoint_grid
80 USE kpoint_types, ONLY: &
81 find_kpoint_rotation_slot, get_kpoint_info, kind_rotmat_type, kpoint_env_create, &
82 kpoint_env_p_type, kpoint_env_type, kpoint_sym_create, kpoint_sym_type, kpoint_type
83 USE local_gemm_api, ONLY: local_gemm_pu_gpu,&
84 local_gemm_pu_host,&
85 local_gemm_ctxt_type
86 USE mathconstants, ONLY: twopi,&
87 z_one,&
88 z_zero
89 USE mathlib, ONLY: inv_3x3
90 USE memory_utilities, ONLY: reallocate
91 USE message_passing, ONLY: mp_cart_type,&
92 mp_para_env_type
93 USE parallel_gemm_api, ONLY: parallel_gemm
94 USE particle_types, ONLY: particle_type
95 USE qs_matrix_pools, ONLY: mpools_create,&
96 mpools_get,&
97 mpools_rebuild_fm_pools,&
98 qs_matrix_pools_type
99 USE qs_mo_types, ONLY: allocate_mo_set,&
100 get_mo_set,&
101 init_mo_set,&
102 mo_set_type,&
103 set_mo_set
104 USE qs_neighbor_list_types, ONLY: get_iterator_info,&
105 get_neighbor_list_set_p,&
106 neighbor_list_iterate,&
107 neighbor_list_iterator_create,&
108 neighbor_list_iterator_p_type,&
109 neighbor_list_iterator_release,&
110 neighbor_list_set_p_type
111 USE scf_control_types, ONLY: smear_type
112 USE smearing_utils, ONLY: smearkp,&
113 smearkp2
114 USE util, ONLY: get_limit,&
115 sort
116
117!$ USE OMP_LIB, ONLY: omp_get_max_threads, &
118!$ omp_get_thread_num
119#include "./base/base_uses.f90"
120
121 IMPLICIT NONE
122
123 PRIVATE
124
125 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpoint_methods'
126
127 TYPE, PUBLIC :: kp_transform_plan_type
128 INTEGER :: nentry = 0, ngroup = 0
129 LOGICAL :: symmetric = .false.
130 INTEGER, ALLOCATABLE, DIMENSION(:) :: col, col_offset, group_start, image, row, row_offset
131 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: cell
132 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: symmetry_sign
133 END TYPE kp_transform_plan_type
134
135 ! Non-owning view of one symmetry contribution; valid while kpoint metadata is unchanged.
136 TYPE :: kp_symmetry_op_type
137 TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kind_rot => null()
138 INTEGER, DIMENSION(:), POINTER :: atom_map => null(), atom_kind => null()
139 INTEGER, DIMENSION(:, :), POINTER :: cell_shift => null()
140 REAL(kind=dp), DIMENSION(3) :: xkp = 0.0_dp
141 LOGICAL :: time_reversal = .false., identity = .false., redistribute = .false.
142 END TYPE kp_symmetry_op_type
143
144 ! Non-owning operation metadata and owned, group-local density contributions.
145 TYPE :: kp_density_group_type
146 TYPE(kp_symmetry_op_type) :: op
147 LOGICAL :: reverse_phase = .false.
148 INTEGER, ALLOCATABLE :: inverse(:), ikp(:)
149 REAL(kind=dp), ALLOCATABLE :: xkp(:, :), weight(:)
150 END TYPE kp_density_group_type
151
152 PUBLIC :: kpoint_initialize, kpoint_env_initialize, kpoint_initialize_mos, kpoint_initialize_mo_set
153 PUBLIC :: kpoint_init_cell_index, kpoint_set_mo_occupation
154 PUBLIC :: kpoint_smearing_edge_status
155 PUBLIC :: kpoint_density_matrices, kpoint_density_transform
156 PUBLIC :: kpoint_ot_energy_weighted_density
157 PUBLIC :: kp_transform_plan_create
158 PUBLIC :: rskp_transform, lowdin_kp_trans, lowdin_kp_mo_coeff
159 PUBLIC :: rskp_grid_type, rskp_transform_grid_prepare, rskp_transform_grid_extract, &
160 rskp_transform_grid_release
161
162! **************************************************************************************************
163
164 TYPE rskp_grid_type
165 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: values
166 INTEGER, ALLOCATABLE, DIMENSION(:) :: block_col, block_nelem, block_offset, block_row
167 LOGICAL :: ready = .false.
168 END TYPE rskp_grid_type
169
170! **************************************************************************************************
171
172CONTAINS
173
174! **************************************************************************************************
175!> \brief Index full-grid links by irreducible k-point without rescanning the complete mesh.
176!> \param representative irreducible representative of each full-grid link
177!> \param nkp number of irreducible k-points
178!> \param start offsets into order for each irreducible k-point
179!> \param order full-grid link indices grouped by irreducible representative
180! **************************************************************************************************
181 SUBROUTINE kp_link_index_create(representative, nkp, start, order)
182
183 INTEGER, INTENT(IN) :: representative(:), nkp
184 INTEGER, ALLOCATABLE, INTENT(OUT) :: start(:), order(:)
185
186 INTEGER :: i, ik
187 INTEGER, ALLOCATABLE :: count(:), next(:)
188
189 ALLOCATE (count(nkp), next(nkp), start(nkp + 1), order(SIZE(representative)))
190 count(:) = 0
191 DO i = 1, SIZE(representative)
192 ik = representative(i)
193 cpassert(ik >= 1 .AND. ik <= nkp)
194 count(ik) = count(ik) + 1
195 END DO
196 start(1) = 1
197 DO ik = 1, nkp
198 start(ik + 1) = start(ik) + count(ik)
199 next(ik) = start(ik)
200 END DO
201 DO i = 1, SIZE(representative)
202 ik = representative(i)
203 order(next(ik)) = i
204 next(ik) = next(ik) + 1
205 END DO
206
207 END SUBROUTINE kp_link_index_create
208
209! **************************************************************************************************
210!> \brief Precompute lattice representations of crystallographic rotations.
211!> \param crys_sym crystallographic symmetry data
212!> \param cell simulation cell
213!> \param srot scaled-coordinate rotations
214!> \param frot integer direct-lattice rotations
215!> \param krot integer reciprocal-lattice rotations
216! **************************************************************************************************
217 SUBROUTINE kp_rotation_cache_create(crys_sym, cell, srot, frot, krot)
218
219 TYPE(csym_type), INTENT(IN) :: crys_sym
220 TYPE(cell_type), POINTER :: cell
221 REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: srot(:, :, :)
222 INTEGER, ALLOCATABLE, INTENT(OUT) :: frot(:, :, :), krot(:, :, :)
223
224 INTEGER :: ic
225
226 ALLOCATE (srot(3, 3, crys_sym%nrtot), frot(3, 3, crys_sym%nrtot), krot(3, 3, crys_sym%nrtot))
227 DO ic = 1, crys_sym%nrtot
228 srot(:, :, ic) = matmul(cell%h_inv, matmul(crys_sym%rt(:, :, ic), cell%hmat))
229 frot(:, :, ic) = nint(srot(:, :, ic))
230 krot(:, :, ic) = nint(transpose(inv_3x3(real(frot(:, :, ic), kind=dp))))
231 END DO
232
233 END SUBROUTINE kp_rotation_cache_create
234
235! **************************************************************************************************
236!> \brief Import symmetry metadata for one irreducible k-point.
237!> \param kpoint k-point environment
238!> \param crys_sym crystallographic symmetry data
239!> \param ik irreducible k-point index
240!> \param natom number of atoms
241!> \param scoord scaled atomic coordinates
242!> \param agauge periodic-image gauge of each atom
243!> \param eps_kpoint reciprocal-space matching tolerance
244!> \param link_start offsets into link_order
245!> \param link_order full-grid links grouped by irreducible representative
246!> \param srot_cache scaled-coordinate rotations
247!> \param frot_cache integer direct-lattice rotations
248!> \param krot_cache integer reciprocal-lattice rotations
249! **************************************************************************************************
250 SUBROUTINE kp_symmetry_import(kpoint, crys_sym, ik, natom, scoord, agauge, eps_kpoint, &
251 link_start, link_order, srot_cache, frot_cache, krot_cache)
252
253 TYPE(kpoint_type), POINTER :: kpoint
254 TYPE(csym_type), INTENT(IN) :: crys_sym
255 INTEGER, INTENT(IN) :: ik, natom
256 REAL(KIND=dp), INTENT(IN) :: scoord(:, :)
257 INTEGER, INTENT(IN) :: agauge(:, :)
258 REAL(KIND=dp), INTENT(IN) :: eps_kpoint
259 INTEGER, ALLOCATABLE, INTENT(IN) :: link_start(:), link_order(:)
260 REAL(KIND=dp), ALLOCATABLE, INTENT(IN) :: srot_cache(:, :, :)
261 INTEGER, ALLOCATABLE, INTENT(IN) :: frot_cache(:, :, :), krot_cache(:, :, :)
262
263 INTEGER :: ic, ilink, ir, ira, is, isign, j, nr, ns
264 INTEGER, DIMENSION(3, 3) :: frot, krot
265 REAL(KIND=dp), DIMENSION(3) :: diff, kgvec, srot
266 REAL(KIND=dp), DIMENSION(3, 3) :: srotmat
267 TYPE(kpoint_sym_type), POINTER :: kpsym
268
269 NULLIFY (kpoint%kp_sym(ik)%kpoint_sym)
270 CALL kpoint_sym_create(kpoint%kp_sym(ik)%kpoint_sym)
271 kpsym => kpoint%kp_sym(ik)%kpoint_sym
272 IF (crys_sym%nrtot <= 0 .OR. crys_sym%fullgrid .OR. &
273 crys_sym%istriz /= 1 .OR. crys_sym%inversion_only) RETURN
274
275 cpassert(ALLOCATED(link_start) .AND. ALLOCATED(link_order))
276 cpassert(ALLOCATED(srot_cache) .AND. ALLOCATED(frot_cache) .AND. ALLOCATED(krot_cache))
277 kpsym%nwght = nint(crys_sym%wkpoint(ik))
278 ns = kpsym%nwght
279 IF (ns <= 1) RETURN
280
281 ! Count mappings not already recorded by kpoint_gen.
282 DO ilink = link_start(ik), link_start(ik + 1) - 1
283 is = link_order(ilink)
284 DO ic = 1, crys_sym%nrtot
285 krot(1:3, 1:3) = krot_cache(1:3, 1:3, ic)
286 DO isign = 1, 2
287 ir = merge(crys_sym%ibrot(ic), -crys_sym%ibrot(ic), isign == 1)
288 IF (ir == crys_sym%kpop(is)) cycle
289 kgvec(1:3) = crys_sym%kpmesh(1:3, is) - &
290 matmul(real(merge(krot(1:3, 1:3), -krot(1:3, 1:3), isign == 1), kind=dp), &
291 kpoint%xkp(1:3, ik))
292 diff(1:3) = kgvec(1:3) - anint(kgvec(1:3))
293 IF (all(abs(diff(1:3)) < eps_kpoint)) ns = ns + 1
294 END DO
295 END DO
296 END DO
297
298 kpsym%apply_symmetry = .true.
299 ALLOCATE (kpsym%rot(3, 3, ns), kpsym%xkp(3, ns), kpsym%rotp(ns))
300 ALLOCATE (kpsym%f0(natom, ns), kpsym%fcell(3, natom, ns), kpsym%fcell_gauge(3, natom, ns))
301 ALLOCATE (kpsym%phase_mode(ns), kpsym%kgphase(natom, ns))
302 kpsym%phase_mode(:) = 0
303
304 nr = 0
305 DO ilink = link_start(ik), link_start(ik + 1) - 1
306 is = link_order(ilink)
307 nr = nr + 1
308 ir = crys_sym%kpop(is)
309 ira = abs(ir)
310 DO ic = 1, crys_sym%nrtot
311 IF (crys_sym%ibrot(ic) /= ira) cycle
312 kpsym%rotp(nr) = ir
313 kpsym%rot(1:3, 1:3, nr) = crys_sym%rt(1:3, 1:3, ic)
314 srotmat(1:3, 1:3) = srot_cache(1:3, 1:3, ic)
315 frot(1:3, 1:3) = frot_cache(1:3, 1:3, ic)
316 kpsym%xkp(1:3, nr) = crys_sym%kpmesh(1:3, is)
317 krot(1:3, 1:3) = krot_cache(1:3, 1:3, ic)
318 IF (ir < 0) krot(1:3, 1:3) = -krot(1:3, 1:3)
319 kgvec(1:3) = kpsym%xkp(1:3, nr) - &
320 matmul(real(krot(1:3, 1:3), kind=dp), kpoint%xkp(1:3, ik))
321 kgvec(1:3) = anint(kgvec(1:3))
322 kpsym%f0(1:natom, nr) = crys_sym%f0(1:natom, ic)
323 DO j = 1, natom
324 srot(1:3) = matmul(srotmat, scoord(1:3, j)) + crys_sym%vt(1:3, ic)
325 kpsym%fcell(1:3, j, nr) = nint(srot(1:3) - scoord(1:3, kpsym%f0(j, nr)))
326 kpsym%fcell_gauge(1:3, j, nr) = kpsym%fcell(1:3, j, nr) + &
327 matmul(frot, agauge(1:3, j)) - &
328 agauge(1:3, kpsym%f0(j, nr))
329 kpsym%kgphase(j, nr) = dot_product(kgvec(1:3), &
330 scoord(1:3, j) + real(agauge(1:3, j), kind=dp))
331 END DO
332 EXIT
333 END DO
334 cpassert(ic <= crys_sym%nrtot)
335 END DO
336
337 ! Complete each star with equivalent operations not selected by kpoint_gen.
338 DO ilink = link_start(ik), link_start(ik + 1) - 1
339 is = link_order(ilink)
340 DO ic = 1, crys_sym%nrtot
341 srotmat(1:3, 1:3) = srot_cache(1:3, 1:3, ic)
342 frot(1:3, 1:3) = frot_cache(1:3, 1:3, ic)
343 krot(1:3, 1:3) = krot_cache(1:3, 1:3, ic)
344 DO isign = 1, 2
345 ir = merge(crys_sym%ibrot(ic), -crys_sym%ibrot(ic), isign == 1)
346 IF (ir == crys_sym%kpop(is)) cycle
347 kgvec(1:3) = crys_sym%kpmesh(1:3, is) - &
348 matmul(real(merge(krot(1:3, 1:3), -krot(1:3, 1:3), isign == 1), kind=dp), &
349 kpoint%xkp(1:3, ik))
350 diff(1:3) = kgvec(1:3) - anint(kgvec(1:3))
351 IF (.NOT. all(abs(diff(1:3)) < eps_kpoint)) cycle
352 nr = nr + 1
353 kpsym%rotp(nr) = ir
354 kpsym%rot(1:3, 1:3, nr) = crys_sym%rt(1:3, 1:3, ic)
355 kpsym%xkp(1:3, nr) = crys_sym%kpmesh(1:3, is)
356 kgvec(1:3) = anint(kgvec(1:3))
357 kpsym%f0(1:natom, nr) = crys_sym%f0(1:natom, ic)
358 DO j = 1, natom
359 srot(1:3) = matmul(srotmat, scoord(1:3, j)) + crys_sym%vt(1:3, ic)
360 kpsym%fcell(1:3, j, nr) = nint(srot(1:3) - scoord(1:3, kpsym%f0(j, nr)))
361 kpsym%fcell_gauge(1:3, j, nr) = kpsym%fcell(1:3, j, nr) + &
362 matmul(frot, agauge(1:3, j)) - &
363 agauge(1:3, kpsym%f0(j, nr))
364 kpsym%kgphase(j, nr) = dot_product(kgvec(1:3), &
365 scoord(1:3, j) + real(agauge(1:3, j), kind=dp))
366 END DO
367 END DO
368 END DO
369 END DO
370 cpassert(nr == ns)
371 kpsym%nwred = nr
372
373 END SUBROUTINE kp_symmetry_import
374
375! **************************************************************************************************
376!> \brief Generate the kpoints and initialize the kpoint environment
377!> \param kpoint The kpoint environment
378!> \param particle_set Particle types and coordinates
379!> \param cell Computational cell information
380! **************************************************************************************************
381 SUBROUTINE kpoint_initialize(kpoint, particle_set, cell)
382
383 TYPE(kpoint_type), POINTER :: kpoint
384 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
385 TYPE(cell_type), POINTER :: cell
386
387 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_initialize'
388
389 INTEGER :: handle, i, ik, iounit, j, natom, nkind, &
390 ns
391 INTEGER, ALLOCATABLE, DIMENSION(:) :: atype, kp_link_order, kp_link_start
392 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: agauge
393 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: frot_cache, krot_cache
394 LOGICAL :: spez
395 REAL(kind=dp) :: eps_kpoint, wsum
396 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coord, scoord
397 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: srot_cache
398 REAL(kind=dp), DIMENSION(3) :: r_pbc, scoord_pbc
399 REAL(kind=dp), DIMENSION(:), POINTER :: wkp_full
400 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp_full
401 TYPE(csym_type) :: crys_sym
402
403 CALL timeset(routinen, handle)
404
405 cpassert(ASSOCIATED(kpoint))
406
407 SELECT CASE (kpoint%kp_scheme)
408 CASE ("NONE")
409 ! do nothing
410 CASE ("GAMMA")
411 kpoint%nkp = 1
412 ALLOCATE (kpoint%xkp(3, 1), kpoint%wkp(1))
413 kpoint%xkp(1:3, 1) = 0.0_dp
414 kpoint%wkp(1) = 1.0_dp
415 ALLOCATE (kpoint%kp_sym(1))
416 NULLIFY (kpoint%kp_sym(1)%kpoint_sym)
417 CALL kpoint_sym_create(kpoint%kp_sym(1)%kpoint_sym)
418 CASE ("MONKHORST-PACK", "MACDONALD")
419
420 IF (.NOT. kpoint%symmetry) THEN
421 ! we set up a random molecule to avoid any possible symmetry
422 natom = 10
423 ALLOCATE (coord(3, natom), scoord(3, natom), atype(natom))
424 DO i = 1, natom
425 atype(i) = i
426 coord(1, i) = sin(i*0.12345_dp)
427 coord(2, i) = cos(i*0.23456_dp)
428 coord(3, i) = sin(i*0.34567_dp)
429 CALL real_to_scaled(scoord(1:3, i), coord(1:3, i), cell)
430 END DO
431 ELSE
432 natom = SIZE(particle_set)
433 ALLOCATE (scoord(3, natom), atype(natom))
434 DO i = 1, natom
435 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=atype(i))
436 CALL real_to_scaled(scoord(1:3, i), particle_set(i)%r(1:3), cell)
437 END DO
438 END IF
439 IF (kpoint%verbose) THEN
440 iounit = cp_logger_get_default_io_unit()
441 ELSE
442 iounit = -1
443 END IF
444 ! kind type list
445 ALLOCATE (kpoint%atype(natom))
446 kpoint%atype = atype
447 ! Match the atom images used by CP2K's periodic neighbor lists.
448 ALLOCATE (agauge(3, natom))
449 agauge = 0
450 IF (kpoint%symmetry) THEN
451 DO i = 1, natom
452 r_pbc(1:3) = pbc_stable(particle_set(i)%r(1:3), cell)
453 CALL real_to_scaled(scoord_pbc, r_pbc, cell)
454 agauge(1:3, i) = nint(scoord_pbc(1:3) - scoord(1:3, i))
455 END DO
456 END IF
457
458 CALL crys_sym_gen(crys_sym, scoord, atype, cell%hmat, delta=kpoint%eps_geo, iounit=iounit, &
459 use_spglib=kpoint%symmetry)
460 CALL kpoint_gen(crys_sym, kpoint%nkp_grid, symm=kpoint%symmetry, shift=kpoint%kp_shift, &
461 full_grid=kpoint%full_grid, gamma_centered=kpoint%gamma_centered, &
462 inversion_symmetry_only=kpoint%inversion_symmetry_only, &
463 use_spglib_reduction= &
464 kpoint%symmetry_reduction_method == use_spglib_kpoint_symmetry, &
465 use_spglib_backend=kpoint%symmetry_backend == use_spglib_kpoint_backend)
466 IF (crys_sym%inversion_only) kpoint%inversion_symmetry_only = .true.
467 kpoint%nkp = crys_sym%nkpoint
468 ALLOCATE (kpoint%xkp(3, kpoint%nkp), kpoint%wkp(kpoint%nkp))
469 wsum = sum(crys_sym%wkpoint)
470 DO ik = 1, kpoint%nkp
471 kpoint%xkp(1:3, ik) = crys_sym%xkpoint(1:3, ik)
472 kpoint%wkp(ik) = crys_sym%wkpoint(ik)/wsum
473 END DO
474 IF (crys_sym%nrtot > 0 .AND. .NOT. crys_sym%fullgrid .AND. &
475 crys_sym%istriz == 1 .AND. .NOT. crys_sym%inversion_only) THEN
476 CALL kp_link_index_create(crys_sym%kplink(2, :), kpoint%nkp, kp_link_start, kp_link_order)
477 CALL kp_rotation_cache_create(crys_sym, cell, srot_cache, frot_cache, krot_cache)
478 END IF
479
480 eps_kpoint = max(1.e-12_dp, 10.0_dp*kpoint%eps_geo)
481 ! print output
482 IF (kpoint%symmetry) CALL print_crys_symmetry(crys_sym)
483 IF (kpoint%symmetry) CALL print_kp_symmetry(crys_sym)
484
485 ! transfer symmetry information
486 ALLOCATE (kpoint%kp_sym(kpoint%nkp))
487 DO ik = 1, kpoint%nkp
488 CALL kp_symmetry_import(kpoint, crys_sym, ik, natom, scoord, agauge, eps_kpoint, &
489 kp_link_start, kp_link_order, srot_cache, frot_cache, krot_cache)
490 END DO
491 IF (kpoint%symmetry) THEN
492 nkind = maxval(atype)
493 ns = crys_sym%nrtot
494 ALLOCATE (kpoint%kind_rotmat(ns, nkind))
495 DO i = 1, ns
496 DO j = 1, nkind
497 NULLIFY (kpoint%kind_rotmat(i, j)%rmat)
498 END DO
499 END DO
500 ALLOCATE (kpoint%ibrot(ns))
501 kpoint%ibrot(1:ns) = crys_sym%ibrot(1:ns)
502 END IF
503
504 CALL release_csym_type(crys_sym)
505 DEALLOCATE (scoord, atype)
506 DEALLOCATE (agauge)
507
508 CASE ("GENERAL")
509 NULLIFY (xkp_full, wkp_full)
510 IF (ASSOCIATED(kpoint%xkp_input)) THEN
511 xkp_full => kpoint%xkp_input
512 wkp_full => kpoint%wkp_input
513 ELSE
514 xkp_full => kpoint%xkp
515 wkp_full => kpoint%wkp
516 END IF
517 cpassert(ASSOCIATED(xkp_full))
518 cpassert(ASSOCIATED(wkp_full))
519 IF (.NOT. ASSOCIATED(kpoint%xkp_input)) THEN
520 ALLOCATE (kpoint%xkp_input(3, SIZE(wkp_full)), kpoint%wkp_input(SIZE(wkp_full)))
521 kpoint%xkp_input(1:3, 1:SIZE(wkp_full)) = xkp_full(1:3, 1:SIZE(wkp_full))
522 kpoint%wkp_input(1:SIZE(wkp_full)) = wkp_full(1:SIZE(wkp_full))
523 xkp_full => kpoint%xkp_input
524 wkp_full => kpoint%wkp_input
525 END IF
526 IF (.NOT. kpoint%symmetry) THEN
527 IF (.NOT. ASSOCIATED(kpoint%xkp)) THEN
528 kpoint%nkp = SIZE(wkp_full)
529 ALLOCATE (kpoint%xkp(3, kpoint%nkp), kpoint%wkp(kpoint%nkp))
530 kpoint%xkp(1:3, 1:kpoint%nkp) = xkp_full(1:3, 1:kpoint%nkp)
531 kpoint%wkp(1:kpoint%nkp) = wkp_full(1:kpoint%nkp)
532 END IF
533 ! default: no symmetry settings
534 ALLOCATE (kpoint%kp_sym(kpoint%nkp))
535 DO i = 1, kpoint%nkp
536 NULLIFY (kpoint%kp_sym(i)%kpoint_sym)
537 CALL kpoint_sym_create(kpoint%kp_sym(i)%kpoint_sym)
538 END DO
539 ELSE
540 IF (kpoint%verbose) THEN
541 iounit = cp_logger_get_default_io_unit()
542 ELSE
543 iounit = -1
544 END IF
545 natom = SIZE(particle_set)
546 ALLOCATE (scoord(3, natom), atype(natom))
547 DO i = 1, natom
548 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=atype(i))
549 CALL real_to_scaled(scoord(1:3, i), particle_set(i)%r(1:3), cell)
550 END DO
551 ALLOCATE (kpoint%atype(natom))
552 kpoint%atype = atype
553 ALLOCATE (agauge(3, natom))
554 DO i = 1, natom
555 r_pbc(1:3) = pbc_stable(particle_set(i)%r(1:3), cell)
556 CALL real_to_scaled(scoord_pbc, r_pbc, cell)
557 agauge(1:3, i) = nint(scoord_pbc(1:3) - scoord(1:3, i))
558 END DO
559
560 CALL crys_sym_gen(crys_sym, scoord, atype, cell%hmat, delta=kpoint%eps_geo, iounit=iounit, &
561 use_spglib=(kpoint%symmetry_backend == use_spglib_kpoint_backend .OR. &
562 kpoint%symmetry_reduction_method == use_spglib_kpoint_symmetry))
563 CALL kpoint_gen_general(crys_sym, xkp_full, wkp_full, symm=kpoint%symmetry, &
564 full_grid=kpoint%full_grid, &
565 inversion_symmetry_only=kpoint%inversion_symmetry_only, &
566 use_spglib_reduction= &
567 kpoint%symmetry_reduction_method == use_spglib_kpoint_symmetry, &
568 use_spglib_backend=kpoint%symmetry_backend == use_spglib_kpoint_backend)
569 IF (crys_sym%inversion_only) kpoint%inversion_symmetry_only = .true.
570 IF (ASSOCIATED(kpoint%xkp)) THEN
571 DEALLOCATE (kpoint%xkp)
572 NULLIFY (kpoint%xkp)
573 END IF
574 IF (ASSOCIATED(kpoint%wkp)) THEN
575 DEALLOCATE (kpoint%wkp)
576 NULLIFY (kpoint%wkp)
577 END IF
578 kpoint%nkp = crys_sym%nkpoint
579 ALLOCATE (kpoint%xkp(3, kpoint%nkp), kpoint%wkp(kpoint%nkp))
580 wsum = sum(crys_sym%wkpoint)
581 DO ik = 1, kpoint%nkp
582 kpoint%xkp(1:3, ik) = crys_sym%xkpoint(1:3, ik)
583 kpoint%wkp(ik) = crys_sym%wkpoint(ik)/wsum
584 END DO
585 IF (crys_sym%nrtot > 0 .AND. .NOT. crys_sym%fullgrid .AND. &
586 crys_sym%istriz == 1 .AND. .NOT. crys_sym%inversion_only) THEN
587 CALL kp_link_index_create(crys_sym%kplink(2, :), kpoint%nkp, kp_link_start, kp_link_order)
588 CALL kp_rotation_cache_create(crys_sym, cell, srot_cache, frot_cache, krot_cache)
589 END IF
590
591 eps_kpoint = max(1.e-12_dp, 10.0_dp*kpoint%eps_geo)
592 CALL print_crys_symmetry(crys_sym)
593 CALL print_kp_symmetry(crys_sym)
594
595 ALLOCATE (kpoint%kp_sym(kpoint%nkp))
596 DO ik = 1, kpoint%nkp
597 CALL kp_symmetry_import(kpoint, crys_sym, ik, natom, scoord, agauge, eps_kpoint, &
598 kp_link_start, kp_link_order, srot_cache, frot_cache, krot_cache)
599 END DO
600 nkind = maxval(atype)
601 ns = crys_sym%nrtot
602 ALLOCATE (kpoint%kind_rotmat(ns, nkind))
603 DO i = 1, ns
604 DO j = 1, nkind
605 NULLIFY (kpoint%kind_rotmat(i, j)%rmat)
606 END DO
607 END DO
608 ALLOCATE (kpoint%ibrot(ns))
609 kpoint%ibrot(1:ns) = crys_sym%ibrot(1:ns)
610
611 CALL release_csym_type(crys_sym)
612 DEALLOCATE (scoord, atype)
613 DEALLOCATE (agauge)
614 END IF
615 CASE DEFAULT
616 cpabort("Option invalid or unavailable for kpoint%kp_scheme")
617 END SELECT
618
619 ! check for consistency of options
620 SELECT CASE (kpoint%kp_scheme)
621 CASE ("NONE")
622 ! don't use k-point code
623 CASE ("GAMMA")
624 cpassert(kpoint%nkp == 1)
625 cpassert(sum(abs(kpoint%xkp)) <= 1.e-12_dp)
626 cpassert(kpoint%wkp(1) == 1.0_dp)
627 cpassert(.NOT. kpoint%symmetry)
628 CASE ("GENERAL")
629 cpassert(kpoint%nkp >= 1)
630 CASE ("MONKHORST-PACK", "MACDONALD")
631 cpassert(kpoint%nkp >= 1)
632 END SELECT
633 IF (kpoint%use_real_wfn) THEN
634 ! what about inversion symmetry?
635 ikloop: DO ik = 1, kpoint%nkp
636 DO i = 1, 3
637 spez = (kpoint%xkp(i, ik) == 0.0_dp .OR. kpoint%xkp(i, ik) == 0.5_dp)
638 IF (.NOT. spez) EXIT ikloop
639 END DO
640 END DO ikloop
641 IF (.NOT. spez) THEN
642 ! Warning: real wfn might be wrong for this system
643 CALL cp_warn(__location__, &
644 "A calculation using real wavefunctions is requested. "// &
645 "We could not determine if the symmetry of the system allows real wavefunctions. ")
646 END IF
647 END IF
648
649 CALL timestop(handle)
650
651 END SUBROUTINE kpoint_initialize
652
653! **************************************************************************************************
654!> \brief Initialize the kpoint environment
655!> \param kpoint Kpoint environment
656!> \param para_env ...
657!> \param blacs_env ...
658!> \param with_aux_fit ...
659! **************************************************************************************************
660 SUBROUTINE kpoint_env_initialize(kpoint, para_env, blacs_env, with_aux_fit)
661
662 TYPE(kpoint_type), INTENT(INOUT) :: kpoint
663 TYPE(mp_para_env_type), INTENT(IN), TARGET :: para_env
664 TYPE(cp_blacs_env_type), INTENT(IN), TARGET :: blacs_env
665 LOGICAL, INTENT(IN), OPTIONAL :: with_aux_fit
666
667 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_env_initialize'
668
669 INTEGER :: handle, igr, ngr, niogrp, nkp, &
670 nkp_grp, nkp_loc, npe, unit_nr
671 INTEGER, DIMENSION(2) :: dims, pos
672 LOGICAL :: aux_fit
673 TYPE(mp_cart_type) :: comm_cart
674 TYPE(mp_para_env_type), POINTER :: para_env_inter_kp, para_env_kp
675
676 CALL timeset(routinen, handle)
677
678 IF (PRESENT(with_aux_fit)) THEN
679 aux_fit = with_aux_fit
680 ELSE
681 aux_fit = .false.
682 END IF
683
684 kpoint%para_env => para_env
685 CALL kpoint%para_env%retain()
686 kpoint%blacs_env_all => blacs_env
687 CALL kpoint%blacs_env_all%retain()
688
689 cpassert(.NOT. ASSOCIATED(kpoint%kp_env))
690 IF (aux_fit) THEN
691 cpassert(.NOT. ASSOCIATED(kpoint%kp_aux_env))
692 END IF
693
694 nkp = kpoint%nkp
695 npe = para_env%num_pe
696 IF (npe == 1) THEN
697 ! only one process available -> owns all kpoints
698 ALLOCATE (kpoint%kp_dist(2, 1))
699 kpoint%kp_dist(1, 1) = 1
700 kpoint%kp_dist(2, 1) = nkp
701 kpoint%kp_range(1) = 1
702 kpoint%kp_range(2) = nkp
703
704 ! parallel environments
705 kpoint%para_env_kp => para_env
706 CALL kpoint%para_env_kp%retain()
707 kpoint%para_env_inter_kp => para_env
708 CALL kpoint%para_env_inter_kp%retain()
709 kpoint%iogrp = .true.
710 kpoint%nkp_groups = 1
711 ELSE
712 IF (kpoint%parallel_group_size == -1) THEN
713 ! Maximum parallelization over kpoints with equal-sized MPI groups.
714 ! Each group must own at least one kpoint; their kpoint counts may differ.
715 DO igr = npe, 1, -1
716 IF (mod(npe, igr) /= 0) cycle
717 nkp_grp = npe/igr
718 IF (nkp_grp > nkp) cycle
719 ngr = igr
720 END DO
721 ELSE IF (kpoint%parallel_group_size == 0) THEN
722 ! no parallelization over kpoints
723 ngr = npe
724 ELSE IF (kpoint%parallel_group_size > 0) THEN
725 ngr = min(kpoint%parallel_group_size, npe)
726 ELSE
727 cpabort("kpoint%parallel_group_size cannot be smaller than -1")
728 END IF
729 nkp_grp = npe/ngr
730 ! processor dimensions
731 dims(1) = ngr
732 dims(2) = nkp_grp
733 IF ((dims(1)*dims(2) /= npe)) THEN
734 cpabort("Number of processors is not divisible by the kpoint group size.")
735 END IF
736 IF (nkp_grp > nkp) THEN
737 cpabort("Too many kpoint groups. Increase PARALLEL_GROUP_SIZE.")
738 END IF
739
740 ! Create the subgroups, one for each k-point group and one interconnecting group
741 CALL comm_cart%create(comm_old=para_env, ndims=2, dims=dims)
742 pos = comm_cart%mepos_cart
743 ALLOCATE (para_env_kp)
744 CALL para_env_kp%from_split(comm_cart, pos(2))
745 ALLOCATE (para_env_inter_kp)
746 CALL para_env_inter_kp%from_split(comm_cart, pos(1))
747 CALL comm_cart%free()
748
749 niogrp = 0
750 IF (para_env%is_source()) niogrp = 1
751 CALL para_env_kp%sum(niogrp)
752 kpoint%iogrp = (niogrp == 1)
753
754 ! parallel groups
755 kpoint%para_env_kp => para_env_kp
756 kpoint%para_env_inter_kp => para_env_inter_kp
757
758 ! distribution of kpoints
759 ALLOCATE (kpoint%kp_dist(2, nkp_grp))
760 DO igr = 1, nkp_grp
761 kpoint%kp_dist(1:2, igr) = get_limit(nkp, nkp_grp, igr - 1)
762 END DO
763 ! local kpoints
764 kpoint%kp_range(1:2) = kpoint%kp_dist(1:2, para_env_inter_kp%mepos + 1)
765 nkp_loc = kpoint%kp_range(2) - kpoint%kp_range(1) + 1
766
767 unit_nr = cp_logger_get_default_io_unit()
768
769 IF (unit_nr > 0 .AND. kpoint%verbose) THEN
770 WRITE (unit_nr, *)
771 WRITE (unit_nr, fmt="(T2,A,T71,I10)") "KPOINTS| Total number of kpoints", nkp
772 WRITE (unit_nr, fmt="(T2,A,T71,I10)") "KPOINTS| Number of kpoint groups ", nkp_grp
773 WRITE (unit_nr, fmt="(T2,A,T71,I10)") "KPOINTS| Size of each kpoint group", ngr
774 IF (mod(nkp, nkp_grp) == 0) THEN
775 WRITE (unit_nr, fmt="(T2,A,T71,I10)") "KPOINTS| Number of kpoints per group", nkp_loc
776 ELSE
777 WRITE (unit_nr, fmt="(T2,A,T71,I10)") "KPOINTS| Minimum number of kpoints per group", nkp/nkp_grp
778 WRITE (unit_nr, fmt="(T2,A,T71,I10)") "KPOINTS| Maximum number of kpoints per group", nkp/nkp_grp + 1
779 END IF
780 END IF
781 kpoint%nkp_groups = nkp_grp
782
783 END IF
784
785 nkp_loc = kpoint%kp_range(2) - kpoint%kp_range(1) + 1
786 CALL create_local_environments(kpoint%kp_env)
787 IF (aux_fit) CALL create_local_environments(kpoint%kp_aux_env)
788
789 CALL timestop(handle)
790
791 CONTAINS
792
793! **************************************************************************************************
794!> \brief ...
795!> \param env ...
796! **************************************************************************************************
797 SUBROUTINE create_local_environments(env)
798 TYPE(kpoint_env_p_type), INTENT(OUT), POINTER :: env(:)
799
800 INTEGER :: ik, ikk
801 TYPE(kpoint_env_type), POINTER :: kp
802
803 ALLOCATE (env(nkp_loc))
804 DO ik = 1, nkp_loc
805 ikk = kpoint%kp_range(1) + ik - 1
806 CALL kpoint_env_create(env(ik)%kpoint_env)
807 kp => env(ik)%kpoint_env
808 kp%nkpoint = ikk
809 kp%wkp = kpoint%wkp(ikk)
810 kp%xkp(1:3) = kpoint%xkp(1:3, ikk)
811 kp%is_local = (kpoint%para_env_kp%num_pe == 1)
812 END DO
813 END SUBROUTINE create_local_environments
814
815 END SUBROUTINE kpoint_env_initialize
816
817! **************************************************************************************************
818!> \brief Initialize a set of MOs and density matrix for each kpoint (kpoint group)
819!> \param kpoint Kpoint environment
820!> \param mos Reference MOs (global)
821!> \param added_mos ...
822!> \param for_aux_fit ...
823! **************************************************************************************************
824 SUBROUTINE kpoint_initialize_mos(kpoint, mos, added_mos, for_aux_fit)
825
826 TYPE(kpoint_type), POINTER :: kpoint
827 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
828 INTEGER, INTENT(IN), OPTIONAL :: added_mos
829 LOGICAL, OPTIONAL :: for_aux_fit
830
831 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_initialize_mos'
832
833 INTEGER :: handle, ic, ik, is, nadd, nao, nc, &
834 nelectron, nkp_loc, nmo, nmorig(2), &
835 nspin
836 LOGICAL :: aux_fit
837 REAL(kind=dp) :: flexible_electron_count, maxocc, n_el_f
838 TYPE(cp_blacs_env_type), POINTER :: blacs_env
839 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools
840 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
841 TYPE(cp_fm_type), POINTER :: fmlocal
842 TYPE(kpoint_env_type), POINTER :: kp
843 TYPE(qs_matrix_pools_type), POINTER :: mpools
844
845 CALL timeset(routinen, handle)
846
847 IF (PRESENT(for_aux_fit)) THEN
848 aux_fit = for_aux_fit
849 ELSE
850 aux_fit = .false.
851 END IF
852
853 cpassert(ASSOCIATED(kpoint))
854
855 IF (aux_fit) THEN
856 cpassert(ASSOCIATED(kpoint%kp_aux_env))
857 END IF
858
859 IF (PRESENT(added_mos)) THEN
860 nadd = added_mos
861 ELSE
862 nadd = 0
863 END IF
864
865 IF (kpoint%use_real_wfn) THEN
866 nc = 1
867 ELSE
868 nc = 2
869 END IF
870 nspin = SIZE(mos, 1)
871 nkp_loc = kpoint%kp_range(2) - kpoint%kp_range(1) + 1
872 IF (nkp_loc > 0) THEN
873 IF (aux_fit) THEN
874 cpassert(SIZE(kpoint%kp_aux_env) == nkp_loc)
875 ELSE
876 cpassert(SIZE(kpoint%kp_env) == nkp_loc)
877 END IF
878 ! Allocate one MO set per spin; its coefficient matrix owns the scalar type.
879 DO ik = 1, nkp_loc
880 IF (aux_fit) THEN
881 kp => kpoint%kp_aux_env(ik)%kpoint_env
882 ELSE
883 kp => kpoint%kp_env(ik)%kpoint_env
884 END IF
885 ALLOCATE (kp%mos(nspin))
886 DO is = 1, nspin
887 CALL get_mo_set(mos(is), nao=nao, nmo=nmo, nelectron=nelectron, &
888 n_el_f=n_el_f, maxocc=maxocc, flexible_electron_count=flexible_electron_count)
889 nmo = min(nao, nmo + nadd)
890 CALL allocate_mo_set(kp%mos(is), nao, nmo, nelectron, n_el_f, maxocc, &
891 flexible_electron_count)
892 END DO
893 ! freshly allocated MOS carry no coefficients: solvers must
894 ! cold-start again until an extrapolation refills them
895 kp%mos_prefilled = .false.
896 END DO
897
898 ! generate the blacs environment for the kpoint group
899 ! we generate a blacs env for each kpoint group in parallel
900 ! we assume here that the group para_env_inter_kp will connect
901 ! equivalent parts of fm matrices, i.e. no reshuffeling of processors
902 NULLIFY (blacs_env)
903 IF (ASSOCIATED(kpoint%blacs_env)) THEN
904 blacs_env => kpoint%blacs_env
905 ELSE
906 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=kpoint%para_env_kp)
907 kpoint%blacs_env => blacs_env
908 END IF
909
910 ! set possible new number of MOs
911 DO is = 1, nspin
912 CALL get_mo_set(mos(is), nmo=nmorig(is))
913 nmo = min(nao, nmorig(is) + nadd)
914 CALL set_mo_set(mos(is), nmo=nmo)
915 END DO
916 ! matrix pools for the kpoint group, information on MOs is transferred using
917 ! generic mos structure
918 NULLIFY (mpools)
919 CALL mpools_create(mpools=mpools)
920 CALL mpools_rebuild_fm_pools(mpools=mpools, mos=mos, &
921 blacs_env=blacs_env, para_env=kpoint%para_env_kp)
922
923 IF (aux_fit) THEN
924 kpoint%mpools_aux_fit => mpools
925 ELSE
926 kpoint%mpools => mpools
927 END IF
928
929 ! reset old number of MOs
930 DO is = 1, nspin
931 CALL set_mo_set(mos(is), nmo=nmorig(is))
932 END DO
933
934 ! allocate density matrices
935 CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
936 ALLOCATE (fmlocal)
937 CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
938 CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
939 DO ik = 1, nkp_loc
940 IF (aux_fit) THEN
941 kp => kpoint%kp_aux_env(ik)%kpoint_env
942 ELSE
943 kp => kpoint%kp_env(ik)%kpoint_env
944 END IF
945 ! density matrix
946 CALL cp_fm_release(kp%pmat)
947 ALLOCATE (kp%pmat(nc, nspin))
948 DO is = 1, nspin
949 DO ic = 1, nc
950 CALL cp_fm_create(kp%pmat(ic, is), matrix_struct)
951 END DO
952 END DO
953 ! energy weighted density matrix
954 CALL cp_fm_release(kp%wmat)
955 ALLOCATE (kp%wmat(nc, nspin))
956 DO is = 1, nspin
957 DO ic = 1, nc
958 CALL cp_fm_create(kp%wmat(ic, is), matrix_struct)
959 END DO
960 END DO
961 END DO
962 CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
963 DEALLOCATE (fmlocal)
964
965 END IF
966
967 CALL timestop(handle)
968
969 END SUBROUTINE kpoint_initialize_mos
970
971! **************************************************************************************************
972!> \brief ...
973!> \param kpoint ...
974! **************************************************************************************************
975 SUBROUTINE kpoint_initialize_mo_set(kpoint)
976 TYPE(kpoint_type), POINTER :: kpoint
977
978 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_initialize_mo_set'
979
980 INTEGER :: handle, ik, ispin
981 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
982 TYPE(mo_set_type), DIMENSION(:), POINTER :: moskp
983
984 CALL timeset(routinen, handle)
985
986 DO ik = 1, SIZE(kpoint%kp_env)
987 CALL mpools_get(kpoint%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
988 moskp => kpoint%kp_env(ik)%kpoint_env%mos
989 cpassert(ASSOCIATED(moskp))
990 DO ispin = 1, SIZE(moskp)
991 IF (ASSOCIATED(moskp(ispin)%mo_coeff) .OR. &
992 ASSOCIATED(moskp(ispin)%cmo_coeff)) cycle
993 CALL init_mo_set(moskp(ispin), fm_pool=ao_mo_fm_pools(ispin)%pool, &
994 name="kpoints", complex_coeff=.NOT. kpoint%use_real_wfn)
995 END DO
996 END DO
997
998 CALL timestop(handle)
999
1000 END SUBROUTINE kpoint_initialize_mo_set
1001
1002! **************************************************************************************************
1003!> \brief Generates the mapping of cell indices and linear RS index
1004!> CELL (0,0,0) is always mapped to index 1
1005!> \param kpoint Kpoint environment
1006!> \param sab_nl Defining neighbour list
1007!> \param para_env Parallel environment
1008!> \param nimages [output]
1009! **************************************************************************************************
1010 SUBROUTINE kpoint_init_cell_index(kpoint, sab_nl, para_env, nimages)
1011
1012 TYPE(kpoint_type), POINTER :: kpoint
1013 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1014 POINTER :: sab_nl
1015 TYPE(mp_para_env_type), POINTER :: para_env
1016 INTEGER, INTENT(OUT) :: nimages
1017
1018 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_init_cell_index'
1019
1020 INTEGER :: handle, i1, i2, i3, ic, icount, it, &
1021 ncount
1022 INTEGER, DIMENSION(3) :: cell, itm
1023 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell, list
1024 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index, cti
1025 LOGICAL :: new
1026 TYPE(neighbor_list_iterator_p_type), &
1027 DIMENSION(:), POINTER :: nl_iterator
1028
1029 NULLIFY (cell_to_index, index_to_cell)
1030
1031 CALL timeset(routinen, handle)
1032
1033 cpassert(ASSOCIATED(kpoint))
1034
1035 ALLOCATE (list(3, 125))
1036 list = 0
1037 icount = 1
1038
1039 CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
1040 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1041 CALL get_iterator_info(nl_iterator, cell=cell)
1042
1043 new = .true.
1044 DO ic = 1, icount
1045 IF (cell(1) == list(1, ic) .AND. cell(2) == list(2, ic) .AND. &
1046 cell(3) == list(3, ic)) THEN
1047 new = .false.
1048 EXIT
1049 END IF
1050 END DO
1051 IF (new) THEN
1052 icount = icount + 1
1053 IF (icount > SIZE(list, 2)) THEN
1054 CALL reallocate(list, 1, 3, 1, 2*SIZE(list, 2))
1055 END IF
1056 list(1:3, icount) = cell(1:3)
1057 END IF
1058
1059 END DO
1060 CALL neighbor_list_iterator_release(nl_iterator)
1061
1062 itm(1) = maxval(abs(list(1, 1:icount)))
1063 itm(2) = maxval(abs(list(2, 1:icount)))
1064 itm(3) = maxval(abs(list(3, 1:icount)))
1065 CALL para_env%max(itm)
1066 it = maxval(itm(1:3))
1067 IF (ASSOCIATED(kpoint%cell_to_index)) THEN
1068 DEALLOCATE (kpoint%cell_to_index)
1069 END IF
1070 ALLOCATE (kpoint%cell_to_index(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
1071 cell_to_index => kpoint%cell_to_index
1072 cti => cell_to_index
1073 cti(:, :, :) = 0
1074 DO ic = 1, icount
1075 i1 = list(1, ic)
1076 i2 = list(2, ic)
1077 i3 = list(3, ic)
1078 cti(i1, i2, i3) = ic
1079 END DO
1080 CALL para_env%sum(cti)
1081 ncount = 0
1082 DO i1 = -itm(1), itm(1)
1083 DO i2 = -itm(2), itm(2)
1084 DO i3 = -itm(3), itm(3)
1085 IF (cti(i1, i2, i3) == 0) THEN
1086 cti(i1, i2, i3) = 1000000
1087 ELSE
1088 ncount = ncount + 1
1089 cti(i1, i2, i3) = (abs(i1) + abs(i2) + abs(i3))*1000 + abs(i3)*100 + abs(i2)*10 + abs(i1)
1090 cti(i1, i2, i3) = cti(i1, i2, i3) + (i1 + i2 + i3)
1091 END IF
1092 END DO
1093 END DO
1094 END DO
1095
1096 IF (ASSOCIATED(kpoint%index_to_cell)) THEN
1097 DEALLOCATE (kpoint%index_to_cell)
1098 END IF
1099 ALLOCATE (kpoint%index_to_cell(3, ncount))
1100 index_to_cell => kpoint%index_to_cell
1101 DO ic = 1, ncount
1102 cell = minloc(cti)
1103 i1 = cell(1) - 1 - itm(1)
1104 i2 = cell(2) - 1 - itm(2)
1105 i3 = cell(3) - 1 - itm(3)
1106 cti(i1, i2, i3) = 1000000
1107 index_to_cell(1, ic) = i1
1108 index_to_cell(2, ic) = i2
1109 index_to_cell(3, ic) = i3
1110 END DO
1111 cti(:, :, :) = 0
1112 DO ic = 1, ncount
1113 i1 = index_to_cell(1, ic)
1114 i2 = index_to_cell(2, ic)
1115 i3 = index_to_cell(3, ic)
1116 cti(i1, i2, i3) = ic
1117 END DO
1118
1119 ! keep pointer to this neighborlist
1120 kpoint%sab_nl => sab_nl
1121
1122 ! set number of images
1123 nimages = SIZE(index_to_cell, 2)
1124
1125 DEALLOCATE (list)
1126
1127 CALL timestop(handle)
1128
1129 END SUBROUTINE kpoint_init_cell_index
1130
1131! **************************************************************************************************
1132!> \brief Transformation of real space matrices to a kpoint
1133!> \param rmatrix Real part of kpoint matrix
1134!> \param cmatrix Complex part of kpoint matrix (optional)
1135!> \param rsmat Real space matrices
1136!> \param ispin Spin index
1137!> \param xkp Kpoint coordinates
1138!> \param cell_to_index mapping of cell indices to RS index
1139!> \param sab_nl Defining neighbor list
1140!> \param is_complex Matrix to be transformed is imaginary
1141!> \param rs_sign Matrix to be transformed is csaled by rs_sign
1142! **************************************************************************************************
1143 SUBROUTINE rskp_transform(rmatrix, cmatrix, rsmat, ispin, &
1144 xkp, cell_to_index, sab_nl, is_complex, rs_sign)
1145
1146 TYPE(dbcsr_type) :: rmatrix
1147 TYPE(dbcsr_type), OPTIONAL :: cmatrix
1148 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
1149 INTEGER, INTENT(IN) :: ispin
1150 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: xkp
1151 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1152 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1153 POINTER :: sab_nl
1154 LOGICAL, INTENT(IN), OPTIONAL :: is_complex
1155 REAL(kind=dp), INTENT(IN), OPTIONAL :: rs_sign
1156
1157 CHARACTER(LEN=*), PARAMETER :: routinen = 'rskp_transform'
1158
1159 INTEGER :: handle, iatom, ic, icol, irow, jatom, &
1160 nimg
1161 INTEGER, DIMENSION(3) :: cell
1162 LOGICAL :: do_symmetric, found, my_complex, &
1163 wfn_real_only
1164 REAL(kind=dp) :: arg, coskl, fsign, fsym, sinkl
1165 REAL(kind=dp), DIMENSION(:, :), POINTER :: cblock, rblock, rsblock
1166 TYPE(neighbor_list_iterator_p_type), &
1167 DIMENSION(:), POINTER :: nl_iterator
1168
1169 CALL timeset(routinen, handle)
1170
1171 my_complex = .false.
1172 IF (PRESENT(is_complex)) my_complex = is_complex
1173
1174 fsign = 1.0_dp
1175 IF (PRESENT(rs_sign)) fsign = rs_sign
1176
1177 wfn_real_only = .true.
1178 IF (PRESENT(cmatrix)) wfn_real_only = .false.
1179
1180 nimg = SIZE(rsmat, 2)
1181
1182 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1183
1184 CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
1185 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1186 CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
1187
1188 ! fsym = +- 1 is due to real space matrices being non-symmetric (although in a symmtric type)
1189 ! with the link S_mu^0,nu^b = S_nu^0,mu^-b, and the KP matrices beeing Hermitian
1190 fsym = 1.0_dp
1191 irow = iatom
1192 icol = jatom
1193 IF (do_symmetric .AND. (iatom > jatom)) THEN
1194 irow = jatom
1195 icol = iatom
1196 fsym = -1.0_dp
1197 END IF
1198
1199 ic = cell_to_index(cell(1), cell(2), cell(3))
1200 IF (ic < 1 .OR. ic > nimg) cycle
1201
1202 arg = real(cell(1), dp)*xkp(1) + real(cell(2), dp)*xkp(2) + real(cell(3), dp)*xkp(3)
1203 IF (my_complex) THEN
1204 coskl = fsign*fsym*cos(twopi*arg)
1205 sinkl = fsign*sin(twopi*arg)
1206 ELSE
1207 coskl = fsign*cos(twopi*arg)
1208 sinkl = fsign*fsym*sin(twopi*arg)
1209 END IF
1210
1211 CALL dbcsr_get_block_p(matrix=rsmat(ispin, ic)%matrix, row=irow, col=icol, &
1212 block=rsblock, found=found)
1213 IF (.NOT. found) cycle
1214
1215 IF (wfn_real_only) THEN
1216 CALL dbcsr_get_block_p(matrix=rmatrix, row=irow, col=icol, &
1217 block=rblock, found=found)
1218 IF (.NOT. found) cycle
1219 rblock = rblock + coskl*rsblock
1220 ELSE
1221 CALL dbcsr_get_block_p(matrix=rmatrix, row=irow, col=icol, &
1222 block=rblock, found=found)
1223 IF (.NOT. found) cycle
1224 CALL dbcsr_get_block_p(matrix=cmatrix, row=irow, col=icol, &
1225 block=cblock, found=found)
1226 IF (.NOT. found) cycle
1227 rblock = rblock + coskl*rsblock
1228 cblock = cblock + sinkl*rsblock
1229 END IF
1230
1231 END DO
1232 CALL neighbor_list_iterator_release(nl_iterator)
1233
1234 CALL timestop(handle)
1235
1236 END SUBROUTINE rskp_transform
1237
1238! **************************************************************************************************
1239!> \brief Prepare a batched real-cell to complete reciprocal-grid transform for local DBCSR blocks.
1240!> The stored blocks remain MPI local. Irregular or too-large grids return used_fft=.FALSE.
1241!> \param grid cached transformed local blocks and their DBCSR coordinates
1242!> \param rmatrix allocated output template defining the local block layout
1243!> \param rsmat real-space matrix set
1244!> \param ispin spin component of rsmat
1245!> \param xkp complete reciprocal grid in arbitrary order
1246!> \param nkp_grid reciprocal grid dimensions
1247!> \param cell_to_index real-cell coordinate mapping
1248!> \param sab_nl neighbor list defining the stored block orientation
1249!> \param used_fft whether the cache was prepared
1250!> \param is_complex whether the real-space operator is imaginary
1251!> \param rs_sign optional overall sign
1252!> \param max_storage_bytes optional conservative per-rank memory limit
1253! **************************************************************************************************
1254 SUBROUTINE rskp_transform_grid_prepare(grid, rmatrix, rsmat, ispin, xkp, nkp_grid, &
1255 cell_to_index, sab_nl, used_fft, is_complex, rs_sign, &
1256 max_storage_bytes)
1257
1258 TYPE(rskp_grid_type), INTENT(INOUT) :: grid
1259 TYPE(dbcsr_type), INTENT(INOUT) :: rmatrix
1260 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
1261 INTEGER, INTENT(IN) :: ispin
1262 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp
1263 INTEGER, DIMENSION(3), INTENT(IN) :: nkp_grid
1264 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1265 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1266 POINTER :: sab_nl
1267 LOGICAL, INTENT(OUT) :: used_fft
1268 LOGICAL, INTENT(IN), OPTIONAL :: is_complex
1269 REAL(kind=dp), INTENT(IN), OPTIONAL :: rs_sign
1270 INTEGER(KIND=int_8), INTENT(IN), OPTIONAL :: max_storage_bytes
1271
1272 CHARACTER(LEN=*), PARAMETER :: routinen = 'rskp_transform_grid_prepare'
1273
1274 INTEGER :: handle, i1, i2, i3, iatom, iblock, ic, &
1275 icol, irow, jatom, nblkcols, nblkrows, &
1276 nblocks, ncell, nelem, nimg, nkp, &
1277 nvalues
1278 INTEGER(KIND=int_8) :: memory_limit, storage_bytes
1279 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: block_map, index_to_cell
1280 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: fft_cell_to_index
1281 INTEGER, DIMENSION(3) :: cell, fft_cell
1282 LOGICAL :: do_symmetric, found, my_complex
1283 REAL(kind=dp) :: fsign, fsym, value_sign
1284 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: values_rs
1285 REAL(kind=dp), DIMENSION(:, :), POINTER :: rblock, rsblock
1286 TYPE(dbcsr_iterator_type) :: iter
1287 TYPE(neighbor_list_iterator_p_type), &
1288 DIMENSION(:), POINTER :: nl_iterator
1289
1290 CALL timeset(routinen, handle)
1291 CALL rskp_transform_grid_release(grid)
1292 used_fft = .false.
1293
1294 nkp = SIZE(xkp, 2)
1295 nimg = SIZE(rsmat, 2)
1296 IF (nkp < 2 .OR. product(nkp_grid) /= nkp) THEN
1297 CALL timestop(handle)
1298 RETURN
1299 END IF
1300 IF (.NOT. regular_kpoint_grid(xkp, nkp_grid)) THEN
1301 CALL timestop(handle)
1302 RETURN
1303 END IF
1304
1305 nblocks = 0
1306 nvalues = 0
1307 CALL dbcsr_iterator_start(iter, rmatrix)
1308 DO WHILE (dbcsr_iterator_blocks_left(iter))
1309 CALL dbcsr_iterator_next_block(iter, irow, icol, rblock)
1310 nblocks = nblocks + 1
1311 nvalues = nvalues + SIZE(rblock)
1312 END DO
1313 CALL dbcsr_iterator_stop(iter)
1314 IF (nblocks == 0 .OR. nvalues == 0) THEN
1315 CALL timestop(handle)
1316 RETURN
1317 END IF
1318
1319 CALL dbcsr_get_info(rmatrix, nblkrows_total=nblkrows, nblkcols_total=nblkcols)
1320 ALLOCATE (block_map(nblkrows, nblkcols), source=0)
1321 ALLOCATE (grid%block_row(nblocks), grid%block_col(nblocks), &
1322 grid%block_offset(nblocks), grid%block_nelem(nblocks))
1323 iblock = 0
1324 nelem = 0
1325 CALL dbcsr_iterator_start(iter, rmatrix)
1326 DO WHILE (dbcsr_iterator_blocks_left(iter))
1327 CALL dbcsr_iterator_next_block(iter, irow, icol, rblock)
1328 iblock = iblock + 1
1329 grid%block_row(iblock) = irow
1330 grid%block_col(iblock) = icol
1331 grid%block_offset(iblock) = nelem + 1
1332 grid%block_nelem(iblock) = SIZE(rblock)
1333 block_map(irow, icol) = iblock
1334 nelem = nelem + SIZE(rblock)
1335 END DO
1336 CALL dbcsr_iterator_stop(iter)
1337
1338 ALLOCATE (fft_cell_to_index(lbound(cell_to_index, 1):ubound(cell_to_index, 1), &
1339 lbound(cell_to_index, 2):ubound(cell_to_index, 2), &
1340 lbound(cell_to_index, 3):ubound(cell_to_index, 3)), source=0)
1341 ncell = 0
1342 DO i3 = lbound(cell_to_index, 3), ubound(cell_to_index, 3)
1343 DO i2 = lbound(cell_to_index, 2), ubound(cell_to_index, 2)
1344 DO i1 = lbound(cell_to_index, 1), ubound(cell_to_index, 1)
1345 IF (cell_to_index(i1, i2, i3) < 1 .OR. cell_to_index(i1, i2, i3) > nimg) cycle
1346 IF (fft_cell_to_index(i1, i2, i3) == 0) THEN
1347 ncell = ncell + 1
1348 fft_cell_to_index(i1, i2, i3) = ncell
1349 END IF
1350 IF (fft_cell_to_index(-i1, -i2, -i3) == 0) THEN
1351 ncell = ncell + 1
1352 fft_cell_to_index(-i1, -i2, -i3) = ncell
1353 END IF
1354 END DO
1355 END DO
1356 END DO
1357 IF (ncell == 0) THEN
1358 CALL rskp_transform_grid_release(grid)
1359 CALL timestop(handle)
1360 RETURN
1361 END IF
1362
1363 memory_limit = 512_int_8*1024_int_8**2
1364 IF (PRESENT(max_storage_bytes)) memory_limit = max_storage_bytes
1365 ! Include the real-cell input, complex result and conservative FFT work/padding allowance.
1366 storage_bytes = int(nvalues, int_8)*(8_int_8*int(ncell, int_8) + &
1367 80_int_8*int(nkp, int_8))
1368 IF (storage_bytes > memory_limit) THEN
1369 CALL rskp_transform_grid_release(grid)
1370 CALL timestop(handle)
1371 RETURN
1372 END IF
1373
1374 ALLOCATE (index_to_cell(3, ncell), source=0)
1375 DO i3 = lbound(fft_cell_to_index, 3), ubound(fft_cell_to_index, 3)
1376 DO i2 = lbound(fft_cell_to_index, 2), ubound(fft_cell_to_index, 2)
1377 DO i1 = lbound(fft_cell_to_index, 1), ubound(fft_cell_to_index, 1)
1378 ic = fft_cell_to_index(i1, i2, i3)
1379 IF (ic > 0) index_to_cell(:, ic) = [i1, i2, i3]
1380 END DO
1381 END DO
1382 END DO
1383
1384 ALLOCATE (values_rs(nvalues, 1, ncell), source=0.0_dp)
1385 ALLOCATE (grid%values(nvalues, 1, nkp))
1386 my_complex = .false.
1387 IF (PRESENT(is_complex)) my_complex = is_complex
1388 fsign = 1.0_dp
1389 IF (PRESENT(rs_sign)) fsign = rs_sign
1390
1391 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1392 CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
1393 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1394 CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
1395
1396 fsym = 1.0_dp
1397 irow = iatom
1398 icol = jatom
1399 IF (do_symmetric .AND. iatom > jatom) THEN
1400 irow = jatom
1401 icol = iatom
1402 fsym = -1.0_dp
1403 END IF
1404 IF (irow < 1 .OR. irow > nblkrows .OR. icol < 1 .OR. icol > nblkcols) cycle
1405 iblock = block_map(irow, icol)
1406 IF (iblock == 0) cycle
1407
1408 ic = cell_to_index(cell(1), cell(2), cell(3))
1409 IF (ic < 1 .OR. ic > nimg) cycle
1410 CALL dbcsr_get_block_p(matrix=rsmat(ispin, ic)%matrix, row=irow, col=icol, &
1411 block=rsblock, found=found)
1412 IF (.NOT. found) cycle
1413
1414 fft_cell = cell
1415 value_sign = fsign
1416 IF (fsym < 0.0_dp) fft_cell = -cell
1417 IF (my_complex) value_sign = value_sign*fsym
1418 ic = fft_cell_to_index(fft_cell(1), fft_cell(2), fft_cell(3))
1419 nelem = grid%block_nelem(iblock)
1420 i1 = grid%block_offset(iblock)
1421 values_rs(i1:i1 + nelem - 1, 1, ic) = values_rs(i1:i1 + nelem - 1, 1, ic) + &
1422 value_sign*reshape(rsblock, [nelem])
1423 END DO
1424 CALL neighbor_list_iterator_release(nl_iterator)
1425
1426 CALL cell_to_k_grid_fft(values_rs, index_to_cell, xkp, nkp_grid, grid%values, used_fft)
1427 IF (used_fft) THEN
1428 grid%ready = .true.
1429 ELSE
1430 CALL rskp_transform_grid_release(grid)
1431 END IF
1432
1433 CALL timestop(handle)
1434
1435 END SUBROUTINE rskp_transform_grid_prepare
1436
1437! **************************************************************************************************
1438!> \brief Extract one reciprocal-grid matrix from a prepared local DBCSR block cache.
1439!> \param grid ...
1440!> \param ikp ...
1441!> \param rmatrix ...
1442!> \param cmatrix ...
1443! **************************************************************************************************
1444 SUBROUTINE rskp_transform_grid_extract(grid, ikp, rmatrix, cmatrix)
1445
1446 TYPE(rskp_grid_type), INTENT(IN) :: grid
1447 INTEGER, INTENT(IN) :: ikp
1448 TYPE(dbcsr_type), INTENT(INOUT) :: rmatrix
1449 TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL :: cmatrix
1450
1451 INTEGER :: iblock, ioffset, nelem
1452 LOGICAL :: found
1453 REAL(kind=dp), DIMENSION(:, :), POINTER :: cblock, rblock
1454
1455 cpassert(grid%ready)
1456 cpassert(ikp >= 1 .AND. ikp <= SIZE(grid%values, 3))
1457 CALL dbcsr_set(rmatrix, 0.0_dp)
1458 IF (PRESENT(cmatrix)) CALL dbcsr_set(cmatrix, 0.0_dp)
1459
1460 DO iblock = 1, SIZE(grid%block_row)
1461 ioffset = grid%block_offset(iblock)
1462 nelem = grid%block_nelem(iblock)
1463 CALL dbcsr_get_block_p(rmatrix, grid%block_row(iblock), grid%block_col(iblock), &
1464 rblock, found=found)
1465 cpassert(found)
1466 rblock = reshape(real(grid%values(ioffset:ioffset + nelem - 1, 1, ikp), kind=dp), &
1467 shape(rblock))
1468 IF (PRESENT(cmatrix)) THEN
1469 CALL dbcsr_get_block_p(cmatrix, grid%block_row(iblock), grid%block_col(iblock), &
1470 cblock, found=found)
1471 cpassert(found)
1472 cblock = reshape(aimag(grid%values(ioffset:ioffset + nelem - 1, 1, ikp)), shape(cblock))
1473 END IF
1474 END DO
1475
1476 END SUBROUTINE rskp_transform_grid_extract
1477
1478! **************************************************************************************************
1479!> \brief Release a reciprocal-grid transformation cache.
1480!> \param grid ...
1481! **************************************************************************************************
1482 SUBROUTINE rskp_transform_grid_release(grid)
1483
1484 TYPE(rskp_grid_type), INTENT(INOUT) :: grid
1485
1486 IF (ALLOCATED(grid%values)) DEALLOCATE (grid%values)
1487 IF (ALLOCATED(grid%block_row)) DEALLOCATE (grid%block_row)
1488 IF (ALLOCATED(grid%block_col)) DEALLOCATE (grid%block_col)
1489 IF (ALLOCATED(grid%block_offset)) DEALLOCATE (grid%block_offset)
1490 IF (ALLOCATED(grid%block_nelem)) DEALLOCATE (grid%block_nelem)
1491 grid%ready = .false.
1492
1493 END SUBROUTINE rskp_transform_grid_release
1494
1495! **************************************************************************************************
1496!> \brief Given the eigenvalues of all kpoints, calculates the occupation numbers
1497!> \param kpoint Kpoint environment
1498!> \param smear Smearing information
1499!> \param probe ...
1500!> \param added_mos_auto ...
1501!> \param added_mos_auto_grow ...
1502!> \param separate_spin_occupations enforce the electron count of each spin channel
1503! **************************************************************************************************
1504 SUBROUTINE kpoint_set_mo_occupation( &
1505 kpoint, smear, probe, added_mos_auto, added_mos_auto_grow, separate_spin_occupations)
1506
1507 TYPE(kpoint_type), POINTER :: kpoint
1508 TYPE(smear_type) :: smear
1509 TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
1510 POINTER :: probe
1511 LOGICAL, DIMENSION(:), INTENT(IN), OPTIONAL :: added_mos_auto
1512 LOGICAL, INTENT(OUT), OPTIONAL :: added_mos_auto_grow
1513 LOGICAL, INTENT(IN), OPTIONAL :: separate_spin_occupations
1514
1515 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_set_mo_occupation'
1516
1517 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: csmatrix
1518 INTEGER :: handle, ik, ikpgr, ispin, kplocal, nao, &
1519 nb, ncol_global, ne_a, ne_b, &
1520 nelectron, nkp, nmo, nrow_global, nspin
1521 INTEGER, DIMENSION(2) :: kp_range
1522 LOGICAL :: my_added_mos_auto_grow, &
1523 my_separate_spin_occupations
1524 REAL(kind=dp) :: kts, kts_spin(2), mu, mus(2), nel
1525 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: smatrix
1526 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: weig, wocc
1527 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: icoeff, rcoeff
1528 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues, occupation, wkp
1529 TYPE(cp_cfm_type), POINTER :: cmo_coeff
1530 TYPE(cp_fm_type), POINTER :: mo_coeff
1531 TYPE(kpoint_env_type), POINTER :: kp
1532 TYPE(mo_set_type), POINTER :: mo_set
1533 TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
1534
1535 CALL timeset(routinen, handle)
1536
1537 my_added_mos_auto_grow = .false.
1538 IF (PRESENT(added_mos_auto_grow)) added_mos_auto_grow = .false.
1539
1540 ! first collect all the eigenvalues
1541 CALL get_kpoint_info(kpoint, nkp=nkp)
1542 kp => kpoint%kp_env(1)%kpoint_env
1543 nspin = SIZE(kp%mos)
1544 my_separate_spin_occupations = smear%fixed_mag_mom > 0.0_dp
1545 IF (PRESENT(separate_spin_occupations)) THEN
1546 my_separate_spin_occupations = my_separate_spin_occupations .OR. separate_spin_occupations
1547 END IF
1548 mo_set => kp%mos(1)
1549 CALL get_mo_set(mo_set, nmo=nmo, nao=nao, nelectron=nelectron)
1550 ne_a = nelectron
1551 IF (nspin == 2) THEN
1552 CALL get_mo_set(kp%mos(2), nmo=nb, nelectron=ne_b)
1553 cpassert(nmo == nb)
1554 END IF
1555 ALLOCATE (weig(nmo, nkp, nspin), wocc(nmo, nkp, nspin))
1556 weig = 0.0_dp
1557 wocc = 0.0_dp
1558 IF (PRESENT(probe)) THEN
1559 ALLOCATE (rcoeff(nao, nmo, nkp, nspin), icoeff(nao, nmo, nkp, nspin))
1560 rcoeff = 0.0_dp !coeff, real part
1561 icoeff = 0.0_dp !coeff, imaginary part
1562 END IF
1563 CALL get_kpoint_info(kpoint, kp_range=kp_range)
1564 kplocal = kp_range(2) - kp_range(1) + 1
1565 DO ikpgr = 1, kplocal
1566 ik = kp_range(1) + ikpgr - 1
1567 kp => kpoint%kp_env(ikpgr)%kpoint_env
1568 DO ispin = 1, nspin
1569 mo_set => kp%mos(ispin)
1570 CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
1571 weig(1:nmo, ik, ispin) = eigenvalues(1:nmo)
1572 IF (PRESENT(probe)) THEN
1573 IF (kpoint%use_real_wfn) THEN
1574 CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
1575 CALL cp_fm_get_info(mo_coeff, nrow_global=nrow_global, ncol_global=ncol_global)
1576 ALLOCATE (smatrix(nrow_global, ncol_global))
1577 CALL cp_fm_get_submatrix(mo_coeff, smatrix)
1578 rcoeff(1:nao, 1:nmo, ik, ispin) = smatrix(1:nao, 1:nmo)
1579 DEALLOCATE (smatrix)
1580 ELSE
1581 CALL get_mo_set(mo_set, cmo_coeff=cmo_coeff)
1582 CALL cp_cfm_get_info(cmo_coeff, nrow_global=nrow_global, ncol_global=ncol_global)
1583 ALLOCATE (csmatrix(nrow_global, ncol_global))
1584 CALL cp_cfm_get_submatrix(cmo_coeff, csmatrix)
1585 rcoeff(1:nao, 1:nmo, ik, ispin) = real(csmatrix(1:nao, 1:nmo), kind=dp)
1586 icoeff(1:nao, 1:nmo, ik, ispin) = aimag(csmatrix(1:nao, 1:nmo))
1587 DEALLOCATE (csmatrix)
1588 END IF
1589 END IF
1590 END DO
1591 END DO
1592 CALL get_kpoint_info(kpoint, para_env_inter_kp=para_env_inter_kp)
1593 CALL para_env_inter_kp%sum(weig)
1594
1595 IF (PRESENT(probe)) THEN
1596 CALL para_env_inter_kp%sum(rcoeff)
1597 CALL para_env_inter_kp%sum(icoeff)
1598 END IF
1599
1600 CALL get_kpoint_info(kpoint, wkp=wkp)
1601 kts_spin = 0.0_dp
1602
1603!calling of HP module HERE, before smear
1604 IF (PRESENT(probe)) THEN
1605 smear%do_smear = .false. !ensures smearing is switched off
1606
1607 IF (nspin == 1) THEN
1608 nel = real(nelectron, kind=dp)
1609 CALL probe_occupancy_kp(wocc(:, :, :), mus(1), kts, weig(:, :, :), rcoeff(:, :, :, :), icoeff(:, :, :, :), 2.0d0, &
1610 probe, nel, wkp)
1611 ELSE
1612 nel = real(ne_a, kind=dp) + real(ne_b, kind=dp)
1613 CALL probe_occupancy_kp(wocc(:, :, :), mu, kts, weig(:, :, :), rcoeff(:, :, :, :), icoeff(:, :, :, :), 1.0d0, &
1614 probe, nel, wkp)
1615 kts = kts/2._dp
1616 mus(1:2) = mu
1617 END IF
1618
1619 DO ikpgr = 1, kplocal
1620 ik = kp_range(1) + ikpgr - 1
1621 kp => kpoint%kp_env(ikpgr)%kpoint_env
1622 DO ispin = 1, nspin
1623 mo_set => kp%mos(ispin)
1624 CALL get_mo_set(mo_set, eigenvalues=eigenvalues, occupation_numbers=occupation)
1625 eigenvalues(1:nmo) = weig(1:nmo, ik, ispin)
1626 occupation(1:nmo) = wocc(1:nmo, ik, ispin)
1627 mo_set%kTS = kts
1628 mo_set%mu = mus(ispin)
1629
1630 END DO
1631 END DO
1632
1633 DEALLOCATE (weig, wocc, rcoeff, icoeff)
1634
1635 END IF
1636
1637 IF (PRESENT(probe) .EQV. .false.) THEN
1638 IF (smear%do_smear) THEN
1639 SELECT CASE (smear%method)
1640 CASE (smear_fermi_dirac)
1641 ! finite electronic temperature
1642 IF (nspin == 1) THEN
1643 nel = real(nelectron, kind=dp)
1644 CALL smearkp(wocc(:, :, 1), mus(1), kts, weig(:, :, 1), nel, wkp, &
1645 smear%electronic_temperature, 2.0_dp, smear_fermi_dirac)
1646 kts_spin(1) = kts
1647 ELSE IF (my_separate_spin_occupations) THEN
1648 nel = real(ne_a, kind=dp)
1649 CALL smearkp(wocc(:, :, 1), mus(1), kts, weig(:, :, 1), nel, wkp, &
1650 smear%electronic_temperature, 1.0_dp, smear_fermi_dirac)
1651 kts_spin(1) = kts
1652 nel = real(ne_b, kind=dp)
1653 CALL smearkp(wocc(:, :, 2), mus(2), kts, weig(:, :, 2), nel, wkp, &
1654 smear%electronic_temperature, 1.0_dp, smear_fermi_dirac)
1655 kts_spin(2) = kts
1656 ELSE
1657 nel = real(ne_a, kind=dp) + real(ne_b, kind=dp)
1658 CALL smearkp2(wocc(:, :, :), mu, kts, weig(:, :, :), nel, wkp, &
1659 smear%electronic_temperature, smear_fermi_dirac)
1660 kts = kts/2._dp
1661 kts_spin(1:2) = kts
1662 mus(1:2) = mu
1663 END IF
1664 CASE (smear_gaussian, smear_mp, smear_mv)
1665 IF (nspin == 1) THEN
1666 nel = real(nelectron, kind=dp)
1667 CALL smearkp(wocc(:, :, 1), mus(1), kts, weig(:, :, 1), nel, wkp, &
1668 smear%smearing_width, 2.0_dp, smear%method)
1669 kts_spin(1) = kts
1670 ELSE IF (my_separate_spin_occupations) THEN
1671 nel = real(ne_a, kind=dp)
1672 CALL smearkp(wocc(:, :, 1), mus(1), kts, weig(:, :, 1), nel, wkp, &
1673 smear%smearing_width, 1.0_dp, smear%method)
1674 kts_spin(1) = kts
1675 nel = real(ne_b, kind=dp)
1676 CALL smearkp(wocc(:, :, 2), mus(2), kts, weig(:, :, 2), nel, wkp, &
1677 smear%smearing_width, 1.0_dp, smear%method)
1678 kts_spin(2) = kts
1679 ELSE
1680 nel = real(ne_a, kind=dp) + real(ne_b, kind=dp)
1681 CALL smearkp2(wocc(:, :, :), mu, kts, weig(:, :, :), nel, wkp, &
1682 smear%smearing_width, smear%method)
1683 kts = kts/2._dp
1684 kts_spin(1:2) = kts
1685 mus(1:2) = mu
1686 END IF
1687 CASE DEFAULT
1688 cpabort("kpoints: Selected smearing not (yet) supported")
1689 END SELECT
1690 ELSE
1691 ! fixed occupations (2/1)
1692 IF (nspin == 1) THEN
1693 nel = real(nelectron, kind=dp)
1694 CALL smearkp(wocc(:, :, 1), mus(1), kts, weig(:, :, 1), nel, wkp, &
1695 0.0_dp, 2.0_dp, smear_gaussian)
1696 kts_spin(1) = kts
1697 ELSE
1698 nel = real(ne_a, kind=dp)
1699 CALL smearkp(wocc(:, :, 1), mus(1), kts, weig(:, :, 1), nel, wkp, &
1700 0.0_dp, 1.0_dp, smear_gaussian)
1701 kts_spin(1) = kts
1702 nel = real(ne_b, kind=dp)
1703 CALL smearkp(wocc(:, :, 2), mus(2), kts, weig(:, :, 2), nel, wkp, &
1704 0.0_dp, 1.0_dp, smear_gaussian)
1705 kts_spin(2) = kts
1706 END IF
1707 END IF
1708 IF (smear%do_smear .AND. PRESENT(added_mos_auto)) THEN
1709 CALL kpoint_check_added_mos_auto_occupation( &
1710 wocc, wkp, smear, nspin, added_mos_auto, nao, my_added_mos_auto_grow)
1711 END IF
1712 IF (my_added_mos_auto_grow) THEN
1713 IF (PRESENT(added_mos_auto_grow)) THEN
1714 added_mos_auto_grow = .true.
1715 DEALLOCATE (weig, wocc)
1716 CALL timestop(handle)
1717 RETURN
1718 ELSE
1719 CALL cp_abort(__location__, &
1720 "K-point ADDED_MOS AUTO needs a larger virtual-space buffer, "// &
1721 "but this call path cannot grow it.")
1722 END IF
1723 END IF
1724 DO ikpgr = 1, kplocal
1725 ik = kp_range(1) + ikpgr - 1
1726 kp => kpoint%kp_env(ikpgr)%kpoint_env
1727 DO ispin = 1, nspin
1728 mo_set => kp%mos(ispin)
1729 CALL get_mo_set(mo_set, eigenvalues=eigenvalues, occupation_numbers=occupation)
1730 eigenvalues(1:nmo) = weig(1:nmo, ik, ispin)
1731 occupation(1:nmo) = wocc(1:nmo, ik, ispin)
1732 mo_set%kTS = kts_spin(ispin)
1733 mo_set%mu = mus(ispin)
1734 IF (.NOT. smear%do_smear) THEN
1735 ! Zero-width filling can give different occupied ranks at different k points.
1736 kp%mos(ispin)%homo = count(occupation(1:nmo) > 0.0_dp)
1737 END IF
1738 END DO
1739 END DO
1740
1741 DEALLOCATE (weig, wocc)
1742
1743 END IF
1744
1745 CALL timestop(handle)
1746
1747 END SUBROUTINE kpoint_set_mo_occupation
1748
1749! **************************************************************************************************
1750!> \brief summarize whether weighted k-point smearing reaches the available band edges
1751!> \param wocc ...
1752!> \param wkp ...
1753!> \param smear ...
1754!> \param nspin ...
1755!> \param has_weight ...
1756!> \param first_fractional ...
1757!> \param last_occupied ...
1758!> \param last_occupied_spin ...
1759! **************************************************************************************************
1760 SUBROUTINE kpoint_smearing_edge_status(wocc, wkp, smear, nspin, has_weight, &
1761 first_fractional, last_occupied, last_occupied_spin)
1762
1763 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: wocc
1764 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wkp
1765 TYPE(smear_type), INTENT(IN) :: smear
1766 INTEGER, INTENT(IN) :: nspin
1767 LOGICAL, INTENT(OUT) :: has_weight, first_fractional, &
1768 last_occupied
1769 LOGICAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: last_occupied_spin
1770
1771 INTEGER :: ik, ispin, nmo
1772 LOGICAL :: band_occupied
1773 REAL(kind=dp) :: eps_occ, maxocc, weight, weight_threshold
1774
1775 has_weight = .false.
1776 first_fractional = .false.
1777 last_occupied = .false.
1778 IF (PRESENT(last_occupied_spin)) last_occupied_spin(:) = .false.
1779 IF (.NOT. smear%do_smear) RETURN
1780
1781 nmo = SIZE(wocc, 1)
1782 IF (nmo < 1) RETURN
1783 eps_occ = max(smear%eps_fermi_dirac, 10.0_dp*epsilon(1.0_dp))
1784 weight_threshold = 10.0_dp*epsilon(1.0_dp)
1785
1786 DO ispin = 1, min(nspin, SIZE(wocc, 3))
1787 maxocc = merge(2.0_dp, 1.0_dp, nspin == 1)
1788 DO ik = 1, min(SIZE(wkp), SIZE(wocc, 2))
1789 weight = abs(wkp(ik))
1790 IF (weight <= weight_threshold) cycle
1791 has_weight = .true.
1792 first_fractional = first_fractional .OR. &
1793 abs(wocc(1, ik, ispin) - maxocc) > eps_occ
1794 band_occupied = abs(wocc(nmo, ik, ispin)) > eps_occ
1795 last_occupied = last_occupied .OR. band_occupied
1796 IF (band_occupied .AND. PRESENT(last_occupied_spin)) THEN
1797 IF (ispin <= SIZE(last_occupied_spin)) last_occupied_spin(ispin) = .true.
1798 END IF
1799 END DO
1800 END DO
1801
1802 END SUBROUTINE kpoint_smearing_edge_status
1803
1804! **************************************************************************************************
1805!> \brief request ADDED_MOS AUTO growth if weighted k-point smearing occupies the last band
1806!> \param wocc ...
1807!> \param wkp ...
1808!> \param smear ...
1809!> \param nspin ...
1810!> \param added_mos_auto ...
1811!> \param nao ...
1812!> \param added_mos_auto_grow ...
1813! **************************************************************************************************
1814 SUBROUTINE kpoint_check_added_mos_auto_occupation(wocc, wkp, smear, nspin, &
1815 added_mos_auto, nao, added_mos_auto_grow)
1816
1817 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: wocc
1818 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wkp
1819 TYPE(smear_type), INTENT(IN) :: smear
1820 INTEGER, INTENT(IN) :: nspin
1821 LOGICAL, DIMENSION(:), INTENT(IN) :: added_mos_auto
1822 INTEGER, INTENT(IN) :: nao
1823 LOGICAL, INTENT(OUT) :: added_mos_auto_grow
1824
1825 INTEGER :: ispin, nmo
1826 LOGICAL :: first_fractional, has_weight, &
1827 last_occupied
1828 LOGICAL, DIMENSION(nspin) :: last_occupied_spin
1829
1830 added_mos_auto_grow = .false.
1831 IF (.NOT. any(added_mos_auto) .OR. .NOT. smear%do_smear) RETURN
1832
1833 CALL kpoint_smearing_edge_status(wocc, wkp, smear, nspin, has_weight, &
1834 first_fractional, last_occupied, &
1835 last_occupied_spin=last_occupied_spin)
1836 IF (.NOT. has_weight) RETURN
1837
1838 nmo = SIZE(wocc, 1)
1839 DO ispin = 1, min(nspin, SIZE(added_mos_auto))
1840 IF (.NOT. added_mos_auto(ispin) .OR. .NOT. last_occupied_spin(ispin)) cycle
1841 IF (nmo >= nao) THEN
1842 CALL cp_abort(__location__, &
1843 "K-point ADDED_MOS AUTO exhausted the AO basis but the highest band "// &
1844 "is still occupied. Use a larger basis or reduce the smearing width.")
1845 END IF
1846 added_mos_auto_grow = .true.
1847 RETURN
1848 END DO
1849
1850 END SUBROUTINE kpoint_check_added_mos_auto_occupation
1851
1852! **************************************************************************************************
1853!> \brief Calculate kpoint density matrices (rho(k), owned by kpoint groups)
1854!> \param kpoint kpoint environment
1855!> \param energy_weighted calculate energy weighted density matrix
1856!> \param for_aux_fit ...
1857! **************************************************************************************************
1858 SUBROUTINE kpoint_density_matrices(kpoint, energy_weighted, for_aux_fit)
1859
1860 TYPE(kpoint_type), POINTER :: kpoint
1861 LOGICAL, OPTIONAL :: energy_weighted, for_aux_fit
1862
1863 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_density_matrices'
1864
1865 INTEGER :: handle, ikpgr, ispin, kplocal, nao, nmo, &
1866 nspin, nworkers, thread
1867 LOGICAL :: aux_fit, local, wtype
1868 TYPE(cp_cfm_type), ALLOCATABLE :: cdensity(:), cwork(:)
1869 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1870 TYPE(cp_fm_type), ALLOCATABLE :: work(:)
1871 TYPE(kpoint_env_type), POINTER :: kp
1872 TYPE(local_gemm_ctxt_type), ALLOCATABLE :: gemm_ctx(:)
1873
1874 CALL timeset(routinen, handle)
1875 wtype = .false.
1876 IF (PRESENT(energy_weighted)) wtype = energy_weighted
1877 aux_fit = .false.
1878 IF (PRESENT(for_aux_fit)) aux_fit = for_aux_fit
1879 IF (aux_fit) THEN
1880 cpassert(ASSOCIATED(kpoint%kp_aux_env))
1881 kp => kpoint%kp_aux_env(1)%kpoint_env
1882 ELSE
1883 kp => kpoint%kp_env(1)%kpoint_env
1884 END IF
1885 CALL get_mo_set(kp%mos(1), nao=nao, nmo=nmo)
1886 IF (kpoint%use_real_wfn) THEN
1887 CALL cp_fm_get_info(kp%mos(1)%mo_coeff, matrix_struct=matrix_struct)
1888 ELSE
1889 CALL cp_cfm_get_info(kp%mos(1)%cmo_coeff, matrix_struct=matrix_struct)
1890 END IF
1891 kplocal = kpoint%kp_range(2) - kpoint%kp_range(1) + 1
1892 nspin = SIZE(kp%mos)
1893 ! Only the group layout matters. Distributed GEMM stays on the calling thread.
1894 local = product(matrix_struct%context%num_pe) == 1
1895 nworkers = 1
1896!$ IF (local) nworkers = MIN(omp_get_max_threads(), kplocal*nspin)
1897 ALLOCATE (gemm_ctx(nworkers))
1898 IF (kpoint%use_real_wfn) THEN
1899 ALLOCATE (work(nworkers))
1900 ELSE
1901 ALLOCATE (cwork(nworkers), cdensity(nworkers))
1902 END IF
1903 DO thread = 1, nworkers
1904 IF (kpoint%use_real_wfn) THEN
1905 CALL cp_fm_create(work(thread), matrix_struct)
1906 ELSE
1907 CALL cp_cfm_create(cwork(thread), matrix_struct)
1908 CALL cp_cfm_create(cdensity(thread), kp%pmat(1, 1)%matrix_struct)
1909 END IF
1910 IF (local) CALL gemm_ctx(thread)%create(local_gemm_pu_gpu, timing=.false.)
1911 END DO
1912
1913!$OMP PARALLEL DO COLLAPSE(2) DEFAULT(NONE) SCHEDULE(DYNAMIC, 1) NUM_THREADS(nworkers) IF(nworkers > 1) &
1914!$OMP SHARED(kpoint,wtype,aux_fit,nao,nmo,kplocal,nspin,nworkers,work,cwork,cdensity,local,gemm_ctx) PRIVATE(ikpgr,ispin,thread)
1915 DO ikpgr = 1, kplocal
1916 DO ispin = 1, nspin
1917 thread = 1
1918!$ thread = omp_get_thread_num() + 1
1919 IF (kpoint%use_real_wfn) THEN
1920 CALL kpoint_density_matrix_job(kpoint, ikpgr, ispin, wtype, aux_fit, nao, nmo, &
1921 local, gemm_ctx(thread), work=work(thread))
1922 ELSE
1923 CALL kpoint_density_matrix_job(kpoint, ikpgr, ispin, wtype, aux_fit, nao, nmo, &
1924 local, gemm_ctx(thread), cwork=cwork(thread), cdensity=cdensity(thread))
1925 END IF
1926 END DO
1927 END DO
1928!$OMP END PARALLEL DO
1929 DO thread = 1, nworkers
1930 CALL gemm_ctx(thread)%destroy()
1931 IF (kpoint%use_real_wfn) THEN
1932 CALL cp_fm_release(work(thread))
1933 ELSE
1934 CALL cp_cfm_release(cwork(thread))
1935 CALL cp_cfm_release(cdensity(thread))
1936 END IF
1937 END DO
1938 IF (.NOT. wtype .AND. .NOT. aux_fit) THEN
1939 kpoint%lowdin_density_ready = .true.
1940 kpoint%lowdin_population_ready = .false.
1941 END IF
1942 CALL timestop(handle)
1943
1944 END SUBROUTINE kpoint_density_matrices
1945
1946! **************************************************************************************************
1947!> \brief Build W(k) for noncanonical complex OT orbitals from H(k) C(k).
1948!> The occupied-space Lagrange multiplier is
1949!> Hermitian[(C^H H C) f], which reduces to the usual eigenvalue-weighted
1950!> expression for canonical orbitals.
1951!> \param coeff C(k)
1952!> \param hc H(k) C(k)
1953!> \param occupation orbital occupations
1954!> \param wmat W(k)
1955! **************************************************************************************************
1956 SUBROUTINE kpoint_ot_energy_weighted_density(coeff, hc, occupation, wmat)
1957
1958 TYPE(cp_cfm_type), INTENT(IN) :: coeff, hc
1959 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: occupation
1960 TYPE(cp_cfm_type), INTENT(INOUT) :: wmat
1961
1962 COMPLEX(KIND=dp), PARAMETER :: zhalf = (0.5_dp, 0.0_dp), &
1963 zone = (1.0_dp, 0.0_dp), &
1964 zzero = (0.0_dp, 0.0_dp)
1965
1966 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: occupation_complex
1967 INTEGER :: nao, nmo
1968 TYPE(cp_cfm_type) :: hblock, hblock_h, weighted
1969
1970 CALL cp_cfm_get_info(coeff, nrow_global=nao, ncol_global=nmo)
1971 cpassert(nmo >= 1)
1972 cpassert(SIZE(occupation) >= nmo)
1973
1974 CALL cp_cfm_create(weighted, coeff%matrix_struct)
1975 CALL cp_cfm_create(hblock, coeff%matrix_struct, nrow=nmo, ncol=nmo)
1976 CALL cp_cfm_create(hblock_h, coeff%matrix_struct, nrow=nmo, ncol=nmo)
1977 CALL parallel_gemm("C", "N", nmo, nmo, nao, zone, coeff, hc, zzero, hblock)
1978
1979 ALLOCATE (occupation_complex(nmo))
1980 occupation_complex(:) = cmplx(occupation(1:nmo), 0.0_dp, kind=dp)
1981 CALL cp_cfm_column_scale(hblock, occupation_complex)
1982 DEALLOCATE (occupation_complex)
1983 CALL cp_cfm_transpose(hblock, "C", hblock_h)
1984 CALL cp_cfm_scale_and_add(zhalf, hblock, zhalf, hblock_h)
1985
1986 CALL parallel_gemm("N", "N", nao, nmo, nmo, zone, coeff, hblock, zzero, weighted)
1987 CALL parallel_gemm("N", "C", nao, nao, nmo, zone, weighted, coeff, zzero, wmat)
1988
1989 CALL cp_cfm_release(hblock_h)
1990 CALL cp_cfm_release(hblock)
1991 CALL cp_cfm_release(weighted)
1992
1993 END SUBROUTINE kpoint_ot_energy_weighted_density
1994
1995! **************************************************************************************************
1996!> \brief Common C-to-P algebra for local and distributed k-point groups.
1997!> \param kpoint host environment
1998!> \param ikpgr local k-point index
1999!> \param ispin spin index
2000!> \param wtype construct the energy-weighted density
2001!> \param aux_fit use the auxiliary MO set
2002!> \param nao number of atomic orbitals
2003!> \param nmo number of molecular orbitals
2004!> \param local use local GEMM, without communication
2005!> \param gemm_ctx private local GEMM context, unused for distributed groups
2006!> \param work real column-scaled MO scratch
2007!> \param cwork complex column-scaled MO scratch
2008!> \param cdensity complex density scratch
2009! **************************************************************************************************
2010 SUBROUTINE kpoint_density_matrix_job(kpoint, ikpgr, ispin, wtype, aux_fit, nao, nmo, &
2011 local, gemm_ctx, work, cwork, cdensity)
2012
2013 TYPE(kpoint_type), POINTER :: kpoint
2014 INTEGER, INTENT(IN) :: ikpgr, ispin
2015 LOGICAL, INTENT(IN) :: wtype, aux_fit
2016 INTEGER, INTENT(IN) :: nao, nmo
2017 LOGICAL, INTENT(IN) :: local
2018 TYPE(local_gemm_ctxt_type), INTENT(INOUT) :: gemm_ctx
2019 TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL :: work
2020 TYPE(cp_cfm_type), INTENT(INOUT), OPTIONAL :: cwork, cdensity
2021
2022 COMPLEX(KIND=dp) :: weights(nmo)
2023 REAL(kind=dp), POINTER :: eigenvalues(:), occupation(:)
2024 TYPE(cp_cfm_type), POINTER :: ccoeff
2025 TYPE(cp_fm_type), POINTER :: coeff, cpmat, rpmat
2026 TYPE(kpoint_env_type), POINTER :: kp
2027
2028 IF (aux_fit) THEN
2029 kp => kpoint%kp_aux_env(ikpgr)%kpoint_env
2030 ELSE
2031 kp => kpoint%kp_env(ikpgr)%kpoint_env
2032 END IF
2033 CALL get_mo_set(kp%mos(ispin), occupation_numbers=occupation, eigenvalues=eigenvalues)
2034 IF (wtype) THEN
2035 rpmat => kp%wmat(1, ispin)
2036 ELSE
2037 rpmat => kp%pmat(1, ispin)
2038 END IF
2039 IF (kpoint%use_real_wfn) THEN
2040 cpassert(PRESENT(work))
2041 coeff => kp%mos(ispin)%mo_coeff
2042 CALL cp_fm_to_fm(coeff, work)
2043 CALL cp_fm_column_scale(work, occupation)
2044 IF (wtype) CALL cp_fm_column_scale(work, eigenvalues)
2045 IF (local) THEN
2046 cpassert(product(coeff%matrix_struct%context%num_pe) == 1)
2047 cpassert(product(work%matrix_struct%context%num_pe) == 1)
2048 cpassert(product(rpmat%matrix_struct%context%num_pe) == 1)
2049 CALL gemm_ctx%gemm("N", "T", nao, nao, nmo, 1.0_dp, &
2050 coeff%local_data, SIZE(coeff%local_data, 1), &
2051 work%local_data, SIZE(work%local_data, 1), 0.0_dp, &
2052 rpmat%local_data, SIZE(rpmat%local_data, 1))
2053 ELSE
2054 CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, coeff, work, 0.0_dp, rpmat)
2055 END IF
2056 ELSE
2057 cpassert(PRESENT(cwork) .AND. PRESENT(cdensity))
2058 IF (wtype) THEN
2059 cpmat => kp%wmat(2, ispin)
2060 ELSE
2061 cpmat => kp%pmat(2, ispin)
2062 END IF
2063 ccoeff => kp%mos(ispin)%cmo_coeff
2064 weights = cmplx(occupation(1:nmo), 0.0_dp, kind=dp)
2065 IF (wtype) weights = weights*eigenvalues(1:nmo)
2066 CALL cp_cfm_to_cfm(ccoeff, cwork)
2067 CALL cp_cfm_column_scale(cwork, weights)
2068 IF (local) THEN
2069 cpassert(product(ccoeff%matrix_struct%context%num_pe) == 1)
2070 CALL gemm_ctx%gemm("N", "C", nao, nao, nmo, z_one, &
2071 ccoeff%local_data, SIZE(ccoeff%local_data, 1), &
2072 cwork%local_data, SIZE(cwork%local_data, 1), z_zero, &
2073 cdensity%local_data, SIZE(cdensity%local_data, 1))
2074 ELSE
2075 CALL parallel_gemm("N", "C", nao, nao, nmo, z_one, ccoeff, cwork, z_zero, cdensity)
2076 END IF
2077 CALL cp_cfm_to_fm(cdensity, rpmat, cpmat)
2078 END IF
2079
2080 END SUBROUTINE kpoint_density_matrix_job
2081
2082! **************************************************************************************************
2083!> \brief Calculate Lowdin transformation of density matrix S^1/2 P S^1/2
2084!> Integrate diagonal elements over k-points to get Lowdin charges
2085!> \param kpoint kpoint environment
2086!> \param pmat_diag Sum over kpoints of diagonal elements
2087!> \par History
2088!> 04.2026 created [JGH]
2089! **************************************************************************************************
2090 SUBROUTINE lowdin_kp_trans(kpoint, pmat_diag)
2091
2092 TYPE(kpoint_type), POINTER :: kpoint
2093 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: pmat_diag
2094
2095 CHARACTER(LEN=*), PARAMETER :: routinen = 'lowdin_kp_trans'
2096 COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp, 0.0_dp), &
2097 czero = (0.0_dp, 0.0_dp)
2098
2099 INTEGER :: handle, ikpgr, ispin, kplocal, nao, nspin
2100 INTEGER, DIMENSION(2) :: kp_range
2101 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dele
2102 TYPE(cp_cfm_type) :: cf1work, cf2work
2103 TYPE(cp_cfm_type), POINTER :: cshalf
2104 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
2105 TYPE(cp_fm_type) :: f1work, f2work
2106 TYPE(cp_fm_type), POINTER :: cpmat, pmat, rpmat, shalf
2107 TYPE(kpoint_env_type), POINTER :: kp
2108 TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
2109
2110 CALL timeset(routinen, handle)
2111
2112 nspin = SIZE(pmat_diag, 2)
2113 pmat_diag = 0.0_dp
2114
2115 ! work matrix
2116 CALL cp_fm_get_info(kpoint%kp_env(1)%kpoint_env%pmat(1, 1), &
2117 matrix_struct=matrix_struct, nrow_global=nao)
2118 IF (kpoint%use_real_wfn) THEN
2119 CALL cp_fm_create(f1work, matrix_struct, nrow=nao, ncol=nao)
2120 CALL cp_fm_create(f2work, matrix_struct, nrow=nao, ncol=nao)
2121 ELSE
2122 CALL cp_fm_create(f2work, matrix_struct, nrow=nao, ncol=nao)
2123 CALL cp_cfm_create(cf1work, matrix_struct, nrow=nao, ncol=nao)
2124 CALL cp_cfm_create(cf2work, matrix_struct, nrow=nao, ncol=nao)
2125 END IF
2126 ALLOCATE (dele(nao))
2127
2128 CALL get_kpoint_info(kpoint, kp_range=kp_range)
2129 kplocal = kp_range(2) - kp_range(1) + 1
2130 DO ikpgr = 1, kplocal
2131 kp => kpoint%kp_env(ikpgr)%kpoint_env
2132 DO ispin = 1, nspin
2133 IF (kpoint%use_real_wfn) THEN
2134 pmat => kp%pmat(1, ispin)
2135 shalf => kp%shalf
2136 CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, pmat, shalf, 0.0_dp, f1work)
2137 CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, shalf, f1work, 0.0_dp, f2work)
2138 ELSE
2139 rpmat => kp%pmat(1, ispin)
2140 cpmat => kp%pmat(2, ispin)
2141 cshalf => kp%cshalf
2142 CALL cp_fm_to_cfm(rpmat, cpmat, cf1work)
2143 CALL parallel_gemm("N", "N", nao, nao, nao, cone, cf1work, cshalf, czero, cf2work)
2144 CALL parallel_gemm("N", "N", nao, nao, nao, cone, cshalf, cf2work, czero, cf1work)
2145 CALL cp_cfm_to_fm(cf1work, mtargetr=f2work)
2146 END IF
2147 CALL cp_fm_get_diag(f2work, dele)
2148 pmat_diag(1:nao, ispin) = pmat_diag(1:nao, ispin) + kp%wkp*dele(1:nao)
2149 END DO
2150 END DO
2151
2152 CALL get_kpoint_info(kpoint, para_env_inter_kp=para_env_inter_kp)
2153 CALL para_env_inter_kp%sum(pmat_diag)
2154
2155 IF (kpoint%use_real_wfn) THEN
2156 CALL cp_fm_release(f1work)
2157 CALL cp_fm_release(f2work)
2158 ELSE
2159 CALL cp_fm_release(f2work)
2160 CALL cp_cfm_release(cf1work)
2161 CALL cp_cfm_release(cf2work)
2162 END IF
2163 DEALLOCATE (dele)
2164
2165 CALL timestop(handle)
2166
2167 END SUBROUTINE lowdin_kp_trans
2168
2169! **************************************************************************************************
2170!> \brief Calculate S(k)^1/2 C(k) for real or complex k-point wavefunctions
2171!> \param kp K-point environment for one local k point
2172!> \param ispin Spin index
2173!> \param use_real_wfn Use real k-point wavefunctions
2174!> \param shalfc Output matrix containing S(k)^1/2 C(k) for real wavefunctions
2175!> \param cshalfc Output matrix containing S(k)^1/2 C(k) for complex wavefunctions
2176! **************************************************************************************************
2177 SUBROUTINE lowdin_kp_mo_coeff(kp, ispin, use_real_wfn, shalfc, cshalfc)
2178
2179 TYPE(kpoint_env_type), POINTER :: kp
2180 INTEGER, INTENT(IN) :: ispin
2181 LOGICAL, INTENT(IN) :: use_real_wfn
2182 TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL :: shalfc
2183 TYPE(cp_cfm_type), INTENT(INOUT), OPTIONAL :: cshalfc
2184
2185 INTEGER :: nao, nmo
2186 TYPE(mo_set_type), POINTER :: mo_set
2187
2188 IF (use_real_wfn) THEN
2189 cpassert(PRESENT(shalfc))
2190 mo_set => kp%mos(ispin)
2191 CALL get_mo_set(mo_set, nao=nao, nmo=nmo)
2192
2193 CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, kp%shalf, &
2194 mo_set%mo_coeff, 0.0_dp, shalfc)
2195 ELSE
2196 cpassert(PRESENT(cshalfc))
2197 mo_set => kp%mos(ispin)
2198 CALL get_mo_set(mo_set, nao=nao, nmo=nmo)
2199 CALL parallel_gemm("N", "N", nao, nmo, nao, z_one, kp%cshalf, &
2200 mo_set%cmo_coeff, z_zero, cshalfc)
2201 END IF
2202
2203 END SUBROUTINE lowdin_kp_mo_coeff
2204
2205! **************************************************************************************************
2206!> \brief generate real space density matrices in DBCSR format
2207!> \param kpoint Kpoint environment
2208!> \param denmat Real space (DBCSR) density matrices
2209!> \param wtype True = energy weighted density matrix
2210!> False = normal density matrix
2211!> \param tempmat DBCSR matrix to be used as template
2212!> \param sab_nl ...
2213!> \param fmwork FM work matrices (kpoint group)
2214!> \param for_aux_fit ...
2215!> \param pmat_ext ...
2216!> \param overlap_rs ...
2217! **************************************************************************************************
2218 SUBROUTINE kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwork, for_aux_fit, &
2219 pmat_ext, overlap_rs)
2220
2221 TYPE(kpoint_type), POINTER :: kpoint
2222 TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
2223 LOGICAL, INTENT(IN) :: wtype
2224 TYPE(dbcsr_type), POINTER :: tempmat
2225 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2226 POINTER :: sab_nl
2227 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN), TARGET :: fmwork
2228 LOGICAL, OPTIONAL :: for_aux_fit
2229 TYPE(cp_fm_type), DIMENSION(:, :, :), INTENT(IN), &
2230 OPTIONAL, TARGET :: pmat_ext
2231 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
2232 POINTER :: overlap_rs
2233
2234 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_density_transform'
2235
2236 INTEGER :: handle, ic, ik, ikp, ispin, kplocal, nc, &
2237 nimg, nkp, nspin
2238 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2239 LOGICAL :: aux_fit, do_symmetric, local, &
2240 regular_grid_done, &
2241 regular_grid_eligible
2242 TYPE(cp_fm_p_type), ALLOCATABLE, &
2243 DIMENSION(:, :, :) :: source
2244 TYPE(dbcsr_p_type), DIMENSION(2) :: components
2245 TYPE(dbcsr_type), POINTER :: cpmat, rpmat
2246 TYPE(kp_transform_plan_type) :: transform_plan
2247 TYPE(kpoint_env_type), POINTER :: kp
2248 TYPE(mp_para_env_type), POINTER :: para_env
2249
2250 CALL timeset(routinen, handle)
2251
2252 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
2253
2254 IF (PRESENT(for_aux_fit)) THEN
2255 aux_fit = for_aux_fit
2256 ELSE
2257 aux_fit = .false.
2258 END IF
2259
2260 IF (aux_fit) THEN
2261 cpassert(ASSOCIATED(kpoint%kp_aux_env))
2262 END IF
2263
2264 ! work storage
2265 ALLOCATE (rpmat)
2266 CALL dbcsr_create(rpmat, template=tempmat, &
2267 matrix_type=merge(dbcsr_type_symmetric, dbcsr_type_no_symmetry, do_symmetric))
2268 CALL cp_dbcsr_alloc_block_from_nbl(rpmat, sab_nl)
2269 CALL dbcsr_set(rpmat, 0.0_dp)
2270 NULLIFY (cpmat)
2271 components(1)%matrix => rpmat
2272
2273 CALL get_kpoint_info(kpoint, nkp=nkp, cell_to_index=cell_to_index)
2274 IF (PRESENT(overlap_rs)) THEN
2275 CALL calibrate_symmetry_phases(kpoint, overlap_rs, tempmat, sab_nl, cell_to_index)
2276 END IF
2277 ! initialize real space density matrices
2278 IF (aux_fit) THEN
2279 kp => kpoint%kp_aux_env(1)%kpoint_env
2280 ELSE
2281 kp => kpoint%kp_env(1)%kpoint_env
2282 END IF
2283 nspin = SIZE(kp%mos)
2284 nc = SIZE(kp%pmat, 1)
2285 nimg = SIZE(denmat, 2)
2286 CALL kp_transform_plan_create(transform_plan, sab_nl, cell_to_index, nimg, &
2287 group_entries=.true.)
2288
2289 ! Borrow one common view for ordinary P, W, auxiliary and external matrices.
2290 ! Pointers remain within this call, including the TARGET lifetime of pmat_ext.
2291 kplocal = kpoint%kp_range(2) - kpoint%kp_range(1) + 1
2292 ALLOCATE (source(kplocal, nc, nspin))
2293 DO ispin = 1, nspin
2294 DO ikp = 1, kplocal
2295 IF (aux_fit) THEN
2296 kp => kpoint%kp_aux_env(ikp)%kpoint_env
2297 ELSE
2298 kp => kpoint%kp_env(ikp)%kpoint_env
2299 END IF
2300 DO ic = 1, nc
2301 IF (PRESENT(pmat_ext)) THEN
2302 source(ikp, ic, ispin)%matrix => pmat_ext(ikp, ic, ispin)
2303 ELSE IF (wtype) THEN
2304 source(ikp, ic, ispin)%matrix => kp%wmat(ic, ispin)
2305 ELSE
2306 source(ikp, ic, ispin)%matrix => kp%pmat(ic, ispin)
2307 END IF
2308 END DO
2309 END DO
2310 END DO
2311
2312 para_env => kpoint%blacs_env_all%para_env
2313 local = para_env%num_pe == 1
2314 IF (local) THEN
2315 cpassert(kpoint%kp_range(1) == 1 .AND. kplocal == nkp)
2316 END IF
2317 ! The regular-grid route is an optional specialization without spatial rotations.
2318 ! The batched block transform handles all remaining cases.
2319 regular_grid_eligible = nimg < nkp
2320 IF (ASSOCIATED(kpoint%index_to_cell)) THEN
2321 regular_grid_eligible = regular_grid_eligible .AND. nimg <= SIZE(kpoint%index_to_cell, 2)
2322 ELSE
2323 regular_grid_eligible = .false.
2324 END IF
2325 IF (regular_grid_eligible) THEN
2326 DO ik = 1, nkp
2327 IF (kpoint%kp_sym(ik)%kpoint_sym%apply_symmetry) THEN
2328 regular_grid_eligible = .false.
2329 EXIT
2330 END IF
2331 END DO
2332 END IF
2333
2334 regular_grid_done = .false.
2335 IF (regular_grid_eligible) THEN
2336 ALLOCATE (cpmat)
2337 CALL dbcsr_create(cpmat, template=rpmat, &
2338 matrix_type=merge(dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, do_symmetric))
2339 CALL cp_dbcsr_alloc_block_from_nbl(cpmat, sab_nl)
2340 CALL dbcsr_set(cpmat, 0.0_dp)
2341 components(2)%matrix => cpmat
2342 CALL kpoint_density_transform_regular_grid(kpoint, source, denmat, components, fmwork, &
2343 transform_plan, regular_grid_done)
2344 END IF
2345 IF (.NOT. regular_grid_done) THEN
2346 CALL kpoint_density_transform_batched(kpoint, source, denmat, rpmat, transform_plan)
2347 END IF
2348
2349 CALL dbcsr_deallocate_matrix(rpmat)
2350 IF (ASSOCIATED(cpmat)) CALL dbcsr_deallocate_matrix(cpmat)
2351
2352 CALL timestop(handle)
2353
2354 END SUBROUTINE kpoint_density_transform
2355
2356! **************************************************************************************************
2357!> \brief Contract group-local densities into bounded batches of real-space AO blocks.
2358!> Small tiles of images sharing an AO block are independent OpenMP work items. Rotations follow the
2359!> k-point contraction; only the resulting blocks are reduced to their DBCSR
2360!> owner. This also works when a source AO block spans several MPI ranks.
2361!> \param kpoint ...
2362!> \param source borrowed P, W or external density matrices
2363!> \param denmat real-space output matrices
2364!> \param template union of the source AO sparsities
2365!> \param plan neighbor-list traversal and storage conventions
2366! **************************************************************************************************
2367 SUBROUTINE kpoint_density_transform_batched(kpoint, source, denmat, template, plan)
2368 TYPE(kpoint_type), POINTER :: kpoint
2369 TYPE(cp_fm_p_type), DIMENSION(:, :, :), INTENT(IN) :: source
2370 TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
2371 TYPE(dbcsr_type), POINTER :: template
2372 TYPE(kp_transform_plan_type), INTENT(IN) :: plan
2373
2374 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_density_transform_batched'
2375 INTEGER(KIND=int_8), PARAMETER :: max_bytes = 64_int_8*1024*1024
2376 INTEGER, PARAMETER :: images_per_tile = 8, max_batch = 256
2377
2378 INTEGER :: a, b, batch, first, handle, i, ic, igroup, ikp, image, ispin, j, last, nblock, &
2379 ncol, njob, nlocal, nrow, nthreads, ntile, nvalue, orientation, owner, slot, thread, tile
2380 INTEGER(KIND=int_8) :: block_size
2381 INTEGER, ALLOCATABLE :: available(:, :, :), counts(:), &
2382 jobs(:, :), local_jobs(:, :), &
2383 offsets(:), tile_start(:), &
2384 tiles(:, :, :, :, :)
2385 INTEGER, DIMENSION(:), POINTER :: col_offset, col_size, row_offset, &
2386 row_size
2387 LOGICAL :: found, trans
2388 REAL(kind=dp), ALLOCATABLE :: packed(:, :), phases(:, :), &
2389 rotation_work(:, :)
2390 REAL(kind=dp), ALLOCATABLE, TARGET :: sums(:, :), values(:)
2391 REAL(kind=dp), POINTER :: block(:, :)
2392 TYPE(kp_density_group_type), ALLOCATABLE :: groups(:)
2393 TYPE(local_gemm_ctxt_type), ALLOCATABLE :: gemm_ctx(:)
2394 TYPE(mp_para_env_type), POINTER :: world
2395
2396 CALL timeset(routinen, handle)
2397 world => kpoint%blacs_env_all%para_env
2398 CALL dbcsr_get_info(template, row_blk_size=row_size, col_blk_size=col_size, &
2399 row_blk_offset=row_offset, col_blk_offset=col_offset)
2400 CALL kp_density_groups_create(kpoint, SIZE(row_size), groups)
2401 cpassert(SIZE(row_size) == SIZE(col_size))
2402 ALLOCATE (tiles(4, SIZE(row_size), SIZE(source, 1), SIZE(source, 2), SIZE(source, 3)))
2403 DO ispin = 1, SIZE(source, 3)
2404 DO ic = 1, SIZE(source, 2)
2405 DO ikp = 1, SIZE(source, 1)
2406 CALL kp_density_layout_create(source(ikp, ic, ispin)%matrix, row_size, col_size, &
2407 row_offset, col_offset, tiles(:, :, ikp, ic, ispin))
2408 END DO
2409 END DO
2410 DO image = 1, SIZE(denmat, 2)
2411 CALL dbcsr_set(denmat(ispin, image)%matrix, 0.0_dp)
2412 END DO
2413 END DO
2414
2415 ! Only owners request blocks. A group contains one (image,row,col), hence
2416 ! all its entries have the same cell; retain both multiplicity and sign.
2417 ALLOCATE (local_jobs(8, plan%ngroup), counts(0:world%num_pe - 1))
2418 nlocal = 0
2419 DO igroup = 1, plan%ngroup
2420 first = plan%group_start(igroup)
2421 last = plan%group_start(igroup + 1) - 1
2422 CALL dbcsr_get_readonly_block_p(denmat(1, plan%image(first))%matrix, &
2423 plan%row(first), plan%col(first), block, found=found)
2424 IF (.NOT. found) cycle
2425 nlocal = nlocal + 1
2426 local_jobs(:, nlocal) = [plan%row(first), plan%col(first), plan%image(first), &
2427 plan%cell(:, first), last - first + 1, &
2428 nint(sum(plan%symmetry_sign(first:last)))]
2429 END DO
2430 CALL world%allgather(nlocal, counts)
2431 IF (maxval(counts) == 0) THEN
2432 CALL timestop(handle)
2433 RETURN
2434 END IF
2435 block_size = int(maxval(row_size), int_8)*maxval(col_size)
2436 batch = max(1, min(max_batch, maxval(counts), int(max_bytes/(8_int_8*block_size))))
2437 ALLOCATE (jobs(8, batch), offsets(batch + 1), tile_start(batch + 1), values(int(block_size)*batch), &
2438 available(2, 0:ubound(groups, 1), batch))
2439 nthreads = 1
2440!$ nthreads = MIN(omp_get_max_threads(), batch)
2441 ALLOCATE (gemm_ctx(nthreads))
2442 DO thread = 1, nthreads
2443 CALL gemm_ctx(thread)%create(local_gemm_pu_host, timing=.false.)
2444 END DO
2445
2446 DO owner = 0, world%num_pe - 1
2447 DO first = 1, counts(owner), batch
2448 njob = min(batch, counts(owner) - first + 1)
2449 IF (world%mepos == owner) jobs(:, 1:njob) = local_jobs(:, first:first + njob - 1)
2450 IF (world%num_pe > 1) CALL world%bcast(jobs(:, 1:njob), owner)
2451 offsets(1) = 1
2452 DO j = 1, njob
2453 offsets(j + 1) = offsets(j) + row_size(jobs(1, j))*col_size(jobs(2, j))
2454 END DO
2455 ! Neighbor-list groups are ordered by AO pair and image. Limit image
2456 ! tiles to retain enough independent work items for small unit cells.
2457 ntile = 0
2458 i = 1
2459 DO WHILE (i <= njob)
2460 ntile = ntile + 1
2461 tile_start(ntile) = i
2462 j = i + 1
2463 DO WHILE (j <= min(njob, i + images_per_tile - 1))
2464 IF (any(jobs(1:2, j) /= jobs(1:2, i))) EXIT
2465 j = j + 1
2466 END DO
2467 i = j
2468 END DO
2469 tile_start(ntile + 1) = njob + 1
2470 available(:, :, :) = 0
2471 DO tile = 1, ntile
2472 j = tile_start(tile)
2473 DO slot = 0, ubound(groups, 1)
2474 IF (.NOT. ALLOCATED(groups(slot)%inverse)) cycle
2475 DO orientation = 1, 2
2476 CALL kp_density_source_block(groups(slot), jobs(1:2, j), plan%symmetric, &
2477 orientation, a, b, trans)
2478 IF (a == 0) cycle
2479 CALL dbcsr_get_readonly_block_p(template, a, b, block, found=found)
2480 available(orientation, slot, tile) = merge(1, 0, found)
2481 END DO
2482 END DO
2483 END DO
2484 ! Source sparsity is distributed independently of the FM layout.
2485 ! Reduce only presence flags, never the k-dependent source matrices.
2486 IF (world%num_pe > 1) CALL world%sum(available(:, :, 1:ntile))
2487 nvalue = offsets(njob + 1) - 1
2488 nthreads = min(SIZE(gemm_ctx), ntile)
2489 DO ispin = 1, SIZE(source, 3)
2490 values(1:nvalue) = 0.0_dp
2491!$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads) &
2492!$OMP SHARED(source,tiles,groups,jobs,available,plan,row_size,col_size,row_offset,col_offset, &
2493!$OMP values,offsets,tile_start,ntile,ispin,gemm_ctx) &
2494!$OMP PRIVATE(tile,j,last,nrow,ncol,block,packed,phases,sums,rotation_work,thread)
2495 thread = 1
2496!$ thread = omp_get_thread_num() + 1
2497!$OMP DO SCHEDULE(DYNAMIC, 1)
2498 DO tile = 1, ntile
2499 j = tile_start(tile)
2500 last = tile_start(tile + 1) - 1
2501 nrow = row_size(jobs(1, j))
2502 ncol = col_size(jobs(2, j))
2503 block(1:nrow*ncol, 1:last - j + 1) => values(offsets(j):offsets(last + 1) - 1)
2504 CALL kp_density_contract_tile(source(:, :, ispin), tiles(:, :, :, :, ispin), groups, &
2505 jobs(:, j:last), available(:, :, tile), plan%symmetric, &
2506 row_size, col_size, row_offset, col_offset, block, &
2507 packed, phases, sums, rotation_work, gemm_ctx(thread))
2508 END DO
2509!$OMP END DO
2510!$OMP END PARALLEL
2511 IF (world%num_pe > 1) CALL world%sum(values(1:nvalue), root=owner)
2512 IF (world%mepos /= owner) cycle
2513 DO j = 1, njob
2514 CALL dbcsr_get_block_p(denmat(ispin, jobs(3, j))%matrix, jobs(1, j), jobs(2, j), block, found=found)
2515 cpassert(found)
2516 nrow = SIZE(block, 1)
2517 DO i = 1, SIZE(block, 2)
2518 nblock = offsets(j) + (i - 1)*nrow
2519 block(:, i) = values(nblock:nblock + nrow - 1)
2520 END DO
2521 END DO
2522 END DO
2523 END DO
2524 END DO
2525 DO thread = 1, SIZE(gemm_ctx)
2526 CALL gemm_ctx(thread)%destroy()
2527 END DO
2528 CALL timestop(handle)
2529 END SUBROUTINE kpoint_density_transform_batched
2530
2531! **************************************************************************************************
2532!> \brief Group the existing density contributions by operation and phase direction.
2533!> Only nwght entries contribute; additional little-group mappings are not
2534!> extra density terms. Metadata is borrowed, local k-point lists are owned.
2535!> \param kpoint ...
2536!> \param natom number of AO blocks
2537!> \param groups operation views and local contribution lists; zero is the identity
2538! **************************************************************************************************
2539 SUBROUTINE kp_density_groups_create(kpoint, natom, groups)
2540 TYPE(kpoint_type), POINTER :: kpoint
2541 INTEGER, INTENT(IN) :: natom
2542 TYPE(kp_density_group_type), ALLOCATABLE, &
2543 INTENT(OUT) :: groups(:)
2544
2545 INTEGER :: i, ik, ikp, is, ncontrib, nrot, pass, &
2546 slot
2547 INTEGER, ALLOCATABLE :: counts(:)
2548 TYPE(kpoint_sym_type), POINTER :: kpsym
2549
2550 nrot = 0
2551 IF (ASSOCIATED(kpoint%ibrot)) nrot = SIZE(kpoint%ibrot)
2552 ALLOCATE (groups(0:4*nrot), counts(0:4*nrot))
2553 DO pass = 1, 2
2554 counts(:) = 0
2555 DO ik = 1, kpoint%nkp
2556 ikp = ik - kpoint%kp_range(1) + 1
2557 kpsym => kpoint%kp_sym(ik)%kpoint_sym
2558 ncontrib = 1
2559 IF (kpsym%apply_symmetry) ncontrib = kpsym%nwght
2560 DO is = 1, ncontrib
2561 slot = 0
2562 IF (kpsym%apply_symmetry) THEN
2563 cpassert(kpsym%phase_mode(is) >= 1 .AND. kpsym%phase_mode(is) <= 2)
2564 slot = find_kpoint_rotation_slot(kpoint, kpsym%rotp(is))
2565 cpassert(slot > 0)
2566 slot = slot + nrot*merge(1, 0, kpsym%rotp(is) < 0) + 2*nrot*(kpsym%phase_mode(is) - 1)
2567 END IF
2568 IF (pass == 1) THEN
2569 IF (.NOT. ALLOCATED(groups(slot)%inverse)) THEN
2570 ALLOCATE (groups(slot)%inverse(natom))
2571 groups(slot)%op%identity = .true.
2572 DO i = 1, natom
2573 groups(slot)%inverse(i) = i
2574 END DO
2575 IF (slot > 0) THEN
2576 CALL kp_symmetry_op_init(groups(slot)%op, kpoint, kpsym, is)
2577 groups(slot)%reverse_phase = kpsym%phase_mode(is) == 2
2578 DO i = 1, natom
2579 groups(slot)%inverse(groups(slot)%op%atom_map(i)) = i
2580 END DO
2581 END IF
2582 ELSE IF (slot > 0) THEN
2583 ! One crystallographic operation has one atom/gauge mapping.
2584 cpassert(all(groups(slot)%op%atom_map == kpsym%f0(:, is)))
2585 cpassert(all(groups(slot)%op%cell_shift == kpsym%fcell_gauge(:, :, is)))
2586 END IF
2587 END IF
2588 IF (ik < kpoint%kp_range(1) .OR. ik > kpoint%kp_range(2)) cycle
2589 counts(slot) = counts(slot) + 1
2590 IF (pass == 2) THEN
2591 i = counts(slot)
2592 groups(slot)%ikp(i) = ikp
2593 groups(slot)%weight(i) = kpoint%wkp(ik)/real(ncontrib, dp)
2594 IF (slot == 0) THEN
2595 groups(slot)%xkp(:, i) = kpoint%xkp(:, ik)
2596 ELSE
2597 groups(slot)%xkp(:, i) = kpsym%xkp(:, is)
2598 END IF
2599 END IF
2600 END DO
2601 END DO
2602 IF (pass == 1) THEN
2603 DO slot = 0, 4*nrot
2604 ALLOCATE (groups(slot)%ikp(counts(slot)), groups(slot)%weight(counts(slot)), &
2605 groups(slot)%xkp(3, counts(slot)))
2606 END DO
2607 END IF
2608 END DO
2609 END SUBROUTINE kp_density_groups_create
2610
2611! **************************************************************************************************
2612!> \brief Map each AO block to its local FM row and column intervals, excluding padding.
2613!> \param matrix ...
2614!> \param row_size ...
2615!> \param col_size ...
2616!> \param row_offset ...
2617!> \param col_offset ...
2618!> \param tile local first/last row and column indices for each atom
2619! **************************************************************************************************
2620 SUBROUTINE kp_density_layout_create(matrix, row_size, col_size, row_offset, col_offset, tile)
2621 TYPE(cp_fm_type), INTENT(IN) :: matrix
2622 INTEGER, INTENT(IN) :: row_size(:), col_size(:), row_offset(:), &
2623 col_offset(:)
2624 INTEGER, INTENT(OUT) :: tile(:, :)
2625
2626 INTEGER :: a, i, ncol, nrow
2627 INTEGER, ALLOCATABLE :: cols(:), rows(:)
2628 INTEGER, POINTER :: col_indices(:), row_indices(:)
2629
2630 cpassert(matrix%matrix_struct%nrow_global == sum(row_size))
2631 cpassert(matrix%matrix_struct%ncol_global == sum(col_size))
2632 CALL cp_fm_get_info(matrix, nrow_local=nrow, ncol_local=ncol, &
2633 row_indices=row_indices, col_indices=col_indices)
2634 ALLOCATE (rows(0:sum(row_size)), cols(0:sum(col_size)))
2635 rows(:) = 0
2636 cols(:) = 0
2637 DO i = 1, nrow
2638 rows(row_indices(i)) = 1
2639 END DO
2640 DO i = 1, ncol
2641 cols(col_indices(i)) = 1
2642 END DO
2643 DO i = 1, ubound(rows, 1)
2644 rows(i) = rows(i) + rows(i - 1)
2645 END DO
2646 DO i = 1, ubound(cols, 1)
2647 cols(i) = cols(i) + cols(i - 1)
2648 END DO
2649 DO a = 1, SIZE(row_size)
2650 tile(1:2, a) = [rows(row_offset(a) - 1) + 1, rows(row_offset(a) + row_size(a) - 1)]
2651 END DO
2652 DO a = 1, SIZE(col_size)
2653 tile(3:4, a) = [cols(col_offset(a) - 1) + 1, cols(col_offset(a) + col_size(a) - 1)]
2654 END DO
2655 END SUBROUTINE kp_density_layout_create
2656
2657! **************************************************************************************************
2658!> \brief Find the source orientation that maps to a stored output block.
2659!> A nonsymmetric source can contribute both orientations to an upper block,
2660!> as in the DBCSR symmetry adapter. An absent orientation returns a=0.
2661!> \param group ...
2662!> \param TARGET output atom pair
2663!> \param symmetric source storage convention
2664!> \param orientation ...
2665!> \param a ...
2666!> \param b ...
2667!> \param trans whether rotation transposes the source block
2668! **************************************************************************************************
2669 SUBROUTINE kp_density_source_block(group, TARGET, symmetric, orientation, a, b, trans)
2670 TYPE(kp_density_group_type), INTENT(IN) :: group
2671 INTEGER, INTENT(IN) :: target(2)
2672 LOGICAL, INTENT(IN) :: symmetric
2673 INTEGER, INTENT(IN) :: orientation
2674 INTEGER, INTENT(OUT) :: a, b
2675 LOGICAL, INTENT(OUT) :: trans
2676
2677 INTEGER :: i, j
2678
2679 a = 0
2680 b = 0
2681 trans = .false.
2682 IF (group%op%identity) THEN
2683 IF (orientation == 1) THEN
2684 a = target(1)
2685 b = target(2)
2686 END IF
2687 RETURN
2688 END IF
2689 IF (target(1) > target(2)) RETURN
2690 i = group%inverse(target(1))
2691 j = group%inverse(target(2))
2692 IF (orientation == 2 .AND. (symmetric .OR. i == j)) RETURN
2693 a = min(i, j)
2694 b = max(i, j)
2695 IF (orientation == 2) THEN
2696 a = max(i, j)
2697 b = min(i, j)
2698 END IF
2699 trans = group%op%atom_map(a) > group%op%atom_map(b)
2700 END SUBROUTINE kp_density_source_block
2701
2702! **************************************************************************************************
2703!> \brief Contract a small tile of images sharing an output AO block.
2704!> Pack bounded k-point batches once and reuse them across images in GEMM.
2705!> Gauge/Fourier phases remain inside the contraction; AO rotations follow
2706!> the completed local sum. MPI can then reduce the rotated partial blocks.
2707!> \param source ...
2708!> \param tiles local first/last row and column indices by atom, k-point and component
2709!> \param groups ...
2710!> \param jobs row, col, image, cell, real multiplicity and imaginary sign sum
2711!> \param available global source-block presence by orientation and operation
2712!> \param symmetric ...
2713!> \param row_size ...
2714!> \param col_size ...
2715!> \param row_offset ...
2716!> \param col_offset ...
2717!> \param TARGET flattened output blocks, one column per image
2718!> \param packed reusable thread-local source batch
2719!> \param phases reusable thread-local contraction coefficients
2720!> \param sums reusable thread-local contracted blocks
2721!> \param rotation_work reusable thread-local GEMM buffer
2722!> \param gemm_ctx thread-owned local GEMM context
2723! **************************************************************************************************
2724 SUBROUTINE kp_density_contract_tile(source, tiles, groups, jobs, available, symmetric, &
2725 row_size, col_size, row_offset, col_offset, TARGET, &
2726 packed, phases, sums, rotation_work, gemm_ctx)
2727 TYPE(cp_fm_p_type), INTENT(IN) :: source(:, :)
2728 INTEGER, INTENT(IN) :: tiles(:, :, :, :)
2729 TYPE(kp_density_group_type), INTENT(IN) :: groups(0:)
2730 INTEGER, INTENT(IN) :: jobs(:, :), available(:, 0:)
2731 LOGICAL, INTENT(IN) :: symmetric
2732 INTEGER, INTENT(IN) :: row_size(:), col_size(:), row_offset(:), &
2733 col_offset(:)
2734 REAL(kind=dp), CONTIGUOUS, INTENT(INOUT), TARGET :: target(:, :)
2735 REAL(kind=dp), ALLOCATABLE, INTENT(INOUT) :: packed(:, :), phases(:, :)
2736 REAL(kind=dp), ALLOCATABLE, INTENT(INOUT), TARGET :: sums(:, :)
2737 REAL(kind=dp), ALLOCATABLE, INTENT(INOUT) :: rotation_work(:, :)
2738 TYPE(local_gemm_ctxt_type), INTENT(INOUT) :: gemm_ctx
2739
2740 INTEGER(KIND=int_8), PARAMETER :: max_elements = 65536_int_8
2741
2742 INTEGER :: a, b, batch, column, first, i, ic, &
2743 ifirst, ii, ikp, ilast, image, j, jj, &
2744 k, last, nc, nelement, nimage, nrow, &
2745 orientation, shift(3), slot
2746 LOGICAL :: contributes, real_only, trans
2747 REAL(kind=dp) :: arg, cosine, factors(2, 2), sine
2748 REAL(kind=dp), POINTER :: left_rot(:, :), right_rot(:, :), &
2749 source_block(:, :), target_block(:, :)
2750 TYPE(cp_fm_type), POINTER :: matrix
2751
2752! Up to 512 KiB of packed input per worker, except for a single larger AO block.
2753
2754 nc = SIZE(source, 2)
2755 real_only = nc == 1
2756 nimage = SIZE(jobs, 2)
2757 DO slot = 0, ubound(groups, 1)
2758 IF (SIZE(groups(slot)%ikp) == 0) cycle
2759 DO orientation = 1, 2
2760 IF (available(orientation, slot) == 0) cycle
2761 CALL kp_density_source_block(groups(slot), jobs(1:2, 1), symmetric, orientation, a, b, trans)
2762 cpassert(a > 0)
2763 nelement = row_size(a)*col_size(b)
2764 batch = max(1, min(128, SIZE(groups(slot)%ikp), int(max_elements/(int(nelement, int_8)*nc))))
2765 CALL ensure_work_matrix(packed, nelement, nc*batch)
2766 CALL ensure_work_matrix(phases, nc*batch, nimage)
2767 CALL ensure_work_matrix(sums, nelement, nimage)
2768 sums(:, :) = 0.0_dp
2769 shift(:) = 0
2770 IF (.NOT. groups(slot)%op%identity) THEN
2771 shift(:) = groups(slot)%op%cell_shift(:, b) - groups(slot)%op%cell_shift(:, a)
2772 IF (groups(slot)%reverse_phase) shift(:) = -shift
2773 END IF
2774 DO first = 1, SIZE(groups(slot)%ikp), batch
2775 last = min(first + batch - 1, SIZE(groups(slot)%ikp))
2776 packed(:, :) = 0.0_dp
2777 contributes = .false.
2778 DO k = first, last
2779 ikp = groups(slot)%ikp(k)
2780 CALL kp_symmetry_phase_factors(shift, groups(slot)%xkp(:, k), &
2781 groups(slot)%op%time_reversal, real_only, factors)
2782 DO image = 1, nimage
2783 arg = twopi*dot_product(real(jobs(4:6, image), dp), groups(slot)%xkp(:, k))
2784 cosine = real(jobs(7, image), dp)*cos(arg)
2785 sine = real(jobs(8, image), dp)*sin(arg)*merge(-1.0_dp, 1.0_dp, trans)
2786 IF (real_only) sine = 0.0_dp
2787 DO ic = 1, nc
2788 column = (k - first)*nc + ic
2789 phases(column, image) = groups(slot)%weight(k)*(cosine*factors(1, ic) + sine*factors(2, ic))
2790 END DO
2791 END DO
2792 DO ic = 1, nc
2793 ifirst = tiles(1, a, ikp, ic)
2794 ilast = tiles(2, a, ikp, ic)
2795 IF (ifirst > ilast .OR. tiles(3, b, ikp, ic) > tiles(4, b, ikp, ic)) cycle
2796 contributes = .true.
2797 column = (k - first)*nc + ic
2798 matrix => source(ikp, ic)%matrix
2799 nrow = ilast - ifirst + 1
2800 DO j = tiles(3, b, ikp, ic), tiles(4, b, ikp, ic)
2801 jj = matrix%matrix_struct%col_indices(j) - col_offset(b)
2802 ii = matrix%matrix_struct%row_indices(ifirst) - row_offset(a) + 1 + jj*row_size(a)
2803 IF (matrix%matrix_struct%row_indices(ilast) - &
2804 matrix%matrix_struct%row_indices(ifirst) == nrow - 1) THEN
2805 packed(ii:ii + nrow - 1, column) = matrix%local_data(ifirst:ilast, j)
2806 ELSE
2807 DO i = ifirst, ilast
2808 ii = matrix%matrix_struct%row_indices(i) - row_offset(a) + 1 + jj*row_size(a)
2809 packed(ii, column) = matrix%local_data(i, j)
2810 END DO
2811 END IF
2812 END DO
2813 END DO
2814 END DO
2815 IF (.NOT. contributes) cycle
2816 CALL gemm_ctx%gemm('N', 'N', nelement, nimage, nc*(last - first + 1), 1.0_dp, &
2817 packed, SIZE(packed, 1), phases, SIZE(phases, 1), &
2818 1.0_dp, sums, SIZE(sums, 1))
2819 END DO
2820 IF (groups(slot)%op%identity) THEN
2821 target(:, :) = TARGET + sums
2822 ELSE
2823 left_rot => groups(slot)%op%kind_rot(groups(slot)%op%atom_kind(groups(slot)%inverse(jobs(1, 1))))%rmat
2824 right_rot => groups(slot)%op%kind_rot(groups(slot)%op%atom_kind(groups(slot)%inverse(jobs(2, 1))))%rmat
2825 DO image = 1, nimage
2826 source_block(1:row_size(a), 1:col_size(b)) => sums(:, image)
2827 target_block(1:row_size(jobs(1, image)), 1:col_size(jobs(2, image))) => target(:, image)
2828 CALL kp_symmetry_rotate_block(source_block, left_rot, right_rot, trans, &
2829 1.0_dp, target_block, rotation_work)
2830 END DO
2831 END IF
2832 END DO
2833 END DO
2834 END SUBROUTINE kp_density_contract_tile
2835
2836! **************************************************************************************************
2837!> \brief Sum compatible group-local FM tiles before returning each image to DBCSR.
2838!> The single-rank case stores directly. Incompatible layouts and cheaper
2839!> local sparse sums return to the common batched transform in the caller.
2840!> \param kpoint ...
2841!> \param source borrowed density matrices for all local k-points, components and spins
2842!> \param denmat ...
2843!> \param components real/imaginary DBCSR work matrices
2844!> \param fmwork global FM work matrices, unused on one rank
2845!> \param plan neighbor-list traversal and storage conventions
2846!> \param completed whether all spin/image results were produced
2847! **************************************************************************************************
2848 SUBROUTINE kpoint_density_transform_regular_grid(kpoint, source, denmat, components, fmwork, plan, completed)
2849
2850 TYPE(kpoint_type), POINTER :: kpoint
2851 TYPE(cp_fm_p_type), DIMENSION(:, :, :), INTENT(IN) :: source
2852 TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
2853 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: components
2854 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN), TARGET :: fmwork
2855 TYPE(kp_transform_plan_type), INTENT(IN) :: plan
2856 LOGICAL, INTENT(OUT) :: completed
2857
2858 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_density_transform_regular_grid'
2859
2860 INTEGER :: compatible, element, first, handle, i, &
2861 ic, igroup, ik, ikp, image, ispin, j, &
2862 kplocal, nc, ncol, nrow
2863 INTEGER, DIMENSION(14) :: layout, layout_max, layout_min
2864 INTEGER, DIMENSION(2) :: fft_range
2865 LOGICAL :: found, local, same_layout
2866 REAL(kind=dp) :: arg, dense_work, sparse_work
2867 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: phase
2868 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: fft_density
2869 REAL(kind=dp), DIMENSION(:, :), POINTER :: block
2870 TYPE(copy_info_type), DIMENSION(2) :: info
2871 TYPE(cp_fm_struct_type), POINTER :: fms
2872 TYPE(cp_fm_type) :: dummy
2873 TYPE(cp_fm_type), DIMENSION(2), TARGET :: partial
2874 TYPE(cp_fm_type), DIMENSION(:), POINTER :: image_source
2875 TYPE(mp_para_env_type), POINTER :: inter, para_env
2876
2877 completed = .false.
2878 para_env => kpoint%blacs_env_all%para_env
2879 inter => kpoint%para_env_inter_kp
2880 local = para_env%num_pe == 1
2881 kplocal = SIZE(source, 1)
2882 nc = SIZE(source, 2)
2883 cpassert(nc >= 1 .AND. nc <= 2)
2884 fms => source(1, 1, 1)%matrix%matrix_struct
2885 same_layout = .true.
2886 DO ispin = 1, SIZE(source, 3)
2887 DO ic = 1, nc
2888 DO ikp = 1, kplocal
2889 same_layout = same_layout .AND. cp_fm_struct_equivalent(fms, source(ikp, ic, ispin)%matrix%matrix_struct)
2890 same_layout = same_layout .AND. &
2891 all(fms%first_p_pos == source(ikp, ic, ispin)%matrix%matrix_struct%first_p_pos)
2892 END DO
2893 END DO
2894 END DO
2895 CALL cp_fm_get_info(source(1, 1, 1)%matrix, nrow_local=nrow, ncol_local=ncol)
2896 layout = [fms%nrow_global, fms%ncol_global, fms%nrow_block, fms%ncol_block, &
2897 fms%first_p_pos, fms%context%num_pe, fms%context%mepos, &
2898 shape(source(1, 1, 1)%matrix%local_data), nrow, ncol]
2899 layout_min = layout
2900 layout_max = layout
2901 IF (inter%num_pe > 1) THEN
2902 CALL inter%min(layout_min)
2903 CALL inter%max(layout_max)
2904 END IF
2905 compatible = merge(1, 0, same_layout .AND. all(layout_min == layout_max))
2906 IF (.NOT. local) CALL para_env%min(compatible)
2907 IF (compatible == 0) RETURN
2908
2909 ! Count the actual stored block work, so a sparse local transform is not
2910 ! replaced by a dense sum merely because the MPI image route is eligible.
2911 sparse_work = -1.0_dp
2912 dense_work = real(nrow, dp)*real(ncol, dp)*SIZE(denmat, 2)
2913 IF (local) THEN
2914 sparse_work = 0.0_dp
2915 DO igroup = 1, plan%ngroup
2916 first = plan%group_start(igroup)
2917 CALL dbcsr_get_readonly_block_p(components(1)%matrix, plan%row(first), plan%col(first), block, found=found)
2918 IF (found) sparse_work = sparse_work + real(SIZE(block), dp)
2919 END DO
2920 END IF
2921
2922 CALL timeset(routinen, handle)
2923 ALLOCATE (phase(kplocal, nc))
2924 DO ispin = 1, SIZE(source, 3)
2925 CALL kp_density_fft(kpoint, source(:, :, ispin), nrow, ncol, SIZE(denmat, 2), &
2926 fft_range, fft_density, local)
2927 IF (ispin == 1) THEN
2928 IF (local .AND. .NOT. ALLOCATED(fft_density) .AND. dense_work > sparse_work) EXIT
2929 DO ic = 1, nc
2930 CALL cp_fm_create(partial(ic), fms)
2931 END DO
2932 image_source => partial
2933 IF (.NOT. local) THEN
2934 cpassert(SIZE(fmwork) >= nc)
2935 image_source => fmwork
2936 END IF
2937 END IF
2938 DO image = 1, SIZE(denmat, 2)
2939 DO ic = 1, nc
2940 partial(ic)%local_data(:, :) = 0.0_dp
2941 END DO
2942 IF (ALLOCATED(fft_density)) THEN
2943!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) &
2944!$OMP SHARED(partial, fft_density, fft_range, image, nc, nrow) PRIVATE(element, i, j, ic)
2945 DO element = fft_range(1), fft_range(2)
2946 i = modulo(element - 1, nrow) + 1
2947 j = (element - 1)/nrow + 1
2948 DO ic = 1, nc
2949 partial(ic)%local_data(i, j) = fft_density(element - fft_range(1) + 1, ic, image)
2950 END DO
2951 END DO
2952!$OMP END PARALLEL DO
2953 ELSE
2954 DO ikp = 1, kplocal
2955 ik = kpoint%kp_range(1) + ikp - 1
2956 arg = twopi*dot_product(real(kpoint%index_to_cell(:, image), dp), kpoint%xkp(:, ik))
2957 phase(ikp, 1) = kpoint%wkp(ik)*cos(arg)
2958 IF (nc == 2) phase(ikp, 2) = kpoint%wkp(ik)*sin(arg)
2959 END DO
2960 ! Keep real/imaginary sums separate for the neighbor-list signs.
2961!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) &
2962!$OMP SHARED(partial, source, phase, nc, nrow, ncol, kplocal, ispin) PRIVATE(j, ikp, ic)
2963 DO j = 1, ncol
2964 DO ikp = 1, kplocal
2965 DO ic = 1, nc
2966 partial(ic)%local_data(1:nrow, j) = partial(ic)%local_data(1:nrow, j) + &
2967 phase(ikp, ic)*source(ikp, ic, ispin)%matrix%local_data(1:nrow, j)
2968 END DO
2969 END DO
2970 END DO
2971!$OMP END PARALLEL DO
2972 END IF
2973
2974 IF (.NOT. local) THEN
2975 DO ic = 1, nc
2976 IF (inter%num_pe > 1) CALL inter%sum(partial(ic)%local_data)
2977 IF (kpoint%iogrp) THEN
2978 CALL cp_fm_start_copy_general(partial(ic), fmwork(ic), para_env, info(ic))
2979 ELSE
2980 CALL cp_fm_start_copy_general(dummy, fmwork(ic), para_env, info(ic))
2981 END IF
2982 END DO
2983 DO ic = 1, nc
2984 CALL cp_fm_finish_copy_general(fmwork(ic), info(ic))
2985 IF (kpoint%iogrp) CALL cp_fm_cleanup_copy_general(info(ic))
2986 END DO
2987 END IF
2988 DO ic = 1, nc
2989 CALL kp_copy_fm_to_dbcsr(image_source(ic), components(ic)%matrix, local)
2990 END DO
2991 CALL kp_accumulate_density_image(denmat(ispin, image)%matrix, components(1)%matrix, &
2992 components(2)%matrix, image, nc == 1, plan)
2993 END DO
2994 IF (ALLOCATED(fft_density)) DEALLOCATE (fft_density)
2995 completed = .true.
2996 END DO
2997 IF (completed) THEN
2998 DO ic = 1, nc
2999 CALL cp_fm_release(partial(ic))
3000 END DO
3001 END IF
3002 CALL timestop(handle)
3003
3004 END SUBROUTINE kpoint_density_transform_regular_grid
3005
3006! **************************************************************************************************
3007!> \brief Exchange k-distributed FM entries into batched FFT pencils.
3008!> Each inter-group rank transforms disjoint FM entries, not a full replicated
3009!> set of matrices. Real F(R) = A(R)+B(R), real F(-R) = A(R)-B(R), where A and B
3010!> are the weighted cosine/real and sine/imaginary sums needed by the DBCSR
3011!> storage convention. Original reduced-grid weights are used without doubling.
3012!> \param kpoint ...
3013!> \param source current group-local real/imaginary density matrices
3014!> \param nrow number of valid local FM rows (excluding padding)
3015!> \param ncol number of valid local FM columns
3016!> \param nimg number of requested real-space images
3017!> \param owned local FM element range assigned to this inter-group rank
3018!> \param density allocated on success; unallocated selects the regular-grid direct sum
3019!> \param local all k-points and matrix entries are on this rank
3020! **************************************************************************************************
3021 SUBROUTINE kp_density_fft(kpoint, source, nrow, ncol, nimg, owned, density, local)
3022
3023 TYPE(kpoint_type), POINTER :: kpoint
3024 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(IN) :: source
3025 INTEGER, INTENT(IN) :: nrow, ncol, nimg
3026 INTEGER, DIMENSION(2), INTENT(OUT) :: owned
3027 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
3028 INTENT(OUT) :: density
3029 LOGICAL, INTENT(IN) :: local
3030
3031 INTEGER(KIND=int_8), PARAMETER :: max_bytes = 64_int_8*1024*1024, &
3032 min_fft_peer_bytes = 32_int_8*1024
3033 INTEGER, PARAMETER :: max_batch = 32
3034
3035 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:), &
3036 TARGET :: recvbuf
3037 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: values_rs
3038 COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:), &
3039 POINTER :: sendbuf
3040 COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
3041 POINTER :: values_k
3042 INTEGER :: batch, count, element, first, group, handle, i, ikp, j, last, local_k, &
3043 min_local_k, mine, nc, nel, ngroup, nkp, nowned, offset, pos, width
3044 INTEGER(KIND=int_8) :: fixed_bytes, peer_bytes, per_entry, &
3045 tile_size
3046 INTEGER, ALLOCATABLE, DIMENSION(:) :: rcount, rdispl, scount, sdispl
3047 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: cells, ranges
3048 INTEGER, DIMENSION(2) :: choice
3049 INTEGER, DIMENSION(3) :: nfft
3050 LOGICAL :: compatible
3051 TYPE(k_grid_to_cell_work_type) :: fft_work_data
3052 TYPE(mp_para_env_type), POINTER :: inter
3053
3054 owned = [1, 0]
3055 IF (kpoint%lattice_fft == lattice_fft_off) RETURN
3056 nkp = SIZE(kpoint%xkp, 2)
3057 IF (kpoint%lattice_fft == lattice_fft_auto .AND. nkp < 27) RETURN
3058 ! Preserve AUTO's local direct sum until a batched FFT executor can amortize
3059 ! per-transform threading. ON uses FFT on any eligible rank/thread layout.
3060 IF (kpoint%lattice_fft == lattice_fft_auto .AND. local) RETURN
3061 IF (any(kpoint%nkp_grid <= 0)) RETURN
3062 ! Bound map/scratch allocations before attempting a candidate-grid mapping.
3063 IF (32_int_8*product(int(kpoint%nkp_grid, int_8)) > max_bytes) RETURN
3064 compatible = lattice_fft_shape(kpoint%nkp_grid, nfft, allow_arbitrary=.true.)
3065 IF (.NOT. compatible) RETURN
3066 IF (32_int_8*product(int(nfft, int_8)) > max_bytes) RETURN
3067 compatible = regular_kpoint_grid(kpoint%xkp, kpoint%nkp_grid, allow_incomplete=.true.)
3068 IF (.NOT. compatible) RETURN
3069
3070 inter => kpoint%para_env_inter_kp
3071 ngroup = inter%num_pe
3072 mine = inter%mepos + 1
3073 local_k = SIZE(source, 1)
3074 nc = SIZE(source, 2)
3075 tile_size = int(nrow, int_8)*ncol
3076 ! Avoid default-integer overflow, including float-rounded partition limits.
3077 nel = 0
3078 IF (tile_size <= huge(nel)) nel = int(tile_size)
3079 ALLOCATE (ranges(2, ngroup))
3080 DO group = 1, ngroup
3081 ranges(1, group) = int(int(nel, int_8)*(group - 1)/ngroup) + 1
3082 ranges(2, group) = int(int(nel, int_8)*group/ngroup)
3083 END DO
3084 owned = ranges(:, mine)
3085 nowned = max(0, owned(2) - owned(1) + 1)
3086 fixed_bytes = 8_int_8*nowned*nc*nimg + 32_int_8*product(int(nfft, int_8)) + &
3087 4_int_8*product(int(kpoint%nkp_grid, int_8)) + 44_int_8*nkp + 52_int_8*nc*nimg + &
3088 24_int_8*ngroup + 64_int_8
3089 per_entry = 16_int_8*(int(ngroup, int_8)*local_k + nkp + int(nc, int_8)*nimg)
3090 width = int(max(0_int_8, min(int(max_batch, int_8), (max_bytes - fixed_bytes)/per_entry)))
3091 choice = [merge(1, 0, tile_size <= huge(nel)), width]
3092 ! All global ranks must either use FFT pencils or retain the direct path.
3093 IF (kpoint%blacs_env_all%para_env%num_pe > 1) CALL kpoint%blacs_env_all%para_env%min(choice)
3094 IF (choice(1) == 0 .OR. choice(2) == 0) RETURN
3095 width = choice(2)
3096 IF (kpoint%lattice_fft == lattice_fft_auto .AND. ngroup > 1) THEN
3097 ! With several k-point groups, small per-peer payloads do not amortize
3098 ! the all-to-all latency and the regular-grid direct sum is faster.
3099 min_local_k = minval(kpoint%kp_dist(2, :) - kpoint%kp_dist(1, :) + 1)
3100 peer_bytes = 16_int_8*int(width, int_8)*int(min_local_k, int_8)
3101 IF (peer_bytes < min_fft_peer_bytes) RETURN
3102 END IF
3103
3104 CALL timeset("kp_density_fft", handle)
3105 ALLOCATE (density(nowned, nc, nimg))
3106 ALLOCATE (cells(3, nc*nimg))
3107 cells(:, 1:nimg) = kpoint%index_to_cell(:, 1:nimg)
3108 IF (nc == 2) cells(:, nimg + 1:) = -cells(:, 1:nimg)
3109 ALLOCATE (scount(ngroup), sdispl(ngroup), rcount(ngroup), rdispl(ngroup))
3110 ALLOCATE (recvbuf(max(1, width*nkp)))
3111 IF (ngroup == 1) THEN
3112 sendbuf => recvbuf
3113 ELSE
3114 ALLOCATE (sendbuf(max(1, width*ngroup*local_k)))
3115 END IF
3116 ALLOCATE (values_rs(width, 1, nc*nimg))
3117 IF (nowned > 0) THEN
3118 CALL k_grid_to_cell_prepare(fft_work_data, kpoint%xkp, kpoint%nkp_grid, cells, &
3119 weights=kpoint%wkp, allow_incomplete=.true.)
3120 END IF
3121 DO batch = 1, maxval(ranges(2, :) - ranges(1, :) + 1), width
3122 count = max(0, min(width, nowned - batch + 1))
3123 offset = 0
3124 DO group = 1, ngroup
3125 first = ranges(1, group) + batch - 1
3126 last = min(first + width - 1, ranges(2, group))
3127 scount(group) = max(0, last - first + 1)*local_k
3128 sdispl(group) = offset
3129 offset = offset + scount(group)
3130 rcount(group) = count*(kpoint%kp_dist(2, group) - kpoint%kp_dist(1, group) + 1)
3131 rdispl(group) = count*(kpoint%kp_dist(1, group) - 1)
3132 ! Pack elements in k-major order; displacements already encode the
3133 ! unequal k ranges, so neither counts nor indices need an exchange.
3134!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) &
3135!$OMP SHARED(source, sendbuf, sdispl, group, first, last, nc, nrow, local_k) &
3136!$OMP PRIVATE(ikp, element, i, j, pos)
3137 DO ikp = 1, local_k
3138 DO element = first, last
3139 i = modulo(element - 1, nrow) + 1
3140 j = (element - 1)/nrow + 1
3141 pos = sdispl(group) + (ikp - 1)*(last - first + 1) + element - first + 1
3142 sendbuf(pos) = cmplx(source(ikp, 1)%matrix%local_data(i, j), 0.0_dp, dp)
3143 IF (nc == 2) sendbuf(pos) = cmplx(real(sendbuf(pos), dp), &
3144 source(ikp, 2)%matrix%local_data(i, j), dp)
3145 END DO
3146 END DO
3147!$OMP END PARALLEL DO
3148 END DO
3149 IF (ngroup > 1) CALL inter%alltoall(sendbuf, scount, sdispl, recvbuf, rcount, rdispl)
3150 IF (count == 0) cycle
3151 values_k(1:count, 1:1, 1:nkp) => recvbuf(1:count*nkp)
3152 ! Keep fft_tools pool/planner access outside OpenMP worker regions.
3153 CALL k_grid_to_cell_execute(fft_work_data, values_k, values_rs(1:count, :, :))
3154 DO i = 1, nimg
3155 IF (nc == 1) THEN
3156 density(batch:batch + count - 1, 1, i) = real(values_rs(1:count, 1, i), dp)
3157 ELSE
3158 density(batch:batch + count - 1, 1, i) = 0.5_dp* &
3159 REAL(values_rs(1:count, 1, i) + values_rs(1:count, 1, nimg + i), dp)
3160 density(batch:batch + count - 1, 2, i) = 0.5_dp* &
3161 REAL(values_rs(1:count, 1, i) - values_rs(1:count, 1, nimg + i), dp)
3162 END IF
3163 END DO
3164 END DO
3165 CALL k_grid_to_cell_release(fft_work_data)
3166 IF (ngroup > 1) DEALLOCATE (sendbuf)
3167 NULLIFY (sendbuf)
3168 CALL timestop(handle)
3169
3170 END SUBROUTINE kp_density_fft
3171
3172! **************************************************************************************************
3173!> \brief Apply the neighbor-list storage convention to an already k-summed image.
3174!> Cell phases and k weights have been applied in the FM distribution;
3175!> retain the multiplicity and imaginary sign of each neighbor-list entry.
3176!> \param denmat ...
3177!> \param rpmat ...
3178!> \param cpmat ...
3179!> \param image ...
3180!> \param real_only ...
3181!> \param plan ...
3182! **************************************************************************************************
3183 SUBROUTINE kp_accumulate_density_image(denmat, rpmat, cpmat, image, real_only, plan)
3184
3185 TYPE(dbcsr_type), POINTER :: denmat, rpmat, cpmat
3186 INTEGER, INTENT(IN) :: image
3187 LOGICAL, INTENT(IN) :: real_only
3188 TYPE(kp_transform_plan_type), INTENT(IN) :: plan
3189
3190 INTEGER :: first, igroup, last
3191 LOGICAL :: found
3192 REAL(kind=dp), DIMENSION(:, :), POINTER :: cblock, dblock, rblock
3193
3194 CALL dbcsr_set(denmat, 0.0_dp)
3195!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) &
3196!$OMP SHARED(denmat, rpmat, cpmat, image, real_only, plan) &
3197!$OMP PRIVATE(igroup, first, last, found, dblock, rblock, cblock)
3198 DO igroup = 1, plan%ngroup
3199 first = plan%group_start(igroup)
3200 IF (plan%image(first) /= image) cycle
3201 last = plan%group_start(igroup + 1) - 1
3202 CALL dbcsr_get_block_p(denmat, plan%row(first), plan%col(first), dblock, found=found)
3203 IF (.NOT. found) cycle
3204 CALL dbcsr_get_readonly_block_p(rpmat, plan%row(first), plan%col(first), rblock, found=found)
3205 IF (.NOT. found) cycle
3206 IF (.NOT. real_only) THEN
3207 CALL dbcsr_get_readonly_block_p(cpmat, plan%row(first), plan%col(first), cblock, found=found)
3208 IF (.NOT. found) cycle
3209 END IF
3210 dblock = real(last - first + 1, dp)*rblock
3211 IF (.NOT. real_only) dblock = dblock + sum(plan%symmetry_sign(first:last))*cblock
3212 END DO
3213!$OMP END PARALLEL DO
3214
3215 END SUBROUTINE kp_accumulate_density_image
3216
3217! **************************************************************************************************
3218!> \brief Copy density FM data to existing DBCSR blocks, with no redistribution on one rank.
3219!> \param fm ...
3220!> \param matrix ...
3221!> \param local all matrix data and DBCSR blocks belong to this rank
3222! **************************************************************************************************
3223 SUBROUTINE kp_copy_fm_to_dbcsr(fm, matrix, local)
3224
3225 TYPE(cp_fm_type), INTENT(IN) :: fm
3226 TYPE(dbcsr_type), POINTER :: matrix
3227 LOGICAL, INTENT(IN) :: local
3228
3229 INTEGER :: col_offset, row_offset
3230 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
3231 POINTER :: full
3232 REAL(kind=dp), DIMENSION(:, :), POINTER :: block
3233 TYPE(dbcsr_iterator_type) :: iterator
3234
3235 IF (.NOT. local) THEN
3236 CALL copy_fm_to_dbcsr(fm, matrix, keep_sparsity=.true.)
3237 RETURN
3238 END IF
3239 CALL cp_fm_get_info(fm, local_data=full)
3240 CALL dbcsr_iterator_start(iterator, matrix, shared=.false.)
3241 DO WHILE (dbcsr_iterator_blocks_left(iterator))
3242 CALL dbcsr_iterator_next_block(iterator, block=block, row_offset=row_offset, &
3243 col_offset=col_offset)
3244 block(:, :) = full(row_offset:row_offset + SIZE(block, 1) - 1, &
3245 col_offset:col_offset + SIZE(block, 2) - 1)
3246 END DO
3247 CALL dbcsr_iterator_stop(iterator)
3248
3249 END SUBROUTINE kp_copy_fm_to_dbcsr
3250
3251! **************************************************************************************************
3252!> \brief Build the immutable neighbor-list traversal shared by R-to-K and K-to-R transforms.
3253!> \param plan ...
3254!> \param sab_nl ...
3255!> \param cell_to_index ...
3256!> \param nimg ...
3257!> \param block_template ...
3258!> \param group_entries ...
3259! **************************************************************************************************
3260 SUBROUTINE kp_transform_plan_create(plan, sab_nl, cell_to_index, nimg, block_template, group_entries)
3261
3262 TYPE(kp_transform_plan_type), INTENT(OUT) :: plan
3263 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3264 POINTER :: sab_nl
3265 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3266 INTEGER, INTENT(IN) :: nimg
3267 TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: block_template
3268 LOGICAL, INTENT(IN), OPTIONAL :: group_entries
3269
3270 INTEGER :: i, iatom, icell, icol, igroup, irow, &
3271 jatom, nblock
3272 INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: key
3273 INTEGER, ALLOCATABLE, DIMENSION(:) :: order
3274 INTEGER, DIMENSION(3) :: cell
3275 INTEGER, DIMENSION(:), POINTER :: col_offsets, row_offsets
3276 LOGICAL :: do_grouping, do_symmetric, store_offsets
3277 TYPE(neighbor_list_iterator_p_type), &
3278 DIMENSION(:), POINTER :: iterator
3279
3280 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
3281 plan%symmetric = do_symmetric
3282 do_grouping = .false.
3283 IF (PRESENT(group_entries)) do_grouping = group_entries
3284 store_offsets = PRESENT(block_template)
3285 IF (store_offsets) THEN
3286 CALL dbcsr_get_info(block_template, row_blk_offset=row_offsets, col_blk_offset=col_offsets)
3287 END IF
3288 CALL neighbor_list_iterator_create(iterator, sab_nl)
3289 DO WHILE (neighbor_list_iterate(iterator) == 0)
3290 CALL get_iterator_info(iterator, cell=cell)
3291 icell = cell_to_index(cell(1), cell(2), cell(3))
3292 IF (icell >= 1 .AND. icell <= nimg) plan%nentry = plan%nentry + 1
3293 END DO
3294 CALL neighbor_list_iterator_release(iterator)
3295
3296 ALLOCATE (plan%row(plan%nentry), plan%col(plan%nentry), plan%image(plan%nentry), &
3297 plan%cell(3, plan%nentry), plan%symmetry_sign(plan%nentry))
3298 IF (store_offsets) THEN
3299 ALLOCATE (plan%row_offset(plan%nentry), plan%col_offset(plan%nentry))
3300 END IF
3301
3302 i = 0
3303 CALL neighbor_list_iterator_create(iterator, sab_nl)
3304 DO WHILE (neighbor_list_iterate(iterator) == 0)
3305 CALL get_iterator_info(iterator, iatom=iatom, jatom=jatom, cell=cell)
3306 icell = cell_to_index(cell(1), cell(2), cell(3))
3307 IF (icell < 1 .OR. icell > nimg) cycle
3308 i = i + 1
3309 irow = iatom
3310 icol = jatom
3311 plan%symmetry_sign(i) = 1.0_dp
3312 IF (do_symmetric .AND. iatom > jatom) THEN
3313 irow = jatom
3314 icol = iatom
3315 plan%symmetry_sign(i) = -1.0_dp
3316 END IF
3317 plan%row(i) = irow
3318 plan%col(i) = icol
3319 IF (store_offsets) THEN
3320 plan%row_offset(i) = row_offsets(irow)
3321 plan%col_offset(i) = col_offsets(icol)
3322 END IF
3323 plan%image(i) = icell
3324 plan%cell(:, i) = cell
3325 END DO
3326 CALL neighbor_list_iterator_release(iterator)
3327 cpassert(i == plan%nentry)
3328
3329 IF (do_grouping .AND. plan%nentry > 0) THEN
3330 nblock = max(maxval(plan%row), maxval(plan%col))
3331 ALLOCATE (key(plan%nentry), order(plan%nentry))
3332 DO i = 1, plan%nentry
3333 key(i) = int(plan%image(i), kind=int_8) + int(nimg, kind=int_8)* &
3334 (int(plan%row(i) - 1, kind=int_8) + int(nblock, kind=int_8)* &
3335 int(plan%col(i) - 1, kind=int_8))
3336 END DO
3337 CALL sort(key, plan%nentry, order)
3338 plan%row(:) = plan%row(order)
3339 plan%col(:) = plan%col(order)
3340 IF (store_offsets) THEN
3341 plan%row_offset(:) = plan%row_offset(order)
3342 plan%col_offset(:) = plan%col_offset(order)
3343 END IF
3344 plan%image(:) = plan%image(order)
3345 plan%cell(:, :) = plan%cell(:, order)
3346 plan%symmetry_sign(:) = plan%symmetry_sign(order)
3347
3348 plan%ngroup = 1
3349 DO i = 2, plan%nentry
3350 IF (key(i) /= key(i - 1)) plan%ngroup = plan%ngroup + 1
3351 END DO
3352 ALLOCATE (plan%group_start(plan%ngroup + 1))
3353 igroup = 1
3354 plan%group_start(igroup) = 1
3355 DO i = 2, plan%nentry
3356 IF (key(i) /= key(i - 1)) THEN
3357 igroup = igroup + 1
3358 plan%group_start(igroup) = i
3359 END IF
3360 END DO
3361 plan%group_start(plan%ngroup + 1) = plan%nentry + 1
3362 DEALLOCATE (key, order)
3363 END IF
3364
3365 END SUBROUTINE kp_transform_plan_create
3366
3367! **************************************************************************************************
3368!> \brief Allocate a dense work matrix with the requested shape
3369!> \param work dense work matrix
3370!> \param nrow number of rows
3371!> \param ncol number of columns
3372! **************************************************************************************************
3373 SUBROUTINE ensure_work_matrix(work, nrow, ncol)
3374
3375 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
3376 INTENT(INOUT) :: work
3377 INTEGER, INTENT(IN) :: nrow, ncol
3378
3379 IF (ALLOCATED(work)) THEN
3380 IF (SIZE(work, 1) == nrow .AND. SIZE(work, 2) == ncol) RETURN
3381 DEALLOCATE (work)
3382 END IF
3383 ALLOCATE (work(nrow, ncol))
3384
3385 END SUBROUTINE ensure_work_matrix
3386
3387! **************************************************************************************************
3388!> \brief Select the Bloch-phase convention that preserves overlap covariance.
3389!> \param kpoint ...
3390!> \param overlap_rs ...
3391!> \param tempmat ...
3392!> \param sab_nl ...
3393!> \param cell_to_index ...
3394! **************************************************************************************************
3395 SUBROUTINE calibrate_symmetry_phases(kpoint, overlap_rs, tempmat, sab_nl, cell_to_index)
3396
3397 TYPE(kpoint_type), POINTER :: kpoint
3398 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: overlap_rs
3399 TYPE(dbcsr_type), POINTER :: tempmat
3400 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3401 POINTER :: sab_nl
3402 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3403
3404 CHARACTER(LEN=*), PARAMETER :: routinen = 'calibrate_symmetry_phases'
3405
3406 CHARACTER(LEN=256) :: phase_error
3407 INTEGER :: best_mode, handle, iatom, ik, is, key, &
3408 mode, nrot, rot_slot
3409 INTEGER, ALLOCATABLE :: best_ik(:), best_is(:), mode_by_op(:)
3410 LOGICAL :: needs_calibration, reverse
3411 REAL(kind=dp) :: arg, best_residual, candidate_norm, &
3412 direct_norm, overlap_dot, phase_score, &
3413 phase_tolerance, relative_residual
3414 REAL(kind=dp), ALLOCATABLE :: best_score(:)
3415 TYPE(dbcsr_type), POINTER :: direct_c, direct_r, source_c, source_r, &
3416 sym_c, sym_r
3417 TYPE(kp_symmetry_op_type) :: op
3418 TYPE(kpoint_sym_type), POINTER :: kpsym
3419
3420 needs_calibration = .false.
3421 DO ik = 1, kpoint%nkp
3422 kpsym => kpoint%kp_sym(ik)%kpoint_sym
3423 IF (kpsym%apply_symmetry) THEN
3424 ! Only the first nwght operations enter the density transform.
3425 IF (any(kpsym%phase_mode(1:kpsym%nwght) == 0)) THEN
3426 needs_calibration = .true.
3427 EXIT
3428 END IF
3429 END IF
3430 END DO
3431 IF (.NOT. needs_calibration) RETURN
3432
3433 CALL timeset(routinen, handle)
3434 cpassert(ASSOCIATED(kpoint%ibrot))
3435 nrot = SIZE(kpoint%ibrot)
3436 ALLOCATE (best_ik(2*nrot), best_is(2*nrot), mode_by_op(2*nrot), best_score(2*nrot))
3437 best_ik(:) = 0
3438 best_is(:) = 0
3439 mode_by_op(:) = 0
3440 best_score(:) = -1.0_dp
3441
3442 ! The Bloch-phase direction is a convention of the signed crystallographic
3443 ! operation, not of an individual member of its k-point star. Select one
3444 ! numerically discriminating representative per signed operation and reuse
3445 ! the result for all equivalent contributions.
3446 DO ik = 1, kpoint%nkp
3447 kpsym => kpoint%kp_sym(ik)%kpoint_sym
3448 IF (.NOT. kpsym%apply_symmetry) cycle
3449 DO is = 1, kpsym%nwght
3450 rot_slot = find_kpoint_rotation_slot(kpoint, kpsym%rotp(is))
3451 cpassert(rot_slot > 0)
3452 key = rot_slot + nrot*merge(1, 0, kpsym%rotp(is) < 0)
3453 IF (kpsym%phase_mode(is) > 0) THEN
3454 IF (mode_by_op(key) == 0) THEN
3455 mode_by_op(key) = kpsym%phase_mode(is)
3456 ELSE
3457 cpassert(mode_by_op(key) == kpsym%phase_mode(is))
3458 END IF
3459 cycle
3460 END IF
3461
3462 phase_score = 0.0_dp
3463 DO iatom = 2, SIZE(kpsym%fcell_gauge, 2)
3464 arg = dot_product(real(kpsym%fcell_gauge(:, iatom, is) - &
3465 kpsym%fcell_gauge(:, 1, is), kind=dp), &
3466 kpsym%xkp(:, is))
3467 phase_score = max(phase_score, abs(sin(twopi*arg)))
3468 END DO
3469 IF (phase_score > best_score(key)) THEN
3470 best_score(key) = phase_score
3471 best_ik(key) = ik
3472 best_is(key) = is
3473 END IF
3474 END DO
3475 END DO
3476
3477 ALLOCATE (source_r, source_c, direct_r, direct_c, sym_r, sym_c)
3478 CALL dbcsr_create(source_r, template=tempmat, matrix_type=dbcsr_type_symmetric)
3479 CALL dbcsr_create(source_c, template=tempmat, matrix_type=dbcsr_type_antisymmetric)
3480 CALL dbcsr_create(direct_r, template=tempmat, matrix_type=dbcsr_type_symmetric)
3481 CALL dbcsr_create(direct_c, template=tempmat, matrix_type=dbcsr_type_antisymmetric)
3482 CALL dbcsr_create(sym_r, template=tempmat, matrix_type=dbcsr_type_symmetric)
3483 CALL dbcsr_create(sym_c, template=tempmat, matrix_type=dbcsr_type_antisymmetric)
3484 CALL cp_dbcsr_alloc_block_from_nbl(source_r, sab_nl)
3485 CALL cp_dbcsr_alloc_block_from_nbl(source_c, sab_nl)
3486 CALL cp_dbcsr_alloc_block_from_nbl(direct_r, sab_nl)
3487 CALL cp_dbcsr_alloc_block_from_nbl(direct_c, sab_nl)
3488 CALL cp_dbcsr_alloc_block_from_nbl(sym_r, sab_nl)
3489 CALL cp_dbcsr_alloc_block_from_nbl(sym_c, sab_nl)
3490
3491 phase_tolerance = max(1.0e-6_dp, 100.0_dp*kpoint%eps_geo)
3492 DO key = 1, 2*nrot
3493 IF (mode_by_op(key) > 0 .OR. best_ik(key) == 0) cycle
3494 ! If both phase directions are algebraically indistinguishable for every
3495 ! occurrence of this operation, avoid an unnecessary overlap transform.
3496 IF (best_score(key) <= 1000.0_dp*epsilon(1.0_dp)) THEN
3497 mode_by_op(key) = 1
3498 cycle
3499 END IF
3500
3501 ik = best_ik(key)
3502 is = best_is(key)
3503 kpsym => kpoint%kp_sym(ik)%kpoint_sym
3504 CALL dbcsr_set(source_r, 0.0_dp)
3505 CALL dbcsr_set(source_c, 0.0_dp)
3506 CALL rskp_transform(source_r, source_c, overlap_rs, 1, kpoint%xkp(1:3, ik), &
3507 cell_to_index, sab_nl)
3508 CALL dbcsr_set(direct_r, 0.0_dp)
3509 CALL dbcsr_set(direct_c, 0.0_dp)
3510 CALL rskp_transform(direct_r, direct_c, overlap_rs, 1, kpsym%xkp(1:3, is), &
3511 cell_to_index, sab_nl)
3512 CALL dbcsr_dot(direct_r, direct_r, direct_norm)
3513 CALL dbcsr_dot(direct_c, direct_c, candidate_norm)
3514 direct_norm = direct_norm + candidate_norm
3515
3516 CALL kp_symmetry_op_init(op, kpoint, kpsym, is)
3517 best_mode = 0
3518 best_residual = huge(1.0_dp)
3519 DO mode = 1, 2
3520 reverse = mode == 2
3521 CALL symtrans_phase(sym_r, sym_c, source_r, source_c, .false., op, reverse)
3522 CALL dbcsr_dot(sym_r, sym_r, candidate_norm)
3523 CALL dbcsr_dot(sym_c, sym_c, relative_residual)
3524 candidate_norm = candidate_norm + relative_residual
3525 CALL dbcsr_dot(sym_r, direct_r, overlap_dot)
3526 CALL dbcsr_dot(sym_c, direct_c, relative_residual)
3527 overlap_dot = overlap_dot + relative_residual
3528 relative_residual = sqrt(max(0.0_dp, candidate_norm + direct_norm - &
3529 2.0_dp*overlap_dot)/max(direct_norm, tiny(1.0_dp)))
3530 IF (relative_residual < best_residual) THEN
3531 best_residual = relative_residual
3532 best_mode = mode
3533 END IF
3534 END DO
3535 IF (best_residual > phase_tolerance) THEN
3536 WRITE (phase_error, '(A,ES12.4,A,I0,A,I0)') &
3537 "No Bloch-phase direction preserves overlap covariance; residual=", &
3538 best_residual, ", irreducible k-point=", ik, ", operation=", is
3539 CALL cp_abort(__location__, trim(phase_error))
3540 END IF
3541 mode_by_op(key) = best_mode
3542 END DO
3543
3544 DO ik = 1, kpoint%nkp
3545 kpsym => kpoint%kp_sym(ik)%kpoint_sym
3546 IF (.NOT. kpsym%apply_symmetry) cycle
3547 DO is = 1, kpsym%nwght
3548 IF (kpsym%phase_mode(is) > 0) cycle
3549 rot_slot = find_kpoint_rotation_slot(kpoint, kpsym%rotp(is))
3550 cpassert(rot_slot > 0)
3551 key = rot_slot + nrot*merge(1, 0, kpsym%rotp(is) < 0)
3552 cpassert(mode_by_op(key) > 0)
3553 kpsym%phase_mode(is) = mode_by_op(key)
3554 END DO
3555 END DO
3556
3557 CALL dbcsr_deallocate_matrix(source_r)
3558 CALL dbcsr_deallocate_matrix(source_c)
3559 CALL dbcsr_deallocate_matrix(direct_r)
3560 CALL dbcsr_deallocate_matrix(direct_c)
3561 CALL dbcsr_deallocate_matrix(sym_r)
3562 CALL dbcsr_deallocate_matrix(sym_c)
3563 CALL timestop(handle)
3564
3565 END SUBROUTINE calibrate_symmetry_phases
3566
3567! **************************************************************************************************
3568!> \brief Borrow the atom mapping and basis rotations for one symmetry contribution.
3569!> \param op non-owning operation view
3570!> \param kpoint k-point environment
3571!> \param kpsym symmetry data for the representative k-point
3572!> \param is index in the existing symmetry contribution list
3573! **************************************************************************************************
3574 SUBROUTINE kp_symmetry_op_init(op, kpoint, kpsym, is)
3575 TYPE(kp_symmetry_op_type), INTENT(OUT) :: op
3576 TYPE(kpoint_type), POINTER :: kpoint
3577 TYPE(kpoint_sym_type), POINTER :: kpsym
3578 INTEGER, INTENT(IN) :: is
3579
3580 INTEGER :: iatom, rot_slot
3581 LOGICAL :: has_phase, perm, rotates
3582 REAL(kind=dp) :: dr
3583 REAL(kind=dp), DIMENSION(:, :), POINTER :: rot
3584
3585 rot_slot = find_kpoint_rotation_slot(kpoint, kpsym%rotp(is))
3586 cpassert(rot_slot > 0)
3587 op%kind_rot => kpoint%kind_rotmat(rot_slot, :)
3588 op%atom_map => kpsym%f0(:, is)
3589 op%atom_kind => kpoint%atype
3590 op%cell_shift => kpsym%fcell_gauge(:, :, is)
3591 op%xkp(:) = kpsym%xkp(:, is)
3592 op%time_reversal = kpsym%rotp(is) < 0
3593
3594 perm = .false.
3595 DO iatom = 1, SIZE(op%atom_map)
3596 IF (op%atom_map(iatom) == iatom) cycle
3597 perm = .true.
3598 EXIT
3599 END DO
3600 rot => kpsym%rot(:, :, is)
3601 dr = abs(rot(1, 1) - 1.0_dp) + abs(rot(2, 2) - 1.0_dp) + abs(rot(3, 3) - 1.0_dp)
3602 rotates = abs(sum(abs(rot)) - 3.0_dp) > 1.e-12_dp .OR. abs(dr) > 1.e-12_dp
3603 has_phase = any(op%cell_shift /= 0) .OR. op%time_reversal
3604 op%identity = .NOT. (rotates .OR. perm .OR. has_phase)
3605 op%redistribute = perm .OR. has_phase
3606
3607 END SUBROUTINE kp_symmetry_op_init
3608
3609! **************************************************************************************************
3610!> \brief DBCSR adapter for a complex k-point symmetry operation.
3611!> \param srpmat real part of transformed matrix
3612!> \param scpmat imaginary part of transformed matrix
3613!> \param rpmat real part of reference matrix
3614!> \param cpmat imaginary part of reference matrix
3615!> \param real_only whether only the real part is present
3616!> \param op symmetry operation and atom mapping
3617!> \param reverse_phase use the reverse Bloch-phase convention
3618! **************************************************************************************************
3619 SUBROUTINE symtrans_phase(srpmat, scpmat, rpmat, cpmat, real_only, op, reverse_phase)
3620 TYPE(dbcsr_type), POINTER :: srpmat, scpmat, rpmat, cpmat
3621 LOGICAL, INTENT(IN) :: real_only
3622 TYPE(kp_symmetry_op_type), INTENT(IN) :: op
3623 LOGICAL, INTENT(IN) :: reverse_phase
3624
3625 CHARACTER(LEN=*), PARAMETER :: routinen = 'symtrans_phase'
3626
3627 INTEGER :: handle, icol, ip, irow, jcol, jp, jrow, &
3628 mynode, nthreads, numnodes, owner
3629 INTEGER, DIMENSION(3) :: shift
3630 LOGICAL :: found, trans
3631 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cwork, rwork, twork
3632 REAL(kind=dp), DIMENSION(:, :), POINTER :: cblock, left_rot, rblock, right_rot, &
3633 scblock, srblock
3634 TYPE(dbcsr_distribution_type) :: dist
3635 TYPE(dbcsr_iterator_type) :: iter
3636
3637 CALL timeset(routinen, handle)
3638
3639 IF (op%identity) THEN
3640 CALL dbcsr_copy(srpmat, rpmat)
3641 IF (.NOT. real_only) CALL dbcsr_copy(scpmat, cpmat)
3642 CALL timestop(handle)
3643 RETURN
3644 END IF
3645
3646 CALL dbcsr_get_info(rpmat, distribution=dist)
3647 CALL dbcsr_distribution_get(dist, mynode=mynode, numnodes=numnodes)
3648 IF (numnodes /= 1 .AND. op%redistribute) THEN
3649 CALL dbcsr_replicate_all(rpmat)
3650 IF (.NOT. real_only) CALL dbcsr_replicate_all(cpmat)
3651 END IF
3652
3653 CALL dbcsr_set(srpmat, 0.0_dp)
3654 IF (.NOT. real_only) CALL dbcsr_set(scpmat, 0.0_dp)
3655
3656 nthreads = 1
3657!$ nthreads = omp_get_max_threads()
3658 ! Atom permutations preserve row ownership, while block scheduling can race on a target row.
3659!$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads) &
3660!$OMP SHARED(rpmat,cpmat,srpmat,scpmat,real_only,op,reverse_phase,mynode) &
3661!$OMP PRIVATE(iter,irow,icol,rblock,cblock,rwork,cwork,twork,left_rot,right_rot,shift, &
3662!$OMP found,ip,jp,jrow,jcol,trans,srblock,scblock,owner)
3663 CALL dbcsr_iterator_readonly_start(iter, rpmat, dynamic=.true., dynamic_byrows=.true.)
3664 DO WHILE (dbcsr_iterator_blocks_left(iter))
3665 CALL dbcsr_iterator_next_block(iter, irow, icol, rblock)
3666 shift(:) = op%cell_shift(:, icol) - op%cell_shift(:, irow)
3667 IF (reverse_phase) shift(:) = -shift
3668 IF (real_only) THEN
3669 CALL kp_symmetry_phase_block(rblock, shift, op%xkp, op%time_reversal, rwork)
3670 ELSE
3671 CALL dbcsr_get_readonly_block_p(matrix=cpmat, row=irow, col=icol, block=cblock, found=found)
3672 IF (.NOT. found) NULLIFY (cblock)
3673 CALL kp_symmetry_phase_block(rblock, shift, op%xkp, op%time_reversal, rwork, cwork, cblock)
3674 END IF
3675
3676 ip = op%atom_map(irow)
3677 jp = op%atom_map(icol)
3678 jrow = min(ip, jp)
3679 jcol = max(ip, jp)
3680 trans = ip > jp
3681 IF (trans) THEN
3682 left_rot => op%kind_rot(op%atom_kind(icol))%rmat
3683 right_rot => op%kind_rot(op%atom_kind(irow))%rmat
3684 ELSE
3685 left_rot => op%kind_rot(op%atom_kind(irow))%rmat
3686 right_rot => op%kind_rot(op%atom_kind(icol))%rmat
3687 END IF
3688
3689 CALL dbcsr_get_block_p(matrix=srpmat, row=jrow, col=jcol, block=srblock, found=found)
3690 IF (.NOT. found) THEN
3691 CALL dbcsr_get_stored_coordinates(srpmat, jrow, jcol, owner)
3692 cpassert(owner /= mynode)
3693 cycle
3694 END IF
3695 CALL kp_symmetry_rotate_block(rwork, left_rot, right_rot, trans, 1.0_dp, srblock, twork)
3696 IF (.NOT. real_only) THEN
3697 CALL dbcsr_get_block_p(matrix=scpmat, row=jrow, col=jcol, block=scblock, found=found)
3698 cpassert(found)
3699 ! Returning to the stored upper triangle conjugates the imaginary part.
3700 CALL kp_symmetry_rotate_block(cwork, left_rot, right_rot, trans, &
3701 merge(-1.0_dp, 1.0_dp, trans), scblock, twork)
3702 END IF
3703 END DO
3704 CALL dbcsr_iterator_stop(iter)
3705!$OMP END PARALLEL
3706 IF (numnodes /= 1 .AND. op%redistribute) THEN
3707 CALL dbcsr_distribute(rpmat)
3708 IF (.NOT. real_only) CALL dbcsr_distribute(cpmat)
3709 END IF
3710
3711 CALL timestop(handle)
3712
3713 END SUBROUTINE symtrans_phase
3714
3715! **************************************************************************************************
3716!> \brief Real two-component representation of the Bloch phase and time reversal.
3717!> \param shift column-minus-row cell shift with the selected phase direction
3718!> \param xkp target k-point coordinates
3719!> \param time_reversal conjugate the source before applying the phase
3720!> \param real_only reject phases not representable by real wavefunctions
3721!> \param phase coefficients mapping source [Re,Im] to transformed [Re,Im]
3722! **************************************************************************************************
3723 SUBROUTINE kp_symmetry_phase_factors(shift, xkp, time_reversal, real_only, phase)
3724 INTEGER, INTENT(IN) :: shift(3)
3725 REAL(kind=dp), INTENT(IN) :: xkp(3)
3726 LOGICAL, INTENT(IN) :: time_reversal, real_only
3727 REAL(kind=dp), INTENT(OUT) :: phase(2, 2)
3728
3729 REAL(kind=dp) :: arg, cosine, sine
3730
3731 arg = real(shift(1), dp)*xkp(1) + real(shift(2), dp)*xkp(2) + real(shift(3), dp)*xkp(3)
3732 cosine = cos(twopi*arg)
3733 sine = sin(twopi*arg)
3734 IF (real_only .AND. abs(sine) > 1.e-12_dp) THEN
3735 CALL cp_abort(__location__, "Real k-point wavefunctions cannot represent symmetry phases")
3736 END IF
3737 phase(:, 1) = [cosine, -sine]
3738 phase(:, 2) = merge(-1.0_dp, 1.0_dp, time_reversal)*[sine, cosine]
3739 END SUBROUTINE kp_symmetry_phase_factors
3740
3741! **************************************************************************************************
3742!> \brief Apply the Bloch phase and optional time reversal to one AO block.
3743!> \param rblock real part of source block
3744!> \param shift column-minus-row cell shift, with the selected phase direction
3745!> \param xkp target k-point coordinates
3746!> \param time_reversal conjugate the source before applying the phase
3747!> \param rwork real part of the phased block; caller-owned reusable workspace
3748!> \param cwork imaginary workspace; absent for real-only wavefunctions
3749!> \param cblock imaginary source block; absent when the sparse block is missing
3750! **************************************************************************************************
3751 SUBROUTINE kp_symmetry_phase_block(rblock, shift, xkp, time_reversal, rwork, cwork, cblock)
3752 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: rblock
3753 INTEGER, DIMENSION(3), INTENT(IN) :: shift
3754 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: xkp
3755 LOGICAL, INTENT(IN) :: time_reversal
3756 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
3757 INTENT(INOUT) :: rwork
3758 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
3759 INTENT(INOUT), OPTIONAL :: cwork
3760 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
3761 OPTIONAL :: cblock
3762
3763 REAL(kind=dp) :: phase(2, 2)
3764
3765 CALL kp_symmetry_phase_factors(shift, xkp, time_reversal,.NOT. PRESENT(cwork), phase)
3766 CALL ensure_work_matrix(rwork, SIZE(rblock, 1), SIZE(rblock, 2))
3767 rwork(:, :) = phase(1, 1)*rblock
3768 IF (.NOT. PRESENT(cwork)) RETURN
3769 CALL ensure_work_matrix(cwork, SIZE(rblock, 1), SIZE(rblock, 2))
3770 cwork(:, :) = phase(2, 1)*rblock
3771 IF (PRESENT(cblock)) THEN
3772 rwork(:, :) = rwork + phase(1, 2)*cblock
3773 cwork(:, :) = cwork + phase(2, 2)*cblock
3774 END IF
3775
3776 END SUBROUTINE kp_symmetry_phase_block
3777
3778! **************************************************************************************************
3779!> \brief Accumulate left_rot * op(block) * right_rot^T into one AO block.
3780!> \param block phased source block
3781!> \param left_rot rotation for the target row
3782!> \param right_rot rotation for the target column
3783!> \param trans transpose the source when the atom mapping reverses block order
3784!> \param alpha sign for real symmetric or imaginary antisymmetric storage
3785!> \param TARGET destination block, updated in place
3786!> \param work caller-owned reusable GEMM workspace
3787! **************************************************************************************************
3788 SUBROUTINE kp_symmetry_rotate_block(block, left_rot, right_rot, trans, alpha, TARGET, work)
3789 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: block, left_rot, right_rot
3790 LOGICAL, INTENT(IN) :: trans
3791 REAL(kind=dp), INTENT(IN) :: alpha
3792 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: target
3793 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
3794 INTENT(INOUT) :: work
3795
3796 INTEGER :: ncol
3797
3798 ncol = merge(SIZE(block, 1), SIZE(block, 2), trans)
3799 CALL ensure_work_matrix(work, SIZE(left_rot, 1), ncol)
3800 CALL dgemm('N', merge('T', 'N', trans), SIZE(left_rot, 1), ncol, SIZE(left_rot, 2), &
3801 1.0_dp, left_rot, SIZE(left_rot, 1), block, SIZE(block, 1), &
3802 0.0_dp, work, SIZE(work, 1))
3803 CALL dgemm('N', 'T', SIZE(work, 1), SIZE(right_rot, 1), SIZE(work, 2), &
3804 alpha, work, SIZE(work, 1), right_rot, SIZE(right_rot, 1), &
3805 1.0_dp, TARGET, SIZE(TARGET, 1))
3806
3807 END SUBROUTINE kp_symmetry_rotate_block
3808
3809END MODULE kpoint_methods
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
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:57
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
memory_utilities::reallocate
Definition memory_utilities.F:27
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
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
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::real_to_scaled
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
Definition cell_types.F:595
cell_types::pbc_stable
real(kind=dp) function, dimension(3), public pbc_stable(r, cell)
Apply a stable periodic-image convention for k-point Bloch gauges.
Definition cell_types.F:422
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_blacs_env::cp_blacs_env_create
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Definition cp_blacs_env.F:123
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_scale_and_add
subroutine, public cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b).
Definition cp_cfm_basic_linalg.F:243
cp_cfm_basic_linalg::cp_cfm_transpose
subroutine, public cp_cfm_transpose(matrix, trans, matrixt)
Transposes a BLACS distributed complex matrix.
Definition cp_cfm_basic_linalg.F:1110
cp_cfm_basic_linalg::cp_cfm_column_scale
subroutine, public cp_cfm_column_scale(matrix_a, scaling)
Scales columns of the full matrix by corresponding factors.
Definition cp_cfm_basic_linalg.F:632
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_get_submatrix
subroutine, public cp_cfm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
Extract a sub-matrix from the full matrix: op(target_m)(1:n_rows,1:n_cols) = fm(start_row:start_row+n...
Definition cp_cfm_types.F:357
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:183
cp_cfm_types::cp_fm_to_cfm
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
Definition cp_cfm_types.F:845
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:126
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:632
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition cp_cfm_types.F:793
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_get_readonly_block_p
subroutine, public dbcsr_get_readonly_block_p(matrix, row, col, block, found, row_size, col_size)
Like dbcsr_get_block_p() but with matrix being INTENT(IN). When invoking this routine,...
Definition cp_dbcsr_api.F:718
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
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_api::dbcsr_replicate_all
subroutine, public dbcsr_replicate_all(matrix)
...
Definition cp_dbcsr_api.F:1156
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_get_stored_coordinates
subroutine, public dbcsr_get_stored_coordinates(matrix, row, column, processor)
...
Definition cp_dbcsr_api.F:911
cp_dbcsr_api::dbcsr_distribute
subroutine, public dbcsr_distribute(matrix)
...
Definition cp_dbcsr_api.F:540
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_iterator_readonly_start
subroutine, public dbcsr_iterator_readonly_start(iterator, matrix, shared, dynamic, dynamic_byrows)
Like dbcsr_iterator_start() but with matrix being INTENT(IN). When invoking this routine,...
Definition cp_dbcsr_api.F:1033
cp_dbcsr_api::dbcsr_distribution_get
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
Definition cp_dbcsr_api.F:572
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:129
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition cp_fm_basic_linalg.F:1690
cp_fm_pool_types
pool for for elements that are retained and released
Definition cp_fm_pool_types.F:14
cp_fm_pool_types::fm_pool_create_fm
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
Definition cp_fm_pool_types.F:185
cp_fm_pool_types::fm_pool_give_back_fm
subroutine, public fm_pool_give_back_fm(pool, element)
returns the element to the pool
Definition cp_fm_pool_types.F:224
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_equivalent
logical function, public cp_fm_struct_equivalent(fmstruct1, fmstruct2)
returns true if the two matrix structures are equivalent, false otherwise.
Definition cp_fm_struct.F:392
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
Definition cp_fm_types.F:564
cp_fm_types::cp_fm_start_copy_general
subroutine, public cp_fm_start_copy_general(source, destination, para_env, info)
Initiates the copy operation: get distribution data, post MPI isend and irecvs.
Definition cp_fm_types.F:1554
cp_fm_types::cp_fm_cleanup_copy_general
subroutine, public cp_fm_cleanup_copy_general(info)
Completes the copy operation: wait for comms clean up MPI state.
Definition cp_fm_types.F:1929
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_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_fm_types::cp_fm_finish_copy_general
subroutine, public cp_fm_finish_copy_general(destination, info)
Completes the copy operation: wait for comms, unpack, clean up MPI state.
Definition cp_fm_types.F:1865
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:925
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cryssym
K-points and crystal symmetry routines.
Definition cryssym.F:12
cryssym::print_crys_symmetry
subroutine, public print_crys_symmetry(csym)
...
Definition cryssym.F:1856
cryssym::kpoint_gen
subroutine, public kpoint_gen(csym, nk, symm, shift, full_grid, gamma_centered, inversion_symmetry_only, use_spglib_reduction, use_spglib_backend)
...
Definition cryssym.F:252
cryssym::release_csym_type
subroutine, public release_csym_type(csym)
Release the CSYM type.
Definition cryssym.F:89
cryssym::kpoint_gen_general
subroutine, public kpoint_gen_general(csym, xkp_in, wkp_in, symm, full_grid, inversion_symmetry_only, use_spglib_reduction, use_spglib_backend)
Reduce an explicitly supplied GENERAL k-point set.
Definition cryssym.F:459
cryssym::print_kp_symmetry
subroutine, public print_kp_symmetry(csym)
...
Definition cryssym.F:1893
cryssym::crys_sym_gen
subroutine, public crys_sym_gen(csym, scoor, types, hmat, delta, iounit, use_spglib)
...
Definition cryssym.F:144
hairy_probes
Definition hairy_probes.F:8
hairy_probes::probe_occupancy_kp
subroutine, public probe_occupancy_kp(occ, fermi, kts, energies, rcoeff, icoeff, maxocc, probe, n, wk)
subroutine to calculate occupation number and 'Fermi' level using the
Definition hairy_probes.F:194
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::smear_fermi_dirac
integer, parameter, public smear_fermi_dirac
Definition input_constants.F:950
input_constants::smear_gaussian
integer, parameter, public smear_gaussian
Definition input_constants.F:950
input_constants::smear_mv
integer, parameter, public smear_mv
Definition input_constants.F:950
input_constants::smear_mp
integer, parameter, public smear_mp
Definition input_constants.F:950
input_cp2k_kpoints
function that build the kpoints section of the input
Definition input_cp2k_kpoints.F:14
input_cp2k_kpoints::use_spglib_kpoint_symmetry
integer, parameter, public use_spglib_kpoint_symmetry
Definition input_cp2k_kpoints.F:37
input_cp2k_kpoints::lattice_fft_auto
integer, parameter, public lattice_fft_auto
Definition input_cp2k_kpoints.F:37
input_cp2k_kpoints::lattice_fft_off
integer, parameter, public lattice_fft_off
Definition input_cp2k_kpoints.F:37
input_cp2k_kpoints::use_real_wfn
integer, parameter, public use_real_wfn
Definition input_cp2k_kpoints.F:37
input_cp2k_kpoints::use_spglib_kpoint_backend
integer, parameter, public use_spglib_kpoint_backend
Definition input_cp2k_kpoints.F:37
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kpoint_lattice_fft
Batched lattice Fourier transforms on regular k-point grids.
Definition kpoint_lattice_fft.F:16
kpoint_lattice_fft::k_grid_to_cell_prepare
subroutine, public k_grid_to_cell_prepare(work, xkp, nkp_grid, index_to_cell, weights, allow_incomplete)
Prepare one inverse lattice transform for repeated matrix batches.
Definition kpoint_lattice_fft.F:326
kpoint_lattice_fft::regular_kpoint_grid
logical function, public regular_kpoint_grid(xkp, nkp_grid, grid_index, k_offset, allow_incomplete)
Test and map a uniformly shifted reciprocal grid or an explicitly allowed subset.
Definition kpoint_lattice_fft.F:67
kpoint_lattice_fft::k_grid_to_cell_execute
subroutine, public k_grid_to_cell_execute(work, values_k, values_rs, used_fft)
Execute a matrix batch using the prepared inverse lattice transform. Call outside OpenMP worker regio...
Definition kpoint_lattice_fft.F:377
kpoint_lattice_fft::k_grid_to_cell_release
subroutine, public k_grid_to_cell_release(work)
Release call-scoped inverse lattice maps and FFT buffers.
Definition kpoint_lattice_fft.F:429
kpoint_lattice_fft::cell_to_k_grid_fft
subroutine, public cell_to_k_grid_fft(values_rs, index_to_cell, xkp, nkp_grid, values_k, used_fft, deriv_direction, hmat, selected_kpoints)
Transform a batch of real matrices from image cells to every supplied k point.
Definition kpoint_lattice_fft.F:171
kpoint_lattice_fft::lattice_fft_shape
logical function, public lattice_fft_shape(n, nfft, allow_arbitrary)
Find backend-supported FFT dimensions that preserve the original lattice grid. Sample the padded tran...
Definition kpoint_lattice_fft.F:126
kpoint_methods
Routines needed for kpoint calculation.
Definition kpoint_methods.F:15
kpoint_methods::kpoint_smearing_edge_status
subroutine, public kpoint_smearing_edge_status(wocc, wkp, smear, nspin, has_weight, first_fractional, last_occupied, last_occupied_spin)
summarize whether weighted k-point smearing reaches the available band edges
Definition kpoint_methods.F:1762
kpoint_methods::kpoint_density_transform
subroutine, public kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwork, for_aux_fit, pmat_ext, overlap_rs)
generate real space density matrices in DBCSR format
Definition kpoint_methods.F:2220
kpoint_methods::lowdin_kp_mo_coeff
subroutine, public lowdin_kp_mo_coeff(kp, ispin, use_real_wfn, shalfc, cshalfc)
Calculate S(k)^1/2 C(k) for real or complex k-point wavefunctions.
Definition kpoint_methods.F:2178
kpoint_methods::kpoint_initialize_mo_set
subroutine, public kpoint_initialize_mo_set(kpoint)
...
Definition kpoint_methods.F:976
kpoint_methods::rskp_transform
subroutine, public rskp_transform(rmatrix, cmatrix, rsmat, ispin, xkp, cell_to_index, sab_nl, is_complex, rs_sign)
Transformation of real space matrices to a kpoint.
Definition kpoint_methods.F:1145
kpoint_methods::kpoint_init_cell_index
subroutine, public kpoint_init_cell_index(kpoint, sab_nl, para_env, nimages)
Generates the mapping of cell indices and linear RS index CELL (0,0,0) is always mapped to index 1.
Definition kpoint_methods.F:1011
kpoint_methods::kp_link_index_create
subroutine kp_link_index_create(representative, nkp, start, order)
Index full-grid links by irreducible k-point without rescanning the complete mesh.
Definition kpoint_methods.F:182
kpoint_methods::kpoint_set_mo_occupation
subroutine, public kpoint_set_mo_occupation(kpoint, smear, probe, added_mos_auto, added_mos_auto_grow, separate_spin_occupations)
Given the eigenvalues of all kpoints, calculates the occupation numbers.
Definition kpoint_methods.F:1506
kpoint_methods::rskp_transform_grid_extract
subroutine, public rskp_transform_grid_extract(grid, ikp, rmatrix, cmatrix)
Extract one reciprocal-grid matrix from a prepared local DBCSR block cache.
Definition kpoint_methods.F:1445
kpoint_methods::rskp_transform_grid_release
subroutine, public rskp_transform_grid_release(grid)
Release a reciprocal-grid transformation cache.
Definition kpoint_methods.F:1483
kpoint_methods::kpoint_initialize_mos
subroutine, public kpoint_initialize_mos(kpoint, mos, added_mos, for_aux_fit)
Initialize a set of MOs and density matrix for each kpoint (kpoint group).
Definition kpoint_methods.F:825
kpoint_methods::kp_transform_plan_create
subroutine, public kp_transform_plan_create(plan, sab_nl, cell_to_index, nimg, block_template, group_entries)
Build the immutable neighbor-list traversal shared by R-to-K and K-to-R transforms.
Definition kpoint_methods.F:3261
kpoint_methods::kpoint_initialize
subroutine, public kpoint_initialize(kpoint, particle_set, cell)
Generate the kpoints and initialize the kpoint environment.
Definition kpoint_methods.F:382
kpoint_methods::kpoint_ot_energy_weighted_density
subroutine, public kpoint_ot_energy_weighted_density(coeff, hc, occupation, wmat)
Build W(k) for noncanonical complex OT orbitals from H(k) C(k). The occupied-space Lagrange multiplie...
Definition kpoint_methods.F:1957
kpoint_methods::rskp_transform_grid_prepare
subroutine, public rskp_transform_grid_prepare(grid, rmatrix, rsmat, ispin, xkp, nkp_grid, cell_to_index, sab_nl, used_fft, is_complex, rs_sign, max_storage_bytes)
Prepare a batched real-cell to complete reciprocal-grid transform for local DBCSR blocks....
Definition kpoint_methods.F:1257
kpoint_methods::kpoint_density_matrices
subroutine, public kpoint_density_matrices(kpoint, energy_weighted, for_aux_fit)
Calculate kpoint density matrices (rho(k), owned by kpoint groups).
Definition kpoint_methods.F:1859
kpoint_methods::lowdin_kp_trans
subroutine, public lowdin_kp_trans(kpoint, pmat_diag)
Calculate Lowdin transformation of density matrix S^1/2 P S^1/2 Integrate diagonal elements over k-po...
Definition kpoint_methods.F:2091
kpoint_methods::kpoint_env_initialize
subroutine, public kpoint_env_initialize(kpoint, para_env, blacs_env, with_aux_fit)
Initialize the kpoint environment.
Definition kpoint_methods.F:661
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::kpoint_sym_create
subroutine, public kpoint_sym_create(kp_sym)
Create a single kpoint symmetry environment.
Definition kpoint_types.F:1197
kpoint_types::find_kpoint_rotation_slot
integer function, public find_kpoint_rotation_slot(qs_kpoint, rotp)
Locate the basis-rotation slot corresponding to a signed k-point symmetry operation.
Definition kpoint_types.F:1266
kpoint_types::kpoint_env_create
subroutine, public kpoint_env_create(kp_env)
Create a single kpoint environment.
Definition kpoint_types.F:1101
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:587
kpsym
Crystal-symmetry routines originating from the K290/ACMI code.
Definition kpsym.F:14
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
local_gemm_api
Definition local_gemm_api.F:8
local_gemm_api::local_gemm_pu_gpu
integer, parameter, public local_gemm_pu_gpu
Definition local_gemm_api.F:41
local_gemm_api::local_gemm_pu_host
integer, parameter, public local_gemm_pu_host
Definition local_gemm_api.F:41
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:144
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:144
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::inv_3x3
pure real(kind=dp) function, dimension(3, 3), public inv_3x3(a)
Returns the inverse of the 3 x 3 matrix a.
Definition mathlib.F:524
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:27
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_matrix_pools
wrapper for the pools of matrixes
Definition qs_matrix_pools.F:14
qs_matrix_pools::mpools_create
subroutine, public mpools_create(mpools)
creates a mpools
Definition qs_matrix_pools.F:184
qs_matrix_pools::mpools_rebuild_fm_pools
subroutine, public mpools_rebuild_fm_pools(mpools, mos, blacs_env, para_env, nmosub)
rebuilds the pools of the (ao x mo, ao x ao , mo x mo) full matrixes
Definition qs_matrix_pools.F:209
qs_matrix_pools::mpools_get
subroutine, public mpools_get(mpools, ao_mo_fm_pools, ao_ao_fm_pools, mo_mo_fm_pools, ao_mosub_fm_pools, mosub_mosub_fm_pools, maxao_maxmo_fm_pool, maxao_maxao_fm_pool, maxmo_maxmo_fm_pool)
returns various attributes of the mpools (notably the pools contained in it)
Definition qs_matrix_pools.F:139
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::set_mo_set
subroutine, public set_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, uniform_occupation, kts, mu, flexible_electron_count)
Set the components of a MO set data structure.
Definition qs_mo_types.F:502
qs_mo_types::init_mo_set
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name, counter, complex_coeff)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
Definition qs_mo_types.F:266
qs_mo_types::allocate_mo_set
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
Definition qs_mo_types.F:224
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count, cmo_coeff)
Get the components of a MO set data structure.
Definition qs_mo_types.F:445
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::get_neighbor_list_set_p
subroutine, public get_neighbor_list_set_p(neighbor_list_sets, nlist, symmetric)
Return the components of the first neighbor list set.
Definition qs_neighbor_list_types.F:996
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:550
scf_control_types
parameters that control an scf iteration
Definition scf_control_types.F:17
smearing_utils
Unified smearing module supporting four methods: smear_fermi_dirac — Fermi-Dirac distribution smear_g...
Definition smearing_utils.F:24
smearing_utils::smearkp
subroutine, public smearkp(f, mu, kts, e, nel, wk, sigma, maxocc, method)
Bisection search for mu given a target electron count (k-point case, single spin channel or spin-dege...
Definition smearing_utils.F:499
smearing_utils::smearkp2
subroutine, public smearkp2(f, mu, kts, e, nel, wk, sigma, method)
Bisection search for mu (k-point, spin-polarised with a shared chemical potential across both spin ch...
Definition smearing_utils.F:643
util
All kind of helpful little routines.
Definition util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition util.F:333
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:70
cp_control_types::hairy_probes_type
Definition cp_control_types.F:49
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
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_pool_types::cp_fm_pool_p_type
to create arrays of pools
Definition cp_fm_pool_types.F:70
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::copy_info_type
Stores the state of a copy between cp_fm_start_copy_general and cp_fm_finish_copy_general.
Definition cp_fm_types.F:138
cp_fm_types::cp_fm_p_type
just to build arrays of pointers to matrices
Definition cp_fm_types.F:128
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cryssym::csym_type
CSM type.
Definition cryssym.F:43
kpoint_lattice_fft::k_grid_to_cell_work_type
Definition kpoint_lattice_fft.F:45
kpoint_methods::kp_transform_plan_type
Definition kpoint_methods.F:127
kpoint_methods::rskp_grid_type
Definition kpoint_methods.F:164
kpoint_types::kind_rotmat_type
Rotation matrices for basis sets.
Definition kpoint_types.F:128
kpoint_types::kpoint_env_p_type
Definition kpoint_types.F:119
kpoint_types::kpoint_env_type
Keeps information about a specific k-point.
Definition kpoint_types.F:94
kpoint_types::kpoint_sym_type
Keeps symmetry information about a specific k-point.
Definition kpoint_types.F:145
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:198
local_gemm_api::local_gemm_ctxt_type
Definition local_gemm_api.F:47
message_passing::mp_cart_type
Definition message_passing.F:745
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
particle_types::particle_type
Definition particle_types.F:35
qs_matrix_pools::qs_matrix_pools_type
container for the pools of matrixes used by qs
Definition qs_matrix_pools.F:66
qs_mo_types::mo_set_type
Definition qs_mo_types.F:53
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
scf_control_types::smear_type
contains the parameters needed by a scf run
Definition scf_control_types.F:98