(git:0622d44)
Loading...
Searching...
No Matches
kpoint_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Routines needed for kpoint calculation
10!> \par History
11!> 2014.07 created [JGH]
12!> 2014.11 unified k-point and gamma-point code [Ole Schuett]
13!> \author JGH
14! **************************************************************************************************
15MODULE kpoint_methods
16 USE atomic_kind_types, ONLY: get_atomic_kind
17 USE cell_types, ONLY: cell_type,&
18 real_to_scaled
19 USE cp_blacs_env, ONLY: cp_blacs_env_create,&
20 cp_blacs_env_type
21 USE cp_cfm_types, ONLY: cp_cfm_create,&
22 cp_cfm_release,&
23 cp_cfm_to_fm,&
24 cp_cfm_type,&
25 cp_fm_to_cfm
26 USE cp_control_types, ONLY: hairy_probes_type
27 USE cp_dbcsr_api, ONLY: &
28 dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribute, &
29 dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_get_block_p, dbcsr_get_info, &
30 dbcsr_get_stored_coordinates, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
31 dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
32 dbcsr_replicate_all, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
33 dbcsr_type_no_symmetry, dbcsr_type_symmetric
34 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
35 USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr
36 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
37 USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
38 fm_pool_create_fm,&
39 fm_pool_give_back_fm
40 USE cp_fm_struct, ONLY: cp_fm_struct_type
41 USE cp_fm_types, ONLY: &
42 copy_info_type, cp_fm_cleanup_copy_general, cp_fm_create, cp_fm_finish_copy_general, &
43 cp_fm_get_diag, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_release, &
44 cp_fm_start_copy_general, cp_fm_to_fm, cp_fm_type
45 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
46 USE cryssym, ONLY: crys_sym_gen,&
47 csym_type,&
48 kpoint_gen,&
49 kpoint_gen_general,&
50 print_crys_symmetry,&
51 print_kp_symmetry,&
52 release_csym_type
53 USE hairy_probes, ONLY: probe_occupancy_kp
54 USE input_constants, ONLY: smear_fermi_dirac,&
55 smear_gaussian,&
56 smear_mp,&
57 smear_mv
58 USE input_cp2k_kpoints, ONLY: use_spglib_kpoint_backend,&
59 use_spglib_kpoint_symmetry
60 USE kinds, ONLY: dp
61 USE kpoint_types, ONLY: get_kpoint_info,&
62 kind_rotmat_type,&
63 kpoint_env_create,&
64 kpoint_env_p_type,&
65 kpoint_env_type,&
66 kpoint_sym_create,&
67 kpoint_sym_type,&
68 kpoint_type
69 USE mathconstants, ONLY: twopi
70 USE mathlib, ONLY: inv_3x3
71 USE memory_utilities, ONLY: reallocate
72 USE message_passing, ONLY: mp_cart_type,&
73 mp_para_env_type
74 USE parallel_gemm_api, ONLY: parallel_gemm
75 USE particle_types, ONLY: particle_type
76 USE qs_matrix_pools, ONLY: mpools_create,&
77 mpools_get,&
78 mpools_rebuild_fm_pools,&
79 qs_matrix_pools_type
80 USE qs_mo_types, ONLY: allocate_mo_set,&
81 get_mo_set,&
82 init_mo_set,&
83 mo_set_type,&
84 set_mo_set
85 USE qs_neighbor_list_types, ONLY: get_iterator_info,&
86 get_neighbor_list_set_p,&
87 neighbor_list_iterate,&
88 neighbor_list_iterator_create,&
89 neighbor_list_iterator_p_type,&
90 neighbor_list_iterator_release,&
91 neighbor_list_set_p_type
92 USE scf_control_types, ONLY: smear_type
93 USE smearing_utils, ONLY: smearkp,&
94 smearkp2
95 USE util, ONLY: get_limit
96#include "./base/base_uses.f90"
97
98 IMPLICIT NONE
99
100 PRIVATE
101
102 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpoint_methods'
103
104 PUBLIC :: kpoint_initialize, kpoint_env_initialize, kpoint_initialize_mos, kpoint_initialize_mo_set
105 PUBLIC :: kpoint_init_cell_index, kpoint_set_mo_occupation
106 PUBLIC :: kpoint_density_matrices, kpoint_density_transform
107 PUBLIC :: rskp_transform, lowdin_kp_trans
108
109! **************************************************************************************************
110
111CONTAINS
112
113! **************************************************************************************************
114!> \brief Generate the kpoints and initialize the kpoint environment
115!> \param kpoint The kpoint environment
116!> \param particle_set Particle types and coordinates
117!> \param cell Computational cell information
118! **************************************************************************************************
119 SUBROUTINE kpoint_initialize(kpoint, particle_set, cell)
120
121 TYPE(kpoint_type), POINTER :: kpoint
122 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
123 TYPE(cell_type), POINTER :: cell
124
125 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_initialize'
126
127 INTEGER :: handle, i, ic, ik, iounit, ir, ira, is, &
128 isign, j, natom, nkind, nr, ns
129 INTEGER, ALLOCATABLE, DIMENSION(:) :: atype
130 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: agauge
131 INTEGER, DIMENSION(3, 3) :: frot, krot
132 LOGICAL :: spez
133 REAL(kind=dp) :: eps_kpoint, wsum
134 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coord, scoord
135 REAL(kind=dp), DIMENSION(3) :: diff, kgvec, srot
136 REAL(kind=dp), DIMENSION(3, 3) :: srotmat
137 REAL(kind=dp), DIMENSION(:), POINTER :: wkp_full
138 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp_full
139 TYPE(csym_type) :: crys_sym
140 TYPE(kpoint_sym_type), POINTER :: kpsym
141
142 CALL timeset(routinen, handle)
143
144 cpassert(ASSOCIATED(kpoint))
145
146 SELECT CASE (kpoint%kp_scheme)
147 CASE ("NONE")
148 ! do nothing
149 CASE ("GAMMA")
150 kpoint%nkp = 1
151 ALLOCATE (kpoint%xkp(3, 1), kpoint%wkp(1))
152 kpoint%xkp(1:3, 1) = 0.0_dp
153 kpoint%wkp(1) = 1.0_dp
154 ALLOCATE (kpoint%kp_sym(1))
155 NULLIFY (kpoint%kp_sym(1)%kpoint_sym)
156 CALL kpoint_sym_create(kpoint%kp_sym(1)%kpoint_sym)
157 CASE ("MONKHORST-PACK", "MACDONALD")
158
159 IF (.NOT. kpoint%symmetry) THEN
160 ! we set up a random molecule to avoid any possible symmetry
161 natom = 10
162 ALLOCATE (coord(3, natom), scoord(3, natom), atype(natom))
163 DO i = 1, natom
164 atype(i) = i
165 coord(1, i) = sin(i*0.12345_dp)
166 coord(2, i) = cos(i*0.23456_dp)
167 coord(3, i) = sin(i*0.34567_dp)
168 CALL real_to_scaled(scoord(1:3, i), coord(1:3, i), cell)
169 END DO
170 ELSE
171 natom = SIZE(particle_set)
172 ALLOCATE (scoord(3, natom), atype(natom))
173 DO i = 1, natom
174 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=atype(i))
175 CALL real_to_scaled(scoord(1:3, i), particle_set(i)%r(1:3), cell)
176 END DO
177 END IF
178 IF (kpoint%verbose) THEN
179 iounit = cp_logger_get_default_io_unit()
180 ELSE
181 iounit = -1
182 END IF
183 ! kind type list
184 ALLOCATE (kpoint%atype(natom))
185 kpoint%atype = atype
186 ! Match the periodic gauge used by the k-space matrices.
187 ALLOCATE (agauge(3, natom))
188 DO i = 1, natom
189 agauge(1:3, i) = -floor(scoord(1:3, i) + 0.5_dp - kpoint%eps_geo)
190 END DO
191
192 CALL crys_sym_gen(crys_sym, scoord, atype, cell%hmat, delta=kpoint%eps_geo, iounit=iounit, &
193 use_spglib=kpoint%symmetry)
194 CALL kpoint_gen(crys_sym, kpoint%nkp_grid, symm=kpoint%symmetry, shift=kpoint%kp_shift, &
195 full_grid=kpoint%full_grid, gamma_centered=kpoint%gamma_centered, &
196 inversion_symmetry_only=kpoint%inversion_symmetry_only, &
197 use_spglib_reduction= &
198 kpoint%symmetry_reduction_method == use_spglib_kpoint_symmetry, &
199 use_spglib_backend=kpoint%symmetry_backend == use_spglib_kpoint_backend)
200 kpoint%nkp = crys_sym%nkpoint
201 ALLOCATE (kpoint%xkp(3, kpoint%nkp), kpoint%wkp(kpoint%nkp))
202 wsum = sum(crys_sym%wkpoint)
203 DO ik = 1, kpoint%nkp
204 kpoint%xkp(1:3, ik) = crys_sym%xkpoint(1:3, ik)
205 kpoint%wkp(ik) = crys_sym%wkpoint(ik)/wsum
206 END DO
207
208 eps_kpoint = max(1.e-12_dp, 10.0_dp*kpoint%eps_geo)
209 ! print output
210 IF (kpoint%symmetry) CALL print_crys_symmetry(crys_sym)
211 IF (kpoint%symmetry) CALL print_kp_symmetry(crys_sym)
212
213 ! transfer symmetry information
214 ALLOCATE (kpoint%kp_sym(kpoint%nkp))
215 DO ik = 1, kpoint%nkp
216 NULLIFY (kpoint%kp_sym(ik)%kpoint_sym)
217 CALL kpoint_sym_create(kpoint%kp_sym(ik)%kpoint_sym)
218 kpsym => kpoint%kp_sym(ik)%kpoint_sym
219 IF (crys_sym%nrtot > 0 .AND. .NOT. crys_sym%fullgrid .AND. &
220 crys_sym%istriz == 1 .AND. .NOT. crys_sym%inversion_only) THEN
221 ! set up the symmetrization information
222 kpsym%nwght = nint(crys_sym%wkpoint(ik))
223 ns = kpsym%nwght
224 !
225 IF (ns > 1) THEN
226 DO is = 1, SIZE(crys_sym%kplink, 2)
227 IF (crys_sym%kplink(2, is) == ik) THEN
228 DO ic = 1, crys_sym%nrtot
229 srotmat = matmul(cell%h_inv, matmul(crys_sym%rt(1:3, 1:3, ic), cell%hmat))
230 frot(1:3, 1:3) = nint(srotmat(1:3, 1:3))
231 krot(1:3, 1:3) = nint(transpose(inv_3x3(real(frot(1:3, 1:3), kind=dp))))
232 DO isign = 1, 2
233 ir = merge(crys_sym%ibrot(ic), -crys_sym%ibrot(ic), isign == 1)
234 IF (ir == crys_sym%kpop(is)) cycle
235 kgvec(1:3) = crys_sym%kpmesh(1:3, is) - &
236 matmul(real(merge(krot(1:3, 1:3), -krot(1:3, 1:3), &
237 isign == 1), kind=dp), &
238 kpoint%xkp(1:3, ik))
239 diff(1:3) = kgvec(1:3) - anint(kgvec(1:3))
240 IF (all(abs(diff(1:3)) < eps_kpoint)) ns = ns + 1
241 END DO
242 END DO
243 END IF
244 END DO
245 kpsym%apply_symmetry = .true.
246 natom = SIZE(particle_set)
247 ALLOCATE (kpsym%rot(3, 3, ns))
248 ALLOCATE (kpsym%xkp(3, ns))
249 ALLOCATE (kpsym%rotp(ns))
250 ALLOCATE (kpsym%f0(natom, ns))
251 ALLOCATE (kpsym%fcell(3, natom, ns))
252 ALLOCATE (kpsym%kgphase(natom, ns))
253 nr = 0
254 DO is = 1, SIZE(crys_sym%kplink, 2)
255 IF (crys_sym%kplink(2, is) == ik) THEN
256 nr = nr + 1
257 ir = crys_sym%kpop(is)
258 ira = abs(ir)
259 DO ic = 1, crys_sym%nrtot
260 IF (crys_sym%ibrot(ic) == ira) THEN
261 kpsym%rotp(nr) = ir
262 kpsym%rot(1:3, 1:3, nr) = crys_sym%rt(1:3, 1:3, ic)
263 srotmat = matmul(cell%h_inv, matmul(crys_sym%rt(1:3, 1:3, ic), cell%hmat))
264 frot(1:3, 1:3) = nint(srotmat(1:3, 1:3))
265 kpsym%xkp(1:3, nr) = crys_sym%kpmesh(1:3, is)
266 krot(1:3, 1:3) = nint(transpose(inv_3x3(real(frot(1:3, 1:3), kind=dp))))
267 IF (ir < 0) krot(1:3, 1:3) = -krot(1:3, 1:3)
268 kgvec(1:3) = kpsym%xkp(1:3, nr) - &
269 matmul(real(krot(1:3, 1:3), kind=dp), &
270 kpoint%xkp(1:3, ik))
271 kgvec(1:3) = anint(kgvec(1:3))
272 kpsym%f0(1:natom, nr) = crys_sym%f0(1:natom, ic)
273 DO j = 1, natom
274 srot(1:3) = matmul(srotmat, scoord(1:3, j)) + crys_sym%vt(1:3, ic)
275 kpsym%fcell(1:3, j, nr) = &
276 nint(srot(1:3) - scoord(1:3, kpsym%f0(j, nr))) + &
277 matmul(frot(1:3, 1:3), agauge(1:3, j)) - &
278 agauge(1:3, kpsym%f0(j, nr))
279 kpsym%kgphase(j, nr) = dot_product(kgvec(1:3), &
280 scoord(1:3, j) + &
281 REAL(agauge(1:3, j), kind=dp))
282 END DO
283 EXIT
284 END IF
285 END DO
286 cpassert(ic <= crys_sym%nrtot)
287 END IF
288 END DO
289 DO is = 1, SIZE(crys_sym%kplink, 2)
290 IF (crys_sym%kplink(2, is) == ik) THEN
291 DO ic = 1, crys_sym%nrtot
292 srotmat = matmul(cell%h_inv, matmul(crys_sym%rt(1:3, 1:3, ic), cell%hmat))
293 frot(1:3, 1:3) = nint(srotmat(1:3, 1:3))
294 krot(1:3, 1:3) = nint(transpose(inv_3x3(real(frot(1:3, 1:3), kind=dp))))
295 DO isign = 1, 2
296 ir = merge(crys_sym%ibrot(ic), -crys_sym%ibrot(ic), isign == 1)
297 IF (ir == crys_sym%kpop(is)) cycle
298 kgvec(1:3) = crys_sym%kpmesh(1:3, is) - &
299 matmul(real(merge(krot(1:3, 1:3), -krot(1:3, 1:3), &
300 isign == 1), kind=dp), &
301 kpoint%xkp(1:3, ik))
302 diff(1:3) = kgvec(1:3) - anint(kgvec(1:3))
303 IF (all(abs(diff(1:3)) < eps_kpoint)) THEN
304 nr = nr + 1
305 kpsym%rotp(nr) = ir
306 kpsym%rot(1:3, 1:3, nr) = crys_sym%rt(1:3, 1:3, ic)
307 kpsym%xkp(1:3, nr) = crys_sym%kpmesh(1:3, is)
308 kgvec(1:3) = anint(kgvec(1:3))
309 kpsym%f0(1:natom, nr) = crys_sym%f0(1:natom, ic)
310 DO j = 1, natom
311 srot(1:3) = matmul(srotmat, scoord(1:3, j)) + crys_sym%vt(1:3, ic)
312 kpsym%fcell(1:3, j, nr) = &
313 nint(srot(1:3) - scoord(1:3, kpsym%f0(j, nr))) + &
314 matmul(frot(1:3, 1:3), agauge(1:3, j)) - &
315 agauge(1:3, kpsym%f0(j, nr))
316 kpsym%kgphase(j, nr) = dot_product(kgvec(1:3), &
317 scoord(1:3, j) + &
318 REAL(agauge(1:3, j), kind=dp))
319 END DO
320 END IF
321 END DO
322 END DO
323 END IF
324 END DO
325 kpsym%nwred = nr
326 END IF
327 END IF
328 END DO
329 IF (kpoint%symmetry) THEN
330 nkind = maxval(atype)
331 ns = crys_sym%nrtot
332 ALLOCATE (kpoint%kind_rotmat(ns, nkind))
333 DO i = 1, ns
334 DO j = 1, nkind
335 NULLIFY (kpoint%kind_rotmat(i, j)%rmat)
336 END DO
337 END DO
338 ALLOCATE (kpoint%ibrot(ns))
339 kpoint%ibrot(1:ns) = crys_sym%ibrot(1:ns)
340 END IF
341
342 CALL release_csym_type(crys_sym)
343 DEALLOCATE (scoord, atype)
344 DEALLOCATE (agauge)
345
346 CASE ("GENERAL")
347 NULLIFY (xkp_full, wkp_full)
348 IF (ASSOCIATED(kpoint%xkp_input)) THEN
349 xkp_full => kpoint%xkp_input
350 wkp_full => kpoint%wkp_input
351 ELSE
352 xkp_full => kpoint%xkp
353 wkp_full => kpoint%wkp
354 END IF
355 cpassert(ASSOCIATED(xkp_full))
356 cpassert(ASSOCIATED(wkp_full))
357 IF (.NOT. ASSOCIATED(kpoint%xkp_input)) THEN
358 ALLOCATE (kpoint%xkp_input(3, SIZE(wkp_full)), kpoint%wkp_input(SIZE(wkp_full)))
359 kpoint%xkp_input(1:3, 1:SIZE(wkp_full)) = xkp_full(1:3, 1:SIZE(wkp_full))
360 kpoint%wkp_input(1:SIZE(wkp_full)) = wkp_full(1:SIZE(wkp_full))
361 xkp_full => kpoint%xkp_input
362 wkp_full => kpoint%wkp_input
363 END IF
364 IF (.NOT. kpoint%symmetry) THEN
365 IF (.NOT. ASSOCIATED(kpoint%xkp)) THEN
366 kpoint%nkp = SIZE(wkp_full)
367 ALLOCATE (kpoint%xkp(3, kpoint%nkp), kpoint%wkp(kpoint%nkp))
368 kpoint%xkp(1:3, 1:kpoint%nkp) = xkp_full(1:3, 1:kpoint%nkp)
369 kpoint%wkp(1:kpoint%nkp) = wkp_full(1:kpoint%nkp)
370 END IF
371 ! default: no symmetry settings
372 ALLOCATE (kpoint%kp_sym(kpoint%nkp))
373 DO i = 1, kpoint%nkp
374 NULLIFY (kpoint%kp_sym(i)%kpoint_sym)
375 CALL kpoint_sym_create(kpoint%kp_sym(i)%kpoint_sym)
376 END DO
377 ELSE
378 IF (kpoint%verbose) THEN
379 iounit = cp_logger_get_default_io_unit()
380 ELSE
381 iounit = -1
382 END IF
383 natom = SIZE(particle_set)
384 ALLOCATE (scoord(3, natom), atype(natom))
385 DO i = 1, natom
386 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=atype(i))
387 CALL real_to_scaled(scoord(1:3, i), particle_set(i)%r(1:3), cell)
388 END DO
389 ALLOCATE (kpoint%atype(natom))
390 kpoint%atype = atype
391 ALLOCATE (agauge(3, natom))
392 DO i = 1, natom
393 agauge(1:3, i) = -floor(scoord(1:3, i) + 0.5_dp - kpoint%eps_geo)
394 END DO
395
396 CALL crys_sym_gen(crys_sym, scoord, atype, cell%hmat, delta=kpoint%eps_geo, iounit=iounit, &
397 use_spglib=(kpoint%symmetry_backend == use_spglib_kpoint_backend .OR. &
398 kpoint%symmetry_reduction_method == use_spglib_kpoint_symmetry))
399 CALL kpoint_gen_general(crys_sym, xkp_full, wkp_full, symm=kpoint%symmetry, &
400 full_grid=kpoint%full_grid, &
401 inversion_symmetry_only=kpoint%inversion_symmetry_only, &
402 use_spglib_reduction= &
403 kpoint%symmetry_reduction_method == use_spglib_kpoint_symmetry, &
404 use_spglib_backend=kpoint%symmetry_backend == use_spglib_kpoint_backend)
405 IF (ASSOCIATED(kpoint%xkp)) THEN
406 DEALLOCATE (kpoint%xkp)
407 NULLIFY (kpoint%xkp)
408 END IF
409 IF (ASSOCIATED(kpoint%wkp)) THEN
410 DEALLOCATE (kpoint%wkp)
411 NULLIFY (kpoint%wkp)
412 END IF
413 kpoint%nkp = crys_sym%nkpoint
414 ALLOCATE (kpoint%xkp(3, kpoint%nkp), kpoint%wkp(kpoint%nkp))
415 wsum = sum(crys_sym%wkpoint)
416 DO ik = 1, kpoint%nkp
417 kpoint%xkp(1:3, ik) = crys_sym%xkpoint(1:3, ik)
418 kpoint%wkp(ik) = crys_sym%wkpoint(ik)/wsum
419 END DO
420
421 eps_kpoint = max(1.e-12_dp, 10.0_dp*kpoint%eps_geo)
422 CALL print_crys_symmetry(crys_sym)
423 CALL print_kp_symmetry(crys_sym)
424
425 ALLOCATE (kpoint%kp_sym(kpoint%nkp))
426 DO ik = 1, kpoint%nkp
427 NULLIFY (kpoint%kp_sym(ik)%kpoint_sym)
428 CALL kpoint_sym_create(kpoint%kp_sym(ik)%kpoint_sym)
429 kpsym => kpoint%kp_sym(ik)%kpoint_sym
430 IF (crys_sym%nrtot > 0 .AND. .NOT. crys_sym%fullgrid .AND. &
431 crys_sym%istriz == 1 .AND. .NOT. crys_sym%inversion_only) THEN
432 kpsym%nwght = nint(crys_sym%wkpoint(ik))
433 ns = kpsym%nwght
434 IF (ns > 1) THEN
435 DO is = 1, SIZE(crys_sym%kplink, 2)
436 IF (crys_sym%kplink(2, is) == ik) THEN
437 DO ic = 1, crys_sym%nrtot
438 srotmat = matmul(cell%h_inv, matmul(crys_sym%rt(1:3, 1:3, ic), cell%hmat))
439 frot(1:3, 1:3) = nint(srotmat(1:3, 1:3))
440 krot(1:3, 1:3) = nint(transpose(inv_3x3(real(frot(1:3, 1:3), kind=dp))))
441 DO isign = 1, 2
442 ir = merge(crys_sym%ibrot(ic), -crys_sym%ibrot(ic), isign == 1)
443 IF (ir == crys_sym%kpop(is)) cycle
444 kgvec(1:3) = crys_sym%kpmesh(1:3, is) - &
445 matmul(real(merge(krot(1:3, 1:3), -krot(1:3, 1:3), &
446 isign == 1), kind=dp), &
447 kpoint%xkp(1:3, ik))
448 diff(1:3) = kgvec(1:3) - anint(kgvec(1:3))
449 IF (all(abs(diff(1:3)) < eps_kpoint)) ns = ns + 1
450 END DO
451 END DO
452 END IF
453 END DO
454 kpsym%apply_symmetry = .true.
455 ALLOCATE (kpsym%rot(3, 3, ns))
456 ALLOCATE (kpsym%xkp(3, ns))
457 ALLOCATE (kpsym%rotp(ns))
458 ALLOCATE (kpsym%f0(natom, ns))
459 ALLOCATE (kpsym%fcell(3, natom, ns))
460 ALLOCATE (kpsym%kgphase(natom, ns))
461 nr = 0
462 DO is = 1, SIZE(crys_sym%kplink, 2)
463 IF (crys_sym%kplink(2, is) == ik) THEN
464 nr = nr + 1
465 ir = crys_sym%kpop(is)
466 ira = abs(ir)
467 DO ic = 1, crys_sym%nrtot
468 IF (crys_sym%ibrot(ic) == ira) THEN
469 kpsym%rotp(nr) = ir
470 kpsym%rot(1:3, 1:3, nr) = crys_sym%rt(1:3, 1:3, ic)
471 srotmat = matmul(cell%h_inv, matmul(crys_sym%rt(1:3, 1:3, ic), cell%hmat))
472 frot(1:3, 1:3) = nint(srotmat(1:3, 1:3))
473 kpsym%xkp(1:3, nr) = crys_sym%kpmesh(1:3, is)
474 krot(1:3, 1:3) = nint(transpose(inv_3x3(real(frot(1:3, 1:3), kind=dp))))
475 IF (ir < 0) krot(1:3, 1:3) = -krot(1:3, 1:3)
476 kgvec(1:3) = kpsym%xkp(1:3, nr) - &
477 matmul(real(krot(1:3, 1:3), kind=dp), &
478 kpoint%xkp(1:3, ik))
479 kgvec(1:3) = anint(kgvec(1:3))
480 kpsym%f0(1:natom, nr) = crys_sym%f0(1:natom, ic)
481 DO j = 1, natom
482 srot(1:3) = matmul(srotmat, scoord(1:3, j)) + crys_sym%vt(1:3, ic)
483 kpsym%fcell(1:3, j, nr) = &
484 nint(srot(1:3) - scoord(1:3, kpsym%f0(j, nr))) + &
485 matmul(frot(1:3, 1:3), agauge(1:3, j)) - &
486 agauge(1:3, kpsym%f0(j, nr))
487 kpsym%kgphase(j, nr) = dot_product(kgvec(1:3), &
488 scoord(1:3, j) + &
489 REAL(agauge(1:3, j), kind=dp))
490 END DO
491 EXIT
492 END IF
493 END DO
494 cpassert(ic <= crys_sym%nrtot)
495 END IF
496 END DO
497 DO is = 1, SIZE(crys_sym%kplink, 2)
498 IF (crys_sym%kplink(2, is) == ik) THEN
499 DO ic = 1, crys_sym%nrtot
500 srotmat = matmul(cell%h_inv, matmul(crys_sym%rt(1:3, 1:3, ic), cell%hmat))
501 frot(1:3, 1:3) = nint(srotmat(1:3, 1:3))
502 krot(1:3, 1:3) = nint(transpose(inv_3x3(real(frot(1:3, 1:3), kind=dp))))
503 DO isign = 1, 2
504 ir = merge(crys_sym%ibrot(ic), -crys_sym%ibrot(ic), isign == 1)
505 IF (ir == crys_sym%kpop(is)) cycle
506 kgvec(1:3) = crys_sym%kpmesh(1:3, is) - &
507 matmul(real(merge(krot(1:3, 1:3), -krot(1:3, 1:3), &
508 isign == 1), kind=dp), &
509 kpoint%xkp(1:3, ik))
510 diff(1:3) = kgvec(1:3) - anint(kgvec(1:3))
511 IF (all(abs(diff(1:3)) < eps_kpoint)) THEN
512 nr = nr + 1
513 kpsym%rotp(nr) = ir
514 kpsym%rot(1:3, 1:3, nr) = crys_sym%rt(1:3, 1:3, ic)
515 kpsym%xkp(1:3, nr) = crys_sym%kpmesh(1:3, is)
516 kgvec(1:3) = anint(kgvec(1:3))
517 kpsym%f0(1:natom, nr) = crys_sym%f0(1:natom, ic)
518 DO j = 1, natom
519 srot(1:3) = matmul(srotmat, scoord(1:3, j)) + crys_sym%vt(1:3, ic)
520 kpsym%fcell(1:3, j, nr) = &
521 nint(srot(1:3) - scoord(1:3, kpsym%f0(j, nr))) + &
522 matmul(frot(1:3, 1:3), agauge(1:3, j)) - &
523 agauge(1:3, kpsym%f0(j, nr))
524 kpsym%kgphase(j, nr) = dot_product(kgvec(1:3), &
525 scoord(1:3, j) + &
526 REAL(agauge(1:3, j), kind=dp))
527 END DO
528 END IF
529 END DO
530 END DO
531 END IF
532 END DO
533 kpsym%nwred = nr
534 END IF
535 END IF
536 END DO
537 nkind = maxval(atype)
538 ns = crys_sym%nrtot
539 ALLOCATE (kpoint%kind_rotmat(ns, nkind))
540 DO i = 1, ns
541 DO j = 1, nkind
542 NULLIFY (kpoint%kind_rotmat(i, j)%rmat)
543 END DO
544 END DO
545 ALLOCATE (kpoint%ibrot(ns))
546 kpoint%ibrot(1:ns) = crys_sym%ibrot(1:ns)
547
548 CALL release_csym_type(crys_sym)
549 DEALLOCATE (scoord, atype)
550 DEALLOCATE (agauge)
551 END IF
552 CASE DEFAULT
553 cpassert(.false.)
554 END SELECT
555
556 ! check for consistency of options
557 SELECT CASE (kpoint%kp_scheme)
558 CASE ("NONE")
559 ! don't use k-point code
560 CASE ("GAMMA")
561 cpassert(kpoint%nkp == 1)
562 cpassert(sum(abs(kpoint%xkp)) <= 1.e-12_dp)
563 cpassert(kpoint%wkp(1) == 1.0_dp)
564 cpassert(.NOT. kpoint%symmetry)
565 CASE ("GENERAL")
566 cpassert(kpoint%nkp >= 1)
567 CASE ("MONKHORST-PACK", "MACDONALD")
568 cpassert(kpoint%nkp >= 1)
569 END SELECT
570 IF (kpoint%use_real_wfn) THEN
571 ! what about inversion symmetry?
572 ikloop: DO ik = 1, kpoint%nkp
573 DO i = 1, 3
574 spez = (kpoint%xkp(i, ik) == 0.0_dp .OR. kpoint%xkp(i, ik) == 0.5_dp)
575 IF (.NOT. spez) EXIT ikloop
576 END DO
577 END DO ikloop
578 IF (.NOT. spez) THEN
579 ! Warning: real wfn might be wrong for this system
580 CALL cp_warn(__location__, &
581 "A calculation using real wavefunctions is requested. "// &
582 "We could not determine if the symmetry of the system allows real wavefunctions. ")
583 END IF
584 END IF
585
586 CALL timestop(handle)
587
588 END SUBROUTINE kpoint_initialize
589
590! **************************************************************************************************
591!> \brief Initialize the kpoint environment
592!> \param kpoint Kpoint environment
593!> \param para_env ...
594!> \param blacs_env ...
595!> \param with_aux_fit ...
596! **************************************************************************************************
597 SUBROUTINE kpoint_env_initialize(kpoint, para_env, blacs_env, with_aux_fit)
598
599 TYPE(kpoint_type), INTENT(INOUT) :: kpoint
600 TYPE(mp_para_env_type), INTENT(IN), TARGET :: para_env
601 TYPE(cp_blacs_env_type), INTENT(IN), TARGET :: blacs_env
602 LOGICAL, INTENT(IN), OPTIONAL :: with_aux_fit
603
604 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_env_initialize'
605
606 INTEGER :: handle, igr, ik, ikk, ngr, niogrp, nkp, &
607 nkp_grp, nkp_loc, npe, unit_nr
608 INTEGER, DIMENSION(2) :: dims, pos
609 LOGICAL :: aux_fit
610 TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_aux_env, kp_env
611 TYPE(kpoint_env_type), POINTER :: kp
612 TYPE(mp_cart_type) :: comm_cart
613 TYPE(mp_para_env_type), POINTER :: para_env_inter_kp, para_env_kp
614
615 CALL timeset(routinen, handle)
616
617 IF (PRESENT(with_aux_fit)) THEN
618 aux_fit = with_aux_fit
619 ELSE
620 aux_fit = .false.
621 END IF
622
623 kpoint%para_env => para_env
624 CALL kpoint%para_env%retain()
625 kpoint%blacs_env_all => blacs_env
626 CALL kpoint%blacs_env_all%retain()
627
628 cpassert(.NOT. ASSOCIATED(kpoint%kp_env))
629 IF (aux_fit) THEN
630 cpassert(.NOT. ASSOCIATED(kpoint%kp_aux_env))
631 END IF
632
633 NULLIFY (kp_env, kp_aux_env)
634 nkp = kpoint%nkp
635 npe = para_env%num_pe
636 IF (npe == 1) THEN
637 ! only one process available -> owns all kpoints
638 ALLOCATE (kp_env(nkp))
639 DO ik = 1, nkp
640 NULLIFY (kp_env(ik)%kpoint_env)
641 CALL kpoint_env_create(kp_env(ik)%kpoint_env)
642 kp => kp_env(ik)%kpoint_env
643 kp%nkpoint = ik
644 kp%wkp = kpoint%wkp(ik)
645 kp%xkp(1:3) = kpoint%xkp(1:3, ik)
646 kp%is_local = .true.
647 END DO
648 kpoint%kp_env => kp_env
649
650 IF (aux_fit) THEN
651 ALLOCATE (kp_aux_env(nkp))
652 DO ik = 1, nkp
653 NULLIFY (kp_aux_env(ik)%kpoint_env)
654 CALL kpoint_env_create(kp_aux_env(ik)%kpoint_env)
655 kp => kp_aux_env(ik)%kpoint_env
656 kp%nkpoint = ik
657 kp%wkp = kpoint%wkp(ik)
658 kp%xkp(1:3) = kpoint%xkp(1:3, ik)
659 kp%is_local = .true.
660 END DO
661
662 kpoint%kp_aux_env => kp_aux_env
663 END IF
664
665 ALLOCATE (kpoint%kp_dist(2, 1))
666 kpoint%kp_dist(1, 1) = 1
667 kpoint%kp_dist(2, 1) = nkp
668 kpoint%kp_range(1) = 1
669 kpoint%kp_range(2) = nkp
670
671 ! parallel environments
672 kpoint%para_env_kp => para_env
673 CALL kpoint%para_env_kp%retain()
674 kpoint%para_env_inter_kp => para_env
675 CALL kpoint%para_env_inter_kp%retain()
676 kpoint%iogrp = .true.
677 kpoint%nkp_groups = 1
678 ELSE
679 IF (kpoint%parallel_group_size == -1) THEN
680 ! maximum parallelization over kpoints
681 ! making sure that the group size divides the npe and the nkp_grp the nkp
682 ! in the worst case, there will be no parallelism over kpoints.
683 DO igr = npe, 1, -1
684 IF (mod(npe, igr) /= 0) cycle
685 nkp_grp = npe/igr
686 IF (mod(nkp, nkp_grp) /= 0) cycle
687 ngr = igr
688 END DO
689 ELSE IF (kpoint%parallel_group_size == 0) THEN
690 ! no parallelization over kpoints
691 ngr = npe
692 ELSE IF (kpoint%parallel_group_size > 0) THEN
693 ngr = min(kpoint%parallel_group_size, npe)
694 ELSE
695 cpassert(.false.)
696 END IF
697 nkp_grp = npe/ngr
698 ! processor dimensions
699 dims(1) = ngr
700 dims(2) = nkp_grp
701 cpassert(mod(nkp, nkp_grp) == 0)
702 nkp_loc = nkp/nkp_grp
703
704 IF ((dims(1)*dims(2) /= npe)) THEN
705 cpabort("Number of processors is not divisible by the kpoint group size.")
706 END IF
707
708 ! Create the subgroups, one for each k-point group and one interconnecting group
709 CALL comm_cart%create(comm_old=para_env, ndims=2, dims=dims)
710 pos = comm_cart%mepos_cart
711 ALLOCATE (para_env_kp)
712 CALL para_env_kp%from_split(comm_cart, pos(2))
713 ALLOCATE (para_env_inter_kp)
714 CALL para_env_inter_kp%from_split(comm_cart, pos(1))
715 CALL comm_cart%free()
716
717 niogrp = 0
718 IF (para_env%is_source()) niogrp = 1
719 CALL para_env_kp%sum(niogrp)
720 kpoint%iogrp = (niogrp == 1)
721
722 ! parallel groups
723 kpoint%para_env_kp => para_env_kp
724 kpoint%para_env_inter_kp => para_env_inter_kp
725
726 ! distribution of kpoints
727 ALLOCATE (kpoint%kp_dist(2, nkp_grp))
728 DO igr = 1, nkp_grp
729 kpoint%kp_dist(1:2, igr) = get_limit(nkp, nkp_grp, igr - 1)
730 END DO
731 ! local kpoints
732 kpoint%kp_range(1:2) = kpoint%kp_dist(1:2, para_env_inter_kp%mepos + 1)
733
734 ALLOCATE (kp_env(nkp_loc))
735 DO ik = 1, nkp_loc
736 NULLIFY (kp_env(ik)%kpoint_env)
737 ikk = kpoint%kp_range(1) + ik - 1
738 CALL kpoint_env_create(kp_env(ik)%kpoint_env)
739 kp => kp_env(ik)%kpoint_env
740 kp%nkpoint = ikk
741 kp%wkp = kpoint%wkp(ikk)
742 kp%xkp(1:3) = kpoint%xkp(1:3, ikk)
743 kp%is_local = (ngr == 1)
744 END DO
745 kpoint%kp_env => kp_env
746
747 IF (aux_fit) THEN
748 ALLOCATE (kp_aux_env(nkp_loc))
749 DO ik = 1, nkp_loc
750 NULLIFY (kp_aux_env(ik)%kpoint_env)
751 ikk = kpoint%kp_range(1) + ik - 1
752 CALL kpoint_env_create(kp_aux_env(ik)%kpoint_env)
753 kp => kp_aux_env(ik)%kpoint_env
754 kp%nkpoint = ikk
755 kp%wkp = kpoint%wkp(ikk)
756 kp%xkp(1:3) = kpoint%xkp(1:3, ikk)
757 kp%is_local = (ngr == 1)
758 END DO
759 kpoint%kp_aux_env => kp_aux_env
760 END IF
761
762 unit_nr = cp_logger_get_default_io_unit()
763
764 IF (unit_nr > 0 .AND. kpoint%verbose) THEN
765 WRITE (unit_nr, *)
766 WRITE (unit_nr, fmt="(T2,A,T71,I10)") "KPOINTS| Number of kpoint groups ", nkp_grp
767 WRITE (unit_nr, fmt="(T2,A,T71,I10)") "KPOINTS| Size of each kpoint group", ngr
768 WRITE (unit_nr, fmt="(T2,A,T71,I10)") "KPOINTS| Number of kpoints per group", nkp_loc
769 END IF
770 kpoint%nkp_groups = nkp_grp
771
772 END IF
773
774 CALL timestop(handle)
775
776 END SUBROUTINE kpoint_env_initialize
777
778! **************************************************************************************************
779!> \brief Initialize a set of MOs and density matrix for each kpoint (kpoint group)
780!> \param kpoint Kpoint environment
781!> \param mos Reference MOs (global)
782!> \param added_mos ...
783!> \param for_aux_fit ...
784! **************************************************************************************************
785 SUBROUTINE kpoint_initialize_mos(kpoint, mos, added_mos, for_aux_fit)
786
787 TYPE(kpoint_type), POINTER :: kpoint
788 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
789 INTEGER, INTENT(IN), OPTIONAL :: added_mos
790 LOGICAL, OPTIONAL :: for_aux_fit
791
792 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_initialize_mos'
793
794 INTEGER :: handle, ic, ik, is, nadd, nao, nc, &
795 nelectron, nkp_loc, nmo, nmorig(2), &
796 nspin
797 LOGICAL :: aux_fit
798 REAL(kind=dp) :: flexible_electron_count, maxocc, n_el_f
799 TYPE(cp_blacs_env_type), POINTER :: blacs_env
800 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools
801 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
802 TYPE(cp_fm_type), POINTER :: fmlocal
803 TYPE(kpoint_env_type), POINTER :: kp
804 TYPE(qs_matrix_pools_type), POINTER :: mpools
805
806 CALL timeset(routinen, handle)
807
808 IF (PRESENT(for_aux_fit)) THEN
809 aux_fit = for_aux_fit
810 ELSE
811 aux_fit = .false.
812 END IF
813
814 cpassert(ASSOCIATED(kpoint))
815
816 IF (.true. .OR. ASSOCIATED(mos(1)%mo_coeff)) THEN
817 IF (aux_fit) THEN
818 cpassert(ASSOCIATED(kpoint%kp_aux_env))
819 END IF
820
821 IF (PRESENT(added_mos)) THEN
822 nadd = added_mos
823 ELSE
824 nadd = 0
825 END IF
826
827 IF (kpoint%use_real_wfn) THEN
828 nc = 1
829 ELSE
830 nc = 2
831 END IF
832 nspin = SIZE(mos, 1)
833 nkp_loc = kpoint%kp_range(2) - kpoint%kp_range(1) + 1
834 IF (nkp_loc > 0) THEN
835 IF (aux_fit) THEN
836 cpassert(SIZE(kpoint%kp_aux_env) == nkp_loc)
837 ELSE
838 cpassert(SIZE(kpoint%kp_env) == nkp_loc)
839 END IF
840 ! allocate the mo sets, correct number of kpoints (local), real/complex, spin
841 DO ik = 1, nkp_loc
842 IF (aux_fit) THEN
843 kp => kpoint%kp_aux_env(ik)%kpoint_env
844 ELSE
845 kp => kpoint%kp_env(ik)%kpoint_env
846 END IF
847 ALLOCATE (kp%mos(nc, nspin))
848 DO is = 1, nspin
849 CALL get_mo_set(mos(is), nao=nao, nmo=nmo, nelectron=nelectron, &
850 n_el_f=n_el_f, maxocc=maxocc, flexible_electron_count=flexible_electron_count)
851 nmo = min(nao, nmo + nadd)
852 DO ic = 1, nc
853 CALL allocate_mo_set(kp%mos(ic, is), nao, nmo, nelectron, n_el_f, maxocc, &
854 flexible_electron_count)
855 END DO
856 END DO
857 END DO
858
859 ! generate the blacs environment for the kpoint group
860 ! we generate a blacs env for each kpoint group in parallel
861 ! we assume here that the group para_env_inter_kp will connect
862 ! equivalent parts of fm matrices, i.e. no reshuffeling of processors
863 NULLIFY (blacs_env)
864 IF (ASSOCIATED(kpoint%blacs_env)) THEN
865 blacs_env => kpoint%blacs_env
866 ELSE
867 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=kpoint%para_env_kp)
868 kpoint%blacs_env => blacs_env
869 END IF
870
871 ! set possible new number of MOs
872 DO is = 1, nspin
873 CALL get_mo_set(mos(is), nmo=nmorig(is))
874 nmo = min(nao, nmorig(is) + nadd)
875 CALL set_mo_set(mos(is), nmo=nmo)
876 END DO
877 ! matrix pools for the kpoint group, information on MOs is transferred using
878 ! generic mos structure
879 NULLIFY (mpools)
880 CALL mpools_create(mpools=mpools)
881 CALL mpools_rebuild_fm_pools(mpools=mpools, mos=mos, &
882 blacs_env=blacs_env, para_env=kpoint%para_env_kp)
883
884 IF (aux_fit) THEN
885 kpoint%mpools_aux_fit => mpools
886 ELSE
887 kpoint%mpools => mpools
888 END IF
889
890 ! reset old number of MOs
891 DO is = 1, nspin
892 CALL set_mo_set(mos(is), nmo=nmorig(is))
893 END DO
894
895 ! allocate density matrices
896 CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
897 ALLOCATE (fmlocal)
898 CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
899 CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
900 DO ik = 1, nkp_loc
901 IF (aux_fit) THEN
902 kp => kpoint%kp_aux_env(ik)%kpoint_env
903 ELSE
904 kp => kpoint%kp_env(ik)%kpoint_env
905 END IF
906 ! density matrix
907 CALL cp_fm_release(kp%pmat)
908 ALLOCATE (kp%pmat(nc, nspin))
909 DO is = 1, nspin
910 DO ic = 1, nc
911 CALL cp_fm_create(kp%pmat(ic, is), matrix_struct)
912 END DO
913 END DO
914 ! energy weighted density matrix
915 CALL cp_fm_release(kp%wmat)
916 ALLOCATE (kp%wmat(nc, nspin))
917 DO is = 1, nspin
918 DO ic = 1, nc
919 CALL cp_fm_create(kp%wmat(ic, is), matrix_struct)
920 END DO
921 END DO
922 END DO
923 CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
924 DEALLOCATE (fmlocal)
925
926 END IF
927
928 END IF
929
930 CALL timestop(handle)
931
932 END SUBROUTINE kpoint_initialize_mos
933
934! **************************************************************************************************
935!> \brief ...
936!> \param kpoint ...
937! **************************************************************************************************
938 SUBROUTINE kpoint_initialize_mo_set(kpoint)
939 TYPE(kpoint_type), POINTER :: kpoint
940
941 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_initialize_mo_set'
942
943 INTEGER :: handle, ic, ik, ikk, ispin
944 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
945 TYPE(cp_fm_type), POINTER :: mo_coeff
946 TYPE(mo_set_type), DIMENSION(:, :), POINTER :: moskp
947
948 CALL timeset(routinen, handle)
949
950 DO ik = 1, SIZE(kpoint%kp_env)
951 CALL mpools_get(kpoint%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
952 moskp => kpoint%kp_env(ik)%kpoint_env%mos
953 ikk = kpoint%kp_range(1) + ik - 1
954 cpassert(ASSOCIATED(moskp))
955 DO ispin = 1, SIZE(moskp, 2)
956 DO ic = 1, SIZE(moskp, 1)
957 CALL get_mo_set(moskp(ic, ispin), mo_coeff=mo_coeff)
958 IF (.NOT. ASSOCIATED(mo_coeff)) THEN
959 CALL init_mo_set(moskp(ic, ispin), &
960 fm_pool=ao_mo_fm_pools(ispin)%pool, name="kpoints")
961 END IF
962 END DO
963 END DO
964 END DO
965
966 CALL timestop(handle)
967
968 END SUBROUTINE kpoint_initialize_mo_set
969
970! **************************************************************************************************
971!> \brief Generates the mapping of cell indices and linear RS index
972!> CELL (0,0,0) is always mapped to index 1
973!> \param kpoint Kpoint environment
974!> \param sab_nl Defining neighbour list
975!> \param para_env Parallel environment
976!> \param nimages [output]
977! **************************************************************************************************
978 SUBROUTINE kpoint_init_cell_index(kpoint, sab_nl, para_env, nimages)
979
980 TYPE(kpoint_type), POINTER :: kpoint
981 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
982 POINTER :: sab_nl
983 TYPE(mp_para_env_type), POINTER :: para_env
984 INTEGER, INTENT(OUT) :: nimages
985
986 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_init_cell_index'
987
988 INTEGER :: handle, i1, i2, i3, ic, icount, it, &
989 ncount
990 INTEGER, DIMENSION(3) :: cell, itm
991 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell, list
992 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index, cti
993 LOGICAL :: new
994 TYPE(neighbor_list_iterator_p_type), &
995 DIMENSION(:), POINTER :: nl_iterator
996
997 NULLIFY (cell_to_index, index_to_cell)
998
999 CALL timeset(routinen, handle)
1000
1001 cpassert(ASSOCIATED(kpoint))
1002
1003 ALLOCATE (list(3, 125))
1004 list = 0
1005 icount = 1
1006
1007 CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
1008 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1009 CALL get_iterator_info(nl_iterator, cell=cell)
1010
1011 new = .true.
1012 DO ic = 1, icount
1013 IF (cell(1) == list(1, ic) .AND. cell(2) == list(2, ic) .AND. &
1014 cell(3) == list(3, ic)) THEN
1015 new = .false.
1016 EXIT
1017 END IF
1018 END DO
1019 IF (new) THEN
1020 icount = icount + 1
1021 IF (icount > SIZE(list, 2)) THEN
1022 CALL reallocate(list, 1, 3, 1, 2*SIZE(list, 2))
1023 END IF
1024 list(1:3, icount) = cell(1:3)
1025 END IF
1026
1027 END DO
1028 CALL neighbor_list_iterator_release(nl_iterator)
1029
1030 itm(1) = maxval(abs(list(1, 1:icount)))
1031 itm(2) = maxval(abs(list(2, 1:icount)))
1032 itm(3) = maxval(abs(list(3, 1:icount)))
1033 CALL para_env%max(itm)
1034 it = maxval(itm(1:3))
1035 IF (ASSOCIATED(kpoint%cell_to_index)) THEN
1036 DEALLOCATE (kpoint%cell_to_index)
1037 END IF
1038 ALLOCATE (kpoint%cell_to_index(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
1039 cell_to_index => kpoint%cell_to_index
1040 cti => cell_to_index
1041 cti(:, :, :) = 0
1042 DO ic = 1, icount
1043 i1 = list(1, ic)
1044 i2 = list(2, ic)
1045 i3 = list(3, ic)
1046 cti(i1, i2, i3) = ic
1047 END DO
1048 CALL para_env%sum(cti)
1049 ncount = 0
1050 DO i1 = -itm(1), itm(1)
1051 DO i2 = -itm(2), itm(2)
1052 DO i3 = -itm(3), itm(3)
1053 IF (cti(i1, i2, i3) == 0) THEN
1054 cti(i1, i2, i3) = 1000000
1055 ELSE
1056 ncount = ncount + 1
1057 cti(i1, i2, i3) = (abs(i1) + abs(i2) + abs(i3))*1000 + abs(i3)*100 + abs(i2)*10 + abs(i1)
1058 cti(i1, i2, i3) = cti(i1, i2, i3) + (i1 + i2 + i3)
1059 END IF
1060 END DO
1061 END DO
1062 END DO
1063
1064 IF (ASSOCIATED(kpoint%index_to_cell)) THEN
1065 DEALLOCATE (kpoint%index_to_cell)
1066 END IF
1067 ALLOCATE (kpoint%index_to_cell(3, ncount))
1068 index_to_cell => kpoint%index_to_cell
1069 DO ic = 1, ncount
1070 cell = minloc(cti)
1071 i1 = cell(1) - 1 - itm(1)
1072 i2 = cell(2) - 1 - itm(2)
1073 i3 = cell(3) - 1 - itm(3)
1074 cti(i1, i2, i3) = 1000000
1075 index_to_cell(1, ic) = i1
1076 index_to_cell(2, ic) = i2
1077 index_to_cell(3, ic) = i3
1078 END DO
1079 cti(:, :, :) = 0
1080 DO ic = 1, ncount
1081 i1 = index_to_cell(1, ic)
1082 i2 = index_to_cell(2, ic)
1083 i3 = index_to_cell(3, ic)
1084 cti(i1, i2, i3) = ic
1085 END DO
1086
1087 ! keep pointer to this neighborlist
1088 kpoint%sab_nl => sab_nl
1089
1090 ! set number of images
1091 nimages = SIZE(index_to_cell, 2)
1092
1093 DEALLOCATE (list)
1094
1095 CALL timestop(handle)
1096
1097 END SUBROUTINE kpoint_init_cell_index
1098
1099! **************************************************************************************************
1100!> \brief Transformation of real space matrices to a kpoint
1101!> \param rmatrix Real part of kpoint matrix
1102!> \param cmatrix Complex part of kpoint matrix (optional)
1103!> \param rsmat Real space matrices
1104!> \param ispin Spin index
1105!> \param xkp Kpoint coordinates
1106!> \param cell_to_index mapping of cell indices to RS index
1107!> \param sab_nl Defining neighbor list
1108!> \param is_complex Matrix to be transformed is imaginary
1109!> \param rs_sign Matrix to be transformed is csaled by rs_sign
1110! **************************************************************************************************
1111 SUBROUTINE rskp_transform(rmatrix, cmatrix, rsmat, ispin, &
1112 xkp, cell_to_index, sab_nl, is_complex, rs_sign)
1113
1114 TYPE(dbcsr_type) :: rmatrix
1115 TYPE(dbcsr_type), OPTIONAL :: cmatrix
1116 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
1117 INTEGER, INTENT(IN) :: ispin
1118 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: xkp
1119 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1120 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1121 POINTER :: sab_nl
1122 LOGICAL, INTENT(IN), OPTIONAL :: is_complex
1123 REAL(kind=dp), INTENT(IN), OPTIONAL :: rs_sign
1124
1125 CHARACTER(LEN=*), PARAMETER :: routinen = 'rskp_transform'
1126
1127 INTEGER :: handle, iatom, ic, icol, irow, jatom, &
1128 nimg
1129 INTEGER, DIMENSION(3) :: cell
1130 LOGICAL :: do_symmetric, found, my_complex, &
1131 wfn_real_only
1132 REAL(kind=dp) :: arg, coskl, fsign, fsym, sinkl
1133 REAL(kind=dp), DIMENSION(:, :), POINTER :: cblock, rblock, rsblock
1134 TYPE(neighbor_list_iterator_p_type), &
1135 DIMENSION(:), POINTER :: nl_iterator
1136
1137 CALL timeset(routinen, handle)
1138
1139 my_complex = .false.
1140 IF (PRESENT(is_complex)) my_complex = is_complex
1141
1142 fsign = 1.0_dp
1143 IF (PRESENT(rs_sign)) fsign = rs_sign
1144
1145 wfn_real_only = .true.
1146 IF (PRESENT(cmatrix)) wfn_real_only = .false.
1147
1148 nimg = SIZE(rsmat, 2)
1149
1150 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1151
1152 CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
1153 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1154 CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
1155
1156 ! fsym = +- 1 is due to real space matrices being non-symmetric (although in a symmtric type)
1157 ! with the link S_mu^0,nu^b = S_nu^0,mu^-b, and the KP matrices beeing Hermitian
1158 fsym = 1.0_dp
1159 irow = iatom
1160 icol = jatom
1161 IF (do_symmetric .AND. (iatom > jatom)) THEN
1162 irow = jatom
1163 icol = iatom
1164 fsym = -1.0_dp
1165 END IF
1166
1167 ic = cell_to_index(cell(1), cell(2), cell(3))
1168 IF (ic < 1 .OR. ic > nimg) cycle
1169
1170 arg = real(cell(1), dp)*xkp(1) + real(cell(2), dp)*xkp(2) + real(cell(3), dp)*xkp(3)
1171 IF (my_complex) THEN
1172 coskl = fsign*fsym*cos(twopi*arg)
1173 sinkl = fsign*sin(twopi*arg)
1174 ELSE
1175 coskl = fsign*cos(twopi*arg)
1176 sinkl = fsign*fsym*sin(twopi*arg)
1177 END IF
1178
1179 CALL dbcsr_get_block_p(matrix=rsmat(ispin, ic)%matrix, row=irow, col=icol, &
1180 block=rsblock, found=found)
1181 IF (.NOT. found) cycle
1182
1183 IF (wfn_real_only) THEN
1184 CALL dbcsr_get_block_p(matrix=rmatrix, row=irow, col=icol, &
1185 block=rblock, found=found)
1186 IF (.NOT. found) cycle
1187 rblock = rblock + coskl*rsblock
1188 ELSE
1189 CALL dbcsr_get_block_p(matrix=rmatrix, row=irow, col=icol, &
1190 block=rblock, found=found)
1191 IF (.NOT. found) cycle
1192 CALL dbcsr_get_block_p(matrix=cmatrix, row=irow, col=icol, &
1193 block=cblock, found=found)
1194 IF (.NOT. found) cycle
1195 rblock = rblock + coskl*rsblock
1196 cblock = cblock + sinkl*rsblock
1197 END IF
1198
1199 END DO
1200 CALL neighbor_list_iterator_release(nl_iterator)
1201
1202 CALL timestop(handle)
1203
1204 END SUBROUTINE rskp_transform
1205
1206! **************************************************************************************************
1207!> \brief Given the eigenvalues of all kpoints, calculates the occupation numbers
1208!> \param kpoint Kpoint environment
1209!> \param smear Smearing information
1210!> \param probe ...
1211! **************************************************************************************************
1212 SUBROUTINE kpoint_set_mo_occupation(kpoint, smear, probe)
1213
1214 TYPE(kpoint_type), POINTER :: kpoint
1215 TYPE(smear_type) :: smear
1216 TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
1217 POINTER :: probe
1218
1219 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_set_mo_occupation'
1220
1221 INTEGER :: handle, ik, ikpgr, ispin, kplocal, nao, &
1222 nb, ncol_global, ne_a, ne_b, &
1223 nelectron, nkp, nmo, nrow_global, nspin
1224 INTEGER, DIMENSION(2) :: kp_range
1225 REAL(kind=dp) :: kts, mu, mus(2), nel
1226 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: smatrix
1227 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: weig, wocc
1228 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: icoeff, rcoeff
1229 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues, occupation, wkp
1230 TYPE(cp_fm_type), POINTER :: mo_coeff
1231 TYPE(kpoint_env_type), POINTER :: kp
1232 TYPE(mo_set_type), POINTER :: mo_set
1233 TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
1234
1235 CALL timeset(routinen, handle)
1236
1237 ! first collect all the eigenvalues
1238 CALL get_kpoint_info(kpoint, nkp=nkp)
1239 kp => kpoint%kp_env(1)%kpoint_env
1240 nspin = SIZE(kp%mos, 2)
1241 mo_set => kp%mos(1, 1)
1242 CALL get_mo_set(mo_set, nmo=nmo, nao=nao, nelectron=nelectron)
1243 ne_a = nelectron
1244 IF (nspin == 2) THEN
1245 CALL get_mo_set(kp%mos(1, 2), nmo=nb, nelectron=ne_b)
1246 cpassert(nmo == nb)
1247 END IF
1248 ALLOCATE (weig(nmo, nkp, nspin), wocc(nmo, nkp, nspin))
1249 weig = 0.0_dp
1250 wocc = 0.0_dp
1251 IF (PRESENT(probe)) THEN
1252 ALLOCATE (rcoeff(nao, nmo, nkp, nspin), icoeff(nao, nmo, nkp, nspin))
1253 rcoeff = 0.0_dp !coeff, real part
1254 icoeff = 0.0_dp !coeff, imaginary part
1255 END IF
1256 CALL get_kpoint_info(kpoint, kp_range=kp_range)
1257 kplocal = kp_range(2) - kp_range(1) + 1
1258 DO ikpgr = 1, kplocal
1259 ik = kp_range(1) + ikpgr - 1
1260 kp => kpoint%kp_env(ikpgr)%kpoint_env
1261 DO ispin = 1, nspin
1262 mo_set => kp%mos(1, ispin)
1263 CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
1264 weig(1:nmo, ik, ispin) = eigenvalues(1:nmo)
1265 IF (PRESENT(probe)) THEN
1266 CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
1267 CALL cp_fm_get_info(mo_coeff, &
1268 nrow_global=nrow_global, &
1269 ncol_global=ncol_global)
1270 ALLOCATE (smatrix(nrow_global, ncol_global))
1271 CALL cp_fm_get_submatrix(mo_coeff, smatrix)
1272
1273 rcoeff(1:nao, 1:nmo, ik, ispin) = smatrix(1:nrow_global, 1:ncol_global)
1274
1275 DEALLOCATE (smatrix)
1276
1277 mo_set => kp%mos(2, ispin)
1278
1279 CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
1280 CALL cp_fm_get_info(mo_coeff, &
1281 nrow_global=nrow_global, &
1282 ncol_global=ncol_global)
1283 ALLOCATE (smatrix(nrow_global, ncol_global))
1284 CALL cp_fm_get_submatrix(mo_coeff, smatrix)
1285
1286 icoeff(1:nao, 1:nmo, ik, ispin) = smatrix(1:nrow_global, 1:ncol_global)
1287
1288 mo_set => kp%mos(1, ispin)
1289
1290 DEALLOCATE (smatrix)
1291 END IF
1292 END DO
1293 END DO
1294 CALL get_kpoint_info(kpoint, para_env_inter_kp=para_env_inter_kp)
1295 CALL para_env_inter_kp%sum(weig)
1296
1297 IF (PRESENT(probe)) THEN
1298 CALL para_env_inter_kp%sum(rcoeff)
1299 CALL para_env_inter_kp%sum(icoeff)
1300 END IF
1301
1302 CALL get_kpoint_info(kpoint, wkp=wkp)
1303
1304!calling of HP module HERE, before smear
1305 IF (PRESENT(probe)) THEN
1306 smear%do_smear = .false. !ensures smearing is switched off
1307
1308 IF (nspin == 1) THEN
1309 nel = real(nelectron, kind=dp)
1310 CALL probe_occupancy_kp(wocc(:, :, :), mus(1), kts, weig(:, :, :), rcoeff(:, :, :, :), icoeff(:, :, :, :), 2.0d0, &
1311 probe, nel, wkp)
1312 ELSE
1313 nel = real(ne_a, kind=dp) + real(ne_b, kind=dp)
1314 CALL probe_occupancy_kp(wocc(:, :, :), mu, kts, weig(:, :, :), rcoeff(:, :, :, :), icoeff(:, :, :, :), 1.0d0, &
1315 probe, nel, wkp)
1316 kts = kts/2._dp
1317 mus(1:2) = mu
1318 END IF
1319
1320 DO ikpgr = 1, kplocal
1321 ik = kp_range(1) + ikpgr - 1
1322 kp => kpoint%kp_env(ikpgr)%kpoint_env
1323 DO ispin = 1, nspin
1324 mo_set => kp%mos(1, ispin)
1325 CALL get_mo_set(mo_set, eigenvalues=eigenvalues, occupation_numbers=occupation)
1326 eigenvalues(1:nmo) = weig(1:nmo, ik, ispin)
1327 occupation(1:nmo) = wocc(1:nmo, ik, ispin)
1328 mo_set%kTS = kts
1329 mo_set%mu = mus(ispin)
1330
1331 CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
1332 !get smatrix for kpoint_env ikp
1333 CALL cp_fm_get_info(mo_coeff, &
1334 nrow_global=nrow_global, &
1335 ncol_global=ncol_global)
1336 ALLOCATE (smatrix(nrow_global, ncol_global))
1337 CALL cp_fm_get_submatrix(mo_coeff, smatrix)
1338
1339 smatrix(1:nrow_global, 1:ncol_global) = rcoeff(1:nao, 1:nmo, ik, ispin)
1340 DEALLOCATE (smatrix)
1341
1342 mo_set => kp%mos(2, ispin)
1343
1344 CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
1345 !get smatrix for kpoint_env ikp
1346 CALL cp_fm_get_info(mo_coeff, &
1347 nrow_global=nrow_global, &
1348 ncol_global=ncol_global)
1349 ALLOCATE (smatrix(nrow_global, ncol_global))
1350 CALL cp_fm_get_submatrix(mo_coeff, smatrix)
1351
1352 smatrix(1:nrow_global, 1:ncol_global) = icoeff(1:nao, 1:nmo, ik, ispin)
1353 DEALLOCATE (smatrix)
1354
1355 mo_set => kp%mos(1, ispin)
1356
1357 END DO
1358 END DO
1359
1360 DEALLOCATE (weig, wocc, rcoeff, icoeff)
1361
1362 END IF
1363
1364 IF (PRESENT(probe) .EQV. .false.) THEN
1365 IF (smear%do_smear) THEN
1366 SELECT CASE (smear%method)
1367 CASE (smear_fermi_dirac)
1368 ! finite electronic temperature
1369 IF (nspin == 1) THEN
1370 nel = real(nelectron, kind=dp)
1371 CALL smearkp(wocc(:, :, 1), mus(1), kts, weig(:, :, 1), nel, wkp, &
1372 smear%electronic_temperature, 2.0_dp, smear_fermi_dirac)
1373 ELSE IF (smear%fixed_mag_mom > 0.0_dp) THEN
1374 nel = real(ne_a, kind=dp)
1375 CALL smearkp(wocc(:, :, 1), mus(1), kts, weig(:, :, 1), nel, wkp, &
1376 smear%electronic_temperature, 1.0_dp, smear_fermi_dirac)
1377 nel = real(ne_b, kind=dp)
1378 CALL smearkp(wocc(:, :, 2), mus(2), kts, weig(:, :, 2), nel, wkp, &
1379 smear%electronic_temperature, 1.0_dp, smear_fermi_dirac)
1380 ELSE
1381 nel = real(ne_a, kind=dp) + real(ne_b, kind=dp)
1382 CALL smearkp2(wocc(:, :, :), mu, kts, weig(:, :, :), nel, wkp, &
1383 smear%electronic_temperature, smear_fermi_dirac)
1384 kts = kts/2._dp
1385 mus(1:2) = mu
1386 END IF
1387 CASE (smear_gaussian, smear_mp, smear_mv)
1388 IF (nspin == 1) THEN
1389 nel = real(nelectron, kind=dp)
1390 CALL smearkp(wocc(:, :, 1), mus(1), kts, weig(:, :, 1), nel, wkp, &
1391 smear%smearing_width, 2.0_dp, smear%method)
1392 ELSE IF (smear%fixed_mag_mom > 0.0_dp) THEN
1393 nel = real(ne_a, kind=dp)
1394 CALL smearkp(wocc(:, :, 1), mus(1), kts, weig(:, :, 1), nel, wkp, &
1395 smear%smearing_width, 1.0_dp, smear%method)
1396 nel = real(ne_b, kind=dp)
1397 CALL smearkp(wocc(:, :, 2), mus(2), kts, weig(:, :, 2), nel, wkp, &
1398 smear%smearing_width, 1.0_dp, smear%method)
1399 ELSE
1400 nel = real(ne_a, kind=dp) + real(ne_b, kind=dp)
1401 CALL smearkp2(wocc(:, :, :), mu, kts, weig(:, :, :), nel, wkp, &
1402 smear%smearing_width, smear%method)
1403 kts = kts/2._dp
1404 mus(1:2) = mu
1405 END IF
1406 CASE DEFAULT
1407 cpabort("kpoints: Selected smearing not (yet) supported")
1408 END SELECT
1409 ELSE
1410 ! fixed occupations (2/1)
1411 IF (nspin == 1) THEN
1412 nel = real(nelectron, kind=dp)
1413 CALL smearkp(wocc(:, :, 1), mus(1), kts, weig(:, :, 1), nel, wkp, &
1414 0.0_dp, 2.0_dp, smear_gaussian)
1415 ELSE
1416 nel = real(ne_a, kind=dp)
1417 CALL smearkp(wocc(:, :, 1), mus(1), kts, weig(:, :, 1), nel, wkp, &
1418 0.0_dp, 1.0_dp, smear_gaussian)
1419 nel = real(ne_b, kind=dp)
1420 CALL smearkp(wocc(:, :, 2), mus(2), kts, weig(:, :, 2), nel, wkp, &
1421 0.0_dp, 1.0_dp, smear_gaussian)
1422 END IF
1423 END IF
1424 DO ikpgr = 1, kplocal
1425 ik = kp_range(1) + ikpgr - 1
1426 kp => kpoint%kp_env(ikpgr)%kpoint_env
1427 DO ispin = 1, nspin
1428 mo_set => kp%mos(1, ispin)
1429 CALL get_mo_set(mo_set, eigenvalues=eigenvalues, occupation_numbers=occupation)
1430 eigenvalues(1:nmo) = weig(1:nmo, ik, ispin)
1431 occupation(1:nmo) = wocc(1:nmo, ik, ispin)
1432 mo_set%kTS = kts
1433 mo_set%mu = mus(ispin)
1434 END DO
1435 END DO
1436
1437 DEALLOCATE (weig, wocc)
1438
1439 END IF
1440
1441 CALL timestop(handle)
1442
1443 END SUBROUTINE kpoint_set_mo_occupation
1444
1445! **************************************************************************************************
1446!> \brief Calculate kpoint density matrices (rho(k), owned by kpoint groups)
1447!> \param kpoint kpoint environment
1448!> \param energy_weighted calculate energy weighted density matrix
1449!> \param for_aux_fit ...
1450! **************************************************************************************************
1451 SUBROUTINE kpoint_density_matrices(kpoint, energy_weighted, for_aux_fit)
1452
1453 TYPE(kpoint_type), POINTER :: kpoint
1454 LOGICAL, OPTIONAL :: energy_weighted, for_aux_fit
1455
1456 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_density_matrices'
1457
1458 INTEGER :: handle, ikpgr, ispin, kplocal, nao, nmo, &
1459 nspin
1460 INTEGER, DIMENSION(2) :: kp_range
1461 LOGICAL :: aux_fit, wtype
1462 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues, occupation
1463 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1464 TYPE(cp_fm_type) :: fwork
1465 TYPE(cp_fm_type), POINTER :: cpmat, pmat, rpmat
1466 TYPE(kpoint_env_type), POINTER :: kp
1467 TYPE(mo_set_type), POINTER :: mo_set
1468
1469 CALL timeset(routinen, handle)
1470
1471 IF (PRESENT(energy_weighted)) THEN
1472 wtype = energy_weighted
1473 ELSE
1474 ! default is normal density matrix
1475 wtype = .false.
1476 END IF
1477
1478 IF (PRESENT(for_aux_fit)) THEN
1479 aux_fit = for_aux_fit
1480 ELSE
1481 aux_fit = .false.
1482 END IF
1483
1484 IF (aux_fit) THEN
1485 cpassert(ASSOCIATED(kpoint%kp_aux_env))
1486 END IF
1487
1488 ! work matrix
1489 IF (aux_fit) THEN
1490 mo_set => kpoint%kp_aux_env(1)%kpoint_env%mos(1, 1)
1491 ELSE
1492 mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)
1493 END IF
1494 CALL get_mo_set(mo_set, nao=nao, nmo=nmo)
1495 CALL cp_fm_get_info(mo_set%mo_coeff, matrix_struct=matrix_struct)
1496 CALL cp_fm_create(fwork, matrix_struct)
1497
1498 CALL get_kpoint_info(kpoint, kp_range=kp_range)
1499 kplocal = kp_range(2) - kp_range(1) + 1
1500 DO ikpgr = 1, kplocal
1501 IF (aux_fit) THEN
1502 kp => kpoint%kp_aux_env(ikpgr)%kpoint_env
1503 ELSE
1504 kp => kpoint%kp_env(ikpgr)%kpoint_env
1505 END IF
1506 nspin = SIZE(kp%mos, 2)
1507 DO ispin = 1, nspin
1508 mo_set => kp%mos(1, ispin)
1509 IF (wtype) THEN
1510 CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
1511 END IF
1512 IF (kpoint%use_real_wfn) THEN
1513 IF (wtype) THEN
1514 pmat => kp%wmat(1, ispin)
1515 ELSE
1516 pmat => kp%pmat(1, ispin)
1517 END IF
1518 CALL get_mo_set(mo_set, occupation_numbers=occupation)
1519 CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
1520 CALL cp_fm_column_scale(fwork, occupation)
1521 IF (wtype) THEN
1522 CALL cp_fm_column_scale(fwork, eigenvalues)
1523 END IF
1524 CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, pmat)
1525 ELSE
1526 IF (wtype) THEN
1527 rpmat => kp%wmat(1, ispin)
1528 cpmat => kp%wmat(2, ispin)
1529 ELSE
1530 rpmat => kp%pmat(1, ispin)
1531 cpmat => kp%pmat(2, ispin)
1532 END IF
1533 CALL get_mo_set(mo_set, occupation_numbers=occupation)
1534 CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
1535 CALL cp_fm_column_scale(fwork, occupation)
1536 IF (wtype) THEN
1537 CALL cp_fm_column_scale(fwork, eigenvalues)
1538 END IF
1539 ! Re(c)*Re(c)
1540 CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, rpmat)
1541 mo_set => kp%mos(2, ispin)
1542 ! Im(c)*Re(c)
1543 CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, cpmat)
1544 ! Re(c)*Im(c)
1545 CALL parallel_gemm("N", "T", nao, nao, nmo, -1.0_dp, fwork, mo_set%mo_coeff, 1.0_dp, cpmat)
1546 CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
1547 CALL cp_fm_column_scale(fwork, occupation)
1548 IF (wtype) THEN
1549 CALL cp_fm_column_scale(fwork, eigenvalues)
1550 END IF
1551 ! Im(c)*Im(c)
1552 CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 1.0_dp, rpmat)
1553 END IF
1554 END DO
1555 END DO
1556
1557 CALL cp_fm_release(fwork)
1558
1559 CALL timestop(handle)
1560
1561 END SUBROUTINE kpoint_density_matrices
1562
1563! **************************************************************************************************
1564!> \brief Calculate Lowdin transformation of density matrix S^1/2 P S^1/2
1565!> Integrate diagonal elements over k-points to get Lowdin charges
1566!> \param kpoint kpoint environment
1567!> \param pmat_diag Sum over kpoints of diagonal elements
1568!> \par History
1569!> 04.2026 created [JGH]
1570! **************************************************************************************************
1571 SUBROUTINE lowdin_kp_trans(kpoint, pmat_diag)
1572
1573 TYPE(kpoint_type), POINTER :: kpoint
1574 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: pmat_diag
1575
1576 CHARACTER(LEN=*), PARAMETER :: routinen = 'lowdin_kp_trans'
1577 COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp, 0.0_dp), &
1578 czero = (0.0_dp, 0.0_dp)
1579
1580 INTEGER :: handle, ikpgr, ispin, kplocal, nao, nspin
1581 INTEGER, DIMENSION(2) :: kp_range
1582 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dele
1583 TYPE(cp_cfm_type) :: cf1work, cf2work
1584 TYPE(cp_cfm_type), POINTER :: cshalf
1585 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1586 TYPE(cp_fm_type) :: f1work, f2work
1587 TYPE(cp_fm_type), POINTER :: cpmat, pmat, rpmat, shalf
1588 TYPE(kpoint_env_type), POINTER :: kp
1589 TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
1590
1591 CALL timeset(routinen, handle)
1592
1593 nspin = SIZE(pmat_diag, 2)
1594 pmat_diag = 0.0_dp
1595
1596 ! work matrix
1597 CALL cp_fm_get_info(kpoint%kp_env(1)%kpoint_env%pmat(1, 1), &
1598 matrix_struct=matrix_struct, nrow_global=nao)
1599 IF (kpoint%use_real_wfn) THEN
1600 CALL cp_fm_create(f1work, matrix_struct, nrow=nao, ncol=nao)
1601 CALL cp_fm_create(f2work, matrix_struct, nrow=nao, ncol=nao)
1602 ELSE
1603 CALL cp_fm_create(f2work, matrix_struct, nrow=nao, ncol=nao)
1604 CALL cp_cfm_create(cf1work, matrix_struct, nrow=nao, ncol=nao)
1605 CALL cp_cfm_create(cf2work, matrix_struct, nrow=nao, ncol=nao)
1606 END IF
1607 ALLOCATE (dele(nao))
1608
1609 CALL get_kpoint_info(kpoint, kp_range=kp_range)
1610 kplocal = kp_range(2) - kp_range(1) + 1
1611 DO ikpgr = 1, kplocal
1612 kp => kpoint%kp_env(ikpgr)%kpoint_env
1613 DO ispin = 1, nspin
1614 IF (kpoint%use_real_wfn) THEN
1615 pmat => kp%pmat(1, ispin)
1616 shalf => kp%shalf
1617 CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, pmat, shalf, 0.0_dp, f1work)
1618 CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, shalf, f1work, 0.0_dp, f2work)
1619 ELSE
1620 rpmat => kp%pmat(1, ispin)
1621 cpmat => kp%pmat(2, ispin)
1622 cshalf => kp%cshalf
1623 CALL cp_fm_to_cfm(rpmat, cpmat, cf1work)
1624 CALL parallel_gemm("N", "N", nao, nao, nao, cone, cf1work, cshalf, czero, cf2work)
1625 CALL parallel_gemm("N", "N", nao, nao, nao, cone, cshalf, cf2work, czero, cf1work)
1626 CALL cp_cfm_to_fm(cf1work, mtargetr=f2work)
1627 END IF
1628 CALL cp_fm_get_diag(f2work, dele)
1629 pmat_diag(1:nao, ispin) = pmat_diag(1:nao, ispin) + kp%wkp*dele(1:nao)
1630 END DO
1631 END DO
1632
1633 CALL get_kpoint_info(kpoint, para_env_inter_kp=para_env_inter_kp)
1634 CALL para_env_inter_kp%sum(pmat_diag)
1635
1636 IF (kpoint%use_real_wfn) THEN
1637 CALL cp_fm_release(f1work)
1638 CALL cp_fm_release(f2work)
1639 ELSE
1640 CALL cp_fm_release(f2work)
1641 CALL cp_cfm_release(cf1work)
1642 CALL cp_cfm_release(cf2work)
1643 END IF
1644 DEALLOCATE (dele)
1645
1646 CALL timestop(handle)
1647
1648 END SUBROUTINE lowdin_kp_trans
1649
1650! **************************************************************************************************
1651!> \brief generate real space density matrices in DBCSR format
1652!> \param kpoint Kpoint environment
1653!> \param denmat Real space (DBCSR) density matrices
1654!> \param wtype True = energy weighted density matrix
1655!> False = normal density matrix
1656!> \param tempmat DBCSR matrix to be used as template
1657!> \param sab_nl ...
1658!> \param fmwork FM work matrices (kpoint group)
1659!> \param for_aux_fit ...
1660!> \param pmat_ext ...
1661! **************************************************************************************************
1662 SUBROUTINE kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwork, for_aux_fit, pmat_ext)
1663
1664 TYPE(kpoint_type), POINTER :: kpoint
1665 TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
1666 LOGICAL, INTENT(IN) :: wtype
1667 TYPE(dbcsr_type), POINTER :: tempmat
1668 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1669 POINTER :: sab_nl
1670 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fmwork
1671 LOGICAL, OPTIONAL :: for_aux_fit
1672 TYPE(cp_fm_type), DIMENSION(:, :, :), INTENT(IN), &
1673 OPTIONAL :: pmat_ext
1674
1675 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_density_transform'
1676
1677 INTEGER :: handle, ic, ik, ikk, indx, ir, ira, is, &
1678 ispin, jr, nc, nimg, nkp, nspin
1679 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1680 LOGICAL :: aux_fit, do_ext, do_symmetric, my_kpgrp, &
1681 real_only, reverse_phase
1682 REAL(kind=dp) :: wkpx
1683 REAL(kind=dp), DIMENSION(:), POINTER :: wkp
1684 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
1685 TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:) :: info
1686 TYPE(cp_fm_type) :: fmdummy
1687 TYPE(dbcsr_type), POINTER :: cpmat, rpmat, scpmat, srpmat
1688 TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kind_rot
1689 TYPE(kpoint_env_type), POINTER :: kp
1690 TYPE(kpoint_sym_type), POINTER :: kpsym
1691 TYPE(mp_para_env_type), POINTER :: para_env
1692
1693 CALL timeset(routinen, handle)
1694
1695 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1696
1697 IF (PRESENT(for_aux_fit)) THEN
1698 aux_fit = for_aux_fit
1699 ELSE
1700 aux_fit = .false.
1701 END IF
1702
1703 do_ext = .false.
1704 IF (PRESENT(pmat_ext)) do_ext = .true.
1705
1706 IF (aux_fit) THEN
1707 cpassert(ASSOCIATED(kpoint%kp_aux_env))
1708 END IF
1709
1710 ! work storage
1711 ALLOCATE (rpmat)
1712 CALL dbcsr_create(rpmat, template=tempmat, &
1713 matrix_type=merge(dbcsr_type_symmetric, dbcsr_type_no_symmetry, do_symmetric))
1714 CALL cp_dbcsr_alloc_block_from_nbl(rpmat, sab_nl)
1715 CALL dbcsr_set(rpmat, 0.0_dp)
1716 ALLOCATE (cpmat)
1717 CALL dbcsr_create(cpmat, template=tempmat, &
1718 matrix_type=merge(dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, do_symmetric))
1719 CALL cp_dbcsr_alloc_block_from_nbl(cpmat, sab_nl)
1720 CALL dbcsr_set(cpmat, 0.0_dp)
1721 IF (.NOT. kpoint%full_grid) THEN
1722 ALLOCATE (srpmat)
1723 CALL dbcsr_create(srpmat, template=rpmat)
1724 CALL cp_dbcsr_alloc_block_from_nbl(srpmat, sab_nl)
1725 CALL dbcsr_set(srpmat, 0.0_dp)
1726 ALLOCATE (scpmat)
1727 CALL dbcsr_create(scpmat, template=cpmat)
1728 CALL cp_dbcsr_alloc_block_from_nbl(scpmat, sab_nl)
1729 CALL dbcsr_set(scpmat, 0.0_dp)
1730 END IF
1731
1732 CALL get_kpoint_info(kpoint, nkp=nkp, xkp=xkp, wkp=wkp, &
1733 cell_to_index=cell_to_index)
1734 ! initialize real space density matrices
1735 IF (aux_fit) THEN
1736 kp => kpoint%kp_aux_env(1)%kpoint_env
1737 ELSE
1738 kp => kpoint%kp_env(1)%kpoint_env
1739 END IF
1740 nspin = SIZE(kp%mos, 2)
1741 nc = SIZE(kp%mos, 1)
1742 nimg = SIZE(denmat, 2)
1743 real_only = (nc == 1)
1744 ! Gamma-centered even grids store atom-cell shifts in the opposite Bloch-phase convention.
1745 reverse_phase = kpoint%gamma_centered .AND. any(mod(kpoint%nkp_grid(1:3), 2) == 0)
1746
1747 para_env => kpoint%blacs_env_all%para_env
1748 ALLOCATE (info(nspin*nkp*nc))
1749
1750 ! Start all the communication
1751 indx = 0
1752 DO ispin = 1, nspin
1753 DO ic = 1, nimg
1754 CALL dbcsr_set(denmat(ispin, ic)%matrix, 0.0_dp)
1755 END DO
1756 !
1757 DO ik = 1, nkp
1758 my_kpgrp = (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2))
1759 IF (my_kpgrp) THEN
1760 ikk = ik - kpoint%kp_range(1) + 1
1761 IF (aux_fit) THEN
1762 kp => kpoint%kp_aux_env(ikk)%kpoint_env
1763 ELSE
1764 kp => kpoint%kp_env(ikk)%kpoint_env
1765 END IF
1766 ELSE
1767 NULLIFY (kp)
1768 END IF
1769 ! collect this density matrix on all processors
1770 cpassert(SIZE(fmwork) >= nc)
1771
1772 IF (my_kpgrp) THEN
1773 DO ic = 1, nc
1774 indx = indx + 1
1775 IF (do_ext) THEN
1776 CALL cp_fm_start_copy_general(pmat_ext(ikk, ic, ispin), fmwork(ic), para_env, info(indx))
1777 ELSE
1778 IF (wtype) THEN
1779 CALL cp_fm_start_copy_general(kp%wmat(ic, ispin), fmwork(ic), para_env, info(indx))
1780 ELSE
1781 CALL cp_fm_start_copy_general(kp%pmat(ic, ispin), fmwork(ic), para_env, info(indx))
1782 END IF
1783 END IF
1784 END DO
1785 ELSE
1786 DO ic = 1, nc
1787 indx = indx + 1
1788 CALL cp_fm_start_copy_general(fmdummy, fmwork(ic), para_env, info(indx))
1789 END DO
1790 END IF
1791 END DO
1792 END DO
1793
1794 ! Finish communication and transform the received matrices
1795 indx = 0
1796 DO ispin = 1, nspin
1797 DO ik = 1, nkp
1798 DO ic = 1, nc
1799 indx = indx + 1
1800 CALL cp_fm_finish_copy_general(fmwork(ic), info(indx))
1801 END DO
1802
1803 ! reduce to dbcsr storage
1804 IF (real_only) THEN
1805 CALL copy_fm_to_dbcsr(fmwork(1), rpmat, keep_sparsity=.true.)
1806 ELSE
1807 CALL copy_fm_to_dbcsr(fmwork(1), rpmat, keep_sparsity=.true.)
1808 CALL copy_fm_to_dbcsr(fmwork(2), cpmat, keep_sparsity=.true.)
1809 END IF
1810
1811 ! symmetrization
1812 kpsym => kpoint%kp_sym(ik)%kpoint_sym
1813 cpassert(ASSOCIATED(kpsym))
1814
1815 IF (kpsym%apply_symmetry) THEN
1816 wkpx = wkp(ik)/real(kpsym%nwght, kind=dp)
1817 DO is = 1, kpsym%nwght
1818 ir = abs(kpsym%rotp(is))
1819 ira = 0
1820 DO jr = 1, SIZE(kpoint%ibrot)
1821 IF (ir == kpoint%ibrot(jr)) ira = jr
1822 END DO
1823 cpassert(ira > 0)
1824 kind_rot => kpoint%kind_rotmat(ira, :)
1825 CALL symtrans_phase(srpmat, scpmat, rpmat, cpmat, real_only, kind_rot, &
1826 kpsym%rot(1:3, 1:3, is), kpsym%f0(:, is), &
1827 kpsym%fcell(:, :, is), kpoint%atype, kpsym%xkp(1:3, is), &
1828 kpsym%rotp(is) < 0, reverse_phase)
1829 CALL transform_dmat(denmat, srpmat, scpmat, ispin, real_only, sab_nl, &
1830 cell_to_index, kpsym%xkp(1:3, is), wkpx)
1831 END DO
1832 ELSE
1833 ! transformation
1834 CALL transform_dmat(denmat, rpmat, cpmat, ispin, real_only, sab_nl, &
1835 cell_to_index, xkp(1:3, ik), wkp(ik))
1836 END IF
1837 END DO
1838 END DO
1839
1840 ! Clean up communication
1841 indx = 0
1842 DO ispin = 1, nspin
1843 DO ik = 1, nkp
1844 my_kpgrp = (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2))
1845 IF (my_kpgrp) THEN
1846 ikk = ik - kpoint%kp_range(1) + 1
1847 IF (aux_fit) THEN
1848 kp => kpoint%kp_aux_env(ikk)%kpoint_env
1849 ELSE
1850 kp => kpoint%kp_env(ikk)%kpoint_env
1851 END IF
1852
1853 DO ic = 1, nc
1854 indx = indx + 1
1855 CALL cp_fm_cleanup_copy_general(info(indx))
1856 END DO
1857 ELSE
1858 ! calls with dummy arguments, so not included
1859 ! therefore just increment counter by trip count
1860 indx = indx + nc
1861 END IF
1862 END DO
1863 END DO
1864
1865 ! All done
1866 DEALLOCATE (info)
1867
1868 CALL dbcsr_deallocate_matrix(rpmat)
1869 CALL dbcsr_deallocate_matrix(cpmat)
1870 IF (.NOT. kpoint%full_grid) THEN
1871 CALL dbcsr_deallocate_matrix(srpmat)
1872 CALL dbcsr_deallocate_matrix(scpmat)
1873 END IF
1874
1875 CALL timestop(handle)
1876
1877 END SUBROUTINE kpoint_density_transform
1878
1879! **************************************************************************************************
1880!> \brief real space density matrices in DBCSR format
1881!> \param denmat Real space (DBCSR) density matrix
1882!> \param rpmat ...
1883!> \param cpmat ...
1884!> \param ispin ...
1885!> \param real_only ...
1886!> \param sab_nl ...
1887!> \param cell_to_index ...
1888!> \param xkp ...
1889!> \param wkp ...
1890! **************************************************************************************************
1891 SUBROUTINE transform_dmat(denmat, rpmat, cpmat, ispin, real_only, sab_nl, cell_to_index, xkp, wkp)
1892
1893 TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
1894 TYPE(dbcsr_type), POINTER :: rpmat, cpmat
1895 INTEGER, INTENT(IN) :: ispin
1896 LOGICAL, INTENT(IN) :: real_only
1897 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1898 POINTER :: sab_nl
1899 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1900 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: xkp
1901 REAL(kind=dp), INTENT(IN) :: wkp
1902
1903 CHARACTER(LEN=*), PARAMETER :: routinen = 'transform_dmat'
1904
1905 INTEGER :: handle, iatom, icell, icol, irow, jatom, &
1906 nimg
1907 INTEGER, DIMENSION(3) :: cell
1908 LOGICAL :: do_symmetric, found
1909 REAL(kind=dp) :: arg, coskl, fc, sinkl
1910 REAL(kind=dp), DIMENSION(:, :), POINTER :: cblock, dblock, rblock
1911 TYPE(neighbor_list_iterator_p_type), &
1912 DIMENSION(:), POINTER :: nl_iterator
1913
1914 CALL timeset(routinen, handle)
1915
1916 nimg = SIZE(denmat, 2)
1917
1918 ! transformation
1919 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1920 CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
1921 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1922 CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
1923
1924 !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
1925 !Therefore, we have: S(R) = sum_k Re(S(k))*cos(k*R) -i^2*Im(S(k))*sin(k*R)
1926 ! = sum_k Re(S(k))*cos(k*R) + Im(S(k))*sin(k*R)
1927 !fc = +- 1 is due to the usual non-symmetric real-space matrices stored as symmetric ones
1928
1929 irow = iatom
1930 icol = jatom
1931 fc = 1.0_dp
1932 IF (do_symmetric .AND. iatom > jatom) THEN
1933 irow = jatom
1934 icol = iatom
1935 fc = -1.0_dp
1936 END IF
1937
1938 icell = cell_to_index(cell(1), cell(2), cell(3))
1939 IF (icell < 1 .OR. icell > nimg) cycle
1940
1941 arg = real(cell(1), dp)*xkp(1) + real(cell(2), dp)*xkp(2) + real(cell(3), dp)*xkp(3)
1942 coskl = wkp*cos(twopi*arg)
1943 sinkl = wkp*fc*sin(twopi*arg)
1944
1945 CALL dbcsr_get_block_p(matrix=denmat(ispin, icell)%matrix, row=irow, col=icol, &
1946 block=dblock, found=found)
1947 IF (.NOT. found) cycle
1948
1949 IF (real_only) THEN
1950 CALL dbcsr_get_block_p(matrix=rpmat, row=irow, col=icol, block=rblock, found=found)
1951 IF (.NOT. found) cycle
1952 dblock = dblock + coskl*rblock
1953 ELSE
1954 CALL dbcsr_get_block_p(matrix=rpmat, row=irow, col=icol, block=rblock, found=found)
1955 IF (.NOT. found) cycle
1956 CALL dbcsr_get_block_p(matrix=cpmat, row=irow, col=icol, block=cblock, found=found)
1957 IF (.NOT. found) cycle
1958 dblock = dblock + coskl*rblock
1959 dblock = dblock + sinkl*cblock
1960 END IF
1961 END DO
1962 CALL neighbor_list_iterator_release(nl_iterator)
1963
1964 CALL timestop(handle)
1965
1966 END SUBROUTINE transform_dmat
1967
1968! **************************************************************************************************
1969!> \brief Allocate a dense work matrix with the requested shape
1970!> \param work dense work matrix
1971!> \param nrow number of rows
1972!> \param ncol number of columns
1973! **************************************************************************************************
1974 SUBROUTINE ensure_work_matrix(work, nrow, ncol)
1975
1976 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
1977 INTENT(INOUT) :: work
1978 INTEGER, INTENT(IN) :: nrow, ncol
1979
1980 IF (ALLOCATED(work)) THEN
1981 IF (SIZE(work, 1) == nrow .AND. SIZE(work, 2) == ncol) RETURN
1982 DEALLOCATE (work)
1983 END IF
1984 ALLOCATE (work(nrow, ncol))
1985
1986 END SUBROUTINE ensure_work_matrix
1987
1988! **************************************************************************************************
1989!> \brief Symmetrize a complex k-point matrix including Bloch phase shifts
1990!> \param srpmat real part of transformed matrix
1991!> \param scpmat imaginary part of transformed matrix
1992!> \param rpmat real part of reference matrix
1993!> \param cpmat imaginary part of reference matrix
1994!> \param real_only ...
1995!> \param kmat kind type rotation matrix
1996!> \param rot rotation matrix
1997!> \param f0 atom permutation
1998!> \param fcell atom cell shifts generated by the symmetry operation
1999!> \param atype atom to kind pointer
2000!> \param xkp target k-point coordinates
2001!> \param time_reversal ...
2002!> \param reverse_phase ...
2003! **************************************************************************************************
2004 SUBROUTINE symtrans_phase(srpmat, scpmat, rpmat, cpmat, real_only, kmat, rot, f0, fcell, atype, &
2005 xkp, time_reversal, reverse_phase)
2006
2007 TYPE(dbcsr_type), POINTER :: srpmat, scpmat, rpmat, cpmat
2008 LOGICAL, INTENT(IN) :: real_only
2009 TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kmat
2010 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: rot
2011 INTEGER, DIMENSION(:), INTENT(IN) :: f0
2012 INTEGER, DIMENSION(:, :), INTENT(IN) :: fcell
2013 INTEGER, DIMENSION(:), INTENT(IN) :: atype
2014 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: xkp
2015 LOGICAL, INTENT(IN) :: time_reversal, reverse_phase
2016
2017 CHARACTER(LEN=*), PARAMETER :: routinen = 'symtrans_phase'
2018
2019 INTEGER :: handle, iatom, icol, ikind, ip, irow, &
2020 jcol, jkind, jp, jrow, mynode, natom, &
2021 numnodes, owner
2022 INTEGER, DIMENSION(3) :: shift
2023 LOGICAL :: dorot, found, has_phase, perm, trans
2024 REAL(kind=dp) :: arg, coskl, dr, sinkl
2025 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cwork, rwork, twork
2026 REAL(kind=dp), DIMENSION(:, :), POINTER :: cblock, kroti, krotj, rblock, scblock, &
2027 srblock
2028 TYPE(dbcsr_distribution_type) :: dist
2029 TYPE(dbcsr_iterator_type) :: iter
2030
2031 CALL timeset(routinen, handle)
2032
2033 natom = SIZE(f0)
2034 perm = .false.
2035 DO iatom = 1, natom
2036 IF (f0(iatom) == iatom) cycle
2037 perm = .true.
2038 EXIT
2039 END DO
2040
2041 dorot = .false.
2042 IF (abs(sum(abs(rot)) - 3.0_dp) > 1.e-12_dp) dorot = .true.
2043 dr = abs(rot(1, 1) - 1.0_dp) + abs(rot(2, 2) - 1.0_dp) + abs(rot(3, 3) - 1.0_dp)
2044 IF (abs(dr) > 1.e-12_dp) dorot = .true.
2045 has_phase = any(fcell /= 0) .OR. time_reversal
2046
2047 IF (.NOT. (dorot .OR. perm .OR. has_phase)) THEN
2048 CALL dbcsr_copy(srpmat, rpmat)
2049 IF (.NOT. real_only) CALL dbcsr_copy(scpmat, cpmat)
2050 CALL timestop(handle)
2051 RETURN
2052 END IF
2053
2054 CALL dbcsr_get_info(rpmat, distribution=dist)
2055 CALL dbcsr_distribution_get(dist, mynode=mynode, numnodes=numnodes)
2056 IF (numnodes /= 1 .AND. (perm .OR. has_phase)) THEN
2057 CALL dbcsr_replicate_all(rpmat)
2058 IF (.NOT. real_only) CALL dbcsr_replicate_all(cpmat)
2059 END IF
2060
2061 CALL dbcsr_set(srpmat, 0.0_dp)
2062 IF (.NOT. real_only) CALL dbcsr_set(scpmat, 0.0_dp)
2063
2064 CALL dbcsr_iterator_start(iter, rpmat)
2065 DO WHILE (dbcsr_iterator_blocks_left(iter))
2066 CALL dbcsr_iterator_next_block(iter, irow, icol, rblock)
2067 IF (.NOT. ALLOCATED(rwork)) THEN
2068 ALLOCATE (rwork(SIZE(rblock, 1), SIZE(rblock, 2)))
2069 ELSEIF (SIZE(rwork, 1) /= SIZE(rblock, 1) .OR. SIZE(rwork, 2) /= SIZE(rblock, 2)) THEN
2070 DEALLOCATE (rwork)
2071 ALLOCATE (rwork(SIZE(rblock, 1), SIZE(rblock, 2)))
2072 END IF
2073 IF (.NOT. real_only) THEN
2074 IF (.NOT. ALLOCATED(cwork)) THEN
2075 ALLOCATE (cwork(SIZE(rblock, 1), SIZE(rblock, 2)))
2076 ELSEIF (SIZE(cwork, 1) /= SIZE(rblock, 1) .OR. SIZE(cwork, 2) /= SIZE(rblock, 2)) THEN
2077 DEALLOCATE (cwork)
2078 ALLOCATE (cwork(SIZE(rblock, 1), SIZE(rblock, 2)))
2079 END IF
2080 END IF
2081
2082 ikind = atype(irow)
2083 jkind = atype(icol)
2084 kroti => kmat(ikind)%rmat
2085 krotj => kmat(jkind)%rmat
2086
2087 IF (reverse_phase) THEN
2088 shift = fcell(1:3, irow) - fcell(1:3, icol)
2089 ELSE
2090 shift = fcell(1:3, icol) - fcell(1:3, irow)
2091 END IF
2092 arg = real(shift(1), dp)*xkp(1) + real(shift(2), dp)*xkp(2) + real(shift(3), dp)*xkp(3)
2093 ! The Bloch phase depends only on the mapped atom-cell shifts and the target k-point.
2094 coskl = cos(twopi*arg)
2095 sinkl = sin(twopi*arg)
2096 IF (real_only) THEN
2097 IF (abs(sinkl) > 1.e-12_dp) THEN
2098 CALL cp_abort(__location__, "Real k-point wavefunctions cannot represent symmetry phases")
2099 END IF
2100 rwork(:, :) = coskl*rblock
2101 ELSE
2102 CALL dbcsr_get_block_p(matrix=cpmat, row=irow, col=icol, block=cblock, found=found)
2103 rwork(:, :) = coskl*rblock
2104 IF (time_reversal) THEN
2105 cwork(:, :) = -sinkl*rblock
2106 IF (found) THEN
2107 rwork(:, :) = rwork - sinkl*cblock
2108 cwork(:, :) = cwork - coskl*cblock
2109 END IF
2110 ELSE
2111 cwork(:, :) = -sinkl*rblock
2112 IF (found) THEN
2113 rwork(:, :) = rwork + sinkl*cblock
2114 cwork(:, :) = cwork + coskl*cblock
2115 END IF
2116 END IF
2117 END IF
2118
2119 ip = f0(irow)
2120 jp = f0(icol)
2121 IF (ip <= jp) THEN
2122 jrow = ip
2123 jcol = jp
2124 trans = .false.
2125 ELSE
2126 jrow = jp
2127 jcol = ip
2128 trans = .true.
2129 END IF
2130
2131 CALL dbcsr_get_block_p(matrix=srpmat, row=jrow, col=jcol, block=srblock, found=found)
2132 IF (.NOT. found) THEN
2133 CALL dbcsr_get_stored_coordinates(srpmat, jrow, jcol, owner)
2134 cpassert(owner /= mynode)
2135 cycle
2136 END IF
2137 IF (trans) THEN
2138 CALL ensure_work_matrix(twork, SIZE(krotj, 1), SIZE(rwork, 1))
2139 CALL dgemm('N', 'T', SIZE(krotj, 1), SIZE(rwork, 1), SIZE(krotj, 2), &
2140 1.0_dp, krotj, SIZE(krotj, 1), rwork, SIZE(rwork, 1), &
2141 0.0_dp, twork, SIZE(twork, 1))
2142 CALL dgemm('N', 'T', SIZE(twork, 1), SIZE(kroti, 1), SIZE(twork, 2), &
2143 1.0_dp, twork, SIZE(twork, 1), kroti, SIZE(kroti, 1), &
2144 1.0_dp, srblock, SIZE(srblock, 1))
2145 ELSE
2146 CALL ensure_work_matrix(twork, SIZE(kroti, 1), SIZE(rwork, 2))
2147 CALL dgemm('N', 'N', SIZE(kroti, 1), SIZE(rwork, 2), SIZE(kroti, 2), &
2148 1.0_dp, kroti, SIZE(kroti, 1), rwork, SIZE(rwork, 1), &
2149 0.0_dp, twork, SIZE(twork, 1))
2150 CALL dgemm('N', 'T', SIZE(twork, 1), SIZE(krotj, 1), SIZE(twork, 2), &
2151 1.0_dp, twork, SIZE(twork, 1), krotj, SIZE(krotj, 1), &
2152 1.0_dp, srblock, SIZE(srblock, 1))
2153 END IF
2154
2155 IF (.NOT. real_only) THEN
2156 CALL dbcsr_get_block_p(matrix=scpmat, row=jrow, col=jcol, block=scblock, found=found)
2157 cpassert(found)
2158 IF (trans) THEN
2159 CALL ensure_work_matrix(twork, SIZE(krotj, 1), SIZE(cwork, 1))
2160 CALL dgemm('N', 'T', SIZE(krotj, 1), SIZE(cwork, 1), SIZE(krotj, 2), &
2161 1.0_dp, krotj, SIZE(krotj, 1), cwork, SIZE(cwork, 1), &
2162 0.0_dp, twork, SIZE(twork, 1))
2163 CALL dgemm('N', 'T', SIZE(twork, 1), SIZE(kroti, 1), SIZE(twork, 2), &
2164 -1.0_dp, twork, SIZE(twork, 1), kroti, SIZE(kroti, 1), &
2165 1.0_dp, scblock, SIZE(scblock, 1))
2166 ELSE
2167 CALL ensure_work_matrix(twork, SIZE(kroti, 1), SIZE(cwork, 2))
2168 CALL dgemm('N', 'N', SIZE(kroti, 1), SIZE(cwork, 2), SIZE(kroti, 2), &
2169 1.0_dp, kroti, SIZE(kroti, 1), cwork, SIZE(cwork, 1), &
2170 0.0_dp, twork, SIZE(twork, 1))
2171 CALL dgemm('N', 'T', SIZE(twork, 1), SIZE(krotj, 1), SIZE(twork, 2), &
2172 1.0_dp, twork, SIZE(twork, 1), krotj, SIZE(krotj, 1), &
2173 1.0_dp, scblock, SIZE(scblock, 1))
2174 END IF
2175 END IF
2176 END DO
2177 CALL dbcsr_iterator_stop(iter)
2178 IF (numnodes /= 1 .AND. (perm .OR. has_phase)) THEN
2179 CALL dbcsr_distribute(rpmat)
2180 IF (.NOT. real_only) CALL dbcsr_distribute(cpmat)
2181 END IF
2182
2183 CALL timestop(handle)
2184
2185 END SUBROUTINE symtrans_phase
2186
2187! **************************************************************************************************
2188!> \brief Symmetrization of density matrix - transform to new k-point
2189!> \param smat density matrix at new kpoint
2190!> \param pmat reference density matrix
2191!> \param kmat Kind type rotation matrix
2192!> \param rot Rotation matrix
2193!> \param f0 Permutation of atoms under transformation
2194!> \param atype Atom to kind pointer
2195!> \param symmetric Symmetric matrix
2196!> \param antisymmetric Anti-Symmetric matrix
2197! **************************************************************************************************
2198 SUBROUTINE symtrans(smat, pmat, kmat, rot, f0, atype, symmetric, antisymmetric)
2199 TYPE(dbcsr_type), POINTER :: smat, pmat
2200 TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kmat
2201 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: rot
2202 INTEGER, DIMENSION(:), INTENT(IN) :: f0, atype
2203 LOGICAL, INTENT(IN), OPTIONAL :: symmetric, antisymmetric
2204
2205 CHARACTER(LEN=*), PARAMETER :: routinen = 'symtrans'
2206
2207 INTEGER :: handle, iatom, icol, ikind, ip, irow, &
2208 jcol, jkind, jp, jrow, natom, numnodes
2209 LOGICAL :: asym, dorot, found, perm, sym, trans
2210 REAL(kind=dp) :: dr, fsign
2211 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
2212 REAL(kind=dp), DIMENSION(:, :), POINTER :: kroti, krotj, pblock, sblock
2213 TYPE(dbcsr_distribution_type) :: dist
2214 TYPE(dbcsr_iterator_type) :: iter
2215
2216 CALL timeset(routinen, handle)
2217
2218 ! check symmetry options
2219 sym = .false.
2220 IF (PRESENT(symmetric)) sym = symmetric
2221 asym = .false.
2222 IF (PRESENT(antisymmetric)) asym = antisymmetric
2223
2224 cpassert(.NOT. (sym .AND. asym))
2225 cpassert((sym .OR. asym))
2226
2227 ! do we have permutation of atoms
2228 natom = SIZE(f0)
2229 perm = .false.
2230 DO iatom = 1, natom
2231 IF (f0(iatom) == iatom) cycle
2232 perm = .true.
2233 EXIT
2234 END DO
2235
2236 ! do we have a real rotation
2237 dorot = .false.
2238 IF (abs(sum(abs(rot)) - 3.0_dp) > 1.e-12_dp) dorot = .true.
2239 dr = abs(rot(1, 1) - 1.0_dp) + abs(rot(2, 2) - 1.0_dp) + abs(rot(3, 3) - 1.0_dp)
2240 IF (abs(dr) > 1.e-12_dp) dorot = .true.
2241
2242 fsign = 1.0_dp
2243 IF (asym) fsign = -1.0_dp
2244
2245 IF (dorot .OR. perm) THEN
2246 CALL cp_abort(__location__, "k-points need FULL_GRID currently. "// &
2247 "Reduced grids not yet working correctly")
2248 CALL dbcsr_set(smat, 0.0_dp)
2249 IF (perm) THEN
2250 CALL dbcsr_get_info(pmat, distribution=dist)
2251 CALL dbcsr_distribution_get(dist, numnodes=numnodes)
2252 IF (numnodes == 1) THEN
2253 ! the matrices are local to this process
2254 CALL dbcsr_iterator_start(iter, pmat)
2255 DO WHILE (dbcsr_iterator_blocks_left(iter))
2256 CALL dbcsr_iterator_next_block(iter, irow, icol, pblock)
2257 ip = f0(irow)
2258 jp = f0(icol)
2259 IF (ip <= jp) THEN
2260 jrow = ip
2261 jcol = jp
2262 trans = .false.
2263 ELSE
2264 jrow = jp
2265 jcol = ip
2266 trans = .true.
2267 END IF
2268 CALL dbcsr_get_block_p(matrix=smat, row=jrow, col=jcol, block=sblock, found=found)
2269 cpassert(found)
2270 ikind = atype(irow)
2271 jkind = atype(icol)
2272 kroti => kmat(ikind)%rmat
2273 krotj => kmat(jkind)%rmat
2274 ! rotation
2275 IF (trans) THEN
2276 CALL ensure_work_matrix(work, SIZE(krotj, 2), SIZE(pblock, 1))
2277 CALL dgemm('T', 'T', SIZE(krotj, 2), SIZE(pblock, 1), SIZE(krotj, 1), &
2278 1.0_dp, krotj, SIZE(krotj, 1), pblock, SIZE(pblock, 1), &
2279 0.0_dp, work, SIZE(work, 1))
2280 CALL dgemm('N', 'N', SIZE(work, 1), SIZE(kroti, 2), SIZE(work, 2), &
2281 fsign, work, SIZE(work, 1), kroti, SIZE(kroti, 1), &
2282 0.0_dp, sblock, SIZE(sblock, 1))
2283 ELSE
2284 CALL ensure_work_matrix(work, SIZE(kroti, 2), SIZE(pblock, 2))
2285 CALL dgemm('T', 'N', SIZE(kroti, 2), SIZE(pblock, 2), SIZE(kroti, 1), &
2286 1.0_dp, kroti, SIZE(kroti, 1), pblock, SIZE(pblock, 1), &
2287 0.0_dp, work, SIZE(work, 1))
2288 CALL dgemm('N', 'N', SIZE(work, 1), SIZE(krotj, 2), SIZE(work, 2), &
2289 fsign, work, SIZE(work, 1), krotj, SIZE(krotj, 1), &
2290 0.0_dp, sblock, SIZE(sblock, 1))
2291 END IF
2292 END DO
2293 CALL dbcsr_iterator_stop(iter)
2294 !
2295 ELSE
2296 ! distributed matrices, most general code needed
2297 CALL cp_abort(__location__, "k-points need FULL_GRID currently. "// &
2298 "Reduced grids not yet working correctly")
2299 END IF
2300 ELSE
2301 ! no atom permutations, this is always local
2302 CALL dbcsr_copy(smat, pmat)
2303 CALL dbcsr_iterator_start(iter, smat)
2304 DO WHILE (dbcsr_iterator_blocks_left(iter))
2305 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
2306 ip = f0(irow)
2307 jp = f0(icol)
2308 IF (ip <= jp) THEN
2309 jrow = ip
2310 jcol = jp
2311 trans = .false.
2312 ELSE
2313 jrow = jp
2314 jcol = ip
2315 trans = .true.
2316 END IF
2317 ikind = atype(irow)
2318 jkind = atype(icol)
2319 kroti => kmat(ikind)%rmat
2320 krotj => kmat(jkind)%rmat
2321 ! rotation
2322 IF (trans) THEN
2323 CALL ensure_work_matrix(work, SIZE(krotj, 2), SIZE(sblock, 1))
2324 CALL dgemm('T', 'T', SIZE(krotj, 2), SIZE(sblock, 1), SIZE(krotj, 1), &
2325 1.0_dp, krotj, SIZE(krotj, 1), sblock, SIZE(sblock, 1), &
2326 0.0_dp, work, SIZE(work, 1))
2327 CALL dgemm('N', 'N', SIZE(work, 1), SIZE(kroti, 2), SIZE(work, 2), &
2328 fsign, work, SIZE(work, 1), kroti, SIZE(kroti, 1), &
2329 0.0_dp, sblock, SIZE(sblock, 1))
2330 ELSE
2331 CALL ensure_work_matrix(work, SIZE(kroti, 2), SIZE(sblock, 2))
2332 CALL dgemm('T', 'N', SIZE(kroti, 2), SIZE(sblock, 2), SIZE(kroti, 1), &
2333 1.0_dp, kroti, SIZE(kroti, 1), sblock, SIZE(sblock, 1), &
2334 0.0_dp, work, SIZE(work, 1))
2335 CALL dgemm('N', 'N', SIZE(work, 1), SIZE(krotj, 2), SIZE(work, 2), &
2336 fsign, work, SIZE(work, 1), krotj, SIZE(krotj, 1), &
2337 0.0_dp, sblock, SIZE(sblock, 1))
2338 END IF
2339 END DO
2340 CALL dbcsr_iterator_stop(iter)
2341 !
2342 END IF
2343 ELSE
2344 ! this is the identity operation, just copy the matrix
2345 CALL dbcsr_copy(smat, pmat)
2346 END IF
2347
2348 CALL timestop(handle)
2349
2350 END SUBROUTINE symtrans
2351
2352! **************************************************************************************************
2353!> \brief ...
2354!> \param mat ...
2355! **************************************************************************************************
2356 SUBROUTINE matprint(mat)
2357 TYPE(dbcsr_type), POINTER :: mat
2358
2359 INTEGER :: i, icol, iounit, irow
2360 REAL(kind=dp), DIMENSION(:, :), POINTER :: mblock
2361 TYPE(dbcsr_iterator_type) :: iter
2362
2363 iounit = cp_logger_get_default_io_unit()
2364 CALL dbcsr_iterator_start(iter, mat)
2365 DO WHILE (dbcsr_iterator_blocks_left(iter))
2366 CALL dbcsr_iterator_next_block(iter, irow, icol, mblock)
2367 !
2368 IF (iounit > 0) THEN
2369 WRITE (iounit, '(A,2I4)') 'BLOCK ', irow, icol
2370 DO i = 1, SIZE(mblock, 1)
2371 WRITE (iounit, '(8F12.6)') mblock(i, :)
2372 END DO
2373 END IF
2374 !
2375 END DO
2376 CALL dbcsr_iterator_stop(iter)
2377
2378 END SUBROUTINE matprint
2379! **************************************************************************************************
2380
2381END MODULE kpoint_methods
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:218
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:89
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:84
memory_utilities::reallocate
Definition memory_utilities.F:27
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::real_to_scaled
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
Definition cell_types.F:491
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_blacs_env::cp_blacs_env_create
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Definition cp_blacs_env.F:123
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:182
cp_cfm_types::cp_fm_to_cfm
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
Definition cp_cfm_types.F:840
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:125
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition cp_cfm_types.F:788
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:943
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1049
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:370
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:692
cp_dbcsr_api::dbcsr_replicate_all
subroutine, public dbcsr_replicate_all(matrix)
...
Definition cp_dbcsr_api.F:1155
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:807
cp_dbcsr_api::dbcsr_get_stored_coordinates
subroutine, public dbcsr_get_stored_coordinates(matrix, row, column, processor)
...
Definition cp_dbcsr_api.F:910
cp_dbcsr_api::dbcsr_distribute
subroutine, public dbcsr_distribute(matrix)
...
Definition cp_dbcsr_api.F:539
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:970
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1011
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1203
cp_dbcsr_api::dbcsr_distribution_get
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
Definition cp_dbcsr_api.F:571
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:111
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition cp_fm_basic_linalg.F:1948
cp_fm_pool_types
pool for for elements that are retained and released
Definition cp_fm_pool_types.F:14
cp_fm_pool_types::fm_pool_create_fm
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
Definition cp_fm_pool_types.F:185
cp_fm_pool_types::fm_pool_give_back_fm
subroutine, public fm_pool_give_back_fm(pool, element)
returns the element to the pool
Definition cp_fm_pool_types.F:224
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
Definition cp_fm_types.F:591
cp_fm_types::cp_fm_start_copy_general
subroutine, public cp_fm_start_copy_general(source, destination, para_env, info)
Initiates the copy operation: get distribution data, post MPI isend and irecvs.
Definition cp_fm_types.F:1645
cp_fm_types::cp_fm_cleanup_copy_general
subroutine, public cp_fm_cleanup_copy_general(info)
Completes the copy operation: wait for comms clean up MPI state.
Definition cp_fm_types.F:2023
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1087
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:169
cp_fm_types::cp_fm_finish_copy_general
subroutine, public cp_fm_finish_copy_general(destination, info)
Completes the copy operation: wait for comms, unpack, clean up MPI state.
Definition cp_fm_types.F:1959
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:972
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:515
cryssym
K-points and crystal symmetry routines.
Definition cryssym.F:12
cryssym::print_crys_symmetry
subroutine, public print_crys_symmetry(csym)
...
Definition cryssym.F:1886
cryssym::kpoint_gen
subroutine, public kpoint_gen(csym, nk, symm, shift, full_grid, gamma_centered, inversion_symmetry_only, use_spglib_reduction, use_spglib_backend)
...
Definition cryssym.F:253
cryssym::release_csym_type
subroutine, public release_csym_type(csym)
Release the CSYM type.
Definition cryssym.F:90
cryssym::kpoint_gen_general
subroutine, public kpoint_gen_general(csym, xkp_in, wkp_in, symm, full_grid, inversion_symmetry_only, use_spglib_reduction, use_spglib_backend)
Reduce an explicitly supplied GENERAL k-point set.
Definition cryssym.F:462
cryssym::print_kp_symmetry
subroutine, public print_kp_symmetry(csym)
...
Definition cryssym.F:1923
cryssym::crys_sym_gen
subroutine, public crys_sym_gen(csym, scoor, types, hmat, delta, iounit, use_spglib)
...
Definition cryssym.F:145
hairy_probes
Definition hairy_probes.F:8
hairy_probes::probe_occupancy_kp
subroutine, public probe_occupancy_kp(occ, fermi, kts, energies, rcoeff, icoeff, maxocc, probe, n, wk)
subroutine to calculate occupation number and 'Fermi' level using the
Definition hairy_probes.F:194
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::smear_fermi_dirac
integer, parameter, public smear_fermi_dirac
Definition input_constants.F:910
input_constants::smear_gaussian
integer, parameter, public smear_gaussian
Definition input_constants.F:910
input_constants::smear_mv
integer, parameter, public smear_mv
Definition input_constants.F:910
input_constants::smear_mp
integer, parameter, public smear_mp
Definition input_constants.F:910
input_cp2k_kpoints
function that build the kpoints section of the input
Definition input_cp2k_kpoints.F:14
input_cp2k_kpoints::use_spglib_kpoint_symmetry
integer, parameter, public use_spglib_kpoint_symmetry
Definition input_cp2k_kpoints.F:37
input_cp2k_kpoints::use_spglib_kpoint_backend
integer, parameter, public use_spglib_kpoint_backend
Definition input_cp2k_kpoints.F:37
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kpoint_methods
Routines needed for kpoint calculation.
Definition kpoint_methods.F:15
kpoint_methods::kpoint_initialize_mo_set
subroutine, public kpoint_initialize_mo_set(kpoint)
...
Definition kpoint_methods.F:939
kpoint_methods::rskp_transform
subroutine, public rskp_transform(rmatrix, cmatrix, rsmat, ispin, xkp, cell_to_index, sab_nl, is_complex, rs_sign)
Transformation of real space matrices to a kpoint.
Definition kpoint_methods.F:1113
kpoint_methods::kpoint_init_cell_index
subroutine, public kpoint_init_cell_index(kpoint, sab_nl, para_env, nimages)
Generates the mapping of cell indices and linear RS index CELL (0,0,0) is always mapped to index 1.
Definition kpoint_methods.F:979
kpoint_methods::kpoint_initialize_mos
subroutine, public kpoint_initialize_mos(kpoint, mos, added_mos, for_aux_fit)
Initialize a set of MOs and density matrix for each kpoint (kpoint group)
Definition kpoint_methods.F:786
kpoint_methods::kpoint_initialize
subroutine, public kpoint_initialize(kpoint, particle_set, cell)
Generate the kpoints and initialize the kpoint environment.
Definition kpoint_methods.F:120
kpoint_methods::kpoint_density_transform
subroutine, public kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwork, for_aux_fit, pmat_ext)
generate real space density matrices in DBCSR format
Definition kpoint_methods.F:1663
kpoint_methods::kpoint_density_matrices
subroutine, public kpoint_density_matrices(kpoint, energy_weighted, for_aux_fit)
Calculate kpoint density matrices (rho(k), owned by kpoint groups)
Definition kpoint_methods.F:1452
kpoint_methods::lowdin_kp_trans
subroutine, public lowdin_kp_trans(kpoint, pmat_diag)
Calculate Lowdin transformation of density matrix S^1/2 P S^1/2 Integrate diagonal elements over k-po...
Definition kpoint_methods.F:1572
kpoint_methods::kpoint_env_initialize
subroutine, public kpoint_env_initialize(kpoint, para_env, blacs_env, with_aux_fit)
Initialize the kpoint environment.
Definition kpoint_methods.F:598
kpoint_methods::kpoint_set_mo_occupation
subroutine, public kpoint_set_mo_occupation(kpoint, smear, probe)
Given the eigenvalues of all kpoints, calculates the occupation numbers.
Definition kpoint_methods.F:1213
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::kpoint_sym_create
subroutine, public kpoint_sym_create(kp_sym)
Create a single kpoint symmetry environment.
Definition kpoint_types.F:1085
kpoint_types::kpoint_env_create
subroutine, public kpoint_env_create(kp_env)
Create a single kpoint environment.
Definition kpoint_types.F:998
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:524
kpsym
K-points and crystal symmetry routines based on.
Definition kpsym.F:28
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::inv_3x3
pure real(kind=dp) function, dimension(3, 3), public inv_3x3(a)
Returns the inverse of the 3 x 3 matrix a.
Definition mathlib.F:524
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_matrix_pools
wrapper for the pools of matrixes
Definition qs_matrix_pools.F:14
qs_matrix_pools::mpools_create
subroutine, public mpools_create(mpools)
creates a mpools
Definition qs_matrix_pools.F:184
qs_matrix_pools::mpools_rebuild_fm_pools
subroutine, public mpools_rebuild_fm_pools(mpools, mos, blacs_env, para_env, nmosub)
rebuilds the pools of the (ao x mo, ao x ao , mo x mo) full matrixes
Definition qs_matrix_pools.F:209
qs_matrix_pools::mpools_get
subroutine, public mpools_get(mpools, ao_mo_fm_pools, ao_ao_fm_pools, mo_mo_fm_pools, ao_mosub_fm_pools, mosub_mosub_fm_pools, maxao_maxmo_fm_pool, maxao_maxao_fm_pool, maxmo_maxmo_fm_pool)
returns various attributes of the mpools (notably the pools contained in it)
Definition qs_matrix_pools.F:139
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::set_mo_set
subroutine, public set_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, uniform_occupation, kts, mu, flexible_electron_count)
Set the components of a MO set data structure.
Definition qs_mo_types.F:454
qs_mo_types::allocate_mo_set
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
Definition qs_mo_types.F:207
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:399
qs_mo_types::init_mo_set
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
Definition qs_mo_types.F:246
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::get_neighbor_list_set_p
subroutine, public get_neighbor_list_set_p(neighbor_list_sets, nlist, symmetric)
Return the components of the first neighbor list set.
Definition qs_neighbor_list_types.F:995
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:549
scf_control_types
parameters that control an scf iteration
Definition scf_control_types.F:17
smearing_utils
Unified smearing module supporting four methods: smear_fermi_dirac — Fermi-Dirac distribution smear_g...
Definition smearing_utils.F:24
smearing_utils::smearkp
subroutine, public smearkp(f, mu, kts, e, nel, wk, sigma, maxocc, method)
Bisection search for mu given a target electron count (k-point case, single spin channel or spin-dege...
Definition smearing_utils.F:464
smearing_utils::smearkp2
subroutine, public smearkp2(f, mu, kts, e, nel, wk, sigma, method)
Bisection search for mu (k-point, spin-polarised with a shared chemical potential across both spin ch...
Definition smearing_utils.F:576
util
All kind of helpful little routines.
Definition util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition util.F:333
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:69
cp_control_types::hairy_probes_type
Definition cp_control_types.F:48
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_pool_types::cp_fm_pool_p_type
to create arrays of pools
Definition cp_fm_pool_types.F:70
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::copy_info_type
Stores the state of a copy between cp_fm_start_copy_general and cp_fm_finish_copy_general.
Definition cp_fm_types.F:141
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:115
cryssym::csym_type
CSM type.
Definition cryssym.F:44
kpoint_types::kind_rotmat_type
Rotation matrices for basis sets.
Definition kpoint_types.F:105
kpoint_types::kpoint_env_p_type
Definition kpoint_types.F:96
kpoint_types::kpoint_env_type
Keeps information about a specific k-point.
Definition kpoint_types.F:82
kpoint_types::kpoint_sym_type
Keeps symmetry information about a specific k-point.
Definition kpoint_types.F:120
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:170
message_passing::mp_cart_type
Definition message_passing.F:720
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
particle_types::particle_type
Definition particle_types.F:35
qs_matrix_pools::qs_matrix_pools_type
container for the pools of matrixes used by qs
Definition qs_matrix_pools.F:66
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
scf_control_types::smear_type
contains the parameters needed by a scf run
Definition scf_control_types.F:91