Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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_control_types, ONLY: dft_control_type
22 : USE cp_dbcsr_api, ONLY: &
23 : dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribution_get, &
24 : dbcsr_distribution_type, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, &
25 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
26 : dbcsr_p_type, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
27 : dbcsr_type_symmetric
28 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
29 : USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr
30 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
31 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
32 : fm_pool_create_fm,&
33 : fm_pool_give_back_fm
34 : USE cp_fm_struct, ONLY: cp_fm_struct_type
35 : USE cp_fm_types, ONLY: &
36 : copy_info_type, cp_fm_cleanup_copy_general, cp_fm_create, cp_fm_finish_copy_general, &
37 : cp_fm_get_info, cp_fm_release, cp_fm_start_copy_general, cp_fm_to_fm, cp_fm_type
38 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
39 : USE cryssym, ONLY: crys_sym_gen,&
40 : csym_type,&
41 : kpoint_gen,&
42 : print_crys_symmetry,&
43 : print_kp_symmetry,&
44 : release_csym_type
45 : USE fermi_utils, ONLY: fermikp,&
46 : fermikp2
47 : USE input_constants, ONLY: smear_fermi_dirac
48 : USE kinds, ONLY: dp
49 : USE kpoint_types, ONLY: get_kpoint_info,&
50 : kind_rotmat_type,&
51 : kpoint_env_create,&
52 : kpoint_env_p_type,&
53 : kpoint_env_type,&
54 : kpoint_sym_create,&
55 : kpoint_sym_type,&
56 : kpoint_type
57 : USE mathconstants, ONLY: twopi
58 : USE memory_utilities, ONLY: reallocate
59 : USE message_passing, ONLY: mp_cart_type,&
60 : mp_para_env_type
61 : USE parallel_gemm_api, ONLY: parallel_gemm
62 : USE particle_types, ONLY: particle_type
63 : USE qs_matrix_pools, ONLY: mpools_create,&
64 : mpools_get,&
65 : mpools_rebuild_fm_pools,&
66 : qs_matrix_pools_type
67 : USE qs_mo_types, ONLY: allocate_mo_set,&
68 : get_mo_set,&
69 : init_mo_set,&
70 : mo_set_type,&
71 : set_mo_set
72 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
73 : get_neighbor_list_set_p,&
74 : neighbor_list_iterate,&
75 : neighbor_list_iterator_create,&
76 : neighbor_list_iterator_p_type,&
77 : neighbor_list_iterator_release,&
78 : neighbor_list_set_p_type
79 : USE scf_control_types, ONLY: smear_type
80 : USE util, ONLY: get_limit
81 : #include "./base/base_uses.f90"
82 :
83 : IMPLICIT NONE
84 :
85 : PRIVATE
86 :
87 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpoint_methods'
88 :
89 : PUBLIC :: kpoint_initialize, kpoint_env_initialize, kpoint_initialize_mos, kpoint_initialize_mo_set
90 : PUBLIC :: kpoint_init_cell_index, kpoint_set_mo_occupation
91 : PUBLIC :: kpoint_density_matrices, kpoint_density_transform
92 : PUBLIC :: rskp_transform
93 :
94 : ! **************************************************************************************************
95 :
96 : CONTAINS
97 :
98 : ! **************************************************************************************************
99 : !> \brief Generate the kpoints and initialize the kpoint environment
100 : !> \param kpoint The kpoint environment
101 : !> \param particle_set Particle types and coordinates
102 : !> \param cell Computational cell information
103 : ! **************************************************************************************************
104 7324 : SUBROUTINE kpoint_initialize(kpoint, particle_set, cell)
105 :
106 : TYPE(kpoint_type), POINTER :: kpoint
107 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
108 : TYPE(cell_type), POINTER :: cell
109 :
110 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_initialize'
111 :
112 : INTEGER :: handle, i, ic, ik, iounit, ir, ira, is, &
113 : j, natom, nkind, nr, ns
114 7324 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atype
115 : LOGICAL :: spez
116 : REAL(KIND=dp) :: wsum
117 7324 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coord, scoord
118 6467092 : TYPE(csym_type) :: crys_sym
119 : TYPE(kpoint_sym_type), POINTER :: kpsym
120 :
121 7324 : CALL timeset(routineN, handle)
122 :
123 7324 : CPASSERT(ASSOCIATED(kpoint))
124 :
125 7340 : SELECT CASE (kpoint%kp_scheme)
126 : CASE ("NONE")
127 : ! do nothing
128 : CASE ("GAMMA")
129 16 : kpoint%nkp = 1
130 16 : ALLOCATE (kpoint%xkp(3, 1), kpoint%wkp(1))
131 64 : kpoint%xkp(1:3, 1) = 0.0_dp
132 16 : kpoint%wkp(1) = 1.0_dp
133 32 : ALLOCATE (kpoint%kp_sym(1))
134 16 : NULLIFY (kpoint%kp_sym(1)%kpoint_sym)
135 16 : CALL kpoint_sym_create(kpoint%kp_sym(1)%kpoint_sym)
136 : CASE ("MONKHORST-PACK", "MACDONALD")
137 :
138 148 : IF (.NOT. kpoint%symmetry) THEN
139 : ! we set up a random molecule to avoid any possible symmetry
140 100 : natom = 10
141 100 : ALLOCATE (coord(3, natom), scoord(3, natom), atype(natom))
142 1100 : DO i = 1, natom
143 1000 : atype(i) = i
144 1000 : coord(1, i) = SIN(i*0.12345_dp)
145 1000 : coord(2, i) = COS(i*0.23456_dp)
146 1000 : coord(3, i) = SIN(i*0.34567_dp)
147 1100 : CALL real_to_scaled(scoord(1:3, i), coord(1:3, i), cell)
148 : END DO
149 : ELSE
150 48 : natom = SIZE(particle_set)
151 240 : ALLOCATE (scoord(3, natom), atype(natom))
152 360 : DO i = 1, natom
153 312 : CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=atype(i))
154 360 : CALL real_to_scaled(scoord(1:3, i), particle_set(i)%r(1:3), cell)
155 : END DO
156 : END IF
157 148 : IF (kpoint%verbose) THEN
158 10 : iounit = cp_logger_get_default_io_unit()
159 : ELSE
160 138 : iounit = -1
161 : END IF
162 : ! kind type list
163 444 : ALLOCATE (kpoint%atype(natom))
164 1460 : kpoint%atype = atype
165 :
166 148 : CALL crys_sym_gen(crys_sym, scoord, atype, cell%hmat, delta=kpoint%eps_geo, iounit=iounit)
167 : CALL kpoint_gen(crys_sym, kpoint%nkp_grid, symm=kpoint%symmetry, shift=kpoint%kp_shift, &
168 148 : full_grid=kpoint%full_grid)
169 148 : kpoint%nkp = crys_sym%nkpoint
170 740 : ALLOCATE (kpoint%xkp(3, kpoint%nkp), kpoint%wkp(kpoint%nkp))
171 1790 : wsum = SUM(crys_sym%wkpoint)
172 1790 : DO ik = 1, kpoint%nkp
173 6568 : kpoint%xkp(1:3, ik) = crys_sym%xkpoint(1:3, ik)
174 1790 : kpoint%wkp(ik) = crys_sym%wkpoint(ik)/wsum
175 : END DO
176 :
177 : ! print output
178 148 : IF (kpoint%symmetry) CALL print_crys_symmetry(crys_sym)
179 148 : IF (kpoint%symmetry) CALL print_kp_symmetry(crys_sym)
180 :
181 : ! transfer symmetry information
182 2086 : ALLOCATE (kpoint%kp_sym(kpoint%nkp))
183 1790 : DO ik = 1, kpoint%nkp
184 1642 : NULLIFY (kpoint%kp_sym(ik)%kpoint_sym)
185 1642 : CALL kpoint_sym_create(kpoint%kp_sym(ik)%kpoint_sym)
186 1642 : kpsym => kpoint%kp_sym(ik)%kpoint_sym
187 1790 : IF (crys_sym%symlib .AND. .NOT. crys_sym%fullgrid .AND. crys_sym%istriz == 1) THEN
188 : ! set up the symmetrization information
189 0 : kpsym%nwght = NINT(crys_sym%wkpoint(ik))
190 0 : ns = kpsym%nwght
191 : !
192 0 : IF (ns > 1) THEN
193 0 : kpsym%apply_symmetry = .TRUE.
194 0 : natom = SIZE(particle_set)
195 0 : ALLOCATE (kpsym%rot(3, 3, ns))
196 0 : ALLOCATE (kpsym%xkp(3, ns))
197 0 : ALLOCATE (kpsym%rotp(ns))
198 0 : ALLOCATE (kpsym%f0(natom, ns))
199 0 : nr = 0
200 0 : DO is = 1, SIZE(crys_sym%kplink, 2)
201 0 : IF (crys_sym%kplink(2, is) == ik) THEN
202 0 : nr = nr + 1
203 0 : ir = crys_sym%kpop(is)
204 0 : ira = ABS(ir)
205 0 : kpsym%rotp(nr) = ir
206 0 : IF (ir > 0) THEN
207 0 : kpsym%rot(1:3, 1:3, nr) = crys_sym%rt(1:3, 1:3, ira)
208 : ELSE
209 0 : kpsym%rot(1:3, 1:3, nr) = -crys_sym%rt(1:3, 1:3, ira)
210 : END IF
211 0 : kpsym%xkp(1:3, nr) = crys_sym%kpmesh(1:3, is)
212 0 : DO ic = 1, crys_sym%nrtot
213 0 : IF (crys_sym%ibrot(ic) == ira) THEN
214 0 : kpsym%f0(1:natom, nr) = crys_sym%f0(1:natom, ic)
215 : EXIT
216 : END IF
217 : END DO
218 : END IF
219 : END DO
220 : ! Reduce inversion?
221 0 : kpsym%nwred = nr
222 : END IF
223 : END IF
224 : END DO
225 148 : IF (kpoint%symmetry) THEN
226 360 : nkind = MAXVAL(atype)
227 48 : ns = crys_sym%nrtot
228 198 : ALLOCATE (kpoint%kind_rotmat(ns, nkind))
229 48 : DO i = 1, ns
230 48 : DO j = 1, nkind
231 0 : NULLIFY (kpoint%kind_rotmat(i, j)%rmat)
232 : END DO
233 : END DO
234 96 : ALLOCATE (kpoint%ibrot(ns))
235 48 : kpoint%ibrot(1:ns) = crys_sym%ibrot(1:ns)
236 : END IF
237 :
238 148 : CALL release_csym_type(crys_sym)
239 148 : DEALLOCATE (scoord, atype)
240 :
241 : CASE ("GENERAL")
242 : ! default: no symmetry settings
243 12 : ALLOCATE (kpoint%kp_sym(kpoint%nkp))
244 8 : DO i = 1, kpoint%nkp
245 6 : NULLIFY (kpoint%kp_sym(i)%kpoint_sym)
246 8 : CALL kpoint_sym_create(kpoint%kp_sym(i)%kpoint_sym)
247 : END DO
248 : CASE DEFAULT
249 7324 : CPASSERT(.FALSE.)
250 : END SELECT
251 :
252 : ! check for consistency of options
253 7340 : SELECT CASE (kpoint%kp_scheme)
254 : CASE ("NONE")
255 : ! don't use k-point code
256 : CASE ("GAMMA")
257 16 : CPASSERT(kpoint%nkp == 1)
258 80 : CPASSERT(SUM(ABS(kpoint%xkp)) <= 1.e-12_dp)
259 16 : CPASSERT(kpoint%wkp(1) == 1.0_dp)
260 16 : CPASSERT(.NOT. kpoint%symmetry)
261 : CASE ("GENERAL")
262 2 : CPASSERT(.NOT. kpoint%symmetry)
263 2 : CPASSERT(kpoint%nkp >= 1)
264 : CASE ("MONKHORST-PACK", "MACDONALD")
265 7324 : CPASSERT(kpoint%nkp >= 1)
266 : END SELECT
267 7324 : IF (kpoint%use_real_wfn) THEN
268 : ! what about inversion symmetry?
269 24 : ikloop: DO ik = 1, kpoint%nkp
270 60 : DO i = 1, 3
271 36 : spez = (kpoint%xkp(i, ik) == 0.0_dp .OR. kpoint%xkp(i, ik) == 0.5_dp)
272 12 : IF (.NOT. spez) EXIT ikloop
273 : END DO
274 : END DO ikloop
275 12 : IF (.NOT. spez) THEN
276 : ! Warning: real wfn might be wrong for this system
277 : CALL cp_warn(__LOCATION__, &
278 : "A calculation using real wavefunctions is requested. "// &
279 0 : "We could not determine if the symmetry of the system allows real wavefunctions. ")
280 : END IF
281 : END IF
282 :
283 7324 : CALL timestop(handle)
284 :
285 7324 : END SUBROUTINE kpoint_initialize
286 :
287 : ! **************************************************************************************************
288 : !> \brief Initialize the kpoint environment
289 : !> \param kpoint Kpoint environment
290 : !> \param para_env ...
291 : !> \param blacs_env ...
292 : !> \param with_aux_fit ...
293 : ! **************************************************************************************************
294 224 : SUBROUTINE kpoint_env_initialize(kpoint, para_env, blacs_env, with_aux_fit)
295 :
296 : TYPE(kpoint_type), INTENT(INOUT) :: kpoint
297 : TYPE(mp_para_env_type), INTENT(IN), TARGET :: para_env
298 : TYPE(cp_blacs_env_type), INTENT(IN), TARGET :: blacs_env
299 : LOGICAL, INTENT(IN), OPTIONAL :: with_aux_fit
300 :
301 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_env_initialize'
302 :
303 : INTEGER :: handle, igr, ik, ikk, ngr, niogrp, nkp, &
304 : nkp_grp, nkp_loc, npe, unit_nr
305 : INTEGER, DIMENSION(2) :: dims, pos
306 : LOGICAL :: aux_fit
307 224 : TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_aux_env, kp_env
308 : TYPE(kpoint_env_type), POINTER :: kp
309 224 : TYPE(mp_cart_type) :: comm_cart
310 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp, para_env_kp
311 :
312 224 : CALL timeset(routineN, handle)
313 :
314 224 : IF (PRESENT(with_aux_fit)) THEN
315 166 : aux_fit = with_aux_fit
316 : ELSE
317 : aux_fit = .FALSE.
318 : END IF
319 :
320 224 : kpoint%para_env => para_env
321 224 : CALL kpoint%para_env%retain()
322 224 : kpoint%blacs_env_all => blacs_env
323 224 : CALL kpoint%blacs_env_all%retain()
324 :
325 224 : CPASSERT(.NOT. ASSOCIATED(kpoint%kp_env))
326 224 : IF (aux_fit) THEN
327 28 : CPASSERT(.NOT. ASSOCIATED(kpoint%kp_aux_env))
328 : END IF
329 :
330 224 : NULLIFY (kp_env, kp_aux_env)
331 224 : nkp = kpoint%nkp
332 224 : npe = para_env%num_pe
333 224 : IF (npe == 1) THEN
334 : ! only one process available -> owns all kpoints
335 0 : ALLOCATE (kp_env(nkp))
336 0 : DO ik = 1, nkp
337 0 : NULLIFY (kp_env(ik)%kpoint_env)
338 0 : CALL kpoint_env_create(kp_env(ik)%kpoint_env)
339 0 : kp => kp_env(ik)%kpoint_env
340 0 : kp%nkpoint = ik
341 0 : kp%wkp = kpoint%wkp(ik)
342 0 : kp%xkp(1:3) = kpoint%xkp(1:3, ik)
343 0 : kp%is_local = .TRUE.
344 : END DO
345 0 : kpoint%kp_env => kp_env
346 :
347 0 : IF (aux_fit) THEN
348 0 : ALLOCATE (kp_aux_env(nkp))
349 0 : DO ik = 1, nkp
350 0 : NULLIFY (kp_aux_env(ik)%kpoint_env)
351 0 : CALL kpoint_env_create(kp_aux_env(ik)%kpoint_env)
352 0 : kp => kp_aux_env(ik)%kpoint_env
353 0 : kp%nkpoint = ik
354 0 : kp%wkp = kpoint%wkp(ik)
355 0 : kp%xkp(1:3) = kpoint%xkp(1:3, ik)
356 0 : kp%is_local = .TRUE.
357 : END DO
358 :
359 0 : kpoint%kp_aux_env => kp_aux_env
360 : END IF
361 :
362 0 : ALLOCATE (kpoint%kp_dist(2, 1))
363 0 : kpoint%kp_dist(1, 1) = 1
364 0 : kpoint%kp_dist(2, 1) = nkp
365 0 : kpoint%kp_range(1) = 1
366 0 : kpoint%kp_range(2) = nkp
367 :
368 : ! parallel environments
369 0 : kpoint%para_env_kp => para_env
370 0 : CALL kpoint%para_env_kp%retain()
371 0 : kpoint%para_env_inter_kp => para_env
372 0 : CALL kpoint%para_env_inter_kp%retain()
373 0 : kpoint%iogrp = .TRUE.
374 0 : kpoint%nkp_groups = 1
375 : ELSE
376 224 : IF (kpoint%parallel_group_size == -1) THEN
377 : ! maximum parallelization over kpoints
378 : ! making sure that the group size divides the npe and the nkp_grp the nkp
379 : ! in the worst case, there will be no parallelism over kpoints.
380 342 : DO igr = npe, 1, -1
381 228 : IF (MOD(npe, igr) .NE. 0) CYCLE
382 228 : nkp_grp = npe/igr
383 228 : IF (MOD(nkp, nkp_grp) .NE. 0) CYCLE
384 342 : ngr = igr
385 : END DO
386 110 : ELSE IF (kpoint%parallel_group_size == 0) THEN
387 : ! no parallelization over kpoints
388 62 : ngr = npe
389 48 : ELSE IF (kpoint%parallel_group_size > 0) THEN
390 48 : ngr = MIN(kpoint%parallel_group_size, npe)
391 : ELSE
392 0 : CPASSERT(.FALSE.)
393 : END IF
394 224 : nkp_grp = npe/ngr
395 : ! processor dimensions
396 224 : dims(1) = ngr
397 224 : dims(2) = nkp_grp
398 224 : CPASSERT(MOD(nkp, nkp_grp) == 0)
399 224 : nkp_loc = nkp/nkp_grp
400 :
401 224 : IF ((dims(1)*dims(2) /= npe)) THEN
402 0 : CPABORT("Number of processors is not divisible by the kpoint group size.")
403 : END IF
404 :
405 : ! Create the subgroups, one for each k-point group and one interconnecting group
406 224 : CALL comm_cart%create(comm_old=para_env, ndims=2, dims=dims)
407 672 : pos = comm_cart%mepos_cart
408 224 : ALLOCATE (para_env_kp)
409 224 : CALL para_env_kp%from_split(comm_cart, pos(2))
410 224 : ALLOCATE (para_env_inter_kp)
411 224 : CALL para_env_inter_kp%from_split(comm_cart, pos(1))
412 224 : CALL comm_cart%free()
413 :
414 224 : niogrp = 0
415 224 : IF (para_env%is_source()) niogrp = 1
416 224 : CALL para_env_kp%sum(niogrp)
417 224 : kpoint%iogrp = (niogrp == 1)
418 :
419 : ! parallel groups
420 224 : kpoint%para_env_kp => para_env_kp
421 224 : kpoint%para_env_inter_kp => para_env_inter_kp
422 :
423 : ! distribution of kpoints
424 672 : ALLOCATE (kpoint%kp_dist(2, nkp_grp))
425 528 : DO igr = 1, nkp_grp
426 1136 : kpoint%kp_dist(1:2, igr) = get_limit(nkp, nkp_grp, igr - 1)
427 : END DO
428 : ! local kpoints
429 672 : kpoint%kp_range(1:2) = kpoint%kp_dist(1:2, para_env_inter_kp%mepos + 1)
430 :
431 2192 : ALLOCATE (kp_env(nkp_loc))
432 1744 : DO ik = 1, nkp_loc
433 1520 : NULLIFY (kp_env(ik)%kpoint_env)
434 1520 : ikk = kpoint%kp_range(1) + ik - 1
435 1520 : CALL kpoint_env_create(kp_env(ik)%kpoint_env)
436 1520 : kp => kp_env(ik)%kpoint_env
437 1520 : kp%nkpoint = ikk
438 1520 : kp%wkp = kpoint%wkp(ikk)
439 6080 : kp%xkp(1:3) = kpoint%xkp(1:3, ikk)
440 1744 : kp%is_local = (ngr == 1)
441 : END DO
442 224 : kpoint%kp_env => kp_env
443 :
444 224 : IF (aux_fit) THEN
445 216 : ALLOCATE (kp_aux_env(nkp_loc))
446 188 : DO ik = 1, nkp_loc
447 160 : NULLIFY (kp_aux_env(ik)%kpoint_env)
448 160 : ikk = kpoint%kp_range(1) + ik - 1
449 160 : CALL kpoint_env_create(kp_aux_env(ik)%kpoint_env)
450 160 : kp => kp_aux_env(ik)%kpoint_env
451 160 : kp%nkpoint = ikk
452 160 : kp%wkp = kpoint%wkp(ikk)
453 640 : kp%xkp(1:3) = kpoint%xkp(1:3, ikk)
454 188 : kp%is_local = (ngr == 1)
455 : END DO
456 28 : kpoint%kp_aux_env => kp_aux_env
457 : END IF
458 :
459 224 : unit_nr = cp_logger_get_default_io_unit()
460 :
461 224 : IF (unit_nr > 0 .AND. kpoint%verbose) THEN
462 5 : WRITE (unit_nr, *)
463 5 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Number of kpoint groups ", nkp_grp
464 5 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Size of each kpoint group", ngr
465 5 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Number of kpoints per group", nkp_loc
466 : END IF
467 224 : kpoint%nkp_groups = nkp_grp
468 :
469 : END IF
470 :
471 224 : CALL timestop(handle)
472 :
473 448 : END SUBROUTINE kpoint_env_initialize
474 :
475 : ! **************************************************************************************************
476 : !> \brief Initialize a set of MOs and density matrix for each kpoint (kpoint group)
477 : !> \param kpoint Kpoint environment
478 : !> \param mos Reference MOs (global)
479 : !> \param added_mos ...
480 : !> \param for_aux_fit ...
481 : ! **************************************************************************************************
482 252 : SUBROUTINE kpoint_initialize_mos(kpoint, mos, added_mos, for_aux_fit)
483 :
484 : TYPE(kpoint_type), POINTER :: kpoint
485 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
486 : INTEGER, INTENT(IN), OPTIONAL :: added_mos
487 : LOGICAL, OPTIONAL :: for_aux_fit
488 :
489 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_initialize_mos'
490 :
491 : INTEGER :: handle, ic, ik, is, nadd, nao, nc, &
492 : nelectron, nkp_loc, nmo, nmorig(2), &
493 : nspin
494 : LOGICAL :: aux_fit
495 : REAL(KIND=dp) :: flexible_electron_count, maxocc, n_el_f
496 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
497 252 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools
498 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
499 : TYPE(cp_fm_type), POINTER :: fmlocal
500 : TYPE(kpoint_env_type), POINTER :: kp
501 : TYPE(qs_matrix_pools_type), POINTER :: mpools
502 :
503 252 : CALL timeset(routineN, handle)
504 :
505 252 : IF (PRESENT(for_aux_fit)) THEN
506 28 : aux_fit = for_aux_fit
507 : ELSE
508 : aux_fit = .FALSE.
509 : END IF
510 :
511 252 : CPASSERT(ASSOCIATED(kpoint))
512 :
513 : IF (.TRUE. .OR. ASSOCIATED(mos(1)%mo_coeff)) THEN
514 252 : IF (aux_fit) THEN
515 28 : CPASSERT(ASSOCIATED(kpoint%kp_aux_env))
516 : END IF
517 :
518 252 : IF (PRESENT(added_mos)) THEN
519 58 : nadd = added_mos
520 : ELSE
521 : nadd = 0
522 : END IF
523 :
524 252 : IF (kpoint%use_real_wfn) THEN
525 : nc = 1
526 : ELSE
527 240 : nc = 2
528 : END IF
529 252 : nspin = SIZE(mos, 1)
530 252 : nkp_loc = kpoint%kp_range(2) - kpoint%kp_range(1) + 1
531 252 : IF (nkp_loc > 0) THEN
532 252 : IF (aux_fit) THEN
533 28 : CPASSERT(SIZE(kpoint%kp_aux_env) == nkp_loc)
534 : ELSE
535 224 : CPASSERT(SIZE(kpoint%kp_env) == nkp_loc)
536 : END IF
537 : ! allocate the mo sets, correct number of kpoints (local), real/complex, spin
538 1932 : DO ik = 1, nkp_loc
539 1680 : IF (aux_fit) THEN
540 160 : kp => kpoint%kp_aux_env(ik)%kpoint_env
541 : ELSE
542 1520 : kp => kpoint%kp_env(ik)%kpoint_env
543 : END IF
544 12490 : ALLOCATE (kp%mos(nc, nspin))
545 3860 : DO is = 1, nspin
546 : CALL get_mo_set(mos(is), nao=nao, nmo=nmo, nelectron=nelectron, &
547 1928 : n_el_f=n_el_f, maxocc=maxocc, flexible_electron_count=flexible_electron_count)
548 1928 : nmo = MIN(nao, nmo + nadd)
549 7450 : DO ic = 1, nc
550 : CALL allocate_mo_set(kp%mos(ic, is), nao, nmo, nelectron, n_el_f, maxocc, &
551 5770 : flexible_electron_count)
552 : END DO
553 : END DO
554 : END DO
555 :
556 : ! generate the blacs environment for the kpoint group
557 : ! we generate a blacs env for each kpoint group in parallel
558 : ! we assume here that the group para_env_inter_kp will connect
559 : ! equivalent parts of fm matrices, i.e. no reshuffeling of processors
560 252 : NULLIFY (blacs_env)
561 252 : IF (ASSOCIATED(kpoint%blacs_env)) THEN
562 28 : blacs_env => kpoint%blacs_env
563 : ELSE
564 224 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=kpoint%para_env_kp)
565 224 : kpoint%blacs_env => blacs_env
566 : END IF
567 :
568 : ! set possible new number of MOs
569 540 : DO is = 1, nspin
570 288 : CALL get_mo_set(mos(is), nmo=nmorig(is))
571 288 : nmo = MIN(nao, nmorig(is) + nadd)
572 540 : CALL set_mo_set(mos(is), nmo=nmo)
573 : END DO
574 : ! matrix pools for the kpoint group, information on MOs is transferred using
575 : ! generic mos structure
576 252 : NULLIFY (mpools)
577 252 : CALL mpools_create(mpools=mpools)
578 : CALL mpools_rebuild_fm_pools(mpools=mpools, mos=mos, &
579 252 : blacs_env=blacs_env, para_env=kpoint%para_env_kp)
580 :
581 252 : IF (aux_fit) THEN
582 28 : kpoint%mpools_aux_fit => mpools
583 : ELSE
584 224 : kpoint%mpools => mpools
585 : END IF
586 :
587 : ! reset old number of MOs
588 540 : DO is = 1, nspin
589 540 : CALL set_mo_set(mos(is), nmo=nmorig(is))
590 : END DO
591 :
592 : ! allocate density matrices
593 252 : CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
594 252 : ALLOCATE (fmlocal)
595 252 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
596 252 : CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
597 1932 : DO ik = 1, nkp_loc
598 1680 : IF (aux_fit) THEN
599 160 : kp => kpoint%kp_aux_env(ik)%kpoint_env
600 : ELSE
601 1520 : kp => kpoint%kp_env(ik)%kpoint_env
602 : END IF
603 : ! density matrix
604 1680 : CALL cp_fm_release(kp%pmat)
605 12490 : ALLOCATE (kp%pmat(nc, nspin))
606 3608 : DO is = 1, nspin
607 7450 : DO ic = 1, nc
608 5770 : CALL cp_fm_create(kp%pmat(ic, is), matrix_struct)
609 : END DO
610 : END DO
611 : ! energy weighted density matrix
612 1680 : CALL cp_fm_release(kp%wmat)
613 10810 : ALLOCATE (kp%wmat(nc, nspin))
614 3860 : DO is = 1, nspin
615 7450 : DO ic = 1, nc
616 5770 : CALL cp_fm_create(kp%wmat(ic, is), matrix_struct)
617 : END DO
618 : END DO
619 : END DO
620 252 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
621 252 : DEALLOCATE (fmlocal)
622 :
623 : END IF
624 :
625 : END IF
626 :
627 252 : CALL timestop(handle)
628 :
629 252 : END SUBROUTINE kpoint_initialize_mos
630 :
631 : ! **************************************************************************************************
632 : !> \brief ...
633 : !> \param kpoint ...
634 : ! **************************************************************************************************
635 58 : SUBROUTINE kpoint_initialize_mo_set(kpoint)
636 : TYPE(kpoint_type), POINTER :: kpoint
637 :
638 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_initialize_mo_set'
639 :
640 : INTEGER :: handle, ic, ik, ikk, ispin
641 58 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
642 : TYPE(cp_fm_type), POINTER :: mo_coeff
643 58 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: moskp
644 :
645 58 : CALL timeset(routineN, handle)
646 :
647 404 : DO ik = 1, SIZE(kpoint%kp_env)
648 346 : CALL mpools_get(kpoint%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
649 346 : moskp => kpoint%kp_env(ik)%kpoint_env%mos
650 346 : ikk = kpoint%kp_range(1) + ik - 1
651 346 : CPASSERT(ASSOCIATED(moskp))
652 784 : DO ispin = 1, SIZE(moskp, 2)
653 1486 : DO ic = 1, SIZE(moskp, 1)
654 760 : CALL get_mo_set(moskp(ic, ispin), mo_coeff=mo_coeff)
655 1140 : IF (.NOT. ASSOCIATED(mo_coeff)) THEN
656 : CALL init_mo_set(moskp(ic, ispin), &
657 760 : fm_pool=ao_mo_fm_pools(ispin)%pool, name="kpoints")
658 : END IF
659 : END DO
660 : END DO
661 : END DO
662 :
663 58 : CALL timestop(handle)
664 :
665 58 : END SUBROUTINE kpoint_initialize_mo_set
666 :
667 : ! **************************************************************************************************
668 : !> \brief Generates the mapping of cell indices and linear RS index
669 : !> CELL (0,0,0) is always mapped to index 1
670 : !> \param kpoint Kpoint environment
671 : !> \param sab_nl Defining neighbour list
672 : !> \param para_env Parallel environment
673 : !> \param dft_control ...
674 : ! **************************************************************************************************
675 988 : SUBROUTINE kpoint_init_cell_index(kpoint, sab_nl, para_env, dft_control)
676 :
677 : TYPE(kpoint_type), POINTER :: kpoint
678 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
679 : POINTER :: sab_nl
680 : TYPE(mp_para_env_type), POINTER :: para_env
681 : TYPE(dft_control_type), POINTER :: dft_control
682 :
683 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_init_cell_index'
684 :
685 : INTEGER :: handle, i1, i2, i3, ic, icount, it, &
686 : ncount
687 : INTEGER, DIMENSION(3) :: cell, itm
688 988 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell, list
689 988 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index, cti
690 : LOGICAL :: new
691 : TYPE(neighbor_list_iterator_p_type), &
692 988 : DIMENSION(:), POINTER :: nl_iterator
693 :
694 988 : NULLIFY (cell_to_index, index_to_cell)
695 :
696 988 : CALL timeset(routineN, handle)
697 :
698 988 : CPASSERT(ASSOCIATED(kpoint))
699 :
700 988 : ALLOCATE (list(3, 125))
701 494988 : list = 0
702 988 : icount = 1
703 :
704 988 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
705 607386 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
706 606398 : CALL get_iterator_info(nl_iterator, cell=cell)
707 :
708 606398 : new = .TRUE.
709 35871597 : DO ic = 1, icount
710 35790766 : IF (cell(1) == list(1, ic) .AND. cell(2) == list(2, ic) .AND. &
711 80831 : cell(3) == list(3, ic)) THEN
712 : new = .FALSE.
713 : EXIT
714 : END IF
715 : END DO
716 607386 : IF (new) THEN
717 80831 : icount = icount + 1
718 80831 : IF (icount > SIZE(list, 2)) THEN
719 273 : CALL reallocate(list, 1, 3, 1, 2*SIZE(list, 2))
720 : END IF
721 323324 : list(1:3, icount) = cell(1:3)
722 : END IF
723 :
724 : END DO
725 988 : CALL neighbor_list_iterator_release(nl_iterator)
726 :
727 82807 : itm(1) = MAXVAL(ABS(list(1, 1:icount)))
728 82807 : itm(2) = MAXVAL(ABS(list(2, 1:icount)))
729 82807 : itm(3) = MAXVAL(ABS(list(3, 1:icount)))
730 988 : CALL para_env%max(itm)
731 3952 : it = MAXVAL(itm(1:3))
732 988 : IF (ASSOCIATED(kpoint%cell_to_index)) THEN
733 984 : DEALLOCATE (kpoint%cell_to_index)
734 : END IF
735 4940 : ALLOCATE (kpoint%cell_to_index(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
736 988 : cell_to_index => kpoint%cell_to_index
737 988 : cti => cell_to_index
738 199900 : cti(:, :, :) = 0
739 82807 : DO ic = 1, icount
740 81819 : i1 = list(1, ic)
741 81819 : i2 = list(2, ic)
742 81819 : i3 = list(3, ic)
743 82807 : cti(i1, i2, i3) = ic
744 : END DO
745 398812 : CALL para_env%sum(cti)
746 988 : ncount = 0
747 5748 : DO i1 = -itm(1), itm(1)
748 34832 : DO i2 = -itm(2), itm(2)
749 200472 : DO i3 = -itm(3), itm(3)
750 195712 : IF (cti(i1, i2, i3) == 0) THEN
751 76028 : cti(i1, i2, i3) = 1000000
752 : ELSE
753 90600 : ncount = ncount + 1
754 90600 : cti(i1, i2, i3) = (ABS(i1) + ABS(i2) + ABS(i3))*1000 + ABS(i3)*100 + ABS(i2)*10 + ABS(i1)
755 90600 : cti(i1, i2, i3) = cti(i1, i2, i3) + (i1 + i2 + i3)
756 : END IF
757 : END DO
758 : END DO
759 : END DO
760 :
761 988 : IF (ASSOCIATED(kpoint%index_to_cell)) THEN
762 988 : DEALLOCATE (kpoint%index_to_cell)
763 : END IF
764 2964 : ALLOCATE (kpoint%index_to_cell(3, ncount))
765 988 : index_to_cell => kpoint%index_to_cell
766 91588 : DO ic = 1, ncount
767 90600 : cell = MINLOC(cti)
768 90600 : i1 = cell(1) - 1 - itm(1)
769 90600 : i2 = cell(2) - 1 - itm(2)
770 90600 : i3 = cell(3) - 1 - itm(3)
771 90600 : cti(i1, i2, i3) = 1000000
772 90600 : index_to_cell(1, ic) = i1
773 90600 : index_to_cell(2, ic) = i2
774 91588 : index_to_cell(3, ic) = i3
775 : END DO
776 199900 : cti(:, :, :) = 0
777 91588 : DO ic = 1, ncount
778 90600 : i1 = index_to_cell(1, ic)
779 90600 : i2 = index_to_cell(2, ic)
780 90600 : i3 = index_to_cell(3, ic)
781 91588 : cti(i1, i2, i3) = ic
782 : END DO
783 :
784 : ! keep pointer to this neighborlist
785 988 : kpoint%sab_nl => sab_nl
786 :
787 : ! set number of images
788 988 : dft_control%nimages = SIZE(index_to_cell, 2)
789 :
790 988 : DEALLOCATE (list)
791 :
792 988 : CALL timestop(handle)
793 988 : END SUBROUTINE kpoint_init_cell_index
794 :
795 : ! **************************************************************************************************
796 : !> \brief Transformation of real space matrices to a kpoint
797 : !> \param rmatrix Real part of kpoint matrix
798 : !> \param cmatrix Complex part of kpoint matrix (optional)
799 : !> \param rsmat Real space matrices
800 : !> \param ispin Spin index
801 : !> \param xkp Kpoint coordinates
802 : !> \param cell_to_index mapping of cell indices to RS index
803 : !> \param sab_nl Defining neighbor list
804 : !> \param is_complex Matrix to be transformed is imaginary
805 : !> \param rs_sign Matrix to be transformed is csaled by rs_sign
806 : ! **************************************************************************************************
807 100260 : SUBROUTINE rskp_transform(rmatrix, cmatrix, rsmat, ispin, &
808 : xkp, cell_to_index, sab_nl, is_complex, rs_sign)
809 :
810 : TYPE(dbcsr_type) :: rmatrix
811 : TYPE(dbcsr_type), OPTIONAL :: cmatrix
812 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
813 : INTEGER, INTENT(IN) :: ispin
814 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
815 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
816 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
817 : POINTER :: sab_nl
818 : LOGICAL, INTENT(IN), OPTIONAL :: is_complex
819 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: rs_sign
820 :
821 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rskp_transform'
822 :
823 : INTEGER :: handle, iatom, ic, icol, irow, jatom, &
824 : nimg
825 : INTEGER, DIMENSION(3) :: cell
826 : LOGICAL :: do_symmetric, found, my_complex, &
827 : wfn_real_only
828 : REAL(KIND=dp) :: arg, coskl, fsign, fsym, sinkl
829 50130 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, rblock, rsblock
830 : TYPE(neighbor_list_iterator_p_type), &
831 50130 : DIMENSION(:), POINTER :: nl_iterator
832 :
833 50130 : CALL timeset(routineN, handle)
834 :
835 50130 : my_complex = .FALSE.
836 50130 : IF (PRESENT(is_complex)) my_complex = is_complex
837 :
838 50130 : fsign = 1.0_dp
839 50130 : IF (PRESENT(rs_sign)) fsign = rs_sign
840 :
841 50130 : wfn_real_only = .TRUE.
842 50130 : IF (PRESENT(cmatrix)) wfn_real_only = .FALSE.
843 :
844 50130 : nimg = SIZE(rsmat, 2)
845 :
846 50130 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
847 :
848 50130 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
849 22192942 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
850 22142812 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
851 :
852 : ! fsym = +- 1 is due to real space matrices being non-symmetric (although in a symmtric type)
853 : ! with the link S_mu^0,nu^b = S_nu^0,mu^-b, and the KP matrices beeing Hermitian
854 22142812 : fsym = 1.0_dp
855 22142812 : irow = iatom
856 22142812 : icol = jatom
857 22142812 : IF (do_symmetric .AND. (iatom > jatom)) THEN
858 9377774 : irow = jatom
859 9377774 : icol = iatom
860 9377774 : fsym = -1.0_dp
861 : END IF
862 :
863 22142812 : ic = cell_to_index(cell(1), cell(2), cell(3))
864 22142812 : IF (ic < 1 .OR. ic > nimg) CYCLE
865 :
866 22142068 : arg = REAL(cell(1), dp)*xkp(1) + REAL(cell(2), dp)*xkp(2) + REAL(cell(3), dp)*xkp(3)
867 22142068 : IF (my_complex) THEN
868 0 : coskl = fsign*fsym*COS(twopi*arg)
869 0 : sinkl = fsign*SIN(twopi*arg)
870 : ELSE
871 22142068 : coskl = fsign*COS(twopi*arg)
872 22142068 : sinkl = fsign*fsym*SIN(twopi*arg)
873 : END IF
874 :
875 : CALL dbcsr_get_block_p(matrix=rsmat(ispin, ic)%matrix, row=irow, col=icol, &
876 22142068 : block=rsblock, found=found)
877 22142068 : IF (.NOT. found) CYCLE
878 :
879 22192198 : IF (wfn_real_only) THEN
880 : CALL dbcsr_get_block_p(matrix=rmatrix, row=irow, col=icol, &
881 303864 : block=rblock, found=found)
882 303864 : IF (.NOT. found) CYCLE
883 153367320 : rblock = rblock + coskl*rsblock
884 : ELSE
885 : CALL dbcsr_get_block_p(matrix=rmatrix, row=irow, col=icol, &
886 21838204 : block=rblock, found=found)
887 21838204 : IF (.NOT. found) CYCLE
888 : CALL dbcsr_get_block_p(matrix=cmatrix, row=irow, col=icol, &
889 21838204 : block=cblock, found=found)
890 21838204 : IF (.NOT. found) CYCLE
891 4280031184 : rblock = rblock + coskl*rsblock
892 4280031184 : cblock = cblock + sinkl*rsblock
893 : END IF
894 :
895 : END DO
896 50130 : CALL neighbor_list_iterator_release(nl_iterator)
897 :
898 50130 : CALL timestop(handle)
899 :
900 50130 : END SUBROUTINE rskp_transform
901 :
902 : ! **************************************************************************************************
903 : !> \brief Given the eigenvalues of all kpoints, calculates the occupation numbers
904 : !> \param kpoint Kpoint environment
905 : !> \param smear Smearing information
906 : ! **************************************************************************************************
907 5464 : SUBROUTINE kpoint_set_mo_occupation(kpoint, smear)
908 :
909 : TYPE(kpoint_type), POINTER :: kpoint
910 : TYPE(smear_type), POINTER :: smear
911 :
912 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_set_mo_occupation'
913 :
914 : INTEGER :: handle, ik, ikpgr, ispin, kplocal, nb, &
915 : ne_a, ne_b, nelectron, nkp, nmo, nspin
916 : INTEGER, DIMENSION(2) :: kp_range
917 : REAL(KIND=dp) :: kTS, mu, mus(2), nel
918 5464 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: weig, wocc
919 5464 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation, wkp
920 : TYPE(kpoint_env_type), POINTER :: kp
921 : TYPE(mo_set_type), POINTER :: mo_set
922 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
923 :
924 5464 : CALL timeset(routineN, handle)
925 :
926 : ! first collect all the eigenvalues
927 5464 : CALL get_kpoint_info(kpoint, nkp=nkp)
928 5464 : kp => kpoint%kp_env(1)%kpoint_env
929 5464 : nspin = SIZE(kp%mos, 2)
930 5464 : mo_set => kp%mos(1, 1)
931 5464 : CALL get_mo_set(mo_set, nmo=nmo, nelectron=nelectron)
932 5464 : ne_a = nelectron
933 5464 : IF (nspin == 2) THEN
934 542 : CALL get_mo_set(kp%mos(1, 2), nmo=nb, nelectron=ne_b)
935 542 : CPASSERT(nmo == nb)
936 : END IF
937 43712 : ALLOCATE (weig(nmo, nkp, nspin), wocc(nmo, nkp, nspin))
938 462772 : weig = 0.0_dp
939 462772 : wocc = 0.0_dp
940 5464 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
941 5464 : kplocal = kp_range(2) - kp_range(1) + 1
942 19024 : DO ikpgr = 1, kplocal
943 13560 : ik = kp_range(1) + ikpgr - 1
944 13560 : kp => kpoint%kp_env(ikpgr)%kpoint_env
945 34452 : DO ispin = 1, nspin
946 15428 : mo_set => kp%mos(1, ispin)
947 15428 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
948 343256 : weig(1:nmo, ik, ispin) = eigenvalues(1:nmo)
949 : END DO
950 : END DO
951 5464 : CALL get_kpoint_info(kpoint, para_env_inter_kp=para_env_inter_kp)
952 5464 : CALL para_env_inter_kp%sum(weig)
953 :
954 5464 : CALL get_kpoint_info(kpoint, wkp=wkp)
955 5464 : IF (smear%do_smear) THEN
956 : ! finite electronic temperature
957 1412 : SELECT CASE (smear%method)
958 : CASE (smear_fermi_dirac)
959 706 : IF (nspin == 1) THEN
960 694 : nel = REAL(nelectron, KIND=dp)
961 : CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
962 694 : smear%electronic_temperature, 2.0_dp)
963 12 : ELSE IF (smear%fixed_mag_mom > 0.0_dp) THEN
964 2 : nel = REAL(ne_a, KIND=dp)
965 : CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
966 2 : smear%electronic_temperature, 1.0_dp)
967 2 : nel = REAL(ne_b, KIND=dp)
968 : CALL Fermikp(wocc(:, :, 2), mus(2), kTS, weig(:, :, 2), nel, wkp, &
969 2 : smear%electronic_temperature, 1.0_dp)
970 : ELSE
971 10 : nel = REAL(ne_a, KIND=dp) + REAL(ne_b, KIND=dp)
972 : CALL Fermikp2(wocc(:, :, :), mu, kTS, weig(:, :, :), nel, wkp, &
973 10 : smear%electronic_temperature)
974 10 : kTS = kTS/2._dp
975 30 : mus(1:2) = mu
976 : END IF
977 : CASE DEFAULT
978 706 : CPABORT("kpoints: Selected smearing not (yet) supported")
979 : END SELECT
980 : ELSE
981 : ! fixed occupations (2/1)
982 4758 : IF (nspin == 1) THEN
983 4228 : nel = REAL(nelectron, KIND=dp)
984 4228 : CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, 0.0_dp, 2.0_dp)
985 : ELSE
986 530 : nel = REAL(ne_a, KIND=dp)
987 530 : CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, 0.0_dp, 1.0_dp)
988 530 : nel = REAL(ne_b, KIND=dp)
989 530 : CALL Fermikp(wocc(:, :, 2), mus(2), kTS, weig(:, :, 2), nel, wkp, 0.0_dp, 1.0_dp)
990 : END IF
991 : END IF
992 19024 : DO ikpgr = 1, kplocal
993 13560 : ik = kp_range(1) + ikpgr - 1
994 13560 : kp => kpoint%kp_env(ikpgr)%kpoint_env
995 34452 : DO ispin = 1, nspin
996 15428 : mo_set => kp%mos(1, ispin)
997 15428 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues, occupation_numbers=occupation)
998 329696 : eigenvalues(1:nmo) = weig(1:nmo, ik, ispin)
999 329696 : occupation(1:nmo) = wocc(1:nmo, ik, ispin)
1000 15428 : mo_set%kTS = kTS
1001 28988 : mo_set%mu = mus(ispin)
1002 : END DO
1003 : END DO
1004 :
1005 5464 : DEALLOCATE (weig, wocc)
1006 :
1007 5464 : CALL timestop(handle)
1008 :
1009 10928 : END SUBROUTINE kpoint_set_mo_occupation
1010 :
1011 : ! **************************************************************************************************
1012 : !> \brief Calculate kpoint density matrices (rho(k), owned by kpoint groups)
1013 : !> \param kpoint kpoint environment
1014 : !> \param energy_weighted calculate energy weighted density matrix
1015 : !> \param for_aux_fit ...
1016 : ! **************************************************************************************************
1017 17082 : SUBROUTINE kpoint_density_matrices(kpoint, energy_weighted, for_aux_fit)
1018 :
1019 : TYPE(kpoint_type), POINTER :: kpoint
1020 : LOGICAL, OPTIONAL :: energy_weighted, for_aux_fit
1021 :
1022 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_density_matrices'
1023 :
1024 : INTEGER :: handle, ikpgr, ispin, kplocal, nao, nmo, &
1025 : nspin
1026 : INTEGER, DIMENSION(2) :: kp_range
1027 : LOGICAL :: aux_fit, wtype
1028 5694 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation
1029 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1030 : TYPE(cp_fm_type) :: fwork
1031 : TYPE(cp_fm_type), POINTER :: cpmat, pmat, rpmat
1032 : TYPE(kpoint_env_type), POINTER :: kp
1033 : TYPE(mo_set_type), POINTER :: mo_set
1034 :
1035 5694 : CALL timeset(routineN, handle)
1036 :
1037 5694 : IF (PRESENT(energy_weighted)) THEN
1038 152 : wtype = energy_weighted
1039 : ELSE
1040 : ! default is normal density matrix
1041 : wtype = .FALSE.
1042 : END IF
1043 :
1044 5694 : IF (PRESENT(for_aux_fit)) THEN
1045 62 : aux_fit = for_aux_fit
1046 : ELSE
1047 : aux_fit = .FALSE.
1048 : END IF
1049 :
1050 62 : IF (aux_fit) THEN
1051 62 : CPASSERT(ASSOCIATED(kpoint%kp_aux_env))
1052 : END IF
1053 :
1054 : ! work matrix
1055 5694 : IF (aux_fit) THEN
1056 62 : mo_set => kpoint%kp_aux_env(1)%kpoint_env%mos(1, 1)
1057 : ELSE
1058 5632 : mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)
1059 : END IF
1060 5694 : CALL get_mo_set(mo_set, nao=nao, nmo=nmo)
1061 5694 : CALL cp_fm_get_info(mo_set%mo_coeff, matrix_struct=matrix_struct)
1062 5694 : CALL cp_fm_create(fwork, matrix_struct)
1063 :
1064 5694 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
1065 5694 : kplocal = kp_range(2) - kp_range(1) + 1
1066 20412 : DO ikpgr = 1, kplocal
1067 14718 : IF (aux_fit) THEN
1068 352 : kp => kpoint%kp_aux_env(ikpgr)%kpoint_env
1069 : ELSE
1070 14366 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1071 : END IF
1072 14718 : nspin = SIZE(kp%mos, 2)
1073 37246 : DO ispin = 1, nspin
1074 16834 : mo_set => kp%mos(1, ispin)
1075 16834 : IF (wtype) THEN
1076 826 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
1077 : END IF
1078 31552 : IF (kpoint%use_real_wfn) THEN
1079 72 : IF (wtype) THEN
1080 10 : pmat => kp%wmat(1, ispin)
1081 : ELSE
1082 62 : pmat => kp%pmat(1, ispin)
1083 : END IF
1084 72 : CALL get_mo_set(mo_set, occupation_numbers=occupation)
1085 72 : CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
1086 72 : CALL cp_fm_column_scale(fwork, occupation)
1087 72 : IF (wtype) THEN
1088 10 : CALL cp_fm_column_scale(fwork, eigenvalues)
1089 : END IF
1090 72 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, pmat)
1091 : ELSE
1092 16762 : IF (wtype) THEN
1093 816 : rpmat => kp%wmat(1, ispin)
1094 816 : cpmat => kp%wmat(2, ispin)
1095 : ELSE
1096 15946 : rpmat => kp%pmat(1, ispin)
1097 15946 : cpmat => kp%pmat(2, ispin)
1098 : END IF
1099 16762 : CALL get_mo_set(mo_set, occupation_numbers=occupation)
1100 16762 : CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
1101 16762 : CALL cp_fm_column_scale(fwork, occupation)
1102 16762 : IF (wtype) THEN
1103 816 : CALL cp_fm_column_scale(fwork, eigenvalues)
1104 : END IF
1105 : ! Re(c)*Re(c)
1106 16762 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, rpmat)
1107 16762 : mo_set => kp%mos(2, ispin)
1108 : ! Im(c)*Re(c)
1109 16762 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, cpmat)
1110 : ! Re(c)*Im(c)
1111 16762 : CALL parallel_gemm("N", "T", nao, nao, nmo, -1.0_dp, fwork, mo_set%mo_coeff, 1.0_dp, cpmat)
1112 16762 : CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
1113 16762 : CALL cp_fm_column_scale(fwork, occupation)
1114 16762 : IF (wtype) THEN
1115 816 : CALL cp_fm_column_scale(fwork, eigenvalues)
1116 : END IF
1117 : ! Im(c)*Im(c)
1118 16762 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 1.0_dp, rpmat)
1119 : END IF
1120 : END DO
1121 : END DO
1122 :
1123 5694 : CALL cp_fm_release(fwork)
1124 :
1125 5694 : CALL timestop(handle)
1126 :
1127 5694 : END SUBROUTINE kpoint_density_matrices
1128 :
1129 : ! **************************************************************************************************
1130 : !> \brief generate real space density matrices in DBCSR format
1131 : !> \param kpoint Kpoint environment
1132 : !> \param denmat Real space (DBCSR) density matrices
1133 : !> \param wtype True = energy weighted density matrix
1134 : !> False = normal density matrix
1135 : !> \param tempmat DBCSR matrix to be used as template
1136 : !> \param sab_nl ...
1137 : !> \param fmwork FM work matrices (kpoint group)
1138 : !> \param for_aux_fit ...
1139 : !> \param pmat_ext ...
1140 : ! **************************************************************************************************
1141 5860 : SUBROUTINE kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwork, for_aux_fit, pmat_ext)
1142 :
1143 : TYPE(kpoint_type), POINTER :: kpoint
1144 : TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
1145 : LOGICAL, INTENT(IN) :: wtype
1146 : TYPE(dbcsr_type), POINTER :: tempmat
1147 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1148 : POINTER :: sab_nl
1149 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fmwork
1150 : LOGICAL, OPTIONAL :: for_aux_fit
1151 : TYPE(cp_fm_type), DIMENSION(:, :, :), INTENT(IN), &
1152 : OPTIONAL :: pmat_ext
1153 :
1154 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_density_transform'
1155 :
1156 : INTEGER :: handle, ic, ik, ikk, indx, ir, ira, is, &
1157 : ispin, jr, nc, nimg, nkp, nspin
1158 5860 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1159 : LOGICAL :: aux_fit, do_ext, do_symmetric, my_kpgrp, &
1160 : real_only
1161 : REAL(KIND=dp) :: wkpx
1162 5860 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
1163 5860 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1164 5860 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:) :: info
1165 : TYPE(cp_fm_type) :: fmdummy
1166 : TYPE(dbcsr_type), POINTER :: cpmat, rpmat, scpmat, srpmat
1167 5860 : TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kind_rot
1168 : TYPE(kpoint_env_type), POINTER :: kp
1169 : TYPE(kpoint_sym_type), POINTER :: kpsym
1170 : TYPE(mp_para_env_type), POINTER :: para_env
1171 :
1172 5860 : CALL timeset(routineN, handle)
1173 :
1174 5860 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1175 :
1176 5860 : IF (PRESENT(for_aux_fit)) THEN
1177 228 : aux_fit = for_aux_fit
1178 : ELSE
1179 : aux_fit = .FALSE.
1180 : END IF
1181 :
1182 5860 : do_ext = .FALSE.
1183 5860 : IF (PRESENT(pmat_ext)) do_ext = .TRUE.
1184 :
1185 5860 : IF (aux_fit) THEN
1186 138 : CPASSERT(ASSOCIATED(kpoint%kp_aux_env))
1187 : END IF
1188 :
1189 : ! work storage
1190 5860 : ALLOCATE (rpmat)
1191 : CALL dbcsr_create(rpmat, template=tempmat, &
1192 5908 : matrix_type=MERGE(dbcsr_type_symmetric, dbcsr_type_no_symmetry, do_symmetric))
1193 5860 : CALL cp_dbcsr_alloc_block_from_nbl(rpmat, sab_nl)
1194 5860 : CALL dbcsr_set(rpmat, 0.0_dp)
1195 5860 : ALLOCATE (cpmat)
1196 : CALL dbcsr_create(cpmat, template=tempmat, &
1197 5908 : matrix_type=MERGE(dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, do_symmetric))
1198 5860 : CALL cp_dbcsr_alloc_block_from_nbl(cpmat, sab_nl)
1199 5860 : CALL dbcsr_set(cpmat, 0.0_dp)
1200 5860 : IF (.NOT. kpoint%full_grid) THEN
1201 4784 : ALLOCATE (srpmat)
1202 4784 : CALL dbcsr_create(srpmat, template=rpmat)
1203 4784 : CALL cp_dbcsr_alloc_block_from_nbl(srpmat, sab_nl)
1204 4784 : CALL dbcsr_set(srpmat, 0.0_dp)
1205 4784 : ALLOCATE (scpmat)
1206 4784 : CALL dbcsr_create(scpmat, template=cpmat)
1207 4784 : CALL cp_dbcsr_alloc_block_from_nbl(scpmat, sab_nl)
1208 4784 : CALL dbcsr_set(scpmat, 0.0_dp)
1209 : END IF
1210 :
1211 : CALL get_kpoint_info(kpoint, nkp=nkp, xkp=xkp, wkp=wkp, &
1212 5860 : cell_to_index=cell_to_index)
1213 : ! initialize real space density matrices
1214 5860 : IF (aux_fit) THEN
1215 138 : kp => kpoint%kp_aux_env(1)%kpoint_env
1216 : ELSE
1217 5722 : kp => kpoint%kp_env(1)%kpoint_env
1218 : END IF
1219 5860 : nspin = SIZE(kp%mos, 2)
1220 5860 : nc = SIZE(kp%mos, 1)
1221 5860 : nimg = SIZE(denmat, 2)
1222 5860 : real_only = (nc == 1)
1223 :
1224 5860 : para_env => kpoint%blacs_env_all%para_env
1225 118840 : ALLOCATE (info(nspin*nkp*nc))
1226 :
1227 : ! Start all the communication
1228 5860 : indx = 0
1229 12342 : DO ispin = 1, nspin
1230 535528 : DO ic = 1, nimg
1231 535528 : CALL dbcsr_set(denmat(ispin, ic)%matrix, 0.0_dp)
1232 : END DO
1233 : !
1234 39568 : DO ik = 1, nkp
1235 27226 : my_kpgrp = (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2))
1236 : IF (my_kpgrp) THEN
1237 18050 : ikk = ik - kpoint%kp_range(1) + 1
1238 18050 : IF (aux_fit) THEN
1239 1008 : kp => kpoint%kp_aux_env(ikk)%kpoint_env
1240 : ELSE
1241 17042 : kp => kpoint%kp_env(ikk)%kpoint_env
1242 : END IF
1243 : ELSE
1244 : NULLIFY (kp)
1245 : END IF
1246 : ! collect this density matrix on all processors
1247 27226 : CPASSERT(SIZE(fmwork) >= nc)
1248 :
1249 33708 : IF (my_kpgrp) THEN
1250 54078 : DO ic = 1, nc
1251 36028 : indx = indx + 1
1252 54078 : IF (do_ext) THEN
1253 2016 : CALL cp_fm_start_copy_general(pmat_ext(ikk, ic, ispin), fmwork(ic), para_env, info(indx))
1254 : ELSE
1255 34012 : IF (wtype) THEN
1256 1642 : CALL cp_fm_start_copy_general(kp%wmat(ic, ispin), fmwork(ic), para_env, info(indx))
1257 : ELSE
1258 32370 : CALL cp_fm_start_copy_general(kp%pmat(ic, ispin), fmwork(ic), para_env, info(indx))
1259 : END IF
1260 : END IF
1261 : END DO
1262 : ELSE
1263 27528 : DO ic = 1, nc
1264 18352 : indx = indx + 1
1265 27528 : CALL cp_fm_start_copy_general(fmdummy, fmwork(ic), para_env, info(indx))
1266 : END DO
1267 : END IF
1268 : END DO
1269 : END DO
1270 :
1271 : ! Finish communication and transform the received matrices
1272 5860 : indx = 0
1273 12342 : DO ispin = 1, nspin
1274 39568 : DO ik = 1, nkp
1275 81606 : DO ic = 1, nc
1276 54380 : indx = indx + 1
1277 81606 : CALL cp_fm_finish_copy_general(fmwork(ic), info(indx))
1278 : END DO
1279 :
1280 : ! reduce to dbcsr storage
1281 27226 : IF (real_only) THEN
1282 72 : CALL copy_fm_to_dbcsr(fmwork(1), rpmat, keep_sparsity=.TRUE.)
1283 : ELSE
1284 27154 : CALL copy_fm_to_dbcsr(fmwork(1), rpmat, keep_sparsity=.TRUE.)
1285 27154 : CALL copy_fm_to_dbcsr(fmwork(2), cpmat, keep_sparsity=.TRUE.)
1286 : END IF
1287 :
1288 : ! symmetrization
1289 27226 : kpsym => kpoint%kp_sym(ik)%kpoint_sym
1290 27226 : CPASSERT(ASSOCIATED(kpsym))
1291 :
1292 33708 : IF (kpsym%apply_symmetry) THEN
1293 0 : wkpx = wkp(ik)/REAL(kpsym%nwght, KIND=dp)
1294 0 : DO is = 1, kpsym%nwght
1295 0 : ir = ABS(kpsym%rotp(is))
1296 0 : ira = 0
1297 0 : DO jr = 1, SIZE(kpoint%ibrot)
1298 0 : IF (ir == kpoint%ibrot(jr)) ira = jr
1299 : END DO
1300 0 : CPASSERT(ira > 0)
1301 0 : kind_rot => kpoint%kind_rotmat(ira, :)
1302 0 : IF (real_only) THEN
1303 : CALL symtrans(srpmat, rpmat, kind_rot, kpsym%rot(1:3, 1:3, is), &
1304 0 : kpsym%f0(:, is), kpoint%atype, symmetric=.TRUE.)
1305 : ELSE
1306 : CALL symtrans(srpmat, rpmat, kind_rot, kpsym%rot(1:3, 1:3, is), &
1307 0 : kpsym%f0(:, is), kpoint%atype, symmetric=.TRUE.)
1308 : CALL symtrans(scpmat, cpmat, kind_rot, kpsym%rot(1:3, 1:3, is), &
1309 0 : kpsym%f0(:, is), kpoint%atype, antisymmetric=.TRUE.)
1310 : END IF
1311 : CALL transform_dmat(denmat, srpmat, scpmat, ispin, real_only, sab_nl, &
1312 0 : cell_to_index, kpsym%xkp(1:3, is), wkpx)
1313 : END DO
1314 : ELSE
1315 : ! transformation
1316 : CALL transform_dmat(denmat, rpmat, cpmat, ispin, real_only, sab_nl, &
1317 27226 : cell_to_index, xkp(1:3, ik), wkp(ik))
1318 : END IF
1319 : END DO
1320 : END DO
1321 :
1322 : ! Clean up communication
1323 5860 : indx = 0
1324 12342 : DO ispin = 1, nspin
1325 39568 : DO ik = 1, nkp
1326 27226 : my_kpgrp = (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2))
1327 6482 : IF (my_kpgrp) THEN
1328 18050 : ikk = ik - kpoint%kp_range(1) + 1
1329 : IF (aux_fit) THEN
1330 18050 : kp => kpoint%kp_aux_env(ikk)%kpoint_env
1331 : ELSE
1332 18050 : kp => kpoint%kp_env(ikk)%kpoint_env
1333 : END IF
1334 :
1335 54078 : DO ic = 1, nc
1336 36028 : indx = indx + 1
1337 54078 : CALL cp_fm_cleanup_copy_general(info(indx))
1338 : END DO
1339 : ELSE
1340 : ! calls with dummy arguments, so not included
1341 : ! therefore just increment counter by trip count
1342 9176 : indx = indx + nc
1343 : END IF
1344 : END DO
1345 : END DO
1346 :
1347 : ! All done
1348 60240 : DEALLOCATE (info)
1349 :
1350 5860 : CALL dbcsr_deallocate_matrix(rpmat)
1351 5860 : CALL dbcsr_deallocate_matrix(cpmat)
1352 5860 : IF (.NOT. kpoint%full_grid) THEN
1353 4784 : CALL dbcsr_deallocate_matrix(srpmat)
1354 4784 : CALL dbcsr_deallocate_matrix(scpmat)
1355 : END IF
1356 :
1357 5860 : CALL timestop(handle)
1358 :
1359 5860 : END SUBROUTINE kpoint_density_transform
1360 :
1361 : ! **************************************************************************************************
1362 : !> \brief real space density matrices in DBCSR format
1363 : !> \param denmat Real space (DBCSR) density matrix
1364 : !> \param rpmat ...
1365 : !> \param cpmat ...
1366 : !> \param ispin ...
1367 : !> \param real_only ...
1368 : !> \param sab_nl ...
1369 : !> \param cell_to_index ...
1370 : !> \param xkp ...
1371 : !> \param wkp ...
1372 : ! **************************************************************************************************
1373 27226 : SUBROUTINE transform_dmat(denmat, rpmat, cpmat, ispin, real_only, sab_nl, cell_to_index, xkp, wkp)
1374 :
1375 : TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
1376 : TYPE(dbcsr_type), POINTER :: rpmat, cpmat
1377 : INTEGER, INTENT(IN) :: ispin
1378 : LOGICAL, INTENT(IN) :: real_only
1379 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1380 : POINTER :: sab_nl
1381 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1382 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
1383 : REAL(KIND=dp), INTENT(IN) :: wkp
1384 :
1385 : CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_dmat'
1386 :
1387 : INTEGER :: handle, iatom, icell, icol, irow, jatom, &
1388 : nimg
1389 : INTEGER, DIMENSION(3) :: cell
1390 : LOGICAL :: do_symmetric, found
1391 : REAL(KIND=dp) :: arg, coskl, fc, sinkl
1392 27226 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, dblock, rblock
1393 : TYPE(neighbor_list_iterator_p_type), &
1394 27226 : DIMENSION(:), POINTER :: nl_iterator
1395 :
1396 27226 : CALL timeset(routineN, handle)
1397 :
1398 27226 : nimg = SIZE(denmat, 2)
1399 :
1400 : ! transformation
1401 27226 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1402 27226 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
1403 11833608 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1404 11806382 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
1405 :
1406 : !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
1407 : !Therefore, we have: S(R) = sum_k Re(S(k))*cos(k*R) -i^2*Im(S(k))*sin(k*R)
1408 : ! = sum_k Re(S(k))*cos(k*R) + Im(S(k))*sin(k*R)
1409 : !fc = +- 1 is due to the usual non-symmetric real-space matrices stored as symmetric ones
1410 :
1411 11806382 : irow = iatom
1412 11806382 : icol = jatom
1413 11806382 : fc = 1.0_dp
1414 11806382 : IF (do_symmetric .AND. iatom > jatom) THEN
1415 4982009 : irow = jatom
1416 4982009 : icol = iatom
1417 4982009 : fc = -1.0_dp
1418 : END IF
1419 :
1420 11806382 : icell = cell_to_index(cell(1), cell(2), cell(3))
1421 11806382 : IF (icell < 1 .OR. icell > nimg) CYCLE
1422 :
1423 11805320 : arg = REAL(cell(1), dp)*xkp(1) + REAL(cell(2), dp)*xkp(2) + REAL(cell(3), dp)*xkp(3)
1424 11805320 : coskl = wkp*COS(twopi*arg)
1425 11805320 : sinkl = wkp*fc*SIN(twopi*arg)
1426 :
1427 : CALL dbcsr_get_block_p(matrix=denmat(ispin, icell)%matrix, row=irow, col=icol, &
1428 11805320 : block=dblock, found=found)
1429 11805320 : IF (.NOT. found) CYCLE
1430 :
1431 11832546 : IF (real_only) THEN
1432 176136 : CALL dbcsr_get_block_p(matrix=rpmat, row=irow, col=icol, block=rblock, found=found)
1433 176136 : IF (.NOT. found) CYCLE
1434 92594280 : dblock = dblock + coskl*rblock
1435 : ELSE
1436 11629184 : CALL dbcsr_get_block_p(matrix=rpmat, row=irow, col=icol, block=rblock, found=found)
1437 11629184 : IF (.NOT. found) CYCLE
1438 11629184 : CALL dbcsr_get_block_p(matrix=cpmat, row=irow, col=icol, block=cblock, found=found)
1439 11629184 : IF (.NOT. found) CYCLE
1440 2340218674 : dblock = dblock + coskl*rblock
1441 2340218674 : dblock = dblock + sinkl*cblock
1442 : END IF
1443 : END DO
1444 27226 : CALL neighbor_list_iterator_release(nl_iterator)
1445 :
1446 27226 : CALL timestop(handle)
1447 :
1448 27226 : END SUBROUTINE transform_dmat
1449 :
1450 : ! **************************************************************************************************
1451 : !> \brief Symmetrization of density matrix - transform to new k-point
1452 : !> \param smat density matrix at new kpoint
1453 : !> \param pmat reference density matrix
1454 : !> \param kmat Kind type rotation matrix
1455 : !> \param rot Rotation matrix
1456 : !> \param f0 Permutation of atoms under transformation
1457 : !> \param atype Atom to kind pointer
1458 : !> \param symmetric Symmetric matrix
1459 : !> \param antisymmetric Anti-Symmetric matrix
1460 : ! **************************************************************************************************
1461 0 : SUBROUTINE symtrans(smat, pmat, kmat, rot, f0, atype, symmetric, antisymmetric)
1462 : TYPE(dbcsr_type), POINTER :: smat, pmat
1463 : TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kmat
1464 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: rot
1465 : INTEGER, DIMENSION(:), INTENT(IN) :: f0, atype
1466 : LOGICAL, INTENT(IN), OPTIONAL :: symmetric, antisymmetric
1467 :
1468 : CHARACTER(LEN=*), PARAMETER :: routineN = 'symtrans'
1469 :
1470 : INTEGER :: handle, iatom, icol, ikind, ip, irow, &
1471 : jcol, jkind, jp, jrow, natom, numnodes
1472 : LOGICAL :: asym, dorot, found, perm, sym, trans
1473 : REAL(KIND=dp) :: dr, fsign
1474 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: kroti, krotj, pblock, sblock
1475 : TYPE(dbcsr_distribution_type) :: dist
1476 : TYPE(dbcsr_iterator_type) :: iter
1477 :
1478 0 : CALL timeset(routineN, handle)
1479 :
1480 : ! check symmetry options
1481 0 : sym = .FALSE.
1482 0 : IF (PRESENT(symmetric)) sym = symmetric
1483 0 : asym = .FALSE.
1484 0 : IF (PRESENT(antisymmetric)) asym = antisymmetric
1485 :
1486 0 : CPASSERT(.NOT. (sym .AND. asym))
1487 0 : CPASSERT((sym .OR. asym))
1488 :
1489 : ! do we have permutation of atoms
1490 0 : natom = SIZE(f0)
1491 0 : perm = .FALSE.
1492 0 : DO iatom = 1, natom
1493 0 : IF (f0(iatom) == iatom) CYCLE
1494 : perm = .TRUE.
1495 0 : EXIT
1496 : END DO
1497 :
1498 : ! do we have a real rotation
1499 0 : dorot = .FALSE.
1500 0 : IF (ABS(SUM(ABS(rot)) - 3.0_dp) > 1.e-12_dp) dorot = .TRUE.
1501 0 : dr = ABS(rot(1, 1) - 1.0_dp) + ABS(rot(2, 2) - 1.0_dp) + ABS(rot(3, 3) - 1.0_dp)
1502 0 : IF (ABS(dr) > 1.e-12_dp) dorot = .TRUE.
1503 :
1504 0 : fsign = 1.0_dp
1505 0 : IF (asym) fsign = -1.0_dp
1506 :
1507 0 : IF (dorot .OR. perm) THEN
1508 : CALL cp_abort(__LOCATION__, "k-points need FULL_GRID currently. "// &
1509 0 : "Reduced grids not yet working correctly")
1510 0 : CALL dbcsr_set(smat, 0.0_dp)
1511 0 : IF (perm) THEN
1512 0 : CALL dbcsr_get_info(pmat, distribution=dist)
1513 0 : CALL dbcsr_distribution_get(dist, numnodes=numnodes)
1514 0 : IF (numnodes == 1) THEN
1515 : ! the matrices are local to this process
1516 0 : CALL dbcsr_iterator_start(iter, pmat)
1517 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1518 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, pblock)
1519 0 : ip = f0(irow)
1520 0 : jp = f0(icol)
1521 0 : IF (ip <= jp) THEN
1522 0 : jrow = ip
1523 0 : jcol = jp
1524 0 : trans = .FALSE.
1525 : ELSE
1526 0 : jrow = jp
1527 0 : jcol = ip
1528 0 : trans = .TRUE.
1529 : END IF
1530 0 : CALL dbcsr_get_block_p(matrix=smat, row=jrow, col=jcol, BLOCK=sblock, found=found)
1531 0 : CPASSERT(found)
1532 0 : ikind = atype(irow)
1533 0 : jkind = atype(icol)
1534 0 : kroti => kmat(ikind)%rmat
1535 0 : krotj => kmat(jkind)%rmat
1536 : ! rotation
1537 0 : IF (trans) THEN
1538 0 : sblock = fsign*MATMUL(TRANSPOSE(krotj), MATMUL(TRANSPOSE(pblock), kroti))
1539 : ELSE
1540 0 : sblock = fsign*MATMUL(TRANSPOSE(kroti), MATMUL(pblock, krotj))
1541 : END IF
1542 : END DO
1543 0 : CALL dbcsr_iterator_stop(iter)
1544 : !
1545 : ELSE
1546 : ! distributed matrices, most general code needed
1547 : CALL cp_abort(__LOCATION__, "k-points need FULL_GRID currently. "// &
1548 0 : "Reduced grids not yet working correctly")
1549 : END IF
1550 : ELSE
1551 : ! no atom permutations, this is always local
1552 0 : CALL dbcsr_copy(smat, pmat)
1553 0 : CALL dbcsr_iterator_start(iter, smat)
1554 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1555 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
1556 0 : ip = f0(irow)
1557 0 : jp = f0(icol)
1558 0 : IF (ip <= jp) THEN
1559 0 : jrow = ip
1560 0 : jcol = jp
1561 0 : trans = .FALSE.
1562 : ELSE
1563 0 : jrow = jp
1564 0 : jcol = ip
1565 0 : trans = .TRUE.
1566 : END IF
1567 0 : ikind = atype(irow)
1568 0 : jkind = atype(icol)
1569 0 : kroti => kmat(ikind)%rmat
1570 0 : krotj => kmat(jkind)%rmat
1571 : ! rotation
1572 0 : IF (trans) THEN
1573 0 : sblock = fsign*MATMUL(TRANSPOSE(krotj), MATMUL(TRANSPOSE(sblock), kroti))
1574 : ELSE
1575 0 : sblock = fsign*MATMUL(TRANSPOSE(kroti), MATMUL(sblock, krotj))
1576 : END IF
1577 : END DO
1578 0 : CALL dbcsr_iterator_stop(iter)
1579 : !
1580 : END IF
1581 : ELSE
1582 : ! this is the identity operation, just copy the matrix
1583 0 : CALL dbcsr_copy(smat, pmat)
1584 : END IF
1585 :
1586 0 : CALL timestop(handle)
1587 :
1588 0 : END SUBROUTINE symtrans
1589 :
1590 : ! **************************************************************************************************
1591 : !> \brief ...
1592 : !> \param mat ...
1593 : ! **************************************************************************************************
1594 0 : SUBROUTINE matprint(mat)
1595 : TYPE(dbcsr_type), POINTER :: mat
1596 :
1597 : INTEGER :: i, icol, iounit, irow
1598 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: mblock
1599 : TYPE(dbcsr_iterator_type) :: iter
1600 :
1601 0 : iounit = cp_logger_get_default_io_unit()
1602 0 : CALL dbcsr_iterator_start(iter, mat)
1603 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1604 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, mblock)
1605 : !
1606 0 : IF (iounit > 0) THEN
1607 0 : WRITE (iounit, '(A,2I4)') 'BLOCK ', irow, icol
1608 0 : DO i = 1, SIZE(mblock, 1)
1609 0 : WRITE (iounit, '(8F12.6)') mblock(i, :)
1610 : END DO
1611 : END IF
1612 : !
1613 : END DO
1614 0 : CALL dbcsr_iterator_stop(iter)
1615 :
1616 0 : END SUBROUTINE matprint
1617 : ! **************************************************************************************************
1618 :
1619 0 : END MODULE kpoint_methods
|
