(git:c641773)
Loading...
Searching...
No Matches
qs_scf_diagonalization.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 Different diagonalization schemes that can be used
10!> for the iterative solution of the eigenvalue problem
11!> \par History
12!> started from routines previously located in the qs_scf module
13!> 05.2009
14! **************************************************************************************************
15MODULE qs_scf_diagonalization
16 USE cp_array_utils, ONLY: cp_1d_r_p_type
17 USE cp_blacs_env, ONLY: cp_blacs_env_type
18 USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,&
19 cp_cfm_scale_and_add,&
20 cp_cfm_scale_and_add_fm,&
21 cp_cfm_transpose,&
22 cp_cfm_triangular_multiply
23 USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose
24 USE cp_cfm_diag, ONLY: cp_cfm_diag_workspace_type,&
25 cp_cfm_geeig,&
26 cp_cfm_geeig_canon,&
27 cp_cfm_geeig_canon_local,&
28 cp_cfm_geeig_local,&
29 cp_cfm_heevd,&
30 cp_cfm_local_workspace_prepare
31 USE cp_cfm_types, ONLY: cp_cfm_create,&
32 cp_cfm_release,&
33 cp_cfm_set_all,&
34 cp_cfm_to_cfm,&
35 cp_cfm_to_fm,&
36 cp_cfm_type,&
37 cp_fm_to_cfm
38 USE cp_control_types, ONLY: dft_control_type,&
39 hairy_probes_type
40 USE cp_dbcsr_api, ONLY: &
41 dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_init_p, dbcsr_p_type, &
42 dbcsr_release_p, dbcsr_set, dbcsr_type, dbcsr_type_symmetric
43 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
44 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
45 copy_fm_to_dbcsr,&
46 cp_dbcsr_sm_fm_multiply,&
47 dbcsr_allocate_matrix_set
48 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
49 cp_fm_scale_and_add,&
50 cp_fm_symm,&
51 cp_fm_uplo_to_full
52 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_reduce,&
53 cp_fm_cholesky_restore
54 USE cp_fm_diag, ONLY: choose_eigv_solver,&
55 cp_fm_geeig,&
56 cp_fm_geeig_canon,&
57 diag_check_requested,&
58 diag_lib_explicit,&
59 direct_generalized_diagonalization
60 USE cp_fm_pool_types, ONLY: cp_fm_pool_type,&
61 fm_pool_create_fm,&
62 fm_pool_get_el_struct,&
63 fm_pool_give_back_fm
64 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
65 cp_fm_struct_release,&
66 cp_fm_struct_type
67 USE cp_fm_types, ONLY: cp_fm_add_to_element,&
68 cp_fm_create,&
69 cp_fm_get_info,&
70 cp_fm_release,&
71 cp_fm_to_fm,&
72 cp_fm_type
73 USE cp_log_handling, ONLY: cp_get_default_logger,&
74 cp_logger_get_default_unit_nr,&
75 cp_logger_type
76 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
77 cp_print_key_unit_nr,&
78 debug_print_level
79 USE dft_plus_u, ONLY: lowdin_kp_apply
80 USE ieee_arithmetic, ONLY: ieee_is_finite
81 USE input_constants, ONLY: &
82 cholesky_dbcsr, cholesky_inverse, cholesky_off, cholesky_reduce, cholesky_restore, &
83 core_guess, diag_update_method_adiis, general_roks, high_spin_roks, ot_precond_full_all, &
84 ot_precond_full_kinetic, ot_precond_full_single_inverse, ot_precond_none, &
85 ot_precond_s_inverse, restart_guess
86 USE input_cp2k_kpoints, ONLY: lattice_fft_auto,&
87 lattice_fft_off,&
88 lattice_fft_on
89 USE input_section_types, ONLY: section_vals_get_subs_vals,&
90 section_vals_type
91 USE kinds, ONLY: dp,&
92 int_8
93 USE kpoint_methods, ONLY: kpoint_density_matrices,&
94 kpoint_density_transform,&
95 kpoint_set_mo_occupation,&
96 rskp_grid_type,&
97 rskp_transform_grid_prepare,&
98 rskp_transform_grid_release
99 USE kpoint_types, ONLY: get_kpoint_info,&
100 kpoint_env_type,&
101 kpoint_type
102 USE local_gemm_api, ONLY: local_gemm_pu_gpu,&
103 local_gemm_ctxt_type
104 USE machine, ONLY: m_flush,&
105 m_walltime
106 USE mathconstants, ONLY: z_one,&
107 z_zero
108 USE message_passing, ONLY: mp_para_env_type
109 USE parallel_gemm_api, ONLY: parallel_gemm
110 USE preconditioner, ONLY: dbcsr_pair_to_cfm,&
111 make_preconditioner_complex_full_kinetic,&
112 make_preconditioner_complex_full_s_inverse,&
113 prepare_preconditioner,&
114 restart_preconditioner
115 USE preconditioner_makes, ONLY: make_complex_full_all,&
116 make_complex_full_single_inverse
117 USE preconditioner_types, ONLY: destroy_preconditioner,&
118 init_preconditioner
119 USE qs_density_matrices, ONLY: calculate_density_matrix
120 USE qs_density_mixing_types, ONLY: direct_mixing_nr,&
121 gspace_mixing_nr
122 USE qs_diis, ONLY: qs_diis_b_calc_err_kp,&
123 qs_diis_b_check_i_alloc_kp,&
124 qs_diis_b_clear,&
125 qs_diis_b_clear_kp,&
126 qs_diis_b_info_kp,&
127 qs_diis_b_step,&
128 qs_diis_b_step_kp
129 USE qs_energy_types, ONLY: qs_energy_type
130 USE qs_environment_types, ONLY: get_qs_env,&
131 qs_environment_type
132 USE qs_gspace_mixing, ONLY: gspace_mixing
133 USE qs_kpoint_operators, ONLY: &
134 kpoint_operator_batch_size, kpoint_operator_context_create, &
135 kpoint_operator_context_release, kpoint_operator_context_type, kpoint_operator_evaluate, &
136 kpoint_operator_evaluate_batch, kpoint_operator_finish, kpoint_operator_get, &
137 kpoint_operator_get_local, kpoint_operator_prepare, kpoint_operator_start, &
138 kpoint_operator_store, kpoint_operator_work_size, kpoint_slot_ks, kpoint_slot_s, &
139 kpoint_slot_t, kpoint_spin_free
140 USE qs_ks_methods, ONLY: qs_ks_update_qs_env
141 USE qs_ks_types, ONLY: qs_ks_did_change,&
142 qs_ks_env_type
143 USE qs_matrix_pools, ONLY: mpools_get,&
144 qs_matrix_pools_type
145 USE qs_mixing_utils, ONLY: charge_mixing_init,&
146 mixing_allocate,&
147 mixing_init,&
148 self_consistency_check
149 USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
150 USE qs_mo_occupation, ONLY: set_mo_occupation
151 USE qs_mo_types, ONLY: get_mo_set,&
152 mo_set_type
153 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
154 USE qs_ot_eigensolver, ONLY: ot_eigensolver,&
155 ot_eigensolver_complex
156 USE qs_rho_atom_types, ONLY: rho_atom_type
157 USE qs_rho_methods, ONLY: qs_rho_update_rho
158 USE qs_rho_types, ONLY: qs_rho_get,&
159 qs_rho_type
160 USE qs_scf_block_davidson, ONLY: generate_extended_space,&
161 generate_extended_space_c,&
162 generate_extended_space_sparse
163 USE qs_scf_lanczos, ONLY: lanczos_refinement,&
164 lanczos_refinement_2v
165 USE qs_scf_lanczos_c, ONLY: lanczos_channel_begin,&
166 lanczos_channel_end,&
167 lanczos_sweep_suffix,&
168 lanczos_workspace_create,&
169 lanczos_workspace_release,&
170 lanczos_workspace_type
171 USE qs_scf_methods, ONLY: combine_ks_matrices,&
172 eigensolver,&
173 eigensolver_dbcsr,&
174 eigensolver_generalized,&
175 eigensolver_simple,&
176 eigensolver_symm,&
177 scf_env_density_mixing
178 USE qs_scf_types, ONLY: qs_scf_env_type,&
179 subspace_env_type
180 USE scf_control_types, ONLY: scf_control_type
181
182!$ USE OMP_LIB, ONLY: omp_get_max_threads, &
183!$ omp_get_thread_num
184#include "./base/base_uses.f90"
185
186 IMPLICIT NONE
187
188 PRIVATE
189
190 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_diagonalization'
191
192 ! occupation below this weight sits beyond the thermal smearing tail
193 ! and converges at the relaxed empty threshold. The Davidson k-point
194 ! driver and both Lanczos k-point paths below use this rule
195 REAL(KIND=dp), PARAMETER, PRIVATE :: occ_tol = 1.0e-3_dp
196
197 TYPE, PRIVATE :: kp_diag_workspace_type
198 TYPE(cp_cfm_type) :: cksmat, cmos, csmat, csmat_base, cwork
199 TYPE(cp_cfm_diag_workspace_type) :: solver
200 TYPE(lanczos_workspace_type) :: lan_ws
201 TYPE(local_gemm_ctxt_type) :: gemm_ctx
202 COMPLEX(KIND=dp), ALLOCATABLE :: phase(:, :)
203 END TYPE kp_diag_workspace_type
204
205 PUBLIC :: do_general_diag, do_general_diag_kp, do_roks_diag, &
206 do_special_diag, do_ot_diag, do_ot_diag_kp, do_block_davidson_diag, &
207 do_block_davidson_diag_kp, &
208 do_block_krylov_diag, do_block_krylov_diag_kp, do_scf_diag_subspace, diag_subspace_allocate, &
209 general_eigenproblem, diag_kp_smat, diag_kp_basic
210
211CONTAINS
212
213! **************************************************************************************************
214!> \brief Fade candidate regularization through handover; never regularize a raw-map check.
215!> \param scf_env Current SCF state.
216!> \return ...
217! **************************************************************************************************
218 FUNCTION adiis_candidate_shift(scf_env) RESULT(shift)
219 TYPE(qs_scf_env_type), INTENT(IN) :: scf_env
220 REAL(kind=dp) :: shift
221
222 shift = 0.0_dp
223 IF (scf_env%iter_count <= 1 .OR. scf_env%adiis_check_next) RETURN
224 IF (.NOT. ASSOCIATED(scf_env%scf_subspace_buffer)) RETURN
225 shift = scf_env%adiis_shift*(1.0_dp - scf_env%scf_subspace_buffer%diis_weight)
226 END FUNCTION adiis_candidate_shift
227
228! **************************************************************************************************
229!> \brief Add -shift*S*C*(occupation/maxocc)*C^T*S to a candidate Fock matrix.
230!> The occupied projector does not require any stored virtual orbitals.
231!> \param fock Candidate matrix; physical F[P] and all history matrices remain unchanged.
232!> \param overlap AO overlap.
233!> \param mos Previous physical orbitals and occupations.
234!> \param shift Nonnegative regularization parameter in hartree.
235! **************************************************************************************************
236 SUBROUTINE shift_candidate_real(fock, overlap, mos, shift)
237 TYPE(cp_fm_type), INTENT(INOUT) :: fock
238 TYPE(cp_fm_type), INTENT(IN) :: overlap
239 TYPE(mo_set_type), INTENT(IN) :: mos
240 REAL(kind=dp), INTENT(IN) :: shift
241
242 INTEGER :: nao, nmo
243 REAL(kind=dp) :: maxocc
244 REAL(kind=dp), DIMENSION(:), POINTER :: occupation
245 TYPE(cp_fm_type) :: sc
246 TYPE(cp_fm_type), POINTER :: coeff
247
248 IF (shift == 0.0_dp) RETURN
249 CALL get_mo_set(mos, nao=nao, nmo=nmo, maxocc=maxocc, occupation_numbers=occupation, mo_coeff=coeff)
250 CALL cp_fm_create(sc, coeff%matrix_struct)
251 CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, overlap, coeff, 0.0_dp, sc)
252 CALL cp_fm_column_scale(sc, sqrt(max(0.0_dp, occupation/maxocc)))
253 CALL parallel_gemm("N", "T", nao, nao, nmo, -shift, sc, sc, 1.0_dp, fock)
254 CALL cp_fm_release(sc)
255 END SUBROUTINE shift_candidate_real
256
257! **************************************************************************************************
258!> \brief Complex occupied-space shift, using conjugate transpose and existing eigensolver workspaces.
259!> \param fock Candidate Fock matrix, modified in place.
260!> \param overlap AO overlap, unchanged.
261!> \param mos Previous complex MOs.
262!> \param coeff AO/MO eigensolver workspace.
263!> \param work AO/AO eigensolver workspace.
264!> \param shift Nonnegative regularization parameter in hartree.
265!> \param local Whether the matrices are complete on this rank.
266!> \param gemm_ctx Private local GEMM context, unused for distributed groups.
267! **************************************************************************************************
268 SUBROUTINE shift_candidate_complex(fock, overlap, mos, coeff, work, shift, local, gemm_ctx)
269 TYPE(cp_cfm_type), INTENT(INOUT) :: fock
270 TYPE(cp_cfm_type), INTENT(IN) :: overlap
271 TYPE(mo_set_type), INTENT(IN) :: mos
272 TYPE(cp_cfm_type), INTENT(INOUT) :: coeff, work
273 REAL(kind=dp), INTENT(IN) :: shift
274 LOGICAL, INTENT(IN) :: local
275 TYPE(local_gemm_ctxt_type), INTENT(INOUT) :: gemm_ctx
276
277 INTEGER :: nao, nmo
278 REAL(kind=dp) :: maxocc
279 REAL(kind=dp), DIMENSION(:), POINTER :: occupation
280 TYPE(cp_cfm_type), POINTER :: cmos
281
282 IF (shift == 0.0_dp) RETURN
283 CALL get_mo_set(mos, nao=nao, nmo=nmo, maxocc=maxocc, occupation_numbers=occupation, cmo_coeff=cmos)
284 CALL cp_cfm_to_cfm(cmos, coeff)
285 CALL cp_cfm_column_scale(coeff, cmplx(sqrt(max(0.0_dp, occupation/maxocc)), kind=dp))
286 IF (local) THEN
287 cpassert(product(overlap%matrix_struct%context%num_pe) == 1)
288 cpassert(product(coeff%matrix_struct%context%num_pe) == 1)
289 cpassert(product(work%matrix_struct%context%num_pe) == 1)
290 cpassert(product(fock%matrix_struct%context%num_pe) == 1)
291 CALL gemm_ctx%gemm("N", "N", nao, nmo, nao, z_one, &
292 overlap%local_data, SIZE(overlap%local_data, 1), &
293 coeff%local_data, SIZE(coeff%local_data, 1), z_zero, &
294 work%local_data, SIZE(work%local_data, 1))
295 CALL gemm_ctx%gemm("N", "C", nao, nao, nmo, cmplx(-shift, kind=dp), &
296 work%local_data, SIZE(work%local_data, 1), &
297 work%local_data, SIZE(work%local_data, 1), z_one, &
298 fock%local_data, SIZE(fock%local_data, 1))
299 ELSE
300 CALL parallel_gemm("N", "N", nao, nmo, nao, z_one, overlap, coeff, z_zero, work)
301 CALL parallel_gemm("N", "C", nao, nao, nmo, cmplx(-shift, kind=dp), work, work, z_one, fock)
302 END IF
303 END SUBROUTINE shift_candidate_complex
304
305! **************************************************************************************************
306!> \brief the inner loop of scf, specific to diagonalization with S matrix
307!> basically, in goes the ks matrix out goes a new p matrix
308!> \param scf_env ...
309!> \param mos ...
310!> \param matrix_ks ...
311!> \param matrix_s ...
312!> \param scf_control ...
313!> \param scf_section ...
314!> \param diis_step ...
315!> \par History
316!> 03.2006 created [Joost VandeVondele]
317! **************************************************************************************************
318
319 SUBROUTINE general_eigenproblem(scf_env, mos, matrix_ks, &
320 matrix_s, scf_control, scf_section, &
321 diis_step)
322
323 TYPE(qs_scf_env_type), POINTER :: scf_env
324 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
325 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
326 TYPE(scf_control_type), POINTER :: scf_control
327 TYPE(section_vals_type), POINTER :: scf_section
328 LOGICAL, INTENT(INOUT) :: diis_step
329
330 INTEGER :: cdiis_cholesky_method, cdiis_ncall, &
331 ispin, nspin
332 LOGICAL :: cdiis_step, do_level_shift, owns_ortho, &
333 use_adiis, use_cdiis, use_jacobi, &
334 use_subspace_fock
335 REAL(kind=dp) :: candidate_shift, cdiis_error, &
336 diis_error, diis_weight, eps_diis
337 TYPE(cp_fm_type), POINTER :: ortho
338 TYPE(dbcsr_type), POINTER :: effective_ks_dbcsr, ortho_dbcsr
339
340 nspin = SIZE(matrix_ks)
341 NULLIFY (effective_ks_dbcsr, ortho, ortho_dbcsr)
342
343 use_adiis = scf_control%diagonalization%update_method == diag_update_method_adiis
344 use_subspace_fock = .false.
345 IF (use_adiis) THEN
346 IF (.NOT. ASSOCIATED(scf_env%scf_subspace_buffer)) THEN
347 cpabort("ADIIS is only available in the main Quickstep SCF workflow")
348 END IF
349 use_subspace_fock = scf_env%scf_subspace_buffer%use_combined_fock
350 END IF
351
352 eps_diis = scf_control%eps_diis
353 use_cdiis = .false.
354 IF (use_adiis .AND. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
355 IF (.NOT. scf_env%skip_diis .AND. scf_env%iter_count > 1 .AND. &
356 .NOT. scf_env%adiis_check_next .AND. &
357 .NOT. scf_env%scf_subspace_buffer%last_restart .AND. &
358 scf_env%scf_subspace_buffer%nstored > 0) THEN
359 use_cdiis = scf_env%iter_delta < eps_diis .OR. scf_env%scf_diis_buffer%ncall > 0
360 END IF
361 END IF
362
363 DO ispin = 1, nspin
364 IF (use_subspace_fock .AND. .NOT. use_cdiis) THEN
365 CALL copy_dbcsr_to_fm(scf_env%scf_subspace_buffer%combined_fock(ispin, 1)%matrix, &
366 scf_env%scf_work1(ispin), scf_env%fm_copy_plan(ispin), check_plan=.true.)
367 ELSE
368 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, scf_env%scf_work1(ispin), &
369 scf_env%fm_copy_plan(ispin), check_plan=.true.)
370 END IF
371 END DO
372
373 IF (use_adiis) THEN
374 diis_step = .false.
375 diis_error = scf_env%scf_subspace_buffer%last_old_fock_weight
376 diis_weight = scf_env%scf_subspace_buffer%diis_weight
377 IF (use_cdiis) THEN
378 cdiis_ncall = scf_env%scf_diis_buffer%ncall
379 CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
380 scf_env%scf_work2, scf_env%iter_delta, cdiis_error, cdiis_step, &
381 eps_diis, 2, s_matrix=matrix_s, scf_section=scf_section)
382 IF (.NOT. ieee_is_finite(cdiis_error) .OR. cdiis_error >= eps_diis) THEN
383 cdiis_step = .false.
384 IF (cdiis_ncall < scf_env%scf_diis_buffer%nbuffer) THEN
385 ! The rejected state has not overwritten a full circular buffer.
386 scf_env%scf_diis_buffer%ncall = cdiis_ncall
387 ELSE
388 CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
389 diis_weight = 0.0_dp
390 END IF
391 END IF
392 IF (cdiis_step) THEN
393 ! Use one half-weight handover step; MAX_DIIS only controls the buffer capacity.
394 diis_weight = min(1.0_dp, diis_weight + 0.5_dp)
395 ELSE
396 ! A rejected or deferred CDIIS step must not leave a partial raw-Fock blend active.
397 diis_weight = 0.0_dp
398 END IF
399 IF (diis_weight < 1.0_dp) THEN
400 DO ispin = 1, nspin
401 IF (use_subspace_fock) THEN
402 CALL copy_dbcsr_to_fm( &
403 scf_env%scf_subspace_buffer%combined_fock(ispin, 1)%matrix, scf_env%scf_work2, &
404 scf_env%fm_copy_plan(ispin), check_plan=.true.)
405 ELSE
406 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, scf_env%scf_work2, &
407 scf_env%fm_copy_plan(ispin), check_plan=.true.)
408 END IF
409 CALL cp_fm_scale_and_add(diis_weight, scf_env%scf_work1(ispin), &
410 1.0_dp - diis_weight, scf_env%scf_work2)
411 END DO
412 END IF
413 diis_step = cdiis_step .AND. diis_weight >= 1.0_dp
414 IF (diis_step) diis_error = cdiis_error
415 ELSE IF (ASSOCIATED(scf_env%scf_diis_buffer)) THEN
416 IF (scf_env%scf_diis_buffer%ncall > 0) CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
417 diis_weight = 0.0_dp
418 END IF
419 scf_env%scf_subspace_buffer%diis_weight = diis_weight
420 ELSE
421 IF (scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis) THEN
422 CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
423 scf_env%scf_work2, scf_env%iter_delta, diis_error, diis_step, &
424 eps_diis, scf_control%nmixing, &
425 s_matrix=matrix_s, &
426 scf_section=scf_section)
427 ELSE
428 diis_step = .false.
429 END IF
430 END IF
431
432 candidate_shift = 0.0_dp
433 IF (use_adiis) candidate_shift = adiis_candidate_shift(scf_env)
434 IF (candidate_shift > 0.0_dp) THEN
435 CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, scf_env%scf_work2, scf_env%fm_copy_plan(0), check_plan=.true.)
436 DO ispin = 1, nspin
437 CALL shift_candidate_real(scf_env%scf_work1(ispin), scf_env%scf_work2, mos(ispin), candidate_shift)
438 END DO
439 END IF
440
441 do_level_shift = ((scf_control%level_shift /= 0.0_dp) .AND. &
442 ((scf_control%density_guess == core_guess) .OR. &
443 (scf_env%iter_count > 1)))
444
445 IF (.NOT. use_adiis .AND. (scf_env%iter_count > 1) .AND. &
446 (scf_env%iter_delta < scf_control%diagonalization%eps_jacobi)) THEN
447 use_jacobi = .true.
448 ELSE
449 use_jacobi = .false.
450 END IF
451
452 IF (use_adiis .AND. .NOT. diis_step) THEN
453 scf_env%iter_param = diis_error
454 IF (scf_env%scf_subspace_buffer%last_restart) THEN
455 scf_env%iter_method = "ADIIS/Rst."
456 ELSE
457 scf_env%iter_method = "ADIIS/Diag."
458 END IF
459 ELSE IF (diis_step) THEN
460 scf_env%iter_param = diis_error
461 IF (use_jacobi) THEN
462 scf_env%iter_method = "DIIS/Jacobi"
463 ELSE
464 scf_env%iter_method = "DIIS/Diag."
465 END IF
466 ELSE
467 IF (scf_env%mixing_method == 0) THEN
468 scf_env%iter_method = "NoMix/Diag."
469 ELSE IF (scf_env%mixing_method == 1) THEN
470 scf_env%iter_param = scf_env%p_mix_alpha
471 IF (use_jacobi) THEN
472 scf_env%iter_method = "P_Mix/Jacobi"
473 ELSE
474 scf_env%iter_method = "P_Mix/Diag."
475 END IF
476 ELSE IF (scf_env%mixing_method > 1) THEN
477 scf_env%iter_param = scf_env%mixing_store%alpha
478 IF (use_jacobi) THEN
479 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Jacobi"
480 ELSE
481 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Diag."
482 END IF
483 END IF
484 END IF
485
486 IF (scf_env%cholesky_method == cholesky_dbcsr .AND. .NOT. use_cdiis) THEN
487 ortho_dbcsr => scf_env%ortho_dbcsr
488 DO ispin = 1, nspin
489 effective_ks_dbcsr => matrix_ks(ispin)%matrix
490 IF (use_subspace_fock) THEN
491 effective_ks_dbcsr => scf_env%scf_subspace_buffer%combined_fock(ispin, 1)%matrix
492 END IF
493 CALL eigensolver_dbcsr(matrix_ks=effective_ks_dbcsr, matrix_ks_fm=scf_env%scf_work1(ispin), &
494 mo_set=mos(ispin), &
495 ortho_dbcsr=ortho_dbcsr, &
496 ksbuf1=scf_env%buf1_dbcsr, ksbuf2=scf_env%buf2_dbcsr)
497 END DO
498
499 ELSE IF (scf_env%cholesky_method > cholesky_off) THEN
500 cdiis_cholesky_method = scf_env%cholesky_method
501 IF (cdiis_cholesky_method == cholesky_dbcsr) THEN
502 cdiis_cholesky_method = cholesky_inverse
503 END IF
504 IF (scf_env%cholesky_method == cholesky_inverse) THEN
505 ortho => scf_env%ortho_m1
506 ELSE
507 ortho => scf_env%ortho
508 END IF
509
510 owns_ortho = .false.
511 IF (.NOT. ASSOCIATED(ortho)) THEN
512 ALLOCATE (ortho)
513 owns_ortho = .true.
514 END IF
515
516 DO ispin = 1, nspin
517 IF (direct_generalized_diagonalization .AND. .NOT. do_level_shift) THEN
518 CALL eigensolver_generalized(matrix_ks_fm=scf_env%scf_work1(ispin), &
519 matrix_s=matrix_s(ispin)%matrix, &
520 mo_set=mos(ispin), &
521 work=scf_env%scf_work2)
522 ELSE
523 IF (do_level_shift) THEN
524 CALL eigensolver(matrix_ks_fm=scf_env%scf_work1(ispin), &
525 mo_set=mos(ispin), &
526 ortho=ortho, &
527 work=scf_env%scf_work2, &
528 cholesky_method=cdiis_cholesky_method, &
529 do_level_shift=do_level_shift, &
530 level_shift=scf_control%level_shift, &
531 matrix_u_fm=scf_env%ortho, &
532 use_jacobi=use_jacobi)
533 ELSE
534 CALL eigensolver(matrix_ks_fm=scf_env%scf_work1(ispin), &
535 mo_set=mos(ispin), &
536 ortho=ortho, &
537 work=scf_env%scf_work2, &
538 cholesky_method=cdiis_cholesky_method, &
539 do_level_shift=do_level_shift, &
540 level_shift=scf_control%level_shift, &
541 use_jacobi=use_jacobi)
542 END IF
543 END IF
544 END DO
545
546 IF (scf_env%cholesky_method /= cholesky_dbcsr) THEN
547 scf_env%cholesky_method = cdiis_cholesky_method
548 END IF
549
550 IF (owns_ortho) DEALLOCATE (ortho)
551 ELSE
552 ortho => scf_env%ortho
553
554 owns_ortho = .false.
555 IF (.NOT. ASSOCIATED(ortho)) THEN
556 ALLOCATE (ortho)
557 owns_ortho = .true.
558 END IF
559
560 IF (do_level_shift) THEN
561 DO ispin = 1, nspin
562 IF (ASSOCIATED(scf_env%scf_work1_red) .AND. ASSOCIATED(scf_env%scf_work2_red) &
563 .AND. ASSOCIATED(scf_env%ortho_red) .AND. ASSOCIATED(scf_env%ortho_m1_red)) THEN
564 CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
565 mo_set=mos(ispin), &
566 ortho=ortho, &
567 work=scf_env%scf_work2, &
568 do_level_shift=do_level_shift, &
569 level_shift=scf_control%level_shift, &
570 matrix_u_fm=scf_env%ortho_m1, &
571 use_jacobi=use_jacobi, &
572 jacobi_threshold=scf_control%diagonalization%jacobi_threshold, &
573 matrix_ks_fm_red=scf_env%scf_work1_red(ispin), &
574 ortho_red=scf_env%ortho_red, &
575 work_red=scf_env%scf_work2_red, &
576 matrix_u_fm_red=scf_env%ortho_m1_red)
577 ELSE
578 CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
579 mo_set=mos(ispin), &
580 ortho=ortho, &
581 work=scf_env%scf_work2, &
582 do_level_shift=do_level_shift, &
583 level_shift=scf_control%level_shift, &
584 matrix_u_fm=scf_env%ortho_m1, &
585 use_jacobi=use_jacobi, &
586 jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
587 END IF
588 END DO
589 ELSE
590 DO ispin = 1, nspin
591 IF (ASSOCIATED(scf_env%scf_work1_red) .AND. ASSOCIATED(scf_env%scf_work2_red) &
592 .AND. ASSOCIATED(scf_env%ortho_red)) THEN
593 CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
594 mo_set=mos(ispin), &
595 ortho=ortho, &
596 work=scf_env%scf_work2, &
597 do_level_shift=do_level_shift, &
598 level_shift=scf_control%level_shift, &
599 use_jacobi=use_jacobi, &
600 jacobi_threshold=scf_control%diagonalization%jacobi_threshold, &
601 matrix_ks_fm_red=scf_env%scf_work1_red(ispin), &
602 ortho_red=scf_env%ortho_red, &
603 work_red=scf_env%scf_work2_red)
604 ELSE
605 CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
606 mo_set=mos(ispin), &
607 ortho=ortho, &
608 work=scf_env%scf_work2, &
609 do_level_shift=do_level_shift, &
610 level_shift=scf_control%level_shift, &
611 use_jacobi=use_jacobi, &
612 jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
613 END IF
614 END DO
615 END IF
616
617 IF (owns_ortho) DEALLOCATE (ortho)
618 END IF
619
620 END SUBROUTINE general_eigenproblem
621
622! **************************************************************************************************
623!> \brief ...
624!> \param scf_env ...
625!> \param mos ...
626!> \param matrix_ks ...
627!> \param matrix_s ...
628!> \param scf_control ...
629!> \param scf_section ...
630!> \param diis_step ...
631!> \param probe ...
632! **************************************************************************************************
633 SUBROUTINE do_general_diag(scf_env, mos, matrix_ks, &
634 matrix_s, scf_control, scf_section, &
635 diis_step, probe)
636
637 TYPE(qs_scf_env_type), POINTER :: scf_env
638 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
639 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
640 TYPE(scf_control_type), POINTER :: scf_control
641 TYPE(section_vals_type), POINTER :: scf_section
642 LOGICAL, INTENT(INOUT) :: diis_step
643 TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
644 POINTER :: probe
645
646 INTEGER :: ispin, nspin
647 REAL(kind=dp) :: total_zeff_corr
648
649 nspin = SIZE(matrix_ks)
650
651 CALL general_eigenproblem(scf_env, mos, matrix_ks, &
652 matrix_s, scf_control, scf_section, diis_step)
653
654 total_zeff_corr = 0.0_dp
655 total_zeff_corr = scf_env%sum_zeff_corr
656
657 IF (abs(total_zeff_corr) > 0.0_dp) THEN
658 CALL set_mo_occupation(mo_array=mos, &
659 smear=scf_control%smear, tot_zeff_corr=total_zeff_corr)
660 ELSE
661 IF (PRESENT(probe) .EQV. .true.) THEN
662 scf_control%smear%do_smear = .false.
663 CALL set_mo_occupation(mo_array=mos, &
664 smear=scf_control%smear, &
665 probe=probe)
666 ELSE
667 IF (.NOT. scf_control%gce%do_gce) THEN
668 CALL set_mo_occupation(mo_array=mos, &
669 smear=scf_control%smear)
670 ELSE
671 CALL set_mo_occupation(mo_array=mos, &
672 smear=scf_control%smear, &
673 gce=scf_control%gce)
674 END IF
675 END IF
676 END IF
677
678 DO ispin = 1, nspin
679 CALL calculate_density_matrix(mos(ispin), &
680 scf_env%p_mix_new(ispin, 1)%matrix)
681 END DO
682
683 END SUBROUTINE do_general_diag
684
685! **************************************************************************************************
686!> \brief Kpoint diagonalization routine
687!> Transforms matrices to kpoint, distributes kpoint groups, performs
688!> general diagonalization (no storgae of overlap decomposition), stores
689!> MOs, calculates occupation numbers, calculates density matrices
690!> in kpoint representation, transforms density matrices to real space
691!> \param matrix_ks Kohn-sham matrices (RS indices, global)
692!> \param matrix_s Overlap matrices (RS indices, global)
693!> \param kpoints Kpoint environment
694!> \param scf_env SCF environment
695!> \param scf_control SCF control variables
696!> \param update_p ...
697!> \param diis_step ...
698!> \param diis_error ...
699!> \param qs_env ...
700!> \param probe ...
701!> \param added_mos_auto_grow true when the adaptive virtual-space buffer must grow
702!> \param potential Fixed Lowdin orbital potential for non-self-consistent k points.
703!> \par History
704!> 08.2014 created [JGH]
705! **************************************************************************************************
706 SUBROUTINE do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, update_p, &
707 diis_step, diis_error, qs_env, probe, added_mos_auto_grow, potential)
708
709 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
710 TYPE(kpoint_type), POINTER :: kpoints
711 TYPE(qs_scf_env_type), POINTER :: scf_env
712 TYPE(scf_control_type), POINTER :: scf_control
713 LOGICAL, INTENT(IN) :: update_p
714 LOGICAL, INTENT(INOUT) :: diis_step
715 REAL(dp), INTENT(INOUT), OPTIONAL :: diis_error
716 TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
717 TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
718 POINTER :: probe
719 LOGICAL, INTENT(OUT), OPTIONAL :: added_mos_auto_grow
720 REAL(dp), INTENT(IN), OPTIONAL, TARGET :: potential(:, :)
721
722 CHARACTER(len=*), PARAMETER :: routinen = 'do_general_diag_kp'
723
724 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coeffs
725 REAL(dp), POINTER :: effective_u(:, :), adiis_u(:, :)
726 INTEGER :: cdiis_ncall, handle, ib, ikp, ispin, kplocal, kpmax, &
727 lattice_fft, nb, nkp, nspin, omp_threads, nworkers, fft_ready
728 INTEGER, DIMENSION(2) :: kp_range
729 INTEGER, DIMENSION(3) :: nkp_grid
730 INTEGER, DIMENSION(:, :), POINTER :: kp_dist
731 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
732 LOGICAL :: adiis_cdiis, cache_smat, cdiis_step, do_diis, my_added_mos_auto_grow, prepared_ops, prepared_part, &
733 reuse_smat, spin_grid_fft, use_adiis, use_grid_fft, use_local_kpoints, use_real_wfn
734 REAL(kind=dp) :: diis_weight
735 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
736 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
737 TYPE(cp_cfm_type), POINTER :: cksmat, cmos, csmat, cwork
738 TYPE(kp_diag_workspace_type), ALLOCATABLE, TARGET :: workspace(:)
739 TYPE(cp_fm_pool_type), POINTER :: ao_ao_fm_pool
740 TYPE(cp_fm_struct_type), POINTER :: matrix_struct, mo_struct
741 TYPE(cp_fm_type) :: fmlocal, rksmat, rsmat
742 TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
743 TYPE(cp_fm_type), POINTER :: mo_coeff
744 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: adiis_ks, effective_ks
745 TYPE(dbcsr_type), POINTER :: rmatrix
746 TYPE(kpoint_env_type), POINTER :: kp
747 TYPE(kpoint_operator_context_type) :: op_ctx
748 TYPE(mp_para_env_type), POINTER :: para_env
749 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
750 POINTER :: sab_nl
751 TYPE(qs_matrix_pools_type), POINTER :: mpools
752 TYPE(rskp_grid_type) :: s_grid
753 TYPE(rskp_grid_type), ALLOCATABLE, DIMENSION(:) :: ks_grid
754 TYPE(section_vals_type), POINTER :: scf_section
755
756 CALL timeset(routinen, handle)
757
758 my_added_mos_auto_grow = .false.
759 IF (PRESENT(added_mos_auto_grow)) added_mos_auto_grow = .false.
760
761 NULLIFY (adiis_ks, effective_ks, sab_nl)
762 CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
763 kp_dist=kp_dist, sab_nl=sab_nl, &
764 cell_to_index=cell_to_index, nkp_grid=nkp_grid, lattice_fft=lattice_fft)
765 cpassert(ASSOCIATED(sab_nl))
766 kplocal = kp_range(2) - kp_range(1) + 1
767 kpmax = maxval(kp_dist(2, :) - kp_dist(1, :) + 1)
768
769 !Whether we use DIIS for k-points
770 do_diis = .false.
771 use_adiis = update_p .AND. scf_control%diagonalization%update_method == diag_update_method_adiis
772 IF (use_adiis) THEN
773 cpassert(ASSOCIATED(scf_env%scf_subspace_buffer))
774 END IF
775 adiis_cdiis = .false.
776 IF (use_adiis .AND. ASSOCIATED(kpoints%scf_diis_buffer)) THEN
777 IF (scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis .AND. &
778 .NOT. scf_env%adiis_check_next .AND. &
779 .NOT. scf_env%scf_subspace_buffer%last_restart .AND. &
780 scf_env%scf_subspace_buffer%nstored > 0 .AND. &
781 PRESENT(diis_error) .AND. PRESENT(qs_env)) THEN
782 adiis_cdiis = scf_env%iter_delta < scf_control%eps_diis .OR. &
783 kpoints%scf_diis_buffer%ncall > 0
784 END IF
785 END IF
786 IF ((.NOT. use_adiis .AND. scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis .AND. &
787 .NOT. use_real_wfn .AND. PRESENT(diis_error) .AND. PRESENT(qs_env)) .OR. &
788 adiis_cdiis) do_diis = .true.
789
790 adiis_ks => matrix_ks
791 NULLIFY (adiis_u)
792 IF (ALLOCATED(kpoints%lowdin_v)) adiis_u => kpoints%lowdin_v
793 IF (use_adiis) THEN
794 IF (scf_env%scf_subspace_buffer%use_combined_fock) THEN
795 adiis_ks => scf_env%scf_subspace_buffer%combined_fock
796 IF (ASSOCIATED(adiis_u)) adiis_u => scf_env%scf_subspace_buffer%combined_aux
797 END IF
798 END IF
799 effective_ks => matrix_ks
800 NULLIFY (effective_u)
801 IF (ALLOCATED(kpoints%lowdin_v)) effective_u => kpoints%lowdin_v
802 IF (PRESENT(potential)) THEN
803 cpassert(.NOT. update_p)
804 effective_u => potential
805 END IF
806 IF (use_adiis .AND. .NOT. adiis_cdiis) THEN
807 effective_ks => adiis_ks
808 effective_u => adiis_u
809 END IF
810 IF (use_adiis .AND. .NOT. do_diis) THEN
811 scf_env%scf_subspace_buffer%diis_weight = 0.0_dp
812 IF (ASSOCIATED(kpoints%scf_diis_buffer)) THEN
813 IF (kpoints%scf_diis_buffer%ncall > 0) CALL qs_diis_b_clear_kp(kpoints%scf_diis_buffer)
814 END IF
815 END IF
816
817 para_env => kpoints%blacs_env_all%para_env
818 nspin = SIZE(matrix_ks, 1)
819 omp_threads = 1
820!$ omp_threads = omp_get_max_threads()
821 ! One-rank groups use the same worker path for MPI, OpenMP and hybrid jobs.
822 use_local_kpoints = kpoints%para_env_kp%num_pe == 1 .AND. .NOT. use_real_wfn
823 IF (diag_check_requested()) use_local_kpoints = .false.
824 ib = 0
825 nb = 0
826 NULLIFY (scf_section)
827
828 NULLIFY (rmatrix)
829
830 fmwork => scf_env%scf_work1
831 ! Both real and complex operators share the same assembly and transfer context.
832 CALL kpoint_operator_context_create(op_ctx, kpoints, fmwork(1), nspin)
833
834 ! fm pools to be used within a kpoint group
835 CALL get_kpoint_info(kpoints, mpools=mpools)
836 CALL mpools_get(mpools, ao_ao_fm_pool=ao_ao_fm_pool)
837
838 CALL fm_pool_create_fm(ao_ao_fm_pool, fmlocal)
839 CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
840
841 IF (use_real_wfn) THEN
842 CALL cp_fm_create(rksmat, matrix_struct)
843 CALL cp_fm_create(rsmat, matrix_struct)
844 END IF
845 ! Retain the unmodified S(k) for this SCF run, outside the driver-local
846 ! operator context. Solvers receive work copies because they overwrite S.
847 cache_smat = update_p .AND. .NOT. use_real_wfn
848 reuse_smat = .false.
849 IF (cache_smat) THEN
850 reuse_smat = ALLOCATED(kpoints%scf_smat)
851 IF (.NOT. reuse_smat) THEN
852 ALLOCATE (kpoints%scf_smat(kplocal))
853 DO ikp = 1, kplocal
854 CALL cp_cfm_create(kpoints%scf_smat(ikp), matrix_struct)
855 END DO
856 END IF
857 cpassert(SIZE(kpoints%scf_smat) == kplocal)
858 END IF
859
860 ! Prepare H/S before starting any round. ADIIS can replace only the idle KS slot.
861 prepared_ops = .false.
862 ! A single distributed group can also assemble its prepared local tiles.
863 CALL kpoint_operator_prepare(op_ctx, effective_ks, kpoint_slot_ks, prepared_ops)
864 IF (.NOT. reuse_smat) THEN
865 CALL kpoint_operator_prepare(op_ctx, matrix_s, kpoint_slot_s, prepared_part, matrix_row=kpoint_spin_free)
866 prepared_ops = prepared_ops .AND. prepared_part
867 END IF
868
869 use_local_kpoints = use_local_kpoints .AND. prepared_ops
870 nworkers = 1
871 IF (use_local_kpoints) nworkers = min(omp_threads, kplocal)
872 IF (.NOT. use_real_wfn .OR. do_diis) THEN
873 kp => kpoints%kp_env(1)%kpoint_env
874 IF (use_real_wfn) THEN
875 CALL get_mo_set(kp%mos(1), mo_coeff=mo_coeff)
876 CALL cp_fm_get_info(mo_coeff, matrix_struct=mo_struct)
877 ELSE
878 mo_struct => kp%mos(1)%cmo_coeff%matrix_struct
879 END IF
880 CALL kp_diag_workspaces_create(workspace, nworkers, matrix_struct, mo_struct, 0, &
881 use_local_kpoints, scf_env%cholesky_method == cholesky_off)
882 cksmat => workspace(1)%cksmat
883 csmat => workspace(1)%csmat
884 cwork => workspace(1)%cwork
885 cmos => workspace(1)%cmos
886 END IF
887
888 ! Prepared operators select FFT or direct assembly inside their common route.
889 ! Only an unprepared operator needs the sparse-grid/transfer fallback below.
890 ! Divide its conservative memory budget between H spins and S.
891 ALLOCATE (ks_grid(nspin))
892 SELECT CASE (lattice_fft)
893 CASE (lattice_fft_auto)
894 use_grid_fft = nkp >= 27
895 CASE (lattice_fft_on)
896 use_grid_fft = .true.
897 CASE (lattice_fft_off)
898 use_grid_fft = .false.
899 CASE DEFAULT
900 cpabort("Unknown lattice FFT selection mode")
901 END SELECT
902 use_grid_fft = use_grid_fft .AND. .NOT. prepared_ops
903 ! The driver needs sparse scratch only for forward lattice-FFT preparation.
904 IF (use_grid_fft) THEN
905 ALLOCATE (rmatrix)
906 CALL dbcsr_create(rmatrix, template=matrix_ks(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
907 CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
908 END IF
909 IF (use_grid_fft .AND. .NOT. reuse_smat) THEN
910 CALL rskp_transform_grid_prepare(s_grid, rmatrix, matrix_s, 1, xkp, nkp_grid, &
911 cell_to_index, sab_nl, use_grid_fft, &
912 max_storage_bytes=512_int_8*1024_int_8**2/int(nspin + 1, int_8))
913 fft_ready = merge(1, 0, use_grid_fft)
914 CALL para_env%min(fft_ready)
915 use_grid_fft = fft_ready == 1
916 END IF
917 IF (use_grid_fft) THEN
918 DO ispin = 1, nspin
919 CALL rskp_transform_grid_prepare(ks_grid(ispin), rmatrix, effective_ks, ispin, xkp, nkp_grid, &
920 cell_to_index, sab_nl, spin_grid_fft, &
921 max_storage_bytes=512_int_8*1024_int_8**2/int(nspin + 1, int_8))
922 use_grid_fft = use_grid_fft .AND. spin_grid_fft
923 END DO
924 fft_ready = merge(1, 0, use_grid_fft)
925 CALL para_env%min(fft_ready)
926 use_grid_fft = fft_ready == 1
927 END IF
928 IF (.NOT. use_grid_fft) THEN
929 CALL rskp_transform_grid_release(s_grid)
930 DO ispin = 1, nspin
931 CALL rskp_transform_grid_release(ks_grid(ispin))
932 END DO
933 END IF
934
935 ! Setup and start all the communication.
936 IF (.NOT. use_local_kpoints) THEN
937 DO ikp = 1, kpmax
938 DO ispin = 1, nspin
939 CALL kpoint_operator_start(op_ctx, ikp, ispin, effective_ks, kpoint_slot_ks, &
940 ks_grid(ispin), use_grid_fft, matrix_row=ispin)
941 IF (.NOT. cache_smat .OR. (.NOT. reuse_smat .AND. ispin == 1)) THEN
942 CALL kpoint_operator_start(op_ctx, ikp, ispin, matrix_s, kpoint_slot_s, s_grid, &
943 use_grid_fft, matrix_row=kpoint_spin_free)
944 END IF
945 END DO
946 END DO
947 END IF
948
949 ! Finish communication then diagonalise in each group
950 diis_weight = 1.0_dp
951 IF (do_diis) THEN
952 scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
953 cdiis_ncall = kpoints%scf_diis_buffer%ncall
954 CALL qs_diis_b_info_kp(kpoints%scf_diis_buffer, ib, nb)
955 CALL qs_diis_b_check_i_alloc_kp(kpoints%scf_diis_buffer, matrix_struct, nspin, kplocal, scf_section)
956 IF (use_local_kpoints) THEN
957 CALL kp_diag_local_phase(op_ctx, kpoints, workspace, scf_env, scf_control, nspin, &
958 cache_smat, reuse_smat, use_adiis, ib, .true., 1.0_dp, potential=effective_u)
959 ELSE
960 DO ikp = 1, kpmax
961 DO ispin = 1, nspin
962 CALL harvest_operators(ikp, ispin)
963 IF (ikp > kplocal) cycle
964 CALL qs_diis_b_calc_err_kp(kpoints%scf_diis_buffer, ib, kp%mos, cksmat, csmat, &
965 ispin, ikp, cmos, real_wfn=use_real_wfn, &
966 local=.false., gemm_ctx=workspace(1)%gemm_ctx)
967
968 END DO !ispin
969 END DO !ikp
970
971 END IF
972
973 ALLOCATE (coeffs(nb))
974 CALL qs_diis_b_step_kp(kpoints%scf_diis_buffer, coeffs, ib, nb, scf_env%iter_delta, diis_error, &
975 cdiis_step, scf_control%eps_diis, nspin, nkp, kplocal, &
976 merge(2, scf_control%nmixing, use_adiis), &
977 scf_section, kpoints%para_env_inter_kp)
978 diis_weight = 1.0_dp
979 IF (use_adiis) THEN
980 diis_weight = scf_env%scf_subspace_buffer%diis_weight
981 IF (.NOT. ieee_is_finite(diis_error) .OR. diis_error >= scf_control%eps_diis) THEN
982 cdiis_step = .false.
983 IF (cdiis_ncall < kpoints%scf_diis_buffer%nbuffer) THEN
984 ! The rejected state has not overwritten a full circular buffer.
985 kpoints%scf_diis_buffer%ncall = cdiis_ncall
986 ELSE
987 CALL qs_diis_b_clear_kp(kpoints%scf_diis_buffer)
988 diis_weight = 0.0_dp
989 END IF
990 END IF
991 IF (cdiis_step) THEN
992 ! Use one half-weight handover step; MAX_DIIS only controls the buffer capacity.
993 diis_weight = min(1.0_dp, diis_weight + 0.5_dp)
994 ELSE
995 ! A rejected or deferred CDIIS step must not leave a partial raw-Fock blend active.
996 diis_weight = 0.0_dp
997 END IF
998 scf_env%scf_subspace_buffer%diis_weight = diis_weight
999 diis_step = cdiis_step .AND. diis_weight >= 1.0_dp
1000 IF (.NOT. diis_step) diis_error = scf_env%scf_subspace_buffer%last_old_fock_weight
1001 ELSE
1002 diis_step = cdiis_step
1003 END IF
1004
1005 ! All first-pass instances have finished. Capture the selected
1006 ! ADIIS operator once instead of retransferring it in every round.
1007 IF (use_adiis .AND. diis_weight < 1.0_dp) THEN
1008 CALL kpoint_operator_prepare(op_ctx, adiis_ks, kpoint_slot_ks, prepared_part)
1009 ! A collective prepare rejection must return all ranks to the transfer path.
1010 use_local_kpoints = use_local_kpoints .AND. prepared_part
1011 END IF
1012
1013 ELSE
1014 IF (.NOT. update_p .OR. use_adiis) diis_step = .false.
1015 END IF
1016
1017 IF (use_local_kpoints) THEN
1018 ! Unallocated arrays and disassociated pointers are absent at nonpointer optional dummies.
1019 IF (do_diis) effective_u => adiis_u
1020 CALL kp_diag_local_phase(op_ctx, kpoints, workspace, scf_env, scf_control, nspin, &
1021 cache_smat, reuse_smat, use_adiis, ib, .false., diis_weight, coeffs, effective_u)
1022 ELSE
1023 DO ikp = 1, kpmax
1024 DO ispin = 1, nspin
1025 IF (do_diis) THEN
1026 IF (use_adiis .AND. diis_weight < 1.0_dp) THEN
1027 CALL kpoint_operator_start(op_ctx, ikp, ispin, adiis_ks, kpoint_slot_ks, &
1028 use_grid=.false., matrix_row=ispin)
1029 CALL kpoint_operator_finish(op_ctx, ikp, ispin, kpoint_slot_ks, cwork)
1030 IF (ikp <= kplocal .AND. ASSOCIATED(adiis_u)) THEN
1031 CALL lowdin_kp_apply(kpoints%kp_env(ikp)%kpoint_env, adiis_u(:, ispin), cmat=cwork)
1032 END IF
1033 END IF
1034 IF (ikp > kplocal) cycle
1035 CALL cp_cfm_to_cfm(kpoints%scf_diis_buffer%smat(ikp), csmat)
1036 CALL kp_diis_candidate(kpoints, ikp, ispin, coeffs, diis_weight, cksmat, cwork)
1037 IF (use_real_wfn) THEN
1038 CALL cp_cfm_to_fm(cksmat, rksmat, fmlocal)
1039 CALL cp_cfm_to_fm(csmat, rsmat, fmlocal)
1040 END IF
1041 ELSE
1042 CALL harvest_operators(ikp, ispin)
1043 IF (ikp > kplocal) cycle
1044 END IF
1045 kp => kpoints%kp_env(ikp)%kpoint_env
1046 IF (use_real_wfn) THEN
1047 CALL get_mo_set(kp%mos(ispin), mo_coeff=mo_coeff, eigenvalues=eigenvalues)
1048 IF (use_adiis) CALL shift_candidate_real(rksmat, rsmat, kp%mos(ispin), &
1049 adiis_candidate_shift(scf_env))
1050 IF (scf_env%cholesky_method == cholesky_off) THEN
1051 CALL cp_fm_geeig_canon(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal, scf_control%eps_eigval)
1052 ELSE
1053 CALL cp_fm_geeig(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal)
1054 END IF
1055 ELSE
1056 CALL kp_solve_complex(kp%mos(ispin), workspace(1), scf_env, scf_control, use_adiis, .false.)
1057 END IF
1058 END DO
1059 END DO
1060 END IF
1061 IF (ALLOCATED(coeffs)) DEALLOCATE (coeffs)
1062
1063 CALL kpoint_operator_context_release(op_ctx)
1064
1065 IF (update_p) THEN
1066 ! MO occupations
1067 IF (PRESENT(probe) .EQV. .true.) THEN
1068 scf_control%smear%do_smear = .false.
1069 CALL kpoint_set_mo_occupation(kpoints, scf_control%smear, &
1070 probe=probe)
1071 ELSE
1072 CALL kpoint_set_mo_occupation(kpoints, scf_control%smear, &
1073 added_mos_auto=scf_control%added_mos_auto, &
1074 added_mos_auto_grow=my_added_mos_auto_grow)
1075 END IF
1076 IF (my_added_mos_auto_grow) THEN
1077 IF (PRESENT(added_mos_auto_grow)) added_mos_auto_grow = .true.
1078 ELSE
1079 ! density matrices
1080 CALL kpoint_density_matrices(kpoints)
1081 ! density matrices in real space
1082 CALL kpoint_density_transform(kpoints, scf_env%p_mix_new, .false., &
1083 matrix_s(1, 1)%matrix, sab_nl, fmwork, overlap_rs=matrix_s)
1084 END IF
1085 END IF
1086
1087 IF (ASSOCIATED(rmatrix)) CALL dbcsr_deallocate_matrix(rmatrix)
1088 CALL rskp_transform_grid_release(s_grid)
1089 DO ispin = 1, nspin
1090 CALL rskp_transform_grid_release(ks_grid(ispin))
1091 END DO
1092 DEALLOCATE (ks_grid)
1093
1094 IF (use_real_wfn) THEN
1095 CALL cp_fm_release(rksmat)
1096 CALL cp_fm_release(rsmat)
1097 END IF
1098 IF (ALLOCATED(workspace)) CALL kp_diag_workspaces_release(workspace)
1099 CALL fm_pool_give_back_fm(ao_ao_fm_pool, fmlocal)
1100
1101 CALL timestop(handle)
1102
1103 CONTAINS
1104
1105 ! Serialized harvest only. Inactive tail rounds participate in every transfer,
1106 ! then return before reading a local matrix or publishing the overlap cache.
1107! **************************************************************************************************
1108!> \brief ...
1109!> \param round ...
1110!> \param spin ...
1111! **************************************************************************************************
1112 SUBROUTINE harvest_operators(round, spin)
1113 INTEGER, INTENT(IN) :: round, spin
1114
1115 IF (use_real_wfn) THEN
1116 CALL kpoint_operator_finish(op_ctx, round, spin, kpoint_slot_ks, fm_re=rksmat)
1117 CALL kpoint_operator_finish(op_ctx, round, spin, kpoint_slot_s, fm_re=rsmat)
1118 IF (round <= kplocal .AND. ASSOCIATED(effective_u)) THEN
1119 CALL lowdin_kp_apply(kpoints%kp_env(round)%kpoint_env, effective_u(:, spin), rmat=rksmat)
1120 END IF
1121 IF (round <= kplocal .AND. do_diis) THEN
1122 CALL cp_cfm_scale_and_add_fm(z_zero, cksmat, z_one, rksmat)
1123 CALL cp_cfm_scale_and_add_fm(z_zero, csmat, z_one, rsmat)
1124 END IF
1125 ELSE
1126 CALL kpoint_operator_finish(op_ctx, round, spin, kpoint_slot_ks, cksmat)
1127 IF (round <= kplocal .AND. ASSOCIATED(effective_u)) THEN
1128 CALL lowdin_kp_apply(kpoints%kp_env(round)%kpoint_env, effective_u(:, spin), cmat=cksmat)
1129 END IF
1130 IF (.NOT. cache_smat .OR. (.NOT. reuse_smat .AND. spin == 1)) THEN
1131 CALL kpoint_operator_finish(op_ctx, round, spin, kpoint_slot_s, csmat)
1132 END IF
1133 END IF
1134
1135 IF (round > kplocal) RETURN
1136 kp => kpoints%kp_env(round)%kpoint_env
1137 IF (cache_smat) THEN
1138 IF (.NOT. reuse_smat .AND. spin == 1) THEN
1139 CALL cp_cfm_to_cfm(csmat, kpoints%scf_smat(round))
1140 ELSE
1141 CALL cp_cfm_to_cfm(kpoints%scf_smat(round), csmat)
1142 END IF
1143 END IF
1144 END SUBROUTINE harvest_operators
1145
1146 END SUBROUTINE do_general_diag_kp
1147
1148! **************************************************************************************************
1149!> \brief Execute a local numerical phase. All SCF decisions and collectives belong to the caller.
1150!> \param ctx immutable prepared operators
1151!> \param kpoints host environment
1152!> \param workspace private dense/phase/LAPACK storage per worker
1153!> \param scf_env SCF settings (read only during this phase)
1154!> \param scf_control SCF control
1155!> \param nspin spin count
1156!> \param cache_smat retain the unmodified overlap
1157!> \param reuse_smat overlap cache is already valid
1158!> \param use_adiis shift an ADIIS candidate only
1159!> \param ib DIIS insertion index
1160!> \param store_error store raw F and its residual instead of solving
1161!> \param diis_weight weight of the CDIIS candidate
1162!> \param coeffs optional DIIS coefficients for the second phase
1163!> \param potential optional Lowdin orbital potential paired with the prepared Hamiltonian
1164! **************************************************************************************************
1165 SUBROUTINE kp_diag_local_phase(ctx, kpoints, workspace, scf_env, scf_control, nspin, &
1166 cache_smat, reuse_smat, use_adiis, ib, store_error, diis_weight, coeffs, potential)
1167 TYPE(kpoint_operator_context_type), INTENT(IN) :: ctx
1168 TYPE(kpoint_type), POINTER :: kpoints
1169 TYPE(kp_diag_workspace_type), TARGET, INTENT(INOUT) :: workspace(:)
1170 TYPE(qs_scf_env_type), POINTER :: scf_env
1171 TYPE(scf_control_type), POINTER :: scf_control
1172 INTEGER, INTENT(IN) :: nspin, ib
1173 LOGICAL, INTENT(IN) :: cache_smat, reuse_smat, use_adiis, store_error
1174 REAL(kind=dp), INTENT(IN) :: diis_weight
1175 COMPLEX(KIND=dp), OPTIONAL, INTENT(IN) :: coeffs(:)
1176 REAL(dp), OPTIONAL, INTENT(IN) :: potential(:, :)
1177
1178 COMPLEX(KIND=dp), ALLOCATABLE :: operators(:, :, :, :)
1179 REAL(kind=dp), ALLOCATABLE :: block_work(:)
1180 INTEGER :: first, ik, ikp, ispin, j, kplocal, nbatch, nmatrices, nworkers, thread, width
1181 INTEGER(KIND=int_8) :: work_size
1182 LOGICAL :: build_h, build_s
1183 TYPE(kpoint_env_type), POINTER :: kp
1184 TYPE(kp_diag_workspace_type), POINTER :: work
1185
1186 kplocal = kpoints%kp_range(2) - kpoints%kp_range(1) + 1
1187 IF (kplocal == 0) RETURN
1188 nworkers = SIZE(workspace)
1189 build_s = .NOT. PRESENT(coeffs) .AND. .NOT. reuse_smat
1190 build_h = .NOT. PRESENT(coeffs) .OR. diis_weight < 1.0_dp
1191 nmatrices = merge(1, 0, build_h) + merge(1, 0, build_s)
1192 width = 1
1193 work_size = 0_int_8
1194 IF (nmatrices > 0) CALL kpoint_operator_batch_size(ctx, nworkers, nmatrices, width, work_size)
1195 ! Each worker assembles a bounded k tile, then consumes it using its private solver.
1196 ! The second CDIIS phase combines history without allocating operator tiles.
1197!$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nworkers) IF(nworkers > 1) &
1198!$OMP SHARED(ctx,kpoints,workspace,scf_env,scf_control,nspin,cache_smat,reuse_smat,use_adiis,ib, &
1199!$OMP store_error,diis_weight,coeffs,potential,kplocal,nworkers,build_h,build_s,nmatrices,width,work_size) &
1200!$OMP PRIVATE(first,ik,ikp,ispin,j,nbatch,thread,kp,work,operators,block_work)
1201 thread = 1
1202!$ thread = omp_get_thread_num() + 1
1203 work => workspace(thread)
1204 ALLOCATE (operators(SIZE(work%cksmat%local_data, 1), SIZE(work%cksmat%local_data, 2), width, nmatrices), &
1205 block_work(work_size))
1206!$OMP DO SCHEDULE(DYNAMIC, 1)
1207 DO first = 1, kplocal, width
1208 nbatch = min(width, kplocal - first + 1)
1209 ik = kpoints%kp_range(1) + first - 1
1210 IF (build_s) THEN
1211 CALL kpoint_operator_evaluate_batch(ctx, 1, kpoint_slot_s, ik, operators(:, :, 1:nbatch, nmatrices), block_work)
1212 END IF
1213 DO ispin = 1, nspin
1214 IF (build_h) THEN
1215 CALL kpoint_operator_evaluate_batch(ctx, ispin, kpoint_slot_ks, ik, operators(:, :, 1:nbatch, 1), block_work)
1216 END IF
1217 DO j = 1, nbatch
1218 ikp = first + j - 1
1219 kp => kpoints%kp_env(ikp)%kpoint_env
1220 IF (build_h) THEN
1221 IF (PRESENT(coeffs)) THEN
1222 work%csmat_base%local_data(:, :) = operators(:, :, j, 1)
1223 IF (PRESENT(potential)) CALL lowdin_kp_apply(kp, potential(:, ispin), cmat=work%csmat_base)
1224 ELSE
1225 work%cksmat%local_data(:, :) = operators(:, :, j, 1)
1226 IF (PRESENT(potential)) CALL lowdin_kp_apply(kp, potential(:, ispin), cmat=work%cksmat)
1227 END IF
1228 END IF
1229 IF (PRESENT(coeffs)) THEN
1230 CALL cp_cfm_to_cfm(kpoints%scf_diis_buffer%smat(ikp), work%csmat)
1231 CALL kp_diis_candidate(kpoints, ikp, ispin, coeffs, diis_weight, work%cksmat, work%csmat_base)
1232 ELSE
1233 IF (reuse_smat) THEN
1234 CALL cp_cfm_to_cfm(kpoints%scf_smat(ikp), work%csmat)
1235 ELSE
1236 work%csmat%local_data(:, :) = operators(:, :, j, nmatrices)
1237 END IF
1238 IF (cache_smat .AND. .NOT. reuse_smat .AND. ispin == 1) THEN
1239 CALL cp_cfm_to_cfm(work%csmat, kpoints%scf_smat(ikp))
1240 END IF
1241 END IF
1242 IF (store_error) THEN
1243 CALL qs_diis_b_calc_err_kp(kpoints%scf_diis_buffer, ib, kp%mos, work%cksmat, work%csmat, &
1244 ispin, ikp, work%cmos, local=.true., gemm_ctx=work%gemm_ctx)
1245 ELSE
1246 CALL kp_solve_complex(kp%mos(ispin), work, scf_env, scf_control, use_adiis, .true.)
1247 END IF
1248 END DO
1249 END DO
1250 END DO
1251!$OMP END DO
1252 DEALLOCATE (operators, block_work)
1253!$OMP END PARALLEL
1254
1255 END SUBROUTINE kp_diag_local_phase
1256
1257! **************************************************************************************************
1258!> \brief Form the same candidate on local and distributed layouts, preserving raw DIIS history.
1259!> \param kpoints host DIIS buffer
1260!> \param ikp local k-point index
1261!> \param ispin spin index
1262!> \param coeffs DIIS coefficients
1263!> \param weight CDIIS blend weight
1264!> \param fock output candidate
1265!> \param adiis_fock ADIIS candidate, read only when weight is less than one
1266! **************************************************************************************************
1267 SUBROUTINE kp_diis_candidate(kpoints, ikp, ispin, coeffs, weight, fock, adiis_fock)
1268 TYPE(kpoint_type), POINTER :: kpoints
1269 INTEGER, INTENT(IN) :: ikp, ispin
1270 COMPLEX(KIND=dp), INTENT(IN) :: coeffs(:)
1271 REAL(kind=dp), INTENT(IN) :: weight
1272 TYPE(cp_cfm_type), INTENT(INOUT) :: fock
1273 TYPE(cp_cfm_type), INTENT(IN) :: adiis_fock
1274
1275 INTEGER :: jb
1276
1277 CALL cp_cfm_set_all(fock, z_zero)
1278 DO jb = 1, SIZE(coeffs)
1279 CALL cp_cfm_scale_and_add(z_one, fock, coeffs(jb), kpoints%scf_diis_buffer%param(jb, ispin, ikp))
1280 END DO
1281 IF (weight < 1.0_dp) CALL cp_cfm_scale_and_add(cmplx(weight, kind=dp), fock, &
1282 cmplx(1.0_dp - weight, kind=dp), adiis_fock)
1283 END SUBROUTINE kp_diis_candidate
1284
1285! **************************************************************************************************
1286!> \brief Shift, solve and publish a complex SCF candidate using private scratch.
1287!> \param mos complex MO set
1288!> \param work private workspace
1289!> \param scf_env SCF settings
1290!> \param scf_control SCF control
1291!> \param use_adiis apply the candidate shift
1292!> \param local select the one-rank local solver; otherwise honor the configured backend
1293! **************************************************************************************************
1294 SUBROUTINE kp_solve_complex(mos, work, scf_env, scf_control, use_adiis, local)
1295 TYPE(mo_set_type), INTENT(INOUT) :: mos
1296 TYPE(kp_diag_workspace_type), INTENT(INOUT) :: work
1297 TYPE(qs_scf_env_type), POINTER :: scf_env
1298 TYPE(scf_control_type), POINTER :: scf_control
1299 LOGICAL, INTENT(IN) :: use_adiis, local
1300
1301 REAL(kind=dp), POINTER :: eigenvalues(:)
1302
1303 CALL get_mo_set(mos, eigenvalues=eigenvalues)
1304 IF (use_adiis) CALL shift_candidate_complex(work%cksmat, work%csmat, mos, work%cmos, work%cwork, &
1305 adiis_candidate_shift(scf_env), local, work%gemm_ctx)
1306 IF (scf_env%cholesky_method == cholesky_off) THEN
1307 IF (local) THEN
1308 CALL cp_cfm_geeig_canon_local(work%cksmat, work%csmat, mos%cmo_coeff, eigenvalues, work%cwork, &
1309 scf_control%eps_eigval, ws=work%solver)
1310 ELSE
1311 CALL cp_cfm_geeig_canon(work%cksmat, work%csmat, mos%cmo_coeff, eigenvalues, work%cwork, &
1312 scf_control%eps_eigval)
1313 END IF
1314 ELSE
1315 IF (local) THEN
1316 CALL cp_cfm_geeig_local(work%cksmat, work%csmat, mos%cmo_coeff, eigenvalues, ws=work%solver)
1317 ELSE
1318 CALL cp_cfm_geeig(work%cksmat, work%csmat, mos%cmo_coeff, eigenvalues, work%cwork)
1319 END IF
1320 END IF
1321 END SUBROUTINE kp_solve_complex
1322
1323! **************************************************************************************************
1324!> \brief Allocate all worker storage before launching numerical phases.
1325!> \param workspace returned private workspaces
1326!> \param nworkers number of workers
1327!> \param matrix_struct AO layout of the k-point group
1328!> \param mo_struct MO layout of the k-point group
1329!> \param nentry number of phase-sum entries
1330!> \param local allocate local phase and solver scratch
1331!> \param canonical use the canonical solver
1332!> \param lanczos allocate the block-Lanczos kernel scratch
1333!> \param nkrylov number of block expansion steps, sizes the kernel scratch
1334!> \param nblock number of columns refined jointly, sizes the kernel scratch
1335! **************************************************************************************************
1336 SUBROUTINE kp_diag_workspaces_create(workspace, nworkers, matrix_struct, mo_struct, nentry, local, canonical, &
1337 lanczos, nkrylov, nblock)
1338 TYPE(kp_diag_workspace_type), ALLOCATABLE, &
1339 INTENT(OUT) :: workspace(:)
1340 INTEGER, INTENT(IN) :: nworkers
1341 TYPE(cp_fm_struct_type), POINTER :: matrix_struct, mo_struct
1342 INTEGER, INTENT(IN) :: nentry
1343 LOGICAL, INTENT(IN) :: local, canonical
1344 LOGICAL, INTENT(IN), OPTIONAL :: lanczos
1345 INTEGER, INTENT(IN), OPTIONAL :: nkrylov, nblock
1346
1347 INTEGER :: i, n
1348 LOGICAL :: my_lanczos
1349
1350 my_lanczos = .false.
1351 IF (PRESENT(lanczos)) my_lanczos = lanczos
1352
1353 ALLOCATE (workspace(nworkers))
1354 DO i = 1, nworkers
1355 CALL cp_cfm_create(workspace(i)%cksmat, matrix_struct)
1356 CALL cp_cfm_create(workspace(i)%csmat, matrix_struct)
1357 CALL cp_cfm_create(workspace(i)%cwork, matrix_struct)
1358 CALL cp_cfm_create(workspace(i)%cmos, mo_struct)
1359 ! the transfer rounds of the Lanczos driver also need the stored
1360 ! overlap copy, although they never run the local solvers
1361 IF (local .OR. my_lanczos) THEN
1362 CALL cp_cfm_create(workspace(i)%csmat_base, matrix_struct)
1363 END IF
1364 IF (local) THEN
1365 CALL workspace(i)%gemm_ctx%create(local_gemm_pu_gpu, timing=.false.)
1366 ALLOCATE (workspace(i)%phase(1, nentry))
1367 n = matrix_struct%nrow_global
1368 IF (canonical) THEN
1369 CALL cp_cfm_local_workspace_prepare(workspace(i)%solver, n, workspace(i)%cksmat%local_data)
1370 ELSE
1371 CALL cp_cfm_local_workspace_prepare(workspace(i)%solver, n, workspace(i)%cksmat%local_data, &
1372 workspace(i)%csmat%local_data)
1373 END IF
1374 END IF
1375 IF (my_lanczos) THEN
1376 CALL lanczos_workspace_create(workspace(i)%lan_ws, mo_struct, matrix_struct, nkrylov, nblock)
1377 END IF
1378 END DO
1379 END SUBROUTINE kp_diag_workspaces_create
1380
1381! **************************************************************************************************
1382!> \brief Release dense matrices and automatic phase/LAPACK allocations after all workers finish.
1383!> \param workspace private workspaces
1384! **************************************************************************************************
1385 SUBROUTINE kp_diag_workspaces_release(workspace)
1386 TYPE(kp_diag_workspace_type), ALLOCATABLE, &
1387 INTENT(INOUT) :: workspace(:)
1388
1389 INTEGER :: i
1390
1391 DO i = 1, SIZE(workspace)
1392 CALL workspace(i)%gemm_ctx%destroy()
1393 CALL cp_cfm_release(workspace(i)%cksmat)
1394 CALL cp_cfm_release(workspace(i)%csmat)
1395 CALL cp_cfm_release(workspace(i)%csmat_base)
1396 CALL cp_cfm_release(workspace(i)%cwork)
1397 CALL cp_cfm_release(workspace(i)%cmos)
1398 CALL lanczos_workspace_release(workspace(i)%lan_ws)
1399 END DO
1400 DEALLOCATE (workspace)
1401 END SUBROUTINE kp_diag_workspaces_release
1402
1403! **************************************************************************************************
1404!> \brief Adaptive inner tolerance of an iterative solver from EPS_ADAPT:
1405!> loose while the density error sits far from self-consistency,
1406!> tracking that error at a fixed fraction and capped in the early
1407!> phase. A nonpositive factor or error disables the loosening and
1408!> returns the user threshold
1409!> \param eps_iter user threshold
1410!> \param eps_adapt loosening factor of the tracking
1411!> \param iter_error outer error the threshold tracks
1412!> \param noise_scale fraction of the outer error the solver noise keeps
1413!> \param cap bound of the early loose phase
1414!> \return ...
1415! **************************************************************************************************
1416 PURE FUNCTION kp_adaptive_eps_iter(eps_iter, eps_adapt, iter_error, noise_scale, cap) RESULT(eps)
1417 REAL(kind=dp), INTENT(IN) :: eps_iter, eps_adapt, iter_error, &
1418 noise_scale, cap
1419 REAL(kind=dp) :: eps
1420
1421 IF (eps_adapt <= 0.0_dp .OR. iter_error <= 0.0_dp) THEN
1422 eps = eps_iter
1423 ELSE
1424 eps = max(eps_iter, min(eps_adapt*iter_error/noise_scale, cap))
1425 END IF
1426 END FUNCTION kp_adaptive_eps_iter
1427
1428! **************************************************************************************************
1429!> \brief Occupations, ADDED_MOS_AUTO growth check and the density exit
1430!> chain shared by the k-point diagonalization drivers. The
1431!> routine keeps its own bookkeeping: a growing step abandons the
1432!> density of the abandoned SCF step
1433!> \param kpoints ...
1434!> \param scf_env ...
1435!> \param scf_control ...
1436!> \param matrix_s_kp real-space overlap matrix, dbcsr image set
1437!> \param sab_nl ...
1438!> \param fmwork global scratch fms of the driver
1439!> \param grow TRUE when the virtual-space buffer is too small for the
1440!> smearing tail: the outer loop grows the MOS and redoes this
1441!> SCF step
1442! **************************************************************************************************
1443 SUBROUTINE kp_diag_exit_chain(kpoints, scf_env, scf_control, matrix_s_kp, sab_nl, fmwork, grow)
1444 TYPE(kpoint_type), POINTER :: kpoints
1445 TYPE(qs_scf_env_type), POINTER :: scf_env
1446 TYPE(scf_control_type), POINTER :: scf_control
1447 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp
1448 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1449 POINTER :: sab_nl
1450 TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
1451 LOGICAL, INTENT(OUT) :: grow
1452
1453 CALL kpoint_set_mo_occupation(kpoints, scf_control%smear, &
1454 added_mos_auto=scf_control%added_mos_auto, &
1455 added_mos_auto_grow=grow)
1456 IF (grow) RETURN
1457 CALL kpoint_density_matrices(kpoints)
1458 CALL kpoint_density_transform(kpoints, scf_env%p_mix_new, .false., &
1459 matrix_s_kp(1, 1)%matrix, sab_nl, fmwork, overlap_rs=matrix_s_kp)
1460 END SUBROUTINE kp_diag_exit_chain
1461
1462! **************************************************************************************************
1463!> \brief Per-column convergence thresholds of an iterative k-point
1464!> solver: the occupation decides which columns sit beyond the
1465!> thermal smearing tail and converge at the relaxed empty
1466!> threshold
1467!> \param occupation occupation numbers of the channel
1468!> \param eps_iter threshold of the occupied columns
1469!> \param eps_iter_empty relaxed threshold of the unoccupied columns
1470!> \param eps_col per-column thresholds, at least as long as the
1471!> occupation numbers
1472! **************************************************************************************************
1473 SUBROUTINE kp_eps_col_by_occupation(occupation, eps_iter, eps_iter_empty, eps_col)
1474 REAL(kind=dp), DIMENSION(:), INTENT(IN), POINTER :: occupation
1475 REAL(kind=dp), INTENT(IN) :: eps_iter, eps_iter_empty
1476 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eps_col
1477
1478 INTEGER :: imo
1479
1480 DO imo = 1, SIZE(occupation)
1481 IF (occupation(imo) < occ_tol) THEN
1482 eps_col(imo) = eps_iter_empty
1483 ELSE
1484 eps_col(imo) = eps_iter
1485 END IF
1486 END DO
1487 END SUBROUTINE kp_eps_col_by_occupation
1488
1489! **************************************************************************************************
1490!> \brief Common setup of the iterative k-point drivers: the global
1491!> scratch fms, the operator context and the group matrix
1492!> structure
1493!> \param qs_env ...
1494!> \param kpoints ...
1495!> \param nspin number of spin components
1496!> \param matrix_ks_kp real-space Kohn-Sham matrices, dbcsr image set
1497!> \param matrix_s_kp real-space overlap matrix, dbcsr image set
1498!> \param fmwork returned global scratch fms
1499!> \param op_ctx returned operator context
1500!> \param matrix_struct returned AO layout of the k-point group
1501!> \param prepare prepare the operator snapshot right away. The
1502!> transfer rounds of the Lanczos driver assemble lazily through
1503!> the get path and skip the snapshot
1504!> \param prepared optional returned status of the snapshot preparation
1505! **************************************************************************************************
1506 SUBROUTINE kp_diag_setup(qs_env, kpoints, nspin, matrix_ks_kp, matrix_s_kp, fmwork, op_ctx, &
1507 matrix_struct, prepare, prepared)
1508 TYPE(qs_environment_type), POINTER :: qs_env
1509 TYPE(kpoint_type), POINTER :: kpoints
1510 INTEGER, INTENT(IN) :: nspin
1511 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
1512 TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
1513 TYPE(kpoint_operator_context_type) :: op_ctx
1514 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1515 LOGICAL, INTENT(IN) :: prepare
1516 LOGICAL, INTENT(OUT), OPTIONAL :: prepared
1517
1518 INTEGER :: i
1519 LOGICAL :: ks_ready, s_ready
1520 TYPE(cp_fm_pool_type), POINTER :: ao_ao_fm_pool, kp_ao_ao_fm_pool
1521 TYPE(qs_matrix_pools_type), POINTER :: kp_mpools, mpools
1522
1523 IF (PRESENT(prepared)) prepared = .false.
1524 CALL get_qs_env(qs_env, mpools=mpools)
1525 CALL mpools_get(mpools, ao_ao_fm_pool=ao_ao_fm_pool)
1526 ALLOCATE (fmwork(4))
1527 DO i = 1, 4
1528 CALL fm_pool_create_fm(ao_ao_fm_pool, fmwork(i))
1529 END DO
1530 CALL kpoint_operator_context_create(op_ctx, kpoints, fmwork(1), nspin)
1531 IF (prepare) THEN
1532 CALL kpoint_operator_prepare(op_ctx, matrix_ks_kp, kpoint_slot_ks, ks_ready)
1533 CALL kpoint_operator_prepare(op_ctx, matrix_s_kp, kpoint_slot_s, s_ready, &
1534 matrix_row=kpoint_spin_free)
1535 IF (PRESENT(prepared)) prepared = ks_ready .AND. s_ready
1536 END IF
1537 CALL get_kpoint_info(kpoints, mpools=kp_mpools)
1538 CALL mpools_get(kp_mpools, ao_ao_fm_pool=kp_ao_ao_fm_pool)
1539 matrix_struct => fm_pool_get_el_struct(kp_ao_ao_fm_pool)
1540 END SUBROUTINE kp_diag_setup
1541
1542! **************************************************************************************************
1543!> \brief Common finalization of the iterative k-point drivers: release
1544!> the operator context and return the global scratch fms
1545!> \param qs_env ...
1546!> \param fmwork global scratch fms of the driver
1547!> \param op_ctx operator context of the driver
1548! **************************************************************************************************
1549 SUBROUTINE kp_diag_finalize(qs_env, fmwork, op_ctx)
1550 TYPE(qs_environment_type), POINTER :: qs_env
1551 TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
1552 TYPE(kpoint_operator_context_type) :: op_ctx
1553
1554 INTEGER :: i
1555 TYPE(cp_fm_pool_type), POINTER :: ao_ao_fm_pool
1556 TYPE(qs_matrix_pools_type), POINTER :: mpools
1557
1558 CALL kpoint_operator_context_release(op_ctx)
1559 CALL get_qs_env(qs_env, mpools=mpools)
1560 CALL mpools_get(mpools, ao_ao_fm_pool=ao_ao_fm_pool)
1561 DO i = 1, 4
1562 CALL fm_pool_give_back_fm(ao_ao_fm_pool, fmwork(i))
1563 END DO
1564 DEALLOCATE (fmwork)
1565 END SUBROUTINE kp_diag_finalize
1566
1567! **************************************************************************************************
1568!> \brief Kpoint diagonalization routine
1569!> Transforms matrices to kpoint, distributes kpoint groups, performs
1570!> general diagonalization (no storgae of overlap decomposition), stores
1571!> MOs, calculates occupation numbers, calculates density matrices
1572!> in kpoint representation, transforms density matrices to real space
1573!> \param matrix_ks Kohn-sham matrices (RS indices, global)
1574!> \param matrix_s Overlap matrices (RS indices, global)
1575!> \param kpoints Kpoint environment
1576!> \param fmwork global-layout FM template in the first array entry
1577!> \par History
1578!> 08.2014 created [JGH]
1579! **************************************************************************************************
1580 SUBROUTINE diag_kp_basic(matrix_ks, matrix_s, kpoints, fmwork)
1581
1582 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
1583 TYPE(kpoint_type), POINTER :: kpoints
1584 TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
1585
1586 CHARACTER(len=*), PARAMETER :: routinen = 'diag_kp_basic'
1587
1588 INTEGER :: handle, ikp, ispin, kplocal, kpmax, nspin
1589 INTEGER, DIMENSION(2) :: kp_range
1590 INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1591 LOGICAL :: use_real_wfn
1592 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
1593 TYPE(cp_cfm_type) :: cksmat, cmos, csmat, cwork
1594 TYPE(cp_fm_pool_type), POINTER :: ao_ao_fm_pool
1595 TYPE(cp_fm_struct_type), POINTER :: matrix_struct, mo_struct
1596 TYPE(cp_fm_type) :: fmlocal, rksmat, rsmat
1597 TYPE(cp_fm_type), POINTER :: mo_coeff
1598 TYPE(kpoint_env_type), POINTER :: kp
1599 TYPE(kpoint_operator_context_type) :: op_ctx
1600 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1601 POINTER :: sab_nl
1602 TYPE(qs_matrix_pools_type), POINTER :: mpools
1603
1604 CALL timeset(routinen, handle)
1605
1606 NULLIFY (sab_nl)
1607 CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn, kp_range=kp_range, &
1608 kp_dist=kp_dist, sab_nl=sab_nl)
1609 cpassert(ASSOCIATED(sab_nl))
1610 kplocal = kp_range(2) - kp_range(1) + 1
1611 kpmax = maxval(kp_dist(2, :) - kp_dist(1, :) + 1)
1612
1613 ! fm pools to be used within a kpoint group
1614 CALL get_kpoint_info(kpoints, mpools=mpools)
1615 CALL mpools_get(mpools, ao_ao_fm_pool=ao_ao_fm_pool)
1616
1617 CALL fm_pool_create_fm(ao_ao_fm_pool, fmlocal)
1618 CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
1619
1620 IF (use_real_wfn) THEN
1621 CALL cp_fm_create(rksmat, matrix_struct)
1622 CALL cp_fm_create(rsmat, matrix_struct)
1623 ELSE
1624 CALL cp_cfm_create(cksmat, matrix_struct)
1625 CALL cp_cfm_create(csmat, matrix_struct)
1626 CALL cp_cfm_create(cwork, matrix_struct)
1627 kp => kpoints%kp_env(1)%kpoint_env
1628 mo_struct => kp%mos(1)%cmo_coeff%matrix_struct
1629 CALL cp_cfm_create(cmos, mo_struct)
1630 END IF
1631
1632 nspin = SIZE(matrix_ks, 1)
1633 CALL kpoint_operator_context_create(op_ctx, kpoints, fmwork(1), nspin)
1634 CALL kpoint_operator_prepare(op_ctx, matrix_ks, kpoint_slot_ks)
1635 CALL kpoint_operator_prepare(op_ctx, matrix_s, kpoint_slot_s, matrix_row=kpoint_spin_free)
1636 DO ikp = 1, kpmax
1637 DO ispin = 1, nspin
1638 CALL kpoint_operator_start(op_ctx, ikp, ispin, matrix_ks, kpoint_slot_ks, &
1639 use_grid=.false., matrix_row=ispin)
1640 CALL kpoint_operator_start(op_ctx, ikp, ispin, matrix_s, kpoint_slot_s, &
1641 use_grid=.false., matrix_row=kpoint_spin_free)
1642 END DO
1643 END DO
1644
1645 ! Finish communication then diagonalise in each group
1646 DO ikp = 1, kpmax
1647 DO ispin = 1, nspin
1648 IF (use_real_wfn) THEN
1649 CALL kpoint_operator_finish(op_ctx, ikp, ispin, kpoint_slot_ks, fm_re=rksmat)
1650 CALL kpoint_operator_finish(op_ctx, ikp, ispin, kpoint_slot_s, fm_re=rsmat)
1651 ELSE
1652 CALL kpoint_operator_finish(op_ctx, ikp, ispin, kpoint_slot_ks, cksmat)
1653 CALL kpoint_operator_finish(op_ctx, ikp, ispin, kpoint_slot_s, csmat)
1654 END IF
1655
1656 ! Each kpoint group has now information on a kpoint to be diagonalized
1657 ! General eigensolver Hermite or Symmetric
1658 IF (ikp > kplocal) cycle
1659 kp => kpoints%kp_env(ikp)%kpoint_env
1660 IF (use_real_wfn) THEN
1661 CALL get_mo_set(kp%mos(ispin), mo_coeff=mo_coeff, eigenvalues=eigenvalues)
1662 CALL cp_fm_geeig(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal)
1663 ELSE
1664 CALL get_mo_set(kp%mos(ispin), eigenvalues=eigenvalues)
1665 CALL cp_cfm_geeig(cksmat, csmat, cmos, eigenvalues, cwork)
1666 CALL cp_cfm_to_cfm(cmos, kp%mos(ispin)%cmo_coeff)
1667 END IF
1668 END DO
1669 END DO
1670
1671 CALL kpoint_operator_context_release(op_ctx)
1672
1673 IF (use_real_wfn) THEN
1674 CALL cp_fm_release(rksmat)
1675 CALL cp_fm_release(rsmat)
1676 ELSE
1677 CALL cp_cfm_release(cksmat)
1678 CALL cp_cfm_release(csmat)
1679 CALL cp_cfm_release(cwork)
1680 CALL cp_cfm_release(cmos)
1681 END IF
1682 CALL fm_pool_give_back_fm(ao_ao_fm_pool, fmlocal)
1683
1684 CALL timestop(handle)
1685
1686 END SUBROUTINE diag_kp_basic
1687
1688! **************************************************************************************************
1689!> \brief inner loop within MOS subspace, to refine occupation and density,
1690!> before next diagonalization of the Hamiltonian
1691!> \param qs_env ...
1692!> \param scf_env ...
1693!> \param subspace_env ...
1694!> \param mos ...
1695!> \param rho ...
1696!> \param ks_env ...
1697!> \param scf_section ...
1698!> \param scf_control ...
1699!> \par History
1700!> 09.2009 created [MI]
1701!> \note it is assumed that when diagonalization is used, also some mixing procedure is active
1702! **************************************************************************************************
1703 SUBROUTINE do_scf_diag_subspace(qs_env, scf_env, subspace_env, mos, rho, &
1704 ks_env, scf_section, scf_control)
1705
1706 TYPE(qs_environment_type), POINTER :: qs_env
1707 TYPE(qs_scf_env_type), POINTER :: scf_env
1708 TYPE(subspace_env_type), POINTER :: subspace_env
1709 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
1710 TYPE(qs_rho_type), POINTER :: rho
1711 TYPE(qs_ks_env_type), POINTER :: ks_env
1712 TYPE(section_vals_type), POINTER :: scf_section
1713 TYPE(scf_control_type), POINTER :: scf_control
1714
1715 CHARACTER(LEN=*), PARAMETER :: routinen = 'do_scf_diag_subspace'
1716 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1717
1718 INTEGER :: handle, i, iloop, ispin, nao, nmo, &
1719 nspin, output_unit
1720 LOGICAL :: converged
1721 REAL(dp) :: ene_diff, ene_old, iter_delta, max_val, &
1722 sum_band, sum_val, t1, t2
1723 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues, mo_occupations
1724 TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: eval_first, occ_first
1725 TYPE(cp_fm_type) :: work
1726 TYPE(cp_fm_type), POINTER :: c0, chc, evec, mo_coeff
1727 TYPE(cp_logger_type), POINTER :: logger
1728 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
1729 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
1730 TYPE(dft_control_type), POINTER :: dft_control
1731 TYPE(mp_para_env_type), POINTER :: para_env
1732 TYPE(qs_energy_type), POINTER :: energy
1733 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
1734
1735 CALL timeset(routinen, handle)
1736 NULLIFY (c0, chc, energy, evec, matrix_ks, mo_coeff, mo_eigenvalues, &
1737 mo_occupations, dft_control, rho_ao, rho_ao_kp)
1738
1739 logger => cp_get_default_logger()
1740 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIAG_SUB_SCF", &
1741 extension=".scfLog")
1742
1743 !Extra loop keeping mos unchanged and refining the subspace occupation
1744 nspin = SIZE(mos)
1745 CALL qs_rho_get(rho, rho_ao=rho_ao, rho_ao_kp=rho_ao_kp)
1746
1747 ALLOCATE (eval_first(nspin))
1748 ALLOCATE (occ_first(nspin))
1749 DO ispin = 1, nspin
1750 CALL get_mo_set(mo_set=mos(ispin), &
1751 nmo=nmo, &
1752 eigenvalues=mo_eigenvalues, &
1753 occupation_numbers=mo_occupations)
1754 ALLOCATE (eval_first(ispin)%array(nmo))
1755 ALLOCATE (occ_first(ispin)%array(nmo))
1756 eval_first(ispin)%array(1:nmo) = mo_eigenvalues(1:nmo)
1757 occ_first(ispin)%array(1:nmo) = mo_occupations(1:nmo)
1758 END DO
1759
1760 DO ispin = 1, nspin
1761 ! does not yet handle k-points
1762 CALL dbcsr_copy(subspace_env%p_matrix_store(ispin)%matrix, rho_ao(ispin)%matrix)
1763 CALL dbcsr_copy(rho_ao(ispin)%matrix, scf_env%p_mix_new(ispin, 1)%matrix)
1764 END DO
1765
1766 subspace_env%p_matrix_mix => scf_env%p_mix_new
1767
1768 NULLIFY (matrix_ks, energy, para_env, matrix_s)
1769 CALL get_qs_env(qs_env, &
1770 matrix_ks=matrix_ks, &
1771 energy=energy, &
1772 matrix_s=matrix_s, &
1773 para_env=para_env, &
1774 dft_control=dft_control)
1775
1776 ! mixing storage allocation
1777 IF (subspace_env%mixing_method >= gspace_mixing_nr) THEN
1778 CALL mixing_allocate(qs_env, subspace_env%mixing_method, scf_env%p_mix_new, &
1779 scf_env%p_delta, nspin, subspace_env%mixing_store)
1780 IF (dft_control%qs_control%gapw) THEN
1781 CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
1782 CALL mixing_init(subspace_env%mixing_method, rho, subspace_env%mixing_store, &
1783 para_env, rho_atom=rho_atom)
1784 ELSE IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
1785 CALL charge_mixing_init(subspace_env%mixing_store)
1786 ELSE IF (dft_control%qs_control%semi_empirical) THEN
1787 cpabort('SE Code not possible')
1788 ELSE
1789 CALL mixing_init(subspace_env%mixing_method, rho, subspace_env%mixing_store, para_env)
1790 END IF
1791 END IF
1792
1793 ene_old = 0.0_dp
1794 ene_diff = 0.0_dp
1795 IF (output_unit > 0) THEN
1796 WRITE (output_unit, "(/T19,A)") '<<<<<<<<< SUBSPACE ROTATION <<<<<<<<<<'
1797 WRITE (output_unit, "(T4,A,T13,A,T21,A,T38,A,T51,A,T65,A/,T4,A)") &
1798 "In-step", "Time", "Convergence", "Band ene.", "Total ene.", "Energy diff.", repeat("-", 74)
1799 END IF
1800
1801 ! recalculate density matrix here
1802
1803 ! update of density
1804 CALL qs_rho_update_rho(rho, qs_env=qs_env)
1805
1806 DO iloop = 1, subspace_env%max_iter
1807 t1 = m_walltime()
1808 converged = .false.
1809 ene_old = energy%total
1810
1811 CALL qs_ks_did_change(ks_env, rho_changed=.true.)
1812 CALL qs_ks_update_qs_env(qs_env, calculate_forces=.false., &
1813 just_energy=.false., print_active=.false.)
1814
1815 max_val = 0.0_dp
1816 sum_val = 0.0_dp
1817 sum_band = 0.0_dp
1818 DO ispin = 1, SIZE(matrix_ks)
1819 CALL get_mo_set(mo_set=mos(ispin), &
1820 nao=nao, &
1821 nmo=nmo, &
1822 eigenvalues=mo_eigenvalues, &
1823 occupation_numbers=mo_occupations, &
1824 mo_coeff=mo_coeff)
1825
1826 !compute C'HC
1827 chc => subspace_env%chc_mat(ispin)
1828 evec => subspace_env%c_vec(ispin)
1829 c0 => subspace_env%c0(ispin)
1830 CALL cp_fm_to_fm(mo_coeff, c0)
1831 CALL cp_fm_create(work, c0%matrix_struct)
1832 CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, c0, work, nmo)
1833 CALL parallel_gemm('T', 'N', nmo, nmo, nao, rone, c0, work, rzero, chc)
1834 CALL cp_fm_release(work)
1835 !diagonalize C'HC
1836 CALL choose_eigv_solver(chc, evec, mo_eigenvalues)
1837
1838 !rotate the mos by the eigenvectors of C'HC
1839 CALL parallel_gemm('N', 'N', nao, nmo, nmo, rone, c0, evec, rzero, mo_coeff)
1840
1841 IF (.NOT. scf_control%gce%do_gce) THEN
1842 CALL set_mo_occupation(mo_set=mos(ispin), &
1843 smear=scf_control%smear)
1844 ELSE
1845 CALL set_mo_occupation(mo_set=mos(ispin), &
1846 smear=scf_control%smear, &
1847 gce=scf_control%gce)
1848 END IF
1849
1850 ! does not yet handle k-points
1851 CALL calculate_density_matrix(mos(ispin), &
1852 subspace_env%p_matrix_mix(ispin, 1)%matrix)
1853
1854 DO i = 1, nmo
1855 sum_band = sum_band + mo_eigenvalues(i)*mo_occupations(i)
1856 END DO
1857
1858 !check for self consistency
1859 END DO
1860
1861 IF (subspace_env%mixing_method == direct_mixing_nr) THEN
1862 CALL scf_env_density_mixing(subspace_env%p_matrix_mix, &
1863 scf_env%mixing_store, rho_ao_kp, para_env, iter_delta, iloop)
1864 ELSE
1865 CALL self_consistency_check(rho_ao_kp, scf_env%p_delta, para_env, &
1866 subspace_env%p_matrix_mix, delta=iter_delta)
1867 END IF
1868
1869 DO ispin = 1, nspin
1870 ! does not yet handle k-points
1871 CALL dbcsr_copy(rho_ao(ispin)%matrix, subspace_env%p_matrix_mix(ispin, 1)%matrix)
1872 END DO
1873 ! update of density
1874 CALL qs_rho_update_rho(rho, qs_env=qs_env)
1875 ! Mixing in reciprocal space
1876 IF (subspace_env%mixing_method >= gspace_mixing_nr) THEN
1877 CALL gspace_mixing(qs_env, scf_env%mixing_method, subspace_env%mixing_store, &
1878 rho, para_env, scf_env%iter_count)
1879 END IF
1880
1881 ene_diff = energy%total - ene_old
1882 converged = (abs(ene_diff) < subspace_env%eps_ene .AND. &
1883 iter_delta < subspace_env%eps_adapt*scf_env%iter_delta)
1884 t2 = m_walltime()
1885 IF (output_unit > 0) THEN
1886 WRITE (output_unit, "(T4,I5,T11,F8.3,T18,E14.4,T34,F12.5,T46,F16.8,T62,E14.4)") &
1887 iloop, t2 - t1, iter_delta, sum_band, energy%total, ene_diff
1888 CALL m_flush(output_unit)
1889 END IF
1890 IF (converged) THEN
1891 IF (output_unit > 0) WRITE (output_unit, "(T10,A,I6,A,/)") &
1892 " Reached convergence in ", iloop, " iterations "
1893 EXIT
1894 END IF
1895
1896 END DO ! iloop
1897
1898 NULLIFY (subspace_env%p_matrix_mix)
1899 DO ispin = 1, nspin
1900 ! does not yet handle k-points
1901 CALL dbcsr_copy(scf_env%p_mix_new(ispin, 1)%matrix, rho_ao(ispin)%matrix)
1902 CALL dbcsr_copy(rho_ao(ispin)%matrix, subspace_env%p_matrix_store(ispin)%matrix)
1903
1904 DEALLOCATE (eval_first(ispin)%array, occ_first(ispin)%array)
1905 END DO
1906 DEALLOCATE (eval_first, occ_first)
1907
1908 CALL timestop(handle)
1909
1910 END SUBROUTINE do_scf_diag_subspace
1911
1912! **************************************************************************************************
1913!> \brief ...
1914!> \param subspace_env ...
1915!> \param qs_env ...
1916!> \param mos ...
1917! **************************************************************************************************
1918 SUBROUTINE diag_subspace_allocate(subspace_env, qs_env, mos)
1919
1920 TYPE(subspace_env_type), POINTER :: subspace_env
1921 TYPE(qs_environment_type), POINTER :: qs_env
1922 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1923
1924 CHARACTER(LEN=*), PARAMETER :: routinen = 'diag_subspace_allocate'
1925
1926 INTEGER :: handle, i, ispin, nmo, nspin
1927 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
1928 TYPE(cp_fm_type), POINTER :: mo_coeff
1929 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1930 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1931 POINTER :: sab_orb
1932
1933 CALL timeset(routinen, handle)
1934
1935 NULLIFY (sab_orb, matrix_s)
1936 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb, &
1937 matrix_s=matrix_s)
1938
1939 nspin = SIZE(mos)
1940! *** allocate p_atrix_store ***
1941 IF (.NOT. ASSOCIATED(subspace_env%p_matrix_store)) THEN
1942 CALL dbcsr_allocate_matrix_set(subspace_env%p_matrix_store, nspin)
1943
1944 DO i = 1, nspin
1945 ALLOCATE (subspace_env%p_matrix_store(i)%matrix)
1946 CALL dbcsr_create(matrix=subspace_env%p_matrix_store(i)%matrix, template=matrix_s(1)%matrix, &
1947 name="DENSITY_STORE", matrix_type=dbcsr_type_symmetric)
1948 CALL cp_dbcsr_alloc_block_from_nbl(subspace_env%p_matrix_store(i)%matrix, &
1949 sab_orb)
1950 CALL dbcsr_set(subspace_env%p_matrix_store(i)%matrix, 0.0_dp)
1951 END DO
1952
1953 END IF
1954
1955 ALLOCATE (subspace_env%chc_mat(nspin))
1956 ALLOCATE (subspace_env%c_vec(nspin))
1957 ALLOCATE (subspace_env%c0(nspin))
1958
1959 DO ispin = 1, nspin
1960 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1961 CALL cp_fm_create(subspace_env%c0(ispin), mo_coeff%matrix_struct)
1962 NULLIFY (fm_struct_tmp)
1963 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo, ncol_global=nmo, &
1964 para_env=mo_coeff%matrix_struct%para_env, &
1965 context=mo_coeff%matrix_struct%context)
1966 CALL cp_fm_create(subspace_env%chc_mat(ispin), fm_struct_tmp, "chc")
1967 CALL cp_fm_create(subspace_env%c_vec(ispin), fm_struct_tmp, "vec")
1968 CALL cp_fm_struct_release(fm_struct_tmp)
1969 END DO
1970
1971 CALL timestop(handle)
1972
1973 END SUBROUTINE diag_subspace_allocate
1974
1975! **************************************************************************************************
1976!> \brief the inner loop of scf, specific to diagonalization without S matrix
1977!> basically, in goes the ks matrix out goes a new p matrix
1978!> \param scf_env ...
1979!> \param mos ...
1980!> \param matrix_ks ...
1981!> \param scf_control ...
1982!> \param scf_section ...
1983!> \param diis_step ...
1984!> \par History
1985!> 03.2006 created [Joost VandeVondele]
1986! **************************************************************************************************
1987 SUBROUTINE do_special_diag(scf_env, mos, matrix_ks, scf_control, &
1988 scf_section, diis_step)
1989
1990 TYPE(qs_scf_env_type), POINTER :: scf_env
1991 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
1992 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1993 TYPE(scf_control_type), POINTER :: scf_control
1994 TYPE(section_vals_type), POINTER :: scf_section
1995 LOGICAL, INTENT(INOUT) :: diis_step
1996
1997 INTEGER :: ispin, nspin
1998 LOGICAL :: do_level_shift, use_jacobi
1999 REAL(kind=dp) :: diis_error
2000
2001 nspin = SIZE(matrix_ks)
2002
2003 DO ispin = 1, nspin
2004 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, scf_env%scf_work1(ispin), &
2005 scf_env%fm_copy_plan(ispin), check_plan=.true.)
2006 END DO
2007 IF (scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis) THEN
2008 CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
2009 scf_env%scf_work2, scf_env%iter_delta, diis_error, diis_step, &
2010 scf_control%eps_diis, scf_control%nmixing, &
2011 scf_section=scf_section)
2012 ELSE
2013 diis_step = .false.
2014 END IF
2015
2016 IF ((scf_env%iter_count > 1) .AND. (scf_env%iter_delta < scf_control%diagonalization%eps_jacobi)) THEN
2017 use_jacobi = .true.
2018 ELSE
2019 use_jacobi = .false.
2020 END IF
2021
2022 do_level_shift = ((scf_control%level_shift /= 0.0_dp) .AND. &
2023 ((scf_control%density_guess == core_guess) .OR. (scf_env%iter_count > 1)))
2024 IF (diis_step) THEN
2025 scf_env%iter_param = diis_error
2026 IF (use_jacobi) THEN
2027 scf_env%iter_method = "DIIS/Jacobi"
2028 ELSE
2029 scf_env%iter_method = "DIIS/Diag."
2030 END IF
2031 ELSE
2032 IF (scf_env%mixing_method == 1) THEN
2033 scf_env%iter_param = scf_env%p_mix_alpha
2034 IF (use_jacobi) THEN
2035 scf_env%iter_method = "P_Mix/Jacobi"
2036 ELSE
2037 scf_env%iter_method = "P_Mix/Diag."
2038 END IF
2039 ELSE IF (scf_env%mixing_method > 1) THEN
2040 scf_env%iter_param = scf_env%mixing_store%alpha
2041 IF (use_jacobi) THEN
2042 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Jacobi"
2043 ELSE
2044 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Diag."
2045 END IF
2046 END IF
2047 END IF
2048 scf_env%iter_delta = 0.0_dp
2049
2050 DO ispin = 1, nspin
2051 CALL eigensolver_simple(matrix_ks=scf_env%scf_work1(ispin), &
2052 mo_set=mos(ispin), &
2053 work=scf_env%scf_work2, &
2054 do_level_shift=do_level_shift, &
2055 level_shift=scf_control%level_shift, &
2056 use_jacobi=use_jacobi, &
2057 jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
2058 END DO
2059
2060 IF (.NOT. scf_control%gce%do_gce) THEN
2061 CALL set_mo_occupation(mo_array=mos, &
2062 smear=scf_control%smear)
2063 ELSE
2064 CALL set_mo_occupation(mo_array=mos, &
2065 smear=scf_control%smear, &
2066 gce=scf_control%gce)
2067 END IF
2068
2069 DO ispin = 1, nspin
2070 ! does not yet handle k-points
2071 CALL calculate_density_matrix(mos(ispin), &
2072 scf_env%p_mix_new(ispin, 1)%matrix)
2073 END DO
2074
2075 END SUBROUTINE do_special_diag
2076
2077! **************************************************************************************************
2078!> \brief the inner loop of scf, specific to iterative diagonalization using OT
2079!> with S matrix; basically, in goes the ks matrix out goes a new p matrix
2080!> \param scf_env ...
2081!> \param mos ...
2082!> \param matrix_ks ...
2083!> \param matrix_s ...
2084!> \param scf_control ...
2085!> \param scf_section ...
2086!> \param diis_step ...
2087!> \par History
2088!> 10.2008 created [JGH]
2089! **************************************************************************************************
2090 SUBROUTINE do_ot_diag(scf_env, mos, matrix_ks, matrix_s, &
2091 scf_control, scf_section, diis_step)
2092
2093 TYPE(qs_scf_env_type), POINTER :: scf_env
2094 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
2095 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
2096 TYPE(scf_control_type), POINTER :: scf_control
2097 TYPE(section_vals_type), POINTER :: scf_section
2098 LOGICAL, INTENT(INOUT) :: diis_step
2099
2100 INTEGER :: homo, ispin, nmo, nspin
2101 REAL(kind=dp) :: diis_error, eps_iter
2102 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
2103 TYPE(cp_fm_type), POINTER :: mo_coeff
2104
2105 NULLIFY (eigenvalues)
2106
2107 nspin = SIZE(matrix_ks)
2108
2109 DO ispin = 1, nspin
2110 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, &
2111 scf_env%scf_work1(ispin))
2112 END DO
2113
2114 IF ((scf_env%iter_count > 1) .AND. (.NOT. scf_env%skip_diis)) THEN
2115 CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
2116 scf_env%scf_work2, scf_env%iter_delta, diis_error, diis_step, &
2117 scf_control%eps_diis, scf_control%nmixing, &
2118 s_matrix=matrix_s, &
2119 scf_section=scf_section)
2120 ELSE
2121 diis_step = .false.
2122 END IF
2123
2124 eps_iter = scf_control%diagonalization%eps_iter
2125 IF (diis_step) THEN
2126 scf_env%iter_param = diis_error
2127 scf_env%iter_method = "DIIS/OTdiag"
2128 DO ispin = 1, nspin
2129 CALL copy_fm_to_dbcsr(scf_env%scf_work1(ispin), &
2130 matrix_ks(ispin)%matrix, keep_sparsity=.true.)
2131 END DO
2132 eps_iter = max(eps_iter, scf_control%diagonalization%eps_adapt*diis_error)
2133 ELSE
2134 IF (scf_env%mixing_method == 1) THEN
2135 scf_env%iter_param = scf_env%p_mix_alpha
2136 scf_env%iter_method = "P_Mix/OTdiag."
2137 ELSE IF (scf_env%mixing_method > 1) THEN
2138 scf_env%iter_param = scf_env%mixing_store%alpha
2139 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/OTdiag."
2140 END IF
2141 END IF
2142
2143 scf_env%iter_delta = 0.0_dp
2144
2145 DO ispin = 1, nspin
2146 CALL get_mo_set(mos(ispin), &
2147 mo_coeff=mo_coeff, &
2148 eigenvalues=eigenvalues, &
2149 nmo=nmo, &
2150 homo=homo)
2151 CALL ot_eigensolver(matrix_h=matrix_ks(ispin)%matrix, &
2152 matrix_s=matrix_s(1)%matrix, &
2153 matrix_c_fm=mo_coeff, &
2154 preconditioner=scf_env%ot_preconditioner( &
2155 min(ispin, SIZE(scf_env%ot_preconditioner)))%preconditioner, &
2156 eps_gradient=eps_iter, &
2157 iter_max=scf_control%diagonalization%max_iter, &
2158 silent=.true., &
2159 ot_settings=scf_control%diagonalization%ot_settings)
2160 CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin)%matrix, &
2161 evals_arg=eigenvalues, &
2162 do_rotation=.true.)
2163 CALL copy_fm_to_dbcsr(mos(ispin)%mo_coeff, &
2164 mos(ispin)%mo_coeff_b)
2165 !fm->dbcsr
2166 END DO
2167
2168 CALL set_mo_occupation(mo_array=mos, &
2169 smear=scf_control%smear)
2170
2171 DO ispin = 1, nspin
2172 ! does not yet handle k-points
2173 CALL calculate_density_matrix(mos(ispin), &
2174 scf_env%p_mix_new(ispin, 1)%matrix)
2175 END DO
2176
2177 END SUBROUTINE do_ot_diag
2178
2179! **************************************************************************************************
2180!> \brief fixed-H OT diagonalization for real or complex K-point channels
2181!> \param matrix_ks real-space Kohn-Sham matrices
2182!> \param matrix_s real-space overlap matrices
2183!> \param matrix_t ...
2184!> \param kpoints K-point environment
2185!> \param scf_env SCF environment
2186!> \param scf_control SCF controls
2187!> \param diis_step set false; density mixing remains the outer accelerator
2188!> \param added_mos_auto_grow request a larger adaptive virtual-space buffer
2189! **************************************************************************************************
2190 SUBROUTINE do_ot_diag_kp(matrix_ks, matrix_s, matrix_t, kpoints, scf_env, scf_control, &
2191 diis_step, added_mos_auto_grow)
2192 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s, matrix_t
2193 TYPE(kpoint_type), POINTER :: kpoints
2194 TYPE(qs_scf_env_type), POINTER :: scf_env
2195 TYPE(scf_control_type), POINTER :: scf_control
2196 LOGICAL, INTENT(INOUT) :: diis_step
2197 LOGICAL, INTENT(OUT) :: added_mos_auto_grow
2198
2199 INTEGER :: ikp, ispin, ks_spin, nspin
2200 LOGICAL :: store_kinetic, use_real_wfn
2201 REAL(kind=dp) :: eps_iter
2202 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
2203 TYPE(cp_cfm_type) :: cmo_work
2204 TYPE(cp_fm_type) :: coeff_im_work
2205 TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
2206 TYPE(cp_fm_type), POINTER :: coeff_re
2207 TYPE(dbcsr_type), POINTER :: matrix_h_im, matrix_h_re, matrix_s_im, &
2208 matrix_s_re, matrix_t_im, matrix_t_re
2209 TYPE(kpoint_env_type), POINTER :: kp
2210 TYPE(kpoint_operator_context_type) :: op_ctx
2211 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2212 POINTER :: sab_nl
2213
2214 NULLIFY (coeff_re, eigenvalues, fmwork, kp, matrix_h_im, matrix_h_re, &
2215 matrix_s_im, matrix_s_re, matrix_t_im, matrix_t_re, sab_nl)
2216 added_mos_auto_grow = .false.
2217 diis_step = .false.
2218 eps_iter = scf_control%diagonalization%eps_iter
2219 nspin = SIZE(matrix_ks, 1)
2220 store_kinetic = scf_control%diagonalization%ot_settings%preconditioner_type == &
2221 ot_precond_full_kinetic
2222 IF (store_kinetic) THEN
2223 cpassert(ASSOCIATED(matrix_t))
2224 END IF
2225 CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn, sab_nl=sab_nl)
2226 cpassert(ASSOCIATED(sab_nl))
2227
2228 ! A full diagonalization supplies an orthonormal start after initialization or an
2229 ! ADDED_MOS AUTO resize. Subsequent SCF steps only redistribute H(k) and S(k).
2230 IF (scf_env%iter_count <= 1) THEN
2231 CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, &
2232 .false., diis_step)
2233 END IF
2234 CALL kpoint_operator_context_create(op_ctx, kpoints, scf_env%scf_work1(1), nspin)
2235 IF (store_kinetic) THEN
2236 CALL kpoint_operator_store(op_ctx, matrix_ks, matrix_s, matrix_t)
2237 ELSE
2238 CALL kpoint_operator_store(op_ctx, matrix_ks, matrix_s)
2239 END IF
2240 CALL kpoint_operator_context_release(op_ctx)
2241
2242 scf_env%iter_delta = 0.0_dp
2243 DO ikp = 1, SIZE(kpoints%kp_env)
2244 kp => kpoints%kp_env(ikp)%kpoint_env
2245 cpassert(ASSOCIATED(kp%ot_hmat) .AND. ASSOCIATED(kp%ot_smat))
2246 IF (SIZE(kp%ot_smat) >= 2) THEN
2247 CALL kpoint_operator_get_local(matrix_s, kpoints, kp, 1, &
2248 kp%ot_smat(1), kp%ot_smat(2), &
2249 matrix_s_re, matrix_s_im)
2250 ELSE
2251 CALL kpoint_operator_get_local(matrix_s, kpoints, kp, 1, &
2252 kp%ot_smat(1), kp%ot_smat(1), &
2253 matrix_s_re, matrix_s_im)
2254 CALL dbcsr_set(matrix_s_im, 0.0_dp)
2255 END IF
2256 IF (store_kinetic) THEN
2257 cpassert(ASSOCIATED(kp%ot_tmat))
2258 IF (SIZE(kp%ot_tmat) >= 2) THEN
2259 CALL kpoint_operator_get_local(matrix_t, kpoints, kp, 1, &
2260 kp%ot_tmat(1), kp%ot_tmat(2), &
2261 matrix_t_re, matrix_t_im)
2262 ELSE
2263 CALL kpoint_operator_get_local(matrix_t, kpoints, kp, 1, &
2264 kp%ot_tmat(1), kp%ot_tmat(1), &
2265 matrix_t_re, matrix_t_im)
2266 CALL dbcsr_set(matrix_t_im, 0.0_dp)
2267 END IF
2268 END IF
2269 DO ispin = 1, nspin
2270 ks_spin = min(ispin, SIZE(kp%ot_hmat, 2))
2271 IF (SIZE(kp%ot_hmat, 1) >= 2) THEN
2272 CALL kpoint_operator_get_local(matrix_ks, kpoints, kp, ks_spin, &
2273 kp%ot_hmat(1, ks_spin), kp%ot_hmat(2, ks_spin), &
2274 matrix_h_re, matrix_h_im)
2275 ELSE
2276 CALL kpoint_operator_get_local(matrix_ks, kpoints, kp, ks_spin, &
2277 kp%ot_hmat(1, ks_spin), kp%ot_hmat(1, ks_spin), &
2278 matrix_h_re, matrix_h_im)
2279 CALL dbcsr_set(matrix_h_im, 0.0_dp)
2280 END IF
2281 CALL get_mo_set(kp%mos(ispin), eigenvalues=eigenvalues)
2282 IF (use_real_wfn) THEN
2283 CALL get_mo_set(kp%mos(ispin), mo_coeff=coeff_re)
2284 CALL cp_fm_create(coeff_im_work, coeff_re%matrix_struct, &
2285 name='real k-point OT imaginary work')
2286 CALL cp_cfm_create(cmo_work, coeff_re%matrix_struct)
2287 CALL cp_cfm_scale_and_add_fm(z_zero, cmo_work, z_one, coeff_re)
2288 CALL ot_eigensolver_complex(matrix_h=matrix_h_re, matrix_h_im=matrix_h_im, &
2289 matrix_s=matrix_s_re, matrix_s_im=matrix_s_im, &
2290 matrix_c_cfm=cmo_work, &
2291 eps_gradient=eps_iter, &
2292 iter_max=scf_control%diagonalization%max_iter, &
2293 eigenvalues=eigenvalues, silent=.true., &
2294 ot_settings=scf_control%diagonalization%ot_settings, &
2295 matrix_t=matrix_t_re, matrix_t_im=matrix_t_im, &
2296 mo_set=kp%mos(ispin))
2297 CALL cp_cfm_to_fm(cmo_work, coeff_re, coeff_im_work)
2298 CALL cp_cfm_release(cmo_work)
2299 CALL cp_fm_release(coeff_im_work)
2300 ELSE
2301 CALL ot_eigensolver_complex(matrix_h=matrix_h_re, matrix_h_im=matrix_h_im, &
2302 matrix_s=matrix_s_re, matrix_s_im=matrix_s_im, &
2303 matrix_c_cfm=kp%mos(ispin)%cmo_coeff, &
2304 eps_gradient=eps_iter, &
2305 iter_max=scf_control%diagonalization%max_iter, &
2306 eigenvalues=eigenvalues, silent=.true., &
2307 ot_settings=scf_control%diagonalization%ot_settings, &
2308 matrix_t=matrix_t_re, matrix_t_im=matrix_t_im, &
2309 mo_set=kp%mos(ispin))
2310 END IF
2311 CALL dbcsr_release_p(matrix_h_re)
2312 CALL dbcsr_release_p(matrix_h_im)
2313 END DO
2314 CALL dbcsr_release_p(matrix_s_re)
2315 CALL dbcsr_release_p(matrix_s_im)
2316 IF (ASSOCIATED(matrix_t_re)) CALL dbcsr_release_p(matrix_t_re)
2317 IF (ASSOCIATED(matrix_t_im)) CALL dbcsr_release_p(matrix_t_im)
2318 END DO
2319
2320 fmwork => scf_env%scf_work1
2321 CALL kp_diag_exit_chain(kpoints, scf_env, scf_control, matrix_s, sab_nl, fmwork, added_mos_auto_grow)
2322 IF (added_mos_auto_grow) RETURN
2323 IF (scf_env%mixing_method == 0) THEN
2324 scf_env%iter_method = 'NoMix/OTdiag.'
2325 ELSE IF (scf_env%mixing_method == 1) THEN
2326 scf_env%iter_param = scf_env%p_mix_alpha
2327 scf_env%iter_method = 'P_Mix/OTdiag.'
2328 ELSE
2329 scf_env%iter_param = scf_env%mixing_store%alpha
2330 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//'/OTdiag.'
2331 END IF
2332 END SUBROUTINE do_ot_diag_kp
2333
2334! **************************************************************************************************
2335!> \brief Solve a set restricted open Kohn-Sham (ROKS) equations based on the
2336!> alpha and beta Kohn-Sham matrices from unrestricted Kohn-Sham.
2337!> \param scf_env ...
2338!> \param mos ...
2339!> \param matrix_ks ...
2340!> \param matrix_s ...
2341!> \param scf_control ...
2342!> \param scf_section ...
2343!> \param diis_step ...
2344!> \param orthogonal_basis ...
2345!> \par History
2346!> 04.2006 created [MK]
2347!> Revised (01.05.06,MK)
2348!> \note
2349!> this is only a high-spin ROKS.
2350! **************************************************************************************************
2351 SUBROUTINE do_roks_diag(scf_env, mos, matrix_ks, matrix_s, &
2352 scf_control, scf_section, diis_step, &
2353 orthogonal_basis)
2354
2355 ! Literature: - C. C. J. Roothaan, Rev. Mod. Phys. 32, 179 (1960)
2356 ! - M. F. Guest and V. R. Saunders, Mol. Phys. 28(3), 819 (1974)
2357 ! - M. Filatov and S. Shaik, Chem. Phys. Lett. 288, 689 (1998)
2358
2359 TYPE(qs_scf_env_type), POINTER :: scf_env
2360 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
2361 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
2362 TYPE(scf_control_type), POINTER :: scf_control
2363 TYPE(section_vals_type), POINTER :: scf_section
2364 LOGICAL, INTENT(INOUT) :: diis_step
2365 LOGICAL, INTENT(IN) :: orthogonal_basis
2366
2367 CHARACTER(LEN=*), PARAMETER :: routinen = 'do_roks_diag'
2368
2369 INTEGER :: handle, homoa, homob, imo, nalpha, nao, &
2370 nbeta, nmo
2371 REAL(kind=dp) :: diis_error, level_shift_loc
2372 REAL(kind=dp), DIMENSION(:), POINTER :: eiga, eigb, occa, occb
2373 TYPE(cp_fm_type), POINTER :: ksa, ksb, mo2ao, moa, mob, ortho, work
2374
2375! -------------------------------------------------------------------------
2376
2377 CALL timeset(routinen, handle)
2378
2379 IF (scf_env%cholesky_method == cholesky_inverse) THEN
2380 ortho => scf_env%ortho_m1
2381 ELSE
2382 ortho => scf_env%ortho
2383 END IF
2384 work => scf_env%scf_work2
2385
2386 ksa => scf_env%scf_work1(1)
2387 ksb => scf_env%scf_work1(2)
2388
2389 CALL copy_dbcsr_to_fm(matrix_ks(1)%matrix, ksa, scf_env%fm_copy_plan(1), check_plan=.true.)
2390 CALL copy_dbcsr_to_fm(matrix_ks(2)%matrix, ksb, scf_env%fm_copy_plan(2), check_plan=.true.)
2391
2392 ! Get MO information
2393
2394 CALL get_mo_set(mo_set=mos(1), &
2395 nao=nao, &
2396 nmo=nmo, &
2397 nelectron=nalpha, &
2398 homo=homoa, &
2399 eigenvalues=eiga, &
2400 occupation_numbers=occa, &
2401 mo_coeff=moa)
2402
2403 CALL get_mo_set(mo_set=mos(2), &
2404 nelectron=nbeta, &
2405 homo=homob, &
2406 eigenvalues=eigb, &
2407 occupation_numbers=occb, &
2408 mo_coeff=mob)
2409
2410 ! Define the amount of level-shifting
2411
2412 IF ((scf_control%level_shift /= 0.0_dp) .AND. &
2413 ((scf_control%density_guess == core_guess) .OR. &
2414 (scf_control%density_guess == restart_guess) .OR. &
2415 (scf_env%iter_count > 1))) THEN
2416 level_shift_loc = scf_control%level_shift
2417 ELSE
2418 level_shift_loc = 0.0_dp
2419 END IF
2420
2421 IF ((scf_env%iter_count > 1) .OR. &
2422 (scf_control%density_guess == core_guess) .OR. &
2423 (scf_control%density_guess == restart_guess)) THEN
2424
2425 ! Transform the spin unrestricted alpha and beta Kohn-Sham matrices
2426 ! from AO basis to MO basis: K(MO) = C(T)*K(AO)*C
2427
2428 CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, moa, 0.0_dp, work)
2429 CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, moa, work, 0.0_dp, ksa)
2430
2431 CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksb, moa, 0.0_dp, work)
2432 CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, moa, work, 0.0_dp, ksb)
2433
2434 ! Combine the spin unrestricted alpha and beta Kohn-Sham matrices
2435 ! in the MO basis
2436
2437 IF (scf_control%roks_scheme == general_roks) THEN
2438 CALL combine_ks_matrices(ksa, ksb, occa, occb, scf_control%roks_f, &
2439 nalpha, nbeta)
2440 ELSE IF (scf_control%roks_scheme == high_spin_roks) THEN
2441 CALL combine_ks_matrices(ksa, ksb, occa, occb, scf_control%roks_parameter)
2442 ELSE
2443 cpabort("Unknown ROKS scheme requested")
2444 END IF
2445
2446 ! Back-transform the restricted open Kohn-Sham matrix from MO basis
2447 ! to AO basis
2448
2449 IF (orthogonal_basis) THEN
2450 ! Q = C
2451 mo2ao => moa
2452 ELSE
2453 ! Q = S*C
2454 mo2ao => mob
2455!MK CALL copy_sm_to_fm(matrix_s(1)%matrix,work)
2456!MK CALL cp_fm_symm("L", "U",nao, nao, 1.0_dp, work, moa, 0.0_dp, mo2ao)
2457 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, moa, mo2ao, nao)
2458 END IF
2459
2460 ! K(AO) = Q*K(MO)*Q(T)
2461
2462 CALL parallel_gemm("N", "T", nao, nao, nao, 1.0_dp, ksa, mo2ao, 0.0_dp, work)
2463 CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, mo2ao, work, 0.0_dp, ksa)
2464
2465 ELSE
2466
2467 ! No transformation matrix available, yet. The closed shell part,
2468 ! i.e. the beta Kohn-Sham matrix in AO basis, is taken.
2469 ! There might be better choices, anyhow.
2470
2471 CALL cp_fm_to_fm(ksb, ksa)
2472
2473 END IF
2474
2475 ! Update DIIS buffer and possibly perform DIIS extrapolation step
2476
2477 IF (scf_env%iter_count > 1) THEN
2478 IF (orthogonal_basis) THEN
2479 CALL qs_diis_b_step(diis_buffer=scf_env%scf_diis_buffer, &
2480 mo_array=mos, &
2481 kc=scf_env%scf_work1, &
2482 sc=work, &
2483 delta=scf_env%iter_delta, &
2484 error_max=diis_error, &
2485 diis_step=diis_step, &
2486 eps_diis=scf_control%eps_diis, &
2487 scf_section=scf_section, &
2488 roks=.true.)
2489 cpassert(scf_env%iter_delta == scf_env%iter_delta)
2490 ELSE
2491 CALL qs_diis_b_step(diis_buffer=scf_env%scf_diis_buffer, &
2492 mo_array=mos, &
2493 kc=scf_env%scf_work1, &
2494 sc=work, &
2495 delta=scf_env%iter_delta, &
2496 error_max=diis_error, &
2497 diis_step=diis_step, &
2498 eps_diis=scf_control%eps_diis, &
2499 scf_section=scf_section, &
2500 s_matrix=matrix_s, &
2501 roks=.true.)
2502 END IF
2503 END IF
2504
2505 IF (diis_step) THEN
2506 scf_env%iter_param = diis_error
2507 scf_env%iter_method = "DIIS/Diag."
2508 ELSE
2509 IF (scf_env%mixing_method == 1) THEN
2510 scf_env%iter_param = scf_env%p_mix_alpha
2511 scf_env%iter_method = "P_Mix/Diag."
2512 ELSE IF (scf_env%mixing_method > 1) THEN
2513 scf_env%iter_param = scf_env%mixing_store%alpha
2514 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Diag."
2515 END IF
2516 END IF
2517
2518 scf_env%iter_delta = 0.0_dp
2519
2520 IF (level_shift_loc /= 0.0_dp) THEN
2521
2522 ! Transform the current Kohn-Sham matrix from AO to MO basis
2523 ! for level-shifting using the current MO set
2524
2525 CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, moa, 0.0_dp, work)
2526 CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, moa, work, 0.0_dp, ksa)
2527
2528 ! Apply level-shifting using 50:50 split of the shift (could be relaxed)
2529
2530 DO imo = homob + 1, homoa
2531 CALL cp_fm_add_to_element(ksa, imo, imo, 0.5_dp*level_shift_loc)
2532 END DO
2533 DO imo = homoa + 1, nmo
2534 CALL cp_fm_add_to_element(ksa, imo, imo, level_shift_loc)
2535 END DO
2536
2537 ELSE IF (.NOT. orthogonal_basis) THEN
2538
2539 ! Transform the current Kohn-Sham matrix to an orthogonal basis
2540 SELECT CASE (scf_env%cholesky_method)
2541 CASE (cholesky_reduce)
2542 CALL cp_fm_cholesky_reduce(ksa, ortho)
2543 CASE (cholesky_restore)
2544 CALL cp_fm_uplo_to_full(ksa, work)
2545 CALL cp_fm_cholesky_restore(ksa, nao, ortho, work, &
2546 "SOLVE", pos="RIGHT")
2547 CALL cp_fm_cholesky_restore(work, nao, ortho, ksa, &
2548 "SOLVE", pos="LEFT", transa="T")
2549 CASE (cholesky_inverse)
2550 CALL cp_fm_uplo_to_full(ksa, work)
2551 CALL cp_fm_cholesky_restore(ksa, nao, ortho, work, &
2552 "MULTIPLY", pos="RIGHT")
2553 CALL cp_fm_cholesky_restore(work, nao, ortho, ksa, &
2554 "MULTIPLY", pos="LEFT", transa="T")
2555 CASE (cholesky_off)
2556 CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, ortho, 0.0_dp, work)
2557 CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, ortho, work, 0.0_dp, ksa)
2558 END SELECT
2559
2560 END IF
2561
2562 ! Diagonalization of the ROKS operator matrix
2563
2564 CALL choose_eigv_solver(ksa, work, eiga)
2565
2566 ! Back-transformation of the orthonormal eigenvectors if needed
2567
2568 IF (level_shift_loc /= 0.0_dp) THEN
2569 ! Use old MO set for back-transformation if level-shifting was applied
2570 CALL cp_fm_to_fm(moa, ortho)
2571 CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, ortho, work, 0.0_dp, moa)
2572 ELSE
2573 IF (orthogonal_basis) THEN
2574 CALL cp_fm_to_fm(work, moa)
2575 ELSE
2576 SELECT CASE (scf_env%cholesky_method)
2577 CASE (cholesky_reduce, cholesky_restore)
2578 CALL cp_fm_cholesky_restore(work, nmo, ortho, moa, "SOLVE")
2579 CASE (cholesky_inverse)
2580 CALL cp_fm_cholesky_restore(work, nmo, ortho, moa, "MULTIPLY")
2581 CASE (cholesky_off)
2582 CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, ortho, work, 0.0_dp, moa)
2583 END SELECT
2584 END IF
2585 END IF
2586
2587 ! Correct MO eigenvalues, if level-shifting was applied
2588
2589 IF (level_shift_loc /= 0.0_dp) THEN
2590 DO imo = homob + 1, homoa
2591 eiga(imo) = eiga(imo) - 0.5_dp*level_shift_loc
2592 END DO
2593 DO imo = homoa + 1, nmo
2594 eiga(imo) = eiga(imo) - level_shift_loc
2595 END DO
2596 END IF
2597
2598 ! Update also the beta MO set
2599
2600 eigb(:) = eiga(:)
2601 CALL cp_fm_to_fm(moa, mob)
2602
2603 ! Calculate the new alpha and beta density matrix
2604
2605 ! does not yet handle k-points
2606 CALL calculate_density_matrix(mos(1), scf_env%p_mix_new(1, 1)%matrix)
2607 CALL calculate_density_matrix(mos(2), scf_env%p_mix_new(2, 1)%matrix)
2608
2609 CALL timestop(handle)
2610
2611 END SUBROUTINE do_roks_diag
2612
2613! **************************************************************************************************
2614!> \brief iterative diagonalization using the block Krylov-space approach
2615!> \param scf_env ...
2616!> \param mos ...
2617!> \param matrix_ks ...
2618!> \param scf_control ...
2619!> \param scf_section ...
2620!> \param check_moconv_only ...
2621!> \param
2622!> \par History
2623!> 05.2009 created [MI]
2624! **************************************************************************************************
2625
2626 SUBROUTINE do_block_krylov_diag(scf_env, mos, matrix_ks, &
2627 scf_control, scf_section, check_moconv_only)
2628
2629 TYPE(qs_scf_env_type), POINTER :: scf_env
2630 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
2631 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
2632 TYPE(scf_control_type), POINTER :: scf_control
2633 TYPE(section_vals_type), POINTER :: scf_section
2634 LOGICAL, INTENT(IN), OPTIONAL :: check_moconv_only
2635
2636 CHARACTER(LEN=*), PARAMETER :: routinen = 'do_block_krylov_diag'
2637 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
2638
2639 INTEGER :: handle, homo, ispin, iter, nao, nmo, &
2640 output_unit
2641 LOGICAL :: converged, my_check_moconv_only
2642 REAL(dp) :: eps_iter, t1, t2
2643 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues
2644 TYPE(cp_fm_type), POINTER :: c0, c1, chc, evec, ks, mo_coeff, ortho, &
2645 work
2646 TYPE(cp_logger_type), POINTER :: logger
2647
2648 logger => cp_get_default_logger()
2649 CALL timeset(routinen, handle)
2650
2651 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%LANCZOS", &
2652 extension=".scfLog")
2653
2654 my_check_moconv_only = .false.
2655 IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
2656
2657 NULLIFY (mo_coeff, ortho, work, ks)
2658 NULLIFY (mo_eigenvalues)
2659 NULLIFY (c0, c1)
2660
2661 IF (scf_env%cholesky_method == cholesky_inverse) THEN
2662 ortho => scf_env%ortho_m1
2663 ELSE
2664 ortho => scf_env%ortho
2665 END IF
2666 work => scf_env%scf_work2
2667
2668 DO ispin = 1, SIZE(matrix_ks)
2669 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, &
2670 scf_env%scf_work1(ispin), scf_env%fm_copy_plan(ispin), check_plan=.true.)
2671 END DO
2672
2673 IF (scf_env%mixing_method == 1) THEN
2674 scf_env%iter_param = scf_env%p_mix_alpha
2675 scf_env%iter_method = "P_Mix/Lanczos"
2676 ELSE
2677 scf_env%iter_param = scf_env%mixing_store%alpha
2678 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Lanc."
2679 END IF
2680
2681 DO ispin = 1, SIZE(matrix_ks)
2682
2683 ks => scf_env%scf_work1(ispin)
2684 CALL cp_fm_uplo_to_full(ks, work)
2685
2686 CALL get_mo_set(mo_set=mos(ispin), &
2687 nao=nao, &
2688 nmo=nmo, &
2689 homo=homo, &
2690 eigenvalues=mo_eigenvalues, &
2691 mo_coeff=mo_coeff)
2692
2693 NULLIFY (c0, c1)
2694 c0 => scf_env%krylov_space%mo_conv(ispin)
2695 c1 => scf_env%krylov_space%mo_refine(ispin)
2696 SELECT CASE (scf_env%cholesky_method)
2697 CASE (cholesky_reduce)
2698 CALL cp_fm_cholesky_reduce(ks, ortho)
2699 CALL cp_fm_uplo_to_full(ks, work)
2700 CALL cp_fm_cholesky_restore(mo_coeff, nmo, ortho, c0, "MULTIPLY")
2701 CASE (cholesky_restore)
2702 CALL cp_fm_cholesky_restore(ks, nao, ortho, work, &
2703 "SOLVE", pos="RIGHT")
2704 CALL cp_fm_cholesky_restore(work, nao, ortho, ks, &
2705 "SOLVE", pos="LEFT", transa="T")
2706 CALL cp_fm_cholesky_restore(mo_coeff, nmo, ortho, c0, "MULTIPLY")
2707 CASE (cholesky_inverse)
2708 CALL cp_fm_cholesky_restore(ks, nao, ortho, work, &
2709 "MULTIPLY", pos="RIGHT")
2710 CALL cp_fm_cholesky_restore(work, nao, ortho, ks, &
2711 "MULTIPLY", pos="LEFT", transa="T")
2712 CALL cp_fm_cholesky_restore(mo_coeff, nmo, ortho, c0, "SOLVE")
2713 END SELECT
2714
2715 scf_env%krylov_space%nmo_nc = nmo
2716 scf_env%krylov_space%nmo_conv = 0
2717
2718 t1 = m_walltime()
2719 IF (output_unit > 0) THEN
2720 WRITE (output_unit, "(/T15,A)") '<<<<<<<<< LANCZOS REFINEMENT <<<<<<<<<<'
2721 WRITE (output_unit, "(T8,A,T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") &
2722 " Spin ", " Cycle ", &
2723 " conv. MOS ", " B2MAX ", " B2MIN ", " Time", repeat("-", 60)
2724 END IF
2725 eps_iter = kp_adaptive_eps_iter(scf_env%krylov_settings%eps_conv, scf_env%krylov_settings%eps_adapt, &
2726 scf_env%iter_delta, 1.0_dp, huge(1.0_dp))
2727 iter = 0
2728 converged = .false.
2729 !Check convergence of MOS
2730 IF (my_check_moconv_only) THEN
2731
2732 CALL lanczos_refinement(scf_env%krylov_space, ks, c0, c1, mo_eigenvalues, &
2733 nao, eps_iter, ispin, check_moconv_only=my_check_moconv_only)
2734 t2 = m_walltime()
2735 IF (output_unit > 0) THEN
2736 WRITE (output_unit, '(T8,I3,T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
2737 ispin, iter, scf_env%krylov_space%nmo_conv, &
2738 scf_env%krylov_space%max_res_norm, scf_env%krylov_space%min_res_norm, t2 - t1
2739 END IF
2740
2741 cycle
2742 ELSE
2743 !Block Lanczos refinement
2744 DO iter = 1, scf_env%krylov_settings%max_iter
2745 CALL lanczos_refinement_2v(scf_env%krylov_space, ks, c0, c1, mo_eigenvalues, &
2746 nao, eps_iter, ispin)
2747 t2 = m_walltime()
2748 IF (output_unit > 0) THEN
2749 WRITE (output_unit, '(T8,I3,T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
2750 ispin, iter, scf_env%krylov_space%nmo_conv, &
2751 scf_env%krylov_space%max_res_norm, scf_env%krylov_space%min_res_norm, t2 - t1
2752 END IF
2753 t1 = m_walltime()
2754 IF (scf_env%krylov_space%max_res_norm < eps_iter) THEN
2755 converged = .true.
2756 IF (output_unit > 0) WRITE (output_unit, *) &
2757 " Reached convergence in ", iter, " iterations "
2758 EXIT
2759 END IF
2760 END DO
2761
2762 IF (.NOT. converged .AND. output_unit > 0) THEN
2763 WRITE (output_unit, "(T4, A)") " WARNING Lanczos refinement could "// &
2764 "not converge all the mos:"
2765 WRITE (output_unit, "(T40,A,T70,I10)") " number of not converged mos ", &
2766 scf_env%krylov_space%nmo_nc
2767 WRITE (output_unit, "(T40,A,T70,E10.2)") " max norm of the residual ", &
2768 scf_env%krylov_space%max_res_norm
2769
2770 END IF
2771
2772 ! For the moment skip the re-orthogonalization
2773 IF (.false.) THEN
2774 !Re-orthogonalization
2775 NULLIFY (chc, evec)
2776 chc => scf_env%krylov_space%chc_mat(ispin)
2777 evec => scf_env%krylov_space%c_vec(ispin)
2778 CALL parallel_gemm('N', 'N', nao, nmo, nao, rone, ks, c0, rzero, work)
2779 CALL parallel_gemm('T', 'N', nmo, nmo, nao, rone, c0, work, rzero, chc)
2780 !Diagonalize (C^t)HC
2781 CALL choose_eigv_solver(chc, evec, mo_eigenvalues)
2782 !Rotate the C vectors
2783 CALL parallel_gemm('N', 'N', nao, nmo, nmo, rone, c0, evec, rzero, c1)
2784 c0 => scf_env%krylov_space%mo_refine(ispin)
2785 END IF
2786
2787 IF (scf_env%cholesky_method == cholesky_inverse) THEN
2788 CALL cp_fm_cholesky_restore(c0, nmo, ortho, mo_coeff, "MULTIPLY")
2789 ELSE
2790 CALL cp_fm_cholesky_restore(c0, nmo, ortho, mo_coeff, "SOLVE")
2791 END IF
2792
2793 IF (.NOT. scf_control%gce%do_gce) THEN
2794 CALL set_mo_occupation(mo_set=mos(ispin), &
2795 smear=scf_control%smear)
2796 ELSE
2797 CALL set_mo_occupation(mo_set=mos(ispin), &
2798 smear=scf_control%smear, &
2799 gce=scf_control%gce)
2800 END IF
2801
2802 ! does not yet handle k-points
2803 CALL calculate_density_matrix(mos(ispin), &
2804 scf_env%p_mix_new(ispin, 1)%matrix)
2805 END IF
2806 END DO ! ispin
2807
2808 IF (output_unit > 0) THEN
2809 WRITE (output_unit, "(T15,A/)") '<<<<<<<<< END LANCZOS REFINEMENT <<<<<<<<<<'
2810 END IF
2811
2812 CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
2813 "PRINT%LANCZOS")
2814
2815 CALL timestop(handle)
2816
2817 END SUBROUTINE do_block_krylov_diag
2818
2819! **************************************************************************************************
2820!> \brief iterative diagonalization using the block davidson space approach
2821!> \param qs_env ...
2822!> \param scf_env ...
2823!> \param mos ...
2824!> \param matrix_ks ...
2825!> \param matrix_s ...
2826!> \param scf_control ...
2827!> \param scf_section ...
2828!> \param check_moconv_only ...
2829!> \param
2830!> \par History
2831!> 05.2011 created [MI]
2832! **************************************************************************************************
2833
2834 SUBROUTINE do_block_davidson_diag(qs_env, scf_env, mos, matrix_ks, matrix_s, &
2835 scf_control, scf_section, check_moconv_only)
2836
2837 TYPE(qs_environment_type), POINTER :: qs_env
2838 TYPE(qs_scf_env_type), POINTER :: scf_env
2839 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
2840 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
2841 TYPE(scf_control_type), POINTER :: scf_control
2842 TYPE(section_vals_type), POINTER :: scf_section
2843 LOGICAL, INTENT(IN), OPTIONAL :: check_moconv_only
2844
2845 CHARACTER(LEN=*), PARAMETER :: routinen = 'do_block_davidson_diag'
2846
2847 INTEGER :: handle, ispin, nspins, output_unit
2848 LOGICAL :: do_prec, my_check_moconv_only
2849 TYPE(cp_logger_type), POINTER :: logger
2850
2851 logger => cp_get_default_logger()
2852 CALL timeset(routinen, handle)
2853
2854 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DAVIDSON", &
2855 extension=".scfLog")
2856
2857 IF (output_unit > 0) THEN
2858 WRITE (output_unit, "(/T15,A)") '<<<<<<<<< DAVIDSON ITERATIONS <<<<<<<<<<'
2859 END IF
2860
2861 IF (scf_env%mixing_method == 1) THEN
2862 scf_env%iter_param = scf_env%p_mix_alpha
2863 scf_env%iter_method = "P_Mix/Dav."
2864 ELSE
2865 scf_env%iter_param = scf_env%mixing_store%alpha
2866 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Dav."
2867 END IF
2868
2869 my_check_moconv_only = .false.
2870 IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
2871 do_prec = .false.
2872 IF (scf_env%block_davidson_env(1)%prec_type /= 0 .AND. &
2873 scf_env%iter_count >= scf_env%block_davidson_env(1)%first_prec) THEN
2874 do_prec = .true.
2875 END IF
2876
2877 nspins = SIZE(matrix_ks)
2878
2879 IF (do_prec .AND. (scf_env%iter_count == scf_env%block_davidson_env(1)%first_prec .OR. &
2880 modulo(scf_env%iter_count, scf_env%block_davidson_env(1)%niter_new_prec) == 0)) THEN
2881 CALL restart_preconditioner(qs_env, scf_env%ot_preconditioner, &
2882 prec_type=scf_env%block_davidson_env(1)%prec_type, nspins=nspins)
2883 CALL prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, scf_env%ot_preconditioner, &
2884 scf_env%block_davidson_env(1)%prec_type, &
2885 scf_env%block_davidson_env(1)%solver_type, &
2886 scf_env%block_davidson_env(1)%energy_gap, nspins, &
2887 convert_to_dbcsr=scf_env%block_davidson_env(1)%use_sparse_mos, &
2888 full_mo_set=.true.)
2889 END IF
2890
2891 DO ispin = 1, nspins
2892 IF (scf_env%block_davidson_env(ispin)%use_sparse_mos) THEN
2893 IF (.NOT. do_prec) THEN
2894 CALL generate_extended_space_sparse(scf_env%block_davidson_env(ispin), mos(ispin), &
2895 matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit)
2896 ELSE
2897 CALL generate_extended_space_sparse(scf_env%block_davidson_env(ispin), mos(ispin), &
2898 matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit, &
2899 scf_env%ot_preconditioner(ispin)%preconditioner)
2900 END IF
2901
2902 ELSE
2903 IF (.NOT. do_prec) THEN
2904 CALL generate_extended_space(scf_env%block_davidson_env(ispin), mos(ispin), &
2905 matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit)
2906 ELSE
2907 CALL generate_extended_space(scf_env%block_davidson_env(ispin), mos(ispin), &
2908 matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit, &
2909 scf_env%ot_preconditioner(ispin)%preconditioner)
2910 END IF
2911 END IF
2912 END DO !ispin
2913
2914 IF (.NOT. scf_control%gce%do_gce) THEN
2915 CALL set_mo_occupation(mo_array=mos, &
2916 smear=scf_control%smear)
2917 ELSE
2918 CALL set_mo_occupation(mo_array=mos, &
2919 smear=scf_control%smear, &
2920 gce=scf_control%gce)
2921 END IF
2922
2923 DO ispin = 1, nspins
2924 ! does not yet handle k-points
2925 CALL calculate_density_matrix(mos(ispin), &
2926 scf_env%p_mix_new(ispin, 1)%matrix)
2927 END DO
2928
2929 IF (output_unit > 0) THEN
2930 WRITE (output_unit, "(T15,A/)") '<<<<<<<<< END DAVIDSON ITERATION <<<<<<<<<<'
2931 END IF
2932
2933 CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
2934 "PRINT%DAVIDSON")
2935
2936 CALL timestop(handle)
2937
2938 END SUBROUTINE do_block_davidson_diag
2939
2940! **************************************************************************************************
2941!> \brief block-Davidson SCF step for complex K points, one channel per
2942!> (local kpoint, spin); mirrors do_block_davidson_diag for the
2943!> preconditioner scheduling and do_general_diag_kp for the H(k)/S(k)
2944!> assembly and the density exit chain
2945!> \param qs_env ...
2946!> \param scf_env ...
2947!> \param matrix_ks_kp real-space Kohn-Sham matrices, dbcsr image set
2948!> \param matrix_s_kp real-space overlap matrix, dbcsr image set
2949!> \param kpoints ...
2950!> \param scf_control ...
2951!> \param added_mos_auto_grow ...
2952! **************************************************************************************************
2953 SUBROUTINE do_block_davidson_diag_kp(qs_env, scf_env, matrix_ks_kp, matrix_s_kp, kpoints, &
2954 scf_control, added_mos_auto_grow)
2955
2956 TYPE(qs_environment_type), POINTER :: qs_env
2957 TYPE(qs_scf_env_type), POINTER :: scf_env
2958 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
2959 TYPE(kpoint_type), POINTER :: kpoints
2960 TYPE(scf_control_type), POINTER :: scf_control
2961 LOGICAL, INTENT(OUT), OPTIONAL :: added_mos_auto_grow
2962
2963 CHARACTER(LEN=*), PARAMETER :: routinen = 'do_block_davidson_diag_kp'
2964 REAL(kind=dp), PARAMETER :: empty_factor = 5.0_dp, &
2965 eps_iter_cap = 1.0e-4_dp
2966
2967 INTEGER :: channel, first_prec, handle, ispin, iw, &
2968 kplocal, kpmax, local_kpoint, &
2969 nchannels, nspin, output_unit
2970 INTEGER, DIMENSION(2) :: kp_range
2971 INTEGER, DIMENSION(:, :), POINTER :: kp_dist
2972 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2973 LOGICAL :: do_prec, my_grow, use_real_wfn
2974 REAL(kind=dp) :: eps_iter, eps_iter_empty
2975 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
2976 TYPE(cp_cfm_type) :: cksmat, csmat, csmat_base
2977 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
2978 TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
2979 TYPE(cp_logger_type), POINTER :: logger
2980 TYPE(kpoint_env_type), POINTER :: kp
2981 TYPE(kpoint_operator_context_type) :: op_ctx
2982 TYPE(mp_para_env_type), POINTER :: para_env
2983 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2984 POINTER :: sab_nl
2985 TYPE(section_vals_type), POINTER :: scf_section
2986
2987! bounds of the adaptive inner tolerance. The loose phase is capped for
2988! the residual-norm semantics of EPS_ITER. A residual of order 1/k would
2989! leave an eigenvector error of order 10 percent at a typical metallic
2990! gap. Unoccupied MOS converge empty_factor times looser. The tracking
2991! divides the mixing delta by noise_scale: the SCF convergence measure
2992! is the LARGEST ELEMENT of the density residual, in which the
2993! per-column residuals of a whole block accumulate. The measured
2994! amplification is ~20 on a smeared Bi2S3/MoS2 slab, so tracking at
2995! delta/100 keeps the injected noise below the outer convergence level.
2996
2997 logger => cp_get_default_logger()
2998 CALL timeset(routinen, handle)
2999
3000 scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
3001 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DAVIDSON", &
3002 extension=".scfLog")
3003 IF (output_unit > 0) THEN
3004 WRITE (output_unit, "(/T15,A)") '<<<<<<<<< DAVIDSON ITERATIONS <<<<<<<<<<'
3005 END IF
3006
3007 CALL get_kpoint_info(kpoints, kp_range=kp_range, kp_dist=kp_dist, xkp=xkp, &
3008 use_real_wfn=use_real_wfn, sab_nl=sab_nl, &
3009 cell_to_index=cell_to_index)
3010 IF (use_real_wfn) THEN
3011 cpabort('Davidson Kp requires complex K-point MOS')
3012 END IF
3013 para_env => kpoints%blacs_env_all%para_env
3014 kplocal = kp_range(2) - kp_range(1) + 1
3015 kpmax = maxval(kp_dist(2, :) - kp_dist(1, :) + 1)
3016 nspin = SIZE(matrix_ks_kp, 1)
3017 nchannels = kplocal*nspin
3018
3019 IF (scf_env%mixing_method == 1) THEN
3020 scf_env%iter_param = scf_env%p_mix_alpha
3021 scf_env%iter_method = "P_Mix/Dav."
3022 ELSE IF (scf_env%mixing_method > 1) THEN
3023 scf_env%iter_param = scf_env%mixing_store%alpha
3024 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Dav."
3025 END IF
3026
3027 ! preconditioner scheduling, mirrors do_block_davidson_diag. The first
3028 ! K-point SCF step cold-starts unless extrapolation prefilled the MOS
3029 ! (mos_prefilled). The FULL_ALL maker consumes the MO coefficients and
3030 ! eigenvalues, so building before the first solve completed would feed it
3031 ! uninitialized data and stall the Davidson convergence. Keep first_prec
3032 ! clamped past that step in both cases. Warming it up under prefilled MOS
3033 ! would mix two variables at once.
3034 first_prec = max(scf_env%block_davidson_env(1)%first_prec, 2)
3035 do_prec = .false.
3036 IF (scf_env%block_davidson_env(1)%prec_type /= ot_precond_none .AND. &
3037 scf_env%iter_count >= first_prec) THEN
3038 do_prec = .true.
3039 END IF
3040 ! adaptive inner tolerance from EPS_ADAPT: loose while the density is
3041 ! still far from self-consistency, tightening with the outer error
3042 ! scaled by noise_scale (mirrors the block Krylov scheduling in
3043 ! do_block_krylov_diag). The
3044 ! ratchet is monotone across SCF steps: the threshold may tighten but
3045 ! never relax again. A rebounding outer error would otherwise re-loosen
3046 ! the solver, and the run could stop on a near-zero mixing delta of a
3047 ! pair of identical, sloppily converged densities. The cap bounds the
3048 ! early loose phase. Unoccupied MOS converge empty_factor times
3049 ! looser (the occupancy decides which columns qualify, see
3050 ! generate_extended_space_c). There is no restoration of the
3051 ! user threshold at the end. The tolerance stays a margin below the
3052 ! outer error. The MOS the run exits with satisfy that same margin
3053 ! rather than the nominal EPS_ITER.
3054 eps_iter = scf_env%block_davidson_env(1)%eps_iter
3055 eps_iter_empty = eps_iter
3056 IF (scf_control%diagonalization%eps_adapt > 0.0_dp .AND. scf_env%iter_delta > 0.0_dp) THEN
3057 ! stall escape: when the mixing delta sank below the tolerance the
3058 ! last solve used, that solve's own sloppiness may be the noise floor
3059 ! pinning the delta. Tighten one decade below what that solve used.
3060 ! Tightening relative to the OPERATING tolerance instead of the
3061 ! dipped delta bounds the reaction to transient mixing dips (a
3062 ! momentary overshoot must not collapse the threshold by orders of
3063 ! magnitude and destabilise the mixer), while a genuine stall moves
3064 ! the threshold down decade by decade. Checking on the stall alone
3065 ! keeps the loose phase as loose as EPS_ADAPT allows
3066 IF (scf_env%iter_delta < scf_env%block_davidson_env(1)%eps_iter_used) THEN
3067 scf_env%block_davidson_env(1)%eps_iter_prev = &
3068 min(scf_env%block_davidson_env(1)%eps_iter_prev, &
3069 max(scf_env%block_davidson_env(1)%eps_iter, &
3070 0.1_dp*scf_env%block_davidson_env(1)%eps_iter_used))
3071 END IF
3072 ! candidate: looser than the user value, capped; ratchet: only
3073 ! tighter than the previous adaptive value, which lives in the
3074 ! persistent channel settings (a local variable would reset on
3075 ! every call into this driver)
3076 eps_iter = min(max(eps_iter, min(scf_control%diagonalization%eps_adapt*scf_env%iter_delta/ &
3077 scf_env%block_davidson_env(1)%noise_scale, &
3078 eps_iter_cap)), scf_env%block_davidson_env(1)%eps_iter_prev)
3079 eps_iter_empty = empty_factor*eps_iter
3080 ! record what this solve uses, after every override above: the outer
3081 ! loop compares the mixing delta against it (qs_scf_check_inner_exit)
3082 scf_env%block_davidson_env(1)%eps_iter_prev = eps_iter
3083 scf_env%block_davidson_env(1)%eps_iter_used = eps_iter
3084 END IF
3085 ! schedule trace of the adaptive tolerance, one line per SCF step at
3086 ! DEBUG print level
3087 IF (scf_control%diagonalization%eps_adapt > 0.0_dp .AND. para_env%is_source() .AND. &
3088 logger%iter_info%print_level >= debug_print_level) THEN
3089 iw = cp_logger_get_default_unit_nr()
3090 IF (iw > 0) THEN
3091 WRITE (unit=iw, fmt="(T3,'EPS_ADAPT step',I5,' delta=',ES9.2,' eps_occ=',ES9.2,' eps_empty=',ES9.2)") &
3092 scf_env%iter_count, scf_env%iter_delta, eps_iter, eps_iter_empty
3093 END IF
3094 END IF
3095
3096 ! work matrices: global scratch fms (also the operator-context
3097 ! template) and the complex H(k)/S(k)
3098 CALL kp_diag_setup(qs_env, kpoints, nspin, matrix_ks_kp, matrix_s_kp, fmwork, op_ctx, &
3099 matrix_struct, .true.)
3100 CALL cp_cfm_create(cksmat, matrix_struct, name="dav_hk")
3101 CALL cp_cfm_create(csmat, matrix_struct, name="dav_sk")
3102 CALL cp_cfm_create(csmat_base, matrix_struct, name="dav_sk_base")
3103
3104 IF (do_prec .AND. (scf_env%iter_count == first_prec .OR. &
3105 modulo(scf_env%iter_count, &
3106 scf_env%block_davidson_env(1)%niter_new_prec) == 0)) THEN
3107 CALL davidson_kp_prepare_preconditioners( &
3108 qs_env, scf_env, matrix_ks_kp, matrix_s_kp, kpoints, op_ctx)
3109 END IF
3110
3111 DO local_kpoint = 1, kpmax
3112 ! Every k-point group executes the same sequence of global transforms.
3113 ! Groups without a k point in the last round still participate
3114 ! in the service calls, but skip their local solver work.
3115 CALL kpoint_operator_get(op_ctx, local_kpoint, kpoint_spin_free, matrix_s_kp, &
3116 cmat=csmat_base, matrix_row=kpoint_spin_free, slot=kpoint_slot_s)
3117
3118 DO ispin = 1, nspin
3119 CALL kpoint_operator_get(op_ctx, local_kpoint, ispin, matrix_ks_kp, &
3120 cmat=cksmat, matrix_row=ispin, slot=kpoint_slot_ks)
3121 IF (local_kpoint > kplocal) cycle
3122 kp => kpoints%kp_env(local_kpoint)%kpoint_env
3123 CALL cp_cfm_to_cfm(csmat_base, csmat)
3124
3125 IF (scf_env%iter_count == 1 .AND. .NOT. kp%mos_prefilled) THEN
3126 ! K-point MOS enter the first SCF step allocated but uninitialized.
3127 ! The stored S(k) copy is restored for every spin because
3128 ! cp_cfm_geeig overwrites its matrix arguments.
3129 CALL cold_start_mo_init(kp, ispin, cksmat, csmat)
3130 ELSE
3131 channel = (local_kpoint - 1)*nspin + ispin
3132 IF (do_prec) THEN
3133 CALL generate_extended_space_c(scf_env%block_davidson_env(channel), &
3134 kp%mos(ispin), cksmat, csmat, output_unit, &
3135 eps_iter, eps_iter_empty, &
3136 scf_env%ot_preconditioner(channel)%preconditioner)
3137 ELSE
3138 CALL generate_extended_space_c(scf_env%block_davidson_env(channel), &
3139 kp%mos(ispin), cksmat, csmat, output_unit, &
3140 eps_iter, eps_iter_empty)
3141 END IF
3142 END IF
3143 END DO
3144 END DO
3145
3146 ! occupations and density: the shared exit chain, same as every
3147 ! other k-point diagonalization driver
3148 CALL kp_diag_exit_chain(kpoints, scf_env, scf_control, matrix_s_kp, sab_nl, fmwork, my_grow)
3149 IF (my_grow .AND. PRESENT(added_mos_auto_grow)) added_mos_auto_grow = .true.
3150
3151 CALL cp_cfm_release(cksmat)
3152 CALL cp_cfm_release(csmat)
3153 CALL cp_cfm_release(csmat_base)
3154 CALL kp_diag_finalize(qs_env, fmwork, op_ctx)
3155
3156 IF (output_unit > 0) THEN
3157 WRITE (output_unit, "(T15,A/)") '<<<<<<<<< END DAVIDSON ITERATION <<<<<<<<<<'
3158 END IF
3159 CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
3160 "PRINT%DAVIDSON")
3161
3162 CALL timestop(handle)
3163
3164 END SUBROUTINE do_block_davidson_diag_kp
3165
3166! **************************************************************************************************
3167!> \brief Block Krylov (Lanczos) refinement of the k-point MOS, the complex
3168!> counterpart of do_block_krylov_diag: every channel (local kpoint,
3169!> spin) is standardized with the Cholesky factor of S(k), refined by
3170!> the complex block-Krylov kernel, and rotated back into the
3171!> generalized problem
3172!> \param qs_env ...
3173!> \param scf_env ...
3174!> \param matrix_ks_kp real-space Kohn-Sham matrix, dbcsr image set
3175!> \param matrix_s_kp real-space overlap matrix, dbcsr image set
3176!> \param kpoints ...
3177!> \param scf_control ...
3178!> \param check_moconv_only only report the convergence statistics
3179!> \param added_mos_auto_grow signals the outer loop that the virtual-space
3180!> buffer must grow before this SCF step is redone
3181! **************************************************************************************************
3182 SUBROUTINE do_block_krylov_diag_kp(qs_env, scf_env, matrix_ks_kp, matrix_s_kp, kpoints, &
3183 scf_control, check_moconv_only, added_mos_auto_grow)
3184
3185 TYPE(qs_environment_type), POINTER :: qs_env
3186 TYPE(qs_scf_env_type), POINTER :: scf_env
3187 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
3188 TYPE(kpoint_type), POINTER :: kpoints
3189 TYPE(scf_control_type), POINTER :: scf_control
3190 LOGICAL, INTENT(IN), OPTIONAL :: check_moconv_only
3191 LOGICAL, INTENT(OUT), OPTIONAL :: added_mos_auto_grow
3192
3193 CHARACTER(LEN=*), PARAMETER :: routinen = 'do_block_krylov_diag_kp'
3194 REAL(kind=dp), PARAMETER :: empty_factor = 5.0_dp, &
3195 eps_iter_cap = 1.0e-4_dp, &
3196 noise_scale = 100.0_dp
3197
3198 INTEGER :: handle, ik, ikp, ispin, kplocal, kpmax, &
3199 local_kpoint, max_iter, nentry, nspin, &
3200 output_unit
3201 INTEGER, DIMENSION(2) :: kp_range
3202 INTEGER, DIMENSION(:, :), POINTER :: kp_dist
3203 LOGICAL :: my_check_moconv_only, my_grow, &
3204 prepared_ops, use_local_kpoints, &
3205 use_real_wfn
3206 REAL(kind=dp) :: eps_iter, eps_iter_empty
3207 TYPE(cp_fm_struct_type), POINTER :: matrix_struct, mo_struct
3208 TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
3209 TYPE(cp_logger_type), POINTER :: logger
3210 TYPE(kp_diag_workspace_type), ALLOCATABLE, TARGET :: workspace(:)
3211 TYPE(kpoint_env_type), POINTER :: kp
3212 TYPE(kpoint_operator_context_type) :: op_ctx
3213 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3214 POINTER :: sab_nl
3215 TYPE(section_vals_type), POINTER :: scf_section
3216
3217! empty_factor and noise_scale mirror the defaults of the
3218! Davidson k-point driver. eps_iter_cap bounds the adaptive tolerance
3219
3220 logger => cp_get_default_logger()
3221 CALL timeset(routinen, handle)
3222
3223 my_check_moconv_only = .false.
3224 IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
3225
3226 scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
3227 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%LANCZOS", &
3228 extension=".scfLog")
3229
3230 CALL get_kpoint_info(kpoints, kp_range=kp_range, kp_dist=kp_dist, &
3231 use_real_wfn=use_real_wfn, sab_nl=sab_nl)
3232 IF (use_real_wfn) THEN
3233 cpabort('k-point Lanczos refinement requires complex k-point MOS')
3234 END IF
3235 kplocal = kp_range(2) - kp_range(1) + 1
3236 ! the largest local k-point count across the groups: kp_dist can be
3237 ! uneven, and every group must join the service calls for every round
3238 kpmax = maxval(kp_dist(2, :) - kp_dist(1, :) + 1)
3239 nspin = SIZE(matrix_ks_kp, 1)
3240
3241 IF (scf_env%mixing_method == 1) THEN
3242 scf_env%iter_param = scf_env%p_mix_alpha
3243 scf_env%iter_method = "P_Mix/Lanczos"
3244 ELSE
3245 scf_env%iter_param = scf_env%mixing_store%alpha
3246 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Lanc."
3247 END IF
3248
3249 ! work matrices: the driver workspace below holds the complex
3250 ! H(k)/S(k) of the channels
3251
3252 ! adaptive inner tolerance, mirrors do_block_krylov_diag; it carries
3253 ! no ratchet state, the threshold follows the sinking mixing delta
3254 eps_iter = kp_adaptive_eps_iter(scf_env%krylov_settings%eps_conv, scf_env%krylov_settings%eps_adapt, &
3255 scf_env%iter_delta, noise_scale, eps_iter_cap)
3256 ! unoccupied MOS relax like the Davidson k-point driver: the
3257 ! occupation decides which columns use the relaxed empty threshold,
3258 ! not the HOMO index (see generate_extended_space_c)
3259 eps_iter_empty = eps_iter
3260 IF (scf_env%krylov_settings%eps_adapt > 0.0_dp) eps_iter_empty = empty_factor*eps_iter
3261 max_iter = scf_env%krylov_settings%max_iter
3262
3263 ! one-rank k-point groups refine their k points locally, without a
3264 ! transfer round; a diag check, a convergence-only probe and an
3265 ! explicitly selected diagonalization library keep the rounds
3266 use_local_kpoints = kpoints%para_env_kp%num_pe == 1
3267 IF (diag_check_requested()) use_local_kpoints = .false.
3268 IF (my_check_moconv_only) use_local_kpoints = .false.
3269 IF (diag_lib_explicit) use_local_kpoints = .false.
3270 ! the transfer rounds assemble their operators lazily through the
3271 ! get path and skip the snapshot
3272 CALL kp_diag_setup(qs_env, kpoints, nspin, matrix_ks_kp, matrix_s_kp, fmwork, op_ctx, &
3273 matrix_struct, use_local_kpoints, prepared_ops)
3274 use_local_kpoints = use_local_kpoints .AND. prepared_ops
3275
3276 nentry = 0
3277 IF (use_local_kpoints) CALL kpoint_operator_work_size(op_ctx, nentry)
3278 kp => kpoints%kp_env(1)%kpoint_env
3279 mo_struct => kp%mos(1)%cmo_coeff%matrix_struct
3280 CALL kp_diag_workspaces_create(workspace, 1, matrix_struct, mo_struct, nentry, &
3281 use_local_kpoints, scf_env%cholesky_method == cholesky_off, &
3282 lanczos=.true., nkrylov=scf_env%krylov_settings%nkrylov, &
3283 nblock=scf_env%krylov_settings%nblock)
3284
3285 IF (output_unit > 0 .AND. .NOT. my_check_moconv_only) THEN
3286 WRITE (output_unit, "(/T15,A)") '<<<<<<<<< LANCZOS REFINEMENT <<<<<<<<<<'
3287 WRITE (output_unit, "(T8,A,T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") &
3288 " Spin ", " Cycle ", &
3289 " conv. MOS ", " B2MAX ", " B2MIN ", " Time", repeat("-", 60)
3290 END IF
3291
3292 IF (use_local_kpoints) THEN
3293 ! serial pass over the locally owned k points. The first SCF step
3294 ! cold-starts every unprefilled channel in place, like the rounds
3295 ! path below. Later steps refine
3296 DO ikp = 1, kplocal
3297 kp => kpoints%kp_env(ikp)%kpoint_env
3298 ik = kp_range(1) + ikp - 1
3299 CALL kpoint_operator_evaluate(op_ctx, 1, kpoint_slot_s, ik, workspace(1)%phase, &
3300 workspace(1)%csmat_base)
3301 DO ispin = 1, nspin
3302 CALL kpoint_operator_evaluate(op_ctx, ispin, kpoint_slot_ks, ik, workspace(1)%phase, &
3303 workspace(1)%cksmat)
3304 IF (scf_env%iter_count == 1 .AND. .NOT. kp%mos_prefilled) THEN
3305 ! the K-point MOS enter the first SCF step allocated but
3306 ! uninitialized
3307 CALL cold_start_mo_init(kp, ispin, workspace(1)%cksmat, workspace(1)%csmat_base)
3308 cycle
3309 END IF
3310 CALL lanczos_refine_channel(kp, ispin, scf_env, my_check_moconv_only, eps_iter, &
3311 eps_iter_empty, max_iter, workspace(1), output_unit, .true.)
3312 END DO
3313 END DO
3314 ELSE
3315
3316 DO local_kpoint = 1, kpmax
3317 ! Every k-point group executes the same sequence of global transforms
3318 ! up to the largest local k-point count (kp_dist can be uneven after
3319 ! #5993). Only the group owning this local_kpoint receives the
3320 ! redistributed matrix. The others participate through the dummy
3321 ! destination. Groups without a k point in the last rounds still
3322 ! join the service calls and skip their local solver work.
3323 CALL kpoint_operator_get(op_ctx, local_kpoint, kpoint_spin_free, matrix_s_kp, &
3324 cmat=workspace(1)%csmat_base, matrix_row=kpoint_spin_free)
3325
3326 DO ispin = 1, nspin
3327 CALL kpoint_operator_get(op_ctx, local_kpoint, ispin, matrix_ks_kp, &
3328 cmat=workspace(1)%cksmat, matrix_row=ispin)
3329 IF (local_kpoint > kplocal) cycle
3330 kp => kpoints%kp_env(local_kpoint)%kpoint_env
3331
3332 IF (scf_env%iter_count == 1 .AND. .NOT. kp%mos_prefilled) THEN
3333 ! first SCF step cold-starts by a direct diagonalization,
3334 ! same entry as the Davidson k-point path
3335 CALL cold_start_mo_init(kp, ispin, workspace(1)%cksmat, workspace(1)%csmat_base)
3336 cycle
3337 END IF
3338
3339 CALL lanczos_refine_channel(kp, ispin, scf_env, my_check_moconv_only, eps_iter, &
3340 eps_iter_empty, max_iter, workspace(1), output_unit, .true.)
3341 END DO
3342 END DO
3343 END IF
3344
3345 ! occupations and density: the shared exit chain. With
3346 ! check_moconv_only the caller falls back to the standard k-point
3347 ! diagonalization. That path supplies MOS and density, and the
3348 ! statistics above are printed for information only
3349 IF (.NOT. my_check_moconv_only) THEN
3350 CALL kp_diag_exit_chain(kpoints, scf_env, scf_control, matrix_s_kp, sab_nl, fmwork, my_grow)
3351 IF (my_grow .AND. PRESENT(added_mos_auto_grow)) added_mos_auto_grow = .true.
3352 END IF
3353
3354 CALL kp_diag_workspaces_release(workspace)
3355 CALL kp_diag_finalize(qs_env, fmwork, op_ctx)
3356
3357 IF (output_unit > 0 .AND. .NOT. my_check_moconv_only) THEN
3358 WRITE (output_unit, "(T15,A/)") '<<<<<<<<< END LANCZOS REFINEMENT <<<<<<<<<<'
3359 END IF
3360 CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
3361 "PRINT%LANCZOS")
3362
3363 CALL timestop(handle)
3364
3365 END SUBROUTINE do_block_krylov_diag_kp
3366
3367! **************************************************************************************************
3368!> \brief Refine the MOS of one (k point, spin) channel. The routine
3369!> projects H(k) onto its Hermitian part, standardizes the channel
3370!> by the Cholesky factor of S(k), sweeps the single-pass kernel,
3371!> and writes the MOS back to the MO sets. The caller owns every
3372!> matrix. The routine allocates no CP2K matrix.
3373!> \param kp k-point environment holding the MOS pair of the channel
3374!> \param ispin spin component
3375!> \param scf_env SCF environment holding the Krylov settings
3376!> \param check_moconv_only only run the entry Rayleigh-Ritz of the kernel
3377!> and report the convergence count
3378!> \param eps_iter convergence threshold of the occupied columns
3379!> \param eps_iter_empty relaxed threshold of the unoccupied columns
3380!> \param max_iter sweeps allowed per channel
3381!> \param work driver workspace of the channel. cksmat holds the
3382!> complex H(k), standardized in place. csmat_base holds the
3383!> untouched S(k). csmat holds its work copy and the Cholesky
3384!> factor after the standardization. cwork holds the
3385!> Hermitian-projection scratch. cmos holds the MOS in the
3386!> standardized basis. lan_ws holds the kernel scratch
3387!> \param output_unit unit of the refinement table. A nonpositive unit
3388!> disables every report
3389!> \param print_row report this channel in the refinement table and in
3390!> the not-converged warnings. The local-kpoint path passes
3391!> .FALSE. and reports nothing, like the local phase of the
3392!> general k-point diagonalization
3393! **************************************************************************************************
3394 SUBROUTINE lanczos_refine_channel(kp, ispin, scf_env, check_moconv_only, eps_iter, eps_iter_empty, &
3395 max_iter, work, output_unit, print_row)
3396 TYPE(kpoint_env_type), POINTER :: kp
3397 INTEGER, INTENT(IN) :: ispin
3398 TYPE(qs_scf_env_type), POINTER :: scf_env
3399 LOGICAL, INTENT(IN) :: check_moconv_only
3400 REAL(kind=dp), INTENT(IN) :: eps_iter, eps_iter_empty
3401 INTEGER, INTENT(IN) :: max_iter
3402 TYPE(kp_diag_workspace_type), INTENT(INOUT) :: work
3403 INTEGER, INTENT(IN) :: output_unit
3404 LOGICAL, INTENT(IN) :: print_row
3405
3406 CHARACTER(LEN=*), PARAMETER :: routinen = 'lanczos_refine_channel'
3407
3408 INTEGER :: handle, iter, nmo, nmo_conv, nmo_nc
3409 LOGICAL :: any_refined, converged, need_suffix_rr, &
3410 next_need_suffix_rr, step_refined
3411 REAL(kind=dp) :: max_res, min_res, t1, t2
3412 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eps_col
3413 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues, occupation
3414 TYPE(cp_cfm_type), POINTER :: cmo_coeff
3415
3416 CALL timeset(routinen, handle)
3417
3418 CALL get_mo_set(mo_set=kp%mos(ispin), nmo=nmo, &
3419 eigenvalues=mo_eigenvalues, occupation_numbers=occupation, &
3420 cmo_coeff=cmo_coeff)
3421
3422 CALL cp_cfm_to_cfm(work%csmat_base, work%csmat)
3423
3424 ! H(k) is Hermitian only to grid accuracy and the kernel forms its
3425 ! residuals with the full matrix: project onto the Hermitian part
3426 ! once per channel
3427 CALL cp_cfm_transpose(work%cksmat, "C", work%cwork)
3428 CALL cp_cfm_scale_and_add(cmplx(0.5_dp, 0.0_dp, kind=dp), work%cksmat, &
3429 cmplx(0.5_dp, 0.0_dp, kind=dp), work%cwork)
3430
3431 ! complex MOS of the channel
3432 CALL cp_cfm_to_cfm(cmo_coeff, work%cmos)
3433
3434 ! standardization S(k) = U^dag U: operator and MOS are pulled
3435 ! through the same transform, X = U C and A = U^-dag H U^-1, so
3436 ! the kernel faces a standard eigenproblem. CHOLESKY OFF is not
3437 ! supported on this path
3438 CALL cp_cfm_cholesky_decompose(work%csmat)
3439 CALL cp_cfm_triangular_multiply(work%csmat, work%cksmat, side="R", invert_tr=.true.)
3440 CALL cp_cfm_triangular_multiply(work%csmat, work%cksmat, side="L", transa_tr="C", invert_tr=.true.)
3441 CALL cp_cfm_triangular_multiply(work%csmat, work%cmos, side="L", invert_tr=.false.)
3442
3443 ! per-column thresholds: the occupation decides which columns
3444 ! qualify for the relaxed empty threshold
3445 ALLOCATE (eps_col(nmo))
3446 CALL kp_eps_col_by_occupation(occupation, eps_iter, eps_iter_empty, eps_col)
3447
3448 converged = .false.
3449 any_refined = .false.
3450 need_suffix_rr = .false.
3451
3452 t1 = m_walltime()
3453 CALL lanczos_channel_begin(work%cksmat, work%cmos, mo_eigenvalues, eps_col, &
3454 work%lan_ws, max_res, min_res, nmo_conv, nmo_nc, &
3455 check_moconv_only=check_moconv_only)
3456 IF (check_moconv_only) THEN
3457 DEALLOCATE (eps_col)
3458 IF (print_row .AND. output_unit > 0) THEN
3459 WRITE (output_unit, '(T8,I3,T16,I5,T24,I6,T33,E12.4,2x,E12.4)') &
3460 ispin, 0, nmo_conv, max_res, min_res
3461 END IF
3462 CALL timestop(handle)
3463 RETURN
3464 END IF
3465
3466 IF (nmo_nc == 0) THEN
3467 t2 = m_walltime()
3468 IF (print_row .AND. output_unit > 0 .AND. max_iter >= 1) THEN
3469 WRITE (output_unit, '(T8,I3,T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
3470 ispin, 1, nmo_conv, max_res, min_res, t2 - t1
3471 END IF
3472 converged = .true.
3473 ELSE
3474 DO iter = 1, max_iter
3475 IF (iter > 1) t1 = m_walltime()
3476 CALL lanczos_sweep_suffix(work%cksmat, work%cmos, mo_eigenvalues, eps_col, &
3477 scf_env%krylov_settings%nkrylov, scf_env%krylov_settings%nblock, &
3478 iter == 1, need_suffix_rr, work%lan_ws, &
3479 max_res, min_res, nmo_conv, nmo_nc, step_refined, &
3480 next_need_suffix_rr)
3481 need_suffix_rr = next_need_suffix_rr
3482 any_refined = any_refined .OR. step_refined
3483 t2 = m_walltime()
3484 IF (print_row .AND. output_unit > 0) THEN
3485 WRITE (output_unit, '(T8,I3,T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
3486 ispin, iter, nmo_conv, max_res, min_res, t2 - t1
3487 END IF
3488 IF (nmo_nc == 0) THEN
3489 converged = .true.
3490 EXIT
3491 END IF
3492 END DO
3493 END IF
3494 DEALLOCATE (eps_col)
3495
3496 IF (any_refined) THEN
3497 CALL lanczos_channel_end(work%cmos, work%lan_ws)
3498 END IF
3499 IF (.NOT. converged .AND. print_row .AND. output_unit > 0) THEN
3500 WRITE (output_unit, "(T4, A)") " WARNING Lanczos refinement could "// &
3501 "not converge all the mos:"
3502 WRITE (output_unit, "(T40,A,T70,I10)") " number of not converged mos ", nmo_nc
3503 WRITE (output_unit, "(T40,A,T70,E10.2)") " max norm of the residual ", max_res
3504 END IF
3505
3506 ! rotate the refined MOS back into the generalized problem, C = U^-1 X
3507 CALL cp_cfm_triangular_multiply(work%csmat, work%cmos, side="L", invert_tr=.true.)
3508 CALL cp_cfm_to_cfm(work%cmos, cmo_coeff)
3509
3510 CALL timestop(handle)
3511 END SUBROUTINE lanczos_refine_channel
3512
3513! **************************************************************************************************
3514!> \brief first-SCF-step initialization of the K-point MOS of one spin by one
3515!> direct diagonalization of H(k); both arguments are assembled by the
3516!> caller because cp_cfm_geeig overwrites them
3517!> \param kp kpoint environment holding the complex MO set
3518!> \param ispin ...
3519!> \param cksmat complex H(k), assembled by the caller
3520!> \param csmat complex S(k), assembled by the caller
3521! **************************************************************************************************
3522 SUBROUTINE cold_start_mo_init(kp, ispin, cksmat, csmat)
3523
3524 TYPE(kpoint_env_type), POINTER :: kp
3525 INTEGER, INTENT(IN) :: ispin
3526 TYPE(cp_cfm_type) :: cksmat, csmat
3527
3528 CHARACTER(LEN=*), PARAMETER :: routinen = 'cold_start_mo_init'
3529
3530 INTEGER :: handle
3531 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
3532 TYPE(cp_cfm_type) :: chwork, cswork, cwork
3533
3534 CALL timeset(routinen, handle)
3535
3536 CALL get_mo_set(kp%mos(ispin), eigenvalues=eigenvalues)
3537 ! cp_cfm_geeig overwrites its matrix arguments: solve on scratch
3538 ! copies so the caller H(k)/S(k) stay valid
3539 CALL cp_cfm_create(chwork, cksmat%matrix_struct, name="kp_cold_h")
3540 CALL cp_cfm_create(cswork, csmat%matrix_struct, name="kp_cold_s")
3541 CALL cp_cfm_to_cfm(cksmat, chwork)
3542 CALL cp_cfm_to_cfm(csmat, cswork)
3543 CALL cp_cfm_create(cwork, cksmat%matrix_struct, name="kp_cold_work")
3544 CALL cp_cfm_geeig(chwork, cswork, kp%mos(ispin)%cmo_coeff, eigenvalues, cwork)
3545 CALL cp_cfm_release(chwork)
3546 CALL cp_cfm_release(cswork)
3547 CALL cp_cfm_release(cwork)
3548
3549 CALL timestop(handle)
3550
3551 END SUBROUTINE cold_start_mo_init
3552
3553! **************************************************************************************************
3554!> \brief builds one complex K-point preconditioner per (local kpoint, spin)
3555!> channel; mirrors prepare_qs_ot_kpoint_preconditioners
3556!> \param qs_env ...
3557!> \param scf_env ...
3558!> \param matrix_ks_kp ...
3559!> \param matrix_s_kp ...
3560!> \param kpoints ...
3561!> \param op_ctx ...
3562! **************************************************************************************************
3563 SUBROUTINE davidson_kp_prepare_preconditioners( &
3564 qs_env, scf_env, matrix_ks_kp, matrix_s_kp, kpoints, op_ctx)
3565
3566 TYPE(qs_environment_type), POINTER :: qs_env
3567 TYPE(qs_scf_env_type), POINTER :: scf_env
3568 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
3569 TYPE(kpoint_type), POINTER :: kpoints
3570 TYPE(kpoint_operator_context_type) :: op_ctx
3571
3572 CHARACTER(LEN=*), PARAMETER :: routinen = 'davidson_kp_prepare_preconditioners'
3573
3574 INTEGER :: channel, handle, ispin, kplocal, kpmax, &
3575 local_kpoint, nchannels, nmo, nspin, &
3576 prec_type, solver_type
3577 INTEGER, DIMENSION(2) :: kp_range
3578 INTEGER, DIMENSION(:, :), POINTER :: kp_dist
3579 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3580 REAL(kind=dp) :: energy_gap
3581 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
3582 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
3583 TYPE(cp_blacs_env_type), POINTER :: blacs_env
3584 TYPE(cp_cfm_type) :: matrix_c, matrix_h, matrix_s
3585 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_t_kp
3586 TYPE(dbcsr_type), POINTER :: matrix_h_im, matrix_h_re, matrix_s_im, &
3587 matrix_s_re, matrix_t_im, matrix_t_re
3588 TYPE(kpoint_env_type), POINTER :: kp
3589 TYPE(mp_para_env_type), POINTER :: para_env
3590 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3591 POINTER :: sab_nl
3592
3593 CALL timeset(routinen, handle)
3594
3595 CALL get_kpoint_info(kpoints, kp_range=kp_range, kp_dist=kp_dist, xkp=xkp, sab_nl=sab_nl, &
3596 cell_to_index=cell_to_index)
3597 CALL get_qs_env(qs_env, kinetic_kp=matrix_t_kp)
3598 ! the preconditioner fms must live on the kpoint-group BLACS context,
3599 ! the context of the MOS they are applied with
3600 blacs_env => kpoints%blacs_env
3601 para_env => kpoints%para_env_kp
3602 kplocal = kp_range(2) - kp_range(1) + 1
3603 kpmax = maxval(kp_dist(2, :) - kp_dist(1, :) + 1)
3604 nspin = SIZE(matrix_ks_kp, 1)
3605 nchannels = kplocal*nspin
3606 prec_type = scf_env%block_davidson_env(1)%prec_type
3607 solver_type = scf_env%block_davidson_env(1)%solver_type
3608 energy_gap = scf_env%block_davidson_env(1)%energy_gap
3609
3610 ! channel storage: every (kpoint, spin) channel carries its own
3611 ! preconditioner, because S(k) and H(k) differ per kpoint. The Gamma
3612 ! reuse pattern of restart_preconditioner (a single entry for the
3613 ! spin-independent types) must not be used here.
3614 IF (ASSOCIATED(scf_env%ot_preconditioner)) THEN
3615 DO channel = 1, SIZE(scf_env%ot_preconditioner)
3616 IF (ASSOCIATED(scf_env%ot_preconditioner(channel)%preconditioner)) THEN
3617 CALL destroy_preconditioner(scf_env%ot_preconditioner(channel)%preconditioner)
3618 DEALLOCATE (scf_env%ot_preconditioner(channel)%preconditioner)
3619 END IF
3620 END DO
3621 DEALLOCATE (scf_env%ot_preconditioner)
3622 NULLIFY (scf_env%ot_preconditioner)
3623 END IF
3624 ALLOCATE (scf_env%ot_preconditioner(nchannels))
3625 DO channel = 1, nchannels
3626 ALLOCATE (scf_env%ot_preconditioner(channel)%preconditioner)
3627 CALL init_preconditioner(scf_env%ot_preconditioner(channel)%preconditioner, &
3628 para_env, blacs_env)
3629 END DO
3630
3631 NULLIFY (eigenvalues, matrix_h_im, matrix_h_re, matrix_s_im, matrix_s_re, &
3632 matrix_t_im, matrix_t_re)
3633 IF (prec_type == ot_precond_full_kinetic) THEN
3634 cpassert(ASSOCIATED(matrix_t_kp))
3635 CALL kpoint_operator_prepare(op_ctx, matrix_t_kp, kpoint_slot_t, matrix_row=kpoint_spin_free)
3636 END IF
3637
3638 DO local_kpoint = 1, kpmax
3639 ! S(k) is spin independent. Build it with the same global-group
3640 ! redistribution protocol as the Davidson operator itself, directly
3641 ! as the group-local block-cyclic DBCSR pair of the preconditioner
3642 ! makers.
3643 CALL dbcsr_init_p(matrix_s_re)
3644 CALL dbcsr_init_p(matrix_s_im)
3645 CALL kpoint_operator_get(op_ctx, local_kpoint, kpoint_spin_free, matrix_s_kp, &
3646 mat_re=matrix_s_re, mat_im=matrix_s_im, &
3647 matrix_row=kpoint_spin_free, slot=kpoint_slot_s)
3648
3649 IF (prec_type == ot_precond_full_kinetic) THEN
3650 CALL dbcsr_init_p(matrix_t_re)
3651 CALL dbcsr_init_p(matrix_t_im)
3652 CALL kpoint_operator_get(op_ctx, local_kpoint, kpoint_spin_free, matrix_t_kp, &
3653 mat_re=matrix_t_re, mat_im=matrix_t_im, &
3654 matrix_row=kpoint_spin_free, slot=kpoint_slot_t)
3655 END IF
3656
3657 DO ispin = 1, nspin
3658 IF (prec_type == ot_precond_full_all .OR. prec_type == ot_precond_full_single_inverse) THEN
3659 CALL dbcsr_init_p(matrix_h_re)
3660 CALL dbcsr_init_p(matrix_h_im)
3661 CALL kpoint_operator_get(op_ctx, local_kpoint, ispin, matrix_ks_kp, &
3662 mat_re=matrix_h_re, mat_im=matrix_h_im, matrix_row=ispin, slot=kpoint_slot_ks)
3663 END IF
3664 IF (local_kpoint <= kplocal) THEN
3665 kp => kpoints%kp_env(local_kpoint)%kpoint_env
3666 channel = (local_kpoint - 1)*nspin + ispin
3667 IF (prec_type == ot_precond_full_all .OR. prec_type == ot_precond_full_single_inverse) THEN
3668 ! Davidson iterates the full MO set, so its active space is all nmo
3669 ! MOs of the channel: assemble the complex channel matrices on the
3670 ! kpoint-group context and call the complex H-based cores directly
3671 ! (mirrors the complex preconditioners of the K-point OT case; the
3672 ! FULL_SINGLE_INVERSE variant trades the state selectivity of
3673 ! FULL_ALL for one shifted factorization and a gemm application)
3674 CALL get_mo_set(kp%mos(ispin), eigenvalues=eigenvalues, &
3675 nmo=nmo)
3676 CALL cp_cfm_create(matrix_c, kp%mos(ispin)%cmo_coeff%matrix_struct, name="dav pc C")
3677 CALL cp_cfm_to_cfm(kp%mos(ispin)%cmo_coeff, matrix_c)
3678 CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, &
3679 scf_env%ot_preconditioner(channel)%preconditioner, &
3680 'dav pc H', matrix_h)
3681 CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, &
3682 scf_env%ot_preconditioner(channel)%preconditioner, &
3683 'dav pc S', matrix_s)
3684 IF (prec_type == ot_precond_full_all) THEN
3685 CALL make_complex_full_all(scf_env%ot_preconditioner(channel)%preconditioner, &
3686 matrix_c, matrix_h, matrix_s, eigenvalues(1:nmo), &
3687 energy_gap)
3688 ELSE
3689 CALL make_complex_full_single_inverse( &
3690 scf_env%ot_preconditioner(channel)%preconditioner, &
3691 matrix_c, matrix_h, matrix_s, energy_gap)
3692 END IF
3693 CALL cp_cfm_release(matrix_c)
3694 CALL cp_cfm_release(matrix_h)
3695 CALL cp_cfm_release(matrix_s)
3696 ELSE IF (prec_type == ot_precond_full_kinetic) THEN
3697 CALL make_preconditioner_complex_full_kinetic( &
3698 scf_env%ot_preconditioner(channel)%preconditioner, &
3699 matrix_t_re, matrix_t_im, matrix_s_re, matrix_s_im, energy_gap, solver_type)
3700 ELSE IF (prec_type == ot_precond_s_inverse) THEN
3701 CALL make_preconditioner_complex_full_s_inverse( &
3702 scf_env%ot_preconditioner(channel)%preconditioner, &
3703 matrix_s_re, matrix_s_im, solver_type)
3704 END IF
3705 END IF
3706 IF (prec_type == ot_precond_full_all .OR. prec_type == ot_precond_full_single_inverse) THEN
3707 CALL dbcsr_release_p(matrix_h_re)
3708 CALL dbcsr_release_p(matrix_h_im)
3709 END IF
3710 END DO
3711
3712 CALL dbcsr_release_p(matrix_s_re)
3713 CALL dbcsr_release_p(matrix_s_im)
3714 IF (prec_type == ot_precond_full_kinetic) THEN
3715 CALL dbcsr_release_p(matrix_t_re)
3716 CALL dbcsr_release_p(matrix_t_im)
3717 END IF
3718 END DO
3719
3720 CALL timestop(handle)
3721
3722 END SUBROUTINE davidson_kp_prepare_preconditioners
3723
3724! **************************************************************************************************
3725!> \brief Kpoint diagonalization routine
3726!> Transforms matrices to kpoint, distributes kpoint groups, performs diagonalization
3727!> \param matrix_s Overlap matrices (RS indices, global)
3728!> \param kpoints Kpoint environment
3729!> \param fmwork full matrices distributed over all groups
3730!> \param lowdin_forces retain the overlap spectrum for Lowdin force/stress derivatives
3731!> \par History
3732!> 02.2026 created [JGH]
3733! **************************************************************************************************
3734 SUBROUTINE diag_kp_smat(matrix_s, kpoints, fmwork, lowdin_forces)
3735
3736 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
3737 TYPE(kpoint_type), POINTER :: kpoints
3738 TYPE(cp_fm_type), DIMENSION(:) :: fmwork
3739 LOGICAL, INTENT(IN), OPTIONAL :: lowdin_forces
3740
3741 CHARACTER(len=*), PARAMETER :: routinen = 'diag_kp_smat'
3742 COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp, 0.0_dp), &
3743 czero = (0.0_dp, 0.0_dp)
3744
3745 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ceig
3746 INTEGER :: handle, ikp, kplocal, kpmax, nao
3747 INTEGER, DIMENSION(2) :: kp_range
3748 INTEGER, DIMENSION(:, :), POINTER :: kp_dist
3749 LOGICAL :: save_spectrum, use_real_wfn
3750 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
3751 TYPE(cp_cfm_type) :: csmat, cwork
3752 TYPE(cp_fm_pool_type), POINTER :: ao_ao_fm_pool
3753 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
3754 TYPE(cp_fm_type) :: fmlocal, rsmat
3755 TYPE(kpoint_env_type), POINTER :: kp
3756 TYPE(kpoint_operator_context_type) :: op_ctx
3757 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3758 POINTER :: sab_nl
3759 TYPE(qs_matrix_pools_type), POINTER :: mpools
3760
3761 CALL timeset(routinen, handle)
3762
3763 NULLIFY (sab_nl)
3764 CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn, kp_range=kp_range, &
3765 kp_dist=kp_dist, sab_nl=sab_nl)
3766 cpassert(ASSOCIATED(sab_nl))
3767 kplocal = kp_range(2) - kp_range(1) + 1
3768 kpmax = maxval(kp_dist(2, :) - kp_dist(1, :) + 1)
3769
3770 ! fm pools to be used within a kpoint group
3771 CALL get_kpoint_info(kpoints, mpools=mpools)
3772 CALL mpools_get(mpools, ao_ao_fm_pool=ao_ao_fm_pool)
3773
3774 CALL fm_pool_create_fm(ao_ao_fm_pool, fmlocal)
3775 CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
3776
3777 IF (use_real_wfn) THEN
3778 CALL cp_fm_create(rsmat, matrix_struct)
3779 ELSE
3780 CALL cp_cfm_create(csmat, matrix_struct)
3781 CALL cp_cfm_create(cwork, matrix_struct)
3782 END IF
3783
3784 CALL cp_fm_get_info(fmwork(1), nrow_global=nao)
3785 ALLOCATE (eigenvalues(nao), ceig(nao))
3786
3787 ! The overlap is spin free: the routine drives every instance under the
3788 ! single spin-free key.
3789 CALL kpoint_operator_context_create(op_ctx, kpoints, fmwork(1), 1)
3790 CALL kpoint_operator_prepare(op_ctx, matrix_s, kpoint_slot_s, matrix_row=kpoint_spin_free)
3791 DO ikp = 1, kpmax
3792 CALL kpoint_operator_start(op_ctx, ikp, kpoint_spin_free, matrix_s, kpoint_slot_s, &
3793 use_grid=.false., matrix_row=kpoint_spin_free)
3794 END DO
3795
3796 ! Finish communication then diagonalise in each group
3797 DO ikp = 1, kpmax
3798 IF (use_real_wfn) THEN
3799 CALL kpoint_operator_finish(op_ctx, ikp, kpoint_spin_free, kpoint_slot_s, fm_re=rsmat)
3800 ELSE
3801 CALL kpoint_operator_finish(op_ctx, ikp, kpoint_spin_free, kpoint_slot_s, csmat)
3802 END IF
3803
3804 ! Each kpoint group has now information on a kpoint to be diagonalized
3805 ! Eigensolver Hermite or Symmetric
3806 IF (ikp > kplocal) cycle
3807 kp => kpoints%kp_env(ikp)%kpoint_env
3808 save_spectrum = ALLOCATED(kp%lowdin_sqrt)
3809 IF (PRESENT(lowdin_forces)) save_spectrum = lowdin_forces
3810 CALL cp_cfm_release(kp%lowdin_evec)
3811 IF (ALLOCATED(kp%lowdin_sqrt)) DEALLOCATE (kp%lowdin_sqrt)
3812 IF (use_real_wfn) THEN
3813 CALL choose_eigv_solver(rsmat, fmlocal, eigenvalues)
3814 ELSE
3815 CALL cp_cfm_heevd(csmat, cwork, eigenvalues)
3816 END IF
3817 cpassert(all(eigenvalues(1:nao) >= 0.0_dp))
3818 IF (save_spectrum) THEN
3819 IF (any(eigenvalues <= 0.0_dp)) THEN
3820 cpabort("Lowdin forces require a positive definite overlap")
3821 END IF
3822 ALLOCATE (kp%lowdin_sqrt(nao))
3823 kp%lowdin_sqrt(:) = sqrt(eigenvalues)
3824 CALL cp_cfm_create(kp%lowdin_evec, matrix_struct)
3825 IF (use_real_wfn) THEN
3826 CALL cp_fm_to_cfm(fmlocal, mtarget=kp%lowdin_evec)
3827 ELSE
3828 CALL cp_cfm_to_cfm(cwork, kp%lowdin_evec)
3829 END IF
3830 END IF
3831 IF (use_real_wfn) THEN
3832 CALL cp_fm_release(kp%shalf)
3833 CALL cp_fm_create(kp%shalf, matrix_struct)
3834 eigenvalues(1:nao) = sqrt(eigenvalues(1:nao))
3835 CALL cp_fm_to_fm(fmlocal, rsmat)
3836 CALL cp_fm_column_scale(rsmat, eigenvalues)
3837 CALL parallel_gemm("N", "T", nao, nao, nao, 1.0_dp, rsmat, fmlocal, &
3838 0.0_dp, kp%shalf)
3839 ELSE
3840 CALL cp_cfm_release(kp%cshalf)
3841 CALL cp_cfm_create(kp%cshalf, matrix_struct)
3842 ceig(1:nao) = sqrt(eigenvalues(1:nao))
3843 CALL cp_cfm_to_cfm(cwork, csmat)
3844 CALL cp_cfm_column_scale(csmat, ceig)
3845 CALL parallel_gemm("N", "C", nao, nao, nao, cone, csmat, cwork, &
3846 czero, kp%cshalf)
3847 END IF
3848 END DO
3849
3850 CALL kpoint_operator_context_release(op_ctx)
3851
3852 ! All done
3853 DEALLOCATE (eigenvalues, ceig)
3854
3855 IF (use_real_wfn) THEN
3856 CALL cp_fm_release(rsmat)
3857 ELSE
3858 CALL cp_cfm_release(csmat)
3859 CALL cp_cfm_release(cwork)
3860 END IF
3861 CALL fm_pool_give_back_fm(ao_ao_fm_pool, fmlocal)
3862
3863 CALL timestop(handle)
3864
3865 END SUBROUTINE diag_kp_smat
3866
3867END MODULE qs_scf_diagonalization
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:126
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:57
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:96
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
qs_density_matrices::calculate_density_matrix
Definition qs_density_matrices.F:58
qs_mo_methods::calculate_subspace_eigenvalues
Definition qs_mo_methods.F:70
qs_mo_occupation::set_mo_occupation
Definition qs_mo_occupation.F:50
qs_scf_methods::combine_ks_matrices
Definition qs_scf_methods.F:72
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition cp_array_utils.F:18
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_scale_and_add
subroutine, public cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b).
Definition cp_cfm_basic_linalg.F:244
cp_cfm_basic_linalg::cp_cfm_transpose
subroutine, public cp_cfm_transpose(matrix, trans, matrixt)
Transposes a BLACS distributed complex matrix.
Definition cp_cfm_basic_linalg.F:1117
cp_cfm_basic_linalg::cp_cfm_scale_and_add_fm
subroutine, public cp_cfm_scale_and_add_fm(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b). where b is a real matrix (adapted from cp_cf...
Definition cp_cfm_basic_linalg.F:353
cp_cfm_basic_linalg::cp_cfm_triangular_multiply
subroutine, public cp_cfm_triangular_multiply(triangular_matrix, matrix_b, side, transa_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, alpha)
Multiplies in place by a triangular matrix: matrix_b = alpha op(triangular_matrix) matrix_b or (if si...
Definition cp_cfm_basic_linalg.F:985
cp_cfm_basic_linalg::cp_cfm_column_scale
subroutine, public cp_cfm_column_scale(matrix_a, scaling, first_col, source)
Scales columns of the full matrix by corresponding factors.
Definition cp_cfm_basic_linalg.F:635
cp_cfm_cholesky
various cholesky decomposition related routines
Definition cp_cfm_cholesky.F:13
cp_cfm_cholesky::cp_cfm_cholesky_decompose
subroutine, public cp_cfm_cholesky_decompose(matrix, n, info_out)
Used to replace a symmetric positive definite matrix M with its Cholesky decomposition U: M = U^T * U...
Definition cp_cfm_cholesky.F:53
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_geeig_canon
subroutine, public cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig, nmo_retained)
General Eigenvalue Problem AX = BXE Use canonical orthogonalization.
Definition cp_cfm_diag.F:686
cp_cfm_diag::cp_cfm_geeig_local
subroutine, public cp_cfm_geeig_local(amatrix, bmatrix, eigenvectors, eigenvalues, ws)
Solve a generalized complex eigenproblem using the local LAPACK backend. This routine is restricted t...
Definition cp_cfm_diag.F:777
cp_cfm_diag::cp_cfm_local_workspace_prepare
subroutine, public cp_cfm_local_workspace_prepare(ws, n, a, b)
Query LAPACK once per matrix order/solver and retain private scratch across k points.
Definition cp_cfm_diag.F:943
cp_cfm_diag::cp_cfm_geeig_canon_local
subroutine, public cp_cfm_geeig_canon_local(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig, ws)
Canonical generalized complex diagonalization on a one-rank BLACS grid.
Definition cp_cfm_diag.F:832
cp_cfm_diag::cp_cfm_heevd
subroutine, public cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
Perform a diagonalisation of a complex matrix.
Definition cp_cfm_diag.F:94
cp_cfm_diag::cp_cfm_geeig
subroutine, public cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work, lowest_subset)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition cp_cfm_diag.F:439
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:183
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:845
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:126
cp_cfm_types::cp_cfm_set_all
subroutine, public cp_cfm_set_all(matrix, alpha, beta)
Set all elements of the full matrix to alpha. Besides, set all diagonal matrix elements to beta (if g...
Definition cp_cfm_types.F:203
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:793
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_release_p
subroutine, public dbcsr_release_p(matrix)
...
Definition cp_dbcsr_api.F:222
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:940
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan, check_plan)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:344
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:240
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:171
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1585
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling, first_col, source)
Scales columns of the full matrix by corresponding factors.
Definition cp_fm_basic_linalg.F:1694
cp_fm_basic_linalg::cp_fm_symm
subroutine, public cp_fm_symm(side, uplo, m, n, alpha, matrix_a, matrix_b, beta, matrix_c)
computes matrix_c = beta * matrix_c + alpha * matrix_a * matrix_b computes matrix_c = beta * matrix_c...
Definition cp_fm_basic_linalg.F:569
cp_fm_cholesky
various cholesky decomposition related routines
Definition cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_restore
subroutine, public cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, pos, transa)
apply Cholesky decomposition op can be "SOLVE" (out = U^-1 * in) or "MULTIPLY" (out = U * in) pos can...
Definition cp_fm_cholesky.F:319
cp_fm_cholesky::cp_fm_cholesky_reduce
subroutine, public cp_fm_cholesky_reduce(matrix, matrixb, itype)
reduce a matrix pencil A,B to normal form B has to be cholesky decomposed with cp_fm_cholesky_decompo...
Definition cp_fm_cholesky.F:261
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::cp_fm_geeig_canon
subroutine, public cp_fm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig, nmo_retained)
General Eigenvalue Problem AX = BXE Use canonical diagonalization : U*s**(-1/2).
Definition cp_fm_diag.F:1718
cp_fm_diag::direct_generalized_diagonalization
logical, save, public direct_generalized_diagonalization
Definition cp_fm_diag.F:106
cp_fm_diag::cp_fm_geeig
subroutine, public cp_fm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE. Use cuSOLVERMp directly when requested and large enough; otherwi...
Definition cp_fm_diag.F:1503
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:262
cp_fm_diag::diag_check_requested
logical function, public diag_check_requested()
Return whether diagonalization checks should be performed.
Definition cp_fm_diag.F:316
cp_fm_diag::diag_lib_explicit
logical, save, public diag_lib_explicit
Definition cp_fm_diag.F:96
cp_fm_pool_types
pool for for elements that are retained and released
Definition cp_fm_pool_types.F:14
cp_fm_pool_types::fm_pool_create_fm
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
Definition cp_fm_pool_types.F:160
cp_fm_pool_types::fm_pool_give_back_fm
subroutine, public fm_pool_give_back_fm(pool, element)
returns the element to the pool
Definition cp_fm_pool_types.F:199
cp_fm_pool_types::fm_pool_get_el_struct
type(cp_fm_struct_type) function, pointer, public fm_pool_get_el_struct(pool)
returns the structure of the elements in this pool
Definition cp_fm_pool_types.F:224
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_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:1046
cp_fm_types::cp_fm_add_to_element
subroutine, public cp_fm_add_to_element(matrix, irow_global, icol_global, alpha)
...
Definition cp_fm_types.F:2141
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:167
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_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:583
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::debug_print_level
integer, parameter, public debug_print_level
Definition cp_output_handling.F:77
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1156
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout, file_access)
...
Definition cp_output_handling.F:935
dft_plus_u
Add the DFT+U contribution to the Hamiltonian matrix.
Definition dft_plus_u.F:18
dft_plus_u::lowdin_kp_apply
subroutine, public lowdin_kp_apply(kp, v, rmat, cmat)
Add B diag(v) B**H before any solver, DIIS residual, or history storage uses H(k).
Definition dft_plus_u.F:1267
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::core_guess
integer, parameter, public core_guess
Definition input_constants.F:224
input_constants::cholesky_restore
integer, parameter, public cholesky_restore
Definition input_constants.F:218
input_constants::cholesky_dbcsr
integer, parameter, public cholesky_dbcsr
Definition input_constants.F:218
input_constants::cholesky_off
integer, parameter, public cholesky_off
Definition input_constants.F:218
input_constants::ot_precond_full_kinetic
integer, parameter, public ot_precond_full_kinetic
Definition input_constants.F:551
input_constants::cholesky_reduce
integer, parameter, public cholesky_reduce
Definition input_constants.F:218
input_constants::high_spin_roks
integer, parameter, public high_spin_roks
Definition input_constants.F:844
input_constants::cholesky_inverse
integer, parameter, public cholesky_inverse
Definition input_constants.F:218
input_constants::ot_precond_none
integer, parameter, public ot_precond_none
Definition input_constants.F:551
input_constants::ot_precond_full_single_inverse
integer, parameter, public ot_precond_full_single_inverse
Definition input_constants.F:551
input_constants::diag_update_method_adiis
integer, parameter, public diag_update_method_adiis
Definition input_constants.F:982
input_constants::general_roks
integer, parameter, public general_roks
Definition input_constants.F:844
input_constants::restart_guess
integer, parameter, public restart_guess
Definition input_constants.F:224
input_constants::ot_precond_s_inverse
integer, parameter, public ot_precond_s_inverse
Definition input_constants.F:551
input_constants::ot_precond_full_all
integer, parameter, public ot_precond_full_all
Definition input_constants.F:551
input_cp2k_kpoints
function that build the kpoints section of the input
Definition input_cp2k_kpoints.F:14
input_cp2k_kpoints::lattice_fft_auto
integer, parameter, public lattice_fft_auto
Definition input_cp2k_kpoints.F:37
input_cp2k_kpoints::lattice_fft_off
integer, parameter, public lattice_fft_off
Definition input_cp2k_kpoints.F:37
input_cp2k_kpoints::use_real_wfn
integer, parameter, public use_real_wfn
Definition input_cp2k_kpoints.F:37
input_cp2k_kpoints::lattice_fft_on
integer, parameter, public lattice_fft_on
Definition input_cp2k_kpoints.F:37
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
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kpoint_methods
Routines needed for kpoint calculation.
Definition kpoint_methods.F:15
kpoint_methods::kpoint_density_transform
subroutine, public kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwork, for_aux_fit, pmat_ext, overlap_rs, mixed_aux_orb)
generate real space density matrices in DBCSR format
Definition kpoint_methods.F:2253
kpoint_methods::kpoint_set_mo_occupation
subroutine, public kpoint_set_mo_occupation(kpoint, smear, probe, added_mos_auto, added_mos_auto_grow, separate_spin_occupations)
Given the eigenvalues of all kpoints, calculates the occupation numbers.
Definition kpoint_methods.F:1535
kpoint_methods::rskp_transform_grid_release
subroutine, public rskp_transform_grid_release(grid)
Release a reciprocal-grid transformation cache.
Definition kpoint_methods.F:1512
kpoint_methods::rskp_transform_grid_prepare
subroutine, public rskp_transform_grid_prepare(grid, rmatrix, rsmat, ispin, xkp, nkp_grid, cell_to_index, sab_nl, used_fft, is_complex, rs_sign, max_storage_bytes)
Prepare a batched real-cell transform to a regular grid or subset for local DBCSR blocks....
Definition kpoint_methods.F:1271
kpoint_methods::kpoint_density_matrices
subroutine, public kpoint_density_matrices(kpoint, energy_weighted, for_aux_fit)
Calculate kpoint density matrices (rho(k), owned by kpoint groups).
Definition kpoint_methods.F:1891
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:613
local_gemm_api
Definition local_gemm_api.F:8
local_gemm_api::local_gemm_pu_gpu
integer, parameter, public local_gemm_pu_gpu
Definition local_gemm_api.F:41
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:127
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:144
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:144
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:144
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
preconditioner_makes
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition preconditioner_makes.F:15
preconditioner_makes::make_complex_full_single_inverse
subroutine, public make_complex_full_single_inverse(preconditioner_env, matrix_c0, matrix_h, matrix_s, energy_gap)
Build a gauge-covariant FULL_SINGLE_INVERSE operator for a complex k-point channel.
Definition preconditioner_makes.F:1172
preconditioner_makes::make_complex_full_all
subroutine, public make_complex_full_all(preconditioner_env, matrix_c0, matrix_h, matrix_s, c0_evals, energy_gap)
Build the state-selective FULL_ALL operator for a complex k-point channel. The occupied/reference sub...
Definition preconditioner_makes.F:106
preconditioner_types
types of preconditioners
Definition preconditioner_types.F:14
preconditioner_types::init_preconditioner
subroutine, public init_preconditioner(preconditioner_env, para_env, blacs_env)
...
Definition preconditioner_types.F:137
preconditioner_types::destroy_preconditioner
subroutine, public destroy_preconditioner(preconditioner_env)
...
Definition preconditioner_types.F:193
preconditioner
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition preconditioner.F:15
preconditioner::make_preconditioner_complex_full_s_inverse
subroutine, public make_preconditioner_complex_full_s_inverse(preconditioner_env, matrix_s_re, matrix_s_im, solver_type)
Construct a complex FULL_S_INVERSE preconditioner.
Definition preconditioner.F:1334
preconditioner::dbcsr_pair_to_cfm
subroutine, public dbcsr_pair_to_cfm(matrix_re, matrix_im, preconditioner_env, matrix_name, matrix)
Copy a real/imaginary DBCSR pair to one distributed complex full matrix.
Definition preconditioner.F:1400
preconditioner::restart_preconditioner
subroutine, public restart_preconditioner(qs_env, preconditioner, prec_type, nspins)
Allows for a restart of the preconditioner depending on the method it purges all arrays or keeps them...
Definition preconditioner.F:1450
preconditioner::prepare_preconditioner
subroutine, public prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, ot_preconditioner, prec_type, solver_type, energy_gap, nspins, has_unit_metric, convert_to_dbcsr, chol_type, full_mo_set, chebyshev_degree, low_rank_base, fermi_low_rank_max_rank, lattice_fft, lattice_fft_local_cells)
...
Definition preconditioner.F:1526
preconditioner::make_preconditioner_complex_full_kinetic
subroutine, public make_preconditioner_complex_full_kinetic(preconditioner_env, matrix_t_re, matrix_t_im, matrix_s_re, matrix_s_im, energy_gap, solver_type)
Construct a complex FULL_KINETIC preconditioner.
Definition preconditioner.F:1367
qs_density_matrices
collects routines that calculate density matrices
Definition qs_density_matrices.F:14
qs_density_mixing_types
module that contains the definitions of the scf types
Definition qs_density_mixing_types.F:14
qs_density_mixing_types::direct_mixing_nr
integer, parameter, public direct_mixing_nr
Definition qs_density_mixing_types.F:50
qs_density_mixing_types::gspace_mixing_nr
integer, parameter, public gspace_mixing_nr
Definition qs_density_mixing_types.F:50
qs_diis
Apply the direct inversion in the iterative subspace (DIIS) of Pulay in the framework of an SCF itera...
Definition qs_diis.F:21
qs_diis::qs_diis_b_info_kp
subroutine, public qs_diis_b_info_kp(diis_buffer, ib, nb)
Update info about the current buffer step ib and the current number of buffers nb.
Definition qs_diis.F:998
qs_diis::qs_diis_b_calc_err_kp
subroutine, public qs_diis_b_calc_err_kp(diis_buffer, ib, mos, kc, sc, ispin, ikp, cmos, real_wfn, local, gemm_ctx)
Calculate and store the error for a given k-point.
Definition qs_diis.F:1024
qs_diis::qs_diis_b_clear
pure subroutine, public qs_diis_b_clear(diis_buffer)
clears the buffer
Definition qs_diis.F:527
qs_diis::qs_diis_b_step_kp
subroutine, public qs_diis_b_step_kp(diis_buffer, coeffs, ib, nb, delta, error_max, diis_step, eps_diis, nspin, nkp, nkp_local, nmixing, scf_section, para_env_inter_kp)
Update the SCF DIIS buffer, and if appropriate does a diis step, for k-points.
Definition qs_diis.F:1128
qs_diis::qs_diis_b_step
subroutine, public qs_diis_b_step(diis_buffer, mo_array, kc, sc, delta, error_max, diis_step, eps_diis, nmixing, s_matrix, scf_section, roks)
Update the SCF DIIS buffer, and if appropriate does a diis step.
Definition qs_diis.F:228
qs_diis::qs_diis_b_clear_kp
pure subroutine, public qs_diis_b_clear_kp(diis_buffer)
clears the buffer
Definition qs_diis.F:984
qs_diis::qs_diis_b_check_i_alloc_kp
subroutine, public qs_diis_b_check_i_alloc_kp(diis_buffer, matrix_struct, nspin, nkp, scf_section)
Allocate and initialize a DIIS buffer for nao*nao parameter variables and with a buffer size of nbuff...
Definition qs_diis.F:901
qs_energy_types
Definition qs_energy_types.F:14
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:544
qs_gspace_mixing
Definition qs_gspace_mixing.F:9
qs_gspace_mixing::gspace_mixing
subroutine, public gspace_mixing(qs_env, mixing_method, mixing_store, rho, para_env, iter_count, auxiliary)
Driver for the g-space mixing, calls the proper routine given the requested method.
Definition qs_gspace_mixing.F:68
qs_kpoint_operators
Assembly of real and complex k-point operators from real-space DBCSR matrices. Outputs use the distri...
Definition qs_kpoint_operators.F:26
qs_kpoint_operators::kpoint_operator_context_create
subroutine, public kpoint_operator_context_create(ctx, kpoints, ao_ao_fm, nspin)
Create an empty operator context. The context owns no buffers after this call; they are created on fi...
Definition qs_kpoint_operators.F:390
qs_kpoint_operators::kpoint_operator_prepare
subroutine, public kpoint_operator_prepare(ctx, rsmat, slot, prepared, matrix_row, target)
Prepare an immutable operator slot for group-local assembly. Refresh after the source values change,...
Definition qs_kpoint_operators.F:574
qs_kpoint_operators::kpoint_operator_finish
subroutine, public kpoint_operator_finish(ctx, ikp, ispin, slot, cmat, fm_re)
Finish one started instance: build its prepared group-local operator, or merge received parts into cm...
Definition qs_kpoint_operators.F:1770
qs_kpoint_operators::kpoint_operator_store
subroutine, public kpoint_operator_store(ctx, matrix_ks, matrix_s, matrix_t)
Prepare owned split-FM operators for OT without entering an SCF eigensolver. H (or the initialization...
Definition qs_kpoint_operators.F:2394
qs_kpoint_operators::kpoint_operator_work_size
subroutine, public kpoint_operator_work_size(ctx, nentry)
Size per-worker phase scratch for direct assembly into a solver matrix.
Definition qs_kpoint_operators.F:1606
qs_kpoint_operators::kpoint_operator_get
subroutine, public kpoint_operator_get(ctx, ikp, ispin, rsmat, cmat, fm_re, fm_im, mat_re, mat_im, matrix_row, slot, grid)
Assemble one operator for one global round and spin in one call: run the single-group transfer of the...
Definition qs_kpoint_operators.F:1868
qs_kpoint_operators::kpoint_spin_free
integer, parameter, public kpoint_spin_free
Row of rsmat that carries the spin-free operator image. S and T carry no spin copies: drivers launch ...
Definition qs_kpoint_operators.F:101
qs_kpoint_operators::kpoint_operator_batch_size
subroutine, public kpoint_operator_batch_size(ctx, nworkers, nmatrices, width, work_size)
Size local k tiles, budgeting operator outputs and block scratch over all workers....
Definition qs_kpoint_operators.F:2214
qs_kpoint_operators::kpoint_operator_evaluate
subroutine, public kpoint_operator_evaluate(ctx, ispin, slot, ik, phase, target, fm_re, fm_im)
Evaluate an immutable prepared operator with exclusively caller-owned output and scratch....
Definition qs_kpoint_operators.F:1627
qs_kpoint_operators::kpoint_slot_s
integer, parameter, public kpoint_slot_s
Definition qs_kpoint_operators.F:94
qs_kpoint_operators::kpoint_slot_t
integer, parameter, public kpoint_slot_t
Definition qs_kpoint_operators.F:94
qs_kpoint_operators::kpoint_slot_ks
integer, parameter, public kpoint_slot_ks
Definition qs_kpoint_operators.F:94
qs_kpoint_operators::kpoint_operator_get_local
subroutine, public kpoint_operator_get_local(matrix_rs, kpoints, kp, ispin, cache_re, cache_im, matrix_re, matrix_im)
Read an operator inside a group-local OT loop, after collective cache preparation....
Definition qs_kpoint_operators.F:1965
qs_kpoint_operators::kpoint_operator_evaluate_batch
subroutine, public kpoint_operator_evaluate_batch(ctx, ispin, slot, first_k, batch, work)
Assemble a caller-owned complex k tile using real block/image contractions. The prepared route and sn...
Definition qs_kpoint_operators.F:2261
qs_kpoint_operators::kpoint_operator_start
subroutine, public kpoint_operator_start(ctx, ikp, ispin, rsmat, slot, grid, use_grid, matrix_row)
Start one operator for one global round and spin. The service walks all groups with a k point in this...
Definition qs_kpoint_operators.F:1684
qs_kpoint_operators::kpoint_operator_context_release
subroutine, public kpoint_operator_context_release(ctx)
Release the context together with every buffer it created. Call once, at the end of the fixed-layout ...
Definition qs_kpoint_operators.F:427
qs_ks_methods
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition qs_ks_methods.F:22
qs_ks_methods::qs_ks_update_qs_env
subroutine, public qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, print_active)
updates the Kohn Sham matrix of the given qs_env (facility method)
Definition qs_ks_methods.F:1657
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::qs_ks_did_change
subroutine, public qs_ks_did_change(ks_env, s_mstruct_changed, rho_changed, potential_changed, full_reset)
tells that some of the things relevant to the ks calculation did change. has to be called when change...
Definition qs_ks_types.F:977
qs_matrix_pools
wrapper for the pools of matrixes
Definition qs_matrix_pools.F:14
qs_matrix_pools::mpools_get
subroutine, public mpools_get(mpools, ao_mo_fm_pools, ao_ao_fm_pool, mo_mo_fm_pools, ao_mosub_fm_pools, mosub_mosub_fm_pools, maxao_maxmo_fm_pool, maxao_maxao_fm_pool, maxmo_maxmo_fm_pool)
returns various attributes of the mpools (notably the pools contained in it)
Definition qs_matrix_pools.F:118
qs_mixing_utils
Definition qs_mixing_utils.F:9
qs_mixing_utils::charge_mixing_init
elemental subroutine, public charge_mixing_init(mixing_store)
initialiation needed when charge mixing is used
Definition qs_mixing_utils.F:756
qs_mixing_utils::mixing_init
subroutine, public mixing_init(mixing_method, rho, mixing_store, para_env, rho_atom, auxiliary)
initialiation needed when gspace mixing is used
Definition qs_mixing_utils.F:462
qs_mixing_utils::mixing_allocate
subroutine, public mixing_allocate(qs_env, mixing_method, p_mix_new, p_delta, nspins, mixing_store)
allocation needed when density mixing is used
Definition qs_mixing_utils.F:117
qs_mixing_utils::self_consistency_check
subroutine, public self_consistency_check(rho_ao, p_delta, para_env, p_out, delta)
...
Definition qs_mixing_utils.F:64
qs_mo_methods
collects routines that perform operations directly related to MOs
Definition qs_mo_methods.F:14
qs_mo_occupation
Set occupation of molecular orbitals.
Definition qs_mo_occupation.F:15
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, cmo_coeff)
Get the components of a MO set data structure.
Definition qs_mo_types.F:440
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_ot_eigensolver
an eigen-space solver for the generalised symmetric eigenvalue problem for sparse matrices,...
Definition qs_ot_eigensolver.F:13
qs_ot_eigensolver::ot_eigensolver_complex
subroutine, public ot_eigensolver_complex(matrix_h, matrix_h_im, matrix_s, matrix_s_im, matrix_c_cfm, preconditioner, eps_gradient, iter_max, eigenvalues, silent, ot_settings, matrix_t, matrix_t_im, mo_set)
solve a fixed complex Hermitian generalized eigenproblem by OT
Definition qs_ot_eigensolver.F:403
qs_ot_eigensolver::ot_eigensolver
subroutine, public ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, matrix_c_fm, preconditioner, eps_gradient, iter_max, size_ortho_space, silent, ot_settings)
...
Definition qs_ot_eigensolver.F:111
qs_rho_atom_types
Definition qs_rho_atom_types.F:8
qs_rho_methods
methods of the rho structure (defined in qs_rho_types)
Definition qs_rho_methods.F:15
qs_rho_methods::qs_rho_update_rho
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
Definition qs_rho_methods.F:381
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_scf_block_davidson
module that contains the algorithms to perform an iterative diagonalization by the block-Davidson app...
Definition qs_scf_block_davidson.F:18
qs_scf_block_davidson::generate_extended_space_c
subroutine, public generate_extended_space_c(bdav_env, mos, matrix_h, matrix_s, output_unit, eps_iter, eps_iter_empty, preconditioner)
iterative diagonalization by the block-Davidson approach for one complex K point; complex counterpart...
Definition qs_scf_block_davidson.F:501
qs_scf_block_davidson::generate_extended_space_sparse
subroutine, public generate_extended_space_sparse(bdav_env, mo_set, matrix_h, matrix_s, output_unit, preconditioner)
...
Definition qs_scf_block_davidson.F:1069
qs_scf_block_davidson::generate_extended_space
subroutine, public generate_extended_space(bdav_env, mo_set, matrix_h, matrix_s, output_unit, preconditioner)
...
Definition qs_scf_block_davidson.F:98
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, added_mos_auto_grow, potential)
Kpoint diagonalization routine Transforms matrices to kpoint, distributes kpoint groups,...
Definition qs_scf_diagonalization.F:708
qs_scf_diagonalization::general_eigenproblem
subroutine, public general_eigenproblem(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step)
the inner loop of scf, specific to diagonalization with S matrix basically, in goes the ks matrix out...
Definition qs_scf_diagonalization.F:322
qs_scf_diagonalization::diag_kp_smat
subroutine, public diag_kp_smat(matrix_s, kpoints, fmwork, lowdin_forces)
Kpoint diagonalization routine Transforms matrices to kpoint, distributes kpoint groups,...
Definition qs_scf_diagonalization.F:3735
qs_scf_diagonalization::do_block_davidson_diag_kp
subroutine, public do_block_davidson_diag_kp(qs_env, scf_env, matrix_ks_kp, matrix_s_kp, kpoints, scf_control, added_mos_auto_grow)
block-Davidson SCF step for complex K points, one channel per (local kpoint, spin); mirrors do_block_...
Definition qs_scf_diagonalization.F:2955
qs_scf_diagonalization::diag_subspace_allocate
subroutine, public diag_subspace_allocate(subspace_env, qs_env, mos)
...
Definition qs_scf_diagonalization.F:1919
qs_scf_diagonalization::do_ot_diag
subroutine, public do_ot_diag(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step)
the inner loop of scf, specific to iterative diagonalization using OT with S matrix; basically,...
Definition qs_scf_diagonalization.F:2092
qs_scf_diagonalization::do_block_davidson_diag
subroutine, public do_block_davidson_diag(qs_env, scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, check_moconv_only)
iterative diagonalization using the block davidson space approach
Definition qs_scf_diagonalization.F:2836
qs_scf_diagonalization::do_roks_diag
subroutine, public do_roks_diag(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step, orthogonal_basis)
Solve a set restricted open Kohn-Sham (ROKS) equations based on the alpha and beta Kohn-Sham matrices...
Definition qs_scf_diagonalization.F:2354
qs_scf_diagonalization::diag_kp_basic
subroutine, public diag_kp_basic(matrix_ks, matrix_s, kpoints, fmwork)
Kpoint diagonalization routine Transforms matrices to kpoint, distributes kpoint groups,...
Definition qs_scf_diagonalization.F:1581
qs_scf_diagonalization::do_scf_diag_subspace
subroutine, public do_scf_diag_subspace(qs_env, scf_env, subspace_env, mos, rho, ks_env, scf_section, scf_control)
inner loop within MOS subspace, to refine occupation and density, before next diagonalization of the ...
Definition qs_scf_diagonalization.F:1705
qs_scf_diagonalization::do_block_krylov_diag_kp
subroutine, public do_block_krylov_diag_kp(qs_env, scf_env, matrix_ks_kp, matrix_s_kp, kpoints, scf_control, check_moconv_only, added_mos_auto_grow)
Block Krylov (Lanczos) refinement of the k-point MOS, the complex counterpart of do_block_krylov_diag...
Definition qs_scf_diagonalization.F:3184
qs_scf_diagonalization::adiis_candidate_shift
real(kind=dp) function adiis_candidate_shift(scf_env)
Fade candidate regularization through handover; never regularize a raw-map check.
Definition qs_scf_diagonalization.F:219
qs_scf_diagonalization::do_block_krylov_diag
subroutine, public do_block_krylov_diag(scf_env, mos, matrix_ks, scf_control, scf_section, check_moconv_only)
iterative diagonalization using the block Krylov-space approach
Definition qs_scf_diagonalization.F:2628
qs_scf_diagonalization::do_special_diag
subroutine, public do_special_diag(scf_env, mos, matrix_ks, scf_control, scf_section, diis_step)
the inner loop of scf, specific to diagonalization without S matrix basically, in goes the ks matrix ...
Definition qs_scf_diagonalization.F:1989
qs_scf_diagonalization::do_ot_diag_kp
subroutine, public do_ot_diag_kp(matrix_ks, matrix_s, matrix_t, kpoints, scf_env, scf_control, diis_step, added_mos_auto_grow)
fixed-H OT diagonalization for real or complex K-point channels
Definition qs_scf_diagonalization.F:2192
qs_scf_diagonalization::do_general_diag
subroutine, public do_general_diag(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step, probe)
...
Definition qs_scf_diagonalization.F:636
qs_scf_lanczos_c
Complex block-Lanczos refinement kernel for k-point SCF diagonalization. Processes unconverged column...
Definition qs_scf_lanczos_c.F:18
qs_scf_lanczos_c::lanczos_channel_end
subroutine, public lanczos_channel_end(c0, ws)
Defensive full-width Cholesky orthonormalization at the end of a channel's refinement sweeps,...
Definition qs_scf_lanczos_c.F:310
qs_scf_lanczos_c::lanczos_workspace_release
subroutine, public lanczos_workspace_release(ws)
Release the kernel workspace. Safe on a workspace that was never created: the matrix components defau...
Definition qs_scf_lanczos_c.F:897
qs_scf_lanczos_c::lanczos_channel_begin
subroutine, public lanczos_channel_begin(ks, c0, eval, eps_col, ws, max_res_norm, min_res_norm, nmo_conv, nmo_nc, check_moconv_only)
Entry preparation of one (k-point, spin) channel before its inner Lanczos sweeps: orthonormalizes the...
Definition qs_scf_lanczos_c.F:169
qs_scf_lanczos_c::lanczos_sweep_suffix
subroutine, public lanczos_sweep_suffix(ks, c0, eval, eps_col, nkrylov, nblock, first_sweep, need_suffix_rr, ws, max_res_norm, min_res_norm, nmo_conv, nmo_nc, step_refined, next_need_suffix_rr)
Run one inner Lanczos sweep on the unconverged suffix of c0. When the previous sweep ran in multi-blo...
Definition qs_scf_lanczos_c.F:246
qs_scf_lanczos_c::lanczos_workspace_create
subroutine, public lanczos_workspace_create(ws, mo_struct, ao_struct, nkrylov, nblock)
Create the kernel workspace at the worst-case widths of the run. The basis budget holds the larger of...
Definition qs_scf_lanczos_c.F:820
qs_scf_lanczos
module that contains the algorithms to perform an iterative diagonalization by the block-Lanczos appr...
Definition qs_scf_lanczos.F:15
qs_scf_lanczos::lanczos_refinement
subroutine, public lanczos_refinement(krylov_space, ks, c0, c1, eval, nao, eps_iter, ispin, check_moconv_only)
lanczos refinement by blocks of non-converged MOs
Definition qs_scf_lanczos.F:184
qs_scf_lanczos::lanczos_refinement_2v
subroutine, public lanczos_refinement_2v(krylov_space, ks, c0, c1, eval, nao, eps_iter, ispin, check_moconv_only)
...
Definition qs_scf_lanczos.F:521
qs_scf_methods
groups fairly general SCF methods, so that modules other than qs_scf can use them too split off from ...
Definition qs_scf_methods.F:19
qs_scf_methods::eigensolver_simple
subroutine, public eigensolver_simple(matrix_ks, mo_set, work, do_level_shift, level_shift, use_jacobi, jacobi_threshold)
...
Definition qs_scf_methods.F:480
qs_scf_methods::eigensolver_dbcsr
subroutine, public eigensolver_dbcsr(matrix_ks, matrix_ks_fm, mo_set, ortho_dbcsr, ksbuf1, ksbuf2)
...
Definition qs_scf_methods.F:315
qs_scf_methods::scf_env_density_mixing
subroutine, public scf_env_density_mixing(p_mix_new, mixing_store, rho_ao, para_env, iter_delta, iter_count, diis, invert)
perform (if requested) a density mixing
Definition qs_scf_methods.F:96
qs_scf_methods::eigensolver
subroutine, public eigensolver(matrix_ks_fm, mo_set, ortho, work, cholesky_method, do_level_shift, level_shift, matrix_u_fm, use_jacobi)
Diagonalise the Kohn-Sham matrix to get a new set of MO eigen- vectors and MO eigenvalues....
Definition qs_scf_methods.F:167
qs_scf_methods::eigensolver_symm
subroutine, public eigensolver_symm(matrix_ks_fm, mo_set, ortho, work, do_level_shift, level_shift, matrix_u_fm, use_jacobi, jacobi_threshold, ortho_red, work_red, matrix_ks_fm_red, matrix_u_fm_red)
...
Definition qs_scf_methods.F:381
qs_scf_methods::eigensolver_generalized
subroutine, public eigensolver_generalized(matrix_ks_fm, matrix_s, mo_set, work)
Solve the generalized eigenvalue problem.
Definition qs_scf_methods.F:271
qs_scf_types
module that contains the definitions of the scf types
Definition qs_scf_types.F:14
scf_control_types
parameters that control an scf iteration
Definition scf_control_types.F:17
cp_array_utils::cp_1d_r_p_type
represent a pointer to a 1d array
Definition cp_array_utils.F:99
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:55
cp_cfm_diag::cp_cfm_diag_workspace_type
Definition cp_cfm_diag.F:72
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:70
cp_control_types::dft_control_type
Definition cp_control_types.F:786
cp_control_types::hairy_probes_type
Definition cp_control_types.F:49
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_pool_types::cp_fm_pool_type
represent a pool of elements with the same structure
Definition cp_fm_pool_types.F:55
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:116
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_methods::rskp_grid_type
Definition kpoint_methods.F:170
kpoint_types::kpoint_env_type
Keeps information about a specific k-point.
Definition kpoint_types.F:94
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:196
local_gemm_api::local_gemm_ctxt_type
Definition local_gemm_api.F:47
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:774
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:223
qs_kpoint_operators::kpoint_operator_context_type
Caller-owned assembly context for a fixed geometry and parallel layout. OT retains it for one inner S...
Definition qs_kpoint_operators.F:167
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_matrix_pools::qs_matrix_pools_type
container for the pools of matrixes used by qs
Definition qs_matrix_pools.F:64
qs_mo_types::mo_set_type
Definition qs_mo_types.F:53
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_atom_types::rho_atom_type
Definition qs_rho_atom_types.F:23
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_scf_lanczos_c::lanczos_workspace_type
Caller-owned scratch of the refinement kernel: every matrix and array that the sweeps of one channel ...
Definition qs_scf_lanczos_c.F:62
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:120
qs_scf_types::subspace_env_type
Definition qs_scf_types.F:107
scf_control_types::scf_control_type
Definition scf_control_types.F:135