Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Interface to Wannier90 code
10 : !> \par History
11 : !> 06.2016 created [JGH]
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE qs_wannier90
15 : USE atomic_kind_types, ONLY: get_atomic_kind
16 : USE basis_set_types, ONLY: get_gto_basis_set,&
17 : gto_basis_set_type
18 : USE bibliography, ONLY: Gresch2017,&
19 : Soluyanov2011,&
20 : cite_reference
21 : USE cell_types, ONLY: cell_type,&
22 : get_cell
23 : USE cp_blacs_env, ONLY: cp_blacs_env_type
24 : USE cp_cfm_types, ONLY: &
25 : cp_cfm_copy_general, cp_cfm_create, cp_cfm_get_element, cp_cfm_get_info, &
26 : cp_cfm_get_submatrix, cp_cfm_release, cp_cfm_set_submatrix, cp_cfm_to_cfm, cp_cfm_type, &
27 : cp_fm_to_cfm
28 : USE cp_control_types, ONLY: dft_control_type
29 : USE cp_dbcsr_api, ONLY: &
30 : dbcsr_add, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_get_block_p, &
31 : dbcsr_p_type, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
32 : dbcsr_type_no_symmetry, dbcsr_type_symmetric
33 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
34 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
35 : dbcsr_deallocate_matrix_set
36 : USE cp_files, ONLY: close_file,&
37 : open_file
38 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
39 : cp_fm_struct_release,&
40 : cp_fm_struct_type
41 : USE cp_fm_types, ONLY: cp_fm_create,&
42 : cp_fm_release,&
43 : cp_fm_type
44 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit,&
45 : cp_logger_type
46 : USE input_section_types, ONLY: section_vals_get,&
47 : section_vals_get_subs_vals,&
48 : section_vals_type,&
49 : section_vals_val_get
50 : USE kinds, ONLY: default_path_length,&
51 : default_string_length,&
52 : dp
53 : USE kpoint_methods, ONLY: kpoint_env_initialize,&
54 : kpoint_init_cell_index,&
55 : kpoint_initialize,&
56 : kpoint_initialize_mo_set,&
57 : kpoint_initialize_mos,&
58 : rskp_transform
59 : USE kpoint_mo_symmetry_methods, ONLY: kpoint_same_periodic,&
60 : kpoint_transform_scf_mo
61 : USE kpoint_types, ONLY: get_kpoint_info,&
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 parallel_gemm_api, ONLY: parallel_gemm
73 : USE particle_types, ONLY: particle_type
74 : USE physcon, ONLY: angstrom,&
75 : evolt
76 : USE qs_environment_types, ONLY: get_qs_env,&
77 : qs_env_release,&
78 : qs_environment_type
79 : USE qs_gamma2kp, ONLY: create_kp_from_gamma
80 : USE qs_kind_types, ONLY: get_qs_kind,&
81 : qs_kind_type
82 : USE qs_mo_types, ONLY: get_mo_set,&
83 : mo_set_type
84 : USE qs_moments, ONLY: build_berry_kpoint_matrix
85 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
86 : get_neighbor_list_set_p,&
87 : neighbor_list_iterate,&
88 : neighbor_list_iterator_create,&
89 : neighbor_list_iterator_p_type,&
90 : neighbor_list_iterator_release,&
91 : neighbor_list_set_p_type
92 : USE qs_scf_diagonalization, ONLY: do_general_diag_kp
93 : USE qs_scf_types, ONLY: qs_scf_env_type
94 : USE scf_control_types, ONLY: scf_control_type
95 : USE soc_pseudopotential_methods, ONLY: V_SOC_xyz_from_pseudopotential
96 : USE topology_inversion, ONLY: gaussian_inversion_action
97 : USE topology_state_io, ONLY: topology_state_begin,&
98 : topology_state_point
99 : USE topology_symmetry, ONLY: inversion_representation
100 : USE topology_tqc, ONLY: inversion_ebr_signature,&
101 : inversion_indicators
102 : USE topology_wilson, ONLY: chern_from_wcc,&
103 : surface_resolved,&
104 : wcc_distance,&
105 : wilson_spectrum,&
106 : wilson_step,&
107 : z2_from_wcc
108 : USE wannier90, ONLY: wannier_setup
109 : USE wannier90_nnkpts, ONLY: read_wannier90_nnkpts
110 : USE wannier90_projections, ONLY: wannier90_hybrid_candidates,&
111 : wannier90_select_projections,&
112 : wannier90_window_mask
113 : #include "./base/base_uses.f90"
114 :
115 : IMPLICIT NONE
116 : PRIVATE
117 :
118 : CHARACTER(LEN=12), PARAMETER :: window_keywords(4) = &
119 : [CHARACTER(LEN=12) :: 'DIS_WIN_MIN', 'DIS_WIN_MAX', 'DIS_FROZ_MIN', 'DIS_FROZ_MAX']
120 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_wannier90'
121 : INTEGER, PARAMETER, PRIVATE :: w90_kpoints_mp_grid = 0, &
122 : w90_kpoints_scf = 1, w90_kpoints_nnkp = 2, w90_kpoints_wilson = 3, &
123 : w90_kpoints_trim = 4
124 :
125 : TYPE berry_matrix_type
126 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: sinmat => NULL(), cosmat => NULL()
127 : END TYPE berry_matrix_type
128 :
129 : PUBLIC :: wannier90_interface, prepare_wannier90_scf_mos
130 :
131 : ! **************************************************************************************************
132 :
133 : CONTAINS
134 :
135 : ! **************************************************************************************************
136 : !> \brief ...
137 : !> \param input ...
138 : !> \param logger ...
139 : !> \param qs_env ...
140 : ! **************************************************************************************************
141 12967 : SUBROUTINE wannier90_interface(input, logger, qs_env)
142 : TYPE(section_vals_type), POINTER :: input
143 : TYPE(cp_logger_type), POINTER :: logger
144 : TYPE(qs_environment_type), POINTER :: qs_env
145 :
146 : CHARACTER(len=*), PARAMETER :: routineN = 'wannier90_interface'
147 :
148 : INTEGER :: first_spin, handle, i, ichern, ispin, &
149 : iw, iz2, last_spin, max_refinement, &
150 : old_chern, old_z2, refinement, &
151 : selection, source
152 : LOGICAL :: converged, explicit, require_chern, &
153 : require_z2
154 : REAL(KIND=dp) :: error, movement, tolerance
155 12967 : REAL(KIND=dp), ALLOCATABLE :: centres(:, :), previous(:, :)
156 : TYPE(dft_control_type), POINTER :: dft_control
157 : TYPE(section_vals_type), POINTER :: w_input
158 :
159 : !--------------------------------------------------------------------------------------------!
160 :
161 12967 : CALL timeset(routineN, handle)
162 : w_input => section_vals_get_subs_vals(section_vals=input, &
163 12967 : subsection_name="DFT%PRINT%WANNIER90")
164 12967 : CALL section_vals_get(w_input, explicit=explicit)
165 12967 : IF (explicit) THEN
166 :
167 58 : iw = cp_logger_get_default_io_unit(logger)
168 :
169 58 : IF (iw > 0) THEN
170 : WRITE (iw, '(/,T2,A)') &
171 29 : '!-----------------------------------------------------------------------------!'
172 29 : WRITE (iw, '(T32,A)') "Interface to Wannier90"
173 : WRITE (iw, '(T2,A)') &
174 29 : '!-----------------------------------------------------------------------------!'
175 : END IF
176 :
177 58 : CALL get_qs_env(qs_env, dft_control=dft_control)
178 58 : first_spin = 1
179 58 : last_spin = dft_control%nspins
180 58 : IF (last_spin /= 1 .AND. last_spin /= 2) CPABORT('Wannier90 requires collinear spin channels.')
181 58 : CALL section_vals_val_get(w_input, 'SPIN_CHANNEL', i_val=selection)
182 58 : IF (selection /= 0) THEN
183 4 : IF (selection > last_spin) CPABORT('The requested Wannier90 spin channel is not present.')
184 : first_spin = selection
185 : last_spin = selection
186 : END IF
187 118 : DO ispin = first_spin, last_spin
188 60 : CALL section_vals_val_get(w_input, "KPOINTS_SOURCE", i_val=source)
189 60 : CALL section_vals_val_get(w_input, "WILSON_MAX_REFINEMENT", i_val=max_refinement)
190 60 : CALL section_vals_val_get(w_input, "WILSON_TOL", r_val=tolerance)
191 60 : CALL section_vals_val_get(w_input, "Z2", l_val=require_z2)
192 60 : CALL section_vals_val_get(w_input, "CHERN", l_val=require_chern)
193 60 : IF (tolerance <= 0.0_dp) CPABORT("WILSON_TOL must be positive.")
194 60 : IF (source == w90_kpoints_wilson) THEN
195 6 : IF (max_refinement < 1 .OR. max_refinement > 10) THEN
196 0 : CPABORT("WILSON_MAX_REFINEMENT must be between 1 and 10.")
197 : END IF
198 6 : converged = .FALSE.
199 6 : old_z2 = -1
200 6 : old_chern = HUGE(0)
201 12 : DO refinement = 0, max_refinement
202 12 : CALL wannier90_files(qs_env, w_input, iw, ispin, refinement, centres, iz2, ichern)
203 12 : IF (refinement > 0) THEN
204 6 : error = 0.0_dp
205 24 : DO i = 1, SIZE(previous, 2)
206 24 : error = MAX(error, wcc_distance(previous(:, i), centres(:, 2*i - 1)))
207 : END DO
208 6 : movement = 0.0_dp
209 30 : DO i = 2, SIZE(centres, 2)
210 30 : movement = MAX(movement, wcc_distance(centres(:, i - 1), centres(:, i)))
211 : END DO
212 6 : IF (iw > 0) WRITE (iw, '(T2,A,I0,A,ES12.4,A,ES12.4)') &
213 3 : "TOPOLOGY| Refinement ", refinement, ": WCC change ", error, ", transverse step ", movement
214 6 : converged = error < tolerance .AND. movement < 0.1_dp .AND. iz2 == old_z2
215 6 : IF (require_z2) converged = converged .AND. iz2 >= 0 .AND. surface_resolved(centres)
216 6 : IF (require_chern) converged = converged .AND. ichern /= HUGE(0) .AND. ichern == old_chern
217 4 : IF (converged) EXIT
218 : END IF
219 6 : CALL MOVE_ALLOC(centres, previous)
220 6 : old_z2 = iz2
221 10 : old_chern = ichern
222 : END DO
223 6 : IF (.NOT. converged) THEN
224 0 : CPABORT("Wilson surface unconverged; increase WILSON_MAX_REFINEMENT or mesh.")
225 : END IF
226 6 : IF (iw > 0) THEN
227 3 : WRITE (iw, '(T2,A)') "TOPOLOGY| Wilson surface sampling converged."
228 3 : IF (iz2 >= 0) WRITE (iw, '(T2,A,I0)') "TOPOLOGY| Converged Z2 invariant: ", iz2
229 3 : IF (require_chern) WRITE (iw, '(T2,A,I0)') "TOPOLOGY| Converged first Chern number: ", ichern
230 : END IF
231 : ELSE
232 54 : CALL wannier90_files(qs_env, w_input, iw, ispin, 0, centres, iz2, ichern)
233 : END IF
234 178 : IF (ALLOCATED(previous)) DEALLOCATE (previous)
235 : END DO
236 :
237 58 : IF (iw > 0) THEN
238 : WRITE (iw, '(/,T2,A)') &
239 29 : '!--------------------------------End of Wannier90-----------------------------!'
240 : END IF
241 : END IF
242 12967 : CALL timestop(handle)
243 :
244 12967 : END SUBROUTINE wannier90_interface
245 :
246 : ! **************************************************************************************************
247 : !> \brief Copies a k-point CFM into the all-k-point Wannier work distribution.
248 : !> \param source source CFM; it may be unassociated on ranks outside its k-point group
249 : !> \param destination destination CFM
250 : !> \param para_env communicator spanning the source and destination distributions
251 : ! **************************************************************************************************
252 6468 : SUBROUTINE wannier90_copy_cfm_general(source, destination, para_env)
253 : TYPE(cp_cfm_type), INTENT(IN), POINTER :: source
254 : TYPE(cp_cfm_type), INTENT(IN) :: destination
255 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
256 :
257 : TYPE(cp_cfm_type) :: empty_source
258 :
259 6468 : IF (ASSOCIATED(source)) THEN
260 6468 : CALL cp_cfm_copy_general(source, destination, para_env)
261 : ELSE
262 0 : CALL cp_cfm_copy_general(empty_source, destination, para_env)
263 : END IF
264 :
265 6468 : END SUBROUTINE wannier90_copy_cfm_general
266 :
267 : ! **************************************************************************************************
268 : !> \brief ...
269 : !> \param qs_env ...
270 : !> \param input ...
271 : !> \param iw ...
272 : !> \param ispin selected spin channel
273 : !> \param refinement number of joint mesh doublings
274 : !> \param wcc_out Wilson centres for each closed loop
275 : !> \param z2_value candidate Z2 parity, or -1 when not requested
276 : !> \param chern_value candidate first Chern number, or HUGE(0) when unavailable
277 : ! **************************************************************************************************
278 66 : SUBROUTINE wannier90_files(qs_env, input, iw, ispin, refinement, wcc_out, z2_value, chern_value)
279 : TYPE(qs_environment_type), POINTER :: qs_env
280 : TYPE(section_vals_type), POINTER :: input
281 : INTEGER, INTENT(IN) :: iw, ispin, refinement
282 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: wcc_out(:, :)
283 : INTEGER, INTENT(OUT) :: z2_value, chern_value
284 :
285 : INTEGER, PARAMETER :: num_nnmax = 12
286 :
287 : CHARACTER(len=2) :: asym
288 66 : CHARACTER(len=20), ALLOCATABLE, DIMENSION(:) :: atom_symbols
289 : CHARACTER(len=default_path_length) :: nnkp_file
290 : CHARACTER(len=default_string_length) :: filename, input_kp_scheme, reuse_reason, &
291 : seed_name
292 : CHARACTER(LEN=timestamp_length) :: timestamp
293 : COMPLEX(KIND=dp) :: zmmn
294 66 : COMPLEX(KIND=dp), ALLOCATABLE :: ao_projections(:, :, :), export_coeff(:, :), &
295 66 : export_scalar(:, :), hybrid_projections(:, :, :), initial_amn(:, :, :), &
296 66 : link_matrix(:, :), parity_metric(:, :), parity_phase(:), reference_mo(:, :, :), &
297 66 : scalar_overlap(:, :), soc_h(:, :), soc_u(:, :), soc_xyz(:, :, :), spinor_coeff(:, :, :), &
298 66 : wilson_product(:, :, :)
299 : INTEGER :: aligned_degenerate_blocks, aligned_degenerate_max_size, axis, base_mesh(2), &
300 : counts(2), first_point, i, i_rep, ib, ib1, ib2, ibs, ik, ik2, ikk, ikpgr, iloop, &
301 : initial_projections, ipoint, iunit, ix, iy, iz, jpar, k, kept_index, kpoints_source, &
302 : loop_direction(3), n_rep, nadd, nao, nberry_images, nbs, nelectron, nexcl, nkp, nloop, &
303 : nmo, nntot, npoint, nscalar, num_atoms, num_bands, num_bands_tot, num_kpts, num_wann, &
304 : spin_channel, state_components, state_unit, status, strong, tqc_dimension, trim_id, &
305 : weak(3), z4
306 66 : INTEGER, ALLOCATABLE :: ebr_coefficients(:, :), hybrid_groups(:, :), loop_index(:), &
307 66 : parity_map(:), parity_odd(:), projection_pivots(:)
308 66 : INTEGER, ALLOCATABLE, DIMENSION(:) :: band_indices, exclude_bands
309 66 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nblist, nnlist
310 66 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: nncell
311 : INTEGER, DIMENSION(2) :: kp_range
312 : INTEGER, DIMENSION(3) :: input_nkp_grid, mp_grid
313 66 : INTEGER, DIMENSION(:), POINTER :: invals
314 66 : INTEGER, DIMENSION(:, :, :), POINTER :: berry_cell_index, cell_to_index
315 : LOGICAL :: diis_step, do_chern, do_kpoints, do_parity, do_soc, do_tqc, do_wilson, do_z2, &
316 : export_state, full_mesh_diagonalized, input_full_grid, input_gamma_centered, &
317 : input_kpoint_symmetry, lowest_bands, mp_grid_explicit, mp_grid_valid, my_kpgrp, mygrp, &
318 : nonnegative_atomic, ordered_berry, require_global_gap, reuse_scf_mos, reused_scf_mos, &
319 : signed_atomic, time_reversal, use_bloch_phases, validate_reuse_ok, &
320 : validate_reuse_scf_mos, window_set(4)
321 66 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: allowed_band, keep_band
322 : REAL(KIND=dp) :: aligned_degenerate_min_svalue, berry_phase, chern_winding, cmmn, &
323 : conduction_min, direct_gap, energy_window(4), gap_tol, gauge_arg, gauge_imag, gauge_real, &
324 : gauge_tmp, ksign, link_sv, loop_origin(3), pair_gap, parity_checks(4), &
325 : parity_energy_error, parity_energy_tolerance, parity_error, parity_gap, parity_origin(3), &
326 : parity_tolerance, projection_min_svalue, reuse_candidate_deviation, &
327 : reuse_candidate_metric_deviation, reuse_candidate_min_svalue, reuse_candidate_residual, &
328 : rmmn, transverse(3), valence_max, validation_eigenvalue_deviation, validation_min_svalue
329 : REAL(KIND=dp) :: validation_subspace_deviation, wkp_ref
330 66 : REAL(KIND=dp), ALLOCATABLE :: loop_sv(:), parity_energies(:), &
331 66 : projection_eigenvalues(:, :), &
332 66 : scalar_values(:), spinor_values(:, :)
333 66 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigval
334 66 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: atoms_cart, b_latt, kpt_latt, &
335 66 : reference_eigenvalues
336 : REAL(KIND=dp), DIMENSION(3) :: bvec, input_kp_shift, phase_center
337 : REAL(KIND=dp), DIMENSION(3, 3) :: h_inv, real_lattice, recip_lattice
338 66 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, wkp, wkp_source
339 66 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp, xkp_source
340 66 : REAL(KIND=dp), POINTER :: rvals(:)
341 66 : TYPE(berry_matrix_type), DIMENSION(:), POINTER :: berry_matrix
342 : TYPE(cell_type), POINTER :: cell
343 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
344 : TYPE(cp_cfm_type) :: fmk1_cfm, fmk2_cfm, mmn_cfm, omat_cfm, &
345 : tmp_cfm
346 : TYPE(cp_cfm_type), POINTER :: cmo_coeff_source
347 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_ao, matrix_struct_mmn, &
348 : matrix_struct_work
349 : TYPE(cp_fm_type) :: mat_imag, mat_real
350 66 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s, soc_matrices
351 : TYPE(dbcsr_type), POINTER :: cmatrix, cmatrix_full, loop_imag, &
352 : loop_real, rmatrix, rmatrix_full
353 : TYPE(dft_control_type), POINTER :: dft_control
354 : TYPE(kpoint_env_type), POINTER :: kp
355 : TYPE(kpoint_type), POINTER :: berry_kpoint, kpoint, qs_kpoint
356 66 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
357 : TYPE(mp_para_env_type), POINTER :: para_env
358 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
359 66 : POINTER :: overlap_nl, sab_nl
360 66 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
361 : TYPE(qs_environment_type), POINTER :: qs_env_kp
362 : TYPE(qs_scf_env_type), POINTER :: scf_env
363 : TYPE(scf_control_type), POINTER :: scf_control
364 :
365 : !--------------------------------------------------------------------------------------------!
366 :
367 : ! generate all arrays needed for the setup call
368 66 : CALL section_vals_val_get(input, "SEED_NAME", c_val=seed_name)
369 66 : CALL get_qs_env(qs_env, dft_control=dft_control)
370 66 : IF (dft_control%nspins == 2) THEN
371 6 : IF (LEN_TRIM(seed_name) + 5 > LEN(seed_name)) THEN
372 0 : CPABORT('Wannier90 seed is too long for the spin suffix.')
373 : END IF
374 6 : IF (ispin == 1) seed_name = TRIM(seed_name)//'_up'
375 6 : IF (ispin == 2) seed_name = TRIM(seed_name)//'_down'
376 6 : IF (iw > 0) WRITE (iw, '(T2,A,I0,2A)') 'WANNIER90| Spin channel ', ispin, ' file prefix ', TRIM(seed_name)
377 : END IF
378 66 : CALL section_vals_val_get(input, "MP_GRID", i_vals=invals, explicit=mp_grid_explicit)
379 66 : CALL section_vals_val_get(input, "KPOINTS_SOURCE", i_val=kpoints_source)
380 66 : ordered_berry = kpoints_source >= w90_kpoints_nnkp
381 66 : CALL section_vals_val_get(input, "NNKP_FILE", c_val=nnkp_file)
382 66 : spin_channel = ispin
383 66 : CALL section_vals_val_get(input, "WILSON_LOOP", l_val=do_wilson)
384 66 : CALL section_vals_val_get(input, "Z2", l_val=do_z2)
385 66 : CALL section_vals_val_get(input, "CHERN", l_val=do_chern)
386 66 : CALL section_vals_val_get(input, "REQUIRE_GLOBAL_GAP", l_val=require_global_gap)
387 66 : CALL section_vals_val_get(input, "STATE_EXPORT", l_val=export_state)
388 66 : CALL section_vals_val_get(input, "PARITY", l_val=do_parity)
389 66 : CALL section_vals_val_get(input, "INVERSION_TQC", l_val=do_tqc)
390 66 : CALL section_vals_val_get(input, "TQC_DIMENSION", i_val=tqc_dimension)
391 66 : CALL section_vals_val_get(input, "PARITY_ORIGIN", r_vals=rvals)
392 264 : parity_origin = rvals
393 66 : CALL section_vals_val_get(input, "PARITY_TOLERANCE", r_val=parity_tolerance)
394 66 : CALL section_vals_val_get(input, "PARITY_ENERGY_TOL", r_val=parity_energy_tolerance)
395 66 : do_parity = do_parity .OR. do_tqc .OR. kpoints_source == w90_kpoints_trim
396 66 : IF (do_parity .AND. .NOT. ordered_berry) THEN
397 0 : CPABORT("PARITY requires explicit NNKP, WILSON or TRIM points.")
398 : END IF
399 66 : IF (tqc_dimension < 2 .OR. tqc_dimension > 3) CPABORT("TQC_DIMENSION must be 2 or 3.")
400 66 : IF (parity_tolerance <= 0.0_dp .OR. parity_energy_tolerance <= 0.0_dp) THEN
401 0 : CPABORT("Invalid parity tolerance.")
402 : END IF
403 66 : IF (export_state .AND. .NOT. ordered_berry) THEN
404 0 : CPABORT("STATE_EXPORT requires explicit NNKP, WILSON or TRIM points.")
405 : END IF
406 66 : CALL section_vals_val_get(input, "TIME_REVERSAL", l_val=time_reversal)
407 66 : CALL section_vals_val_get(input, "SOC", l_val=do_soc)
408 66 : IF (do_tqc .AND. (.NOT. do_soc .OR. .NOT. time_reversal)) THEN
409 0 : CPABORT("INVERSION_TQC requires SOC T and TIME_REVERSAL T.")
410 : END IF
411 66 : CALL section_vals_val_get(input, "WILSON_GAP_TOL", r_val=gap_tol)
412 66 : IF (gap_tol <= 0.0_dp) CPABORT("WILSON_GAP_TOL must be positive.")
413 66 : z2_value = -1
414 66 : chern_value = HUGE(0)
415 66 : do_wilson = do_wilson .OR. kpoints_source == w90_kpoints_wilson
416 66 : IF (do_wilson .AND. kpoints_source == w90_kpoints_trim) THEN
417 0 : CPABORT("TRIM points are not Wilson loops; use KPOINTS_SOURCE WILSON.")
418 : END IF
419 66 : IF (do_wilson) CALL cite_reference(Gresch2017)
420 66 : IF (do_z2) CALL cite_reference(Soluyanov2011)
421 66 : IF ((do_wilson .OR. do_soc) .AND. kpoints_source < w90_kpoints_nnkp) THEN
422 0 : CPABORT("Native Wilson/SOC requires KPOINTS_SOURCE NNKP or WILSON.")
423 : END IF
424 66 : IF (do_z2 .AND. (kpoints_source /= w90_kpoints_wilson .OR. .NOT. do_soc .OR. .NOT. time_reversal)) THEN
425 0 : CPABORT("Z2 requires KPOINTS_SOURCE WILSON, SOC T, and TIME_REVERSAL T.")
426 : END IF
427 66 : IF (do_chern .AND. (kpoints_source /= w90_kpoints_wilson .OR. do_z2)) THEN
428 0 : CPABORT("CHERN requires a full WILSON surface and cannot be combined with Z2.")
429 : END IF
430 66 : IF (require_global_gap .AND. .NOT. do_wilson) THEN
431 0 : CPABORT("REQUIRE_GLOBAL_GAP requires Wilson analysis.")
432 : END IF
433 66 : CALL get_qs_env(qs_env, dft_control=dft_control)
434 66 : IF (spin_channel < 1 .OR. spin_channel > dft_control%nspins) THEN
435 0 : CPABORT("WANNIER90%SPIN_CHANNEL is not available in this calculation.")
436 : END IF
437 66 : IF (do_soc .AND. dft_control%nspins /= 1) THEN
438 0 : CPABORT("WANNIER90 SOC currently requires a restricted SCF.")
439 : END IF
440 66 : CALL section_vals_val_get(input, "WANNIER_FUNCTIONS", n_rep_val=n_rep)
441 66 : IF (n_rep > 1 .AND. n_rep /= dft_control%nspins) THEN
442 0 : CPABORT('Specify WANNIER_FUNCTIONS once, or once per spin channel.')
443 : END IF
444 66 : i_rep = 1
445 66 : IF (n_rep > 1) i_rep = ispin
446 66 : CALL section_vals_val_get(input, "WANNIER_FUNCTIONS", i_rep_val=i_rep, i_val=num_wann)
447 66 : CALL section_vals_val_get(input, "ADDED_MOS", i_val=nadd)
448 66 : CALL section_vals_val_get(input, "REUSE_SCF_MOS", l_val=reuse_scf_mos)
449 66 : CALL section_vals_val_get(input, "VALIDATE_REUSE_SCF_MOS", l_val=validate_reuse_scf_mos)
450 66 : CALL section_vals_val_get(input, "USE_BLOCH_PHASES", l_val=use_bloch_phases)
451 66 : CALL section_vals_val_get(input, "INITIAL_PROJECTIONS", i_val=initial_projections)
452 66 : IF (initial_projections == -1 .AND. use_bloch_phases) initial_projections = 0
453 66 : energy_window = [-HUGE(0.0_dp)/100, HUGE(0.0_dp)/100, -HUGE(0.0_dp)/100, HUGE(0.0_dp)/100]
454 330 : DO i = 1, 4
455 264 : CALL section_vals_val_get(input, TRIM(window_keywords(i)), explicit=window_set(i))
456 330 : IF (window_set(i)) CALL section_vals_val_get(input, TRIM(window_keywords(i)), r_val=energy_window(i))
457 : END DO
458 322 : IF (ordered_berry .AND. (initial_projections > 0 .OR. ANY(window_set))) THEN
459 0 : CPABORT("AO projections and energy windows require KPOINTS_SOURCE MP_GRID or SCF.")
460 : END IF
461 66 : IF (window_set(3) .AND. .NOT. window_set(4)) CPABORT("DIS_FROZ_MIN requires DIS_FROZ_MAX.")
462 66 : IF (energy_window(1) > energy_window(2) .OR. energy_window(3) > energy_window(4)) THEN
463 0 : CPABORT("Invalid Wannier90 energy window.")
464 : END IF
465 66 : reuse_scf_mos = reuse_scf_mos .AND. kpoints_source == w90_kpoints_scf
466 66 : validate_reuse_scf_mos = validate_reuse_scf_mos .AND. reuse_scf_mos
467 264 : mp_grid(1:3) = invals(1:3)
468 : ! excluded bands
469 66 : CALL section_vals_val_get(input, "EXCLUDE_BANDS", n_rep_val=n_rep)
470 66 : nexcl = 0
471 92 : DO i_rep = 1, n_rep
472 26 : CALL section_vals_val_get(input, "EXCLUDE_BANDS", i_rep_val=i_rep, i_vals=invals)
473 92 : nexcl = nexcl + SIZE(invals)
474 : END DO
475 66 : IF (nexcl > 0) THEN
476 72 : ALLOCATE (exclude_bands(nexcl))
477 24 : nexcl = 0
478 50 : DO i_rep = 1, n_rep
479 26 : CALL section_vals_val_get(input, "EXCLUDE_BANDS", i_rep_val=i_rep, i_vals=invals)
480 144 : exclude_bands(nexcl + 1:nexcl + SIZE(invals)) = invals(:)
481 50 : nexcl = nexcl + SIZE(invals)
482 : END DO
483 : END IF
484 : !
485 : ! lattice -> Angstrom
486 66 : CALL get_qs_env(qs_env, cell=cell)
487 66 : CALL get_cell(cell, h=real_lattice, h_inv=h_inv)
488 : ! k-points
489 66 : CALL get_qs_env(qs_env, particle_set=particle_set)
490 66 : CALL get_qs_env(qs_env, para_env=para_env)
491 66 : phase_center = 0.0_dp
492 302 : DO i = 1, SIZE(particle_set)
493 3842 : phase_center(1:3) = phase_center(1:3) + MATMUL(h_inv, particle_set(i)%r)
494 : END DO
495 264 : phase_center(1:3) = phase_center(1:3)/REAL(SIZE(particle_set), KIND=dp)
496 264 : phase_center(1:3) = phase_center(1:3) - FLOOR(phase_center(1:3))
497 66 : recip_lattice(1:3, 1:3) = h_inv(1:3, 1:3)
498 858 : real_lattice(1:3, 1:3) = angstrom*real_lattice(1:3, 1:3)
499 1650 : recip_lattice(1:3, 1:3) = (twopi/angstrom)*TRANSPOSE(recip_lattice(1:3, 1:3))
500 66 : NULLIFY (kpoint, qs_kpoint, xkp, wkp, xkp_source, wkp_source)
501 66 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints, kpoints=qs_kpoint)
502 66 : input_kpoint_symmetry = .FALSE.
503 66 : input_full_grid = .FALSE.
504 66 : input_kp_scheme = ""
505 66 : IF (do_kpoints .AND. ASSOCIATED(qs_kpoint)) THEN
506 : CALL get_kpoint_info(qs_kpoint, kp_scheme=input_kp_scheme, nkp_grid=input_nkp_grid, &
507 : kp_shift=input_kp_shift, symmetry=input_kpoint_symmetry, &
508 : full_grid=input_full_grid, gamma_centered=input_gamma_centered, &
509 66 : nkp=nkp, xkp=xkp, wkp=wkp)
510 : END IF
511 66 : CALL kpoint_create(kpoint)
512 :
513 22 : SELECT CASE (kpoints_source)
514 : CASE (w90_kpoints_nnkp, w90_kpoints_wilson, w90_kpoints_trim)
515 22 : IF (kpoints_source == w90_kpoints_nnkp) THEN
516 4 : CALL read_wannier90_nnkpts(nnkp_file, real_lattice, recip_lattice, kpt_latt, nnlist, nncell)
517 18 : ELSE IF (kpoints_source == w90_kpoints_trim) THEN
518 6 : num_kpts = 2**tqc_dimension
519 42 : ALLOCATE (kpt_latt(3, num_kpts), nnlist(num_kpts, 1), nncell(3, num_kpts, 1))
520 6 : kpt_latt = 0.0_dp
521 6 : nncell = 0
522 38 : DO i = 1, num_kpts
523 112 : DO axis = 1, tqc_dimension
524 112 : IF (BTEST(i - 1, axis - 1)) kpt_latt(axis, i) = 0.5_dp
525 : END DO
526 38 : nnlist(i, 1) = i
527 : END DO
528 : ELSE
529 12 : CALL section_vals_val_get(input, "WILSON_MESH", i_vals=invals)
530 36 : base_mesh = invals
531 12 : IF (base_mesh(1) < 2 .OR. base_mesh(2) < 2) THEN
532 0 : CPABORT("WILSON_MESH entries must be at least two.")
533 : END IF
534 12 : npoint = base_mesh(1)*2**refinement
535 12 : nloop = (base_mesh(2) - 1)*2**refinement + 1
536 12 : CALL section_vals_val_get(input, "WILSON_DIRECTION", i_vals=invals)
537 48 : loop_direction = invals
538 12 : IF (ALL(loop_direction == 0)) CPABORT("WILSON_DIRECTION must be nonzero.")
539 12 : CALL section_vals_val_get(input, "WILSON_ORIGIN", r_vals=rvals)
540 48 : loop_origin = rvals
541 12 : CALL section_vals_val_get(input, "WILSON_TRANSVERSE", r_vals=rvals)
542 48 : transverse = rvals
543 : bvec = [loop_direction(2)*transverse(3) - loop_direction(3)*transverse(2), &
544 : loop_direction(3)*transverse(1) - loop_direction(1)*transverse(3), &
545 48 : loop_direction(1)*transverse(2) - loop_direction(2)*transverse(1)]
546 48 : IF (SUM(bvec**2) < 1.e-20_dp) THEN
547 0 : CPABORT("Wilson loop and transverse vectors must be linearly independent.")
548 : END IF
549 48 : IF (do_chern .AND. MAXVAL(ABS(transverse - NINT(transverse))) > 1.e-10_dp) THEN
550 0 : CPABORT("CHERN requires integer WILSON_TRANSVERSE: a closed full surface, not a half-plane.")
551 : END IF
552 12 : IF (do_z2) THEN
553 28 : IF (MAXVAL(ABS(2*loop_origin - NINT(2*loop_origin))) > 1.e-10_dp .OR. &
554 : MAXVAL(ABS(2*transverse - NINT(2*transverse))) > 1.e-10_dp) THEN
555 0 : CPABORT("Z2 surface boundaries must pass through time-reversal-invariant momenta.")
556 : END IF
557 : ! Initially restrict native Z2 to standard half-planes, avoiding multiple coverings.
558 : IF (SUM(ABS(loop_direction)) /= 1 .OR. &
559 40 : ABS(SUM(ABS(transverse)) - 0.5_dp) > 1.e-10_dp .OR. &
560 : ABS(DOT_PRODUCT(REAL(loop_direction, dp), transverse)) > 1.e-10_dp) THEN
561 0 : CPABORT("Native Z2 requires distinct coordinate axes with windings one and one half.")
562 : END IF
563 : END IF
564 12 : num_kpts = npoint*nloop
565 84 : ALLOCATE (kpt_latt(3, num_kpts), nnlist(num_kpts, 1), nncell(3, num_kpts, 1))
566 12 : nncell = 0
567 60 : DO iloop = 1, nloop
568 48 : first_point = (iloop - 1)*npoint + 1
569 360 : DO ipoint = 1, npoint
570 312 : i = first_point + ipoint - 1
571 : kpt_latt(:, i) = loop_origin + REAL(iloop - 1, dp)/REAL(nloop - 1, dp)*transverse + &
572 1248 : REAL(ipoint - 1, dp)/REAL(npoint, dp)*loop_direction
573 360 : nnlist(i, 1) = i + 1
574 : END DO
575 48 : nnlist(i, 1) = first_point
576 204 : nncell(:, i, 1) = loop_direction
577 : END DO
578 : END IF
579 22 : num_kpts = SIZE(kpt_latt, 2)
580 22 : nntot = SIZE(nnlist, 2)
581 22 : kpoint%kp_scheme = "GENERAL"
582 22 : kpoint%symmetry = .FALSE.
583 22 : kpoint%verbose = .FALSE.
584 22 : kpoint%full_grid = .TRUE.
585 22 : kpoint%eps_geo = 1.0e-6_dp
586 22 : kpoint%use_real_wfn = .FALSE.
587 22 : kpoint%parallel_group_size = para_env%num_pe
588 22 : kpoint%nkp = num_kpts
589 110 : ALLOCATE (kpoint%xkp(3, num_kpts), kpoint%wkp(num_kpts))
590 1462 : kpoint%xkp = kpt_latt
591 : ! Export weights do not enter the converged SCF density.
592 382 : kpoint%wkp = 1.0_dp/REAL(num_kpts, KIND=dp)
593 22 : IF (iw > 0) WRITE (iw, '(T2,A,I0,A,I0)') &
594 11 : "WANNIER90| Explicit points: ", num_kpts, ", neighbours per point: ", nntot
595 : CASE (w90_kpoints_mp_grid)
596 0 : num_kpts = mp_grid(1)*mp_grid(2)*mp_grid(3)
597 0 : ALLOCATE (kpt_latt(3, num_kpts))
598 0 : kpoint%kp_scheme = "MONKHORST-PACK"
599 0 : kpoint%symmetry = .FALSE.
600 0 : kpoint%nkp_grid(1:3) = mp_grid(1:3)
601 0 : kpoint%verbose = .FALSE.
602 0 : kpoint%full_grid = .TRUE.
603 0 : kpoint%eps_geo = 1.0e-6_dp
604 0 : kpoint%use_real_wfn = .FALSE.
605 0 : kpoint%parallel_group_size = para_env%num_pe
606 0 : i = 0
607 0 : DO ix = 0, mp_grid(1) - 1
608 0 : DO iy = 0, mp_grid(2) - 1
609 0 : DO iz = 0, mp_grid(3) - 1
610 0 : i = i + 1
611 0 : kpt_latt(1, i) = REAL(ix, KIND=dp)/REAL(mp_grid(1), KIND=dp)
612 0 : kpt_latt(2, i) = REAL(iy, KIND=dp)/REAL(mp_grid(2), KIND=dp)
613 0 : kpt_latt(3, i) = REAL(iz, KIND=dp)/REAL(mp_grid(3), KIND=dp)
614 : END DO
615 : END DO
616 : END DO
617 0 : kpoint%nkp = num_kpts
618 0 : ALLOCATE (kpoint%xkp(3, num_kpts), kpoint%wkp(num_kpts))
619 0 : kpoint%wkp(:) = 1._dp/REAL(num_kpts, KIND=dp)
620 : ! rskp_transform uses reduced coordinates in exp(2*pi*i*k.R), with integer cell shifts R.
621 0 : kpoint%xkp(:, :) = kpt_latt(:, :)
622 :
623 : CASE (w90_kpoints_scf)
624 44 : IF (.NOT. do_kpoints .OR. .NOT. ASSOCIATED(qs_kpoint)) THEN
625 0 : CPABORT("WANNIER90%KPOINTS_SOURCE SCF requires an active DFT%KPOINTS section.")
626 : END IF
627 44 : SELECT CASE (TRIM(input_kp_scheme))
628 : CASE ("GAMMA")
629 0 : mp_grid(:) = 1
630 0 : num_kpts = 1
631 0 : ALLOCATE (kpt_latt(3, num_kpts))
632 0 : kpt_latt(1:3, 1) = 0.0_dp
633 0 : kpoint%kp_scheme = "GAMMA"
634 0 : kpoint%symmetry = .FALSE.
635 0 : kpoint%verbose = .FALSE.
636 0 : kpoint%full_grid = .TRUE.
637 0 : kpoint%eps_geo = 1.0e-6_dp
638 0 : kpoint%use_real_wfn = .FALSE.
639 0 : kpoint%parallel_group_size = para_env%num_pe
640 0 : kpoint%nkp = num_kpts
641 0 : ALLOCATE (kpoint%xkp(3, num_kpts), kpoint%wkp(num_kpts))
642 0 : kpoint%xkp(1:3, 1) = 0.0_dp
643 0 : kpoint%wkp(1) = 1.0_dp
644 :
645 : CASE ("MONKHORST-PACK", "MACDONALD")
646 38 : mp_grid(1:3) = input_nkp_grid(1:3)
647 38 : kpoint%kp_scheme = input_kp_scheme
648 38 : kpoint%symmetry = .FALSE.
649 152 : kpoint%nkp_grid(1:3) = input_nkp_grid(1:3)
650 152 : kpoint%kp_shift(1:3) = input_kp_shift(1:3)
651 38 : kpoint%gamma_centered = input_gamma_centered
652 38 : kpoint%verbose = .FALSE.
653 38 : kpoint%full_grid = .TRUE.
654 38 : kpoint%eps_geo = 1.0e-6_dp
655 38 : kpoint%use_real_wfn = .FALSE.
656 38 : kpoint%parallel_group_size = para_env%num_pe
657 38 : CALL kpoint_initialize(kpoint, particle_set, cell)
658 38 : num_kpts = kpoint%nkp
659 114 : ALLOCATE (kpt_latt(3, num_kpts))
660 2150 : kpt_latt(1:3, 1:num_kpts) = kpoint%xkp(1:3, 1:num_kpts)
661 38 : IF (input_kpoint_symmetry .AND. .NOT. input_full_grid .AND. iw > 0) THEN
662 : WRITE (iw, '(T2,A)') &
663 18 : "WANNIER90| SCF k-points are symmetry-reduced; regenerating the full SCF mesh."
664 18 : IF (reuse_scf_mos) THEN
665 : WRITE (iw, '(T2,A)') &
666 17 : "WANNIER90| CP2K will try to reconstruct the full-mesh MOs from the SCF orbitals."
667 : ELSE
668 : WRITE (iw, '(T2,A)') &
669 1 : "WANNIER90| The full exported mesh is diagonalized for the Wannier90 files."
670 : END IF
671 : END IF
672 :
673 : CASE ("GENERAL")
674 6 : IF (ASSOCIATED(qs_kpoint%xkp_input)) THEN
675 6 : xkp_source => qs_kpoint%xkp_input
676 6 : wkp_source => qs_kpoint%wkp_input
677 : ELSE
678 0 : xkp_source => xkp
679 0 : wkp_source => wkp
680 : END IF
681 6 : IF (.NOT. ASSOCIATED(xkp_source) .OR. .NOT. ASSOCIATED(wkp_source)) THEN
682 0 : CPABORT("Could not access the SCF GENERAL k-point set for the Wannier90 export.")
683 : END IF
684 6 : num_kpts = SIZE(wkp_source)
685 18 : ALLOCATE (kpt_latt(3, num_kpts))
686 198 : kpt_latt(1:3, 1:num_kpts) = xkp_source(1:3, 1:num_kpts)
687 6 : IF (mp_grid_explicit) THEN
688 0 : IF (mp_grid(1)*mp_grid(2)*mp_grid(3) /= num_kpts) THEN
689 0 : CPABORT("WANNIER90%MP_GRID must contain exactly as many points as the SCF GENERAL mesh.")
690 : END IF
691 : ELSE
692 6 : CALL infer_wannier_mp_grid(kpt_latt, mp_grid, mp_grid_valid)
693 6 : IF (.NOT. mp_grid_valid) THEN
694 0 : CPABORT("Could not infer WANNIER90%MP_GRID from the SCF GENERAL mesh.")
695 : END IF
696 : END IF
697 6 : wkp_ref = 1.0_dp/REAL(num_kpts, KIND=dp)
698 54 : DO i = 1, num_kpts
699 54 : IF (ABS(wkp_source(i) - wkp_ref) > 1.0e-10_dp) THEN
700 0 : CPABORT("WANNIER90%KPOINTS_SOURCE SCF requires equally weighted GENERAL k-points.")
701 : END IF
702 : END DO
703 6 : kpoint%kp_scheme = "GENERAL"
704 6 : kpoint%symmetry = .FALSE.
705 24 : kpoint%nkp_grid(1:3) = mp_grid(1:3)
706 6 : kpoint%verbose = .FALSE.
707 6 : kpoint%full_grid = .TRUE.
708 6 : kpoint%eps_geo = 1.0e-6_dp
709 6 : kpoint%use_real_wfn = .FALSE.
710 6 : kpoint%parallel_group_size = para_env%num_pe
711 6 : kpoint%nkp = num_kpts
712 30 : ALLOCATE (kpoint%xkp(3, num_kpts), kpoint%wkp(num_kpts))
713 390 : kpoint%xkp(1:3, 1:num_kpts) = xkp_source(1:3, 1:num_kpts)
714 54 : kpoint%wkp(1:num_kpts) = wkp_ref
715 6 : IF (input_kpoint_symmetry .AND. .NOT. input_full_grid .AND. iw > 0) THEN
716 : WRITE (iw, '(T2,A)') &
717 2 : "WANNIER90| SCF k-points are symmetry-reduced; using the full input GENERAL mesh."
718 2 : IF (reuse_scf_mos) THEN
719 : WRITE (iw, '(T2,A)') &
720 2 : "WANNIER90| CP2K will try to reconstruct the full-mesh MOs from the SCF orbitals."
721 : ELSE
722 : WRITE (iw, '(T2,A)') &
723 0 : "WANNIER90| The full exported mesh is diagonalized for the Wannier90 files."
724 : END IF
725 : END IF
726 :
727 : CASE DEFAULT
728 44 : CPABORT("WANNIER90%KPOINTS_SOURCE SCF does not support this DFT%KPOINTS scheme.")
729 : END SELECT
730 : CASE DEFAULT
731 66 : CPABORT("Unknown WANNIER90%KPOINTS_SOURCE setting.")
732 : END SELECT
733 : ! number of bands in calculation
734 66 : CALL get_qs_env(qs_env, mos=mos)
735 66 : CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=num_bands_tot, nelectron=nelectron)
736 66 : num_bands_tot = MIN(nao, num_bands_tot + nadd)
737 66 : nscalar = num_bands_tot
738 66 : IF (do_soc) num_bands_tot = 2*nscalar
739 198 : ALLOCATE (keep_band(num_bands_tot))
740 418 : keep_band = .TRUE.
741 184 : DO i = 1, nexcl
742 118 : ib = exclude_bands(i)
743 118 : IF (ib < 1 .OR. ib > num_bands_tot) CPABORT("WANNIER90%EXCLUDE_BANDS: index out of range.")
744 118 : IF (ordered_berry .AND. .NOT. keep_band(ib)) THEN
745 0 : CPABORT("WANNIER90%EXCLUDE_BANDS: duplicate band index.")
746 : END IF
747 184 : keep_band(ib) = .FALSE.
748 : END DO
749 418 : num_bands = COUNT(keep_band)
750 66 : IF (num_bands < 1 .OR. num_wann > num_bands) THEN
751 0 : CPABORT("WANNIER90 EXCLUDE_BANDS leaves too few bands for WANNIER_FUNCTIONS.")
752 : END IF
753 198 : ALLOCATE (band_indices(num_bands))
754 770 : band_indices(:) = PACK([(i, i=1, num_bands_tot)], keep_band)
755 530 : lowest_bands = num_bands < num_bands_tot .AND. ALL(band_indices == [(i, i=1, num_bands)])
756 66 : IF (require_global_gap .AND. .NOT. lowest_bands) THEN
757 0 : CPABORT("REQUIRE_GLOBAL_GAP needs a lowest-band prefix and an excluded band above it.")
758 : END IF
759 66 : IF (do_z2) THEN
760 4 : IF (num_bands /= nelectron) THEN
761 0 : CPABORT("Z2 requires one occupied spinor per electron; adjust EXCLUDE_BANDS.")
762 : END IF
763 4 : IF (num_bands == num_bands_tot .OR. MOD(num_bands, 2) /= 0) THEN
764 0 : CPABORT("Z2 needs an even occupied spinor subspace and at least one excluded conduction band.")
765 : END IF
766 68 : IF (ANY(band_indices /= [(i, i=1, num_bands)])) THEN
767 0 : CPABORT("Z2 requires the lowest occupied spinor bands.")
768 : END IF
769 : END IF
770 66 : IF (do_wilson) THEN
771 16 : IF (nntot /= 1) CPABORT("WILSON_LOOP requires exactly one directed neighbour per point.")
772 48 : ALLOCATE (loop_index(num_kpts))
773 68 : i = 1
774 68 : nloop = 0
775 68 : DO WHILE (i <= num_kpts)
776 52 : first_point = i
777 52 : nloop = nloop + 1
778 : DO
779 328 : loop_index(i) = nloop
780 328 : IF (nnlist(i, 1) == first_point) EXIT
781 276 : IF (nnlist(i, 1) /= i + 1 .OR. i == num_kpts) THEN
782 0 : CPABORT("WILSON_LOOP requires contiguous ordered closed loops.")
783 : END IF
784 52 : i = i + 1
785 : END DO
786 52 : i = i + 1
787 : END DO
788 112 : ALLOCATE (wilson_product(num_bands, num_bands, nloop), loop_sv(nloop))
789 16 : wilson_product = CMPLX(0.0_dp, 0.0_dp, dp)
790 60 : DO i = 1, num_bands
791 224 : wilson_product(i, i, :) = CMPLX(1.0_dp, 0.0_dp, dp)
792 : END DO
793 68 : loop_sv = 1.0_dp
794 64 : ALLOCATE (wcc_out(num_bands, nloop))
795 : ELSE
796 50 : ALLOCATE (wcc_out(0, 0))
797 : END IF
798 264 : ALLOCATE (link_matrix(num_bands, num_bands))
799 66 : IF (initial_projections == 0 .AND. num_wann /= num_bands .AND. .NOT. ordered_berry) THEN
800 0 : CPABORT("IDENTITY projections require WANNIER_FUNCTIONS to match the number of bands.")
801 : END IF
802 66 : IF (initial_projections >= 0 .AND. .NOT. ordered_berry .AND. num_wann < 1) THEN
803 0 : CPABORT("Initial projections require a positive WANNIER_FUNCTIONS count.")
804 : END IF
805 322 : IF (ANY(window_set) .AND. num_wann == num_bands) THEN
806 0 : CPABORT("Energy windows require more bands than Wannier functions.")
807 : END IF
808 322 : IF (ANY(window_set) .AND. num_wann < 1) THEN
809 0 : CPABORT("Energy windows require a positive WANNIER_FUNCTIONS count.")
810 : END IF
811 198 : ALLOCATE (allowed_band(num_bands))
812 328 : IF (ANY(window_set)) ALLOCATE (projection_eigenvalues(num_bands, num_kpts))
813 66 : num_atoms = SIZE(particle_set)
814 198 : ALLOCATE (atoms_cart(3, num_atoms))
815 198 : ALLOCATE (atom_symbols(num_atoms))
816 302 : DO i = 1, num_atoms
817 944 : atoms_cart(1:3, i) = particle_set(i)%r(1:3)
818 236 : CALL get_atomic_kind(particle_set(i)%atomic_kind, element_symbol=asym)
819 302 : atom_symbols(i) = asym
820 : END DO
821 66 : IF (kpoints_source < w90_kpoints_nnkp) THEN
822 264 : ALLOCATE (nnlist(num_kpts, num_nnmax), nncell(3, num_kpts, num_nnmax))
823 44 : nnlist(:, :) = 0
824 44 : nncell(:, :, :) = 0
825 44 : nntot = 0
826 44 : IF (iw > 0) THEN
827 : CALL wannier_setup(mp_grid, num_kpts, real_lattice, recip_lattice, &
828 22 : kpt_latt, nntot, nnlist, nncell, iw)
829 : END IF
830 44 : CALL para_env%sum(nntot)
831 44 : CALL para_env%sum(nnlist)
832 44 : CALL para_env%sum(nncell)
833 : END IF
834 :
835 66 : IF (para_env%is_source() .AND. kpoints_source < w90_kpoints_nnkp) THEN
836 : ! Write the ordinary file-export input "seed_name.win".
837 22 : WRITE (filename, '(A,A)') TRIM(seed_name), ".win"
838 22 : CALL open_file(filename, unit_number=iunit, file_status="REPLACE", file_action="WRITE")
839 : !
840 22 : CALL m_timestamp(timestamp)
841 22 : WRITE (iunit, "(A)") "! Wannier90 input file generated by CP2K "
842 22 : WRITE (iunit, "(A,/)") "! Creation date "//timestamp
843 : !
844 22 : WRITE (iunit, "(A,I5)") "num_wann = ", num_wann
845 22 : IF (num_bands /= num_wann .OR. use_bloch_phases) THEN
846 20 : WRITE (iunit, "(A,I5)") "num_bands = ", num_bands
847 : END IF
848 110 : DO i = 1, 4
849 110 : IF (window_set(i)) WRITE (iunit, '(2A,ES24.16)') TRIM(window_keywords(i)), ' = ', energy_window(i)*evolt
850 : END DO
851 22 : WRITE (iunit, "(/,A,/)") "length_unit = bohr"
852 22 : WRITE (iunit, "(/,A,/)") "! System"
853 22 : WRITE (iunit, "(/,A)") "begin unit_cell_cart"
854 22 : WRITE (iunit, "(A)") "bohr"
855 88 : DO i = 1, 3
856 286 : WRITE (iunit, "(3ES24.16)") cell%hmat(1:3, i)
857 : END DO
858 22 : WRITE (iunit, "(A,/)") "end unit_cell_cart"
859 22 : WRITE (iunit, "(/,A)") "begin atoms_cart"
860 22 : WRITE (iunit, "(A)") "bohr"
861 128 : DO i = 1, num_atoms
862 128 : WRITE (iunit, "(A,3F15.10)") atom_symbols(i), atoms_cart(1:3, i)
863 : END DO
864 22 : WRITE (iunit, "(A,/)") "end atoms_cart"
865 22 : WRITE (iunit, "(/,A,/)") "! Kpoints"
866 22 : WRITE (iunit, "(/,A,3I6/)") "mp_grid = ", mp_grid(1:3)
867 22 : WRITE (iunit, "(A)") "begin kpoints"
868 310 : DO i = 1, num_kpts
869 310 : WRITE (iunit, "(3ES24.16)") kpt_latt(1:3, i)
870 : END DO
871 22 : WRITE (iunit, "(A)") "end kpoints"
872 22 : CALL close_file(iunit)
873 : ELSE
874 44 : iunit = -1
875 : END IF
876 :
877 : ! calculate bands
878 66 : NULLIFY (qs_env_kp)
879 66 : IF (kpoints_source == w90_kpoints_mp_grid .AND. input_kpoint_symmetry .AND. iw > 0) THEN
880 : WRITE (iw, '(T2,A)') &
881 0 : "WANNIER90| Atomic k-point symmetry from the SCF calculation is not reused."
882 : WRITE (iw, '(T2,A)') &
883 0 : "WANNIER90| A full Monkhorst-Pack grid is generated for the Wannier90 interface."
884 : END IF
885 66 : IF (do_kpoints) THEN
886 : ! we already do kpoints
887 66 : qs_env_kp => qs_env
888 : ELSE
889 : ! we start from gamma point only
890 0 : CALL create_kp_from_gamma(qs_env, qs_env_kp)
891 : END IF
892 66 : IF (iw > 0) THEN
893 33 : WRITE (unit=iw, FMT="(/,T2,A)") "Start K-Point Calculation ..."
894 : END IF
895 : ! The Gamma-to-k-point conversion normalizes the spin dimensions for diagonalization.
896 : ! Keep the original band selection, but use the converted MO sizes for its workspace.
897 66 : CALL get_qs_env(qs_env=qs_env_kp, para_env=para_env, blacs_env=blacs_env, mos=mos)
898 66 : CALL kpoint_env_initialize(kpoint, para_env, blacs_env)
899 66 : CALL kpoint_initialize_mos(kpoint, mos, nadd)
900 66 : CALL kpoint_initialize_mo_set(kpoint)
901 : !
902 66 : CALL get_qs_env(qs_env=qs_env_kp, sab_orb=sab_nl, dft_control=dft_control)
903 66 : CALL kpoint_init_cell_index(kpoint, sab_nl, para_env, dft_control%nimages)
904 : !
905 : CALL get_qs_env(qs_env=qs_env_kp, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s, &
906 66 : scf_env=scf_env, scf_control=scf_control)
907 66 : full_mesh_diagonalized = .FALSE.
908 66 : reused_scf_mos = .FALSE.
909 66 : reuse_reason = ""
910 66 : aligned_degenerate_blocks = 0
911 66 : aligned_degenerate_max_size = 0
912 66 : aligned_degenerate_min_svalue = 0.0_dp
913 66 : IF (reuse_scf_mos) THEN
914 42 : CALL get_kpoint_info(kpoint=kpoint, cell_to_index=cell_to_index)
915 : CALL do_general_diag_kp(matrix_ks, matrix_s, qs_kpoint, scf_env, scf_control, .FALSE., &
916 42 : diis_step)
917 42 : IF (validate_reuse_scf_mos) THEN
918 14 : IF (iw > 0) THEN
919 : WRITE (iw, '(T2,A)') &
920 7 : "WANNIER90| Validating SCF MO reuse against a full-mesh diagonalization reference."
921 : END IF
922 14 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoint, scf_env, scf_control, .FALSE., diis_step)
923 14 : full_mesh_diagonalized = .TRUE.
924 14 : CALL save_wannier90_mo_snapshot(kpoint, ispin, para_env, reference_mo, reference_eigenvalues)
925 : CALL diagnose_wannier90_scf_reuse_candidates(kpoint, qs_kpoint, matrix_s, matrix_ks, &
926 : cell_to_index, sab_nl, para_env, iw, ispin, &
927 : reuse_candidate_deviation, &
928 : reuse_candidate_min_svalue, &
929 : reuse_candidate_metric_deviation, &
930 14 : reuse_candidate_residual)
931 14 : IF (iw > 0 .AND. reuse_candidate_deviation < 1.0e100_dp) THEN
932 : WRITE (iw, '(T2,A,ES10.3,A,ES10.3,A,ES10.3,A,ES10.3)') &
933 7 : "WANNIER90| Best atom/AO candidate subspace deviation ", &
934 7 : reuse_candidate_deviation, ", minimum singular value ", &
935 7 : reuse_candidate_min_svalue, ", max metric deviation ", &
936 14 : reuse_candidate_metric_deviation, ", max residual ", reuse_candidate_residual
937 : END IF
938 : END IF
939 : CALL prepare_wannier90_scf_mos(kpoint, qs_kpoint, matrix_s, matrix_ks, cell_to_index, &
940 : sab_nl, para_env, reused_scf_mos, reuse_reason, &
941 : aligned_degenerate_blocks, aligned_degenerate_max_size, &
942 42 : aligned_degenerate_min_svalue, spin_channel=ispin)
943 42 : IF (validate_reuse_scf_mos) THEN
944 14 : IF (reused_scf_mos) THEN
945 : CALL validate_wannier90_reused_mos(kpoint, matrix_s, cell_to_index, sab_nl, &
946 : para_env, ispin, reference_mo, &
947 : reference_eigenvalues, validate_reuse_ok, &
948 : validation_subspace_deviation, validation_min_svalue, &
949 14 : validation_eigenvalue_deviation)
950 14 : IF (iw > 0) THEN
951 : WRITE (iw, '(T2,A,ES10.3,A,ES10.3,A,ES10.3)') &
952 7 : "WANNIER90| Reused MO validation: subspace deviation ", &
953 7 : validation_subspace_deviation, ", minimum singular value ", &
954 7 : validation_min_svalue, ", eigenvalue deviation ", &
955 14 : validation_eigenvalue_deviation
956 : END IF
957 14 : IF (.NOT. validate_reuse_ok) THEN
958 2 : reused_scf_mos = .FALSE.
959 : WRITE (reuse_reason, "(A,ES10.3,A,ES10.3)") &
960 2 : "validation failed: dS=", &
961 4 : validation_subspace_deviation, ", dE=", validation_eigenvalue_deviation
962 : END IF
963 : END IF
964 14 : IF (.NOT. reused_scf_mos) THEN
965 : CALL restore_wannier90_mo_snapshot(kpoint, ispin, reference_mo, &
966 2 : reference_eigenvalues)
967 : END IF
968 : END IF
969 42 : IF (iw > 0) THEN
970 21 : IF (reused_scf_mos) THEN
971 : WRITE (iw, '(T2,A)') &
972 17 : "WANNIER90| Reused SCF MO coefficients for the Wannier90 full k-point mesh."
973 17 : IF (use_bloch_phases) THEN
974 10 : IF (initial_projections == 0) WRITE (iw, '(T2,A)') &
975 6 : "WANNIER90| Wrote identity projections for Bloch-phase complete band subspaces."
976 : WRITE (iw, '(T2,A,3F10.6)') &
977 10 : "WANNIER90| Applied Bloch phase gauge to reused overlaps around fractional center", &
978 20 : phase_center(1:3)
979 : END IF
980 17 : IF (aligned_degenerate_blocks > 0) THEN
981 : WRITE (iw, '(T2,A,I0,A,I0,A,ES10.3)') &
982 10 : "WANNIER90| Ritz-stabilized ", aligned_degenerate_blocks, &
983 10 : " degenerate SCF MO subspace(s) with S(k),H(k); largest block has ", &
984 10 : aligned_degenerate_max_size, " band(s), min metric eigenvalue ", &
985 20 : aligned_degenerate_min_svalue
986 : END IF
987 : ELSE
988 : WRITE (iw, '(T2,A,A)') &
989 4 : "WANNIER90| Could not reuse SCF MOs: ", TRIM(reuse_reason)
990 : WRITE (iw, '(T2,A)') &
991 4 : "WANNIER90| Falling back to full-mesh diagonalization for the Wannier90 files."
992 : END IF
993 : END IF
994 : END IF
995 66 : IF (.NOT. reused_scf_mos .AND. .NOT. full_mesh_diagonalized) THEN
996 30 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoint, scf_env, scf_control, .FALSE., diis_step)
997 : END IF
998 66 : IF (ALLOCATED(reference_mo)) DEALLOCATE (reference_mo)
999 66 : IF (ALLOCATED(reference_eigenvalues)) DEALLOCATE (reference_eigenvalues)
1000 : !
1001 66 : IF (iw > 0) THEN
1002 33 : WRITE (iw, '(T69,A)') "... Finished"
1003 : END IF
1004 : !
1005 : ! Calculate and print Overlaps
1006 : !
1007 : ! create a list of unique b vectors and a table of pointers
1008 : ! nblist(ik,i) -> +/- b_latt(1:3,x)
1009 264 : ALLOCATE (nblist(num_kpts, nntot))
1010 198 : ALLOCATE (b_latt(3, num_kpts*nntot))
1011 66 : nblist(:, :) = 0
1012 66 : nbs = 0
1013 1002 : DO ik = 1, num_kpts
1014 4914 : DO i = 1, nntot
1015 15648 : bvec(1:3) = kpt_latt(1:3, nnlist(ik, i)) - kpt_latt(1:3, ik) + nncell(1:3, ik, i)
1016 7656 : ibs = 0
1017 7656 : DO k = 1, nbs
1018 29984 : IF (SUM(ABS(bvec(1:3) - b_latt(1:3, k))) < 1.e-6_dp) THEN
1019 : ibs = k
1020 : EXIT
1021 : END IF
1022 22240 : IF (SUM(ABS(bvec(1:3) + b_latt(1:3, k))) < 1.e-6_dp) THEN
1023 1776 : ibs = -k
1024 1776 : EXIT
1025 : END IF
1026 : END DO
1027 4848 : IF (ibs /= 0) THEN
1028 : ! old lattice vector
1029 3752 : nblist(ik, i) = ibs
1030 : ELSE
1031 : ! new lattice vector
1032 160 : nbs = nbs + 1
1033 640 : b_latt(1:3, nbs) = bvec(1:3)
1034 160 : nblist(ik, i) = nbs
1035 : END IF
1036 : END DO
1037 : END DO
1038 : ! calculate all the operator matrices (a|bvec|b)
1039 66 : overlap_nl => sab_nl
1040 66 : NULLIFY (berry_kpoint)
1041 66 : IF (ordered_berry) THEN
1042 22 : CALL get_qs_env(qs_env_kp, sab_all=overlap_nl)
1043 22 : CALL kpoint_create(berry_kpoint)
1044 22 : CALL kpoint_init_cell_index(berry_kpoint, overlap_nl, para_env, nberry_images)
1045 22 : CALL get_kpoint_info(berry_kpoint, cell_to_index=berry_cell_index)
1046 : END IF
1047 66 : IF (.NOT. ASSOCIATED(overlap_nl)) CPABORT("Explicit overlaps require k-point neighbour lists.")
1048 358 : ALLOCATE (berry_matrix(nbs))
1049 226 : DO i = 1, nbs
1050 160 : NULLIFY (berry_matrix(i)%cosmat)
1051 160 : NULLIFY (berry_matrix(i)%sinmat)
1052 800 : bvec(1:3) = -twopi*MATMUL(TRANSPOSE(cell%h_inv(1:3, 1:3)), b_latt(1:3, i))
1053 : CALL build_berry_kpoint_matrix(qs_env_kp, berry_matrix(i)%cosmat, &
1054 226 : berry_matrix(i)%sinmat, bvec, ordered=ordered_berry, ordered_kpoints=berry_kpoint)
1055 : END DO
1056 : ! work matrices for MOs (all group)
1057 66 : kp => kpoint%kp_env(1)%kpoint_env
1058 66 : CALL get_mo_set(kp%mos(ispin), nmo=nmo)
1059 66 : IF (ordered_berry .AND. nmo /= nscalar) THEN
1060 0 : CPABORT("WANNIER90: unexpected orbital count in export.")
1061 : END IF
1062 66 : NULLIFY (matrix_struct_ao, matrix_struct_work)
1063 : CALL cp_fm_struct_create(matrix_struct_work, nrow_global=nao, &
1064 : ncol_global=nmo, &
1065 : para_env=para_env, &
1066 66 : context=blacs_env)
1067 66 : CALL cp_cfm_create(fmk1_cfm, matrix_struct_work)
1068 66 : CALL cp_cfm_create(fmk2_cfm, matrix_struct_work)
1069 66 : CALL cp_cfm_create(tmp_cfm, matrix_struct_work)
1070 : CALL cp_fm_struct_create(matrix_struct_ao, nrow_global=nao, &
1071 : ncol_global=nao, &
1072 : para_env=para_env, &
1073 66 : context=blacs_env)
1074 66 : CALL cp_fm_create(mat_real, matrix_struct_ao)
1075 66 : CALL cp_fm_create(mat_imag, matrix_struct_ao)
1076 66 : CALL cp_cfm_create(omat_cfm, matrix_struct_ao)
1077 : ! work matrices for Mmn(k,b) integrals
1078 66 : NULLIFY (matrix_struct_mmn)
1079 : CALL cp_fm_struct_create(matrix_struct_mmn, nrow_global=nmo, &
1080 : ncol_global=nmo, &
1081 : para_env=para_env, &
1082 66 : context=blacs_env)
1083 66 : CALL cp_cfm_create(mmn_cfm, matrix_struct_mmn)
1084 : ! allocate some work matrices
1085 66 : ALLOCATE (rmatrix_full, cmatrix_full)
1086 : CALL dbcsr_create(rmatrix_full, template=matrix_s(1, 1)%matrix, &
1087 66 : matrix_type=dbcsr_type_no_symmetry)
1088 : CALL dbcsr_create(cmatrix_full, template=matrix_s(1, 1)%matrix, &
1089 66 : matrix_type=dbcsr_type_no_symmetry)
1090 66 : IF (ordered_berry) THEN
1091 22 : ALLOCATE (loop_real, loop_imag)
1092 22 : CALL dbcsr_create(loop_real, template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1093 22 : CALL dbcsr_create(loop_imag, template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1094 22 : CALL cp_dbcsr_alloc_block_from_nbl(loop_real, overlap_nl)
1095 22 : CALL cp_dbcsr_alloc_block_from_nbl(loop_imag, overlap_nl)
1096 22 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix_full, overlap_nl)
1097 22 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix_full, overlap_nl)
1098 : ELSE
1099 44 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix_full, sab_nl, desymmetrize=.TRUE.)
1100 44 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix_full, sab_nl, desymmetrize=.TRUE.)
1101 : END IF
1102 : !
1103 66 : CALL get_kpoint_info(kpoint=kpoint, cell_to_index=cell_to_index)
1104 66 : NULLIFY (cmo_coeff_source)
1105 66 : IF (do_soc .OR. do_parity) THEN
1106 : ! Parity needs the overlap work matrices even without SOC, and after SOC is done.
1107 12 : ALLOCATE (rmatrix, cmatrix)
1108 12 : CALL dbcsr_create(rmatrix, template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
1109 12 : CALL dbcsr_create(cmatrix, template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
1110 12 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
1111 12 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
1112 : END IF
1113 66 : IF (do_soc) THEN
1114 : ! Second variation in the scalar KS eigenbasis. Retain selected spinors at each k.
1115 30 : CPASSERT(ALL(kpoint%kp_range == [1, num_kpts]))
1116 10 : NULLIFY (soc_matrices)
1117 10 : CALL V_SOC_xyz_from_pseudopotential(qs_env_kp, soc_matrices)
1118 100 : ALLOCATE (soc_xyz(nmo, nmo, 3), soc_h(2*nmo, 2*nmo), soc_u(2*nmo, 2*nmo))
1119 60 : ALLOCATE (scalar_values(2*nmo), spinor_values(2*nmo, num_kpts))
1120 80 : ALLOCATE (spinor_coeff(2*nmo, num_bands, num_kpts), scalar_overlap(nmo, nmo))
1121 146 : DO ik = 1, num_kpts
1122 136 : kp => kpoint%kp_env(ik)%kpoint_env
1123 136 : CALL cp_cfm_copy_general(kp%mos(1)%cmo_coeff, fmk1_cfm, para_env)
1124 136 : CALL get_mo_set(kp%mos(1), eigenvalues=eigenvalues)
1125 984 : scalar_values(1:nmo) = eigenvalues(1:nmo)
1126 984 : scalar_values(nmo + 1:) = eigenvalues(1:nmo)
1127 544 : DO axis = 1, 3
1128 408 : CALL dbcsr_set(rmatrix, 0.0_dp)
1129 408 : CALL dbcsr_set(cmatrix, 0.0_dp)
1130 : CALL rskp_transform(cmatrix, rmatrix, rsmat=soc_matrices, ispin=axis, &
1131 408 : xkp=kpoint%xkp(:, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1132 408 : CALL dbcsr_scale(rmatrix, -1.0_dp)
1133 408 : CALL dbcsr_desymmetrize(rmatrix, rmatrix_full)
1134 408 : CALL dbcsr_desymmetrize(cmatrix, cmatrix_full)
1135 408 : CALL copy_dbcsr_to_fm(rmatrix_full, mat_real)
1136 408 : CALL copy_dbcsr_to_fm(cmatrix_full, mat_imag)
1137 408 : CALL cp_fm_to_cfm(mat_real, mat_imag, omat_cfm)
1138 : CALL parallel_gemm("N", "N", nao, nmo, nao, CMPLX(1.0_dp, 0.0_dp, dp), &
1139 408 : omat_cfm, fmk1_cfm, CMPLX(0.0_dp, 0.0_dp, dp), tmp_cfm)
1140 : CALL parallel_gemm("C", "N", nmo, nmo, nao, CMPLX(1.0_dp, 0.0_dp, dp), &
1141 408 : fmk1_cfm, tmp_cfm, CMPLX(0.0_dp, 0.0_dp, dp), mmn_cfm)
1142 544 : CALL cp_cfm_get_submatrix(mmn_cfm, soc_xyz(:, :, axis))
1143 : END DO
1144 9592 : soc_h(1:nmo, 1:nmo) = soc_xyz(:, :, 3)
1145 9592 : soc_h(nmo + 1:, nmo + 1:) = -soc_xyz(:, :, 3)
1146 : ! Match H_KS_spinor_kp and CP2K's real-spherical angular-momentum convention.
1147 9592 : soc_h(1:nmo, nmo + 1:) = soc_xyz(:, :, 1) + CMPLX(0.0_dp, 1.0_dp, dp)*soc_xyz(:, :, 2)
1148 9592 : soc_h(nmo + 1:, 1:nmo) = soc_xyz(:, :, 1) - CMPLX(0.0_dp, 1.0_dp, dp)*soc_xyz(:, :, 2)
1149 1832 : DO ib = 1, 2*nmo
1150 1832 : soc_h(ib, ib) = soc_h(ib, ib) + scalar_values(ib)
1151 : END DO
1152 36264 : IF (MAXVAL(ABS(soc_h - CONJG(TRANSPOSE(soc_h)))) > 1.e-8_dp) THEN
1153 0 : IF (iw > 0) WRITE (iw, '(T2,A,I0,A,ES18.10)') &
1154 0 : "TOPOLOGY| SOC Hermiticity error at point ", ik, ": ", &
1155 0 : MAXVAL(ABS(soc_h - CONJG(TRANSPOSE(soc_h))))
1156 0 : CPABORT("WANNIER90 SOC: non-Hermitian second-variational Hamiltonian.")
1157 : END IF
1158 136 : CALL diag_complex(soc_h, soc_u, spinor_values(:, ik))
1159 14802 : spinor_coeff(:, :, ik) = soc_u(:, band_indices)
1160 : END DO
1161 10 : CALL dbcsr_deallocate_matrix_set(soc_matrices)
1162 10 : DEALLOCATE (soc_xyz, soc_h, soc_u, scalar_values)
1163 10 : IF (iw > 0) WRITE (iw, '(T2,A,I0)') "TOPOLOGY| Scalar bands in SOC second variation: ", nmo
1164 : END IF
1165 : ! Eigenvalues determine the projection windows and the sampled spectral gaps.
1166 66 : kp => kpoint%kp_env(1)%kpoint_env
1167 66 : CALL get_mo_set(kp%mos(ispin), nmo=nmo)
1168 66 : CPASSERT(nmo >= nscalar)
1169 254 : ALLOCATE (eigval(MERGE(2*nmo, nmo, do_soc)))
1170 66 : direct_gap = HUGE(1.0_dp)
1171 66 : valence_max = -HUGE(1.0_dp)
1172 66 : conduction_min = HUGE(1.0_dp)
1173 66 : CALL get_kpoint_info(kpoint, nkp=nkp, kp_range=kp_range, xkp=xkp)
1174 66 : IF (para_env%is_source()) THEN
1175 33 : WRITE (filename, '(A,A)') TRIM(seed_name), ".eig"
1176 33 : CALL open_file(filename, unit_number=iunit, file_status="REPLACE", file_action="WRITE")
1177 : ELSE
1178 33 : iunit = -1
1179 : END IF
1180 : !
1181 1002 : DO ik = 1, nkp
1182 936 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
1183 936 : IF (do_soc) THEN
1184 1832 : eigval(:) = spinor_values(:, ik)
1185 800 : ELSE IF (my_kpgrp) THEN
1186 800 : ikpgr = ik - kp_range(1) + 1
1187 800 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1188 800 : CALL get_mo_set(kp%mos(ispin), eigenvalues=eigenvalues)
1189 3296 : eigval(1:nmo) = eigenvalues(1:nmo)
1190 : ELSE
1191 0 : eigval(1:nmo) = 0.0_dp
1192 : END IF
1193 936 : IF (.NOT. do_soc) CALL kpoint%para_env_inter_kp%sum(eigval)
1194 936 : IF (.NOT. ordered_berry .AND. nexcl > 0) THEN
1195 80 : DO ib = 1, num_bands_tot - 1
1196 64 : IF (keep_band(ib) .EQV. keep_band(ib + 1)) CYCLE
1197 : ! Degenerate states outside the outer window cannot enter the Wannier subspace.
1198 32 : kept_index = MERGE(ib, ib + 1, keep_band(ib))
1199 32 : IF (eigval(kept_index) < energy_window(1) .OR. eigval(kept_index) > energy_window(2)) CYCLE
1200 16 : IF (ABS(eigval(ib + 1) - eigval(ib)) < 1.0e-8_dp) THEN
1201 0 : CPABORT("WANNIER90 EXCLUDE_BANDS must retain or remove complete degenerate subspaces.")
1202 : END IF
1203 : END DO
1204 : END IF
1205 936 : IF (do_wilson) THEN
1206 328 : IF (lowest_bands) THEN
1207 328 : valence_max = MAX(valence_max, eigval(num_bands))
1208 328 : conduction_min = MIN(conduction_min, eigval(num_bands + 1))
1209 : END IF
1210 1488 : DO ib = 1, num_bands_tot - 1
1211 10032 : IF (ANY(band_indices == ib) .NEQV. ANY(band_indices == ib + 1)) THEN
1212 328 : pair_gap = eigval(ib + 1) - eigval(ib)
1213 328 : direct_gap = MIN(direct_gap, pair_gap)
1214 : END IF
1215 : END DO
1216 : END IF
1217 5128 : eigval(:) = eigval*evolt
1218 4616 : IF (ANY(window_set)) THEN
1219 80 : projection_eigenvalues(:, ik) = eigval(band_indices)
1220 : CALL wannier90_window_mask(projection_eigenvalues(:, ik), num_wann, energy_window*evolt, &
1221 80 : window_set, allowed_band)
1222 : END IF
1223 : ! output
1224 1002 : IF (iunit > 0) THEN
1225 2140 : DO ib = 1, num_bands
1226 2140 : WRITE (iunit, "(2I8,F24.14)") ib, ik, eigval(band_indices(ib))
1227 : END DO
1228 : END IF
1229 : END DO
1230 66 : IF (iunit > 0) THEN
1231 33 : CALL close_file(iunit)
1232 : END IF
1233 66 : IF (.NOT. ordered_berry .AND. initial_projections >= 0) THEN
1234 120 : ALLOCATE (initial_amn(num_bands, num_wann, num_kpts))
1235 24 : IF (initial_projections == 0) THEN
1236 14 : initial_amn(:, :, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1237 62 : DO i = 1, num_wann
1238 1342 : initial_amn(i, i, :) = CMPLX(1.0_dp, 0.0_dp, KIND=dp)
1239 : END DO
1240 : ELSE
1241 : CALL build_wannier90_ao_projections(kpoint, matrix_s, sab_nl, para_env, ispin, band_indices, &
1242 12 : reused_scf_mos .AND. use_bloch_phases, phase_center, ao_projections)
1243 42 : IF (para_env%is_source() .AND. ANY(window_set)) THEN
1244 9 : DO ik = 1, num_kpts
1245 : CALL wannier90_window_mask(projection_eigenvalues(:, ik), num_wann, energy_window*evolt, &
1246 40 : window_set, allowed_band)
1247 41 : DO ib = 1, num_bands
1248 160 : IF (.NOT. allowed_band(ib)) ao_projections(ib, :, ik) = CMPLX(0.0_dp, 0.0_dp, dp)
1249 : END DO
1250 : END DO
1251 : END IF
1252 10 : IF (initial_projections == 2) THEN
1253 4 : CALL build_wannier90_hybrid_groups(qs_env, hybrid_groups)
1254 4 : CALL wannier90_hybrid_candidates(ao_projections, hybrid_groups, hybrid_projections)
1255 4 : DEALLOCATE (ao_projections, hybrid_groups)
1256 4 : CALL MOVE_ALLOC(hybrid_projections, ao_projections)
1257 4 : IF (iw > 0) WRITE (iw, '(T2,A)') "WANNIER90| Selecting from the full radial s/p hybrid candidate bank."
1258 : END IF
1259 30 : ALLOCATE (projection_pivots(num_wann))
1260 : CALL wannier90_select_projections(ao_projections, para_env, initial_amn, projection_pivots, &
1261 10 : projection_min_svalue)
1262 10 : IF (iw > 0) THEN
1263 5 : WRITE (iw, '(T2,A,ES12.4)') "WANNIER90| Minimum projection singular value ", projection_min_svalue
1264 5 : WRITE (iw, '(T2,A,*(I6))') "WANNIER90| Selected projection candidates ", projection_pivots
1265 : END IF
1266 10 : IF (projection_min_svalue < 1.0e-10_dp) THEN
1267 0 : CPABORT("No full-rank Wannier90 AO trial set was found on the entire k-point mesh.")
1268 : END IF
1269 10 : DEALLOCATE (projection_pivots)
1270 20 : DEALLOCATE (ao_projections)
1271 : END IF
1272 24 : CALL write_wannier90_amn(initial_amn, para_env, seed_name)
1273 24 : DEALLOCATE (initial_amn)
1274 : END IF
1275 66 : IF (para_env%is_source()) THEN
1276 33 : WRITE (filename, '(A,A)') TRIM(seed_name), ".mmn"
1277 33 : CALL open_file(filename, unit_number=iunit, file_status="REPLACE", file_action="WRITE")
1278 33 : CALL m_timestamp(timestamp)
1279 33 : WRITE (iunit, "(A)") "! Wannier90 file generated by CP2K "//timestamp
1280 33 : WRITE (iunit, "(3I8)") num_bands, num_kpts, nntot
1281 : ELSE
1282 33 : iunit = -1
1283 : END IF
1284 : ! loop over all k-points
1285 1002 : DO ik = 1, num_kpts
1286 : ! get the MO coefficients for this k-point
1287 936 : my_kpgrp = (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2))
1288 : IF (my_kpgrp) THEN
1289 936 : ikk = ik - kpoint%kp_range(1) + 1
1290 936 : kp => kpoint%kp_env(ikk)%kpoint_env
1291 936 : cmo_coeff_source => kp%mos(ispin)%cmo_coeff
1292 : ELSE
1293 0 : NULLIFY (cmo_coeff_source, kp)
1294 : END IF
1295 936 : CALL wannier90_copy_cfm_general(cmo_coeff_source, fmk1_cfm, para_env)
1296 : ! loop over all connected neighbors
1297 4914 : DO i = 1, nntot
1298 : ! get the MO coefficients for the connected k-point
1299 3912 : ik2 = nnlist(ik, i)
1300 3912 : mygrp = (ik2 >= kpoint%kp_range(1) .AND. ik2 <= kpoint%kp_range(2))
1301 : IF (mygrp) THEN
1302 3912 : ikk = ik2 - kpoint%kp_range(1) + 1
1303 3912 : kp => kpoint%kp_env(ikk)%kpoint_env
1304 3912 : cmo_coeff_source => kp%mos(ispin)%cmo_coeff
1305 : ELSE
1306 0 : NULLIFY (cmo_coeff_source, kp)
1307 : END IF
1308 3912 : CALL wannier90_copy_cfm_general(cmo_coeff_source, fmk2_cfm, para_env)
1309 : !
1310 : ! transfer realspace overlaps to connected k-point
1311 3912 : ibs = nblist(ik, i)
1312 3912 : ksign = SIGN(1.0_dp, REAL(ibs, KIND=dp))
1313 3912 : ibs = ABS(ibs)
1314 3912 : IF (ordered_berry) THEN
1315 : ! The cross-k AO operator is not Hermitian. Retain ordered atom pairs
1316 : ! and use exp(i*k'.R) [cos(b.r) - i*sin(b.r)] without symmetry shortcuts.
1317 360 : CALL dbcsr_set(rmatrix_full, 0.0_dp)
1318 360 : CALL dbcsr_set(cmatrix_full, 0.0_dp)
1319 360 : CALL dbcsr_set(loop_real, 0.0_dp)
1320 360 : CALL dbcsr_set(loop_imag, 0.0_dp)
1321 : CALL rskp_transform(rmatrix_full, cmatrix_full, berry_matrix(ibs)%cosmat, 1, &
1322 360 : kpoint%xkp(:, ik2), berry_cell_index, overlap_nl)
1323 : CALL rskp_transform(loop_real, loop_imag, berry_matrix(ibs)%sinmat, 1, &
1324 360 : kpoint%xkp(:, ik2), berry_cell_index, overlap_nl, rs_sign=ksign)
1325 360 : CALL dbcsr_add(rmatrix_full, loop_imag, 1.0_dp, -1.0_dp)
1326 360 : CALL dbcsr_add(cmatrix_full, loop_real, 1.0_dp, 1.0_dp)
1327 : ELSE
1328 : CALL transform_wannier90_overlap(berry_matrix(ibs)%cosmat, berry_matrix(ibs)%sinmat, &
1329 : ksign, kpoint%xkp(:, ik), kpoint%xkp(:, ik2), &
1330 3552 : cell_to_index, sab_nl, rmatrix_full, cmatrix_full)
1331 : END IF
1332 : !
1333 : ! calculate M_(mn)^(k,b) = C(k)^H O(k,b) C(k+b)
1334 3912 : CALL copy_dbcsr_to_fm(rmatrix_full, mat_real)
1335 3912 : CALL copy_dbcsr_to_fm(cmatrix_full, mat_imag)
1336 3912 : CALL cp_fm_to_cfm(mat_real, mat_imag, omat_cfm)
1337 : CALL parallel_gemm("N", "N", nao, nmo, nao, CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
1338 3912 : omat_cfm, fmk2_cfm, CMPLX(0.0_dp, 0.0_dp, KIND=dp), tmp_cfm)
1339 : CALL parallel_gemm("C", "N", nmo, nmo, nao, CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
1340 3912 : fmk1_cfm, tmp_cfm, CMPLX(0.0_dp, 0.0_dp, KIND=dp), mmn_cfm)
1341 3912 : IF (do_soc) THEN
1342 136 : CALL cp_cfm_get_submatrix(mmn_cfm, scalar_overlap)
1343 136 : link_matrix(:, :) = MATMUL(CONJG(TRANSPOSE(spinor_coeff(1:nmo, :, ik))), &
1344 544 : MATMUL(scalar_overlap, spinor_coeff(1:nmo, :, ik2))) + &
1345 544 : MATMUL(CONJG(TRANSPOSE(spinor_coeff(nmo + 1:, :, ik))), &
1346 922480 : MATMUL(scalar_overlap, spinor_coeff(nmo + 1:, :, ik2)))
1347 : END IF
1348 : !
1349 : ! write to output file
1350 3912 : IF (reused_scf_mos .AND. use_bloch_phases) THEN
1351 : ! Reused SCF MOs need the same global Bloch gauge in every overlap block.
1352 : gauge_arg = twopi*DOT_PRODUCT(kpoint%xkp(1:3, ik2) - kpoint%xkp(1:3, ik), &
1353 9344 : phase_center(1:3))
1354 2336 : gauge_real = COS(gauge_arg)
1355 2336 : gauge_imag = SIN(gauge_arg)
1356 : ELSE
1357 : gauge_real = 1.0_dp
1358 : gauge_imag = 0.0_dp
1359 : END IF
1360 3912 : IF (iunit > 0) THEN
1361 1956 : WRITE (iunit, "(2I8,3I5)") ik, ik2, nncell(1:3, ik, i)
1362 : END IF
1363 17800 : DO ib2 = 1, num_bands
1364 77704 : DO ib1 = 1, num_bands
1365 59904 : IF (do_soc) THEN
1366 8704 : rmmn = REAL(link_matrix(ib1, ib2), dp)
1367 8704 : cmmn = AIMAG(link_matrix(ib1, ib2))
1368 : ELSE
1369 51200 : CALL cp_cfm_get_element(mmn_cfm, band_indices(ib1), band_indices(ib2), zmmn)
1370 51200 : rmmn = REAL(zmmn, KIND=dp)
1371 51200 : cmmn = AIMAG(zmmn)
1372 : END IF
1373 59904 : gauge_tmp = gauge_real*rmmn - gauge_imag*cmmn
1374 59904 : cmmn = gauge_imag*rmmn + gauge_real*cmmn
1375 59904 : rmmn = gauge_tmp
1376 59904 : link_matrix(ib1, ib2) = CMPLX(rmmn, cmmn, dp)
1377 73792 : IF (iunit > 0) THEN
1378 29952 : WRITE (iunit, "(2E30.14)") rmmn, cmmn
1379 : END IF
1380 : END DO
1381 : END DO
1382 4848 : IF (do_wilson) THEN
1383 328 : iloop = loop_index(ik)
1384 328 : CALL wilson_step(wilson_product(:, :, iloop), link_matrix, link_sv, status, 1.e-10_dp)
1385 328 : IF (status /= 0) THEN
1386 0 : CPABORT("Wilson link singular or SVD failed; refine sampling and check subspace.")
1387 : END IF
1388 328 : loop_sv(iloop) = MIN(loop_sv(iloop), link_sv)
1389 : END IF
1390 : !
1391 : END DO
1392 : END DO
1393 : ! Optional full-state snapshots, before the distributed work matrices are released.
1394 66 : IF (export_state .OR. do_parity) THEN
1395 12 : state_components = 1
1396 12 : IF (do_soc) state_components = 2
1397 84 : ALLOCATE (export_scalar(nao, nscalar), export_coeff(nao*state_components, num_bands))
1398 12 : state_unit = -1
1399 12 : IF (export_state) CALL topology_state_begin(qs_env_kp, TRIM(seed_name)//".topology", nao, num_bands, num_kpts, &
1400 4 : state_components, num_bands_tot, spin_channel, band_indices, state_unit)
1401 12 : IF (do_parity) THEN
1402 24 : CPASSERT(ALL(kpoint%kp_range == [1, num_kpts]))
1403 0 : ALLOCATE (parity_metric(nao, nao), parity_phase(nao), parity_map(nao), &
1404 0 : parity_odd(2**tqc_dimension), ebr_coefficients(2, 2**tqc_dimension), &
1405 112 : parity_energies(num_bands_tot))
1406 56 : parity_odd = -1
1407 8 : IF (do_tqc) THEN
1408 102 : IF (num_bands /= nelectron .OR. MOD(num_bands, 2) /= 0 .OR. &
1409 : ANY(band_indices /= [(i, i=1, num_bands)])) THEN
1410 0 : CPABORT("INVERSION_TQC requires the full occupied spinor band prefix.")
1411 : END IF
1412 : END IF
1413 : END IF
1414 156 : DO ik = 1, num_kpts
1415 144 : kp => kpoint%kp_env(ik)%kpoint_env
1416 144 : CALL cp_cfm_copy_general(kp%mos(spin_channel)%cmo_coeff, fmk1_cfm, para_env)
1417 144 : CALL cp_cfm_get_submatrix(fmk1_cfm, export_scalar)
1418 144 : IF (do_soc) THEN
1419 230304 : export_coeff(:nao, :) = MATMUL(export_scalar, spinor_coeff(:nscalar, :, ik))
1420 174560 : export_coeff(nao + 1:, :) = MATMUL(export_scalar, spinor_coeff(nscalar + 1:, :, ik))
1421 32 : CALL topology_state_point(state_unit, ik, kpt_latt(:, ik), spinor_values(:, ik), export_coeff)
1422 688 : IF (do_parity) parity_energies(:) = spinor_values(:, ik)
1423 : ELSE
1424 784 : export_coeff(:, :) = export_scalar(:, band_indices)
1425 112 : CALL get_mo_set(kp%mos(spin_channel), eigenvalues=eigenvalues)
1426 112 : CALL topology_state_point(state_unit, ik, kpt_latt(:, ik), eigenvalues(:nscalar), export_coeff)
1427 128 : IF (do_parity) parity_energies(:) = eigenvalues(:nscalar)
1428 : END IF
1429 144 : IF (.NOT. do_parity) CYCLE
1430 40 : parity_gap = HUGE(1.0_dp)
1431 712 : DO jpar = 1, num_bands_tot
1432 4776 : IF (ANY(band_indices == jpar)) CYCLE
1433 3920 : parity_gap = MIN(parity_gap, MINVAL(ABS(parity_energies(band_indices) - parity_energies(jpar))))
1434 : END DO
1435 40 : IF (num_bands == num_bands_tot .OR. parity_gap <= gap_tol) THEN
1436 0 : CPABORT("PARITY requires an isolated selected subspace and computed excluded bands.")
1437 : END IF
1438 160 : IF (MAXVAL(ABS(2*kpt_latt(:, ik) - NINT(2*kpt_latt(:, ik)))) > 1.e-8_dp) CYCLE
1439 : CALL gaussian_inversion_action(qs_env_kp, kpt_latt(:, ik), parity_origin, parity_map, &
1440 36 : parity_phase, parity_tolerance, status)
1441 36 : IF (status /= 0) THEN
1442 0 : CPABORT("PARITY: inversion does not preserve geometry and atomic kinds.")
1443 : END IF
1444 36 : CALL dbcsr_set(rmatrix, 0.0_dp)
1445 36 : CALL dbcsr_set(cmatrix, 0.0_dp)
1446 : CALL rskp_transform(rmatrix, cmatrix, rsmat=matrix_s, ispin=1, xkp=kpoint%xkp(:, ik), &
1447 36 : cell_to_index=cell_to_index, sab_nl=sab_nl)
1448 36 : CALL dbcsr_desymmetrize(rmatrix, rmatrix_full)
1449 36 : CALL dbcsr_desymmetrize(cmatrix, cmatrix_full)
1450 36 : CALL copy_dbcsr_to_fm(rmatrix_full, mat_real)
1451 36 : CALL copy_dbcsr_to_fm(cmatrix_full, mat_imag)
1452 36 : CALL cp_fm_to_cfm(mat_real, mat_imag, omat_cfm)
1453 36 : CALL cp_cfm_get_submatrix(omat_cfm, parity_metric)
1454 : CALL inversion_representation(parity_metric, export_coeff, parity_map, parity_phase, &
1455 : parity_energies(band_indices), do_soc .AND. time_reversal, parity_tolerance, &
1456 300 : parity_energy_tolerance, counts, parity_error, parity_energy_error, status, parity_checks)
1457 36 : IF (iw > 0) THEN
1458 18 : WRITE (iw, '(T2,A,3F10.5,A,2I8,A,3ES14.5)') 'TQC| TRIM ', kpt_latt(:, ik), &
1459 36 : ' even/odd states ', counts, ' residuals/gap ', parity_error, parity_energy_error, parity_gap
1460 18 : WRITE (iw, '(T2,A,4ES14.5)') 'TQC| Metric/normalization/inversion/TR residuals: ', parity_checks
1461 : END IF
1462 36 : IF (status /= 0) THEN
1463 0 : CPABORT("PARITY: metric, subspace or symmetry representation check failed.")
1464 : END IF
1465 36 : IF (tqc_dimension == 2 .AND. MODULO(NINT(2*kpt_latt(3, ik)), 2) /= 0) CYCLE
1466 36 : trim_id = 1
1467 128 : DO axis = 1, tqc_dimension
1468 128 : trim_id = trim_id + 2**(axis - 1)*MODULO(NINT(2*kpt_latt(axis, ik)), 2)
1469 : END DO
1470 36 : IF (parity_odd(trim_id) >= 0 .AND. parity_odd(trim_id) /= counts(2)) THEN
1471 0 : CPABORT("PARITY: repeated equivalent TRIM have inconsistent multiplicities.")
1472 : END IF
1473 156 : parity_odd(trim_id) = counts(2)
1474 : END DO
1475 12 : IF (state_unit > 0) CALL close_file(state_unit)
1476 12 : IF (do_parity) THEN
1477 8 : IF (ALL(parity_odd < 0)) THEN
1478 0 : CPABORT("PARITY requires a TRIM in the requested dimension.")
1479 : END IF
1480 : END IF
1481 12 : IF (do_tqc) THEN
1482 38 : IF (ANY(parity_odd < 0)) THEN
1483 0 : CPABORT("INVERSION_TQC requires all four (2D) or eight (3D) TRIM.")
1484 : END IF
1485 38 : IF (ANY(MOD(parity_odd, 2) /= 0)) THEN
1486 0 : CPABORT("INVERSION_TQC requires whole odd-parity Kramers pairs.")
1487 : END IF
1488 38 : parity_odd(:) = parity_odd/2
1489 6 : CALL inversion_indicators(parity_odd, num_bands/2, tqc_dimension, strong, weak, z4, status)
1490 6 : IF (status /= 0) CPABORT("Invalid inversion indicator data.")
1491 : CALL inversion_ebr_signature(parity_odd, num_bands/2, signed_atomic, nonnegative_atomic, &
1492 6 : ebr_coefficients, status)
1493 6 : IF (status /= 0) CPABORT("Invalid inversion EBR data.")
1494 6 : IF (iw > 0) THEN
1495 3 : WRITE (iw, '(T2,A,I0)') 'TQC| Inversion-subgroup dimension: ', tqc_dimension
1496 3 : WRITE (iw, '(T2,A,I0)') 'TQC| Fu-Kane parity index: ', strong
1497 3 : IF (tqc_dimension == 3) THEN
1498 1 : WRITE (iw, '(T2,A,3I3)') 'TQC| Weak parity indices: ', weak
1499 1 : WRITE (iw, '(T2,A,I0)') 'TQC| Inversion Z4 (sum odd pairs mod 4): ', z4
1500 : END IF
1501 3 : WRITE (iw, '(T2,A,L1)') 'TQC| Signed atomic signature: ', signed_atomic
1502 3 : WRITE (iw, '(T2,A,L1)') 'TQC| Nonnegative atomic signature: ', nonnegative_atomic
1503 3 : IF (nonnegative_atomic) THEN
1504 14 : DO i = 1, SIZE(parity_odd)
1505 12 : WRITE (iw, '(T2,A,I0,A,2I8)') 'TQC| Center bit index ', i - 1, &
1506 26 : ' even/odd EBR pairs ', ebr_coefficients(:, i)
1507 : END DO
1508 2 : WRITE (iw, '(T2,A)') 'TQC| Atomic-compatible symmetry signature; not a proof of trivial topology.'
1509 1 : ELSE IF (signed_atomic) THEN
1510 : WRITE (iw, '(T2,A)') 'TQC| Nonnegative EBR obstruction; '// &
1511 0 : 'exclude hidden stable topology before calling it fragile.'
1512 : ELSE
1513 1 : WRITE (iw, '(T2,A)') 'TQC| Stable inversion-symmetry indicator obstruction.'
1514 : END IF
1515 3 : WRITE (iw, '(T2,A)') 'TQC| TRIM checks do not establish a bulk gap; verify band isolation over the BZ.'
1516 3 : WRITE (iw, '(T2,A)') 'TQC| Converge the scalar-state basis used for second-variational SOC.'
1517 : END IF
1518 : END IF
1519 : END IF
1520 226 : DO i = 1, nbs
1521 160 : CALL dbcsr_deallocate_matrix_set(berry_matrix(i)%cosmat)
1522 226 : CALL dbcsr_deallocate_matrix_set(berry_matrix(i)%sinmat)
1523 : END DO
1524 66 : DEALLOCATE (berry_matrix)
1525 66 : CALL cp_fm_struct_release(matrix_struct_work)
1526 66 : CALL cp_cfm_release(fmk1_cfm)
1527 66 : CALL cp_cfm_release(fmk2_cfm)
1528 66 : CALL cp_cfm_release(tmp_cfm)
1529 66 : CALL cp_fm_struct_release(matrix_struct_ao)
1530 66 : CALL cp_fm_release(mat_real)
1531 66 : CALL cp_fm_release(mat_imag)
1532 66 : CALL cp_cfm_release(omat_cfm)
1533 66 : CALL cp_fm_struct_release(matrix_struct_mmn)
1534 66 : CALL cp_cfm_release(mmn_cfm)
1535 66 : CALL dbcsr_deallocate_matrix(rmatrix_full)
1536 66 : CALL dbcsr_deallocate_matrix(cmatrix_full)
1537 66 : IF (do_soc .OR. do_parity) THEN
1538 12 : CALL dbcsr_deallocate_matrix(rmatrix)
1539 12 : CALL dbcsr_deallocate_matrix(cmatrix)
1540 : END IF
1541 66 : IF (ordered_berry) THEN
1542 22 : CALL dbcsr_deallocate_matrix(loop_real)
1543 22 : CALL dbcsr_deallocate_matrix(loop_imag)
1544 : END IF
1545 : !
1546 66 : IF (iunit > 0) THEN
1547 33 : CALL close_file(iunit)
1548 : END IF
1549 : !
1550 66 : IF (do_wilson) THEN
1551 16 : IF (iw > 0 .AND. num_bands < num_bands_tot) WRITE (iw, '(T2,A,ES18.10)') &
1552 8 : "TOPOLOGY| Minimum sampled subspace gap [eV]: ", direct_gap*evolt
1553 16 : IF (direct_gap < gap_tol) CPABORT("Wilson subspace is not isolated on the sampled points.")
1554 16 : IF (lowest_bands) THEN
1555 16 : IF (iw > 0) THEN
1556 : WRITE (iw, '(T2,A,ES18.10)') &
1557 8 : "TOPOLOGY| Maximum selected-band energy [eV]: ", valence_max*evolt
1558 : WRITE (iw, '(T2,A,ES18.10)') &
1559 8 : "TOPOLOGY| Minimum excluded-band energy [eV]: ", conduction_min*evolt
1560 : WRITE (iw, '(T2,A,ES18.10)') &
1561 8 : "TOPOLOGY| Sampled indirect gap [eV]: ", (conduction_min - valence_max)*evolt
1562 8 : IF (conduction_min - valence_max <= gap_tol) WRITE (iw, '(T2,A)') &
1563 0 : "TOPOLOGY| Isolated band subspace, but no resolved common spectral gap on these samples."
1564 : END IF
1565 16 : IF (require_global_gap .AND. conduction_min - valence_max <= gap_tol) THEN
1566 0 : CPABORT("No positive sampled indirect gap above the selected bands.")
1567 : END IF
1568 : END IF
1569 68 : DO iloop = 1, nloop
1570 52 : CALL wilson_spectrum(wilson_product(:, :, iloop), wcc_out(:, iloop), berry_phase, status)
1571 52 : IF (status /= 0) CPABORT("Wilson eigenvalue calculation failed.")
1572 52 : IF (iw > 0) WRITE (iw, '(T2,A,I0,A,F18.12,A,ES12.4)') &
1573 42 : "TOPOLOGY| Loop ", iloop, " Berry phase [rad]: ", berry_phase, " min singular value: ", loop_sv(iloop)
1574 : END DO
1575 16 : IF (para_env%is_source()) THEN
1576 8 : CALL open_file(TRIM(seed_name)//".wilson", unit_number=iunit, file_status="UNKNOWN", file_action="WRITE")
1577 8 : WRITE (iunit, '(A)') "# loop, minimum link singular value, sorted WCC (arg(lambda)/(2*pi))"
1578 34 : DO iloop = 1, nloop
1579 34 : WRITE (iunit, '(I8,*(1X,ES24.16))') iloop, loop_sv(iloop), wcc_out(:, iloop)
1580 : END DO
1581 8 : CALL close_file(iunit)
1582 : END IF
1583 16 : IF (do_z2) THEN
1584 4 : CALL z2_from_wcc(wcc_out, z2_value, status, 1.e-5_dp)
1585 4 : IF (status /= 0 .AND. iw > 0) WRITE (iw, '(T2,A,I0)') &
1586 0 : "TOPOLOGY| Kramers/crossing check requests refinement, status: ", status
1587 4 : IF (iw > 0) WRITE (iw, '(T2,A,I0)') "TOPOLOGY| Z2 candidate (before sampling convergence): ", z2_value
1588 : END IF
1589 16 : IF (do_chern) THEN
1590 4 : CALL chern_from_wcc(wcc_out, chern_value, chern_winding, status, 1.e-5_dp)
1591 4 : IF (status /= 0) THEN
1592 0 : chern_value = HUGE(0)
1593 0 : IF (iw > 0) WRITE (iw, '(T2,A,I0)') &
1594 0 : "TOPOLOGY| Chern closure/winding check requests refinement, status: ", status
1595 4 : ELSE IF (iw > 0) THEN
1596 : WRITE (iw, '(T2,A,I0,A,F18.12)') &
1597 2 : "TOPOLOGY| First Chern candidate: ", chern_value, ", winding: ", chern_winding
1598 : END IF
1599 : END IF
1600 : END IF
1601 : ! clean up
1602 66 : IF (ordered_berry) CALL kpoint_release(berry_kpoint)
1603 66 : DEALLOCATE (kpt_latt, atoms_cart, atom_symbols, eigval)
1604 66 : DEALLOCATE (nnlist, nncell)
1605 66 : DEALLOCATE (nblist, b_latt)
1606 66 : IF (para_env%is_source()) THEN
1607 : ! Also replace an earlier exclusion map when the seed is reused without exclusions.
1608 : CALL open_file(TRIM(seed_name)//'_band_indices.dat', unit_number=iunit, &
1609 33 : file_status='REPLACE', file_action='WRITE')
1610 33 : WRITE (iunit, '(A)') '# Exported band, original MO index (before EXCLUDE_BANDS)'
1611 150 : DO ib = 1, num_bands
1612 150 : WRITE (iunit, '(2I10)') ib, band_indices(ib)
1613 : END DO
1614 33 : CALL close_file(iunit)
1615 : END IF
1616 66 : IF (nexcl > 0) DEALLOCATE (exclude_bands)
1617 66 : DEALLOCATE (band_indices, keep_band, allowed_band)
1618 66 : IF (ALLOCATED(projection_eigenvalues)) DEALLOCATE (projection_eigenvalues)
1619 66 : IF (do_kpoints) THEN
1620 66 : NULLIFY (qs_env_kp)
1621 : ELSE
1622 0 : CALL qs_env_release(qs_env_kp)
1623 0 : DEALLOCATE (qs_env_kp)
1624 : NULLIFY (qs_env_kp)
1625 : END IF
1626 :
1627 66 : CALL kpoint_release(kpoint)
1628 :
1629 1056 : END SUBROUTINE wannier90_files
1630 :
1631 : ! **************************************************************************************************
1632 : !> \brief Reconstruct a full Wannier90 k-point MO set from the SCF k-point MOs.
1633 : !> \param kpoint full Wannier90 export k-point object
1634 : !> \param qs_kpoint SCF k-point object
1635 : !> \param matrix_s real-space overlap matrix
1636 : !> \param matrix_ks real-space Kohn-Sham matrix
1637 : !> \param cell_to_index real-space cell index table
1638 : !> \param sab_nl overlap neighbor list
1639 : !> \param para_env global parallel environment
1640 : !> \param success true if all full-mesh MOs were reconstructed
1641 : !> \param reason diagnostic message when reconstruction is not possible
1642 : !> \param aligned_degenerate_blocks number of aligned degenerate MO blocks
1643 : !> \param aligned_degenerate_max_size largest aligned degenerate MO block
1644 : !> \param aligned_degenerate_min_svalue smallest S(k)-metric subspace singular value
1645 : !> \param spin_channel Optional single spin channel; otherwise reconstruct all channels.
1646 : ! **************************************************************************************************
1647 46 : SUBROUTINE prepare_wannier90_scf_mos(kpoint, qs_kpoint, matrix_s, matrix_ks, cell_to_index, &
1648 : sab_nl, para_env, success, reason, aligned_degenerate_blocks, &
1649 : aligned_degenerate_max_size, aligned_degenerate_min_svalue, spin_channel)
1650 : TYPE(kpoint_type), POINTER :: kpoint, qs_kpoint
1651 : TYPE(dbcsr_p_type), POINTER :: matrix_s(:, :), matrix_ks(:, :)
1652 : INTEGER, POINTER :: cell_to_index(:, :, :)
1653 : TYPE(neighbor_list_set_p_type), POINTER :: sab_nl(:)
1654 : TYPE(mp_para_env_type), POINTER :: para_env
1655 : LOGICAL, INTENT(OUT) :: success
1656 : CHARACTER(LEN=*), INTENT(OUT) :: reason
1657 : INTEGER, INTENT(OUT) :: aligned_degenerate_blocks, &
1658 : aligned_degenerate_max_size
1659 : REAL(KIND=dp), INTENT(OUT) :: aligned_degenerate_min_svalue
1660 : INTEGER, INTENT(IN), OPTIONAL :: spin_channel
1661 :
1662 : INTEGER :: blocks, first_spin, ispin, last_spin, &
1663 : max_size
1664 : REAL(KIND=dp) :: min_svalue
1665 :
1666 46 : success = .FALSE.
1667 46 : reason = 'Wannier90 export k-point MOs are not initialized'
1668 46 : aligned_degenerate_blocks = 0
1669 46 : aligned_degenerate_max_size = 0
1670 46 : aligned_degenerate_min_svalue = 0.0_dp
1671 52 : IF (.NOT. ASSOCIATED(kpoint)) RETURN
1672 46 : IF (.NOT. ASSOCIATED(kpoint%kp_env)) RETURN
1673 46 : IF (SIZE(kpoint%kp_env) == 0) RETURN
1674 46 : first_spin = 1
1675 46 : last_spin = SIZE(kpoint%kp_env(1)%kpoint_env%mos)
1676 46 : IF (PRESENT(spin_channel)) THEN
1677 42 : CPASSERT(spin_channel >= 1 .AND. spin_channel <= last_spin)
1678 : first_spin = spin_channel
1679 : last_spin = spin_channel
1680 : END IF
1681 86 : DO ispin = first_spin, last_spin
1682 : CALL prepare_wannier90_spin_mos(kpoint, qs_kpoint, matrix_s, matrix_ks, cell_to_index, &
1683 46 : sab_nl, para_env, ispin, success, reason, blocks, max_size, min_svalue)
1684 46 : IF (.NOT. success) RETURN
1685 126 : IF (blocks > 0) THEN
1686 24 : IF (aligned_degenerate_blocks == 0) aligned_degenerate_min_svalue = min_svalue
1687 24 : aligned_degenerate_min_svalue = MIN(aligned_degenerate_min_svalue, min_svalue)
1688 24 : aligned_degenerate_blocks = aligned_degenerate_blocks + blocks
1689 24 : aligned_degenerate_max_size = MAX(aligned_degenerate_max_size, max_size)
1690 : END IF
1691 : END DO
1692 46 : END SUBROUTINE prepare_wannier90_scf_mos
1693 :
1694 : ! **************************************************************************************************
1695 : !> \brief Reconstruct one spin channel with its own MO dimensions and Hamiltonian.
1696 : !> \param kpoint Full Wannier90 export k-point object.
1697 : !> \param qs_kpoint Source SCF k-point object.
1698 : !> \param matrix_s Real-space overlap matrices.
1699 : !> \param matrix_ks Spin-dependent real-space Hamiltonians.
1700 : !> \param cell_to_index Real-space cell index table.
1701 : !> \param sab_nl Overlap neighbour list.
1702 : !> \param para_env Global parallel environment.
1703 : !> \param ispin Selected spin channel.
1704 : !> \param success Whether all MOs in this channel were reconstructed.
1705 : !> \param reason Explanation if reconstruction is unavailable.
1706 : !> \param aligned_degenerate_blocks Number of aligned blocks.
1707 : !> \param aligned_degenerate_max_size Largest aligned block.
1708 : !> \param aligned_degenerate_min_svalue Smallest metric singular value.
1709 : ! **************************************************************************************************
1710 46 : SUBROUTINE prepare_wannier90_spin_mos(kpoint, qs_kpoint, matrix_s, matrix_ks, cell_to_index, &
1711 : sab_nl, para_env, ispin, success, reason, aligned_degenerate_blocks, &
1712 : aligned_degenerate_max_size, &
1713 : aligned_degenerate_min_svalue)
1714 : TYPE(kpoint_type), POINTER :: kpoint, qs_kpoint
1715 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_ks
1716 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1717 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1718 : POINTER :: sab_nl
1719 : TYPE(mp_para_env_type), POINTER :: para_env
1720 : INTEGER, INTENT(IN) :: ispin
1721 : LOGICAL, INTENT(OUT) :: success
1722 : CHARACTER(LEN=*), INTENT(OUT) :: reason
1723 : INTEGER, INTENT(OUT) :: aligned_degenerate_blocks, &
1724 : aligned_degenerate_max_size
1725 : REAL(KIND=dp), INTENT(OUT) :: aligned_degenerate_min_svalue
1726 :
1727 : CHARACTER(LEN=default_string_length) :: best_reason, candidate_reason
1728 : INTEGER :: aligned_blocks, aligned_max_size, candidate_aligned_blocks, &
1729 : candidate_aligned_max_size, ik, ikpgr, ikred, isym_try, min_gap_band, min_gap_kpoint, &
1730 : min_gap_spin, nao, nao_src, nmo, nmo_src, nsymmetry, num_candidates
1731 46 : INTEGER, ALLOCATABLE, DIMENSION(:) :: source_kpoint, sym_index
1732 : INTEGER, DIMENSION(2) :: kp_range, source_kp_range
1733 : LOGICAL :: my_kpgrp, my_source_kpgrp, ok, &
1734 : source_window
1735 : REAL(KIND=dp) :: aligned_min_svalue, band_gap, best_residual, candidate_residual, &
1736 : degenerate_band_tol, local_min_band_gap, min_band_gap, source_owner_count, &
1737 : source_window_min_svalue
1738 46 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_buffer, occupation_buffer, &
1739 46 : source_eigenvalues_buffer, &
1740 46 : source_occupation_buffer
1741 46 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation
1742 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1743 : TYPE(cp_cfm_type) :: dst_cfm_full, src_cfm, src_cfm_full
1744 : TYPE(cp_cfm_type), POINTER :: destination_cfm, source_cfm
1745 : TYPE(cp_cfm_type), TARGET :: dst_cfm, empty_destination
1746 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_source, matrix_struct_work
1747 : TYPE(kpoint_env_type), POINTER :: kp, kp_source
1748 : TYPE(kpoint_sym_type), POINTER :: kpsym
1749 :
1750 46 : success = .FALSE.
1751 46 : reason = ""
1752 46 : aligned_degenerate_blocks = 0
1753 46 : aligned_degenerate_max_size = 0
1754 46 : aligned_degenerate_min_svalue = HUGE(1.0_dp)
1755 46 : NULLIFY (matrix_struct_source, matrix_struct_work)
1756 :
1757 46 : IF (.NOT. ASSOCIATED(kpoint)) THEN
1758 0 : reason = "internal Wannier90 k-point object is not available"
1759 0 : RETURN
1760 : END IF
1761 46 : IF (.NOT. ASSOCIATED(qs_kpoint)) THEN
1762 0 : reason = "SCF k-point object is not available"
1763 0 : RETURN
1764 : END IF
1765 46 : IF (.NOT. ASSOCIATED(kpoint%kp_env) .OR. .NOT. ASSOCIATED(qs_kpoint%kp_env)) THEN
1766 0 : reason = "k-point MO environments are not initialized"
1767 0 : RETURN
1768 : END IF
1769 46 : IF (.NOT. ASSOCIATED(kpoint%blacs_env)) THEN
1770 0 : reason = "Wannier90 k-point BLACS environment is not initialized"
1771 0 : RETURN
1772 : END IF
1773 :
1774 46 : CALL build_wannier90_scf_mapping(kpoint, qs_kpoint, source_kpoint, sym_index, ok, reason)
1775 46 : IF (.NOT. ok) RETURN
1776 626 : nsymmetry = COUNT(sym_index > 0)
1777 :
1778 46 : kp => kpoint%kp_env(1)%kpoint_env
1779 46 : IF (.NOT. ASSOCIATED(kp%mos(ispin)%cmo_coeff)) THEN
1780 0 : reason = "Wannier90 export k-point MOs are not complex-valued"
1781 0 : DEALLOCATE (source_kpoint, sym_index)
1782 0 : RETURN
1783 : END IF
1784 46 : CALL get_mo_set(kp%mos(ispin), nao=nao, nmo=nmo)
1785 :
1786 46 : kp_source => qs_kpoint%kp_env(1)%kpoint_env
1787 46 : IF (.NOT. ASSOCIATED(kp_source%mos(ispin)%cmo_coeff)) THEN
1788 0 : reason = "SCF MOs are real-valued; complex symmetry phases cannot be reconstructed"
1789 0 : DEALLOCATE (source_kpoint, sym_index)
1790 0 : RETURN
1791 : END IF
1792 46 : CALL get_mo_set(kp_source%mos(ispin), nao=nao_src, nmo=nmo_src)
1793 46 : CALL para_env%max(nao_src)
1794 46 : CALL para_env%max(nmo_src)
1795 46 : IF (nao_src /= nao) THEN
1796 0 : reason = "SCF and Wannier90 MO bases have different AO dimensions"
1797 0 : DEALLOCATE (source_kpoint, sym_index)
1798 0 : RETURN
1799 : END IF
1800 46 : IF (nmo_src < nmo) THEN
1801 0 : reason = "SCF MO set has fewer bands than the Wannier90 export"
1802 0 : DEALLOCATE (source_kpoint, sym_index)
1803 0 : RETURN
1804 : END IF
1805 46 : source_window = nmo_src > nmo
1806 46 : degenerate_band_tol = 1.0e-8_dp
1807 46 : CALL get_kpoint_info(qs_kpoint, kp_range=source_kp_range)
1808 46 : IF (source_kp_range(1) /= 1 .OR. source_kp_range(2) /= qs_kpoint%nkp) THEN
1809 6 : reason = "SCF k-point symmetry data are distributed over k-point parallel groups"
1810 6 : DEALLOCATE (source_kpoint, sym_index)
1811 6 : RETURN
1812 : END IF
1813 : ! Positive symmetry entries require atom/AO rotations and Bloch phases. Degenerate subspaces
1814 : ! fully contained in the exported band window are aligned below; only guard when the Wannier90
1815 : ! window cuts through a degenerate SCF manifold at the upper band edge.
1816 40 : IF (nsymmetry > 0 .AND. nmo_src > nmo) THEN
1817 0 : local_min_band_gap = HUGE(1.0_dp)
1818 0 : min_gap_band = nmo
1819 0 : min_gap_kpoint = 0
1820 0 : min_gap_spin = 0
1821 0 : DO ikred = source_kp_range(1), source_kp_range(2)
1822 0 : ikpgr = ikred - source_kp_range(1) + 1
1823 0 : kp_source => qs_kpoint%kp_env(ikpgr)%kpoint_env
1824 0 : CALL get_mo_set(kp_source%mos(ispin), eigenvalues=eigenvalues)
1825 0 : band_gap = ABS(eigenvalues(nmo + 1) - eigenvalues(nmo))
1826 0 : IF (band_gap < local_min_band_gap) THEN
1827 0 : local_min_band_gap = band_gap
1828 0 : min_gap_band = nmo
1829 0 : min_gap_kpoint = ikred
1830 0 : min_gap_spin = ispin
1831 : END IF
1832 : END DO
1833 0 : min_band_gap = local_min_band_gap
1834 0 : CALL para_env%min(min_band_gap)
1835 0 : IF (ABS(local_min_band_gap - min_band_gap) > degenerate_band_tol*EPSILON(1.0_dp)) THEN
1836 0 : min_gap_kpoint = 0
1837 0 : min_gap_spin = 0
1838 : END IF
1839 0 : CALL para_env%max(min_gap_kpoint)
1840 0 : CALL para_env%max(min_gap_spin)
1841 0 : CALL para_env%max(min_gap_band)
1842 0 : IF (min_band_gap < degenerate_band_tol) THEN
1843 : WRITE (reason, "(A,ES9.2,A,I0,A,I0,A,I0)") &
1844 0 : "degenerate atom/AO W90 reuse guarded: edge gap=", min_band_gap, ", k=", &
1845 0 : min_gap_kpoint, ", s=", min_gap_spin, ", nband=", min_gap_band
1846 0 : DEALLOCATE (source_kpoint, sym_index)
1847 0 : RETURN
1848 : END IF
1849 : END IF
1850 40 : blacs_env => kpoint%blacs_env
1851 : CALL cp_fm_struct_create(matrix_struct_work, nrow_global=nao, ncol_global=nmo, &
1852 40 : para_env=para_env, context=blacs_env)
1853 40 : CALL cp_cfm_create(src_cfm, matrix_struct_work)
1854 40 : CALL cp_cfm_create(dst_cfm, matrix_struct_work)
1855 40 : IF (source_window) THEN
1856 : CALL cp_fm_struct_create(matrix_struct_source, nrow_global=nao, ncol_global=nmo_src, &
1857 0 : para_env=para_env, context=blacs_env)
1858 0 : CALL cp_cfm_create(src_cfm_full, matrix_struct_source)
1859 0 : CALL cp_cfm_create(dst_cfm_full, matrix_struct_source)
1860 : END IF
1861 0 : ALLOCATE (eigenvalues_buffer(nmo), occupation_buffer(nmo), &
1862 280 : source_eigenvalues_buffer(nmo_src), source_occupation_buffer(nmo_src))
1863 :
1864 40 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
1865 40 : CALL get_kpoint_info(qs_kpoint, kp_range=source_kp_range)
1866 :
1867 572 : DO ik = 1, kpoint%nkp
1868 532 : ikred = source_kpoint(ik)
1869 532 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
1870 532 : my_source_kpgrp = (ikred >= source_kp_range(1) .AND. ikred <= source_kp_range(2))
1871 2628 : source_eigenvalues_buffer(1:nmo_src) = 0.0_dp
1872 2628 : source_occupation_buffer(1:nmo_src) = 0.0_dp
1873 532 : NULLIFY (source_cfm)
1874 532 : IF (my_source_kpgrp) THEN
1875 532 : ikpgr = ikred - source_kp_range(1) + 1
1876 532 : kp_source => qs_kpoint%kp_env(ikpgr)%kpoint_env
1877 532 : source_cfm => kp_source%mos(ispin)%cmo_coeff
1878 : CALL get_mo_set(kp_source%mos(ispin), eigenvalues=eigenvalues, &
1879 532 : occupation_numbers=occupation)
1880 2628 : source_eigenvalues_buffer(1:nmo_src) = eigenvalues(1:nmo_src)
1881 2628 : source_occupation_buffer(1:nmo_src) = occupation(1:nmo_src)
1882 : END IF
1883 532 : IF (source_window) THEN
1884 0 : CALL wannier90_copy_cfm_general(source_cfm, src_cfm_full, para_env)
1885 : ELSE
1886 532 : CALL wannier90_copy_cfm_general(source_cfm, src_cfm, para_env)
1887 : END IF
1888 532 : IF (my_source_kpgrp) THEN
1889 532 : source_owner_count = 1.0_dp
1890 : ELSE
1891 0 : source_owner_count = 0.0_dp
1892 : END IF
1893 532 : CALL para_env%sum(source_owner_count)
1894 532 : CALL para_env%sum(source_eigenvalues_buffer)
1895 532 : CALL para_env%sum(source_occupation_buffer)
1896 532 : IF (source_owner_count > 0.0_dp) THEN
1897 : source_eigenvalues_buffer(1:nmo_src) = &
1898 2628 : source_eigenvalues_buffer(1:nmo_src)/source_owner_count
1899 : source_occupation_buffer(1:nmo_src) = &
1900 2628 : source_occupation_buffer(1:nmo_src)/source_owner_count
1901 : END IF
1902 2628 : eigenvalues_buffer(1:nmo) = source_eigenvalues_buffer(1:nmo)
1903 2628 : occupation_buffer(1:nmo) = source_occupation_buffer(1:nmo)
1904 532 : IF (source_window) THEN
1905 0 : CALL cp_cfm_to_cfm(src_cfm_full, src_cfm, nmo)
1906 : END IF
1907 :
1908 532 : ok = .FALSE.
1909 532 : reason = ""
1910 532 : aligned_blocks = 0
1911 532 : aligned_max_size = 0
1912 532 : aligned_min_svalue = 0.0_dp
1913 532 : IF (sym_index(ik) > 0) THEN
1914 412 : kpsym => qs_kpoint%kp_sym(ikred)%kpoint_sym
1915 412 : IF (ASSOCIATED(kpsym)) THEN
1916 412 : best_reason = ""
1917 412 : best_residual = HUGE(1.0_dp)
1918 412 : num_candidates = 0
1919 : ! Little-group operations can reach the same target k-point; keep the first valid eigenspace.
1920 3198 : DO isym_try = 1, kpsym%nwred
1921 3198 : IF (.NOT. kpoint_same_periodic(kpoint%xkp(1:3, ik), &
1922 : kpsym%xkp(1:3, isym_try))) CYCLE
1923 412 : num_candidates = num_candidates + 1
1924 : CALL kpoint_transform_scf_mo(src_cfm, dst_cfm, &
1925 : qs_kpoint, ikred, isym_try, ok, &
1926 412 : candidate_reason)
1927 412 : IF (.NOT. ok) THEN
1928 0 : reason = candidate_reason
1929 : CYCLE
1930 : END IF
1931 : CALL ritz_stabilize_wannier90_subspace(dst_cfm, matrix_s, matrix_ks, &
1932 : kpoint%xkp(1:3, ik), cell_to_index, &
1933 : sab_nl, ispin, eigenvalues_buffer, &
1934 : degenerate_band_tol, ok, candidate_reason, &
1935 : candidate_aligned_blocks, &
1936 : candidate_aligned_max_size, &
1937 412 : aligned_min_svalue, candidate_residual)
1938 412 : IF (candidate_residual < best_residual) THEN
1939 412 : best_residual = candidate_residual
1940 412 : best_reason = candidate_reason
1941 : END IF
1942 412 : IF (.NOT. ok) THEN
1943 0 : IF (source_window) THEN
1944 : CALL kpoint_transform_scf_mo(src_cfm_full, dst_cfm_full, qs_kpoint, &
1945 0 : ikred, isym_try, ok, candidate_reason)
1946 0 : IF (ok) THEN
1947 : CALL ritz_reconstruct_wannier90_window(dst_cfm_full, dst_cfm, matrix_s, &
1948 : matrix_ks, kpoint%xkp(1:3, ik), &
1949 : cell_to_index, sab_nl, ispin, &
1950 : eigenvalues_buffer, nmo, ok, &
1951 : candidate_reason, source_window_min_svalue, &
1952 0 : candidate_residual)
1953 0 : IF (candidate_residual < best_residual) THEN
1954 0 : best_residual = candidate_residual
1955 0 : best_reason = candidate_reason
1956 : END IF
1957 : END IF
1958 : ELSE
1959 : CALL ritz_reconstruct_wannier90_window(dst_cfm, dst_cfm, matrix_s, matrix_ks, &
1960 : kpoint%xkp(1:3, ik), cell_to_index, &
1961 : sab_nl, ispin, eigenvalues_buffer, nmo, &
1962 : ok, candidate_reason, source_window_min_svalue, &
1963 0 : candidate_residual)
1964 0 : IF (candidate_residual < best_residual) THEN
1965 0 : best_residual = candidate_residual
1966 0 : best_reason = candidate_reason
1967 : END IF
1968 : END IF
1969 : END IF
1970 412 : IF (ok) THEN
1971 412 : aligned_blocks = candidate_aligned_blocks
1972 412 : aligned_max_size = candidate_aligned_max_size
1973 412 : sym_index(ik) = isym_try
1974 412 : EXIT
1975 : END IF
1976 412 : reason = candidate_reason
1977 : END DO
1978 412 : IF (.NOT. ok .AND. num_candidates > 0 .AND. best_residual < HUGE(1.0_dp)) THEN
1979 : WRITE (reason, "(A,I0,A,ES9.2,A,I0,A,A32)") &
1980 0 : "atom/AO W90 guarded: best/", num_candidates, "=", best_residual, &
1981 0 : " k=", ik, " ", TRIM(best_reason)
1982 412 : ELSE IF (.NOT. ok .AND. num_candidates == 0) THEN
1983 0 : reason = "no matching SCF symmetry operation candidate"
1984 : END IF
1985 : ELSE
1986 0 : reason = "SCF k-point symmetry operation is not available"
1987 : END IF
1988 : ELSE
1989 : CALL kpoint_transform_scf_mo(src_cfm, dst_cfm, qs_kpoint, &
1990 120 : ikred, sym_index(ik), ok, reason)
1991 : END IF
1992 532 : IF (ok .AND. sym_index(ik) <= 0) THEN
1993 : ! Even a direct k-point copy must be a closed H(k),S(k) subspace. This catches
1994 : ! incomplete degenerate band windows before they can be exported to Wannier90.
1995 : CALL ritz_stabilize_wannier90_subspace(dst_cfm, matrix_s, matrix_ks, &
1996 : kpoint%xkp(1:3, ik), cell_to_index, sab_nl, &
1997 : ispin, eigenvalues_buffer, degenerate_band_tol, &
1998 : ok, reason, aligned_blocks, aligned_max_size, &
1999 120 : aligned_min_svalue, candidate_residual)
2000 120 : IF (.NOT. ok .AND. source_window) THEN
2001 : CALL kpoint_transform_scf_mo(src_cfm_full, dst_cfm_full, qs_kpoint, &
2002 0 : ikred, sym_index(ik), ok, reason)
2003 0 : IF (ok) THEN
2004 : CALL ritz_reconstruct_wannier90_window(dst_cfm_full, dst_cfm, matrix_s, matrix_ks, &
2005 : kpoint%xkp(1:3, ik), cell_to_index, &
2006 : sab_nl, ispin, eigenvalues_buffer, nmo, ok, &
2007 : reason, source_window_min_svalue, &
2008 0 : candidate_residual)
2009 : END IF
2010 120 : ELSE IF (.NOT. ok) THEN
2011 : CALL ritz_reconstruct_wannier90_window(dst_cfm, dst_cfm, matrix_s, matrix_ks, &
2012 : kpoint%xkp(1:3, ik), &
2013 : cell_to_index, sab_nl, ispin, &
2014 : eigenvalues_buffer, nmo, ok, reason, &
2015 0 : source_window_min_svalue, candidate_residual)
2016 : END IF
2017 : END IF
2018 532 : IF (.NOT. ok) THEN
2019 0 : CALL cp_cfm_release(src_cfm)
2020 0 : CALL cp_cfm_release(dst_cfm)
2021 0 : CALL cp_fm_struct_release(matrix_struct_work)
2022 0 : IF (source_window) THEN
2023 0 : CALL cp_cfm_release(src_cfm_full)
2024 0 : CALL cp_cfm_release(dst_cfm_full)
2025 0 : CALL cp_fm_struct_release(matrix_struct_source)
2026 : END IF
2027 0 : DEALLOCATE (source_kpoint, sym_index, eigenvalues_buffer, occupation_buffer, &
2028 0 : source_eigenvalues_buffer, source_occupation_buffer)
2029 0 : RETURN
2030 : END IF
2031 532 : IF (sym_index(ik) /= 0) THEN
2032 456 : aligned_degenerate_blocks = aligned_degenerate_blocks + aligned_blocks
2033 456 : aligned_degenerate_max_size = MAX(aligned_degenerate_max_size, aligned_max_size)
2034 456 : IF (aligned_blocks > 0) THEN
2035 360 : aligned_degenerate_min_svalue = MIN(aligned_degenerate_min_svalue, aligned_min_svalue)
2036 : END IF
2037 : END IF
2038 :
2039 532 : IF (my_kpgrp) THEN
2040 532 : ikpgr = ik - kp_range(1) + 1
2041 532 : kp => kpoint%kp_env(ikpgr)%kpoint_env
2042 : CALL get_mo_set(kp%mos(ispin), eigenvalues=eigenvalues, &
2043 532 : occupation_numbers=occupation)
2044 2628 : eigenvalues(1:nmo) = eigenvalues_buffer(1:nmo)
2045 2628 : occupation(1:nmo) = occupation_buffer(1:nmo)
2046 532 : destination_cfm => kp%mos(ispin)%cmo_coeff
2047 : ELSE
2048 : destination_cfm => empty_destination
2049 : END IF
2050 572 : CALL cp_cfm_copy_general(dst_cfm, destination_cfm, para_env)
2051 : END DO
2052 :
2053 40 : CALL cp_cfm_release(src_cfm)
2054 40 : CALL cp_cfm_release(dst_cfm)
2055 40 : CALL cp_fm_struct_release(matrix_struct_work)
2056 40 : IF (source_window) THEN
2057 0 : CALL cp_cfm_release(src_cfm_full)
2058 0 : CALL cp_cfm_release(dst_cfm_full)
2059 0 : CALL cp_fm_struct_release(matrix_struct_source)
2060 : END IF
2061 0 : DEALLOCATE (source_kpoint, sym_index, eigenvalues_buffer, occupation_buffer, &
2062 40 : source_eigenvalues_buffer, source_occupation_buffer)
2063 40 : IF (aligned_degenerate_blocks == 0) aligned_degenerate_min_svalue = 0.0_dp
2064 40 : success = .TRUE.
2065 :
2066 224 : END SUBROUTINE prepare_wannier90_spin_mos
2067 :
2068 : ! **************************************************************************************************
2069 : !> \brief Save a full-mesh Wannier90 MO reference on all ranks for diagnostic validation.
2070 : !> \param kpoint full Wannier90 export k-point object
2071 : !> \param ispin Selected spin channel.
2072 : !> \param para_env global parallel environment
2073 : !> \param mo MO coefficients for the selected spin, indexed as AO, MO, k-point
2074 : !> \param eigenvalue_snapshot MO eigenvalues for the selected spin, indexed as MO, k-point
2075 : ! **************************************************************************************************
2076 14 : SUBROUTINE save_wannier90_mo_snapshot(kpoint, ispin, para_env, mo, &
2077 : eigenvalue_snapshot)
2078 : TYPE(kpoint_type), POINTER :: kpoint
2079 : INTEGER, INTENT(IN) :: ispin
2080 : TYPE(mp_para_env_type), POINTER :: para_env
2081 : COMPLEX(KIND=dp), ALLOCATABLE, &
2082 : DIMENSION(:, :, :), INTENT(OUT) :: mo
2083 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
2084 : INTENT(OUT) :: eigenvalue_snapshot
2085 :
2086 : INTEGER :: ik, ikpgr, nao, nkp, nmo
2087 : INTEGER, DIMENSION(2) :: kp_range
2088 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: owner_weight
2089 14 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
2090 : TYPE(kpoint_env_type), POINTER :: kp
2091 :
2092 14 : CALL get_kpoint_info(kpoint, nkp=nkp, kp_range=kp_range)
2093 14 : kp => kpoint%kp_env(1)%kpoint_env
2094 14 : CALL get_mo_set(kp%mos(ispin), nao=nao, nmo=nmo)
2095 0 : ALLOCATE (mo(nao, nmo, nkp), &
2096 140 : eigenvalue_snapshot(nmo, nkp), owner_weight(nkp))
2097 14 : mo(:, :, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2098 14 : eigenvalue_snapshot(:, :) = 0.0_dp
2099 14 : owner_weight(:) = 0.0_dp
2100 350 : DO ik = kp_range(1), kp_range(2)
2101 336 : ikpgr = ik - kp_range(1) + 1
2102 336 : kp => kpoint%kp_env(ikpgr)%kpoint_env
2103 336 : CALL cp_cfm_get_submatrix(kp%mos(ispin)%cmo_coeff, mo(:, :, ik))
2104 336 : CALL get_mo_set(kp%mos(ispin), eigenvalues=eigenvalues)
2105 1600 : eigenvalue_snapshot(1:nmo, ik) = eigenvalues(1:nmo)
2106 350 : owner_weight(ik) = 1.0_dp
2107 : END DO
2108 14 : CALL para_env%sum(mo)
2109 14 : CALL para_env%sum(eigenvalue_snapshot)
2110 14 : CALL para_env%sum(owner_weight)
2111 350 : DO ik = 1, nkp
2112 350 : IF (owner_weight(ik) > 0.0_dp) THEN
2113 50000 : mo(:, :, ik) = mo(:, :, ik)/owner_weight(ik)
2114 : eigenvalue_snapshot(:, ik) = &
2115 1600 : eigenvalue_snapshot(:, ik)/owner_weight(ik)
2116 : END IF
2117 : END DO
2118 14 : DEALLOCATE (owner_weight)
2119 :
2120 14 : END SUBROUTINE save_wannier90_mo_snapshot
2121 :
2122 : ! **************************************************************************************************
2123 : !> \brief Restore a full-mesh Wannier90 MO reference after a failed diagnostic reuse attempt.
2124 : !> \param kpoint full Wannier90 export k-point object
2125 : !> \param ispin Selected spin channel.
2126 : !> \param mo complex MO coefficient snapshot
2127 : !> \param eigenvalue_snapshot MO eigenvalue snapshot
2128 : ! **************************************************************************************************
2129 2 : SUBROUTINE restore_wannier90_mo_snapshot(kpoint, ispin, mo, eigenvalue_snapshot)
2130 : TYPE(kpoint_type), POINTER :: kpoint
2131 : INTEGER, INTENT(IN) :: ispin
2132 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: mo
2133 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eigenvalue_snapshot
2134 :
2135 : INTEGER :: ik, ikpgr, nmo
2136 : INTEGER, DIMENSION(2) :: kp_range
2137 2 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
2138 : TYPE(kpoint_env_type), POINTER :: kp
2139 :
2140 2 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
2141 2 : nmo = SIZE(eigenvalue_snapshot, 1)
2142 18 : DO ik = kp_range(1), kp_range(2)
2143 16 : ikpgr = ik - kp_range(1) + 1
2144 16 : kp => kpoint%kp_env(ikpgr)%kpoint_env
2145 16 : CALL cp_cfm_set_submatrix(kp%mos(ispin)%cmo_coeff, mo(:, :, ik))
2146 16 : CALL get_mo_set(kp%mos(ispin), eigenvalues=eigenvalues)
2147 146 : eigenvalues(1:nmo) = eigenvalue_snapshot(1:nmo, ik)
2148 : END DO
2149 :
2150 2 : END SUBROUTINE restore_wannier90_mo_snapshot
2151 :
2152 : ! **************************************************************************************************
2153 : !> \brief Validate current Wannier90 MOs against a saved full-mesh diagonalization reference.
2154 : !> \param kpoint full Wannier90 export k-point object
2155 : !> \param matrix_s real-space overlap matrix
2156 : !> \param cell_to_index real-space cell index table
2157 : !> \param sab_nl overlap neighbor list
2158 : !> \param para_env global parallel environment
2159 : !> \param ispin Selected spin channel.
2160 : !> \param reference_mo complex MO coefficient reference
2161 : !> \param reference_eigenvalues MO eigenvalue reference
2162 : !> \param success true if the reconstructed MOs match the reference subspaces
2163 : !> \param max_subspace_deviation largest deviation of S(k)-metric singular values from one
2164 : !> \param min_svalue smallest S(k)-metric singular value
2165 : !> \param max_eigenvalue_deviation largest eigenvalue deviation
2166 : ! **************************************************************************************************
2167 14 : SUBROUTINE validate_wannier90_reused_mos(kpoint, matrix_s, cell_to_index, sab_nl, para_env, ispin, &
2168 14 : reference_mo, reference_eigenvalues, &
2169 : success, max_subspace_deviation, min_svalue, &
2170 : max_eigenvalue_deviation)
2171 : TYPE(kpoint_type), POINTER :: kpoint
2172 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
2173 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2174 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2175 : POINTER :: sab_nl
2176 : TYPE(mp_para_env_type), POINTER :: para_env
2177 : INTEGER, INTENT(IN) :: ispin
2178 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: reference_mo
2179 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: reference_eigenvalues
2180 : LOGICAL, INTENT(OUT) :: success
2181 : REAL(KIND=dp), INTENT(OUT) :: max_subspace_deviation, min_svalue, &
2182 : max_eigenvalue_deviation
2183 :
2184 : REAL(KIND=dp), PARAMETER :: eigenvalue_tol = 1.0e-8_dp, &
2185 : subspace_tol = 1.0e-4_dp
2186 :
2187 : INTEGER :: ik, ikpgr, nao, nkp, nmo
2188 : INTEGER, DIMENSION(2) :: kp_range
2189 : LOGICAL :: my_kpgrp, ok
2190 : REAL(KIND=dp) :: candidate_deviation, candidate_svalue, &
2191 : owner_count
2192 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalue_buffer
2193 14 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
2194 : TYPE(cp_cfm_type) :: cand_cfm, ref_cfm
2195 : TYPE(cp_cfm_type), POINTER :: candidate_source
2196 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_work
2197 : TYPE(kpoint_env_type), POINTER :: kp
2198 :
2199 14 : success = .FALSE.
2200 14 : max_subspace_deviation = 0.0_dp
2201 14 : min_svalue = HUGE(1.0_dp)
2202 14 : max_eigenvalue_deviation = 0.0_dp
2203 14 : NULLIFY (matrix_struct_work)
2204 :
2205 14 : CALL get_kpoint_info(kpoint, nkp=nkp, kp_range=kp_range)
2206 14 : kp => kpoint%kp_env(1)%kpoint_env
2207 14 : CALL get_mo_set(kp%mos(ispin), nao=nao, nmo=nmo)
2208 14 : CPASSERT(SIZE(reference_mo, 1) == nao)
2209 14 : CPASSERT(SIZE(reference_mo, 2) == nmo)
2210 14 : CPASSERT(SIZE(reference_mo, 3) == nkp)
2211 :
2212 14 : matrix_struct_work => kp%mos(ispin)%cmo_coeff%matrix_struct
2213 14 : CALL cp_cfm_create(ref_cfm, matrix_struct_work)
2214 14 : CALL cp_cfm_create(cand_cfm, matrix_struct_work)
2215 42 : ALLOCATE (eigenvalue_buffer(nmo))
2216 :
2217 350 : DO ik = 1, nkp
2218 336 : CALL cp_cfm_set_submatrix(ref_cfm, reference_mo(:, :, ik))
2219 336 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
2220 0 : NULLIFY (candidate_source)
2221 : IF (my_kpgrp) THEN
2222 336 : ikpgr = ik - kp_range(1) + 1
2223 336 : kp => kpoint%kp_env(ikpgr)%kpoint_env
2224 336 : candidate_source => kp%mos(ispin)%cmo_coeff
2225 336 : CALL get_mo_set(kp%mos(ispin), eigenvalues=eigenvalues)
2226 1600 : eigenvalue_buffer(1:nmo) = eigenvalues(1:nmo)
2227 : ELSE
2228 0 : eigenvalue_buffer(1:nmo) = 0.0_dp
2229 : END IF
2230 336 : CALL wannier90_copy_cfm_general(candidate_source, cand_cfm, para_env)
2231 336 : IF (my_kpgrp) THEN
2232 336 : owner_count = 1.0_dp
2233 : ELSE
2234 0 : owner_count = 0.0_dp
2235 : END IF
2236 336 : CALL para_env%sum(owner_count)
2237 336 : CALL para_env%sum(eigenvalue_buffer)
2238 1600 : IF (owner_count > 0.0_dp) eigenvalue_buffer(1:nmo) = eigenvalue_buffer(1:nmo)/owner_count
2239 : max_eigenvalue_deviation = MAX(max_eigenvalue_deviation, &
2240 : MAXVAL(ABS(eigenvalue_buffer(1:nmo) - &
2241 1600 : reference_eigenvalues(1:nmo, ik))))
2242 : CALL measure_wannier90_subspace_error(ref_cfm, cand_cfm, matrix_s, &
2243 : kpoint%xkp(1:3, ik), cell_to_index, sab_nl, para_env, &
2244 336 : ok, candidate_deviation, candidate_svalue)
2245 350 : IF (.NOT. ok) THEN
2246 0 : max_subspace_deviation = HUGE(1.0_dp)
2247 : ELSE
2248 336 : max_subspace_deviation = MAX(max_subspace_deviation, candidate_deviation)
2249 336 : min_svalue = MIN(min_svalue, candidate_svalue)
2250 : END IF
2251 : END DO
2252 14 : CALL para_env%max(max_subspace_deviation)
2253 14 : CALL para_env%min(min_svalue)
2254 14 : CALL para_env%max(max_eigenvalue_deviation)
2255 14 : success = max_subspace_deviation < subspace_tol .AND. max_eigenvalue_deviation < eigenvalue_tol
2256 :
2257 14 : DEALLOCATE (eigenvalue_buffer)
2258 14 : CALL cp_cfm_release(ref_cfm)
2259 14 : CALL cp_cfm_release(cand_cfm)
2260 :
2261 28 : END SUBROUTINE validate_wannier90_reused_mos
2262 :
2263 : ! **************************************************************************************************
2264 : !> \brief Compare atom/AO reuse candidates directly to the full-mesh reference MOs.
2265 : !> \param kpoint full Wannier90 export k-point object holding reference MOs
2266 : !> \param qs_kpoint SCF k-point object
2267 : !> \param matrix_s real-space overlap matrix
2268 : !> \param matrix_ks real-space Kohn-Sham matrix
2269 : !> \param cell_to_index real-space cell index table
2270 : !> \param sab_nl overlap neighbor list
2271 : !> \param para_env global parallel environment
2272 : !> \param iw output unit
2273 : !> \param ispin Selected spin channel.
2274 : !> \param max_subspace_deviation largest best-candidate subspace deviation
2275 : !> \param min_svalue smallest best-candidate singular value
2276 : !> \param max_metric_deviation largest S(k)-metric deviation of a candidate
2277 : !> \param max_residual largest H(k),S(k) eigen-residual of a candidate
2278 : ! **************************************************************************************************
2279 14 : SUBROUTINE diagnose_wannier90_scf_reuse_candidates(kpoint, qs_kpoint, matrix_s, matrix_ks, &
2280 : cell_to_index, sab_nl, para_env, iw, ispin, &
2281 : max_subspace_deviation, min_svalue, &
2282 : max_metric_deviation, max_residual)
2283 : TYPE(kpoint_type), POINTER :: kpoint, qs_kpoint
2284 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_ks
2285 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2286 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2287 : POINTER :: sab_nl
2288 : TYPE(mp_para_env_type), POINTER :: para_env
2289 : INTEGER, INTENT(IN) :: iw, ispin
2290 : REAL(KIND=dp), INTENT(OUT) :: max_subspace_deviation, min_svalue, &
2291 : max_metric_deviation, max_residual
2292 :
2293 : REAL(KIND=dp), PARAMETER :: print_tol = 1.0e-4_dp, &
2294 : residual_print_tol = 1.0e-3_dp
2295 :
2296 : CHARACTER(LEN=default_string_length) :: reason
2297 : INTEGER :: ik, ikpgr, ikred, isym_try, nao, &
2298 : nao_src, nkp, nmo, nmo_src
2299 14 : INTEGER, ALLOCATABLE, DIMENSION(:) :: source_kpoint, sym_index
2300 : INTEGER, DIMENSION(2) :: kp_range, source_kp_range
2301 : LOGICAL :: my_kpgrp, ok, source_window
2302 : REAL(KIND=dp) :: best_deviation, best_metric_deviation, best_residual, best_svalue, &
2303 : candidate_deviation, candidate_metric_deviation, candidate_metric_min, &
2304 : candidate_residual, candidate_svalue, owner_count, ref_metric_deviation, ref_metric_min, &
2305 : ref_residual
2306 14 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_buffer
2307 14 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
2308 : TYPE(cp_cfm_type) :: dst_cfm, ref_cfm, src_cfm, src_cfm_full
2309 : TYPE(cp_cfm_type), POINTER :: source_cfm
2310 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_source, matrix_struct_work
2311 : TYPE(kpoint_env_type), POINTER :: kp, kp_source
2312 : TYPE(kpoint_sym_type), POINTER :: kpsym
2313 :
2314 14 : max_subspace_deviation = HUGE(1.0_dp)
2315 14 : min_svalue = 0.0_dp
2316 14 : max_metric_deviation = HUGE(1.0_dp)
2317 14 : max_residual = HUGE(1.0_dp)
2318 14 : NULLIFY (matrix_struct_source, matrix_struct_work)
2319 :
2320 14 : CALL build_wannier90_scf_mapping(kpoint, qs_kpoint, source_kpoint, sym_index, ok, reason)
2321 14 : IF (.NOT. ok) RETURN
2322 14 : kp => kpoint%kp_env(1)%kpoint_env
2323 14 : kp_source => qs_kpoint%kp_env(1)%kpoint_env
2324 14 : IF (.NOT. ASSOCIATED(kp%mos(ispin)%cmo_coeff) .OR. &
2325 : .NOT. ASSOCIATED(kp_source%mos(ispin)%cmo_coeff)) THEN
2326 0 : DEALLOCATE (source_kpoint, sym_index)
2327 0 : RETURN
2328 : END IF
2329 14 : CALL get_mo_set(kp%mos(ispin), nao=nao, nmo=nmo)
2330 14 : CALL get_mo_set(kp_source%mos(ispin), nao=nao_src, nmo=nmo_src)
2331 14 : CALL para_env%max(nao_src)
2332 14 : CALL para_env%max(nmo_src)
2333 14 : IF (nao_src /= nao .OR. nmo_src < nmo) THEN
2334 0 : DEALLOCATE (source_kpoint, sym_index)
2335 0 : RETURN
2336 : END IF
2337 14 : source_window = nmo_src > nmo
2338 14 : CALL get_kpoint_info(kpoint, nkp=nkp, kp_range=kp_range)
2339 14 : CALL get_kpoint_info(qs_kpoint, kp_range=source_kp_range)
2340 14 : IF (source_kp_range(1) /= 1 .OR. source_kp_range(2) /= qs_kpoint%nkp) THEN
2341 0 : DEALLOCATE (source_kpoint, sym_index)
2342 0 : RETURN
2343 : END IF
2344 :
2345 : CALL cp_fm_struct_create(matrix_struct_work, nrow_global=nao, ncol_global=nmo, &
2346 14 : para_env=para_env, context=kpoint%blacs_env)
2347 14 : CALL cp_cfm_create(ref_cfm, matrix_struct_work)
2348 14 : CALL cp_cfm_create(src_cfm, matrix_struct_work)
2349 14 : CALL cp_cfm_create(dst_cfm, matrix_struct_work)
2350 42 : ALLOCATE (eigenvalues_buffer(nmo))
2351 14 : IF (source_window) THEN
2352 : CALL cp_fm_struct_create(matrix_struct_source, nrow_global=nao, ncol_global=nmo_src, &
2353 0 : para_env=para_env, context=kpoint%blacs_env)
2354 0 : CALL cp_cfm_create(src_cfm_full, matrix_struct_source)
2355 : END IF
2356 :
2357 14 : max_subspace_deviation = 0.0_dp
2358 14 : min_svalue = HUGE(1.0_dp)
2359 14 : max_metric_deviation = 0.0_dp
2360 14 : max_residual = 0.0_dp
2361 350 : DO ik = 1, nkp
2362 336 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
2363 0 : NULLIFY (source_cfm)
2364 : IF (my_kpgrp) THEN
2365 336 : ikpgr = ik - kp_range(1) + 1
2366 336 : kp => kpoint%kp_env(ikpgr)%kpoint_env
2367 336 : source_cfm => kp%mos(ispin)%cmo_coeff
2368 : END IF
2369 336 : CALL wannier90_copy_cfm_general(source_cfm, ref_cfm, para_env)
2370 :
2371 336 : ikred = source_kpoint(ik)
2372 336 : my_kpgrp = (ikred >= source_kp_range(1) .AND. ikred <= source_kp_range(2))
2373 0 : NULLIFY (source_cfm)
2374 : IF (my_kpgrp) THEN
2375 336 : ikpgr = ikred - source_kp_range(1) + 1
2376 336 : kp_source => qs_kpoint%kp_env(ikpgr)%kpoint_env
2377 336 : source_cfm => kp_source%mos(ispin)%cmo_coeff
2378 336 : CALL get_mo_set(kp_source%mos(ispin), eigenvalues=eigenvalues)
2379 1600 : eigenvalues_buffer(1:nmo) = eigenvalues(1:nmo)
2380 : ELSE
2381 0 : eigenvalues_buffer(1:nmo) = 0.0_dp
2382 : END IF
2383 336 : IF (source_window) THEN
2384 0 : CALL wannier90_copy_cfm_general(source_cfm, src_cfm_full, para_env)
2385 : ELSE
2386 336 : CALL wannier90_copy_cfm_general(source_cfm, src_cfm, para_env)
2387 : END IF
2388 336 : IF (my_kpgrp) THEN
2389 336 : owner_count = 1.0_dp
2390 : ELSE
2391 0 : owner_count = 0.0_dp
2392 : END IF
2393 336 : CALL para_env%sum(owner_count)
2394 336 : CALL para_env%sum(eigenvalues_buffer)
2395 1600 : IF (owner_count > 0.0_dp) eigenvalues_buffer(1:nmo) = eigenvalues_buffer(1:nmo)/owner_count
2396 : CALL measure_wannier90_eigenspace_quality(ref_cfm, matrix_s, matrix_ks, &
2397 : kpoint%xkp(1:3, ik), cell_to_index, sab_nl, &
2398 : para_env, ispin, eigenvalues_buffer, ok, &
2399 336 : ref_metric_deviation, ref_metric_min, ref_residual)
2400 336 : IF (ok .AND. para_env%is_source() .AND. iw > 0 .AND. &
2401 : (ref_metric_deviation > print_tol .OR. ref_residual > residual_print_tol)) THEN
2402 : WRITE (iw, '(T2,A,I0,A,ES10.3,A,ES10.3,A,ES10.3)') &
2403 0 : "WANNIER90| reference k=", ik, " dM=", ref_metric_deviation, &
2404 0 : " smin=", ref_metric_min, " resid=", ref_residual
2405 : END IF
2406 336 : IF (source_window) THEN
2407 0 : CALL cp_cfm_to_cfm(src_cfm_full, src_cfm, nmo)
2408 : END IF
2409 :
2410 336 : best_deviation = HUGE(1.0_dp)
2411 336 : best_metric_deviation = 0.0_dp
2412 336 : best_residual = 0.0_dp
2413 336 : best_svalue = 0.0_dp
2414 336 : IF (sym_index(ik) <= 0) THEN
2415 : CALL kpoint_transform_scf_mo(src_cfm, dst_cfm, qs_kpoint, &
2416 68 : ikred, sym_index(ik), ok, reason)
2417 68 : IF (ok) THEN
2418 : CALL measure_wannier90_subspace_error(ref_cfm, dst_cfm, matrix_s, &
2419 : kpoint%xkp(1:3, ik), cell_to_index, sab_nl, &
2420 68 : para_env, ok, candidate_deviation, candidate_svalue)
2421 68 : IF (ok) THEN
2422 : CALL measure_wannier90_eigenspace_quality(dst_cfm, matrix_s, matrix_ks, &
2423 : kpoint%xkp(1:3, ik), cell_to_index, sab_nl, &
2424 : para_env, ispin, eigenvalues_buffer, ok, &
2425 : candidate_metric_deviation, &
2426 68 : candidate_metric_min, candidate_residual)
2427 : END IF
2428 68 : IF (ok) THEN
2429 68 : best_deviation = candidate_deviation
2430 68 : best_metric_deviation = candidate_metric_deviation
2431 68 : best_residual = candidate_residual
2432 68 : best_svalue = candidate_svalue
2433 68 : IF (para_env%is_source() .AND. iw > 0 .AND. &
2434 : (candidate_deviation > print_tol .OR. candidate_metric_deviation > print_tol .OR. &
2435 : candidate_residual > residual_print_tol)) THEN
2436 : WRITE (iw, '(T2,A,I0,A,I0,A,I0,A,ES10.3,A,ES10.3,A,ES10.3,A,ES10.3)') &
2437 0 : "WANNIER90| reuse candidate k=", ik, " src=", ikred, " sym=", &
2438 0 : sym_index(ik), " dRef=", candidate_deviation, " dM=", &
2439 0 : candidate_metric_deviation, " smin=", candidate_metric_min, &
2440 0 : " resid=", candidate_residual
2441 : END IF
2442 : END IF
2443 : END IF
2444 268 : ELSE IF (ASSOCIATED(qs_kpoint%kp_sym)) THEN
2445 268 : kpsym => qs_kpoint%kp_sym(ikred)%kpoint_sym
2446 268 : IF (ASSOCIATED(kpsym)) THEN
2447 25036 : DO isym_try = 1, kpsym%nwred
2448 24768 : IF (.NOT. kpoint_same_periodic(kpoint%xkp(1:3, ik), &
2449 : kpsym%xkp(1:3, isym_try))) CYCLE
2450 : CALL kpoint_transform_scf_mo(src_cfm, dst_cfm, &
2451 2128 : qs_kpoint, ikred, isym_try, ok, reason)
2452 2128 : IF (.NOT. ok) CYCLE
2453 : CALL measure_wannier90_subspace_error(ref_cfm, dst_cfm, &
2454 : matrix_s, kpoint%xkp(1:3, ik), cell_to_index, &
2455 : sab_nl, para_env, ok, candidate_deviation, &
2456 2128 : candidate_svalue)
2457 2128 : IF (ok) THEN
2458 : CALL measure_wannier90_eigenspace_quality(dst_cfm, matrix_s, matrix_ks, &
2459 : kpoint%xkp(1:3, ik), cell_to_index, &
2460 : sab_nl, para_env, ispin, &
2461 : eigenvalues_buffer, ok, &
2462 : candidate_metric_deviation, &
2463 2128 : candidate_metric_min, candidate_residual)
2464 : END IF
2465 4524 : IF (ok .AND. candidate_deviation < best_deviation) THEN
2466 314 : best_deviation = candidate_deviation
2467 314 : best_metric_deviation = candidate_metric_deviation
2468 314 : best_residual = candidate_residual
2469 314 : best_svalue = candidate_svalue
2470 : END IF
2471 : END DO
2472 : END IF
2473 : END IF
2474 350 : IF (best_deviation < HUGE(1.0_dp)) THEN
2475 336 : max_subspace_deviation = MAX(max_subspace_deviation, best_deviation)
2476 336 : min_svalue = MIN(min_svalue, best_svalue)
2477 336 : max_metric_deviation = MAX(max_metric_deviation, best_metric_deviation)
2478 336 : max_residual = MAX(max_residual, best_residual)
2479 : END IF
2480 : END DO
2481 14 : CALL para_env%max(max_subspace_deviation)
2482 14 : CALL para_env%min(min_svalue)
2483 14 : CALL para_env%max(max_metric_deviation)
2484 14 : CALL para_env%max(max_residual)
2485 :
2486 14 : IF (source_window) THEN
2487 0 : CALL cp_cfm_release(src_cfm_full)
2488 0 : CALL cp_fm_struct_release(matrix_struct_source)
2489 : END IF
2490 14 : CALL cp_cfm_release(ref_cfm)
2491 14 : CALL cp_cfm_release(src_cfm)
2492 14 : CALL cp_cfm_release(dst_cfm)
2493 14 : CALL cp_fm_struct_release(matrix_struct_work)
2494 14 : DEALLOCATE (eigenvalues_buffer)
2495 14 : DEALLOCATE (source_kpoint, sym_index)
2496 :
2497 56 : END SUBROUTINE diagnose_wannier90_scf_reuse_candidates
2498 :
2499 : ! **************************************************************************************************
2500 : !> \brief Measure the S(k)-metric distance between two Wannier90 MO subspaces.
2501 : !> \param ref_cfm reference MO coefficients
2502 : !> \param cand_cfm candidate MO coefficients
2503 : !> \param matrix_s real-space overlap matrix
2504 : !> \param xkp target k-point coordinate
2505 : !> \param cell_to_index real-space cell index table
2506 : !> \param sab_nl overlap neighbor list
2507 : !> \param para_env global parallel environment
2508 : !> \param success true if the metric comparison was performed
2509 : !> \param max_subspace_deviation largest deviation of singular values from one
2510 : !> \param min_svalue smallest singular value of C_ref^+ S(k) C_candidate
2511 : ! **************************************************************************************************
2512 2532 : SUBROUTINE measure_wannier90_subspace_error(ref_cfm, cand_cfm, matrix_s, &
2513 : xkp, cell_to_index, sab_nl, para_env, success, &
2514 : max_subspace_deviation, min_svalue)
2515 : TYPE(cp_cfm_type), INTENT(IN) :: ref_cfm, cand_cfm
2516 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
2517 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
2518 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2519 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2520 : POINTER :: sab_nl
2521 : TYPE(mp_para_env_type), POINTER :: para_env
2522 : LOGICAL, INTENT(OUT) :: success
2523 : REAL(KIND=dp), INTENT(OUT) :: max_subspace_deviation, min_svalue
2524 :
2525 2532 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: metric_projected, metric_vectors, &
2526 2532 : overlap, ref_coeff, s_cand_coeff
2527 : INTEGER :: ib, nao, nmo, nmo_candidate
2528 : REAL(KIND=dp) :: singular_value
2529 2532 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: metric_values
2530 : TYPE(cp_cfm_type) :: s_cand
2531 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_metric
2532 :
2533 2532 : success = .FALSE.
2534 2532 : max_subspace_deviation = HUGE(1.0_dp)
2535 2532 : min_svalue = 0.0_dp
2536 2532 : NULLIFY (matrix_struct_metric)
2537 :
2538 : CALL cp_cfm_get_info(ref_cfm, nrow_global=nao, ncol_global=nmo, &
2539 2532 : matrix_struct=matrix_struct_metric)
2540 2532 : CALL cp_cfm_get_info(cand_cfm, ncol_global=nmo_candidate)
2541 2532 : IF (nmo_candidate /= nmo) RETURN
2542 :
2543 2532 : CALL cp_cfm_create(s_cand, matrix_struct_metric)
2544 : CALL apply_wannier90_kp_matrix(matrix_s, 1, xkp, cell_to_index, sab_nl, &
2545 2532 : cand_cfm, s_cand)
2546 :
2547 15192 : ALLOCATE (ref_coeff(nao, nmo), s_cand_coeff(nao, nmo))
2548 2532 : CALL cp_cfm_get_submatrix(ref_cfm, ref_coeff)
2549 2532 : CALL cp_cfm_get_submatrix(s_cand, s_cand_coeff)
2550 :
2551 : ALLOCATE (overlap(nmo, nmo), &
2552 25320 : metric_projected(nmo, nmo), metric_vectors(nmo, nmo), metric_values(nmo))
2553 5311724 : overlap(:, :) = MATMUL(CONJG(TRANSPOSE(ref_coeff)), s_cand_coeff)
2554 196096 : metric_projected(:, :) = MATMUL(CONJG(TRANSPOSE(overlap)), overlap)
2555 87124 : metric_projected(:, :) = 0.5_dp*(metric_projected + CONJG(TRANSPOSE(metric_projected)))
2556 2532 : CALL diag_complex(metric_projected, metric_vectors, metric_values)
2557 :
2558 2532 : min_svalue = HUGE(1.0_dp)
2559 2532 : max_subspace_deviation = 0.0_dp
2560 11096 : DO ib = 1, nmo
2561 8564 : singular_value = SQRT(MAX(metric_values(ib), 0.0_dp))
2562 8564 : min_svalue = MIN(min_svalue, singular_value)
2563 11096 : max_subspace_deviation = MAX(max_subspace_deviation, ABS(singular_value - 1.0_dp))
2564 : END DO
2565 2532 : CALL para_env%max(max_subspace_deviation)
2566 2532 : CALL para_env%min(min_svalue)
2567 2532 : success = .TRUE.
2568 :
2569 2532 : DEALLOCATE (ref_coeff, s_cand_coeff, overlap, metric_projected, metric_vectors, metric_values)
2570 2532 : CALL cp_cfm_release(s_cand)
2571 :
2572 7596 : END SUBROUTINE measure_wannier90_subspace_error
2573 :
2574 : ! **************************************************************************************************
2575 : !> \brief Measure whether transformed MOs are an H(k),S(k) invariant eigenspace.
2576 : !> \param cand_cfm candidate MO coefficients
2577 : !> \param matrix_s real-space overlap matrix
2578 : !> \param matrix_ks real-space Kohn-Sham matrix
2579 : !> \param xkp target k-point coordinate
2580 : !> \param cell_to_index real-space cell index table
2581 : !> \param sab_nl overlap neighbor list
2582 : !> \param para_env global parallel environment
2583 : !> \param ispin spin index
2584 : !> \param eigenvalues source MO eigenvalues corresponding to the candidate columns
2585 : !> \param success true if the metric and residual checks were performed
2586 : !> \param metric_deviation largest deviation of eigenvalues of C^+ S(k) C from one
2587 : !> \param min_metric_eigenvalue smallest eigenvalue of C^+ S(k) C
2588 : !> \param residual_norm largest element of H(k) C - S(k) C eps
2589 : ! **************************************************************************************************
2590 2532 : SUBROUTINE measure_wannier90_eigenspace_quality(cand_cfm, matrix_s, matrix_ks, xkp, &
2591 2532 : cell_to_index, sab_nl, para_env, ispin, eigenvalues, &
2592 : success, metric_deviation, min_metric_eigenvalue, &
2593 : residual_norm)
2594 : TYPE(cp_cfm_type), INTENT(IN) :: cand_cfm
2595 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_ks
2596 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
2597 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2598 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2599 : POINTER :: sab_nl
2600 : TYPE(mp_para_env_type), POINTER :: para_env
2601 : INTEGER, INTENT(IN) :: ispin
2602 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: eigenvalues
2603 : LOGICAL, INTENT(OUT) :: success
2604 : REAL(KIND=dp), INTENT(OUT) :: metric_deviation, min_metric_eigenvalue, &
2605 : residual_norm
2606 :
2607 2532 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: h_coeff, metric_vectors, residual_block, &
2608 2532 : s_coeff, s_projected
2609 : INTEGER :: ib, nao, nmo
2610 2532 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: metric_values
2611 : TYPE(cp_cfm_type) :: h_cand, metric_cfm, s_cand
2612 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_metric, &
2613 : matrix_struct_projected
2614 :
2615 2532 : success = .FALSE.
2616 2532 : metric_deviation = HUGE(1.0_dp)
2617 2532 : min_metric_eigenvalue = 0.0_dp
2618 2532 : residual_norm = HUGE(1.0_dp)
2619 2532 : NULLIFY (matrix_struct_metric, matrix_struct_projected)
2620 :
2621 : CALL cp_cfm_get_info(cand_cfm, nrow_global=nao, ncol_global=nmo, &
2622 2532 : matrix_struct=matrix_struct_metric)
2623 2532 : IF (SIZE(eigenvalues) < nmo) RETURN
2624 :
2625 2532 : CALL cp_cfm_create(s_cand, matrix_struct_metric)
2626 2532 : CALL cp_cfm_create(h_cand, matrix_struct_metric)
2627 :
2628 : CALL apply_wannier90_kp_matrix(matrix_s, 1, xkp, cell_to_index, sab_nl, &
2629 2532 : cand_cfm, s_cand)
2630 : CALL apply_wannier90_kp_matrix(matrix_ks, ispin, xkp, cell_to_index, sab_nl, &
2631 2532 : cand_cfm, h_cand)
2632 :
2633 : ALLOCATE (h_coeff(nao, nmo), metric_vectors(nmo, nmo), &
2634 : residual_block(nao, nmo), s_coeff(nao, nmo), s_projected(nmo, nmo), &
2635 37980 : metric_values(nmo))
2636 2532 : CALL cp_cfm_get_submatrix(s_cand, s_coeff)
2637 2532 : CALL cp_cfm_get_submatrix(h_cand, h_coeff)
2638 : CALL cp_fm_struct_create(matrix_struct_projected, nrow_global=nmo, ncol_global=nmo, &
2639 : para_env=matrix_struct_metric%para_env, &
2640 2532 : context=matrix_struct_metric%context)
2641 2532 : CALL cp_cfm_create(metric_cfm, matrix_struct_projected)
2642 : CALL parallel_gemm("C", "N", nmo, nmo, nao, CMPLX(1.0_dp, 0.0_dp, KIND=dp), cand_cfm, &
2643 2532 : s_cand, CMPLX(0.0_dp, 0.0_dp, KIND=dp), metric_cfm)
2644 2532 : CALL cp_cfm_get_submatrix(metric_cfm, s_projected)
2645 87124 : s_projected(:, :) = 0.5_dp*(s_projected + CONJG(TRANSPOSE(s_projected)))
2646 2532 : CALL diag_complex(s_projected, metric_vectors, metric_values)
2647 11096 : metric_deviation = MAXVAL(ABS(metric_values - 1.0_dp))
2648 11096 : min_metric_eigenvalue = MINVAL(metric_values)
2649 :
2650 309164 : residual_block(:, :) = h_coeff
2651 11096 : DO ib = 1, nmo
2652 309164 : residual_block(:, ib) = residual_block(:, ib) - eigenvalues(ib)*s_coeff(:, ib)
2653 : END DO
2654 309164 : residual_norm = MAXVAL(ABS(residual_block))
2655 2532 : CALL para_env%max(metric_deviation)
2656 2532 : CALL para_env%min(min_metric_eigenvalue)
2657 2532 : CALL para_env%max(residual_norm)
2658 2532 : success = .TRUE.
2659 :
2660 2532 : DEALLOCATE (h_coeff, metric_vectors, residual_block, s_coeff, s_projected, metric_values)
2661 2532 : CALL cp_cfm_release(s_cand)
2662 2532 : CALL cp_cfm_release(h_cand)
2663 2532 : CALL cp_cfm_release(metric_cfm)
2664 2532 : CALL cp_fm_struct_release(matrix_struct_projected)
2665 :
2666 10128 : END SUBROUTINE measure_wannier90_eigenspace_quality
2667 :
2668 : ! **************************************************************************************************
2669 : !> \brief Write fixed trial projections in the exported MO gauge for external Wannier90.
2670 : !> \param amn Band/trial/k-point projection matrices.
2671 : !> \param para_env Global parallel environment.
2672 : !> \param seed_name Output filename prefix.
2673 : ! **************************************************************************************************
2674 24 : SUBROUTINE write_wannier90_amn(amn, para_env, seed_name)
2675 : COMPLEX(KIND=dp), INTENT(IN) :: amn(:, :, :)
2676 : TYPE(mp_para_env_type), POINTER :: para_env
2677 : CHARACTER(LEN=*), INTENT(IN) :: seed_name
2678 :
2679 : INTEGER :: ib, ik, iunit, iwann
2680 :
2681 24 : IF (.NOT. para_env%is_source()) RETURN
2682 12 : CALL open_file(TRIM(seed_name)//'.amn', unit_number=iunit, file_status='REPLACE', file_action='WRITE')
2683 12 : WRITE (iunit, '(A)') '! Initial Wannier90 projections generated by CP2K'
2684 12 : WRITE (iunit, '(3I10)') SIZE(amn, 1), SIZE(amn, 3), SIZE(amn, 2)
2685 220 : DO ik = 1, SIZE(amn, 3)
2686 924 : DO iwann = 1, SIZE(amn, 2)
2687 3792 : DO ib = 1, SIZE(amn, 1)
2688 2880 : WRITE (iunit, '(3I10,2ES24.16)') ib, iwann, ik, REAL(amn(ib, iwann, ik), dp), &
2689 6464 : AIMAG(amn(ib, iwann, ik))
2690 : END DO
2691 : END DO
2692 : END DO
2693 12 : CALL close_file(iunit)
2694 : END SUBROUTINE write_wannier90_amn
2695 :
2696 : ! **************************************************************************************************
2697 : !> \brief Enumerate all on-atom radial s/p shell pairs in global spherical-AO ordering.
2698 : !> \param qs_env Electronic environment providing the basis and atom ordering.
2699 : !> \param groups Global indices of an s function and the three functions of a p shell.
2700 : ! **************************************************************************************************
2701 4 : SUBROUTINE build_wannier90_hybrid_groups(qs_env, groups)
2702 : TYPE(qs_environment_type), POINTER :: qs_env
2703 : INTEGER, ALLOCATABLE, INTENT(OUT) :: groups(:, :)
2704 :
2705 : INTEGER :: iatom, ig, ikind, ip, is, iset, ishell, &
2706 : n, ng, np, ns, nset, nsgf, offset, pass
2707 4 : INTEGER, ALLOCATABLE :: p_shells(:), s_shells(:)
2708 4 : INTEGER, POINTER :: first_sgf(:, :), l(:, :), nshell(:)
2709 : TYPE(gto_basis_set_type), POINTER :: basis
2710 4 : TYPE(particle_type), POINTER :: particles(:)
2711 4 : TYPE(qs_kind_type), POINTER :: qs_kinds(:)
2712 :
2713 4 : CALL get_qs_env(qs_env, particle_set=particles, qs_kind_set=qs_kinds)
2714 4 : ng = 0
2715 12 : DO pass = 1, 2
2716 8 : offset = 0
2717 8 : ig = 0
2718 24 : DO iatom = 1, SIZE(particles)
2719 16 : CALL get_atomic_kind(particles(iatom)%atomic_kind, kind_number=ikind)
2720 16 : CALL get_qs_kind(qs_kinds(ikind), basis_set=basis, nsgf=nsgf)
2721 16 : IF (ASSOCIATED(basis)) THEN
2722 16 : CALL get_gto_basis_set(basis, nset=nset, nshell=nshell, l=l, first_sgf=first_sgf)
2723 160 : ALLOCATE (s_shells(SUM(nshell)), p_shells(SUM(nshell)))
2724 16 : ns = 0
2725 16 : np = 0
2726 48 : DO iset = 1, nset
2727 96 : DO ishell = 1, nshell(iset)
2728 48 : n = offset + first_sgf(ishell, iset)
2729 32 : SELECT CASE (l(ishell, iset))
2730 : CASE (0)
2731 32 : ns = ns + 1
2732 32 : s_shells(ns) = n
2733 : CASE (1)
2734 16 : np = np + 1
2735 48 : p_shells(np) = n
2736 : END SELECT
2737 : END DO
2738 : END DO
2739 16 : IF (pass == 1) THEN
2740 8 : ng = ng + ns*np
2741 : ELSE
2742 24 : DO is = 1, ns
2743 40 : DO ip = 1, np
2744 16 : ig = ig + 1
2745 16 : n = p_shells(ip)
2746 96 : groups(:, ig) = [s_shells(is), n, n + 1, n + 2]
2747 : END DO
2748 : END DO
2749 : END IF
2750 32 : DEALLOCATE (s_shells, p_shells)
2751 : END IF
2752 40 : offset = offset + nsgf
2753 : END DO
2754 12 : IF (pass == 1) THEN
2755 12 : ALLOCATE (groups(4, ng))
2756 : END IF
2757 : END DO
2758 8 : END SUBROUTINE build_wannier90_hybrid_groups
2759 :
2760 : ! **************************************************************************************************
2761 : !> \brief Form C(k)^H S(k) projections on all AO trials using the exported MO coefficients.
2762 : !> \param kpoint Full export mesh and its distributed MO coefficients.
2763 : !> \param matrix_s Real-space overlap matrices.
2764 : !> \param sab_nl Overlap neighbour lists.
2765 : !> \param para_env Global parallel environment.
2766 : !> \param ispin Selected spin channel.
2767 : !> \param bands Original MO indices retained in the export.
2768 : !> \param apply_phase Whether the overlap export uses the additional global Bloch gauge.
2769 : !> \param phase_center Fractional origin for that gauge.
2770 : !> \param projections Band/AO/k-point projections, stored only on the source rank.
2771 : ! **************************************************************************************************
2772 10 : SUBROUTINE build_wannier90_ao_projections(kpoint, matrix_s, sab_nl, para_env, ispin, bands, &
2773 : apply_phase, phase_center, projections)
2774 : TYPE(kpoint_type), POINTER :: kpoint
2775 : TYPE(dbcsr_p_type), POINTER :: matrix_s(:, :)
2776 : TYPE(neighbor_list_set_p_type), POINTER :: sab_nl(:)
2777 : TYPE(mp_para_env_type), POINTER :: para_env
2778 : INTEGER, INTENT(IN) :: ispin, bands(:)
2779 : LOGICAL, INTENT(IN) :: apply_phase
2780 : REAL(KIND=dp), INTENT(IN) :: phase_center(3)
2781 : COMPLEX(KIND=dp), ALLOCATABLE, INTENT(OUT) :: projections(:, :, :)
2782 :
2783 : COMPLEX(KIND=dp), ALLOCATABLE :: buffer(:, :)
2784 : INTEGER :: ik, iklocal, nao, nk, nmo
2785 10 : INTEGER, POINTER :: cell_to_index(:, :, :)
2786 : REAL(KIND=dp) :: phase
2787 : TYPE(cp_cfm_type) :: mo_cfm, sc_cfm
2788 : TYPE(cp_cfm_type), POINTER :: source_cfm
2789 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
2790 :
2791 10 : NULLIFY (matrix_struct, cell_to_index, source_cfm)
2792 10 : CALL get_mo_set(kpoint%kp_env(1)%kpoint_env%mos(ispin), nao=nao, nmo=nmo)
2793 10 : nk = SIZE(kpoint%xkp, 2)
2794 10 : CALL get_kpoint_info(kpoint, cell_to_index=cell_to_index)
2795 : CALL cp_fm_struct_create(matrix_struct, nrow_global=nao, ncol_global=nmo, &
2796 10 : para_env=para_env, context=kpoint%blacs_env)
2797 10 : CALL cp_cfm_create(mo_cfm, matrix_struct)
2798 10 : CALL cp_cfm_create(sc_cfm, matrix_struct)
2799 40 : ALLOCATE (buffer(nmo, nao))
2800 10 : IF (para_env%is_source()) THEN
2801 25 : ALLOCATE (projections(SIZE(bands), nao, nk))
2802 : ELSE
2803 5 : ALLOCATE (projections(0, 0, 0))
2804 : END IF
2805 90 : DO ik = 1, nk
2806 80 : IF (ik >= kpoint%kp_range(1) .AND. ik <= kpoint%kp_range(2)) THEN
2807 80 : iklocal = ik - kpoint%kp_range(1) + 1
2808 80 : source_cfm => kpoint%kp_env(iklocal)%kpoint_env%mos(ispin)%cmo_coeff
2809 : ELSE
2810 0 : NULLIFY (source_cfm)
2811 : END IF
2812 80 : CALL wannier90_copy_cfm_general(source_cfm, mo_cfm, para_env)
2813 : CALL apply_wannier90_kp_matrix(matrix_s, 1, kpoint%xkp(:, ik), cell_to_index, sab_nl, &
2814 80 : mo_cfm, sc_cfm)
2815 80 : CALL cp_cfm_get_submatrix(sc_cfm, buffer, transpose=.TRUE.)
2816 90 : IF (para_env%is_source()) THEN
2817 1008 : projections(:, :, ik) = CONJG(buffer(bands, :))
2818 40 : IF (apply_phase) THEN
2819 : ! M'(k,b) = exp(i[phi(k+b)-phi(k)]) M(k,b) requires A'(k) = exp(-i phi(k)) A(k).
2820 128 : phase = -twopi*DOT_PRODUCT(kpoint%xkp(:, ik), phase_center)
2821 760 : projections(:, :, ik) = projections(:, :, ik)*CMPLX(COS(phase), SIN(phase), KIND=dp)
2822 : END IF
2823 : END IF
2824 : END DO
2825 10 : DEALLOCATE (buffer)
2826 10 : CALL cp_cfm_release(mo_cfm)
2827 10 : CALL cp_cfm_release(sc_cfm)
2828 10 : CALL cp_fm_struct_release(matrix_struct)
2829 20 : END SUBROUTINE build_wannier90_ao_projections
2830 :
2831 : ! **************************************************************************************************
2832 : !> \brief Assemble the non-Hermitian inter-k-point AO overlap from its real-space images.
2833 : !> \param cosmat Real part of <mu,0|exp(-i b.r)|nu,R> for the stored direction of b.
2834 : !> \param sinmat Imaginary part of the same integrals.
2835 : !> \param bsign Sign of b relative to the stored integrals; cosine is even and sine is odd.
2836 : !> \param kleft Reduced left k-point.
2837 : !> \param kright Reduced right k-point, equal to kleft+b modulo a reciprocal lattice vector.
2838 : !> \param cell_to_index Real-space image index table.
2839 : !> \param sab_nl Symmetric orbital neighbour list.
2840 : !> \param real_matrix Real part of the general AO overlap matrix.
2841 : !> \param imag_matrix Imaginary part of the general AO overlap matrix.
2842 : ! **************************************************************************************************
2843 3552 : SUBROUTINE transform_wannier90_overlap(cosmat, sinmat, bsign, kleft, kright, cell_to_index, &
2844 : sab_nl, real_matrix, imag_matrix)
2845 : TYPE(dbcsr_p_type), POINTER :: cosmat(:, :), sinmat(:, :)
2846 : REAL(KIND=dp), INTENT(IN) :: bsign, kleft(3), kright(3)
2847 : INTEGER, POINTER :: cell_to_index(:, :, :)
2848 : TYPE(neighbor_list_set_p_type), POINTER :: sab_nl(:)
2849 : TYPE(dbcsr_type), INTENT(INOUT) :: real_matrix, imag_matrix
2850 :
2851 : INTEGER :: direction, iatom, ic, icol, &
2852 : image_cell(3), irow, jatom
2853 : LOGICAL :: found, symmetric, transposed
2854 : REAL(KIND=dp) :: ci, cr, phase
2855 3552 : REAL(KIND=dp), POINTER :: cb(:, :), ib(:, :), rb(:, :), sb(:, :)
2856 3552 : TYPE(neighbor_list_iterator_p_type), POINTER :: iterator(:)
2857 :
2858 3552 : CALL get_neighbor_list_set_p(sab_nl, symmetric=symmetric)
2859 3552 : CPASSERT(symmetric)
2860 3552 : CALL dbcsr_set(real_matrix, 0.0_dp)
2861 3552 : CALL dbcsr_set(imag_matrix, 0.0_dp)
2862 3552 : CALL neighbor_list_iterator_create(iterator, sab_nl)
2863 3729408 : DO WHILE (neighbor_list_iterate(iterator) == 0)
2864 3725856 : CALL get_iterator_info(iterator, iatom=iatom, jatom=jatom, cell=image_cell)
2865 3725856 : ic = cell_to_index(image_cell(1), image_cell(2), image_cell(3))
2866 3725856 : CPASSERT(ic > 0)
2867 3725856 : irow = MIN(iatom, jatom)
2868 3725856 : icol = MAX(iatom, jatom)
2869 3725856 : CALL dbcsr_get_block_p(cosmat(1, ic)%matrix, irow, icol, cb, found=found)
2870 3725856 : CPASSERT(found)
2871 3725856 : CALL dbcsr_get_block_p(sinmat(1, ic)%matrix, irow, icol, sb, found=found)
2872 3725856 : CPASSERT(found)
2873 11181120 : DO direction = 1, 2
2874 7451712 : IF (direction == 1) THEN
2875 3725856 : irow = iatom
2876 3725856 : icol = jatom
2877 3725856 : transposed = iatom > jatom
2878 14903424 : phase = twopi*DOT_PRODUCT(kright, REAL(image_cell, KIND=dp))
2879 : ELSE
2880 : ! B_nu,mu(-R,b) = exp(i b.R) B_mu,nu(R,b)^T, not its Hermitian transpose.
2881 : ! Self-pair images already occur in both directions in the neighbour list.
2882 3725856 : IF (iatom == jatom) CYCLE
2883 2786304 : irow = jatom
2884 2786304 : icol = iatom
2885 2786304 : transposed = iatom < jatom
2886 11145216 : phase = -twopi*DOT_PRODUCT(kleft, REAL(image_cell, KIND=dp))
2887 : END IF
2888 6512160 : CALL dbcsr_get_block_p(real_matrix, irow, icol, rb, found=found)
2889 6512160 : CPASSERT(found)
2890 6512160 : CALL dbcsr_get_block_p(imag_matrix, irow, icol, ib, found=found)
2891 6512160 : CPASSERT(found)
2892 6512160 : cr = COS(phase)
2893 6512160 : ci = SIN(phase)
2894 23262336 : IF (transposed) THEN
2895 205685760 : rb(:, :) = rb(:, :) + cr*TRANSPOSE(cb) - bsign*ci*TRANSPOSE(sb)
2896 205685760 : ib(:, :) = ib(:, :) + ci*TRANSPOSE(cb) + bsign*cr*TRANSPOSE(sb)
2897 : ELSE
2898 290364576 : rb(:, :) = rb(:, :) + cr*cb(:, :) - bsign*ci*sb(:, :)
2899 290364576 : ib(:, :) = ib(:, :) + ci*cb(:, :) + bsign*cr*sb(:, :)
2900 : END IF
2901 : END DO
2902 : END DO
2903 3552 : CALL neighbor_list_iterator_release(iterator)
2904 3552 : END SUBROUTINE transform_wannier90_overlap
2905 :
2906 : ! **************************************************************************************************
2907 : !> \brief Apply a complex k-point matrix to a complex MO coefficient matrix.
2908 : !> \param rsmat real-space matrix images
2909 : !> \param ispin spin index for rsmat
2910 : !> \param xkp target k-point coordinate
2911 : !> \param cell_to_index real-space cell index table
2912 : !> \param sab_nl overlap neighbor list
2913 : !> \param coeff input MO coefficients
2914 : !> \param RESULT output of applying the k-point matrix to the MO coefficients
2915 : ! **************************************************************************************************
2916 34960 : SUBROUTINE apply_wannier90_kp_matrix(rsmat, ispin, xkp, cell_to_index, sab_nl, &
2917 : coeff, RESULT)
2918 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rsmat
2919 : INTEGER, INTENT(IN) :: ispin
2920 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
2921 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2922 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2923 : POINTER :: sab_nl
2924 : TYPE(cp_cfm_type), INTENT(IN) :: coeff
2925 : TYPE(cp_cfm_type), INTENT(INOUT) :: result
2926 :
2927 : INTEGER :: nao, ncol
2928 : TYPE(cp_cfm_type) :: kmat_cfm
2929 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_ao, matrix_struct_coeff
2930 : TYPE(cp_fm_type) :: mat_imag, mat_real
2931 : TYPE(dbcsr_type), POINTER :: kmat_imag, kmat_imag_full, kmat_real, &
2932 : kmat_real_full
2933 :
2934 8740 : NULLIFY (matrix_struct_ao, matrix_struct_coeff, kmat_imag, kmat_imag_full, kmat_real, &
2935 8740 : kmat_real_full)
2936 :
2937 : CALL cp_cfm_get_info(coeff, nrow_global=nao, ncol_global=ncol, &
2938 8740 : matrix_struct=matrix_struct_coeff)
2939 :
2940 8740 : ALLOCATE (kmat_real, kmat_imag, kmat_real_full, kmat_imag_full)
2941 : CALL dbcsr_create(kmat_real, template=rsmat(ispin, 1)%matrix, &
2942 8740 : matrix_type=dbcsr_type_symmetric)
2943 : CALL dbcsr_create(kmat_imag, template=rsmat(ispin, 1)%matrix, &
2944 8740 : matrix_type=dbcsr_type_antisymmetric)
2945 : CALL dbcsr_create(kmat_real_full, template=rsmat(ispin, 1)%matrix, &
2946 8740 : matrix_type=dbcsr_type_no_symmetry)
2947 : CALL dbcsr_create(kmat_imag_full, template=rsmat(ispin, 1)%matrix, &
2948 8740 : matrix_type=dbcsr_type_no_symmetry)
2949 8740 : CALL cp_dbcsr_alloc_block_from_nbl(kmat_real, sab_nl)
2950 8740 : CALL cp_dbcsr_alloc_block_from_nbl(kmat_imag, sab_nl)
2951 8740 : CALL dbcsr_set(kmat_real, 0.0_dp)
2952 8740 : CALL dbcsr_set(kmat_imag, 0.0_dp)
2953 : CALL rskp_transform(kmat_real, kmat_imag, rsmat=rsmat, ispin=ispin, &
2954 8740 : xkp=xkp, cell_to_index=cell_to_index, sab_nl=sab_nl)
2955 8740 : CALL dbcsr_desymmetrize(kmat_real, kmat_real_full)
2956 8740 : CALL dbcsr_desymmetrize(kmat_imag, kmat_imag_full)
2957 :
2958 : CALL cp_fm_struct_create(matrix_struct_ao, nrow_global=nao, ncol_global=nao, &
2959 : para_env=matrix_struct_coeff%para_env, &
2960 8740 : context=matrix_struct_coeff%context)
2961 8740 : CALL cp_fm_create(mat_real, matrix_struct_ao)
2962 8740 : CALL cp_fm_create(mat_imag, matrix_struct_ao)
2963 8740 : CALL copy_dbcsr_to_fm(kmat_real_full, mat_real)
2964 8740 : CALL copy_dbcsr_to_fm(kmat_imag_full, mat_imag)
2965 :
2966 8740 : CALL cp_cfm_create(kmat_cfm, matrix_struct_ao)
2967 8740 : CALL cp_fm_to_cfm(mat_real, mat_imag, kmat_cfm)
2968 : CALL parallel_gemm("N", "N", nao, ncol, nao, CMPLX(1.0_dp, 0.0_dp, KIND=dp), kmat_cfm, &
2969 8740 : coeff, CMPLX(0.0_dp, 0.0_dp, KIND=dp), RESULT)
2970 :
2971 8740 : CALL cp_fm_release(mat_real)
2972 8740 : CALL cp_fm_release(mat_imag)
2973 8740 : CALL cp_cfm_release(kmat_cfm)
2974 8740 : CALL cp_fm_struct_release(matrix_struct_ao)
2975 8740 : CALL dbcsr_deallocate_matrix(kmat_real)
2976 8740 : CALL dbcsr_deallocate_matrix(kmat_imag)
2977 8740 : CALL dbcsr_deallocate_matrix(kmat_real_full)
2978 8740 : CALL dbcsr_deallocate_matrix(kmat_imag_full)
2979 :
2980 8740 : END SUBROUTINE apply_wannier90_kp_matrix
2981 :
2982 : ! **************************************************************************************************
2983 : !> \brief Rayleigh-Ritz stabilize a symmetry-reconstructed Wannier90 MO subspace.
2984 : !> \param dst_cfm transformed MO coefficients
2985 : !> \param matrix_s real-space overlap matrix
2986 : !> \param matrix_ks real-space Kohn-Sham matrix
2987 : !> \param xkp target k-point coordinate
2988 : !> \param cell_to_index real-space cell index table
2989 : !> \param sab_nl overlap neighbor list
2990 : !> \param ispin spin index
2991 : !> \param eigenvalues Ritz eigenvalues of the stabilized subspace
2992 : !> \param degenerate_band_tol degeneracy threshold
2993 : !> \param success true if the subspace was stabilized
2994 : !> \param reason diagnostic message
2995 : !> \param aligned_blocks number of stabilized subspaces
2996 : !> \param aligned_max_size largest stabilized subspace
2997 : !> \param aligned_min_svalue smallest S(k)-metric eigenvalue
2998 : !> \param max_residual largest Ritz residual
2999 : ! **************************************************************************************************
3000 532 : SUBROUTINE ritz_stabilize_wannier90_subspace(dst_cfm, matrix_s, matrix_ks, &
3001 1064 : xkp, cell_to_index, sab_nl, ispin, eigenvalues, &
3002 : degenerate_band_tol, success, reason, aligned_blocks, &
3003 : aligned_max_size, aligned_min_svalue, max_residual)
3004 : TYPE(cp_cfm_type), INTENT(IN) :: dst_cfm
3005 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_ks
3006 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
3007 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3008 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3009 : POINTER :: sab_nl
3010 : INTEGER, INTENT(IN) :: ispin
3011 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: eigenvalues
3012 : REAL(KIND=dp), INTENT(IN) :: degenerate_band_tol
3013 : LOGICAL, INTENT(OUT) :: success
3014 : CHARACTER(LEN=*), INTENT(OUT) :: reason
3015 : INTEGER, INTENT(OUT) :: aligned_blocks, aligned_max_size
3016 : REAL(KIND=dp), INTENT(OUT) :: aligned_min_svalue, max_residual
3017 :
3018 : REAL(KIND=dp), PARAMETER :: residual_tol = 1.0e-2_dp
3019 :
3020 532 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: block_coeff, dst_coeff, h_block, h_coeff, &
3021 532 : h_projected, h_projected_work, metric_vectors, residual_block, ritz_vectors, s_block, &
3022 532 : s_coeff, s_projected, stabilized
3023 : INTEGER :: block_first, block_last, block_size, ib, &
3024 : nao, nmo
3025 : REAL(KIND=dp) :: metric_deviation, norm_value, &
3026 : residual_norm
3027 532 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: metric_values, ritz_values
3028 : TYPE(cp_cfm_type) :: h_dst, s_dst
3029 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_metric
3030 :
3031 532 : success = .FALSE.
3032 532 : reason = ""
3033 532 : aligned_blocks = 0
3034 532 : aligned_max_size = 0
3035 532 : aligned_min_svalue = HUGE(1.0_dp)
3036 532 : max_residual = 0.0_dp
3037 :
3038 532 : NULLIFY (matrix_struct_metric)
3039 : CALL cp_cfm_get_info(dst_cfm, nrow_global=nao, ncol_global=nmo, &
3040 532 : matrix_struct=matrix_struct_metric)
3041 532 : IF (SIZE(eigenvalues) < nmo) THEN
3042 0 : reason = "not enough eigenvalues for Wannier90 Ritz subspace stabilization"
3043 : RETURN
3044 : END IF
3045 :
3046 532 : CALL cp_cfm_create(s_dst, matrix_struct_metric)
3047 532 : CALL cp_cfm_create(h_dst, matrix_struct_metric)
3048 :
3049 : stabilize_subspace: BLOCK
3050 : CALL apply_wannier90_kp_matrix(matrix_s, 1, xkp, cell_to_index, sab_nl, &
3051 532 : dst_cfm, s_dst)
3052 : CALL apply_wannier90_kp_matrix(matrix_ks, ispin, xkp, cell_to_index, sab_nl, &
3053 532 : dst_cfm, h_dst)
3054 :
3055 4256 : ALLOCATE (dst_coeff(nao, nmo), s_coeff(nao, nmo), h_coeff(nao, nmo))
3056 532 : CALL cp_cfm_get_submatrix(dst_cfm, dst_coeff)
3057 532 : CALL cp_cfm_get_submatrix(s_dst, s_coeff)
3058 532 : CALL cp_cfm_get_submatrix(h_dst, h_coeff)
3059 :
3060 532 : block_first = 1
3061 1916 : DO WHILE (block_first <= nmo)
3062 : block_last = block_first
3063 2096 : DO WHILE (block_last < nmo)
3064 1564 : IF (ABS(eigenvalues(block_last + 1) - eigenvalues(block_last)) >= degenerate_band_tol) EXIT
3065 1384 : block_last = block_last + 1
3066 : END DO
3067 1384 : block_size = block_last - block_first + 1
3068 1384 : IF (block_size > 1) THEN
3069 : ! The atom/AO operation fixes the subspace, while the little-group gauge inside an
3070 : ! exactly degenerate manifold is arbitrary. Stabilize only that manifold and verify
3071 : ! that it is an invariant H(k),S(k) subspace before exporting it to Wannier90.
3072 0 : ALLOCATE (block_coeff(nao, block_size), h_block(nao, block_size), &
3073 0 : h_projected(block_size, block_size), h_projected_work(block_size, block_size), &
3074 0 : metric_vectors(block_size, block_size), residual_block(nao, block_size), &
3075 0 : ritz_vectors(block_size, block_size), s_block(nao, block_size), &
3076 0 : s_projected(block_size, block_size), stabilized(nao, block_size), &
3077 12688 : metric_values(block_size), ritz_values(block_size))
3078 63192 : block_coeff(:, :) = dst_coeff(:, block_first:block_last)
3079 63192 : s_block(:, :) = s_coeff(:, block_first:block_last)
3080 63192 : h_block(:, :) = h_coeff(:, block_first:block_last)
3081 977816 : s_projected(:, :) = MATMUL(CONJG(TRANSPOSE(block_coeff)), s_block)
3082 977816 : h_projected(:, :) = MATMUL(CONJG(TRANSPOSE(block_coeff)), h_block)
3083 9032 : s_projected(:, :) = 0.5_dp*(s_projected + CONJG(TRANSPOSE(s_projected)))
3084 9032 : h_projected(:, :) = 0.5_dp*(h_projected + CONJG(TRANSPOSE(h_projected)))
3085 :
3086 488 : CALL diag_complex(s_projected, metric_vectors, metric_values)
3087 1688 : aligned_min_svalue = MIN(aligned_min_svalue, MINVAL(metric_values))
3088 1688 : metric_deviation = MAXVAL(ABS(metric_values - 1.0_dp))
3089 1688 : IF (MINVAL(metric_values) < 1.0e-10_dp) THEN
3090 : WRITE (reason, "(A,I0,A,ES9.2,A,ES9.2)") &
3091 0 : "singular metric blk=", block_first, " smin=", MINVAL(metric_values), &
3092 0 : " dS=", metric_deviation
3093 0 : max_residual = HUGE(1.0_dp)
3094 0 : EXIT stabilize_subspace
3095 : END IF
3096 :
3097 1688 : DO ib = 1, block_size
3098 4760 : metric_vectors(:, ib) = metric_vectors(:, ib)/SQRT(metric_values(ib))
3099 : END DO
3100 53720 : h_projected_work(:, :) = MATMUL(h_projected, metric_vectors)
3101 53720 : h_projected(:, :) = MATMUL(CONJG(TRANSPOSE(metric_vectors)), h_projected_work)
3102 9032 : h_projected(:, :) = 0.5_dp*(h_projected + CONJG(TRANSPOSE(h_projected)))
3103 488 : CALL diag_complex(h_projected, ritz_vectors, ritz_values)
3104 53720 : h_projected_work(:, :) = MATMUL(metric_vectors, ritz_vectors)
3105 4760 : ritz_vectors(:, :) = h_projected_work
3106 653464 : stabilized(:, :) = MATMUL(block_coeff, ritz_vectors)
3107 653464 : residual_block(:, :) = MATMUL(h_block, ritz_vectors)
3108 63192 : h_block(:, :) = residual_block
3109 653464 : residual_block(:, :) = MATMUL(s_block, ritz_vectors)
3110 63192 : s_block(:, :) = residual_block
3111 1688 : DO ib = 1, block_size
3112 62704 : norm_value = SQRT(ABS(REAL(DOT_PRODUCT(stabilized(:, ib), s_block(:, ib)), KIND=dp)))
3113 1688 : IF (norm_value > EPSILON(1.0_dp)) THEN
3114 62704 : stabilized(:, ib) = stabilized(:, ib)/norm_value
3115 62704 : h_block(:, ib) = h_block(:, ib)/norm_value
3116 62704 : s_block(:, ib) = s_block(:, ib)/norm_value
3117 : END IF
3118 : END DO
3119 63192 : residual_block(:, :) = h_block
3120 1688 : DO ib = 1, block_size
3121 : residual_block(:, ib) = residual_block(:, ib) - &
3122 63192 : eigenvalues(block_first + ib - 1)*s_block(:, ib)
3123 : END DO
3124 63192 : residual_norm = MAXVAL(ABS(residual_block))
3125 488 : max_residual = MAX(max_residual, residual_norm)
3126 488 : IF (residual_norm > residual_tol) THEN
3127 : WRITE (reason, "(A,I0,A,ES9.2)") &
3128 0 : "blk=", block_first, " dS=", metric_deviation
3129 0 : EXIT stabilize_subspace
3130 : END IF
3131 :
3132 63192 : dst_coeff(:, block_first:block_last) = stabilized
3133 63192 : h_coeff(:, block_first:block_last) = h_block
3134 63192 : s_coeff(:, block_first:block_last) = s_block
3135 488 : aligned_blocks = aligned_blocks + 1
3136 488 : aligned_max_size = MAX(aligned_max_size, block_size)
3137 0 : DEALLOCATE (block_coeff, h_block, h_projected, h_projected_work, metric_vectors, &
3138 0 : residual_block, ritz_vectors, s_block, s_projected, stabilized, metric_values, &
3139 488 : ritz_values)
3140 : END IF
3141 1384 : block_first = block_last + 1
3142 : END DO
3143 :
3144 2628 : DO ib = 1, nmo
3145 93744 : residual_norm = MAXVAL(ABS(h_coeff(:, ib) - eigenvalues(ib)*s_coeff(:, ib)))
3146 2628 : max_residual = MAX(max_residual, residual_norm)
3147 : END DO
3148 532 : IF (max_residual > residual_tol) THEN
3149 : WRITE (reason, "(A,ES10.3)") &
3150 0 : "atom/AO W90 reuse guarded: Ritz residual=", max_residual
3151 0 : EXIT stabilize_subspace
3152 : END IF
3153 :
3154 532 : CALL cp_cfm_set_submatrix(dst_cfm, dst_coeff)
3155 :
3156 532 : IF (aligned_blocks == 0) aligned_min_svalue = 0.0_dp
3157 532 : success = .TRUE.
3158 : END BLOCK stabilize_subspace
3159 :
3160 532 : CALL cp_cfm_release(s_dst)
3161 532 : CALL cp_cfm_release(h_dst)
3162 :
3163 1064 : END SUBROUTINE ritz_stabilize_wannier90_subspace
3164 :
3165 : ! **************************************************************************************************
3166 : !> \brief Reconstruct a Wannier90 export window from a larger symmetry-transformed SCF MO space.
3167 : !> \param src_cfm transformed source MO window
3168 : !> \param dst_cfm exported reconstructed MO coefficients
3169 : !> \param matrix_s real-space overlap matrix
3170 : !> \param matrix_ks real-space Kohn-Sham matrix
3171 : !> \param xkp target k-point coordinate
3172 : !> \param cell_to_index real-space cell index table
3173 : !> \param sab_nl overlap neighbor list
3174 : !> \param ispin spin index
3175 : !> \param eigenvalues reconstructed target eigenvalues for the exported window
3176 : !> \param nmo_export number of MOs to export
3177 : !> \param success true if the reconstructed window is an invariant H(k),S(k) subspace
3178 : !> \param reason diagnostic message
3179 : !> \param min_svalue smallest S(k)-metric eigenvalue in the source window
3180 : !> \param max_residual largest target Ritz residual
3181 : ! **************************************************************************************************
3182 0 : SUBROUTINE ritz_reconstruct_wannier90_window(src_cfm, dst_cfm, matrix_s, &
3183 : matrix_ks, xkp, cell_to_index, sab_nl, ispin, &
3184 0 : eigenvalues, nmo_export, success, reason, min_svalue, &
3185 : max_residual)
3186 : TYPE(cp_cfm_type), INTENT(IN) :: src_cfm, dst_cfm
3187 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_ks
3188 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
3189 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3190 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3191 : POINTER :: sab_nl
3192 : INTEGER, INTENT(IN) :: ispin
3193 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: eigenvalues
3194 : INTEGER, INTENT(IN) :: nmo_export
3195 : LOGICAL, INTENT(OUT) :: success
3196 : CHARACTER(LEN=*), INTENT(OUT) :: reason
3197 : REAL(KIND=dp), INTENT(OUT) :: min_svalue, max_residual
3198 :
3199 : REAL(KIND=dp), PARAMETER :: eigenvalue_tol = 1.0e-6_dp, &
3200 : residual_tol = 1.0e-7_dp
3201 :
3202 0 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coeff_work, h_coeff, h_projected, &
3203 0 : h_projected_work, metric_vectors, residual_block, ritz_vectors, s_coeff, s_projected, &
3204 0 : source_coeff, stabilized
3205 : INTEGER :: ib, nao, nmo_source
3206 : REAL(KIND=dp) :: max_eigenvalue_shift, metric_deviation, &
3207 : norm_value
3208 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: metric_values, ritz_values, &
3209 0 : source_eigenvalues
3210 : TYPE(cp_cfm_type) :: h_src, s_src
3211 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_metric
3212 :
3213 0 : success = .FALSE.
3214 0 : reason = ""
3215 0 : min_svalue = HUGE(1.0_dp)
3216 0 : max_residual = HUGE(1.0_dp)
3217 :
3218 0 : NULLIFY (matrix_struct_metric)
3219 : CALL cp_cfm_get_info(src_cfm, nrow_global=nao, ncol_global=nmo_source, &
3220 0 : matrix_struct=matrix_struct_metric)
3221 0 : IF (nmo_export > nmo_source) THEN
3222 0 : reason = "Wannier90 export window is larger than the transformed SCF MO space"
3223 0 : RETURN
3224 : END IF
3225 0 : IF (SIZE(eigenvalues) < nmo_export) THEN
3226 0 : reason = "not enough eigenvalue storage for Wannier90 source-window reconstruction"
3227 : RETURN
3228 : END IF
3229 :
3230 0 : CALL cp_cfm_create(s_src, matrix_struct_metric)
3231 0 : CALL cp_cfm_create(h_src, matrix_struct_metric)
3232 :
3233 : CALL apply_wannier90_kp_matrix(matrix_s, 1, xkp, cell_to_index, sab_nl, &
3234 0 : src_cfm, s_src)
3235 : CALL apply_wannier90_kp_matrix(matrix_ks, ispin, xkp, cell_to_index, sab_nl, &
3236 0 : src_cfm, h_src)
3237 :
3238 0 : ALLOCATE (source_coeff(nao, nmo_source), s_coeff(nao, nmo_source), &
3239 0 : h_coeff(nao, nmo_source), h_projected(nmo_source, nmo_source), &
3240 0 : h_projected_work(nmo_source, nmo_source), metric_vectors(nmo_source, nmo_source), &
3241 0 : residual_block(nao, nmo_export), ritz_vectors(nmo_source, nmo_source), &
3242 0 : s_projected(nmo_source, nmo_source), stabilized(nao, nmo_export), &
3243 0 : coeff_work(nao, nmo_source), metric_values(nmo_source), ritz_values(nmo_source), &
3244 0 : source_eigenvalues(nmo_export))
3245 0 : source_eigenvalues(1:nmo_export) = eigenvalues(1:nmo_export)
3246 0 : CALL cp_cfm_get_submatrix(src_cfm, source_coeff)
3247 0 : CALL cp_cfm_get_submatrix(s_src, s_coeff)
3248 0 : CALL cp_cfm_get_submatrix(h_src, h_coeff)
3249 0 : s_projected(:, :) = MATMUL(CONJG(TRANSPOSE(source_coeff)), s_coeff)
3250 0 : h_projected(:, :) = MATMUL(CONJG(TRANSPOSE(source_coeff)), h_coeff)
3251 0 : s_projected(:, :) = 0.5_dp*(s_projected + CONJG(TRANSPOSE(s_projected)))
3252 0 : h_projected(:, :) = 0.5_dp*(h_projected + CONJG(TRANSPOSE(h_projected)))
3253 :
3254 : reconstruct_window: BLOCK
3255 0 : CALL diag_complex(s_projected, metric_vectors, metric_values)
3256 0 : min_svalue = MINVAL(metric_values)
3257 0 : metric_deviation = MAXVAL(ABS(metric_values - 1.0_dp))
3258 0 : IF (min_svalue < 1.0e-10_dp) THEN
3259 : WRITE (reason, "(A,ES9.2,A,ES9.2)") &
3260 0 : "singular expanded metric smin=", min_svalue, " dS=", metric_deviation
3261 0 : EXIT reconstruct_window
3262 : END IF
3263 :
3264 0 : DO ib = 1, nmo_source
3265 0 : metric_vectors(:, ib) = metric_vectors(:, ib)/SQRT(metric_values(ib))
3266 : END DO
3267 0 : h_projected_work(:, :) = MATMUL(h_projected, metric_vectors)
3268 0 : h_projected(:, :) = MATMUL(CONJG(TRANSPOSE(metric_vectors)), h_projected_work)
3269 0 : h_projected(:, :) = 0.5_dp*(h_projected + CONJG(TRANSPOSE(h_projected)))
3270 0 : CALL diag_complex(h_projected, ritz_vectors, ritz_values)
3271 0 : h_projected_work(:, :) = MATMUL(metric_vectors, ritz_vectors)
3272 0 : ritz_vectors(:, :) = h_projected_work
3273 0 : stabilized(:, :) = MATMUL(source_coeff, ritz_vectors(:, 1:nmo_export))
3274 0 : coeff_work(:, :) = MATMUL(h_coeff, ritz_vectors)
3275 0 : h_coeff(:, :) = coeff_work
3276 0 : coeff_work(:, :) = MATMUL(s_coeff, ritz_vectors)
3277 0 : s_coeff(:, :) = coeff_work
3278 0 : DO ib = 1, nmo_export
3279 0 : norm_value = SQRT(ABS(REAL(DOT_PRODUCT(stabilized(:, ib), s_coeff(:, ib)), KIND=dp)))
3280 0 : IF (norm_value > EPSILON(1.0_dp)) THEN
3281 0 : stabilized(:, ib) = stabilized(:, ib)/norm_value
3282 0 : h_coeff(:, ib) = h_coeff(:, ib)/norm_value
3283 0 : s_coeff(:, ib) = s_coeff(:, ib)/norm_value
3284 : END IF
3285 : END DO
3286 0 : residual_block(:, :) = h_coeff(:, 1:nmo_export)
3287 0 : DO ib = 1, nmo_export
3288 0 : residual_block(:, ib) = residual_block(:, ib) - ritz_values(ib)*s_coeff(:, ib)
3289 : END DO
3290 0 : max_residual = MAXVAL(ABS(residual_block))
3291 0 : IF (max_residual > residual_tol) THEN
3292 : WRITE (reason, "(A,ES9.2)") &
3293 0 : "expanded dS=", metric_deviation
3294 0 : EXIT reconstruct_window
3295 : END IF
3296 0 : max_eigenvalue_shift = MAXVAL(ABS(ritz_values(1:nmo_export) - source_eigenvalues(1:nmo_export)))
3297 0 : IF (max_eigenvalue_shift > eigenvalue_tol) THEN
3298 : WRITE (reason, "(A,ES9.2)") &
3299 0 : "expanded dS=", metric_deviation
3300 0 : EXIT reconstruct_window
3301 : END IF
3302 :
3303 0 : CALL cp_cfm_set_submatrix(dst_cfm, stabilized)
3304 0 : success = .TRUE.
3305 :
3306 : END BLOCK reconstruct_window
3307 :
3308 0 : DEALLOCATE (source_coeff, s_coeff, h_coeff, h_projected, h_projected_work, metric_vectors, &
3309 0 : residual_block, ritz_vectors, s_projected, stabilized, coeff_work, metric_values, &
3310 0 : ritz_values, source_eigenvalues)
3311 0 : CALL cp_cfm_release(s_src)
3312 0 : CALL cp_cfm_release(h_src)
3313 :
3314 0 : END SUBROUTINE ritz_reconstruct_wannier90_window
3315 :
3316 : ! **************************************************************************************************
3317 : !> \brief Map the full Wannier90 mesh to SCF representative k-points and symmetry operations.
3318 : !> \param kpoint full Wannier90 export k-point object
3319 : !> \param qs_kpoint SCF k-point object
3320 : !> \param source_kpoint source representative index for each full k-point
3321 : !> \param sym_index symmetry entry in source kp_sym; 0 direct, -1 time reversal only
3322 : !> \param success true if every full k-point was mapped
3323 : !> \param reason diagnostic message
3324 : ! **************************************************************************************************
3325 60 : SUBROUTINE build_wannier90_scf_mapping(kpoint, qs_kpoint, source_kpoint, sym_index, success, &
3326 : reason)
3327 : TYPE(kpoint_type), POINTER :: kpoint, qs_kpoint
3328 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: source_kpoint, sym_index
3329 : LOGICAL, INTENT(OUT) :: success
3330 : CHARACTER(LEN=*), INTENT(OUT) :: reason
3331 :
3332 : INTEGER :: ik, ikred, imatch, isym, nfull
3333 : TYPE(kpoint_sym_type), POINTER :: kpsym
3334 :
3335 60 : success = .FALSE.
3336 60 : reason = ""
3337 60 : nfull = kpoint%nkp
3338 240 : ALLOCATE (source_kpoint(nfull), sym_index(nfull))
3339 60 : source_kpoint(:) = 0
3340 60 : sym_index(:) = 0
3341 :
3342 226 : DO ikred = 1, qs_kpoint%nkp
3343 166 : imatch = find_matching_kpoint(kpoint%xkp, qs_kpoint%xkp(1:3, ikred))
3344 226 : IF (imatch > 0 .AND. source_kpoint(imatch) == 0) THEN
3345 166 : source_kpoint(imatch) = ikred
3346 166 : sym_index(imatch) = 0
3347 : END IF
3348 : END DO
3349 :
3350 : ! Prefer pure time-reversal partners before general atom/AO symmetry operations.
3351 226 : DO ikred = 1, qs_kpoint%nkp
3352 664 : imatch = find_matching_kpoint(kpoint%xkp, -qs_kpoint%xkp(1:3, ikred))
3353 226 : IF (imatch > 0 .AND. source_kpoint(imatch) == 0) THEN
3354 70 : source_kpoint(imatch) = ikred
3355 70 : sym_index(imatch) = -1
3356 : END IF
3357 : END DO
3358 :
3359 60 : IF (ASSOCIATED(qs_kpoint%kp_sym)) THEN
3360 226 : DO ikred = 1, qs_kpoint%nkp
3361 166 : kpsym => qs_kpoint%kp_sym(ikred)%kpoint_sym
3362 166 : IF (.NOT. ASSOCIATED(kpsym)) CYCLE
3363 166 : IF (.NOT. kpsym%apply_symmetry) CYCLE
3364 8058 : DO isym = 1, kpsym%nwred
3365 7904 : imatch = find_matching_kpoint(kpoint%xkp, kpsym%xkp(1:3, isym))
3366 8070 : IF (imatch > 0 .AND. source_kpoint(imatch) == 0) THEN
3367 680 : source_kpoint(imatch) = ikred
3368 680 : sym_index(imatch) = isym
3369 : END IF
3370 : END DO
3371 : END DO
3372 : END IF
3373 :
3374 976 : DO ik = 1, nfull
3375 976 : IF (source_kpoint(ik) == 0) THEN
3376 0 : reason = "not all full-mesh k-points are represented by the SCF symmetry orbits"
3377 0 : RETURN
3378 : END IF
3379 : END DO
3380 60 : success = .TRUE.
3381 :
3382 60 : END SUBROUTINE build_wannier90_scf_mapping
3383 :
3384 : ! **************************************************************************************************
3385 : !> \brief Find a fractional k-point in a periodic mesh.
3386 : !> \param xkp_mesh mesh coordinates
3387 : !> \param xkp_search coordinate to find
3388 : !> \return matching index, or zero when no match is found
3389 : ! **************************************************************************************************
3390 8236 : INTEGER FUNCTION find_matching_kpoint(xkp_mesh, xkp_search) RESULT(ik_match)
3391 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp_mesh
3392 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp_search
3393 :
3394 : INTEGER :: ik
3395 :
3396 8236 : ik_match = 0
3397 125162 : DO ik = 1, SIZE(xkp_mesh, 2)
3398 125162 : IF (kpoint_same_periodic(xkp_mesh(1:3, ik), xkp_search)) THEN
3399 8236 : ik_match = ik
3400 8236 : RETURN
3401 : END IF
3402 : END DO
3403 :
3404 : END FUNCTION find_matching_kpoint
3405 :
3406 : ! **************************************************************************************************
3407 : !> \brief Infer a tensor-product Wannier90 mesh from explicit fractional k-point coordinates.
3408 : !> \param kpt_latt explicit k-point coordinates in reciprocal-lattice units
3409 : !> \param mp_grid inferred mesh dimensions
3410 : !> \param valid true if the coordinate set is compatible with a tensor-product mesh
3411 : ! **************************************************************************************************
3412 6 : SUBROUTINE infer_wannier_mp_grid(kpt_latt, mp_grid, valid)
3413 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: kpt_latt
3414 : INTEGER, DIMENSION(3), INTENT(OUT) :: mp_grid
3415 : LOGICAL, INTENT(OUT) :: valid
3416 :
3417 : INTEGER :: coord_id, i, idim, idx, n_unique, &
3418 : num_kpts, stride, unique_id
3419 : LOGICAL :: known
3420 6 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: seen
3421 : REAL(KIND=dp) :: coord
3422 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: unique_coord
3423 :
3424 6 : num_kpts = SIZE(kpt_latt, 2)
3425 6 : mp_grid(:) = 0
3426 18 : ALLOCATE (unique_coord(3, num_kpts))
3427 24 : DO idim = 1, 3
3428 : n_unique = 0
3429 162 : DO i = 1, num_kpts
3430 144 : coord = kpt_latt(idim, i) - FLOOR(kpt_latt(idim, i))
3431 144 : IF (ABS(coord - 1.0_dp) < 1.0e-8_dp) coord = 0.0_dp
3432 144 : known = .FALSE.
3433 216 : DO unique_id = 1, n_unique
3434 216 : IF (ABS(unique_coord(idim, unique_id) - coord) < 1.0e-8_dp) THEN
3435 : known = .TRUE.
3436 : EXIT
3437 : END IF
3438 : END DO
3439 162 : IF (.NOT. known) THEN
3440 36 : n_unique = n_unique + 1
3441 36 : unique_coord(idim, n_unique) = coord
3442 : END IF
3443 : END DO
3444 24 : mp_grid(idim) = n_unique
3445 : END DO
3446 6 : valid = (mp_grid(1)*mp_grid(2)*mp_grid(3) == num_kpts)
3447 6 : IF (valid) THEN
3448 18 : ALLOCATE (seen(num_kpts))
3449 6 : seen(:) = .FALSE.
3450 54 : DO i = 1, num_kpts
3451 : idx = 1
3452 : stride = 1
3453 192 : DO idim = 1, 3
3454 144 : coord = kpt_latt(idim, i) - FLOOR(kpt_latt(idim, i))
3455 144 : IF (ABS(coord - 1.0_dp) < 1.0e-8_dp) coord = 0.0_dp
3456 144 : coord_id = 0
3457 216 : DO unique_id = 1, mp_grid(idim)
3458 216 : IF (ABS(unique_coord(idim, unique_id) - coord) < 1.0e-8_dp) THEN
3459 : coord_id = unique_id
3460 : EXIT
3461 : END IF
3462 : END DO
3463 144 : CPASSERT(coord_id > 0)
3464 144 : idx = idx + (coord_id - 1)*stride
3465 192 : stride = stride*mp_grid(idim)
3466 : END DO
3467 48 : IF (seen(idx)) valid = .FALSE.
3468 54 : seen(idx) = .TRUE.
3469 : END DO
3470 54 : valid = valid .AND. ALL(seen)
3471 6 : DEALLOCATE (seen)
3472 : END IF
3473 6 : DEALLOCATE (unique_coord)
3474 :
3475 6 : END SUBROUTINE infer_wannier_mp_grid
3476 :
3477 272 : END MODULE qs_wannier90
|
