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