Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 7334 : 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 7334 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atype
115 : LOGICAL :: spez
116 : REAL(KIND=dp) :: wsum
117 7334 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coord, scoord
118 6475922 : TYPE(csym_type) :: crys_sym
119 : TYPE(kpoint_sym_type), POINTER :: kpsym
120 :
121 7334 : CALL timeset(routineN, handle)
122 :
123 7334 : CPASSERT(ASSOCIATED(kpoint))
124 :
125 7350 : 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 144 : IF (.NOT. kpoint%symmetry) THEN
139 : ! we set up a random molecule to avoid any possible symmetry
140 96 : natom = 10
141 96 : ALLOCATE (coord(3, natom), scoord(3, natom), atype(natom))
142 1056 : DO i = 1, natom
143 960 : atype(i) = i
144 960 : coord(1, i) = SIN(i*0.12345_dp)
145 960 : coord(2, i) = COS(i*0.23456_dp)
146 960 : coord(3, i) = SIN(i*0.34567_dp)
147 1056 : 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 144 : IF (kpoint%verbose) THEN
158 8 : iounit = cp_logger_get_default_io_unit()
159 : ELSE
160 136 : iounit = -1
161 : END IF
162 : ! kind type list
163 432 : ALLOCATE (kpoint%atype(natom))
164 1416 : kpoint%atype = atype
165 :
166 144 : 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 144 : full_grid=kpoint%full_grid)
169 144 : kpoint%nkp = crys_sym%nkpoint
170 720 : ALLOCATE (kpoint%xkp(3, kpoint%nkp), kpoint%wkp(kpoint%nkp))
171 1768 : wsum = SUM(crys_sym%wkpoint)
172 1768 : DO ik = 1, kpoint%nkp
173 6496 : kpoint%xkp(1:3, ik) = crys_sym%xkpoint(1:3, ik)
174 1768 : kpoint%wkp(ik) = crys_sym%wkpoint(ik)/wsum
175 : END DO
176 :
177 : ! print output
178 144 : IF (kpoint%symmetry) CALL print_crys_symmetry(crys_sym)
179 144 : IF (kpoint%symmetry) CALL print_kp_symmetry(crys_sym)
180 :
181 : ! transfer symmetry information
182 2056 : ALLOCATE (kpoint%kp_sym(kpoint%nkp))
183 1768 : DO ik = 1, kpoint%nkp
184 1624 : NULLIFY (kpoint%kp_sym(ik)%kpoint_sym)
185 1624 : CALL kpoint_sym_create(kpoint%kp_sym(ik)%kpoint_sym)
186 1624 : kpsym => kpoint%kp_sym(ik)%kpoint_sym
187 1768 : 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 144 : 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 144 : CALL release_csym_type(crys_sym)
239 144 : 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 7334 : CPASSERT(.FALSE.)
250 : END SELECT
251 :
252 : ! check for consistency of options
253 7350 : 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 7334 : CPASSERT(kpoint%nkp >= 1)
266 : END SELECT
267 7334 : 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 7334 : CALL timestop(handle)
284 :
285 7334 : 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 226 : 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 226 : TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_aux_env, kp_env
308 : TYPE(kpoint_env_type), POINTER :: kp
309 226 : TYPE(mp_cart_type) :: comm_cart
310 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp, para_env_kp
311 :
312 226 : CALL timeset(routineN, handle)
313 :
314 226 : IF (PRESENT(with_aux_fit)) THEN
315 162 : aux_fit = with_aux_fit
316 : ELSE
317 : aux_fit = .FALSE.
318 : END IF
319 :
320 226 : kpoint%para_env => para_env
321 226 : CALL kpoint%para_env%retain()
322 226 : kpoint%blacs_env_all => blacs_env
323 226 : CALL kpoint%blacs_env_all%retain()
324 :
325 226 : CPASSERT(.NOT. ASSOCIATED(kpoint%kp_env))
326 226 : IF (aux_fit) THEN
327 28 : CPASSERT(.NOT. ASSOCIATED(kpoint%kp_aux_env))
328 : END IF
329 :
330 226 : NULLIFY (kp_env, kp_aux_env)
331 226 : nkp = kpoint%nkp
332 226 : npe = para_env%num_pe
333 226 : 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 226 : 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 336 : DO igr = npe, 1, -1
381 224 : IF (MOD(npe, igr) .NE. 0) CYCLE
382 224 : nkp_grp = npe/igr
383 224 : IF (MOD(nkp, nkp_grp) .NE. 0) CYCLE
384 336 : ngr = igr
385 : END DO
386 114 : ELSE IF (kpoint%parallel_group_size == 0) THEN
387 : ! no parallelization over kpoints
388 60 : ngr = npe
389 54 : ELSE IF (kpoint%parallel_group_size > 0) THEN
390 54 : ngr = MIN(kpoint%parallel_group_size, npe)
391 : ELSE
392 0 : CPASSERT(.FALSE.)
393 : END IF
394 226 : nkp_grp = npe/ngr
395 : ! processor dimensions
396 226 : dims(1) = ngr
397 226 : dims(2) = nkp_grp
398 226 : CPASSERT(MOD(nkp, nkp_grp) == 0)
399 226 : nkp_loc = nkp/nkp_grp
400 :
401 226 : 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 226 : CALL comm_cart%create(comm_old=para_env, ndims=2, dims=dims)
407 678 : pos = comm_cart%mepos_cart
408 226 : ALLOCATE (para_env_kp)
409 226 : CALL para_env_kp%from_split(comm_cart, pos(2))
410 226 : ALLOCATE (para_env_inter_kp)
411 226 : CALL para_env_inter_kp%from_split(comm_cart, pos(1))
412 226 : CALL comm_cart%free()
413 :
414 226 : niogrp = 0
415 226 : IF (para_env%is_source()) niogrp = 1
416 226 : CALL para_env_kp%sum(niogrp)
417 226 : kpoint%iogrp = (niogrp == 1)
418 :
419 : ! parallel groups
420 226 : kpoint%para_env_kp => para_env_kp
421 226 : kpoint%para_env_inter_kp => para_env_inter_kp
422 :
423 : ! distribution of kpoints
424 678 : ALLOCATE (kpoint%kp_dist(2, nkp_grp))
425 532 : DO igr = 1, nkp_grp
426 1144 : kpoint%kp_dist(1:2, igr) = get_limit(nkp, nkp_grp, igr - 1)
427 : END DO
428 : ! local kpoints
429 678 : kpoint%kp_range(1:2) = kpoint%kp_dist(1:2, para_env_inter_kp%mepos + 1)
430 :
431 2214 : ALLOCATE (kp_env(nkp_loc))
432 1762 : DO ik = 1, nkp_loc
433 1536 : NULLIFY (kp_env(ik)%kpoint_env)
434 1536 : ikk = kpoint%kp_range(1) + ik - 1
435 1536 : CALL kpoint_env_create(kp_env(ik)%kpoint_env)
436 1536 : kp => kp_env(ik)%kpoint_env
437 1536 : kp%nkpoint = ikk
438 1536 : kp%wkp = kpoint%wkp(ikk)
439 6144 : kp%xkp(1:3) = kpoint%xkp(1:3, ikk)
440 1762 : kp%is_local = (ngr == 1)
441 : END DO
442 226 : kpoint%kp_env => kp_env
443 :
444 226 : 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 226 : unit_nr = cp_logger_get_default_io_unit()
460 :
461 226 : IF (unit_nr > 0 .AND. kpoint%verbose) THEN
462 4 : WRITE (unit_nr, *)
463 4 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Number of kpoint groups ", nkp_grp
464 4 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Size of each kpoint group", ngr
465 4 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Number of kpoints per group", nkp_loc
466 : END IF
467 226 : kpoint%nkp_groups = nkp_grp
468 :
469 : END IF
470 :
471 226 : CALL timestop(handle)
472 :
473 452 : 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 254 : 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 254 : 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 254 : CALL timeset(routineN, handle)
504 :
505 254 : IF (PRESENT(for_aux_fit)) THEN
506 28 : aux_fit = for_aux_fit
507 : ELSE
508 : aux_fit = .FALSE.
509 : END IF
510 :
511 254 : CPASSERT(ASSOCIATED(kpoint))
512 :
513 : IF (.TRUE. .OR. ASSOCIATED(mos(1)%mo_coeff)) THEN
514 254 : IF (aux_fit) THEN
515 28 : CPASSERT(ASSOCIATED(kpoint%kp_aux_env))
516 : END IF
517 :
518 254 : IF (PRESENT(added_mos)) THEN
519 64 : nadd = added_mos
520 : ELSE
521 : nadd = 0
522 : END IF
523 :
524 254 : IF (kpoint%use_real_wfn) THEN
525 : nc = 1
526 : ELSE
527 242 : nc = 2
528 : END IF
529 254 : nspin = SIZE(mos, 1)
530 254 : nkp_loc = kpoint%kp_range(2) - kpoint%kp_range(1) + 1
531 254 : IF (nkp_loc > 0) THEN
532 254 : IF (aux_fit) THEN
533 28 : CPASSERT(SIZE(kpoint%kp_aux_env) == nkp_loc)
534 : ELSE
535 226 : CPASSERT(SIZE(kpoint%kp_env) == nkp_loc)
536 : END IF
537 : ! allocate the mo sets, correct number of kpoints (local), real/complex, spin
538 1950 : DO ik = 1, nkp_loc
539 1696 : IF (aux_fit) THEN
540 160 : kp => kpoint%kp_aux_env(ik)%kpoint_env
541 : ELSE
542 1536 : kp => kpoint%kp_env(ik)%kpoint_env
543 : END IF
544 12674 : ALLOCATE (kp%mos(nc, nspin))
545 3918 : DO is = 1, nspin
546 : CALL get_mo_set(mos(is), nao=nao, nmo=nmo, nelectron=nelectron, &
547 1968 : n_el_f=n_el_f, maxocc=maxocc, flexible_electron_count=flexible_electron_count)
548 1968 : nmo = MIN(nao, nmo + nadd)
549 7586 : DO ic = 1, nc
550 : CALL allocate_mo_set(kp%mos(ic, is), nao, nmo, nelectron, n_el_f, maxocc, &
551 5890 : 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 254 : NULLIFY (blacs_env)
561 254 : IF (ASSOCIATED(kpoint%blacs_env)) THEN
562 28 : blacs_env => kpoint%blacs_env
563 : ELSE
564 226 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=kpoint%para_env_kp)
565 226 : kpoint%blacs_env => blacs_env
566 : END IF
567 :
568 : ! set possible new number of MOs
569 548 : DO is = 1, nspin
570 294 : CALL get_mo_set(mos(is), nmo=nmorig(is))
571 294 : nmo = MIN(nao, nmorig(is) + nadd)
572 548 : 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 254 : NULLIFY (mpools)
577 254 : CALL mpools_create(mpools=mpools)
578 : CALL mpools_rebuild_fm_pools(mpools=mpools, mos=mos, &
579 254 : blacs_env=blacs_env, para_env=kpoint%para_env_kp)
580 :
581 254 : IF (aux_fit) THEN
582 28 : kpoint%mpools_aux_fit => mpools
583 : ELSE
584 226 : kpoint%mpools => mpools
585 : END IF
586 :
587 : ! reset old number of MOs
588 548 : DO is = 1, nspin
589 548 : CALL set_mo_set(mos(is), nmo=nmorig(is))
590 : END DO
591 :
592 : ! allocate density matrices
593 254 : CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
594 254 : ALLOCATE (fmlocal)
595 254 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
596 254 : CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
597 1950 : DO ik = 1, nkp_loc
598 1696 : IF (aux_fit) THEN
599 160 : kp => kpoint%kp_aux_env(ik)%kpoint_env
600 : ELSE
601 1536 : kp => kpoint%kp_env(ik)%kpoint_env
602 : END IF
603 : ! density matrix
604 1696 : CALL cp_fm_release(kp%pmat)
605 12674 : ALLOCATE (kp%pmat(nc, nspin))
606 3664 : DO is = 1, nspin
607 7586 : DO ic = 1, nc
608 5890 : CALL cp_fm_create(kp%pmat(ic, is), matrix_struct)
609 : END DO
610 : END DO
611 : ! energy weighted density matrix
612 1696 : CALL cp_fm_release(kp%wmat)
613 10978 : ALLOCATE (kp%wmat(nc, nspin))
614 3918 : DO is = 1, nspin
615 7586 : DO ic = 1, nc
616 5890 : CALL cp_fm_create(kp%wmat(ic, is), matrix_struct)
617 : END DO
618 : END DO
619 : END DO
620 254 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
621 254 : DEALLOCATE (fmlocal)
622 :
623 : END IF
624 :
625 : END IF
626 :
627 254 : CALL timestop(handle)
628 :
629 254 : END SUBROUTINE kpoint_initialize_mos
630 :
631 : ! **************************************************************************************************
632 : !> \brief ...
633 : !> \param kpoint ...
634 : ! **************************************************************************************************
635 64 : 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 64 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
642 : TYPE(cp_fm_type), POINTER :: mo_coeff
643 64 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: moskp
644 :
645 64 : CALL timeset(routineN, handle)
646 :
647 444 : DO ik = 1, SIZE(kpoint%kp_env)
648 380 : CALL mpools_get(kpoint%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
649 380 : moskp => kpoint%kp_env(ik)%kpoint_env%mos
650 380 : ikk = kpoint%kp_range(1) + ik - 1
651 380 : CPASSERT(ASSOCIATED(moskp))
652 892 : DO ispin = 1, SIZE(moskp, 2)
653 1724 : DO ic = 1, SIZE(moskp, 1)
654 896 : CALL get_mo_set(moskp(ic, ispin), mo_coeff=mo_coeff)
655 1344 : IF (.NOT. ASSOCIATED(mo_coeff)) THEN
656 : CALL init_mo_set(moskp(ic, ispin), &
657 896 : fm_pool=ao_mo_fm_pools(ispin)%pool, name="kpoints")
658 : END IF
659 : END DO
660 : END DO
661 : END DO
662 :
663 64 : CALL timestop(handle)
664 :
665 64 : 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 992 : 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 992 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell, list
689 992 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index, cti
690 : LOGICAL :: new
691 : TYPE(neighbor_list_iterator_p_type), &
692 992 : DIMENSION(:), POINTER :: nl_iterator
693 :
694 992 : NULLIFY (cell_to_index, index_to_cell)
695 :
696 992 : CALL timeset(routineN, handle)
697 :
698 992 : CPASSERT(ASSOCIATED(kpoint))
699 :
700 992 : ALLOCATE (list(3, 125))
701 496992 : list = 0
702 992 : icount = 1
703 :
704 992 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
705 599262 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
706 598270 : CALL get_iterator_info(nl_iterator, cell=cell)
707 :
708 598270 : new = .TRUE.
709 35771204 : DO ic = 1, icount
710 35690474 : IF (cell(1) == list(1, ic) .AND. cell(2) == list(2, ic) .AND. &
711 80730 : cell(3) == list(3, ic)) THEN
712 : new = .FALSE.
713 : EXIT
714 : END IF
715 : END DO
716 599262 : IF (new) THEN
717 80730 : icount = icount + 1
718 80730 : IF (icount > SIZE(list, 2)) THEN
719 273 : CALL reallocate(list, 1, 3, 1, 2*SIZE(list, 2))
720 : END IF
721 322920 : list(1:3, icount) = cell(1:3)
722 : END IF
723 :
724 : END DO
725 992 : CALL neighbor_list_iterator_release(nl_iterator)
726 :
727 82714 : itm(1) = MAXVAL(ABS(list(1, 1:icount)))
728 82714 : itm(2) = MAXVAL(ABS(list(2, 1:icount)))
729 82714 : itm(3) = MAXVAL(ABS(list(3, 1:icount)))
730 992 : CALL para_env%max(itm)
731 3968 : it = MAXVAL(itm(1:3))
732 992 : IF (ASSOCIATED(kpoint%cell_to_index)) THEN
733 988 : DEALLOCATE (kpoint%cell_to_index)
734 : END IF
735 4960 : ALLOCATE (kpoint%cell_to_index(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
736 992 : cell_to_index => kpoint%cell_to_index
737 992 : cti => cell_to_index
738 199480 : cti(:, :, :) = 0
739 82714 : DO ic = 1, icount
740 81722 : i1 = list(1, ic)
741 81722 : i2 = list(2, ic)
742 81722 : i3 = list(3, ic)
743 82714 : cti(i1, i2, i3) = ic
744 : END DO
745 397968 : CALL para_env%sum(cti)
746 992 : ncount = 0
747 5768 : DO i1 = -itm(1), itm(1)
748 34852 : DO i2 = -itm(2), itm(2)
749 200040 : DO i3 = -itm(3), itm(3)
750 195264 : IF (cti(i1, i2, i3) == 0) THEN
751 75688 : cti(i1, i2, i3) = 1000000
752 : ELSE
753 90492 : ncount = ncount + 1
754 90492 : cti(i1, i2, i3) = (ABS(i1) + ABS(i2) + ABS(i3))*1000 + ABS(i3)*100 + ABS(i2)*10 + ABS(i1)
755 90492 : 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 992 : IF (ASSOCIATED(kpoint%index_to_cell)) THEN
762 992 : DEALLOCATE (kpoint%index_to_cell)
763 : END IF
764 2976 : ALLOCATE (kpoint%index_to_cell(3, ncount))
765 992 : index_to_cell => kpoint%index_to_cell
766 91484 : DO ic = 1, ncount
767 90492 : cell = MINLOC(cti)
768 90492 : i1 = cell(1) - 1 - itm(1)
769 90492 : i2 = cell(2) - 1 - itm(2)
770 90492 : i3 = cell(3) - 1 - itm(3)
771 90492 : cti(i1, i2, i3) = 1000000
772 90492 : index_to_cell(1, ic) = i1
773 90492 : index_to_cell(2, ic) = i2
774 91484 : index_to_cell(3, ic) = i3
775 : END DO
776 199480 : cti(:, :, :) = 0
777 91484 : DO ic = 1, ncount
778 90492 : i1 = index_to_cell(1, ic)
779 90492 : i2 = index_to_cell(2, ic)
780 90492 : i3 = index_to_cell(3, ic)
781 91484 : cti(i1, i2, i3) = ic
782 : END DO
783 :
784 : ! keep pointer to this neighborlist
785 992 : kpoint%sab_nl => sab_nl
786 :
787 : ! set number of images
788 992 : dft_control%nimages = SIZE(index_to_cell, 2)
789 :
790 992 : DEALLOCATE (list)
791 :
792 992 : CALL timestop(handle)
793 992 : 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 100752 : 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 50376 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, rblock, rsblock
830 : TYPE(neighbor_list_iterator_p_type), &
831 50376 : DIMENSION(:), POINTER :: nl_iterator
832 :
833 50376 : CALL timeset(routineN, handle)
834 :
835 50376 : my_complex = .FALSE.
836 50376 : IF (PRESENT(is_complex)) my_complex = is_complex
837 :
838 50376 : fsign = 1.0_dp
839 50376 : IF (PRESENT(rs_sign)) fsign = rs_sign
840 :
841 50376 : wfn_real_only = .TRUE.
842 50376 : IF (PRESENT(cmatrix)) wfn_real_only = .FALSE.
843 :
844 50376 : nimg = SIZE(rsmat, 2)
845 :
846 50376 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
847 :
848 50376 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
849 21459347 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
850 21408971 : 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 21408971 : fsym = 1.0_dp
855 21408971 : irow = iatom
856 21408971 : icol = jatom
857 21408971 : IF (do_symmetric .AND. (iatom > jatom)) THEN
858 9010083 : irow = jatom
859 9010083 : icol = iatom
860 9010083 : fsym = -1.0_dp
861 : END IF
862 :
863 21408971 : ic = cell_to_index(cell(1), cell(2), cell(3))
864 21408971 : IF (ic < 1 .OR. ic > nimg) CYCLE
865 :
866 21408227 : arg = REAL(cell(1), dp)*xkp(1) + REAL(cell(2), dp)*xkp(2) + REAL(cell(3), dp)*xkp(3)
867 21408227 : IF (my_complex) THEN
868 0 : coskl = fsign*fsym*COS(twopi*arg)
869 0 : sinkl = fsign*SIN(twopi*arg)
870 : ELSE
871 21408227 : coskl = fsign*COS(twopi*arg)
872 21408227 : 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 21408227 : block=rsblock, found=found)
877 21408227 : IF (.NOT. found) CYCLE
878 :
879 21458603 : 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 21104363 : block=rblock, found=found)
887 21104363 : IF (.NOT. found) CYCLE
888 : CALL dbcsr_get_block_p(matrix=cmatrix, row=irow, col=icol, &
889 21104363 : block=cblock, found=found)
890 21104363 : IF (.NOT. found) CYCLE
891 4150841687 : rblock = rblock + coskl*rsblock
892 4150841687 : cblock = cblock + sinkl*rsblock
893 : END IF
894 :
895 : END DO
896 50376 : CALL neighbor_list_iterator_release(nl_iterator)
897 :
898 50376 : CALL timestop(handle)
899 :
900 50376 : 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 5430 : 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 5430 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: weig, wocc
919 5430 : 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 5430 : CALL timeset(routineN, handle)
925 :
926 : ! first collect all the eigenvalues
927 5430 : CALL get_kpoint_info(kpoint, nkp=nkp)
928 5430 : kp => kpoint%kp_env(1)%kpoint_env
929 5430 : nspin = SIZE(kp%mos, 2)
930 5430 : mo_set => kp%mos(1, 1)
931 5430 : CALL get_mo_set(mo_set, nmo=nmo, nelectron=nelectron)
932 5430 : ne_a = nelectron
933 5430 : IF (nspin == 2) THEN
934 530 : CALL get_mo_set(kp%mos(1, 2), nmo=nb, nelectron=ne_b)
935 530 : CPASSERT(nmo == nb)
936 : END IF
937 43440 : ALLOCATE (weig(nmo, nkp, nspin), wocc(nmo, nkp, nspin))
938 435306 : weig = 0.0_dp
939 435306 : wocc = 0.0_dp
940 5430 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
941 5430 : kplocal = kp_range(2) - kp_range(1) + 1
942 18888 : DO ikpgr = 1, kplocal
943 13458 : ik = kp_range(1) + ikpgr - 1
944 13458 : kp => kpoint%kp_env(ikpgr)%kpoint_env
945 34194 : DO ispin = 1, nspin
946 15306 : mo_set => kp%mos(1, ispin)
947 15306 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
948 316128 : weig(1:nmo, ik, ispin) = eigenvalues(1:nmo)
949 : END DO
950 : END DO
951 5430 : CALL get_kpoint_info(kpoint, para_env_inter_kp=para_env_inter_kp)
952 5430 : CALL para_env_inter_kp%sum(weig)
953 :
954 5430 : CALL get_kpoint_info(kpoint, wkp=wkp)
955 5430 : IF (smear%do_smear) THEN
956 : ! finite electronic temperature
957 1364 : SELECT CASE (smear%method)
958 : CASE (smear_fermi_dirac)
959 682 : IF (nspin == 1) THEN
960 672 : nel = REAL(nelectron, KIND=dp)
961 : CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
962 672 : smear%electronic_temperature, 2.0_dp)
963 10 : ELSE IF (smear%fixed_mag_mom > 0.0_dp) THEN
964 0 : CPABORT("kpoints: Smearing with fixed magnetic moments not (yet) supported")
965 0 : nel = REAL(ne_a, KIND=dp)
966 : CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
967 0 : smear%electronic_temperature, 1.0_dp)
968 0 : nel = REAL(ne_b, KIND=dp)
969 : CALL Fermikp(wocc(:, :, 2), mus(2), kTS, weig(:, :, 2), nel, wkp, &
970 0 : smear%electronic_temperature, 1.0_dp)
971 : ELSE
972 10 : nel = REAL(ne_a, KIND=dp) + REAL(ne_b, KIND=dp)
973 : CALL Fermikp2(wocc(:, :, :), mu, kTS, weig(:, :, :), nel, wkp, &
974 10 : smear%electronic_temperature)
975 10 : kTS = kTS/2._dp
976 30 : mus(1:2) = mu
977 : END IF
978 : CASE DEFAULT
979 682 : CPABORT("kpoints: Selected smearing not (yet) supported")
980 : END SELECT
981 : ELSE
982 : ! fixed occupations (2/1)
983 4748 : IF (nspin == 1) THEN
984 4228 : nel = REAL(nelectron, KIND=dp)
985 4228 : CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, 0.0_dp, 2.0_dp)
986 : ELSE
987 520 : nel = REAL(ne_a, KIND=dp)
988 520 : CALL Fermikp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, 0.0_dp, 1.0_dp)
989 520 : nel = REAL(ne_b, KIND=dp)
990 520 : CALL Fermikp(wocc(:, :, 2), mus(2), kTS, weig(:, :, 2), nel, wkp, 0.0_dp, 1.0_dp)
991 : END IF
992 : END IF
993 18888 : DO ikpgr = 1, kplocal
994 13458 : ik = kp_range(1) + ikpgr - 1
995 13458 : kp => kpoint%kp_env(ikpgr)%kpoint_env
996 34194 : DO ispin = 1, nspin
997 15306 : mo_set => kp%mos(1, ispin)
998 15306 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues, occupation_numbers=occupation)
999 302670 : eigenvalues(1:nmo) = weig(1:nmo, ik, ispin)
1000 302670 : occupation(1:nmo) = wocc(1:nmo, ik, ispin)
1001 15306 : mo_set%kTS = kTS
1002 28764 : mo_set%mu = mus(ispin)
1003 : END DO
1004 : END DO
1005 :
1006 5430 : DEALLOCATE (weig, wocc)
1007 :
1008 5430 : CALL timestop(handle)
1009 :
1010 10860 : END SUBROUTINE kpoint_set_mo_occupation
1011 :
1012 : ! **************************************************************************************************
1013 : !> \brief Calculate kpoint density matrices (rho(k), owned by kpoint groups)
1014 : !> \param kpoint kpoint environment
1015 : !> \param energy_weighted calculate energy weighted density matrix
1016 : !> \param for_aux_fit ...
1017 : ! **************************************************************************************************
1018 16974 : SUBROUTINE kpoint_density_matrices(kpoint, energy_weighted, for_aux_fit)
1019 :
1020 : TYPE(kpoint_type), POINTER :: kpoint
1021 : LOGICAL, OPTIONAL :: energy_weighted, for_aux_fit
1022 :
1023 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_density_matrices'
1024 :
1025 : INTEGER :: handle, ikpgr, ispin, kplocal, nao, nmo, &
1026 : nspin
1027 : INTEGER, DIMENSION(2) :: kp_range
1028 : LOGICAL :: aux_fit, wtype
1029 5658 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation
1030 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1031 : TYPE(cp_fm_type) :: fwork
1032 : TYPE(cp_fm_type), POINTER :: cpmat, pmat, rpmat
1033 : TYPE(kpoint_env_type), POINTER :: kp
1034 : TYPE(mo_set_type), POINTER :: mo_set
1035 :
1036 5658 : CALL timeset(routineN, handle)
1037 :
1038 5658 : IF (PRESENT(energy_weighted)) THEN
1039 150 : wtype = energy_weighted
1040 : ELSE
1041 : ! default is normal density matrix
1042 : wtype = .FALSE.
1043 : END IF
1044 :
1045 5658 : IF (PRESENT(for_aux_fit)) THEN
1046 62 : aux_fit = for_aux_fit
1047 : ELSE
1048 : aux_fit = .FALSE.
1049 : END IF
1050 :
1051 62 : IF (aux_fit) THEN
1052 62 : CPASSERT(ASSOCIATED(kpoint%kp_aux_env))
1053 : END IF
1054 :
1055 : ! work matrix
1056 5658 : IF (aux_fit) THEN
1057 62 : mo_set => kpoint%kp_aux_env(1)%kpoint_env%mos(1, 1)
1058 : ELSE
1059 5596 : mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)
1060 : END IF
1061 5658 : CALL get_mo_set(mo_set, nao=nao, nmo=nmo)
1062 5658 : CALL cp_fm_get_info(mo_set%mo_coeff, matrix_struct=matrix_struct)
1063 5658 : CALL cp_fm_create(fwork, matrix_struct)
1064 :
1065 5658 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
1066 5658 : kplocal = kp_range(2) - kp_range(1) + 1
1067 20266 : DO ikpgr = 1, kplocal
1068 14608 : IF (aux_fit) THEN
1069 352 : kp => kpoint%kp_aux_env(ikpgr)%kpoint_env
1070 : ELSE
1071 14256 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1072 : END IF
1073 14608 : nspin = SIZE(kp%mos, 2)
1074 36970 : DO ispin = 1, nspin
1075 16704 : mo_set => kp%mos(1, ispin)
1076 16704 : IF (wtype) THEN
1077 818 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
1078 : END IF
1079 31312 : IF (kpoint%use_real_wfn) THEN
1080 72 : IF (wtype) THEN
1081 10 : pmat => kp%wmat(1, ispin)
1082 : ELSE
1083 62 : pmat => kp%pmat(1, ispin)
1084 : END IF
1085 72 : CALL get_mo_set(mo_set, occupation_numbers=occupation)
1086 72 : CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
1087 72 : CALL cp_fm_column_scale(fwork, occupation)
1088 72 : IF (wtype) THEN
1089 10 : CALL cp_fm_column_scale(fwork, eigenvalues)
1090 : END IF
1091 72 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, pmat)
1092 : ELSE
1093 16632 : IF (wtype) THEN
1094 808 : rpmat => kp%wmat(1, ispin)
1095 808 : cpmat => kp%wmat(2, ispin)
1096 : ELSE
1097 15824 : rpmat => kp%pmat(1, ispin)
1098 15824 : cpmat => kp%pmat(2, ispin)
1099 : END IF
1100 16632 : CALL get_mo_set(mo_set, occupation_numbers=occupation)
1101 16632 : CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
1102 16632 : CALL cp_fm_column_scale(fwork, occupation)
1103 16632 : IF (wtype) THEN
1104 808 : CALL cp_fm_column_scale(fwork, eigenvalues)
1105 : END IF
1106 : ! Re(c)*Re(c)
1107 16632 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, rpmat)
1108 16632 : mo_set => kp%mos(2, ispin)
1109 : ! Im(c)*Re(c)
1110 16632 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, cpmat)
1111 : ! Re(c)*Im(c)
1112 16632 : CALL parallel_gemm("N", "T", nao, nao, nmo, -1.0_dp, fwork, mo_set%mo_coeff, 1.0_dp, cpmat)
1113 16632 : CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
1114 16632 : CALL cp_fm_column_scale(fwork, occupation)
1115 16632 : IF (wtype) THEN
1116 808 : CALL cp_fm_column_scale(fwork, eigenvalues)
1117 : END IF
1118 : ! Im(c)*Im(c)
1119 16632 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 1.0_dp, rpmat)
1120 : END IF
1121 : END DO
1122 : END DO
1123 :
1124 5658 : CALL cp_fm_release(fwork)
1125 :
1126 5658 : CALL timestop(handle)
1127 :
1128 5658 : END SUBROUTINE kpoint_density_matrices
1129 :
1130 : ! **************************************************************************************************
1131 : !> \brief generate real space density matrices in DBCSR format
1132 : !> \param kpoint Kpoint environment
1133 : !> \param denmat Real space (DBCSR) density matrices
1134 : !> \param wtype True = energy weighted density matrix
1135 : !> False = normal density matrix
1136 : !> \param tempmat DBCSR matrix to be used as template
1137 : !> \param sab_nl ...
1138 : !> \param fmwork FM work matrices (kpoint group)
1139 : !> \param for_aux_fit ...
1140 : !> \param pmat_ext ...
1141 : ! **************************************************************************************************
1142 5824 : SUBROUTINE kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwork, for_aux_fit, pmat_ext)
1143 :
1144 : TYPE(kpoint_type), POINTER :: kpoint
1145 : TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
1146 : LOGICAL, INTENT(IN) :: wtype
1147 : TYPE(dbcsr_type), POINTER :: tempmat
1148 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1149 : POINTER :: sab_nl
1150 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fmwork
1151 : LOGICAL, OPTIONAL :: for_aux_fit
1152 : TYPE(cp_fm_type), DIMENSION(:, :, :), INTENT(IN), &
1153 : OPTIONAL :: pmat_ext
1154 :
1155 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_density_transform'
1156 :
1157 : INTEGER :: handle, ic, ik, ikk, indx, ir, ira, is, &
1158 : ispin, jr, nc, nimg, nkp, nspin
1159 5824 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1160 : LOGICAL :: aux_fit, do_ext, do_symmetric, my_kpgrp, &
1161 : real_only
1162 : REAL(KIND=dp) :: wkpx
1163 5824 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
1164 5824 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1165 5824 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:) :: info
1166 : TYPE(cp_fm_type) :: fmdummy
1167 : TYPE(dbcsr_type), POINTER :: cpmat, rpmat, scpmat, srpmat
1168 5824 : TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kind_rot
1169 : TYPE(kpoint_env_type), POINTER :: kp
1170 : TYPE(kpoint_sym_type), POINTER :: kpsym
1171 : TYPE(mp_para_env_type), POINTER :: para_env
1172 :
1173 5824 : CALL timeset(routineN, handle)
1174 :
1175 5824 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1176 :
1177 5824 : IF (PRESENT(for_aux_fit)) THEN
1178 228 : aux_fit = for_aux_fit
1179 : ELSE
1180 : aux_fit = .FALSE.
1181 : END IF
1182 :
1183 5824 : do_ext = .FALSE.
1184 5824 : IF (PRESENT(pmat_ext)) do_ext = .TRUE.
1185 :
1186 5824 : IF (aux_fit) THEN
1187 138 : CPASSERT(ASSOCIATED(kpoint%kp_aux_env))
1188 : END IF
1189 :
1190 : ! work storage
1191 5824 : ALLOCATE (rpmat)
1192 : CALL dbcsr_create(rpmat, template=tempmat, &
1193 5872 : matrix_type=MERGE(dbcsr_type_symmetric, dbcsr_type_no_symmetry, do_symmetric))
1194 5824 : CALL cp_dbcsr_alloc_block_from_nbl(rpmat, sab_nl)
1195 5824 : CALL dbcsr_set(rpmat, 0.0_dp)
1196 5824 : ALLOCATE (cpmat)
1197 : CALL dbcsr_create(cpmat, template=tempmat, &
1198 5872 : matrix_type=MERGE(dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, do_symmetric))
1199 5824 : CALL cp_dbcsr_alloc_block_from_nbl(cpmat, sab_nl)
1200 5824 : CALL dbcsr_set(cpmat, 0.0_dp)
1201 5824 : IF (.NOT. kpoint%full_grid) THEN
1202 4772 : ALLOCATE (srpmat)
1203 4772 : CALL dbcsr_create(srpmat, template=rpmat)
1204 4772 : CALL cp_dbcsr_alloc_block_from_nbl(srpmat, sab_nl)
1205 4772 : CALL dbcsr_set(srpmat, 0.0_dp)
1206 4772 : ALLOCATE (scpmat)
1207 4772 : CALL dbcsr_create(scpmat, template=cpmat)
1208 4772 : CALL cp_dbcsr_alloc_block_from_nbl(scpmat, sab_nl)
1209 4772 : CALL dbcsr_set(scpmat, 0.0_dp)
1210 : END IF
1211 :
1212 : CALL get_kpoint_info(kpoint, nkp=nkp, xkp=xkp, wkp=wkp, &
1213 5824 : cell_to_index=cell_to_index)
1214 : ! initialize real space density matrices
1215 5824 : IF (aux_fit) THEN
1216 138 : kp => kpoint%kp_aux_env(1)%kpoint_env
1217 : ELSE
1218 5686 : kp => kpoint%kp_env(1)%kpoint_env
1219 : END IF
1220 5824 : nspin = SIZE(kp%mos, 2)
1221 5824 : nc = SIZE(kp%mos, 1)
1222 5824 : nimg = SIZE(denmat, 2)
1223 5824 : real_only = (nc == 1)
1224 :
1225 5824 : para_env => kpoint%blacs_env_all%para_env
1226 118200 : ALLOCATE (info(nspin*nkp*nc))
1227 :
1228 : ! Start all the communication
1229 5824 : indx = 0
1230 12258 : DO ispin = 1, nspin
1231 531026 : DO ic = 1, nimg
1232 531026 : CALL dbcsr_set(denmat(ispin, ic)%matrix, 0.0_dp)
1233 : END DO
1234 : !
1235 39362 : DO ik = 1, nkp
1236 27104 : my_kpgrp = (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2))
1237 : IF (my_kpgrp) THEN
1238 17920 : ikk = ik - kpoint%kp_range(1) + 1
1239 17920 : IF (aux_fit) THEN
1240 1008 : kp => kpoint%kp_aux_env(ikk)%kpoint_env
1241 : ELSE
1242 16912 : kp => kpoint%kp_env(ikk)%kpoint_env
1243 : END IF
1244 : ELSE
1245 : NULLIFY (kp)
1246 : END IF
1247 : ! collect this density matrix on all processors
1248 27104 : CPASSERT(SIZE(fmwork) >= nc)
1249 :
1250 33538 : IF (my_kpgrp) THEN
1251 53688 : DO ic = 1, nc
1252 35768 : indx = indx + 1
1253 53688 : IF (do_ext) THEN
1254 2016 : CALL cp_fm_start_copy_general(pmat_ext(ikk, ic, ispin), fmwork(ic), para_env, info(indx))
1255 : ELSE
1256 33752 : IF (wtype) THEN
1257 1626 : CALL cp_fm_start_copy_general(kp%wmat(ic, ispin), fmwork(ic), para_env, info(indx))
1258 : ELSE
1259 32126 : CALL cp_fm_start_copy_general(kp%pmat(ic, ispin), fmwork(ic), para_env, info(indx))
1260 : END IF
1261 : END IF
1262 : END DO
1263 : ELSE
1264 27552 : DO ic = 1, nc
1265 18368 : indx = indx + 1
1266 27552 : CALL cp_fm_start_copy_general(fmdummy, fmwork(ic), para_env, info(indx))
1267 : END DO
1268 : END IF
1269 : END DO
1270 : END DO
1271 :
1272 : ! Finish communication and transform the received matrices
1273 5824 : indx = 0
1274 12258 : DO ispin = 1, nspin
1275 39362 : DO ik = 1, nkp
1276 81240 : DO ic = 1, nc
1277 54136 : indx = indx + 1
1278 81240 : CALL cp_fm_finish_copy_general(fmwork(ic), info(indx))
1279 : END DO
1280 :
1281 : ! reduce to dbcsr storage
1282 27104 : IF (real_only) THEN
1283 72 : CALL copy_fm_to_dbcsr(fmwork(1), rpmat, keep_sparsity=.TRUE.)
1284 : ELSE
1285 27032 : CALL copy_fm_to_dbcsr(fmwork(1), rpmat, keep_sparsity=.TRUE.)
1286 27032 : CALL copy_fm_to_dbcsr(fmwork(2), cpmat, keep_sparsity=.TRUE.)
1287 : END IF
1288 :
1289 : ! symmetrization
1290 27104 : kpsym => kpoint%kp_sym(ik)%kpoint_sym
1291 27104 : CPASSERT(ASSOCIATED(kpsym))
1292 :
1293 33538 : IF (kpsym%apply_symmetry) THEN
1294 0 : wkpx = wkp(ik)/REAL(kpsym%nwght, KIND=dp)
1295 0 : DO is = 1, kpsym%nwght
1296 0 : ir = ABS(kpsym%rotp(is))
1297 0 : ira = 0
1298 0 : DO jr = 1, SIZE(kpoint%ibrot)
1299 0 : IF (ir == kpoint%ibrot(jr)) ira = jr
1300 : END DO
1301 0 : CPASSERT(ira > 0)
1302 0 : kind_rot => kpoint%kind_rotmat(ira, :)
1303 0 : IF (real_only) THEN
1304 : CALL symtrans(srpmat, rpmat, kind_rot, kpsym%rot(1:3, 1:3, is), &
1305 0 : kpsym%f0(:, is), kpoint%atype, symmetric=.TRUE.)
1306 : ELSE
1307 : CALL symtrans(srpmat, rpmat, kind_rot, kpsym%rot(1:3, 1:3, is), &
1308 0 : kpsym%f0(:, is), kpoint%atype, symmetric=.TRUE.)
1309 : CALL symtrans(scpmat, cpmat, kind_rot, kpsym%rot(1:3, 1:3, is), &
1310 0 : kpsym%f0(:, is), kpoint%atype, antisymmetric=.TRUE.)
1311 : END IF
1312 : CALL transform_dmat(denmat, srpmat, scpmat, ispin, real_only, sab_nl, &
1313 0 : cell_to_index, kpsym%xkp(1:3, is), wkpx)
1314 : END DO
1315 : ELSE
1316 : ! transformation
1317 : CALL transform_dmat(denmat, rpmat, cpmat, ispin, real_only, sab_nl, &
1318 27104 : cell_to_index, xkp(1:3, ik), wkp(ik))
1319 : END IF
1320 : END DO
1321 : END DO
1322 :
1323 : ! Clean up communication
1324 5824 : indx = 0
1325 12258 : DO ispin = 1, nspin
1326 39362 : DO ik = 1, nkp
1327 27104 : my_kpgrp = (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2))
1328 6434 : IF (my_kpgrp) THEN
1329 17920 : ikk = ik - kpoint%kp_range(1) + 1
1330 : IF (aux_fit) THEN
1331 17920 : kp => kpoint%kp_aux_env(ikk)%kpoint_env
1332 : ELSE
1333 17920 : kp => kpoint%kp_env(ikk)%kpoint_env
1334 : END IF
1335 :
1336 53688 : DO ic = 1, nc
1337 35768 : indx = indx + 1
1338 53688 : CALL cp_fm_cleanup_copy_general(info(indx))
1339 : END DO
1340 : ELSE
1341 : ! calls with dummy arguments, so not included
1342 : ! therefore just increment counter by trip count
1343 9184 : indx = indx + nc
1344 : END IF
1345 : END DO
1346 : END DO
1347 :
1348 : ! All done
1349 59960 : DEALLOCATE (info)
1350 :
1351 5824 : CALL dbcsr_deallocate_matrix(rpmat)
1352 5824 : CALL dbcsr_deallocate_matrix(cpmat)
1353 5824 : IF (.NOT. kpoint%full_grid) THEN
1354 4772 : CALL dbcsr_deallocate_matrix(srpmat)
1355 4772 : CALL dbcsr_deallocate_matrix(scpmat)
1356 : END IF
1357 :
1358 5824 : CALL timestop(handle)
1359 :
1360 5824 : END SUBROUTINE kpoint_density_transform
1361 :
1362 : ! **************************************************************************************************
1363 : !> \brief real space density matrices in DBCSR format
1364 : !> \param denmat Real space (DBCSR) density matrix
1365 : !> \param rpmat ...
1366 : !> \param cpmat ...
1367 : !> \param ispin ...
1368 : !> \param real_only ...
1369 : !> \param sab_nl ...
1370 : !> \param cell_to_index ...
1371 : !> \param xkp ...
1372 : !> \param wkp ...
1373 : ! **************************************************************************************************
1374 27104 : SUBROUTINE transform_dmat(denmat, rpmat, cpmat, ispin, real_only, sab_nl, cell_to_index, xkp, wkp)
1375 :
1376 : TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
1377 : TYPE(dbcsr_type), POINTER :: rpmat, cpmat
1378 : INTEGER, INTENT(IN) :: ispin
1379 : LOGICAL, INTENT(IN) :: real_only
1380 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1381 : POINTER :: sab_nl
1382 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1383 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
1384 : REAL(KIND=dp), INTENT(IN) :: wkp
1385 :
1386 : CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_dmat'
1387 :
1388 : INTEGER :: handle, iatom, icell, icol, irow, jatom, &
1389 : nimg
1390 : INTEGER, DIMENSION(3) :: cell
1391 : LOGICAL :: do_symmetric, found
1392 : REAL(KIND=dp) :: arg, coskl, fc, sinkl
1393 27104 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, dblock, rblock
1394 : TYPE(neighbor_list_iterator_p_type), &
1395 27104 : DIMENSION(:), POINTER :: nl_iterator
1396 :
1397 27104 : CALL timeset(routineN, handle)
1398 :
1399 27104 : nimg = SIZE(denmat, 2)
1400 :
1401 : ! transformation
1402 27104 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1403 27104 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
1404 11433220 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1405 11406116 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
1406 :
1407 : !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
1408 : !Therefore, we have: S(R) = sum_k Re(S(k))*cos(k*R) -i^2*Im(S(k))*sin(k*R)
1409 : ! = sum_k Re(S(k))*cos(k*R) + Im(S(k))*sin(k*R)
1410 : !fc = +- 1 is due to the usual non-symmetric real-space matrices stored as symmetric ones
1411 :
1412 11406116 : irow = iatom
1413 11406116 : icol = jatom
1414 11406116 : fc = 1.0_dp
1415 11406116 : IF (do_symmetric .AND. iatom > jatom) THEN
1416 4781933 : irow = jatom
1417 4781933 : icol = iatom
1418 4781933 : fc = -1.0_dp
1419 : END IF
1420 :
1421 11406116 : icell = cell_to_index(cell(1), cell(2), cell(3))
1422 11406116 : IF (icell < 1 .OR. icell > nimg) CYCLE
1423 :
1424 11405054 : arg = REAL(cell(1), dp)*xkp(1) + REAL(cell(2), dp)*xkp(2) + REAL(cell(3), dp)*xkp(3)
1425 11405054 : coskl = wkp*COS(twopi*arg)
1426 11405054 : sinkl = wkp*fc*SIN(twopi*arg)
1427 :
1428 : CALL dbcsr_get_block_p(matrix=denmat(ispin, icell)%matrix, row=irow, col=icol, &
1429 11405054 : block=dblock, found=found)
1430 11405054 : IF (.NOT. found) CYCLE
1431 :
1432 11432158 : IF (real_only) THEN
1433 176136 : CALL dbcsr_get_block_p(matrix=rpmat, row=irow, col=icol, block=rblock, found=found)
1434 176136 : IF (.NOT. found) CYCLE
1435 92594280 : dblock = dblock + coskl*rblock
1436 : ELSE
1437 11228918 : CALL dbcsr_get_block_p(matrix=rpmat, row=irow, col=icol, block=rblock, found=found)
1438 11228918 : IF (.NOT. found) CYCLE
1439 11228918 : CALL dbcsr_get_block_p(matrix=cpmat, row=irow, col=icol, block=cblock, found=found)
1440 11228918 : IF (.NOT. found) CYCLE
1441 2269532456 : dblock = dblock + coskl*rblock
1442 2269532456 : dblock = dblock + sinkl*cblock
1443 : END IF
1444 : END DO
1445 27104 : CALL neighbor_list_iterator_release(nl_iterator)
1446 :
1447 27104 : CALL timestop(handle)
1448 :
1449 27104 : END SUBROUTINE transform_dmat
1450 :
1451 : ! **************************************************************************************************
1452 : !> \brief Symmetrization of density matrix - transform to new k-point
1453 : !> \param smat density matrix at new kpoint
1454 : !> \param pmat reference density matrix
1455 : !> \param kmat Kind type rotation matrix
1456 : !> \param rot Rotation matrix
1457 : !> \param f0 Permutation of atoms under transformation
1458 : !> \param atype Atom to kind pointer
1459 : !> \param symmetric Symmetric matrix
1460 : !> \param antisymmetric Anti-Symmetric matrix
1461 : ! **************************************************************************************************
1462 0 : SUBROUTINE symtrans(smat, pmat, kmat, rot, f0, atype, symmetric, antisymmetric)
1463 : TYPE(dbcsr_type), POINTER :: smat, pmat
1464 : TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kmat
1465 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: rot
1466 : INTEGER, DIMENSION(:), INTENT(IN) :: f0, atype
1467 : LOGICAL, INTENT(IN), OPTIONAL :: symmetric, antisymmetric
1468 :
1469 : CHARACTER(LEN=*), PARAMETER :: routineN = 'symtrans'
1470 :
1471 : INTEGER :: handle, iatom, icol, ikind, ip, irow, &
1472 : jcol, jkind, jp, jrow, natom, numnodes
1473 : LOGICAL :: asym, dorot, found, perm, sym, trans
1474 : REAL(KIND=dp) :: dr, fsign
1475 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: kroti, krotj, pblock, sblock
1476 : TYPE(dbcsr_distribution_type) :: dist
1477 : TYPE(dbcsr_iterator_type) :: iter
1478 :
1479 0 : CALL timeset(routineN, handle)
1480 :
1481 : ! check symmetry options
1482 0 : sym = .FALSE.
1483 0 : IF (PRESENT(symmetric)) sym = symmetric
1484 0 : asym = .FALSE.
1485 0 : IF (PRESENT(antisymmetric)) asym = antisymmetric
1486 :
1487 0 : CPASSERT(.NOT. (sym .AND. asym))
1488 0 : CPASSERT((sym .OR. asym))
1489 :
1490 : ! do we have permutation of atoms
1491 0 : natom = SIZE(f0)
1492 0 : perm = .FALSE.
1493 0 : DO iatom = 1, natom
1494 0 : IF (f0(iatom) == iatom) CYCLE
1495 : perm = .TRUE.
1496 0 : EXIT
1497 : END DO
1498 :
1499 : ! do we have a real rotation
1500 0 : dorot = .FALSE.
1501 0 : IF (ABS(SUM(ABS(rot)) - 3.0_dp) > 1.e-12_dp) dorot = .TRUE.
1502 0 : dr = ABS(rot(1, 1) - 1.0_dp) + ABS(rot(2, 2) - 1.0_dp) + ABS(rot(3, 3) - 1.0_dp)
1503 0 : IF (ABS(dr) > 1.e-12_dp) dorot = .TRUE.
1504 :
1505 0 : fsign = 1.0_dp
1506 0 : IF (asym) fsign = -1.0_dp
1507 :
1508 0 : IF (dorot .OR. perm) THEN
1509 : CALL cp_abort(__LOCATION__, "k-points need FULL_GRID currently. "// &
1510 0 : "Reduced grids not yet working correctly")
1511 0 : CALL dbcsr_set(smat, 0.0_dp)
1512 0 : IF (perm) THEN
1513 0 : CALL dbcsr_get_info(pmat, distribution=dist)
1514 0 : CALL dbcsr_distribution_get(dist, numnodes=numnodes)
1515 0 : IF (numnodes == 1) THEN
1516 : ! the matrices are local to this process
1517 0 : CALL dbcsr_iterator_start(iter, pmat)
1518 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1519 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, pblock)
1520 0 : ip = f0(irow)
1521 0 : jp = f0(icol)
1522 0 : IF (ip <= jp) THEN
1523 0 : jrow = ip
1524 0 : jcol = jp
1525 0 : trans = .FALSE.
1526 : ELSE
1527 0 : jrow = jp
1528 0 : jcol = ip
1529 0 : trans = .TRUE.
1530 : END IF
1531 0 : CALL dbcsr_get_block_p(matrix=smat, row=jrow, col=jcol, BLOCK=sblock, found=found)
1532 0 : CPASSERT(found)
1533 0 : ikind = atype(irow)
1534 0 : jkind = atype(icol)
1535 0 : kroti => kmat(ikind)%rmat
1536 0 : krotj => kmat(jkind)%rmat
1537 : ! rotation
1538 0 : IF (trans) THEN
1539 0 : sblock = fsign*MATMUL(TRANSPOSE(krotj), MATMUL(TRANSPOSE(pblock), kroti))
1540 : ELSE
1541 0 : sblock = fsign*MATMUL(TRANSPOSE(kroti), MATMUL(pblock, krotj))
1542 : END IF
1543 : END DO
1544 0 : CALL dbcsr_iterator_stop(iter)
1545 : !
1546 : ELSE
1547 : ! distributed matrices, most general code needed
1548 : CALL cp_abort(__LOCATION__, "k-points need FULL_GRID currently. "// &
1549 0 : "Reduced grids not yet working correctly")
1550 : END IF
1551 : ELSE
1552 : ! no atom permutations, this is always local
1553 0 : CALL dbcsr_copy(smat, pmat)
1554 0 : CALL dbcsr_iterator_start(iter, smat)
1555 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1556 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
1557 0 : ip = f0(irow)
1558 0 : jp = f0(icol)
1559 0 : IF (ip <= jp) THEN
1560 0 : jrow = ip
1561 0 : jcol = jp
1562 0 : trans = .FALSE.
1563 : ELSE
1564 0 : jrow = jp
1565 0 : jcol = ip
1566 0 : trans = .TRUE.
1567 : END IF
1568 0 : ikind = atype(irow)
1569 0 : jkind = atype(icol)
1570 0 : kroti => kmat(ikind)%rmat
1571 0 : krotj => kmat(jkind)%rmat
1572 : ! rotation
1573 0 : IF (trans) THEN
1574 0 : sblock = fsign*MATMUL(TRANSPOSE(krotj), MATMUL(TRANSPOSE(sblock), kroti))
1575 : ELSE
1576 0 : sblock = fsign*MATMUL(TRANSPOSE(kroti), MATMUL(sblock, krotj))
1577 : END IF
1578 : END DO
1579 0 : CALL dbcsr_iterator_stop(iter)
1580 : !
1581 : END IF
1582 : ELSE
1583 : ! this is the identity operation, just copy the matrix
1584 0 : CALL dbcsr_copy(smat, pmat)
1585 : END IF
1586 :
1587 0 : CALL timestop(handle)
1588 :
1589 0 : END SUBROUTINE symtrans
1590 :
1591 : ! **************************************************************************************************
1592 : !> \brief ...
1593 : !> \param mat ...
1594 : ! **************************************************************************************************
1595 0 : SUBROUTINE matprint(mat)
1596 : TYPE(dbcsr_type), POINTER :: mat
1597 :
1598 : INTEGER :: i, icol, iounit, irow
1599 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: mblock
1600 : TYPE(dbcsr_iterator_type) :: iter
1601 :
1602 0 : iounit = cp_logger_get_default_io_unit()
1603 0 : CALL dbcsr_iterator_start(iter, mat)
1604 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1605 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, mblock)
1606 : !
1607 0 : IF (iounit > 0) THEN
1608 0 : WRITE (iounit, '(A,2I4)') 'BLOCK ', irow, icol
1609 0 : DO i = 1, SIZE(mblock, 1)
1610 0 : WRITE (iounit, '(8F12.6)') mblock(i, :)
1611 : END DO
1612 : END IF
1613 : !
1614 : END DO
1615 0 : CALL dbcsr_iterator_stop(iter)
1616 :
1617 0 : END SUBROUTINE matprint
1618 : ! **************************************************************************************************
1619 :
1620 0 : END MODULE kpoint_methods
|
