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