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