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 Interface to Wannier90 code
10 : !> \par History
11 : !> 06.2016 created [JGH]
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE qs_wannier90
15 : USE atomic_kind_types, ONLY: get_atomic_kind
16 : USE cell_types, ONLY: cell_type,&
17 : get_cell
18 : USE cp_blacs_env, ONLY: cp_blacs_env_type
19 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm
20 : USE cp_cfm_types, ONLY: cp_cfm_create,&
21 : cp_cfm_get_submatrix,&
22 : cp_cfm_release,&
23 : cp_cfm_to_fm,&
24 : cp_cfm_type,&
25 : cp_fm_to_cfm
26 : USE cp_control_types, ONLY: dft_control_type
27 : USE cp_dbcsr_api, ONLY: &
28 : dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_p_type, dbcsr_set, &
29 : dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
30 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
31 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
32 : dbcsr_deallocate_matrix_set
33 : USE cp_files, ONLY: close_file,&
34 : open_file
35 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
36 : cp_fm_struct_release,&
37 : cp_fm_struct_type
38 : USE cp_fm_types, ONLY: cp_fm_copy_general,&
39 : cp_fm_create,&
40 : cp_fm_get_element,&
41 : cp_fm_get_info,&
42 : cp_fm_get_submatrix,&
43 : cp_fm_release,&
44 : cp_fm_set_submatrix,&
45 : cp_fm_type
46 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit,&
47 : cp_logger_type
48 : USE input_section_types, ONLY: section_vals_get,&
49 : section_vals_get_subs_vals,&
50 : section_vals_type,&
51 : section_vals_val_get
52 : USE kinds, ONLY: default_string_length,&
53 : dp
54 : USE kpoint_methods, ONLY: kpoint_env_initialize,&
55 : kpoint_init_cell_index,&
56 : kpoint_initialize,&
57 : kpoint_initialize_mo_set,&
58 : kpoint_initialize_mos,&
59 : rskp_transform
60 : USE kpoint_types, ONLY: get_kpoint_info,&
61 : kind_rotmat_type,&
62 : kpoint_create,&
63 : kpoint_env_type,&
64 : kpoint_release,&
65 : kpoint_sym_type,&
66 : kpoint_type
67 : USE machine, ONLY: m_timestamp,&
68 : timestamp_length
69 : USE mathconstants, ONLY: twopi
70 : USE mathlib, ONLY: diag_complex
71 : USE message_passing, ONLY: mp_para_env_type
72 : USE particle_types, ONLY: particle_type
73 : USE physcon, ONLY: angstrom,&
74 : evolt
75 : USE qs_environment_types, ONLY: get_qs_env,&
76 : qs_env_release,&
77 : qs_environment_type
78 : USE qs_gamma2kp, ONLY: create_kp_from_gamma
79 : USE qs_mo_types, ONLY: get_mo_set,&
80 : mo_set_type
81 : USE qs_moments, ONLY: build_berry_kpoint_matrix
82 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
83 : USE qs_scf_diagonalization, ONLY: do_general_diag_kp
84 : USE qs_scf_types, ONLY: qs_scf_env_type
85 : USE scf_control_types, ONLY: scf_control_type
86 : USE wannier90, ONLY: wannier_setup
87 : #include "./base/base_uses.f90"
88 :
89 : IMPLICIT NONE
90 : PRIVATE
91 :
92 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_wannier90'
93 : INTEGER, PARAMETER, PRIVATE :: w90_kpoints_mp_grid = 0, &
94 : w90_kpoints_scf = 1
95 :
96 : TYPE berry_matrix_type
97 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: sinmat => NULL(), cosmat => NULL()
98 : END TYPE berry_matrix_type
99 :
100 : PUBLIC :: wannier90_interface, prepare_wannier90_scf_mos
101 :
102 : ! **************************************************************************************************
103 :
104 : CONTAINS
105 :
106 : ! **************************************************************************************************
107 : !> \brief ...
108 : !> \param input ...
109 : !> \param logger ...
110 : !> \param qs_env ...
111 : ! **************************************************************************************************
112 22654 : SUBROUTINE wannier90_interface(input, logger, qs_env)
113 : TYPE(section_vals_type), POINTER :: input
114 : TYPE(cp_logger_type), POINTER :: logger
115 : TYPE(qs_environment_type), POINTER :: qs_env
116 :
117 : CHARACTER(len=*), PARAMETER :: routineN = 'wannier90_interface'
118 :
119 : INTEGER :: handle, iw
120 : LOGICAL :: explicit
121 : TYPE(section_vals_type), POINTER :: w_input
122 :
123 : !--------------------------------------------------------------------------------------------!
124 :
125 11327 : CALL timeset(routineN, handle)
126 : w_input => section_vals_get_subs_vals(section_vals=input, &
127 11327 : subsection_name="DFT%PRINT%WANNIER90")
128 11327 : CALL section_vals_get(w_input, explicit=explicit)
129 11327 : IF (explicit) THEN
130 :
131 32 : iw = cp_logger_get_default_io_unit(logger)
132 :
133 32 : IF (iw > 0) THEN
134 : WRITE (iw, '(/,T2,A)') &
135 16 : '!-----------------------------------------------------------------------------!'
136 16 : WRITE (iw, '(T32,A)') "Interface to Wannier90"
137 : WRITE (iw, '(T2,A)') &
138 16 : '!-----------------------------------------------------------------------------!'
139 : END IF
140 :
141 32 : CALL wannier90_files(qs_env, w_input, iw)
142 :
143 32 : IF (iw > 0) THEN
144 : WRITE (iw, '(/,T2,A)') &
145 16 : '!--------------------------------End of Wannier90-----------------------------!'
146 : END IF
147 : END IF
148 11327 : CALL timestop(handle)
149 :
150 11327 : END SUBROUTINE wannier90_interface
151 :
152 : ! **************************************************************************************************
153 : !> \brief ...
154 : !> \param qs_env ...
155 : !> \param input ...
156 : !> \param iw ...
157 : ! **************************************************************************************************
158 32 : SUBROUTINE wannier90_files(qs_env, input, iw)
159 : TYPE(qs_environment_type), POINTER :: qs_env
160 : TYPE(section_vals_type), POINTER :: input
161 : INTEGER, INTENT(IN) :: iw
162 :
163 : INTEGER, PARAMETER :: num_nnmax = 12
164 :
165 : CHARACTER(len=2) :: asym
166 32 : CHARACTER(len=20), ALLOCATABLE, DIMENSION(:) :: atom_symbols
167 : CHARACTER(len=default_string_length) :: filename, input_kp_scheme, reuse_reason, &
168 : seed_name
169 : CHARACTER(LEN=timestamp_length) :: timestamp
170 : INTEGER :: aligned_degenerate_blocks, aligned_degenerate_max_size, i, i_rep, ib, ib1, ib2, &
171 : ibs, ik, ik2, ikk, ikpgr, ispin, iunit, ix, iy, iz, k, kpoints_source, n_rep, nadd, nao, &
172 : nbs, nexcl, nkp, nmo, nntot, nspins, num_atoms, num_bands, num_bands_tot, num_kpts, &
173 : num_wann
174 32 : INTEGER, ALLOCATABLE, DIMENSION(:) :: exclude_bands
175 32 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nblist, nnlist
176 32 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: nncell
177 : INTEGER, DIMENSION(2) :: kp_range
178 : INTEGER, DIMENSION(3) :: input_nkp_grid, mp_grid
179 32 : INTEGER, DIMENSION(:), POINTER :: invals
180 32 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
181 : LOGICAL :: diis_step, do_kpoints, full_mesh_diagonalized, gamma_only, input_full_grid, &
182 : input_gamma_centered, input_kpoint_symmetry, mp_grid_explicit, mp_grid_valid, my_kpgrp, &
183 : mygrp, reuse_scf_mos, reused_scf_mos, spinors, use_bloch_phases, validate_reuse_ok, &
184 : validate_reuse_scf_mos
185 : REAL(KIND=dp) :: aligned_degenerate_min_svalue, cmmn, gauge_arg, gauge_imag, gauge_real, &
186 : gauge_tmp, ksign, reuse_candidate_deviation, reuse_candidate_metric_deviation, &
187 : reuse_candidate_min_svalue, reuse_candidate_residual, rmmn, &
188 : validation_eigenvalue_deviation, validation_min_svalue, validation_subspace_deviation, &
189 : wkp_ref
190 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigval
191 64 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: atoms_cart, b_latt, kpt_latt
192 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: reference_eigenvalues
193 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: reference_mo_imag, reference_mo_real
194 : REAL(KIND=dp), DIMENSION(3) :: bvec, input_kp_shift, phase_center
195 : REAL(KIND=dp), DIMENSION(3, 3) :: h_inv, real_lattice, recip_lattice
196 64 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, wkp, wkp_source
197 32 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp, xkp_source
198 32 : TYPE(berry_matrix_type), DIMENSION(:), POINTER :: berry_matrix
199 : TYPE(cell_type), POINTER :: cell
200 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
201 : TYPE(cp_cfm_type) :: fmk1_cfm, fmk2_cfm, mmn_cfm, omat_cfm, &
202 : tmp_cfm
203 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_ao, matrix_struct_mmn, &
204 : matrix_struct_work
205 : TYPE(cp_fm_type) :: mat_imag, mat_real, mmn_imag, mmn_real
206 192 : TYPE(cp_fm_type), DIMENSION(2) :: fmk1, fmk2
207 : TYPE(cp_fm_type), POINTER :: fmdummy, fmi, fmr
208 32 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
209 : TYPE(dbcsr_type), POINTER :: cmatrix, cmatrix_full, rmatrix, &
210 : rmatrix_full
211 : TYPE(dft_control_type), POINTER :: dft_control
212 : TYPE(kpoint_env_type), POINTER :: kp
213 : TYPE(kpoint_type), POINTER :: kpoint, qs_kpoint
214 32 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
215 : TYPE(mp_para_env_type), POINTER :: para_env
216 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
217 32 : POINTER :: sab_nl
218 32 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
219 : TYPE(qs_environment_type), POINTER :: qs_env_kp
220 : TYPE(qs_scf_env_type), POINTER :: scf_env
221 : TYPE(scf_control_type), POINTER :: scf_control
222 :
223 : !--------------------------------------------------------------------------------------------!
224 :
225 : ! add code for exclude_bands and projectors
226 :
227 : ! generate all arrays needed for the setup call
228 32 : CALL section_vals_val_get(input, "SEED_NAME", c_val=seed_name)
229 32 : CALL section_vals_val_get(input, "MP_GRID", i_vals=invals, explicit=mp_grid_explicit)
230 32 : CALL section_vals_val_get(input, "KPOINTS_SOURCE", i_val=kpoints_source)
231 32 : CALL section_vals_val_get(input, "WANNIER_FUNCTIONS", i_val=num_wann)
232 32 : CALL section_vals_val_get(input, "ADDED_MOS", i_val=nadd)
233 32 : CALL section_vals_val_get(input, "REUSE_SCF_MOS", l_val=reuse_scf_mos)
234 32 : CALL section_vals_val_get(input, "VALIDATE_REUSE_SCF_MOS", l_val=validate_reuse_scf_mos)
235 32 : CALL section_vals_val_get(input, "USE_BLOCH_PHASES", l_val=use_bloch_phases)
236 32 : reuse_scf_mos = reuse_scf_mos .AND. kpoints_source == w90_kpoints_scf
237 32 : validate_reuse_scf_mos = validate_reuse_scf_mos .AND. reuse_scf_mos
238 128 : mp_grid(1:3) = invals(1:3)
239 : ! excluded bands
240 32 : CALL section_vals_val_get(input, "EXCLUDE_BANDS", n_rep_val=n_rep)
241 32 : nexcl = 0
242 32 : DO i_rep = 1, n_rep
243 0 : CALL section_vals_val_get(input, "EXCLUDE_BANDS", i_rep_val=i_rep, i_vals=invals)
244 32 : nexcl = nexcl + SIZE(invals)
245 : END DO
246 32 : IF (nexcl > 0) THEN
247 0 : ALLOCATE (exclude_bands(nexcl))
248 0 : nexcl = 0
249 0 : DO i_rep = 1, n_rep
250 0 : CALL section_vals_val_get(input, "EXCLUDE_BANDS", i_rep_val=i_rep, i_vals=invals)
251 0 : exclude_bands(nexcl + 1:nexcl + SIZE(invals)) = invals(:)
252 0 : nexcl = nexcl + SIZE(invals)
253 : END DO
254 : END IF
255 : !
256 : ! lattice -> Angstrom
257 32 : CALL get_qs_env(qs_env, cell=cell)
258 32 : CALL get_cell(cell, h=real_lattice, h_inv=h_inv)
259 : ! k-points
260 32 : CALL get_qs_env(qs_env, particle_set=particle_set)
261 32 : CALL get_qs_env(qs_env, para_env=para_env)
262 32 : phase_center = 0.0_dp
263 222 : DO i = 1, SIZE(particle_set)
264 3072 : phase_center(1:3) = phase_center(1:3) + MATMUL(h_inv, particle_set(i)%r)
265 : END DO
266 128 : phase_center(1:3) = phase_center(1:3)/REAL(SIZE(particle_set), KIND=dp)
267 128 : phase_center(1:3) = phase_center(1:3) - FLOOR(phase_center(1:3))
268 32 : recip_lattice(1:3, 1:3) = h_inv(1:3, 1:3)
269 416 : real_lattice(1:3, 1:3) = angstrom*real_lattice(1:3, 1:3)
270 800 : recip_lattice(1:3, 1:3) = (twopi/angstrom)*TRANSPOSE(recip_lattice(1:3, 1:3))
271 32 : NULLIFY (kpoint, qs_kpoint, xkp, wkp, xkp_source, wkp_source)
272 32 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints, kpoints=qs_kpoint)
273 32 : input_kpoint_symmetry = .FALSE.
274 32 : input_full_grid = .FALSE.
275 32 : input_kp_scheme = ""
276 32 : IF (do_kpoints .AND. ASSOCIATED(qs_kpoint)) THEN
277 : CALL get_kpoint_info(qs_kpoint, kp_scheme=input_kp_scheme, nkp_grid=input_nkp_grid, &
278 : kp_shift=input_kp_shift, symmetry=input_kpoint_symmetry, &
279 : full_grid=input_full_grid, gamma_centered=input_gamma_centered, &
280 32 : nkp=nkp, xkp=xkp, wkp=wkp)
281 : END IF
282 32 : CALL kpoint_create(kpoint)
283 :
284 0 : SELECT CASE (kpoints_source)
285 : CASE (w90_kpoints_mp_grid)
286 0 : num_kpts = mp_grid(1)*mp_grid(2)*mp_grid(3)
287 0 : ALLOCATE (kpt_latt(3, num_kpts))
288 0 : kpoint%kp_scheme = "MONKHORST-PACK"
289 0 : kpoint%symmetry = .FALSE.
290 0 : kpoint%nkp_grid(1:3) = mp_grid(1:3)
291 0 : kpoint%verbose = .FALSE.
292 0 : kpoint%full_grid = .TRUE.
293 0 : kpoint%eps_geo = 1.0e-6_dp
294 0 : kpoint%use_real_wfn = .FALSE.
295 0 : kpoint%parallel_group_size = para_env%num_pe
296 0 : i = 0
297 0 : DO ix = 0, mp_grid(1) - 1
298 0 : DO iy = 0, mp_grid(2) - 1
299 0 : DO iz = 0, mp_grid(3) - 1
300 0 : i = i + 1
301 0 : kpt_latt(1, i) = REAL(ix, KIND=dp)/REAL(mp_grid(1), KIND=dp)
302 0 : kpt_latt(2, i) = REAL(iy, KIND=dp)/REAL(mp_grid(2), KIND=dp)
303 0 : kpt_latt(3, i) = REAL(iz, KIND=dp)/REAL(mp_grid(3), KIND=dp)
304 : END DO
305 : END DO
306 : END DO
307 0 : kpoint%nkp = num_kpts
308 0 : ALLOCATE (kpoint%xkp(3, num_kpts), kpoint%wkp(num_kpts))
309 0 : kpoint%wkp(:) = 1._dp/REAL(num_kpts, KIND=dp)
310 0 : DO i = 1, num_kpts
311 0 : kpoint%xkp(1:3, i) = (angstrom/twopi)*MATMUL(recip_lattice, kpt_latt(:, i))
312 : END DO
313 :
314 : CASE (w90_kpoints_scf)
315 32 : IF (.NOT. do_kpoints .OR. .NOT. ASSOCIATED(qs_kpoint)) THEN
316 0 : CPABORT("WANNIER90%KPOINTS_SOURCE SCF requires an active DFT%KPOINTS section.")
317 : END IF
318 32 : SELECT CASE (TRIM(input_kp_scheme))
319 : CASE ("GAMMA")
320 0 : mp_grid(:) = 1
321 0 : num_kpts = 1
322 0 : ALLOCATE (kpt_latt(3, num_kpts))
323 0 : kpt_latt(1:3, 1) = 0.0_dp
324 0 : kpoint%kp_scheme = "GAMMA"
325 0 : kpoint%symmetry = .FALSE.
326 0 : kpoint%verbose = .FALSE.
327 0 : kpoint%full_grid = .TRUE.
328 0 : kpoint%eps_geo = 1.0e-6_dp
329 0 : kpoint%use_real_wfn = .FALSE.
330 0 : kpoint%parallel_group_size = para_env%num_pe
331 0 : kpoint%nkp = num_kpts
332 0 : ALLOCATE (kpoint%xkp(3, num_kpts), kpoint%wkp(num_kpts))
333 0 : kpoint%xkp(1:3, 1) = 0.0_dp
334 0 : kpoint%wkp(1) = 1.0_dp
335 :
336 : CASE ("MONKHORST-PACK", "MACDONALD")
337 26 : mp_grid(1:3) = input_nkp_grid(1:3)
338 26 : kpoint%kp_scheme = input_kp_scheme
339 26 : kpoint%symmetry = .FALSE.
340 104 : kpoint%nkp_grid(1:3) = input_nkp_grid(1:3)
341 104 : kpoint%kp_shift(1:3) = input_kp_shift(1:3)
342 26 : kpoint%gamma_centered = input_gamma_centered
343 26 : kpoint%verbose = .FALSE.
344 26 : kpoint%full_grid = .TRUE.
345 26 : kpoint%eps_geo = 1.0e-6_dp
346 26 : kpoint%use_real_wfn = .FALSE.
347 26 : kpoint%parallel_group_size = para_env%num_pe
348 26 : CALL kpoint_initialize(kpoint, particle_set, cell)
349 26 : num_kpts = kpoint%nkp
350 78 : ALLOCATE (kpt_latt(3, num_kpts))
351 1754 : kpt_latt(1:3, 1:num_kpts) = kpoint%xkp(1:3, 1:num_kpts)
352 26 : IF (input_kpoint_symmetry .AND. .NOT. input_full_grid .AND. iw > 0) THEN
353 : WRITE (iw, '(T2,A)') &
354 12 : "WANNIER90| SCF k-points are symmetry-reduced; regenerating the full SCF mesh."
355 12 : IF (reuse_scf_mos) THEN
356 : WRITE (iw, '(T2,A)') &
357 12 : "WANNIER90| CP2K will try to reconstruct the full-mesh MOs from the SCF orbitals."
358 : ELSE
359 : WRITE (iw, '(T2,A)') &
360 0 : "WANNIER90| The full exported mesh is diagonalized for the Wannier90 files."
361 : END IF
362 : END IF
363 :
364 : CASE ("GENERAL")
365 6 : IF (ASSOCIATED(qs_kpoint%xkp_input)) THEN
366 6 : xkp_source => qs_kpoint%xkp_input
367 6 : wkp_source => qs_kpoint%wkp_input
368 : ELSE
369 0 : xkp_source => xkp
370 0 : wkp_source => wkp
371 : END IF
372 6 : IF (.NOT. ASSOCIATED(xkp_source) .OR. .NOT. ASSOCIATED(wkp_source)) THEN
373 0 : CPABORT("Could not access the SCF GENERAL k-point set for the Wannier90 export.")
374 : END IF
375 6 : num_kpts = SIZE(wkp_source)
376 18 : ALLOCATE (kpt_latt(3, num_kpts))
377 198 : kpt_latt(1:3, 1:num_kpts) = xkp_source(1:3, 1:num_kpts)
378 6 : IF (mp_grid_explicit) THEN
379 0 : IF (mp_grid(1)*mp_grid(2)*mp_grid(3) /= num_kpts) THEN
380 0 : CPABORT("WANNIER90%MP_GRID must contain exactly as many points as the SCF GENERAL mesh.")
381 : END IF
382 : ELSE
383 6 : CALL infer_wannier_mp_grid(kpt_latt, mp_grid, mp_grid_valid)
384 6 : IF (.NOT. mp_grid_valid) THEN
385 0 : CPABORT("Could not infer WANNIER90%MP_GRID from the SCF GENERAL mesh.")
386 : END IF
387 : END IF
388 6 : wkp_ref = 1.0_dp/REAL(num_kpts, KIND=dp)
389 54 : DO i = 1, num_kpts
390 54 : IF (ABS(wkp_source(i) - wkp_ref) > 1.0e-10_dp) THEN
391 0 : CPABORT("WANNIER90%KPOINTS_SOURCE SCF requires equally weighted GENERAL k-points.")
392 : END IF
393 : END DO
394 6 : kpoint%kp_scheme = "GENERAL"
395 6 : kpoint%symmetry = .FALSE.
396 24 : kpoint%nkp_grid(1:3) = mp_grid(1:3)
397 6 : kpoint%verbose = .FALSE.
398 6 : kpoint%full_grid = .TRUE.
399 6 : kpoint%eps_geo = 1.0e-6_dp
400 6 : kpoint%use_real_wfn = .FALSE.
401 6 : kpoint%parallel_group_size = para_env%num_pe
402 6 : kpoint%nkp = num_kpts
403 30 : ALLOCATE (kpoint%xkp(3, num_kpts), kpoint%wkp(num_kpts))
404 390 : kpoint%xkp(1:3, 1:num_kpts) = xkp_source(1:3, 1:num_kpts)
405 54 : kpoint%wkp(1:num_kpts) = wkp_ref
406 6 : IF (input_kpoint_symmetry .AND. .NOT. input_full_grid .AND. iw > 0) THEN
407 : WRITE (iw, '(T2,A)') &
408 2 : "WANNIER90| SCF k-points are symmetry-reduced; using the full input GENERAL mesh."
409 2 : IF (reuse_scf_mos) THEN
410 : WRITE (iw, '(T2,A)') &
411 2 : "WANNIER90| CP2K will try to reconstruct the full-mesh MOs from the SCF orbitals."
412 : ELSE
413 : WRITE (iw, '(T2,A)') &
414 0 : "WANNIER90| The full exported mesh is diagonalized for the Wannier90 files."
415 : END IF
416 : END IF
417 :
418 : CASE DEFAULT
419 32 : CPABORT("WANNIER90%KPOINTS_SOURCE SCF does not support this DFT%KPOINTS scheme.")
420 : END SELECT
421 : CASE DEFAULT
422 32 : CPABORT("Unknown WANNIER90%KPOINTS_SOURCE setting.")
423 : END SELECT
424 : ! number of bands in calculation
425 32 : CALL get_qs_env(qs_env, mos=mos)
426 32 : CALL get_mo_set(mo_set=mos(1), nao=nao, nmo=num_bands_tot)
427 32 : num_bands_tot = MIN(nao, num_bands_tot + nadd)
428 32 : num_bands = num_bands_tot
429 32 : IF (use_bloch_phases .AND. num_wann /= num_bands) THEN
430 0 : CPABORT("WANNIER90%USE_BLOCH_PHASES requires WANNIER_FUNCTIONS to match the number of bands.")
431 : END IF
432 32 : num_atoms = SIZE(particle_set)
433 96 : ALLOCATE (atoms_cart(3, num_atoms))
434 96 : ALLOCATE (atom_symbols(num_atoms))
435 222 : DO i = 1, num_atoms
436 760 : atoms_cart(1:3, i) = particle_set(i)%r(1:3)
437 190 : CALL get_atomic_kind(particle_set(i)%atomic_kind, element_symbol=asym)
438 222 : atom_symbols(i) = asym
439 : END DO
440 32 : gamma_only = .FALSE.
441 32 : spinors = .FALSE.
442 : ! output
443 96 : ALLOCATE (nnlist(num_kpts, num_nnmax))
444 128 : ALLOCATE (nncell(3, num_kpts, num_nnmax))
445 32 : nnlist(:, :) = 0
446 32 : nncell(:, :, :) = 0
447 32 : nntot = 0
448 :
449 32 : IF (iw > 0) THEN
450 : ! setup
451 : CALL wannier_setup(mp_grid, num_kpts, real_lattice, recip_lattice, &
452 16 : kpt_latt, nntot, nnlist, nncell, iw)
453 : END IF
454 :
455 32 : CALL get_qs_env(qs_env, para_env=para_env)
456 32 : CALL para_env%sum(nntot)
457 32 : CALL para_env%sum(nnlist)
458 32 : CALL para_env%sum(nncell)
459 :
460 32 : IF (para_env%is_source()) THEN
461 : ! Write the Wannier90 input file "seed_name.win"
462 16 : WRITE (filename, '(A,A)') TRIM(seed_name), ".win"
463 16 : CALL open_file(filename, unit_number=iunit, file_status="UNKNOWN", file_action="WRITE")
464 : !
465 16 : CALL m_timestamp(timestamp)
466 16 : WRITE (iunit, "(A)") "! Wannier90 input file generated by CP2K "
467 16 : WRITE (iunit, "(A,/)") "! Creation date "//timestamp
468 : !
469 16 : WRITE (iunit, "(A,I5)") "num_wann = ", num_wann
470 16 : IF (num_bands /= num_wann .OR. use_bloch_phases) THEN
471 14 : WRITE (iunit, "(A,I5)") "num_bands = ", num_bands
472 : END IF
473 16 : IF (use_bloch_phases) THEN
474 : ! Keep the external Wannier90 projection matrix fully defined for
475 : ! complete-band Bloch-phase subspaces by writing explicit identity projections.
476 6 : WRITE (iunit, "(A)") "! CP2K writes identity projections for Bloch-phase complete subspaces."
477 : END IF
478 16 : WRITE (iunit, "(/,A,/)") "length_unit = bohr "
479 16 : WRITE (iunit, "(/,A,/)") "! System"
480 16 : WRITE (iunit, "(/,A)") "begin unit_cell_cart"
481 16 : WRITE (iunit, "(A)") "bohr"
482 64 : DO i = 1, 3
483 208 : WRITE (iunit, "(3F12.6)") cell%hmat(i, 1:3)
484 : END DO
485 16 : WRITE (iunit, "(A,/)") "end unit_cell_cart"
486 16 : WRITE (iunit, "(/,A)") "begin atoms_cart"
487 16 : WRITE (iunit, "(A)") "bohr"
488 111 : DO i = 1, num_atoms
489 111 : WRITE (iunit, "(A,3F15.10)") atom_symbols(i), atoms_cart(1:3, i)
490 : END DO
491 16 : WRITE (iunit, "(A,/)") "end atoms_cart"
492 16 : WRITE (iunit, "(/,A,/)") "! Kpoints"
493 16 : WRITE (iunit, "(/,A,3I6/)") "mp_grid = ", mp_grid(1:3)
494 16 : WRITE (iunit, "(A)") "begin kpoints"
495 256 : DO i = 1, num_kpts
496 256 : WRITE (iunit, "(3F12.6)") kpt_latt(1:3, i)
497 : END DO
498 16 : WRITE (iunit, "(A)") "end kpoints"
499 16 : CALL close_file(iunit)
500 16 : IF (use_bloch_phases) THEN
501 6 : WRITE (filename, '(A,A)') TRIM(seed_name), ".amn"
502 6 : CALL open_file(filename, unit_number=iunit, file_status="UNKNOWN", file_action="WRITE")
503 6 : WRITE (iunit, "(A)") "! Wannier90 identity projections generated by CP2K"
504 6 : WRITE (iunit, "(3I8)") num_bands, num_kpts, num_wann
505 166 : DO ik = 1, num_kpts
506 742 : DO ib2 = 1, num_wann
507 2960 : DO ib1 = 1, num_bands
508 2800 : IF (ib1 == ib2) THEN
509 576 : WRITE (iunit, "(3I8,2E30.14)") ib1, ib2, ik, 1.0_dp, 0.0_dp
510 : ELSE
511 1648 : WRITE (iunit, "(3I8,2E30.14)") ib1, ib2, ik, 0.0_dp, 0.0_dp
512 : END IF
513 : END DO
514 : END DO
515 : END DO
516 6 : CALL close_file(iunit)
517 : END IF
518 : ELSE
519 16 : iunit = -1
520 : END IF
521 :
522 : ! calculate bands
523 32 : NULLIFY (qs_env_kp)
524 32 : IF (kpoints_source == w90_kpoints_mp_grid .AND. input_kpoint_symmetry .AND. iw > 0) THEN
525 : WRITE (iw, '(T2,A)') &
526 0 : "WANNIER90| Atomic k-point symmetry from the SCF calculation is not reused."
527 : WRITE (iw, '(T2,A)') &
528 0 : "WANNIER90| A full Monkhorst-Pack grid is generated for the Wannier90 interface."
529 : END IF
530 32 : IF (do_kpoints) THEN
531 : ! we already do kpoints
532 32 : qs_env_kp => qs_env
533 : ELSE
534 : ! we start from gamma point only
535 0 : ALLOCATE (qs_env_kp)
536 0 : CALL create_kp_from_gamma(qs_env, qs_env_kp)
537 : END IF
538 32 : IF (iw > 0) THEN
539 16 : WRITE (unit=iw, FMT="(/,T2,A)") "Start K-Point Calculation ..."
540 : END IF
541 32 : CALL get_qs_env(qs_env=qs_env_kp, para_env=para_env, blacs_env=blacs_env)
542 32 : CALL kpoint_env_initialize(kpoint, para_env, blacs_env)
543 32 : CALL kpoint_initialize_mos(kpoint, mos, nadd)
544 32 : CALL kpoint_initialize_mo_set(kpoint)
545 : !
546 32 : CALL get_qs_env(qs_env=qs_env_kp, sab_orb=sab_nl, dft_control=dft_control)
547 32 : CALL kpoint_init_cell_index(kpoint, sab_nl, para_env, dft_control%nimages)
548 : !
549 : CALL get_qs_env(qs_env=qs_env_kp, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s, &
550 32 : scf_env=scf_env, scf_control=scf_control)
551 32 : full_mesh_diagonalized = .FALSE.
552 32 : reused_scf_mos = .FALSE.
553 32 : reuse_reason = ""
554 32 : aligned_degenerate_blocks = 0
555 32 : aligned_degenerate_max_size = 0
556 32 : aligned_degenerate_min_svalue = 0.0_dp
557 32 : IF (reuse_scf_mos) THEN
558 32 : CALL get_kpoint_info(kpoint=kpoint, cell_to_index=cell_to_index)
559 : CALL do_general_diag_kp(matrix_ks, matrix_s, qs_kpoint, scf_env, scf_control, .FALSE., &
560 32 : diis_step)
561 32 : IF (validate_reuse_scf_mos) THEN
562 6 : IF (iw > 0) THEN
563 : WRITE (iw, '(T2,A)') &
564 3 : "WANNIER90| Validating SCF MO reuse against a full-mesh diagonalization reference."
565 : END IF
566 6 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoint, scf_env, scf_control, .FALSE., diis_step)
567 6 : full_mesh_diagonalized = .TRUE.
568 6 : nspins = dft_control%nspins
569 : CALL save_wannier90_mo_snapshot(kpoint, nspins, para_env, reference_mo_real, &
570 6 : reference_mo_imag, reference_eigenvalues)
571 : CALL diagnose_wannier90_scf_reuse_candidates(kpoint, qs_kpoint, matrix_s, matrix_ks, &
572 : cell_to_index, sab_nl, para_env, iw, &
573 : reuse_candidate_deviation, &
574 : reuse_candidate_min_svalue, &
575 : reuse_candidate_metric_deviation, &
576 6 : reuse_candidate_residual)
577 6 : IF (iw > 0 .AND. reuse_candidate_deviation < 1.0e100_dp) THEN
578 : WRITE (iw, '(T2,A,ES10.3,A,ES10.3,A,ES10.3,A,ES10.3)') &
579 3 : "WANNIER90| Best atom/AO candidate subspace deviation ", &
580 3 : reuse_candidate_deviation, ", minimum singular value ", &
581 3 : reuse_candidate_min_svalue, ", max metric deviation ", &
582 6 : reuse_candidate_metric_deviation, ", max residual ", reuse_candidate_residual
583 : END IF
584 : END IF
585 : CALL prepare_wannier90_scf_mos(kpoint, qs_kpoint, matrix_s, matrix_ks, cell_to_index, &
586 : sab_nl, para_env, reused_scf_mos, reuse_reason, &
587 : aligned_degenerate_blocks, aligned_degenerate_max_size, &
588 32 : aligned_degenerate_min_svalue)
589 32 : IF (validate_reuse_scf_mos) THEN
590 6 : IF (reused_scf_mos) THEN
591 : CALL validate_wannier90_reused_mos(kpoint, matrix_s, cell_to_index, sab_nl, &
592 : para_env, reference_mo_real, reference_mo_imag, &
593 : reference_eigenvalues, validate_reuse_ok, &
594 : validation_subspace_deviation, validation_min_svalue, &
595 6 : validation_eigenvalue_deviation)
596 6 : IF (iw > 0) THEN
597 : WRITE (iw, '(T2,A,ES10.3,A,ES10.3,A,ES10.3)') &
598 3 : "WANNIER90| Reused MO validation: subspace deviation ", &
599 3 : validation_subspace_deviation, ", minimum singular value ", &
600 3 : validation_min_svalue, ", eigenvalue deviation ", &
601 6 : validation_eigenvalue_deviation
602 : END IF
603 6 : IF (.NOT. validate_reuse_ok) THEN
604 0 : reused_scf_mos = .FALSE.
605 : WRITE (reuse_reason, "(A,ES10.3,A,ES10.3)") &
606 0 : "validation failed: dS=", &
607 0 : validation_subspace_deviation, ", dE=", validation_eigenvalue_deviation
608 : END IF
609 : END IF
610 6 : IF (.NOT. reused_scf_mos) THEN
611 : CALL restore_wannier90_mo_snapshot(kpoint, reference_mo_real, reference_mo_imag, &
612 0 : reference_eigenvalues)
613 : END IF
614 : END IF
615 32 : IF (iw > 0) THEN
616 16 : IF (reused_scf_mos) THEN
617 : WRITE (iw, '(T2,A)') &
618 13 : "WANNIER90| Reused SCF MO coefficients for the Wannier90 full k-point mesh."
619 13 : IF (use_bloch_phases) THEN
620 : WRITE (iw, '(T2,A)') &
621 6 : "WANNIER90| Wrote identity projections for Bloch-phase complete band subspaces."
622 : WRITE (iw, '(T2,A,3F10.6)') &
623 6 : "WANNIER90| Applied Bloch phase gauge to reused overlaps around fractional center", &
624 12 : phase_center(1:3)
625 : END IF
626 13 : IF (aligned_degenerate_blocks > 0) THEN
627 : WRITE (iw, '(T2,A,I0,A,I0,A,ES10.3)') &
628 8 : "WANNIER90| Ritz-stabilized ", aligned_degenerate_blocks, &
629 8 : " degenerate SCF MO subspace(s) with S(k),H(k); largest block has ", &
630 8 : aligned_degenerate_max_size, " band(s), min metric eigenvalue ", &
631 16 : aligned_degenerate_min_svalue
632 : END IF
633 : ELSE
634 : WRITE (iw, '(T2,A,A)') &
635 3 : "WANNIER90| Could not reuse SCF MOs: ", TRIM(reuse_reason)
636 : WRITE (iw, '(T2,A)') &
637 3 : "WANNIER90| Falling back to full-mesh diagonalization for the Wannier90 files."
638 : END IF
639 : END IF
640 : END IF
641 32 : IF (.NOT. reused_scf_mos .AND. .NOT. full_mesh_diagonalized) THEN
642 6 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoint, scf_env, scf_control, .FALSE., diis_step)
643 : END IF
644 32 : IF (ALLOCATED(reference_mo_real)) DEALLOCATE (reference_mo_real)
645 32 : IF (ALLOCATED(reference_mo_imag)) DEALLOCATE (reference_mo_imag)
646 32 : IF (ALLOCATED(reference_eigenvalues)) DEALLOCATE (reference_eigenvalues)
647 : !
648 32 : IF (iw > 0) THEN
649 16 : WRITE (iw, '(T69,A)') "... Finished"
650 : END IF
651 : !
652 : ! Calculate and print Overlaps
653 : !
654 32 : IF (para_env%is_source()) THEN
655 16 : WRITE (filename, '(A,A)') TRIM(seed_name), ".mmn"
656 16 : CALL open_file(filename, unit_number=iunit, file_status="UNKNOWN", file_action="WRITE")
657 16 : CALL m_timestamp(timestamp)
658 16 : WRITE (iunit, "(A)") "! Wannier90 file generated by CP2K "//timestamp
659 16 : WRITE (iunit, "(3I8)") num_bands, num_kpts, nntot
660 : ELSE
661 16 : iunit = -1
662 : END IF
663 : ! create a list of unique b vectors and a table of pointers
664 : ! nblist(ik,i) -> +/- b_latt(1:3,x)
665 128 : ALLOCATE (nblist(num_kpts, nntot))
666 96 : ALLOCATE (b_latt(3, num_kpts*nntot))
667 32 : nblist(:, :) = 0
668 32 : nbs = 0
669 512 : DO ik = 1, num_kpts
670 3392 : DO i = 1, nntot
671 11520 : bvec(1:3) = kpt_latt(1:3, nnlist(ik, i)) - kpt_latt(1:3, ik) + nncell(1:3, ik, i)
672 5760 : ibs = 0
673 5760 : DO k = 1, nbs
674 22656 : IF (SUM(ABS(bvec(1:3) - b_latt(1:3, k))) < 1.e-6_dp) THEN
675 : ibs = k
676 : EXIT
677 : END IF
678 17376 : IF (SUM(ABS(bvec(1:3) + b_latt(1:3, k))) < 1.e-6_dp) THEN
679 1440 : ibs = -k
680 1440 : EXIT
681 : END IF
682 : END DO
683 3360 : IF (ibs /= 0) THEN
684 : ! old lattice vector
685 2784 : nblist(ik, i) = ibs
686 : ELSE
687 : ! new lattice vector
688 96 : nbs = nbs + 1
689 384 : b_latt(1:3, nbs) = bvec(1:3)
690 96 : nblist(ik, i) = nbs
691 : END IF
692 : END DO
693 : END DO
694 : ! calculate all the operator matrices (a|bvec|b)
695 192 : ALLOCATE (berry_matrix(nbs))
696 128 : DO i = 1, nbs
697 96 : NULLIFY (berry_matrix(i)%cosmat)
698 96 : NULLIFY (berry_matrix(i)%sinmat)
699 480 : bvec(1:3) = twopi*MATMUL(TRANSPOSE(cell%h_inv(1:3, 1:3)), b_latt(1:3, i))
700 : CALL build_berry_kpoint_matrix(qs_env_kp, berry_matrix(i)%cosmat, &
701 128 : berry_matrix(i)%sinmat, bvec)
702 : END DO
703 : ! work matrices for MOs (all group)
704 32 : kp => kpoint%kp_env(1)%kpoint_env
705 32 : CALL get_mo_set(kp%mos(1, 1), nmo=nmo)
706 32 : NULLIFY (matrix_struct_ao, matrix_struct_work)
707 : CALL cp_fm_struct_create(matrix_struct_work, nrow_global=nao, &
708 : ncol_global=nmo, &
709 : para_env=para_env, &
710 32 : context=blacs_env)
711 96 : DO i = 1, 2
712 64 : CALL cp_fm_create(fmk1(i), matrix_struct_work)
713 96 : CALL cp_fm_create(fmk2(i), matrix_struct_work)
714 : END DO
715 32 : CALL cp_cfm_create(fmk1_cfm, matrix_struct_work)
716 32 : CALL cp_cfm_create(fmk2_cfm, matrix_struct_work)
717 32 : CALL cp_cfm_create(tmp_cfm, matrix_struct_work)
718 : CALL cp_fm_struct_create(matrix_struct_ao, nrow_global=nao, &
719 : ncol_global=nao, &
720 : para_env=para_env, &
721 32 : context=blacs_env)
722 32 : CALL cp_fm_create(mat_real, matrix_struct_ao)
723 32 : CALL cp_fm_create(mat_imag, matrix_struct_ao)
724 32 : CALL cp_cfm_create(omat_cfm, matrix_struct_ao)
725 : ! work matrices for Mmn(k,b) integrals
726 32 : NULLIFY (matrix_struct_mmn)
727 : CALL cp_fm_struct_create(matrix_struct_mmn, nrow_global=nmo, &
728 : ncol_global=nmo, &
729 : para_env=para_env, &
730 32 : context=blacs_env)
731 32 : CALL cp_fm_create(mmn_real, matrix_struct_mmn)
732 32 : CALL cp_fm_create(mmn_imag, matrix_struct_mmn)
733 32 : CALL cp_cfm_create(mmn_cfm, matrix_struct_mmn)
734 : ! allocate some work matrices
735 32 : ALLOCATE (rmatrix, cmatrix, rmatrix_full, cmatrix_full)
736 : CALL dbcsr_create(rmatrix, template=matrix_s(1, 1)%matrix, &
737 32 : matrix_type=dbcsr_type_symmetric)
738 : CALL dbcsr_create(cmatrix, template=matrix_s(1, 1)%matrix, &
739 32 : matrix_type=dbcsr_type_antisymmetric)
740 : CALL dbcsr_create(rmatrix_full, template=matrix_s(1, 1)%matrix, &
741 32 : matrix_type=dbcsr_type_no_symmetry)
742 : CALL dbcsr_create(cmatrix_full, template=matrix_s(1, 1)%matrix, &
743 32 : matrix_type=dbcsr_type_no_symmetry)
744 32 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
745 32 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
746 : !
747 32 : CALL get_kpoint_info(kpoint=kpoint, cell_to_index=cell_to_index)
748 32 : NULLIFY (fmdummy)
749 32 : nspins = dft_control%nspins
750 64 : DO ispin = 1, nspins
751 : ! loop over all k-points
752 544 : DO ik = 1, num_kpts
753 : ! get the MO coefficients for this k-point
754 480 : my_kpgrp = (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2))
755 : IF (my_kpgrp) THEN
756 480 : ikk = ik - kpoint%kp_range(1) + 1
757 480 : kp => kpoint%kp_env(ikk)%kpoint_env
758 480 : CPASSERT(SIZE(kp%mos, 1) == 2)
759 480 : fmr => kp%mos(1, ispin)%mo_coeff
760 480 : fmi => kp%mos(2, ispin)%mo_coeff
761 480 : CALL cp_fm_copy_general(fmr, fmk1(1), para_env)
762 480 : CALL cp_fm_copy_general(fmi, fmk1(2), para_env)
763 : ELSE
764 0 : NULLIFY (fmr, fmi, kp)
765 0 : CALL cp_fm_copy_general(fmdummy, fmk1(1), para_env)
766 0 : CALL cp_fm_copy_general(fmdummy, fmk1(2), para_env)
767 : END IF
768 480 : CALL cp_fm_to_cfm(fmk1(1), fmk1(2), fmk1_cfm)
769 : ! loop over all connected neighbors
770 3392 : DO i = 1, nntot
771 : ! get the MO coefficients for the connected k-point
772 2880 : ik2 = nnlist(ik, i)
773 2880 : mygrp = (ik2 >= kpoint%kp_range(1) .AND. ik2 <= kpoint%kp_range(2))
774 : IF (mygrp) THEN
775 2880 : ikk = ik2 - kpoint%kp_range(1) + 1
776 2880 : kp => kpoint%kp_env(ikk)%kpoint_env
777 2880 : CPASSERT(SIZE(kp%mos, 1) == 2)
778 2880 : fmr => kp%mos(1, ispin)%mo_coeff
779 2880 : fmi => kp%mos(2, ispin)%mo_coeff
780 2880 : CALL cp_fm_copy_general(fmr, fmk2(1), para_env)
781 2880 : CALL cp_fm_copy_general(fmi, fmk2(2), para_env)
782 : ELSE
783 0 : NULLIFY (fmr, fmi, kp)
784 0 : CALL cp_fm_copy_general(fmdummy, fmk2(1), para_env)
785 0 : CALL cp_fm_copy_general(fmdummy, fmk2(2), para_env)
786 : END IF
787 2880 : CALL cp_fm_to_cfm(fmk2(1), fmk2(2), fmk2_cfm)
788 : !
789 : ! transfer realspace overlaps to connected k-point
790 2880 : ibs = nblist(ik, i)
791 2880 : ksign = SIGN(1.0_dp, REAL(ibs, KIND=dp))
792 2880 : ibs = ABS(ibs)
793 2880 : CALL dbcsr_set(rmatrix, 0.0_dp)
794 2880 : CALL dbcsr_set(cmatrix, 0.0_dp)
795 : CALL rskp_transform(rmatrix, cmatrix, rsmat=berry_matrix(ibs)%cosmat, ispin=1, &
796 : xkp=kpoint%xkp(1:3, ik2), cell_to_index=cell_to_index, sab_nl=sab_nl, &
797 2880 : is_complex=.FALSE., rs_sign=ksign)
798 : CALL rskp_transform(cmatrix, rmatrix, rsmat=berry_matrix(ibs)%sinmat, ispin=1, &
799 : xkp=kpoint%xkp(1:3, ik2), cell_to_index=cell_to_index, sab_nl=sab_nl, &
800 2880 : is_complex=.TRUE., rs_sign=ksign)
801 : !
802 : ! calculate M_(mn)^(k,b) = C(k)^H O(k,b) C(k+b)
803 2880 : CALL dbcsr_desymmetrize(rmatrix, rmatrix_full)
804 2880 : CALL dbcsr_desymmetrize(cmatrix, cmatrix_full)
805 2880 : CALL copy_dbcsr_to_fm(rmatrix_full, mat_real)
806 2880 : CALL copy_dbcsr_to_fm(cmatrix_full, mat_imag)
807 2880 : CALL cp_fm_to_cfm(mat_real, mat_imag, omat_cfm)
808 : CALL cp_cfm_gemm("N", "N", nao, nmo, nao, CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
809 2880 : omat_cfm, fmk2_cfm, CMPLX(0.0_dp, 0.0_dp, KIND=dp), tmp_cfm)
810 : CALL cp_cfm_gemm("C", "N", nmo, nmo, nao, CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
811 2880 : fmk1_cfm, tmp_cfm, CMPLX(0.0_dp, 0.0_dp, KIND=dp), mmn_cfm)
812 2880 : CALL cp_cfm_to_fm(mmn_cfm, mmn_real, mmn_imag)
813 : !
814 : ! write to output file
815 2880 : IF (reused_scf_mos .AND. use_bloch_phases) THEN
816 : ! Reused SCF MOs need the same global Bloch gauge in every overlap block.
817 : gauge_arg = twopi*DOT_PRODUCT(kpoint%xkp(1:3, ik2) - kpoint%xkp(1:3, ik), &
818 7680 : phase_center(1:3))
819 1920 : gauge_real = COS(gauge_arg)
820 1920 : gauge_imag = SIN(gauge_arg)
821 : ELSE
822 : gauge_real = 1.0_dp
823 : gauge_imag = 0.0_dp
824 : END IF
825 2880 : IF (para_env%is_source()) THEN
826 1440 : WRITE (iunit, "(2I8,3I5)") ik, ik2, nncell(1:3, ik, i)
827 : END IF
828 13632 : DO ib2 = 1, nmo
829 52608 : DO ib1 = 1, nmo
830 39456 : CALL cp_fm_get_element(mmn_real, ib1, ib2, rmmn)
831 39456 : CALL cp_fm_get_element(mmn_imag, ib1, ib2, cmmn)
832 39456 : gauge_tmp = gauge_real*rmmn - gauge_imag*cmmn
833 39456 : cmmn = gauge_imag*rmmn + gauge_real*cmmn
834 39456 : rmmn = gauge_tmp
835 49728 : IF (para_env%is_source()) THEN
836 19728 : WRITE (iunit, "(2E30.14)") rmmn, cmmn
837 : END IF
838 : END DO
839 : END DO
840 : !
841 : END DO
842 : END DO
843 : END DO
844 128 : DO i = 1, nbs
845 96 : CALL dbcsr_deallocate_matrix_set(berry_matrix(i)%cosmat)
846 128 : CALL dbcsr_deallocate_matrix_set(berry_matrix(i)%sinmat)
847 : END DO
848 32 : DEALLOCATE (berry_matrix)
849 32 : CALL cp_fm_struct_release(matrix_struct_work)
850 96 : DO i = 1, 2
851 64 : CALL cp_fm_release(fmk1(i))
852 96 : CALL cp_fm_release(fmk2(i))
853 : END DO
854 32 : CALL cp_cfm_release(fmk1_cfm)
855 32 : CALL cp_cfm_release(fmk2_cfm)
856 32 : CALL cp_cfm_release(tmp_cfm)
857 32 : CALL cp_fm_struct_release(matrix_struct_ao)
858 32 : CALL cp_fm_release(mat_real)
859 32 : CALL cp_fm_release(mat_imag)
860 32 : CALL cp_cfm_release(omat_cfm)
861 32 : CALL cp_fm_struct_release(matrix_struct_mmn)
862 32 : CALL cp_fm_release(mmn_real)
863 32 : CALL cp_fm_release(mmn_imag)
864 32 : CALL cp_cfm_release(mmn_cfm)
865 32 : CALL dbcsr_deallocate_matrix(rmatrix)
866 32 : CALL dbcsr_deallocate_matrix(cmatrix)
867 32 : CALL dbcsr_deallocate_matrix(rmatrix_full)
868 32 : CALL dbcsr_deallocate_matrix(cmatrix_full)
869 : !
870 32 : IF (para_env%is_source()) THEN
871 16 : CALL close_file(iunit)
872 : END IF
873 : !
874 : ! Calculate and print Projections
875 : !
876 : ! Print eigenvalues
877 32 : nspins = dft_control%nspins
878 32 : kp => kpoint%kp_env(1)%kpoint_env
879 32 : CALL get_mo_set(kp%mos(1, 1), nmo=nmo)
880 96 : ALLOCATE (eigval(nmo))
881 32 : CALL get_kpoint_info(kpoint, nkp=nkp, kp_range=kp_range, xkp=xkp)
882 32 : IF (para_env%is_source()) THEN
883 16 : WRITE (filename, '(A,A)') TRIM(seed_name), ".eig"
884 16 : CALL open_file(filename, unit_number=iunit, file_status="UNKNOWN", file_action="WRITE")
885 : ELSE
886 16 : iunit = -1
887 : END IF
888 : !
889 512 : DO ik = 1, nkp
890 480 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
891 992 : DO ispin = 1, nspins
892 480 : IF (my_kpgrp) THEN
893 480 : ikpgr = ik - kp_range(1) + 1
894 480 : kp => kpoint%kp_env(ikpgr)%kpoint_env
895 480 : CALL get_mo_set(kp%mos(1, ispin), eigenvalues=eigenvalues)
896 2192 : eigval(1:nmo) = eigenvalues(1:nmo)
897 : ELSE
898 0 : eigval(1:nmo) = 0.0_dp
899 : END IF
900 480 : CALL kpoint%para_env_inter_kp%sum(eigval)
901 2192 : eigval(1:nmo) = eigval(1:nmo)*evolt
902 : ! output
903 960 : IF (iunit > 0) THEN
904 1096 : DO ib = 1, nmo
905 1096 : WRITE (iunit, "(2I8,F24.14)") ib, ik, eigval(ib)
906 : END DO
907 : END IF
908 : END DO
909 : END DO
910 32 : IF (para_env%is_source()) THEN
911 16 : CALL close_file(iunit)
912 : END IF
913 : !
914 : ! clean up
915 32 : DEALLOCATE (kpt_latt, atoms_cart, atom_symbols, eigval)
916 32 : DEALLOCATE (nnlist, nncell)
917 32 : DEALLOCATE (nblist, b_latt)
918 32 : IF (nexcl > 0) THEN
919 0 : DEALLOCATE (exclude_bands)
920 : END IF
921 32 : IF (do_kpoints) THEN
922 32 : NULLIFY (qs_env_kp)
923 : ELSE
924 0 : CALL qs_env_release(qs_env_kp)
925 0 : DEALLOCATE (qs_env_kp)
926 : NULLIFY (qs_env_kp)
927 : END IF
928 :
929 32 : CALL kpoint_release(kpoint)
930 :
931 352 : END SUBROUTINE wannier90_files
932 :
933 : ! **************************************************************************************************
934 : !> \brief Reconstruct a full Wannier90 k-point MO set from the SCF k-point MOs.
935 : !> \param kpoint full Wannier90 export k-point object
936 : !> \param qs_kpoint SCF k-point object
937 : !> \param matrix_s real-space overlap matrix
938 : !> \param matrix_ks real-space Kohn-Sham matrix
939 : !> \param cell_to_index real-space cell index table
940 : !> \param sab_nl overlap neighbor list
941 : !> \param para_env global parallel environment
942 : !> \param success true if all full-mesh MOs were reconstructed
943 : !> \param reason diagnostic message when reconstruction is not possible
944 : !> \param aligned_degenerate_blocks number of aligned degenerate MO blocks
945 : !> \param aligned_degenerate_max_size largest aligned degenerate MO block
946 : !> \param aligned_degenerate_min_svalue smallest S(k)-metric subspace singular value
947 : ! **************************************************************************************************
948 36 : SUBROUTINE prepare_wannier90_scf_mos(kpoint, qs_kpoint, matrix_s, matrix_ks, cell_to_index, &
949 : sab_nl, para_env, success, reason, aligned_degenerate_blocks, &
950 : aligned_degenerate_max_size, &
951 : aligned_degenerate_min_svalue)
952 : TYPE(kpoint_type), POINTER :: kpoint, qs_kpoint
953 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_ks
954 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
955 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
956 : POINTER :: sab_nl
957 : TYPE(mp_para_env_type), POINTER :: para_env
958 : LOGICAL, INTENT(OUT) :: success
959 : CHARACTER(LEN=*), INTENT(OUT) :: reason
960 : INTEGER, INTENT(OUT) :: aligned_degenerate_blocks, &
961 : aligned_degenerate_max_size
962 : REAL(KIND=dp), INTENT(OUT) :: aligned_degenerate_min_svalue
963 :
964 : CHARACTER(LEN=default_string_length) :: best_reason, candidate_reason
965 : INTEGER :: aligned_blocks, aligned_max_size, candidate_aligned_blocks, &
966 : candidate_aligned_max_size, ik, ikpgr, ikred, ispin, isym_try, min_gap_band, &
967 : min_gap_kpoint, min_gap_spin, nao, nao_src, nmo, nmo_src, nspins, nsymmetry, &
968 : num_candidates
969 36 : INTEGER, ALLOCATABLE, DIMENSION(:) :: source_kpoint, sym_index
970 : INTEGER, DIMENSION(2) :: kp_range, source_kp_range
971 : LOGICAL :: my_kpgrp, my_source_kpgrp, ok, &
972 : source_window
973 : REAL(KIND=dp) :: aligned_min_svalue, band_gap, best_residual, candidate_residual, &
974 : degenerate_band_tol, local_min_band_gap, min_band_gap, source_owner_count, &
975 : source_window_min_svalue
976 36 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_buffer, occupation_buffer, &
977 36 : source_eigenvalues_buffer, &
978 36 : source_occupation_buffer
979 36 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation
980 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
981 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_source, matrix_struct_work
982 : TYPE(cp_fm_type) :: dst_imag, dst_imag_full, dst_real, &
983 : dst_real_full, src_imag, &
984 : src_imag_full, src_real, src_real_full
985 : TYPE(cp_fm_type), POINTER :: dst_fmi, dst_fmr, src_fmi, src_fmr
986 : TYPE(kpoint_env_type), POINTER :: kp, kp_source
987 : TYPE(kpoint_sym_type), POINTER :: kpsym
988 :
989 36 : success = .FALSE.
990 36 : reason = ""
991 36 : aligned_degenerate_blocks = 0
992 36 : aligned_degenerate_max_size = 0
993 36 : aligned_degenerate_min_svalue = HUGE(1.0_dp)
994 36 : NULLIFY (matrix_struct_source, matrix_struct_work, src_fmr, src_fmi, dst_fmr, dst_fmi)
995 :
996 36 : IF (.NOT. ASSOCIATED(kpoint)) THEN
997 0 : reason = "internal Wannier90 k-point object is not available"
998 0 : RETURN
999 : END IF
1000 36 : IF (.NOT. ASSOCIATED(qs_kpoint)) THEN
1001 0 : reason = "SCF k-point object is not available"
1002 0 : RETURN
1003 : END IF
1004 36 : IF (.NOT. ASSOCIATED(kpoint%kp_env) .OR. .NOT. ASSOCIATED(qs_kpoint%kp_env)) THEN
1005 0 : reason = "k-point MO environments are not initialized"
1006 0 : RETURN
1007 : END IF
1008 36 : IF (.NOT. ASSOCIATED(kpoint%blacs_env)) THEN
1009 0 : reason = "Wannier90 k-point BLACS environment is not initialized"
1010 0 : RETURN
1011 : END IF
1012 :
1013 36 : CALL build_wannier90_scf_mapping(kpoint, qs_kpoint, source_kpoint, sym_index, ok, reason)
1014 36 : IF (.NOT. ok) RETURN
1015 536 : nsymmetry = COUNT(sym_index > 0)
1016 :
1017 36 : kp => kpoint%kp_env(1)%kpoint_env
1018 36 : nspins = SIZE(kp%mos, 2)
1019 36 : IF (SIZE(kp%mos, 1) < 2) THEN
1020 0 : reason = "Wannier90 export k-point MOs are not complex-valued"
1021 0 : DEALLOCATE (source_kpoint, sym_index)
1022 0 : RETURN
1023 : END IF
1024 36 : CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo)
1025 :
1026 36 : kp_source => qs_kpoint%kp_env(1)%kpoint_env
1027 36 : IF (SIZE(kp_source%mos, 1) < 2) THEN
1028 0 : reason = "SCF MOs are real-valued; complex symmetry phases cannot be reconstructed"
1029 0 : DEALLOCATE (source_kpoint, sym_index)
1030 0 : RETURN
1031 : END IF
1032 36 : CALL get_mo_set(kp_source%mos(1, 1), nao=nao_src, nmo=nmo_src)
1033 36 : CALL para_env%max(nao_src)
1034 36 : CALL para_env%max(nmo_src)
1035 36 : IF (nao_src /= nao) THEN
1036 0 : reason = "SCF and Wannier90 MO bases have different AO dimensions"
1037 0 : DEALLOCATE (source_kpoint, sym_index)
1038 0 : RETURN
1039 : END IF
1040 36 : IF (nmo_src < nmo) THEN
1041 0 : reason = "SCF MO set has fewer bands than the Wannier90 export"
1042 0 : DEALLOCATE (source_kpoint, sym_index)
1043 0 : RETURN
1044 : END IF
1045 36 : source_window = nmo_src > nmo
1046 36 : degenerate_band_tol = 1.0e-8_dp
1047 36 : CALL get_kpoint_info(qs_kpoint, kp_range=source_kp_range)
1048 36 : IF (source_kp_range(1) /= 1 .OR. source_kp_range(2) /= qs_kpoint%nkp) THEN
1049 6 : reason = "SCF k-point symmetry data are distributed over k-point parallel groups"
1050 6 : DEALLOCATE (source_kpoint, sym_index)
1051 6 : RETURN
1052 : END IF
1053 : ! Positive symmetry entries require atom/AO rotations and Bloch phases. Degenerate subspaces
1054 : ! fully contained in the exported band window are aligned below; only guard when the Wannier90
1055 : ! window cuts through a degenerate SCF manifold at the upper band edge.
1056 30 : IF (nsymmetry > 0 .AND. nmo_src > nmo) THEN
1057 0 : local_min_band_gap = HUGE(1.0_dp)
1058 0 : min_gap_band = nmo
1059 0 : min_gap_kpoint = 0
1060 0 : min_gap_spin = 0
1061 0 : DO ikred = source_kp_range(1), source_kp_range(2)
1062 0 : ikpgr = ikred - source_kp_range(1) + 1
1063 0 : kp_source => qs_kpoint%kp_env(ikpgr)%kpoint_env
1064 0 : DO ispin = 1, nspins
1065 0 : CALL get_mo_set(kp_source%mos(1, ispin), eigenvalues=eigenvalues)
1066 0 : band_gap = ABS(eigenvalues(nmo + 1) - eigenvalues(nmo))
1067 0 : IF (band_gap < local_min_band_gap) THEN
1068 0 : local_min_band_gap = band_gap
1069 0 : min_gap_band = nmo
1070 0 : min_gap_kpoint = ikred
1071 0 : min_gap_spin = ispin
1072 : END IF
1073 : END DO
1074 : END DO
1075 0 : min_band_gap = local_min_band_gap
1076 0 : CALL para_env%min(min_band_gap)
1077 0 : IF (ABS(local_min_band_gap - min_band_gap) > degenerate_band_tol*EPSILON(1.0_dp)) THEN
1078 0 : min_gap_kpoint = 0
1079 0 : min_gap_spin = 0
1080 : END IF
1081 0 : CALL para_env%max(min_gap_kpoint)
1082 0 : CALL para_env%max(min_gap_spin)
1083 0 : CALL para_env%max(min_gap_band)
1084 0 : IF (min_band_gap < degenerate_band_tol) THEN
1085 : WRITE (reason, "(A,ES9.2,A,I0,A,I0,A,I0)") &
1086 0 : "degenerate atom/AO W90 reuse guarded: edge gap=", min_band_gap, ", k=", &
1087 0 : min_gap_kpoint, ", s=", min_gap_spin, ", nband=", min_gap_band
1088 0 : DEALLOCATE (source_kpoint, sym_index)
1089 0 : RETURN
1090 : END IF
1091 : END IF
1092 30 : blacs_env => kpoint%blacs_env
1093 : CALL cp_fm_struct_create(matrix_struct_work, nrow_global=nao, ncol_global=nmo, &
1094 30 : para_env=para_env, context=blacs_env)
1095 30 : CALL cp_fm_create(src_real, matrix_struct_work)
1096 30 : CALL cp_fm_create(src_imag, matrix_struct_work)
1097 30 : CALL cp_fm_create(dst_real, matrix_struct_work)
1098 30 : CALL cp_fm_create(dst_imag, matrix_struct_work)
1099 30 : IF (source_window) THEN
1100 : CALL cp_fm_struct_create(matrix_struct_source, nrow_global=nao, ncol_global=nmo_src, &
1101 0 : para_env=para_env, context=blacs_env)
1102 0 : CALL cp_fm_create(src_real_full, matrix_struct_source)
1103 0 : CALL cp_fm_create(src_imag_full, matrix_struct_source)
1104 0 : CALL cp_fm_create(dst_real_full, matrix_struct_source)
1105 0 : CALL cp_fm_create(dst_imag_full, matrix_struct_source)
1106 : END IF
1107 0 : ALLOCATE (eigenvalues_buffer(nmo), occupation_buffer(nmo), &
1108 210 : source_eigenvalues_buffer(nmo_src), source_occupation_buffer(nmo_src))
1109 :
1110 30 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
1111 30 : CALL get_kpoint_info(qs_kpoint, kp_range=source_kp_range)
1112 :
1113 482 : DO ik = 1, kpoint%nkp
1114 452 : ikred = source_kpoint(ik)
1115 452 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
1116 452 : my_source_kpgrp = (ikred >= source_kp_range(1) .AND. ikred <= source_kp_range(2))
1117 934 : DO ispin = 1, nspins
1118 2244 : source_eigenvalues_buffer(1:nmo_src) = 0.0_dp
1119 2244 : source_occupation_buffer(1:nmo_src) = 0.0_dp
1120 452 : IF (my_source_kpgrp) THEN
1121 452 : ikpgr = ikred - source_kp_range(1) + 1
1122 452 : kp_source => qs_kpoint%kp_env(ikpgr)%kpoint_env
1123 452 : src_fmr => kp_source%mos(1, ispin)%mo_coeff
1124 452 : src_fmi => kp_source%mos(2, ispin)%mo_coeff
1125 : CALL get_mo_set(kp_source%mos(1, ispin), eigenvalues=eigenvalues, &
1126 452 : occupation_numbers=occupation)
1127 2244 : source_eigenvalues_buffer(1:nmo_src) = eigenvalues(1:nmo_src)
1128 2244 : source_occupation_buffer(1:nmo_src) = occupation(1:nmo_src)
1129 : ELSE
1130 : NULLIFY (src_fmr, src_fmi)
1131 : END IF
1132 : IF (my_source_kpgrp) THEN
1133 452 : source_owner_count = 1.0_dp
1134 : ELSE
1135 0 : source_owner_count = 0.0_dp
1136 : END IF
1137 452 : CALL para_env%sum(source_owner_count)
1138 452 : CALL para_env%sum(source_eigenvalues_buffer)
1139 452 : CALL para_env%sum(source_occupation_buffer)
1140 452 : IF (source_owner_count > 0.0_dp) THEN
1141 : source_eigenvalues_buffer(1:nmo_src) = &
1142 2244 : source_eigenvalues_buffer(1:nmo_src)/source_owner_count
1143 : source_occupation_buffer(1:nmo_src) = &
1144 2244 : source_occupation_buffer(1:nmo_src)/source_owner_count
1145 : END IF
1146 2244 : eigenvalues_buffer(1:nmo) = source_eigenvalues_buffer(1:nmo)
1147 2244 : occupation_buffer(1:nmo) = source_occupation_buffer(1:nmo)
1148 452 : IF (source_window) THEN
1149 0 : CALL cp_fm_copy_general(src_fmr, src_real_full, para_env)
1150 0 : CALL cp_fm_copy_general(src_fmi, src_imag_full, para_env)
1151 0 : CALL copy_wannier90_mo_window(src_real_full, src_real, nmo)
1152 0 : CALL copy_wannier90_mo_window(src_imag_full, src_imag, nmo)
1153 : ELSE
1154 452 : CALL cp_fm_copy_general(src_fmr, src_real, para_env)
1155 452 : CALL cp_fm_copy_general(src_fmi, src_imag, para_env)
1156 : END IF
1157 :
1158 452 : ok = .FALSE.
1159 452 : reason = ""
1160 452 : aligned_blocks = 0
1161 452 : aligned_max_size = 0
1162 452 : aligned_min_svalue = 0.0_dp
1163 452 : IF (sym_index(ik) > 0) THEN
1164 368 : kpsym => qs_kpoint%kp_sym(ikred)%kpoint_sym
1165 368 : IF (ASSOCIATED(kpsym)) THEN
1166 368 : best_reason = ""
1167 368 : best_residual = HUGE(1.0_dp)
1168 368 : num_candidates = 0
1169 : ! Little-group operations can reach the same target k-point; keep the first valid eigenspace.
1170 3064 : DO isym_try = 1, kpsym%nwred
1171 3064 : IF (.NOT. same_periodic_kpoint(kpoint%xkp(1:3, ik), &
1172 : kpsym%xkp(1:3, isym_try))) CYCLE
1173 368 : num_candidates = num_candidates + 1
1174 : CALL transform_wannier90_scf_mo(src_real, src_imag, dst_real, dst_imag, &
1175 : qs_kpoint, ikred, isym_try, para_env, ok, &
1176 368 : candidate_reason)
1177 368 : IF (.NOT. ok) THEN
1178 0 : reason = candidate_reason
1179 : CYCLE
1180 : END IF
1181 : CALL ritz_stabilize_wannier90_subspace(dst_real, dst_imag, matrix_s, matrix_ks, &
1182 : kpoint%xkp(1:3, ik), cell_to_index, &
1183 : sab_nl, ispin, eigenvalues_buffer, &
1184 : degenerate_band_tol, ok, candidate_reason, &
1185 : candidate_aligned_blocks, &
1186 : candidate_aligned_max_size, &
1187 368 : aligned_min_svalue, candidate_residual)
1188 368 : IF (candidate_residual < best_residual) THEN
1189 368 : best_residual = candidate_residual
1190 368 : best_reason = candidate_reason
1191 : END IF
1192 368 : IF (.NOT. ok) THEN
1193 0 : IF (source_window) THEN
1194 : CALL transform_wannier90_scf_mo(src_real_full, src_imag_full, &
1195 : dst_real_full, dst_imag_full, qs_kpoint, &
1196 0 : ikred, isym_try, para_env, ok, candidate_reason)
1197 0 : IF (ok) THEN
1198 : CALL ritz_reconstruct_wannier90_window(dst_real_full, dst_imag_full, &
1199 : dst_real, dst_imag, matrix_s, &
1200 : matrix_ks, kpoint%xkp(1:3, ik), &
1201 : cell_to_index, sab_nl, ispin, &
1202 : eigenvalues_buffer, nmo, ok, &
1203 : candidate_reason, source_window_min_svalue, &
1204 0 : candidate_residual)
1205 0 : IF (candidate_residual < best_residual) THEN
1206 0 : best_residual = candidate_residual
1207 0 : best_reason = candidate_reason
1208 : END IF
1209 : END IF
1210 : ELSE
1211 : CALL ritz_reconstruct_wannier90_window(dst_real, dst_imag, dst_real, &
1212 : dst_imag, matrix_s, matrix_ks, &
1213 : kpoint%xkp(1:3, ik), cell_to_index, &
1214 : sab_nl, ispin, eigenvalues_buffer, nmo, &
1215 : ok, candidate_reason, source_window_min_svalue, &
1216 0 : candidate_residual)
1217 0 : IF (candidate_residual < best_residual) THEN
1218 0 : best_residual = candidate_residual
1219 0 : best_reason = candidate_reason
1220 : END IF
1221 : END IF
1222 : END IF
1223 368 : IF (ok) THEN
1224 368 : aligned_blocks = candidate_aligned_blocks
1225 368 : aligned_max_size = candidate_aligned_max_size
1226 368 : sym_index(ik) = isym_try
1227 368 : EXIT
1228 : END IF
1229 368 : reason = candidate_reason
1230 : END DO
1231 368 : IF (.NOT. ok .AND. num_candidates > 0 .AND. best_residual < HUGE(1.0_dp)) THEN
1232 : WRITE (reason, "(A,I0,A,ES9.2,A,I0,A,A32)") &
1233 0 : "atom/AO W90 guarded: best/", num_candidates, "=", best_residual, &
1234 0 : " k=", ik, " ", TRIM(best_reason)
1235 368 : ELSEIF (.NOT. ok .AND. num_candidates == 0) THEN
1236 0 : reason = "no matching SCF symmetry operation candidate"
1237 : END IF
1238 : ELSE
1239 0 : reason = "SCF k-point symmetry operation is not available"
1240 : END IF
1241 : ELSE
1242 : CALL transform_wannier90_scf_mo(src_real, src_imag, dst_real, dst_imag, qs_kpoint, &
1243 84 : ikred, sym_index(ik), para_env, ok, reason)
1244 : END IF
1245 452 : IF (ok .AND. sym_index(ik) <= 0) THEN
1246 : ! Even a direct k-point copy must be a closed H(k),S(k) subspace. This catches
1247 : ! incomplete degenerate band windows before they can be exported to Wannier90.
1248 : CALL ritz_stabilize_wannier90_subspace(dst_real, dst_imag, matrix_s, matrix_ks, &
1249 : kpoint%xkp(1:3, ik), cell_to_index, sab_nl, &
1250 : ispin, eigenvalues_buffer, degenerate_band_tol, &
1251 : ok, reason, aligned_blocks, aligned_max_size, &
1252 84 : aligned_min_svalue, candidate_residual)
1253 84 : IF (.NOT. ok .AND. source_window) THEN
1254 : CALL transform_wannier90_scf_mo(src_real_full, src_imag_full, dst_real_full, &
1255 : dst_imag_full, qs_kpoint, ikred, sym_index(ik), &
1256 0 : para_env, ok, reason)
1257 0 : IF (ok) THEN
1258 : CALL ritz_reconstruct_wannier90_window(dst_real_full, dst_imag_full, dst_real, &
1259 : dst_imag, matrix_s, matrix_ks, &
1260 : kpoint%xkp(1:3, ik), cell_to_index, &
1261 : sab_nl, ispin, eigenvalues_buffer, nmo, ok, &
1262 : reason, source_window_min_svalue, &
1263 0 : candidate_residual)
1264 : END IF
1265 84 : ELSEIF (.NOT. ok) THEN
1266 : CALL ritz_reconstruct_wannier90_window(dst_real, dst_imag, dst_real, dst_imag, &
1267 : matrix_s, matrix_ks, kpoint%xkp(1:3, ik), &
1268 : cell_to_index, sab_nl, ispin, &
1269 : eigenvalues_buffer, nmo, ok, reason, &
1270 0 : source_window_min_svalue, candidate_residual)
1271 : END IF
1272 : END IF
1273 452 : IF (.NOT. ok) THEN
1274 0 : CALL cp_fm_release(src_real)
1275 0 : CALL cp_fm_release(src_imag)
1276 0 : CALL cp_fm_release(dst_real)
1277 0 : CALL cp_fm_release(dst_imag)
1278 0 : CALL cp_fm_struct_release(matrix_struct_work)
1279 0 : IF (source_window) THEN
1280 0 : CALL cp_fm_release(src_real_full)
1281 0 : CALL cp_fm_release(src_imag_full)
1282 0 : CALL cp_fm_release(dst_real_full)
1283 0 : CALL cp_fm_release(dst_imag_full)
1284 0 : CALL cp_fm_struct_release(matrix_struct_source)
1285 : END IF
1286 0 : DEALLOCATE (source_kpoint, sym_index, eigenvalues_buffer, occupation_buffer, &
1287 0 : source_eigenvalues_buffer, source_occupation_buffer)
1288 0 : RETURN
1289 : END IF
1290 452 : IF (sym_index(ik) /= 0) THEN
1291 410 : aligned_degenerate_blocks = aligned_degenerate_blocks + aligned_blocks
1292 410 : aligned_degenerate_max_size = MAX(aligned_degenerate_max_size, aligned_max_size)
1293 410 : IF (aligned_blocks > 0) THEN
1294 338 : aligned_degenerate_min_svalue = MIN(aligned_degenerate_min_svalue, aligned_min_svalue)
1295 : END IF
1296 : END IF
1297 :
1298 452 : IF (my_kpgrp) THEN
1299 452 : ikpgr = ik - kp_range(1) + 1
1300 452 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1301 452 : dst_fmr => kp%mos(1, ispin)%mo_coeff
1302 452 : dst_fmi => kp%mos(2, ispin)%mo_coeff
1303 : CALL get_mo_set(kp%mos(1, ispin), eigenvalues=eigenvalues, &
1304 452 : occupation_numbers=occupation)
1305 2244 : eigenvalues(1:nmo) = eigenvalues_buffer(1:nmo)
1306 2244 : occupation(1:nmo) = occupation_buffer(1:nmo)
1307 : CALL get_mo_set(kp%mos(2, ispin), eigenvalues=eigenvalues, &
1308 452 : occupation_numbers=occupation)
1309 2244 : IF (ASSOCIATED(eigenvalues)) eigenvalues(1:nmo) = eigenvalues_buffer(1:nmo)
1310 2244 : IF (ASSOCIATED(occupation)) occupation(1:nmo) = occupation_buffer(1:nmo)
1311 : ELSE
1312 : NULLIFY (dst_fmr, dst_fmi)
1313 : END IF
1314 452 : CALL cp_fm_copy_general(dst_real, dst_fmr, para_env)
1315 904 : CALL cp_fm_copy_general(dst_imag, dst_fmi, para_env)
1316 : END DO
1317 : END DO
1318 :
1319 30 : CALL cp_fm_release(src_real)
1320 30 : CALL cp_fm_release(src_imag)
1321 30 : CALL cp_fm_release(dst_real)
1322 30 : CALL cp_fm_release(dst_imag)
1323 30 : CALL cp_fm_struct_release(matrix_struct_work)
1324 30 : IF (source_window) THEN
1325 0 : CALL cp_fm_release(src_real_full)
1326 0 : CALL cp_fm_release(src_imag_full)
1327 0 : CALL cp_fm_release(dst_real_full)
1328 0 : CALL cp_fm_release(dst_imag_full)
1329 0 : CALL cp_fm_struct_release(matrix_struct_source)
1330 : END IF
1331 0 : DEALLOCATE (source_kpoint, sym_index, eigenvalues_buffer, occupation_buffer, &
1332 30 : source_eigenvalues_buffer, source_occupation_buffer)
1333 30 : IF (aligned_degenerate_blocks == 0) aligned_degenerate_min_svalue = 0.0_dp
1334 30 : success = .TRUE.
1335 :
1336 174 : END SUBROUTINE prepare_wannier90_scf_mos
1337 :
1338 : ! **************************************************************************************************
1339 : !> \brief Save a full-mesh Wannier90 MO reference on all ranks for diagnostic validation.
1340 : !> \param kpoint full Wannier90 export k-point object
1341 : !> \param nspins number of spin channels
1342 : !> \param para_env global parallel environment
1343 : !> \param mo_real real MO coefficients, indexed as AO, MO, k-point, spin
1344 : !> \param mo_imag imaginary MO coefficients, indexed as AO, MO, k-point, spin
1345 : !> \param eigenvalue_snapshot MO eigenvalues, indexed as MO, k-point, spin
1346 : ! **************************************************************************************************
1347 6 : SUBROUTINE save_wannier90_mo_snapshot(kpoint, nspins, para_env, mo_real, mo_imag, &
1348 : eigenvalue_snapshot)
1349 : TYPE(kpoint_type), POINTER :: kpoint
1350 : INTEGER, INTENT(IN) :: nspins
1351 : TYPE(mp_para_env_type), POINTER :: para_env
1352 : REAL(KIND=dp), ALLOCATABLE, &
1353 : DIMENSION(:, :, :, :), INTENT(OUT) :: mo_real, mo_imag
1354 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1355 : INTENT(OUT) :: eigenvalue_snapshot
1356 :
1357 : INTEGER :: ik, ikpgr, ispin, nao, nkp, nmo
1358 : INTEGER, DIMENSION(2) :: kp_range
1359 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: owner_weight
1360 6 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1361 : TYPE(cp_fm_type), POINTER :: fmi, fmr
1362 : TYPE(kpoint_env_type), POINTER :: kp
1363 :
1364 6 : CALL get_kpoint_info(kpoint, nkp=nkp, kp_range=kp_range)
1365 6 : kp => kpoint%kp_env(1)%kpoint_env
1366 6 : CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo)
1367 0 : ALLOCATE (mo_real(nao, nmo, nkp, nspins), mo_imag(nao, nmo, nkp, nspins), &
1368 102 : eigenvalue_snapshot(nmo, nkp, nspins), owner_weight(nkp, nspins))
1369 6 : mo_real(:, :, :, :) = 0.0_dp
1370 6 : mo_imag(:, :, :, :) = 0.0_dp
1371 6 : eigenvalue_snapshot(:, :, :) = 0.0_dp
1372 6 : owner_weight(:, :) = 0.0_dp
1373 278 : DO ik = kp_range(1), kp_range(2)
1374 272 : ikpgr = ik - kp_range(1) + 1
1375 272 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1376 550 : DO ispin = 1, nspins
1377 272 : fmr => kp%mos(1, ispin)%mo_coeff
1378 272 : fmi => kp%mos(2, ispin)%mo_coeff
1379 272 : CALL cp_fm_get_submatrix(fmr, mo_real(:, :, ik, ispin))
1380 272 : CALL cp_fm_get_submatrix(fmi, mo_imag(:, :, ik, ispin))
1381 272 : CALL get_mo_set(kp%mos(1, ispin), eigenvalues=eigenvalues)
1382 1312 : eigenvalue_snapshot(1:nmo, ik, ispin) = eigenvalues(1:nmo)
1383 544 : owner_weight(ik, ispin) = 1.0_dp
1384 : END DO
1385 : END DO
1386 6 : CALL para_env%sum(mo_real)
1387 6 : CALL para_env%sum(mo_imag)
1388 6 : CALL para_env%sum(eigenvalue_snapshot)
1389 6 : CALL para_env%sum(owner_weight)
1390 278 : DO ik = 1, nkp
1391 550 : DO ispin = 1, nspins
1392 544 : IF (owner_weight(ik, ispin) > 0.0_dp) THEN
1393 42352 : mo_real(:, :, ik, ispin) = mo_real(:, :, ik, ispin)/owner_weight(ik, ispin)
1394 42352 : mo_imag(:, :, ik, ispin) = mo_imag(:, :, ik, ispin)/owner_weight(ik, ispin)
1395 : eigenvalue_snapshot(:, ik, ispin) = &
1396 1312 : eigenvalue_snapshot(:, ik, ispin)/owner_weight(ik, ispin)
1397 : END IF
1398 : END DO
1399 : END DO
1400 6 : DEALLOCATE (owner_weight)
1401 :
1402 6 : END SUBROUTINE save_wannier90_mo_snapshot
1403 :
1404 : ! **************************************************************************************************
1405 : !> \brief Restore a full-mesh Wannier90 MO reference after a failed diagnostic reuse attempt.
1406 : !> \param kpoint full Wannier90 export k-point object
1407 : !> \param mo_real real MO coefficient snapshot
1408 : !> \param mo_imag imaginary MO coefficient snapshot
1409 : !> \param eigenvalue_snapshot MO eigenvalue snapshot
1410 : ! **************************************************************************************************
1411 0 : SUBROUTINE restore_wannier90_mo_snapshot(kpoint, mo_real, mo_imag, eigenvalue_snapshot)
1412 : TYPE(kpoint_type), POINTER :: kpoint
1413 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: mo_real, mo_imag
1414 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: eigenvalue_snapshot
1415 :
1416 : INTEGER :: ik, ikpgr, ispin, nmo, nspins
1417 : INTEGER, DIMENSION(2) :: kp_range
1418 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1419 : TYPE(cp_fm_type), POINTER :: fmi, fmr
1420 : TYPE(kpoint_env_type), POINTER :: kp
1421 :
1422 0 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
1423 0 : nmo = SIZE(eigenvalue_snapshot, 1)
1424 0 : nspins = SIZE(eigenvalue_snapshot, 3)
1425 0 : DO ik = kp_range(1), kp_range(2)
1426 0 : ikpgr = ik - kp_range(1) + 1
1427 0 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1428 0 : DO ispin = 1, nspins
1429 0 : fmr => kp%mos(1, ispin)%mo_coeff
1430 0 : fmi => kp%mos(2, ispin)%mo_coeff
1431 0 : CALL cp_fm_set_submatrix(fmr, mo_real(:, :, ik, ispin))
1432 0 : CALL cp_fm_set_submatrix(fmi, mo_imag(:, :, ik, ispin))
1433 0 : CALL get_mo_set(kp%mos(1, ispin), eigenvalues=eigenvalues)
1434 0 : eigenvalues(1:nmo) = eigenvalue_snapshot(1:nmo, ik, ispin)
1435 0 : CALL get_mo_set(kp%mos(2, ispin), eigenvalues=eigenvalues)
1436 0 : IF (ASSOCIATED(eigenvalues)) eigenvalues(1:nmo) = eigenvalue_snapshot(1:nmo, ik, ispin)
1437 : END DO
1438 : END DO
1439 :
1440 0 : END SUBROUTINE restore_wannier90_mo_snapshot
1441 :
1442 : ! **************************************************************************************************
1443 : !> \brief Validate current Wannier90 MOs against a saved full-mesh diagonalization reference.
1444 : !> \param kpoint full Wannier90 export k-point object
1445 : !> \param matrix_s real-space overlap matrix
1446 : !> \param cell_to_index real-space cell index table
1447 : !> \param sab_nl overlap neighbor list
1448 : !> \param para_env global parallel environment
1449 : !> \param reference_mo_real real MO coefficient reference
1450 : !> \param reference_mo_imag imaginary MO coefficient reference
1451 : !> \param reference_eigenvalues MO eigenvalue reference
1452 : !> \param success true if the reconstructed MOs match the reference subspaces
1453 : !> \param max_subspace_deviation largest deviation of S(k)-metric singular values from one
1454 : !> \param min_svalue smallest S(k)-metric singular value
1455 : !> \param max_eigenvalue_deviation largest eigenvalue deviation
1456 : ! **************************************************************************************************
1457 6 : SUBROUTINE validate_wannier90_reused_mos(kpoint, matrix_s, cell_to_index, sab_nl, para_env, &
1458 6 : reference_mo_real, reference_mo_imag, reference_eigenvalues, &
1459 : success, max_subspace_deviation, min_svalue, &
1460 : max_eigenvalue_deviation)
1461 : TYPE(kpoint_type), POINTER :: kpoint
1462 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
1463 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1464 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1465 : POINTER :: sab_nl
1466 : TYPE(mp_para_env_type), POINTER :: para_env
1467 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: reference_mo_real, reference_mo_imag
1468 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: reference_eigenvalues
1469 : LOGICAL, INTENT(OUT) :: success
1470 : REAL(KIND=dp), INTENT(OUT) :: max_subspace_deviation, min_svalue, &
1471 : max_eigenvalue_deviation
1472 :
1473 : REAL(KIND=dp), PARAMETER :: eigenvalue_tol = 1.0e-8_dp, &
1474 : subspace_tol = 1.0e-4_dp
1475 :
1476 : INTEGER :: ik, ikpgr, ispin, nao, nkp, nmo, nspins
1477 : INTEGER, DIMENSION(2) :: kp_range
1478 : LOGICAL :: my_kpgrp, ok
1479 : REAL(KIND=dp) :: candidate_deviation, candidate_svalue, &
1480 : owner_count
1481 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalue_buffer
1482 6 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1483 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_work
1484 : TYPE(cp_fm_type) :: cand_imag, cand_real, ref_imag, ref_real
1485 : TYPE(cp_fm_type), POINTER :: fmi, fmr
1486 : TYPE(kpoint_env_type), POINTER :: kp
1487 :
1488 6 : success = .FALSE.
1489 6 : max_subspace_deviation = 0.0_dp
1490 6 : min_svalue = HUGE(1.0_dp)
1491 6 : max_eigenvalue_deviation = 0.0_dp
1492 6 : NULLIFY (matrix_struct_work, fmr, fmi)
1493 :
1494 6 : CALL get_kpoint_info(kpoint, nkp=nkp, kp_range=kp_range)
1495 6 : kp => kpoint%kp_env(1)%kpoint_env
1496 6 : nspins = SIZE(kp%mos, 2)
1497 6 : CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo)
1498 6 : CPASSERT(SIZE(reference_mo_real, 1) == nao)
1499 6 : CPASSERT(SIZE(reference_mo_real, 2) == nmo)
1500 6 : CPASSERT(SIZE(reference_mo_real, 3) == nkp)
1501 6 : CPASSERT(SIZE(reference_mo_real, 4) == nspins)
1502 :
1503 6 : CALL cp_fm_get_info(kp%mos(1, 1)%mo_coeff, matrix_struct=matrix_struct_work)
1504 6 : CALL cp_fm_create(ref_real, matrix_struct_work)
1505 6 : CALL cp_fm_create(ref_imag, matrix_struct_work)
1506 6 : CALL cp_fm_create(cand_real, matrix_struct_work)
1507 6 : CALL cp_fm_create(cand_imag, matrix_struct_work)
1508 18 : ALLOCATE (eigenvalue_buffer(nmo))
1509 :
1510 12 : DO ispin = 1, nspins
1511 284 : DO ik = 1, nkp
1512 272 : CALL cp_fm_set_submatrix(ref_real, reference_mo_real(:, :, ik, ispin))
1513 272 : CALL cp_fm_set_submatrix(ref_imag, reference_mo_imag(:, :, ik, ispin))
1514 272 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
1515 : IF (my_kpgrp) THEN
1516 272 : ikpgr = ik - kp_range(1) + 1
1517 272 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1518 272 : fmr => kp%mos(1, ispin)%mo_coeff
1519 272 : fmi => kp%mos(2, ispin)%mo_coeff
1520 272 : CALL get_mo_set(kp%mos(1, ispin), eigenvalues=eigenvalues)
1521 1312 : eigenvalue_buffer(1:nmo) = eigenvalues(1:nmo)
1522 : ELSE
1523 0 : NULLIFY (fmr, fmi)
1524 0 : eigenvalue_buffer(1:nmo) = 0.0_dp
1525 : END IF
1526 272 : CALL cp_fm_copy_general(fmr, cand_real, para_env)
1527 272 : CALL cp_fm_copy_general(fmi, cand_imag, para_env)
1528 272 : IF (my_kpgrp) THEN
1529 272 : owner_count = 1.0_dp
1530 : ELSE
1531 0 : owner_count = 0.0_dp
1532 : END IF
1533 272 : CALL para_env%sum(owner_count)
1534 272 : CALL para_env%sum(eigenvalue_buffer)
1535 1312 : IF (owner_count > 0.0_dp) eigenvalue_buffer(1:nmo) = eigenvalue_buffer(1:nmo)/owner_count
1536 : max_eigenvalue_deviation = MAX(max_eigenvalue_deviation, &
1537 : MAXVAL(ABS(eigenvalue_buffer(1:nmo) - &
1538 1312 : reference_eigenvalues(1:nmo, ik, ispin))))
1539 : CALL measure_wannier90_subspace_error(ref_real, ref_imag, cand_real, cand_imag, matrix_s, &
1540 : kpoint%xkp(1:3, ik), cell_to_index, sab_nl, para_env, &
1541 272 : ok, candidate_deviation, candidate_svalue)
1542 278 : IF (.NOT. ok) THEN
1543 0 : max_subspace_deviation = HUGE(1.0_dp)
1544 : ELSE
1545 272 : max_subspace_deviation = MAX(max_subspace_deviation, candidate_deviation)
1546 272 : min_svalue = MIN(min_svalue, candidate_svalue)
1547 : END IF
1548 : END DO
1549 : END DO
1550 6 : CALL para_env%max(max_subspace_deviation)
1551 6 : CALL para_env%min(min_svalue)
1552 6 : CALL para_env%max(max_eigenvalue_deviation)
1553 6 : success = max_subspace_deviation < subspace_tol .AND. max_eigenvalue_deviation < eigenvalue_tol
1554 :
1555 6 : DEALLOCATE (eigenvalue_buffer)
1556 6 : CALL cp_fm_release(ref_real)
1557 6 : CALL cp_fm_release(ref_imag)
1558 6 : CALL cp_fm_release(cand_real)
1559 6 : CALL cp_fm_release(cand_imag)
1560 :
1561 12 : END SUBROUTINE validate_wannier90_reused_mos
1562 :
1563 : ! **************************************************************************************************
1564 : !> \brief Compare atom/AO reuse candidates directly to the full-mesh reference MOs.
1565 : !> \param kpoint full Wannier90 export k-point object holding reference MOs
1566 : !> \param qs_kpoint SCF k-point object
1567 : !> \param matrix_s real-space overlap matrix
1568 : !> \param matrix_ks real-space Kohn-Sham matrix
1569 : !> \param cell_to_index real-space cell index table
1570 : !> \param sab_nl overlap neighbor list
1571 : !> \param para_env global parallel environment
1572 : !> \param iw output unit
1573 : !> \param max_subspace_deviation largest best-candidate subspace deviation
1574 : !> \param min_svalue smallest best-candidate singular value
1575 : !> \param max_metric_deviation largest S(k)-metric deviation of a candidate
1576 : !> \param max_residual largest H(k),S(k) eigen-residual of a candidate
1577 : ! **************************************************************************************************
1578 6 : SUBROUTINE diagnose_wannier90_scf_reuse_candidates(kpoint, qs_kpoint, matrix_s, matrix_ks, &
1579 : cell_to_index, sab_nl, para_env, iw, &
1580 : max_subspace_deviation, min_svalue, &
1581 : max_metric_deviation, max_residual)
1582 : TYPE(kpoint_type), POINTER :: kpoint, qs_kpoint
1583 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_ks
1584 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1585 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1586 : POINTER :: sab_nl
1587 : TYPE(mp_para_env_type), POINTER :: para_env
1588 : INTEGER, INTENT(IN) :: iw
1589 : REAL(KIND=dp), INTENT(OUT) :: max_subspace_deviation, min_svalue, &
1590 : max_metric_deviation, max_residual
1591 :
1592 : REAL(KIND=dp), PARAMETER :: print_tol = 1.0e-4_dp, &
1593 : residual_print_tol = 1.0e-3_dp
1594 :
1595 : CHARACTER(LEN=default_string_length) :: reason
1596 : INTEGER :: ik, ikpgr, ikred, ispin, isym_try, nao, &
1597 : nao_src, nkp, nmo, nmo_src, nspins
1598 6 : INTEGER, ALLOCATABLE, DIMENSION(:) :: source_kpoint, sym_index
1599 : INTEGER, DIMENSION(2) :: kp_range, source_kp_range
1600 : LOGICAL :: my_kpgrp, ok, source_window
1601 : REAL(KIND=dp) :: best_deviation, best_metric_deviation, best_residual, best_svalue, &
1602 : candidate_deviation, candidate_metric_deviation, candidate_metric_min, &
1603 : candidate_residual, candidate_svalue, owner_count, ref_metric_deviation, ref_metric_min, &
1604 : ref_residual
1605 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_buffer
1606 6 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1607 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_source, matrix_struct_work
1608 : TYPE(cp_fm_type) :: dst_imag, dst_real, ref_imag, ref_real, &
1609 : src_imag, src_imag_full, src_real, &
1610 : src_real_full
1611 : TYPE(cp_fm_type), POINTER :: fmi, fmr, src_fmi, src_fmr
1612 : TYPE(kpoint_env_type), POINTER :: kp, kp_source
1613 : TYPE(kpoint_sym_type), POINTER :: kpsym
1614 :
1615 6 : max_subspace_deviation = HUGE(1.0_dp)
1616 6 : min_svalue = 0.0_dp
1617 6 : max_metric_deviation = HUGE(1.0_dp)
1618 6 : max_residual = HUGE(1.0_dp)
1619 6 : NULLIFY (matrix_struct_source, matrix_struct_work, fmi, fmr, src_fmi, src_fmr)
1620 :
1621 6 : CALL build_wannier90_scf_mapping(kpoint, qs_kpoint, source_kpoint, sym_index, ok, reason)
1622 6 : IF (.NOT. ok) RETURN
1623 6 : kp => kpoint%kp_env(1)%kpoint_env
1624 6 : kp_source => qs_kpoint%kp_env(1)%kpoint_env
1625 6 : IF (SIZE(kp%mos, 1) < 2 .OR. SIZE(kp_source%mos, 1) < 2) THEN
1626 0 : DEALLOCATE (source_kpoint, sym_index)
1627 0 : RETURN
1628 : END IF
1629 6 : nspins = SIZE(kp%mos, 2)
1630 6 : CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo)
1631 6 : CALL get_mo_set(kp_source%mos(1, 1), nao=nao_src, nmo=nmo_src)
1632 6 : CALL para_env%max(nao_src)
1633 6 : CALL para_env%max(nmo_src)
1634 6 : IF (nao_src /= nao .OR. nmo_src < nmo) THEN
1635 0 : DEALLOCATE (source_kpoint, sym_index)
1636 0 : RETURN
1637 : END IF
1638 6 : source_window = nmo_src > nmo
1639 6 : CALL get_kpoint_info(kpoint, nkp=nkp, kp_range=kp_range)
1640 6 : CALL get_kpoint_info(qs_kpoint, kp_range=source_kp_range)
1641 6 : IF (source_kp_range(1) /= 1 .OR. source_kp_range(2) /= qs_kpoint%nkp) THEN
1642 0 : DEALLOCATE (source_kpoint, sym_index)
1643 0 : RETURN
1644 : END IF
1645 :
1646 : CALL cp_fm_struct_create(matrix_struct_work, nrow_global=nao, ncol_global=nmo, &
1647 6 : para_env=para_env, context=kpoint%blacs_env)
1648 6 : CALL cp_fm_create(ref_real, matrix_struct_work)
1649 6 : CALL cp_fm_create(ref_imag, matrix_struct_work)
1650 6 : CALL cp_fm_create(src_real, matrix_struct_work)
1651 6 : CALL cp_fm_create(src_imag, matrix_struct_work)
1652 6 : CALL cp_fm_create(dst_real, matrix_struct_work)
1653 6 : CALL cp_fm_create(dst_imag, matrix_struct_work)
1654 18 : ALLOCATE (eigenvalues_buffer(nmo))
1655 6 : IF (source_window) THEN
1656 : CALL cp_fm_struct_create(matrix_struct_source, nrow_global=nao, ncol_global=nmo_src, &
1657 0 : para_env=para_env, context=kpoint%blacs_env)
1658 0 : CALL cp_fm_create(src_real_full, matrix_struct_source)
1659 0 : CALL cp_fm_create(src_imag_full, matrix_struct_source)
1660 : END IF
1661 :
1662 6 : max_subspace_deviation = 0.0_dp
1663 6 : min_svalue = HUGE(1.0_dp)
1664 6 : max_metric_deviation = 0.0_dp
1665 6 : max_residual = 0.0_dp
1666 12 : DO ispin = 1, nspins
1667 284 : DO ik = 1, nkp
1668 272 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
1669 : IF (my_kpgrp) THEN
1670 272 : ikpgr = ik - kp_range(1) + 1
1671 272 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1672 272 : fmr => kp%mos(1, ispin)%mo_coeff
1673 272 : fmi => kp%mos(2, ispin)%mo_coeff
1674 : ELSE
1675 : NULLIFY (fmr, fmi)
1676 : END IF
1677 272 : CALL cp_fm_copy_general(fmr, ref_real, para_env)
1678 272 : CALL cp_fm_copy_general(fmi, ref_imag, para_env)
1679 :
1680 272 : ikred = source_kpoint(ik)
1681 272 : my_kpgrp = (ikred >= source_kp_range(1) .AND. ikred <= source_kp_range(2))
1682 : IF (my_kpgrp) THEN
1683 272 : ikpgr = ikred - source_kp_range(1) + 1
1684 272 : kp_source => qs_kpoint%kp_env(ikpgr)%kpoint_env
1685 272 : src_fmr => kp_source%mos(1, ispin)%mo_coeff
1686 272 : src_fmi => kp_source%mos(2, ispin)%mo_coeff
1687 272 : CALL get_mo_set(kp_source%mos(1, ispin), eigenvalues=eigenvalues)
1688 1312 : eigenvalues_buffer(1:nmo) = eigenvalues(1:nmo)
1689 : ELSE
1690 0 : NULLIFY (src_fmr, src_fmi)
1691 0 : eigenvalues_buffer(1:nmo) = 0.0_dp
1692 : END IF
1693 272 : IF (my_kpgrp) THEN
1694 272 : owner_count = 1.0_dp
1695 : ELSE
1696 0 : owner_count = 0.0_dp
1697 : END IF
1698 272 : CALL para_env%sum(owner_count)
1699 272 : CALL para_env%sum(eigenvalues_buffer)
1700 1312 : IF (owner_count > 0.0_dp) eigenvalues_buffer(1:nmo) = eigenvalues_buffer(1:nmo)/owner_count
1701 : CALL measure_wannier90_eigenspace_quality(ref_real, ref_imag, matrix_s, matrix_ks, &
1702 : kpoint%xkp(1:3, ik), cell_to_index, sab_nl, &
1703 : para_env, ispin, eigenvalues_buffer, ok, &
1704 272 : ref_metric_deviation, ref_metric_min, ref_residual)
1705 272 : IF (ok .AND. para_env%is_source() .AND. iw > 0 .AND. &
1706 : (ref_metric_deviation > print_tol .OR. ref_residual > residual_print_tol)) THEN
1707 : WRITE (iw, '(T2,A,I0,A,ES10.3,A,ES10.3,A,ES10.3)') &
1708 0 : "WANNIER90| reference k=", ik, " dM=", ref_metric_deviation, &
1709 0 : " smin=", ref_metric_min, " resid=", ref_residual
1710 : END IF
1711 272 : IF (source_window) THEN
1712 0 : CALL cp_fm_copy_general(src_fmr, src_real_full, para_env)
1713 0 : CALL cp_fm_copy_general(src_fmi, src_imag_full, para_env)
1714 0 : CALL copy_wannier90_mo_window(src_real_full, src_real, nmo)
1715 0 : CALL copy_wannier90_mo_window(src_imag_full, src_imag, nmo)
1716 : ELSE
1717 272 : CALL cp_fm_copy_general(src_fmr, src_real, para_env)
1718 272 : CALL cp_fm_copy_general(src_fmi, src_imag, para_env)
1719 : END IF
1720 :
1721 272 : best_deviation = HUGE(1.0_dp)
1722 272 : best_metric_deviation = 0.0_dp
1723 272 : best_residual = 0.0_dp
1724 272 : best_svalue = 0.0_dp
1725 272 : IF (sym_index(ik) <= 0) THEN
1726 : CALL transform_wannier90_scf_mo(src_real, src_imag, dst_real, dst_imag, qs_kpoint, &
1727 36 : ikred, sym_index(ik), para_env, ok, reason)
1728 36 : IF (ok) THEN
1729 : CALL measure_wannier90_subspace_error(ref_real, ref_imag, dst_real, dst_imag, matrix_s, &
1730 : kpoint%xkp(1:3, ik), cell_to_index, sab_nl, &
1731 36 : para_env, ok, candidate_deviation, candidate_svalue)
1732 36 : IF (ok) THEN
1733 : CALL measure_wannier90_eigenspace_quality(dst_real, dst_imag, matrix_s, matrix_ks, &
1734 : kpoint%xkp(1:3, ik), cell_to_index, sab_nl, &
1735 : para_env, ispin, eigenvalues_buffer, ok, &
1736 : candidate_metric_deviation, &
1737 36 : candidate_metric_min, candidate_residual)
1738 : END IF
1739 36 : IF (ok) THEN
1740 36 : best_deviation = candidate_deviation
1741 36 : best_metric_deviation = candidate_metric_deviation
1742 36 : best_residual = candidate_residual
1743 36 : best_svalue = candidate_svalue
1744 36 : IF (para_env%is_source() .AND. iw > 0 .AND. &
1745 : (candidate_deviation > print_tol .OR. candidate_metric_deviation > print_tol .OR. &
1746 : candidate_residual > residual_print_tol)) THEN
1747 : WRITE (iw, '(T2,A,I0,A,I0,A,I0,A,ES10.3,A,ES10.3,A,ES10.3,A,ES10.3)') &
1748 0 : "WANNIER90| reuse candidate k=", ik, " src=", ikred, " sym=", &
1749 0 : sym_index(ik), " dRef=", candidate_deviation, " dM=", &
1750 0 : candidate_metric_deviation, " smin=", candidate_metric_min, &
1751 0 : " resid=", candidate_residual
1752 : END IF
1753 : END IF
1754 : END IF
1755 236 : ELSEIF (ASSOCIATED(qs_kpoint%kp_sym)) THEN
1756 236 : kpsym => qs_kpoint%kp_sym(ikred)%kpoint_sym
1757 236 : IF (ASSOCIATED(kpsym)) THEN
1758 87404 : DO isym_try = 1, kpsym%nwred
1759 87168 : IF (.NOT. same_periodic_kpoint(kpoint%xkp(1:3, ik), &
1760 : kpsym%xkp(1:3, isym_try))) CYCLE
1761 : CALL transform_wannier90_scf_mo(src_real, src_imag, dst_real, dst_imag, &
1762 5264 : qs_kpoint, ikred, isym_try, para_env, ok, reason)
1763 5264 : IF (.NOT. ok) CYCLE
1764 : CALL measure_wannier90_subspace_error(ref_real, ref_imag, dst_real, dst_imag, &
1765 : matrix_s, kpoint%xkp(1:3, ik), cell_to_index, &
1766 : sab_nl, para_env, ok, candidate_deviation, &
1767 5264 : candidate_svalue)
1768 5264 : IF (ok) THEN
1769 : CALL measure_wannier90_eigenspace_quality(dst_real, dst_imag, matrix_s, matrix_ks, &
1770 : kpoint%xkp(1:3, ik), cell_to_index, &
1771 : sab_nl, para_env, ispin, &
1772 : eigenvalues_buffer, ok, &
1773 : candidate_metric_deviation, &
1774 5264 : candidate_metric_min, candidate_residual)
1775 : END IF
1776 10764 : IF (ok .AND. candidate_deviation < best_deviation) THEN
1777 388 : best_deviation = candidate_deviation
1778 388 : best_metric_deviation = candidate_metric_deviation
1779 388 : best_residual = candidate_residual
1780 388 : best_svalue = candidate_svalue
1781 : END IF
1782 : END DO
1783 : END IF
1784 : END IF
1785 278 : IF (best_deviation < HUGE(1.0_dp)) THEN
1786 272 : max_subspace_deviation = MAX(max_subspace_deviation, best_deviation)
1787 272 : min_svalue = MIN(min_svalue, best_svalue)
1788 272 : max_metric_deviation = MAX(max_metric_deviation, best_metric_deviation)
1789 272 : max_residual = MAX(max_residual, best_residual)
1790 : END IF
1791 : END DO
1792 : END DO
1793 6 : CALL para_env%max(max_subspace_deviation)
1794 6 : CALL para_env%min(min_svalue)
1795 6 : CALL para_env%max(max_metric_deviation)
1796 6 : CALL para_env%max(max_residual)
1797 :
1798 6 : IF (source_window) THEN
1799 0 : CALL cp_fm_release(src_real_full)
1800 0 : CALL cp_fm_release(src_imag_full)
1801 0 : CALL cp_fm_struct_release(matrix_struct_source)
1802 : END IF
1803 6 : CALL cp_fm_release(ref_real)
1804 6 : CALL cp_fm_release(ref_imag)
1805 6 : CALL cp_fm_release(src_real)
1806 6 : CALL cp_fm_release(src_imag)
1807 6 : CALL cp_fm_release(dst_real)
1808 6 : CALL cp_fm_release(dst_imag)
1809 6 : CALL cp_fm_struct_release(matrix_struct_work)
1810 6 : DEALLOCATE (eigenvalues_buffer)
1811 6 : DEALLOCATE (source_kpoint, sym_index)
1812 :
1813 24 : END SUBROUTINE diagnose_wannier90_scf_reuse_candidates
1814 :
1815 : ! **************************************************************************************************
1816 : !> \brief Measure the S(k)-metric distance between two Wannier90 MO subspaces.
1817 : !> \param ref_real real part of reference MO coefficients
1818 : !> \param ref_imag imaginary part of reference MO coefficients
1819 : !> \param cand_real real part of candidate MO coefficients
1820 : !> \param cand_imag imaginary part of candidate MO coefficients
1821 : !> \param matrix_s real-space overlap matrix
1822 : !> \param xkp target k-point coordinate
1823 : !> \param cell_to_index real-space cell index table
1824 : !> \param sab_nl overlap neighbor list
1825 : !> \param para_env global parallel environment
1826 : !> \param success true if the metric comparison was performed
1827 : !> \param max_subspace_deviation largest deviation of singular values from one
1828 : !> \param min_svalue smallest singular value of C_ref^+ S(k) C_candidate
1829 : ! **************************************************************************************************
1830 5572 : SUBROUTINE measure_wannier90_subspace_error(ref_real, ref_imag, cand_real, cand_imag, matrix_s, &
1831 : xkp, cell_to_index, sab_nl, para_env, success, &
1832 : max_subspace_deviation, min_svalue)
1833 : TYPE(cp_fm_type), INTENT(IN) :: ref_real, ref_imag, cand_real, cand_imag
1834 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
1835 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
1836 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1837 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1838 : POINTER :: sab_nl
1839 : TYPE(mp_para_env_type), POINTER :: para_env
1840 : LOGICAL, INTENT(OUT) :: success
1841 : REAL(KIND=dp), INTENT(OUT) :: max_subspace_deviation, min_svalue
1842 :
1843 5572 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: metric_projected, metric_vectors, &
1844 5572 : overlap, ref_coeff, s_cand
1845 : INTEGER :: ib, nao, nmo, nmo_candidate
1846 : REAL(KIND=dp) :: singular_value
1847 5572 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: metric_values
1848 5572 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ref_i, ref_r, s_cand_i, s_cand_r
1849 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_metric
1850 : TYPE(cp_fm_type) :: s_cand_imag, s_cand_real
1851 :
1852 5572 : success = .FALSE.
1853 5572 : max_subspace_deviation = HUGE(1.0_dp)
1854 5572 : min_svalue = 0.0_dp
1855 5572 : NULLIFY (matrix_struct_metric)
1856 :
1857 : CALL cp_fm_get_info(ref_real, nrow_global=nao, ncol_global=nmo, &
1858 5572 : matrix_struct=matrix_struct_metric)
1859 5572 : CALL cp_fm_get_info(cand_real, ncol_global=nmo_candidate)
1860 5572 : IF (nmo_candidate /= nmo) RETURN
1861 :
1862 5572 : CALL cp_fm_create(s_cand_real, matrix_struct_metric)
1863 5572 : CALL cp_fm_create(s_cand_imag, matrix_struct_metric)
1864 : CALL apply_wannier90_kp_matrix(matrix_s, 1, xkp, cell_to_index, sab_nl, &
1865 5572 : cand_real, cand_imag, s_cand_real, s_cand_imag)
1866 :
1867 55720 : ALLOCATE (ref_r(nao, nmo), ref_i(nao, nmo), s_cand_r(nao, nmo), s_cand_i(nao, nmo))
1868 5572 : CALL cp_fm_get_submatrix(ref_real, ref_r)
1869 5572 : CALL cp_fm_get_submatrix(ref_imag, ref_i)
1870 5572 : CALL cp_fm_get_submatrix(s_cand_real, s_cand_r)
1871 5572 : CALL cp_fm_get_submatrix(s_cand_imag, s_cand_i)
1872 :
1873 : ALLOCATE (ref_coeff(nao, nmo), s_cand(nao, nmo), overlap(nmo, nmo), &
1874 83580 : metric_projected(nmo, nmo), metric_vectors(nmo, nmo), metric_values(nmo))
1875 893468 : ref_coeff(:, :) = CMPLX(ref_r, ref_i, KIND=dp)
1876 893468 : s_cand(:, :) = CMPLX(s_cand_r, s_cand_i, KIND=dp)
1877 3576000 : overlap(:, :) = MATMUL(CONJG(TRANSPOSE(ref_coeff)), s_cand)
1878 460336 : metric_projected(:, :) = MATMUL(CONJG(TRANSPOSE(overlap)), overlap)
1879 222548 : metric_projected(:, :) = 0.5_dp*(metric_projected + CONJG(TRANSPOSE(metric_projected)))
1880 5572 : CALL diag_complex(metric_projected, metric_vectors, metric_values)
1881 :
1882 5572 : min_svalue = HUGE(1.0_dp)
1883 5572 : max_subspace_deviation = 0.0_dp
1884 27368 : DO ib = 1, nmo
1885 21796 : singular_value = SQRT(MAX(metric_values(ib), 0.0_dp))
1886 21796 : min_svalue = MIN(min_svalue, singular_value)
1887 27368 : max_subspace_deviation = MAX(max_subspace_deviation, ABS(singular_value - 1.0_dp))
1888 : END DO
1889 5572 : CALL para_env%max(max_subspace_deviation)
1890 5572 : CALL para_env%min(min_svalue)
1891 5572 : success = .TRUE.
1892 :
1893 5572 : DEALLOCATE (ref_coeff, s_cand, overlap, metric_projected, metric_vectors, metric_values)
1894 5572 : DEALLOCATE (ref_r, ref_i, s_cand_r, s_cand_i)
1895 5572 : CALL cp_fm_release(s_cand_real)
1896 5572 : CALL cp_fm_release(s_cand_imag)
1897 :
1898 16716 : END SUBROUTINE measure_wannier90_subspace_error
1899 :
1900 : ! **************************************************************************************************
1901 : !> \brief Measure whether transformed MOs are an H(k),S(k) invariant eigenspace.
1902 : !> \param cand_real real part of candidate MO coefficients
1903 : !> \param cand_imag imaginary part of candidate MO coefficients
1904 : !> \param matrix_s real-space overlap matrix
1905 : !> \param matrix_ks real-space Kohn-Sham matrix
1906 : !> \param xkp target k-point coordinate
1907 : !> \param cell_to_index real-space cell index table
1908 : !> \param sab_nl overlap neighbor list
1909 : !> \param para_env global parallel environment
1910 : !> \param ispin spin index
1911 : !> \param eigenvalues source MO eigenvalues corresponding to the candidate columns
1912 : !> \param success true if the metric and residual checks were performed
1913 : !> \param metric_deviation largest deviation of eigenvalues of C^+ S(k) C from one
1914 : !> \param min_metric_eigenvalue smallest eigenvalue of C^+ S(k) C
1915 : !> \param residual_norm largest element of H(k) C - S(k) C eps
1916 : ! **************************************************************************************************
1917 5572 : SUBROUTINE measure_wannier90_eigenspace_quality(cand_real, cand_imag, matrix_s, matrix_ks, xkp, &
1918 5572 : cell_to_index, sab_nl, para_env, ispin, eigenvalues, &
1919 : success, metric_deviation, min_metric_eigenvalue, &
1920 : residual_norm)
1921 : TYPE(cp_fm_type), INTENT(IN) :: cand_real, cand_imag
1922 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_ks
1923 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
1924 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1925 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1926 : POINTER :: sab_nl
1927 : TYPE(mp_para_env_type), POINTER :: para_env
1928 : INTEGER, INTENT(IN) :: ispin
1929 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: eigenvalues
1930 : LOGICAL, INTENT(OUT) :: success
1931 : REAL(KIND=dp), INTENT(OUT) :: metric_deviation, min_metric_eigenvalue, &
1932 : residual_norm
1933 :
1934 5572 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cand_coeff, h_coeff, metric_vectors, &
1935 5572 : residual_block, s_coeff, s_projected
1936 : INTEGER :: ib, nao, nmo
1937 5572 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: metric_values
1938 5572 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cand_i, cand_r, h_coeff_i, h_coeff_r, &
1939 5572 : s_coeff_i, s_coeff_r
1940 : TYPE(cp_cfm_type) :: cand_cfm, metric_cfm, s_cfm
1941 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_metric, &
1942 : matrix_struct_projected
1943 : TYPE(cp_fm_type) :: h_cand_imag, h_cand_real, s_cand_imag, &
1944 : s_cand_real, tmp_fm
1945 :
1946 5572 : success = .FALSE.
1947 5572 : metric_deviation = HUGE(1.0_dp)
1948 5572 : min_metric_eigenvalue = 0.0_dp
1949 5572 : residual_norm = HUGE(1.0_dp)
1950 5572 : NULLIFY (matrix_struct_metric, matrix_struct_projected)
1951 :
1952 : CALL cp_fm_get_info(cand_real, nrow_global=nao, ncol_global=nmo, &
1953 5572 : matrix_struct=matrix_struct_metric)
1954 5572 : IF (SIZE(eigenvalues) < nmo) RETURN
1955 :
1956 5572 : CALL cp_fm_create(s_cand_real, matrix_struct_metric)
1957 5572 : CALL cp_fm_create(s_cand_imag, matrix_struct_metric)
1958 5572 : CALL cp_fm_create(h_cand_real, matrix_struct_metric)
1959 5572 : CALL cp_fm_create(h_cand_imag, matrix_struct_metric)
1960 5572 : CALL cp_fm_create(tmp_fm, matrix_struct_metric)
1961 5572 : CALL cp_cfm_create(cand_cfm, matrix_struct_metric)
1962 5572 : CALL cp_cfm_create(s_cfm, matrix_struct_metric)
1963 :
1964 : CALL apply_wannier90_kp_matrix(matrix_s, 1, xkp, cell_to_index, sab_nl, &
1965 5572 : cand_real, cand_imag, s_cand_real, s_cand_imag)
1966 : CALL apply_wannier90_kp_matrix(matrix_ks, ispin, xkp, cell_to_index, sab_nl, &
1967 5572 : cand_real, cand_imag, h_cand_real, h_cand_imag)
1968 :
1969 : ALLOCATE (cand_r(nao, nmo), cand_i(nao, nmo), s_coeff_r(nao, nmo), &
1970 78008 : s_coeff_i(nao, nmo), h_coeff_r(nao, nmo), h_coeff_i(nao, nmo))
1971 5572 : CALL cp_fm_get_submatrix(cand_real, cand_r)
1972 5572 : CALL cp_fm_get_submatrix(cand_imag, cand_i)
1973 5572 : CALL cp_fm_get_submatrix(s_cand_real, s_coeff_r)
1974 5572 : CALL cp_fm_get_submatrix(s_cand_imag, s_coeff_i)
1975 5572 : CALL cp_fm_get_submatrix(h_cand_real, h_coeff_r)
1976 5572 : CALL cp_fm_get_submatrix(h_cand_imag, h_coeff_i)
1977 :
1978 : ALLOCATE (cand_coeff(nao, nmo), h_coeff(nao, nmo), metric_vectors(nmo, nmo), &
1979 : residual_block(nao, nmo), s_coeff(nao, nmo), s_projected(nmo, nmo), &
1980 94724 : metric_values(nmo))
1981 893468 : cand_coeff(:, :) = CMPLX(cand_r, cand_i, KIND=dp)
1982 893468 : s_coeff(:, :) = CMPLX(s_coeff_r, s_coeff_i, KIND=dp)
1983 893468 : h_coeff(:, :) = CMPLX(h_coeff_r, h_coeff_i, KIND=dp)
1984 : CALL cp_fm_struct_create(matrix_struct_projected, nrow_global=nmo, ncol_global=nmo, &
1985 : para_env=matrix_struct_metric%para_env, &
1986 5572 : context=matrix_struct_metric%context)
1987 5572 : CALL cp_cfm_create(metric_cfm, matrix_struct_projected)
1988 5572 : CALL cp_fm_to_cfm(cand_real, cand_imag, cand_cfm)
1989 5572 : CALL cp_fm_to_cfm(s_cand_real, s_cand_imag, s_cfm)
1990 : CALL cp_cfm_gemm("C", "N", nmo, nmo, nao, CMPLX(1.0_dp, 0.0_dp, KIND=dp), cand_cfm, &
1991 5572 : s_cfm, CMPLX(0.0_dp, 0.0_dp, KIND=dp), metric_cfm)
1992 5572 : CALL cp_cfm_get_submatrix(metric_cfm, s_projected)
1993 222548 : s_projected(:, :) = 0.5_dp*(s_projected + CONJG(TRANSPOSE(s_projected)))
1994 5572 : CALL diag_complex(s_projected, metric_vectors, metric_values)
1995 27368 : metric_deviation = MAXVAL(ABS(metric_values - 1.0_dp))
1996 27368 : min_metric_eigenvalue = MINVAL(metric_values)
1997 :
1998 893468 : residual_block(:, :) = h_coeff
1999 27368 : DO ib = 1, nmo
2000 893468 : residual_block(:, ib) = residual_block(:, ib) - eigenvalues(ib)*s_coeff(:, ib)
2001 : END DO
2002 893468 : residual_norm = MAXVAL(ABS(residual_block))
2003 5572 : CALL para_env%max(metric_deviation)
2004 5572 : CALL para_env%min(min_metric_eigenvalue)
2005 5572 : CALL para_env%max(residual_norm)
2006 5572 : success = .TRUE.
2007 :
2008 0 : DEALLOCATE (cand_coeff, h_coeff, metric_vectors, residual_block, s_coeff, s_projected, &
2009 5572 : metric_values)
2010 5572 : DEALLOCATE (cand_r, cand_i, s_coeff_r, s_coeff_i, h_coeff_r, h_coeff_i)
2011 5572 : CALL cp_fm_release(s_cand_real)
2012 5572 : CALL cp_fm_release(s_cand_imag)
2013 5572 : CALL cp_fm_release(h_cand_real)
2014 5572 : CALL cp_fm_release(h_cand_imag)
2015 5572 : CALL cp_fm_release(tmp_fm)
2016 5572 : CALL cp_cfm_release(cand_cfm)
2017 5572 : CALL cp_cfm_release(s_cfm)
2018 5572 : CALL cp_cfm_release(metric_cfm)
2019 5572 : CALL cp_fm_struct_release(matrix_struct_projected)
2020 :
2021 22288 : END SUBROUTINE measure_wannier90_eigenspace_quality
2022 :
2023 : ! **************************************************************************************************
2024 : !> \brief Copy the leading MO columns from a larger SCF MO matrix into the Wannier90 export window.
2025 : !> \param source source MO coefficient matrix
2026 : !> \param destination destination MO coefficient matrix
2027 : !> \param ncol number of columns to copy
2028 : ! **************************************************************************************************
2029 0 : SUBROUTINE copy_wannier90_mo_window(source, destination, ncol)
2030 : TYPE(cp_fm_type), INTENT(IN) :: source, destination
2031 : INTEGER, INTENT(IN) :: ncol
2032 :
2033 : INTEGER :: ncol_source, nrow
2034 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: destination_buffer, source_buffer
2035 :
2036 0 : CALL cp_fm_get_info(source, nrow_global=nrow, ncol_global=ncol_source)
2037 0 : CPASSERT(ncol_source >= ncol)
2038 0 : ALLOCATE (source_buffer(nrow, ncol_source), destination_buffer(nrow, ncol))
2039 0 : CALL cp_fm_get_submatrix(source, source_buffer)
2040 0 : destination_buffer(1:nrow, 1:ncol) = source_buffer(1:nrow, 1:ncol)
2041 0 : CALL cp_fm_set_submatrix(destination, destination_buffer)
2042 0 : DEALLOCATE (source_buffer, destination_buffer)
2043 :
2044 0 : END SUBROUTINE copy_wannier90_mo_window
2045 :
2046 : ! **************************************************************************************************
2047 : !> \brief Apply a complex k-point matrix to a complex MO coefficient matrix.
2048 : !> \param rsmat real-space matrix images
2049 : !> \param ispin spin index for rsmat
2050 : !> \param xkp target k-point coordinate
2051 : !> \param cell_to_index real-space cell index table
2052 : !> \param sab_nl overlap neighbor list
2053 : !> \param coeff_real real part of input MO coefficients
2054 : !> \param coeff_imag imaginary part of input MO coefficients
2055 : !> \param result_real real part of matrix-vector product
2056 : !> \param result_imag imaginary part of matrix-vector product
2057 : ! **************************************************************************************************
2058 105720 : SUBROUTINE apply_wannier90_kp_matrix(rsmat, ispin, xkp, cell_to_index, sab_nl, &
2059 : coeff_real, coeff_imag, result_real, result_imag)
2060 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
2061 : INTEGER, INTENT(IN) :: ispin
2062 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
2063 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2064 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2065 : POINTER :: sab_nl
2066 : TYPE(cp_fm_type), INTENT(IN) :: coeff_real, coeff_imag, result_real, &
2067 : result_imag
2068 :
2069 : INTEGER :: nao, ncol
2070 : TYPE(cp_cfm_type) :: coeff_cfm, kmat_cfm, result_cfm
2071 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_ao, matrix_struct_coeff
2072 : TYPE(cp_fm_type) :: mat_imag, mat_real
2073 : TYPE(dbcsr_type), POINTER :: kmat_imag, kmat_imag_full, kmat_real, &
2074 : kmat_real_full
2075 :
2076 17620 : NULLIFY (matrix_struct_ao, matrix_struct_coeff, kmat_imag, kmat_imag_full, kmat_real, &
2077 17620 : kmat_real_full)
2078 :
2079 : CALL cp_fm_get_info(coeff_real, nrow_global=nao, ncol_global=ncol, &
2080 17620 : matrix_struct=matrix_struct_coeff)
2081 :
2082 17620 : ALLOCATE (kmat_real, kmat_imag, kmat_real_full, kmat_imag_full)
2083 : CALL dbcsr_create(kmat_real, template=rsmat(ispin, 1)%matrix, &
2084 17620 : matrix_type=dbcsr_type_symmetric)
2085 : CALL dbcsr_create(kmat_imag, template=rsmat(ispin, 1)%matrix, &
2086 17620 : matrix_type=dbcsr_type_antisymmetric)
2087 : CALL dbcsr_create(kmat_real_full, template=rsmat(ispin, 1)%matrix, &
2088 17620 : matrix_type=dbcsr_type_no_symmetry)
2089 : CALL dbcsr_create(kmat_imag_full, template=rsmat(ispin, 1)%matrix, &
2090 17620 : matrix_type=dbcsr_type_no_symmetry)
2091 17620 : CALL cp_dbcsr_alloc_block_from_nbl(kmat_real, sab_nl)
2092 17620 : CALL cp_dbcsr_alloc_block_from_nbl(kmat_imag, sab_nl)
2093 17620 : CALL dbcsr_set(kmat_real, 0.0_dp)
2094 17620 : CALL dbcsr_set(kmat_imag, 0.0_dp)
2095 : CALL rskp_transform(kmat_real, kmat_imag, rsmat=rsmat, ispin=ispin, &
2096 17620 : xkp=xkp, cell_to_index=cell_to_index, sab_nl=sab_nl)
2097 17620 : CALL dbcsr_desymmetrize(kmat_real, kmat_real_full)
2098 17620 : CALL dbcsr_desymmetrize(kmat_imag, kmat_imag_full)
2099 :
2100 : CALL cp_fm_struct_create(matrix_struct_ao, nrow_global=nao, ncol_global=nao, &
2101 : para_env=matrix_struct_coeff%para_env, &
2102 17620 : context=matrix_struct_coeff%context)
2103 17620 : CALL cp_fm_create(mat_real, matrix_struct_ao)
2104 17620 : CALL cp_fm_create(mat_imag, matrix_struct_ao)
2105 17620 : CALL copy_dbcsr_to_fm(kmat_real_full, mat_real)
2106 17620 : CALL copy_dbcsr_to_fm(kmat_imag_full, mat_imag)
2107 :
2108 17620 : CALL cp_cfm_create(kmat_cfm, matrix_struct_ao)
2109 17620 : CALL cp_cfm_create(coeff_cfm, matrix_struct_coeff)
2110 17620 : CALL cp_cfm_create(result_cfm, matrix_struct_coeff)
2111 17620 : CALL cp_fm_to_cfm(mat_real, mat_imag, kmat_cfm)
2112 17620 : CALL cp_fm_to_cfm(coeff_real, coeff_imag, coeff_cfm)
2113 : CALL cp_cfm_gemm("N", "N", nao, ncol, nao, CMPLX(1.0_dp, 0.0_dp, KIND=dp), kmat_cfm, &
2114 17620 : coeff_cfm, CMPLX(0.0_dp, 0.0_dp, KIND=dp), result_cfm)
2115 17620 : CALL cp_cfm_to_fm(result_cfm, result_real, result_imag)
2116 :
2117 17620 : CALL cp_fm_release(mat_real)
2118 17620 : CALL cp_fm_release(mat_imag)
2119 17620 : CALL cp_cfm_release(kmat_cfm)
2120 17620 : CALL cp_cfm_release(coeff_cfm)
2121 17620 : CALL cp_cfm_release(result_cfm)
2122 17620 : CALL cp_fm_struct_release(matrix_struct_ao)
2123 17620 : CALL dbcsr_deallocate_matrix(kmat_real)
2124 17620 : CALL dbcsr_deallocate_matrix(kmat_imag)
2125 17620 : CALL dbcsr_deallocate_matrix(kmat_real_full)
2126 17620 : CALL dbcsr_deallocate_matrix(kmat_imag_full)
2127 :
2128 17620 : END SUBROUTINE apply_wannier90_kp_matrix
2129 :
2130 : ! **************************************************************************************************
2131 : !> \brief Rayleigh-Ritz stabilize a symmetry-reconstructed Wannier90 MO subspace.
2132 : !> \param dst_real real part of transformed MO coefficients
2133 : !> \param dst_imag imaginary part of transformed MO coefficients
2134 : !> \param matrix_s real-space overlap matrix
2135 : !> \param matrix_ks real-space Kohn-Sham matrix
2136 : !> \param xkp target k-point coordinate
2137 : !> \param cell_to_index real-space cell index table
2138 : !> \param sab_nl overlap neighbor list
2139 : !> \param ispin spin index
2140 : !> \param eigenvalues Ritz eigenvalues of the stabilized subspace
2141 : !> \param degenerate_band_tol degeneracy threshold
2142 : !> \param success true if the subspace was stabilized
2143 : !> \param reason diagnostic message
2144 : !> \param aligned_blocks number of stabilized subspaces
2145 : !> \param aligned_max_size largest stabilized subspace
2146 : !> \param aligned_min_svalue smallest S(k)-metric eigenvalue
2147 : !> \param max_residual largest Ritz residual
2148 : ! **************************************************************************************************
2149 452 : SUBROUTINE ritz_stabilize_wannier90_subspace(dst_real, dst_imag, matrix_s, matrix_ks, &
2150 904 : xkp, cell_to_index, sab_nl, ispin, eigenvalues, &
2151 : degenerate_band_tol, success, reason, aligned_blocks, &
2152 : aligned_max_size, aligned_min_svalue, max_residual)
2153 : TYPE(cp_fm_type), INTENT(IN) :: dst_real, dst_imag
2154 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_ks
2155 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
2156 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2157 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2158 : POINTER :: sab_nl
2159 : INTEGER, INTENT(IN) :: ispin
2160 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: eigenvalues
2161 : REAL(KIND=dp), INTENT(IN) :: degenerate_band_tol
2162 : LOGICAL, INTENT(OUT) :: success
2163 : CHARACTER(LEN=*), INTENT(OUT) :: reason
2164 : INTEGER, INTENT(OUT) :: aligned_blocks, aligned_max_size
2165 : REAL(KIND=dp), INTENT(OUT) :: aligned_min_svalue, max_residual
2166 :
2167 : REAL(KIND=dp), PARAMETER :: residual_tol = 1.0e-2_dp
2168 :
2169 452 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: block_coeff, h_block, h_coeff, &
2170 452 : h_projected, h_projected_work, metric_vectors, residual_block, ritz_vectors, s_block, &
2171 452 : s_coeff, s_projected, stabilized
2172 : INTEGER :: block_first, block_last, block_size, ib, &
2173 : nao, nmo
2174 : REAL(KIND=dp) :: metric_deviation, norm_value, &
2175 : residual_norm
2176 452 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: metric_values, ritz_values
2177 452 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dst_i, dst_r, h_coeff_i, h_coeff_r, &
2178 452 : s_coeff_i, s_coeff_r
2179 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_metric
2180 : TYPE(cp_fm_type) :: h_dst_imag, h_dst_real, s_dst_imag, &
2181 : s_dst_real, tmp_fm
2182 :
2183 452 : success = .FALSE.
2184 452 : reason = ""
2185 452 : aligned_blocks = 0
2186 452 : aligned_max_size = 0
2187 452 : aligned_min_svalue = HUGE(1.0_dp)
2188 452 : max_residual = 0.0_dp
2189 :
2190 452 : NULLIFY (matrix_struct_metric)
2191 : CALL cp_fm_get_info(dst_real, nrow_global=nao, ncol_global=nmo, &
2192 452 : matrix_struct=matrix_struct_metric)
2193 452 : IF (SIZE(eigenvalues) < nmo) THEN
2194 0 : reason = "not enough eigenvalues for Wannier90 Ritz subspace stabilization"
2195 : RETURN
2196 : END IF
2197 :
2198 452 : CALL cp_fm_create(s_dst_real, matrix_struct_metric)
2199 452 : CALL cp_fm_create(s_dst_imag, matrix_struct_metric)
2200 452 : CALL cp_fm_create(h_dst_real, matrix_struct_metric)
2201 452 : CALL cp_fm_create(h_dst_imag, matrix_struct_metric)
2202 452 : CALL cp_fm_create(tmp_fm, matrix_struct_metric)
2203 :
2204 : CALL apply_wannier90_kp_matrix(matrix_s, 1, xkp, cell_to_index, sab_nl, &
2205 452 : dst_real, dst_imag, s_dst_real, s_dst_imag)
2206 : CALL apply_wannier90_kp_matrix(matrix_ks, ispin, xkp, cell_to_index, sab_nl, &
2207 452 : dst_real, dst_imag, h_dst_real, h_dst_imag)
2208 :
2209 0 : ALLOCATE (dst_r(nao, nmo), dst_i(nao, nmo), s_coeff_r(nao, nmo), s_coeff_i(nao, nmo), &
2210 6328 : h_coeff_r(nao, nmo), h_coeff_i(nao, nmo))
2211 452 : CALL cp_fm_get_submatrix(dst_real, dst_r)
2212 452 : CALL cp_fm_get_submatrix(dst_imag, dst_i)
2213 452 : CALL cp_fm_get_submatrix(s_dst_real, s_coeff_r)
2214 452 : CALL cp_fm_get_submatrix(s_dst_imag, s_coeff_i)
2215 452 : CALL cp_fm_get_submatrix(h_dst_real, h_coeff_r)
2216 452 : CALL cp_fm_get_submatrix(h_dst_imag, h_coeff_i)
2217 :
2218 2712 : ALLOCATE (s_coeff(nao, nmo), h_coeff(nao, nmo))
2219 86132 : s_coeff(:, :) = CMPLX(s_coeff_r, s_coeff_i, KIND=dp)
2220 86132 : h_coeff(:, :) = CMPLX(h_coeff_r, h_coeff_i, KIND=dp)
2221 :
2222 452 : block_first = 1
2223 1588 : DO WHILE (block_first <= nmo)
2224 : block_last = block_first
2225 1792 : DO WHILE (block_last < nmo)
2226 1340 : IF (ABS(eigenvalues(block_last + 1) - eigenvalues(block_last)) >= degenerate_band_tol) EXIT
2227 1136 : block_last = block_last + 1
2228 : END DO
2229 1136 : block_size = block_last - block_first + 1
2230 1136 : IF (block_size > 1) THEN
2231 : ! The atom/AO operation fixes the subspace, while the little-group gauge inside an
2232 : ! exactly degenerate manifold is arbitrary. Stabilize only that manifold and verify
2233 : ! that it is an invariant H(k),S(k) subspace before exporting it to Wannier90.
2234 0 : ALLOCATE (block_coeff(nao, block_size), h_block(nao, block_size), &
2235 0 : h_projected(block_size, block_size), h_projected_work(block_size, block_size), &
2236 0 : metric_vectors(block_size, block_size), residual_block(nao, block_size), &
2237 0 : ritz_vectors(block_size, block_size), s_block(nao, block_size), &
2238 0 : s_projected(block_size, block_size), stabilized(nao, block_size), &
2239 11232 : metric_values(block_size), ritz_values(block_size))
2240 : block_coeff(:, :) = CMPLX(dst_r(:, block_first:block_last), &
2241 59376 : dst_i(:, block_first:block_last), KIND=dp)
2242 59376 : s_block(:, :) = s_coeff(:, block_first:block_last)
2243 59376 : h_block(:, :) = h_coeff(:, block_first:block_last)
2244 933648 : s_projected(:, :) = MATMUL(CONJG(TRANSPOSE(block_coeff)), s_block)
2245 933648 : h_projected(:, :) = MATMUL(CONJG(TRANSPOSE(block_coeff)), h_block)
2246 8304 : s_projected(:, :) = 0.5_dp*(s_projected + CONJG(TRANSPOSE(s_projected)))
2247 8304 : h_projected(:, :) = 0.5_dp*(h_projected + CONJG(TRANSPOSE(h_projected)))
2248 :
2249 432 : CALL diag_complex(s_projected, metric_vectors, metric_values)
2250 1520 : aligned_min_svalue = MIN(aligned_min_svalue, MINVAL(metric_values))
2251 1520 : metric_deviation = MAXVAL(ABS(metric_values - 1.0_dp))
2252 1520 : IF (MINVAL(metric_values) < 1.0e-10_dp) THEN
2253 : WRITE (reason, "(A,I0,A,ES9.2,A,ES9.2)") &
2254 0 : "singular metric blk=", block_first, " smin=", MINVAL(metric_values), &
2255 0 : " dS=", metric_deviation
2256 0 : max_residual = HUGE(1.0_dp)
2257 0 : DEALLOCATE (block_coeff, h_block, h_projected, h_projected_work, metric_vectors, &
2258 0 : residual_block, ritz_vectors, s_block, s_projected, stabilized, &
2259 0 : metric_values, ritz_values)
2260 0 : DEALLOCATE (s_coeff, h_coeff, dst_r, dst_i, s_coeff_r, s_coeff_i, h_coeff_r, &
2261 0 : h_coeff_i)
2262 0 : CALL cp_fm_release(s_dst_real)
2263 0 : CALL cp_fm_release(s_dst_imag)
2264 0 : CALL cp_fm_release(h_dst_real)
2265 0 : CALL cp_fm_release(h_dst_imag)
2266 0 : CALL cp_fm_release(tmp_fm)
2267 0 : RETURN
2268 : END IF
2269 :
2270 1520 : DO ib = 1, block_size
2271 4368 : metric_vectors(:, ib) = metric_vectors(:, ib)/SQRT(metric_values(ib))
2272 : END DO
2273 50640 : h_projected_work(:, :) = MATMUL(h_projected, metric_vectors)
2274 50640 : h_projected(:, :) = MATMUL(CONJG(TRANSPOSE(metric_vectors)), h_projected_work)
2275 8304 : h_projected(:, :) = 0.5_dp*(h_projected + CONJG(TRANSPOSE(h_projected)))
2276 432 : CALL diag_complex(h_projected, ritz_vectors, ritz_values)
2277 50640 : h_projected_work(:, :) = MATMUL(metric_vectors, ritz_vectors)
2278 4368 : ritz_vectors(:, :) = h_projected_work
2279 623888 : stabilized(:, :) = MATMUL(block_coeff, ritz_vectors)
2280 623888 : residual_block(:, :) = MATMUL(h_block, ritz_vectors)
2281 59376 : h_block(:, :) = residual_block
2282 623888 : residual_block(:, :) = MATMUL(s_block, ritz_vectors)
2283 59376 : s_block(:, :) = residual_block
2284 1520 : DO ib = 1, block_size
2285 58944 : norm_value = SQRT(ABS(REAL(DOT_PRODUCT(stabilized(:, ib), s_block(:, ib)), KIND=dp)))
2286 1520 : IF (norm_value > EPSILON(1.0_dp)) THEN
2287 58944 : stabilized(:, ib) = stabilized(:, ib)/norm_value
2288 58944 : h_block(:, ib) = h_block(:, ib)/norm_value
2289 58944 : s_block(:, ib) = s_block(:, ib)/norm_value
2290 : END IF
2291 : END DO
2292 59376 : residual_block(:, :) = h_block
2293 1520 : DO ib = 1, block_size
2294 : residual_block(:, ib) = residual_block(:, ib) - &
2295 59376 : eigenvalues(block_first + ib - 1)*s_block(:, ib)
2296 : END DO
2297 59376 : residual_norm = MAXVAL(ABS(residual_block))
2298 432 : max_residual = MAX(max_residual, residual_norm)
2299 432 : IF (residual_norm > residual_tol) THEN
2300 : WRITE (reason, "(A,I0,A,ES9.2)") &
2301 0 : "blk=", block_first, " dS=", metric_deviation
2302 0 : DEALLOCATE (block_coeff, h_block, h_projected, h_projected_work, metric_vectors, &
2303 0 : residual_block, ritz_vectors, s_block, s_projected, stabilized, &
2304 0 : metric_values, ritz_values)
2305 0 : DEALLOCATE (s_coeff, h_coeff, dst_r, dst_i, s_coeff_r, s_coeff_i, h_coeff_r, &
2306 0 : h_coeff_i)
2307 0 : CALL cp_fm_release(s_dst_real)
2308 0 : CALL cp_fm_release(s_dst_imag)
2309 0 : CALL cp_fm_release(h_dst_real)
2310 0 : CALL cp_fm_release(h_dst_imag)
2311 0 : CALL cp_fm_release(tmp_fm)
2312 0 : RETURN
2313 : END IF
2314 :
2315 59376 : dst_r(:, block_first:block_last) = REAL(stabilized, KIND=dp)
2316 59376 : dst_i(:, block_first:block_last) = AIMAG(stabilized)
2317 59376 : h_coeff(:, block_first:block_last) = h_block
2318 59376 : s_coeff(:, block_first:block_last) = s_block
2319 432 : aligned_blocks = aligned_blocks + 1
2320 432 : aligned_max_size = MAX(aligned_max_size, block_size)
2321 0 : DEALLOCATE (block_coeff, h_block, h_projected, h_projected_work, metric_vectors, &
2322 0 : residual_block, ritz_vectors, s_block, s_projected, stabilized, metric_values, &
2323 432 : ritz_values)
2324 : END IF
2325 1136 : block_first = block_last + 1
2326 : END DO
2327 :
2328 2244 : DO ib = 1, nmo
2329 85680 : residual_norm = MAXVAL(ABS(h_coeff(:, ib) - eigenvalues(ib)*s_coeff(:, ib)))
2330 2244 : max_residual = MAX(max_residual, residual_norm)
2331 : END DO
2332 452 : IF (max_residual > residual_tol) THEN
2333 : WRITE (reason, "(A,ES10.3)") &
2334 0 : "atom/AO W90 reuse guarded: Ritz residual=", max_residual
2335 0 : DEALLOCATE (s_coeff, h_coeff, dst_r, dst_i, s_coeff_r, s_coeff_i, h_coeff_r, h_coeff_i)
2336 0 : CALL cp_fm_release(s_dst_real)
2337 0 : CALL cp_fm_release(s_dst_imag)
2338 0 : CALL cp_fm_release(h_dst_real)
2339 0 : CALL cp_fm_release(h_dst_imag)
2340 0 : CALL cp_fm_release(tmp_fm)
2341 0 : RETURN
2342 : END IF
2343 :
2344 452 : CALL cp_fm_set_submatrix(dst_real, dst_r)
2345 452 : CALL cp_fm_set_submatrix(dst_imag, dst_i)
2346 :
2347 452 : IF (aligned_blocks == 0) aligned_min_svalue = 0.0_dp
2348 452 : success = .TRUE.
2349 :
2350 452 : DEALLOCATE (s_coeff, h_coeff, dst_r, dst_i, s_coeff_r, s_coeff_i, h_coeff_r, h_coeff_i)
2351 452 : CALL cp_fm_release(s_dst_real)
2352 452 : CALL cp_fm_release(s_dst_imag)
2353 452 : CALL cp_fm_release(h_dst_real)
2354 452 : CALL cp_fm_release(h_dst_imag)
2355 452 : CALL cp_fm_release(tmp_fm)
2356 :
2357 1808 : END SUBROUTINE ritz_stabilize_wannier90_subspace
2358 :
2359 : ! **************************************************************************************************
2360 : !> \brief Reconstruct a Wannier90 export window from a larger symmetry-transformed SCF MO space.
2361 : !> \param src_real real part of the transformed source MO window
2362 : !> \param src_imag imaginary part of the transformed source MO window
2363 : !> \param dst_real real part of the exported reconstructed MO coefficients
2364 : !> \param dst_imag imaginary part of the exported reconstructed MO coefficients
2365 : !> \param matrix_s real-space overlap matrix
2366 : !> \param matrix_ks real-space Kohn-Sham matrix
2367 : !> \param xkp target k-point coordinate
2368 : !> \param cell_to_index real-space cell index table
2369 : !> \param sab_nl overlap neighbor list
2370 : !> \param ispin spin index
2371 : !> \param eigenvalues reconstructed target eigenvalues for the exported window
2372 : !> \param nmo_export number of MOs to export
2373 : !> \param success true if the reconstructed window is an invariant H(k),S(k) subspace
2374 : !> \param reason diagnostic message
2375 : !> \param min_svalue smallest S(k)-metric eigenvalue in the source window
2376 : !> \param max_residual largest target Ritz residual
2377 : ! **************************************************************************************************
2378 0 : SUBROUTINE ritz_reconstruct_wannier90_window(src_real, src_imag, dst_real, dst_imag, matrix_s, &
2379 : matrix_ks, xkp, cell_to_index, sab_nl, ispin, &
2380 0 : eigenvalues, nmo_export, success, reason, min_svalue, &
2381 : max_residual)
2382 : TYPE(cp_fm_type), INTENT(IN) :: src_real, src_imag, dst_real, dst_imag
2383 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_ks
2384 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
2385 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2386 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2387 : POINTER :: sab_nl
2388 : INTEGER, INTENT(IN) :: ispin
2389 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: eigenvalues
2390 : INTEGER, INTENT(IN) :: nmo_export
2391 : LOGICAL, INTENT(OUT) :: success
2392 : CHARACTER(LEN=*), INTENT(OUT) :: reason
2393 : REAL(KIND=dp), INTENT(OUT) :: min_svalue, max_residual
2394 :
2395 : REAL(KIND=dp), PARAMETER :: eigenvalue_tol = 1.0e-6_dp, &
2396 : residual_tol = 1.0e-7_dp
2397 :
2398 0 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coeff_work, h_coeff, h_projected, &
2399 0 : h_projected_work, metric_vectors, residual_block, ritz_vectors, s_coeff, s_projected, &
2400 0 : source_coeff, stabilized
2401 : INTEGER :: ib, nao, nmo_source
2402 : REAL(KIND=dp) :: max_eigenvalue_shift, metric_deviation, &
2403 : norm_value
2404 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: metric_values, ritz_values, &
2405 0 : source_eigenvalues
2406 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dst_i, dst_r, h_coeff_i, h_coeff_r, &
2407 0 : s_coeff_i, s_coeff_r, src_i, src_r
2408 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_metric
2409 : TYPE(cp_fm_type) :: h_src_imag, h_src_real, s_src_imag, &
2410 : s_src_real, tmp_fm
2411 :
2412 0 : success = .FALSE.
2413 0 : reason = ""
2414 0 : min_svalue = HUGE(1.0_dp)
2415 0 : max_residual = HUGE(1.0_dp)
2416 :
2417 0 : NULLIFY (matrix_struct_metric)
2418 : CALL cp_fm_get_info(src_real, nrow_global=nao, ncol_global=nmo_source, &
2419 0 : matrix_struct=matrix_struct_metric)
2420 0 : IF (nmo_export > nmo_source) THEN
2421 0 : reason = "Wannier90 export window is larger than the transformed SCF MO space"
2422 0 : RETURN
2423 : END IF
2424 0 : IF (SIZE(eigenvalues) < nmo_export) THEN
2425 0 : reason = "not enough eigenvalue storage for Wannier90 source-window reconstruction"
2426 : RETURN
2427 : END IF
2428 :
2429 0 : CALL cp_fm_create(s_src_real, matrix_struct_metric)
2430 0 : CALL cp_fm_create(s_src_imag, matrix_struct_metric)
2431 0 : CALL cp_fm_create(h_src_real, matrix_struct_metric)
2432 0 : CALL cp_fm_create(h_src_imag, matrix_struct_metric)
2433 0 : CALL cp_fm_create(tmp_fm, matrix_struct_metric)
2434 :
2435 : CALL apply_wannier90_kp_matrix(matrix_s, 1, xkp, cell_to_index, sab_nl, &
2436 0 : src_real, src_imag, s_src_real, s_src_imag)
2437 : CALL apply_wannier90_kp_matrix(matrix_ks, ispin, xkp, cell_to_index, sab_nl, &
2438 0 : src_real, src_imag, h_src_real, h_src_imag)
2439 :
2440 0 : ALLOCATE (src_r(nao, nmo_source), src_i(nao, nmo_source), &
2441 0 : s_coeff_r(nao, nmo_source), s_coeff_i(nao, nmo_source), &
2442 0 : h_coeff_r(nao, nmo_source), h_coeff_i(nao, nmo_source), &
2443 0 : dst_r(nao, nmo_export), dst_i(nao, nmo_export))
2444 0 : CALL cp_fm_get_submatrix(src_real, src_r)
2445 0 : CALL cp_fm_get_submatrix(src_imag, src_i)
2446 0 : CALL cp_fm_get_submatrix(s_src_real, s_coeff_r)
2447 0 : CALL cp_fm_get_submatrix(s_src_imag, s_coeff_i)
2448 0 : CALL cp_fm_get_submatrix(h_src_real, h_coeff_r)
2449 0 : CALL cp_fm_get_submatrix(h_src_imag, h_coeff_i)
2450 :
2451 0 : ALLOCATE (source_coeff(nao, nmo_source), s_coeff(nao, nmo_source), &
2452 0 : h_coeff(nao, nmo_source), h_projected(nmo_source, nmo_source), &
2453 0 : h_projected_work(nmo_source, nmo_source), metric_vectors(nmo_source, nmo_source), &
2454 0 : residual_block(nao, nmo_export), ritz_vectors(nmo_source, nmo_source), &
2455 0 : s_projected(nmo_source, nmo_source), stabilized(nao, nmo_export), &
2456 0 : coeff_work(nao, nmo_source), metric_values(nmo_source), ritz_values(nmo_source), &
2457 0 : source_eigenvalues(nmo_export))
2458 0 : source_eigenvalues(1:nmo_export) = eigenvalues(1:nmo_export)
2459 0 : source_coeff(:, :) = CMPLX(src_r, src_i, KIND=dp)
2460 0 : s_coeff(:, :) = CMPLX(s_coeff_r, s_coeff_i, KIND=dp)
2461 0 : h_coeff(:, :) = CMPLX(h_coeff_r, h_coeff_i, KIND=dp)
2462 0 : s_projected(:, :) = MATMUL(CONJG(TRANSPOSE(source_coeff)), s_coeff)
2463 0 : h_projected(:, :) = MATMUL(CONJG(TRANSPOSE(source_coeff)), h_coeff)
2464 0 : s_projected(:, :) = 0.5_dp*(s_projected + CONJG(TRANSPOSE(s_projected)))
2465 0 : h_projected(:, :) = 0.5_dp*(h_projected + CONJG(TRANSPOSE(h_projected)))
2466 :
2467 : reconstruct_window: BLOCK
2468 0 : CALL diag_complex(s_projected, metric_vectors, metric_values)
2469 0 : min_svalue = MINVAL(metric_values)
2470 0 : metric_deviation = MAXVAL(ABS(metric_values - 1.0_dp))
2471 0 : IF (min_svalue < 1.0e-10_dp) THEN
2472 : WRITE (reason, "(A,ES9.2,A,ES9.2)") &
2473 0 : "singular expanded metric smin=", min_svalue, " dS=", metric_deviation
2474 0 : EXIT reconstruct_window
2475 : END IF
2476 :
2477 0 : DO ib = 1, nmo_source
2478 0 : metric_vectors(:, ib) = metric_vectors(:, ib)/SQRT(metric_values(ib))
2479 : END DO
2480 0 : h_projected_work(:, :) = MATMUL(h_projected, metric_vectors)
2481 0 : h_projected(:, :) = MATMUL(CONJG(TRANSPOSE(metric_vectors)), h_projected_work)
2482 0 : h_projected(:, :) = 0.5_dp*(h_projected + CONJG(TRANSPOSE(h_projected)))
2483 0 : CALL diag_complex(h_projected, ritz_vectors, ritz_values)
2484 0 : h_projected_work(:, :) = MATMUL(metric_vectors, ritz_vectors)
2485 0 : ritz_vectors(:, :) = h_projected_work
2486 0 : stabilized(:, :) = MATMUL(source_coeff, ritz_vectors(:, 1:nmo_export))
2487 0 : coeff_work(:, :) = MATMUL(h_coeff, ritz_vectors)
2488 0 : h_coeff(:, :) = coeff_work
2489 0 : coeff_work(:, :) = MATMUL(s_coeff, ritz_vectors)
2490 0 : s_coeff(:, :) = coeff_work
2491 0 : DO ib = 1, nmo_export
2492 0 : norm_value = SQRT(ABS(REAL(DOT_PRODUCT(stabilized(:, ib), s_coeff(:, ib)), KIND=dp)))
2493 0 : IF (norm_value > EPSILON(1.0_dp)) THEN
2494 0 : stabilized(:, ib) = stabilized(:, ib)/norm_value
2495 0 : h_coeff(:, ib) = h_coeff(:, ib)/norm_value
2496 0 : s_coeff(:, ib) = s_coeff(:, ib)/norm_value
2497 : END IF
2498 : END DO
2499 0 : residual_block(:, :) = h_coeff(:, 1:nmo_export)
2500 0 : DO ib = 1, nmo_export
2501 0 : residual_block(:, ib) = residual_block(:, ib) - ritz_values(ib)*s_coeff(:, ib)
2502 : END DO
2503 0 : max_residual = MAXVAL(ABS(residual_block))
2504 0 : IF (max_residual > residual_tol) THEN
2505 : WRITE (reason, "(A,ES9.2)") &
2506 0 : "expanded dS=", metric_deviation
2507 0 : EXIT reconstruct_window
2508 : END IF
2509 0 : max_eigenvalue_shift = MAXVAL(ABS(ritz_values(1:nmo_export) - source_eigenvalues(1:nmo_export)))
2510 0 : IF (max_eigenvalue_shift > eigenvalue_tol) THEN
2511 : WRITE (reason, "(A,ES9.2)") &
2512 0 : "expanded dS=", metric_deviation
2513 0 : EXIT reconstruct_window
2514 : END IF
2515 :
2516 0 : dst_r(:, :) = REAL(stabilized, KIND=dp)
2517 0 : dst_i(:, :) = AIMAG(stabilized)
2518 0 : CALL cp_fm_set_submatrix(dst_real, dst_r)
2519 0 : CALL cp_fm_set_submatrix(dst_imag, dst_i)
2520 0 : success = .TRUE.
2521 :
2522 : END BLOCK reconstruct_window
2523 :
2524 0 : DEALLOCATE (source_coeff, s_coeff, h_coeff, h_projected, h_projected_work, metric_vectors, &
2525 0 : residual_block, ritz_vectors, s_projected, stabilized, coeff_work, metric_values, &
2526 0 : ritz_values, source_eigenvalues)
2527 0 : DEALLOCATE (src_r, src_i, s_coeff_r, s_coeff_i, h_coeff_r, h_coeff_i, dst_r, dst_i)
2528 0 : CALL cp_fm_release(s_src_real)
2529 0 : CALL cp_fm_release(s_src_imag)
2530 0 : CALL cp_fm_release(h_src_real)
2531 0 : CALL cp_fm_release(h_src_imag)
2532 0 : CALL cp_fm_release(tmp_fm)
2533 :
2534 0 : END SUBROUTINE ritz_reconstruct_wannier90_window
2535 :
2536 : ! **************************************************************************************************
2537 : !> \brief Map the full Wannier90 mesh to SCF representative k-points and symmetry operations.
2538 : !> \param kpoint full Wannier90 export k-point object
2539 : !> \param qs_kpoint SCF k-point object
2540 : !> \param source_kpoint source representative index for each full k-point
2541 : !> \param sym_index symmetry entry in source kp_sym; 0 direct, -1 time reversal only
2542 : !> \param success true if every full k-point was mapped
2543 : !> \param reason diagnostic message
2544 : ! **************************************************************************************************
2545 42 : SUBROUTINE build_wannier90_scf_mapping(kpoint, qs_kpoint, source_kpoint, sym_index, success, &
2546 : reason)
2547 : TYPE(kpoint_type), POINTER :: kpoint, qs_kpoint
2548 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: source_kpoint, sym_index
2549 : LOGICAL, INTENT(OUT) :: success
2550 : CHARACTER(LEN=*), INTENT(OUT) :: reason
2551 :
2552 : INTEGER :: ik, ikred, imatch, isym, nfull
2553 : TYPE(kpoint_sym_type), POINTER :: kpsym
2554 :
2555 42 : success = .FALSE.
2556 42 : reason = ""
2557 42 : nfull = kpoint%nkp
2558 168 : ALLOCATE (source_kpoint(nfull), sym_index(nfull))
2559 42 : source_kpoint(:) = 0
2560 42 : sym_index(:) = 0
2561 :
2562 142 : DO ikred = 1, qs_kpoint%nkp
2563 100 : imatch = find_matching_kpoint(kpoint%xkp, qs_kpoint%xkp(1:3, ikred))
2564 142 : IF (imatch > 0 .AND. source_kpoint(imatch) == 0) THEN
2565 100 : source_kpoint(imatch) = ikred
2566 100 : sym_index(imatch) = 0
2567 : END IF
2568 : END DO
2569 :
2570 : ! Prefer pure time-reversal partners before general atom/AO symmetry operations.
2571 142 : DO ikred = 1, qs_kpoint%nkp
2572 400 : imatch = find_matching_kpoint(kpoint%xkp, -qs_kpoint%xkp(1:3, ikred))
2573 142 : IF (imatch > 0 .AND. source_kpoint(imatch) == 0) THEN
2574 68 : source_kpoint(imatch) = ikred
2575 68 : sym_index(imatch) = -1
2576 : END IF
2577 : END DO
2578 :
2579 42 : IF (ASSOCIATED(qs_kpoint%kp_sym)) THEN
2580 142 : DO ikred = 1, qs_kpoint%nkp
2581 100 : kpsym => qs_kpoint%kp_sym(ikred)%kpoint_sym
2582 100 : IF (.NOT. ASSOCIATED(kpsym)) CYCLE
2583 100 : IF (.NOT. kpsym%apply_symmetry) CYCLE
2584 18758 : DO isym = 1, kpsym%nwred
2585 18656 : imatch = find_matching_kpoint(kpoint%xkp, kpsym%xkp(1:3, isym))
2586 18756 : IF (imatch > 0 .AND. source_kpoint(imatch) == 0) THEN
2587 604 : source_kpoint(imatch) = ikred
2588 604 : sym_index(imatch) = isym
2589 : END IF
2590 : END DO
2591 : END DO
2592 : END IF
2593 :
2594 814 : DO ik = 1, nfull
2595 814 : IF (source_kpoint(ik) == 0) THEN
2596 0 : reason = "not all full-mesh k-points are represented by the SCF symmetry orbits"
2597 0 : RETURN
2598 : END IF
2599 : END DO
2600 42 : success = .TRUE.
2601 :
2602 42 : END SUBROUTINE build_wannier90_scf_mapping
2603 :
2604 : ! **************************************************************************************************
2605 : !> \brief Find a fractional k-point in a periodic mesh.
2606 : !> \param xkp_mesh mesh coordinates
2607 : !> \param xkp_search coordinate to find
2608 : !> \return matching index, or zero when no match is found
2609 : ! **************************************************************************************************
2610 18856 : INTEGER FUNCTION find_matching_kpoint(xkp_mesh, xkp_search) RESULT(ik_match)
2611 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp_mesh
2612 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp_search
2613 :
2614 : INTEGER :: ik
2615 :
2616 18856 : ik_match = 0
2617 430600 : DO ik = 1, SIZE(xkp_mesh, 2)
2618 430600 : IF (same_periodic_kpoint(xkp_mesh(1:3, ik), xkp_search)) THEN
2619 18856 : ik_match = ik
2620 18856 : RETURN
2621 : END IF
2622 : END DO
2623 :
2624 : END FUNCTION find_matching_kpoint
2625 :
2626 : ! **************************************************************************************************
2627 : !> \brief Compare two fractional k-points modulo reciprocal lattice vectors.
2628 : !> \param xkp_a first k-point
2629 : !> \param xkp_b second k-point
2630 : !> \return true if the k-points are equivalent
2631 : ! **************************************************************************************************
2632 520832 : LOGICAL FUNCTION same_periodic_kpoint(xkp_a, xkp_b) RESULT(same)
2633 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp_a, xkp_b
2634 :
2635 : REAL(KIND=dp), DIMENSION(3) :: diff
2636 :
2637 2083328 : diff(1:3) = xkp_a(1:3) - xkp_b(1:3)
2638 2083328 : diff(1:3) = diff(1:3) - REAL(NINT(diff(1:3)), KIND=dp)
2639 2083328 : same = SUM(ABS(diff(1:3))) < 1.0e-8_dp
2640 :
2641 520832 : END FUNCTION same_periodic_kpoint
2642 :
2643 : ! **************************************************************************************************
2644 : !> \brief Transform one SCF MO coefficient matrix to an equivalent full-mesh k-point.
2645 : !> \param src_real real part of source MO coefficients
2646 : !> \param src_imag imaginary part of source MO coefficients
2647 : !> \param dst_real real part of transformed MO coefficients
2648 : !> \param dst_imag imaginary part of transformed MO coefficients
2649 : !> \param qs_kpoint SCF k-point object containing symmetry operations
2650 : !> \param ikred representative k-point index
2651 : !> \param isym symmetry operation index; 0 direct copy, -1 time reversal
2652 : !> \param para_env global parallel environment
2653 : !> \param success true if the transformation was performed
2654 : !> \param reason diagnostic message
2655 : ! **************************************************************************************************
2656 5752 : SUBROUTINE transform_wannier90_scf_mo(src_real, src_imag, dst_real, dst_imag, qs_kpoint, ikred, &
2657 : isym, para_env, success, reason)
2658 : TYPE(cp_fm_type), INTENT(IN) :: src_real, src_imag, dst_real, dst_imag
2659 : TYPE(kpoint_type), POINTER :: qs_kpoint
2660 : INTEGER, INTENT(IN) :: ikred, isym
2661 : TYPE(mp_para_env_type), POINTER :: para_env
2662 : LOGICAL, INTENT(OUT) :: success
2663 : CHARACTER(LEN=*), INTENT(OUT) :: reason
2664 :
2665 : INTEGER :: iao, iatom, ikind, imo, irow, irow_source, irow_target, irow_work, nao, natom, &
2666 : nmo, rot_slot, source_atom, source_dim, target_atom, target_dim
2667 5752 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ao_size, ao_start
2668 : LOGICAL :: reverse_phase, time_reversal
2669 : REAL(KIND=dp) :: arg, coeff_imag, coeff_real, coskl, &
2670 : phase_imag, phase_real, sinkl
2671 5752 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dst_i, dst_r, src_i, src_r
2672 : REAL(KIND=dp), DIMENSION(3) :: xkp_phase
2673 5752 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rotmat
2674 5752 : TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kind_rot
2675 : TYPE(kpoint_sym_type), POINTER :: kpsym
2676 :
2677 5752 : success = .FALSE.
2678 5752 : reason = ""
2679 5752 : IF (isym == 0) THEN
2680 60 : CALL cp_fm_copy_general(src_real, dst_real, para_env)
2681 60 : CALL cp_fm_copy_general(src_imag, dst_imag, para_env)
2682 60 : success = .TRUE.
2683 60 : RETURN
2684 : END IF
2685 :
2686 5692 : CALL cp_fm_get_info(src_real, nrow_global=nao, ncol_global=nmo)
2687 56920 : ALLOCATE (src_r(nao, nmo), src_i(nao, nmo), dst_r(nao, nmo), dst_i(nao, nmo))
2688 5692 : CALL cp_fm_get_submatrix(src_real, src_r)
2689 5692 : CALL cp_fm_get_submatrix(src_imag, src_i)
2690 5692 : dst_r(:, :) = 0.0_dp
2691 5692 : dst_i(:, :) = 0.0_dp
2692 :
2693 5692 : IF (isym == -1) THEN
2694 10812 : dst_r(1:nao, 1:nmo) = src_r(1:nao, 1:nmo)
2695 10812 : dst_i(1:nao, 1:nmo) = -src_i(1:nao, 1:nmo)
2696 60 : CALL cp_fm_set_submatrix(dst_real, dst_r)
2697 60 : CALL cp_fm_set_submatrix(dst_imag, dst_i)
2698 60 : DEALLOCATE (src_r, src_i, dst_r, dst_i)
2699 60 : success = .TRUE.
2700 60 : RETURN
2701 : END IF
2702 :
2703 5632 : IF (.NOT. ASSOCIATED(qs_kpoint%kp_sym)) THEN
2704 0 : reason = "SCF symmetry operation data are not available"
2705 0 : DEALLOCATE (src_r, src_i, dst_r, dst_i)
2706 0 : RETURN
2707 : END IF
2708 5632 : kpsym => qs_kpoint%kp_sym(ikred)%kpoint_sym
2709 5632 : IF (.NOT. ASSOCIATED(kpsym)) THEN
2710 0 : reason = "SCF k-point symmetry operation is not available"
2711 0 : DEALLOCATE (src_r, src_i, dst_r, dst_i)
2712 0 : RETURN
2713 : END IF
2714 5632 : IF (isym > kpsym%nwred) THEN
2715 0 : reason = "SCF k-point symmetry operation index is out of range"
2716 0 : DEALLOCATE (src_r, src_i, dst_r, dst_i)
2717 0 : RETURN
2718 : END IF
2719 5632 : IF (.NOT. ASSOCIATED(qs_kpoint%atype) .OR. .NOT. ASSOCIATED(qs_kpoint%kind_rotmat)) THEN
2720 0 : reason = "SCF atom mappings or basis rotation matrices are not available"
2721 0 : DEALLOCATE (src_r, src_i, dst_r, dst_i)
2722 0 : RETURN
2723 : END IF
2724 :
2725 5632 : rot_slot = find_wannier90_rotation_slot(qs_kpoint, kpsym%rotp(isym))
2726 5632 : IF (rot_slot == 0) THEN
2727 0 : reason = "could not match the SCF symmetry operation to a basis rotation"
2728 0 : DEALLOCATE (src_r, src_i, dst_r, dst_i)
2729 0 : RETURN
2730 : END IF
2731 5632 : kind_rot => qs_kpoint%kind_rotmat(rot_slot, :)
2732 5632 : natom = SIZE(qs_kpoint%atype)
2733 22528 : ALLOCATE (ao_start(natom), ao_size(natom))
2734 49284 : irow = 1
2735 49284 : DO iatom = 1, natom
2736 43652 : ikind = qs_kpoint%atype(iatom)
2737 43652 : IF (.NOT. ASSOCIATED(kind_rot(ikind)%rmat)) THEN
2738 0 : reason = "a required basis rotation matrix is not available"
2739 0 : DEALLOCATE (src_r, src_i, dst_r, dst_i, ao_start, ao_size)
2740 0 : RETURN
2741 : END IF
2742 43652 : ao_start(iatom) = irow
2743 43652 : ao_size(iatom) = SIZE(kind_rot(ikind)%rmat, 2)
2744 49284 : irow = irow + ao_size(iatom)
2745 : END DO
2746 5632 : IF (irow - 1 /= nao) THEN
2747 0 : reason = "atom-resolved AO dimensions do not match the MO coefficient matrix"
2748 0 : DEALLOCATE (src_r, src_i, dst_r, dst_i, ao_start, ao_size)
2749 0 : RETURN
2750 : END IF
2751 :
2752 5632 : time_reversal = kpsym%rotp(isym) < 0
2753 5632 : reverse_phase = qs_kpoint%gamma_centered .AND. ANY(MOD(qs_kpoint%nkp_grid, 2) == 0)
2754 22528 : xkp_phase(1:3) = kpsym%xkp(1:3, isym)
2755 49284 : DO iatom = 1, natom
2756 43652 : source_atom = iatom
2757 43652 : target_atom = kpsym%f0(iatom, isym)
2758 43652 : ikind = qs_kpoint%atype(source_atom)
2759 43652 : rotmat => kind_rot(ikind)%rmat
2760 43652 : source_dim = ao_size(source_atom)
2761 43652 : target_dim = ao_size(target_atom)
2762 43652 : IF (SIZE(rotmat, 1) /= target_dim .OR. SIZE(rotmat, 2) /= source_dim) THEN
2763 0 : reason = "basis rotation dimensions do not match the atom/AO symmetry transform"
2764 0 : DEALLOCATE (src_r, src_i, dst_r, dst_i, ao_start, ao_size)
2765 0 : RETURN
2766 : END IF
2767 : arg = REAL(kpsym%fcell(1, source_atom, isym), KIND=dp)*xkp_phase(1) + &
2768 : REAL(kpsym%fcell(2, source_atom, isym), KIND=dp)*xkp_phase(2) + &
2769 43652 : REAL(kpsym%fcell(3, source_atom, isym), KIND=dp)*xkp_phase(3)
2770 43652 : IF (ASSOCIATED(kpsym%kgphase)) THEN
2771 43652 : arg = arg + kpsym%kgphase(source_atom, isym)
2772 : END IF
2773 43652 : IF (reverse_phase) arg = -arg
2774 43652 : coskl = COS(twopi*arg)
2775 43652 : sinkl = SIN(twopi*arg)
2776 224408 : DO imo = 1, nmo
2777 1106684 : DO irow_work = 1, source_dim
2778 887908 : irow_source = ao_start(source_atom) + irow_work - 1
2779 887908 : coeff_real = src_r(irow_source, imo)
2780 887908 : coeff_imag = src_i(irow_source, imo)
2781 887908 : IF (time_reversal) coeff_imag = -coeff_imag
2782 887908 : phase_real = coskl*coeff_real - sinkl*coeff_imag
2783 887908 : phase_imag = sinkl*coeff_real + coskl*coeff_imag
2784 5662316 : DO iao = 1, target_dim
2785 4599284 : irow_target = ao_start(target_atom) + iao - 1
2786 : dst_r(irow_target, imo) = dst_r(irow_target, imo) + &
2787 4599284 : rotmat(iao, irow_work)*phase_real
2788 : dst_i(irow_target, imo) = dst_i(irow_target, imo) + &
2789 5487192 : rotmat(iao, irow_work)*phase_imag
2790 : END DO
2791 : END DO
2792 : END DO
2793 : END DO
2794 :
2795 5632 : CALL cp_fm_set_submatrix(dst_real, dst_r)
2796 5632 : CALL cp_fm_set_submatrix(dst_imag, dst_i)
2797 5632 : DEALLOCATE (src_r, src_i, dst_r, dst_i, ao_start, ao_size)
2798 5632 : success = .TRUE.
2799 :
2800 17196 : END SUBROUTINE transform_wannier90_scf_mo
2801 :
2802 : ! **************************************************************************************************
2803 : !> \brief Locate the basis-rotation slot corresponding to a signed k-point symmetry operation.
2804 : !> \param qs_kpoint SCF k-point object
2805 : !> \param rotp signed operation identifier
2806 : !> \return rotation slot, or zero when no slot matches
2807 : ! **************************************************************************************************
2808 5632 : INTEGER FUNCTION find_wannier90_rotation_slot(qs_kpoint, rotp) RESULT(rot_slot)
2809 : TYPE(kpoint_type), POINTER :: qs_kpoint
2810 : INTEGER, INTENT(IN) :: rotp
2811 :
2812 : INTEGER :: irot, rot_abs
2813 :
2814 5632 : rot_slot = 0
2815 5632 : rot_abs = ABS(rotp)
2816 5632 : IF (.NOT. ASSOCIATED(qs_kpoint%ibrot)) RETURN
2817 508008 : DO irot = 1, SIZE(qs_kpoint%ibrot)
2818 508008 : IF (rot_abs == qs_kpoint%ibrot(irot)) THEN
2819 5632 : rot_slot = irot
2820 : RETURN
2821 : END IF
2822 : END DO
2823 :
2824 : END FUNCTION find_wannier90_rotation_slot
2825 :
2826 : ! **************************************************************************************************
2827 : !> \brief Infer a tensor-product Wannier90 mesh from explicit fractional k-point coordinates.
2828 : !> \param kpt_latt explicit k-point coordinates in reciprocal-lattice units
2829 : !> \param mp_grid inferred mesh dimensions
2830 : !> \param valid true if the coordinate set is compatible with a tensor-product mesh
2831 : ! **************************************************************************************************
2832 6 : SUBROUTINE infer_wannier_mp_grid(kpt_latt, mp_grid, valid)
2833 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: kpt_latt
2834 : INTEGER, DIMENSION(3), INTENT(OUT) :: mp_grid
2835 : LOGICAL, INTENT(OUT) :: valid
2836 :
2837 : INTEGER :: coord_id, i, idim, idx, n_unique, &
2838 : num_kpts, stride, unique_id
2839 : LOGICAL :: known
2840 6 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: seen
2841 : REAL(KIND=dp) :: coord
2842 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: unique_coord
2843 :
2844 6 : num_kpts = SIZE(kpt_latt, 2)
2845 6 : mp_grid(:) = 0
2846 18 : ALLOCATE (unique_coord(3, num_kpts))
2847 24 : DO idim = 1, 3
2848 : n_unique = 0
2849 162 : DO i = 1, num_kpts
2850 144 : coord = kpt_latt(idim, i) - FLOOR(kpt_latt(idim, i))
2851 144 : IF (ABS(coord - 1.0_dp) < 1.0e-8_dp) coord = 0.0_dp
2852 144 : known = .FALSE.
2853 216 : DO unique_id = 1, n_unique
2854 216 : IF (ABS(unique_coord(idim, unique_id) - coord) < 1.0e-8_dp) THEN
2855 : known = .TRUE.
2856 : EXIT
2857 : END IF
2858 : END DO
2859 162 : IF (.NOT. known) THEN
2860 36 : n_unique = n_unique + 1
2861 36 : unique_coord(idim, n_unique) = coord
2862 : END IF
2863 : END DO
2864 24 : mp_grid(idim) = n_unique
2865 : END DO
2866 6 : valid = (mp_grid(1)*mp_grid(2)*mp_grid(3) == num_kpts)
2867 6 : IF (valid) THEN
2868 18 : ALLOCATE (seen(num_kpts))
2869 6 : seen(:) = .FALSE.
2870 54 : DO i = 1, num_kpts
2871 : idx = 1
2872 : stride = 1
2873 192 : DO idim = 1, 3
2874 144 : coord = kpt_latt(idim, i) - FLOOR(kpt_latt(idim, i))
2875 144 : IF (ABS(coord - 1.0_dp) < 1.0e-8_dp) coord = 0.0_dp
2876 144 : coord_id = 0
2877 216 : DO unique_id = 1, mp_grid(idim)
2878 216 : IF (ABS(unique_coord(idim, unique_id) - coord) < 1.0e-8_dp) THEN
2879 : coord_id = unique_id
2880 : EXIT
2881 : END IF
2882 : END DO
2883 144 : CPASSERT(coord_id > 0)
2884 144 : idx = idx + (coord_id - 1)*stride
2885 192 : stride = stride*mp_grid(idim)
2886 : END DO
2887 48 : IF (seen(idx)) valid = .FALSE.
2888 54 : seen(idx) = .TRUE.
2889 : END DO
2890 54 : valid = valid .AND. ALL(seen)
2891 6 : DEALLOCATE (seen)
2892 : END IF
2893 6 : DEALLOCATE (unique_coord)
2894 :
2895 6 : END SUBROUTINE infer_wannier_mp_grid
2896 :
2897 0 : END MODULE qs_wannier90
|
