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