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_types, ONLY: cp_cfm_create,&
22 : cp_cfm_release,&
23 : cp_cfm_to_fm,&
24 : cp_cfm_type,&
25 : cp_fm_to_cfm
26 : USE cp_control_types, ONLY: hairy_probes_type
27 : USE cp_dbcsr_api, ONLY: &
28 : dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribution_get, &
29 : dbcsr_distribution_type, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, &
30 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
31 : dbcsr_p_type, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
32 : dbcsr_type_symmetric
33 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
34 : USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr
35 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
36 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
37 : fm_pool_create_fm,&
38 : fm_pool_give_back_fm
39 : USE cp_fm_struct, ONLY: cp_fm_struct_type
40 : USE cp_fm_types, ONLY: &
41 : copy_info_type, cp_fm_cleanup_copy_general, cp_fm_create, cp_fm_finish_copy_general, &
42 : cp_fm_get_diag, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_release, &
43 : cp_fm_start_copy_general, cp_fm_to_fm, cp_fm_type
44 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
45 : USE cryssym, ONLY: crys_sym_gen,&
46 : csym_type,&
47 : kpoint_gen,&
48 : print_crys_symmetry,&
49 : print_kp_symmetry,&
50 : release_csym_type
51 : USE hairy_probes, ONLY: probe_occupancy_kp
52 : USE input_constants, ONLY: smear_fermi_dirac,&
53 : smear_gaussian,&
54 : smear_mp,&
55 : smear_mv
56 : USE kinds, ONLY: dp
57 : USE kpoint_types, ONLY: get_kpoint_info,&
58 : kind_rotmat_type,&
59 : kpoint_env_create,&
60 : kpoint_env_p_type,&
61 : kpoint_env_type,&
62 : kpoint_sym_create,&
63 : kpoint_sym_type,&
64 : kpoint_type
65 : USE mathconstants, ONLY: twopi
66 : USE memory_utilities, ONLY: reallocate
67 : USE message_passing, ONLY: mp_cart_type,&
68 : mp_para_env_type
69 : USE parallel_gemm_api, ONLY: parallel_gemm
70 : USE particle_types, ONLY: particle_type
71 : USE qs_matrix_pools, ONLY: mpools_create,&
72 : mpools_get,&
73 : mpools_rebuild_fm_pools,&
74 : qs_matrix_pools_type
75 : USE qs_mo_types, ONLY: allocate_mo_set,&
76 : get_mo_set,&
77 : init_mo_set,&
78 : mo_set_type,&
79 : set_mo_set
80 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
81 : get_neighbor_list_set_p,&
82 : neighbor_list_iterate,&
83 : neighbor_list_iterator_create,&
84 : neighbor_list_iterator_p_type,&
85 : neighbor_list_iterator_release,&
86 : neighbor_list_set_p_type
87 : USE scf_control_types, ONLY: smear_type
88 : USE smearing_utils, ONLY: Smearkp,&
89 : Smearkp2
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 7750 : 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 7750 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atype
125 : LOGICAL :: spez
126 : REAL(KIND=dp) :: wsum
127 7750 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coord, scoord
128 6843250 : TYPE(csym_type) :: crys_sym
129 : TYPE(kpoint_sym_type), POINTER :: kpsym
130 :
131 7750 : CALL timeset(routineN, handle)
132 :
133 7750 : CPASSERT(ASSOCIATED(kpoint))
134 :
135 7766 : 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 186 : 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 74 : natom = SIZE(particle_set)
161 370 : ALLOCATE (scoord(3, natom), atype(natom))
162 558 : DO i = 1, natom
163 484 : CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=atype(i))
164 558 : CALL real_to_scaled(scoord(1:3, i), particle_set(i)%r(1:3), cell)
165 : END DO
166 : END IF
167 186 : IF (kpoint%verbose) THEN
168 14 : iounit = cp_logger_get_default_io_unit()
169 : ELSE
170 172 : iounit = -1
171 : END IF
172 : ! kind type list
173 558 : ALLOCATE (kpoint%atype(natom))
174 1790 : kpoint%atype = atype
175 :
176 186 : 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 186 : full_grid=kpoint%full_grid)
179 186 : kpoint%nkp = crys_sym%nkpoint
180 930 : ALLOCATE (kpoint%xkp(3, kpoint%nkp), kpoint%wkp(kpoint%nkp))
181 2002 : wsum = SUM(crys_sym%wkpoint)
182 2002 : DO ik = 1, kpoint%nkp
183 7264 : kpoint%xkp(1:3, ik) = crys_sym%xkpoint(1:3, ik)
184 2002 : kpoint%wkp(ik) = crys_sym%wkpoint(ik)/wsum
185 : END DO
186 :
187 : ! print output
188 186 : IF (kpoint%symmetry) CALL print_crys_symmetry(crys_sym)
189 186 : IF (kpoint%symmetry) CALL print_kp_symmetry(crys_sym)
190 :
191 : ! transfer symmetry information
192 2374 : ALLOCATE (kpoint%kp_sym(kpoint%nkp))
193 2002 : DO ik = 1, kpoint%nkp
194 1816 : NULLIFY (kpoint%kp_sym(ik)%kpoint_sym)
195 1816 : CALL kpoint_sym_create(kpoint%kp_sym(ik)%kpoint_sym)
196 1816 : kpsym => kpoint%kp_sym(ik)%kpoint_sym
197 2002 : 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 186 : IF (kpoint%symmetry) THEN
236 558 : nkind = MAXVAL(atype)
237 74 : ns = crys_sym%nrtot
238 306 : ALLOCATE (kpoint%kind_rotmat(ns, nkind))
239 74 : DO i = 1, ns
240 74 : DO j = 1, nkind
241 0 : NULLIFY (kpoint%kind_rotmat(i, j)%rmat)
242 : END DO
243 : END DO
244 148 : ALLOCATE (kpoint%ibrot(ns))
245 74 : kpoint%ibrot(1:ns) = crys_sym%ibrot(1:ns)
246 : END IF
247 :
248 186 : CALL release_csym_type(crys_sym)
249 186 : 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 7750 : CPASSERT(.FALSE.)
260 : END SELECT
261 :
262 : ! check for consistency of options
263 7766 : 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 7750 : CPASSERT(kpoint%nkp >= 1)
276 : END SELECT
277 7750 : 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 7750 : CALL timestop(handle)
294 :
295 7750 : 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 262 : 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 262 : TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_aux_env, kp_env
318 : TYPE(kpoint_env_type), POINTER :: kp
319 262 : TYPE(mp_cart_type) :: comm_cart
320 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp, para_env_kp
321 :
322 262 : CALL timeset(routineN, handle)
323 :
324 262 : IF (PRESENT(with_aux_fit)) THEN
325 198 : aux_fit = with_aux_fit
326 : ELSE
327 : aux_fit = .FALSE.
328 : END IF
329 :
330 262 : kpoint%para_env => para_env
331 262 : CALL kpoint%para_env%retain()
332 262 : kpoint%blacs_env_all => blacs_env
333 262 : CALL kpoint%blacs_env_all%retain()
334 :
335 262 : CPASSERT(.NOT. ASSOCIATED(kpoint%kp_env))
336 262 : IF (aux_fit) THEN
337 30 : CPASSERT(.NOT. ASSOCIATED(kpoint%kp_aux_env))
338 : END IF
339 :
340 262 : NULLIFY (kp_env, kp_aux_env)
341 262 : nkp = kpoint%nkp
342 262 : npe = para_env%num_pe
343 262 : 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 262 : 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 136 : ELSE IF (kpoint%parallel_group_size == 0) THEN
397 : ! no parallelization over kpoints
398 88 : 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 262 : nkp_grp = npe/ngr
405 : ! processor dimensions
406 262 : dims(1) = ngr
407 262 : dims(2) = nkp_grp
408 262 : CPASSERT(MOD(nkp, nkp_grp) == 0)
409 262 : nkp_loc = nkp/nkp_grp
410 :
411 262 : 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 262 : CALL comm_cart%create(comm_old=para_env, ndims=2, dims=dims)
417 786 : pos = comm_cart%mepos_cart
418 262 : ALLOCATE (para_env_kp)
419 262 : CALL para_env_kp%from_split(comm_cart, pos(2))
420 262 : ALLOCATE (para_env_inter_kp)
421 262 : CALL para_env_inter_kp%from_split(comm_cart, pos(1))
422 262 : CALL comm_cart%free()
423 :
424 262 : niogrp = 0
425 262 : IF (para_env%is_source()) niogrp = 1
426 262 : CALL para_env_kp%sum(niogrp)
427 262 : kpoint%iogrp = (niogrp == 1)
428 :
429 : ! parallel groups
430 262 : kpoint%para_env_kp => para_env_kp
431 262 : kpoint%para_env_inter_kp => para_env_inter_kp
432 :
433 : ! distribution of kpoints
434 786 : ALLOCATE (kpoint%kp_dist(2, nkp_grp))
435 610 : DO igr = 1, nkp_grp
436 1306 : kpoint%kp_dist(1:2, igr) = get_limit(nkp, nkp_grp, igr - 1)
437 : END DO
438 : ! local kpoints
439 786 : kpoint%kp_range(1:2) = kpoint%kp_dist(1:2, para_env_inter_kp%mepos + 1)
440 :
441 2542 : ALLOCATE (kp_env(nkp_loc))
442 2018 : DO ik = 1, nkp_loc
443 1756 : NULLIFY (kp_env(ik)%kpoint_env)
444 1756 : ikk = kpoint%kp_range(1) + ik - 1
445 1756 : CALL kpoint_env_create(kp_env(ik)%kpoint_env)
446 1756 : kp => kp_env(ik)%kpoint_env
447 1756 : kp%nkpoint = ikk
448 1756 : kp%wkp = kpoint%wkp(ikk)
449 7024 : kp%xkp(1:3) = kpoint%xkp(1:3, ikk)
450 2018 : kp%is_local = (ngr == 1)
451 : END DO
452 262 : kpoint%kp_env => kp_env
453 :
454 262 : 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 262 : unit_nr = cp_logger_get_default_io_unit()
470 :
471 262 : IF (unit_nr > 0 .AND. kpoint%verbose) THEN
472 7 : WRITE (unit_nr, *)
473 7 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Number of kpoint groups ", nkp_grp
474 7 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Size of each kpoint group", ngr
475 7 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Number of kpoints per group", nkp_loc
476 : END IF
477 262 : kpoint%nkp_groups = nkp_grp
478 :
479 : END IF
480 :
481 262 : CALL timestop(handle)
482 :
483 524 : 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 292 : 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 292 : 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 292 : CALL timeset(routineN, handle)
514 :
515 292 : IF (PRESENT(for_aux_fit)) THEN
516 30 : aux_fit = for_aux_fit
517 : ELSE
518 : aux_fit = .FALSE.
519 : END IF
520 :
521 292 : CPASSERT(ASSOCIATED(kpoint))
522 :
523 : IF (.TRUE. .OR. ASSOCIATED(mos(1)%mo_coeff)) THEN
524 292 : IF (aux_fit) THEN
525 30 : CPASSERT(ASSOCIATED(kpoint%kp_aux_env))
526 : END IF
527 :
528 292 : IF (PRESENT(added_mos)) THEN
529 58 : nadd = added_mos
530 : ELSE
531 : nadd = 0
532 : END IF
533 :
534 292 : IF (kpoint%use_real_wfn) THEN
535 : nc = 1
536 : ELSE
537 280 : nc = 2
538 : END IF
539 292 : nspin = SIZE(mos, 1)
540 292 : nkp_loc = kpoint%kp_range(2) - kpoint%kp_range(1) + 1
541 292 : IF (nkp_loc > 0) THEN
542 292 : IF (aux_fit) THEN
543 30 : CPASSERT(SIZE(kpoint%kp_aux_env) == nkp_loc)
544 : ELSE
545 262 : CPASSERT(SIZE(kpoint%kp_env) == nkp_loc)
546 : END IF
547 : ! allocate the mo sets, correct number of kpoints (local), real/complex, spin
548 2262 : DO ik = 1, nkp_loc
549 1970 : IF (aux_fit) THEN
550 214 : kp => kpoint%kp_aux_env(ik)%kpoint_env
551 : ELSE
552 1756 : kp => kpoint%kp_env(ik)%kpoint_env
553 : END IF
554 14592 : ALLOCATE (kp%mos(nc, nspin))
555 4504 : DO is = 1, nspin
556 : CALL get_mo_set(mos(is), nao=nao, nmo=nmo, nelectron=nelectron, &
557 2242 : n_el_f=n_el_f, maxocc=maxocc, flexible_electron_count=flexible_electron_count)
558 2242 : nmo = MIN(nao, nmo + nadd)
559 8682 : DO ic = 1, nc
560 : CALL allocate_mo_set(kp%mos(ic, is), nao, nmo, nelectron, n_el_f, maxocc, &
561 6712 : 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 292 : NULLIFY (blacs_env)
571 292 : IF (ASSOCIATED(kpoint%blacs_env)) THEN
572 30 : blacs_env => kpoint%blacs_env
573 : ELSE
574 262 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=kpoint%para_env_kp)
575 262 : kpoint%blacs_env => blacs_env
576 : END IF
577 :
578 : ! set possible new number of MOs
579 624 : DO is = 1, nspin
580 332 : CALL get_mo_set(mos(is), nmo=nmorig(is))
581 332 : nmo = MIN(nao, nmorig(is) + nadd)
582 624 : 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 292 : NULLIFY (mpools)
587 292 : CALL mpools_create(mpools=mpools)
588 : CALL mpools_rebuild_fm_pools(mpools=mpools, mos=mos, &
589 292 : blacs_env=blacs_env, para_env=kpoint%para_env_kp)
590 :
591 292 : IF (aux_fit) THEN
592 30 : kpoint%mpools_aux_fit => mpools
593 : ELSE
594 262 : kpoint%mpools => mpools
595 : END IF
596 :
597 : ! reset old number of MOs
598 624 : DO is = 1, nspin
599 624 : CALL set_mo_set(mos(is), nmo=nmorig(is))
600 : END DO
601 :
602 : ! allocate density matrices
603 292 : CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
604 292 : ALLOCATE (fmlocal)
605 292 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
606 292 : CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
607 2262 : DO ik = 1, nkp_loc
608 1970 : IF (aux_fit) THEN
609 214 : kp => kpoint%kp_aux_env(ik)%kpoint_env
610 : ELSE
611 1756 : kp => kpoint%kp_env(ik)%kpoint_env
612 : END IF
613 : ! density matrix
614 1970 : CALL cp_fm_release(kp%pmat)
615 14592 : ALLOCATE (kp%pmat(nc, nspin))
616 4212 : DO is = 1, nspin
617 8682 : DO ic = 1, nc
618 6712 : CALL cp_fm_create(kp%pmat(ic, is), matrix_struct)
619 : END DO
620 : END DO
621 : ! energy weighted density matrix
622 1970 : CALL cp_fm_release(kp%wmat)
623 12622 : ALLOCATE (kp%wmat(nc, nspin))
624 4504 : DO is = 1, nspin
625 8682 : DO ic = 1, nc
626 6712 : CALL cp_fm_create(kp%wmat(ic, is), matrix_struct)
627 : END DO
628 : END DO
629 : END DO
630 292 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
631 292 : DEALLOCATE (fmlocal)
632 :
633 : END IF
634 :
635 : END IF
636 :
637 292 : CALL timestop(handle)
638 :
639 292 : END SUBROUTINE kpoint_initialize_mos
640 :
641 : ! **************************************************************************************************
642 : !> \brief ...
643 : !> \param kpoint ...
644 : ! **************************************************************************************************
645 64 : 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 64 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
652 : TYPE(cp_fm_type), POINTER :: mo_coeff
653 64 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: moskp
654 :
655 64 : CALL timeset(routineN, handle)
656 :
657 434 : DO ik = 1, SIZE(kpoint%kp_env)
658 370 : CALL mpools_get(kpoint%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
659 370 : moskp => kpoint%kp_env(ik)%kpoint_env%mos
660 370 : ikk = kpoint%kp_range(1) + ik - 1
661 370 : CPASSERT(ASSOCIATED(moskp))
662 838 : DO ispin = 1, SIZE(moskp, 2)
663 1582 : DO ic = 1, SIZE(moskp, 1)
664 808 : CALL get_mo_set(moskp(ic, ispin), mo_coeff=mo_coeff)
665 1212 : IF (.NOT. ASSOCIATED(mo_coeff)) THEN
666 : CALL init_mo_set(moskp(ic, ispin), &
667 808 : fm_pool=ao_mo_fm_pools(ispin)%pool, name="kpoints")
668 : END IF
669 : END DO
670 : END DO
671 : END DO
672 :
673 64 : CALL timestop(handle)
674 :
675 64 : 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 nimages [output]
684 : ! **************************************************************************************************
685 1080 : SUBROUTINE kpoint_init_cell_index(kpoint, sab_nl, para_env, nimages)
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 : INTEGER, INTENT(OUT) :: nimages
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 1080 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell, list
699 1080 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index, cti
700 : LOGICAL :: new
701 : TYPE(neighbor_list_iterator_p_type), &
702 1080 : DIMENSION(:), POINTER :: nl_iterator
703 :
704 1080 : NULLIFY (cell_to_index, index_to_cell)
705 :
706 1080 : CALL timeset(routineN, handle)
707 :
708 1080 : CPASSERT(ASSOCIATED(kpoint))
709 :
710 1080 : ALLOCATE (list(3, 125))
711 541080 : list = 0
712 1080 : icount = 1
713 :
714 1080 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
715 707931 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
716 706851 : CALL get_iterator_info(nl_iterator, cell=cell)
717 :
718 706851 : new = .TRUE.
719 43240198 : DO ic = 1, icount
720 43146617 : IF (cell(1) == list(1, ic) .AND. cell(2) == list(2, ic) .AND. &
721 93581 : cell(3) == list(3, ic)) THEN
722 : new = .FALSE.
723 : EXIT
724 : END IF
725 : END DO
726 707931 : IF (new) THEN
727 93581 : icount = icount + 1
728 93581 : IF (icount > SIZE(list, 2)) THEN
729 312 : CALL reallocate(list, 1, 3, 1, 2*SIZE(list, 2))
730 : END IF
731 374324 : list(1:3, icount) = cell(1:3)
732 : END IF
733 :
734 : END DO
735 1080 : CALL neighbor_list_iterator_release(nl_iterator)
736 :
737 95741 : itm(1) = MAXVAL(ABS(list(1, 1:icount)))
738 95741 : itm(2) = MAXVAL(ABS(list(2, 1:icount)))
739 95741 : itm(3) = MAXVAL(ABS(list(3, 1:icount)))
740 1080 : CALL para_env%max(itm)
741 4320 : it = MAXVAL(itm(1:3))
742 1080 : IF (ASSOCIATED(kpoint%cell_to_index)) THEN
743 1076 : DEALLOCATE (kpoint%cell_to_index)
744 : END IF
745 5400 : ALLOCATE (kpoint%cell_to_index(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
746 1080 : cell_to_index => kpoint%cell_to_index
747 1080 : cti => cell_to_index
748 227732 : cti(:, :, :) = 0
749 95741 : DO ic = 1, icount
750 94661 : i1 = list(1, ic)
751 94661 : i2 = list(2, ic)
752 94661 : i3 = list(3, ic)
753 95741 : cti(i1, i2, i3) = ic
754 : END DO
755 454384 : CALL para_env%sum(cti)
756 1080 : ncount = 0
757 6524 : DO i1 = -itm(1), itm(1)
758 41672 : DO i2 = -itm(2), itm(2)
759 230980 : DO i3 = -itm(3), itm(3)
760 225536 : IF (cti(i1, i2, i3) == 0) THEN
761 85522 : cti(i1, i2, i3) = 1000000
762 : ELSE
763 104866 : ncount = ncount + 1
764 104866 : cti(i1, i2, i3) = (ABS(i1) + ABS(i2) + ABS(i3))*1000 + ABS(i3)*100 + ABS(i2)*10 + ABS(i1)
765 104866 : 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 1080 : IF (ASSOCIATED(kpoint%index_to_cell)) THEN
772 1080 : DEALLOCATE (kpoint%index_to_cell)
773 : END IF
774 3240 : ALLOCATE (kpoint%index_to_cell(3, ncount))
775 1080 : index_to_cell => kpoint%index_to_cell
776 105946 : DO ic = 1, ncount
777 104866 : cell = MINLOC(cti)
778 104866 : i1 = cell(1) - 1 - itm(1)
779 104866 : i2 = cell(2) - 1 - itm(2)
780 104866 : i3 = cell(3) - 1 - itm(3)
781 104866 : cti(i1, i2, i3) = 1000000
782 104866 : index_to_cell(1, ic) = i1
783 104866 : index_to_cell(2, ic) = i2
784 105946 : index_to_cell(3, ic) = i3
785 : END DO
786 227732 : cti(:, :, :) = 0
787 105946 : DO ic = 1, ncount
788 104866 : i1 = index_to_cell(1, ic)
789 104866 : i2 = index_to_cell(2, ic)
790 104866 : i3 = index_to_cell(3, ic)
791 105946 : cti(i1, i2, i3) = ic
792 : END DO
793 :
794 : ! keep pointer to this neighborlist
795 1080 : kpoint%sab_nl => sab_nl
796 :
797 : ! set number of images
798 1080 : nimages = SIZE(index_to_cell, 2)
799 :
800 1080 : DEALLOCATE (list)
801 :
802 1080 : CALL timestop(handle)
803 :
804 1080 : END SUBROUTINE kpoint_init_cell_index
805 :
806 : ! **************************************************************************************************
807 : !> \brief Transformation of real space matrices to a kpoint
808 : !> \param rmatrix Real part of kpoint matrix
809 : !> \param cmatrix Complex part of kpoint matrix (optional)
810 : !> \param rsmat Real space matrices
811 : !> \param ispin Spin index
812 : !> \param xkp Kpoint coordinates
813 : !> \param cell_to_index mapping of cell indices to RS index
814 : !> \param sab_nl Defining neighbor list
815 : !> \param is_complex Matrix to be transformed is imaginary
816 : !> \param rs_sign Matrix to be transformed is csaled by rs_sign
817 : ! **************************************************************************************************
818 126224 : SUBROUTINE rskp_transform(rmatrix, cmatrix, rsmat, ispin, &
819 : xkp, cell_to_index, sab_nl, is_complex, rs_sign)
820 :
821 : TYPE(dbcsr_type) :: rmatrix
822 : TYPE(dbcsr_type), OPTIONAL :: cmatrix
823 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
824 : INTEGER, INTENT(IN) :: ispin
825 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
826 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
827 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
828 : POINTER :: sab_nl
829 : LOGICAL, INTENT(IN), OPTIONAL :: is_complex
830 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: rs_sign
831 :
832 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rskp_transform'
833 :
834 : INTEGER :: handle, iatom, ic, icol, irow, jatom, &
835 : nimg
836 : INTEGER, DIMENSION(3) :: cell
837 : LOGICAL :: do_symmetric, found, my_complex, &
838 : wfn_real_only
839 : REAL(KIND=dp) :: arg, coskl, fsign, fsym, sinkl
840 63112 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, rblock, rsblock
841 : TYPE(neighbor_list_iterator_p_type), &
842 63112 : DIMENSION(:), POINTER :: nl_iterator
843 :
844 63112 : CALL timeset(routineN, handle)
845 :
846 63112 : my_complex = .FALSE.
847 63112 : IF (PRESENT(is_complex)) my_complex = is_complex
848 :
849 63112 : fsign = 1.0_dp
850 63112 : IF (PRESENT(rs_sign)) fsign = rs_sign
851 :
852 63112 : wfn_real_only = .TRUE.
853 63112 : IF (PRESENT(cmatrix)) wfn_real_only = .FALSE.
854 :
855 63112 : nimg = SIZE(rsmat, 2)
856 :
857 63112 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
858 :
859 63112 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
860 30550931 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
861 30487819 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
862 :
863 : ! fsym = +- 1 is due to real space matrices being non-symmetric (although in a symmtric type)
864 : ! with the link S_mu^0,nu^b = S_nu^0,mu^-b, and the KP matrices beeing Hermitian
865 30487819 : fsym = 1.0_dp
866 30487819 : irow = iatom
867 30487819 : icol = jatom
868 30487819 : IF (do_symmetric .AND. (iatom > jatom)) THEN
869 12842959 : irow = jatom
870 12842959 : icol = iatom
871 12842959 : fsym = -1.0_dp
872 : END IF
873 :
874 30487819 : ic = cell_to_index(cell(1), cell(2), cell(3))
875 30487819 : IF (ic < 1 .OR. ic > nimg) CYCLE
876 :
877 30487075 : arg = REAL(cell(1), dp)*xkp(1) + REAL(cell(2), dp)*xkp(2) + REAL(cell(3), dp)*xkp(3)
878 30487075 : IF (my_complex) THEN
879 0 : coskl = fsign*fsym*COS(twopi*arg)
880 0 : sinkl = fsign*SIN(twopi*arg)
881 : ELSE
882 30487075 : coskl = fsign*COS(twopi*arg)
883 30487075 : sinkl = fsign*fsym*SIN(twopi*arg)
884 : END IF
885 :
886 : CALL dbcsr_get_block_p(matrix=rsmat(ispin, ic)%matrix, row=irow, col=icol, &
887 30487075 : block=rsblock, found=found)
888 30487075 : IF (.NOT. found) CYCLE
889 :
890 30550187 : IF (wfn_real_only) THEN
891 : CALL dbcsr_get_block_p(matrix=rmatrix, row=irow, col=icol, &
892 456120 : block=rblock, found=found)
893 456120 : IF (.NOT. found) CYCLE
894 208940760 : rblock = rblock + coskl*rsblock
895 : ELSE
896 : CALL dbcsr_get_block_p(matrix=rmatrix, row=irow, col=icol, &
897 30030955 : block=rblock, found=found)
898 30030955 : IF (.NOT. found) CYCLE
899 : CALL dbcsr_get_block_p(matrix=cmatrix, row=irow, col=icol, &
900 30030955 : block=cblock, found=found)
901 30030955 : IF (.NOT. found) CYCLE
902 5606550649 : rblock = rblock + coskl*rsblock
903 5606550649 : cblock = cblock + sinkl*rsblock
904 : END IF
905 :
906 : END DO
907 63112 : CALL neighbor_list_iterator_release(nl_iterator)
908 :
909 63112 : CALL timestop(handle)
910 :
911 63112 : END SUBROUTINE rskp_transform
912 :
913 : ! **************************************************************************************************
914 : !> \brief Given the eigenvalues of all kpoints, calculates the occupation numbers
915 : !> \param kpoint Kpoint environment
916 : !> \param smear Smearing information
917 : !> \param probe ...
918 : ! **************************************************************************************************
919 5160 : SUBROUTINE kpoint_set_mo_occupation(kpoint, smear, probe)
920 :
921 : TYPE(kpoint_type), POINTER :: kpoint
922 : TYPE(smear_type) :: smear
923 : TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
924 : POINTER :: probe
925 :
926 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_set_mo_occupation'
927 :
928 : INTEGER :: handle, ik, ikpgr, ispin, kplocal, nao, &
929 : nb, ncol_global, ne_a, ne_b, &
930 : nelectron, nkp, nmo, nrow_global, nspin
931 : INTEGER, DIMENSION(2) :: kp_range
932 : REAL(KIND=dp) :: kTS, mu, mus(2), nel
933 5160 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: smatrix
934 5160 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: weig, wocc
935 5160 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: icoeff, rcoeff
936 5160 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation, wkp
937 : TYPE(cp_fm_type), POINTER :: mo_coeff
938 : TYPE(kpoint_env_type), POINTER :: kp
939 : TYPE(mo_set_type), POINTER :: mo_set
940 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
941 :
942 5160 : CALL timeset(routineN, handle)
943 :
944 : ! first collect all the eigenvalues
945 5160 : CALL get_kpoint_info(kpoint, nkp=nkp)
946 5160 : kp => kpoint%kp_env(1)%kpoint_env
947 5160 : nspin = SIZE(kp%mos, 2)
948 5160 : mo_set => kp%mos(1, 1)
949 5160 : CALL get_mo_set(mo_set, nmo=nmo, nao=nao, nelectron=nelectron)
950 5160 : ne_a = nelectron
951 5160 : IF (nspin == 2) THEN
952 612 : CALL get_mo_set(kp%mos(1, 2), nmo=nb, nelectron=ne_b)
953 612 : CPASSERT(nmo == nb)
954 : END IF
955 41280 : ALLOCATE (weig(nmo, nkp, nspin), wocc(nmo, nkp, nspin))
956 503684 : weig = 0.0_dp
957 503684 : wocc = 0.0_dp
958 5160 : IF (PRESENT(probe)) THEN
959 0 : ALLOCATE (rcoeff(nao, nmo, nkp, nspin), icoeff(nao, nmo, nkp, nspin))
960 0 : rcoeff = 0.0_dp !coeff, real part
961 0 : icoeff = 0.0_dp !coeff, imaginary part
962 : END IF
963 5160 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
964 5160 : kplocal = kp_range(2) - kp_range(1) + 1
965 23566 : DO ikpgr = 1, kplocal
966 18406 : ik = kp_range(1) + ikpgr - 1
967 18406 : kp => kpoint%kp_env(ikpgr)%kpoint_env
968 44248 : DO ispin = 1, nspin
969 20682 : mo_set => kp%mos(1, ispin)
970 20682 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
971 370286 : weig(1:nmo, ik, ispin) = eigenvalues(1:nmo)
972 39088 : IF (PRESENT(probe)) THEN
973 0 : CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
974 : CALL cp_fm_get_info(mo_coeff, &
975 : nrow_global=nrow_global, &
976 0 : ncol_global=ncol_global)
977 0 : ALLOCATE (smatrix(nrow_global, ncol_global))
978 0 : CALL cp_fm_get_submatrix(mo_coeff, smatrix)
979 :
980 0 : rcoeff(1:nao, 1:nmo, ik, ispin) = smatrix(1:nrow_global, 1:ncol_global)
981 :
982 0 : DEALLOCATE (smatrix)
983 :
984 0 : mo_set => kp%mos(2, ispin)
985 :
986 0 : CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
987 : CALL cp_fm_get_info(mo_coeff, &
988 : nrow_global=nrow_global, &
989 0 : ncol_global=ncol_global)
990 0 : ALLOCATE (smatrix(nrow_global, ncol_global))
991 0 : CALL cp_fm_get_submatrix(mo_coeff, smatrix)
992 :
993 0 : icoeff(1:nao, 1:nmo, ik, ispin) = smatrix(1:nrow_global, 1:ncol_global)
994 :
995 0 : mo_set => kp%mos(1, ispin)
996 :
997 0 : DEALLOCATE (smatrix)
998 : END IF
999 : END DO
1000 : END DO
1001 5160 : CALL get_kpoint_info(kpoint, para_env_inter_kp=para_env_inter_kp)
1002 5160 : CALL para_env_inter_kp%sum(weig)
1003 :
1004 5160 : IF (PRESENT(probe)) THEN
1005 0 : CALL para_env_inter_kp%sum(rcoeff)
1006 0 : CALL para_env_inter_kp%sum(icoeff)
1007 : END IF
1008 :
1009 5160 : CALL get_kpoint_info(kpoint, wkp=wkp)
1010 :
1011 : !calling of HP module HERE, before smear
1012 5160 : IF (PRESENT(probe)) THEN
1013 0 : smear%do_smear = .FALSE. !ensures smearing is switched off
1014 :
1015 0 : IF (nspin == 1) THEN
1016 0 : nel = REAL(nelectron, KIND=dp)
1017 : CALL probe_occupancy_kp(wocc(:, :, :), mus(1), kTS, weig(:, :, :), rcoeff(:, :, :, :), icoeff(:, :, :, :), 2.0d0, &
1018 0 : probe, nel, wkp)
1019 : ELSE
1020 0 : nel = REAL(ne_a, KIND=dp) + REAL(ne_b, KIND=dp)
1021 : CALL probe_occupancy_kp(wocc(:, :, :), mu, kTS, weig(:, :, :), rcoeff(:, :, :, :), icoeff(:, :, :, :), 1.0d0, &
1022 0 : probe, nel, wkp)
1023 0 : kTS = kTS/2._dp
1024 0 : mus(1:2) = mu
1025 : END IF
1026 :
1027 0 : DO ikpgr = 1, kplocal
1028 0 : ik = kp_range(1) + ikpgr - 1
1029 0 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1030 0 : DO ispin = 1, nspin
1031 0 : mo_set => kp%mos(1, ispin)
1032 0 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues, occupation_numbers=occupation)
1033 0 : eigenvalues(1:nmo) = weig(1:nmo, ik, ispin)
1034 0 : occupation(1:nmo) = wocc(1:nmo, ik, ispin)
1035 0 : mo_set%kTS = kTS
1036 0 : mo_set%mu = mus(ispin)
1037 :
1038 0 : CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
1039 : !get smatrix for kpoint_env ikp
1040 : CALL cp_fm_get_info(mo_coeff, &
1041 : nrow_global=nrow_global, &
1042 0 : ncol_global=ncol_global)
1043 0 : ALLOCATE (smatrix(nrow_global, ncol_global))
1044 0 : CALL cp_fm_get_submatrix(mo_coeff, smatrix)
1045 :
1046 0 : smatrix(1:nrow_global, 1:ncol_global) = rcoeff(1:nao, 1:nmo, ik, ispin)
1047 0 : DEALLOCATE (smatrix)
1048 :
1049 0 : mo_set => kp%mos(2, ispin)
1050 :
1051 0 : CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
1052 : !get smatrix for kpoint_env ikp
1053 : CALL cp_fm_get_info(mo_coeff, &
1054 : nrow_global=nrow_global, &
1055 0 : ncol_global=ncol_global)
1056 0 : ALLOCATE (smatrix(nrow_global, ncol_global))
1057 0 : CALL cp_fm_get_submatrix(mo_coeff, smatrix)
1058 :
1059 0 : smatrix(1:nrow_global, 1:ncol_global) = icoeff(1:nao, 1:nmo, ik, ispin)
1060 0 : DEALLOCATE (smatrix)
1061 :
1062 0 : mo_set => kp%mos(1, ispin)
1063 :
1064 : END DO
1065 : END DO
1066 :
1067 0 : DEALLOCATE (weig, wocc, rcoeff, icoeff)
1068 :
1069 : END IF
1070 :
1071 : IF (PRESENT(probe) .EQV. .FALSE.) THEN
1072 5160 : IF (smear%do_smear) THEN
1073 1516 : SELECT CASE (smear%method)
1074 : CASE (smear_fermi_dirac)
1075 : ! finite electronic temperature
1076 696 : IF (nspin == 1) THEN
1077 696 : nel = REAL(nelectron, KIND=dp)
1078 : CALL Smearkp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
1079 696 : smear%electronic_temperature, 2.0_dp, smear_fermi_dirac)
1080 0 : ELSE IF (smear%fixed_mag_mom > 0.0_dp) THEN
1081 0 : nel = REAL(ne_a, KIND=dp)
1082 : CALL Smearkp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
1083 0 : smear%electronic_temperature, 1.0_dp, smear_fermi_dirac)
1084 0 : nel = REAL(ne_b, KIND=dp)
1085 : CALL Smearkp(wocc(:, :, 2), mus(2), kTS, weig(:, :, 2), nel, wkp, &
1086 0 : smear%electronic_temperature, 1.0_dp, smear_fermi_dirac)
1087 : ELSE
1088 0 : nel = REAL(ne_a, KIND=dp) + REAL(ne_b, KIND=dp)
1089 : CALL Smearkp2(wocc(:, :, :), mu, kTS, weig(:, :, :), nel, wkp, &
1090 0 : smear%electronic_temperature, smear_fermi_dirac)
1091 0 : kTS = kTS/2._dp
1092 0 : mus(1:2) = mu
1093 : END IF
1094 : CASE (smear_gaussian, smear_mp, smear_mv)
1095 124 : IF (nspin == 1) THEN
1096 96 : nel = REAL(nelectron, KIND=dp)
1097 : CALL Smearkp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
1098 96 : smear%smearing_width, 2.0_dp, smear%method)
1099 28 : ELSE IF (smear%fixed_mag_mom > 0.0_dp) THEN
1100 0 : nel = REAL(ne_a, KIND=dp)
1101 : CALL Smearkp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
1102 0 : smear%smearing_width, 1.0_dp, smear%method)
1103 0 : nel = REAL(ne_b, KIND=dp)
1104 : CALL Smearkp(wocc(:, :, 2), mus(2), kTS, weig(:, :, 2), nel, wkp, &
1105 0 : smear%smearing_width, 1.0_dp, smear%method)
1106 : ELSE
1107 28 : nel = REAL(ne_a, KIND=dp) + REAL(ne_b, KIND=dp)
1108 : CALL Smearkp2(wocc(:, :, :), mu, kTS, weig(:, :, :), nel, wkp, &
1109 28 : smear%smearing_width, smear%method)
1110 28 : kTS = kTS/2._dp
1111 84 : mus(1:2) = mu
1112 : END IF
1113 : CASE DEFAULT
1114 820 : CPABORT("kpoints: Selected smearing not (yet) supported")
1115 : END SELECT
1116 : ELSE
1117 : ! fixed occupations (2/1)
1118 4340 : IF (nspin == 1) THEN
1119 3756 : nel = REAL(nelectron, KIND=dp)
1120 : CALL Smearkp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
1121 3756 : 0.0_dp, 2.0_dp, smear_gaussian)
1122 : ELSE
1123 584 : nel = REAL(ne_a, KIND=dp)
1124 : CALL Smearkp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
1125 584 : 0.0_dp, 1.0_dp, smear_gaussian)
1126 584 : nel = REAL(ne_b, KIND=dp)
1127 : CALL Smearkp(wocc(:, :, 2), mus(2), kTS, weig(:, :, 2), nel, wkp, &
1128 584 : 0.0_dp, 1.0_dp, smear_gaussian)
1129 : END IF
1130 : END IF
1131 23566 : DO ikpgr = 1, kplocal
1132 18406 : ik = kp_range(1) + ikpgr - 1
1133 18406 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1134 44248 : DO ispin = 1, nspin
1135 20682 : mo_set => kp%mos(1, ispin)
1136 20682 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues, occupation_numbers=occupation)
1137 370286 : eigenvalues(1:nmo) = weig(1:nmo, ik, ispin)
1138 370286 : occupation(1:nmo) = wocc(1:nmo, ik, ispin)
1139 20682 : mo_set%kTS = kTS
1140 39088 : mo_set%mu = mus(ispin)
1141 : END DO
1142 : END DO
1143 :
1144 5160 : DEALLOCATE (weig, wocc)
1145 :
1146 : END IF
1147 :
1148 5160 : CALL timestop(handle)
1149 :
1150 15480 : END SUBROUTINE kpoint_set_mo_occupation
1151 :
1152 : ! **************************************************************************************************
1153 : !> \brief Calculate kpoint density matrices (rho(k), owned by kpoint groups)
1154 : !> \param kpoint kpoint environment
1155 : !> \param energy_weighted calculate energy weighted density matrix
1156 : !> \param for_aux_fit ...
1157 : ! **************************************************************************************************
1158 16476 : SUBROUTINE kpoint_density_matrices(kpoint, energy_weighted, for_aux_fit)
1159 :
1160 : TYPE(kpoint_type), POINTER :: kpoint
1161 : LOGICAL, OPTIONAL :: energy_weighted, for_aux_fit
1162 :
1163 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_density_matrices'
1164 :
1165 : INTEGER :: handle, ikpgr, ispin, kplocal, nao, nmo, &
1166 : nspin
1167 : INTEGER, DIMENSION(2) :: kp_range
1168 : LOGICAL :: aux_fit, wtype
1169 5492 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation
1170 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1171 : TYPE(cp_fm_type) :: fwork
1172 : TYPE(cp_fm_type), POINTER :: cpmat, pmat, rpmat
1173 : TYPE(kpoint_env_type), POINTER :: kp
1174 : TYPE(mo_set_type), POINTER :: mo_set
1175 :
1176 5492 : CALL timeset(routineN, handle)
1177 :
1178 5492 : IF (PRESENT(energy_weighted)) THEN
1179 198 : wtype = energy_weighted
1180 : ELSE
1181 : ! default is normal density matrix
1182 : wtype = .FALSE.
1183 : END IF
1184 :
1185 5492 : IF (PRESENT(for_aux_fit)) THEN
1186 118 : aux_fit = for_aux_fit
1187 : ELSE
1188 : aux_fit = .FALSE.
1189 : END IF
1190 :
1191 118 : IF (aux_fit) THEN
1192 118 : CPASSERT(ASSOCIATED(kpoint%kp_aux_env))
1193 : END IF
1194 :
1195 : ! work matrix
1196 5492 : IF (aux_fit) THEN
1197 118 : mo_set => kpoint%kp_aux_env(1)%kpoint_env%mos(1, 1)
1198 : ELSE
1199 5374 : mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)
1200 : END IF
1201 5492 : CALL get_mo_set(mo_set, nao=nao, nmo=nmo)
1202 5492 : CALL cp_fm_get_info(mo_set%mo_coeff, matrix_struct=matrix_struct)
1203 5492 : CALL cp_fm_create(fwork, matrix_struct)
1204 :
1205 5492 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
1206 5492 : kplocal = kp_range(2) - kp_range(1) + 1
1207 26794 : DO ikpgr = 1, kplocal
1208 21302 : IF (aux_fit) THEN
1209 1864 : kp => kpoint%kp_aux_env(ikpgr)%kpoint_env
1210 : ELSE
1211 19438 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1212 : END IF
1213 21302 : nspin = SIZE(kp%mos, 2)
1214 50624 : DO ispin = 1, nspin
1215 23830 : mo_set => kp%mos(1, ispin)
1216 23830 : IF (wtype) THEN
1217 1056 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
1218 : END IF
1219 45132 : IF (kpoint%use_real_wfn) THEN
1220 104 : IF (wtype) THEN
1221 10 : pmat => kp%wmat(1, ispin)
1222 : ELSE
1223 94 : pmat => kp%pmat(1, ispin)
1224 : END IF
1225 104 : CALL get_mo_set(mo_set, occupation_numbers=occupation)
1226 104 : CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
1227 104 : CALL cp_fm_column_scale(fwork, occupation)
1228 104 : IF (wtype) THEN
1229 10 : CALL cp_fm_column_scale(fwork, eigenvalues)
1230 : END IF
1231 104 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, pmat)
1232 : ELSE
1233 23726 : IF (wtype) THEN
1234 1046 : rpmat => kp%wmat(1, ispin)
1235 1046 : cpmat => kp%wmat(2, ispin)
1236 : ELSE
1237 22680 : rpmat => kp%pmat(1, ispin)
1238 22680 : cpmat => kp%pmat(2, ispin)
1239 : END IF
1240 23726 : CALL get_mo_set(mo_set, occupation_numbers=occupation)
1241 23726 : CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
1242 23726 : CALL cp_fm_column_scale(fwork, occupation)
1243 23726 : IF (wtype) THEN
1244 1046 : CALL cp_fm_column_scale(fwork, eigenvalues)
1245 : END IF
1246 : ! Re(c)*Re(c)
1247 23726 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, rpmat)
1248 23726 : mo_set => kp%mos(2, ispin)
1249 : ! Im(c)*Re(c)
1250 23726 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, cpmat)
1251 : ! Re(c)*Im(c)
1252 23726 : CALL parallel_gemm("N", "T", nao, nao, nmo, -1.0_dp, fwork, mo_set%mo_coeff, 1.0_dp, cpmat)
1253 23726 : CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
1254 23726 : CALL cp_fm_column_scale(fwork, occupation)
1255 23726 : IF (wtype) THEN
1256 1046 : CALL cp_fm_column_scale(fwork, eigenvalues)
1257 : END IF
1258 : ! Im(c)*Im(c)
1259 23726 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 1.0_dp, rpmat)
1260 : END IF
1261 : END DO
1262 : END DO
1263 :
1264 5492 : CALL cp_fm_release(fwork)
1265 :
1266 5492 : CALL timestop(handle)
1267 :
1268 5492 : END SUBROUTINE kpoint_density_matrices
1269 :
1270 : ! **************************************************************************************************
1271 : !> \brief Calculate Lowdin transformation of density matrix S^1/2 P S^1/2
1272 : !> Integrate diagonal elements over k-points to get Lowdin charges
1273 : !> \param kpoint kpoint environment
1274 : !> \param pmat_diag Sum over kpoints of diagonal elements
1275 : !> \par History
1276 : !> 04.2026 created [JGH]
1277 : ! **************************************************************************************************
1278 4 : SUBROUTINE lowdin_kp_trans(kpoint, pmat_diag)
1279 :
1280 : TYPE(kpoint_type), POINTER :: kpoint
1281 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: pmat_diag
1282 :
1283 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lowdin_kp_trans'
1284 : COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp, 0.0_dp), &
1285 : czero = (0.0_dp, 0.0_dp)
1286 :
1287 : INTEGER :: handle, ikpgr, ispin, kplocal, nao, nspin
1288 : INTEGER, DIMENSION(2) :: kp_range
1289 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dele
1290 : TYPE(cp_cfm_type) :: cf1work, cf2work
1291 : TYPE(cp_cfm_type), POINTER :: cshalf
1292 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1293 : TYPE(cp_fm_type) :: f1work, f2work
1294 : TYPE(cp_fm_type), POINTER :: cpmat, pmat, rpmat, shalf
1295 : TYPE(kpoint_env_type), POINTER :: kp
1296 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
1297 :
1298 4 : CALL timeset(routineN, handle)
1299 :
1300 4 : nspin = SIZE(pmat_diag, 2)
1301 174 : pmat_diag = 0.0_dp
1302 :
1303 : ! work matrix
1304 : CALL cp_fm_get_info(kpoint%kp_env(1)%kpoint_env%pmat(1, 1), &
1305 4 : matrix_struct=matrix_struct, nrow_global=nao)
1306 4 : IF (kpoint%use_real_wfn) THEN
1307 0 : CALL cp_fm_create(f1work, matrix_struct, nrow=nao, ncol=nao)
1308 0 : CALL cp_fm_create(f2work, matrix_struct, nrow=nao, ncol=nao)
1309 : ELSE
1310 4 : CALL cp_fm_create(f2work, matrix_struct, nrow=nao, ncol=nao)
1311 4 : CALL cp_cfm_create(cf1work, matrix_struct, nrow=nao, ncol=nao)
1312 4 : CALL cp_cfm_create(cf2work, matrix_struct, nrow=nao, ncol=nao)
1313 : END IF
1314 12 : ALLOCATE (dele(nao))
1315 :
1316 4 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
1317 4 : kplocal = kp_range(2) - kp_range(1) + 1
1318 228 : DO ikpgr = 1, kplocal
1319 224 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1320 460 : DO ispin = 1, nspin
1321 232 : IF (kpoint%use_real_wfn) THEN
1322 0 : pmat => kp%pmat(1, ispin)
1323 0 : shalf => kp%shalf
1324 0 : CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, pmat, shalf, 0.0_dp, f1work)
1325 0 : CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, shalf, f1work, 0.0_dp, f2work)
1326 : ELSE
1327 232 : rpmat => kp%pmat(1, ispin)
1328 232 : cpmat => kp%pmat(2, ispin)
1329 232 : cshalf => kp%cshalf
1330 232 : CALL cp_fm_to_cfm(rpmat, cpmat, cf1work)
1331 232 : CALL parallel_gemm("N", "N", nao, nao, nao, cone, cf1work, cshalf, czero, cf2work)
1332 232 : CALL parallel_gemm("N", "N", nao, nao, nao, cone, cshalf, cf2work, czero, cf1work)
1333 232 : CALL cp_cfm_to_fm(cf1work, mtargetr=f2work)
1334 : END IF
1335 232 : CALL cp_fm_get_diag(f2work, dele)
1336 1944 : pmat_diag(1:nao, ispin) = pmat_diag(1:nao, ispin) + kp%wkp*dele(1:nao)
1337 : END DO
1338 : END DO
1339 :
1340 4 : CALL get_kpoint_info(kpoint, para_env_inter_kp=para_env_inter_kp)
1341 344 : CALL para_env_inter_kp%sum(pmat_diag)
1342 :
1343 4 : IF (kpoint%use_real_wfn) THEN
1344 0 : CALL cp_fm_release(f1work)
1345 0 : CALL cp_fm_release(f2work)
1346 : ELSE
1347 4 : CALL cp_fm_release(f2work)
1348 4 : CALL cp_cfm_release(cf1work)
1349 4 : CALL cp_cfm_release(cf2work)
1350 : END IF
1351 4 : DEALLOCATE (dele)
1352 :
1353 4 : CALL timestop(handle)
1354 :
1355 8 : END SUBROUTINE lowdin_kp_trans
1356 :
1357 : ! **************************************************************************************************
1358 : !> \brief generate real space density matrices in DBCSR format
1359 : !> \param kpoint Kpoint environment
1360 : !> \param denmat Real space (DBCSR) density matrices
1361 : !> \param wtype True = energy weighted density matrix
1362 : !> False = normal density matrix
1363 : !> \param tempmat DBCSR matrix to be used as template
1364 : !> \param sab_nl ...
1365 : !> \param fmwork FM work matrices (kpoint group)
1366 : !> \param for_aux_fit ...
1367 : !> \param pmat_ext ...
1368 : ! **************************************************************************************************
1369 5722 : SUBROUTINE kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwork, for_aux_fit, pmat_ext)
1370 :
1371 : TYPE(kpoint_type), POINTER :: kpoint
1372 : TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
1373 : LOGICAL, INTENT(IN) :: wtype
1374 : TYPE(dbcsr_type), POINTER :: tempmat
1375 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1376 : POINTER :: sab_nl
1377 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fmwork
1378 : LOGICAL, OPTIONAL :: for_aux_fit
1379 : TYPE(cp_fm_type), DIMENSION(:, :, :), INTENT(IN), &
1380 : OPTIONAL :: pmat_ext
1381 :
1382 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_density_transform'
1383 :
1384 : INTEGER :: handle, ic, ik, ikk, indx, ir, ira, is, &
1385 : ispin, jr, nc, nimg, nkp, nspin
1386 5722 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1387 : LOGICAL :: aux_fit, do_ext, do_symmetric, my_kpgrp, &
1388 : real_only
1389 : REAL(KIND=dp) :: wkpx
1390 5722 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
1391 5722 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1392 5722 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:) :: info
1393 : TYPE(cp_fm_type) :: fmdummy
1394 : TYPE(dbcsr_type), POINTER :: cpmat, rpmat, scpmat, srpmat
1395 5722 : TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kind_rot
1396 : TYPE(kpoint_env_type), POINTER :: kp
1397 : TYPE(kpoint_sym_type), POINTER :: kpsym
1398 : TYPE(mp_para_env_type), POINTER :: para_env
1399 :
1400 5722 : CALL timeset(routineN, handle)
1401 :
1402 5722 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1403 :
1404 5722 : IF (PRESENT(for_aux_fit)) THEN
1405 348 : aux_fit = for_aux_fit
1406 : ELSE
1407 : aux_fit = .FALSE.
1408 : END IF
1409 :
1410 5722 : do_ext = .FALSE.
1411 5722 : IF (PRESENT(pmat_ext)) do_ext = .TRUE.
1412 :
1413 5722 : IF (aux_fit) THEN
1414 200 : CPASSERT(ASSOCIATED(kpoint%kp_aux_env))
1415 : END IF
1416 :
1417 : ! work storage
1418 5722 : ALLOCATE (rpmat)
1419 : CALL dbcsr_create(rpmat, template=tempmat, &
1420 5774 : matrix_type=MERGE(dbcsr_type_symmetric, dbcsr_type_no_symmetry, do_symmetric))
1421 5722 : CALL cp_dbcsr_alloc_block_from_nbl(rpmat, sab_nl)
1422 5722 : CALL dbcsr_set(rpmat, 0.0_dp)
1423 5722 : ALLOCATE (cpmat)
1424 : CALL dbcsr_create(cpmat, template=tempmat, &
1425 5774 : matrix_type=MERGE(dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, do_symmetric))
1426 5722 : CALL cp_dbcsr_alloc_block_from_nbl(cpmat, sab_nl)
1427 5722 : CALL dbcsr_set(cpmat, 0.0_dp)
1428 5722 : IF (.NOT. kpoint%full_grid) THEN
1429 4066 : ALLOCATE (srpmat)
1430 4066 : CALL dbcsr_create(srpmat, template=rpmat)
1431 4066 : CALL cp_dbcsr_alloc_block_from_nbl(srpmat, sab_nl)
1432 4066 : CALL dbcsr_set(srpmat, 0.0_dp)
1433 4066 : ALLOCATE (scpmat)
1434 4066 : CALL dbcsr_create(scpmat, template=cpmat)
1435 4066 : CALL cp_dbcsr_alloc_block_from_nbl(scpmat, sab_nl)
1436 4066 : CALL dbcsr_set(scpmat, 0.0_dp)
1437 : END IF
1438 :
1439 : CALL get_kpoint_info(kpoint, nkp=nkp, xkp=xkp, wkp=wkp, &
1440 5722 : cell_to_index=cell_to_index)
1441 : ! initialize real space density matrices
1442 5722 : IF (aux_fit) THEN
1443 200 : kp => kpoint%kp_aux_env(1)%kpoint_env
1444 : ELSE
1445 5522 : kp => kpoint%kp_env(1)%kpoint_env
1446 : END IF
1447 5722 : nspin = SIZE(kp%mos, 2)
1448 5722 : nc = SIZE(kp%mos, 1)
1449 5722 : nimg = SIZE(denmat, 2)
1450 5722 : real_only = (nc == 1)
1451 :
1452 5722 : para_env => kpoint%blacs_env_all%para_env
1453 134938 : ALLOCATE (info(nspin*nkp*nc))
1454 :
1455 : ! Start all the communication
1456 5722 : indx = 0
1457 12138 : DO ispin = 1, nspin
1458 634664 : DO ic = 1, nimg
1459 634664 : CALL dbcsr_set(denmat(ispin, ic)%matrix, 0.0_dp)
1460 : END DO
1461 : !
1462 48188 : DO ik = 1, nkp
1463 36050 : my_kpgrp = (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2))
1464 : IF (my_kpgrp) THEN
1465 26774 : ikk = ik - kpoint%kp_range(1) + 1
1466 26774 : IF (aux_fit) THEN
1467 2682 : kp => kpoint%kp_aux_env(ikk)%kpoint_env
1468 : ELSE
1469 24092 : kp => kpoint%kp_env(ikk)%kpoint_env
1470 : END IF
1471 : ELSE
1472 : NULLIFY (kp)
1473 : END IF
1474 : ! collect this density matrix on all processors
1475 36050 : CPASSERT(SIZE(fmwork) >= nc)
1476 :
1477 42466 : IF (my_kpgrp) THEN
1478 80218 : DO ic = 1, nc
1479 53444 : indx = indx + 1
1480 80218 : IF (do_ext) THEN
1481 5364 : CALL cp_fm_start_copy_general(pmat_ext(ikk, ic, ispin), fmwork(ic), para_env, info(indx))
1482 : ELSE
1483 48080 : IF (wtype) THEN
1484 2102 : CALL cp_fm_start_copy_general(kp%wmat(ic, ispin), fmwork(ic), para_env, info(indx))
1485 : ELSE
1486 45978 : CALL cp_fm_start_copy_general(kp%pmat(ic, ispin), fmwork(ic), para_env, info(indx))
1487 : END IF
1488 : END IF
1489 : END DO
1490 : ELSE
1491 27828 : DO ic = 1, nc
1492 18552 : indx = indx + 1
1493 27828 : CALL cp_fm_start_copy_general(fmdummy, fmwork(ic), para_env, info(indx))
1494 : END DO
1495 : END IF
1496 : END DO
1497 : END DO
1498 :
1499 : ! Finish communication and transform the received matrices
1500 5722 : indx = 0
1501 12138 : DO ispin = 1, nspin
1502 48188 : DO ik = 1, nkp
1503 108046 : DO ic = 1, nc
1504 71996 : indx = indx + 1
1505 108046 : CALL cp_fm_finish_copy_general(fmwork(ic), info(indx))
1506 : END DO
1507 :
1508 : ! reduce to dbcsr storage
1509 36050 : IF (real_only) THEN
1510 104 : CALL copy_fm_to_dbcsr(fmwork(1), rpmat, keep_sparsity=.TRUE.)
1511 : ELSE
1512 35946 : CALL copy_fm_to_dbcsr(fmwork(1), rpmat, keep_sparsity=.TRUE.)
1513 35946 : CALL copy_fm_to_dbcsr(fmwork(2), cpmat, keep_sparsity=.TRUE.)
1514 : END IF
1515 :
1516 : ! symmetrization
1517 36050 : kpsym => kpoint%kp_sym(ik)%kpoint_sym
1518 36050 : CPASSERT(ASSOCIATED(kpsym))
1519 :
1520 42466 : IF (kpsym%apply_symmetry) THEN
1521 0 : wkpx = wkp(ik)/REAL(kpsym%nwght, KIND=dp)
1522 0 : DO is = 1, kpsym%nwght
1523 0 : ir = ABS(kpsym%rotp(is))
1524 0 : ira = 0
1525 0 : DO jr = 1, SIZE(kpoint%ibrot)
1526 0 : IF (ir == kpoint%ibrot(jr)) ira = jr
1527 : END DO
1528 0 : CPASSERT(ira > 0)
1529 0 : kind_rot => kpoint%kind_rotmat(ira, :)
1530 0 : IF (real_only) THEN
1531 : CALL symtrans(srpmat, rpmat, kind_rot, kpsym%rot(1:3, 1:3, is), &
1532 0 : kpsym%f0(:, is), kpoint%atype, symmetric=.TRUE.)
1533 : ELSE
1534 : CALL symtrans(srpmat, rpmat, kind_rot, kpsym%rot(1:3, 1:3, is), &
1535 0 : kpsym%f0(:, is), kpoint%atype, symmetric=.TRUE.)
1536 : CALL symtrans(scpmat, cpmat, kind_rot, kpsym%rot(1:3, 1:3, is), &
1537 0 : kpsym%f0(:, is), kpoint%atype, antisymmetric=.TRUE.)
1538 : END IF
1539 : CALL transform_dmat(denmat, srpmat, scpmat, ispin, real_only, sab_nl, &
1540 0 : cell_to_index, kpsym%xkp(1:3, is), wkpx)
1541 : END DO
1542 : ELSE
1543 : ! transformation
1544 : CALL transform_dmat(denmat, rpmat, cpmat, ispin, real_only, sab_nl, &
1545 36050 : cell_to_index, xkp(1:3, ik), wkp(ik))
1546 : END IF
1547 : END DO
1548 : END DO
1549 :
1550 : ! Clean up communication
1551 5722 : indx = 0
1552 12138 : DO ispin = 1, nspin
1553 48188 : DO ik = 1, nkp
1554 36050 : my_kpgrp = (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2))
1555 6416 : IF (my_kpgrp) THEN
1556 26774 : ikk = ik - kpoint%kp_range(1) + 1
1557 : IF (aux_fit) THEN
1558 26774 : kp => kpoint%kp_aux_env(ikk)%kpoint_env
1559 : ELSE
1560 26774 : kp => kpoint%kp_env(ikk)%kpoint_env
1561 : END IF
1562 :
1563 80218 : DO ic = 1, nc
1564 53444 : indx = indx + 1
1565 80218 : CALL cp_fm_cleanup_copy_general(info(indx))
1566 : END DO
1567 : ELSE
1568 : ! calls with dummy arguments, so not included
1569 : ! therefore just increment counter by trip count
1570 9276 : indx = indx + nc
1571 : END IF
1572 : END DO
1573 : END DO
1574 :
1575 : ! All done
1576 77718 : DEALLOCATE (info)
1577 :
1578 5722 : CALL dbcsr_deallocate_matrix(rpmat)
1579 5722 : CALL dbcsr_deallocate_matrix(cpmat)
1580 5722 : IF (.NOT. kpoint%full_grid) THEN
1581 4066 : CALL dbcsr_deallocate_matrix(srpmat)
1582 4066 : CALL dbcsr_deallocate_matrix(scpmat)
1583 : END IF
1584 :
1585 5722 : CALL timestop(handle)
1586 :
1587 5722 : END SUBROUTINE kpoint_density_transform
1588 :
1589 : ! **************************************************************************************************
1590 : !> \brief real space density matrices in DBCSR format
1591 : !> \param denmat Real space (DBCSR) density matrix
1592 : !> \param rpmat ...
1593 : !> \param cpmat ...
1594 : !> \param ispin ...
1595 : !> \param real_only ...
1596 : !> \param sab_nl ...
1597 : !> \param cell_to_index ...
1598 : !> \param xkp ...
1599 : !> \param wkp ...
1600 : ! **************************************************************************************************
1601 36050 : SUBROUTINE transform_dmat(denmat, rpmat, cpmat, ispin, real_only, sab_nl, cell_to_index, xkp, wkp)
1602 :
1603 : TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
1604 : TYPE(dbcsr_type), POINTER :: rpmat, cpmat
1605 : INTEGER, INTENT(IN) :: ispin
1606 : LOGICAL, INTENT(IN) :: real_only
1607 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1608 : POINTER :: sab_nl
1609 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1610 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
1611 : REAL(KIND=dp), INTENT(IN) :: wkp
1612 :
1613 : CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_dmat'
1614 :
1615 : INTEGER :: handle, iatom, icell, icol, irow, jatom, &
1616 : nimg
1617 : INTEGER, DIMENSION(3) :: cell
1618 : LOGICAL :: do_symmetric, found
1619 : REAL(KIND=dp) :: arg, coskl, fc, sinkl
1620 36050 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, dblock, rblock
1621 : TYPE(neighbor_list_iterator_p_type), &
1622 36050 : DIMENSION(:), POINTER :: nl_iterator
1623 :
1624 36050 : CALL timeset(routineN, handle)
1625 :
1626 36050 : nimg = SIZE(denmat, 2)
1627 :
1628 : ! transformation
1629 36050 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1630 36050 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
1631 16335415 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1632 16299365 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
1633 :
1634 : !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
1635 : !Therefore, we have: S(R) = sum_k Re(S(k))*cos(k*R) -i^2*Im(S(k))*sin(k*R)
1636 : ! = sum_k Re(S(k))*cos(k*R) + Im(S(k))*sin(k*R)
1637 : !fc = +- 1 is due to the usual non-symmetric real-space matrices stored as symmetric ones
1638 :
1639 16299365 : irow = iatom
1640 16299365 : icol = jatom
1641 16299365 : fc = 1.0_dp
1642 16299365 : IF (do_symmetric .AND. iatom > jatom) THEN
1643 6852556 : irow = jatom
1644 6852556 : icol = iatom
1645 6852556 : fc = -1.0_dp
1646 : END IF
1647 :
1648 16299365 : icell = cell_to_index(cell(1), cell(2), cell(3))
1649 16299365 : IF (icell < 1 .OR. icell > nimg) CYCLE
1650 :
1651 16298087 : arg = REAL(cell(1), dp)*xkp(1) + REAL(cell(2), dp)*xkp(2) + REAL(cell(3), dp)*xkp(3)
1652 16298087 : coskl = wkp*COS(twopi*arg)
1653 16298087 : sinkl = wkp*fc*SIN(twopi*arg)
1654 :
1655 : CALL dbcsr_get_block_p(matrix=denmat(ispin, icell)%matrix, row=irow, col=icol, &
1656 16298087 : block=dblock, found=found)
1657 16298087 : IF (.NOT. found) CYCLE
1658 :
1659 16334137 : IF (real_only) THEN
1660 252264 : CALL dbcsr_get_block_p(matrix=rpmat, row=irow, col=icol, block=rblock, found=found)
1661 252264 : IF (.NOT. found) CYCLE
1662 120381000 : dblock = dblock + coskl*rblock
1663 : ELSE
1664 16045823 : CALL dbcsr_get_block_p(matrix=rpmat, row=irow, col=icol, block=rblock, found=found)
1665 16045823 : IF (.NOT. found) CYCLE
1666 16045823 : CALL dbcsr_get_block_p(matrix=cpmat, row=irow, col=icol, block=cblock, found=found)
1667 16045823 : IF (.NOT. found) CYCLE
1668 2995013781 : dblock = dblock + coskl*rblock
1669 2995013781 : dblock = dblock + sinkl*cblock
1670 : END IF
1671 : END DO
1672 36050 : CALL neighbor_list_iterator_release(nl_iterator)
1673 :
1674 36050 : CALL timestop(handle)
1675 :
1676 36050 : END SUBROUTINE transform_dmat
1677 :
1678 : ! **************************************************************************************************
1679 : !> \brief Symmetrization of density matrix - transform to new k-point
1680 : !> \param smat density matrix at new kpoint
1681 : !> \param pmat reference density matrix
1682 : !> \param kmat Kind type rotation matrix
1683 : !> \param rot Rotation matrix
1684 : !> \param f0 Permutation of atoms under transformation
1685 : !> \param atype Atom to kind pointer
1686 : !> \param symmetric Symmetric matrix
1687 : !> \param antisymmetric Anti-Symmetric matrix
1688 : ! **************************************************************************************************
1689 0 : SUBROUTINE symtrans(smat, pmat, kmat, rot, f0, atype, symmetric, antisymmetric)
1690 : TYPE(dbcsr_type), POINTER :: smat, pmat
1691 : TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kmat
1692 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: rot
1693 : INTEGER, DIMENSION(:), INTENT(IN) :: f0, atype
1694 : LOGICAL, INTENT(IN), OPTIONAL :: symmetric, antisymmetric
1695 :
1696 : CHARACTER(LEN=*), PARAMETER :: routineN = 'symtrans'
1697 :
1698 : INTEGER :: handle, iatom, icol, ikind, ip, irow, &
1699 : jcol, jkind, jp, jrow, natom, numnodes
1700 : LOGICAL :: asym, dorot, found, perm, sym, trans
1701 : REAL(KIND=dp) :: dr, fsign
1702 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: kroti, krotj, pblock, sblock
1703 : TYPE(dbcsr_distribution_type) :: dist
1704 : TYPE(dbcsr_iterator_type) :: iter
1705 :
1706 0 : CALL timeset(routineN, handle)
1707 :
1708 : ! check symmetry options
1709 0 : sym = .FALSE.
1710 0 : IF (PRESENT(symmetric)) sym = symmetric
1711 0 : asym = .FALSE.
1712 0 : IF (PRESENT(antisymmetric)) asym = antisymmetric
1713 :
1714 0 : CPASSERT(.NOT. (sym .AND. asym))
1715 0 : CPASSERT((sym .OR. asym))
1716 :
1717 : ! do we have permutation of atoms
1718 0 : natom = SIZE(f0)
1719 0 : perm = .FALSE.
1720 0 : DO iatom = 1, natom
1721 0 : IF (f0(iatom) == iatom) CYCLE
1722 : perm = .TRUE.
1723 0 : EXIT
1724 : END DO
1725 :
1726 : ! do we have a real rotation
1727 0 : dorot = .FALSE.
1728 0 : IF (ABS(SUM(ABS(rot)) - 3.0_dp) > 1.e-12_dp) dorot = .TRUE.
1729 0 : dr = ABS(rot(1, 1) - 1.0_dp) + ABS(rot(2, 2) - 1.0_dp) + ABS(rot(3, 3) - 1.0_dp)
1730 0 : IF (ABS(dr) > 1.e-12_dp) dorot = .TRUE.
1731 :
1732 0 : fsign = 1.0_dp
1733 0 : IF (asym) fsign = -1.0_dp
1734 :
1735 0 : IF (dorot .OR. perm) THEN
1736 : CALL cp_abort(__LOCATION__, "k-points need FULL_GRID currently. "// &
1737 0 : "Reduced grids not yet working correctly")
1738 0 : CALL dbcsr_set(smat, 0.0_dp)
1739 0 : IF (perm) THEN
1740 0 : CALL dbcsr_get_info(pmat, distribution=dist)
1741 0 : CALL dbcsr_distribution_get(dist, numnodes=numnodes)
1742 0 : IF (numnodes == 1) THEN
1743 : ! the matrices are local to this process
1744 0 : CALL dbcsr_iterator_start(iter, pmat)
1745 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1746 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, pblock)
1747 0 : ip = f0(irow)
1748 0 : jp = f0(icol)
1749 0 : IF (ip <= jp) THEN
1750 0 : jrow = ip
1751 0 : jcol = jp
1752 0 : trans = .FALSE.
1753 : ELSE
1754 0 : jrow = jp
1755 0 : jcol = ip
1756 0 : trans = .TRUE.
1757 : END IF
1758 0 : CALL dbcsr_get_block_p(matrix=smat, row=jrow, col=jcol, BLOCK=sblock, found=found)
1759 0 : CPASSERT(found)
1760 0 : ikind = atype(irow)
1761 0 : jkind = atype(icol)
1762 0 : kroti => kmat(ikind)%rmat
1763 0 : krotj => kmat(jkind)%rmat
1764 : ! rotation
1765 0 : IF (trans) THEN
1766 0 : sblock = fsign*MATMUL(TRANSPOSE(krotj), MATMUL(TRANSPOSE(pblock), kroti))
1767 : ELSE
1768 0 : sblock = fsign*MATMUL(TRANSPOSE(kroti), MATMUL(pblock, krotj))
1769 : END IF
1770 : END DO
1771 0 : CALL dbcsr_iterator_stop(iter)
1772 : !
1773 : ELSE
1774 : ! distributed matrices, most general code needed
1775 : CALL cp_abort(__LOCATION__, "k-points need FULL_GRID currently. "// &
1776 0 : "Reduced grids not yet working correctly")
1777 : END IF
1778 : ELSE
1779 : ! no atom permutations, this is always local
1780 0 : CALL dbcsr_copy(smat, pmat)
1781 0 : CALL dbcsr_iterator_start(iter, smat)
1782 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1783 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
1784 0 : ip = f0(irow)
1785 0 : jp = f0(icol)
1786 0 : IF (ip <= jp) THEN
1787 0 : jrow = ip
1788 0 : jcol = jp
1789 0 : trans = .FALSE.
1790 : ELSE
1791 0 : jrow = jp
1792 0 : jcol = ip
1793 0 : trans = .TRUE.
1794 : END IF
1795 0 : ikind = atype(irow)
1796 0 : jkind = atype(icol)
1797 0 : kroti => kmat(ikind)%rmat
1798 0 : krotj => kmat(jkind)%rmat
1799 : ! rotation
1800 0 : IF (trans) THEN
1801 0 : sblock = fsign*MATMUL(TRANSPOSE(krotj), MATMUL(TRANSPOSE(sblock), kroti))
1802 : ELSE
1803 0 : sblock = fsign*MATMUL(TRANSPOSE(kroti), MATMUL(sblock, krotj))
1804 : END IF
1805 : END DO
1806 0 : CALL dbcsr_iterator_stop(iter)
1807 : !
1808 : END IF
1809 : ELSE
1810 : ! this is the identity operation, just copy the matrix
1811 0 : CALL dbcsr_copy(smat, pmat)
1812 : END IF
1813 :
1814 0 : CALL timestop(handle)
1815 :
1816 0 : END SUBROUTINE symtrans
1817 :
1818 : ! **************************************************************************************************
1819 : !> \brief ...
1820 : !> \param mat ...
1821 : ! **************************************************************************************************
1822 0 : SUBROUTINE matprint(mat)
1823 : TYPE(dbcsr_type), POINTER :: mat
1824 :
1825 : INTEGER :: i, icol, iounit, irow
1826 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: mblock
1827 : TYPE(dbcsr_iterator_type) :: iter
1828 :
1829 0 : iounit = cp_logger_get_default_io_unit()
1830 0 : CALL dbcsr_iterator_start(iter, mat)
1831 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1832 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, mblock)
1833 : !
1834 0 : IF (iounit > 0) THEN
1835 0 : WRITE (iounit, '(A,2I4)') 'BLOCK ', irow, icol
1836 0 : DO i = 1, SIZE(mblock, 1)
1837 0 : WRITE (iounit, '(8F12.6)') mblock(i, :)
1838 : END DO
1839 : END IF
1840 : !
1841 : END DO
1842 0 : CALL dbcsr_iterator_stop(iter)
1843 :
1844 0 : END SUBROUTINE matprint
1845 : ! **************************************************************************************************
1846 :
1847 0 : END MODULE kpoint_methods
|
