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 : pbc_stable,&
19 : real_to_scaled
20 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
21 : cp_blacs_env_type
22 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,&
23 : cp_cfm_gemm,&
24 : cp_cfm_scale_and_add,&
25 : cp_cfm_transpose
26 : USE cp_cfm_types, ONLY: cp_cfm_create,&
27 : cp_cfm_get_info,&
28 : cp_cfm_release,&
29 : cp_cfm_to_fm,&
30 : cp_cfm_type,&
31 : cp_fm_to_cfm
32 : USE cp_control_types, ONLY: hairy_probes_type
33 : USE cp_dbcsr_api, ONLY: &
34 : dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribute, &
35 : dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_get_block_p, dbcsr_get_info, &
36 : dbcsr_get_readonly_block_p, dbcsr_get_stored_coordinates, dbcsr_iterator_blocks_left, &
37 : dbcsr_iterator_next_block, dbcsr_iterator_readonly_start, dbcsr_iterator_start, &
38 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_replicate_all, dbcsr_set, &
39 : dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
40 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot
41 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
42 : USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,&
43 : cp_dbcsr_sm_fm_multiply
44 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
45 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
46 : fm_pool_create_fm,&
47 : fm_pool_give_back_fm
48 : USE cp_fm_struct, ONLY: cp_fm_struct_type
49 : USE cp_fm_types, ONLY: &
50 : copy_info_type, cp_fm_cleanup_copy_general, cp_fm_create, cp_fm_finish_copy_general, &
51 : cp_fm_get_diag, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_release, &
52 : cp_fm_start_copy_general, cp_fm_to_fm, cp_fm_type
53 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
54 : USE cryssym, ONLY: crys_sym_gen,&
55 : csym_type,&
56 : kpoint_gen,&
57 : kpoint_gen_general,&
58 : print_crys_symmetry,&
59 : print_kp_symmetry,&
60 : release_csym_type
61 : USE hairy_probes, ONLY: probe_occupancy_kp
62 : USE input_constants, ONLY: smear_fermi_dirac,&
63 : smear_gaussian,&
64 : smear_mp,&
65 : smear_mv
66 : USE input_cp2k_kpoints, ONLY: use_spglib_kpoint_backend,&
67 : use_spglib_kpoint_symmetry
68 : USE kinds, ONLY: dp,&
69 : int_8
70 : USE kpoint_lattice_fft, ONLY: cell_to_k_grid_fft,&
71 : regular_kpoint_grid
72 : USE kpoint_types, ONLY: get_kpoint_info,&
73 : kind_rotmat_type,&
74 : kpoint_env_create,&
75 : kpoint_env_p_type,&
76 : kpoint_env_type,&
77 : kpoint_sym_create,&
78 : kpoint_sym_type,&
79 : kpoint_type
80 : USE mathconstants, ONLY: twopi
81 : USE mathlib, ONLY: inv_3x3
82 : USE memory_utilities, ONLY: reallocate
83 : USE message_passing, ONLY: mp_cart_type,&
84 : mp_para_env_type
85 : USE parallel_gemm_api, ONLY: parallel_gemm
86 : USE particle_types, ONLY: particle_type
87 : USE qs_matrix_pools, ONLY: mpools_create,&
88 : mpools_get,&
89 : mpools_rebuild_fm_pools,&
90 : qs_matrix_pools_type
91 : USE qs_mo_types, ONLY: allocate_mo_set,&
92 : get_mo_set,&
93 : init_mo_set,&
94 : mo_set_type,&
95 : set_mo_set
96 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
97 : get_neighbor_list_set_p,&
98 : neighbor_list_iterate,&
99 : neighbor_list_iterator_create,&
100 : neighbor_list_iterator_p_type,&
101 : neighbor_list_iterator_release,&
102 : neighbor_list_set_p_type
103 : USE scf_control_types, ONLY: smear_type
104 : USE smearing_utils, ONLY: Smearkp,&
105 : Smearkp2
106 : USE util, ONLY: get_limit,&
107 : sort
108 :
109 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, &
110 : !$ omp_get_thread_num
111 : #include "./base/base_uses.f90"
112 :
113 : IMPLICIT NONE
114 :
115 : PRIVATE
116 :
117 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpoint_methods'
118 :
119 : TYPE, PUBLIC :: kp_transform_plan_type
120 : INTEGER :: nentry = 0, ngroup = 0
121 : LOGICAL :: symmetric = .FALSE.
122 : INTEGER, ALLOCATABLE, DIMENSION(:) :: col, col_offset, group_start, image, row, row_offset
123 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: cell
124 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: symmetry_sign
125 : END TYPE kp_transform_plan_type
126 :
127 : PUBLIC :: kpoint_initialize, kpoint_env_initialize, kpoint_initialize_mos, kpoint_initialize_mo_set
128 : PUBLIC :: kpoint_init_cell_index, kpoint_set_mo_occupation
129 : PUBLIC :: kpoint_smearing_edge_status
130 : PUBLIC :: kpoint_density_matrices, kpoint_density_transform
131 : PUBLIC :: kpoint_ot_energy_weighted_density, kpoint_ot_energy_weighted_matrices
132 : PUBLIC :: kp_transform_plan_create, kp_transform_plan_release
133 : PUBLIC :: rskp_transform, lowdin_kp_trans, lowdin_kp_mo_coeff
134 : PUBLIC :: rskp_grid_type, rskp_transform_grid_prepare, rskp_transform_grid_extract, &
135 : rskp_transform_grid_release
136 :
137 : ! **************************************************************************************************
138 :
139 : TYPE rskp_grid_type
140 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: values
141 : INTEGER, ALLOCATABLE, DIMENSION(:) :: block_col, block_nelem, block_offset, block_row
142 : LOGICAL :: ready = .FALSE.
143 : END TYPE rskp_grid_type
144 :
145 : ! **************************************************************************************************
146 :
147 : CONTAINS
148 :
149 : ! **************************************************************************************************
150 : !> \brief Generate the kpoints and initialize the kpoint environment
151 : !> \param kpoint The kpoint environment
152 : !> \param particle_set Particle types and coordinates
153 : !> \param cell Computational cell information
154 : ! **************************************************************************************************
155 11750 : SUBROUTINE kpoint_initialize(kpoint, particle_set, cell)
156 :
157 : TYPE(kpoint_type), POINTER :: kpoint
158 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
159 : TYPE(cell_type), POINTER :: cell
160 :
161 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_initialize'
162 :
163 : INTEGER :: handle, i, ic, ik, iounit, ir, ira, is, &
164 : isign, j, natom, nkind, nr, ns
165 11750 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atype
166 11750 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: agauge
167 : INTEGER, DIMENSION(3, 3) :: frot, krot
168 : LOGICAL :: spez
169 : REAL(KIND=dp) :: eps_kpoint, wsum
170 11750 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coord, scoord
171 : REAL(KIND=dp), DIMENSION(3) :: diff, kgvec, r_pbc, scoord_pbc, srot
172 : REAL(KIND=dp), DIMENSION(3, 3) :: srotmat
173 11750 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp_full
174 11750 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp_full
175 223250 : TYPE(csym_type) :: crys_sym
176 : TYPE(kpoint_sym_type), POINTER :: kpsym
177 :
178 11750 : CALL timeset(routineN, handle)
179 :
180 11750 : CPASSERT(ASSOCIATED(kpoint))
181 :
182 11770 : SELECT CASE (kpoint%kp_scheme)
183 : CASE ("NONE")
184 : ! do nothing
185 : CASE ("GAMMA")
186 20 : kpoint%nkp = 1
187 20 : ALLOCATE (kpoint%xkp(3, 1), kpoint%wkp(1))
188 80 : kpoint%xkp(1:3, 1) = 0.0_dp
189 20 : kpoint%wkp(1) = 1.0_dp
190 40 : ALLOCATE (kpoint%kp_sym(1))
191 20 : NULLIFY (kpoint%kp_sym(1)%kpoint_sym)
192 20 : CALL kpoint_sym_create(kpoint%kp_sym(1)%kpoint_sym)
193 : CASE ("MONKHORST-PACK", "MACDONALD")
194 :
195 3316 : IF (.NOT. kpoint%symmetry) THEN
196 : ! we set up a random molecule to avoid any possible symmetry
197 182 : natom = 10
198 182 : ALLOCATE (coord(3, natom), scoord(3, natom), atype(natom))
199 2002 : DO i = 1, natom
200 1820 : atype(i) = i
201 1820 : coord(1, i) = SIN(i*0.12345_dp)
202 1820 : coord(2, i) = COS(i*0.23456_dp)
203 1820 : coord(3, i) = SIN(i*0.34567_dp)
204 2002 : CALL real_to_scaled(scoord(1:3, i), coord(1:3, i), cell)
205 : END DO
206 : ELSE
207 3134 : natom = SIZE(particle_set)
208 15670 : ALLOCATE (scoord(3, natom), atype(natom))
209 17104 : DO i = 1, natom
210 13970 : CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=atype(i))
211 17104 : CALL real_to_scaled(scoord(1:3, i), particle_set(i)%r(1:3), cell)
212 : END DO
213 : END IF
214 3316 : IF (kpoint%verbose) THEN
215 2340 : iounit = cp_logger_get_default_io_unit()
216 : ELSE
217 976 : iounit = -1
218 : END IF
219 : ! kind type list
220 9948 : ALLOCATE (kpoint%atype(natom))
221 19106 : kpoint%atype = atype
222 : ! Match the atom images used by CP2K's periodic neighbor lists.
223 9948 : ALLOCATE (agauge(3, natom))
224 3316 : agauge = 0
225 3316 : IF (kpoint%symmetry) THEN
226 17104 : DO i = 1, natom
227 13970 : r_pbc(1:3) = pbc_stable(particle_set(i)%r(1:3), cell)
228 13970 : CALL real_to_scaled(scoord_pbc, r_pbc, cell)
229 59014 : agauge(1:3, i) = NINT(scoord_pbc(1:3) - scoord(1:3, i))
230 : END DO
231 : END IF
232 :
233 : CALL crys_sym_gen(crys_sym, scoord, atype, cell%hmat, delta=kpoint%eps_geo, iounit=iounit, &
234 3316 : use_spglib=kpoint%symmetry)
235 : CALL kpoint_gen(crys_sym, kpoint%nkp_grid, symm=kpoint%symmetry, shift=kpoint%kp_shift, &
236 : full_grid=kpoint%full_grid, gamma_centered=kpoint%gamma_centered, &
237 : inversion_symmetry_only=kpoint%inversion_symmetry_only, &
238 : use_spglib_reduction= &
239 : kpoint%symmetry_reduction_method == use_spglib_kpoint_symmetry, &
240 3316 : use_spglib_backend=kpoint%symmetry_backend == use_spglib_kpoint_backend)
241 3316 : IF (crys_sym%inversion_only) kpoint%inversion_symmetry_only = .TRUE.
242 3316 : kpoint%nkp = crys_sym%nkpoint
243 16580 : ALLOCATE (kpoint%xkp(3, kpoint%nkp), kpoint%wkp(kpoint%nkp))
244 14984 : wsum = SUM(crys_sym%wkpoint)
245 14984 : DO ik = 1, kpoint%nkp
246 46672 : kpoint%xkp(1:3, ik) = crys_sym%xkpoint(1:3, ik)
247 14984 : kpoint%wkp(ik) = crys_sym%wkpoint(ik)/wsum
248 : END DO
249 :
250 3316 : eps_kpoint = MAX(1.e-12_dp, 10.0_dp*kpoint%eps_geo)
251 : ! print output
252 3316 : IF (kpoint%symmetry) CALL print_crys_symmetry(crys_sym)
253 3316 : IF (kpoint%symmetry) CALL print_kp_symmetry(crys_sym)
254 :
255 : ! transfer symmetry information
256 21616 : ALLOCATE (kpoint%kp_sym(kpoint%nkp))
257 14984 : DO ik = 1, kpoint%nkp
258 11668 : NULLIFY (kpoint%kp_sym(ik)%kpoint_sym)
259 11668 : CALL kpoint_sym_create(kpoint%kp_sym(ik)%kpoint_sym)
260 11668 : kpsym => kpoint%kp_sym(ik)%kpoint_sym
261 : IF (crys_sym%nrtot > 0 .AND. .NOT. crys_sym%fullgrid .AND. &
262 14984 : crys_sym%istriz == 1 .AND. .NOT. crys_sym%inversion_only) THEN
263 : ! set up the symmetrization information
264 3982 : kpsym%nwght = NINT(crys_sym%wkpoint(ik))
265 3982 : ns = kpsym%nwght
266 : !
267 3982 : IF (ns > 1) THEN
268 81538 : DO is = 1, SIZE(crys_sym%kplink, 2)
269 81538 : IF (crys_sym%kplink(2, is) == ik) THEN
270 557368 : DO ic = 1, crys_sym%nrtot
271 42541816 : srotmat = MATMUL(cell%h_inv, MATMUL(crys_sym%rt(1:3, 1:3, ic), cell%hmat))
272 7000552 : frot(1:3, 1:3) = NINT(srotmat(1:3, 1:3))
273 14001104 : krot(1:3, 1:3) = NINT(TRANSPOSE(inv_3x3(REAL(frot(1:3, 1:3), KIND=dp))))
274 1634376 : DO isign = 1, 2
275 1077008 : ir = MERGE(crys_sym%ibrot(ic), -crys_sym%ibrot(ic), isign == 1)
276 1077008 : IF (ir == crys_sym%kpop(is)) CYCLE
277 : kgvec(1:3) = crys_sym%kpmesh(1:3, is) - &
278 : MATMUL(REAL(MERGE(krot(1:3, 1:3), -krot(1:3, 1:3), &
279 : isign == 1), KIND=dp), &
280 29628032 : kpoint%xkp(1:3, ik))
281 4232576 : diff(1:3) = kgvec(1:3) - ANINT(kgvec(1:3))
282 2218356 : IF (ALL(ABS(diff(1:3)) < eps_kpoint)) ns = ns + 1
283 : END DO
284 : END DO
285 : END IF
286 : END DO
287 3572 : kpsym%apply_symmetry = .TRUE.
288 3572 : natom = SIZE(particle_set)
289 10716 : ALLOCATE (kpsym%rot(3, 3, ns))
290 10716 : ALLOCATE (kpsym%xkp(3, ns))
291 10716 : ALLOCATE (kpsym%rotp(ns))
292 14288 : ALLOCATE (kpsym%f0(natom, ns))
293 14288 : ALLOCATE (kpsym%fcell(3, natom, ns))
294 10716 : ALLOCATE (kpsym%fcell_gauge(3, natom, ns))
295 7144 : ALLOCATE (kpsym%phase_mode(ns))
296 101760 : kpsym%phase_mode = 0
297 14288 : ALLOCATE (kpsym%kgphase(natom, ns))
298 3572 : nr = 0
299 81538 : DO is = 1, SIZE(crys_sym%kplink, 2)
300 81538 : IF (crys_sym%kplink(2, is) == ik) THEN
301 18864 : nr = nr + 1
302 18864 : ir = crys_sym%kpop(is)
303 18864 : ira = ABS(ir)
304 100776 : DO ic = 1, crys_sym%nrtot
305 100776 : IF (crys_sym%ibrot(ic) == ira) THEN
306 18864 : kpsym%rotp(nr) = ir
307 245232 : kpsym%rot(1:3, 1:3, nr) = crys_sym%rt(1:3, 1:3, ic)
308 1490256 : srotmat = MATMUL(cell%h_inv, MATMUL(crys_sym%rt(1:3, 1:3, ic), cell%hmat))
309 245232 : frot(1:3, 1:3) = NINT(srotmat(1:3, 1:3))
310 75456 : kpsym%xkp(1:3, nr) = crys_sym%kpmesh(1:3, is)
311 490464 : krot(1:3, 1:3) = NINT(TRANSPOSE(inv_3x3(REAL(frot(1:3, 1:3), KIND=dp))))
312 129960 : IF (ir < 0) krot(1:3, 1:3) = -krot(1:3, 1:3)
313 : kgvec(1:3) = kpsym%xkp(1:3, nr) - &
314 : MATMUL(REAL(krot(1:3, 1:3), KIND=dp), &
315 528192 : kpoint%xkp(1:3, ik))
316 75456 : kgvec(1:3) = ANINT(kgvec(1:3))
317 135952 : kpsym%f0(1:natom, nr) = crys_sym%f0(1:natom, ic)
318 135952 : DO j = 1, natom
319 1873408 : srot(1:3) = MATMUL(srotmat, scoord(1:3, j)) + crys_sym%vt(1:3, ic)
320 : kpsym%fcell(1:3, j, nr) = &
321 468352 : NINT(srot(1:3) - scoord(1:3, kpsym%f0(j, nr)))
322 : kpsym%fcell_gauge(1:3, j, nr) = &
323 : NINT(srot(1:3) - scoord(1:3, kpsym%f0(j, nr))) + &
324 : MATMUL(frot(1:3, 1:3), agauge(1:3, j)) - &
325 2224672 : agauge(1:3, kpsym%f0(j, nr))
326 : kpsym%kgphase(j, nr) = DOT_PRODUCT(kgvec(1:3), &
327 : scoord(1:3, j) + &
328 487216 : REAL(agauge(1:3, j), KIND=dp))
329 : END DO
330 : EXIT
331 : END IF
332 : END DO
333 18864 : CPASSERT(ic <= crys_sym%nrtot)
334 : END IF
335 : END DO
336 81538 : DO is = 1, SIZE(crys_sym%kplink, 2)
337 81538 : IF (crys_sym%kplink(2, is) == ik) THEN
338 557368 : DO ic = 1, crys_sym%nrtot
339 42541816 : srotmat = MATMUL(cell%h_inv, MATMUL(crys_sym%rt(1:3, 1:3, ic), cell%hmat))
340 7000552 : frot(1:3, 1:3) = NINT(srotmat(1:3, 1:3))
341 14001104 : krot(1:3, 1:3) = NINT(TRANSPOSE(inv_3x3(REAL(frot(1:3, 1:3), KIND=dp))))
342 1634376 : DO isign = 1, 2
343 1077008 : ir = MERGE(crys_sym%ibrot(ic), -crys_sym%ibrot(ic), isign == 1)
344 1077008 : IF (ir == crys_sym%kpop(is)) CYCLE
345 : kgvec(1:3) = crys_sym%kpmesh(1:3, is) - &
346 : MATMUL(REAL(MERGE(krot(1:3, 1:3), -krot(1:3, 1:3), &
347 : isign == 1), KIND=dp), &
348 29628032 : kpoint%xkp(1:3, ik))
349 4232576 : diff(1:3) = kgvec(1:3) - ANINT(kgvec(1:3))
350 2218356 : IF (ALL(ABS(diff(1:3)) < eps_kpoint)) THEN
351 79324 : nr = nr + 1
352 79324 : kpsym%rotp(nr) = ir
353 1031212 : kpsym%rot(1:3, 1:3, nr) = crys_sym%rt(1:3, 1:3, ic)
354 317296 : kpsym%xkp(1:3, nr) = crys_sym%kpmesh(1:3, is)
355 317296 : kgvec(1:3) = ANINT(kgvec(1:3))
356 512720 : kpsym%f0(1:natom, nr) = crys_sym%f0(1:natom, ic)
357 512720 : DO j = 1, natom
358 6934336 : srot(1:3) = MATMUL(srotmat, scoord(1:3, j)) + crys_sym%vt(1:3, ic)
359 : kpsym%fcell(1:3, j, nr) = &
360 1733584 : NINT(srot(1:3) - scoord(1:3, kpsym%f0(j, nr)))
361 : kpsym%fcell_gauge(1:3, j, nr) = &
362 : NINT(srot(1:3) - scoord(1:3, kpsym%f0(j, nr))) + &
363 : MATMUL(frot(1:3, 1:3), agauge(1:3, j)) - &
364 8234524 : agauge(1:3, kpsym%f0(j, nr))
365 : kpsym%kgphase(j, nr) = DOT_PRODUCT(kgvec(1:3), &
366 : scoord(1:3, j) + &
367 1812908 : REAL(agauge(1:3, j), KIND=dp))
368 : END DO
369 : END IF
370 : END DO
371 : END DO
372 : END IF
373 : END DO
374 3572 : kpsym%nwred = nr
375 : END IF
376 : END IF
377 : END DO
378 3316 : IF (kpoint%symmetry) THEN
379 17104 : nkind = MAXVAL(atype)
380 3134 : ns = crys_sym%nrtot
381 40762 : ALLOCATE (kpoint%kind_rotmat(ns, nkind))
382 26390 : DO i = 1, ns
383 51690 : DO j = 1, nkind
384 48556 : NULLIFY (kpoint%kind_rotmat(i, j)%rmat)
385 : END DO
386 : END DO
387 7942 : ALLOCATE (kpoint%ibrot(ns))
388 26390 : kpoint%ibrot(1:ns) = crys_sym%ibrot(1:ns)
389 : END IF
390 :
391 3316 : CALL release_csym_type(crys_sym)
392 3316 : DEALLOCATE (scoord, atype)
393 3316 : DEALLOCATE (agauge)
394 :
395 : CASE ("GENERAL")
396 : NULLIFY (xkp_full, wkp_full)
397 36 : IF (ASSOCIATED(kpoint%xkp_input)) THEN
398 36 : xkp_full => kpoint%xkp_input
399 36 : wkp_full => kpoint%wkp_input
400 : ELSE
401 0 : xkp_full => kpoint%xkp
402 0 : wkp_full => kpoint%wkp
403 : END IF
404 36 : CPASSERT(ASSOCIATED(xkp_full))
405 36 : CPASSERT(ASSOCIATED(wkp_full))
406 36 : IF (.NOT. ASSOCIATED(kpoint%xkp_input)) THEN
407 0 : ALLOCATE (kpoint%xkp_input(3, SIZE(wkp_full)), kpoint%wkp_input(SIZE(wkp_full)))
408 0 : kpoint%xkp_input(1:3, 1:SIZE(wkp_full)) = xkp_full(1:3, 1:SIZE(wkp_full))
409 0 : kpoint%wkp_input(1:SIZE(wkp_full)) = wkp_full(1:SIZE(wkp_full))
410 0 : xkp_full => kpoint%xkp_input
411 0 : wkp_full => kpoint%wkp_input
412 : END IF
413 36 : IF (.NOT. kpoint%symmetry) THEN
414 10 : IF (.NOT. ASSOCIATED(kpoint%xkp)) THEN
415 0 : kpoint%nkp = SIZE(wkp_full)
416 0 : ALLOCATE (kpoint%xkp(3, kpoint%nkp), kpoint%wkp(kpoint%nkp))
417 0 : kpoint%xkp(1:3, 1:kpoint%nkp) = xkp_full(1:3, 1:kpoint%nkp)
418 0 : kpoint%wkp(1:kpoint%nkp) = wkp_full(1:kpoint%nkp)
419 : END IF
420 : ! default: no symmetry settings
421 74 : ALLOCATE (kpoint%kp_sym(kpoint%nkp))
422 54 : DO i = 1, kpoint%nkp
423 44 : NULLIFY (kpoint%kp_sym(i)%kpoint_sym)
424 54 : CALL kpoint_sym_create(kpoint%kp_sym(i)%kpoint_sym)
425 : END DO
426 : ELSE
427 26 : IF (kpoint%verbose) THEN
428 16 : iounit = cp_logger_get_default_io_unit()
429 : ELSE
430 10 : iounit = -1
431 : END IF
432 26 : natom = SIZE(particle_set)
433 130 : ALLOCATE (scoord(3, natom), atype(natom))
434 234 : DO i = 1, natom
435 208 : CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=atype(i))
436 234 : CALL real_to_scaled(scoord(1:3, i), particle_set(i)%r(1:3), cell)
437 : END DO
438 52 : ALLOCATE (kpoint%atype(natom))
439 234 : kpoint%atype = atype
440 78 : ALLOCATE (agauge(3, natom))
441 234 : DO i = 1, natom
442 208 : r_pbc(1:3) = pbc_stable(particle_set(i)%r(1:3), cell)
443 208 : CALL real_to_scaled(scoord_pbc, r_pbc, cell)
444 858 : agauge(1:3, i) = NINT(scoord_pbc(1:3) - scoord(1:3, i))
445 : END DO
446 :
447 : CALL crys_sym_gen(crys_sym, scoord, atype, cell%hmat, delta=kpoint%eps_geo, iounit=iounit, &
448 : use_spglib=(kpoint%symmetry_backend == use_spglib_kpoint_backend .OR. &
449 30 : kpoint%symmetry_reduction_method == use_spglib_kpoint_symmetry))
450 : CALL kpoint_gen_general(crys_sym, xkp_full, wkp_full, symm=kpoint%symmetry, &
451 : full_grid=kpoint%full_grid, &
452 : inversion_symmetry_only=kpoint%inversion_symmetry_only, &
453 : use_spglib_reduction= &
454 : kpoint%symmetry_reduction_method == use_spglib_kpoint_symmetry, &
455 26 : use_spglib_backend=kpoint%symmetry_backend == use_spglib_kpoint_backend)
456 26 : IF (crys_sym%inversion_only) kpoint%inversion_symmetry_only = .TRUE.
457 26 : IF (ASSOCIATED(kpoint%xkp)) THEN
458 26 : DEALLOCATE (kpoint%xkp)
459 26 : NULLIFY (kpoint%xkp)
460 : END IF
461 26 : IF (ASSOCIATED(kpoint%wkp)) THEN
462 26 : DEALLOCATE (kpoint%wkp)
463 26 : NULLIFY (kpoint%wkp)
464 : END IF
465 26 : kpoint%nkp = crys_sym%nkpoint
466 130 : ALLOCATE (kpoint%xkp(3, kpoint%nkp), kpoint%wkp(kpoint%nkp))
467 52 : wsum = SUM(crys_sym%wkpoint)
468 52 : DO ik = 1, kpoint%nkp
469 104 : kpoint%xkp(1:3, ik) = crys_sym%xkpoint(1:3, ik)
470 52 : kpoint%wkp(ik) = crys_sym%wkpoint(ik)/wsum
471 : END DO
472 :
473 26 : eps_kpoint = MAX(1.e-12_dp, 10.0_dp*kpoint%eps_geo)
474 26 : CALL print_crys_symmetry(crys_sym)
475 26 : CALL print_kp_symmetry(crys_sym)
476 :
477 104 : ALLOCATE (kpoint%kp_sym(kpoint%nkp))
478 52 : DO ik = 1, kpoint%nkp
479 26 : NULLIFY (kpoint%kp_sym(ik)%kpoint_sym)
480 26 : CALL kpoint_sym_create(kpoint%kp_sym(ik)%kpoint_sym)
481 26 : kpsym => kpoint%kp_sym(ik)%kpoint_sym
482 : IF (crys_sym%nrtot > 0 .AND. .NOT. crys_sym%fullgrid .AND. &
483 52 : crys_sym%istriz == 1 .AND. .NOT. crys_sym%inversion_only) THEN
484 26 : kpsym%nwght = NINT(crys_sym%wkpoint(ik))
485 26 : ns = kpsym%nwght
486 26 : IF (ns > 1) THEN
487 234 : DO is = 1, SIZE(crys_sym%kplink, 2)
488 234 : IF (crys_sym%kplink(2, is) == ik) THEN
489 10192 : DO ic = 1, crys_sym%nrtot
490 788736 : srotmat = MATMUL(cell%h_inv, MATMUL(crys_sym%rt(1:3, 1:3, ic), cell%hmat))
491 129792 : frot(1:3, 1:3) = NINT(srotmat(1:3, 1:3))
492 259584 : krot(1:3, 1:3) = NINT(TRANSPOSE(inv_3x3(REAL(frot(1:3, 1:3), KIND=dp))))
493 30160 : DO isign = 1, 2
494 19968 : ir = MERGE(crys_sym%ibrot(ic), -crys_sym%ibrot(ic), isign == 1)
495 19968 : IF (ir == crys_sym%kpop(is)) CYCLE
496 : kgvec(1:3) = crys_sym%kpmesh(1:3, is) - &
497 : MATMUL(REAL(MERGE(krot(1:3, 1:3), -krot(1:3, 1:3), &
498 : isign == 1), KIND=dp), &
499 553280 : kpoint%xkp(1:3, ik))
500 79040 : diff(1:3) = kgvec(1:3) - ANINT(kgvec(1:3))
501 46592 : IF (ALL(ABS(diff(1:3)) < eps_kpoint)) ns = ns + 1
502 : END DO
503 : END DO
504 : END IF
505 : END DO
506 26 : kpsym%apply_symmetry = .TRUE.
507 78 : ALLOCATE (kpsym%rot(3, 3, ns))
508 78 : ALLOCATE (kpsym%xkp(3, ns))
509 78 : ALLOCATE (kpsym%rotp(ns))
510 104 : ALLOCATE (kpsym%f0(natom, ns))
511 104 : ALLOCATE (kpsym%fcell(3, natom, ns))
512 78 : ALLOCATE (kpsym%fcell_gauge(3, natom, ns))
513 52 : ALLOCATE (kpsym%phase_mode(ns))
514 2522 : kpsym%phase_mode = 0
515 104 : ALLOCATE (kpsym%kgphase(natom, ns))
516 26 : nr = 0
517 234 : DO is = 1, SIZE(crys_sym%kplink, 2)
518 234 : IF (crys_sym%kplink(2, is) == ik) THEN
519 208 : nr = nr + 1
520 208 : ir = crys_sym%kpop(is)
521 208 : ira = ABS(ir)
522 628 : DO ic = 1, crys_sym%nrtot
523 628 : IF (crys_sym%ibrot(ic) == ira) THEN
524 208 : kpsym%rotp(nr) = ir
525 2704 : kpsym%rot(1:3, 1:3, nr) = crys_sym%rt(1:3, 1:3, ic)
526 16432 : srotmat = MATMUL(cell%h_inv, MATMUL(crys_sym%rt(1:3, 1:3, ic), cell%hmat))
527 2704 : frot(1:3, 1:3) = NINT(srotmat(1:3, 1:3))
528 832 : kpsym%xkp(1:3, nr) = crys_sym%kpmesh(1:3, is)
529 5408 : krot(1:3, 1:3) = NINT(TRANSPOSE(inv_3x3(REAL(frot(1:3, 1:3), KIND=dp))))
530 1456 : IF (ir < 0) krot(1:3, 1:3) = -krot(1:3, 1:3)
531 : kgvec(1:3) = kpsym%xkp(1:3, nr) - &
532 : MATMUL(REAL(krot(1:3, 1:3), KIND=dp), &
533 5824 : kpoint%xkp(1:3, ik))
534 832 : kgvec(1:3) = ANINT(kgvec(1:3))
535 1872 : kpsym%f0(1:natom, nr) = crys_sym%f0(1:natom, ic)
536 1872 : DO j = 1, natom
537 26624 : srot(1:3) = MATMUL(srotmat, scoord(1:3, j)) + crys_sym%vt(1:3, ic)
538 : kpsym%fcell(1:3, j, nr) = &
539 6656 : NINT(srot(1:3) - scoord(1:3, kpsym%f0(j, nr)))
540 : kpsym%fcell_gauge(1:3, j, nr) = &
541 : NINT(srot(1:3) - scoord(1:3, kpsym%f0(j, nr))) + &
542 : MATMUL(frot(1:3, 1:3), agauge(1:3, j)) - &
543 31616 : agauge(1:3, kpsym%f0(j, nr))
544 : kpsym%kgphase(j, nr) = DOT_PRODUCT(kgvec(1:3), &
545 : scoord(1:3, j) + &
546 6864 : REAL(agauge(1:3, j), KIND=dp))
547 : END DO
548 : EXIT
549 : END IF
550 : END DO
551 208 : CPASSERT(ic <= crys_sym%nrtot)
552 : END IF
553 : END DO
554 234 : DO is = 1, SIZE(crys_sym%kplink, 2)
555 234 : IF (crys_sym%kplink(2, is) == ik) THEN
556 10192 : DO ic = 1, crys_sym%nrtot
557 788736 : srotmat = MATMUL(cell%h_inv, MATMUL(crys_sym%rt(1:3, 1:3, ic), cell%hmat))
558 129792 : frot(1:3, 1:3) = NINT(srotmat(1:3, 1:3))
559 259584 : krot(1:3, 1:3) = NINT(TRANSPOSE(inv_3x3(REAL(frot(1:3, 1:3), KIND=dp))))
560 30160 : DO isign = 1, 2
561 19968 : ir = MERGE(crys_sym%ibrot(ic), -crys_sym%ibrot(ic), isign == 1)
562 19968 : IF (ir == crys_sym%kpop(is)) CYCLE
563 : kgvec(1:3) = crys_sym%kpmesh(1:3, is) - &
564 : MATMUL(REAL(MERGE(krot(1:3, 1:3), -krot(1:3, 1:3), &
565 : isign == 1), KIND=dp), &
566 553280 : kpoint%xkp(1:3, ik))
567 79040 : diff(1:3) = kgvec(1:3) - ANINT(kgvec(1:3))
568 46592 : IF (ALL(ABS(diff(1:3)) < eps_kpoint)) THEN
569 2288 : nr = nr + 1
570 2288 : kpsym%rotp(nr) = ir
571 29744 : kpsym%rot(1:3, 1:3, nr) = crys_sym%rt(1:3, 1:3, ic)
572 9152 : kpsym%xkp(1:3, nr) = crys_sym%kpmesh(1:3, is)
573 9152 : kgvec(1:3) = ANINT(kgvec(1:3))
574 20592 : kpsym%f0(1:natom, nr) = crys_sym%f0(1:natom, ic)
575 20592 : DO j = 1, natom
576 292864 : srot(1:3) = MATMUL(srotmat, scoord(1:3, j)) + crys_sym%vt(1:3, ic)
577 : kpsym%fcell(1:3, j, nr) = &
578 73216 : NINT(srot(1:3) - scoord(1:3, kpsym%f0(j, nr)))
579 : kpsym%fcell_gauge(1:3, j, nr) = &
580 : NINT(srot(1:3) - scoord(1:3, kpsym%f0(j, nr))) + &
581 : MATMUL(frot(1:3, 1:3), agauge(1:3, j)) - &
582 347776 : agauge(1:3, kpsym%f0(j, nr))
583 : kpsym%kgphase(j, nr) = DOT_PRODUCT(kgvec(1:3), &
584 : scoord(1:3, j) + &
585 75504 : REAL(agauge(1:3, j), KIND=dp))
586 : END DO
587 : END IF
588 : END DO
589 : END DO
590 : END IF
591 : END DO
592 26 : kpsym%nwred = nr
593 : END IF
594 : END IF
595 : END DO
596 234 : nkind = MAXVAL(atype)
597 26 : ns = crys_sym%nrtot
598 1378 : ALLOCATE (kpoint%kind_rotmat(ns, nkind))
599 1274 : DO i = 1, ns
600 2522 : DO j = 1, nkind
601 2496 : NULLIFY (kpoint%kind_rotmat(i, j)%rmat)
602 : END DO
603 : END DO
604 78 : ALLOCATE (kpoint%ibrot(ns))
605 1274 : kpoint%ibrot(1:ns) = crys_sym%ibrot(1:ns)
606 :
607 26 : CALL release_csym_type(crys_sym)
608 26 : DEALLOCATE (scoord, atype)
609 26 : DEALLOCATE (agauge)
610 : END IF
611 : CASE DEFAULT
612 11750 : CPABORT("Option invalid or unavailable for kpoint%kp_scheme")
613 : END SELECT
614 :
615 : ! check for consistency of options
616 11770 : SELECT CASE (kpoint%kp_scheme)
617 : CASE ("NONE")
618 : ! don't use k-point code
619 : CASE ("GAMMA")
620 20 : CPASSERT(kpoint%nkp == 1)
621 100 : CPASSERT(SUM(ABS(kpoint%xkp)) <= 1.e-12_dp)
622 20 : CPASSERT(kpoint%wkp(1) == 1.0_dp)
623 20 : CPASSERT(.NOT. kpoint%symmetry)
624 : CASE ("GENERAL")
625 36 : CPASSERT(kpoint%nkp >= 1)
626 : CASE ("MONKHORST-PACK", "MACDONALD")
627 11750 : CPASSERT(kpoint%nkp >= 1)
628 : END SELECT
629 11750 : IF (kpoint%use_real_wfn) THEN
630 : ! what about inversion symmetry?
631 40 : ikloop: DO ik = 1, kpoint%nkp
632 100 : DO i = 1, 3
633 60 : spez = (kpoint%xkp(i, ik) == 0.0_dp .OR. kpoint%xkp(i, ik) == 0.5_dp)
634 20 : IF (.NOT. spez) EXIT ikloop
635 : END DO
636 : END DO ikloop
637 20 : IF (.NOT. spez) THEN
638 : ! Warning: real wfn might be wrong for this system
639 : CALL cp_warn(__LOCATION__, &
640 : "A calculation using real wavefunctions is requested. "// &
641 0 : "We could not determine if the symmetry of the system allows real wavefunctions. ")
642 : END IF
643 : END IF
644 :
645 11750 : CALL timestop(handle)
646 :
647 23500 : END SUBROUTINE kpoint_initialize
648 :
649 : ! **************************************************************************************************
650 : !> \brief Initialize the kpoint environment
651 : !> \param kpoint Kpoint environment
652 : !> \param para_env ...
653 : !> \param blacs_env ...
654 : !> \param with_aux_fit ...
655 : ! **************************************************************************************************
656 3126 : SUBROUTINE kpoint_env_initialize(kpoint, para_env, blacs_env, with_aux_fit)
657 :
658 : TYPE(kpoint_type), INTENT(INOUT) :: kpoint
659 : TYPE(mp_para_env_type), INTENT(IN), TARGET :: para_env
660 : TYPE(cp_blacs_env_type), INTENT(IN), TARGET :: blacs_env
661 : LOGICAL, INTENT(IN), OPTIONAL :: with_aux_fit
662 :
663 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_env_initialize'
664 :
665 : INTEGER :: handle, igr, ik, ikk, ngr, niogrp, nkp, &
666 : nkp_grp, nkp_loc, npe, unit_nr
667 : INTEGER, DIMENSION(2) :: dims, pos
668 : LOGICAL :: aux_fit
669 3126 : TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_aux_env, kp_env
670 : TYPE(kpoint_env_type), POINTER :: kp
671 3126 : TYPE(mp_cart_type) :: comm_cart
672 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp, para_env_kp
673 :
674 3126 : CALL timeset(routineN, handle)
675 :
676 3126 : IF (PRESENT(with_aux_fit)) THEN
677 3018 : aux_fit = with_aux_fit
678 : ELSE
679 : aux_fit = .FALSE.
680 : END IF
681 :
682 3126 : kpoint%para_env => para_env
683 3126 : CALL kpoint%para_env%retain()
684 3126 : kpoint%blacs_env_all => blacs_env
685 3126 : CALL kpoint%blacs_env_all%retain()
686 :
687 3126 : CPASSERT(.NOT. ASSOCIATED(kpoint%kp_env))
688 3126 : IF (aux_fit) THEN
689 32 : CPASSERT(.NOT. ASSOCIATED(kpoint%kp_aux_env))
690 : END IF
691 :
692 3126 : NULLIFY (kp_env, kp_aux_env)
693 3126 : nkp = kpoint%nkp
694 3126 : npe = para_env%num_pe
695 3126 : IF (npe == 1) THEN
696 : ! only one process available -> owns all kpoints
697 0 : ALLOCATE (kp_env(nkp))
698 0 : DO ik = 1, nkp
699 0 : NULLIFY (kp_env(ik)%kpoint_env)
700 0 : CALL kpoint_env_create(kp_env(ik)%kpoint_env)
701 0 : kp => kp_env(ik)%kpoint_env
702 0 : kp%nkpoint = ik
703 0 : kp%wkp = kpoint%wkp(ik)
704 0 : kp%xkp(1:3) = kpoint%xkp(1:3, ik)
705 0 : kp%is_local = .TRUE.
706 : END DO
707 0 : kpoint%kp_env => kp_env
708 :
709 0 : IF (aux_fit) THEN
710 0 : ALLOCATE (kp_aux_env(nkp))
711 0 : DO ik = 1, nkp
712 0 : NULLIFY (kp_aux_env(ik)%kpoint_env)
713 0 : CALL kpoint_env_create(kp_aux_env(ik)%kpoint_env)
714 0 : kp => kp_aux_env(ik)%kpoint_env
715 0 : kp%nkpoint = ik
716 0 : kp%wkp = kpoint%wkp(ik)
717 0 : kp%xkp(1:3) = kpoint%xkp(1:3, ik)
718 0 : kp%is_local = .TRUE.
719 : END DO
720 :
721 0 : kpoint%kp_aux_env => kp_aux_env
722 : END IF
723 :
724 0 : ALLOCATE (kpoint%kp_dist(2, 1))
725 0 : kpoint%kp_dist(1, 1) = 1
726 0 : kpoint%kp_dist(2, 1) = nkp
727 0 : kpoint%kp_range(1) = 1
728 0 : kpoint%kp_range(2) = nkp
729 :
730 : ! parallel environments
731 0 : kpoint%para_env_kp => para_env
732 0 : CALL kpoint%para_env_kp%retain()
733 0 : kpoint%para_env_inter_kp => para_env
734 0 : CALL kpoint%para_env_inter_kp%retain()
735 0 : kpoint%iogrp = .TRUE.
736 0 : kpoint%nkp_groups = 1
737 : ELSE
738 3126 : IF (kpoint%parallel_group_size == -1) THEN
739 : ! maximum parallelization over kpoints
740 : ! making sure that the group size divides the npe and the nkp_grp the nkp
741 : ! in the worst case, there will be no parallelism over kpoints.
742 8106 : DO igr = npe, 1, -1
743 5404 : IF (MOD(npe, igr) /= 0) CYCLE
744 5404 : nkp_grp = npe/igr
745 5404 : IF (MOD(nkp, nkp_grp) /= 0) CYCLE
746 8106 : ngr = igr
747 : END DO
748 424 : ELSE IF (kpoint%parallel_group_size == 0) THEN
749 : ! no parallelization over kpoints
750 332 : ngr = npe
751 92 : ELSE IF (kpoint%parallel_group_size > 0) THEN
752 92 : ngr = MIN(kpoint%parallel_group_size, npe)
753 : ELSE
754 0 : CPABORT("kpoint%parallel_group_size cannot be smaller than -1")
755 : END IF
756 3126 : nkp_grp = npe/ngr
757 : ! processor dimensions
758 3126 : dims(1) = ngr
759 3126 : dims(2) = nkp_grp
760 3126 : CPASSERT(MOD(nkp, nkp_grp) == 0)
761 3126 : nkp_loc = nkp/nkp_grp
762 :
763 3126 : IF ((dims(1)*dims(2) /= npe)) THEN
764 0 : CPABORT("Number of processors is not divisible by the kpoint group size.")
765 : END IF
766 :
767 : ! Create the subgroups, one for each k-point group and one interconnecting group
768 3126 : CALL comm_cart%create(comm_old=para_env, ndims=2, dims=dims)
769 9378 : pos = comm_cart%mepos_cart
770 3126 : ALLOCATE (para_env_kp)
771 3126 : CALL para_env_kp%from_split(comm_cart, pos(2))
772 3126 : ALLOCATE (para_env_inter_kp)
773 3126 : CALL para_env_inter_kp%from_split(comm_cart, pos(1))
774 3126 : CALL comm_cart%free()
775 :
776 3126 : niogrp = 0
777 3126 : IF (para_env%is_source()) niogrp = 1
778 3126 : CALL para_env_kp%sum(niogrp)
779 3126 : kpoint%iogrp = (niogrp == 1)
780 :
781 : ! parallel groups
782 3126 : kpoint%para_env_kp => para_env_kp
783 3126 : kpoint%para_env_inter_kp => para_env_inter_kp
784 :
785 : ! distribution of kpoints
786 9378 : ALLOCATE (kpoint%kp_dist(2, nkp_grp))
787 8202 : DO igr = 1, nkp_grp
788 18354 : kpoint%kp_dist(1:2, igr) = get_limit(nkp, nkp_grp, igr - 1)
789 : END DO
790 : ! local kpoints
791 9378 : kpoint%kp_range(1:2) = kpoint%kp_dist(1:2, para_env_inter_kp%mepos + 1)
792 :
793 16422 : ALLOCATE (kp_env(nkp_loc))
794 10170 : DO ik = 1, nkp_loc
795 7044 : NULLIFY (kp_env(ik)%kpoint_env)
796 7044 : ikk = kpoint%kp_range(1) + ik - 1
797 7044 : CALL kpoint_env_create(kp_env(ik)%kpoint_env)
798 7044 : kp => kp_env(ik)%kpoint_env
799 7044 : kp%nkpoint = ikk
800 7044 : kp%wkp = kpoint%wkp(ikk)
801 28176 : kp%xkp(1:3) = kpoint%xkp(1:3, ikk)
802 10170 : kp%is_local = (ngr == 1)
803 : END DO
804 3126 : kpoint%kp_env => kp_env
805 :
806 3126 : IF (aux_fit) THEN
807 282 : ALLOCATE (kp_aux_env(nkp_loc))
808 250 : DO ik = 1, nkp_loc
809 218 : NULLIFY (kp_aux_env(ik)%kpoint_env)
810 218 : ikk = kpoint%kp_range(1) + ik - 1
811 218 : CALL kpoint_env_create(kp_aux_env(ik)%kpoint_env)
812 218 : kp => kp_aux_env(ik)%kpoint_env
813 218 : kp%nkpoint = ikk
814 218 : kp%wkp = kpoint%wkp(ikk)
815 872 : kp%xkp(1:3) = kpoint%xkp(1:3, ikk)
816 250 : kp%is_local = (ngr == 1)
817 : END DO
818 32 : kpoint%kp_aux_env => kp_aux_env
819 : END IF
820 :
821 3126 : unit_nr = cp_logger_get_default_io_unit()
822 :
823 3126 : IF (unit_nr > 0 .AND. kpoint%verbose) THEN
824 1179 : WRITE (unit_nr, *)
825 1179 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Number of kpoint groups ", nkp_grp
826 1179 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Size of each kpoint group", ngr
827 1179 : WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Number of kpoints per group", nkp_loc
828 : END IF
829 3126 : kpoint%nkp_groups = nkp_grp
830 :
831 : END IF
832 :
833 3126 : CALL timestop(handle)
834 :
835 6252 : END SUBROUTINE kpoint_env_initialize
836 :
837 : ! **************************************************************************************************
838 : !> \brief Initialize a set of MOs and density matrix for each kpoint (kpoint group)
839 : !> \param kpoint Kpoint environment
840 : !> \param mos Reference MOs (global)
841 : !> \param added_mos ...
842 : !> \param for_aux_fit ...
843 : ! **************************************************************************************************
844 3164 : SUBROUTINE kpoint_initialize_mos(kpoint, mos, added_mos, for_aux_fit)
845 :
846 : TYPE(kpoint_type), POINTER :: kpoint
847 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
848 : INTEGER, INTENT(IN), OPTIONAL :: added_mos
849 : LOGICAL, OPTIONAL :: for_aux_fit
850 :
851 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_initialize_mos'
852 :
853 : INTEGER :: handle, ic, ik, is, nadd, nao, nc, &
854 : nelectron, nkp_loc, nmo, nmorig(2), &
855 : nspin
856 : LOGICAL :: aux_fit
857 : REAL(KIND=dp) :: flexible_electron_count, maxocc, n_el_f
858 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
859 3164 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools
860 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
861 : TYPE(cp_fm_type), POINTER :: fmlocal
862 : TYPE(kpoint_env_type), POINTER :: kp
863 : TYPE(qs_matrix_pools_type), POINTER :: mpools
864 :
865 3164 : CALL timeset(routineN, handle)
866 :
867 3164 : IF (PRESENT(for_aux_fit)) THEN
868 32 : aux_fit = for_aux_fit
869 : ELSE
870 : aux_fit = .FALSE.
871 : END IF
872 :
873 3164 : CPASSERT(ASSOCIATED(kpoint))
874 :
875 : IF (.TRUE. .OR. ASSOCIATED(mos(1)%mo_coeff)) THEN
876 3164 : IF (aux_fit) THEN
877 32 : CPASSERT(ASSOCIATED(kpoint%kp_aux_env))
878 : END IF
879 :
880 3164 : IF (PRESENT(added_mos)) THEN
881 90 : nadd = added_mos
882 : ELSE
883 : nadd = 0
884 : END IF
885 :
886 3164 : IF (kpoint%use_real_wfn) THEN
887 : nc = 1
888 : ELSE
889 3146 : nc = 2
890 : END IF
891 3164 : nspin = SIZE(mos, 1)
892 3164 : nkp_loc = kpoint%kp_range(2) - kpoint%kp_range(1) + 1
893 3164 : IF (nkp_loc > 0) THEN
894 3164 : IF (aux_fit) THEN
895 32 : CPASSERT(SIZE(kpoint%kp_aux_env) == nkp_loc)
896 : ELSE
897 3132 : CPASSERT(SIZE(kpoint%kp_env) == nkp_loc)
898 : END IF
899 : ! allocate the mo sets, correct number of kpoints (local), real/complex, spin
900 10434 : DO ik = 1, nkp_loc
901 7270 : IF (aux_fit) THEN
902 218 : kp => kpoint%kp_aux_env(ik)%kpoint_env
903 : ELSE
904 7052 : kp => kpoint%kp_env(ik)%kpoint_env
905 : END IF
906 52910 : ALLOCATE (kp%mos(nc, nspin))
907 15220 : DO is = 1, nspin
908 : CALL get_mo_set(mos(is), nao=nao, nmo=nmo, nelectron=nelectron, &
909 7950 : n_el_f=n_el_f, maxocc=maxocc, flexible_electron_count=flexible_electron_count)
910 7950 : nmo = MIN(nao, nmo + nadd)
911 31100 : DO ic = 1, nc
912 : CALL allocate_mo_set(kp%mos(ic, is), nao, nmo, nelectron, n_el_f, maxocc, &
913 23830 : flexible_electron_count)
914 : END DO
915 : END DO
916 : ! freshly allocated MOS carry no coefficients: solvers must
917 : ! cold-start again until an extrapolation refills them
918 10434 : kp%mos_prefilled = .FALSE.
919 : END DO
920 :
921 : ! generate the blacs environment for the kpoint group
922 : ! we generate a blacs env for each kpoint group in parallel
923 : ! we assume here that the group para_env_inter_kp will connect
924 : ! equivalent parts of fm matrices, i.e. no reshuffeling of processors
925 3164 : NULLIFY (blacs_env)
926 3164 : IF (ASSOCIATED(kpoint%blacs_env)) THEN
927 38 : blacs_env => kpoint%blacs_env
928 : ELSE
929 3126 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=kpoint%para_env_kp)
930 3126 : kpoint%blacs_env => blacs_env
931 : END IF
932 :
933 : ! set possible new number of MOs
934 6530 : DO is = 1, nspin
935 3366 : CALL get_mo_set(mos(is), nmo=nmorig(is))
936 3366 : nmo = MIN(nao, nmorig(is) + nadd)
937 6530 : CALL set_mo_set(mos(is), nmo=nmo)
938 : END DO
939 : ! matrix pools for the kpoint group, information on MOs is transferred using
940 : ! generic mos structure
941 3164 : NULLIFY (mpools)
942 3164 : CALL mpools_create(mpools=mpools)
943 : CALL mpools_rebuild_fm_pools(mpools=mpools, mos=mos, &
944 3164 : blacs_env=blacs_env, para_env=kpoint%para_env_kp)
945 :
946 3164 : IF (aux_fit) THEN
947 32 : kpoint%mpools_aux_fit => mpools
948 : ELSE
949 3132 : kpoint%mpools => mpools
950 : END IF
951 :
952 : ! reset old number of MOs
953 6530 : DO is = 1, nspin
954 6530 : CALL set_mo_set(mos(is), nmo=nmorig(is))
955 : END DO
956 :
957 : ! allocate density matrices
958 3164 : CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
959 3164 : ALLOCATE (fmlocal)
960 3164 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
961 3164 : CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
962 10434 : DO ik = 1, nkp_loc
963 7270 : IF (aux_fit) THEN
964 218 : kp => kpoint%kp_aux_env(ik)%kpoint_env
965 : ELSE
966 7052 : kp => kpoint%kp_env(ik)%kpoint_env
967 : END IF
968 : ! density matrix
969 7270 : CALL cp_fm_release(kp%pmat)
970 52910 : ALLOCATE (kp%pmat(nc, nspin))
971 15220 : DO is = 1, nspin
972 31100 : DO ic = 1, nc
973 23830 : CALL cp_fm_create(kp%pmat(ic, is), matrix_struct)
974 : END DO
975 : END DO
976 : ! energy weighted density matrix
977 7270 : CALL cp_fm_release(kp%wmat)
978 45640 : ALLOCATE (kp%wmat(nc, nspin))
979 18384 : DO is = 1, nspin
980 31100 : DO ic = 1, nc
981 23830 : CALL cp_fm_create(kp%wmat(ic, is), matrix_struct)
982 : END DO
983 : END DO
984 : END DO
985 3164 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
986 3164 : DEALLOCATE (fmlocal)
987 :
988 : END IF
989 :
990 : END IF
991 :
992 3164 : CALL timestop(handle)
993 :
994 3164 : END SUBROUTINE kpoint_initialize_mos
995 :
996 : ! **************************************************************************************************
997 : !> \brief ...
998 : !> \param kpoint ...
999 : ! **************************************************************************************************
1000 114 : SUBROUTINE kpoint_initialize_mo_set(kpoint)
1001 : TYPE(kpoint_type), POINTER :: kpoint
1002 :
1003 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_initialize_mo_set'
1004 :
1005 : INTEGER :: handle, ic, ik, ikk, ispin
1006 114 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
1007 : TYPE(cp_fm_type), POINTER :: mo_coeff
1008 114 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: moskp
1009 :
1010 114 : CALL timeset(routineN, handle)
1011 :
1012 1008 : DO ik = 1, SIZE(kpoint%kp_env)
1013 894 : CALL mpools_get(kpoint%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
1014 894 : moskp => kpoint%kp_env(ik)%kpoint_env%mos
1015 894 : ikk = kpoint%kp_range(1) + ik - 1
1016 894 : CPASSERT(ASSOCIATED(moskp))
1017 1940 : DO ispin = 1, SIZE(moskp, 2)
1018 3690 : DO ic = 1, SIZE(moskp, 1)
1019 1864 : CALL get_mo_set(moskp(ic, ispin), mo_coeff=mo_coeff)
1020 2796 : IF (.NOT. ASSOCIATED(mo_coeff)) THEN
1021 : CALL init_mo_set(moskp(ic, ispin), &
1022 1864 : fm_pool=ao_mo_fm_pools(ispin)%pool, name="kpoints")
1023 : END IF
1024 : END DO
1025 : END DO
1026 : END DO
1027 :
1028 114 : CALL timestop(handle)
1029 :
1030 114 : END SUBROUTINE kpoint_initialize_mo_set
1031 :
1032 : ! **************************************************************************************************
1033 : !> \brief Generates the mapping of cell indices and linear RS index
1034 : !> CELL (0,0,0) is always mapped to index 1
1035 : !> \param kpoint Kpoint environment
1036 : !> \param sab_nl Defining neighbour list
1037 : !> \param para_env Parallel environment
1038 : !> \param nimages [output]
1039 : ! **************************************************************************************************
1040 3898 : SUBROUTINE kpoint_init_cell_index(kpoint, sab_nl, para_env, nimages)
1041 :
1042 : TYPE(kpoint_type), POINTER :: kpoint
1043 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1044 : POINTER :: sab_nl
1045 : TYPE(mp_para_env_type), POINTER :: para_env
1046 : INTEGER, INTENT(OUT) :: nimages
1047 :
1048 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_init_cell_index'
1049 :
1050 : INTEGER :: handle, i1, i2, i3, ic, icount, it, &
1051 : ncount
1052 : INTEGER, DIMENSION(3) :: cell, itm
1053 3898 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell, list
1054 3898 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index, cti
1055 : LOGICAL :: new
1056 : TYPE(neighbor_list_iterator_p_type), &
1057 3898 : DIMENSION(:), POINTER :: nl_iterator
1058 :
1059 3898 : NULLIFY (cell_to_index, index_to_cell)
1060 :
1061 3898 : CALL timeset(routineN, handle)
1062 :
1063 3898 : CPASSERT(ASSOCIATED(kpoint))
1064 :
1065 3898 : ALLOCATE (list(3, 125))
1066 1952898 : list = 0
1067 3898 : icount = 1
1068 :
1069 3898 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
1070 1326902 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1071 1323004 : CALL get_iterator_info(nl_iterator, cell=cell)
1072 :
1073 1323004 : new = .TRUE.
1074 78922169 : DO ic = 1, icount
1075 78719072 : IF (cell(1) == list(1, ic) .AND. cell(2) == list(2, ic) .AND. &
1076 203097 : cell(3) == list(3, ic)) THEN
1077 : new = .FALSE.
1078 : EXIT
1079 : END IF
1080 : END DO
1081 1326902 : IF (new) THEN
1082 203097 : icount = icount + 1
1083 203097 : IF (icount > SIZE(list, 2)) THEN
1084 607 : CALL reallocate(list, 1, 3, 1, 2*SIZE(list, 2))
1085 : END IF
1086 812388 : list(1:3, icount) = cell(1:3)
1087 : END IF
1088 :
1089 : END DO
1090 3898 : CALL neighbor_list_iterator_release(nl_iterator)
1091 :
1092 210893 : itm(1) = MAXVAL(ABS(list(1, 1:icount)))
1093 210893 : itm(2) = MAXVAL(ABS(list(2, 1:icount)))
1094 210893 : itm(3) = MAXVAL(ABS(list(3, 1:icount)))
1095 3898 : CALL para_env%max(itm)
1096 15592 : it = MAXVAL(itm(1:3))
1097 3898 : IF (ASSOCIATED(kpoint%cell_to_index)) THEN
1098 3894 : DEALLOCATE (kpoint%cell_to_index)
1099 : END IF
1100 19490 : ALLOCATE (kpoint%cell_to_index(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
1101 3898 : cell_to_index => kpoint%cell_to_index
1102 3898 : cti => cell_to_index
1103 570260 : cti(:, :, :) = 0
1104 210893 : DO ic = 1, icount
1105 206995 : i1 = list(1, ic)
1106 206995 : i2 = list(2, ic)
1107 206995 : i3 = list(3, ic)
1108 210893 : cti(i1, i2, i3) = ic
1109 : END DO
1110 1136622 : CALL para_env%sum(cti)
1111 3898 : ncount = 0
1112 20528 : DO i1 = -itm(1), itm(1)
1113 106530 : DO i2 = -itm(2), itm(2)
1114 572290 : DO i3 = -itm(3), itm(3)
1115 555660 : IF (cti(i1, i2, i3) == 0) THEN
1116 221234 : cti(i1, i2, i3) = 1000000
1117 : ELSE
1118 248424 : ncount = ncount + 1
1119 248424 : cti(i1, i2, i3) = (ABS(i1) + ABS(i2) + ABS(i3))*1000 + ABS(i3)*100 + ABS(i2)*10 + ABS(i1)
1120 248424 : cti(i1, i2, i3) = cti(i1, i2, i3) + (i1 + i2 + i3)
1121 : END IF
1122 : END DO
1123 : END DO
1124 : END DO
1125 :
1126 3898 : IF (ASSOCIATED(kpoint%index_to_cell)) THEN
1127 3898 : DEALLOCATE (kpoint%index_to_cell)
1128 : END IF
1129 11694 : ALLOCATE (kpoint%index_to_cell(3, ncount))
1130 3898 : index_to_cell => kpoint%index_to_cell
1131 252322 : DO ic = 1, ncount
1132 160754864 : cell = MINLOC(cti)
1133 248424 : i1 = cell(1) - 1 - itm(1)
1134 248424 : i2 = cell(2) - 1 - itm(2)
1135 248424 : i3 = cell(3) - 1 - itm(3)
1136 248424 : cti(i1, i2, i3) = 1000000
1137 248424 : index_to_cell(1, ic) = i1
1138 248424 : index_to_cell(2, ic) = i2
1139 252322 : index_to_cell(3, ic) = i3
1140 : END DO
1141 570260 : cti(:, :, :) = 0
1142 252322 : DO ic = 1, ncount
1143 248424 : i1 = index_to_cell(1, ic)
1144 248424 : i2 = index_to_cell(2, ic)
1145 248424 : i3 = index_to_cell(3, ic)
1146 252322 : cti(i1, i2, i3) = ic
1147 : END DO
1148 :
1149 : ! keep pointer to this neighborlist
1150 3898 : kpoint%sab_nl => sab_nl
1151 :
1152 : ! set number of images
1153 3898 : nimages = SIZE(index_to_cell, 2)
1154 :
1155 3898 : DEALLOCATE (list)
1156 :
1157 3898 : CALL timestop(handle)
1158 :
1159 3898 : END SUBROUTINE kpoint_init_cell_index
1160 :
1161 : ! **************************************************************************************************
1162 : !> \brief Transformation of real space matrices to a kpoint
1163 : !> \param rmatrix Real part of kpoint matrix
1164 : !> \param cmatrix Complex part of kpoint matrix (optional)
1165 : !> \param rsmat Real space matrices
1166 : !> \param ispin Spin index
1167 : !> \param xkp Kpoint coordinates
1168 : !> \param cell_to_index mapping of cell indices to RS index
1169 : !> \param sab_nl Defining neighbor list
1170 : !> \param is_complex Matrix to be transformed is imaginary
1171 : !> \param rs_sign Matrix to be transformed is csaled by rs_sign
1172 : ! **************************************************************************************************
1173 586752 : SUBROUTINE rskp_transform(rmatrix, cmatrix, rsmat, ispin, &
1174 : xkp, cell_to_index, sab_nl, is_complex, rs_sign)
1175 :
1176 : TYPE(dbcsr_type) :: rmatrix
1177 : TYPE(dbcsr_type), OPTIONAL :: cmatrix
1178 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
1179 : INTEGER, INTENT(IN) :: ispin
1180 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
1181 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1182 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1183 : POINTER :: sab_nl
1184 : LOGICAL, INTENT(IN), OPTIONAL :: is_complex
1185 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: rs_sign
1186 :
1187 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rskp_transform'
1188 :
1189 : INTEGER :: handle, iatom, ic, icol, irow, jatom, &
1190 : nimg
1191 : INTEGER, DIMENSION(3) :: cell
1192 : LOGICAL :: do_symmetric, found, my_complex, &
1193 : wfn_real_only
1194 : REAL(KIND=dp) :: arg, coskl, fsign, fsym, sinkl
1195 293376 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, rblock, rsblock
1196 : TYPE(neighbor_list_iterator_p_type), &
1197 293376 : DIMENSION(:), POINTER :: nl_iterator
1198 :
1199 293376 : CALL timeset(routineN, handle)
1200 :
1201 293376 : my_complex = .FALSE.
1202 293376 : IF (PRESENT(is_complex)) my_complex = is_complex
1203 :
1204 293376 : fsign = 1.0_dp
1205 293376 : IF (PRESENT(rs_sign)) fsign = rs_sign
1206 :
1207 293376 : wfn_real_only = .TRUE.
1208 293376 : IF (PRESENT(cmatrix)) wfn_real_only = .FALSE.
1209 :
1210 293376 : nimg = SIZE(rsmat, 2)
1211 :
1212 293376 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1213 :
1214 293376 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
1215 104270755 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1216 103977379 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
1217 :
1218 : ! fsym = +- 1 is due to real space matrices being non-symmetric (although in a symmtric type)
1219 : ! with the link S_mu^0,nu^b = S_nu^0,mu^-b, and the KP matrices beeing Hermitian
1220 103977379 : fsym = 1.0_dp
1221 103977379 : irow = iatom
1222 103977379 : icol = jatom
1223 103977379 : IF (do_symmetric .AND. (iatom > jatom)) THEN
1224 44285608 : irow = jatom
1225 44285608 : icol = iatom
1226 44285608 : fsym = -1.0_dp
1227 : END IF
1228 :
1229 103977379 : ic = cell_to_index(cell(1), cell(2), cell(3))
1230 103977379 : IF (ic < 1 .OR. ic > nimg) CYCLE
1231 :
1232 103976635 : arg = REAL(cell(1), dp)*xkp(1) + REAL(cell(2), dp)*xkp(2) + REAL(cell(3), dp)*xkp(3)
1233 103976635 : IF (my_complex) THEN
1234 3466896 : coskl = fsign*fsym*COS(twopi*arg)
1235 3466896 : sinkl = fsign*SIN(twopi*arg)
1236 : ELSE
1237 100509739 : coskl = fsign*COS(twopi*arg)
1238 100509739 : sinkl = fsign*fsym*SIN(twopi*arg)
1239 : END IF
1240 :
1241 : CALL dbcsr_get_block_p(matrix=rsmat(ispin, ic)%matrix, row=irow, col=icol, &
1242 103976635 : block=rsblock, found=found)
1243 103976635 : IF (.NOT. found) CYCLE
1244 :
1245 104270011 : IF (wfn_real_only) THEN
1246 : CALL dbcsr_get_block_p(matrix=rmatrix, row=irow, col=icol, &
1247 529850 : block=rblock, found=found)
1248 529850 : IF (.NOT. found) CYCLE
1249 249444630 : rblock = rblock + coskl*rsblock
1250 : ELSE
1251 : CALL dbcsr_get_block_p(matrix=rmatrix, row=irow, col=icol, &
1252 103446785 : block=rblock, found=found)
1253 103446785 : IF (.NOT. found) CYCLE
1254 : CALL dbcsr_get_block_p(matrix=cmatrix, row=irow, col=icol, &
1255 103446785 : block=cblock, found=found)
1256 103446785 : IF (.NOT. found) CYCLE
1257 13238731357 : rblock = rblock + coskl*rsblock
1258 13238731357 : cblock = cblock + sinkl*rsblock
1259 : END IF
1260 :
1261 : END DO
1262 293376 : CALL neighbor_list_iterator_release(nl_iterator)
1263 :
1264 293376 : CALL timestop(handle)
1265 :
1266 293376 : END SUBROUTINE rskp_transform
1267 :
1268 : ! **************************************************************************************************
1269 : !> \brief Prepare a batched real-cell to complete reciprocal-grid transform for local DBCSR blocks.
1270 : !> The stored blocks remain MPI local. Irregular or too-large grids return used_fft=.FALSE.
1271 : !> \param grid cached transformed local blocks and their DBCSR coordinates
1272 : !> \param rmatrix allocated output template defining the local block layout
1273 : !> \param rsmat real-space matrix set
1274 : !> \param ispin spin component of rsmat
1275 : !> \param xkp complete reciprocal grid in arbitrary order
1276 : !> \param nkp_grid reciprocal grid dimensions
1277 : !> \param cell_to_index real-cell coordinate mapping
1278 : !> \param sab_nl neighbor list defining the stored block orientation
1279 : !> \param used_fft whether the cache was prepared
1280 : !> \param is_complex whether the real-space operator is imaginary
1281 : !> \param rs_sign optional overall sign
1282 : !> \param max_storage_bytes optional conservative per-rank memory limit
1283 : ! **************************************************************************************************
1284 216 : SUBROUTINE rskp_transform_grid_prepare(grid, rmatrix, rsmat, ispin, xkp, nkp_grid, &
1285 : cell_to_index, sab_nl, used_fft, is_complex, rs_sign, &
1286 : max_storage_bytes)
1287 :
1288 : TYPE(rskp_grid_type), INTENT(INOUT) :: grid
1289 : TYPE(dbcsr_type), INTENT(INOUT) :: rmatrix
1290 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
1291 : INTEGER, INTENT(IN) :: ispin
1292 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp
1293 : INTEGER, DIMENSION(3), INTENT(IN) :: nkp_grid
1294 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1295 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1296 : POINTER :: sab_nl
1297 : LOGICAL, INTENT(OUT) :: used_fft
1298 : LOGICAL, INTENT(IN), OPTIONAL :: is_complex
1299 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: rs_sign
1300 : INTEGER(KIND=int_8), INTENT(IN), OPTIONAL :: max_storage_bytes
1301 :
1302 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rskp_transform_grid_prepare'
1303 :
1304 : INTEGER :: handle, i1, i2, i3, iatom, iblock, ic, &
1305 : icol, irow, jatom, nblkcols, nblkrows, &
1306 : nblocks, ncell, nelem, nimg, nkp, &
1307 : nvalues
1308 : INTEGER(KIND=int_8) :: memory_limit, storage_bytes
1309 216 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: block_map, index_to_cell
1310 216 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: fft_cell_to_index
1311 : INTEGER, DIMENSION(3) :: cell, fft_cell
1312 : LOGICAL :: do_symmetric, found, my_complex
1313 : REAL(KIND=dp) :: fsign, fsym, value_sign
1314 216 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: values_rs
1315 216 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rblock, rsblock
1316 : TYPE(dbcsr_iterator_type) :: iter
1317 : TYPE(neighbor_list_iterator_p_type), &
1318 216 : DIMENSION(:), POINTER :: nl_iterator
1319 :
1320 216 : CALL timeset(routineN, handle)
1321 216 : CALL rskp_transform_grid_release(grid)
1322 216 : used_fft = .FALSE.
1323 :
1324 216 : nkp = SIZE(xkp, 2)
1325 216 : nimg = SIZE(rsmat, 2)
1326 864 : IF (nkp < 2 .OR. PRODUCT(nkp_grid) /= nkp) THEN
1327 208 : CALL timestop(handle)
1328 208 : RETURN
1329 : END IF
1330 8 : IF (.NOT. regular_kpoint_grid(xkp, nkp_grid)) THEN
1331 0 : CALL timestop(handle)
1332 0 : RETURN
1333 : END IF
1334 :
1335 8 : nblocks = 0
1336 8 : nvalues = 0
1337 8 : CALL dbcsr_iterator_start(iter, rmatrix)
1338 152 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1339 144 : CALL dbcsr_iterator_next_block(iter, irow, icol, rblock)
1340 144 : nblocks = nblocks + 1
1341 440 : nvalues = nvalues + SIZE(rblock)
1342 : END DO
1343 8 : CALL dbcsr_iterator_stop(iter)
1344 8 : IF (nblocks == 0 .OR. nvalues == 0) THEN
1345 0 : CALL timestop(handle)
1346 0 : RETURN
1347 : END IF
1348 :
1349 8 : CALL dbcsr_get_info(rmatrix, nblkrows_total=nblkrows, nblkcols_total=nblkcols)
1350 32 : ALLOCATE (block_map(nblkrows, nblkcols), source=0)
1351 0 : ALLOCATE (grid%block_row(nblocks), grid%block_col(nblocks), &
1352 48 : grid%block_offset(nblocks), grid%block_nelem(nblocks))
1353 8 : iblock = 0
1354 8 : nelem = 0
1355 8 : CALL dbcsr_iterator_start(iter, rmatrix)
1356 152 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1357 144 : CALL dbcsr_iterator_next_block(iter, irow, icol, rblock)
1358 144 : iblock = iblock + 1
1359 144 : grid%block_row(iblock) = irow
1360 144 : grid%block_col(iblock) = icol
1361 144 : grid%block_offset(iblock) = nelem + 1
1362 432 : grid%block_nelem(iblock) = SIZE(rblock)
1363 144 : block_map(irow, icol) = iblock
1364 440 : nelem = nelem + SIZE(rblock)
1365 : END DO
1366 8 : CALL dbcsr_iterator_stop(iter)
1367 :
1368 : ALLOCATE (fft_cell_to_index(LBOUND(cell_to_index, 1):UBOUND(cell_to_index, 1), &
1369 : LBOUND(cell_to_index, 2):UBOUND(cell_to_index, 2), &
1370 40 : LBOUND(cell_to_index, 3):UBOUND(cell_to_index, 3)), source=0)
1371 8 : ncell = 0
1372 80 : DO i3 = LBOUND(cell_to_index, 3), UBOUND(cell_to_index, 3)
1373 568 : DO i2 = LBOUND(cell_to_index, 2), UBOUND(cell_to_index, 2)
1374 3976 : DO i1 = LBOUND(cell_to_index, 1), UBOUND(cell_to_index, 1)
1375 2744 : IF (cell_to_index(i1, i2, i3) < 1 .OR. cell_to_index(i1, i2, i3) > nimg) CYCLE
1376 1360 : IF (fft_cell_to_index(i1, i2, i3) == 0) THEN
1377 720 : ncell = ncell + 1
1378 720 : fft_cell_to_index(i1, i2, i3) = ncell
1379 : END IF
1380 1752 : IF (fft_cell_to_index(-i1, -i2, -i3) == 0) THEN
1381 712 : ncell = ncell + 1
1382 712 : fft_cell_to_index(-i1, -i2, -i3) = ncell
1383 : END IF
1384 : END DO
1385 : END DO
1386 : END DO
1387 8 : IF (ncell == 0) THEN
1388 0 : CALL rskp_transform_grid_release(grid)
1389 0 : CALL timestop(handle)
1390 0 : RETURN
1391 : END IF
1392 :
1393 8 : memory_limit = 512_int_8*1024_int_8**2
1394 8 : IF (PRESENT(max_storage_bytes)) memory_limit = max_storage_bytes
1395 : ! Include the real-cell input, complex result and conservative FFT work/padding allowance.
1396 : storage_bytes = INT(nvalues, int_8)*(8_int_8*INT(ncell, int_8) + &
1397 8 : 80_int_8*INT(nkp, int_8))
1398 8 : IF (storage_bytes > memory_limit) THEN
1399 0 : CALL rskp_transform_grid_release(grid)
1400 0 : CALL timestop(handle)
1401 0 : RETURN
1402 : END IF
1403 :
1404 24 : ALLOCATE (index_to_cell(3, ncell), source=0)
1405 80 : DO i3 = LBOUND(fft_cell_to_index, 3), UBOUND(fft_cell_to_index, 3)
1406 568 : DO i2 = LBOUND(fft_cell_to_index, 2), UBOUND(fft_cell_to_index, 2)
1407 3976 : DO i1 = LBOUND(fft_cell_to_index, 1), UBOUND(fft_cell_to_index, 1)
1408 2744 : ic = fft_cell_to_index(i1, i2, i3)
1409 7432 : IF (ic > 0) index_to_cell(:, ic) = [i1, i2, i3]
1410 : END DO
1411 : END DO
1412 : END DO
1413 :
1414 32 : ALLOCATE (values_rs(nvalues, 1, ncell), source=0.0_dp)
1415 32 : ALLOCATE (grid%values(nvalues, 1, nkp))
1416 8 : my_complex = .FALSE.
1417 8 : IF (PRESENT(is_complex)) my_complex = is_complex
1418 8 : fsign = 1.0_dp
1419 8 : IF (PRESENT(rs_sign)) fsign = rs_sign
1420 :
1421 8 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1422 8 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
1423 11432 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1424 11424 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
1425 :
1426 11424 : fsym = 1.0_dp
1427 11424 : irow = iatom
1428 11424 : icol = jatom
1429 11424 : IF (do_symmetric .AND. iatom > jatom) THEN
1430 5056 : irow = jatom
1431 5056 : icol = iatom
1432 5056 : fsym = -1.0_dp
1433 : END IF
1434 11424 : IF (irow < 1 .OR. irow > nblkrows .OR. icol < 1 .OR. icol > nblkcols) CYCLE
1435 11424 : iblock = block_map(irow, icol)
1436 11424 : IF (iblock == 0) CYCLE
1437 :
1438 11424 : ic = cell_to_index(cell(1), cell(2), cell(3))
1439 11424 : IF (ic < 1 .OR. ic > nimg) CYCLE
1440 : CALL dbcsr_get_block_p(matrix=rsmat(ispin, ic)%matrix, row=irow, col=icol, &
1441 11424 : block=rsblock, found=found)
1442 11424 : IF (.NOT. found) CYCLE
1443 :
1444 11424 : fft_cell = cell
1445 11424 : value_sign = fsign
1446 26592 : IF (fsym < 0.0_dp) fft_cell = -cell
1447 11424 : IF (my_complex) value_sign = value_sign*fsym
1448 11424 : ic = fft_cell_to_index(fft_cell(1), fft_cell(2), fft_cell(3))
1449 11424 : nelem = grid%block_nelem(iblock)
1450 11424 : i1 = grid%block_offset(iblock)
1451 : values_rs(i1:i1 + nelem - 1, 1, ic) = values_rs(i1:i1 + nelem - 1, 1, ic) + &
1452 308456 : value_sign*RESHAPE(rsblock, [nelem])
1453 : END DO
1454 8 : CALL neighbor_list_iterator_release(nl_iterator)
1455 :
1456 8 : CALL cell_to_k_grid_fft(values_rs, index_to_cell, xkp, nkp_grid, grid%values, used_fft)
1457 8 : IF (used_fft) THEN
1458 8 : grid%ready = .TRUE.
1459 : ELSE
1460 0 : CALL rskp_transform_grid_release(grid)
1461 : END IF
1462 :
1463 8 : CALL timestop(handle)
1464 :
1465 464 : END SUBROUTINE rskp_transform_grid_prepare
1466 :
1467 : ! **************************************************************************************************
1468 : !> \brief Extract one reciprocal-grid matrix from a prepared local DBCSR block cache.
1469 : !> \param grid ...
1470 : !> \param ikp ...
1471 : !> \param rmatrix ...
1472 : !> \param cmatrix ...
1473 : ! **************************************************************************************************
1474 512 : SUBROUTINE rskp_transform_grid_extract(grid, ikp, rmatrix, cmatrix)
1475 :
1476 : TYPE(rskp_grid_type), INTENT(IN) :: grid
1477 : INTEGER, INTENT(IN) :: ikp
1478 : TYPE(dbcsr_type), INTENT(INOUT) :: rmatrix
1479 : TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL :: cmatrix
1480 :
1481 : INTEGER :: iblock, ioffset, nelem
1482 : LOGICAL :: found
1483 512 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, rblock
1484 :
1485 0 : CPASSERT(grid%ready)
1486 512 : CPASSERT(ikp >= 1 .AND. ikp <= SIZE(grid%values, 3))
1487 512 : CALL dbcsr_set(rmatrix, 0.0_dp)
1488 512 : IF (PRESENT(cmatrix)) CALL dbcsr_set(cmatrix, 0.0_dp)
1489 :
1490 9728 : DO iblock = 1, SIZE(grid%block_row)
1491 9216 : ioffset = grid%block_offset(iblock)
1492 9216 : nelem = grid%block_nelem(iblock)
1493 : CALL dbcsr_get_block_p(rmatrix, grid%block_row(iblock), grid%block_col(iblock), &
1494 9216 : rblock, found=found)
1495 9216 : CPASSERT(found)
1496 : rblock = RESHAPE(REAL(grid%values(ioffset:ioffset + nelem - 1, 1, ikp), KIND=dp), &
1497 258048 : SHAPE(rblock))
1498 18944 : IF (PRESENT(cmatrix)) THEN
1499 : CALL dbcsr_get_block_p(cmatrix, grid%block_row(iblock), grid%block_col(iblock), &
1500 9216 : cblock, found=found)
1501 9216 : CPASSERT(found)
1502 258048 : cblock = RESHAPE(AIMAG(grid%values(ioffset:ioffset + nelem - 1, 1, ikp)), SHAPE(cblock))
1503 : END IF
1504 : END DO
1505 :
1506 512 : END SUBROUTINE rskp_transform_grid_extract
1507 :
1508 : ! **************************************************************************************************
1509 : !> \brief Release a reciprocal-grid transformation cache.
1510 : !> \param grid ...
1511 : ! **************************************************************************************************
1512 141072 : SUBROUTINE rskp_transform_grid_release(grid)
1513 :
1514 : TYPE(rskp_grid_type), INTENT(INOUT) :: grid
1515 :
1516 141072 : IF (ALLOCATED(grid%values)) DEALLOCATE (grid%values)
1517 141072 : IF (ALLOCATED(grid%block_row)) DEALLOCATE (grid%block_row)
1518 141072 : IF (ALLOCATED(grid%block_col)) DEALLOCATE (grid%block_col)
1519 141072 : IF (ALLOCATED(grid%block_offset)) DEALLOCATE (grid%block_offset)
1520 141072 : IF (ALLOCATED(grid%block_nelem)) DEALLOCATE (grid%block_nelem)
1521 141072 : grid%ready = .FALSE.
1522 :
1523 141072 : END SUBROUTINE rskp_transform_grid_release
1524 :
1525 : ! **************************************************************************************************
1526 : !> \brief Given the eigenvalues of all kpoints, calculates the occupation numbers
1527 : !> \param kpoint Kpoint environment
1528 : !> \param smear Smearing information
1529 : !> \param probe ...
1530 : !> \param added_mos_auto ...
1531 : !> \param added_mos_auto_grow ...
1532 : !> \param separate_spin_occupations enforce the electron count of each spin channel
1533 : ! **************************************************************************************************
1534 36112 : SUBROUTINE kpoint_set_mo_occupation( &
1535 36112 : kpoint, smear, probe, added_mos_auto, added_mos_auto_grow, separate_spin_occupations)
1536 :
1537 : TYPE(kpoint_type), POINTER :: kpoint
1538 : TYPE(smear_type) :: smear
1539 : TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
1540 : POINTER :: probe
1541 : LOGICAL, DIMENSION(:), INTENT(IN), OPTIONAL :: added_mos_auto
1542 : LOGICAL, INTENT(OUT), OPTIONAL :: added_mos_auto_grow
1543 : LOGICAL, INTENT(IN), OPTIONAL :: separate_spin_occupations
1544 :
1545 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_set_mo_occupation'
1546 :
1547 : INTEGER :: handle, ik, ikpgr, ispin, kplocal, nao, &
1548 : nb, ncol_global, ne_a, ne_b, &
1549 : nelectron, nkp, nmo, nrow_global, nspin
1550 : INTEGER, DIMENSION(2) :: kp_range
1551 : LOGICAL :: my_added_mos_auto_grow, &
1552 : my_separate_spin_occupations
1553 : REAL(KIND=dp) :: kTS, kTS_spin(2), mu, mus(2), nel
1554 36112 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: smatrix
1555 36112 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: weig, wocc
1556 36112 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: icoeff, rcoeff
1557 36112 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation, wkp
1558 : TYPE(cp_fm_type), POINTER :: mo_coeff
1559 : TYPE(kpoint_env_type), POINTER :: kp
1560 : TYPE(mo_set_type), POINTER :: mo_set
1561 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
1562 :
1563 36112 : CALL timeset(routineN, handle)
1564 :
1565 36112 : my_added_mos_auto_grow = .FALSE.
1566 36112 : IF (PRESENT(added_mos_auto_grow)) added_mos_auto_grow = .FALSE.
1567 :
1568 : ! first collect all the eigenvalues
1569 36112 : CALL get_kpoint_info(kpoint, nkp=nkp)
1570 36112 : kp => kpoint%kp_env(1)%kpoint_env
1571 36112 : nspin = SIZE(kp%mos, 2)
1572 36112 : my_separate_spin_occupations = smear%fixed_mag_mom > 0.0_dp
1573 36112 : IF (PRESENT(separate_spin_occupations)) THEN
1574 2092 : my_separate_spin_occupations = my_separate_spin_occupations .OR. separate_spin_occupations
1575 : END IF
1576 36112 : mo_set => kp%mos(1, 1)
1577 36112 : CALL get_mo_set(mo_set, nmo=nmo, nao=nao, nelectron=nelectron)
1578 36112 : ne_a = nelectron
1579 36112 : IF (nspin == 2) THEN
1580 2222 : CALL get_mo_set(kp%mos(1, 2), nmo=nb, nelectron=ne_b)
1581 2222 : CPASSERT(nmo == nb)
1582 : END IF
1583 288896 : ALLOCATE (weig(nmo, nkp, nspin), wocc(nmo, nkp, nspin))
1584 36112 : weig = 0.0_dp
1585 36112 : wocc = 0.0_dp
1586 36112 : IF (PRESENT(probe)) THEN
1587 0 : ALLOCATE (rcoeff(nao, nmo, nkp, nspin), icoeff(nao, nmo, nkp, nspin))
1588 0 : rcoeff = 0.0_dp !coeff, real part
1589 0 : icoeff = 0.0_dp !coeff, imaginary part
1590 : END IF
1591 36112 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
1592 36112 : kplocal = kp_range(2) - kp_range(1) + 1
1593 102274 : DO ikpgr = 1, kplocal
1594 66162 : ik = kp_range(1) + ikpgr - 1
1595 66162 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1596 174970 : DO ispin = 1, nspin
1597 72696 : mo_set => kp%mos(1, ispin)
1598 72696 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
1599 1301550 : weig(1:nmo, ik, ispin) = eigenvalues(1:nmo)
1600 138858 : IF (PRESENT(probe)) THEN
1601 0 : CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
1602 : CALL cp_fm_get_info(mo_coeff, &
1603 : nrow_global=nrow_global, &
1604 0 : ncol_global=ncol_global)
1605 0 : ALLOCATE (smatrix(nrow_global, ncol_global))
1606 0 : CALL cp_fm_get_submatrix(mo_coeff, smatrix)
1607 :
1608 0 : rcoeff(1:nao, 1:nmo, ik, ispin) = smatrix(1:nrow_global, 1:ncol_global)
1609 :
1610 0 : DEALLOCATE (smatrix)
1611 :
1612 0 : mo_set => kp%mos(2, ispin)
1613 :
1614 0 : CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
1615 : CALL cp_fm_get_info(mo_coeff, &
1616 : nrow_global=nrow_global, &
1617 0 : ncol_global=ncol_global)
1618 0 : ALLOCATE (smatrix(nrow_global, ncol_global))
1619 0 : CALL cp_fm_get_submatrix(mo_coeff, smatrix)
1620 :
1621 0 : icoeff(1:nao, 1:nmo, ik, ispin) = smatrix(1:nrow_global, 1:ncol_global)
1622 :
1623 0 : mo_set => kp%mos(1, ispin)
1624 :
1625 0 : DEALLOCATE (smatrix)
1626 : END IF
1627 : END DO
1628 : END DO
1629 36112 : CALL get_kpoint_info(kpoint, para_env_inter_kp=para_env_inter_kp)
1630 36112 : CALL para_env_inter_kp%sum(weig)
1631 :
1632 36112 : IF (PRESENT(probe)) THEN
1633 0 : CALL para_env_inter_kp%sum(rcoeff)
1634 0 : CALL para_env_inter_kp%sum(icoeff)
1635 : END IF
1636 :
1637 36112 : CALL get_kpoint_info(kpoint, wkp=wkp)
1638 36112 : kTS_spin = 0.0_dp
1639 :
1640 : !calling of HP module HERE, before smear
1641 36112 : IF (PRESENT(probe)) THEN
1642 0 : smear%do_smear = .FALSE. !ensures smearing is switched off
1643 :
1644 0 : IF (nspin == 1) THEN
1645 0 : nel = REAL(nelectron, KIND=dp)
1646 : CALL probe_occupancy_kp(wocc(:, :, :), mus(1), kTS, weig(:, :, :), rcoeff(:, :, :, :), icoeff(:, :, :, :), 2.0d0, &
1647 0 : probe, nel, wkp)
1648 : ELSE
1649 0 : nel = REAL(ne_a, KIND=dp) + REAL(ne_b, KIND=dp)
1650 : CALL probe_occupancy_kp(wocc(:, :, :), mu, kTS, weig(:, :, :), rcoeff(:, :, :, :), icoeff(:, :, :, :), 1.0d0, &
1651 0 : probe, nel, wkp)
1652 0 : kTS = kTS/2._dp
1653 0 : mus(1:2) = mu
1654 : END IF
1655 :
1656 0 : DO ikpgr = 1, kplocal
1657 0 : ik = kp_range(1) + ikpgr - 1
1658 0 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1659 0 : DO ispin = 1, nspin
1660 0 : mo_set => kp%mos(1, ispin)
1661 0 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues, occupation_numbers=occupation)
1662 0 : eigenvalues(1:nmo) = weig(1:nmo, ik, ispin)
1663 0 : occupation(1:nmo) = wocc(1:nmo, ik, ispin)
1664 0 : mo_set%kTS = kTS
1665 0 : mo_set%mu = mus(ispin)
1666 :
1667 0 : CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
1668 : !get smatrix for kpoint_env ikp
1669 : CALL cp_fm_get_info(mo_coeff, &
1670 : nrow_global=nrow_global, &
1671 0 : ncol_global=ncol_global)
1672 0 : ALLOCATE (smatrix(nrow_global, ncol_global))
1673 0 : CALL cp_fm_get_submatrix(mo_coeff, smatrix)
1674 :
1675 0 : smatrix(1:nrow_global, 1:ncol_global) = rcoeff(1:nao, 1:nmo, ik, ispin)
1676 0 : DEALLOCATE (smatrix)
1677 :
1678 0 : mo_set => kp%mos(2, ispin)
1679 :
1680 0 : CALL get_mo_set(mo_set, mo_coeff=mo_coeff)
1681 : !get smatrix for kpoint_env ikp
1682 : CALL cp_fm_get_info(mo_coeff, &
1683 : nrow_global=nrow_global, &
1684 0 : ncol_global=ncol_global)
1685 0 : ALLOCATE (smatrix(nrow_global, ncol_global))
1686 0 : CALL cp_fm_get_submatrix(mo_coeff, smatrix)
1687 :
1688 0 : smatrix(1:nrow_global, 1:ncol_global) = icoeff(1:nao, 1:nmo, ik, ispin)
1689 0 : DEALLOCATE (smatrix)
1690 :
1691 0 : mo_set => kp%mos(1, ispin)
1692 :
1693 : END DO
1694 : END DO
1695 :
1696 0 : DEALLOCATE (weig, wocc, rcoeff, icoeff)
1697 :
1698 : END IF
1699 :
1700 : IF (PRESENT(probe) .EQV. .FALSE.) THEN
1701 36112 : IF (smear%do_smear) THEN
1702 28356 : SELECT CASE (smear%method)
1703 : CASE (smear_fermi_dirac)
1704 : ! finite electronic temperature
1705 14036 : IF (nspin == 1) THEN
1706 13362 : nel = REAL(nelectron, KIND=dp)
1707 : CALL Smearkp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
1708 13362 : smear%electronic_temperature, 2.0_dp, smear_fermi_dirac)
1709 13362 : kTS_spin(1) = kTS
1710 674 : ELSE IF (my_separate_spin_occupations) THEN
1711 86 : nel = REAL(ne_a, KIND=dp)
1712 : CALL Smearkp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
1713 86 : smear%electronic_temperature, 1.0_dp, smear_fermi_dirac)
1714 86 : kTS_spin(1) = kTS
1715 86 : nel = REAL(ne_b, KIND=dp)
1716 : CALL Smearkp(wocc(:, :, 2), mus(2), kTS, weig(:, :, 2), nel, wkp, &
1717 86 : smear%electronic_temperature, 1.0_dp, smear_fermi_dirac)
1718 86 : kTS_spin(2) = kTS
1719 : ELSE
1720 588 : nel = REAL(ne_a, KIND=dp) + REAL(ne_b, KIND=dp)
1721 : CALL Smearkp2(wocc(:, :, :), mu, kTS, weig(:, :, :), nel, wkp, &
1722 588 : smear%electronic_temperature, smear_fermi_dirac)
1723 588 : kTS = kTS/2._dp
1724 1764 : kTS_spin(1:2) = kTS
1725 1764 : mus(1:2) = mu
1726 : END IF
1727 : CASE (smear_gaussian, smear_mp, smear_mv)
1728 284 : IF (nspin == 1) THEN
1729 252 : nel = REAL(nelectron, KIND=dp)
1730 : CALL Smearkp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
1731 252 : smear%smearing_width, 2.0_dp, smear%method)
1732 252 : kTS_spin(1) = kTS
1733 32 : ELSE IF (my_separate_spin_occupations) THEN
1734 0 : nel = REAL(ne_a, KIND=dp)
1735 : CALL Smearkp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
1736 0 : smear%smearing_width, 1.0_dp, smear%method)
1737 0 : kTS_spin(1) = kTS
1738 0 : nel = REAL(ne_b, KIND=dp)
1739 : CALL Smearkp(wocc(:, :, 2), mus(2), kTS, weig(:, :, 2), nel, wkp, &
1740 0 : smear%smearing_width, 1.0_dp, smear%method)
1741 0 : kTS_spin(2) = kTS
1742 : ELSE
1743 32 : nel = REAL(ne_a, KIND=dp) + REAL(ne_b, KIND=dp)
1744 : CALL Smearkp2(wocc(:, :, :), mu, kTS, weig(:, :, :), nel, wkp, &
1745 32 : smear%smearing_width, smear%method)
1746 32 : kTS = kTS/2._dp
1747 96 : kTS_spin(1:2) = kTS
1748 96 : mus(1:2) = mu
1749 : END IF
1750 : CASE DEFAULT
1751 14320 : CPABORT("kpoints: Selected smearing not (yet) supported")
1752 : END SELECT
1753 : ELSE
1754 : ! fixed occupations (2/1)
1755 21792 : IF (nspin == 1) THEN
1756 20276 : nel = REAL(nelectron, KIND=dp)
1757 : CALL Smearkp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
1758 20276 : 0.0_dp, 2.0_dp, smear_gaussian)
1759 20276 : kTS_spin(1) = kTS
1760 : ELSE
1761 1516 : nel = REAL(ne_a, KIND=dp)
1762 : CALL Smearkp(wocc(:, :, 1), mus(1), kTS, weig(:, :, 1), nel, wkp, &
1763 1516 : 0.0_dp, 1.0_dp, smear_gaussian)
1764 1516 : kTS_spin(1) = kTS
1765 1516 : nel = REAL(ne_b, KIND=dp)
1766 : CALL Smearkp(wocc(:, :, 2), mus(2), kTS, weig(:, :, 2), nel, wkp, &
1767 1516 : 0.0_dp, 1.0_dp, smear_gaussian)
1768 1516 : kTS_spin(2) = kTS
1769 : END IF
1770 : END IF
1771 36112 : IF (smear%do_smear .AND. PRESENT(added_mos_auto)) THEN
1772 : CALL kpoint_check_added_mos_auto_occupation( &
1773 13780 : wocc, wkp, smear, nspin, added_mos_auto, nao, my_added_mos_auto_grow)
1774 : END IF
1775 36112 : IF (my_added_mos_auto_grow) THEN
1776 6 : IF (PRESENT(added_mos_auto_grow)) THEN
1777 6 : added_mos_auto_grow = .TRUE.
1778 6 : DEALLOCATE (weig, wocc)
1779 6 : CALL timestop(handle)
1780 6 : RETURN
1781 : ELSE
1782 : CALL cp_abort(__LOCATION__, &
1783 : "K-point ADDED_MOS AUTO needs a larger virtual-space buffer, "// &
1784 0 : "but this call path cannot grow it.")
1785 : END IF
1786 : END IF
1787 102260 : DO ikpgr = 1, kplocal
1788 66154 : ik = kp_range(1) + ikpgr - 1
1789 66154 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1790 174944 : DO ispin = 1, nspin
1791 72684 : mo_set => kp%mos(1, ispin)
1792 72684 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues, occupation_numbers=occupation)
1793 1301466 : eigenvalues(1:nmo) = weig(1:nmo, ik, ispin)
1794 1301466 : occupation(1:nmo) = wocc(1:nmo, ik, ispin)
1795 72684 : mo_set%kTS = kTS_spin(ispin)
1796 138838 : mo_set%mu = mus(ispin)
1797 : END DO
1798 : END DO
1799 :
1800 36106 : DEALLOCATE (weig, wocc)
1801 :
1802 : END IF
1803 :
1804 36106 : CALL timestop(handle)
1805 :
1806 108336 : END SUBROUTINE kpoint_set_mo_occupation
1807 :
1808 : ! **************************************************************************************************
1809 : !> \brief summarize whether weighted k-point smearing reaches the available band edges
1810 : !> \param wocc ...
1811 : !> \param wkp ...
1812 : !> \param smear ...
1813 : !> \param nspin ...
1814 : !> \param has_weight ...
1815 : !> \param first_fractional ...
1816 : !> \param last_occupied ...
1817 : !> \param last_occupied_spin ...
1818 : ! **************************************************************************************************
1819 788 : SUBROUTINE kpoint_smearing_edge_status(wocc, wkp, smear, nspin, has_weight, &
1820 788 : first_fractional, last_occupied, last_occupied_spin)
1821 :
1822 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: wocc
1823 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wkp
1824 : TYPE(smear_type), INTENT(IN) :: smear
1825 : INTEGER, INTENT(IN) :: nspin
1826 : LOGICAL, INTENT(OUT) :: has_weight, first_fractional, &
1827 : last_occupied
1828 : LOGICAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: last_occupied_spin
1829 :
1830 : INTEGER :: ik, ispin, nmo
1831 : LOGICAL :: band_occupied
1832 : REAL(KIND=dp) :: eps_occ, maxocc, weight, weight_threshold
1833 :
1834 788 : has_weight = .FALSE.
1835 788 : first_fractional = .FALSE.
1836 788 : last_occupied = .FALSE.
1837 1756 : IF (PRESENT(last_occupied_spin)) last_occupied_spin(:) = .FALSE.
1838 788 : IF (.NOT. smear%do_smear) RETURN
1839 :
1840 786 : nmo = SIZE(wocc, 1)
1841 786 : IF (nmo < 1) RETURN
1842 786 : eps_occ = MAX(smear%eps_fermi_dirac, 10.0_dp*EPSILON(1.0_dp))
1843 786 : weight_threshold = 10.0_dp*EPSILON(1.0_dp)
1844 :
1845 1752 : DO ispin = 1, MIN(nspin, SIZE(wocc, 3))
1846 966 : maxocc = MERGE(2.0_dp, 1.0_dp, nspin == 1)
1847 2946 : DO ik = 1, MIN(SIZE(wkp), SIZE(wocc, 2))
1848 1194 : weight = ABS(wkp(ik))
1849 1194 : IF (weight <= weight_threshold) CYCLE
1850 1182 : has_weight = .TRUE.
1851 : first_fractional = first_fractional .OR. &
1852 1182 : ABS(wocc(1, ik, ispin) - maxocc) > eps_occ
1853 1182 : band_occupied = ABS(wocc(nmo, ik, ispin)) > eps_occ
1854 1182 : last_occupied = last_occupied .OR. band_occupied
1855 2148 : IF (band_occupied .AND. PRESENT(last_occupied_spin)) THEN
1856 8 : IF (ispin <= SIZE(last_occupied_spin)) last_occupied_spin(ispin) = .TRUE.
1857 : END IF
1858 : END DO
1859 : END DO
1860 :
1861 : END SUBROUTINE kpoint_smearing_edge_status
1862 :
1863 : ! **************************************************************************************************
1864 : !> \brief request ADDED_MOS AUTO growth if weighted k-point smearing occupies the last band
1865 : !> \param wocc ...
1866 : !> \param wkp ...
1867 : !> \param smear ...
1868 : !> \param nspin ...
1869 : !> \param added_mos_auto ...
1870 : !> \param nao ...
1871 : !> \param added_mos_auto_grow ...
1872 : ! **************************************************************************************************
1873 13780 : SUBROUTINE kpoint_check_added_mos_auto_occupation(wocc, wkp, smear, nspin, &
1874 13780 : added_mos_auto, nao, added_mos_auto_grow)
1875 :
1876 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: wocc
1877 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wkp
1878 : TYPE(smear_type), INTENT(IN) :: smear
1879 : INTEGER, INTENT(IN) :: nspin
1880 : LOGICAL, DIMENSION(:), INTENT(IN) :: added_mos_auto
1881 : INTEGER, INTENT(IN) :: nao
1882 : LOGICAL, INTENT(OUT) :: added_mos_auto_grow
1883 :
1884 : INTEGER :: ispin, nmo
1885 : LOGICAL :: first_fractional, has_weight, &
1886 : last_occupied
1887 13780 : LOGICAL, DIMENSION(nspin) :: last_occupied_spin
1888 :
1889 13780 : added_mos_auto_grow = .FALSE.
1890 39784 : IF (.NOT. ANY(added_mos_auto) .OR. .NOT. smear%do_smear) RETURN
1891 :
1892 : CALL kpoint_smearing_edge_status(wocc, wkp, smear, nspin, has_weight, &
1893 : first_fractional, last_occupied, &
1894 778 : last_occupied_spin=last_occupied_spin)
1895 778 : IF (.NOT. has_weight) RETURN
1896 :
1897 778 : nmo = SIZE(wocc, 1)
1898 1720 : DO ispin = 1, MIN(nspin, SIZE(added_mos_auto))
1899 948 : IF (.NOT. added_mos_auto(ispin) .OR. .NOT. last_occupied_spin(ispin)) CYCLE
1900 6 : IF (nmo >= nao) THEN
1901 : CALL cp_abort(__LOCATION__, &
1902 : "K-point ADDED_MOS AUTO exhausted the AO basis but the highest band "// &
1903 0 : "is still occupied. Use a larger basis or reduce the smearing width.")
1904 : END IF
1905 6 : added_mos_auto_grow = .TRUE.
1906 1714 : RETURN
1907 : END DO
1908 :
1909 : END SUBROUTINE kpoint_check_added_mos_auto_occupation
1910 :
1911 : ! **************************************************************************************************
1912 : !> \brief Calculate kpoint density matrices (rho(k), owned by kpoint groups)
1913 : !> \param kpoint kpoint environment
1914 : !> \param energy_weighted calculate energy weighted density matrix
1915 : !> \param for_aux_fit ...
1916 : ! **************************************************************************************************
1917 113370 : SUBROUTINE kpoint_density_matrices(kpoint, energy_weighted, for_aux_fit)
1918 :
1919 : TYPE(kpoint_type), POINTER :: kpoint
1920 : LOGICAL, OPTIONAL :: energy_weighted, for_aux_fit
1921 :
1922 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_density_matrices'
1923 :
1924 : INTEGER :: handle, ikpgr, ispin, kplocal, nao, nmo, &
1925 : nspin, omp_threads
1926 : INTEGER, DIMENSION(2) :: kp_range
1927 : LOGICAL :: aux_fit, wtype
1928 37790 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation
1929 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1930 : TYPE(cp_fm_type) :: fwork
1931 : TYPE(cp_fm_type), POINTER :: cpmat, pmat, rpmat
1932 : TYPE(kpoint_env_type), POINTER :: kp
1933 : TYPE(mo_set_type), POINTER :: mo_set
1934 :
1935 37790 : CALL timeset(routineN, handle)
1936 :
1937 37790 : IF (PRESENT(energy_weighted)) THEN
1938 454 : wtype = energy_weighted
1939 : ELSE
1940 : ! default is normal density matrix
1941 37336 : wtype = .FALSE.
1942 : END IF
1943 :
1944 37790 : IF (PRESENT(for_aux_fit)) THEN
1945 124 : aux_fit = for_aux_fit
1946 : ELSE
1947 37666 : aux_fit = .FALSE.
1948 : END IF
1949 :
1950 37790 : IF (aux_fit) THEN
1951 124 : CPASSERT(ASSOCIATED(kpoint%kp_aux_env))
1952 : END IF
1953 :
1954 : ! work matrix
1955 37790 : IF (aux_fit) THEN
1956 124 : mo_set => kpoint%kp_aux_env(1)%kpoint_env%mos(1, 1)
1957 : ELSE
1958 37666 : mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)
1959 : END IF
1960 37790 : CALL get_mo_set(mo_set, nao=nao, nmo=nmo)
1961 37790 : CALL cp_fm_get_info(mo_set%mo_coeff, matrix_struct=matrix_struct)
1962 :
1963 37790 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
1964 37790 : kplocal = kp_range(2) - kp_range(1) + 1
1965 37790 : IF (aux_fit) THEN
1966 124 : kp => kpoint%kp_aux_env(1)%kpoint_env
1967 : ELSE
1968 37666 : kp => kpoint%kp_env(1)%kpoint_env
1969 : END IF
1970 37790 : nspin = SIZE(kp%mos, 2)
1971 37790 : omp_threads = 1
1972 37790 : !$ omp_threads = omp_get_max_threads()
1973 37790 : IF (kpoint%blacs_env_all%para_env%num_pe == 1 .AND. omp_threads > 1 .AND. &
1974 : kplocal*nspin > 1) THEN
1975 0 : CALL kpoint_density_matrices_omp(kpoint, wtype, aux_fit, nao, nmo, kplocal, nspin)
1976 0 : CALL timestop(handle)
1977 0 : RETURN
1978 : END IF
1979 :
1980 37790 : CALL cp_fm_create(fwork, matrix_struct)
1981 109502 : DO ikpgr = 1, kplocal
1982 71712 : IF (aux_fit) THEN
1983 1876 : kp => kpoint%kp_aux_env(ikpgr)%kpoint_env
1984 : ELSE
1985 69836 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1986 : END IF
1987 71712 : nspin = SIZE(kp%mos, 2)
1988 188262 : DO ispin = 1, nspin
1989 78760 : mo_set => kp%mos(1, ispin)
1990 78760 : IF (wtype) THEN
1991 1550 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
1992 : END IF
1993 150472 : IF (kpoint%use_real_wfn) THEN
1994 168 : IF (wtype) THEN
1995 12 : pmat => kp%wmat(1, ispin)
1996 : ELSE
1997 156 : pmat => kp%pmat(1, ispin)
1998 : END IF
1999 168 : CALL get_mo_set(mo_set, occupation_numbers=occupation)
2000 168 : CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
2001 168 : CALL cp_fm_column_scale(fwork, occupation)
2002 168 : IF (wtype) THEN
2003 12 : CALL cp_fm_column_scale(fwork, eigenvalues)
2004 : END IF
2005 168 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, pmat)
2006 : ELSE
2007 78592 : IF (wtype) THEN
2008 1538 : rpmat => kp%wmat(1, ispin)
2009 1538 : cpmat => kp%wmat(2, ispin)
2010 : ELSE
2011 77054 : rpmat => kp%pmat(1, ispin)
2012 77054 : cpmat => kp%pmat(2, ispin)
2013 : END IF
2014 78592 : CALL get_mo_set(mo_set, occupation_numbers=occupation)
2015 78592 : CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
2016 78592 : CALL cp_fm_column_scale(fwork, occupation)
2017 78592 : IF (wtype) THEN
2018 1538 : CALL cp_fm_column_scale(fwork, eigenvalues)
2019 : END IF
2020 : ! Re(c)*Re(c)
2021 78592 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, rpmat)
2022 78592 : mo_set => kp%mos(2, ispin)
2023 : ! Im(c)*Re(c)
2024 78592 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, cpmat)
2025 : ! Re(c)*Im(c)
2026 78592 : CALL parallel_gemm("N", "T", nao, nao, nmo, -1.0_dp, fwork, mo_set%mo_coeff, 1.0_dp, cpmat)
2027 78592 : CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
2028 78592 : CALL cp_fm_column_scale(fwork, occupation)
2029 78592 : IF (wtype) THEN
2030 1538 : CALL cp_fm_column_scale(fwork, eigenvalues)
2031 : END IF
2032 : ! Im(c)*Im(c)
2033 78592 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 1.0_dp, rpmat)
2034 : END IF
2035 : END DO
2036 : END DO
2037 :
2038 37790 : CALL cp_fm_release(fwork)
2039 :
2040 37790 : CALL timestop(handle)
2041 :
2042 37790 : END SUBROUTINE kpoint_density_matrices
2043 :
2044 : ! **************************************************************************************************
2045 : !> \brief Build W(k) for noncanonical complex OT orbitals from H(k) C(k).
2046 : !> The occupied-space Lagrange multiplier is
2047 : !> Hermitian[(C^H H C) f], which reduces to the usual eigenvalue-weighted
2048 : !> expression for canonical orbitals.
2049 : !> \param coeff_re real part of C(k)
2050 : !> \param coeff_im imaginary part of C(k)
2051 : !> \param hc_re real part of H(k) C(k)
2052 : !> \param hc_im imaginary part of H(k) C(k)
2053 : !> \param occupation orbital occupations
2054 : !> \param wmat_re real part of W(k)
2055 : !> \param wmat_im imaginary part of W(k)
2056 : ! **************************************************************************************************
2057 34 : SUBROUTINE kpoint_ot_energy_weighted_density(coeff_re, coeff_im, hc_re, hc_im, occupation, &
2058 : wmat_re, wmat_im)
2059 :
2060 : TYPE(cp_fm_type), INTENT(IN) :: coeff_re, coeff_im, hc_re, hc_im
2061 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occupation
2062 : TYPE(cp_fm_type), INTENT(IN) :: wmat_re, wmat_im
2063 :
2064 : COMPLEX(KIND=dp), PARAMETER :: zhalf = (0.5_dp, 0.0_dp), &
2065 : zone = (1.0_dp, 0.0_dp), &
2066 : zzero = (0.0_dp, 0.0_dp)
2067 :
2068 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: occupation_complex
2069 : INTEGER :: nao, nmo
2070 : TYPE(cp_cfm_type) :: coeff, hblock, hblock_h, hc, weighted, &
2071 : wmat
2072 :
2073 34 : CALL cp_fm_get_info(coeff_re, nrow_global=nao, ncol_global=nmo)
2074 34 : CPASSERT(nmo >= 1)
2075 34 : CPASSERT(SIZE(occupation) >= nmo)
2076 :
2077 34 : CALL cp_cfm_create(coeff, coeff_re%matrix_struct)
2078 34 : CALL cp_cfm_create(hc, hc_re%matrix_struct)
2079 34 : CALL cp_cfm_create(weighted, coeff_re%matrix_struct)
2080 34 : CALL cp_cfm_create(hblock, coeff_re%matrix_struct, nrow=nmo, ncol=nmo)
2081 34 : CALL cp_cfm_create(hblock_h, coeff_re%matrix_struct, nrow=nmo, ncol=nmo)
2082 34 : CALL cp_cfm_create(wmat, wmat_re%matrix_struct)
2083 :
2084 34 : CALL cp_fm_to_cfm(coeff_re, coeff_im, coeff)
2085 34 : CALL cp_fm_to_cfm(hc_re, hc_im, hc)
2086 34 : CALL cp_cfm_gemm("C", "N", nmo, nmo, nao, zone, coeff, hc, zzero, hblock)
2087 :
2088 102 : ALLOCATE (occupation_complex(nmo))
2089 158 : occupation_complex(:) = CMPLX(occupation(1:nmo), 0.0_dp, KIND=dp)
2090 34 : CALL cp_cfm_column_scale(hblock, occupation_complex)
2091 34 : DEALLOCATE (occupation_complex)
2092 34 : CALL cp_cfm_transpose(hblock, "C", hblock_h)
2093 34 : CALL cp_cfm_scale_and_add(zhalf, hblock, zhalf, hblock_h)
2094 :
2095 34 : CALL cp_cfm_gemm("N", "N", nao, nmo, nmo, zone, coeff, hblock, zzero, weighted)
2096 34 : CALL cp_cfm_gemm("N", "C", nao, nao, nmo, zone, weighted, coeff, zzero, wmat)
2097 34 : CALL cp_cfm_to_fm(wmat, wmat_re, wmat_im)
2098 :
2099 34 : CALL cp_cfm_release(wmat)
2100 34 : CALL cp_cfm_release(hblock_h)
2101 34 : CALL cp_cfm_release(hblock)
2102 34 : CALL cp_cfm_release(weighted)
2103 34 : CALL cp_cfm_release(hc)
2104 34 : CALL cp_cfm_release(coeff)
2105 :
2106 68 : END SUBROUTINE kpoint_ot_energy_weighted_density
2107 :
2108 : ! **************************************************************************************************
2109 : !> \brief Build energy-weighted density matrices for complex K-point OT.
2110 : !> H(k) is transformed only while W(k) is needed, avoiding a persistent
2111 : !> dense AO-by-AO Hamiltonian for every spin and K point.
2112 : !> \param kpoint K-point environment
2113 : !> \param matrix_ks real-space Kohn-Sham matrices
2114 : ! **************************************************************************************************
2115 18 : SUBROUTINE kpoint_ot_energy_weighted_matrices(kpoint, matrix_ks)
2116 :
2117 : TYPE(kpoint_type), POINTER :: kpoint
2118 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks
2119 :
2120 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_ot_energy_weighted_matrices'
2121 :
2122 : INTEGER :: handle, ikpgr, ikpoint, ispin, nmo, nspin
2123 18 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2124 : LOGICAL :: use_real_wfn
2125 18 : REAL(KIND=dp), DIMENSION(:), POINTER :: occupation
2126 18 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
2127 : TYPE(cp_fm_type) :: hc_im, hc_re
2128 : TYPE(cp_fm_type), POINTER :: coeff_im, coeff_re
2129 : TYPE(dbcsr_type), POINTER :: matrix_k_im, matrix_k_re
2130 : TYPE(kpoint_env_type), POINTER :: kp
2131 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2132 18 : POINTER :: sab_nl
2133 :
2134 18 : CALL timeset(routineN, handle)
2135 18 : NULLIFY (cell_to_index, coeff_im, coeff_re, matrix_k_im, matrix_k_re, occupation, sab_nl, xkp)
2136 :
2137 : CALL get_kpoint_info(kpoint, xkp=xkp, use_real_wfn=use_real_wfn, &
2138 18 : cell_to_index=cell_to_index, sab_nl=sab_nl)
2139 18 : CPASSERT(.NOT. use_real_wfn)
2140 18 : CPASSERT(ASSOCIATED(cell_to_index))
2141 18 : CPASSERT(ASSOCIATED(sab_nl))
2142 18 : CPASSERT(ASSOCIATED(matrix_ks))
2143 :
2144 18 : ALLOCATE (matrix_k_re, matrix_k_im)
2145 : CALL dbcsr_create(matrix_k_re, template=matrix_ks(1, 1)%matrix, &
2146 18 : matrix_type=dbcsr_type_symmetric)
2147 : CALL dbcsr_create(matrix_k_im, template=matrix_ks(1, 1)%matrix, &
2148 18 : matrix_type=dbcsr_type_antisymmetric)
2149 18 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_k_re, sab_nl)
2150 18 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_k_im, sab_nl)
2151 :
2152 42 : DO ikpgr = 1, SIZE(kpoint%kp_env)
2153 24 : kp => kpoint%kp_env(ikpgr)%kpoint_env
2154 24 : ikpoint = kp%nkpoint
2155 24 : nspin = SIZE(kp%mos, 2)
2156 72 : DO ispin = 1, nspin
2157 30 : CALL dbcsr_set(matrix_k_re, 0.0_dp)
2158 30 : CALL dbcsr_set(matrix_k_im, 0.0_dp)
2159 : CALL rskp_transform(rmatrix=matrix_k_re, cmatrix=matrix_k_im, rsmat=matrix_ks, &
2160 : ispin=ispin, xkp=xkp(1:3, ikpoint), &
2161 30 : cell_to_index=cell_to_index, sab_nl=sab_nl)
2162 :
2163 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=coeff_re, nmo=nmo, &
2164 30 : occupation_numbers=occupation)
2165 30 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=coeff_im)
2166 30 : CPASSERT(nmo >= 1)
2167 30 : CPASSERT(SIZE(occupation) >= nmo)
2168 :
2169 : ! Smearing occupations are k-point resolved, whereas mo_set%homo is not.
2170 30 : CALL cp_fm_create(hc_re, coeff_re%matrix_struct, set_zero=.TRUE.)
2171 30 : CALL cp_fm_create(hc_im, coeff_im%matrix_struct, set_zero=.TRUE.)
2172 30 : CALL cp_dbcsr_sm_fm_multiply(matrix_k_re, coeff_re, hc_re, nmo)
2173 : CALL cp_dbcsr_sm_fm_multiply(matrix_k_im, coeff_im, hc_re, nmo, &
2174 30 : alpha=-1.0_dp, beta=1.0_dp)
2175 30 : CALL cp_dbcsr_sm_fm_multiply(matrix_k_re, coeff_im, hc_im, nmo)
2176 : CALL cp_dbcsr_sm_fm_multiply(matrix_k_im, coeff_re, hc_im, nmo, &
2177 30 : alpha=1.0_dp, beta=1.0_dp)
2178 :
2179 : CALL kpoint_ot_energy_weighted_density( &
2180 : coeff_re, coeff_im, hc_re, hc_im, occupation, &
2181 30 : kp%wmat(1, ispin), kp%wmat(2, ispin))
2182 :
2183 30 : CALL cp_fm_release(hc_im)
2184 84 : CALL cp_fm_release(hc_re)
2185 : END DO
2186 : END DO
2187 :
2188 18 : CALL dbcsr_deallocate_matrix(matrix_k_im)
2189 18 : CALL dbcsr_deallocate_matrix(matrix_k_re)
2190 18 : CALL timestop(handle)
2191 :
2192 18 : END SUBROUTINE kpoint_ot_energy_weighted_matrices
2193 :
2194 : ! **************************************************************************************************
2195 : !> \brief Thread k-point density-matrix construction on a one-rank BLACS grid.
2196 : !> The local BLAS path avoids concurrent ScaLAPACK calls and gives every
2197 : !> worker a private column-scaled MO buffer.
2198 : !> \param kpoint ...
2199 : !> \param wtype ...
2200 : !> \param aux_fit ...
2201 : !> \param nao ...
2202 : !> \param nmo ...
2203 : !> \param kplocal ...
2204 : !> \param nspin ...
2205 : ! **************************************************************************************************
2206 0 : SUBROUTINE kpoint_density_matrices_omp(kpoint, wtype, aux_fit, nao, nmo, kplocal, nspin)
2207 :
2208 : TYPE(kpoint_type), POINTER :: kpoint
2209 : LOGICAL, INTENT(IN) :: wtype, aux_fit
2210 : INTEGER, INTENT(IN) :: nao, nmo, kplocal, nspin
2211 :
2212 : INTEGER :: ikpgr, ispin, nworkers, thread
2213 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: work
2214 :
2215 0 : nworkers = 1
2216 0 : !$ nworkers = MIN(omp_get_max_threads(), kplocal*nspin)
2217 0 : ALLOCATE (work(nao, nmo, nworkers))
2218 :
2219 : !$OMP PARALLEL DO COLLAPSE(2) DEFAULT(NONE) SCHEDULE(DYNAMIC, 1) NUM_THREADS(nworkers) &
2220 0 : !$OMP SHARED(kpoint,wtype,aux_fit,nao,nmo,kplocal,nspin,nworkers,work) PRIVATE(ikpgr,ispin,thread)
2221 : DO ikpgr = 1, kplocal
2222 : DO ispin = 1, nspin
2223 : thread = 1
2224 : !$ thread = omp_get_thread_num() + 1
2225 : CALL kpoint_density_matrix_job(kpoint, ikpgr, ispin, wtype, aux_fit, nao, nmo, &
2226 : work(:, :, thread))
2227 : END DO
2228 : END DO
2229 : !$OMP END PARALLEL DO
2230 :
2231 0 : DEALLOCATE (work)
2232 :
2233 0 : END SUBROUTINE kpoint_density_matrices_omp
2234 :
2235 : ! **************************************************************************************************
2236 : !> \brief Construct one (k-point,spin) density matrix using local BLAS.
2237 : !> \param kpoint ...
2238 : !> \param ikpgr ...
2239 : !> \param ispin ...
2240 : !> \param wtype ...
2241 : !> \param aux_fit ...
2242 : !> \param nao ...
2243 : !> \param nmo ...
2244 : !> \param work ...
2245 : ! **************************************************************************************************
2246 0 : SUBROUTINE kpoint_density_matrix_job(kpoint, ikpgr, ispin, wtype, aux_fit, nao, nmo, work)
2247 :
2248 : TYPE(kpoint_type), POINTER :: kpoint
2249 : INTEGER, INTENT(IN) :: ikpgr, ispin
2250 : LOGICAL, INTENT(IN) :: wtype, aux_fit
2251 : INTEGER, INTENT(IN) :: nao, nmo
2252 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: work
2253 :
2254 : INTEGER :: imo
2255 : REAL(KIND=dp) :: scale
2256 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2257 0 : POINTER :: ci, cmat, cr, rmat
2258 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation
2259 : TYPE(cp_fm_type), POINTER :: cpmat, rpmat
2260 : TYPE(kpoint_env_type), POINTER :: kp
2261 : TYPE(mo_set_type), POINTER :: mo_set
2262 :
2263 0 : IF (aux_fit) THEN
2264 0 : kp => kpoint%kp_aux_env(ikpgr)%kpoint_env
2265 : ELSE
2266 0 : kp => kpoint%kp_env(ikpgr)%kpoint_env
2267 : END IF
2268 0 : mo_set => kp%mos(1, ispin)
2269 0 : IF (wtype) THEN
2270 0 : CALL get_mo_set(mo_set, occupation_numbers=occupation, eigenvalues=eigenvalues)
2271 : ELSE
2272 0 : CALL get_mo_set(mo_set, occupation_numbers=occupation)
2273 : END IF
2274 0 : CALL cp_fm_get_info(mo_set%mo_coeff, local_data=cr)
2275 :
2276 0 : IF (wtype) THEN
2277 0 : rpmat => kp%wmat(1, ispin)
2278 : ELSE
2279 0 : rpmat => kp%pmat(1, ispin)
2280 : END IF
2281 0 : CALL cp_fm_get_info(rpmat, local_data=rmat)
2282 :
2283 0 : DO imo = 1, nmo
2284 0 : scale = occupation(imo)
2285 0 : IF (wtype) scale = scale*eigenvalues(imo)
2286 0 : work(1:nao, imo) = scale*cr(1:nao, imo)
2287 : END DO
2288 : CALL dgemm('N', 'T', nao, nao, nmo, 1.0_dp, cr, SIZE(cr, 1), &
2289 0 : work, SIZE(work, 1), 0.0_dp, rmat, SIZE(rmat, 1))
2290 :
2291 0 : IF (.NOT. kpoint%use_real_wfn) THEN
2292 0 : mo_set => kp%mos(2, ispin)
2293 0 : CALL cp_fm_get_info(mo_set%mo_coeff, local_data=ci)
2294 0 : IF (wtype) THEN
2295 0 : cpmat => kp%wmat(2, ispin)
2296 : ELSE
2297 0 : cpmat => kp%pmat(2, ispin)
2298 : END IF
2299 0 : CALL cp_fm_get_info(cpmat, local_data=cmat)
2300 :
2301 : CALL dgemm('N', 'T', nao, nao, nmo, 1.0_dp, ci, SIZE(ci, 1), &
2302 0 : work, SIZE(work, 1), 0.0_dp, cmat, SIZE(cmat, 1))
2303 : CALL dgemm('N', 'T', nao, nao, nmo, -1.0_dp, work, SIZE(work, 1), &
2304 0 : ci, SIZE(ci, 1), 1.0_dp, cmat, SIZE(cmat, 1))
2305 :
2306 0 : DO imo = 1, nmo
2307 0 : scale = occupation(imo)
2308 0 : IF (wtype) scale = scale*eigenvalues(imo)
2309 0 : work(1:nao, imo) = scale*ci(1:nao, imo)
2310 : END DO
2311 : CALL dgemm('N', 'T', nao, nao, nmo, 1.0_dp, ci, SIZE(ci, 1), &
2312 0 : work, SIZE(work, 1), 1.0_dp, rmat, SIZE(rmat, 1))
2313 : END IF
2314 :
2315 0 : END SUBROUTINE kpoint_density_matrix_job
2316 :
2317 : ! **************************************************************************************************
2318 : !> \brief Calculate Lowdin transformation of density matrix S^1/2 P S^1/2
2319 : !> Integrate diagonal elements over k-points to get Lowdin charges
2320 : !> \param kpoint kpoint environment
2321 : !> \param pmat_diag Sum over kpoints of diagonal elements
2322 : !> \par History
2323 : !> 04.2026 created [JGH]
2324 : ! **************************************************************************************************
2325 12 : SUBROUTINE lowdin_kp_trans(kpoint, pmat_diag)
2326 :
2327 : TYPE(kpoint_type), POINTER :: kpoint
2328 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: pmat_diag
2329 :
2330 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lowdin_kp_trans'
2331 : COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp, 0.0_dp), &
2332 : czero = (0.0_dp, 0.0_dp)
2333 :
2334 : INTEGER :: handle, ikpgr, ispin, kplocal, nao, nspin
2335 : INTEGER, DIMENSION(2) :: kp_range
2336 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dele
2337 : TYPE(cp_cfm_type) :: cf1work, cf2work
2338 : TYPE(cp_cfm_type), POINTER :: cshalf
2339 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
2340 : TYPE(cp_fm_type) :: f1work, f2work
2341 : TYPE(cp_fm_type), POINTER :: cpmat, pmat, rpmat, shalf
2342 : TYPE(kpoint_env_type), POINTER :: kp
2343 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
2344 :
2345 12 : CALL timeset(routineN, handle)
2346 :
2347 12 : nspin = SIZE(pmat_diag, 2)
2348 822 : pmat_diag = 0.0_dp
2349 :
2350 : ! work matrix
2351 : CALL cp_fm_get_info(kpoint%kp_env(1)%kpoint_env%pmat(1, 1), &
2352 12 : matrix_struct=matrix_struct, nrow_global=nao)
2353 12 : IF (kpoint%use_real_wfn) THEN
2354 0 : CALL cp_fm_create(f1work, matrix_struct, nrow=nao, ncol=nao)
2355 0 : CALL cp_fm_create(f2work, matrix_struct, nrow=nao, ncol=nao)
2356 : ELSE
2357 12 : CALL cp_fm_create(f2work, matrix_struct, nrow=nao, ncol=nao)
2358 12 : CALL cp_cfm_create(cf1work, matrix_struct, nrow=nao, ncol=nao)
2359 12 : CALL cp_cfm_create(cf2work, matrix_struct, nrow=nao, ncol=nao)
2360 : END IF
2361 36 : ALLOCATE (dele(nao))
2362 :
2363 12 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
2364 12 : kplocal = kp_range(2) - kp_range(1) + 1
2365 268 : DO ikpgr = 1, kplocal
2366 256 : kp => kpoint%kp_env(ikpgr)%kpoint_env
2367 564 : DO ispin = 1, nspin
2368 296 : IF (kpoint%use_real_wfn) THEN
2369 0 : pmat => kp%pmat(1, ispin)
2370 0 : shalf => kp%shalf
2371 0 : CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, pmat, shalf, 0.0_dp, f1work)
2372 0 : CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, shalf, f1work, 0.0_dp, f2work)
2373 : ELSE
2374 296 : rpmat => kp%pmat(1, ispin)
2375 296 : cpmat => kp%pmat(2, ispin)
2376 296 : cshalf => kp%cshalf
2377 296 : CALL cp_fm_to_cfm(rpmat, cpmat, cf1work)
2378 296 : CALL parallel_gemm("N", "N", nao, nao, nao, cone, cf1work, cshalf, czero, cf2work)
2379 296 : CALL parallel_gemm("N", "N", nao, nao, nao, cone, cshalf, cf2work, czero, cf1work)
2380 296 : CALL cp_cfm_to_fm(cf1work, mtargetr=f2work)
2381 : END IF
2382 296 : CALL cp_fm_get_diag(f2work, dele)
2383 4536 : pmat_diag(1:nao, ispin) = pmat_diag(1:nao, ispin) + kp%wkp*dele(1:nao)
2384 : END DO
2385 : END DO
2386 :
2387 12 : CALL get_kpoint_info(kpoint, para_env_inter_kp=para_env_inter_kp)
2388 1632 : CALL para_env_inter_kp%sum(pmat_diag)
2389 :
2390 12 : IF (kpoint%use_real_wfn) THEN
2391 0 : CALL cp_fm_release(f1work)
2392 0 : CALL cp_fm_release(f2work)
2393 : ELSE
2394 12 : CALL cp_fm_release(f2work)
2395 12 : CALL cp_cfm_release(cf1work)
2396 12 : CALL cp_cfm_release(cf2work)
2397 : END IF
2398 12 : DEALLOCATE (dele)
2399 :
2400 12 : CALL timestop(handle)
2401 :
2402 24 : END SUBROUTINE lowdin_kp_trans
2403 :
2404 : ! **************************************************************************************************
2405 : !> \brief Calculate S(k)^1/2 C(k) for real or complex k-point wavefunctions
2406 : !> \param kp K-point environment for one local k point
2407 : !> \param ispin Spin index
2408 : !> \param use_real_wfn Use real k-point wavefunctions
2409 : !> \param shalfc Output matrix containing S(k)^1/2 C(k) for real wavefunctions
2410 : !> \param cshalfc Output matrix containing S(k)^1/2 C(k) for complex wavefunctions
2411 : ! **************************************************************************************************
2412 0 : SUBROUTINE lowdin_kp_mo_coeff(kp, ispin, use_real_wfn, shalfc, cshalfc)
2413 :
2414 : TYPE(kpoint_env_type), POINTER :: kp
2415 : INTEGER, INTENT(IN) :: ispin
2416 : LOGICAL, INTENT(IN) :: use_real_wfn
2417 : TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL :: shalfc
2418 : TYPE(cp_cfm_type), INTENT(INOUT), OPTIONAL :: cshalfc
2419 :
2420 : INTEGER :: nao, nmo
2421 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_mo, matrix_struct_shalf
2422 : TYPE(cp_fm_type) :: cshalf_im, cshalf_re, shalf_im, shalf_re
2423 : TYPE(mo_set_type), POINTER :: mo_set, mo_set_im, mo_set_re
2424 :
2425 0 : IF (use_real_wfn) THEN
2426 0 : CPASSERT(PRESENT(shalfc))
2427 0 : mo_set => kp%mos(1, ispin)
2428 0 : CALL get_mo_set(mo_set, nao=nao, nmo=nmo)
2429 :
2430 : CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, kp%shalf, &
2431 0 : mo_set%mo_coeff, 0.0_dp, shalfc)
2432 : ELSE
2433 0 : CPASSERT(PRESENT(cshalfc))
2434 0 : mo_set_re => kp%mos(1, ispin)
2435 0 : mo_set_im => kp%mos(2, ispin)
2436 0 : CALL get_mo_set(mo_set_re, nao=nao, nmo=nmo)
2437 0 : CALL cp_fm_get_info(mo_set_re%mo_coeff, matrix_struct=matrix_struct_mo)
2438 0 : CALL cp_cfm_get_info(kp%cshalf, matrix_struct=matrix_struct_shalf)
2439 :
2440 0 : CALL cp_fm_create(shalf_re, matrix_struct_shalf, nrow=nao, ncol=nao)
2441 0 : CALL cp_fm_create(shalf_im, matrix_struct_shalf, nrow=nao, ncol=nao)
2442 0 : CALL cp_fm_create(cshalf_re, matrix_struct_mo, nrow=nao, ncol=nmo)
2443 0 : CALL cp_fm_create(cshalf_im, matrix_struct_mo, nrow=nao, ncol=nmo)
2444 :
2445 0 : CALL cp_cfm_to_fm(kp%cshalf, mtargetr=shalf_re, mtargeti=shalf_im)
2446 :
2447 : ! Re[S(k)^1/2 C(k)] = Re[S(k)^1/2] C_re(k) - Im[S(k)^1/2] C_im(k)
2448 : CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, shalf_re, &
2449 0 : mo_set_re%mo_coeff, 0.0_dp, cshalf_re)
2450 : CALL parallel_gemm("N", "N", nao, nmo, nao, -1.0_dp, shalf_im, &
2451 0 : mo_set_im%mo_coeff, 1.0_dp, cshalf_re)
2452 :
2453 : ! Im[S(k)^1/2 C(k)] = Re[S(k)^1/2] C_im(k) + Im[S(k)^1/2] C_re(k)
2454 : CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, shalf_re, &
2455 0 : mo_set_im%mo_coeff, 0.0_dp, cshalf_im)
2456 : CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, shalf_im, &
2457 0 : mo_set_re%mo_coeff, 1.0_dp, cshalf_im)
2458 :
2459 0 : CALL cp_fm_to_cfm(cshalf_re, cshalf_im, cshalfc)
2460 :
2461 0 : CALL cp_fm_release(shalf_re)
2462 0 : CALL cp_fm_release(shalf_im)
2463 0 : CALL cp_fm_release(cshalf_re)
2464 0 : CALL cp_fm_release(cshalf_im)
2465 : END IF
2466 :
2467 0 : END SUBROUTINE lowdin_kp_mo_coeff
2468 :
2469 : ! **************************************************************************************************
2470 : !> \brief generate real space density matrices in DBCSR format
2471 : !> \param kpoint Kpoint environment
2472 : !> \param denmat Real space (DBCSR) density matrices
2473 : !> \param wtype True = energy weighted density matrix
2474 : !> False = normal density matrix
2475 : !> \param tempmat DBCSR matrix to be used as template
2476 : !> \param sab_nl ...
2477 : !> \param fmwork FM work matrices (kpoint group)
2478 : !> \param for_aux_fit ...
2479 : !> \param pmat_ext ...
2480 : !> \param overlap_rs ...
2481 : ! **************************************************************************************************
2482 38254 : SUBROUTINE kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwork, for_aux_fit, &
2483 38254 : pmat_ext, overlap_rs)
2484 :
2485 : TYPE(kpoint_type), POINTER :: kpoint
2486 : TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
2487 : LOGICAL, INTENT(IN) :: wtype
2488 : TYPE(dbcsr_type), POINTER :: tempmat
2489 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2490 : POINTER :: sab_nl
2491 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fmwork
2492 : LOGICAL, OPTIONAL :: for_aux_fit
2493 : TYPE(cp_fm_type), DIMENSION(:, :, :), INTENT(IN), &
2494 : OPTIONAL :: pmat_ext
2495 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
2496 : POINTER :: overlap_rs
2497 :
2498 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_density_transform'
2499 :
2500 : INTEGER :: handle, ic, ik, ikk, indx, ispin, nc, &
2501 : nimg, nkp, nspin, omp_threads
2502 38254 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2503 : LOGICAL :: aux_fit, do_ext, do_symmetric, my_kpgrp, &
2504 : real_only
2505 38254 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
2506 38254 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
2507 38254 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:) :: info
2508 : TYPE(cp_fm_type) :: fmdummy
2509 : TYPE(dbcsr_type), POINTER :: cpmat, rpmat, scpmat, srpmat
2510 38254 : TYPE(kp_transform_plan_type) :: transform_plan
2511 : TYPE(kpoint_env_type), POINTER :: kp
2512 : TYPE(mp_para_env_type), POINTER :: para_env
2513 :
2514 38254 : CALL timeset(routineN, handle)
2515 :
2516 38254 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
2517 :
2518 38254 : IF (PRESENT(for_aux_fit)) THEN
2519 372 : aux_fit = for_aux_fit
2520 : ELSE
2521 37882 : aux_fit = .FALSE.
2522 : END IF
2523 :
2524 38254 : do_ext = .FALSE.
2525 38254 : IF (PRESENT(pmat_ext)) do_ext = .TRUE.
2526 :
2527 38254 : IF (aux_fit) THEN
2528 216 : CPASSERT(ASSOCIATED(kpoint%kp_aux_env))
2529 : END IF
2530 :
2531 : ! work storage
2532 38254 : ALLOCATE (rpmat)
2533 : CALL dbcsr_create(rpmat, template=tempmat, &
2534 38314 : matrix_type=MERGE(dbcsr_type_symmetric, dbcsr_type_no_symmetry, do_symmetric))
2535 38254 : CALL cp_dbcsr_alloc_block_from_nbl(rpmat, sab_nl)
2536 38254 : CALL dbcsr_set(rpmat, 0.0_dp)
2537 38254 : ALLOCATE (cpmat)
2538 : CALL dbcsr_create(cpmat, template=tempmat, &
2539 38314 : matrix_type=MERGE(dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, do_symmetric))
2540 38254 : CALL cp_dbcsr_alloc_block_from_nbl(cpmat, sab_nl)
2541 38254 : CALL dbcsr_set(cpmat, 0.0_dp)
2542 38254 : IF (.NOT. kpoint%full_grid) THEN
2543 35012 : ALLOCATE (srpmat)
2544 35012 : CALL dbcsr_create(srpmat, template=rpmat)
2545 35012 : CALL cp_dbcsr_alloc_block_from_nbl(srpmat, sab_nl)
2546 35012 : CALL dbcsr_set(srpmat, 0.0_dp)
2547 35012 : ALLOCATE (scpmat)
2548 35012 : CALL dbcsr_create(scpmat, template=cpmat)
2549 35012 : CALL cp_dbcsr_alloc_block_from_nbl(scpmat, sab_nl)
2550 35012 : CALL dbcsr_set(scpmat, 0.0_dp)
2551 : END IF
2552 :
2553 : CALL get_kpoint_info(kpoint, nkp=nkp, xkp=xkp, wkp=wkp, &
2554 38254 : cell_to_index=cell_to_index)
2555 38254 : IF (PRESENT(overlap_rs)) THEN
2556 37218 : CALL calibrate_symmetry_phases(kpoint, overlap_rs, tempmat, sab_nl, cell_to_index)
2557 : END IF
2558 : ! initialize real space density matrices
2559 38254 : IF (aux_fit) THEN
2560 216 : kp => kpoint%kp_aux_env(1)%kpoint_env
2561 : ELSE
2562 38038 : kp => kpoint%kp_env(1)%kpoint_env
2563 : END IF
2564 38254 : nspin = SIZE(kp%mos, 2)
2565 38254 : nc = SIZE(kp%mos, 1)
2566 38254 : nimg = SIZE(denmat, 2)
2567 38254 : real_only = (nc == 1)
2568 : CALL kp_transform_plan_create(transform_plan, sab_nl, cell_to_index, nimg, &
2569 38254 : group_entries=.TRUE.)
2570 :
2571 38254 : para_env => kpoint%blacs_env_all%para_env
2572 38254 : omp_threads = 1
2573 38254 : !$ omp_threads = omp_get_max_threads()
2574 38254 : IF (para_env%num_pe == 1 .AND. omp_threads > 1) THEN
2575 : CALL kpoint_density_transform_local(kpoint, denmat, wtype, aux_fit, do_ext, pmat_ext, &
2576 : rpmat, cpmat, srpmat, scpmat, real_only, nspin, nimg, &
2577 0 : nkp, xkp, wkp, transform_plan)
2578 0 : CALL kp_transform_plan_release(transform_plan)
2579 0 : CALL dbcsr_deallocate_matrix(rpmat)
2580 0 : CALL dbcsr_deallocate_matrix(cpmat)
2581 0 : IF (.NOT. kpoint%full_grid) THEN
2582 0 : CALL dbcsr_deallocate_matrix(srpmat)
2583 0 : CALL dbcsr_deallocate_matrix(scpmat)
2584 : END IF
2585 0 : CALL timestop(handle)
2586 0 : RETURN
2587 : END IF
2588 :
2589 598710 : ALLOCATE (info(nspin*nkp*nc))
2590 :
2591 : ! Start all the communication
2592 38254 : indx = 0
2593 78994 : DO ispin = 1, nspin
2594 2408242 : DO ic = 1, nimg
2595 2408242 : CALL dbcsr_set(denmat(ispin, ic)%matrix, 0.0_dp)
2596 : END DO
2597 : !
2598 206290 : DO ik = 1, nkp
2599 127296 : my_kpgrp = (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2))
2600 : IF (my_kpgrp) THEN
2601 82088 : ikk = ik - kpoint%kp_range(1) + 1
2602 82088 : IF (aux_fit) THEN
2603 2714 : kp => kpoint%kp_aux_env(ikk)%kpoint_env
2604 : ELSE
2605 79374 : kp => kpoint%kp_env(ikk)%kpoint_env
2606 : END IF
2607 : ELSE
2608 : NULLIFY (kp)
2609 : END IF
2610 : ! collect this density matrix on all processors
2611 127296 : CPASSERT(SIZE(fmwork) >= nc)
2612 :
2613 168036 : IF (my_kpgrp) THEN
2614 246096 : DO ic = 1, nc
2615 164008 : indx = indx + 1
2616 246096 : IF (do_ext) THEN
2617 6064 : CALL cp_fm_start_copy_general(pmat_ext(ikk, ic, ispin), fmwork(ic), para_env, info(indx))
2618 : ELSE
2619 157944 : IF (wtype) THEN
2620 3148 : CALL cp_fm_start_copy_general(kp%wmat(ic, ispin), fmwork(ic), para_env, info(indx))
2621 : ELSE
2622 154796 : CALL cp_fm_start_copy_general(kp%pmat(ic, ispin), fmwork(ic), para_env, info(indx))
2623 : END IF
2624 : END IF
2625 : END DO
2626 : ELSE
2627 135624 : DO ic = 1, nc
2628 90416 : indx = indx + 1
2629 135624 : CALL cp_fm_start_copy_general(fmdummy, fmwork(ic), para_env, info(indx))
2630 : END DO
2631 : END IF
2632 : END DO
2633 : END DO
2634 :
2635 : ! Finish communication and transform the received matrices
2636 38254 : indx = 0
2637 78994 : DO ispin = 1, nspin
2638 206290 : DO ik = 1, nkp
2639 381720 : DO ic = 1, nc
2640 254424 : indx = indx + 1
2641 381720 : CALL cp_fm_finish_copy_general(fmwork(ic), info(indx))
2642 : END DO
2643 :
2644 : ! reduce to dbcsr storage
2645 127296 : IF (real_only) THEN
2646 168 : CALL copy_fm_to_dbcsr(fmwork(1), rpmat, keep_sparsity=.TRUE.)
2647 : ELSE
2648 127128 : CALL copy_fm_to_dbcsr(fmwork(1), rpmat, keep_sparsity=.TRUE.)
2649 127128 : CALL copy_fm_to_dbcsr(fmwork(2), cpmat, keep_sparsity=.TRUE.)
2650 : END IF
2651 :
2652 : CALL kp_transform_density_matrix(kpoint, denmat, rpmat, cpmat, srpmat, scpmat, &
2653 : ispin, real_only, ik, xkp(1:3, ik), wkp(ik), &
2654 168036 : transform_plan)
2655 : END DO
2656 : END DO
2657 :
2658 : ! Clean up communication
2659 38254 : indx = 0
2660 78994 : DO ispin = 1, nspin
2661 206290 : DO ik = 1, nkp
2662 127296 : my_kpgrp = (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2))
2663 40740 : IF (my_kpgrp) THEN
2664 246096 : ikk = ik - kpoint%kp_range(1) + 1
2665 : IF (aux_fit) THEN
2666 246096 : kp => kpoint%kp_aux_env(ikk)%kpoint_env
2667 : ELSE
2668 246096 : kp => kpoint%kp_env(ikk)%kpoint_env
2669 : END IF
2670 :
2671 246096 : DO ic = 1, nc
2672 164008 : indx = indx + 1
2673 246096 : CALL cp_fm_cleanup_copy_general(info(indx))
2674 : END DO
2675 : ELSE
2676 : ! calls with dummy arguments, so not included
2677 : ! therefore just increment counter by trip count
2678 45208 : indx = indx + nc
2679 : END IF
2680 : END DO
2681 : END DO
2682 :
2683 : ! All done
2684 38254 : CALL kp_transform_plan_release(transform_plan)
2685 292678 : DEALLOCATE (info)
2686 :
2687 38254 : CALL dbcsr_deallocate_matrix(rpmat)
2688 38254 : CALL dbcsr_deallocate_matrix(cpmat)
2689 38254 : IF (.NOT. kpoint%full_grid) THEN
2690 35012 : CALL dbcsr_deallocate_matrix(srpmat)
2691 35012 : CALL dbcsr_deallocate_matrix(scpmat)
2692 : END IF
2693 :
2694 38254 : CALL timestop(handle)
2695 :
2696 76508 : END SUBROUTINE kpoint_density_transform
2697 :
2698 : ! **************************************************************************************************
2699 : !> \brief One-rank K-to-R density transform without full-matrix redistribution.
2700 : !> MPI communication remains on the generic path; this routine is called
2701 : !> only when every k-point is local to the current rank.
2702 : !> \param kpoint ...
2703 : !> \param denmat ...
2704 : !> \param wtype ...
2705 : !> \param aux_fit ...
2706 : !> \param do_ext ...
2707 : !> \param pmat_ext ...
2708 : !> \param rpmat ...
2709 : !> \param cpmat ...
2710 : !> \param srpmat ...
2711 : !> \param scpmat ...
2712 : !> \param real_only ...
2713 : !> \param nspin ...
2714 : !> \param nimg ...
2715 : !> \param nkp ...
2716 : !> \param xkp ...
2717 : !> \param wkp ...
2718 : !> \param plan ...
2719 : ! **************************************************************************************************
2720 0 : SUBROUTINE kpoint_density_transform_local(kpoint, denmat, wtype, aux_fit, do_ext, pmat_ext, &
2721 : rpmat, cpmat, srpmat, scpmat, real_only, nspin, nimg, &
2722 : nkp, xkp, wkp, plan)
2723 :
2724 : TYPE(kpoint_type), POINTER :: kpoint
2725 : TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
2726 : LOGICAL, INTENT(IN) :: wtype, aux_fit, do_ext
2727 : TYPE(cp_fm_type), DIMENSION(:, :, :), INTENT(IN), &
2728 : OPTIONAL :: pmat_ext
2729 : TYPE(dbcsr_type), POINTER :: rpmat, cpmat, srpmat, scpmat
2730 : LOGICAL, INTENT(IN) :: real_only
2731 : INTEGER, INTENT(IN) :: nspin, nimg, nkp
2732 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
2733 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
2734 : TYPE(kp_transform_plan_type), INTENT(IN) :: plan
2735 :
2736 : INTEGER :: ic, ik, ispin
2737 : TYPE(cp_fm_type), POINTER :: source_c, source_r
2738 : TYPE(kpoint_env_type), POINTER :: kp
2739 :
2740 0 : CPASSERT(kpoint%kp_range(1) == 1 .AND. kpoint%kp_range(2) == nkp)
2741 0 : CPASSERT(.NOT. do_ext .OR. PRESENT(pmat_ext))
2742 :
2743 0 : DO ispin = 1, nspin
2744 0 : DO ic = 1, nimg
2745 0 : CALL dbcsr_set(denmat(ispin, ic)%matrix, 0.0_dp)
2746 : END DO
2747 :
2748 0 : DO ik = 1, nkp
2749 0 : IF (aux_fit) THEN
2750 0 : kp => kpoint%kp_aux_env(ik)%kpoint_env
2751 : ELSE
2752 0 : kp => kpoint%kp_env(ik)%kpoint_env
2753 : END IF
2754 :
2755 0 : IF (do_ext) THEN
2756 0 : CALL kp_copy_fm_to_dbcsr_local(pmat_ext(ik, 1, ispin), rpmat)
2757 0 : IF (.NOT. real_only) CALL kp_copy_fm_to_dbcsr_local(pmat_ext(ik, 2, ispin), cpmat)
2758 : ELSE
2759 0 : IF (wtype) THEN
2760 0 : source_r => kp%wmat(1, ispin)
2761 0 : IF (.NOT. real_only) source_c => kp%wmat(2, ispin)
2762 : ELSE
2763 0 : source_r => kp%pmat(1, ispin)
2764 0 : IF (.NOT. real_only) source_c => kp%pmat(2, ispin)
2765 : END IF
2766 0 : CALL kp_copy_fm_to_dbcsr_local(source_r, rpmat)
2767 0 : IF (.NOT. real_only) CALL kp_copy_fm_to_dbcsr_local(source_c, cpmat)
2768 : END IF
2769 :
2770 : CALL kp_transform_density_matrix(kpoint, denmat, rpmat, cpmat, srpmat, scpmat, &
2771 0 : ispin, real_only, ik, xkp(1:3, ik), wkp(ik), plan)
2772 : END DO
2773 : END DO
2774 :
2775 0 : END SUBROUTINE kpoint_density_transform_local
2776 :
2777 : ! **************************************************************************************************
2778 : !> \brief Apply k-point symmetry and accumulate one reciprocal-space density matrix in real space.
2779 : !> \param kpoint ...
2780 : !> \param denmat ...
2781 : !> \param rpmat ...
2782 : !> \param cpmat ...
2783 : !> \param srpmat ...
2784 : !> \param scpmat ...
2785 : !> \param ispin ...
2786 : !> \param real_only ...
2787 : !> \param ik ...
2788 : !> \param xkp ...
2789 : !> \param wkp ...
2790 : !> \param plan ...
2791 : ! **************************************************************************************************
2792 127296 : SUBROUTINE kp_transform_density_matrix(kpoint, denmat, rpmat, cpmat, srpmat, scpmat, &
2793 : ispin, real_only, ik, xkp, wkp, plan)
2794 :
2795 : TYPE(kpoint_type), POINTER :: kpoint
2796 : TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
2797 : TYPE(dbcsr_type), POINTER :: rpmat, cpmat, srpmat, scpmat
2798 : INTEGER, INTENT(IN) :: ispin
2799 : LOGICAL, INTENT(IN) :: real_only
2800 : INTEGER, INTENT(IN) :: ik
2801 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
2802 : REAL(KIND=dp), INTENT(IN) :: wkp
2803 : TYPE(kp_transform_plan_type), INTENT(IN) :: plan
2804 :
2805 : INTEGER :: ir, ira, is, jr
2806 : LOGICAL :: reverse_phase
2807 : REAL(KIND=dp) :: symmetry_weight
2808 127296 : TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kind_rot
2809 : TYPE(kpoint_sym_type), POINTER :: kpsym
2810 :
2811 127296 : kpsym => kpoint%kp_sym(ik)%kpoint_sym
2812 0 : CPASSERT(ASSOCIATED(kpsym))
2813 127296 : IF (kpsym%apply_symmetry) THEN
2814 33070 : symmetry_weight = wkp/REAL(kpsym%nwght, KIND=dp)
2815 170518 : DO is = 1, kpsym%nwght
2816 137448 : ir = ABS(kpsym%rotp(is))
2817 137448 : ira = 0
2818 3538096 : DO jr = 1, SIZE(kpoint%ibrot)
2819 3538096 : IF (ir == kpoint%ibrot(jr)) ira = jr
2820 : END DO
2821 137448 : CPASSERT(ira > 0)
2822 137448 : kind_rot => kpoint%kind_rotmat(ira, :)
2823 137448 : CPASSERT(kpsym%phase_mode(is) > 0)
2824 137448 : reverse_phase = kpsym%phase_mode(is) == 2
2825 : CALL symtrans_phase(srpmat, scpmat, rpmat, cpmat, real_only, kind_rot, &
2826 : kpsym%rot(1:3, 1:3, is), kpsym%f0(:, is), &
2827 : kpsym%fcell_gauge(:, :, is), kpoint%atype, &
2828 137448 : kpsym%xkp(1:3, is), kpsym%rotp(is) < 0, reverse_phase)
2829 : CALL transform_dmat(denmat, srpmat, scpmat, ispin, real_only, plan, &
2830 170518 : kpsym%xkp(1:3, is), symmetry_weight)
2831 : END DO
2832 : ELSE
2833 94226 : CALL transform_dmat(denmat, rpmat, cpmat, ispin, real_only, plan, xkp, wkp)
2834 : END IF
2835 :
2836 127296 : END SUBROUTINE kp_transform_density_matrix
2837 :
2838 : ! **************************************************************************************************
2839 : !> \brief Copy a one-rank full matrix directly into existing DBCSR blocks.
2840 : !> \param fm ...
2841 : !> \param matrix ...
2842 : ! **************************************************************************************************
2843 0 : SUBROUTINE kp_copy_fm_to_dbcsr_local(fm, matrix)
2844 :
2845 : TYPE(cp_fm_type), INTENT(IN) :: fm
2846 : TYPE(dbcsr_type), POINTER :: matrix
2847 :
2848 : INTEGER :: col_offset, row_offset
2849 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2850 0 : POINTER :: full
2851 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
2852 : TYPE(dbcsr_iterator_type) :: iterator
2853 :
2854 0 : CALL cp_fm_get_info(fm, local_data=full)
2855 0 : CALL dbcsr_iterator_start(iterator, matrix, shared=.FALSE.)
2856 0 : DO WHILE (dbcsr_iterator_blocks_left(iterator))
2857 : CALL dbcsr_iterator_next_block(iterator, block=block, row_offset=row_offset, &
2858 0 : col_offset=col_offset)
2859 : block(:, :) = full(row_offset:row_offset + SIZE(block, 1) - 1, &
2860 0 : col_offset:col_offset + SIZE(block, 2) - 1)
2861 : END DO
2862 0 : CALL dbcsr_iterator_stop(iterator)
2863 :
2864 0 : END SUBROUTINE kp_copy_fm_to_dbcsr_local
2865 :
2866 : ! **************************************************************************************************
2867 : !> \brief Build the immutable neighbor-list traversal shared by R-to-K and K-to-R transforms.
2868 : !> \param plan ...
2869 : !> \param sab_nl ...
2870 : !> \param cell_to_index ...
2871 : !> \param nimg ...
2872 : !> \param block_template ...
2873 : !> \param group_entries ...
2874 : ! **************************************************************************************************
2875 38254 : SUBROUTINE kp_transform_plan_create(plan, sab_nl, cell_to_index, nimg, block_template, group_entries)
2876 :
2877 : TYPE(kp_transform_plan_type), INTENT(OUT) :: plan
2878 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2879 : POINTER :: sab_nl
2880 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2881 : INTEGER, INTENT(IN) :: nimg
2882 : TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: block_template
2883 : LOGICAL, INTENT(IN), OPTIONAL :: group_entries
2884 :
2885 : INTEGER :: i, iatom, icell, icol, igroup, irow, &
2886 : jatom, nblock
2887 38254 : INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: key
2888 38254 : INTEGER, ALLOCATABLE, DIMENSION(:) :: order
2889 : INTEGER, DIMENSION(3) :: cell
2890 38254 : INTEGER, DIMENSION(:), POINTER :: col_offsets, row_offsets
2891 : LOGICAL :: do_grouping, do_symmetric, store_offsets
2892 : TYPE(neighbor_list_iterator_p_type), &
2893 38254 : DIMENSION(:), POINTER :: iterator
2894 :
2895 38254 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
2896 38254 : plan%symmetric = do_symmetric
2897 38254 : do_grouping = .FALSE.
2898 38254 : IF (PRESENT(group_entries)) do_grouping = group_entries
2899 38254 : store_offsets = PRESENT(block_template)
2900 38254 : IF (store_offsets) THEN
2901 0 : CALL dbcsr_get_info(block_template, row_blk_offset=row_offsets, col_blk_offset=col_offsets)
2902 : END IF
2903 38254 : plan%nentry = 0
2904 38254 : plan%ngroup = 0
2905 38254 : CALL neighbor_list_iterator_create(iterator, sab_nl)
2906 9130919 : DO WHILE (neighbor_list_iterate(iterator) == 0)
2907 9092665 : CALL get_iterator_info(iterator, cell=cell)
2908 9092665 : icell = cell_to_index(cell(1), cell(2), cell(3))
2909 9130919 : IF (icell >= 1 .AND. icell <= nimg) plan%nentry = plan%nentry + 1
2910 : END DO
2911 38254 : CALL neighbor_list_iterator_release(iterator)
2912 :
2913 0 : ALLOCATE (plan%row(plan%nentry), plan%col(plan%nentry), plan%image(plan%nentry), &
2914 339371 : plan%cell(3, plan%nentry), plan%symmetry_sign(plan%nentry))
2915 38254 : IF (store_offsets) THEN
2916 0 : ALLOCATE (plan%row_offset(plan%nentry), plan%col_offset(plan%nentry))
2917 : END IF
2918 :
2919 38254 : i = 0
2920 38254 : CALL neighbor_list_iterator_create(iterator, sab_nl)
2921 9130919 : DO WHILE (neighbor_list_iterate(iterator) == 0)
2922 9092665 : CALL get_iterator_info(iterator, iatom=iatom, jatom=jatom, cell=cell)
2923 9092665 : icell = cell_to_index(cell(1), cell(2), cell(3))
2924 9092665 : IF (icell < 1 .OR. icell > nimg) CYCLE
2925 9092577 : i = i + 1
2926 9092577 : irow = iatom
2927 9092577 : icol = jatom
2928 9092577 : plan%symmetry_sign(i) = 1.0_dp
2929 9092577 : IF (do_symmetric .AND. iatom > jatom) THEN
2930 3842507 : irow = jatom
2931 3842507 : icol = iatom
2932 3842507 : plan%symmetry_sign(i) = -1.0_dp
2933 : END IF
2934 9092577 : plan%row(i) = irow
2935 9092577 : plan%col(i) = icol
2936 9092577 : IF (store_offsets) THEN
2937 0 : plan%row_offset(i) = row_offsets(irow)
2938 0 : plan%col_offset(i) = col_offsets(icol)
2939 : END IF
2940 9092577 : plan%image(i) = icell
2941 36408650 : plan%cell(:, i) = cell
2942 : END DO
2943 38254 : CALL neighbor_list_iterator_release(iterator)
2944 38254 : CPASSERT(i == plan%nentry)
2945 :
2946 38254 : IF (do_grouping .AND. plan%nentry > 0) THEN
2947 18222425 : nblock = MAX(MAXVAL(plan%row), MAXVAL(plan%col))
2948 186355 : ALLOCATE (key(plan%nentry), order(plan%nentry))
2949 9129848 : DO i = 1, plan%nentry
2950 : key(i) = INT(plan%image(i), KIND=int_8) + INT(nimg, KIND=int_8)* &
2951 : (INT(plan%row(i) - 1, KIND=int_8) + INT(nblock, KIND=int_8)* &
2952 9129848 : INT(plan%col(i) - 1, KIND=int_8))
2953 : END DO
2954 37271 : CALL sort(key, plan%nentry, order)
2955 18259696 : plan%row(:) = plan%row(order)
2956 18259696 : plan%col(:) = plan%col(order)
2957 37271 : IF (store_offsets) THEN
2958 0 : plan%row_offset(:) = plan%row_offset(order)
2959 0 : plan%col_offset(:) = plan%col_offset(order)
2960 : END IF
2961 18259696 : plan%image(:) = plan%image(order)
2962 72815158 : plan%cell(:, :) = plan%cell(:, order)
2963 18259696 : plan%symmetry_sign(:) = plan%symmetry_sign(order)
2964 :
2965 37271 : plan%ngroup = 1
2966 9092577 : DO i = 2, plan%nentry
2967 9092577 : IF (key(i) /= key(i - 1)) plan%ngroup = plan%ngroup + 1
2968 : END DO
2969 111813 : ALLOCATE (plan%group_start(plan%ngroup + 1))
2970 37271 : igroup = 1
2971 37271 : plan%group_start(igroup) = 1
2972 9092577 : DO i = 2, plan%nentry
2973 9092577 : IF (key(i) /= key(i - 1)) THEN
2974 9055306 : igroup = igroup + 1
2975 9055306 : plan%group_start(igroup) = i
2976 : END IF
2977 : END DO
2978 37271 : plan%group_start(plan%ngroup + 1) = plan%nentry + 1
2979 37271 : DEALLOCATE (key, order)
2980 : END IF
2981 :
2982 76508 : END SUBROUTINE kp_transform_plan_create
2983 :
2984 : ! **************************************************************************************************
2985 : !> \brief Release a K-to-R traversal plan.
2986 : !> \param plan ...
2987 : ! **************************************************************************************************
2988 38254 : SUBROUTINE kp_transform_plan_release(plan)
2989 :
2990 : TYPE(kp_transform_plan_type), INTENT(INOUT) :: plan
2991 :
2992 38254 : DEALLOCATE (plan%row, plan%col, plan%image, plan%cell, plan%symmetry_sign)
2993 38254 : IF (ALLOCATED(plan%row_offset)) DEALLOCATE (plan%row_offset, plan%col_offset)
2994 38254 : IF (ALLOCATED(plan%group_start)) DEALLOCATE (plan%group_start)
2995 38254 : plan%nentry = 0
2996 38254 : plan%ngroup = 0
2997 38254 : plan%symmetric = .FALSE.
2998 :
2999 38254 : END SUBROUTINE kp_transform_plan_release
3000 :
3001 : ! **************************************************************************************************
3002 : !> \brief real space density matrices in DBCSR format
3003 : !> \param denmat Real space (DBCSR) density matrix
3004 : !> \param rpmat ...
3005 : !> \param cpmat ...
3006 : !> \param ispin ...
3007 : !> \param real_only ...
3008 : !> \param plan precomputed neighbor-list traversal
3009 : !> \param xkp ...
3010 : !> \param wkp ...
3011 : ! **************************************************************************************************
3012 231674 : SUBROUTINE transform_dmat(denmat, rpmat, cpmat, ispin, real_only, plan, xkp, wkp)
3013 :
3014 : TYPE(dbcsr_p_type), DIMENSION(:, :) :: denmat
3015 : TYPE(dbcsr_type), POINTER :: rpmat, cpmat
3016 : INTEGER, INTENT(IN) :: ispin
3017 : LOGICAL, INTENT(IN) :: real_only
3018 : TYPE(kp_transform_plan_type), INTENT(IN) :: plan
3019 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
3020 : REAL(KIND=dp), INTENT(IN) :: wkp
3021 :
3022 : CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_dmat'
3023 :
3024 : INTEGER :: handle, i, ifirst, igroup, ilast, &
3025 : nthreads
3026 : LOGICAL :: found
3027 : REAL(KIND=dp) :: arg, coskl, sinkl
3028 231674 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, dblock, rblock
3029 :
3030 231674 : CALL timeset(routineN, handle)
3031 :
3032 231674 : nthreads = 1
3033 231674 : !$ nthreads = MAX(1, MIN(omp_get_max_threads(), plan%ngroup))
3034 : ! Entries are grouped by (image,row-block,col-block), hence workers
3035 : ! update disjoint DBCSR blocks without atomics or false sharing.
3036 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(DYNAMIC, 4) NUM_THREADS(nthreads) &
3037 : !$OMP SHARED(plan,denmat,rpmat,cpmat,ispin,real_only,xkp,wkp,nthreads) &
3038 231674 : !$OMP PRIVATE(igroup,ifirst,ilast,i,arg,coskl,sinkl,found,dblock,rblock,cblock)
3039 : DO igroup = 1, plan%ngroup
3040 : ifirst = plan%group_start(igroup)
3041 : ilast = plan%group_start(igroup + 1) - 1
3042 :
3043 : !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
3044 : !Therefore, we have: S(R) = sum_k Re(S(k))*cos(k*R) -i^2*Im(S(k))*sin(k*R)
3045 : ! = sum_k Re(S(k))*cos(k*R) + Im(S(k))*sin(k*R)
3046 : !fc = +- 1 is due to the usual non-symmetric real-space matrices stored as symmetric ones
3047 : coskl = 0.0_dp
3048 : sinkl = 0.0_dp
3049 : DO i = ifirst, ilast
3050 : arg = DOT_PRODUCT(REAL(plan%cell(:, i), KIND=dp), xkp)
3051 : coskl = coskl + wkp*COS(twopi*arg)
3052 : sinkl = sinkl + wkp*plan%symmetry_sign(i)*SIN(twopi*arg)
3053 : END DO
3054 :
3055 : CALL dbcsr_get_block_p(matrix=denmat(ispin, plan%image(ifirst))%matrix, &
3056 : row=plan%row(ifirst), col=plan%col(ifirst), &
3057 : block=dblock, found=found)
3058 : IF (.NOT. found) CYCLE
3059 :
3060 : IF (real_only) THEN
3061 : CALL dbcsr_get_readonly_block_p(matrix=rpmat, row=plan%row(ifirst), col=plan%col(ifirst), &
3062 : block=rblock, found=found)
3063 : IF (.NOT. found) CYCLE
3064 : dblock = dblock + coskl*rblock
3065 : ELSE
3066 : CALL dbcsr_get_readonly_block_p(matrix=rpmat, row=plan%row(ifirst), col=plan%col(ifirst), &
3067 : block=rblock, found=found)
3068 : IF (.NOT. found) CYCLE
3069 : CALL dbcsr_get_readonly_block_p(matrix=cpmat, row=plan%row(ifirst), col=plan%col(ifirst), &
3070 : block=cblock, found=found)
3071 : IF (.NOT. found) CYCLE
3072 : dblock = dblock + coskl*rblock
3073 : dblock = dblock + sinkl*cblock
3074 : END IF
3075 : END DO
3076 : !$OMP END PARALLEL DO
3077 :
3078 231674 : CALL timestop(handle)
3079 :
3080 231674 : END SUBROUTINE transform_dmat
3081 :
3082 : ! **************************************************************************************************
3083 : !> \brief Allocate a dense work matrix with the requested shape
3084 : !> \param work dense work matrix
3085 : !> \param nrow number of rows
3086 : !> \param ncol number of columns
3087 : ! **************************************************************************************************
3088 3066398 : SUBROUTINE ensure_work_matrix(work, nrow, ncol)
3089 :
3090 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
3091 : INTENT(INOUT) :: work
3092 : INTEGER, INTENT(IN) :: nrow, ncol
3093 :
3094 3066398 : IF (ALLOCATED(work)) THEN
3095 2955155 : IF (SIZE(work, 1) == nrow .AND. SIZE(work, 2) == ncol) RETURN
3096 47397 : DEALLOCATE (work)
3097 : END IF
3098 634560 : ALLOCATE (work(nrow, ncol))
3099 :
3100 : END SUBROUTINE ensure_work_matrix
3101 :
3102 : ! **************************************************************************************************
3103 : !> \brief Select the Bloch-phase convention that preserves overlap covariance.
3104 : !> \param kpoint ...
3105 : !> \param overlap_rs ...
3106 : !> \param tempmat ...
3107 : !> \param sab_nl ...
3108 : !> \param cell_to_index ...
3109 : ! **************************************************************************************************
3110 37218 : SUBROUTINE calibrate_symmetry_phases(kpoint, overlap_rs, tempmat, sab_nl, cell_to_index)
3111 :
3112 : TYPE(kpoint_type), POINTER :: kpoint
3113 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: overlap_rs
3114 : TYPE(dbcsr_type), POINTER :: tempmat
3115 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3116 : POINTER :: sab_nl
3117 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3118 :
3119 : CHARACTER(LEN=256) :: phase_error
3120 : INTEGER :: best_mode, ik, ir, ira, is, jr, mode
3121 : LOGICAL :: needs_calibration, reverse
3122 : REAL(KIND=dp) :: best_residual, candidate_norm, &
3123 : direct_norm, overlap_dot, &
3124 : phase_tolerance, relative_residual
3125 : TYPE(dbcsr_type), POINTER :: direct_c, direct_r, source_c, source_r, &
3126 : sym_c, sym_r
3127 37218 : TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kind_rot
3128 : TYPE(kpoint_sym_type), POINTER :: kpsym
3129 :
3130 37218 : needs_calibration = .FALSE.
3131 127178 : DO ik = 1, kpoint%nkp
3132 97858 : kpsym => kpoint%kp_sym(ik)%kpoint_sym
3133 127178 : IF (kpsym%apply_symmetry) THEN
3134 92236 : IF (ANY(kpsym%phase_mode == 0)) THEN
3135 : needs_calibration = .TRUE.
3136 : EXIT
3137 : END IF
3138 : END IF
3139 : END DO
3140 37218 : IF (.NOT. needs_calibration) RETURN
3141 :
3142 7898 : ALLOCATE (source_r, source_c, direct_r, direct_c, sym_r, sym_c)
3143 7898 : CALL dbcsr_create(source_r, template=tempmat, matrix_type=dbcsr_type_symmetric)
3144 7898 : CALL dbcsr_create(source_c, template=tempmat, matrix_type=dbcsr_type_antisymmetric)
3145 7898 : CALL dbcsr_create(direct_r, template=tempmat, matrix_type=dbcsr_type_symmetric)
3146 7898 : CALL dbcsr_create(direct_c, template=tempmat, matrix_type=dbcsr_type_antisymmetric)
3147 7898 : CALL dbcsr_create(sym_r, template=tempmat, matrix_type=dbcsr_type_symmetric)
3148 7898 : CALL dbcsr_create(sym_c, template=tempmat, matrix_type=dbcsr_type_antisymmetric)
3149 7898 : CALL cp_dbcsr_alloc_block_from_nbl(source_r, sab_nl)
3150 7898 : CALL cp_dbcsr_alloc_block_from_nbl(source_c, sab_nl)
3151 7898 : CALL cp_dbcsr_alloc_block_from_nbl(direct_r, sab_nl)
3152 7898 : CALL cp_dbcsr_alloc_block_from_nbl(direct_c, sab_nl)
3153 7898 : CALL cp_dbcsr_alloc_block_from_nbl(sym_r, sab_nl)
3154 7898 : CALL cp_dbcsr_alloc_block_from_nbl(sym_c, sab_nl)
3155 :
3156 7898 : phase_tolerance = MAX(1.0e-6_dp, 100.0_dp*kpoint%eps_geo)
3157 28380 : DO ik = 1, kpoint%nkp
3158 20482 : kpsym => kpoint%kp_sym(ik)%kpoint_sym
3159 20482 : IF (.NOT. kpsym%apply_symmetry) CYCLE
3160 97480 : IF (ALL(kpsym%phase_mode > 0)) CYCLE
3161 :
3162 15652 : CALL dbcsr_set(source_r, 0.0_dp)
3163 15652 : CALL dbcsr_set(source_c, 0.0_dp)
3164 : CALL rskp_transform(source_r, source_c, overlap_rs, 1, kpoint%xkp(1:3, ik), &
3165 15652 : cell_to_index, sab_nl)
3166 :
3167 114114 : DO is = 1, kpsym%nwght
3168 90564 : IF (kpsym%phase_mode(is) > 0) CYCLE
3169 12424 : CALL dbcsr_set(direct_r, 0.0_dp)
3170 12424 : CALL dbcsr_set(direct_c, 0.0_dp)
3171 : CALL rskp_transform(direct_r, direct_c, overlap_rs, 1, kpsym%xkp(1:3, is), &
3172 12424 : cell_to_index, sab_nl)
3173 12424 : CALL dbcsr_dot(direct_r, direct_r, direct_norm)
3174 12424 : CALL dbcsr_dot(direct_c, direct_c, candidate_norm)
3175 12424 : direct_norm = direct_norm + candidate_norm
3176 :
3177 12424 : ir = ABS(kpsym%rotp(is))
3178 12424 : ira = 0
3179 310008 : DO jr = 1, SIZE(kpoint%ibrot)
3180 310008 : IF (ir == kpoint%ibrot(jr)) ira = jr
3181 : END DO
3182 12424 : CPASSERT(ira > 0)
3183 12424 : kind_rot => kpoint%kind_rotmat(ira, :)
3184 :
3185 12424 : best_mode = 0
3186 12424 : best_residual = HUGE(1.0_dp)
3187 37272 : DO mode = 1, 2
3188 24848 : reverse = mode == 2
3189 : CALL symtrans_phase(sym_r, sym_c, source_r, source_c, .FALSE., kind_rot, &
3190 : kpsym%rot(1:3, 1:3, is), kpsym%f0(:, is), &
3191 : kpsym%fcell_gauge(:, :, is), kpoint%atype, &
3192 24848 : kpsym%xkp(1:3, is), kpsym%rotp(is) < 0, reverse)
3193 24848 : CALL dbcsr_dot(sym_r, sym_r, candidate_norm)
3194 24848 : CALL dbcsr_dot(sym_c, sym_c, relative_residual)
3195 24848 : candidate_norm = candidate_norm + relative_residual
3196 24848 : CALL dbcsr_dot(sym_r, direct_r, overlap_dot)
3197 24848 : CALL dbcsr_dot(sym_c, direct_c, relative_residual)
3198 24848 : overlap_dot = overlap_dot + relative_residual
3199 : relative_residual = SQRT(MAX(0.0_dp, candidate_norm + direct_norm - &
3200 24848 : 2.0_dp*overlap_dot)/MAX(direct_norm, TINY(1.0_dp)))
3201 37272 : IF (relative_residual < best_residual) THEN
3202 16928 : best_residual = relative_residual
3203 16928 : best_mode = mode
3204 : END IF
3205 : END DO
3206 12424 : IF (best_residual > phase_tolerance) THEN
3207 : WRITE (phase_error, '(A,ES12.4,A,I0,A,I0)') &
3208 0 : "No Bloch-phase direction preserves overlap covariance; residual=", &
3209 0 : best_residual, ", irreducible k-point=", ik, ", operation=", is
3210 0 : CALL cp_abort(__LOCATION__, TRIM(phase_error))
3211 : END IF
3212 123470 : kpsym%phase_mode(is) = best_mode
3213 : END DO
3214 : END DO
3215 :
3216 7898 : CALL dbcsr_deallocate_matrix(source_r)
3217 7898 : CALL dbcsr_deallocate_matrix(source_c)
3218 7898 : CALL dbcsr_deallocate_matrix(direct_r)
3219 7898 : CALL dbcsr_deallocate_matrix(direct_c)
3220 7898 : CALL dbcsr_deallocate_matrix(sym_r)
3221 7898 : CALL dbcsr_deallocate_matrix(sym_c)
3222 :
3223 37218 : END SUBROUTINE calibrate_symmetry_phases
3224 :
3225 : ! **************************************************************************************************
3226 : !> \brief Symmetrize a complex k-point matrix including Bloch phase shifts
3227 : !> \param srpmat real part of transformed matrix
3228 : !> \param scpmat imaginary part of transformed matrix
3229 : !> \param rpmat real part of reference matrix
3230 : !> \param cpmat imaginary part of reference matrix
3231 : !> \param real_only ...
3232 : !> \param kmat kind type rotation matrix
3233 : !> \param rot rotation matrix
3234 : !> \param f0 atom permutation
3235 : !> \param fcell atom cell shifts generated by the symmetry operation
3236 : !> \param atype atom to kind pointer
3237 : !> \param xkp target k-point coordinates
3238 : !> \param time_reversal ...
3239 : !> \param reverse_phase ...
3240 : ! **************************************************************************************************
3241 162296 : SUBROUTINE symtrans_phase(srpmat, scpmat, rpmat, cpmat, real_only, kmat, rot, f0, fcell, atype, &
3242 : xkp, time_reversal, reverse_phase)
3243 :
3244 : TYPE(dbcsr_type), POINTER :: srpmat, scpmat, rpmat, cpmat
3245 : LOGICAL, INTENT(IN) :: real_only
3246 : TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kmat
3247 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: rot
3248 : INTEGER, DIMENSION(:), INTENT(IN) :: f0
3249 : INTEGER, DIMENSION(:, :), INTENT(IN) :: fcell
3250 : INTEGER, DIMENSION(:), INTENT(IN) :: atype
3251 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
3252 : LOGICAL, INTENT(IN) :: time_reversal, reverse_phase
3253 :
3254 : CHARACTER(LEN=*), PARAMETER :: routineN = 'symtrans_phase'
3255 :
3256 : INTEGER :: handle, iatom, icol, ikind, ip, irow, &
3257 : jcol, jkind, jp, jrow, mynode, natom, &
3258 : nthreads, numnodes, owner
3259 : INTEGER, DIMENSION(3) :: shift
3260 : LOGICAL :: byrows, dorot, found, has_phase, perm, &
3261 : trans
3262 : REAL(KIND=dp) :: arg, coskl, dr, sinkl
3263 162296 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cwork, rwork, twork
3264 162296 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, kroti, krotj, rblock, scblock, &
3265 162296 : srblock
3266 : TYPE(dbcsr_distribution_type) :: dist
3267 : TYPE(dbcsr_iterator_type) :: iter
3268 :
3269 162296 : CALL timeset(routineN, handle)
3270 :
3271 162296 : natom = SIZE(f0)
3272 162296 : perm = .FALSE.
3273 653304 : DO iatom = 1, natom
3274 549400 : IF (f0(iatom) == iatom) CYCLE
3275 : perm = .TRUE.
3276 594912 : EXIT
3277 : END DO
3278 :
3279 162296 : dorot = .FALSE.
3280 2109848 : IF (ABS(SUM(ABS(rot)) - 3.0_dp) > 1.e-12_dp) dorot = .TRUE.
3281 162296 : dr = ABS(rot(1, 1) - 1.0_dp) + ABS(rot(2, 2) - 1.0_dp) + ABS(rot(3, 3) - 1.0_dp)
3282 162296 : IF (ABS(dr) > 1.e-12_dp) dorot = .TRUE.
3283 2218444 : has_phase = ANY(fcell /= 0) .OR. time_reversal
3284 :
3285 162296 : IF (.NOT. (dorot .OR. perm .OR. has_phase)) THEN
3286 38494 : CALL dbcsr_copy(srpmat, rpmat)
3287 38494 : IF (.NOT. real_only) CALL dbcsr_copy(scpmat, cpmat)
3288 38494 : CALL timestop(handle)
3289 38494 : RETURN
3290 : END IF
3291 :
3292 123802 : CALL dbcsr_get_info(rpmat, distribution=dist)
3293 123802 : CALL dbcsr_distribution_get(dist, mynode=mynode, numnodes=numnodes)
3294 123802 : IF (numnodes /= 1 .AND. (perm .OR. has_phase)) THEN
3295 116958 : CALL dbcsr_replicate_all(rpmat)
3296 116958 : IF (.NOT. real_only) CALL dbcsr_replicate_all(cpmat)
3297 : END IF
3298 :
3299 123802 : CALL dbcsr_set(srpmat, 0.0_dp)
3300 123802 : IF (.NOT. real_only) CALL dbcsr_set(scpmat, 0.0_dp)
3301 :
3302 123802 : nthreads = 1
3303 123802 : !$ nthreads = omp_get_max_threads()
3304 : ! Atom permutations preserve row ownership, while block scheduling can race on a target row.
3305 123802 : byrows = .TRUE.
3306 : !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads) &
3307 : !$OMP SHARED(rpmat,cpmat,srpmat,scpmat,real_only,kmat,rot,f0,fcell,atype,xkp,time_reversal, &
3308 : !$OMP reverse_phase,natom,perm,dorot,has_phase,mynode,numnodes,nthreads,byrows) &
3309 : !$OMP PRIVATE(iter,irow,icol,rblock,rwork,cwork,twork,ikind,jkind,kroti,krotj,shift,arg,coskl, &
3310 : !$OMP sinkl,cblock,found,ip,jp,jrow,jcol,trans,srblock,scblock,owner)
3311 : CALL dbcsr_iterator_readonly_start(iter, rpmat, dynamic=.TRUE., dynamic_byrows=byrows)
3312 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3313 : CALL dbcsr_iterator_next_block(iter, irow, icol, rblock)
3314 : IF (.NOT. ALLOCATED(rwork)) THEN
3315 : ALLOCATE (rwork(SIZE(rblock, 1), SIZE(rblock, 2)))
3316 : ELSE IF (SIZE(rwork, 1) /= SIZE(rblock, 1) .OR. SIZE(rwork, 2) /= SIZE(rblock, 2)) THEN
3317 : DEALLOCATE (rwork)
3318 : ALLOCATE (rwork(SIZE(rblock, 1), SIZE(rblock, 2)))
3319 : END IF
3320 : IF (.NOT. real_only) THEN
3321 : IF (.NOT. ALLOCATED(cwork)) THEN
3322 : ALLOCATE (cwork(SIZE(rblock, 1), SIZE(rblock, 2)))
3323 : ELSE IF (SIZE(cwork, 1) /= SIZE(rblock, 1) .OR. SIZE(cwork, 2) /= SIZE(rblock, 2)) THEN
3324 : DEALLOCATE (cwork)
3325 : ALLOCATE (cwork(SIZE(rblock, 1), SIZE(rblock, 2)))
3326 : END IF
3327 : END IF
3328 :
3329 : ikind = atype(irow)
3330 : jkind = atype(icol)
3331 : kroti => kmat(ikind)%rmat
3332 : krotj => kmat(jkind)%rmat
3333 :
3334 : IF (reverse_phase) THEN
3335 : shift = fcell(1:3, irow) - fcell(1:3, icol)
3336 : ELSE
3337 : shift = fcell(1:3, icol) - fcell(1:3, irow)
3338 : END IF
3339 : arg = REAL(shift(1), dp)*xkp(1) + REAL(shift(2), dp)*xkp(2) + &
3340 : REAL(shift(3), dp)*xkp(3)
3341 : coskl = COS(twopi*arg)
3342 : sinkl = SIN(twopi*arg)
3343 : IF (real_only) THEN
3344 : IF (ABS(sinkl) > 1.e-12_dp) THEN
3345 : CALL cp_abort(__LOCATION__, "Real k-point wavefunctions cannot represent symmetry phases")
3346 : END IF
3347 : rwork(:, :) = coskl*rblock
3348 : ELSE
3349 : CALL dbcsr_get_readonly_block_p(matrix=cpmat, row=irow, col=icol, block=cblock, found=found)
3350 : rwork(:, :) = coskl*rblock
3351 : IF (time_reversal) THEN
3352 : cwork(:, :) = -sinkl*rblock
3353 : IF (found) THEN
3354 : rwork(:, :) = rwork - sinkl*cblock
3355 : cwork(:, :) = cwork - coskl*cblock
3356 : END IF
3357 : ELSE
3358 : cwork(:, :) = -sinkl*rblock
3359 : IF (found) THEN
3360 : rwork(:, :) = rwork + sinkl*cblock
3361 : cwork(:, :) = cwork + coskl*cblock
3362 : END IF
3363 : END IF
3364 : END IF
3365 :
3366 : ip = f0(irow)
3367 : jp = f0(icol)
3368 : IF (ip <= jp) THEN
3369 : jrow = ip
3370 : jcol = jp
3371 : trans = .FALSE.
3372 : ELSE
3373 : jrow = jp
3374 : jcol = ip
3375 : trans = .TRUE.
3376 : END IF
3377 :
3378 : CALL dbcsr_get_block_p(matrix=srpmat, row=jrow, col=jcol, block=srblock, found=found)
3379 : IF (.NOT. found) THEN
3380 : CALL dbcsr_get_stored_coordinates(srpmat, jrow, jcol, owner)
3381 : CPASSERT(owner /= mynode)
3382 : CYCLE
3383 : END IF
3384 : IF (trans) THEN
3385 : CALL ensure_work_matrix(twork, SIZE(krotj, 1), SIZE(rwork, 1))
3386 : CALL dgemm('N', 'T', SIZE(krotj, 1), SIZE(rwork, 1), SIZE(krotj, 2), &
3387 : 1.0_dp, krotj, SIZE(krotj, 1), rwork, SIZE(rwork, 1), &
3388 : 0.0_dp, twork, SIZE(twork, 1))
3389 : CALL dgemm('N', 'T', SIZE(twork, 1), SIZE(kroti, 1), SIZE(twork, 2), &
3390 : 1.0_dp, twork, SIZE(twork, 1), kroti, SIZE(kroti, 1), &
3391 : 1.0_dp, srblock, SIZE(srblock, 1))
3392 : ELSE
3393 : CALL ensure_work_matrix(twork, SIZE(kroti, 1), SIZE(rwork, 2))
3394 : CALL dgemm('N', 'N', SIZE(kroti, 1), SIZE(rwork, 2), SIZE(kroti, 2), &
3395 : 1.0_dp, kroti, SIZE(kroti, 1), rwork, SIZE(rwork, 1), &
3396 : 0.0_dp, twork, SIZE(twork, 1))
3397 : CALL dgemm('N', 'T', SIZE(twork, 1), SIZE(krotj, 1), SIZE(twork, 2), &
3398 : 1.0_dp, twork, SIZE(twork, 1), krotj, SIZE(krotj, 1), &
3399 : 1.0_dp, srblock, SIZE(srblock, 1))
3400 : END IF
3401 :
3402 : IF (.NOT. real_only) THEN
3403 : CALL dbcsr_get_block_p(matrix=scpmat, row=jrow, col=jcol, block=scblock, found=found)
3404 : CPASSERT(found)
3405 : IF (trans) THEN
3406 : CALL ensure_work_matrix(twork, SIZE(krotj, 1), SIZE(cwork, 1))
3407 : CALL dgemm('N', 'T', SIZE(krotj, 1), SIZE(cwork, 1), SIZE(krotj, 2), &
3408 : 1.0_dp, krotj, SIZE(krotj, 1), cwork, SIZE(cwork, 1), &
3409 : 0.0_dp, twork, SIZE(twork, 1))
3410 : CALL dgemm('N', 'T', SIZE(twork, 1), SIZE(kroti, 1), SIZE(twork, 2), &
3411 : -1.0_dp, twork, SIZE(twork, 1), kroti, SIZE(kroti, 1), &
3412 : 1.0_dp, scblock, SIZE(scblock, 1))
3413 : ELSE
3414 : CALL ensure_work_matrix(twork, SIZE(kroti, 1), SIZE(cwork, 2))
3415 : CALL dgemm('N', 'N', SIZE(kroti, 1), SIZE(cwork, 2), SIZE(kroti, 2), &
3416 : 1.0_dp, kroti, SIZE(kroti, 1), cwork, SIZE(cwork, 1), &
3417 : 0.0_dp, twork, SIZE(twork, 1))
3418 : CALL dgemm('N', 'T', SIZE(twork, 1), SIZE(krotj, 1), SIZE(twork, 2), &
3419 : 1.0_dp, twork, SIZE(twork, 1), krotj, SIZE(krotj, 1), &
3420 : 1.0_dp, scblock, SIZE(scblock, 1))
3421 : END IF
3422 : END IF
3423 : END DO
3424 123802 : CALL dbcsr_iterator_stop(iter)
3425 : !$OMP END PARALLEL
3426 123802 : IF (numnodes /= 1 .AND. (perm .OR. has_phase)) THEN
3427 116958 : CALL dbcsr_distribute(rpmat)
3428 116958 : IF (.NOT. real_only) CALL dbcsr_distribute(cpmat)
3429 : END IF
3430 :
3431 123802 : CALL timestop(handle)
3432 :
3433 324592 : END SUBROUTINE symtrans_phase
3434 :
3435 : ! **************************************************************************************************
3436 : !> \brief Symmetrization of density matrix - transform to new k-point
3437 : !> \param smat density matrix at new kpoint
3438 : !> \param pmat reference density matrix
3439 : !> \param kmat Kind type rotation matrix
3440 : !> \param rot Rotation matrix
3441 : !> \param f0 Permutation of atoms under transformation
3442 : !> \param atype Atom to kind pointer
3443 : !> \param symmetric Symmetric matrix
3444 : !> \param antisymmetric Anti-Symmetric matrix
3445 : ! **************************************************************************************************
3446 0 : SUBROUTINE symtrans(smat, pmat, kmat, rot, f0, atype, symmetric, antisymmetric)
3447 : TYPE(dbcsr_type), POINTER :: smat, pmat
3448 : TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kmat
3449 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: rot
3450 : INTEGER, DIMENSION(:), INTENT(IN) :: f0, atype
3451 : LOGICAL, INTENT(IN), OPTIONAL :: symmetric, antisymmetric
3452 :
3453 : CHARACTER(LEN=*), PARAMETER :: routineN = 'symtrans'
3454 :
3455 : INTEGER :: handle, iatom, icol, ikind, ip, irow, &
3456 : jcol, jkind, jp, jrow, natom, numnodes
3457 : LOGICAL :: asym, dorot, found, perm, sym, trans
3458 : REAL(KIND=dp) :: dr, fsign
3459 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
3460 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: kroti, krotj, pblock, sblock
3461 : TYPE(dbcsr_distribution_type) :: dist
3462 : TYPE(dbcsr_iterator_type) :: iter
3463 :
3464 0 : CALL timeset(routineN, handle)
3465 :
3466 : ! check symmetry options
3467 0 : sym = .FALSE.
3468 0 : IF (PRESENT(symmetric)) sym = symmetric
3469 0 : asym = .FALSE.
3470 0 : IF (PRESENT(antisymmetric)) asym = antisymmetric
3471 :
3472 0 : CPASSERT(.NOT. (sym .AND. asym))
3473 0 : CPASSERT((sym .OR. asym))
3474 :
3475 : ! do we have permutation of atoms
3476 0 : natom = SIZE(f0)
3477 0 : perm = .FALSE.
3478 0 : DO iatom = 1, natom
3479 0 : IF (f0(iatom) == iatom) CYCLE
3480 : perm = .TRUE.
3481 0 : EXIT
3482 : END DO
3483 :
3484 : ! do we have a real rotation
3485 0 : dorot = .FALSE.
3486 0 : IF (ABS(SUM(ABS(rot)) - 3.0_dp) > 1.e-12_dp) dorot = .TRUE.
3487 0 : dr = ABS(rot(1, 1) - 1.0_dp) + ABS(rot(2, 2) - 1.0_dp) + ABS(rot(3, 3) - 1.0_dp)
3488 0 : IF (ABS(dr) > 1.e-12_dp) dorot = .TRUE.
3489 :
3490 0 : fsign = 1.0_dp
3491 0 : IF (asym) fsign = -1.0_dp
3492 :
3493 0 : IF (dorot .OR. perm) THEN
3494 : CALL cp_abort(__LOCATION__, "k-points need FULL_GRID currently. "// &
3495 0 : "Reduced grids not yet working correctly")
3496 0 : CALL dbcsr_set(smat, 0.0_dp)
3497 0 : IF (perm) THEN
3498 0 : CALL dbcsr_get_info(pmat, distribution=dist)
3499 0 : CALL dbcsr_distribution_get(dist, numnodes=numnodes)
3500 0 : IF (numnodes == 1) THEN
3501 : ! the matrices are local to this process
3502 0 : CALL dbcsr_iterator_start(iter, pmat)
3503 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3504 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, pblock)
3505 0 : ip = f0(irow)
3506 0 : jp = f0(icol)
3507 0 : IF (ip <= jp) THEN
3508 0 : jrow = ip
3509 0 : jcol = jp
3510 0 : trans = .FALSE.
3511 : ELSE
3512 0 : jrow = jp
3513 0 : jcol = ip
3514 0 : trans = .TRUE.
3515 : END IF
3516 0 : CALL dbcsr_get_block_p(matrix=smat, row=jrow, col=jcol, BLOCK=sblock, found=found)
3517 0 : CPASSERT(found)
3518 0 : ikind = atype(irow)
3519 0 : jkind = atype(icol)
3520 0 : kroti => kmat(ikind)%rmat
3521 0 : krotj => kmat(jkind)%rmat
3522 : ! rotation
3523 0 : IF (trans) THEN
3524 0 : CALL ensure_work_matrix(work, SIZE(krotj, 2), SIZE(pblock, 1))
3525 : CALL dgemm('T', 'T', SIZE(krotj, 2), SIZE(pblock, 1), SIZE(krotj, 1), &
3526 : 1.0_dp, krotj, SIZE(krotj, 1), pblock, SIZE(pblock, 1), &
3527 0 : 0.0_dp, work, SIZE(work, 1))
3528 : CALL dgemm('N', 'N', SIZE(work, 1), SIZE(kroti, 2), SIZE(work, 2), &
3529 : fsign, work, SIZE(work, 1), kroti, SIZE(kroti, 1), &
3530 0 : 0.0_dp, sblock, SIZE(sblock, 1))
3531 : ELSE
3532 0 : CALL ensure_work_matrix(work, SIZE(kroti, 2), SIZE(pblock, 2))
3533 : CALL dgemm('T', 'N', SIZE(kroti, 2), SIZE(pblock, 2), SIZE(kroti, 1), &
3534 : 1.0_dp, kroti, SIZE(kroti, 1), pblock, SIZE(pblock, 1), &
3535 0 : 0.0_dp, work, SIZE(work, 1))
3536 : CALL dgemm('N', 'N', SIZE(work, 1), SIZE(krotj, 2), SIZE(work, 2), &
3537 : fsign, work, SIZE(work, 1), krotj, SIZE(krotj, 1), &
3538 0 : 0.0_dp, sblock, SIZE(sblock, 1))
3539 : END IF
3540 : END DO
3541 0 : CALL dbcsr_iterator_stop(iter)
3542 : !
3543 : ELSE
3544 : ! distributed matrices, most general code needed
3545 : CALL cp_abort(__LOCATION__, "k-points need FULL_GRID currently. "// &
3546 0 : "Reduced grids not yet working correctly")
3547 : END IF
3548 : ELSE
3549 : ! no atom permutations, this is always local
3550 0 : CALL dbcsr_copy(smat, pmat)
3551 0 : CALL dbcsr_iterator_start(iter, smat)
3552 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3553 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
3554 0 : ip = f0(irow)
3555 0 : jp = f0(icol)
3556 0 : IF (ip <= jp) THEN
3557 : jrow = ip
3558 : jcol = jp
3559 0 : trans = .FALSE.
3560 : ELSE
3561 : jrow = jp
3562 : jcol = ip
3563 0 : trans = .TRUE.
3564 : END IF
3565 0 : ikind = atype(irow)
3566 0 : jkind = atype(icol)
3567 0 : kroti => kmat(ikind)%rmat
3568 0 : krotj => kmat(jkind)%rmat
3569 : ! rotation
3570 0 : IF (trans) THEN
3571 0 : CALL ensure_work_matrix(work, SIZE(krotj, 2), SIZE(sblock, 1))
3572 : CALL dgemm('T', 'T', SIZE(krotj, 2), SIZE(sblock, 1), SIZE(krotj, 1), &
3573 : 1.0_dp, krotj, SIZE(krotj, 1), sblock, SIZE(sblock, 1), &
3574 0 : 0.0_dp, work, SIZE(work, 1))
3575 : CALL dgemm('N', 'N', SIZE(work, 1), SIZE(kroti, 2), SIZE(work, 2), &
3576 : fsign, work, SIZE(work, 1), kroti, SIZE(kroti, 1), &
3577 0 : 0.0_dp, sblock, SIZE(sblock, 1))
3578 : ELSE
3579 0 : CALL ensure_work_matrix(work, SIZE(kroti, 2), SIZE(sblock, 2))
3580 : CALL dgemm('T', 'N', SIZE(kroti, 2), SIZE(sblock, 2), SIZE(kroti, 1), &
3581 : 1.0_dp, kroti, SIZE(kroti, 1), sblock, SIZE(sblock, 1), &
3582 0 : 0.0_dp, work, SIZE(work, 1))
3583 : CALL dgemm('N', 'N', SIZE(work, 1), SIZE(krotj, 2), SIZE(work, 2), &
3584 : fsign, work, SIZE(work, 1), krotj, SIZE(krotj, 1), &
3585 0 : 0.0_dp, sblock, SIZE(sblock, 1))
3586 : END IF
3587 : END DO
3588 0 : CALL dbcsr_iterator_stop(iter)
3589 : !
3590 : END IF
3591 : ELSE
3592 : ! this is the identity operation, just copy the matrix
3593 0 : CALL dbcsr_copy(smat, pmat)
3594 : END IF
3595 :
3596 0 : CALL timestop(handle)
3597 :
3598 0 : END SUBROUTINE symtrans
3599 :
3600 : ! **************************************************************************************************
3601 : !> \brief ...
3602 : !> \param mat ...
3603 : ! **************************************************************************************************
3604 0 : SUBROUTINE matprint(mat)
3605 : TYPE(dbcsr_type), POINTER :: mat
3606 :
3607 : INTEGER :: i, icol, iounit, irow
3608 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: mblock
3609 : TYPE(dbcsr_iterator_type) :: iter
3610 :
3611 0 : iounit = cp_logger_get_default_io_unit()
3612 0 : CALL dbcsr_iterator_start(iter, mat)
3613 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3614 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, mblock)
3615 : !
3616 0 : IF (iounit > 0) THEN
3617 0 : WRITE (iounit, '(A,2I4)') 'BLOCK ', irow, icol
3618 0 : DO i = 1, SIZE(mblock, 1)
3619 0 : WRITE (iounit, '(8F12.6)') mblock(i, :)
3620 : END DO
3621 : END IF
3622 : !
3623 : END DO
3624 0 : CALL dbcsr_iterator_stop(iter)
3625 :
3626 0 : END SUBROUTINE matprint
3627 : ! **************************************************************************************************
3628 :
3629 0 : END MODULE kpoint_methods
|
