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 : real_to_scaled
19 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
20 : cp_blacs_env_type
21 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale_and_add_fm
22 : USE cp_cfm_types, ONLY: cp_cfm_create,&
23 : cp_cfm_release,&
24 : cp_cfm_to_fm,&
25 : cp_cfm_type
26 : USE cp_control_types, ONLY: dft_control_type,&
27 : hairy_probes_type
28 : USE cp_dbcsr_api, ONLY: &
29 : dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
30 : dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_get_block_p, dbcsr_get_info, &
31 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
32 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_set, dbcsr_type, &
33 : dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
34 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
35 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
36 : copy_fm_to_dbcsr
37 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
38 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
39 : fm_pool_create_fm,&
40 : fm_pool_give_back_fm
41 : USE cp_fm_struct, ONLY: cp_fm_struct_type
42 : USE cp_fm_types, ONLY: &
43 : copy_info_type, cp_fm_cleanup_copy_general, cp_fm_create, cp_fm_finish_copy_general, &
44 : cp_fm_get_info, cp_fm_get_submatrix, cp_fm_release, cp_fm_start_copy_general, cp_fm_to_fm, &
45 : cp_fm_type
46 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
47 : USE cryssym, ONLY: crys_sym_gen,&
48 : csym_type,&
49 : kpoint_gen,&
50 : print_crys_symmetry,&
51 : print_kp_symmetry,&
52 : release_csym_type
53 : USE fermi_utils, ONLY: fermikp,&
54 : fermikp2
55 : USE hairy_probes, ONLY: probe_occupancy_kp
56 : USE input_constants, ONLY: smear_fermi_dirac
57 : USE kinds, ONLY: dp
58 : USE kpoint_types, ONLY: get_kpoint_info,&
59 : kind_rotmat_type,&
60 : kpoint_env_create,&
61 : kpoint_env_p_type,&
62 : kpoint_env_type,&
63 : kpoint_sym_create,&
64 : kpoint_sym_type,&
65 : kpoint_type
66 : USE mathconstants, ONLY: gaussi,&
67 : twopi
68 : USE memory_utilities, ONLY: reallocate
69 : USE message_passing, ONLY: mp_cart_type,&
70 : mp_para_env_type
71 : USE parallel_gemm_api, ONLY: parallel_gemm
72 : USE particle_types, ONLY: particle_type
73 : USE qs_matrix_pools, ONLY: mpools_create,&
74 : mpools_get,&
75 : mpools_rebuild_fm_pools,&
76 : qs_matrix_pools_type
77 : USE qs_mo_types, ONLY: allocate_mo_set,&
78 : get_mo_set,&
79 : init_mo_set,&
80 : mo_set_type,&
81 : set_mo_set
82 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
83 : get_neighbor_list_set_p,&
84 : neighbor_list_iterate,&
85 : neighbor_list_iterator_create,&
86 : neighbor_list_iterator_p_type,&
87 : neighbor_list_iterator_release,&
88 : neighbor_list_set_p_type
89 : USE scf_control_types, ONLY: smear_type
90 : USE util, ONLY: get_limit
91 : #include "./base/base_uses.f90"
92 :
93 : IMPLICIT NONE
94 :
95 : PRIVATE
96 :
97 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpoint_methods'
98 :
99 : PUBLIC :: kpoint_initialize, kpoint_env_initialize, kpoint_initialize_mos, kpoint_initialize_mo_set
100 : PUBLIC :: kpoint_init_cell_index, kpoint_set_mo_occupation
101 : PUBLIC :: kpoint_density_matrices, kpoint_density_transform
102 : PUBLIC :: rskp_transform, lowdin_kp_trans
103 :
104 : ! **************************************************************************************************
105 :
106 : CONTAINS
107 :
108 : ! **************************************************************************************************
109 : !> \brief Generate the kpoints and initialize the kpoint environment
110 : !> \param kpoint The kpoint environment
111 : !> \param particle_set Particle types and coordinates
112 : !> \param cell Computational cell information
113 : ! **************************************************************************************************
114 7682 : SUBROUTINE kpoint_initialize(kpoint, particle_set, cell)
115 :
116 : TYPE(kpoint_type), POINTER :: kpoint
117 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
118 : TYPE(cell_type), POINTER :: cell
119 :
120 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_initialize'
121 :
122 : INTEGER :: handle, i, ic, ik, iounit, ir, ira, is, &
123 : j, natom, nkind, nr, ns
124 7682 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atype
125 : LOGICAL :: spez
126 : REAL(KIND=dp) :: wsum
127 7682 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coord, scoord
128 6783206 : TYPE(csym_type) :: crys_sym
129 : TYPE(kpoint_sym_type), POINTER :: kpsym
130 :
131 7682 : CALL timeset(routineN, handle)
132 :
133 7682 : CPASSERT(ASSOCIATED(kpoint))
134 :
135 7698 : SELECT CASE (kpoint%kp_scheme)
136 : CASE ("NONE")
137 : ! do nothing
138 : CASE ("GAMMA")
139 16 : kpoint%nkp = 1
140 16 : ALLOCATE (kpoint%xkp(3, 1), kpoint%wkp(1))
141 64 : kpoint%xkp(1:3, 1) = 0.0_dp
142 16 : kpoint%wkp(1) = 1.0_dp
143 32 : ALLOCATE (kpoint%kp_sym(1))
144 16 : NULLIFY (kpoint%kp_sym(1)%kpoint_sym)
145 16 : CALL kpoint_sym_create(kpoint%kp_sym(1)%kpoint_sym)
146 : CASE ("MONKHORST-PACK", "MACDONALD")
147 :
148 162 : IF (.NOT. kpoint%symmetry) THEN
149 : ! we set up a random molecule to avoid any possible symmetry
150 112 : natom = 10
151 112 : ALLOCATE (coord(3, natom), scoord(3, natom), atype(natom))
152 1232 : DO i = 1, natom
153 1120 : atype(i) = i
154 1120 : coord(1, i) = SIN(i*0.12345_dp)
155 1120 : coord(2, i) = COS(i*0.23456_dp)
156 1120 : coord(3, i) = SIN(i*0.34567_dp)
157 1232 : CALL real_to_scaled(scoord(1:3, i), coord(1:3, i), cell)
158 : END DO
159 : ELSE
160 50 : natom = SIZE(particle_set)
161 250 : ALLOCATE (scoord(3, natom), atype(natom))
162 368 : DO i = 1, natom
163 318 : CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=atype(i))
164 368 : CALL real_to_scaled(scoord(1:3, i), particle_set(i)%r(1:3), cell)
165 : END DO
166 : END IF
167 162 : IF (kpoint%verbose) THEN
168 10 : iounit = cp_logger_get_default_io_unit()
169 : ELSE
170 152 : iounit = -1
171 : END IF
172 : ! kind type list
173 486 : ALLOCATE (kpoint%atype(natom))
174 1600 : kpoint%atype = atype
175 :
176 162 : CALL crys_sym_gen(crys_sym, scoord, atype, cell%hmat, delta=kpoint%eps_geo, iounit=iounit)
177 : CALL kpoint_gen(crys_sym, kpoint%nkp_grid, symm=kpoint%symmetry, shift=kpoint%kp_shift, &
178 162 : full_grid=kpoint%full_grid)
179 162 : kpoint%nkp = crys_sym%nkpoint
180 810 : ALLOCATE (kpoint%xkp(3, kpoint%nkp), kpoint%wkp(kpoint%nkp))
181 1770 : wsum = SUM(crys_sym%wkpoint)
182 1770 : DO ik = 1, kpoint%nkp
183 6432 : kpoint%xkp(1:3, ik) = crys_sym%xkpoint(1:3, ik)
184 1770 : kpoint%wkp(ik) = crys_sym%wkpoint(ik)/wsum
185 : END DO
186 :
187 : ! print output
188 162 : IF (kpoint%symmetry) CALL print_crys_symmetry(crys_sym)
189 162 : IF (kpoint%symmetry) CALL print_kp_symmetry(crys_sym)
190 :
191 : ! transfer symmetry information
192 2094 : ALLOCATE (kpoint%kp_sym(kpoint%nkp))
193 1770 : DO ik = 1, kpoint%nkp
194 1608 : NULLIFY (kpoint%kp_sym(ik)%kpoint_sym)
195 1608 : CALL kpoint_sym_create(kpoint%kp_sym(ik)%kpoint_sym)
196 1608 : kpsym => kpoint%kp_sym(ik)%kpoint_sym
197 1770 : IF (crys_sym%symlib .AND. .NOT. crys_sym%fullgrid .AND. crys_sym%istriz == 1) THEN
198 : ! set up the symmetrization information
199 0 : kpsym%nwght = NINT(crys_sym%wkpoint(ik))
200 0 : ns = kpsym%nwght
201 : !
202 0 : IF (ns > 1) THEN
203 0 : kpsym%apply_symmetry = .TRUE.
204 0 : natom = SIZE(particle_set)
205 0 : ALLOCATE (kpsym%rot(3, 3, ns))
206 0 : ALLOCATE (kpsym%xkp(3, ns))
207 0 : ALLOCATE (kpsym%rotp(ns))
208 0 : ALLOCATE (kpsym%f0(natom, ns))
209 0 : nr = 0
210 0 : DO is = 1, SIZE(crys_sym%kplink, 2)
211 0 : IF (crys_sym%kplink(2, is) == ik) THEN
212 0 : nr = nr + 1
213 0 : ir = crys_sym%kpop(is)
214 0 : ira = ABS(ir)
215 0 : kpsym%rotp(nr) = ir
216 0 : IF (ir > 0) THEN
217 0 : kpsym%rot(1:3, 1:3, nr) = crys_sym%rt(1:3, 1:3, ira)
218 : ELSE
219 0 : kpsym%rot(1:3, 1:3, nr) = -crys_sym%rt(1:3, 1:3, ira)
220 : END IF
221 0 : kpsym%xkp(1:3, nr) = crys_sym%kpmesh(1:3, is)
222 0 : DO ic = 1, crys_sym%nrtot
223 0 : IF (crys_sym%ibrot(ic) == ira) THEN
224 0 : kpsym%f0(1:natom, nr) = crys_sym%f0(1:natom, ic)
225 : EXIT
226 : END IF
227 : END DO
228 : END IF
229 : END DO
230 : ! Reduce inversion?
231 0 : kpsym%nwred = nr
232 : END IF
233 : END IF
234 : END DO
235 162 : IF (kpoint%symmetry) THEN
236 368 : nkind = MAXVAL(atype)
237 50 : ns = crys_sym%nrtot
238 208 : ALLOCATE (kpoint%kind_rotmat(ns, nkind))
239 50 : DO i = 1, ns
240 50 : DO j = 1, nkind
241 0 : NULLIFY (kpoint%kind_rotmat(i, j)%rmat)
242 : END DO
243 : END DO
244 100 : ALLOCATE (kpoint%ibrot(ns))
245 50 : kpoint%ibrot(1:ns) = crys_sym%ibrot(1:ns)
246 : END IF
247 :
248 162 : CALL release_csym_type(crys_sym)
249 162 : DEALLOCATE (scoord, atype)
250 :
251 : CASE ("GENERAL")
252 : ! default: no symmetry settings
253 20 : ALLOCATE (kpoint%kp_sym(kpoint%nkp))
254 12 : DO i = 1, kpoint%nkp
255 8 : NULLIFY (kpoint%kp_sym(i)%kpoint_sym)
256 12 : CALL kpoint_sym_create(kpoint%kp_sym(i)%kpoint_sym)
257 : END DO
258 : CASE DEFAULT
259 7682 : CPASSERT(.FALSE.)
260 : END SELECT
261 :
262 : ! check for consistency of options
263 7698 : SELECT CASE (kpoint%kp_scheme)
264 : CASE ("NONE")
265 : ! don't use k-point code
266 : CASE ("GAMMA")
267 16 : CPASSERT(kpoint%nkp == 1)
268 80 : CPASSERT(SUM(ABS(kpoint%xkp)) <= 1.e-12_dp)
269 16 : CPASSERT(kpoint%wkp(1) == 1.0_dp)
270 16 : CPASSERT(.NOT. kpoint%symmetry)
271 : CASE ("GENERAL")
272 4 : CPASSERT(.NOT. kpoint%symmetry)
273 4 : CPASSERT(kpoint%nkp >= 1)
274 : CASE ("MONKHORST-PACK", "MACDONALD")
275 7682 : CPASSERT(kpoint%nkp >= 1)
276 : END SELECT
277 7682 : IF (kpoint%use_real_wfn) THEN
278 : ! what about inversion symmetry?
279 24 : ikloop: DO ik = 1, kpoint%nkp
280 60 : DO i = 1, 3
281 36 : spez = (kpoint%xkp(i, ik) == 0.0_dp .OR. kpoint%xkp(i, ik) == 0.5_dp)
282 12 : IF (.NOT. spez) EXIT ikloop
283 : END DO
284 : END DO ikloop
285 12 : IF (.NOT. spez) THEN
286 : ! Warning: real wfn might be wrong for this system
287 : CALL cp_warn(__LOCATION__, &
288 : "A calculation using real wavefunctions is requested. "// &
289 0 : "We could not determine if the symmetry of the system allows real wavefunctions. ")
290 : END IF
291 : END IF
292 :
293 7682 : CALL timestop(handle)
294 :
295 7682 : END SUBROUTINE kpoint_initialize
296 :
297 : ! **************************************************************************************************
298 : !> \brief Initialize the kpoint environment
299 : !> \param kpoint Kpoint environment
300 : !> \param para_env ...
301 : !> \param blacs_env ...
302 : !> \param with_aux_fit ...
303 : ! **************************************************************************************************
304 238 : SUBROUTINE kpoint_env_initialize(kpoint, para_env, blacs_env, with_aux_fit)
305 :
306 : TYPE(kpoint_type), INTENT(INOUT) :: kpoint
307 : TYPE(mp_para_env_type), INTENT(IN), TARGET :: para_env
308 : TYPE(cp_blacs_env_type), INTENT(IN), TARGET :: blacs_env
309 : LOGICAL, INTENT(IN), OPTIONAL :: with_aux_fit
310 :
311 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_env_initialize'
312 :
313 : INTEGER :: handle, igr, ik, ikk, ngr, niogrp, nkp, &
314 : nkp_grp, nkp_loc, npe, unit_nr
315 : INTEGER, DIMENSION(2) :: dims, pos
316 : LOGICAL :: aux_fit
317 238 : TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_aux_env, kp_env
318 : TYPE(kpoint_env_type), POINTER :: kp
319 238 : TYPE(mp_cart_type) :: comm_cart
320 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp, para_env_kp
321 :
322 238 : CALL timeset(routineN, handle)
323 :
324 238 : IF (PRESENT(with_aux_fit)) THEN
325 180 : aux_fit = with_aux_fit
326 : ELSE
327 : aux_fit = .FALSE.
328 : END IF
329 :
330 238 : kpoint%para_env => para_env
331 238 : CALL kpoint%para_env%retain()
332 238 : kpoint%blacs_env_all => blacs_env
333 238 : CALL kpoint%blacs_env_all%retain()
334 :
335 238 : CPASSERT(.NOT. ASSOCIATED(kpoint%kp_env))
336 238 : IF (aux_fit) THEN
337 30 : CPASSERT(.NOT. ASSOCIATED(kpoint%kp_aux_env))
338 : END IF
339 :
340 238 : NULLIFY (kp_env, kp_aux_env)
341 238 : nkp = kpoint%nkp
342 238 : npe = para_env%num_pe
343 238 : IF (npe == 1) THEN
344 : ! only one process available -> owns all kpoints
345 0 : ALLOCATE (kp_env(nkp))
346 0 : DO ik = 1, nkp
347 0 : NULLIFY (kp_env(ik)%kpoint_env)
348 0 : CALL kpoint_env_create(kp_env(ik)%kpoint_env)
349 0 : kp => kp_env(ik)%kpoint_env
350 0 : kp%nkpoint = ik
351 0 : kp%wkp = kpoint%wkp(ik)
352 0 : kp%xkp(1:3) = kpoint%xkp(1:3, ik)
353 0 : kp%is_local = .TRUE.
354 : END DO
355 0 : kpoint%kp_env => kp_env
356 :
357 0 : IF (aux_fit) THEN
358 0 : ALLOCATE (kp_aux_env(nkp))
359 0 : DO ik = 1, nkp
360 0 : NULLIFY (kp_aux_env(ik)%kpoint_env)
361 0 : CALL kpoint_env_create(kp_aux_env(ik)%kpoint_env)
362 0 : kp => kp_aux_env(ik)%kpoint_env
363 0 : kp%nkpoint = ik
364 0 : kp%wkp = kpoint%wkp(ik)
365 0 : kp%xkp(1:3) = kpoint%xkp(1:3, ik)
366 0 : kp%is_local = .TRUE.
367 : END DO
368 :
369 0 : kpoint%kp_aux_env => kp_aux_env
370 : END IF
371 :
372 0 : ALLOCATE (kpoint%kp_dist(2, 1))
373 0 : kpoint%kp_dist(1, 1) = 1
374 0 : kpoint%kp_dist(2, 1) = nkp
375 0 : kpoint%kp_range(1) = 1
376 0 : kpoint%kp_range(2) = nkp
377 :
378 : ! parallel environments
379 0 : kpoint%para_env_kp => para_env
380 0 : CALL kpoint%para_env_kp%retain()
381 0 : kpoint%para_env_inter_kp => para_env
382 0 : CALL kpoint%para_env_inter_kp%retain()
383 0 : kpoint%iogrp = .TRUE.
384 0 : kpoint%nkp_groups = 1
385 : ELSE
386 238 : IF (kpoint%parallel_group_size == -1) THEN
387 : ! maximum parallelization over kpoints
388 : ! making sure that the group size divides the npe and the nkp_grp the nkp
389 : ! in the worst case, there will be no parallelism over kpoints.
390 378 : DO igr = npe, 1, -1
391 252 : IF (MOD(npe, igr) /= 0) CYCLE
392 252 : nkp_grp = npe/igr
393 252 : IF (MOD(nkp, nkp_grp) /= 0) CYCLE
394 378 : ngr = igr
395 : END DO
396 112 : ELSE IF (kpoint%parallel_group_size == 0) THEN
397 : ! no parallelization over kpoints
398 64 : ngr = npe
399 48 : ELSE IF (kpoint%parallel_group_size > 0) THEN
400 48 : ngr = MIN(kpoint%parallel_group_size, npe)
401 : ELSE
402 0 : CPASSERT(.FALSE.)
403 : END IF
404 238 : nkp_grp = npe/ngr
405 : ! processor dimensions
406 238 : dims(1) = ngr
407 238 : dims(2) = nkp_grp
408 238 : CPASSERT(MOD(nkp, nkp_grp) == 0)
409 238 : nkp_loc = nkp/nkp_grp
410 :
411 238 : IF ((dims(1)*dims(2) /= npe)) THEN
412 0 : CPABORT("Number of processors is not divisible by the kpoint group size.")
413 : END IF
414 :
415 : ! Create the subgroups, one for each k-point group and one interconnecting group
416 238 : CALL comm_cart%create(comm_old=para_env, ndims=2, dims=dims)
417 714 : pos = comm_cart%mepos_cart
418 238 : ALLOCATE (para_env_kp)
419 238 : CALL para_env_kp%from_split(comm_cart, pos(2))
420 238 : ALLOCATE (para_env_inter_kp)
421 238 : CALL para_env_inter_kp%from_split(comm_cart, pos(1))
422 238 : CALL comm_cart%free()
423 :
424 238 : niogrp = 0
425 238 : IF (para_env%is_source()) niogrp = 1
426 238 : CALL para_env_kp%sum(niogrp)
427 238 : kpoint%iogrp = (niogrp == 1)
428 :
429 : ! parallel groups
430 238 : kpoint%para_env_kp => para_env_kp
431 238 : kpoint%para_env_inter_kp => para_env_inter_kp
432 :
433 : ! distribution of kpoints
434 714 : ALLOCATE (kpoint%kp_dist(2, nkp_grp))
435 562 : DO igr = 1, nkp_grp
436 1210 : kpoint%kp_dist(1:2, igr) = get_limit(nkp, nkp_grp, igr - 1)
437 : END DO
438 : ! local kpoints
439 714 : kpoint%kp_range(1:2) = kpoint%kp_dist(1:2, para_env_inter_kp%mepos + 1)
440 :
441 2262 : ALLOCATE (kp_env(nkp_loc))
442 1786 : DO ik = 1, nkp_loc
443 1548 : NULLIFY (kp_env(ik)%kpoint_env)
444 1548 : ikk = kpoint%kp_range(1) + ik - 1
445 1548 : CALL kpoint_env_create(kp_env(ik)%kpoint_env)
446 1548 : kp => kp_env(ik)%kpoint_env
447 1548 : kp%nkpoint = ikk
448 1548 : kp%wkp = kpoint%wkp(ikk)
449 6192 : kp%xkp(1:3) = kpoint%xkp(1:3, ikk)
450 1786 : kp%is_local = (ngr == 1)
451 : END DO
452 238 : kpoint%kp_env => kp_env
453 :
454 238 : IF (aux_fit) THEN
455 274 : ALLOCATE (kp_aux_env(nkp_loc))
456 244 : DO ik = 1, nkp_loc
457 214 : NULLIFY (kp_aux_env(ik)%kpoint_env)
458 214 : ikk = kpoint%kp_range(1) + ik - 1
459 214 : CALL kpoint_env_create(kp_aux_env(ik)%kpoint_env)
460 214 : kp => kp_aux_env(ik)%kpoint_env
461 214 : kp%nkpoint = ikk
462 214 : kp%wkp = kpoint%wkp(ikk)
463 856 : kp%xkp(1:3) = kpoint%xkp(1:3, ikk)
464 244 : kp%is_local = (ngr == 1)
465 : END DO
466 30 : kpoint%kp_aux_env => kp_aux_env
467 : END IF
468 :
469 238 : unit_nr = cp_logger_get_default_io_unit()
470 :
471 238 : IF (unit_nr > 0 .AND. kpoint%verbose) THEN
472 5 : WRITE (unit_nr, *)
473 5 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Number of kpoint groups ", nkp_grp
474 5 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Size of each kpoint group", ngr
475 5 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Number of kpoints per group", nkp_loc
476 : END IF
477 238 : kpoint%nkp_groups = nkp_grp
478 :
479 : END IF
480 :
481 238 : CALL timestop(handle)
482 :
483 476 : END SUBROUTINE kpoint_env_initialize
484 :
485 : ! **************************************************************************************************
486 : !> \brief Initialize a set of MOs and density matrix for each kpoint (kpoint group)
487 : !> \param kpoint Kpoint environment
488 : !> \param mos Reference MOs (global)
489 : !> \param added_mos ...
490 : !> \param for_aux_fit ...
491 : ! **************************************************************************************************
492 268 : SUBROUTINE kpoint_initialize_mos(kpoint, mos, added_mos, for_aux_fit)
493 :
494 : TYPE(kpoint_type), POINTER :: kpoint
495 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
496 : INTEGER, INTENT(IN), OPTIONAL :: added_mos
497 : LOGICAL, OPTIONAL :: for_aux_fit
498 :
499 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_initialize_mos'
500 :
501 : INTEGER :: handle, ic, ik, is, nadd, nao, nc, &
502 : nelectron, nkp_loc, nmo, nmorig(2), &
503 : nspin
504 : LOGICAL :: aux_fit
505 : REAL(KIND=dp) :: flexible_electron_count, maxocc, n_el_f
506 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
507 268 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools
508 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
509 : TYPE(cp_fm_type), POINTER :: fmlocal
510 : TYPE(kpoint_env_type), POINTER :: kp
511 : TYPE(qs_matrix_pools_type), POINTER :: mpools
512 :
513 268 : CALL timeset(routineN, handle)
514 :
515 268 : IF (PRESENT(for_aux_fit)) THEN
516 30 : aux_fit = for_aux_fit
517 : ELSE
518 : aux_fit = .FALSE.
519 : END IF
520 :
521 268 : CPASSERT(ASSOCIATED(kpoint))
522 :
523 : IF (.TRUE. .OR. ASSOCIATED(mos(1)%mo_coeff)) THEN
524 268 : IF (aux_fit) THEN
525 30 : CPASSERT(ASSOCIATED(kpoint%kp_aux_env))
526 : END IF
527 :
528 268 : IF (PRESENT(added_mos)) THEN
529 58 : nadd = added_mos
530 : ELSE
531 : nadd = 0
532 : END IF
533 :
534 268 : IF (kpoint%use_real_wfn) THEN
535 : nc = 1
536 : ELSE
537 256 : nc = 2
538 : END IF
539 268 : nspin = SIZE(mos, 1)
540 268 : nkp_loc = kpoint%kp_range(2) - kpoint%kp_range(1) + 1
541 268 : IF (nkp_loc > 0) THEN
542 268 : IF (aux_fit) THEN
543 30 : CPASSERT(SIZE(kpoint%kp_aux_env) == nkp_loc)
544 : ELSE
545 238 : CPASSERT(SIZE(kpoint%kp_env) == nkp_loc)
546 : END IF
547 : ! allocate the mo sets, correct number of kpoints (local), real/complex, spin
548 2030 : DO ik = 1, nkp_loc
549 1762 : IF (aux_fit) THEN
550 214 : kp => kpoint%kp_aux_env(ik)%kpoint_env
551 : ELSE
552 1548 : kp => kpoint%kp_env(ik)%kpoint_env
553 : END IF
554 13112 : ALLOCATE (kp%mos(nc, nspin))
555 4056 : DO is = 1, nspin
556 : CALL get_mo_set(mos(is), nao=nao, nmo=nmo, nelectron=nelectron, &
557 2026 : n_el_f=n_el_f, maxocc=maxocc, flexible_electron_count=flexible_electron_count)
558 2026 : nmo = MIN(nao, nmo + nadd)
559 7826 : DO ic = 1, nc
560 : CALL allocate_mo_set(kp%mos(ic, is), nao, nmo, nelectron, n_el_f, maxocc, &
561 6064 : flexible_electron_count)
562 : END DO
563 : END DO
564 : END DO
565 :
566 : ! generate the blacs environment for the kpoint group
567 : ! we generate a blacs env for each kpoint group in parallel
568 : ! we assume here that the group para_env_inter_kp will connect
569 : ! equivalent parts of fm matrices, i.e. no reshuffeling of processors
570 268 : NULLIFY (blacs_env)
571 268 : IF (ASSOCIATED(kpoint%blacs_env)) THEN
572 30 : blacs_env => kpoint%blacs_env
573 : ELSE
574 238 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=kpoint%para_env_kp)
575 238 : kpoint%blacs_env => blacs_env
576 : END IF
577 :
578 : ! set possible new number of MOs
579 574 : DO is = 1, nspin
580 306 : CALL get_mo_set(mos(is), nmo=nmorig(is))
581 306 : nmo = MIN(nao, nmorig(is) + nadd)
582 574 : CALL set_mo_set(mos(is), nmo=nmo)
583 : END DO
584 : ! matrix pools for the kpoint group, information on MOs is transferred using
585 : ! generic mos structure
586 268 : NULLIFY (mpools)
587 268 : CALL mpools_create(mpools=mpools)
588 : CALL mpools_rebuild_fm_pools(mpools=mpools, mos=mos, &
589 268 : blacs_env=blacs_env, para_env=kpoint%para_env_kp)
590 :
591 268 : IF (aux_fit) THEN
592 30 : kpoint%mpools_aux_fit => mpools
593 : ELSE
594 238 : kpoint%mpools => mpools
595 : END IF
596 :
597 : ! reset old number of MOs
598 574 : DO is = 1, nspin
599 574 : CALL set_mo_set(mos(is), nmo=nmorig(is))
600 : END DO
601 :
602 : ! allocate density matrices
603 268 : CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
604 268 : ALLOCATE (fmlocal)
605 268 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
606 268 : CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
607 2030 : DO ik = 1, nkp_loc
608 1762 : IF (aux_fit) THEN
609 214 : kp => kpoint%kp_aux_env(ik)%kpoint_env
610 : ELSE
611 1548 : kp => kpoint%kp_env(ik)%kpoint_env
612 : END IF
613 : ! density matrix
614 1762 : CALL cp_fm_release(kp%pmat)
615 13112 : ALLOCATE (kp%pmat(nc, nspin))
616 3788 : DO is = 1, nspin
617 7826 : DO ic = 1, nc
618 6064 : CALL cp_fm_create(kp%pmat(ic, is), matrix_struct)
619 : END DO
620 : END DO
621 : ! energy weighted density matrix
622 1762 : CALL cp_fm_release(kp%wmat)
623 11350 : ALLOCATE (kp%wmat(nc, nspin))
624 4056 : DO is = 1, nspin
625 7826 : DO ic = 1, nc
626 6064 : CALL cp_fm_create(kp%wmat(ic, is), matrix_struct)
627 : END DO
628 : END DO
629 : END DO
630 268 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
631 268 : DEALLOCATE (fmlocal)
632 :
633 : END IF
634 :
635 : END IF
636 :
637 268 : CALL timestop(handle)
638 :
639 268 : END SUBROUTINE kpoint_initialize_mos
640 :
641 : ! **************************************************************************************************
642 : !> \brief ...
643 : !> \param kpoint ...
644 : ! **************************************************************************************************
645 58 : SUBROUTINE kpoint_initialize_mo_set(kpoint)
646 : TYPE(kpoint_type), POINTER :: kpoint
647 :
648 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_initialize_mo_set'
649 :
650 : INTEGER :: handle, ic, ik, ikk, ispin
651 58 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
652 : TYPE(cp_fm_type), POINTER :: mo_coeff
653 58 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: moskp
654 :
655 58 : CALL timeset(routineN, handle)
656 :
657 404 : DO ik = 1, SIZE(kpoint%kp_env)
658 346 : CALL mpools_get(kpoint%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
659 346 : moskp => kpoint%kp_env(ik)%kpoint_env%mos
660 346 : ikk = kpoint%kp_range(1) + ik - 1
661 346 : CPASSERT(ASSOCIATED(moskp))
662 784 : DO ispin = 1, SIZE(moskp, 2)
663 1486 : DO ic = 1, SIZE(moskp, 1)
664 760 : CALL get_mo_set(moskp(ic, ispin), mo_coeff=mo_coeff)
665 1140 : IF (.NOT. ASSOCIATED(mo_coeff)) THEN
666 : CALL init_mo_set(moskp(ic, ispin), &
667 760 : fm_pool=ao_mo_fm_pools(ispin)%pool, name="kpoints")
668 : END IF
669 : END DO
670 : END DO
671 : END DO
672 :
673 58 : CALL timestop(handle)
674 :
675 58 : END SUBROUTINE kpoint_initialize_mo_set
676 :
677 : ! **************************************************************************************************
678 : !> \brief Generates the mapping of cell indices and linear RS index
679 : !> CELL (0,0,0) is always mapped to index 1
680 : !> \param kpoint Kpoint environment
681 : !> \param sab_nl Defining neighbour list
682 : !> \param para_env Parallel environment
683 : !> \param dft_control ...
684 : ! **************************************************************************************************
685 1032 : SUBROUTINE kpoint_init_cell_index(kpoint, sab_nl, para_env, dft_control)
686 :
687 : TYPE(kpoint_type), POINTER :: kpoint
688 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
689 : POINTER :: sab_nl
690 : TYPE(mp_para_env_type), POINTER :: para_env
691 : TYPE(dft_control_type), POINTER :: dft_control
692 :
693 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_init_cell_index'
694 :
695 : INTEGER :: handle, i1, i2, i3, ic, icount, it, &
696 : ncount
697 : INTEGER, DIMENSION(3) :: cell, itm
698 1032 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell, list
699 1032 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index, cti
700 : LOGICAL :: new
701 : TYPE(neighbor_list_iterator_p_type), &
702 1032 : DIMENSION(:), POINTER :: nl_iterator
703 :
704 1032 : NULLIFY (cell_to_index, index_to_cell)
705 :
706 1032 : CALL timeset(routineN, handle)
707 :
708 1032 : CPASSERT(ASSOCIATED(kpoint))
709 :
710 1032 : ALLOCATE (list(3, 125))
711 517032 : list = 0
712 1032 : icount = 1
713 :
714 1032 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
715 614395 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
716 613363 : CALL get_iterator_info(nl_iterator, cell=cell)
717 :
718 613363 : new = .TRUE.
719 36116490 : DO ic = 1, icount
720 36032898 : IF (cell(1) == list(1, ic) .AND. cell(2) == list(2, ic) .AND. &
721 83592 : cell(3) == list(3, ic)) THEN
722 : new = .FALSE.
723 : EXIT
724 : END IF
725 : END DO
726 614395 : IF (new) THEN
727 83592 : icount = icount + 1
728 83592 : IF (icount > SIZE(list, 2)) THEN
729 273 : CALL reallocate(list, 1, 3, 1, 2*SIZE(list, 2))
730 : END IF
731 334368 : list(1:3, icount) = cell(1:3)
732 : END IF
733 :
734 : END DO
735 1032 : CALL neighbor_list_iterator_release(nl_iterator)
736 :
737 85656 : itm(1) = MAXVAL(ABS(list(1, 1:icount)))
738 85656 : itm(2) = MAXVAL(ABS(list(2, 1:icount)))
739 85656 : itm(3) = MAXVAL(ABS(list(3, 1:icount)))
740 1032 : CALL para_env%max(itm)
741 4128 : it = MAXVAL(itm(1:3))
742 1032 : IF (ASSOCIATED(kpoint%cell_to_index)) THEN
743 1028 : DEALLOCATE (kpoint%cell_to_index)
744 : END IF
745 5160 : ALLOCATE (kpoint%cell_to_index(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
746 1032 : cell_to_index => kpoint%cell_to_index
747 1032 : cti => cell_to_index
748 205320 : cti(:, :, :) = 0
749 85656 : DO ic = 1, icount
750 84624 : i1 = list(1, ic)
751 84624 : i2 = list(2, ic)
752 84624 : i3 = list(3, ic)
753 85656 : cti(i1, i2, i3) = ic
754 : END DO
755 409608 : CALL para_env%sum(cti)
756 1032 : ncount = 0
757 6148 : DO i1 = -itm(1), itm(1)
758 38876 : DO i2 = -itm(2), itm(2)
759 209112 : DO i3 = -itm(3), itm(3)
760 203996 : IF (cti(i1, i2, i3) == 0) THEN
761 77714 : cti(i1, i2, i3) = 1000000
762 : ELSE
763 93554 : ncount = ncount + 1
764 93554 : cti(i1, i2, i3) = (ABS(i1) + ABS(i2) + ABS(i3))*1000 + ABS(i3)*100 + ABS(i2)*10 + ABS(i1)
765 93554 : cti(i1, i2, i3) = cti(i1, i2, i3) + (i1 + i2 + i3)
766 : END IF
767 : END DO
768 : END DO
769 : END DO
770 :
771 1032 : IF (ASSOCIATED(kpoint%index_to_cell)) THEN
772 1032 : DEALLOCATE (kpoint%index_to_cell)
773 : END IF
774 3096 : ALLOCATE (kpoint%index_to_cell(3, ncount))
775 1032 : index_to_cell => kpoint%index_to_cell
776 94586 : DO ic = 1, ncount
777 93554 : cell = MINLOC(cti)
778 93554 : i1 = cell(1) - 1 - itm(1)
779 93554 : i2 = cell(2) - 1 - itm(2)
780 93554 : i3 = cell(3) - 1 - itm(3)
781 93554 : cti(i1, i2, i3) = 1000000
782 93554 : index_to_cell(1, ic) = i1
783 93554 : index_to_cell(2, ic) = i2
784 94586 : index_to_cell(3, ic) = i3
785 : END DO
786 205320 : cti(:, :, :) = 0
787 94586 : DO ic = 1, ncount
788 93554 : i1 = index_to_cell(1, ic)
789 93554 : i2 = index_to_cell(2, ic)
790 93554 : i3 = index_to_cell(3, ic)
791 94586 : cti(i1, i2, i3) = ic
792 : END DO
793 :
794 : ! keep pointer to this neighborlist
795 1032 : kpoint%sab_nl => sab_nl
796 :
797 : ! set number of images
798 1032 : dft_control%nimages = SIZE(index_to_cell, 2)
799 :
800 1032 : DEALLOCATE (list)
801 :
802 1032 : CALL timestop(handle)
803 1032 : END SUBROUTINE kpoint_init_cell_index
804 :
805 : ! **************************************************************************************************
806 : !> \brief Transformation of real space matrices to a kpoint
807 : !> \param rmatrix Real part of kpoint matrix
808 : !> \param cmatrix Complex part of kpoint matrix (optional)
809 : !> \param rsmat Real space matrices
810 : !> \param ispin Spin index
811 : !> \param xkp Kpoint coordinates
812 : !> \param cell_to_index mapping of cell indices to RS index
813 : !> \param sab_nl Defining neighbor list
814 : !> \param is_complex Matrix to be transformed is imaginary
815 : !> \param rs_sign Matrix to be transformed is csaled by rs_sign
816 : ! **************************************************************************************************
817 112188 : SUBROUTINE rskp_transform(rmatrix, cmatrix, rsmat, ispin, &
818 : xkp, cell_to_index, sab_nl, is_complex, rs_sign)
819 :
820 : TYPE(dbcsr_type) :: rmatrix
821 : TYPE(dbcsr_type), OPTIONAL :: cmatrix
822 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
823 : INTEGER, INTENT(IN) :: ispin
824 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
825 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
826 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
827 : POINTER :: sab_nl
828 : LOGICAL, INTENT(IN), OPTIONAL :: is_complex
829 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: rs_sign
830 :
831 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rskp_transform'
832 :
833 : INTEGER :: handle, iatom, ic, icol, irow, jatom, &
834 : nimg
835 : INTEGER, DIMENSION(3) :: cell
836 : LOGICAL :: do_symmetric, found, my_complex, &
837 : wfn_real_only
838 : REAL(KIND=dp) :: arg, coskl, fsign, fsym, sinkl
839 56094 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, rblock, rsblock
840 : TYPE(neighbor_list_iterator_p_type), &
841 56094 : DIMENSION(:), POINTER :: nl_iterator
842 :
843 56094 : CALL timeset(routineN, handle)
844 :
845 56094 : my_complex = .FALSE.
846 56094 : IF (PRESENT(is_complex)) my_complex = is_complex
847 :
848 56094 : fsign = 1.0_dp
849 56094 : IF (PRESENT(rs_sign)) fsign = rs_sign
850 :
851 56094 : wfn_real_only = .TRUE.
852 56094 : IF (PRESENT(cmatrix)) wfn_real_only = .FALSE.
853 :
854 56094 : nimg = SIZE(rsmat, 2)
855 :
856 56094 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
857 :
858 56094 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
859 22574814 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
860 22518720 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
861 :
862 : ! fsym = +- 1 is due to real space matrices being non-symmetric (although in a symmtric type)
863 : ! with the link S_mu^0,nu^b = S_nu^0,mu^-b, and the KP matrices beeing Hermitian
864 22518720 : fsym = 1.0_dp
865 22518720 : irow = iatom
866 22518720 : icol = jatom
867 22518720 : IF (do_symmetric .AND. (iatom > jatom)) THEN
868 9490958 : irow = jatom
869 9490958 : icol = iatom
870 9490958 : fsym = -1.0_dp
871 : END IF
872 :
873 22518720 : ic = cell_to_index(cell(1), cell(2), cell(3))
874 22518720 : IF (ic < 1 .OR. ic > nimg) CYCLE
875 :
876 22517976 : arg = REAL(cell(1), dp)*xkp(1) + REAL(cell(2), dp)*xkp(2) + REAL(cell(3), dp)*xkp(3)
877 22517976 : IF (my_complex) THEN
878 0 : coskl = fsign*fsym*COS(twopi*arg)
879 0 : sinkl = fsign*SIN(twopi*arg)
880 : ELSE
881 22517976 : coskl = fsign*COS(twopi*arg)
882 22517976 : sinkl = fsign*fsym*SIN(twopi*arg)
883 : END IF
884 :
885 : CALL dbcsr_get_block_p(matrix=rsmat(ispin, ic)%matrix, row=irow, col=icol, &
886 22517976 : block=rsblock, found=found)
887 22517976 : IF (.NOT. found) CYCLE
888 :
889 22574070 : IF (wfn_real_only) THEN
890 : CALL dbcsr_get_block_p(matrix=rmatrix, row=irow, col=icol, &
891 303864 : block=rblock, found=found)
892 303864 : IF (.NOT. found) CYCLE
893 153367320 : rblock = rblock + coskl*rsblock
894 : ELSE
895 : CALL dbcsr_get_block_p(matrix=rmatrix, row=irow, col=icol, &
896 22214112 : block=rblock, found=found)
897 22214112 : IF (.NOT. found) CYCLE
898 : CALL dbcsr_get_block_p(matrix=cmatrix, row=irow, col=icol, &
899 22214112 : block=cblock, found=found)
900 22214112 : IF (.NOT. found) CYCLE
901 4326617316 : rblock = rblock + coskl*rsblock
902 4326617316 : cblock = cblock + sinkl*rsblock
903 : END IF
904 :
905 : END DO
906 56094 : CALL neighbor_list_iterator_release(nl_iterator)
907 :
908 56094 : CALL timestop(handle)
909 :
910 56094 : END SUBROUTINE rskp_transform
911 :
912 : ! **************************************************************************************************
913 : !> \brief Given the eigenvalues of all kpoints, calculates the occupation numbers
914 : !> \param kpoint Kpoint environment
915 : !> \param smear Smearing information
916 : !> \param probe ...
917 : ! **************************************************************************************************
918 5666 : SUBROUTINE kpoint_set_mo_occupation(kpoint, smear, probe)
919 :
920 : TYPE(kpoint_type), POINTER :: kpoint
921 : TYPE(smear_type), POINTER :: smear
922 : TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
923 : POINTER :: probe
924 :
925 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_set_mo_occupation'
926 :
927 : INTEGER :: handle, ik, ikpgr, ispin, kplocal, nao, &
928 : nb, ncol_global, ne_a, ne_b, &
929 : nelectron, nkp, nmo, nrow_global, nspin
930 : INTEGER, DIMENSION(2) :: kp_range
931 : REAL(KIND=dp) :: kTS, mu, mus(2), nel
932 5666 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: smatrix
933 5666 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: weig, wocc
934 5666 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: icoeff, rcoeff
935 5666 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation, wkp
936 : TYPE(cp_fm_type), POINTER :: mo_coeff
937 : TYPE(kpoint_env_type), POINTER :: kp
938 : TYPE(mo_set_type), POINTER :: mo_set
939 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
940 :
941 5666 : CALL timeset(routineN, handle)
942 :
943 : ! first collect all the eigenvalues
944 5666 : CALL get_kpoint_info(kpoint, nkp=nkp)
945 5666 : kp => kpoint%kp_env(1)%kpoint_env
946 5666 : nspin = SIZE(kp%mos, 2)
947 5666 : mo_set => kp%mos(1, 1)
948 5666 : CALL get_mo_set(mo_set, nmo=nmo, nao=nao, nelectron=nelectron)
949 5666 : ne_a = nelectron
950 5666 : IF (nspin == 2) THEN
951 548 : CALL get_mo_set(kp%mos(1, 2), nmo=nb, nelectron=ne_b)
952 548 : CPASSERT(nmo == nb)
953 : END IF
954 45328 : ALLOCATE (weig(nmo, nkp, nspin), wocc(nmo, nkp, nspin))
955 476642 : weig = 0.0_dp
956 476642 : wocc = 0.0_dp
957 5666 : IF (PRESENT(probe)) THEN
958 0 : ALLOCATE (rcoeff(nao, nmo, nkp, nspin), icoeff(nao, nmo, nkp, nspin))
959 0 : rcoeff = 0.0_dp !coeff, real part
960 0 : icoeff = 0.0_dp !coeff, imaginary part
961 : END IF
962 5666 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
963 5666 : kplocal = kp_range(2) - kp_range(1) + 1
964 21036 : DO ikpgr = 1, kplocal
965 15370 : ik = kp_range(1) + ikpgr - 1
966 15370 : kp => kpoint%kp_env(ikpgr)%kpoint_env
967 38298 : DO ispin = 1, nspin
968 17262 : mo_set => kp%mos(1, ispin)
969 17262 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
970 338600 : weig(1:nmo, ik, ispin) = eigenvalues(1:nmo)
971 32632 : IF (PRESENT(probe)) THEN
972 0 : CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
973 : CALL cp_fm_get_info(mo_coeff, &
974 : nrow_global=nrow_global, &
975 0 : ncol_global=ncol_global)
976 0 : ALLOCATE (smatrix(nrow_global, ncol_global))
977 0 : CALL cp_fm_get_submatrix(mo_coeff, smatrix)
978 :
979 0 : rcoeff(1:nao, 1:nmo, ik, ispin) = smatrix(1:nrow_global, 1:ncol_global)
980 :
981 0 : DEALLOCATE (smatrix)
982 :
983 0 : mo_set => kp%mos(2, ispin)
984 :
985 0 : CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
986 : CALL cp_fm_get_info(mo_coeff, &
987 : nrow_global=nrow_global, &
988 0 : ncol_global=ncol_global)
989 0 : ALLOCATE (smatrix(nrow_global, ncol_global))
990 0 : CALL cp_fm_get_submatrix(mo_coeff, smatrix)
991 :
992 0 : icoeff(1:nao, 1:nmo, ik, ispin) = smatrix(1:nrow_global, 1:ncol_global)
993 :
994 0 : mo_set => kp%mos(1, ispin)
995 :
996 0 : DEALLOCATE (smatrix)
997 : END IF
998 : END DO
999 : END DO
1000 5666 : CALL get_kpoint_info(kpoint, para_env_inter_kp=para_env_inter_kp)
1001 5666 : CALL para_env_inter_kp%sum(weig)
1002 :
1003 5666 : IF (PRESENT(probe)) THEN
1004 0 : CALL para_env_inter_kp%sum(rcoeff)
1005 0 : CALL para_env_inter_kp%sum(icoeff)
1006 : END IF
1007 :
1008 5666 : CALL get_kpoint_info(kpoint, wkp=wkp)
1009 :
1010 : !calling of HP module HERE, before smear
1011 5666 : IF (PRESENT(probe)) THEN
1012 0 : smear%do_smear = .FALSE. !ensures smearing is switched off
1013 :
1014 0 : IF (nspin == 1) THEN
1015 0 : nel = REAL(nelectron, KIND=dp)
1016 : CALL probe_occupancy_kp(wocc(:, :, :), mus(1), kTS, weig(:, :, :), rcoeff(:, :, :, :), icoeff(:, :, :, :), 2.0d0, &
1017 0 : probe, nel, wkp)
1018 : ELSE
1019 0 : nel = REAL(ne_a, KIND=dp) + REAL(ne_b, KIND=dp)
1020 : CALL probe_occupancy_kp(wocc(:, :, :), mu, kTS, weig(:, :, :), rcoeff(:, :, :, :), icoeff(:, :, :, :), 1.0d0, &
1021 0 : probe, nel, wkp)
1022 0 : kTS = kTS/2._dp
1023 0 : mus(1:2) = mu
1024 : END IF
1025 :
1026 0 : DO ikpgr = 1, kplocal
1027 0 : ik = kp_range(1) + ikpgr - 1
1028 0 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1029 0 : DO ispin = 1, nspin
1030 0 : mo_set => kp%mos(1, ispin)
1031 0 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues, occupation_numbers=occupation)
1032 0 : eigenvalues(1:nmo) = weig(1:nmo, ik, ispin)
1033 0 : occupation(1:nmo) = wocc(1:nmo, ik, ispin)
1034 0 : mo_set%kTS = kTS
1035 0 : mo_set%mu = mus(ispin)
1036 :
1037 0 : CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
1038 : !get smatrix for kpoint_env ikp
1039 : CALL cp_fm_get_info(mo_coeff, &
1040 : nrow_global=nrow_global, &
1041 0 : ncol_global=ncol_global)
1042 0 : ALLOCATE (smatrix(nrow_global, ncol_global))
1043 0 : CALL cp_fm_get_submatrix(mo_coeff, smatrix)
1044 :
1045 0 : smatrix(1:nrow_global, 1:ncol_global) = rcoeff(1:nao, 1:nmo, ik, ispin)
1046 0 : DEALLOCATE (smatrix)
1047 :
1048 0 : mo_set => kp%mos(2, ispin)
1049 :
1050 0 : CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
1051 : !get smatrix for kpoint_env ikp
1052 : CALL cp_fm_get_info(mo_coeff, &
1053 : nrow_global=nrow_global, &
1054 0 : ncol_global=ncol_global)
1055 0 : ALLOCATE (smatrix(nrow_global, ncol_global))
1056 0 : CALL cp_fm_get_submatrix(mo_coeff, smatrix)
1057 :
1058 0 : smatrix(1:nrow_global, 1:ncol_global) = icoeff(1:nao, 1:nmo, ik, ispin)
1059 0 : DEALLOCATE (smatrix)
1060 :
1061 0 : mo_set => kp%mos(1, ispin)
1062 :
1063 : END DO
1064 : END DO
1065 :
1066 0 : DEALLOCATE (weig, wocc, rcoeff, icoeff)
1067 :
1068 : END IF
1069 :
1070 : IF (PRESENT(probe) .EQV. .FALSE.) THEN
1071 5666 : IF (smear%do_smear) THEN
1072 : ! finite electronic temperature
1073 1412 : SELECT CASE (smear%method)
1074 : CASE (smear_fermi_dirac)
1075 706 : IF (nspin == 1) THEN
1076 694 : nel = REAL(nelectron, KIND=dp)
1077 : CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
1078 694 : smear%electronic_temperature, 2.0_dp)
1079 12 : ELSE IF (smear%fixed_mag_mom > 0.0_dp) THEN
1080 2 : nel = REAL(ne_a, KIND=dp)
1081 : CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
1082 2 : smear%electronic_temperature, 1.0_dp)
1083 2 : nel = REAL(ne_b, KIND=dp)
1084 : CALL Fermikp(wocc(:, :, 2), mus(2), kTS, weig(:, :, 2), nel, wkp, &
1085 2 : smear%electronic_temperature, 1.0_dp)
1086 : ELSE
1087 10 : nel = REAL(ne_a, KIND=dp) + REAL(ne_b, KIND=dp)
1088 : CALL Fermikp2(wocc(:, :, :), mu, kTS, weig(:, :, :), nel, wkp, &
1089 10 : smear%electronic_temperature)
1090 10 : kTS = kTS/2._dp
1091 30 : mus(1:2) = mu
1092 : END IF
1093 : CASE DEFAULT
1094 706 : CPABORT("kpoints: Selected smearing not (yet) supported")
1095 : END SELECT
1096 : ELSE
1097 : ! fixed occupations (2/1)
1098 4960 : IF (nspin == 1) THEN
1099 4424 : nel = REAL(nelectron, KIND=dp)
1100 4424 : CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, 0.0_dp, 2.0_dp)
1101 : ELSE
1102 536 : nel = REAL(ne_a, KIND=dp)
1103 536 : CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, 0.0_dp, 1.0_dp)
1104 536 : nel = REAL(ne_b, KIND=dp)
1105 536 : CALL Fermikp(wocc(:, :, 2), mus(2), kTS, weig(:, :, 2), nel, wkp, 0.0_dp, 1.0_dp)
1106 : END IF
1107 : END IF
1108 21036 : DO ikpgr = 1, kplocal
1109 15370 : ik = kp_range(1) + ikpgr - 1
1110 15370 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1111 38298 : DO ispin = 1, nspin
1112 17262 : mo_set => kp%mos(1, ispin)
1113 17262 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues, occupation_numbers=occupation)
1114 338600 : eigenvalues(1:nmo) = weig(1:nmo, ik, ispin)
1115 338600 : occupation(1:nmo) = wocc(1:nmo, ik, ispin)
1116 17262 : mo_set%kTS = kTS
1117 32632 : mo_set%mu = mus(ispin)
1118 : END DO
1119 : END DO
1120 :
1121 5666 : DEALLOCATE (weig, wocc)
1122 :
1123 : END IF
1124 :
1125 5666 : CALL timestop(handle)
1126 :
1127 16998 : END SUBROUTINE kpoint_set_mo_occupation
1128 :
1129 : ! **************************************************************************************************
1130 : !> \brief Calculate kpoint density matrices (rho(k), owned by kpoint groups)
1131 : !> \param kpoint kpoint environment
1132 : !> \param energy_weighted calculate energy weighted density matrix
1133 : !> \param for_aux_fit ...
1134 : ! **************************************************************************************************
1135 17844 : SUBROUTINE kpoint_density_matrices(kpoint, energy_weighted, for_aux_fit)
1136 :
1137 : TYPE(kpoint_type), POINTER :: kpoint
1138 : LOGICAL, OPTIONAL :: energy_weighted, for_aux_fit
1139 :
1140 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_density_matrices'
1141 :
1142 : INTEGER :: handle, ikpgr, ispin, kplocal, nao, nmo, &
1143 : nspin
1144 : INTEGER, DIMENSION(2) :: kp_range
1145 : LOGICAL :: aux_fit, wtype
1146 5948 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation
1147 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1148 : TYPE(cp_fm_type) :: fwork
1149 : TYPE(cp_fm_type), POINTER :: cpmat, pmat, rpmat
1150 : TYPE(kpoint_env_type), POINTER :: kp
1151 : TYPE(mo_set_type), POINTER :: mo_set
1152 :
1153 5948 : CALL timeset(routineN, handle)
1154 :
1155 5948 : IF (PRESENT(energy_weighted)) THEN
1156 158 : wtype = energy_weighted
1157 : ELSE
1158 : ! default is normal density matrix
1159 : wtype = .FALSE.
1160 : END IF
1161 :
1162 5948 : IF (PRESENT(for_aux_fit)) THEN
1163 108 : aux_fit = for_aux_fit
1164 : ELSE
1165 : aux_fit = .FALSE.
1166 : END IF
1167 :
1168 108 : IF (aux_fit) THEN
1169 108 : CPASSERT(ASSOCIATED(kpoint%kp_aux_env))
1170 : END IF
1171 :
1172 : ! work matrix
1173 5948 : IF (aux_fit) THEN
1174 108 : mo_set => kpoint%kp_aux_env(1)%kpoint_env%mos(1, 1)
1175 : ELSE
1176 5840 : mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)
1177 : END IF
1178 5948 : CALL get_mo_set(mo_set, nao=nao, nmo=nmo)
1179 5948 : CALL cp_fm_get_info(mo_set%mo_coeff, matrix_struct=matrix_struct)
1180 5948 : CALL cp_fm_create(fwork, matrix_struct)
1181 :
1182 5948 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
1183 5948 : kplocal = kp_range(2) - kp_range(1) + 1
1184 23784 : DO ikpgr = 1, kplocal
1185 17836 : IF (aux_fit) THEN
1186 1594 : kp => kpoint%kp_aux_env(ikpgr)%kpoint_env
1187 : ELSE
1188 16242 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1189 : END IF
1190 17836 : nspin = SIZE(kp%mos, 2)
1191 43764 : DO ispin = 1, nspin
1192 19980 : mo_set => kp%mos(1, ispin)
1193 19980 : IF (wtype) THEN
1194 896 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
1195 : END IF
1196 37816 : IF (kpoint%use_real_wfn) THEN
1197 72 : IF (wtype) THEN
1198 10 : pmat => kp%wmat(1, ispin)
1199 : ELSE
1200 62 : pmat => kp%pmat(1, ispin)
1201 : END IF
1202 72 : CALL get_mo_set(mo_set, occupation_numbers=occupation)
1203 72 : CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
1204 72 : CALL cp_fm_column_scale(fwork, occupation)
1205 72 : IF (wtype) THEN
1206 10 : CALL cp_fm_column_scale(fwork, eigenvalues)
1207 : END IF
1208 72 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, pmat)
1209 : ELSE
1210 19908 : IF (wtype) THEN
1211 886 : rpmat => kp%wmat(1, ispin)
1212 886 : cpmat => kp%wmat(2, ispin)
1213 : ELSE
1214 19022 : rpmat => kp%pmat(1, ispin)
1215 19022 : cpmat => kp%pmat(2, ispin)
1216 : END IF
1217 19908 : CALL get_mo_set(mo_set, occupation_numbers=occupation)
1218 19908 : CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
1219 19908 : CALL cp_fm_column_scale(fwork, occupation)
1220 19908 : IF (wtype) THEN
1221 886 : CALL cp_fm_column_scale(fwork, eigenvalues)
1222 : END IF
1223 : ! Re(c)*Re(c)
1224 19908 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, rpmat)
1225 19908 : mo_set => kp%mos(2, ispin)
1226 : ! Im(c)*Re(c)
1227 19908 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, cpmat)
1228 : ! Re(c)*Im(c)
1229 19908 : CALL parallel_gemm("N", "T", nao, nao, nmo, -1.0_dp, fwork, mo_set%mo_coeff, 1.0_dp, cpmat)
1230 19908 : CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
1231 19908 : CALL cp_fm_column_scale(fwork, occupation)
1232 19908 : IF (wtype) THEN
1233 886 : CALL cp_fm_column_scale(fwork, eigenvalues)
1234 : END IF
1235 : ! Im(c)*Im(c)
1236 19908 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 1.0_dp, rpmat)
1237 : END IF
1238 : END DO
1239 : END DO
1240 :
1241 5948 : CALL cp_fm_release(fwork)
1242 :
1243 5948 : CALL timestop(handle)
1244 :
1245 5948 : END SUBROUTINE kpoint_density_matrices
1246 :
1247 : ! **************************************************************************************************
1248 : !> \brief Kpoint transformation of density matrix for Lowdin population analysis
1249 : !> Transforms matrices to kpoint, distributes kpoint groups, performs S^1/2 P S^1/2
1250 : !> \param matrix_p Density matrices (RS indices, global)
1251 : !> \param kpoints Kpoint environment
1252 : !> \param ispin spin component
1253 : !> \param fmwork full matrices distributed over all groups
1254 : !> \par History
1255 : !> 02.2026 created [JGH]
1256 : ! **************************************************************************************************
1257 0 : SUBROUTINE lowdin_kp_trans(matrix_p, kpoints, ispin, fmwork)
1258 :
1259 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
1260 : TYPE(kpoint_type), POINTER :: kpoints
1261 : INTEGER, INTENT(IN) :: ispin
1262 : TYPE(cp_fm_type), DIMENSION(:) :: fmwork
1263 :
1264 : CHARACTER(len=*), PARAMETER :: routineN = 'lowdin_kp_trans'
1265 : COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp, 0.0_dp), &
1266 : czero = (0.0_dp, 0.0_dp)
1267 :
1268 : INTEGER :: handle, igroup, ik, ikp, indx, kplocal, &
1269 : nao, nkp, nkp_groups
1270 : INTEGER, DIMENSION(2) :: kp_range
1271 0 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1272 0 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1273 : LOGICAL :: my_kpgrp, use_real_wfn
1274 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1275 0 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
1276 : TYPE(cp_cfm_type) :: csmat, cwork
1277 0 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools
1278 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1279 : TYPE(cp_fm_type) :: fmdummy, fmlocal, rsmat
1280 : TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix, tmpmat
1281 : TYPE(kpoint_env_type), POINTER :: kp
1282 : TYPE(mp_para_env_type), POINTER :: para_env
1283 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1284 0 : POINTER :: sab_nl
1285 : TYPE(qs_matrix_pools_type), POINTER :: mpools
1286 :
1287 0 : CALL timeset(routineN, handle)
1288 :
1289 0 : NULLIFY (sab_nl)
1290 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
1291 : nkp_groups=nkp_groups, kp_dist=kp_dist, sab_nl=sab_nl, &
1292 0 : cell_to_index=cell_to_index)
1293 0 : CPASSERT(ASSOCIATED(sab_nl))
1294 0 : kplocal = kp_range(2) - kp_range(1) + 1
1295 :
1296 : ! allocate some work matrices
1297 0 : ALLOCATE (rmatrix, cmatrix, tmpmat)
1298 : CALL dbcsr_create(rmatrix, template=matrix_p(1, 1)%matrix, &
1299 0 : matrix_type=dbcsr_type_symmetric)
1300 : CALL dbcsr_create(cmatrix, template=matrix_p(1, 1)%matrix, &
1301 0 : matrix_type=dbcsr_type_antisymmetric)
1302 : CALL dbcsr_create(tmpmat, template=matrix_p(1, 1)%matrix, &
1303 0 : matrix_type=dbcsr_type_no_symmetry)
1304 0 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
1305 0 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
1306 :
1307 : ! fm pools to be used within a kpoint group
1308 0 : CALL get_kpoint_info(kpoints, mpools=mpools)
1309 0 : CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
1310 :
1311 0 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
1312 0 : CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
1313 :
1314 0 : IF (use_real_wfn) THEN
1315 0 : CALL cp_fm_create(rsmat, matrix_struct)
1316 : ELSE
1317 0 : CALL cp_cfm_create(csmat, matrix_struct)
1318 0 : CALL cp_cfm_create(cwork, matrix_struct)
1319 : END IF
1320 :
1321 0 : CALL cp_fm_get_info(fmwork(1), nrow_global=nao)
1322 :
1323 0 : para_env => kpoints%blacs_env_all%para_env
1324 0 : ALLOCATE (info(kplocal*nkp_groups, 2))
1325 :
1326 : ! Setup and start all the communication
1327 0 : indx = 0
1328 0 : DO ikp = 1, kplocal
1329 0 : DO igroup = 1, nkp_groups
1330 : ! number of current kpoint
1331 0 : ik = kp_dist(1, igroup) + ikp - 1
1332 0 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
1333 0 : indx = indx + 1
1334 0 : IF (use_real_wfn) THEN
1335 0 : CALL dbcsr_set(rmatrix, 0.0_dp)
1336 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_p, ispin=ispin, &
1337 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1338 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1339 0 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
1340 : ELSE
1341 0 : CALL dbcsr_set(rmatrix, 0.0_dp)
1342 0 : CALL dbcsr_set(cmatrix, 0.0_dp)
1343 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_p, ispin=ispin, &
1344 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1345 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1346 0 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
1347 0 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
1348 0 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(2))
1349 : END IF
1350 : ! transfer to kpoint group
1351 : ! redistribution of matrices, new blacs environment
1352 : ! fmwork -> fmlocal -> rsmat/csmat
1353 0 : IF (use_real_wfn) THEN
1354 0 : IF (my_kpgrp) THEN
1355 0 : CALL cp_fm_start_copy_general(fmwork(1), rsmat, para_env, info(indx, 1))
1356 : ELSE
1357 0 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
1358 : END IF
1359 : ELSE
1360 0 : IF (my_kpgrp) THEN
1361 0 : CALL cp_fm_start_copy_general(fmwork(1), fmlocal, para_env, info(indx, 1))
1362 0 : CALL cp_fm_start_copy_general(fmwork(2), fmlocal, para_env, info(indx, 2))
1363 : ELSE
1364 0 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
1365 0 : CALL cp_fm_start_copy_general(fmwork(2), fmdummy, para_env, info(indx, 2))
1366 : END IF
1367 : END IF
1368 : END DO
1369 : END DO
1370 :
1371 : ! Finish communication then tranform in each group
1372 : indx = 0
1373 0 : DO ikp = 1, kplocal
1374 0 : DO igroup = 1, nkp_groups
1375 : ! number of current kpoint
1376 0 : ik = kp_dist(1, igroup) + ikp - 1
1377 0 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
1378 0 : indx = indx + 1
1379 0 : IF (my_kpgrp) THEN
1380 0 : IF (use_real_wfn) THEN
1381 0 : CALL cp_fm_finish_copy_general(rsmat, info(indx, 1))
1382 : ELSE
1383 0 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
1384 0 : CALL cp_cfm_scale_and_add_fm(czero, csmat, cone, fmlocal)
1385 0 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
1386 0 : CALL cp_cfm_scale_and_add_fm(cone, csmat, gaussi, fmlocal)
1387 : END IF
1388 : END IF
1389 : END DO
1390 :
1391 : ! Each kpoint group has now information on a kpoint to be transformed
1392 0 : kp => kpoints%kp_env(ikp)%kpoint_env
1393 0 : IF (use_real_wfn) THEN
1394 : CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, rsmat, kp%shalf, &
1395 0 : 0.0_dp, fmlocal)
1396 : CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, kp%shalf, fmlocal, &
1397 0 : 0.0_dp, kp%pmat(1, ispin))
1398 : ELSE
1399 : CALL parallel_gemm("N", "N", nao, nao, nao, cone, csmat, kp%cshalf, &
1400 0 : czero, cwork)
1401 : CALL parallel_gemm("N", "N", nao, nao, nao, cone, kp%cshalf, cwork, &
1402 0 : czero, csmat)
1403 0 : CALL cp_cfm_to_fm(csmat, kp%pmat(1, ispin), kp%pmat(2, ispin))
1404 : END IF
1405 : END DO
1406 :
1407 : ! Clean up communication
1408 : indx = 0
1409 0 : DO ikp = 1, kplocal
1410 0 : DO igroup = 1, nkp_groups
1411 : ! number of current kpoint
1412 0 : ik = kp_dist(1, igroup) + ikp - 1
1413 0 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
1414 0 : indx = indx + 1
1415 0 : IF (use_real_wfn) THEN
1416 0 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
1417 : ELSE
1418 0 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
1419 0 : CALL cp_fm_cleanup_copy_general(info(indx, 2))
1420 : END IF
1421 : END DO
1422 : END DO
1423 :
1424 : ! All done
1425 0 : DEALLOCATE (info)
1426 :
1427 0 : CALL dbcsr_deallocate_matrix(rmatrix)
1428 0 : CALL dbcsr_deallocate_matrix(cmatrix)
1429 0 : CALL dbcsr_deallocate_matrix(tmpmat)
1430 :
1431 0 : IF (use_real_wfn) THEN
1432 0 : CALL cp_fm_release(rsmat)
1433 : ELSE
1434 0 : CALL cp_cfm_release(csmat)
1435 0 : CALL cp_cfm_release(cwork)
1436 : END IF
1437 0 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
1438 :
1439 0 : CALL timestop(handle)
1440 :
1441 0 : END SUBROUTINE lowdin_kp_trans
1442 :
1443 : ! **************************************************************************************************
1444 : !> \brief generate real space density matrices in DBCSR format
1445 : !> \param kpoint Kpoint environment
1446 : !> \param denmat Real space (DBCSR) density matrices
1447 : !> \param wtype True = energy weighted density matrix
1448 : !> False = normal density matrix
1449 : !> \param tempmat DBCSR matrix to be used as template
1450 : !> \param sab_nl ...
1451 : !> \param fmwork FM work matrices (kpoint group)
1452 : !> \param for_aux_fit ...
1453 : !> \param pmat_ext ...
1454 : ! **************************************************************************************************
1455 6168 : SUBROUTINE kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwork, for_aux_fit, pmat_ext)
1456 :
1457 : TYPE(kpoint_type), POINTER :: kpoint
1458 : TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
1459 : LOGICAL, INTENT(IN) :: wtype
1460 : TYPE(dbcsr_type), POINTER :: tempmat
1461 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1462 : POINTER :: sab_nl
1463 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fmwork
1464 : LOGICAL, OPTIONAL :: for_aux_fit
1465 : TYPE(cp_fm_type), DIMENSION(:, :, :), INTENT(IN), &
1466 : OPTIONAL :: pmat_ext
1467 :
1468 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_density_transform'
1469 :
1470 : INTEGER :: handle, ic, ik, ikk, indx, ir, ira, is, &
1471 : ispin, jr, nc, nimg, nkp, nspin
1472 6168 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1473 : LOGICAL :: aux_fit, do_ext, do_symmetric, my_kpgrp, &
1474 : real_only
1475 : REAL(KIND=dp) :: wkpx
1476 6168 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
1477 6168 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1478 6168 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:) :: info
1479 : TYPE(cp_fm_type) :: fmdummy
1480 : TYPE(dbcsr_type), POINTER :: cpmat, rpmat, scpmat, srpmat
1481 6168 : TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kind_rot
1482 : TYPE(kpoint_env_type), POINTER :: kp
1483 : TYPE(kpoint_sym_type), POINTER :: kpsym
1484 : TYPE(mp_para_env_type), POINTER :: para_env
1485 :
1486 6168 : CALL timeset(routineN, handle)
1487 :
1488 6168 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1489 :
1490 6168 : IF (PRESENT(for_aux_fit)) THEN
1491 328 : aux_fit = for_aux_fit
1492 : ELSE
1493 : aux_fit = .FALSE.
1494 : END IF
1495 :
1496 6168 : do_ext = .FALSE.
1497 6168 : IF (PRESENT(pmat_ext)) do_ext = .TRUE.
1498 :
1499 6168 : IF (aux_fit) THEN
1500 190 : CPASSERT(ASSOCIATED(kpoint%kp_aux_env))
1501 : END IF
1502 :
1503 : ! work storage
1504 6168 : ALLOCATE (rpmat)
1505 : CALL dbcsr_create(rpmat, template=tempmat, &
1506 6220 : matrix_type=MERGE(dbcsr_type_symmetric, dbcsr_type_no_symmetry, do_symmetric))
1507 6168 : CALL cp_dbcsr_alloc_block_from_nbl(rpmat, sab_nl)
1508 6168 : CALL dbcsr_set(rpmat, 0.0_dp)
1509 6168 : ALLOCATE (cpmat)
1510 : CALL dbcsr_create(cpmat, template=tempmat, &
1511 6220 : matrix_type=MERGE(dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, do_symmetric))
1512 6168 : CALL cp_dbcsr_alloc_block_from_nbl(cpmat, sab_nl)
1513 6168 : CALL dbcsr_set(cpmat, 0.0_dp)
1514 6168 : IF (.NOT. kpoint%full_grid) THEN
1515 4936 : ALLOCATE (srpmat)
1516 4936 : CALL dbcsr_create(srpmat, template=rpmat)
1517 4936 : CALL cp_dbcsr_alloc_block_from_nbl(srpmat, sab_nl)
1518 4936 : CALL dbcsr_set(srpmat, 0.0_dp)
1519 4936 : ALLOCATE (scpmat)
1520 4936 : CALL dbcsr_create(scpmat, template=cpmat)
1521 4936 : CALL cp_dbcsr_alloc_block_from_nbl(scpmat, sab_nl)
1522 4936 : CALL dbcsr_set(scpmat, 0.0_dp)
1523 : END IF
1524 :
1525 : CALL get_kpoint_info(kpoint, nkp=nkp, xkp=xkp, wkp=wkp, &
1526 6168 : cell_to_index=cell_to_index)
1527 : ! initialize real space density matrices
1528 6168 : IF (aux_fit) THEN
1529 190 : kp => kpoint%kp_aux_env(1)%kpoint_env
1530 : ELSE
1531 5978 : kp => kpoint%kp_env(1)%kpoint_env
1532 : END IF
1533 6168 : nspin = SIZE(kp%mos, 2)
1534 6168 : nc = SIZE(kp%mos, 1)
1535 6168 : nimg = SIZE(denmat, 2)
1536 6168 : real_only = (nc == 1)
1537 :
1538 6168 : para_env => kpoint%blacs_env_all%para_env
1539 132020 : ALLOCATE (info(nspin*nkp*nc))
1540 :
1541 : ! Start all the communication
1542 6168 : indx = 0
1543 12966 : DO ispin = 1, nspin
1544 554120 : DO ic = 1, nimg
1545 554120 : CALL dbcsr_set(denmat(ispin, ic)%matrix, 0.0_dp)
1546 : END DO
1547 : !
1548 45088 : DO ik = 1, nkp
1549 32122 : my_kpgrp = (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2))
1550 : IF (my_kpgrp) THEN
1551 22654 : ikk = ik - kpoint%kp_range(1) + 1
1552 22654 : IF (aux_fit) THEN
1553 2412 : kp => kpoint%kp_aux_env(ikk)%kpoint_env
1554 : ELSE
1555 20242 : kp => kpoint%kp_env(ikk)%kpoint_env
1556 : END IF
1557 : ELSE
1558 : NULLIFY (kp)
1559 : END IF
1560 : ! collect this density matrix on all processors
1561 32122 : CPASSERT(SIZE(fmwork) >= nc)
1562 :
1563 38920 : IF (my_kpgrp) THEN
1564 67890 : DO ic = 1, nc
1565 45236 : indx = indx + 1
1566 67890 : IF (do_ext) THEN
1567 4824 : CALL cp_fm_start_copy_general(pmat_ext(ikk, ic, ispin), fmwork(ic), para_env, info(indx))
1568 : ELSE
1569 40412 : IF (wtype) THEN
1570 1782 : CALL cp_fm_start_copy_general(kp%wmat(ic, ispin), fmwork(ic), para_env, info(indx))
1571 : ELSE
1572 38630 : CALL cp_fm_start_copy_general(kp%pmat(ic, ispin), fmwork(ic), para_env, info(indx))
1573 : END IF
1574 : END IF
1575 : END DO
1576 : ELSE
1577 28404 : DO ic = 1, nc
1578 18936 : indx = indx + 1
1579 28404 : CALL cp_fm_start_copy_general(fmdummy, fmwork(ic), para_env, info(indx))
1580 : END DO
1581 : END IF
1582 : END DO
1583 : END DO
1584 :
1585 : ! Finish communication and transform the received matrices
1586 6168 : indx = 0
1587 12966 : DO ispin = 1, nspin
1588 45088 : DO ik = 1, nkp
1589 96294 : DO ic = 1, nc
1590 64172 : indx = indx + 1
1591 96294 : CALL cp_fm_finish_copy_general(fmwork(ic), info(indx))
1592 : END DO
1593 :
1594 : ! reduce to dbcsr storage
1595 32122 : IF (real_only) THEN
1596 72 : CALL copy_fm_to_dbcsr(fmwork(1), rpmat, keep_sparsity=.TRUE.)
1597 : ELSE
1598 32050 : CALL copy_fm_to_dbcsr(fmwork(1), rpmat, keep_sparsity=.TRUE.)
1599 32050 : CALL copy_fm_to_dbcsr(fmwork(2), cpmat, keep_sparsity=.TRUE.)
1600 : END IF
1601 :
1602 : ! symmetrization
1603 32122 : kpsym => kpoint%kp_sym(ik)%kpoint_sym
1604 32122 : CPASSERT(ASSOCIATED(kpsym))
1605 :
1606 38920 : IF (kpsym%apply_symmetry) THEN
1607 0 : wkpx = wkp(ik)/REAL(kpsym%nwght, KIND=dp)
1608 0 : DO is = 1, kpsym%nwght
1609 0 : ir = ABS(kpsym%rotp(is))
1610 0 : ira = 0
1611 0 : DO jr = 1, SIZE(kpoint%ibrot)
1612 0 : IF (ir == kpoint%ibrot(jr)) ira = jr
1613 : END DO
1614 0 : CPASSERT(ira > 0)
1615 0 : kind_rot => kpoint%kind_rotmat(ira, :)
1616 0 : IF (real_only) THEN
1617 : CALL symtrans(srpmat, rpmat, kind_rot, kpsym%rot(1:3, 1:3, is), &
1618 0 : kpsym%f0(:, is), kpoint%atype, symmetric=.TRUE.)
1619 : ELSE
1620 : CALL symtrans(srpmat, rpmat, kind_rot, kpsym%rot(1:3, 1:3, is), &
1621 0 : kpsym%f0(:, is), kpoint%atype, symmetric=.TRUE.)
1622 : CALL symtrans(scpmat, cpmat, kind_rot, kpsym%rot(1:3, 1:3, is), &
1623 0 : kpsym%f0(:, is), kpoint%atype, antisymmetric=.TRUE.)
1624 : END IF
1625 : CALL transform_dmat(denmat, srpmat, scpmat, ispin, real_only, sab_nl, &
1626 0 : cell_to_index, kpsym%xkp(1:3, is), wkpx)
1627 : END DO
1628 : ELSE
1629 : ! transformation
1630 : CALL transform_dmat(denmat, rpmat, cpmat, ispin, real_only, sab_nl, &
1631 32122 : cell_to_index, xkp(1:3, ik), wkp(ik))
1632 : END IF
1633 : END DO
1634 : END DO
1635 :
1636 : ! Clean up communication
1637 6168 : indx = 0
1638 12966 : DO ispin = 1, nspin
1639 45088 : DO ik = 1, nkp
1640 32122 : my_kpgrp = (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2))
1641 6798 : IF (my_kpgrp) THEN
1642 22654 : ikk = ik - kpoint%kp_range(1) + 1
1643 : IF (aux_fit) THEN
1644 22654 : kp => kpoint%kp_aux_env(ikk)%kpoint_env
1645 : ELSE
1646 22654 : kp => kpoint%kp_env(ikk)%kpoint_env
1647 : END IF
1648 :
1649 67890 : DO ic = 1, nc
1650 45236 : indx = indx + 1
1651 67890 : CALL cp_fm_cleanup_copy_general(info(indx))
1652 : END DO
1653 : ELSE
1654 : ! calls with dummy arguments, so not included
1655 : ! therefore just increment counter by trip count
1656 9468 : indx = indx + nc
1657 : END IF
1658 : END DO
1659 : END DO
1660 :
1661 : ! All done
1662 70340 : DEALLOCATE (info)
1663 :
1664 6168 : CALL dbcsr_deallocate_matrix(rpmat)
1665 6168 : CALL dbcsr_deallocate_matrix(cpmat)
1666 6168 : IF (.NOT. kpoint%full_grid) THEN
1667 4936 : CALL dbcsr_deallocate_matrix(srpmat)
1668 4936 : CALL dbcsr_deallocate_matrix(scpmat)
1669 : END IF
1670 :
1671 6168 : CALL timestop(handle)
1672 :
1673 6168 : END SUBROUTINE kpoint_density_transform
1674 :
1675 : ! **************************************************************************************************
1676 : !> \brief real space density matrices in DBCSR format
1677 : !> \param denmat Real space (DBCSR) density matrix
1678 : !> \param rpmat ...
1679 : !> \param cpmat ...
1680 : !> \param ispin ...
1681 : !> \param real_only ...
1682 : !> \param sab_nl ...
1683 : !> \param cell_to_index ...
1684 : !> \param xkp ...
1685 : !> \param wkp ...
1686 : ! **************************************************************************************************
1687 32122 : SUBROUTINE transform_dmat(denmat, rpmat, cpmat, ispin, real_only, sab_nl, cell_to_index, xkp, wkp)
1688 :
1689 : TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
1690 : TYPE(dbcsr_type), POINTER :: rpmat, cpmat
1691 : INTEGER, INTENT(IN) :: ispin
1692 : LOGICAL, INTENT(IN) :: real_only
1693 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1694 : POINTER :: sab_nl
1695 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1696 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
1697 : REAL(KIND=dp), INTENT(IN) :: wkp
1698 :
1699 : CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_dmat'
1700 :
1701 : INTEGER :: handle, iatom, icell, icol, irow, jatom, &
1702 : nimg
1703 : INTEGER, DIMENSION(3) :: cell
1704 : LOGICAL :: do_symmetric, found
1705 : REAL(KIND=dp) :: arg, coskl, fc, sinkl
1706 32122 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, dblock, rblock
1707 : TYPE(neighbor_list_iterator_p_type), &
1708 32122 : DIMENSION(:), POINTER :: nl_iterator
1709 :
1710 32122 : CALL timeset(routineN, handle)
1711 :
1712 32122 : nimg = SIZE(denmat, 2)
1713 :
1714 : ! transformation
1715 32122 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1716 32122 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
1717 12128844 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1718 12096722 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
1719 :
1720 : !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
1721 : !Therefore, we have: S(R) = sum_k Re(S(k))*cos(k*R) -i^2*Im(S(k))*sin(k*R)
1722 : ! = sum_k Re(S(k))*cos(k*R) + Im(S(k))*sin(k*R)
1723 : !fc = +- 1 is due to the usual non-symmetric real-space matrices stored as symmetric ones
1724 :
1725 12096722 : irow = iatom
1726 12096722 : icol = jatom
1727 12096722 : fc = 1.0_dp
1728 12096722 : IF (do_symmetric .AND. iatom > jatom) THEN
1729 5076237 : irow = jatom
1730 5076237 : icol = iatom
1731 5076237 : fc = -1.0_dp
1732 : END IF
1733 :
1734 12096722 : icell = cell_to_index(cell(1), cell(2), cell(3))
1735 12096722 : IF (icell < 1 .OR. icell > nimg) CYCLE
1736 :
1737 12095444 : arg = REAL(cell(1), dp)*xkp(1) + REAL(cell(2), dp)*xkp(2) + REAL(cell(3), dp)*xkp(3)
1738 12095444 : coskl = wkp*COS(twopi*arg)
1739 12095444 : sinkl = wkp*fc*SIN(twopi*arg)
1740 :
1741 : CALL dbcsr_get_block_p(matrix=denmat(ispin, icell)%matrix, row=irow, col=icol, &
1742 12095444 : block=dblock, found=found)
1743 12095444 : IF (.NOT. found) CYCLE
1744 :
1745 12127566 : IF (real_only) THEN
1746 176136 : CALL dbcsr_get_block_p(matrix=rpmat, row=irow, col=icol, block=rblock, found=found)
1747 176136 : IF (.NOT. found) CYCLE
1748 92594280 : dblock = dblock + coskl*rblock
1749 : ELSE
1750 11919308 : CALL dbcsr_get_block_p(matrix=rpmat, row=irow, col=icol, block=rblock, found=found)
1751 11919308 : IF (.NOT. found) CYCLE
1752 11919308 : CALL dbcsr_get_block_p(matrix=cpmat, row=irow, col=icol, block=cblock, found=found)
1753 11919308 : IF (.NOT. found) CYCLE
1754 2365081102 : dblock = dblock + coskl*rblock
1755 2365081102 : dblock = dblock + sinkl*cblock
1756 : END IF
1757 : END DO
1758 32122 : CALL neighbor_list_iterator_release(nl_iterator)
1759 :
1760 32122 : CALL timestop(handle)
1761 :
1762 32122 : END SUBROUTINE transform_dmat
1763 :
1764 : ! **************************************************************************************************
1765 : !> \brief Symmetrization of density matrix - transform to new k-point
1766 : !> \param smat density matrix at new kpoint
1767 : !> \param pmat reference density matrix
1768 : !> \param kmat Kind type rotation matrix
1769 : !> \param rot Rotation matrix
1770 : !> \param f0 Permutation of atoms under transformation
1771 : !> \param atype Atom to kind pointer
1772 : !> \param symmetric Symmetric matrix
1773 : !> \param antisymmetric Anti-Symmetric matrix
1774 : ! **************************************************************************************************
1775 0 : SUBROUTINE symtrans(smat, pmat, kmat, rot, f0, atype, symmetric, antisymmetric)
1776 : TYPE(dbcsr_type), POINTER :: smat, pmat
1777 : TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kmat
1778 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: rot
1779 : INTEGER, DIMENSION(:), INTENT(IN) :: f0, atype
1780 : LOGICAL, INTENT(IN), OPTIONAL :: symmetric, antisymmetric
1781 :
1782 : CHARACTER(LEN=*), PARAMETER :: routineN = 'symtrans'
1783 :
1784 : INTEGER :: handle, iatom, icol, ikind, ip, irow, &
1785 : jcol, jkind, jp, jrow, natom, numnodes
1786 : LOGICAL :: asym, dorot, found, perm, sym, trans
1787 : REAL(KIND=dp) :: dr, fsign
1788 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: kroti, krotj, pblock, sblock
1789 : TYPE(dbcsr_distribution_type) :: dist
1790 : TYPE(dbcsr_iterator_type) :: iter
1791 :
1792 0 : CALL timeset(routineN, handle)
1793 :
1794 : ! check symmetry options
1795 0 : sym = .FALSE.
1796 0 : IF (PRESENT(symmetric)) sym = symmetric
1797 0 : asym = .FALSE.
1798 0 : IF (PRESENT(antisymmetric)) asym = antisymmetric
1799 :
1800 0 : CPASSERT(.NOT. (sym .AND. asym))
1801 0 : CPASSERT((sym .OR. asym))
1802 :
1803 : ! do we have permutation of atoms
1804 0 : natom = SIZE(f0)
1805 0 : perm = .FALSE.
1806 0 : DO iatom = 1, natom
1807 0 : IF (f0(iatom) == iatom) CYCLE
1808 : perm = .TRUE.
1809 0 : EXIT
1810 : END DO
1811 :
1812 : ! do we have a real rotation
1813 0 : dorot = .FALSE.
1814 0 : IF (ABS(SUM(ABS(rot)) - 3.0_dp) > 1.e-12_dp) dorot = .TRUE.
1815 0 : dr = ABS(rot(1, 1) - 1.0_dp) + ABS(rot(2, 2) - 1.0_dp) + ABS(rot(3, 3) - 1.0_dp)
1816 0 : IF (ABS(dr) > 1.e-12_dp) dorot = .TRUE.
1817 :
1818 0 : fsign = 1.0_dp
1819 0 : IF (asym) fsign = -1.0_dp
1820 :
1821 0 : IF (dorot .OR. perm) THEN
1822 : CALL cp_abort(__LOCATION__, "k-points need FULL_GRID currently. "// &
1823 0 : "Reduced grids not yet working correctly")
1824 0 : CALL dbcsr_set(smat, 0.0_dp)
1825 0 : IF (perm) THEN
1826 0 : CALL dbcsr_get_info(pmat, distribution=dist)
1827 0 : CALL dbcsr_distribution_get(dist, numnodes=numnodes)
1828 0 : IF (numnodes == 1) THEN
1829 : ! the matrices are local to this process
1830 0 : CALL dbcsr_iterator_start(iter, pmat)
1831 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1832 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, pblock)
1833 0 : ip = f0(irow)
1834 0 : jp = f0(icol)
1835 0 : IF (ip <= jp) THEN
1836 0 : jrow = ip
1837 0 : jcol = jp
1838 0 : trans = .FALSE.
1839 : ELSE
1840 0 : jrow = jp
1841 0 : jcol = ip
1842 0 : trans = .TRUE.
1843 : END IF
1844 0 : CALL dbcsr_get_block_p(matrix=smat, row=jrow, col=jcol, BLOCK=sblock, found=found)
1845 0 : CPASSERT(found)
1846 0 : ikind = atype(irow)
1847 0 : jkind = atype(icol)
1848 0 : kroti => kmat(ikind)%rmat
1849 0 : krotj => kmat(jkind)%rmat
1850 : ! rotation
1851 0 : IF (trans) THEN
1852 0 : sblock = fsign*MATMUL(TRANSPOSE(krotj), MATMUL(TRANSPOSE(pblock), kroti))
1853 : ELSE
1854 0 : sblock = fsign*MATMUL(TRANSPOSE(kroti), MATMUL(pblock, krotj))
1855 : END IF
1856 : END DO
1857 0 : CALL dbcsr_iterator_stop(iter)
1858 : !
1859 : ELSE
1860 : ! distributed matrices, most general code needed
1861 : CALL cp_abort(__LOCATION__, "k-points need FULL_GRID currently. "// &
1862 0 : "Reduced grids not yet working correctly")
1863 : END IF
1864 : ELSE
1865 : ! no atom permutations, this is always local
1866 0 : CALL dbcsr_copy(smat, pmat)
1867 0 : CALL dbcsr_iterator_start(iter, smat)
1868 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1869 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
1870 0 : ip = f0(irow)
1871 0 : jp = f0(icol)
1872 0 : IF (ip <= jp) THEN
1873 0 : jrow = ip
1874 0 : jcol = jp
1875 0 : trans = .FALSE.
1876 : ELSE
1877 0 : jrow = jp
1878 0 : jcol = ip
1879 0 : trans = .TRUE.
1880 : END IF
1881 0 : ikind = atype(irow)
1882 0 : jkind = atype(icol)
1883 0 : kroti => kmat(ikind)%rmat
1884 0 : krotj => kmat(jkind)%rmat
1885 : ! rotation
1886 0 : IF (trans) THEN
1887 0 : sblock = fsign*MATMUL(TRANSPOSE(krotj), MATMUL(TRANSPOSE(sblock), kroti))
1888 : ELSE
1889 0 : sblock = fsign*MATMUL(TRANSPOSE(kroti), MATMUL(sblock, krotj))
1890 : END IF
1891 : END DO
1892 0 : CALL dbcsr_iterator_stop(iter)
1893 : !
1894 : END IF
1895 : ELSE
1896 : ! this is the identity operation, just copy the matrix
1897 0 : CALL dbcsr_copy(smat, pmat)
1898 : END IF
1899 :
1900 0 : CALL timestop(handle)
1901 :
1902 0 : END SUBROUTINE symtrans
1903 :
1904 : ! **************************************************************************************************
1905 : !> \brief ...
1906 : !> \param mat ...
1907 : ! **************************************************************************************************
1908 0 : SUBROUTINE matprint(mat)
1909 : TYPE(dbcsr_type), POINTER :: mat
1910 :
1911 : INTEGER :: i, icol, iounit, irow
1912 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: mblock
1913 : TYPE(dbcsr_iterator_type) :: iter
1914 :
1915 0 : iounit = cp_logger_get_default_io_unit()
1916 0 : CALL dbcsr_iterator_start(iter, mat)
1917 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1918 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, mblock)
1919 : !
1920 0 : IF (iounit > 0) THEN
1921 0 : WRITE (iounit, '(A,2I4)') 'BLOCK ', irow, icol
1922 0 : DO i = 1, SIZE(mblock, 1)
1923 0 : WRITE (iounit, '(8F12.6)') mblock(i, :)
1924 : END DO
1925 : END IF
1926 : !
1927 : END DO
1928 0 : CALL dbcsr_iterator_stop(iter)
1929 :
1930 0 : END SUBROUTINE matprint
1931 : ! **************************************************************************************************
1932 :
1933 0 : END MODULE kpoint_methods
|
