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