Line data Source code
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 Assembly of real and complex k-point operators from real-space DBCSR matrices.
10 : !> Outputs use the distribution of the k-point group; real operators transfer one part only.
11 : !> The module splits this work into a start step and a finish step.
12 : !> One process can start many operators before it finishes any of them.
13 : !> The STANDARD k-point driver uses this order.
14 : !> The routines start, finish and get are collective over the
15 : !> environment-wide communicator. Every rank must call the same
16 : !> sequence of service calls.
17 : !> A context references its host kpoint environment without owning it,
18 : !> so it must not outlive that environment.
19 : !> An explicitly prepared, immutable operator batch can instead be built
20 : !> from bounded context-owned real-space snapshots inside each group.
21 : !> The public round and output-ownership contracts stay the same.
22 : !> Start/finish/get/prepare are not thread safe: serialize them, because their collective
23 : !> communication tolerates no concurrent service call. A separate
24 : !> context per thread only keeps the scratch from aliasing.
25 : ! **************************************************************************************************
26 : MODULE qs_kpoint_operators
27 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale_and_add_fm
28 : USE cp_cfm_types, ONLY: cp_cfm_to_fm,&
29 : cp_cfm_type
30 : USE cp_dbcsr_api, ONLY: &
31 : dbcsr_create, dbcsr_deallocate_matrix, dbcsr_get_info, dbcsr_get_readonly_block_p, &
32 : dbcsr_init_p, dbcsr_p_type, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
33 : dbcsr_type_symmetric
34 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
35 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
36 : copy_fm_to_dbcsr_bc,&
37 : dbcsr_fm_copy_plan_type
38 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
39 : fm_pool_create_fm,&
40 : fm_pool_give_back_fm
41 : USE cp_fm_struct, ONLY: cp_fm_struct_equivalent,&
42 : cp_fm_struct_release,&
43 : cp_fm_struct_retain,&
44 : cp_fm_struct_type
45 : USE cp_fm_types, ONLY: copy_info_type,&
46 : cp_fm_cleanup_copy_general,&
47 : cp_fm_create,&
48 : cp_fm_finish_copy_general,&
49 : cp_fm_get_info,&
50 : cp_fm_release,&
51 : cp_fm_start_copy_general,&
52 : cp_fm_type
53 : USE input_cp2k_kpoints, ONLY: lattice_fft_auto,&
54 : lattice_fft_on
55 : USE kinds, ONLY: dp,&
56 : int_8
57 : USE kpoint_methods, ONLY: kp_transform_plan_create,&
58 : kp_transform_plan_type,&
59 : rskp_grid_type,&
60 : rskp_transform,&
61 : rskp_transform_grid_extract,&
62 : rskp_transform_grid_prepare,&
63 : rskp_transform_grid_release
64 : USE kpoint_types, ONLY: kpoint_env_type,&
65 : kpoint_type
66 : USE mathconstants, ONLY: gaussi,&
67 : twopi,&
68 : z_one,&
69 : z_zero
70 : USE message_passing, ONLY: mp_para_env_type
71 : USE qs_matrix_pools, ONLY: mpools_get
72 : USE qs_mo_types, ONLY: get_mo_set,&
73 : mo_set_type
74 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
75 :
76 : !$ USE omp_lib, ONLY: omp_get_max_threads, omp_get_num_threads, omp_get_thread_num
77 : #include "./base/base_uses.f90"
78 :
79 : IMPLICIT NONE
80 :
81 : PRIVATE
82 :
83 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_kpoint_operators'
84 :
85 : PUBLIC :: kpoint_operator_cfm_to_mo
86 : PUBLIC :: kpoint_operator_context_create, kpoint_operator_context_release
87 : PUBLIC :: kpoint_operator_start, kpoint_operator_finish
88 : PUBLIC :: kpoint_operator_get, kpoint_operator_prepare
89 : PUBLIC :: kpoint_operator_evaluate, kpoint_operator_work_size
90 : PUBLIC :: kpoint_operator_store, kpoint_operator_get_local, kpoint_operator_get_all
91 :
92 : INTEGER, PARAMETER, PUBLIC :: kpoint_slot_ks = 1, &
93 : kpoint_slot_s = 2, &
94 : kpoint_slot_t = 3
95 :
96 : !> Row of rsmat that carries the spin-free operator image. S and T carry
97 : !> no spin copies: drivers launch one instance per spin of their loop and
98 : !> read this same row in every instance.
99 : INTEGER, PARAMETER, PUBLIC :: kpoint_spin_free = 1
100 :
101 : !> Number of defined slots; grows only inside this module.
102 : INTEGER, PARAMETER, PRIVATE :: kpoint_num_slots = kpoint_slot_t - kpoint_slot_ks + 1
103 :
104 : !> Part indices of one transfer: its real and its imaginary half.
105 : INTEGER, PARAMETER, PRIVATE :: kpoint_part_re = 1, &
106 : kpoint_part_im = 2
107 :
108 : !> Instance states of one (k point, spin, slot) transfer.
109 : !> State transitions happen only inside the service routines.
110 : INTEGER, PARAMETER, PRIVATE :: kpoint_slot_idle = 0, &
111 : kpoint_slot_in_flight = 1, &
112 : kpoint_slot_local = 2
113 :
114 : ! Structural route shared by all captured rows. Whole blocks go to one
115 : ! representative group, then only equivalent tiles are replicated.
116 : TYPE kpoint_operator_rs_route_type
117 : TYPE(kp_transform_plan_type) :: traversal = kp_transform_plan_type()
118 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct => NULL()
119 : INTEGER :: nimage = 0, nvalue = 0, nrow = 0, ncol = 0
120 : INTEGER(KIND=int_8) :: storage = 0
121 : INTEGER, ALLOCATABLE :: row_size(:), col_size(:), ENTRY(:, :), first(:), tile(:, :, :)
122 : INTEGER, ALLOCATABLE :: send_entry(:), send_first(:), scount(:, :), rcount(:, :), &
123 : sdisp(:, :), rdisp(:, :)
124 : REAL(KIND=dp), ALLOCATABLE :: send_values(:)
125 : END TYPE kpoint_operator_rs_route_type
126 :
127 : ! Immutable values, not per-k matrices. Metadata belong to the shared route;
128 : ! source identifies starts belonging to the explicitly prepared batch only.
129 : TYPE kpoint_operator_block_type
130 : REAL(KIND=dp), POINTER :: data(:, :) => NULL()
131 : END TYPE kpoint_operator_block_type
132 :
133 : TYPE kpoint_operator_rs_type
134 : INTEGER :: matrix_row = 0
135 : TYPE(dbcsr_p_type), POINTER :: source(:, :) => NULL()
136 : ! Values are owned only for distributed sources. Blocks always borrow:
137 : ! either from values below or, on one world rank, directly from DBCSR.
138 : REAL(KIND=dp), CONTIGUOUS, POINTER :: values(:) => NULL()
139 : TYPE(kpoint_operator_block_type), ALLOCATABLE :: blocks(:)
140 : LOGICAL, ALLOCATABLE :: nonzero(:)
141 : END TYPE kpoint_operator_rs_type
142 :
143 : ! Small read-ahead window per prepared row. Keep mutable output separate
144 : ! from the immutable snapshot; prepare invalidates only the selected slot.
145 : TYPE kpoint_operator_rs_batch_type
146 : INTEGER :: first = 0
147 : COMPLEX(KIND=dp), ALLOCATABLE :: data(:, :, :)
148 : END TYPE kpoint_operator_rs_batch_type
149 :
150 : INTEGER, PARAMETER :: rs_batch_size = 4
151 : ! Includes cached tiles and batch phase scratch, divided between all rows/slots.
152 : INTEGER(KIND=int_8), PARAMETER :: rs_batch_bytes = 32_int_8*1024_int_8**2
153 :
154 : ! **************************************************************************************************
155 : !> \brief Caller-owned assembly context for one driver call.
156 : !> Created and released in the same driver routine and passed to helpers
157 : !> as an argument; never stored on the host kpoint_type.
158 : !> Components are private and touched only by this module.
159 : ! **************************************************************************************************
160 : TYPE, PUBLIC :: kpoint_operator_context_type
161 : PRIVATE
162 : ! Instance spin space of the owning driver call: the bound of the spin
163 : ! loop that start and finish key. Declared once at creation. It never
164 : ! comes from the operators: the first matrix dimension of a spin-free
165 : ! operator such as S counts derivatives, and the images live in the
166 : ! second, so callers address operator rows through matrix_row instead.
167 : INTEGER :: nspin = 0
168 : ! Borrowed host/distribution references; the context cannot span an energy evaluation.
169 : TYPE(kpoint_type), POINTER :: kpoints => NULL()
170 : TYPE(cp_fm_struct_type), POINTER :: ao_ao_fmstruct => NULL()
171 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: group_fm_pools => NULL()
172 : ! Identity of the neighbor lists used by scratch and prepared snapshots.
173 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
174 : POINTER :: work_sab => NULL()
175 : ! Owned scratch, created on first use and released with the context.
176 : TYPE(dbcsr_type), POINTER :: rmatrix => NULL()
177 : TYPE(dbcsr_type), POINTER :: cmatrix => NULL()
178 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fmwork
179 : TYPE(cp_fm_type), ALLOCATABLE :: fmlocal, fmimag
180 : ! Plans share the lifetime and fixed layouts of the two densify workspaces.
181 : TYPE(dbcsr_fm_copy_plan_type), ALLOCATABLE :: copy_plan(:)
182 : ! Immutable snapshots and bounded local output batches. Prepared starts
183 : ! need neither per-k transfers nor outstanding MPI requests.
184 : TYPE(kpoint_operator_rs_type), ALLOCATABLE :: rs(:, :)
185 : TYPE(kpoint_operator_rs_route_type), ALLOCATABLE :: rs_route
186 : TYPE(kpoint_operator_rs_batch_type), ALLOCATABLE :: rs_batch(:, :)
187 : ! Bookkeeping. info is keyed (global k point, spin, slot, part),
188 : ! part kpoint_part_re holding the real half and kpoint_part_im holding
189 : ! the imaginary half. Status is keyed (global round, spin, slot)
190 : ! because start and finish treat all groups of one instance
191 : ! atomically. Allocated on first start.
192 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :, :, :) :: info
193 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: status
194 : END TYPE kpoint_operator_context_type
195 :
196 : CONTAINS
197 :
198 : ! **************************************************************************************************
199 : !> \brief Transform and densify one operator using the context-owned workspaces and routing plans.
200 : !> \param ctx initialized transfer context
201 : !> \param rsmat real-space operator images
202 : !> \param matrix_row operator row, independent of the instance spin
203 : !> \param ik global k-point index
204 : !> \param grid caller-owned reciprocal-grid cache
205 : !> \param use_grid extract from grid instead of a direct phase sum
206 : ! **************************************************************************************************
207 40 : SUBROUTINE kpoint_operator_densify(ctx, rsmat, matrix_row, ik, grid, use_grid)
208 :
209 : TYPE(kpoint_operator_context_type), INTENT(INOUT) :: ctx
210 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
211 : INTEGER, INTENT(IN) :: matrix_row, ik
212 : TYPE(rskp_grid_type), INTENT(IN), OPTIONAL :: grid
213 : LOGICAL, INTENT(IN) :: use_grid
214 :
215 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_operator_densify'
216 :
217 : INTEGER :: handle
218 40 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
219 :
220 40 : CALL timeset(routineN, handle)
221 :
222 : ! Avoid nested component array sections in GCC 16 bounds checks.
223 40 : xkp => ctx%kpoints%xkp
224 40 : IF (ctx%kpoints%use_real_wfn) THEN
225 16 : IF (use_grid) THEN
226 0 : CPASSERT(PRESENT(grid))
227 0 : CALL rskp_transform_grid_extract(grid, ik, ctx%rmatrix)
228 : ELSE
229 16 : CALL dbcsr_set(ctx%rmatrix, 0.0_dp)
230 : CALL rskp_transform(rmatrix=ctx%rmatrix, rsmat=rsmat, ispin=matrix_row, &
231 : xkp=xkp(:, ik), cell_to_index=ctx%kpoints%cell_to_index, &
232 16 : sab_nl=ctx%kpoints%sab_nl)
233 : END IF
234 24 : ELSE IF (use_grid) THEN
235 0 : CPASSERT(PRESENT(grid))
236 0 : CALL rskp_transform_grid_extract(grid, ik, ctx%rmatrix, ctx%cmatrix)
237 : ELSE
238 24 : CALL dbcsr_set(ctx%rmatrix, 0.0_dp)
239 24 : CALL dbcsr_set(ctx%cmatrix, 0.0_dp)
240 : CALL rskp_transform(rmatrix=ctx%rmatrix, cmatrix=ctx%cmatrix, rsmat=rsmat, ispin=matrix_row, &
241 : xkp=xkp(:, ik), cell_to_index=ctx%kpoints%cell_to_index, &
242 24 : sab_nl=ctx%kpoints%sab_nl)
243 : END IF
244 40 : CALL copy_dbcsr_to_fm(ctx%rmatrix, ctx%fmwork(kpoint_part_re), ctx%copy_plan(kpoint_part_re))
245 40 : IF (.NOT. ctx%kpoints%use_real_wfn) THEN
246 24 : CALL copy_dbcsr_to_fm(ctx%cmatrix, ctx%fmwork(kpoint_part_im), ctx%copy_plan(kpoint_part_im))
247 : END IF
248 :
249 40 : CALL timestop(handle)
250 :
251 40 : END SUBROUTINE kpoint_operator_densify
252 :
253 : ! **************************************************************************************************
254 : !> \brief Launch a densified operator to its owning group (only the real part for real wavefunctions).
255 : !> Absent explicit targets, receive both parts through fmlocal. Each copy
256 : !> owns its receive buffer until finish, so the targets can alias.
257 : !> Non-owning groups participate with an empty dummy target.
258 : !> \param kpoints host environment, including group ownership and transfer communicator
259 : !> \param ik global k-point index
260 : !> \param fmwork real/imaginary source matrices on the global distribution
261 : !> \param fmlocal default receiver on the k-point group distribution
262 : !> \param info transfer states, real part first; real wavefunctions allocate only the first part
263 : !> \param target_re optional caller-owned real-part target
264 : !> \param target_im optional imaginary-part target; required with target_re for complex wavefunctions
265 : ! **************************************************************************************************
266 40 : SUBROUTINE kpoint_operator_launch(kpoints, ik, fmwork, fmlocal, info, target_re, target_im)
267 :
268 : TYPE(kpoint_type), POINTER :: kpoints
269 : INTEGER, INTENT(IN) :: ik
270 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fmwork
271 : TYPE(cp_fm_type), INTENT(IN) :: fmlocal
272 : TYPE(copy_info_type), DIMENSION(:) :: info
273 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: target_re, target_im
274 :
275 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_operator_launch'
276 :
277 : INTEGER :: handle
278 : TYPE(cp_fm_type) :: fmdummy
279 : TYPE(mp_para_env_type), POINTER :: para_env
280 :
281 40 : CALL timeset(routineN, handle)
282 40 : CPASSERT(.NOT. PRESENT(target_im) .OR. PRESENT(target_re))
283 40 : IF (.NOT. kpoints%use_real_wfn) THEN
284 24 : CPASSERT(PRESENT(target_re) .EQV. PRESENT(target_im))
285 : END IF
286 :
287 40 : para_env => kpoints%blacs_env_all%para_env
288 40 : IF (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2)) THEN
289 40 : IF (PRESENT(target_re)) THEN
290 0 : CALL cp_fm_start_copy_general(fmwork(kpoint_part_re), target_re, para_env, info(kpoint_part_re))
291 0 : IF (.NOT. kpoints%use_real_wfn) THEN
292 0 : CALL cp_fm_start_copy_general(fmwork(kpoint_part_im), target_im, para_env, info(kpoint_part_im))
293 : END IF
294 : ELSE
295 40 : CALL cp_fm_start_copy_general(fmwork(kpoint_part_re), fmlocal, para_env, info(kpoint_part_re))
296 40 : IF (.NOT. kpoints%use_real_wfn) THEN
297 24 : CALL cp_fm_start_copy_general(fmwork(kpoint_part_im), fmlocal, para_env, info(kpoint_part_im))
298 : END IF
299 : END IF
300 : ELSE
301 0 : CALL cp_fm_start_copy_general(fmwork(kpoint_part_re), fmdummy, para_env, info(kpoint_part_re))
302 0 : IF (.NOT. kpoints%use_real_wfn) THEN
303 0 : CALL cp_fm_start_copy_general(fmwork(kpoint_part_im), fmdummy, para_env, info(kpoint_part_im))
304 : END IF
305 : END IF
306 :
307 40 : CALL timestop(handle)
308 :
309 40 : END SUBROUTINE kpoint_operator_launch
310 :
311 : ! **************************************************************************************************
312 : !> \brief Receive one group's operator and reclaim every rank's send buffers.
313 : !> CFM/DBCSR outputs consume the context buffer before receiving the next part;
314 : !> FM outputs receive directly. Inactive groups only clean up their sends.
315 : !> \param kpoints parallel environment
316 : !> \param fmlocal shared receive buffer, consumed before receiving the next part
317 : !> \param ik global k-point index
318 : !> \param info real/imaginary transfers started by kpoint_operator_launch
319 : !> \param cmat complex output, real part overwrites and imaginary part adds
320 : !> \param fm_re real FM output
321 : !> \param fm_im optional imaginary FM output
322 : !> \param mat_re real DBCSR output
323 : !> \param mat_im imaginary DBCSR output
324 : ! **************************************************************************************************
325 40 : SUBROUTINE kpoint_operator_receive(kpoints, fmlocal, ik, info, cmat, fm_re, fm_im, mat_re, mat_im)
326 : TYPE(kpoint_type), INTENT(IN), POINTER :: kpoints
327 : TYPE(cp_fm_type) :: fmlocal
328 : INTEGER, INTENT(IN) :: ik
329 : TYPE(copy_info_type), INTENT(INOUT) :: info(:)
330 : TYPE(cp_cfm_type), OPTIONAL :: cmat
331 : TYPE(cp_fm_type), OPTIONAL :: fm_re, fm_im
332 : TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL :: mat_re, mat_im
333 :
334 : INTEGER :: part
335 :
336 40 : IF (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2)) THEN
337 40 : IF (PRESENT(cmat)) THEN
338 24 : CALL cp_fm_finish_copy_general(fmlocal, info(kpoint_part_re))
339 24 : CALL cp_cfm_scale_and_add_fm(z_zero, cmat, z_one, fmlocal)
340 24 : IF (.NOT. kpoints%use_real_wfn) THEN
341 24 : CALL cp_fm_finish_copy_general(fmlocal, info(kpoint_part_im))
342 24 : CALL cp_cfm_scale_and_add_fm(z_one, cmat, gaussi, fmlocal)
343 : END IF
344 16 : ELSE IF (PRESENT(fm_re)) THEN
345 16 : CALL cp_fm_finish_copy_general(fm_re, info(kpoint_part_re))
346 16 : IF (kpoints%use_real_wfn) THEN
347 16 : IF (PRESENT(fm_im)) fm_im%local_data(:, :) = 0.0_dp
348 : ELSE
349 0 : CALL cp_fm_finish_copy_general(fm_im, info(kpoint_part_im))
350 : END IF
351 : ELSE
352 0 : CALL cp_fm_finish_copy_general(fmlocal, info(kpoint_part_re))
353 0 : CALL copy_fm_to_dbcsr_bc(fmlocal, mat_re)
354 0 : IF (kpoints%use_real_wfn) THEN
355 0 : fmlocal%local_data(:, :) = 0.0_dp
356 : ELSE
357 0 : CALL cp_fm_finish_copy_general(fmlocal, info(kpoint_part_im))
358 : END IF
359 0 : CALL copy_fm_to_dbcsr_bc(fmlocal, mat_im)
360 : END IF
361 : END IF
362 128 : DO part = 1, MERGE(1, 2, kpoints%use_real_wfn)
363 104 : CALL cp_fm_cleanup_copy_general(info(part))
364 : END DO
365 40 : END SUBROUTINE kpoint_operator_receive
366 :
367 : ! **************************************************************************************************
368 : !> \brief Split a complex matrix into the real and imaginary MO sets and copy
369 : !> the eigenvalues to the imaginary set.
370 : !> \param cmat complex MO coefficients
371 : !> \param mo_re real part MO set
372 : !> \param mo_im imaginary part MO set
373 : !> \param eigenvalues eigenvalues of mo_re, copied to mo_im
374 : ! **************************************************************************************************
375 68281 : SUBROUTINE kpoint_operator_cfm_to_mo(cmat, mo_re, mo_im, eigenvalues)
376 :
377 : TYPE(cp_cfm_type) :: cmat
378 : TYPE(mo_set_type) :: mo_re, mo_im
379 : REAL(KIND=dp), DIMENSION(:) :: eigenvalues
380 :
381 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_operator_cfm_to_mo'
382 :
383 : INTEGER :: handle
384 : TYPE(cp_fm_type), POINTER :: imos, rmos
385 :
386 68281 : CALL timeset(routineN, handle)
387 :
388 68281 : CALL get_mo_set(mo_re, mo_coeff=rmos)
389 68281 : CALL get_mo_set(mo_im, mo_coeff=imos)
390 1252635 : mo_im%eigenvalues(:) = eigenvalues
391 68281 : CALL cp_cfm_to_fm(cmat, rmos, imos)
392 :
393 68281 : CALL timestop(handle)
394 :
395 68281 : END SUBROUTINE kpoint_operator_cfm_to_mo
396 :
397 : ! **************************************************************************************************
398 : !> \brief Create an empty operator context.
399 : !> The context owns no buffers after this call; they are created on first
400 : !> use and released with the context.
401 : !> \param ctx context to initialize
402 : !> \param kpoints host k-point environment, referenced, not copied
403 : !> \param ao_ao_fm any square full matrix on the environment-wide communicator;
404 : !> its matrix struct is the template for the context work buffers,
405 : !> which densify on the source side of transfers that cross into the
406 : !> k-point groups
407 : !> \param nspin instance spin space of the owning driver call: the bound of
408 : !> the spin loop its start and finish calls iterate. Fixed for the
409 : !> lifetime of the context
410 : ! **************************************************************************************************
411 36524 : SUBROUTINE kpoint_operator_context_create(ctx, kpoints, ao_ao_fm, nspin)
412 :
413 : TYPE(kpoint_operator_context_type), INTENT(OUT) :: ctx
414 : TYPE(kpoint_type), POINTER :: kpoints
415 : TYPE(cp_fm_type), INTENT(IN) :: ao_ao_fm
416 : INTEGER, INTENT(IN) :: nspin
417 :
418 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_operator_context_create'
419 :
420 : INTEGER :: handle
421 :
422 36524 : CALL timeset(routineN, handle)
423 :
424 36524 : CPASSERT(nspin >= 1)
425 36524 : ctx%nspin = nspin
426 36524 : CPASSERT(ASSOCIATED(kpoints))
427 36524 : ctx%kpoints => kpoints
428 36524 : CALL cp_fm_get_info(ao_ao_fm, matrix_struct=ctx%ao_ao_fmstruct)
429 36524 : CPASSERT(ASSOCIATED(ctx%ao_ao_fmstruct))
430 : ! The template must live on the transfer world. A template from
431 : ! another communicator fails during the transfers, far from this call.
432 36524 : IF (.NOT. ASSOCIATED(ctx%ao_ao_fmstruct%para_env, kpoints%blacs_env_all%para_env)) THEN
433 0 : CALL cp_abort(__LOCATION__, "Operator context requires the global parallel environment")
434 : END IF
435 :
436 36524 : CALL timestop(handle)
437 :
438 36524 : END SUBROUTINE kpoint_operator_context_create
439 :
440 : ! **************************************************************************************************
441 : !> \brief Release the context together with every buffer it created.
442 : !> Call once, at the exit of the routine that created the context.
443 : !> Every started instance must be idle at this point. The call aborts
444 : !> when a transfer is still in flight, because nothing would ever wait
445 : !> for its MPI requests.
446 : !> \param ctx context to release
447 : ! **************************************************************************************************
448 36524 : SUBROUTINE kpoint_operator_context_release(ctx)
449 :
450 : TYPE(kpoint_operator_context_type), INTENT(INOUT) :: ctx
451 :
452 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_operator_context_release'
453 :
454 : INTEGER :: handle, ispin, iwork, slot
455 :
456 36524 : CALL timeset(routineN, handle)
457 :
458 : ! Every started instance must have been finished: an unfinished start
459 : ! holds MPI requests that nothing would ever wait for.
460 36524 : IF (ALLOCATED(ctx%status)) THEN
461 150450 : CPASSERT(ALL(ctx%status == kpoint_slot_idle))
462 : END IF
463 :
464 36524 : IF (ASSOCIATED(ctx%rmatrix)) CALL dbcsr_deallocate_matrix(ctx%rmatrix)
465 36524 : IF (ASSOCIATED(ctx%cmatrix)) CALL dbcsr_deallocate_matrix(ctx%cmatrix)
466 36524 : NULLIFY (ctx%rmatrix, ctx%cmatrix, ctx%work_sab)
467 36556 : IF (ALLOCATED(ctx%copy_plan)) DEALLOCATE (ctx%copy_plan)
468 36524 : IF (ALLOCATED(ctx%rs)) THEN
469 146096 : DO slot = 1, kpoint_num_slots
470 261956 : DO ispin = 1, ctx%nspin
471 225432 : CALL kpoint_operator_rs_clear(ctx%rs(ispin, slot))
472 : END DO
473 : END DO
474 268244 : DEALLOCATE (ctx%rs, ctx%rs_batch)
475 : END IF
476 36524 : IF (ALLOCATED(ctx%rs_route)) THEN
477 36512 : CALL cp_fm_struct_release(ctx%rs_route%matrix_struct)
478 36512 : DEALLOCATE (ctx%rs_route)
479 : END IF
480 36524 : IF (ALLOCATED(ctx%fmwork)) THEN
481 28 : DO iwork = 1, SIZE(ctx%fmwork)
482 28 : CALL cp_fm_release(ctx%fmwork(iwork))
483 : END DO
484 12 : DEALLOCATE (ctx%fmwork)
485 : END IF
486 36524 : IF (ALLOCATED(ctx%fmimag)) THEN
487 82 : CALL cp_fm_release(ctx%fmimag)
488 82 : DEALLOCATE (ctx%fmimag)
489 : END IF
490 36524 : IF (ALLOCATED(ctx%fmlocal)) THEN
491 36524 : CALL fm_pool_give_back_fm(ctx%group_fm_pools(1)%pool, ctx%fmlocal)
492 36524 : DEALLOCATE (ctx%fmlocal)
493 : END IF
494 36644 : IF (ALLOCATED(ctx%info)) DEALLOCATE (ctx%info)
495 36524 : IF (ALLOCATED(ctx%status)) DEALLOCATE (ctx%status)
496 36524 : NULLIFY (ctx%kpoints, ctx%ao_ao_fmstruct, ctx%group_fm_pools)
497 :
498 36524 : CALL timestop(handle)
499 :
500 36524 : END SUBROUTINE kpoint_operator_context_release
501 :
502 : ! **************************************************************************************************
503 : !> \brief Create the context scratch that the first service call needs: the
504 : !> DBCSR work pair with blocks from the current neighbor lists, two
505 : !> densify buffers on the environment-wide struct, and the group
506 : !> transfer target from the group pools.
507 : !> The per-instance slot bookkeeping is not scratch: start creates it
508 : !> on its first call, because the synchronous get call never touches it.
509 : !> The buffers assume the context lifetime of one driver call: the
510 : !> host neighbor lists do not change within it, and every buffer use
511 : !> asserts this identity. Buffers from an earlier call stay in place.
512 : !> Called from the service routines only.
513 : !> \param ctx context that receives the buffers
514 : !> \param rsmat real-space image matrices of one operator. The call uses the
515 : !> first matrix as the DBCSR template of the work pair. It reads no
516 : !> matrix data.
517 : ! **************************************************************************************************
518 40 : SUBROUTINE kpoint_operator_ensure(ctx, rsmat)
519 :
520 : TYPE(kpoint_operator_context_type), INTENT(INOUT) :: ctx
521 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
522 :
523 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_operator_ensure'
524 :
525 : INTEGER :: handle
526 :
527 40 : CALL timeset(routineN, handle)
528 :
529 : ! The pair must still match the host neighbor lists: they are
530 : ! reallocated at every energy evaluation, while the context spans one
531 : ! driver call within it.
532 40 : IF (ASSOCIATED(ctx%work_sab)) THEN
533 28 : CPASSERT(ASSOCIATED(ctx%work_sab, ctx%kpoints%sab_nl))
534 : END IF
535 40 : IF (.NOT. ASSOCIATED(ctx%rmatrix)) THEN
536 12 : ALLOCATE (ctx%rmatrix)
537 12 : CALL dbcsr_create(ctx%rmatrix, template=rsmat(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
538 12 : CALL cp_dbcsr_alloc_block_from_nbl(ctx%rmatrix, ctx%kpoints%sab_nl)
539 12 : IF (.NOT. ctx%kpoints%use_real_wfn) THEN
540 4 : ALLOCATE (ctx%cmatrix)
541 4 : CALL dbcsr_create(ctx%cmatrix, template=rsmat(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
542 4 : CALL cp_dbcsr_alloc_block_from_nbl(ctx%cmatrix, ctx%kpoints%sab_nl)
543 : END IF
544 12 : ctx%work_sab => ctx%kpoints%sab_nl
545 : END IF
546 :
547 40 : IF (.NOT. ALLOCATED(ctx%fmwork)) THEN
548 44 : ALLOCATE (ctx%fmwork(MERGE(1, 2, ctx%kpoints%use_real_wfn)))
549 52 : ALLOCATE (ctx%copy_plan(SIZE(ctx%fmwork)))
550 : CALL cp_fm_create(ctx%fmwork(kpoint_part_re), ctx%ao_ao_fmstruct, &
551 12 : name="KPOINT-OPERATOR-WORK-RE")
552 12 : IF (.NOT. ctx%kpoints%use_real_wfn) THEN
553 : CALL cp_fm_create(ctx%fmwork(kpoint_part_im), ctx%ao_ao_fmstruct, &
554 4 : name="KPOINT-OPERATOR-WORK-IM")
555 : END IF
556 : END IF
557 :
558 40 : CALL kpoint_operator_ensure_group(ctx)
559 :
560 40 : CALL timestop(handle)
561 :
562 40 : END SUBROUTINE kpoint_operator_ensure
563 :
564 : ! **************************************************************************************************
565 : !> \brief Acquire only the group-local output workspace; both backends share it.
566 : !> \param ctx operator context
567 : ! **************************************************************************************************
568 43986 : SUBROUTINE kpoint_operator_ensure_group(ctx)
569 :
570 : TYPE(kpoint_operator_context_type), INTENT(INOUT) :: ctx
571 :
572 43986 : IF (.NOT. ALLOCATED(ctx%fmlocal)) THEN
573 36524 : CALL mpools_get(ctx%kpoints%mpools, ao_ao_fm_pools=ctx%group_fm_pools)
574 36524 : ALLOCATE (ctx%fmlocal)
575 36524 : CALL fm_pool_create_fm(ctx%group_fm_pools(1)%pool, ctx%fmlocal)
576 : END IF
577 :
578 43986 : END SUBROUTINE kpoint_operator_ensure_group
579 :
580 : ! **************************************************************************************************
581 : !> \brief Prepare an immutable operator slot for group-local assembly.
582 : !> Refresh after the source values change, with every instance idle.
583 : !> Other prepared slots and structural routing remain valid. Without
584 : !> matrix_row, capture one row per instance spin; with matrix_row, store
585 : !> one row shared by all instance spins. Failure clears this slot on
586 : !> all ranks and retains the existing transfer fallback.
587 : !> \param ctx context spanning the current initialized k-point environment
588 : !> \param rsmat current real-space operator images
589 : !> \param slot slot to replace (KS, S or T)
590 : !> \param prepared optional collective success flag
591 : !> \param matrix_row optional explicit spin-free operator row
592 : !> \param TARGET optional group-local FM layout for prepared get outputs;
593 : !> otherwise use the host AO/AO pool layout. One layout per context
594 : ! **************************************************************************************************
595 42706 : SUBROUTINE kpoint_operator_prepare(ctx, rsmat, slot, prepared, matrix_row, TARGET)
596 :
597 : TYPE(kpoint_operator_context_type), INTENT(INOUT) :: ctx
598 : TYPE(dbcsr_p_type), POINTER :: rsmat(:, :)
599 : INTEGER, INTENT(IN) :: slot
600 : LOGICAL, INTENT(OUT), OPTIONAL :: prepared
601 : INTEGER, INTENT(IN), OPTIONAL :: matrix_row
602 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: target
603 :
604 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_operator_prepare'
605 : INTEGER(KIND=int_8), PARAMETER :: max_bytes = 64_int_8*1024_int_8**2
606 :
607 : INTEGER :: compatible, handle, ispin, jslot, nrow, &
608 : row
609 : INTEGER(KIND=int_8) :: retained, storage
610 : LOGICAL :: ready
611 : TYPE(cp_fm_struct_type), POINTER :: fms
612 : TYPE(mp_para_env_type), POINTER :: para_env
613 :
614 42706 : CALL timeset(routineN, handle)
615 42706 : IF (PRESENT(prepared)) prepared = .FALSE.
616 42706 : CPASSERT(ASSOCIATED(ctx%kpoints))
617 42706 : CPASSERT(slot >= kpoint_slot_ks .AND. slot <= kpoint_slot_t)
618 42706 : CPASSERT(ASSOCIATED(ctx%kpoints%sab_nl))
619 42706 : IF (ASSOCIATED(ctx%work_sab)) THEN
620 6172 : CPASSERT(ASSOCIATED(ctx%work_sab, ctx%kpoints%sab_nl))
621 : END IF
622 42706 : IF (ALLOCATED(ctx%status)) THEN
623 0 : CPASSERT(ALL(ctx%status == kpoint_slot_idle))
624 : END IF
625 42706 : IF (.NOT. ALLOCATED(ctx%rs)) THEN
626 633484 : ALLOCATE (ctx%rs(ctx%nspin, kpoint_num_slots), ctx%rs_batch(ctx%nspin, kpoint_num_slots))
627 : END IF
628 88148 : DO ispin = 1, ctx%nspin
629 45442 : CALL kpoint_operator_rs_clear(ctx%rs(ispin, slot))
630 45442 : IF (ALLOCATED(ctx%rs_batch(ispin, slot)%data)) DEALLOCATE (ctx%rs_batch(ispin, slot)%data)
631 88148 : ctx%rs_batch(ispin, slot)%first = 0
632 : END DO
633 42706 : nrow = ctx%nspin
634 42706 : IF (PRESENT(matrix_row)) THEN
635 6634 : CPASSERT(matrix_row >= 1 .AND. matrix_row <= SIZE(rsmat, 1))
636 : nrow = 1
637 : ELSE
638 36072 : CPASSERT(SIZE(rsmat, 1) == ctx%nspin)
639 : END IF
640 42706 : para_env => ctx%kpoints%blacs_env_all%para_env
641 42706 : CALL kpoint_operator_ensure_group(ctx)
642 42706 : fms => ctx%fmlocal%matrix_struct
643 42706 : IF (PRESENT(TARGET)) fms => target%matrix_struct
644 42706 : CPASSERT(ASSOCIATED(fms))
645 : ! Routing uses the host group-rank numbering. Different BLACS layouts
646 : ! on that communicator are supported, foreign communicators are not.
647 42706 : compatible = MERGE(1, 0, ASSOCIATED(fms%para_env, ctx%kpoints%para_env_kp))
648 42706 : CALL para_env%min(compatible)
649 42706 : IF (compatible == 0) THEN
650 0 : CALL timestop(handle)
651 22 : RETURN
652 : END IF
653 42706 : IF (.NOT. ALLOCATED(ctx%rs_route)) THEN
654 36534 : ALLOCATE (ctx%rs_route)
655 36534 : CALL kpoint_operator_rs_plan_create(ctx%kpoints, rsmat, fms, max_bytes, ctx%rs_route, ready)
656 36534 : IF (.NOT. ready) THEN
657 22 : DEALLOCATE (ctx%rs_route)
658 22 : CALL timestop(handle)
659 22 : RETURN
660 : END IF
661 36512 : ctx%work_sab => ctx%kpoints%sab_nl
662 : END IF
663 42684 : compatible = MERGE(1, 0, cp_fm_struct_equivalent(fms, ctx%rs_route%matrix_struct))
664 87476 : DO ispin = 1, nrow
665 44792 : row = ispin
666 44792 : IF (PRESENT(matrix_row)) row = matrix_row
667 87476 : IF (.NOT. kpoint_operator_rs_compatible(ctx%rs_route, rsmat, row)) compatible = 0
668 : END DO
669 42684 : CALL para_env%min(compatible)
670 : ! Keep the single-k phase allowance in the snapshot budget for the uncached fallback.
671 42684 : retained = 0_int_8
672 170736 : DO jslot = kpoint_slot_ks, kpoint_slot_t
673 306996 : DO ispin = 1, ctx%nspin
674 136260 : IF (.NOT. ALLOCATED(ctx%rs(ispin, jslot)%blocks)) CYCLE
675 6914 : retained = retained + 128_int_8*SIZE(ctx%rs(ispin, jslot)%blocks, KIND=int_8)
676 134966 : IF (ASSOCIATED(ctx%rs(ispin, jslot)%values)) THEN
677 6914 : retained = retained + 8_int_8*SIZE(ctx%rs(ispin, jslot)%values, KIND=int_8)
678 : END IF
679 : END DO
680 : END DO
681 : storage = ctx%rs_route%storage + retained + INT(nrow, int_8)* &
682 : (MERGE(0_int_8, 8_int_8*ctx%rs_route%nvalue, para_env%num_pe == 1) + &
683 42684 : 128_int_8*SIZE(ctx%rs_route%entry, 2))
684 42684 : CALL para_env%max(storage)
685 42684 : ready = compatible == 1 .AND. storage <= max_bytes
686 42684 : IF (ready) THEN
687 87476 : DO ispin = 1, nrow
688 44792 : row = ispin
689 44792 : IF (PRESENT(matrix_row)) row = matrix_row
690 87476 : CALL kpoint_operator_rs_capture(rsmat, row, ctx%rs_route, ctx%kpoints, ctx%rs(ispin, slot))
691 : END DO
692 : END IF
693 42684 : IF (PRESENT(prepared)) prepared = ready
694 42684 : CALL timestop(handle)
695 :
696 : END SUBROUTINE kpoint_operator_prepare
697 :
698 : ! **************************************************************************************************
699 : !> \brief Release one captured row without invalidating other slots or routing.
700 : !> \param snapshot row to clear
701 : ! **************************************************************************************************
702 161302 : SUBROUTINE kpoint_operator_rs_clear(snapshot)
703 :
704 : TYPE(kpoint_operator_rs_type), INTENT(INOUT) :: snapshot
705 :
706 161302 : IF (ALLOCATED(snapshot%blocks)) DEALLOCATE (snapshot%blocks, snapshot%nonzero)
707 161302 : IF (ASSOCIATED(snapshot%values)) DEALLOCATE (snapshot%values)
708 161302 : NULLIFY (snapshot%source)
709 161302 : snapshot%matrix_row = 0
710 :
711 161302 : END SUBROUTINE kpoint_operator_rs_clear
712 :
713 : ! **************************************************************************************************
714 : !> \brief Match an explicitly prepared row, sharing spin-free rows across spins.
715 : !> \param ctx operator context
716 : !> \param source live array identity identifying the prepared operation
717 : !> \param ispin instance spin
718 : !> \param slot prepared slot
719 : !> \param matrix_row operator row
720 : !> \return snapshot spin index, or zero for the transfer backend
721 : ! **************************************************************************************************
722 33818 : FUNCTION kpoint_operator_rs_match(ctx, source, ispin, slot, matrix_row) RESULT(irs)
723 :
724 : TYPE(kpoint_operator_context_type), INTENT(IN) :: ctx
725 : TYPE(dbcsr_p_type), POINTER :: source(:, :)
726 : INTEGER, INTENT(IN) :: ispin, slot, matrix_row
727 : INTEGER :: irs
728 :
729 : INTEGER :: i
730 :
731 33818 : irs = 0
732 33818 : IF (.NOT. ALLOCATED(ctx%rs)) RETURN
733 33818 : i = ispin
734 33818 : IF (.NOT. ASSOCIATED(ctx%rs(i, slot)%source)) i = 1
735 33818 : IF (ASSOCIATED(ctx%rs(i, slot)%source, source) .AND. ctx%rs(i, slot)%matrix_row == matrix_row) irs = i
736 :
737 : END FUNCTION kpoint_operator_rs_match
738 :
739 : ! **************************************************************************************************
740 : !> \brief Cache whole-block routing to the FM tiles of one representative group.
741 : !> Equivalent ranks in the other groups receive only their tile's blocks.
742 : !> Boundary blocks may overlap several tiles. Missing blocks are packed
743 : !> as zero, permitting different operator sparsities on a fixed AO layout.
744 : !> \param kpoints initialized parallel environment and neighbor lists
745 : !> \param source operator layout template
746 : !> \param matrix_struct target group layout
747 : !> \param max_bytes per-rank budget for routing and capture buffers
748 : !> \param route context-owned structural route
749 : !> \param ready false collectively if unsupported or over budget
750 : ! **************************************************************************************************
751 36534 : SUBROUTINE kpoint_operator_rs_plan_create(kpoints, source, matrix_struct, max_bytes, route, ready)
752 :
753 : TYPE(kpoint_type), POINTER :: kpoints
754 : TYPE(dbcsr_p_type), POINTER :: source(:, :)
755 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
756 : INTEGER(KIND=int_8), INTENT(IN) :: max_bytes
757 : TYPE(kpoint_operator_rs_route_type), INTENT(OUT) :: route
758 : LOGICAL, INTENT(OUT) :: ready
759 :
760 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_operator_rs_plan_create'
761 :
762 : INTEGER :: compatible, dest, handle, i, iblock, &
763 : nblock, ncol, nrow, p, pos, rank
764 : INTEGER(KIND=int_8) :: nr(2), ns(2), storage
765 36534 : INTEGER(KIND=int_8), ALLOCATABLE :: recv_size(:, :), send_size(:, :)
766 36534 : INTEGER, ALLOCATABLE :: identities(:, :), owner(:), &
767 36534 : position(:, :), send_metadata(:, :)
768 36534 : INTEGER, POINTER :: col_size(:), row_size(:)
769 36534 : LOGICAL, ALLOCATABLE :: needed(:)
770 : TYPE(mp_para_env_type), POINTER :: inter, world
771 :
772 36534 : CALL timeset(routineN, handle)
773 36534 : ready = .FALSE.
774 36534 : world => kpoints%blacs_env_all%para_env
775 36534 : inter => kpoints%para_env_inter_kp
776 36534 : CALL dbcsr_get_info(source(1, 1)%matrix, row_blk_size=row_size, col_blk_size=col_size)
777 : CALL kp_transform_plan_create(route%traversal, kpoints%sab_nl, kpoints%cell_to_index, &
778 36534 : SIZE(source, 2), block_template=source(1, 1)%matrix)
779 36534 : compatible = MERGE(1, 0, route%traversal%symmetric)
780 352402 : IF (SUM(row_size) /= matrix_struct%nrow_global .OR. &
781 0 : SUM(col_size) /= matrix_struct%ncol_global) compatible = 0
782 36534 : IF (SIZE(row_size) == SIZE(col_size)) THEN
783 194468 : IF (ANY(row_size /= col_size)) compatible = 0
784 : ELSE
785 0 : compatible = 0
786 : END IF
787 36534 : CALL world%min(compatible)
788 36534 : IF (compatible == 0) THEN
789 0 : CALL timestop(handle)
790 0 : RETURN
791 : END IF
792 36534 : route%nimage = SIZE(source, 2)
793 182670 : ALLOCATE (route%row_size(SIZE(row_size)), route%col_size(SIZE(col_size)))
794 194468 : route%row_size(:) = row_size
795 194468 : route%col_size(:) = col_size
796 : ALLOCATE (identities(2, 0:world%num_pe - 1), owner(0:kpoints%para_env_kp%num_pe - 1), &
797 219204 : needed(0:kpoints%para_env_kp%num_pe - 1))
798 109602 : CALL world%allgather([kpoints%para_env_kp%mepos, inter%mepos], identities)
799 84198 : owner(:) = -1
800 109602 : DO rank = 0, world%num_pe - 1
801 109602 : IF (identities(2, rank) == 0) owner(identities(1, rank)) = rank
802 : END DO
803 84198 : CPASSERT(ALL(owner >= 0))
804 146136 : ALLOCATE (send_size(2, 0:world%num_pe - 1), recv_size(2, 0:world%num_pe - 1))
805 36534 : send_size(:, :) = 0_int_8
806 8454946 : DO i = 1, route%traversal%nentry
807 8418412 : nrow = row_size(route%traversal%row(i))
808 8418412 : ncol = col_size(route%traversal%col(i))
809 8418412 : CALL kpoint_operator_rs_owners(route, matrix_struct, i, needed)
810 20821370 : DO p = 0, SIZE(owner) - 1
811 12366424 : IF (.NOT. needed(p)) CYCLE
812 10253764 : dest = owner(p)
813 10253764 : send_size(1, dest) = send_size(1, dest) + 1_int_8
814 20784836 : send_size(2, dest) = send_size(2, dest) + INT(nrow, int_8)*ncol
815 : END DO
816 : END DO
817 36534 : CALL world%alltoall(send_size, recv_size, 2)
818 : ! Only inter-rank zero receives from the world. The other groups hold
819 : ! equivalent FM tiles, as required by kpoint_env_initialize.
820 255738 : ns = SUM(send_size, DIM=2)
821 255738 : nr = SUM(recv_size, DIM=2)
822 36534 : CALL inter%bcast(nr, 0)
823 : ! Conservative peak: traversal, retained route, temporary metadata and
824 : ! one snapshot. The full batch budget is checked again by prepare.
825 : storage = 40_int_8*route%traversal%nentry + &
826 : 4_int_8*(SIZE(row_size) + SIZE(col_size) + matrix_struct%nrow_global + matrix_struct%ncol_global) + &
827 : 48_int_8*ns(1) + MERGE(0_int_8, 8_int_8*ns(2), world%num_pe == 1) + 68_int_8*nr(1) + &
828 36534 : 128_int_8*world%num_pe + 8_int_8*SIZE(owner) + 8_int_8
829 36534 : route%storage = storage
830 73068 : storage = storage + MERGE(0_int_8, 8_int_8*nr(2), world%num_pe == 1) + 128_int_8*nr(1)
831 : ! Borrowed values do not consume the snapshot budget; still bound integer route offsets.
832 182670 : IF (MAXVAL(ns) >= INT(HUGE(0), int_8) .OR. MAXVAL(nr) >= INT(HUGE(0), int_8)) THEN
833 0 : storage = MAX(storage, max_bytes + 1_int_8)
834 : END IF
835 36534 : CALL world%max(storage)
836 36534 : IF (storage > max_bytes) THEN
837 22 : CALL timestop(handle)
838 22 : RETURN
839 : END IF
840 : ! The byte bound also keeps all counts/displacements within default INTEGER.
841 36512 : nblock = INT(nr(1))
842 36512 : route%nvalue = INT(nr(2))
843 0 : ALLOCATE (route%scount(2, 0:world%num_pe - 1), route%rcount(2, 0:world%num_pe - 1), &
844 219072 : route%sdisp(2, 0:world%num_pe - 1), route%rdisp(2, 0:world%num_pe - 1))
845 255584 : route%scount(:, :) = INT(send_size)
846 255584 : route%rcount(:, :) = INT(recv_size)
847 36512 : IF (SIZE(owner) == 1) CALL inter%bcast(route%rcount, 0)
848 109536 : route%sdisp(:, 0) = 0
849 109536 : route%rdisp(:, 0) = 0
850 73024 : DO rank = 1, world%num_pe - 1
851 109536 : route%sdisp(:, rank) = route%sdisp(:, rank - 1) + route%scount(:, rank - 1)
852 146048 : route%rdisp(:, rank) = route%rdisp(:, rank - 1) + route%rcount(:, rank - 1)
853 : END DO
854 0 : ALLOCATE (route%send_entry(INT(ns(1))), route%send_first(INT(ns(1)) + 1), &
855 0 : route%send_values(MERGE(0, INT(ns(2)), world%num_pe == 1)), send_metadata(8, INT(ns(1))), &
856 0 : position(2, 0:world%num_pe - 1), route%entry(8, nblock), &
857 582260 : route%first(nblock + 1))
858 255584 : position(:, :) = route%sdisp
859 8334804 : DO i = 1, route%traversal%nentry
860 8298292 : nrow = row_size(route%traversal%row(i))
861 8298292 : ncol = col_size(route%traversal%col(i))
862 8298292 : CALL kpoint_operator_rs_owners(route, matrix_struct, i, needed)
863 20460988 : DO p = 0, SIZE(owner) - 1
864 12126184 : IF (.NOT. needed(p)) CYCLE
865 10057029 : dest = owner(p)
866 10057029 : iblock = position(1, dest) + 1
867 10057029 : pos = position(2, dest) + 1
868 10057029 : route%send_entry(iblock) = i
869 10057029 : route%send_first(iblock) = pos
870 50285145 : send_metadata(1:4, iblock) = [route%traversal%row_offset(i), route%traversal%col_offset(i), nrow, ncol]
871 40228116 : send_metadata(5:7, iblock) = NINT(route%traversal%symmetry_sign(i))*route%traversal%cell(:, i)
872 10057029 : send_metadata(8, iblock) = MERGE(1, 0, route%traversal%row(i) /= route%traversal%col(i))
873 40538534 : position(:, dest) = position(:, dest) + [1, nrow*ncol]
874 : END DO
875 : END DO
876 36512 : route%send_first(INT(ns(1)) + 1) = INT(ns(2)) + 1
877 36512 : IF (SIZE(owner) == 1) THEN
878 127020 : CALL world%allgatherv(send_metadata, route%entry, 8*route%rcount(1, :), 8*route%rdisp(1, :))
879 : ELSE
880 : CALL world%alltoall(send_metadata, 8*route%scount(1, :), 8*route%sdisp(1, :), &
881 99972 : route%entry, 8*route%rcount(1, :), 8*route%rdisp(1, :))
882 11108 : CALL inter%bcast(route%entry, 0)
883 : END IF
884 36512 : route%first(1) = 1
885 14563941 : DO i = 1, nblock
886 14563941 : route%first(i + 1) = route%first(i) + route%entry(3, i)*route%entry(4, i)
887 : END DO
888 36512 : CPASSERT(route%first(nblock + 1) == route%nvalue + 1)
889 36512 : route%matrix_struct => matrix_struct
890 36512 : CALL cp_fm_struct_retain(route%matrix_struct)
891 36512 : CALL kpoint_operator_rs_tile_map(route)
892 36512 : ready = .TRUE.
893 36512 : CALL timestop(handle)
894 :
895 109602 : END SUBROUTINE kpoint_operator_rs_plan_create
896 :
897 : ! **************************************************************************************************
898 : !> \brief Cache local intervals intersecting each routed block and its Hermitian mirror.
899 : !> FM indices are ordered globally; dummy indices on empty tiles lie outside the matrix.
900 : !> \param route fixed routing and target distribution
901 : ! **************************************************************************************************
902 36512 : SUBROUTINE kpoint_operator_rs_tile_map(route)
903 : TYPE(kpoint_operator_rs_route_type), INTENT(INOUT) :: route
904 :
905 73024 : INTEGER :: c0, i, nc, nr, r0, rows(0:route%matrix_struct%nrow_global), &
906 36512 : cols(0:route%matrix_struct%ncol_global)
907 :
908 951378 : rows(:) = 0
909 951378 : cols(:) = 0
910 804533 : DO i = 1, SIZE(route%matrix_struct%row_indices)
911 768021 : r0 = route%matrix_struct%row_indices(i)
912 804533 : IF (r0 >= 1 .AND. r0 <= UBOUND(rows, 1)) rows(r0) = 1
913 : END DO
914 914866 : DO i = 1, SIZE(route%matrix_struct%col_indices)
915 878354 : c0 = route%matrix_struct%col_indices(i)
916 914866 : IF (c0 >= 1 .AND. c0 <= UBOUND(cols, 1)) cols(c0) = 1
917 : END DO
918 914866 : DO i = 1, UBOUND(rows, 1)
919 914866 : rows(i) = rows(i) + rows(i - 1)
920 : END DO
921 914866 : DO i = 1, UBOUND(cols, 1)
922 914866 : cols(i) = cols(i) + cols(i - 1)
923 : END DO
924 36512 : route%nrow = rows(UBOUND(rows, 1))
925 36512 : route%ncol = cols(UBOUND(cols, 1))
926 109536 : ALLOCATE (route%tile(4, 2, SIZE(route%entry, 2)))
927 14563941 : DO i = 1, SIZE(route%entry, 2)
928 14527429 : r0 = route%entry(1, i)
929 14527429 : c0 = route%entry(2, i)
930 14527429 : nr = route%entry(3, i)
931 14527429 : nc = route%entry(4, i)
932 72637145 : route%tile(:, 1, i) = [rows(r0 - 1) + 1, rows(r0 + nr - 1), cols(c0 - 1) + 1, cols(c0 + nc - 1)]
933 72673657 : route%tile(:, 2, i) = [rows(c0 - 1) + 1, rows(c0 + nc - 1), cols(r0 - 1) + 1, cols(r0 + nr - 1)]
934 : END DO
935 36512 : END SUBROUTINE kpoint_operator_rs_tile_map
936 :
937 : ! **************************************************************************************************
938 : !> \brief Find group-local FM owners intersecting a stored block or its mirror.
939 : !> \param route source AO block layout and neighbor traversal
940 : !> \param fms target FM layout within each group
941 : !> \param i traversal entry
942 : !> \param needed recipient group-local ranks, without duplicate mirror sends
943 : ! **************************************************************************************************
944 16716704 : SUBROUTINE kpoint_operator_rs_owners(route, fms, i, needed)
945 :
946 : TYPE(kpoint_operator_rs_route_type), INTENT(IN) :: route
947 : TYPE(cp_fm_struct_type), INTENT(IN) :: fms
948 : INTEGER, INTENT(IN) :: i
949 : LOGICAL, INTENT(OUT) :: needed(0:)
950 :
951 : INTEGER :: c0, j, nc, nr, pc, pr, r0
952 33433408 : LOGICAL :: rows(2, 0:fms%context%num_pe(1) - 1), cols(2, 0:fms%context%num_pe(2) - 1)
953 :
954 16716704 : r0 = route%traversal%row_offset(i)
955 16716704 : c0 = route%traversal%col_offset(i)
956 16716704 : nr = route%row_size(route%traversal%row(i))
957 16716704 : nc = route%col_size(route%traversal%col(i))
958 90194528 : rows(:, :) = .FALSE.
959 66866816 : cols(:, :) = .FALSE.
960 123687914 : DO j = r0, r0 + nr - 1
961 106971210 : rows(1, fms%g2p_row(j)) = .TRUE.
962 123687914 : cols(2, fms%g2p_col(j)) = .TRUE.
963 : END DO
964 122127850 : DO j = c0, c0 + nc - 1
965 105411146 : rows(2, fms%g2p_row(j)) = .TRUE.
966 122127850 : cols(1, fms%g2p_col(j)) = .TRUE.
967 : END DO
968 27641048 : IF (route%traversal%row(i) == route%traversal%col(i)) rows(2, :) = .FALSE.
969 41209312 : needed(:) = .FALSE.
970 33433408 : DO pc = 0, fms%context%num_pe(2) - 1
971 57926016 : DO pr = 0, fms%context%num_pe(1) - 1
972 : needed(fms%context%blacs2mpi(pr, pc)) = &
973 45391127 : (rows(1, pr) .AND. cols(1, pc)) .OR. (rows(2, pr) .AND. cols(2, pc))
974 : END DO
975 : END DO
976 :
977 16716704 : END SUBROUTINE kpoint_operator_rs_owners
978 :
979 : ! **************************************************************************************************
980 : !> \brief Check the source AO layout, not its values or sparse block presence.
981 : !> \param route fixed context route
982 : !> \param source current operator images
983 : !> \param matrix_row row to capture
984 : !> \return whether the row has the route's image count and AO layout
985 : ! **************************************************************************************************
986 44792 : FUNCTION kpoint_operator_rs_compatible(route, source, matrix_row) RESULT(compatible)
987 :
988 : TYPE(kpoint_operator_rs_route_type), INTENT(IN) :: route
989 : TYPE(dbcsr_p_type), POINTER :: source(:, :)
990 : INTEGER, INTENT(IN) :: matrix_row
991 : LOGICAL :: compatible
992 :
993 : INTEGER :: image
994 44792 : INTEGER, POINTER :: col_size(:), row_size(:)
995 :
996 44792 : compatible = .FALSE.
997 44792 : IF (SIZE(source, 2) /= route%nimage) RETURN
998 44792 : IF (matrix_row < 1 .OR. matrix_row > SIZE(source, 1)) RETURN
999 2552982 : DO image = 1, SIZE(source, 2)
1000 2508190 : CALL dbcsr_get_info(source(matrix_row, image)%matrix, row_blk_size=row_size, col_blk_size=col_size)
1001 2508190 : IF (SIZE(row_size) /= SIZE(route%row_size) .OR. SIZE(col_size) /= SIZE(route%col_size)) RETURN
1002 26482606 : IF (ANY(row_size /= route%row_size) .OR. ANY(col_size /= route%col_size)) RETURN
1003 : END DO
1004 44792 : compatible = .TRUE.
1005 :
1006 44792 : END FUNCTION kpoint_operator_rs_compatible
1007 :
1008 : ! **************************************************************************************************
1009 : !> \brief Capture current values using cached routing, with no per-k communication.
1010 : !> \param source current real-space operator images
1011 : !> \param matrix_row row to snapshot
1012 : !> \param route context-owned metadata, counts and reusable packing buffer
1013 : !> \param kpoints initialized host parallel environment
1014 : !> \param snapshot immutable tile-local values
1015 : ! **************************************************************************************************
1016 44792 : SUBROUTINE kpoint_operator_rs_capture(source, matrix_row, route, kpoints, snapshot)
1017 :
1018 : TYPE(dbcsr_p_type), POINTER :: source(:, :)
1019 : INTEGER, INTENT(IN) :: matrix_row
1020 : TYPE(kpoint_operator_rs_route_type), INTENT(INOUT) :: route
1021 : TYPE(kpoint_type), POINTER :: kpoints
1022 : TYPE(kpoint_operator_rs_type), INTENT(OUT) :: snapshot
1023 :
1024 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_operator_rs_capture'
1025 :
1026 : INTEGER :: col, handle, i, iblock, ncol, nrow, pos, &
1027 : row
1028 : LOGICAL :: found
1029 44792 : REAL(KIND=dp), POINTER :: block(:, :)
1030 : TYPE(mp_para_env_type), POINTER :: inter, world
1031 :
1032 44792 : CALL timeset(routineN, handle)
1033 44792 : world => kpoints%blacs_env_all%para_env
1034 44792 : inter => kpoints%para_env_inter_kp
1035 19041691 : ALLOCATE (snapshot%blocks(SIZE(route%entry, 2)), snapshot%nonzero(SIZE(route%entry, 2)))
1036 44792 : snapshot%source => source
1037 44792 : snapshot%matrix_row = matrix_row
1038 44792 : IF (world%num_pe == 1) THEN
1039 : ! No value copy on a single world rank. Lookups happen before workers start.
1040 0 : DO iblock = 1, SIZE(route%send_entry)
1041 0 : i = route%send_entry(iblock)
1042 : CALL dbcsr_get_readonly_block_p(source(matrix_row, route%traversal%image(i))%matrix, &
1043 0 : route%traversal%row(i), route%traversal%col(i), block, found)
1044 0 : snapshot%nonzero(iblock) = found
1045 0 : IF (found) THEN
1046 0 : CPASSERT(SIZE(block, 1) == route%entry(3, iblock))
1047 0 : CPASSERT(SIZE(block, 2) == route%entry(4, iblock))
1048 0 : snapshot%blocks(iblock)%data => block
1049 : END IF
1050 : END DO
1051 0 : CALL timestop(handle)
1052 0 : RETURN
1053 : END IF
1054 134376 : ALLOCATE (snapshot%values(route%nvalue))
1055 748700624 : route%send_values(:) = 0.0_dp
1056 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(source, matrix_row, route) &
1057 44792 : !$OMP PRIVATE(iblock, i, block, found, nrow, ncol, pos, row, col) SCHEDULE(STATIC)
1058 : DO iblock = 1, SIZE(route%send_entry)
1059 : i = route%send_entry(iblock)
1060 : CALL dbcsr_get_readonly_block_p(source(matrix_row, route%traversal%image(i))%matrix, &
1061 : route%traversal%row(i), route%traversal%col(i), block, found)
1062 : IF (.NOT. found) CYCLE
1063 : nrow = route%row_size(route%traversal%row(i))
1064 : ncol = route%col_size(route%traversal%col(i))
1065 : CPASSERT(SIZE(block, 1) == nrow .AND. SIZE(block, 2) == ncol)
1066 : pos = route%send_first(iblock)
1067 : DO col = 1, ncol
1068 : DO row = 1, nrow
1069 : route%send_values(pos) = block(row, col)
1070 : pos = pos + 1
1071 : END DO
1072 : END DO
1073 : END DO
1074 : !$OMP END PARALLEL DO
1075 44792 : IF (kpoints%para_env_kp%num_pe == 1) THEN
1076 150360 : CALL world%allgatherv(route%send_values, snapshot%values, route%rcount(2, :), route%rdisp(2, :))
1077 : ELSE
1078 : CALL world%alltoall(route%send_values, route%scount(2, :), route%sdisp(2, :), &
1079 132480 : snapshot%values, route%rcount(2, :), route%rdisp(2, :))
1080 14720 : CALL inter%bcast(snapshot%values, 0)
1081 : END IF
1082 18862523 : DO iblock = 1, SIZE(route%entry, 2)
1083 96583655 : snapshot%nonzero(iblock) = ANY(snapshot%values(route%first(iblock):route%first(iblock + 1) - 1) /= 0.0_dp)
1084 : snapshot%blocks(iblock)%data(1:route%entry(3, iblock), 1:route%entry(4, iblock)) => &
1085 18862523 : snapshot%values(route%first(iblock):route%first(iblock + 1) - 1)
1086 : END DO
1087 44792 : snapshot%source => source
1088 44792 : snapshot%matrix_row = matrix_row
1089 44792 : CALL timestop(handle)
1090 :
1091 89584 : END SUBROUTINE kpoint_operator_rs_capture
1092 :
1093 : ! **************************************************************************************************
1094 : !> \brief Deliver one prepared operator, building a bounded window of local k points on a miss.
1095 : !> The callers retain their round/spin/slot ordering. Each row has its own
1096 : !> window so interleaved spins and H/S/T do not evict one another.
1097 : !> Only prepared values are read ahead; transfer and grid paths are unchanged.
1098 : !> \param ctx serialized operator context
1099 : !> \param irs captured row, possibly shared by spin-free instances
1100 : !> \param slot prepared operator slot
1101 : !> \param ik global k-point index owned by this group
1102 : !> \param TARGET optional complex output on the prepared distribution
1103 : !> \param fm_re optional real output instead of target; fm_im may be omitted for a real operator
1104 : !> \param fm_im optional imaginary output on the same distribution
1105 : ! **************************************************************************************************
1106 33776 : SUBROUTINE kpoint_operator_rs_build(ctx, irs, slot, ik, TARGET, fm_re, fm_im)
1107 :
1108 : TYPE(kpoint_operator_context_type), INTENT(INOUT) :: ctx
1109 : INTEGER, INTENT(IN) :: irs, slot, ik
1110 : TYPE(cp_cfm_type), INTENT(INOUT), OPTIONAL :: target
1111 : TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL :: fm_re, fm_im
1112 :
1113 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_operator_rs_build'
1114 :
1115 : INTEGER :: first, handle, k, last, ncol, nrow, width
1116 : INTEGER(KIND=int_8) :: bytes_per_k
1117 : LOGICAL :: layout_ok
1118 33776 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1119 :
1120 33776 : CALL timeset(routineN, handle)
1121 33776 : IF (PRESENT(TARGET)) THEN
1122 32096 : layout_ok = cp_fm_struct_equivalent(target%matrix_struct, ctx%rs_route%matrix_struct)
1123 : ELSE
1124 1680 : CPASSERT(PRESENT(fm_re))
1125 1680 : layout_ok = cp_fm_struct_equivalent(fm_re%matrix_struct, ctx%rs_route%matrix_struct)
1126 1680 : IF (PRESENT(fm_im)) THEN
1127 1384 : layout_ok = layout_ok .AND. cp_fm_struct_equivalent(fm_im%matrix_struct, ctx%rs_route%matrix_struct)
1128 : END IF
1129 : END IF
1130 32392 : CPASSERT(layout_ok)
1131 : ! Avoid nested component array sections in GCC 16 bounds checks.
1132 33776 : xkp => ctx%kpoints%xkp
1133 33776 : first = ctx%kpoints%kp_range(1)
1134 33776 : last = ctx%kpoints%kp_range(2)
1135 33776 : CPASSERT(ik >= first .AND. ik <= last)
1136 33776 : nrow = ctx%rs_route%nrow
1137 33776 : ncol = ctx%rs_route%ncol
1138 33776 : bytes_per_k = 16_int_8*(INT(nrow, int_8)*ncol + SIZE(ctx%rs_route%entry, 2, KIND=int_8))
1139 : width = MIN(rs_batch_size, last - first + 1, &
1140 33776 : INT(rs_batch_bytes/(INT(ctx%nspin*kpoint_num_slots, int_8)*MAX(1_int_8, bytes_per_k))))
1141 33776 : IF (width < 2 .OR. nrow == 0 .OR. ncol == 0) THEN
1142 : ! No extra dense tile for a large/empty layout or a single local k point.
1143 : CALL kpoint_operator_rs_build_batch(ctx%rs(irs, slot), ctx%rs_route, xkp(:, ik:ik), &
1144 10882 : TARGET=TARGET, fm_re=fm_re, fm_im=fm_im)
1145 : ELSE
1146 : ASSOCIATE (work => ctx%rs_batch(irs, slot))
1147 25750 : IF (.NOT. ALLOCATED(work%data)) ALLOCATE (work%data(width, nrow, ncol))
1148 22894 : IF (work%first == 0 .OR. ik < work%first .OR. ik >= work%first + width) THEN
1149 6568 : work%first = first + ((ik - first)/width)*width
1150 : CALL kpoint_operator_rs_build_batch(ctx%rs(irs, slot), ctx%rs_route, &
1151 6568 : xkp(:, work%first:MIN(work%first + width - 1, last)), batch=work%data)
1152 : END IF
1153 22894 : k = ik - work%first + 1
1154 : ! Deliver a work copy: eigensolvers may overwrite H/S without corrupting the window.
1155 45788 : IF (PRESENT(TARGET)) THEN
1156 29380020 : target%local_data(1:nrow, 1:ncol) = work%data(k, :, :)
1157 707710 : target%local_data(nrow + 1:, :) = z_zero
1158 21950 : target%local_data(:, ncol + 1:) = z_zero
1159 : ELSE
1160 944 : CPASSERT(PRESENT(fm_re))
1161 103934 : fm_re%local_data(1:nrow, 1:ncol) = REAL(work%data(k, :, :), dp)
1162 103934 : IF (PRESENT(fm_im)) fm_im%local_data(1:nrow, 1:ncol) = AIMAG(work%data(k, :, :))
1163 9568 : fm_re%local_data(nrow + 1:, :) = 0.0_dp
1164 944 : fm_re%local_data(:, ncol + 1:) = 0.0_dp
1165 9568 : IF (PRESENT(fm_im)) fm_im%local_data(nrow + 1:, :) = 0.0_dp
1166 944 : IF (PRESENT(fm_im)) fm_im%local_data(:, ncol + 1:) = 0.0_dp
1167 : END IF
1168 : END ASSOCIATE
1169 : END IF
1170 33776 : IF (ctx%kpoints%use_real_wfn) THEN
1171 296 : IF (PRESENT(TARGET)) target%local_data(:, :) = CMPLX(REAL(target%local_data, dp), 0.0_dp, dp)
1172 296 : IF (PRESENT(fm_im)) fm_im%local_data(:, :) = 0.0_dp
1173 : END IF
1174 33776 : CALL timestop(handle)
1175 :
1176 33776 : END SUBROUTINE kpoint_operator_rs_build
1177 :
1178 : ! **************************************************************************************************
1179 : !> \brief Accumulate several k points in one route/value traversal, without MPI or BLACS.
1180 : !> The leading batch dimension makes updates to different k points contiguous.
1181 : !> Columns remain independent OpenMP work items; each target retains the
1182 : !> original normal/mirrored contribution order. Phase scratch is call-local.
1183 : !> The singleton fallback writes directly to the supplied CFM or FM pair.
1184 : !> \param snapshot immutable real-space contributions
1185 : !> \param route structural block routing and target tile map
1186 : !> \param xkp actual coordinates, one column per k point
1187 : !> \param batch optional packed complex output (k point, local row, local column)
1188 : !> \param TARGET optional single complex output instead of batch
1189 : !> \param fm_re optional single real output instead of batch/target
1190 : !> \param fm_im optional single imaginary output
1191 : !> \param phase_work optional private phase storage supplied by a worker
1192 : !> \param parallel_columns split columns across threads, disabled inside k-point workers
1193 : ! **************************************************************************************************
1194 70243 : SUBROUTINE kpoint_operator_rs_build_batch(snapshot, route, xkp, batch, TARGET, fm_re, fm_im, phase_work, parallel_columns)
1195 :
1196 : TYPE(kpoint_operator_rs_type), INTENT(IN) :: snapshot
1197 : TYPE(kpoint_operator_rs_route_type), INTENT(IN) :: route
1198 : REAL(KIND=dp), INTENT(IN) :: xkp(:, :)
1199 : COMPLEX(KIND=dp), CONTIGUOUS, INTENT(INOUT), &
1200 : OPTIONAL, TARGET :: batch(:, :, :)
1201 : TYPE(cp_cfm_type), INTENT(INOUT), OPTIONAL :: target
1202 : TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL :: fm_re, fm_im
1203 : COMPLEX(KIND=dp), INTENT(INOUT), OPTIONAL, TARGET :: phase_work(:, :)
1204 : LOGICAL, INTENT(IN), OPTIONAL :: parallel_columns
1205 :
1206 70243 : COMPLEX(KIND=dp), ALLOCATABLE, TARGET :: owned_phase(:, :)
1207 70243 : COMPLEX(KIND=dp), CONTIGUOUS, POINTER :: output(:, :, :)
1208 70243 : COMPLEX(KIND=dp), POINTER :: phase(:, :)
1209 : INTEGER :: first, iblock, k, last, nbatch, ncol, &
1210 : nthreads, thread
1211 : LOGICAL :: parallel
1212 : REAL(KIND=dp) :: arg
1213 :
1214 70243 : nbatch = SIZE(xkp, 2)
1215 70243 : CPASSERT(nbatch >= 1 .AND. nbatch <= rs_batch_size .AND. SIZE(xkp, 1) == 3)
1216 280972 : CPASSERT(COUNT([PRESENT(batch), PRESENT(TARGET), PRESENT(fm_re)]) == 1)
1217 70243 : NULLIFY (output)
1218 70243 : IF (PRESENT(batch)) THEN
1219 6568 : CPASSERT(SIZE(batch, 1) >= nbatch)
1220 6568 : output => batch
1221 : ELSE
1222 63675 : CPASSERT(nbatch == 1)
1223 63675 : IF (PRESENT(TARGET)) THEN
1224 59833 : output(1:1, 1:SIZE(target%local_data, 1), 1:SIZE(target%local_data, 2)) => target%local_data
1225 : ELSE
1226 3842 : CPASSERT(PRESENT(fm_re))
1227 10829256 : fm_re%local_data(:, :) = 0.0_dp
1228 9732754 : IF (PRESENT(fm_im)) fm_im%local_data(:, :) = 0.0_dp
1229 : END IF
1230 : END IF
1231 70243 : IF (ASSOCIATED(output)) output(:, :, :) = z_zero
1232 70243 : IF (route%nrow == 0 .OR. route%ncol == 0) RETURN
1233 70243 : IF (PRESENT(phase_work)) THEN
1234 52793 : CPASSERT(SIZE(phase_work, 1) >= nbatch)
1235 52793 : CPASSERT(SIZE(phase_work, 2) == SIZE(route%entry, 2))
1236 52793 : phase => phase_work(1:nbatch, :)
1237 : ELSE
1238 69800 : ALLOCATE (owned_phase(nbatch, SIZE(route%entry, 2)))
1239 17450 : phase => owned_phase
1240 : END IF
1241 36285278 : DO iblock = 1, SIZE(route%entry, 2)
1242 36215035 : IF (.NOT. snapshot%nonzero(iblock)) CYCLE
1243 86924269 : DO k = 1, nbatch
1244 203219524 : arg = twopi*DOT_PRODUCT(REAL(route%entry(5:7, iblock), dp), xkp(:, k))
1245 87019916 : phase(k, iblock) = CMPLX(COS(arg), SIN(arg), KIND=dp)
1246 : END DO
1247 : END DO
1248 70243 : ncol = route%ncol
1249 70243 : parallel = .TRUE.
1250 70243 : IF (PRESENT(parallel_columns)) parallel = parallel_columns
1251 70243 : nthreads = 1
1252 70243 : !$ IF (parallel) nthreads = MIN(omp_get_max_threads(), ncol)
1253 17450 : IF (nthreads > 1) THEN
1254 : !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads) SHARED(snapshot,route,phase,output,fm_re,fm_im,ncol) &
1255 : !$OMP PRIVATE(first,last,thread,nthreads)
1256 : thread = 0
1257 : nthreads = 1
1258 : !$ thread = omp_get_thread_num()
1259 : !$ nthreads = omp_get_num_threads()
1260 : first = ncol*thread/nthreads + 1
1261 : last = ncol*(thread + 1)/nthreads
1262 0 : CALL kpoint_operator_rs_accumulate(snapshot, route, phase, first, last, output, fm_re, fm_im)
1263 : !$OMP END PARALLEL
1264 : ELSE
1265 : ! Called inside a k-point worker: no nested parallel region, timer or allocation.
1266 70243 : CALL kpoint_operator_rs_accumulate(snapshot, route, phase, 1, ncol, output, fm_re, fm_im)
1267 : END IF
1268 :
1269 140486 : END SUBROUTINE kpoint_operator_rs_build_batch
1270 :
1271 : ! **************************************************************************************************
1272 : !> \brief Shared borrowed/packed phase-sum kernel; each caller owns a disjoint column interval.
1273 : !> \param snapshot read-only block views
1274 : !> \param route structural layout
1275 : !> \param phase precomputed phases (k, block)
1276 : !> \param first first local column owned by this call
1277 : !> \param last last local column owned by this call
1278 : !> \param output packed complex output, or disassociated for the FM pair
1279 : !> \param fm_re real output
1280 : !> \param fm_im imaginary output
1281 : ! **************************************************************************************************
1282 70243 : SUBROUTINE kpoint_operator_rs_accumulate(snapshot, route, phase, first, last, output, fm_re, fm_im)
1283 : TYPE(kpoint_operator_rs_type), INTENT(IN) :: snapshot
1284 : TYPE(kpoint_operator_rs_route_type), INTENT(IN) :: route
1285 : COMPLEX(KIND=dp), INTENT(IN) :: phase(:, :)
1286 : INTEGER, INTENT(IN) :: first, last
1287 : COMPLEX(KIND=dp), CONTIGUOUS, INTENT(IN), POINTER :: output(:, :, :)
1288 : TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL :: fm_re, fm_im
1289 :
1290 : COMPLEX(KIND=dp) :: factors(rs_batch_size), value
1291 70243 : COMPLEX(KIND=dp), CONTIGUOUS, POINTER :: full(:, :)
1292 : INTEGER :: c0, col, hi, iblock, j, k, lo, nc, nr, &
1293 : part, r0, row, source_col, source_row
1294 : REAL(KIND=dp) :: element
1295 :
1296 0 : IF (first > last) RETURN
1297 70243 : NULLIFY (full)
1298 70243 : IF (ASSOCIATED(output)) THEN
1299 : IF (SIZE(output, 1) == 1 .AND. &
1300 66401 : route%nrow == route%matrix_struct%nrow_global .AND. &
1301 : route%ncol == route%matrix_struct%ncol_global) THEN
1302 49831 : full(1:SIZE(output, 2), 1:SIZE(output, 3)) => output
1303 : END IF
1304 : END IF
1305 36285278 : DO iblock = 1, SIZE(route%entry, 2)
1306 36215035 : IF (.NOT. snapshot%nonzero(iblock)) CYCLE
1307 36049145 : r0 = route%entry(1, iblock)
1308 36049145 : c0 = route%entry(2, iblock)
1309 36049145 : nr = route%entry(3, iblock)
1310 36049145 : nc = route%entry(4, iblock)
1311 36049145 : IF (ASSOCIATED(full)) THEN
1312 : ! Update whole contiguous blocks, clipped to this worker's columns.
1313 22503928 : lo = MAX(first, c0)
1314 22503928 : hi = MIN(last, c0 + nc - 1)
1315 22503928 : IF (lo <= hi) THEN
1316 : full(r0:r0 + nr - 1, lo:hi) = full(r0:r0 + nr - 1, lo:hi) + &
1317 1391994932 : phase(1, iblock)*snapshot%blocks(iblock)%data(:, lo - c0 + 1:hi - c0 + 1)
1318 : END IF
1319 22503928 : IF (route%entry(8, iblock) /= 0) THEN
1320 15634493 : lo = MAX(first, r0)
1321 15634493 : hi = MIN(last, r0 + nr - 1)
1322 15634493 : IF (lo <= hi) THEN
1323 : full(c0:c0 + nc - 1, lo:hi) = full(c0:c0 + nc - 1, lo:hi) + &
1324 : CONJG(phase(1, iblock))* &
1325 905700813 : TRANSPOSE(snapshot%blocks(iblock)%data(lo - r0 + 1:hi - r0 + 1, :))
1326 : END IF
1327 : END IF
1328 : CYCLE
1329 : END IF
1330 40705894 : DO part = 1, 2
1331 27090434 : IF (part == 2 .AND. route%entry(8, iblock) == 0) CYCLE
1332 23912072 : lo = MAX(first, route%tile(3, part, iblock))
1333 23912072 : hi = MIN(last, route%tile(4, part, iblock))
1334 23912072 : IF (lo > hi) CYCLE
1335 73897462 : factors(1:SIZE(phase, 1)) = phase(:, iblock)
1336 45596509 : IF (part == 2) factors(1:SIZE(phase, 1)) = CONJG(factors(1:SIZE(phase, 1)))
1337 219496369 : DO j = lo, hi
1338 181873190 : col = route%matrix_struct%col_indices(j)
1339 1520376409 : DO row = route%tile(1, part, iblock), route%tile(2, part, iblock)
1340 1311412785 : source_row = route%matrix_struct%row_indices(row)
1341 1311412785 : source_col = col
1342 1311412785 : IF (part == 2) THEN
1343 545020995 : source_col = source_row
1344 545020995 : source_row = col
1345 : END IF
1346 1311412785 : element = snapshot%blocks(iblock)%data(source_row - r0 + 1, source_col - c0 + 1)
1347 1493285975 : IF (ASSOCIATED(output)) THEN
1348 : !$OMP SIMD
1349 : DO k = 1, SIZE(phase, 1)
1350 1925489731 : output(k, row, j) = output(k, row, j) + factors(k)*element
1351 : END DO
1352 : ELSE
1353 716397876 : value = factors(1)*element
1354 716397876 : fm_re%local_data(row, j) = fm_re%local_data(row, j) + REAL(value, dp)
1355 716397876 : IF (PRESENT(fm_im)) fm_im%local_data(row, j) = fm_im%local_data(row, j) + AIMAG(value)
1356 : END IF
1357 : END DO
1358 : END DO
1359 : END DO
1360 : END DO
1361 70243 : END SUBROUTINE kpoint_operator_rs_accumulate
1362 :
1363 : ! **************************************************************************************************
1364 : !> \brief Size per-worker phase scratch for direct assembly into a solver matrix.
1365 : !> \param ctx prepared context
1366 : !> \param nentry number of route entries
1367 : ! **************************************************************************************************
1368 25820 : SUBROUTINE kpoint_operator_work_size(ctx, nentry)
1369 : TYPE(kpoint_operator_context_type), INTENT(IN) :: ctx
1370 : INTEGER, INTENT(OUT) :: nentry
1371 :
1372 25820 : CPASSERT(ALLOCATED(ctx%rs_route))
1373 25820 : nentry = SIZE(ctx%rs_route%entry, 2)
1374 25820 : END SUBROUTINE kpoint_operator_work_size
1375 :
1376 : ! **************************************************************************************************
1377 : !> \brief Evaluate an immutable prepared operator with exclusively caller-owned output and scratch.
1378 : !> Prepare/refresh/release remain collective and serialized. Evaluate contains no MPI,
1379 : !> DBCSR lookup, window update or nested OpenMP region. Source blocks must stay unchanged.
1380 : !> \param ctx prepared read-only context
1381 : !> \param ispin instance spin
1382 : !> \param slot operator slot
1383 : !> \param ik global k index owned by this group
1384 : !> \param phase private phase scratch, sized with work_size
1385 : !> \param TARGET optional private complex singleton output
1386 : !> \param fm_re optional private real output instead of TARGET
1387 : !> \param fm_im optional private imaginary output, accompanying fm_re for a complex operator
1388 : ! **************************************************************************************************
1389 52793 : SUBROUTINE kpoint_operator_evaluate(ctx, ispin, slot, ik, phase, TARGET, fm_re, fm_im)
1390 : TYPE(kpoint_operator_context_type), INTENT(IN) :: ctx
1391 : INTEGER, INTENT(IN) :: ispin, slot, ik
1392 : TYPE(cp_cfm_type), INTENT(INOUT), OPTIONAL :: target
1393 : COMPLEX(KIND=dp), INTENT(INOUT), TARGET :: phase(:, :)
1394 : TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL :: fm_re, fm_im
1395 :
1396 52793 : REAL(KIND=dp), POINTER :: xkp(:, :)
1397 : INTEGER :: row
1398 :
1399 0 : CPASSERT(ik >= ctx%kpoints%kp_range(1))
1400 52793 : CPASSERT(ik <= ctx%kpoints%kp_range(2))
1401 52793 : row = ispin
1402 52793 : IF (.NOT. ASSOCIATED(ctx%rs(row, slot)%source)) row = 1
1403 52793 : CPASSERT(ASSOCIATED(ctx%rs(row, slot)%source))
1404 52793 : xkp => ctx%kpoints%xkp
1405 : CALL kpoint_operator_rs_build_batch(ctx%rs(row, slot), ctx%rs_route, xkp(:, ik:ik), &
1406 : TARGET=TARGET, fm_re=fm_re, fm_im=fm_im, &
1407 52793 : phase_work=phase, parallel_columns=.FALSE.)
1408 52793 : IF (ctx%kpoints%use_real_wfn) THEN
1409 0 : IF (PRESENT(TARGET)) target%local_data(:, :) = CMPLX(REAL(target%local_data, dp), 0.0_dp, dp)
1410 0 : IF (PRESENT(fm_im)) fm_im%local_data(:, :) = 0.0_dp
1411 : END IF
1412 52793 : END SUBROUTINE kpoint_operator_evaluate
1413 :
1414 : ! **************************************************************************************************
1415 : !> \brief Start one operator for one global round and spin. The service walks
1416 : !> all groups with a k point in this round; groups without one still
1417 : !> participate in the transfers. Every rank must iterate through the
1418 : !> maximum local k-point count. Finish every started instance with
1419 : !> kpoint_operator_finish. The first start also creates the slot
1420 : !> bookkeeping of the context.
1421 : !> Collective over the environment-wide communicator: every rank must
1422 : !> call the same sequence of service calls.
1423 : !> \param ctx context created by kpoint_operator_context_create
1424 : !> \param ikp global round, up to the maximum k-point count across groups
1425 : !> \param ispin instance spin: keys the bookkeeping
1426 : !> \param rsmat real-space image matrices of the operator. A prepared operator
1427 : !> batch uses the values captured by prepare; otherwise start reads
1428 : !> the current values. The matching finish never reads rsmat
1429 : !> \param slot instance slot, optionally backed by a prepared snapshot
1430 : !> \param grid prepared reciprocal-grid cache of the caller
1431 : !> \param use_grid extract from grid instead of a direct phase sum
1432 : !> \param matrix_row row of rsmat to transform. Spin-free operators such
1433 : !> as S and T must pass kpoint_spin_free: their first matrix
1434 : !> dimension counts derivatives, not spins. A spin index is out of
1435 : !> bounds in plain SCF, but with derivatives stored it reads another
1436 : !> valid row silently
1437 : ! **************************************************************************************************
1438 28232 : SUBROUTINE kpoint_operator_start(ctx, ikp, ispin, rsmat, slot, grid, use_grid, matrix_row)
1439 :
1440 : TYPE(kpoint_operator_context_type), INTENT(INOUT) :: ctx
1441 : INTEGER, INTENT(IN) :: ikp, ispin
1442 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
1443 : INTEGER, INTENT(IN) :: slot
1444 : TYPE(rskp_grid_type), INTENT(IN), OPTIONAL :: grid
1445 : LOGICAL, INTENT(IN) :: use_grid
1446 : INTEGER, INTENT(IN) :: matrix_row
1447 :
1448 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_operator_start'
1449 :
1450 : INTEGER :: handle, igroup, ik, irs, kpmax
1451 28232 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1452 :
1453 28232 : CALL timeset(routineN, handle)
1454 :
1455 28232 : CPASSERT(ASSOCIATED(ctx%kpoints))
1456 28232 : CPASSERT(ASSOCIATED(ctx%kpoints%sab_nl))
1457 : ! The instance spin is bounded by the spin space declared at creation;
1458 : ! the matrix row by this operator's own rows. The two bounds differ for
1459 : ! spin-free operators such as S, whose first dimension counts
1460 : ! real-space blocks.
1461 28232 : CPASSERT(ispin >= 1 .AND. ispin <= ctx%nspin)
1462 28232 : CPASSERT(matrix_row >= 1 .AND. matrix_row <= SIZE(rsmat, 1))
1463 56472 : kpmax = MAXVAL(ctx%kpoints%kp_dist(2, :) - ctx%kpoints%kp_dist(1, :) + 1)
1464 28232 : CPASSERT(ikp >= 1 .AND. ikp <= kpmax)
1465 28232 : CPASSERT(slot >= kpoint_slot_ks .AND. slot <= kpoint_slot_t)
1466 : ! State is per round; allocate expensive copy descriptors only when
1467 : ! an actual distributed transfer is launched.
1468 28232 : IF (.NOT. ALLOCATED(ctx%status)) THEN
1469 51810 : ALLOCATE (ctx%status(kpmax, ctx%nspin, kpoint_num_slots))
1470 150450 : ctx%status(:, :, :) = kpoint_slot_idle
1471 : END IF
1472 28232 : CPASSERT(ctx%status(ikp, ispin, slot) == kpoint_slot_idle)
1473 28232 : IF (.NOT. use_grid) THEN
1474 28232 : irs = kpoint_operator_rs_match(ctx, rsmat, ispin, slot, matrix_row)
1475 28232 : IF (irs > 0) THEN
1476 28192 : CPASSERT(ASSOCIATED(ctx%work_sab, ctx%kpoints%sab_nl))
1477 28192 : ctx%status(ikp, ispin, slot) = kpoint_slot_local
1478 28192 : CALL timestop(handle)
1479 28192 : RETURN
1480 : END IF
1481 : END IF
1482 40 : CALL kpoint_operator_ensure(ctx, rsmat)
1483 :
1484 : ! Only the transfer backend needs copy descriptors. A prepared-only
1485 : ! or get-only context never carries them. Size by the declared spin
1486 : ! space, not by the possibly spin-free operator's row count.
1487 40 : IF (.NOT. ALLOCATED(ctx%info)) THEN
1488 : ! One copy_info per in-flight part: start takes it INTENT(OUT),
1489 : ! so two concurrent starts cannot share one instance.
1490 0 : ALLOCATE (ctx%info(SIZE(ctx%kpoints%xkp, 2), ctx%nspin, &
1491 380 : kpoint_num_slots, MERGE(1, 2, ctx%kpoints%use_real_wfn)))
1492 : END IF
1493 :
1494 40 : kp_dist => ctx%kpoints%kp_dist
1495 80 : DO igroup = 1, ctx%kpoints%nkp_groups
1496 40 : ik = kp_dist(1, igroup) + ikp - 1
1497 40 : IF (ik > kp_dist(2, igroup)) CYCLE
1498 40 : CALL kpoint_operator_densify(ctx, rsmat, matrix_row, ik, grid, use_grid)
1499 : CALL kpoint_operator_launch(ctx%kpoints, ik, ctx%fmwork, ctx%fmlocal, &
1500 80 : ctx%info(ik, ispin, slot, :))
1501 : END DO
1502 40 : ctx%status(ikp, ispin, slot) = kpoint_slot_in_flight
1503 :
1504 40 : CALL timestop(handle)
1505 :
1506 28232 : END SUBROUTINE kpoint_operator_start
1507 :
1508 : ! **************************************************************************************************
1509 : !> \brief Finish one started instance: build its prepared group-local operator,
1510 : !> or merge received parts into cmat (real part first, overwrite then add)
1511 : !> and reclaim the send side on every group.
1512 : !> The slot returns to idle.
1513 : !> Collective over the environment-wide communicator: every rank must
1514 : !> call the same sequence of service calls.
1515 : !> \param ctx context that started the instance
1516 : !> \param ikp global round of the start call
1517 : !> \param ispin spin component of the start call
1518 : !> \param slot slot of the start call; the slot must be in flight
1519 : !> \param cmat caller-allocated complex matrix that receives the operator.
1520 : !> It must share the distribution of the k-point group. Unchanged on
1521 : !> groups without a k point in this round
1522 : !> \param fm_re optional real FM output instead of cmat for real wavefunctions, on the group layout
1523 : ! **************************************************************************************************
1524 28232 : SUBROUTINE kpoint_operator_finish(ctx, ikp, ispin, slot, cmat, fm_re)
1525 :
1526 : TYPE(kpoint_operator_context_type), INTENT(INOUT) :: ctx
1527 : INTEGER, INTENT(IN) :: ikp, ispin, slot
1528 : TYPE(cp_cfm_type), OPTIONAL :: cmat
1529 : TYPE(cp_fm_type), OPTIONAL :: fm_re
1530 :
1531 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_operator_finish'
1532 :
1533 : INTEGER :: handle, igroup, ik, irs, kpmax
1534 : INTEGER, DIMENSION(2) :: kp_range
1535 28232 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1536 : LOGICAL :: layout_ok
1537 :
1538 28232 : CALL timeset(routineN, handle)
1539 :
1540 28232 : CPASSERT(ASSOCIATED(ctx%kpoints))
1541 28232 : CPASSERT(ALLOCATED(ctx%status))
1542 28232 : kpmax = SIZE(ctx%status, 1)
1543 28232 : CPASSERT(ikp >= 1 .AND. ikp <= kpmax)
1544 28232 : CPASSERT(ispin >= 1 .AND. ispin <= ctx%nspin)
1545 28232 : CPASSERT(slot >= kpoint_slot_ks .AND. slot <= kpoint_slot_t)
1546 28232 : CPASSERT(ctx%status(ikp, ispin, slot) /= kpoint_slot_idle)
1547 : ! The merge writes cmat in place. A foreign distribution would fail
1548 : ! inside the BLACS calls, far from this argument.
1549 :
1550 28232 : CPASSERT(PRESENT(cmat) .NEQV. PRESENT(fm_re))
1551 28232 : CPASSERT(.NOT. PRESENT(fm_re) .OR. ctx%kpoints%use_real_wfn)
1552 84696 : kp_range = ctx%kpoints%kp_range
1553 28232 : IF (ctx%status(ikp, ispin, slot) == kpoint_slot_local) THEN
1554 28192 : irs = ispin
1555 28192 : IF (.NOT. ASSOCIATED(ctx%rs(irs, slot)%source)) irs = 1
1556 28192 : ik = kp_range(1) + ikp - 1
1557 28192 : IF (ik <= kp_range(2)) THEN
1558 28192 : CALL kpoint_operator_rs_build(ctx, irs, slot, ik, cmat, fm_re=fm_re)
1559 : END IF
1560 28192 : ctx%status(ikp, ispin, slot) = kpoint_slot_idle
1561 28192 : CALL timestop(handle)
1562 28192 : RETURN
1563 : END IF
1564 40 : CPASSERT(ctx%status(ikp, ispin, slot) == kpoint_slot_in_flight)
1565 40 : IF (PRESENT(cmat)) THEN
1566 24 : layout_ok = cp_fm_struct_equivalent(cmat%matrix_struct, ctx%fmlocal%matrix_struct)
1567 : ELSE
1568 16 : layout_ok = cp_fm_struct_equivalent(fm_re%matrix_struct, ctx%fmlocal%matrix_struct)
1569 : END IF
1570 40 : CPASSERT(layout_ok)
1571 40 : kp_dist => ctx%kpoints%kp_dist
1572 80 : DO igroup = 1, ctx%kpoints%nkp_groups
1573 40 : ik = kp_dist(1, igroup) + ikp - 1
1574 40 : IF (ik > kp_dist(2, igroup)) CYCLE
1575 : CALL kpoint_operator_receive(ctx%kpoints, ctx%fmlocal, ik, ctx%info(ik, ispin, slot, :), &
1576 80 : cmat=cmat, fm_re=fm_re)
1577 : END DO
1578 40 : ctx%status(ikp, ispin, slot) = kpoint_slot_idle
1579 :
1580 40 : CALL timestop(handle)
1581 :
1582 28232 : END SUBROUTINE kpoint_operator_finish
1583 :
1584 : ! **************************************************************************************************
1585 : !> \brief Assemble one operator for one global round and spin in one call:
1586 : !> run the single-group transfer of the service for every group and
1587 : !> finish it immediately, through transfer states local to this call.
1588 : !> The persistent slot machinery stays with start and finish.
1589 : !> Exactly one output flavor must be present: cmat receives the complex
1590 : !> operator (real part overwrites, imaginary part adds) through the
1591 : !> context buffer; the fm pair receives the two parts separately,
1592 : !> transferred directly into fm_re and fm_im on the owning group (any
1593 : !> distribution: the general copy redistributes); or the dbcsr pair
1594 : !> receives the two parts as block-cyclic matrices built by the service,
1595 : !> replacing the caller-side copy_fm_to_dbcsr_bc detour. The dbcsr
1596 : !> matrices must be freshly initialized (dbcsr_init_p) and the cmat
1597 : !> flavor must share the distribution of the k-point group.
1598 : !> Collective over the environment-wide communicator: every rank must
1599 : !> call the same sequence of service calls.
1600 : !> \param ctx context created by kpoint_operator_context_create
1601 : !> \param ikp global round, up to the maximum k-point count across groups
1602 : !> \param ispin instance spin: bounds-checked against the context spin space
1603 : !> \param rsmat real-space image matrices of the operator
1604 : !> \param cmat complex matrix output
1605 : !> \param fm_re real-part matrix output
1606 : !> \param fm_im imaginary-part output; may be omitted for real wavefunctions
1607 : !> \param mat_re real-part block-cyclic dbcsr output
1608 : !> \param mat_im imaginary-part block-cyclic dbcsr output
1609 : !> \param matrix_row row of rsmat to transform. Spin-free operators such
1610 : !> as S and T must pass kpoint_spin_free: their first matrix
1611 : !> dimension counts derivatives, not spins. A spin index is out of
1612 : !> bounds in plain SCF, but with derivatives stored it reads another
1613 : !> valid row silently
1614 : !> \param slot optional prepared slot: opt in to its immutable snapshot.
1615 : !> Without slot, get always reads current source values. A missing or
1616 : !> mismatched snapshot retains the original transfer path. Prepared FM
1617 : !> targets must match the layout selected at prepare; inactive groups
1618 : !> do not inspect or modify their output arguments.
1619 : !> \param grid optional lattice-FFT grid, overriding a prepared snapshot for this request
1620 : ! **************************************************************************************************
1621 5586 : SUBROUTINE kpoint_operator_get(ctx, ikp, ispin, rsmat, cmat, fm_re, fm_im, &
1622 : mat_re, mat_im, matrix_row, slot, grid)
1623 :
1624 : TYPE(kpoint_operator_context_type), INTENT(INOUT) :: ctx
1625 : INTEGER, INTENT(IN) :: ikp, ispin
1626 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
1627 : TYPE(cp_cfm_type), OPTIONAL :: cmat
1628 : TYPE(cp_fm_type), OPTIONAL :: fm_re, fm_im
1629 : TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL :: mat_re, mat_im
1630 : INTEGER, INTENT(IN) :: matrix_row
1631 : INTEGER, INTENT(IN), OPTIONAL :: slot
1632 : TYPE(rskp_grid_type), INTENT(IN), OPTIONAL :: grid
1633 :
1634 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_operator_get'
1635 :
1636 : INTEGER :: handle, igroup, ik, irs, kpmax
1637 5586 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1638 : LOGICAL :: layout_ok
1639 67032 : TYPE(copy_info_type) :: info(2)
1640 :
1641 5586 : CALL timeset(routineN, handle)
1642 :
1643 22344 : CPASSERT(COUNT([PRESENT(cmat), PRESENT(fm_re), PRESENT(mat_re)]) == 1)
1644 5586 : IF (PRESENT(fm_re) .OR. PRESENT(fm_im)) THEN
1645 1114 : CPASSERT(PRESENT(fm_re))
1646 1114 : CPASSERT(PRESENT(fm_im) .OR. ctx%kpoints%use_real_wfn)
1647 : END IF
1648 5586 : IF (PRESENT(mat_re) .OR. PRESENT(mat_im)) THEN
1649 272 : CPASSERT(PRESENT(mat_re) .AND. PRESENT(mat_im))
1650 : END IF
1651 :
1652 5586 : CPASSERT(ASSOCIATED(ctx%kpoints))
1653 5586 : CPASSERT(ASSOCIATED(ctx%kpoints%sab_nl))
1654 5586 : CPASSERT(ispin >= 1 .AND. ispin <= ctx%nspin)
1655 5586 : CPASSERT(matrix_row >= 1 .AND. matrix_row <= SIZE(rsmat, 1))
1656 13528 : kpmax = MAXVAL(ctx%kpoints%kp_dist(2, :) - ctx%kpoints%kp_dist(1, :) + 1)
1657 5586 : CPASSERT(ikp >= 1 .AND. ikp <= kpmax)
1658 5586 : IF (PRESENT(slot) .AND. .NOT. PRESENT(grid)) THEN
1659 5586 : CPASSERT(slot >= kpoint_slot_ks .AND. slot <= kpoint_slot_t)
1660 5586 : irs = kpoint_operator_rs_match(ctx, rsmat, ispin, slot, matrix_row)
1661 5586 : IF (irs > 0) THEN
1662 5586 : ik = ctx%kpoints%kp_range(1) + ikp - 1
1663 5586 : IF (ik <= ctx%kpoints%kp_range(2)) THEN
1664 5584 : IF (PRESENT(mat_re)) THEN
1665 270 : IF (.NOT. ALLOCATED(ctx%fmimag)) THEN
1666 82 : ALLOCATE (ctx%fmimag)
1667 82 : CALL cp_fm_create(ctx%fmimag, ctx%fmlocal%matrix_struct)
1668 : END IF
1669 270 : CALL kpoint_operator_rs_build(ctx, irs, slot, ik, fm_re=ctx%fmlocal, fm_im=ctx%fmimag)
1670 270 : CALL copy_fm_to_dbcsr_bc(ctx%fmlocal, mat_re)
1671 270 : CALL copy_fm_to_dbcsr_bc(ctx%fmimag, mat_im)
1672 : ELSE
1673 5314 : CALL kpoint_operator_rs_build(ctx, irs, slot, ik, TARGET=cmat, fm_re=fm_re, fm_im=fm_im)
1674 : END IF
1675 : END IF
1676 5586 : CALL timestop(handle)
1677 5586 : RETURN
1678 : END IF
1679 : END IF
1680 0 : CALL kpoint_operator_ensure(ctx, rsmat)
1681 0 : IF (PRESENT(cmat)) THEN
1682 : ! The merge writes cmat in place. A foreign distribution would fail
1683 : ! inside the BLACS calls, far from this argument.
1684 0 : layout_ok = cp_fm_struct_equivalent(cmat%matrix_struct, ctx%fmlocal%matrix_struct)
1685 0 : CPASSERT(layout_ok)
1686 : END IF
1687 :
1688 0 : kp_dist => ctx%kpoints%kp_dist
1689 0 : DO igroup = 1, ctx%kpoints%nkp_groups
1690 0 : ik = kp_dist(1, igroup) + ikp - 1
1691 0 : IF (ik > kp_dist(2, igroup)) CYCLE
1692 : ! Absent fm targets are forwarded as absent optional arguments. Both
1693 : ! CFM and DBCSR outputs use the same context-owned receive buffer.
1694 0 : CALL kpoint_operator_densify(ctx, rsmat, matrix_row, ik, grid=grid, use_grid=PRESENT(grid))
1695 : CALL kpoint_operator_launch(ctx%kpoints, ik, ctx%fmwork, ctx%fmlocal, info, &
1696 0 : target_re=fm_re, target_im=fm_im)
1697 : CALL kpoint_operator_receive(ctx%kpoints, ctx%fmlocal, ik, info, cmat=cmat, fm_re=fm_re, fm_im=fm_im, &
1698 0 : mat_re=mat_re, mat_im=mat_im)
1699 : END DO
1700 :
1701 0 : CALL timestop(handle)
1702 :
1703 5586 : END SUBROUTINE kpoint_operator_get
1704 :
1705 : ! **************************************************************************************************
1706 : !> \brief Read an operator inside a group-local OT loop, after collective cache preparation.
1707 : !> One group retains the sparse transform and distribution required by OT's DBCSR kernels.
1708 : !> Multiple groups read their owned cache; no world collective may occur in this loop.
1709 : !> \param matrix_rs current real-space operator
1710 : !> \param kpoints host environment
1711 : !> \param kp local k-point state
1712 : !> \param ispin source row (one for spin-free S/T)
1713 : !> \param cache_re owned real part, used with multiple groups
1714 : !> \param cache_im owned imaginary part, used with multiple groups
1715 : !> \param matrix_re newly initialized sparse real output
1716 : !> \param matrix_im newly initialized sparse imaginary output
1717 : ! **************************************************************************************************
1718 9636 : SUBROUTINE kpoint_operator_get_local(matrix_rs, kpoints, kp, ispin, cache_re, cache_im, &
1719 : matrix_re, matrix_im)
1720 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_rs
1721 : TYPE(kpoint_type), POINTER :: kpoints
1722 : TYPE(kpoint_env_type), POINTER :: kp
1723 : INTEGER, INTENT(IN) :: ispin
1724 : TYPE(cp_fm_type), INTENT(IN) :: cache_re, cache_im
1725 : TYPE(dbcsr_type), POINTER :: matrix_re, matrix_im
1726 :
1727 : INTEGER :: nkp_groups
1728 9636 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1729 9636 : REAL(dp), DIMENSION(:, :), POINTER :: xkp
1730 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1731 9636 : POINTER :: sab_nl
1732 :
1733 9636 : NULLIFY (cell_to_index, sab_nl, xkp)
1734 0 : CPASSERT(ASSOCIATED(matrix_rs))
1735 9636 : CPASSERT(ispin >= 1 .AND. ispin <= SIZE(matrix_rs, 1))
1736 9636 : CPASSERT(ASSOCIATED(kpoints))
1737 9636 : CPASSERT(ASSOCIATED(kp))
1738 :
1739 9636 : CALL dbcsr_init_p(matrix_re)
1740 9636 : CALL dbcsr_init_p(matrix_im)
1741 9636 : nkp_groups = kpoints%nkp_groups
1742 9636 : xkp => kpoints%xkp
1743 9636 : cell_to_index => kpoints%cell_to_index
1744 9636 : sab_nl => kpoints%sab_nl
1745 9636 : IF (nkp_groups == 1) THEN
1746 6482 : CPASSERT(ASSOCIATED(xkp))
1747 6482 : CPASSERT(ASSOCIATED(cell_to_index))
1748 6482 : CPASSERT(ASSOCIATED(sab_nl))
1749 : CALL dbcsr_create(matrix_re, template=matrix_rs(1, 1)%matrix, &
1750 6482 : matrix_type=dbcsr_type_symmetric)
1751 : CALL dbcsr_create(matrix_im, template=matrix_rs(1, 1)%matrix, &
1752 6482 : matrix_type=dbcsr_type_antisymmetric)
1753 6482 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_re, sab_nl)
1754 6482 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_im, sab_nl)
1755 6482 : CALL dbcsr_set(matrix_re, 0.0_dp)
1756 6482 : CALL dbcsr_set(matrix_im, 0.0_dp)
1757 : CALL rskp_transform(rmatrix=matrix_re, cmatrix=matrix_im, rsmat=matrix_rs, ispin=ispin, &
1758 6482 : xkp=xkp(1:3, kp%nkpoint), cell_to_index=cell_to_index, sab_nl=sab_nl)
1759 : ELSE
1760 3154 : CALL copy_fm_to_dbcsr_bc(cache_re, matrix_re)
1761 3154 : CALL copy_fm_to_dbcsr_bc(cache_im, matrix_im)
1762 : END IF
1763 9636 : END SUBROUTINE kpoint_operator_get_local
1764 :
1765 : ! **************************************************************************************************
1766 : !> \brief Fill caller-owned matrices for one operator row, including independent history snapshots.
1767 : !> Prepare and fallback transfers are collective. Prepared one-rank groups evaluate k points
1768 : !> concurrently, with no MPI or matrix allocation/refcounts inside the worker loop.
1769 : !> \param kpoints host environment
1770 : !> \param ao_ao_fm global-layout template
1771 : !> \param rsmat current real-space images
1772 : !> \param slot operator slot
1773 : !> \param cmat allocated group-local outputs, one per locally owned k point
1774 : !> \param matrix_row source row; S/T use kpoint_spin_free
1775 : ! **************************************************************************************************
1776 480 : SUBROUTINE kpoint_operator_get_all(kpoints, ao_ao_fm, rsmat, slot, cmat, matrix_row)
1777 : TYPE(kpoint_type), POINTER :: kpoints
1778 : TYPE(cp_fm_type), INTENT(IN) :: ao_ao_fm
1779 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
1780 : INTEGER, INTENT(IN) :: slot
1781 : TYPE(cp_cfm_type), DIMENSION(:), INTENT(INOUT) :: cmat
1782 : INTEGER, INTENT(IN) :: matrix_row
1783 :
1784 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_operator_get_all'
1785 :
1786 480 : COMPLEX(KIND=dp), ALLOCATABLE :: phase(:, :)
1787 : INTEGER :: handle, ik, ikp, kplocal, kpmax, nentry, &
1788 : nworkers
1789 : LOGICAL :: layout_ok, prepared
1790 480 : TYPE(kpoint_operator_context_type) :: ctx
1791 :
1792 480 : CALL timeset(routineN, handle)
1793 480 : kplocal = SIZE(kpoints%kp_env)
1794 480 : CPASSERT(SIZE(cmat) == kplocal .AND. kplocal > 0)
1795 480 : CALL kpoint_operator_context_create(ctx, kpoints, ao_ao_fm, 1)
1796 480 : CALL kpoint_operator_prepare(ctx, rsmat, slot, prepared, matrix_row=matrix_row)
1797 480 : IF (prepared .AND. kpoints%para_env_kp%num_pe == 1) THEN
1798 1242 : DO ikp = 1, kplocal
1799 916 : layout_ok = cp_fm_struct_equivalent(cmat(ikp)%matrix_struct, ctx%rs_route%matrix_struct)
1800 1242 : CPASSERT(layout_ok)
1801 : END DO
1802 326 : CALL kpoint_operator_work_size(ctx, nentry)
1803 326 : nworkers = 1
1804 326 : !$ nworkers = MIN(omp_get_max_threads(), kplocal)
1805 : !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nworkers) SHARED(ctx,kpoints,cmat,slot,kplocal,nentry) &
1806 326 : !$OMP PRIVATE(phase,ikp,ik)
1807 : ALLOCATE (phase(1, nentry))
1808 : !$OMP DO SCHEDULE(DYNAMIC, 1)
1809 : DO ikp = 1, kplocal
1810 : ik = kpoints%kp_range(1) + ikp - 1
1811 : CALL kpoint_operator_evaluate(ctx, 1, slot, ik, phase, cmat(ikp))
1812 : END DO
1813 : !$OMP END DO
1814 326 : DEALLOCATE (phase)
1815 : !$OMP END PARALLEL
1816 : ELSE
1817 308 : kpmax = MAXVAL(kpoints%kp_dist(2, :) - kpoints%kp_dist(1, :) + 1)
1818 398 : DO ikp = 1, kpmax
1819 : CALL kpoint_operator_get(ctx, ikp, 1, rsmat, cmat=cmat(MIN(ikp, kplocal)), &
1820 398 : matrix_row=matrix_row, slot=slot)
1821 : END DO
1822 : END IF
1823 480 : CALL kpoint_operator_context_release(ctx)
1824 480 : CALL timestop(handle)
1825 480 : END SUBROUTINE kpoint_operator_get_all
1826 :
1827 : ! **************************************************************************************************
1828 : !> \brief Prepare owned split-FM operators for OT without entering an SCF eigensolver.
1829 : !> H (or the initialization density) retains all source spins; S/T are built only once.
1830 : !> Outputs survive the transient context. Compatible cache storage is reused, values refreshed.
1831 : !> \param kpoints host environment and owned OT caches
1832 : !> \param ao_ao_fm global-layout template
1833 : !> \param matrix_ks Hamiltonian or initialization density
1834 : !> \param matrix_s spin-free overlap
1835 : !> \param matrix_t optional spin-free kinetic operator; absent leaves the existing kinetic cache alone
1836 : ! **************************************************************************************************
1837 1240 : SUBROUTINE kpoint_operator_store(kpoints, ao_ao_fm, matrix_ks, matrix_s, matrix_t)
1838 : TYPE(kpoint_type), POINTER :: kpoints
1839 : TYPE(cp_fm_type), INTENT(IN) :: ao_ao_fm
1840 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
1841 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
1842 : POINTER :: matrix_t
1843 :
1844 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_operator_store'
1845 : INTEGER :: handle, ic, ikp, ispin, ncomp, nspin
1846 : LOGICAL :: kinetic, reuse
1847 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1848 : TYPE(kpoint_env_type), POINTER :: kp
1849 1240 : TYPE(kpoint_operator_context_type) :: ctx
1850 :
1851 1240 : CALL timeset(routineN, handle)
1852 1240 : nspin = SIZE(matrix_ks, 1)
1853 1240 : ncomp = MERGE(1, 2, kpoints%use_real_wfn)
1854 1240 : kinetic = PRESENT(matrix_t)
1855 1240 : IF (kinetic) kinetic = ASSOCIATED(matrix_t)
1856 1240 : CALL kpoint_operator_context_create(ctx, kpoints, ao_ao_fm, nspin)
1857 1240 : CALL kpoint_operator_ensure_group(ctx)
1858 1240 : matrix_struct => ctx%fmlocal%matrix_struct
1859 3099 : DO ikp = 1, SIZE(kpoints%kp_env)
1860 1859 : kp => kpoints%kp_env(ikp)%kpoint_env
1861 1859 : reuse = .FALSE.
1862 1859 : IF (ASSOCIATED(kp%ot_hmat)) THEN
1863 4986 : IF (ALL(SHAPE(kp%ot_hmat) == [ncomp, nspin])) THEN
1864 1662 : reuse = cp_fm_struct_equivalent(kp%ot_hmat(1, 1)%matrix_struct, matrix_struct)
1865 : END IF
1866 : END IF
1867 1662 : IF (.NOT. reuse) THEN
1868 197 : CALL cp_fm_release(kp%ot_hmat)
1869 1472 : ALLOCATE (kp%ot_hmat(ncomp, nspin))
1870 425 : DO ispin = 1, nspin
1871 881 : DO ic = 1, ncomp
1872 684 : CALL cp_fm_create(kp%ot_hmat(ic, ispin), matrix_struct)
1873 : END DO
1874 : END DO
1875 : END IF
1876 1859 : CALL allocate_spin_free(kp%ot_smat)
1877 3099 : IF (kinetic) CALL allocate_spin_free(kp%ot_tmat)
1878 : END DO
1879 1240 : CALL store_slot(matrix_ks, kpoint_slot_ks, nspin)
1880 1240 : CALL store_slot(matrix_s, kpoint_slot_s, 1)
1881 1240 : IF (kinetic) CALL store_slot(matrix_t, kpoint_slot_t, 1)
1882 1240 : CALL kpoint_operator_context_release(ctx)
1883 1240 : CALL timestop(handle)
1884 :
1885 : CONTAINS
1886 :
1887 : ! **************************************************************************************************
1888 : !> \brief ...
1889 : !> \param cache ...
1890 : ! **************************************************************************************************
1891 2058 : SUBROUTINE allocate_spin_free(cache)
1892 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: cache
1893 :
1894 : INTEGER :: ic
1895 : LOGICAL :: reuse
1896 :
1897 2058 : reuse = .FALSE.
1898 2058 : IF (ASSOCIATED(cache)) THEN
1899 1702 : IF (SIZE(cache) == ncomp) reuse = cp_fm_struct_equivalent(cache(1)%matrix_struct, matrix_struct)
1900 : END IF
1901 1702 : IF (reuse) RETURN
1902 356 : CALL cp_fm_release(cache)
1903 1780 : ALLOCATE (cache(ncomp))
1904 1068 : DO ic = 1, ncomp
1905 1068 : CALL cp_fm_create(cache(ic), matrix_struct)
1906 : END DO
1907 : END SUBROUTINE allocate_spin_free
1908 :
1909 : ! **************************************************************************************************
1910 : !> \brief ...
1911 : !> \param rsmat ...
1912 : !> \param slot ...
1913 : !> \param nrows ...
1914 : ! **************************************************************************************************
1915 2606 : SUBROUTINE store_slot(rsmat, slot, nrows)
1916 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
1917 : INTEGER, INTENT(IN) :: slot, nrows
1918 :
1919 2606 : COMPLEX(KIND=dp), ALLOCATABLE :: phase(:, :)
1920 : INTEGER :: ikp, ispin, kplocal, kpmax, nentry, &
1921 : nworkers
1922 : LOGICAL :: prepared, ready, use_grid
1923 2606 : TYPE(rskp_grid_type), ALLOCATABLE :: grids(:)
1924 :
1925 2606 : nworkers = 1
1926 2606 : !$ nworkers = omp_get_max_threads()
1927 : use_grid = kpoints%lattice_fft == lattice_fft_on .OR. &
1928 2606 : (kpoints%lattice_fft == lattice_fft_auto .AND. SIZE(kpoints%xkp, 2) >= 27)
1929 2606 : prepared = .FALSE.
1930 : ! Match STANDARD's single-world-rank, single-thread preference for the forward lattice FFT.
1931 2606 : IF (.NOT. use_grid .OR. nworkers > 1 .OR. kpoints%blacs_env_all%para_env%num_pe > 1) THEN
1932 2606 : IF (slot == kpoint_slot_ks) THEN
1933 1240 : CALL kpoint_operator_prepare(ctx, rsmat, slot, prepared)
1934 : ELSE
1935 1366 : CALL kpoint_operator_prepare(ctx, rsmat, slot, prepared, matrix_row=kpoint_spin_free)
1936 : END IF
1937 : END IF
1938 2606 : kplocal = SIZE(kpoints%kp_env)
1939 2606 : IF (prepared .AND. kpoints%para_env_kp%num_pe == 1) THEN
1940 2044 : CALL kpoint_operator_work_size(ctx, nentry)
1941 2044 : nworkers = MIN(nworkers, kplocal*nrows)
1942 : !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nworkers) SHARED(ctx,kpoints,rsmat,slot,nrows,kplocal,nentry) &
1943 2044 : !$OMP PRIVATE(phase,ikp,ispin)
1944 : ALLOCATE (phase(1, nentry))
1945 : !$OMP DO COLLAPSE(2) SCHEDULE(DYNAMIC, 1)
1946 : DO ikp = 1, kplocal
1947 : DO ispin = 1, nrows
1948 : CALL store_round(rsmat, slot, ikp, ispin, phase)
1949 : END DO
1950 : END DO
1951 : !$OMP END DO
1952 2044 : DEALLOCATE (phase)
1953 : !$OMP END PARALLEL
1954 : ELSE
1955 1124 : kpmax = MAXVAL(kpoints%kp_dist(2, :) - kpoints%kp_dist(1, :) + 1)
1956 : ! Keep the existing forward lattice-FFT choice when distributed snapshots cannot be prepared.
1957 562 : use_grid = use_grid .AND. .NOT. prepared
1958 : IF (use_grid) THEN
1959 0 : CALL kpoint_operator_ensure(ctx, rsmat)
1960 0 : ALLOCATE (grids(nrows))
1961 0 : DO ispin = 1, nrows
1962 : CALL rskp_transform_grid_prepare(grids(ispin), ctx%rmatrix, rsmat, ispin, &
1963 : kpoints%xkp, kpoints%nkp_grid, kpoints%cell_to_index, &
1964 : kpoints%sab_nl, ready, &
1965 0 : max_storage_bytes=512_int_8*1024_int_8**2/INT(nspin + 1, int_8))
1966 0 : use_grid = use_grid .AND. ready
1967 : END DO
1968 : END IF
1969 1576 : DO ikp = 1, kpmax
1970 2636 : DO ispin = 1, nrows
1971 2074 : IF (use_grid) THEN
1972 0 : CALL store_round(rsmat, slot, ikp, ispin, grid=grids(ispin))
1973 : ELSE
1974 1060 : CALL store_round(rsmat, slot, ikp, ispin)
1975 : END IF
1976 : END DO
1977 : END DO
1978 562 : IF (ALLOCATED(grids)) THEN
1979 0 : DO ispin = 1, nrows
1980 0 : CALL rskp_transform_grid_release(grids(ispin))
1981 : END DO
1982 0 : DEALLOCATE (grids)
1983 : END IF
1984 : END IF
1985 2606 : END SUBROUTINE store_slot
1986 :
1987 : ! **************************************************************************************************
1988 : !> \brief ...
1989 : !> \param rsmat ...
1990 : !> \param slot ...
1991 : !> \param ikp ...
1992 : !> \param ispin ...
1993 : !> \param phase ...
1994 : !> \param grid ...
1995 : ! **************************************************************************************************
1996 4166 : SUBROUTINE store_round(rsmat, slot, ikp, ispin, phase, grid)
1997 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
1998 : INTEGER, INTENT(IN) :: slot, ikp, ispin
1999 : COMPLEX(KIND=dp), INTENT(INOUT), OPTIONAL :: phase(:, :)
2000 : TYPE(rskp_grid_type), INTENT(IN), OPTIONAL :: grid
2001 :
2002 : INTEGER :: ik
2003 4166 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: cache
2004 : TYPE(cp_fm_type), POINTER :: imag
2005 : TYPE(kpoint_env_type), POINTER :: kp
2006 :
2007 : ! A valid dummy cache is passed on inactive tail groups; get leaves it untouched.
2008 4166 : kp => kpoints%kp_env(MIN(ikp, SIZE(kpoints%kp_env)))%kpoint_env
2009 6274 : SELECT CASE (slot)
2010 : CASE (kpoint_slot_ks)
2011 2108 : cache => kp%ot_hmat(:, ispin)
2012 : CASE (kpoint_slot_s)
2013 1859 : cache => kp%ot_smat
2014 : CASE (kpoint_slot_t)
2015 199 : cache => kp%ot_tmat
2016 : END SELECT
2017 4166 : NULLIFY (imag)
2018 4166 : IF (SIZE(cache) == 2) imag => cache(2)
2019 4166 : IF (PRESENT(phase)) THEN
2020 3106 : ik = kpoints%kp_range(1) + ikp - 1
2021 3106 : CALL kpoint_operator_evaluate(ctx, ispin, slot, ik, phase, fm_re=cache(1), fm_im=imag)
2022 : ELSE
2023 : CALL kpoint_operator_get(ctx, ikp, ispin, rsmat, fm_re=cache(1), fm_im=imag, &
2024 1060 : matrix_row=ispin, slot=slot, grid=grid)
2025 : END IF
2026 4166 : END SUBROUTINE store_round
2027 : END SUBROUTINE kpoint_operator_store
2028 :
2029 0 : END MODULE qs_kpoint_operators
|
