(git:fdbe441)
Loading...
Searching...
No Matches
qs_wannier90.F
Go to the documentation of this file.
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! **************************************************************************************************
14MODULE 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
105CONTAINS
106
107! **************************************************************************************************
108!> \brief ...
109!> \param input ...
110!> \param logger ...
111!> \param qs_env ...
112! **************************************************************************************************
113 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 CALL timeset(routinen, handle)
127 w_input => section_vals_get_subs_vals(section_vals=input, &
128 subsection_name="DFT%PRINT%WANNIER90")
129 CALL section_vals_get(w_input, explicit=explicit)
130 IF (explicit) THEN
131
132 iw = cp_logger_get_default_io_unit(logger)
133
134 IF (iw > 0) THEN
135 WRITE (iw, '(/,T2,A)') &
136 '!-----------------------------------------------------------------------------!'
137 WRITE (iw, '(T32,A)') "Interface to Wannier90"
138 WRITE (iw, '(T2,A)') &
139 '!-----------------------------------------------------------------------------!'
140 END IF
141
142 CALL wannier90_files(qs_env, w_input, iw)
143
144 IF (iw > 0) THEN
145 WRITE (iw, '(/,T2,A)') &
146 '!--------------------------------End of Wannier90-----------------------------!'
147 END IF
148 END IF
149 CALL timestop(handle)
150
151 END SUBROUTINE wannier90_interface
152
153! **************************************************************************************************
154!> \brief ...
155!> \param qs_env ...
156!> \param input ...
157!> \param iw ...
158! **************************************************************************************************
159 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: exclude_bands
176 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nblist, nnlist
177 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: nncell
178 INTEGER, DIMENSION(2) :: kp_range
179 INTEGER, DIMENSION(3) :: input_nkp_grid, mp_grid
180 INTEGER, DIMENSION(:), POINTER :: invals
181 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigval
192 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: atoms_cart, b_latt, kpt_latt
193 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: reference_eigenvalues
194 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 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues, wkp, wkp_source
198 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp, xkp_source
199 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 TYPE(cp_fm_type), DIMENSION(2) :: fmk1, fmk2
208 TYPE(cp_fm_type), POINTER :: fmdummy, fmi, fmr
209 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 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 POINTER :: sab_nl
219 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 CALL section_vals_val_get(input, "SEED_NAME", c_val=seed_name)
230 CALL section_vals_val_get(input, "MP_GRID", i_vals=invals, explicit=mp_grid_explicit)
231 CALL section_vals_val_get(input, "KPOINTS_SOURCE", i_val=kpoints_source)
232 CALL section_vals_val_get(input, "WANNIER_FUNCTIONS", i_val=num_wann)
233 CALL section_vals_val_get(input, "ADDED_MOS", i_val=nadd)
234 CALL section_vals_val_get(input, "REUSE_SCF_MOS", l_val=reuse_scf_mos)
235 CALL section_vals_val_get(input, "VALIDATE_REUSE_SCF_MOS", l_val=validate_reuse_scf_mos)
236 CALL section_vals_val_get(input, "USE_BLOCH_PHASES", l_val=use_bloch_phases)
237 reuse_scf_mos = reuse_scf_mos .AND. kpoints_source == w90_kpoints_scf
238 validate_reuse_scf_mos = validate_reuse_scf_mos .AND. reuse_scf_mos
239 mp_grid(1:3) = invals(1:3)
240 ! excluded bands
241 CALL section_vals_val_get(input, "EXCLUDE_BANDS", n_rep_val=n_rep)
242 nexcl = 0
243 DO i_rep = 1, n_rep
244 CALL section_vals_val_get(input, "EXCLUDE_BANDS", i_rep_val=i_rep, i_vals=invals)
245 nexcl = nexcl + SIZE(invals)
246 END DO
247 IF (nexcl > 0) THEN
248 ALLOCATE (exclude_bands(nexcl))
249 nexcl = 0
250 DO i_rep = 1, n_rep
251 CALL section_vals_val_get(input, "EXCLUDE_BANDS", i_rep_val=i_rep, i_vals=invals)
252 exclude_bands(nexcl + 1:nexcl + SIZE(invals)) = invals(:)
253 nexcl = nexcl + SIZE(invals)
254 END DO
255 END IF
256 !
257 ! lattice -> Angstrom
258 CALL get_qs_env(qs_env, cell=cell)
259 CALL get_cell(cell, h=real_lattice, h_inv=h_inv)
260 ! k-points
261 CALL get_qs_env(qs_env, particle_set=particle_set)
262 CALL get_qs_env(qs_env, para_env=para_env)
263 phase_center = 0.0_dp
264 DO i = 1, SIZE(particle_set)
265 phase_center(1:3) = phase_center(1:3) + matmul(h_inv, particle_set(i)%r)
266 END DO
267 phase_center(1:3) = phase_center(1:3)/real(SIZE(particle_set), kind=dp)
268 phase_center(1:3) = phase_center(1:3) - floor(phase_center(1:3))
269 recip_lattice(1:3, 1:3) = h_inv(1:3, 1:3)
270 real_lattice(1:3, 1:3) = angstrom*real_lattice(1:3, 1:3)
271 recip_lattice(1:3, 1:3) = (twopi/angstrom)*transpose(recip_lattice(1:3, 1:3))
272 NULLIFY (kpoint, qs_kpoint, xkp, wkp, xkp_source, wkp_source)
273 CALL get_qs_env(qs_env, do_kpoints=do_kpoints, kpoints=qs_kpoint)
274 input_kpoint_symmetry = .false.
275 input_full_grid = .false.
276 input_kp_scheme = ""
277 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 nkp=nkp, xkp=xkp, wkp=wkp)
282 END IF
283 CALL kpoint_create(kpoint)
284
285 SELECT CASE (kpoints_source)
286 CASE (w90_kpoints_mp_grid)
287 num_kpts = mp_grid(1)*mp_grid(2)*mp_grid(3)
288 ALLOCATE (kpt_latt(3, num_kpts))
289 kpoint%kp_scheme = "MONKHORST-PACK"
290 kpoint%symmetry = .false.
291 kpoint%nkp_grid(1:3) = mp_grid(1:3)
292 kpoint%verbose = .false.
293 kpoint%full_grid = .true.
294 kpoint%eps_geo = 1.0e-6_dp
295 kpoint%use_real_wfn = .false.
296 kpoint%parallel_group_size = para_env%num_pe
297 i = 0
298 DO ix = 0, mp_grid(1) - 1
299 DO iy = 0, mp_grid(2) - 1
300 DO iz = 0, mp_grid(3) - 1
301 i = i + 1
302 kpt_latt(1, i) = real(ix, kind=dp)/real(mp_grid(1), kind=dp)
303 kpt_latt(2, i) = real(iy, kind=dp)/real(mp_grid(2), kind=dp)
304 kpt_latt(3, i) = real(iz, kind=dp)/real(mp_grid(3), kind=dp)
305 END DO
306 END DO
307 END DO
308 kpoint%nkp = num_kpts
309 ALLOCATE (kpoint%xkp(3, num_kpts), kpoint%wkp(num_kpts))
310 kpoint%wkp(:) = 1._dp/real(num_kpts, kind=dp)
311 DO i = 1, num_kpts
312 kpoint%xkp(1:3, i) = (angstrom/twopi)*matmul(recip_lattice, kpt_latt(:, i))
313 END DO
314
315 CASE (w90_kpoints_scf)
316 IF (.NOT. do_kpoints .OR. .NOT. ASSOCIATED(qs_kpoint)) THEN
317 cpabort("WANNIER90%KPOINTS_SOURCE SCF requires an active DFT%KPOINTS section.")
318 END IF
319 SELECT CASE (trim(input_kp_scheme))
320 CASE ("GAMMA")
321 mp_grid(:) = 1
322 num_kpts = 1
323 ALLOCATE (kpt_latt(3, num_kpts))
324 kpt_latt(1:3, 1) = 0.0_dp
325 kpoint%kp_scheme = "GAMMA"
326 kpoint%symmetry = .false.
327 kpoint%verbose = .false.
328 kpoint%full_grid = .true.
329 kpoint%eps_geo = 1.0e-6_dp
330 kpoint%use_real_wfn = .false.
331 kpoint%parallel_group_size = para_env%num_pe
332 kpoint%nkp = num_kpts
333 ALLOCATE (kpoint%xkp(3, num_kpts), kpoint%wkp(num_kpts))
334 kpoint%xkp(1:3, 1) = 0.0_dp
335 kpoint%wkp(1) = 1.0_dp
336
337 CASE ("MONKHORST-PACK", "MACDONALD")
338 mp_grid(1:3) = input_nkp_grid(1:3)
339 kpoint%kp_scheme = input_kp_scheme
340 kpoint%symmetry = .false.
341 kpoint%nkp_grid(1:3) = input_nkp_grid(1:3)
342 kpoint%kp_shift(1:3) = input_kp_shift(1:3)
343 kpoint%gamma_centered = input_gamma_centered
344 kpoint%verbose = .false.
345 kpoint%full_grid = .true.
346 kpoint%eps_geo = 1.0e-6_dp
347 kpoint%use_real_wfn = .false.
348 kpoint%parallel_group_size = para_env%num_pe
349 CALL kpoint_initialize(kpoint, particle_set, cell)
350 num_kpts = kpoint%nkp
351 ALLOCATE (kpt_latt(3, num_kpts))
352 kpt_latt(1:3, 1:num_kpts) = kpoint%xkp(1:3, 1:num_kpts)
353 IF (input_kpoint_symmetry .AND. .NOT. input_full_grid .AND. iw > 0) THEN
354 WRITE (iw, '(T2,A)') &
355 "WANNIER90| SCF k-points are symmetry-reduced; regenerating the full SCF mesh."
356 IF (reuse_scf_mos) THEN
357 WRITE (iw, '(T2,A)') &
358 "WANNIER90| CP2K will try to reconstruct the full-mesh MOs from the SCF orbitals."
359 ELSE
360 WRITE (iw, '(T2,A)') &
361 "WANNIER90| The full exported mesh is diagonalized for the Wannier90 files."
362 END IF
363 END IF
364
365 CASE ("GENERAL")
366 IF (ASSOCIATED(qs_kpoint%xkp_input)) THEN
367 xkp_source => qs_kpoint%xkp_input
368 wkp_source => qs_kpoint%wkp_input
369 ELSE
370 xkp_source => xkp
371 wkp_source => wkp
372 END IF
373 IF (.NOT. ASSOCIATED(xkp_source) .OR. .NOT. ASSOCIATED(wkp_source)) THEN
374 cpabort("Could not access the SCF GENERAL k-point set for the Wannier90 export.")
375 END IF
376 num_kpts = SIZE(wkp_source)
377 ALLOCATE (kpt_latt(3, num_kpts))
378 kpt_latt(1:3, 1:num_kpts) = xkp_source(1:3, 1:num_kpts)
379 IF (mp_grid_explicit) THEN
380 IF (mp_grid(1)*mp_grid(2)*mp_grid(3) /= num_kpts) THEN
381 cpabort("WANNIER90%MP_GRID must contain exactly as many points as the SCF GENERAL mesh.")
382 END IF
383 ELSE
384 CALL infer_wannier_mp_grid(kpt_latt, mp_grid, mp_grid_valid)
385 IF (.NOT. mp_grid_valid) THEN
386 cpabort("Could not infer WANNIER90%MP_GRID from the SCF GENERAL mesh.")
387 END IF
388 END IF
389 wkp_ref = 1.0_dp/real(num_kpts, kind=dp)
390 DO i = 1, num_kpts
391 IF (abs(wkp_source(i) - wkp_ref) > 1.0e-10_dp) THEN
392 cpabort("WANNIER90%KPOINTS_SOURCE SCF requires equally weighted GENERAL k-points.")
393 END IF
394 END DO
395 kpoint%kp_scheme = "GENERAL"
396 kpoint%symmetry = .false.
397 kpoint%nkp_grid(1:3) = mp_grid(1:3)
398 kpoint%verbose = .false.
399 kpoint%full_grid = .true.
400 kpoint%eps_geo = 1.0e-6_dp
401 kpoint%use_real_wfn = .false.
402 kpoint%parallel_group_size = para_env%num_pe
403 kpoint%nkp = num_kpts
404 ALLOCATE (kpoint%xkp(3, num_kpts), kpoint%wkp(num_kpts))
405 kpoint%xkp(1:3, 1:num_kpts) = xkp_source(1:3, 1:num_kpts)
406 kpoint%wkp(1:num_kpts) = wkp_ref
407 IF (input_kpoint_symmetry .AND. .NOT. input_full_grid .AND. iw > 0) THEN
408 WRITE (iw, '(T2,A)') &
409 "WANNIER90| SCF k-points are symmetry-reduced; using the full input GENERAL mesh."
410 IF (reuse_scf_mos) THEN
411 WRITE (iw, '(T2,A)') &
412 "WANNIER90| CP2K will try to reconstruct the full-mesh MOs from the SCF orbitals."
413 ELSE
414 WRITE (iw, '(T2,A)') &
415 "WANNIER90| The full exported mesh is diagonalized for the Wannier90 files."
416 END IF
417 END IF
418
419 CASE DEFAULT
420 cpabort("WANNIER90%KPOINTS_SOURCE SCF does not support this DFT%KPOINTS scheme.")
421 END SELECT
422 CASE DEFAULT
423 cpabort("Unknown WANNIER90%KPOINTS_SOURCE setting.")
424 END SELECT
425 ! number of bands in calculation
426 CALL get_qs_env(qs_env, mos=mos)
427 CALL get_mo_set(mo_set=mos(1), nao=nao, nmo=num_bands_tot)
428 num_bands_tot = min(nao, num_bands_tot + nadd)
429 num_bands = num_bands_tot
430 IF (use_bloch_phases .AND. num_wann /= num_bands) THEN
431 cpabort("WANNIER90%USE_BLOCH_PHASES requires WANNIER_FUNCTIONS to match the number of bands.")
432 END IF
433 num_atoms = SIZE(particle_set)
434 ALLOCATE (atoms_cart(3, num_atoms))
435 ALLOCATE (atom_symbols(num_atoms))
436 DO i = 1, num_atoms
437 atoms_cart(1:3, i) = particle_set(i)%r(1:3)
438 CALL get_atomic_kind(particle_set(i)%atomic_kind, element_symbol=asym)
439 atom_symbols(i) = asym
440 END DO
441 gamma_only = .false.
442 spinors = .false.
443 ! output
444 ALLOCATE (nnlist(num_kpts, num_nnmax))
445 ALLOCATE (nncell(3, num_kpts, num_nnmax))
446 nnlist(:, :) = 0
447 nncell(:, :, :) = 0
448 nntot = 0
449
450 IF (iw > 0) THEN
451 ! setup
452 CALL wannier_setup(mp_grid, num_kpts, real_lattice, recip_lattice, &
453 kpt_latt, nntot, nnlist, nncell, iw)
454 END IF
455
456 CALL get_qs_env(qs_env, para_env=para_env)
457 CALL para_env%sum(nntot)
458 CALL para_env%sum(nnlist)
459 CALL para_env%sum(nncell)
460
461 IF (para_env%is_source()) THEN
462 ! Write the Wannier90 input file "seed_name.win"
463 WRITE (filename, '(A,A)') trim(seed_name), ".win"
464 CALL open_file(filename, unit_number=iunit, file_status="UNKNOWN", file_action="WRITE")
465 !
466 CALL m_timestamp(timestamp)
467 WRITE (iunit, "(A)") "! Wannier90 input file generated by CP2K "
468 WRITE (iunit, "(A,/)") "! Creation date "//timestamp
469 !
470 WRITE (iunit, "(A,I5)") "num_wann = ", num_wann
471 IF (num_bands /= num_wann .OR. use_bloch_phases) THEN
472 WRITE (iunit, "(A,I5)") "num_bands = ", num_bands
473 END IF
474 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 WRITE (iunit, "(A)") "! CP2K writes identity projections for Bloch-phase complete subspaces."
478 END IF
479 WRITE (iunit, "(/,A,/)") "length_unit = bohr "
480 WRITE (iunit, "(/,A,/)") "! System"
481 WRITE (iunit, "(/,A)") "begin unit_cell_cart"
482 WRITE (iunit, "(A)") "bohr"
483 DO i = 1, 3
484 WRITE (iunit, "(3F12.6)") cell%hmat(i, 1:3)
485 END DO
486 WRITE (iunit, "(A,/)") "end unit_cell_cart"
487 WRITE (iunit, "(/,A)") "begin atoms_cart"
488 WRITE (iunit, "(A)") "bohr"
489 DO i = 1, num_atoms
490 WRITE (iunit, "(A,3F15.10)") atom_symbols(i), atoms_cart(1:3, i)
491 END DO
492 WRITE (iunit, "(A,/)") "end atoms_cart"
493 WRITE (iunit, "(/,A,/)") "! Kpoints"
494 WRITE (iunit, "(/,A,3I6/)") "mp_grid = ", mp_grid(1:3)
495 WRITE (iunit, "(A)") "begin kpoints"
496 DO i = 1, num_kpts
497 WRITE (iunit, "(3F12.6)") kpt_latt(1:3, i)
498 END DO
499 WRITE (iunit, "(A)") "end kpoints"
500 CALL close_file(iunit)
501 IF (use_bloch_phases) THEN
502 WRITE (filename, '(A,A)') trim(seed_name), ".amn"
503 CALL open_file(filename, unit_number=iunit, file_status="UNKNOWN", file_action="WRITE")
504 WRITE (iunit, "(A)") "! Wannier90 identity projections generated by CP2K"
505 WRITE (iunit, "(3I8)") num_bands, num_kpts, num_wann
506 DO ik = 1, num_kpts
507 DO ib2 = 1, num_wann
508 DO ib1 = 1, num_bands
509 IF (ib1 == ib2) THEN
510 WRITE (iunit, "(3I8,2E30.14)") ib1, ib2, ik, 1.0_dp, 0.0_dp
511 ELSE
512 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 CALL close_file(iunit)
518 END IF
519 ELSE
520 iunit = -1
521 END IF
522
523 ! calculate bands
524 NULLIFY (qs_env_kp)
525 IF (kpoints_source == w90_kpoints_mp_grid .AND. input_kpoint_symmetry .AND. iw > 0) THEN
526 WRITE (iw, '(T2,A)') &
527 "WANNIER90| Atomic k-point symmetry from the SCF calculation is not reused."
528 WRITE (iw, '(T2,A)') &
529 "WANNIER90| A full Monkhorst-Pack grid is generated for the Wannier90 interface."
530 END IF
531 IF (do_kpoints) THEN
532 ! we already do kpoints
533 qs_env_kp => qs_env
534 ELSE
535 ! we start from gamma point only
536 ALLOCATE (qs_env_kp)
537 CALL create_kp_from_gamma(qs_env, qs_env_kp)
538 END IF
539 IF (iw > 0) THEN
540 WRITE (unit=iw, fmt="(/,T2,A)") "Start K-Point Calculation ..."
541 END IF
542 CALL get_qs_env(qs_env=qs_env_kp, para_env=para_env, blacs_env=blacs_env)
543 CALL kpoint_env_initialize(kpoint, para_env, blacs_env)
544 CALL kpoint_initialize_mos(kpoint, mos, nadd)
545 CALL kpoint_initialize_mo_set(kpoint)
546 !
547 CALL get_qs_env(qs_env=qs_env_kp, sab_orb=sab_nl, dft_control=dft_control)
548 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 scf_env=scf_env, scf_control=scf_control)
552 full_mesh_diagonalized = .false.
553 reused_scf_mos = .false.
554 reuse_reason = ""
555 aligned_degenerate_blocks = 0
556 aligned_degenerate_max_size = 0
557 aligned_degenerate_min_svalue = 0.0_dp
558 IF (reuse_scf_mos) THEN
559 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 diis_step)
562 IF (validate_reuse_scf_mos) THEN
563 IF (iw > 0) THEN
564 WRITE (iw, '(T2,A)') &
565 "WANNIER90| Validating SCF MO reuse against a full-mesh diagonalization reference."
566 END IF
567 CALL do_general_diag_kp(matrix_ks, matrix_s, kpoint, scf_env, scf_control, .false., diis_step)
568 full_mesh_diagonalized = .true.
569 nspins = dft_control%nspins
570 CALL save_wannier90_mo_snapshot(kpoint, nspins, para_env, reference_mo_real, &
571 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 reuse_candidate_residual)
578 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 "WANNIER90| Best atom/AO candidate subspace deviation ", &
581 reuse_candidate_deviation, ", minimum singular value ", &
582 reuse_candidate_min_svalue, ", max metric deviation ", &
583 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 aligned_degenerate_min_svalue)
590 IF (validate_reuse_scf_mos) THEN
591 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 validation_eigenvalue_deviation)
597 IF (iw > 0) THEN
598 WRITE (iw, '(T2,A,ES10.3,A,ES10.3,A,ES10.3)') &
599 "WANNIER90| Reused MO validation: subspace deviation ", &
600 validation_subspace_deviation, ", minimum singular value ", &
601 validation_min_svalue, ", eigenvalue deviation ", &
602 validation_eigenvalue_deviation
603 END IF
604 IF (.NOT. validate_reuse_ok) THEN
605 reused_scf_mos = .false.
606 WRITE (reuse_reason, "(A,ES10.3,A,ES10.3)") &
607 "validation failed: dS=", &
608 validation_subspace_deviation, ", dE=", validation_eigenvalue_deviation
609 END IF
610 END IF
611 IF (.NOT. reused_scf_mos) THEN
612 CALL restore_wannier90_mo_snapshot(kpoint, reference_mo_real, reference_mo_imag, &
613 reference_eigenvalues)
614 END IF
615 END IF
616 IF (iw > 0) THEN
617 IF (reused_scf_mos) THEN
618 WRITE (iw, '(T2,A)') &
619 "WANNIER90| Reused SCF MO coefficients for the Wannier90 full k-point mesh."
620 IF (use_bloch_phases) THEN
621 WRITE (iw, '(T2,A)') &
622 "WANNIER90| Wrote identity projections for Bloch-phase complete band subspaces."
623 WRITE (iw, '(T2,A,3F10.6)') &
624 "WANNIER90| Applied Bloch phase gauge to reused overlaps around fractional center", &
625 phase_center(1:3)
626 END IF
627 IF (aligned_degenerate_blocks > 0) THEN
628 WRITE (iw, '(T2,A,I0,A,I0,A,ES10.3)') &
629 "WANNIER90| Ritz-stabilized ", aligned_degenerate_blocks, &
630 " degenerate SCF MO subspace(s) with S(k),H(k); largest block has ", &
631 aligned_degenerate_max_size, " band(s), min metric eigenvalue ", &
632 aligned_degenerate_min_svalue
633 END IF
634 ELSE
635 WRITE (iw, '(T2,A,A)') &
636 "WANNIER90| Could not reuse SCF MOs: ", trim(reuse_reason)
637 WRITE (iw, '(T2,A)') &
638 "WANNIER90| Falling back to full-mesh diagonalization for the Wannier90 files."
639 END IF
640 END IF
641 END IF
642 IF (.NOT. reused_scf_mos .AND. .NOT. full_mesh_diagonalized) THEN
643 CALL do_general_diag_kp(matrix_ks, matrix_s, kpoint, scf_env, scf_control, .false., diis_step)
644 END IF
645 IF (ALLOCATED(reference_mo_real)) DEALLOCATE (reference_mo_real)
646 IF (ALLOCATED(reference_mo_imag)) DEALLOCATE (reference_mo_imag)
647 IF (ALLOCATED(reference_eigenvalues)) DEALLOCATE (reference_eigenvalues)
648 !
649 IF (iw > 0) THEN
650 WRITE (iw, '(T69,A)') "... Finished"
651 END IF
652 !
653 ! Calculate and print Overlaps
654 !
655 IF (para_env%is_source()) THEN
656 WRITE (filename, '(A,A)') trim(seed_name), ".mmn"
657 CALL open_file(filename, unit_number=iunit, file_status="UNKNOWN", file_action="WRITE")
658 CALL m_timestamp(timestamp)
659 WRITE (iunit, "(A)") "! Wannier90 file generated by CP2K "//timestamp
660 WRITE (iunit, "(3I8)") num_bands, num_kpts, nntot
661 ELSE
662 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 ALLOCATE (nblist(num_kpts, nntot))
667 ALLOCATE (b_latt(3, num_kpts*nntot))
668 nblist(:, :) = 0
669 nbs = 0
670 DO ik = 1, num_kpts
671 DO i = 1, nntot
672 bvec(1:3) = kpt_latt(1:3, nnlist(ik, i)) - kpt_latt(1:3, ik) + nncell(1:3, ik, i)
673 ibs = 0
674 DO k = 1, nbs
675 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 IF (sum(abs(bvec(1:3) + b_latt(1:3, k))) < 1.e-6_dp) THEN
680 ibs = -k
681 EXIT
682 END IF
683 END DO
684 IF (ibs /= 0) THEN
685 ! old lattice vector
686 nblist(ik, i) = ibs
687 ELSE
688 ! new lattice vector
689 nbs = nbs + 1
690 b_latt(1:3, nbs) = bvec(1:3)
691 nblist(ik, i) = nbs
692 END IF
693 END DO
694 END DO
695 ! calculate all the operator matrices (a|bvec|b)
696 ALLOCATE (berry_matrix(nbs))
697 DO i = 1, nbs
698 NULLIFY (berry_matrix(i)%cosmat)
699 NULLIFY (berry_matrix(i)%sinmat)
700 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 berry_matrix(i)%sinmat, bvec)
703 END DO
704 ! work matrices for MOs (all group)
705 kp => kpoint%kp_env(1)%kpoint_env
706 CALL get_mo_set(kp%mos(1, 1), nmo=nmo)
707 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 context=blacs_env)
712 DO i = 1, 2
713 CALL cp_fm_create(fmk1(i), matrix_struct_work)
714 CALL cp_fm_create(fmk2(i), matrix_struct_work)
715 END DO
716 CALL cp_cfm_create(fmk1_cfm, matrix_struct_work)
717 CALL cp_cfm_create(fmk2_cfm, matrix_struct_work)
718 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 context=blacs_env)
723 CALL cp_fm_create(mat_real, matrix_struct_ao)
724 CALL cp_fm_create(mat_imag, matrix_struct_ao)
725 CALL cp_cfm_create(omat_cfm, matrix_struct_ao)
726 ! work matrices for Mmn(k,b) integrals
727 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 context=blacs_env)
732 CALL cp_fm_create(mmn_real, matrix_struct_mmn)
733 CALL cp_fm_create(mmn_imag, matrix_struct_mmn)
734 CALL cp_cfm_create(mmn_cfm, matrix_struct_mmn)
735 ! allocate some work matrices
736 ALLOCATE (rmatrix, cmatrix, rmatrix_full, cmatrix_full)
737 CALL dbcsr_create(rmatrix, template=matrix_s(1, 1)%matrix, &
738 matrix_type=dbcsr_type_symmetric)
739 CALL dbcsr_create(cmatrix, template=matrix_s(1, 1)%matrix, &
740 matrix_type=dbcsr_type_antisymmetric)
741 CALL dbcsr_create(rmatrix_full, template=matrix_s(1, 1)%matrix, &
742 matrix_type=dbcsr_type_no_symmetry)
743 CALL dbcsr_create(cmatrix_full, template=matrix_s(1, 1)%matrix, &
744 matrix_type=dbcsr_type_no_symmetry)
745 CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
746 CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
747 !
748 CALL get_kpoint_info(kpoint=kpoint, cell_to_index=cell_to_index)
749 NULLIFY (fmdummy)
750 nspins = dft_control%nspins
751 DO ispin = 1, nspins
752 ! loop over all k-points
753 DO ik = 1, num_kpts
754 ! get the MO coefficients for this k-point
755 my_kpgrp = (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2))
756 IF (my_kpgrp) THEN
757 ikk = ik - kpoint%kp_range(1) + 1
758 kp => kpoint%kp_env(ikk)%kpoint_env
759 cpassert(SIZE(kp%mos, 1) == 2)
760 fmr => kp%mos(1, ispin)%mo_coeff
761 fmi => kp%mos(2, ispin)%mo_coeff
762 CALL cp_fm_copy_general(fmr, fmk1(1), para_env)
763 CALL cp_fm_copy_general(fmi, fmk1(2), para_env)
764 ELSE
765 NULLIFY (fmr, fmi, kp)
766 CALL cp_fm_copy_general(fmdummy, fmk1(1), para_env)
767 CALL cp_fm_copy_general(fmdummy, fmk1(2), para_env)
768 END IF
769 CALL cp_fm_to_cfm(fmk1(1), fmk1(2), fmk1_cfm)
770 ! loop over all connected neighbors
771 DO i = 1, nntot
772 ! get the MO coefficients for the connected k-point
773 ik2 = nnlist(ik, i)
774 mygrp = (ik2 >= kpoint%kp_range(1) .AND. ik2 <= kpoint%kp_range(2))
775 IF (mygrp) THEN
776 ikk = ik2 - kpoint%kp_range(1) + 1
777 kp => kpoint%kp_env(ikk)%kpoint_env
778 cpassert(SIZE(kp%mos, 1) == 2)
779 fmr => kp%mos(1, ispin)%mo_coeff
780 fmi => kp%mos(2, ispin)%mo_coeff
781 CALL cp_fm_copy_general(fmr, fmk2(1), para_env)
782 CALL cp_fm_copy_general(fmi, fmk2(2), para_env)
783 ELSE
784 NULLIFY (fmr, fmi, kp)
785 CALL cp_fm_copy_general(fmdummy, fmk2(1), para_env)
786 CALL cp_fm_copy_general(fmdummy, fmk2(2), para_env)
787 END IF
788 CALL cp_fm_to_cfm(fmk2(1), fmk2(2), fmk2_cfm)
789 !
790 ! transfer realspace overlaps to connected k-point
791 ibs = nblist(ik, i)
792 ksign = sign(1.0_dp, real(ibs, kind=dp))
793 ibs = abs(ibs)
794 CALL dbcsr_set(rmatrix, 0.0_dp)
795 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 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 is_complex=.true., rs_sign=ksign)
802 !
803 ! calculate M_(mn)^(k,b) = C(k)^H O(k,b) C(k+b)
804 CALL dbcsr_desymmetrize(rmatrix, rmatrix_full)
805 CALL dbcsr_desymmetrize(cmatrix, cmatrix_full)
806 CALL copy_dbcsr_to_fm(rmatrix_full, mat_real)
807 CALL copy_dbcsr_to_fm(cmatrix_full, mat_imag)
808 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 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 fmk1_cfm, tmp_cfm, cmplx(0.0_dp, 0.0_dp, kind=dp), mmn_cfm)
813 CALL cp_cfm_to_fm(mmn_cfm, mmn_real, mmn_imag)
814 !
815 ! write to output file
816 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 phase_center(1:3))
820 gauge_real = cos(gauge_arg)
821 gauge_imag = sin(gauge_arg)
822 ELSE
823 gauge_real = 1.0_dp
824 gauge_imag = 0.0_dp
825 END IF
826 IF (para_env%is_source()) THEN
827 WRITE (iunit, "(2I8,3I5)") ik, ik2, nncell(1:3, ik, i)
828 END IF
829 DO ib2 = 1, nmo
830 DO ib1 = 1, nmo
831 CALL cp_fm_get_element(mmn_real, ib1, ib2, rmmn)
832 CALL cp_fm_get_element(mmn_imag, ib1, ib2, cmmn)
833 gauge_tmp = gauge_real*rmmn - gauge_imag*cmmn
834 cmmn = gauge_imag*rmmn + gauge_real*cmmn
835 rmmn = gauge_tmp
836 IF (para_env%is_source()) THEN
837 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 DO i = 1, nbs
846 CALL dbcsr_deallocate_matrix_set(berry_matrix(i)%cosmat)
847 CALL dbcsr_deallocate_matrix_set(berry_matrix(i)%sinmat)
848 END DO
849 DEALLOCATE (berry_matrix)
850 CALL cp_fm_struct_release(matrix_struct_work)
851 DO i = 1, 2
852 CALL cp_fm_release(fmk1(i))
853 CALL cp_fm_release(fmk2(i))
854 END DO
855 CALL cp_cfm_release(fmk1_cfm)
856 CALL cp_cfm_release(fmk2_cfm)
857 CALL cp_cfm_release(tmp_cfm)
858 CALL cp_fm_struct_release(matrix_struct_ao)
859 CALL cp_fm_release(mat_real)
860 CALL cp_fm_release(mat_imag)
861 CALL cp_cfm_release(omat_cfm)
862 CALL cp_fm_struct_release(matrix_struct_mmn)
863 CALL cp_fm_release(mmn_real)
864 CALL cp_fm_release(mmn_imag)
865 CALL cp_cfm_release(mmn_cfm)
866 CALL dbcsr_deallocate_matrix(rmatrix)
867 CALL dbcsr_deallocate_matrix(cmatrix)
868 CALL dbcsr_deallocate_matrix(rmatrix_full)
869 CALL dbcsr_deallocate_matrix(cmatrix_full)
870 !
871 IF (para_env%is_source()) THEN
872 CALL close_file(iunit)
873 END IF
874 !
875 ! Calculate and print Projections
876 !
877 ! Print eigenvalues
878 nspins = dft_control%nspins
879 kp => kpoint%kp_env(1)%kpoint_env
880 CALL get_mo_set(kp%mos(1, 1), nmo=nmo)
881 ALLOCATE (eigval(nmo))
882 CALL get_kpoint_info(kpoint, nkp=nkp, kp_range=kp_range, xkp=xkp)
883 IF (para_env%is_source()) THEN
884 WRITE (filename, '(A,A)') trim(seed_name), ".eig"
885 CALL open_file(filename, unit_number=iunit, file_status="UNKNOWN", file_action="WRITE")
886 ELSE
887 iunit = -1
888 END IF
889 !
890 DO ik = 1, nkp
891 my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
892 DO ispin = 1, nspins
893 IF (my_kpgrp) THEN
894 ikpgr = ik - kp_range(1) + 1
895 kp => kpoint%kp_env(ikpgr)%kpoint_env
896 CALL get_mo_set(kp%mos(1, ispin), eigenvalues=eigenvalues)
897 eigval(1:nmo) = eigenvalues(1:nmo)
898 ELSE
899 eigval(1:nmo) = 0.0_dp
900 END IF
901 CALL kpoint%para_env_inter_kp%sum(eigval)
902 eigval(1:nmo) = eigval(1:nmo)*evolt
903 ! output
904 IF (iunit > 0) THEN
905 DO ib = 1, nmo
906 WRITE (iunit, "(2I8,F24.14)") ib, ik, eigval(ib)
907 END DO
908 END IF
909 END DO
910 END DO
911 IF (para_env%is_source()) THEN
912 CALL close_file(iunit)
913 END IF
914 !
915 ! clean up
916 DEALLOCATE (kpt_latt, atoms_cart, atom_symbols, eigval)
917 DEALLOCATE (nnlist, nncell)
918 DEALLOCATE (nblist, b_latt)
919 IF (nexcl > 0) THEN
920 DEALLOCATE (exclude_bands)
921 END IF
922 IF (do_kpoints) THEN
923 NULLIFY (qs_env_kp)
924 ELSE
925 CALL qs_env_release(qs_env_kp)
926 DEALLOCATE (qs_env_kp)
927 NULLIFY (qs_env_kp)
928 END IF
929
930 CALL kpoint_release(kpoint)
931
932 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 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_buffer, occupation_buffer, &
978 source_eigenvalues_buffer, &
979 source_occupation_buffer
980 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 success = .false.
991 reason = ""
992 aligned_degenerate_blocks = 0
993 aligned_degenerate_max_size = 0
994 aligned_degenerate_min_svalue = huge(1.0_dp)
995 NULLIFY (matrix_struct_source, matrix_struct_work, src_fmr, src_fmi, dst_fmr, dst_fmi)
996
997 IF (.NOT. ASSOCIATED(kpoint)) THEN
998 reason = "internal Wannier90 k-point object is not available"
999 RETURN
1000 END IF
1001 IF (.NOT. ASSOCIATED(qs_kpoint)) THEN
1002 reason = "SCF k-point object is not available"
1003 RETURN
1004 END IF
1005 IF (.NOT. ASSOCIATED(kpoint%kp_env) .OR. .NOT. ASSOCIATED(qs_kpoint%kp_env)) THEN
1006 reason = "k-point MO environments are not initialized"
1007 RETURN
1008 END IF
1009 IF (.NOT. ASSOCIATED(kpoint%blacs_env)) THEN
1010 reason = "Wannier90 k-point BLACS environment is not initialized"
1011 RETURN
1012 END IF
1013
1014 CALL build_wannier90_scf_mapping(kpoint, qs_kpoint, source_kpoint, sym_index, ok, reason)
1015 IF (.NOT. ok) RETURN
1016 nsymmetry = count(sym_index > 0)
1017
1018 kp => kpoint%kp_env(1)%kpoint_env
1019 nspins = SIZE(kp%mos, 2)
1020 IF (SIZE(kp%mos, 1) < 2) THEN
1021 reason = "Wannier90 export k-point MOs are not complex-valued"
1022 DEALLOCATE (source_kpoint, sym_index)
1023 RETURN
1024 END IF
1025 CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo)
1026
1027 kp_source => qs_kpoint%kp_env(1)%kpoint_env
1028 IF (SIZE(kp_source%mos, 1) < 2) THEN
1029 reason = "SCF MOs are real-valued; complex symmetry phases cannot be reconstructed"
1030 DEALLOCATE (source_kpoint, sym_index)
1031 RETURN
1032 END IF
1033 CALL get_mo_set(kp_source%mos(1, 1), nao=nao_src, nmo=nmo_src)
1034 CALL para_env%max(nao_src)
1035 CALL para_env%max(nmo_src)
1036 IF (nao_src /= nao) THEN
1037 reason = "SCF and Wannier90 MO bases have different AO dimensions"
1038 DEALLOCATE (source_kpoint, sym_index)
1039 RETURN
1040 END IF
1041 IF (nmo_src < nmo) THEN
1042 reason = "SCF MO set has fewer bands than the Wannier90 export"
1043 DEALLOCATE (source_kpoint, sym_index)
1044 RETURN
1045 END IF
1046 source_window = nmo_src > nmo
1047 degenerate_band_tol = 1.0e-8_dp
1048 CALL get_kpoint_info(qs_kpoint, kp_range=source_kp_range)
1049 IF (source_kp_range(1) /= 1 .OR. source_kp_range(2) /= qs_kpoint%nkp) THEN
1050 reason = "SCF k-point symmetry data are distributed over k-point parallel groups"
1051 DEALLOCATE (source_kpoint, sym_index)
1052 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 IF (nsymmetry > 0 .AND. nmo_src > nmo) THEN
1058 local_min_band_gap = huge(1.0_dp)
1059 min_gap_band = nmo
1060 min_gap_kpoint = 0
1061 min_gap_spin = 0
1062 DO ikred = source_kp_range(1), source_kp_range(2)
1063 ikpgr = ikred - source_kp_range(1) + 1
1064 kp_source => qs_kpoint%kp_env(ikpgr)%kpoint_env
1065 DO ispin = 1, nspins
1066 CALL get_mo_set(kp_source%mos(1, ispin), eigenvalues=eigenvalues)
1067 band_gap = abs(eigenvalues(nmo + 1) - eigenvalues(nmo))
1068 IF (band_gap < local_min_band_gap) THEN
1069 local_min_band_gap = band_gap
1070 min_gap_band = nmo
1071 min_gap_kpoint = ikred
1072 min_gap_spin = ispin
1073 END IF
1074 END DO
1075 END DO
1076 min_band_gap = local_min_band_gap
1077 CALL para_env%min(min_band_gap)
1078 IF (abs(local_min_band_gap - min_band_gap) > degenerate_band_tol*epsilon(1.0_dp)) THEN
1079 min_gap_kpoint = 0
1080 min_gap_spin = 0
1081 END IF
1082 CALL para_env%max(min_gap_kpoint)
1083 CALL para_env%max(min_gap_spin)
1084 CALL para_env%max(min_gap_band)
1085 IF (min_band_gap < degenerate_band_tol) THEN
1086 WRITE (reason, "(A,ES9.2,A,I0,A,I0,A,I0)") &
1087 "degenerate atom/AO W90 reuse guarded: edge gap=", min_band_gap, ", k=", &
1088 min_gap_kpoint, ", s=", min_gap_spin, ", nband=", min_gap_band
1089 DEALLOCATE (source_kpoint, sym_index)
1090 RETURN
1091 END IF
1092 END IF
1093 blacs_env => kpoint%blacs_env
1094 CALL cp_fm_struct_create(matrix_struct_work, nrow_global=nao, ncol_global=nmo, &
1095 para_env=para_env, context=blacs_env)
1096 CALL cp_fm_create(src_real, matrix_struct_work)
1097 CALL cp_fm_create(src_imag, matrix_struct_work)
1098 CALL cp_fm_create(dst_real, matrix_struct_work)
1099 CALL cp_fm_create(dst_imag, matrix_struct_work)
1100 IF (source_window) THEN
1101 CALL cp_fm_struct_create(matrix_struct_source, nrow_global=nao, ncol_global=nmo_src, &
1102 para_env=para_env, context=blacs_env)
1103 CALL cp_fm_create(src_real_full, matrix_struct_source)
1104 CALL cp_fm_create(src_imag_full, matrix_struct_source)
1105 CALL cp_fm_create(dst_real_full, matrix_struct_source)
1106 CALL cp_fm_create(dst_imag_full, matrix_struct_source)
1107 END IF
1108 ALLOCATE (eigenvalues_buffer(nmo), occupation_buffer(nmo), &
1109 source_eigenvalues_buffer(nmo_src), source_occupation_buffer(nmo_src))
1110
1111 CALL get_kpoint_info(kpoint, kp_range=kp_range)
1112 CALL get_kpoint_info(qs_kpoint, kp_range=source_kp_range)
1113
1114 DO ik = 1, kpoint%nkp
1115 ikred = source_kpoint(ik)
1116 my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
1117 my_source_kpgrp = (ikred >= source_kp_range(1) .AND. ikred <= source_kp_range(2))
1118 DO ispin = 1, nspins
1119 source_eigenvalues_buffer(1:nmo_src) = 0.0_dp
1120 source_occupation_buffer(1:nmo_src) = 0.0_dp
1121 IF (my_source_kpgrp) THEN
1122 ikpgr = ikred - source_kp_range(1) + 1
1123 kp_source => qs_kpoint%kp_env(ikpgr)%kpoint_env
1124 src_fmr => kp_source%mos(1, ispin)%mo_coeff
1125 src_fmi => kp_source%mos(2, ispin)%mo_coeff
1126 CALL get_mo_set(kp_source%mos(1, ispin), eigenvalues=eigenvalues, &
1127 occupation_numbers=occupation)
1128 source_eigenvalues_buffer(1:nmo_src) = eigenvalues(1:nmo_src)
1129 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 source_owner_count = 1.0_dp
1135 ELSE
1136 source_owner_count = 0.0_dp
1137 END IF
1138 CALL para_env%sum(source_owner_count)
1139 CALL para_env%sum(source_eigenvalues_buffer)
1140 CALL para_env%sum(source_occupation_buffer)
1141 IF (source_owner_count > 0.0_dp) THEN
1142 source_eigenvalues_buffer(1:nmo_src) = &
1143 source_eigenvalues_buffer(1:nmo_src)/source_owner_count
1144 source_occupation_buffer(1:nmo_src) = &
1145 source_occupation_buffer(1:nmo_src)/source_owner_count
1146 END IF
1147 eigenvalues_buffer(1:nmo) = source_eigenvalues_buffer(1:nmo)
1148 occupation_buffer(1:nmo) = source_occupation_buffer(1:nmo)
1149 IF (source_window) THEN
1150 CALL cp_fm_copy_general(src_fmr, src_real_full, para_env)
1151 CALL cp_fm_copy_general(src_fmi, src_imag_full, para_env)
1152 CALL copy_wannier90_mo_window(src_real_full, src_real, nmo)
1153 CALL copy_wannier90_mo_window(src_imag_full, src_imag, nmo)
1154 ELSE
1155 CALL cp_fm_copy_general(src_fmr, src_real, para_env)
1156 CALL cp_fm_copy_general(src_fmi, src_imag, para_env)
1157 END IF
1158
1159 ok = .false.
1160 reason = ""
1161 aligned_blocks = 0
1162 aligned_max_size = 0
1163 aligned_min_svalue = 0.0_dp
1164 IF (sym_index(ik) > 0) THEN
1165 kpsym => qs_kpoint%kp_sym(ikred)%kpoint_sym
1166 IF (ASSOCIATED(kpsym)) THEN
1167 best_reason = ""
1168 best_residual = huge(1.0_dp)
1169 num_candidates = 0
1170 ! Little-group operations can reach the same target k-point; keep the first valid eigenspace.
1171 DO isym_try = 1, kpsym%nwred
1172 IF (.NOT. kpoint_same_periodic(kpoint%xkp(1:3, ik), &
1173 kpsym%xkp(1:3, isym_try))) cycle
1174 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 candidate_reason)
1178 IF (.NOT. ok) THEN
1179 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 aligned_min_svalue, candidate_residual)
1189 IF (candidate_residual < best_residual) THEN
1190 best_residual = candidate_residual
1191 best_reason = candidate_reason
1192 END IF
1193 IF (.NOT. ok) THEN
1194 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 ikred, isym_try, para_env, ok, candidate_reason)
1198 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 candidate_residual)
1206 IF (candidate_residual < best_residual) THEN
1207 best_residual = candidate_residual
1208 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 candidate_residual)
1218 IF (candidate_residual < best_residual) THEN
1219 best_residual = candidate_residual
1220 best_reason = candidate_reason
1221 END IF
1222 END IF
1223 END IF
1224 IF (ok) THEN
1225 aligned_blocks = candidate_aligned_blocks
1226 aligned_max_size = candidate_aligned_max_size
1227 sym_index(ik) = isym_try
1228 EXIT
1229 END IF
1230 reason = candidate_reason
1231 END DO
1232 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 "atom/AO W90 guarded: best/", num_candidates, "=", best_residual, &
1235 " k=", ik, " ", trim(best_reason)
1236 ELSE IF (.NOT. ok .AND. num_candidates == 0) THEN
1237 reason = "no matching SCF symmetry operation candidate"
1238 END IF
1239 ELSE
1240 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 ikred, sym_index(ik), para_env, ok, reason)
1245 END IF
1246 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 aligned_min_svalue, candidate_residual)
1254 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 para_env, ok, reason)
1258 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 candidate_residual)
1265 END IF
1266 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 source_window_min_svalue, candidate_residual)
1272 END IF
1273 END IF
1274 IF (.NOT. ok) THEN
1275 CALL cp_fm_release(src_real)
1276 CALL cp_fm_release(src_imag)
1277 CALL cp_fm_release(dst_real)
1278 CALL cp_fm_release(dst_imag)
1279 CALL cp_fm_struct_release(matrix_struct_work)
1280 IF (source_window) THEN
1281 CALL cp_fm_release(src_real_full)
1282 CALL cp_fm_release(src_imag_full)
1283 CALL cp_fm_release(dst_real_full)
1284 CALL cp_fm_release(dst_imag_full)
1285 CALL cp_fm_struct_release(matrix_struct_source)
1286 END IF
1287 DEALLOCATE (source_kpoint, sym_index, eigenvalues_buffer, occupation_buffer, &
1288 source_eigenvalues_buffer, source_occupation_buffer)
1289 RETURN
1290 END IF
1291 IF (sym_index(ik) /= 0) THEN
1292 aligned_degenerate_blocks = aligned_degenerate_blocks + aligned_blocks
1293 aligned_degenerate_max_size = max(aligned_degenerate_max_size, aligned_max_size)
1294 IF (aligned_blocks > 0) THEN
1295 aligned_degenerate_min_svalue = min(aligned_degenerate_min_svalue, aligned_min_svalue)
1296 END IF
1297 END IF
1298
1299 IF (my_kpgrp) THEN
1300 ikpgr = ik - kp_range(1) + 1
1301 kp => kpoint%kp_env(ikpgr)%kpoint_env
1302 dst_fmr => kp%mos(1, ispin)%mo_coeff
1303 dst_fmi => kp%mos(2, ispin)%mo_coeff
1304 CALL get_mo_set(kp%mos(1, ispin), eigenvalues=eigenvalues, &
1305 occupation_numbers=occupation)
1306 eigenvalues(1:nmo) = eigenvalues_buffer(1:nmo)
1307 occupation(1:nmo) = occupation_buffer(1:nmo)
1308 CALL get_mo_set(kp%mos(2, ispin), eigenvalues=eigenvalues, &
1309 occupation_numbers=occupation)
1310 IF (ASSOCIATED(eigenvalues)) eigenvalues(1:nmo) = eigenvalues_buffer(1:nmo)
1311 IF (ASSOCIATED(occupation)) occupation(1:nmo) = occupation_buffer(1:nmo)
1312 ELSE
1313 NULLIFY (dst_fmr, dst_fmi)
1314 END IF
1315 CALL cp_fm_copy_general(dst_real, dst_fmr, para_env)
1316 CALL cp_fm_copy_general(dst_imag, dst_fmi, para_env)
1317 END DO
1318 END DO
1319
1320 CALL cp_fm_release(src_real)
1321 CALL cp_fm_release(src_imag)
1322 CALL cp_fm_release(dst_real)
1323 CALL cp_fm_release(dst_imag)
1324 CALL cp_fm_struct_release(matrix_struct_work)
1325 IF (source_window) THEN
1326 CALL cp_fm_release(src_real_full)
1327 CALL cp_fm_release(src_imag_full)
1328 CALL cp_fm_release(dst_real_full)
1329 CALL cp_fm_release(dst_imag_full)
1330 CALL cp_fm_struct_release(matrix_struct_source)
1331 END IF
1332 DEALLOCATE (source_kpoint, sym_index, eigenvalues_buffer, occupation_buffer, &
1333 source_eigenvalues_buffer, source_occupation_buffer)
1334 IF (aligned_degenerate_blocks == 0) aligned_degenerate_min_svalue = 0.0_dp
1335 success = .true.
1336
1337 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
1362 TYPE(cp_fm_type), POINTER :: fmi, fmr
1363 TYPE(kpoint_env_type), POINTER :: kp
1364
1365 CALL get_kpoint_info(kpoint, nkp=nkp, kp_range=kp_range)
1366 kp => kpoint%kp_env(1)%kpoint_env
1367 CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo)
1368 ALLOCATE (mo_real(nao, nmo, nkp, nspins), mo_imag(nao, nmo, nkp, nspins), &
1369 eigenvalue_snapshot(nmo, nkp, nspins), owner_weight(nkp, nspins))
1370 mo_real(:, :, :, :) = 0.0_dp
1371 mo_imag(:, :, :, :) = 0.0_dp
1372 eigenvalue_snapshot(:, :, :) = 0.0_dp
1373 owner_weight(:, :) = 0.0_dp
1374 DO ik = kp_range(1), kp_range(2)
1375 ikpgr = ik - kp_range(1) + 1
1376 kp => kpoint%kp_env(ikpgr)%kpoint_env
1377 DO ispin = 1, nspins
1378 fmr => kp%mos(1, ispin)%mo_coeff
1379 fmi => kp%mos(2, ispin)%mo_coeff
1380 CALL cp_fm_get_submatrix(fmr, mo_real(:, :, ik, ispin))
1381 CALL cp_fm_get_submatrix(fmi, mo_imag(:, :, ik, ispin))
1382 CALL get_mo_set(kp%mos(1, ispin), eigenvalues=eigenvalues)
1383 eigenvalue_snapshot(1:nmo, ik, ispin) = eigenvalues(1:nmo)
1384 owner_weight(ik, ispin) = 1.0_dp
1385 END DO
1386 END DO
1387 CALL para_env%sum(mo_real)
1388 CALL para_env%sum(mo_imag)
1389 CALL para_env%sum(eigenvalue_snapshot)
1390 CALL para_env%sum(owner_weight)
1391 DO ik = 1, nkp
1392 DO ispin = 1, nspins
1393 IF (owner_weight(ik, ispin) > 0.0_dp) THEN
1394 mo_real(:, :, ik, ispin) = mo_real(:, :, ik, ispin)/owner_weight(ik, ispin)
1395 mo_imag(:, :, ik, ispin) = mo_imag(:, :, ik, ispin)/owner_weight(ik, ispin)
1396 eigenvalue_snapshot(:, ik, ispin) = &
1397 eigenvalue_snapshot(:, ik, ispin)/owner_weight(ik, ispin)
1398 END IF
1399 END DO
1400 END DO
1401 DEALLOCATE (owner_weight)
1402
1403 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
1420 TYPE(cp_fm_type), POINTER :: fmi, fmr
1421 TYPE(kpoint_env_type), POINTER :: kp
1422
1423 CALL get_kpoint_info(kpoint, kp_range=kp_range)
1424 nmo = SIZE(eigenvalue_snapshot, 1)
1425 nspins = SIZE(eigenvalue_snapshot, 3)
1426 DO ik = kp_range(1), kp_range(2)
1427 ikpgr = ik - kp_range(1) + 1
1428 kp => kpoint%kp_env(ikpgr)%kpoint_env
1429 DO ispin = 1, nspins
1430 fmr => kp%mos(1, ispin)%mo_coeff
1431 fmi => kp%mos(2, ispin)%mo_coeff
1432 CALL cp_fm_set_submatrix(fmr, mo_real(:, :, ik, ispin))
1433 CALL cp_fm_set_submatrix(fmi, mo_imag(:, :, ik, ispin))
1434 CALL get_mo_set(kp%mos(1, ispin), eigenvalues=eigenvalues)
1435 eigenvalues(1:nmo) = eigenvalue_snapshot(1:nmo, ik, ispin)
1436 CALL get_mo_set(kp%mos(2, ispin), eigenvalues=eigenvalues)
1437 IF (ASSOCIATED(eigenvalues)) eigenvalues(1:nmo) = eigenvalue_snapshot(1:nmo, ik, ispin)
1438 END DO
1439 END DO
1440
1441 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 SUBROUTINE validate_wannier90_reused_mos(kpoint, matrix_s, cell_to_index, sab_nl, para_env, &
1459 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 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 success = .false.
1490 max_subspace_deviation = 0.0_dp
1491 min_svalue = huge(1.0_dp)
1492 max_eigenvalue_deviation = 0.0_dp
1493 NULLIFY (matrix_struct_work, fmr, fmi)
1494
1495 CALL get_kpoint_info(kpoint, nkp=nkp, kp_range=kp_range)
1496 kp => kpoint%kp_env(1)%kpoint_env
1497 nspins = SIZE(kp%mos, 2)
1498 CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo)
1499 cpassert(SIZE(reference_mo_real, 1) == nao)
1500 cpassert(SIZE(reference_mo_real, 2) == nmo)
1501 cpassert(SIZE(reference_mo_real, 3) == nkp)
1502 cpassert(SIZE(reference_mo_real, 4) == nspins)
1503
1504 CALL cp_fm_get_info(kp%mos(1, 1)%mo_coeff, matrix_struct=matrix_struct_work)
1505 CALL cp_fm_create(ref_real, matrix_struct_work)
1506 CALL cp_fm_create(ref_imag, matrix_struct_work)
1507 CALL cp_fm_create(cand_real, matrix_struct_work)
1508 CALL cp_fm_create(cand_imag, matrix_struct_work)
1509 ALLOCATE (eigenvalue_buffer(nmo))
1510
1511 DO ispin = 1, nspins
1512 DO ik = 1, nkp
1513 CALL cp_fm_set_submatrix(ref_real, reference_mo_real(:, :, ik, ispin))
1514 CALL cp_fm_set_submatrix(ref_imag, reference_mo_imag(:, :, ik, ispin))
1515 my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
1516 IF (my_kpgrp) THEN
1517 ikpgr = ik - kp_range(1) + 1
1518 kp => kpoint%kp_env(ikpgr)%kpoint_env
1519 fmr => kp%mos(1, ispin)%mo_coeff
1520 fmi => kp%mos(2, ispin)%mo_coeff
1521 CALL get_mo_set(kp%mos(1, ispin), eigenvalues=eigenvalues)
1522 eigenvalue_buffer(1:nmo) = eigenvalues(1:nmo)
1523 ELSE
1524 NULLIFY (fmr, fmi)
1525 eigenvalue_buffer(1:nmo) = 0.0_dp
1526 END IF
1527 CALL cp_fm_copy_general(fmr, cand_real, para_env)
1528 CALL cp_fm_copy_general(fmi, cand_imag, para_env)
1529 IF (my_kpgrp) THEN
1530 owner_count = 1.0_dp
1531 ELSE
1532 owner_count = 0.0_dp
1533 END IF
1534 CALL para_env%sum(owner_count)
1535 CALL para_env%sum(eigenvalue_buffer)
1536 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 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 ok, candidate_deviation, candidate_svalue)
1543 IF (.NOT. ok) THEN
1544 max_subspace_deviation = huge(1.0_dp)
1545 ELSE
1546 max_subspace_deviation = max(max_subspace_deviation, candidate_deviation)
1547 min_svalue = min(min_svalue, candidate_svalue)
1548 END IF
1549 END DO
1550 END DO
1551 CALL para_env%max(max_subspace_deviation)
1552 CALL para_env%min(min_svalue)
1553 CALL para_env%max(max_eigenvalue_deviation)
1554 success = max_subspace_deviation < subspace_tol .AND. max_eigenvalue_deviation < eigenvalue_tol
1555
1556 DEALLOCATE (eigenvalue_buffer)
1557 CALL cp_fm_release(ref_real)
1558 CALL cp_fm_release(ref_imag)
1559 CALL cp_fm_release(cand_real)
1560 CALL cp_fm_release(cand_imag)
1561
1562 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 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_buffer
1607 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 max_subspace_deviation = huge(1.0_dp)
1617 min_svalue = 0.0_dp
1618 max_metric_deviation = huge(1.0_dp)
1619 max_residual = huge(1.0_dp)
1620 NULLIFY (matrix_struct_source, matrix_struct_work, fmi, fmr, src_fmi, src_fmr)
1621
1622 CALL build_wannier90_scf_mapping(kpoint, qs_kpoint, source_kpoint, sym_index, ok, reason)
1623 IF (.NOT. ok) RETURN
1624 kp => kpoint%kp_env(1)%kpoint_env
1625 kp_source => qs_kpoint%kp_env(1)%kpoint_env
1626 IF (SIZE(kp%mos, 1) < 2 .OR. SIZE(kp_source%mos, 1) < 2) THEN
1627 DEALLOCATE (source_kpoint, sym_index)
1628 RETURN
1629 END IF
1630 nspins = SIZE(kp%mos, 2)
1631 CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo)
1632 CALL get_mo_set(kp_source%mos(1, 1), nao=nao_src, nmo=nmo_src)
1633 CALL para_env%max(nao_src)
1634 CALL para_env%max(nmo_src)
1635 IF (nao_src /= nao .OR. nmo_src < nmo) THEN
1636 DEALLOCATE (source_kpoint, sym_index)
1637 RETURN
1638 END IF
1639 source_window = nmo_src > nmo
1640 CALL get_kpoint_info(kpoint, nkp=nkp, kp_range=kp_range)
1641 CALL get_kpoint_info(qs_kpoint, kp_range=source_kp_range)
1642 IF (source_kp_range(1) /= 1 .OR. source_kp_range(2) /= qs_kpoint%nkp) THEN
1643 DEALLOCATE (source_kpoint, sym_index)
1644 RETURN
1645 END IF
1646
1647 CALL cp_fm_struct_create(matrix_struct_work, nrow_global=nao, ncol_global=nmo, &
1648 para_env=para_env, context=kpoint%blacs_env)
1649 CALL cp_fm_create(ref_real, matrix_struct_work)
1650 CALL cp_fm_create(ref_imag, matrix_struct_work)
1651 CALL cp_fm_create(src_real, matrix_struct_work)
1652 CALL cp_fm_create(src_imag, matrix_struct_work)
1653 CALL cp_fm_create(dst_real, matrix_struct_work)
1654 CALL cp_fm_create(dst_imag, matrix_struct_work)
1655 ALLOCATE (eigenvalues_buffer(nmo))
1656 IF (source_window) THEN
1657 CALL cp_fm_struct_create(matrix_struct_source, nrow_global=nao, ncol_global=nmo_src, &
1658 para_env=para_env, context=kpoint%blacs_env)
1659 CALL cp_fm_create(src_real_full, matrix_struct_source)
1660 CALL cp_fm_create(src_imag_full, matrix_struct_source)
1661 END IF
1662
1663 max_subspace_deviation = 0.0_dp
1664 min_svalue = huge(1.0_dp)
1665 max_metric_deviation = 0.0_dp
1666 max_residual = 0.0_dp
1667 DO ispin = 1, nspins
1668 DO ik = 1, nkp
1669 my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
1670 IF (my_kpgrp) THEN
1671 ikpgr = ik - kp_range(1) + 1
1672 kp => kpoint%kp_env(ikpgr)%kpoint_env
1673 fmr => kp%mos(1, ispin)%mo_coeff
1674 fmi => kp%mos(2, ispin)%mo_coeff
1675 ELSE
1676 NULLIFY (fmr, fmi)
1677 END IF
1678 CALL cp_fm_copy_general(fmr, ref_real, para_env)
1679 CALL cp_fm_copy_general(fmi, ref_imag, para_env)
1680
1681 ikred = source_kpoint(ik)
1682 my_kpgrp = (ikred >= source_kp_range(1) .AND. ikred <= source_kp_range(2))
1683 IF (my_kpgrp) THEN
1684 ikpgr = ikred - source_kp_range(1) + 1
1685 kp_source => qs_kpoint%kp_env(ikpgr)%kpoint_env
1686 src_fmr => kp_source%mos(1, ispin)%mo_coeff
1687 src_fmi => kp_source%mos(2, ispin)%mo_coeff
1688 CALL get_mo_set(kp_source%mos(1, ispin), eigenvalues=eigenvalues)
1689 eigenvalues_buffer(1:nmo) = eigenvalues(1:nmo)
1690 ELSE
1691 NULLIFY (src_fmr, src_fmi)
1692 eigenvalues_buffer(1:nmo) = 0.0_dp
1693 END IF
1694 IF (my_kpgrp) THEN
1695 owner_count = 1.0_dp
1696 ELSE
1697 owner_count = 0.0_dp
1698 END IF
1699 CALL para_env%sum(owner_count)
1700 CALL para_env%sum(eigenvalues_buffer)
1701 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 ref_metric_deviation, ref_metric_min, ref_residual)
1706 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 "WANNIER90| reference k=", ik, " dM=", ref_metric_deviation, &
1710 " smin=", ref_metric_min, " resid=", ref_residual
1711 END IF
1712 IF (source_window) THEN
1713 CALL cp_fm_copy_general(src_fmr, src_real_full, para_env)
1714 CALL cp_fm_copy_general(src_fmi, src_imag_full, para_env)
1715 CALL copy_wannier90_mo_window(src_real_full, src_real, nmo)
1716 CALL copy_wannier90_mo_window(src_imag_full, src_imag, nmo)
1717 ELSE
1718 CALL cp_fm_copy_general(src_fmr, src_real, para_env)
1719 CALL cp_fm_copy_general(src_fmi, src_imag, para_env)
1720 END IF
1721
1722 best_deviation = huge(1.0_dp)
1723 best_metric_deviation = 0.0_dp
1724 best_residual = 0.0_dp
1725 best_svalue = 0.0_dp
1726 IF (sym_index(ik) <= 0) THEN
1727 CALL kpoint_transform_scf_mo(src_real, src_imag, dst_real, dst_imag, qs_kpoint, &
1728 ikred, sym_index(ik), para_env, ok, reason)
1729 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 para_env, ok, candidate_deviation, candidate_svalue)
1733 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 candidate_metric_min, candidate_residual)
1739 END IF
1740 IF (ok) THEN
1741 best_deviation = candidate_deviation
1742 best_metric_deviation = candidate_metric_deviation
1743 best_residual = candidate_residual
1744 best_svalue = candidate_svalue
1745 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 "WANNIER90| reuse candidate k=", ik, " src=", ikred, " sym=", &
1750 sym_index(ik), " dRef=", candidate_deviation, " dM=", &
1751 candidate_metric_deviation, " smin=", candidate_metric_min, &
1752 " resid=", candidate_residual
1753 END IF
1754 END IF
1755 END IF
1756 ELSE IF (ASSOCIATED(qs_kpoint%kp_sym)) THEN
1757 kpsym => qs_kpoint%kp_sym(ikred)%kpoint_sym
1758 IF (ASSOCIATED(kpsym)) THEN
1759 DO isym_try = 1, kpsym%nwred
1760 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 qs_kpoint, ikred, isym_try, para_env, ok, reason)
1764 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 candidate_svalue)
1769 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 candidate_metric_min, candidate_residual)
1776 END IF
1777 IF (ok .AND. candidate_deviation < best_deviation) THEN
1778 best_deviation = candidate_deviation
1779 best_metric_deviation = candidate_metric_deviation
1780 best_residual = candidate_residual
1781 best_svalue = candidate_svalue
1782 END IF
1783 END DO
1784 END IF
1785 END IF
1786 IF (best_deviation < huge(1.0_dp)) THEN
1787 max_subspace_deviation = max(max_subspace_deviation, best_deviation)
1788 min_svalue = min(min_svalue, best_svalue)
1789 max_metric_deviation = max(max_metric_deviation, best_metric_deviation)
1790 max_residual = max(max_residual, best_residual)
1791 END IF
1792 END DO
1793 END DO
1794 CALL para_env%max(max_subspace_deviation)
1795 CALL para_env%min(min_svalue)
1796 CALL para_env%max(max_metric_deviation)
1797 CALL para_env%max(max_residual)
1798
1799 IF (source_window) THEN
1800 CALL cp_fm_release(src_real_full)
1801 CALL cp_fm_release(src_imag_full)
1802 CALL cp_fm_struct_release(matrix_struct_source)
1803 END IF
1804 CALL cp_fm_release(ref_real)
1805 CALL cp_fm_release(ref_imag)
1806 CALL cp_fm_release(src_real)
1807 CALL cp_fm_release(src_imag)
1808 CALL cp_fm_release(dst_real)
1809 CALL cp_fm_release(dst_imag)
1810 CALL cp_fm_struct_release(matrix_struct_work)
1811 DEALLOCATE (eigenvalues_buffer)
1812 DEALLOCATE (source_kpoint, sym_index)
1813
1814 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 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 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: metric_projected, metric_vectors, &
1845 overlap, ref_coeff, s_cand
1846 INTEGER :: ib, nao, nmo, nmo_candidate
1847 REAL(kind=dp) :: singular_value
1848 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: metric_values
1849 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 success = .false.
1854 max_subspace_deviation = huge(1.0_dp)
1855 min_svalue = 0.0_dp
1856 NULLIFY (matrix_struct_metric)
1857
1858 CALL cp_fm_get_info(ref_real, nrow_global=nao, ncol_global=nmo, &
1859 matrix_struct=matrix_struct_metric)
1860 CALL cp_fm_get_info(cand_real, ncol_global=nmo_candidate)
1861 IF (nmo_candidate /= nmo) RETURN
1862
1863 CALL cp_fm_create(s_cand_real, matrix_struct_metric)
1864 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 cand_real, cand_imag, s_cand_real, s_cand_imag)
1867
1868 ALLOCATE (ref_r(nao, nmo), ref_i(nao, nmo), s_cand_r(nao, nmo), s_cand_i(nao, nmo))
1869 CALL cp_fm_get_submatrix(ref_real, ref_r)
1870 CALL cp_fm_get_submatrix(ref_imag, ref_i)
1871 CALL cp_fm_get_submatrix(s_cand_real, s_cand_r)
1872 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 metric_projected(nmo, nmo), metric_vectors(nmo, nmo), metric_values(nmo))
1876 ref_coeff(:, :) = cmplx(ref_r, ref_i, kind=dp)
1877 s_cand(:, :) = cmplx(s_cand_r, s_cand_i, kind=dp)
1878 overlap(:, :) = matmul(conjg(transpose(ref_coeff)), s_cand)
1879 metric_projected(:, :) = matmul(conjg(transpose(overlap)), overlap)
1880 metric_projected(:, :) = 0.5_dp*(metric_projected + conjg(transpose(metric_projected)))
1881 CALL diag_complex(metric_projected, metric_vectors, metric_values)
1882
1883 min_svalue = huge(1.0_dp)
1884 max_subspace_deviation = 0.0_dp
1885 DO ib = 1, nmo
1886 singular_value = sqrt(max(metric_values(ib), 0.0_dp))
1887 min_svalue = min(min_svalue, singular_value)
1888 max_subspace_deviation = max(max_subspace_deviation, abs(singular_value - 1.0_dp))
1889 END DO
1890 CALL para_env%max(max_subspace_deviation)
1891 CALL para_env%min(min_svalue)
1892 success = .true.
1893
1894 DEALLOCATE (ref_coeff, s_cand, overlap, metric_projected, metric_vectors, metric_values)
1895 DEALLOCATE (ref_r, ref_i, s_cand_r, s_cand_i)
1896 CALL cp_fm_release(s_cand_real)
1897 CALL cp_fm_release(s_cand_imag)
1898
1899 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 SUBROUTINE measure_wannier90_eigenspace_quality(cand_real, cand_imag, matrix_s, matrix_ks, xkp, &
1919 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 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cand_coeff, h_coeff, metric_vectors, &
1936 residual_block, s_coeff, s_projected
1937 INTEGER :: ib, nao, nmo
1938 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: metric_values
1939 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cand_i, cand_r, h_coeff_i, h_coeff_r, &
1940 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 success = .false.
1948 metric_deviation = huge(1.0_dp)
1949 min_metric_eigenvalue = 0.0_dp
1950 residual_norm = huge(1.0_dp)
1951 NULLIFY (matrix_struct_metric, matrix_struct_projected)
1952
1953 CALL cp_fm_get_info(cand_real, nrow_global=nao, ncol_global=nmo, &
1954 matrix_struct=matrix_struct_metric)
1955 IF (SIZE(eigenvalues) < nmo) RETURN
1956
1957 CALL cp_fm_create(s_cand_real, matrix_struct_metric)
1958 CALL cp_fm_create(s_cand_imag, matrix_struct_metric)
1959 CALL cp_fm_create(h_cand_real, matrix_struct_metric)
1960 CALL cp_fm_create(h_cand_imag, matrix_struct_metric)
1961 CALL cp_fm_create(tmp_fm, matrix_struct_metric)
1962 CALL cp_cfm_create(cand_cfm, matrix_struct_metric)
1963 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 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 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 s_coeff_i(nao, nmo), h_coeff_r(nao, nmo), h_coeff_i(nao, nmo))
1972 CALL cp_fm_get_submatrix(cand_real, cand_r)
1973 CALL cp_fm_get_submatrix(cand_imag, cand_i)
1974 CALL cp_fm_get_submatrix(s_cand_real, s_coeff_r)
1975 CALL cp_fm_get_submatrix(s_cand_imag, s_coeff_i)
1976 CALL cp_fm_get_submatrix(h_cand_real, h_coeff_r)
1977 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 metric_values(nmo))
1982 cand_coeff(:, :) = cmplx(cand_r, cand_i, kind=dp)
1983 s_coeff(:, :) = cmplx(s_coeff_r, s_coeff_i, kind=dp)
1984 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 context=matrix_struct_metric%context)
1988 CALL cp_cfm_create(metric_cfm, matrix_struct_projected)
1989 CALL cp_fm_to_cfm(cand_real, cand_imag, cand_cfm)
1990 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 s_cfm, cmplx(0.0_dp, 0.0_dp, kind=dp), metric_cfm)
1993 CALL cp_cfm_get_submatrix(metric_cfm, s_projected)
1994 s_projected(:, :) = 0.5_dp*(s_projected + conjg(transpose(s_projected)))
1995 CALL diag_complex(s_projected, metric_vectors, metric_values)
1996 metric_deviation = maxval(abs(metric_values - 1.0_dp))
1997 min_metric_eigenvalue = minval(metric_values)
1998
1999 residual_block(:, :) = h_coeff
2000 DO ib = 1, nmo
2001 residual_block(:, ib) = residual_block(:, ib) - eigenvalues(ib)*s_coeff(:, ib)
2002 END DO
2003 residual_norm = maxval(abs(residual_block))
2004 CALL para_env%max(metric_deviation)
2005 CALL para_env%min(min_metric_eigenvalue)
2006 CALL para_env%max(residual_norm)
2007 success = .true.
2008
2009 DEALLOCATE (cand_coeff, h_coeff, metric_vectors, residual_block, s_coeff, s_projected, &
2010 metric_values)
2011 DEALLOCATE (cand_r, cand_i, s_coeff_r, s_coeff_i, h_coeff_r, h_coeff_i)
2012 CALL cp_fm_release(s_cand_real)
2013 CALL cp_fm_release(s_cand_imag)
2014 CALL cp_fm_release(h_cand_real)
2015 CALL cp_fm_release(h_cand_imag)
2016 CALL cp_fm_release(tmp_fm)
2017 CALL cp_cfm_release(cand_cfm)
2018 CALL cp_cfm_release(s_cfm)
2019 CALL cp_cfm_release(metric_cfm)
2020 CALL cp_fm_struct_release(matrix_struct_projected)
2021
2022 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 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 CALL cp_fm_get_info(source, nrow_global=nrow, ncol_global=ncol_source)
2038 cpassert(ncol_source >= ncol)
2039 ALLOCATE (source_buffer(nrow, ncol_source), destination_buffer(nrow, ncol))
2040 CALL cp_fm_get_submatrix(source, source_buffer)
2041 destination_buffer(1:nrow, 1:ncol) = source_buffer(1:nrow, 1:ncol)
2042 CALL cp_fm_set_submatrix(destination, destination_buffer)
2043 DEALLOCATE (source_buffer, destination_buffer)
2044
2045 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 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 NULLIFY (matrix_struct_ao, matrix_struct_coeff, kmat_imag, kmat_imag_full, kmat_real, &
2078 kmat_real_full)
2079
2080 CALL cp_fm_get_info(coeff_real, nrow_global=nao, ncol_global=ncol, &
2081 matrix_struct=matrix_struct_coeff)
2082
2083 ALLOCATE (kmat_real, kmat_imag, kmat_real_full, kmat_imag_full)
2084 CALL dbcsr_create(kmat_real, template=rsmat(ispin, 1)%matrix, &
2085 matrix_type=dbcsr_type_symmetric)
2086 CALL dbcsr_create(kmat_imag, template=rsmat(ispin, 1)%matrix, &
2087 matrix_type=dbcsr_type_antisymmetric)
2088 CALL dbcsr_create(kmat_real_full, template=rsmat(ispin, 1)%matrix, &
2089 matrix_type=dbcsr_type_no_symmetry)
2090 CALL dbcsr_create(kmat_imag_full, template=rsmat(ispin, 1)%matrix, &
2091 matrix_type=dbcsr_type_no_symmetry)
2092 CALL cp_dbcsr_alloc_block_from_nbl(kmat_real, sab_nl)
2093 CALL cp_dbcsr_alloc_block_from_nbl(kmat_imag, sab_nl)
2094 CALL dbcsr_set(kmat_real, 0.0_dp)
2095 CALL dbcsr_set(kmat_imag, 0.0_dp)
2096 CALL rskp_transform(kmat_real, kmat_imag, rsmat=rsmat, ispin=ispin, &
2097 xkp=xkp, cell_to_index=cell_to_index, sab_nl=sab_nl)
2098 CALL dbcsr_desymmetrize(kmat_real, kmat_real_full)
2099 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 context=matrix_struct_coeff%context)
2104 CALL cp_fm_create(mat_real, matrix_struct_ao)
2105 CALL cp_fm_create(mat_imag, matrix_struct_ao)
2106 CALL copy_dbcsr_to_fm(kmat_real_full, mat_real)
2107 CALL copy_dbcsr_to_fm(kmat_imag_full, mat_imag)
2108
2109 CALL cp_cfm_create(kmat_cfm, matrix_struct_ao)
2110 CALL cp_cfm_create(coeff_cfm, matrix_struct_coeff)
2111 CALL cp_cfm_create(result_cfm, matrix_struct_coeff)
2112 CALL cp_fm_to_cfm(mat_real, mat_imag, kmat_cfm)
2113 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 coeff_cfm, cmplx(0.0_dp, 0.0_dp, kind=dp), result_cfm)
2116 CALL cp_cfm_to_fm(result_cfm, result_real, result_imag)
2117
2118 CALL cp_fm_release(mat_real)
2119 CALL cp_fm_release(mat_imag)
2120 CALL cp_cfm_release(kmat_cfm)
2121 CALL cp_cfm_release(coeff_cfm)
2122 CALL cp_cfm_release(result_cfm)
2123 CALL cp_fm_struct_release(matrix_struct_ao)
2124 CALL dbcsr_deallocate_matrix(kmat_real)
2125 CALL dbcsr_deallocate_matrix(kmat_imag)
2126 CALL dbcsr_deallocate_matrix(kmat_real_full)
2127 CALL dbcsr_deallocate_matrix(kmat_imag_full)
2128
2129 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 SUBROUTINE ritz_stabilize_wannier90_subspace(dst_real, dst_imag, matrix_s, matrix_ks, &
2151 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 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: block_coeff, h_block, h_coeff, &
2171 h_projected, h_projected_work, metric_vectors, residual_block, ritz_vectors, s_block, &
2172 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: metric_values, ritz_values
2178 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dst_i, dst_r, h_coeff_i, h_coeff_r, &
2179 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 success = .false.
2185 reason = ""
2186 aligned_blocks = 0
2187 aligned_max_size = 0
2188 aligned_min_svalue = huge(1.0_dp)
2189 max_residual = 0.0_dp
2190
2191 NULLIFY (matrix_struct_metric)
2192 CALL cp_fm_get_info(dst_real, nrow_global=nao, ncol_global=nmo, &
2193 matrix_struct=matrix_struct_metric)
2194 IF (SIZE(eigenvalues) < nmo) THEN
2195 reason = "not enough eigenvalues for Wannier90 Ritz subspace stabilization"
2196 RETURN
2197 END IF
2198
2199 CALL cp_fm_create(s_dst_real, matrix_struct_metric)
2200 CALL cp_fm_create(s_dst_imag, matrix_struct_metric)
2201 CALL cp_fm_create(h_dst_real, matrix_struct_metric)
2202 CALL cp_fm_create(h_dst_imag, matrix_struct_metric)
2203 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 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 dst_real, dst_imag, h_dst_real, h_dst_imag)
2209
2210 ALLOCATE (dst_r(nao, nmo), dst_i(nao, nmo), s_coeff_r(nao, nmo), s_coeff_i(nao, nmo), &
2211 h_coeff_r(nao, nmo), h_coeff_i(nao, nmo))
2212 CALL cp_fm_get_submatrix(dst_real, dst_r)
2213 CALL cp_fm_get_submatrix(dst_imag, dst_i)
2214 CALL cp_fm_get_submatrix(s_dst_real, s_coeff_r)
2215 CALL cp_fm_get_submatrix(s_dst_imag, s_coeff_i)
2216 CALL cp_fm_get_submatrix(h_dst_real, h_coeff_r)
2217 CALL cp_fm_get_submatrix(h_dst_imag, h_coeff_i)
2218
2219 ALLOCATE (s_coeff(nao, nmo), h_coeff(nao, nmo))
2220 s_coeff(:, :) = cmplx(s_coeff_r, s_coeff_i, kind=dp)
2221 h_coeff(:, :) = cmplx(h_coeff_r, h_coeff_i, kind=dp)
2222
2223 block_first = 1
2224 DO WHILE (block_first <= nmo)
2225 block_last = block_first
2226 DO WHILE (block_last < nmo)
2227 IF (abs(eigenvalues(block_last + 1) - eigenvalues(block_last)) >= degenerate_band_tol) EXIT
2228 block_last = block_last + 1
2229 END DO
2230 block_size = block_last - block_first + 1
2231 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 ALLOCATE (block_coeff(nao, block_size), h_block(nao, block_size), &
2236 h_projected(block_size, block_size), h_projected_work(block_size, block_size), &
2237 metric_vectors(block_size, block_size), residual_block(nao, block_size), &
2238 ritz_vectors(block_size, block_size), s_block(nao, block_size), &
2239 s_projected(block_size, block_size), stabilized(nao, block_size), &
2240 metric_values(block_size), ritz_values(block_size))
2241 block_coeff(:, :) = cmplx(dst_r(:, block_first:block_last), &
2242 dst_i(:, block_first:block_last), kind=dp)
2243 s_block(:, :) = s_coeff(:, block_first:block_last)
2244 h_block(:, :) = h_coeff(:, block_first:block_last)
2245 s_projected(:, :) = matmul(conjg(transpose(block_coeff)), s_block)
2246 h_projected(:, :) = matmul(conjg(transpose(block_coeff)), h_block)
2247 s_projected(:, :) = 0.5_dp*(s_projected + conjg(transpose(s_projected)))
2248 h_projected(:, :) = 0.5_dp*(h_projected + conjg(transpose(h_projected)))
2249
2250 CALL diag_complex(s_projected, metric_vectors, metric_values)
2251 aligned_min_svalue = min(aligned_min_svalue, minval(metric_values))
2252 metric_deviation = maxval(abs(metric_values - 1.0_dp))
2253 IF (minval(metric_values) < 1.0e-10_dp) THEN
2254 WRITE (reason, "(A,I0,A,ES9.2,A,ES9.2)") &
2255 "singular metric blk=", block_first, " smin=", minval(metric_values), &
2256 " dS=", metric_deviation
2257 max_residual = huge(1.0_dp)
2258 DEALLOCATE (block_coeff, h_block, h_projected, h_projected_work, metric_vectors, &
2259 residual_block, ritz_vectors, s_block, s_projected, stabilized, &
2260 metric_values, ritz_values)
2261 DEALLOCATE (s_coeff, h_coeff, dst_r, dst_i, s_coeff_r, s_coeff_i, h_coeff_r, &
2262 h_coeff_i)
2263 CALL cp_fm_release(s_dst_real)
2264 CALL cp_fm_release(s_dst_imag)
2265 CALL cp_fm_release(h_dst_real)
2266 CALL cp_fm_release(h_dst_imag)
2267 CALL cp_fm_release(tmp_fm)
2268 RETURN
2269 END IF
2270
2271 DO ib = 1, block_size
2272 metric_vectors(:, ib) = metric_vectors(:, ib)/sqrt(metric_values(ib))
2273 END DO
2274 h_projected_work(:, :) = matmul(h_projected, metric_vectors)
2275 h_projected(:, :) = matmul(conjg(transpose(metric_vectors)), h_projected_work)
2276 h_projected(:, :) = 0.5_dp*(h_projected + conjg(transpose(h_projected)))
2277 CALL diag_complex(h_projected, ritz_vectors, ritz_values)
2278 h_projected_work(:, :) = matmul(metric_vectors, ritz_vectors)
2279 ritz_vectors(:, :) = h_projected_work
2280 stabilized(:, :) = matmul(block_coeff, ritz_vectors)
2281 residual_block(:, :) = matmul(h_block, ritz_vectors)
2282 h_block(:, :) = residual_block
2283 residual_block(:, :) = matmul(s_block, ritz_vectors)
2284 s_block(:, :) = residual_block
2285 DO ib = 1, block_size
2286 norm_value = sqrt(abs(real(dot_product(stabilized(:, ib), s_block(:, ib)), kind=dp)))
2287 IF (norm_value > epsilon(1.0_dp)) THEN
2288 stabilized(:, ib) = stabilized(:, ib)/norm_value
2289 h_block(:, ib) = h_block(:, ib)/norm_value
2290 s_block(:, ib) = s_block(:, ib)/norm_value
2291 END IF
2292 END DO
2293 residual_block(:, :) = h_block
2294 DO ib = 1, block_size
2295 residual_block(:, ib) = residual_block(:, ib) - &
2296 eigenvalues(block_first + ib - 1)*s_block(:, ib)
2297 END DO
2298 residual_norm = maxval(abs(residual_block))
2299 max_residual = max(max_residual, residual_norm)
2300 IF (residual_norm > residual_tol) THEN
2301 WRITE (reason, "(A,I0,A,ES9.2)") &
2302 "blk=", block_first, " dS=", metric_deviation
2303 DEALLOCATE (block_coeff, h_block, h_projected, h_projected_work, metric_vectors, &
2304 residual_block, ritz_vectors, s_block, s_projected, stabilized, &
2305 metric_values, ritz_values)
2306 DEALLOCATE (s_coeff, h_coeff, dst_r, dst_i, s_coeff_r, s_coeff_i, h_coeff_r, &
2307 h_coeff_i)
2308 CALL cp_fm_release(s_dst_real)
2309 CALL cp_fm_release(s_dst_imag)
2310 CALL cp_fm_release(h_dst_real)
2311 CALL cp_fm_release(h_dst_imag)
2312 CALL cp_fm_release(tmp_fm)
2313 RETURN
2314 END IF
2315
2316 dst_r(:, block_first:block_last) = real(stabilized, kind=dp)
2317 dst_i(:, block_first:block_last) = aimag(stabilized)
2318 h_coeff(:, block_first:block_last) = h_block
2319 s_coeff(:, block_first:block_last) = s_block
2320 aligned_blocks = aligned_blocks + 1
2321 aligned_max_size = max(aligned_max_size, block_size)
2322 DEALLOCATE (block_coeff, h_block, h_projected, h_projected_work, metric_vectors, &
2323 residual_block, ritz_vectors, s_block, s_projected, stabilized, metric_values, &
2324 ritz_values)
2325 END IF
2326 block_first = block_last + 1
2327 END DO
2328
2329 DO ib = 1, nmo
2330 residual_norm = maxval(abs(h_coeff(:, ib) - eigenvalues(ib)*s_coeff(:, ib)))
2331 max_residual = max(max_residual, residual_norm)
2332 END DO
2333 IF (max_residual > residual_tol) THEN
2334 WRITE (reason, "(A,ES10.3)") &
2335 "atom/AO W90 reuse guarded: Ritz residual=", max_residual
2336 DEALLOCATE (s_coeff, h_coeff, dst_r, dst_i, s_coeff_r, s_coeff_i, h_coeff_r, h_coeff_i)
2337 CALL cp_fm_release(s_dst_real)
2338 CALL cp_fm_release(s_dst_imag)
2339 CALL cp_fm_release(h_dst_real)
2340 CALL cp_fm_release(h_dst_imag)
2341 CALL cp_fm_release(tmp_fm)
2342 RETURN
2343 END IF
2344
2345 CALL cp_fm_set_submatrix(dst_real, dst_r)
2346 CALL cp_fm_set_submatrix(dst_imag, dst_i)
2347
2348 IF (aligned_blocks == 0) aligned_min_svalue = 0.0_dp
2349 success = .true.
2350
2351 DEALLOCATE (s_coeff, h_coeff, dst_r, dst_i, s_coeff_r, s_coeff_i, h_coeff_r, h_coeff_i)
2352 CALL cp_fm_release(s_dst_real)
2353 CALL cp_fm_release(s_dst_imag)
2354 CALL cp_fm_release(h_dst_real)
2355 CALL cp_fm_release(h_dst_imag)
2356 CALL cp_fm_release(tmp_fm)
2357
2358 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 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 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 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coeff_work, h_coeff, h_projected, &
2400 h_projected_work, metric_vectors, residual_block, ritz_vectors, s_coeff, s_projected, &
2401 source_coeff, stabilized
2402 INTEGER :: ib, nao, nmo_source
2403 REAL(kind=dp) :: max_eigenvalue_shift, metric_deviation, &
2404 norm_value
2405 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: metric_values, ritz_values, &
2406 source_eigenvalues
2407 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dst_i, dst_r, h_coeff_i, h_coeff_r, &
2408 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 success = .false.
2414 reason = ""
2415 min_svalue = huge(1.0_dp)
2416 max_residual = huge(1.0_dp)
2417
2418 NULLIFY (matrix_struct_metric)
2419 CALL cp_fm_get_info(src_real, nrow_global=nao, ncol_global=nmo_source, &
2420 matrix_struct=matrix_struct_metric)
2421 IF (nmo_export > nmo_source) THEN
2422 reason = "Wannier90 export window is larger than the transformed SCF MO space"
2423 RETURN
2424 END IF
2425 IF (SIZE(eigenvalues) < nmo_export) THEN
2426 reason = "not enough eigenvalue storage for Wannier90 source-window reconstruction"
2427 RETURN
2428 END IF
2429
2430 CALL cp_fm_create(s_src_real, matrix_struct_metric)
2431 CALL cp_fm_create(s_src_imag, matrix_struct_metric)
2432 CALL cp_fm_create(h_src_real, matrix_struct_metric)
2433 CALL cp_fm_create(h_src_imag, matrix_struct_metric)
2434 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 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 src_real, src_imag, h_src_real, h_src_imag)
2440
2441 ALLOCATE (src_r(nao, nmo_source), src_i(nao, nmo_source), &
2442 s_coeff_r(nao, nmo_source), s_coeff_i(nao, nmo_source), &
2443 h_coeff_r(nao, nmo_source), h_coeff_i(nao, nmo_source), &
2444 dst_r(nao, nmo_export), dst_i(nao, nmo_export))
2445 CALL cp_fm_get_submatrix(src_real, src_r)
2446 CALL cp_fm_get_submatrix(src_imag, src_i)
2447 CALL cp_fm_get_submatrix(s_src_real, s_coeff_r)
2448 CALL cp_fm_get_submatrix(s_src_imag, s_coeff_i)
2449 CALL cp_fm_get_submatrix(h_src_real, h_coeff_r)
2450 CALL cp_fm_get_submatrix(h_src_imag, h_coeff_i)
2451
2452 ALLOCATE (source_coeff(nao, nmo_source), s_coeff(nao, nmo_source), &
2453 h_coeff(nao, nmo_source), h_projected(nmo_source, nmo_source), &
2454 h_projected_work(nmo_source, nmo_source), metric_vectors(nmo_source, nmo_source), &
2455 residual_block(nao, nmo_export), ritz_vectors(nmo_source, nmo_source), &
2456 s_projected(nmo_source, nmo_source), stabilized(nao, nmo_export), &
2457 coeff_work(nao, nmo_source), metric_values(nmo_source), ritz_values(nmo_source), &
2458 source_eigenvalues(nmo_export))
2459 source_eigenvalues(1:nmo_export) = eigenvalues(1:nmo_export)
2460 source_coeff(:, :) = cmplx(src_r, src_i, kind=dp)
2461 s_coeff(:, :) = cmplx(s_coeff_r, s_coeff_i, kind=dp)
2462 h_coeff(:, :) = cmplx(h_coeff_r, h_coeff_i, kind=dp)
2463 s_projected(:, :) = matmul(conjg(transpose(source_coeff)), s_coeff)
2464 h_projected(:, :) = matmul(conjg(transpose(source_coeff)), h_coeff)
2465 s_projected(:, :) = 0.5_dp*(s_projected + conjg(transpose(s_projected)))
2466 h_projected(:, :) = 0.5_dp*(h_projected + conjg(transpose(h_projected)))
2467
2468 reconstruct_window: block
2469 CALL diag_complex(s_projected, metric_vectors, metric_values)
2470 min_svalue = minval(metric_values)
2471 metric_deviation = maxval(abs(metric_values - 1.0_dp))
2472 IF (min_svalue < 1.0e-10_dp) THEN
2473 WRITE (reason, "(A,ES9.2,A,ES9.2)") &
2474 "singular expanded metric smin=", min_svalue, " dS=", metric_deviation
2475 EXIT reconstruct_window
2476 END IF
2477
2478 DO ib = 1, nmo_source
2479 metric_vectors(:, ib) = metric_vectors(:, ib)/sqrt(metric_values(ib))
2480 END DO
2481 h_projected_work(:, :) = matmul(h_projected, metric_vectors)
2482 h_projected(:, :) = matmul(conjg(transpose(metric_vectors)), h_projected_work)
2483 h_projected(:, :) = 0.5_dp*(h_projected + conjg(transpose(h_projected)))
2484 CALL diag_complex(h_projected, ritz_vectors, ritz_values)
2485 h_projected_work(:, :) = matmul(metric_vectors, ritz_vectors)
2486 ritz_vectors(:, :) = h_projected_work
2487 stabilized(:, :) = matmul(source_coeff, ritz_vectors(:, 1:nmo_export))
2488 coeff_work(:, :) = matmul(h_coeff, ritz_vectors)
2489 h_coeff(:, :) = coeff_work
2490 coeff_work(:, :) = matmul(s_coeff, ritz_vectors)
2491 s_coeff(:, :) = coeff_work
2492 DO ib = 1, nmo_export
2493 norm_value = sqrt(abs(real(dot_product(stabilized(:, ib), s_coeff(:, ib)), kind=dp)))
2494 IF (norm_value > epsilon(1.0_dp)) THEN
2495 stabilized(:, ib) = stabilized(:, ib)/norm_value
2496 h_coeff(:, ib) = h_coeff(:, ib)/norm_value
2497 s_coeff(:, ib) = s_coeff(:, ib)/norm_value
2498 END IF
2499 END DO
2500 residual_block(:, :) = h_coeff(:, 1:nmo_export)
2501 DO ib = 1, nmo_export
2502 residual_block(:, ib) = residual_block(:, ib) - ritz_values(ib)*s_coeff(:, ib)
2503 END DO
2504 max_residual = maxval(abs(residual_block))
2505 IF (max_residual > residual_tol) THEN
2506 WRITE (reason, "(A,ES9.2)") &
2507 "expanded dS=", metric_deviation
2508 EXIT reconstruct_window
2509 END IF
2510 max_eigenvalue_shift = maxval(abs(ritz_values(1:nmo_export) - source_eigenvalues(1:nmo_export)))
2511 IF (max_eigenvalue_shift > eigenvalue_tol) THEN
2512 WRITE (reason, "(A,ES9.2)") &
2513 "expanded dS=", metric_deviation
2514 EXIT reconstruct_window
2515 END IF
2516
2517 dst_r(:, :) = real(stabilized, kind=dp)
2518 dst_i(:, :) = aimag(stabilized)
2519 CALL cp_fm_set_submatrix(dst_real, dst_r)
2520 CALL cp_fm_set_submatrix(dst_imag, dst_i)
2521 success = .true.
2522
2523 END BLOCK reconstruct_window
2524
2525 DEALLOCATE (source_coeff, s_coeff, h_coeff, h_projected, h_projected_work, metric_vectors, &
2526 residual_block, ritz_vectors, s_projected, stabilized, coeff_work, metric_values, &
2527 ritz_values, source_eigenvalues)
2528 DEALLOCATE (src_r, src_i, s_coeff_r, s_coeff_i, h_coeff_r, h_coeff_i, dst_r, dst_i)
2529 CALL cp_fm_release(s_src_real)
2530 CALL cp_fm_release(s_src_imag)
2531 CALL cp_fm_release(h_src_real)
2532 CALL cp_fm_release(h_src_imag)
2533 CALL cp_fm_release(tmp_fm)
2534
2535 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 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 success = .false.
2557 reason = ""
2558 nfull = kpoint%nkp
2559 ALLOCATE (source_kpoint(nfull), sym_index(nfull))
2560 source_kpoint(:) = 0
2561 sym_index(:) = 0
2562
2563 DO ikred = 1, qs_kpoint%nkp
2564 imatch = find_matching_kpoint(kpoint%xkp, qs_kpoint%xkp(1:3, ikred))
2565 IF (imatch > 0 .AND. source_kpoint(imatch) == 0) THEN
2566 source_kpoint(imatch) = ikred
2567 sym_index(imatch) = 0
2568 END IF
2569 END DO
2570
2571 ! Prefer pure time-reversal partners before general atom/AO symmetry operations.
2572 DO ikred = 1, qs_kpoint%nkp
2573 imatch = find_matching_kpoint(kpoint%xkp, -qs_kpoint%xkp(1:3, ikred))
2574 IF (imatch > 0 .AND. source_kpoint(imatch) == 0) THEN
2575 source_kpoint(imatch) = ikred
2576 sym_index(imatch) = -1
2577 END IF
2578 END DO
2579
2580 IF (ASSOCIATED(qs_kpoint%kp_sym)) THEN
2581 DO ikred = 1, qs_kpoint%nkp
2582 kpsym => qs_kpoint%kp_sym(ikred)%kpoint_sym
2583 IF (.NOT. ASSOCIATED(kpsym)) cycle
2584 IF (.NOT. kpsym%apply_symmetry) cycle
2585 DO isym = 1, kpsym%nwred
2586 imatch = find_matching_kpoint(kpoint%xkp, kpsym%xkp(1:3, isym))
2587 IF (imatch > 0 .AND. source_kpoint(imatch) == 0) THEN
2588 source_kpoint(imatch) = ikred
2589 sym_index(imatch) = isym
2590 END IF
2591 END DO
2592 END DO
2593 END IF
2594
2595 DO ik = 1, nfull
2596 IF (source_kpoint(ik) == 0) THEN
2597 reason = "not all full-mesh k-points are represented by the SCF symmetry orbits"
2598 RETURN
2599 END IF
2600 END DO
2601 success = .true.
2602
2603 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 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 ik_match = 0
2618 DO ik = 1, SIZE(xkp_mesh, 2)
2619 IF (kpoint_same_periodic(xkp_mesh(1:3, ik), xkp_search)) THEN
2620 ik_match = ik
2621 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 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 LOGICAL, ALLOCATABLE, DIMENSION(:) :: seen
2642 REAL(kind=dp) :: coord
2643 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: unique_coord
2644
2645 num_kpts = SIZE(kpt_latt, 2)
2646 mp_grid(:) = 0
2647 ALLOCATE (unique_coord(3, num_kpts))
2648 DO idim = 1, 3
2649 n_unique = 0
2650 DO i = 1, num_kpts
2651 coord = kpt_latt(idim, i) - floor(kpt_latt(idim, i))
2652 IF (abs(coord - 1.0_dp) < 1.0e-8_dp) coord = 0.0_dp
2653 known = .false.
2654 DO unique_id = 1, n_unique
2655 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 IF (.NOT. known) THEN
2661 n_unique = n_unique + 1
2662 unique_coord(idim, n_unique) = coord
2663 END IF
2664 END DO
2665 mp_grid(idim) = n_unique
2666 END DO
2667 valid = (mp_grid(1)*mp_grid(2)*mp_grid(3) == num_kpts)
2668 IF (valid) THEN
2669 ALLOCATE (seen(num_kpts))
2670 seen(:) = .false.
2671 DO i = 1, num_kpts
2672 idx = 1
2673 stride = 1
2674 DO idim = 1, 3
2675 coord = kpt_latt(idim, i) - floor(kpt_latt(idim, i))
2676 IF (abs(coord - 1.0_dp) < 1.0e-8_dp) coord = 0.0_dp
2677 coord_id = 0
2678 DO unique_id = 1, mp_grid(idim)
2679 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 cpassert(coord_id > 0)
2685 idx = idx + (coord_id - 1)*stride
2686 stride = stride*mp_grid(idim)
2687 END DO
2688 IF (seen(idx)) valid = .false.
2689 seen(idx) = .true.
2690 END DO
2691 valid = valid .AND. all(seen)
2692 DEALLOCATE (seen)
2693 END IF
2694 DEALLOCATE (unique_coord)
2695
2696 END SUBROUTINE infer_wannier_mp_grid
2697
2698END MODULE qs_wannier90
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition grid_common.h:161
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:233
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_gemm
subroutine, public cp_cfm_gemm(transa, transb, m, n, k, alpha, matrix_a, matrix_b, beta, matrix_c, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, c_first_row)
Performs one of the matrix-matrix operations: matrix_c = alpha * op1( matrix_a ) * op2( matrix_b ) + ...
Definition cp_cfm_basic_linalg.F:559
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_get_submatrix
subroutine, public cp_cfm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
Extract a sub-matrix from the full matrix: op(target_m)(1:n_rows,1:n_cols) = fm(start_row:start_row+n...
Definition cp_cfm_types.F:356
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:182
cp_cfm_types::cp_fm_to_cfm
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
Definition cp_cfm_types.F:844
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:125
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition cp_cfm_types.F:792
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:525
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:212
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:311
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:122
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_copy_general
subroutine, public cp_fm_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
Definition cp_fm_types.F:1524
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_get_element
subroutine, public cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
returns an element of a fm this value is valid on every cpu using this call is expensive
Definition cp_fm_types.F:627
cp_fm_types::cp_fm_set_submatrix
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition cp_fm_types.F:750
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:925
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:708
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kpoint_methods
Routines needed for kpoint calculation.
Definition kpoint_methods.F:15
kpoint_methods::kpoint_initialize_mo_set
subroutine, public kpoint_initialize_mo_set(kpoint)
...
Definition kpoint_methods.F:965
kpoint_methods::rskp_transform
subroutine, public rskp_transform(rmatrix, cmatrix, rsmat, ispin, xkp, cell_to_index, sab_nl, is_complex, rs_sign)
Transformation of real space matrices to a kpoint.
Definition kpoint_methods.F:1139
kpoint_methods::kpoint_init_cell_index
subroutine, public kpoint_init_cell_index(kpoint, sab_nl, para_env, nimages)
Generates the mapping of cell indices and linear RS index CELL (0,0,0) is always mapped to index 1.
Definition kpoint_methods.F:1005
kpoint_methods::kpoint_initialize_mos
subroutine, public kpoint_initialize_mos(kpoint, mos, added_mos, for_aux_fit)
Initialize a set of MOs and density matrix for each kpoint (kpoint group)
Definition kpoint_methods.F:812
kpoint_methods::kpoint_initialize
subroutine, public kpoint_initialize(kpoint, particle_set, cell)
Generate the kpoints and initialize the kpoint environment.
Definition kpoint_methods.F:123
kpoint_methods::kpoint_env_initialize
subroutine, public kpoint_env_initialize(kpoint, para_env, blacs_env, with_aux_fit)
Initialize the kpoint environment.
Definition kpoint_methods.F:624
kpoint_mo_symmetry_methods
Utilities for transforming k-point MO coefficients under symmetry operations.
Definition kpoint_mo_symmetry_methods.F:11
kpoint_mo_symmetry_methods::kpoint_same_periodic
logical function, public kpoint_same_periodic(xkp_a, xkp_b)
Compare two fractional k-points modulo reciprocal lattice vectors.
Definition kpoint_mo_symmetry_methods.F:42
kpoint_mo_symmetry_methods::kpoint_transform_scf_mo
subroutine, public kpoint_transform_scf_mo(src_real, src_imag, dst_real, dst_imag, qs_kpoint, ikred, isym, para_env, success, reason)
Transform one SCF MO coefficient matrix to an equivalent full-mesh k-point.
Definition kpoint_mo_symmetry_methods.F:67
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:529
kpoint_types::kpoint_release
subroutine, public kpoint_release(kpoint)
Release a kpoint environment, deallocate all data.
Definition kpoint_types.F:288
kpoint_types::kpoint_create
subroutine, public kpoint_create(kpoint)
Create a kpoint environment.
Definition kpoint_types.F:237
kpsym
K-points and crystal symmetry routines based on.
Definition kpsym.F:28
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::timestamp_length
integer, parameter, public timestamp_length
Definition machine.F:46
machine::m_timestamp
subroutine, public m_timestamp(timestamp)
Returns a human readable timestamp.
Definition machine.F:381
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::diag_complex
subroutine, public diag_complex(matrix, eigenvectors, eigenvalues)
Diagonalizes a local complex Hermitian matrix using LAPACK. Based on cp_cfm_heevd.
Definition mathlib.F:1878
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:534
qs_environment_types::qs_env_release
subroutine, public qs_env_release(qs_env)
releases the given qs_env (see doc/ReferenceCounting.html)
Definition qs_environment_types.F:1467
qs_gamma2kp
Initialize a qs_env for kpoint calculations starting from a gamma point qs_env.
Definition qs_gamma2kp.F:14
qs_gamma2kp::create_kp_from_gamma
subroutine, public create_kp_from_gamma(qs_env, qs_env_kp, with_xc_terms)
...
Definition qs_gamma2kp.F:63
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:401
qs_moments
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition qs_moments.F:14
qs_moments::build_berry_kpoint_matrix
subroutine, public build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type)
...
Definition qs_moments.F:1527
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_scf_diagonalization
Different diagonalization schemes that can be used for the iterative solution of the eigenvalue probl...
Definition qs_scf_diagonalization.F:15
qs_scf_diagonalization::do_general_diag_kp
subroutine, public do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, update_p, diis_step, diis_error, qs_env, probe)
Kpoint diagonalization routine Transforms matrices to kpoint, distributes kpoint groups,...
Definition qs_scf_diagonalization.F:446
qs_scf_types
module that contains the definitions of the scf types
Definition qs_scf_types.F:14
qs_wannier90
Interface to Wannier90 code.
Definition qs_wannier90.F:14
qs_wannier90::prepare_wannier90_scf_mos
subroutine, public prepare_wannier90_scf_mos(kpoint, qs_kpoint, matrix_s, matrix_ks, cell_to_index, sab_nl, para_env, success, reason, aligned_degenerate_blocks, aligned_degenerate_max_size, aligned_degenerate_min_svalue)
Reconstruct a full Wannier90 k-point MO set from the SCF k-point MOs.
Definition qs_wannier90.F:953
qs_wannier90::wannier90_interface
subroutine, public wannier90_interface(input, logger, qs_env)
...
Definition qs_wannier90.F:114
scf_control_types
parameters that control an scf iteration
Definition scf_control_types.F:17
wannier90
Outtakes from Wannier90 code.
Definition wannier90.F:66
wannier90::wannier_setup
subroutine, public wannier_setup(mp_grid_loc, num_kpts_loc, real_lattice_loc, recip_lattice_loc, kpt_latt_loc, nntot_loc, nnlist_loc, nncell_loc, iounit)
...
Definition wannier90.F:140
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:69
cp_control_types::dft_control_type
Definition cp_control_types.F:757
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_env_type
Keeps information about a specific k-point.
Definition kpoint_types.F:83
kpoint_types::kpoint_sym_type
Keeps symmetry information about a specific k-point.
Definition kpoint_types.F:123
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:175
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:750
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:97
scf_control_types::scf_control_type
Definition scf_control_types.F:126