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