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