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