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