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