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