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