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