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