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