Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 : ! **************************************************************************************************
15 : MODULE 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 :
189 : CONTAINS
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 8826 : 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 8826 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
216 8826 : REAL(dp), DIMENSION(:, :), POINTER :: xkp
217 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
218 8826 : POINTER :: sab_nl
219 :
220 8826 : NULLIFY (cell_to_index, sab_nl, xkp)
221 0 : CPASSERT(ASSOCIATED(matrix_rs))
222 8826 : CPASSERT(ispin >= 1 .AND. ispin <= SIZE(matrix_rs, 1))
223 8826 : CPASSERT(ASSOCIATED(kpoints))
224 8826 : CPASSERT(ASSOCIATED(kp))
225 :
226 8826 : CALL dbcsr_init_p(matrix_re)
227 8826 : CALL dbcsr_init_p(matrix_im)
228 : CALL get_kpoint_info(kpoints, nkp_groups=nkp_groups, xkp=xkp, &
229 8826 : cell_to_index=cell_to_index, sab_nl=sab_nl)
230 8826 : IF (nkp_groups == 1) THEN
231 6352 : CPASSERT(ASSOCIATED(xkp))
232 6352 : CPASSERT(ASSOCIATED(cell_to_index))
233 6352 : CPASSERT(ASSOCIATED(sab_nl))
234 : CALL dbcsr_create(matrix_re, template=matrix_rs(1, 1)%matrix, &
235 6352 : matrix_type=dbcsr_type_symmetric)
236 : CALL dbcsr_create(matrix_im, template=matrix_rs(1, 1)%matrix, &
237 6352 : matrix_type=dbcsr_type_antisymmetric)
238 6352 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_re, sab_nl)
239 6352 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_im, sab_nl)
240 6352 : CALL dbcsr_set(matrix_re, 0.0_dp)
241 6352 : CALL dbcsr_set(matrix_im, 0.0_dp)
242 : CALL rskp_transform(rmatrix=matrix_re, cmatrix=matrix_im, rsmat=matrix_rs, ispin=ispin, &
243 6352 : xkp=xkp(1:3, kp%nkpoint), cell_to_index=cell_to_index, sab_nl=sab_nl)
244 : ELSE
245 2474 : CALL copy_fm_to_dbcsr_bc(cache_re, matrix_re)
246 2474 : CALL copy_fm_to_dbcsr_bc(cache_im, matrix_im)
247 : END IF
248 8826 : 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 99657 : 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 99657 : nspin = SIZE(matrix_ks)
286 99657 : NULLIFY (effective_ks_dbcsr, ortho, ortho_dbcsr)
287 :
288 99657 : use_adiis = scf_control%diagonalization%update_method == diag_update_method_adiis
289 99657 : use_subspace_fock = .FALSE.
290 99657 : IF (use_adiis) THEN
291 98 : IF (.NOT. ASSOCIATED(scf_env%scf_subspace_buffer)) THEN
292 0 : CPABORT("ADIIS is only available in the main Quickstep SCF workflow")
293 : END IF
294 98 : use_subspace_fock = scf_env%scf_subspace_buffer%use_combined_fock
295 : END IF
296 :
297 99657 : eps_diis = scf_control%eps_diis
298 99657 : use_cdiis = .FALSE.
299 99657 : 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 98 : scf_env%scf_subspace_buffer%diis_state_valid .AND. &
303 : scf_env%scf_subspace_buffer%nstored > 0) THEN
304 70 : use_cdiis = scf_env%iter_delta < eps_diis .OR. scf_env%scf_diis_buffer%ncall > 0
305 : END IF
306 : END IF
307 :
308 216654 : DO ispin = 1, nspin
309 216654 : IF (use_cdiis) THEN
310 96 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, scf_env%scf_work1(ispin))
311 116901 : ELSE IF (use_subspace_fock) THEN
312 : CALL copy_dbcsr_to_fm(scf_env%scf_subspace_buffer%combined_fock(ispin, 1)%matrix, &
313 32 : scf_env%scf_work1(ispin))
314 : ELSE
315 116869 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, scf_env%scf_work1(ispin))
316 : END IF
317 : END DO
318 :
319 99657 : IF (use_adiis) THEN
320 98 : diis_step = .FALSE.
321 98 : diis_error = scf_env%scf_subspace_buffer%last_old_fock_weight
322 98 : diis_weight = scf_env%scf_subspace_buffer%diis_weight
323 98 : IF (use_cdiis) THEN
324 48 : 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 48 : eps_diis, 2, s_matrix=matrix_s, scf_section=scf_section)
328 48 : IF (.NOT. ieee_is_finite(cdiis_error) .OR. cdiis_error >= eps_diis) THEN
329 0 : cdiis_step = .FALSE.
330 0 : IF (cdiis_ncall < scf_env%scf_diis_buffer%nbuffer) THEN
331 : ! The rejected state has not overwritten a full circular buffer.
332 0 : scf_env%scf_diis_buffer%ncall = cdiis_ncall
333 : ELSE
334 0 : CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
335 0 : diis_weight = 0.0_dp
336 : END IF
337 : END IF
338 48 : IF (cdiis_step) THEN
339 : ! Use one half-weight handover step; MAX_DIIS only controls the buffer capacity.
340 34 : 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 14 : diis_weight = 0.0_dp
344 : END IF
345 48 : IF (diis_weight < 1.0_dp) THEN
346 60 : DO ispin = 1, nspin
347 40 : IF (use_subspace_fock) THEN
348 : CALL copy_dbcsr_to_fm( &
349 24 : scf_env%scf_subspace_buffer%combined_fock(ispin, 1)%matrix, scf_env%scf_work2)
350 : ELSE
351 16 : 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 60 : 1.0_dp - diis_weight, scf_env%scf_work2)
355 : END DO
356 : END IF
357 48 : diis_step = cdiis_step .AND. diis_weight >= 1.0_dp
358 48 : IF (diis_step) diis_error = cdiis_error
359 50 : ELSE IF (ASSOCIATED(scf_env%scf_diis_buffer)) THEN
360 50 : IF (scf_env%scf_diis_buffer%ncall > 0) CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
361 50 : diis_weight = 0.0_dp
362 : END IF
363 98 : scf_env%scf_subspace_buffer%diis_weight = diis_weight
364 : ELSE
365 99559 : 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 80973 : scf_section=scf_section)
371 : ELSE
372 18586 : 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 99657 : (scf_env%iter_count > 1)))
379 :
380 99657 : IF (.NOT. use_adiis .AND. (scf_env%iter_count > 1) .AND. &
381 : (scf_env%iter_delta < scf_control%diagonalization%eps_jacobi)) THEN
382 0 : use_jacobi = .TRUE.
383 : ELSE
384 99559 : use_jacobi = .FALSE.
385 : END IF
386 :
387 98 : IF (use_adiis .AND. .NOT. diis_step) THEN
388 70 : scf_env%iter_param = diis_error
389 70 : IF (scf_env%scf_subspace_buffer%last_restart) THEN
390 0 : scf_env%iter_method = "ADIIS/Rst."
391 : ELSE
392 70 : scf_env%iter_method = "ADIIS/Diag."
393 : END IF
394 99587 : ELSE IF (diis_step) THEN
395 49733 : scf_env%iter_param = diis_error
396 49733 : IF (use_jacobi) THEN
397 0 : scf_env%iter_method = "DIIS/Jacobi"
398 : ELSE
399 49733 : scf_env%iter_method = "DIIS/Diag."
400 : END IF
401 : ELSE
402 49854 : IF (scf_env%mixing_method == 0) THEN
403 0 : scf_env%iter_method = "NoMix/Diag."
404 49854 : ELSE IF (scf_env%mixing_method == 1) THEN
405 45026 : scf_env%iter_param = scf_env%p_mix_alpha
406 45026 : IF (use_jacobi) THEN
407 0 : scf_env%iter_method = "P_Mix/Jacobi"
408 : ELSE
409 45026 : scf_env%iter_method = "P_Mix/Diag."
410 : END IF
411 4828 : ELSE IF (scf_env%mixing_method > 1) THEN
412 4828 : scf_env%iter_param = scf_env%mixing_store%alpha
413 4828 : IF (use_jacobi) THEN
414 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Jacobi"
415 : ELSE
416 4828 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Diag."
417 : END IF
418 : END IF
419 : END IF
420 :
421 99657 : IF (scf_env%cholesky_method == cholesky_dbcsr .AND. .NOT. use_cdiis) THEN
422 1072 : ortho_dbcsr => scf_env%ortho_dbcsr
423 3198 : DO ispin = 1, nspin
424 2126 : effective_ks_dbcsr => matrix_ks(ispin)%matrix
425 2126 : IF (use_subspace_fock) THEN
426 0 : 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 3198 : ksbuf1=scf_env%buf1_dbcsr, ksbuf2=scf_env%buf2_dbcsr)
432 : END DO
433 :
434 98585 : ELSE IF (scf_env%cholesky_method > cholesky_off) THEN
435 98151 : cdiis_cholesky_method = scf_env%cholesky_method
436 98151 : IF (cdiis_cholesky_method == cholesky_dbcsr) THEN
437 0 : cdiis_cholesky_method = cholesky_inverse
438 : END IF
439 98151 : IF (scf_env%cholesky_method == cholesky_inverse) THEN
440 264 : ortho => scf_env%ortho_m1
441 : ELSE
442 97887 : ortho => scf_env%ortho
443 : END IF
444 :
445 98151 : owns_ortho = .FALSE.
446 98151 : IF (.NOT. ASSOCIATED(ortho)) THEN
447 0 : ALLOCATE (ortho)
448 0 : owns_ortho = .TRUE.
449 : END IF
450 :
451 212506 : DO ispin = 1, nspin
452 212506 : 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 0 : work=scf_env%scf_work2)
457 : ELSE
458 114355 : 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 124 : 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 114231 : use_jacobi=use_jacobi)
477 : END IF
478 : END IF
479 : END DO
480 :
481 98151 : IF (scf_env%cholesky_method /= cholesky_dbcsr) THEN
482 98151 : scf_env%cholesky_method = cdiis_cholesky_method
483 : END IF
484 :
485 98151 : IF (owns_ortho) DEALLOCATE (ortho)
486 : ELSE
487 434 : ortho => scf_env%ortho
488 :
489 434 : owns_ortho = .FALSE.
490 434 : IF (.NOT. ASSOCIATED(ortho)) THEN
491 0 : ALLOCATE (ortho)
492 0 : owns_ortho = .TRUE.
493 : END IF
494 :
495 434 : IF (do_level_shift) THEN
496 172 : DO ispin = 1, nspin
497 : IF (ASSOCIATED(scf_env%scf_work1_red) .AND. ASSOCIATED(scf_env%scf_work2_red) &
498 172 : .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 86 : 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 0 : jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
522 : END IF
523 : END DO
524 : ELSE
525 778 : DO ispin = 1, nspin
526 : IF (ASSOCIATED(scf_env%scf_work1_red) .AND. ASSOCIATED(scf_env%scf_work2_red) &
527 778 : .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 430 : 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 0 : jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
548 : END IF
549 : END DO
550 : END IF
551 :
552 434 : IF (owns_ortho) DEALLOCATE (ortho)
553 : END IF
554 :
555 99657 : 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 98617 : 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 98617 : nspin = SIZE(matrix_ks)
585 :
586 : CALL general_eigenproblem(scf_env, mos, matrix_ks, &
587 98617 : matrix_s, scf_control, scf_section, diis_step)
588 :
589 : total_zeff_corr = 0.0_dp
590 98617 : total_zeff_corr = scf_env%sum_zeff_corr
591 :
592 98617 : IF (ABS(total_zeff_corr) > 0.0_dp) THEN
593 : CALL set_mo_occupation(mo_array=mos, &
594 40 : smear=scf_control%smear, tot_zeff_corr=total_zeff_corr)
595 : ELSE
596 98577 : IF (PRESENT(probe) .EQV. .TRUE.) THEN
597 14 : scf_control%smear%do_smear = .FALSE.
598 : CALL set_mo_occupation(mo_array=mos, &
599 : smear=scf_control%smear, &
600 14 : probe=probe)
601 : ELSE
602 98563 : IF (.NOT. scf_control%gce%do_gce) THEN
603 : CALL set_mo_occupation(mo_array=mos, &
604 98499 : smear=scf_control%smear)
605 : ELSE
606 : CALL set_mo_occupation(mo_array=mos, &
607 : smear=scf_control%smear, &
608 64 : gce=scf_control%gce)
609 : END IF
610 : END IF
611 : END IF
612 :
613 213534 : DO ispin = 1, nspin
614 : CALL calculate_density_matrix(mos(ispin), &
615 213534 : scf_env%p_mix_new(ispin, 1)%matrix)
616 : END DO
617 :
618 98617 : 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 34124 : 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 34124 : 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 34124 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
670 34124 : 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 34124 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
676 34124 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
677 34124 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
678 : TYPE(cp_cfm_type) :: cksmat, cmos, csmat, ctmat, cwork
679 34124 : 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 34124 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
683 : TYPE(cp_fm_type), POINTER :: mo_coeff
684 34124 : 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 34124 : 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 34124 : POINTER :: sab_nl
691 : TYPE(qs_matrix_pools_type), POINTER :: mpools
692 34124 : TYPE(rskp_grid_type) :: s_grid
693 34124 : TYPE(rskp_grid_type), ALLOCATABLE, DIMENSION(:) :: ks_grid
694 : TYPE(section_vals_type), POINTER :: scf_section
695 :
696 34124 : CALL timeset(routineN, handle)
697 :
698 34124 : my_store_ot_matrices = .FALSE.
699 34124 : IF (PRESENT(store_ot_matrices)) my_store_ot_matrices = store_ot_matrices
700 34124 : my_transfer_only = .FALSE.
701 34124 : IF (PRESENT(transfer_only)) my_transfer_only = transfer_only
702 988 : IF (my_transfer_only) my_store_ot_matrices = .TRUE.
703 34124 : my_store_ot_kinetic = PRESENT(matrix_t)
704 34124 : IF (my_store_ot_kinetic) my_store_ot_kinetic = ASSOCIATED(matrix_t)
705 34124 : CPASSERT(.NOT. my_transfer_only .OR. .NOT. update_p)
706 34124 : CPASSERT(.NOT. my_store_ot_kinetic .OR. my_transfer_only)
707 34124 : my_added_mos_auto_grow = .FALSE.
708 34124 : IF (PRESENT(added_mos_auto_grow)) added_mos_auto_grow = .FALSE.
709 :
710 34124 : 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 34124 : cell_to_index=cell_to_index, nkp_grid=nkp_grid, lattice_fft=lattice_fft)
714 34124 : CPASSERT(ASSOCIATED(sab_nl))
715 34124 : kplocal = kp_range(2) - kp_range(1) + 1
716 :
717 : !Whether we use DIIS for k-points
718 34124 : do_diis = .FALSE.
719 34124 : use_adiis = update_p .AND. scf_control%diagonalization%update_method == diag_update_method_adiis
720 : IF (use_adiis) THEN
721 122 : CPASSERT(ASSOCIATED(scf_env%scf_subspace_buffer))
722 : END IF
723 122 : adiis_cdiis = .FALSE.
724 122 : 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 122 : PRESENT(diis_error) .AND. PRESENT(qs_env)) THEN
730 : adiis_cdiis = scf_env%iter_delta < scf_control%eps_diis .OR. &
731 108 : 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 34124 : .NOT. use_real_wfn .AND. PRESENT(diis_error) .AND. PRESENT(qs_env)) .OR. &
736 27112 : adiis_cdiis) do_diis = .TRUE.
737 34124 : IF (my_transfer_only) do_diis = .FALSE.
738 :
739 34124 : adiis_ks => matrix_ks
740 34124 : IF (use_adiis) THEN
741 122 : IF (scf_env%scf_subspace_buffer%use_combined_fock) THEN
742 50 : adiis_ks => scf_env%scf_subspace_buffer%combined_fock
743 : END IF
744 : END IF
745 34124 : effective_ks => matrix_ks
746 34124 : IF (use_adiis .AND. .NOT. adiis_cdiis) THEN
747 50 : effective_ks => adiis_ks
748 : END IF
749 34124 : IF (use_adiis .AND. .NOT. do_diis) THEN
750 50 : scf_env%scf_subspace_buffer%diis_weight = 0.0_dp
751 50 : IF (ASSOCIATED(kpoints%scf_diis_buffer)) THEN
752 50 : 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 34124 : para_env => kpoints%blacs_env_all%para_env
757 34124 : nspin = SIZE(matrix_ks, 1)
758 34124 : omp_threads = 1
759 34124 : !$ 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 34124 : .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 34124 : IF (use_adiis) THEN
766 122 : use_omp_kpoints = use_omp_kpoints .AND. adiis_cdiis
767 : END IF
768 : ! Keep accelerator and distributed eigensolver selections on their existing paths.
769 34124 : IF (diag_type == FM_DIAG_TYPE_CUSOLVER .OR. diag_type == FM_DIAG_TYPE_DLAF) THEN
770 0 : use_omp_kpoints = .FALSE.
771 : END IF
772 34124 : IF (diag_check_requested()) use_omp_kpoints = .FALSE.
773 :
774 34070 : 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 0 : cell_to_index, xkp)
778 0 : CALL timestop(handle)
779 0 : RETURN
780 : END IF
781 :
782 : ! allocate some work matrices
783 34124 : ALLOCATE (rmatrix, cmatrix, tmpmat)
784 : CALL dbcsr_create(rmatrix, template=matrix_ks(1, 1)%matrix, &
785 34124 : matrix_type=dbcsr_type_symmetric)
786 : CALL dbcsr_create(cmatrix, template=matrix_ks(1, 1)%matrix, &
787 34124 : matrix_type=dbcsr_type_antisymmetric)
788 : CALL dbcsr_create(tmpmat, template=matrix_ks(1, 1)%matrix, &
789 34124 : matrix_type=dbcsr_type_no_symmetry)
790 34124 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
791 34124 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
792 :
793 34124 : 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 34124 : nspin = SIZE(matrix_ks, 1)
798 34124 : CALL kpoint_operator_context_create(op_ctx, kpoints, fmwork(1), nspin)
799 :
800 : ! fm pools to be used within a kpoint group
801 34124 : CALL get_kpoint_info(kpoints, mpools=mpools)
802 34124 : CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
803 :
804 34124 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
805 34124 : CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
806 :
807 34124 : IF (use_real_wfn) THEN
808 130 : CALL cp_fm_create(rksmat, matrix_struct)
809 130 : CALL cp_fm_create(rsmat, matrix_struct)
810 130 : IF (my_store_ot_kinetic) CALL cp_fm_create(rtmat, matrix_struct)
811 : END IF
812 34124 : IF (.NOT. use_real_wfn .OR. do_diis) THEN
813 33994 : CALL cp_cfm_create(cksmat, matrix_struct)
814 33994 : CALL cp_cfm_create(csmat, matrix_struct)
815 33994 : CALL cp_cfm_create(cwork, matrix_struct)
816 33994 : IF (my_store_ot_kinetic) CALL cp_cfm_create(ctmat, matrix_struct)
817 33994 : IF (.NOT. use_real_wfn) THEN
818 33994 : kp => kpoints%kp_env(1)%kpoint_env
819 33994 : CALL get_mo_set(kp%mos(1, 1), mo_coeff=mo_coeff)
820 33994 : CALL cp_fm_get_info(mo_coeff, matrix_struct=mo_struct)
821 33994 : CALL cp_cfm_create(cmos, mo_struct)
822 : END IF
823 : END IF
824 :
825 34124 : 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 138402 : ALLOCATE (ks_grid(nspin))
831 34116 : SELECT CASE (lattice_fft)
832 : CASE (lattice_fft_auto)
833 34116 : use_grid_fft = nkp >= 27
834 : CASE (lattice_fft_on)
835 8 : use_grid_fft = .TRUE.
836 : CASE (lattice_fft_off)
837 0 : use_grid_fft = .FALSE.
838 : CASE DEFAULT
839 34124 : CPABORT("Unknown lattice FFT selection mode")
840 : END SELECT
841 34124 : 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 214 : max_storage_bytes=512_int_8*1024_int_8**2/INT(nspin + 1, int_8))
845 : END IF
846 34124 : IF (use_grid_fft) THEN
847 8 : 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 4 : max_storage_bytes=512_int_8*1024_int_8**2/INT(nspin + 1, int_8))
851 8 : use_grid_fft = use_grid_fft .AND. spin_grid_fft
852 : END DO
853 : END IF
854 34124 : IF (.NOT. use_grid_fft) THEN
855 34120 : CALL rskp_transform_grid_release(s_grid)
856 70146 : DO ispin = 1, nspin
857 70146 : CALL rskp_transform_grid_release(ks_grid(ispin))
858 : END DO
859 : END IF
860 :
861 34124 : IF (my_store_ot_matrices) THEN
862 988 : ot_components = MERGE(1, 2, use_real_wfn)
863 2530 : DO ikp = 1, kplocal
864 1542 : kp => kpoints%kp_env(ikp)%kpoint_env
865 1542 : CALL cp_fm_release(kp%ot_hmat)
866 1542 : CALL cp_fm_release(kp%ot_smat)
867 15594 : ALLOCATE (kp%ot_hmat(ot_components, nspin), kp%ot_smat(ot_components))
868 4626 : DO ic = 1, ot_components
869 3084 : CALL cp_fm_create(kp%ot_smat(ic), matrix_struct)
870 7826 : DO ispin = 1, nspin
871 6284 : CALL cp_fm_create(kp%ot_hmat(ic, ispin), matrix_struct)
872 : END DO
873 : END DO
874 2530 : IF (my_store_ot_kinetic) THEN
875 192 : CALL cp_fm_release(kp%ot_tmat)
876 768 : ALLOCATE (kp%ot_tmat(ot_components))
877 576 : DO ic = 1, ot_components
878 576 : 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 34124 : IF (use_real_wfn) THEN
886 130 : ninfo = 4
887 130 : IF (my_store_ot_kinetic) ninfo = 6
888 2444 : ALLOCATE (info(kplocal*nspin*nkp_groups, ninfo))
889 : END IF
890 :
891 : ! Setup and start all the communication
892 34124 : indx = 0
893 97088 : DO ikp = 1, kplocal
894 166084 : DO ispin = 1, nspin
895 131960 : IF (use_real_wfn) THEN
896 312 : DO igroup = 1, nkp_groups
897 : ! number of current kpoint
898 156 : ik = kp_dist(1, igroup) + ikp - 1
899 156 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
900 156 : indx = indx + 1
901 : ! FT of matrices KS and S, then transfer to FM type
902 156 : IF (use_grid_fft) THEN
903 0 : CALL rskp_transform_grid_extract(ks_grid(ispin), ik, rmatrix)
904 : ELSE
905 156 : CALL dbcsr_set(rmatrix, 0.0_dp)
906 : CALL rskp_transform(rmatrix=rmatrix, rsmat=effective_ks, ispin=ispin, &
907 156 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
908 : END IF
909 156 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
910 156 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
911 : ! s matrix is not spin dependent
912 156 : IF (use_grid_fft) THEN
913 0 : CALL rskp_transform_grid_extract(s_grid, ik, rmatrix)
914 : ELSE
915 156 : CALL dbcsr_set(rmatrix, 0.0_dp)
916 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_s, ispin=1, &
917 156 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
918 : END IF
919 156 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
920 156 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(3))
921 : ! transfer to the kpoint group with a new blacs environment
922 156 : IF (my_kpgrp) THEN
923 156 : CALL cp_fm_start_copy_general(fmwork(1), rksmat, para_env, info(indx, 1))
924 156 : CALL cp_fm_start_copy_general(fmwork(3), rsmat, para_env, info(indx, 2))
925 : ELSE
926 0 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
927 0 : CALL cp_fm_start_copy_general(fmwork(3), fmdummy, para_env, info(indx, 2))
928 : END IF
929 312 : IF (my_store_ot_kinetic) THEN
930 0 : CALL dbcsr_set(rmatrix, 0.0_dp)
931 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_t, ispin=1, &
932 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
933 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
934 0 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(3))
935 0 : IF (my_kpgrp) THEN
936 0 : CALL cp_fm_start_copy_general(fmwork(3), rtmat, para_env, info(indx, 3))
937 : ELSE
938 0 : 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 68840 : 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 68840 : use_grid_fft, matrix_row=kpoint_spin_free)
949 68840 : IF (my_store_ot_kinetic) THEN
950 : CALL kpoint_operator_start(op_ctx, ikp, ispin, matrix_t, kpoint_slot_t, &
951 230 : 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 34124 : IF (do_diis) THEN
959 27112 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
960 27112 : cdiis_ncall = kpoints%scf_diis_buffer%ncall
961 27112 : CALL qs_diis_b_info_kp(kpoints%scf_diis_buffer, ib, nb)
962 27112 : indx = 0
963 71078 : DO ikp = 1, kplocal
964 120094 : DO ispin = 1, nspin
965 49016 : IF (use_real_wfn) THEN
966 0 : DO igroup = 1, nkp_groups
967 : ! number of current kpoint
968 0 : ik = kp_dist(1, igroup) + ikp - 1
969 0 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
970 0 : indx = indx + 1
971 0 : IF (my_kpgrp) THEN
972 0 : CALL cp_fm_finish_copy_general(rksmat, info(indx, 1))
973 0 : CALL cp_fm_finish_copy_general(rsmat, info(indx, 2))
974 0 : IF (my_store_ot_kinetic) THEN
975 0 : CALL cp_fm_finish_copy_general(rtmat, info(indx, 3))
976 : END IF
977 0 : CALL cp_cfm_scale_and_add_fm(z_zero, cksmat, z_one, rksmat)
978 0 : CALL cp_cfm_scale_and_add_fm(z_zero, csmat, z_one, rsmat)
979 : END IF
980 : END DO !igroup
981 : ELSE
982 49016 : CALL kpoint_operator_finish(op_ctx, ikp, ispin, kpoint_slot_ks, cksmat)
983 49016 : CALL kpoint_operator_finish(op_ctx, ikp, ispin, kpoint_slot_s, csmat)
984 49016 : IF (my_store_ot_kinetic) THEN
985 0 : CALL kpoint_operator_finish(op_ctx, ikp, ispin, kpoint_slot_t, ctmat)
986 : END IF
987 : END IF
988 :
989 49016 : 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 92982 : ispin, ikp, kplocal, scf_section, real_wfn=use_real_wfn)
992 :
993 : END DO !ispin
994 : END DO !ikp
995 :
996 81336 : 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 54152 : scf_section, kpoints%para_env_inter_kp)
1001 27112 : diis_weight = 1.0_dp
1002 27112 : IF (use_adiis) THEN
1003 72 : diis_weight = scf_env%scf_subspace_buffer%diis_weight
1004 72 : IF (.NOT. ieee_is_finite(diis_error) .OR. diis_error >= scf_control%eps_diis) THEN
1005 8 : cdiis_step = .FALSE.
1006 8 : IF (cdiis_ncall < kpoints%scf_diis_buffer%nbuffer) THEN
1007 : ! The rejected state has not overwritten a full circular buffer.
1008 8 : kpoints%scf_diis_buffer%ncall = cdiis_ncall
1009 : ELSE
1010 0 : CALL qs_diis_b_clear_kp(kpoints%scf_diis_buffer)
1011 0 : diis_weight = 0.0_dp
1012 : END IF
1013 : END IF
1014 72 : IF (cdiis_step) THEN
1015 : ! Use one half-weight handover step; MAX_DIIS only controls the buffer capacity.
1016 56 : 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 72 : scf_env%scf_subspace_buffer%diis_weight = diis_weight
1022 72 : diis_step = cdiis_step .AND. diis_weight >= 1.0_dp
1023 72 : IF (.NOT. diis_step) diis_error = scf_env%scf_subspace_buffer%last_old_fock_weight
1024 : ELSE
1025 27040 : diis_step = cdiis_step
1026 : END IF
1027 :
1028 : !build the ks matrices and diagonalize
1029 98190 : DO ikp = 1, kplocal
1030 120094 : DO ispin = 1, nspin
1031 49016 : IF (use_adiis) THEN
1032 120 : 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 36 : use_grid=.FALSE., matrix_row=ispin)
1037 36 : CALL kpoint_operator_finish(op_ctx, ikp, ispin, kpoint_slot_ks, cwork)
1038 : END IF
1039 : END IF
1040 49016 : kp => kpoints%kp_env(ikp)%kpoint_env
1041 49016 : CALL cp_cfm_to_cfm(kpoints%scf_diis_buffer%smat(ikp), csmat)
1042 :
1043 49016 : CALL cp_cfm_set_all(cksmat, z_zero)
1044 212684 : DO jb = 1, nb
1045 212684 : CALL cp_cfm_scale_and_add(z_one, cksmat, coeffs(jb), kpoints%scf_diis_buffer%param(jb, ispin, ikp))
1046 : END DO
1047 49016 : IF (use_adiis) THEN
1048 120 : IF (diis_weight < 1.0_dp) THEN
1049 : CALL cp_cfm_scale_and_add(CMPLX(diis_weight, KIND=dp), cksmat, &
1050 36 : CMPLX(1.0_dp - diis_weight, KIND=dp), cwork)
1051 : END IF
1052 : END IF
1053 :
1054 92982 : IF (use_real_wfn) THEN
1055 0 : CALL cp_cfm_to_fm(cksmat, rksmat, fmlocal)
1056 0 : CALL cp_cfm_to_fm(csmat, rsmat, fmlocal)
1057 0 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=mo_coeff, eigenvalues=eigenvalues)
1058 0 : IF (scf_env%cholesky_method == cholesky_off) THEN
1059 : CALL cp_fm_geeig_canon(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal, &
1060 0 : scf_control%eps_eigval)
1061 : ELSE
1062 0 : CALL cp_fm_geeig(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal)
1063 : END IF
1064 : ELSE
1065 49016 : CALL get_mo_set(kp%mos(1, ispin), eigenvalues=eigenvalues)
1066 49016 : IF (scf_env%cholesky_method == cholesky_off) THEN
1067 : CALL cp_cfm_geeig_canon(cksmat, csmat, cmos, eigenvalues, cwork, &
1068 50 : scf_control%eps_eigval)
1069 : ELSE
1070 48966 : 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 49016 : 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 7012 : IF (.NOT. update_p .OR. use_adiis) diis_step = .FALSE.
1080 7012 : indx = 0
1081 26010 : DO ikp = 1, kplocal
1082 45990 : DO ispin = 1, nspin
1083 19980 : IF (use_real_wfn) THEN
1084 312 : DO igroup = 1, nkp_groups
1085 : ! number of current kpoint
1086 156 : ik = kp_dist(1, igroup) + ikp - 1
1087 156 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
1088 0 : indx = indx + 1
1089 156 : IF (my_kpgrp) THEN
1090 156 : CALL cp_fm_finish_copy_general(rksmat, info(indx, 1))
1091 156 : CALL cp_fm_finish_copy_general(rsmat, info(indx, 2))
1092 156 : IF (my_store_ot_kinetic) THEN
1093 0 : CALL cp_fm_finish_copy_general(rtmat, info(indx, 3))
1094 : END IF
1095 : END IF
1096 : END DO
1097 : ELSE
1098 19824 : CALL kpoint_operator_finish(op_ctx, ikp, ispin, kpoint_slot_ks, cksmat)
1099 19824 : CALL kpoint_operator_finish(op_ctx, ikp, ispin, kpoint_slot_s, csmat)
1100 19824 : IF (my_store_ot_kinetic) THEN
1101 230 : 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 19980 : kp => kpoints%kp_env(ikp)%kpoint_env
1108 38978 : IF (use_real_wfn) THEN
1109 156 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=mo_coeff, eigenvalues=eigenvalues)
1110 156 : IF (my_store_ot_matrices) THEN
1111 0 : CALL cp_fm_to_fm(rksmat, kp%ot_hmat(1, ispin))
1112 0 : IF (ispin == 1) CALL cp_fm_to_fm(rsmat, kp%ot_smat(1))
1113 0 : IF (ispin == 1 .AND. my_store_ot_kinetic) THEN
1114 0 : CALL cp_fm_to_fm(rtmat, kp%ot_tmat(1))
1115 : END IF
1116 : END IF
1117 156 : IF (.NOT. my_transfer_only) THEN
1118 156 : IF (scf_env%cholesky_method == cholesky_off) THEN
1119 : CALL cp_fm_geeig_canon(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal, &
1120 40 : scf_control%eps_eigval)
1121 : ELSE
1122 116 : CALL cp_fm_geeig(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal)
1123 : END IF
1124 : END IF
1125 : ELSE
1126 19824 : CALL get_mo_set(kp%mos(1, ispin), eigenvalues=eigenvalues)
1127 19824 : IF (my_store_ot_matrices) THEN
1128 1600 : CALL cp_cfm_to_fm(cksmat, kp%ot_hmat(1, ispin), kp%ot_hmat(2, ispin))
1129 1600 : IF (ispin == 1) THEN
1130 1542 : CALL cp_cfm_to_fm(csmat, kp%ot_smat(1), kp%ot_smat(2))
1131 1542 : IF (my_store_ot_kinetic) THEN
1132 192 : CALL cp_cfm_to_fm(ctmat, kp%ot_tmat(1), kp%ot_tmat(2))
1133 : END IF
1134 : END IF
1135 : END IF
1136 19824 : IF (.NOT. my_transfer_only) THEN
1137 18224 : IF (scf_env%cholesky_method == cholesky_off) THEN
1138 : CALL cp_cfm_geeig_canon(cksmat, csmat, cmos, eigenvalues, cwork, &
1139 3914 : scf_control%eps_eigval)
1140 : ELSE
1141 14310 : 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 18224 : 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 34124 : IF (use_real_wfn) THEN
1155 130 : indx = 0
1156 260 : DO ikp = 1, kplocal
1157 416 : DO ispin = 1, nspin
1158 442 : DO igroup = 1, nkp_groups
1159 : ! number of current kpoint
1160 156 : ik = kp_dist(1, igroup) + ikp - 1
1161 156 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
1162 156 : indx = indx + 1
1163 156 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
1164 156 : CALL cp_fm_cleanup_copy_general(info(indx, 2))
1165 312 : IF (my_store_ot_kinetic) THEN
1166 0 : 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 34124 : CALL kpoint_operator_context_release(op_ctx)
1174 :
1175 : ! All done
1176 34748 : IF (ALLOCATED(info)) DEALLOCATE (info)
1177 :
1178 34124 : IF (update_p .AND. .NOT. my_transfer_only) THEN
1179 : ! MO occupations
1180 32994 : IF (PRESENT(probe) .EQV. .TRUE.) THEN
1181 0 : scf_control%smear%do_smear = .FALSE.
1182 : CALL kpoint_set_mo_occupation(kpoints, scf_control%smear, &
1183 0 : probe=probe)
1184 : ELSE
1185 : CALL kpoint_set_mo_occupation(kpoints, scf_control%smear, &
1186 : added_mos_auto=scf_control%added_mos_auto, &
1187 32994 : added_mos_auto_grow=my_added_mos_auto_grow)
1188 : END IF
1189 32994 : IF (my_added_mos_auto_grow) THEN
1190 0 : IF (PRESENT(added_mos_auto_grow)) added_mos_auto_grow = .TRUE.
1191 : ELSE
1192 : ! density matrices
1193 32994 : CALL kpoint_density_matrices(kpoints)
1194 : ! density matrices in real space
1195 : CALL kpoint_density_transform(kpoints, scf_env%p_mix_new, .FALSE., &
1196 32994 : matrix_s(1, 1)%matrix, sab_nl, fmwork, overlap_rs=matrix_s)
1197 : END IF
1198 : END IF
1199 :
1200 34124 : CALL dbcsr_deallocate_matrix(rmatrix)
1201 34124 : CALL dbcsr_deallocate_matrix(cmatrix)
1202 34124 : CALL dbcsr_deallocate_matrix(tmpmat)
1203 34124 : CALL rskp_transform_grid_release(s_grid)
1204 70154 : DO ispin = 1, nspin
1205 70154 : CALL rskp_transform_grid_release(ks_grid(ispin))
1206 : END DO
1207 70154 : DEALLOCATE (ks_grid)
1208 :
1209 34124 : IF (use_real_wfn) THEN
1210 130 : CALL cp_fm_release(rksmat)
1211 130 : CALL cp_fm_release(rsmat)
1212 130 : IF (my_store_ot_kinetic) CALL cp_fm_release(rtmat)
1213 : END IF
1214 34124 : IF (.NOT. use_real_wfn .OR. do_diis) THEN
1215 33994 : CALL cp_cfm_release(cksmat)
1216 33994 : CALL cp_cfm_release(csmat)
1217 33994 : CALL cp_cfm_release(cwork)
1218 33994 : IF (my_store_ot_kinetic) CALL cp_cfm_release(ctmat)
1219 : END IF
1220 34124 : IF (.NOT. use_real_wfn) THEN
1221 33994 : CALL cp_cfm_release(cmos)
1222 : END IF
1223 34124 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
1224 :
1225 34124 : CALL timestop(handle)
1226 :
1227 102372 : 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 0 : 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 0 : 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 0 : DIMENSION(:) :: workspace
1281 0 : TYPE(kp_transform_plan_type) :: transform_plan
1282 : TYPE(kpoint_env_type), POINTER :: kp
1283 : TYPE(section_vals_type), POINTER :: scf_section
1284 :
1285 0 : CALL timeset(routineN, handle)
1286 0 : NULLIFY (matrix_struct, mo_struct, mo_coeff, kp, scf_section)
1287 :
1288 0 : kplocal = SIZE(xkp, 2)
1289 0 : nspin = SIZE(matrix_ks, 1)
1290 : use_adiis = update_p .AND. &
1291 0 : scf_control%diagonalization%update_method == diag_update_method_adiis
1292 : nworkers = 1
1293 0 : !$ nworkers = MIN(omp_get_max_threads(), kplocal)
1294 0 : CALL cp_fm_get_info(scf_env%scf_work1(1), matrix_struct=matrix_struct)
1295 0 : kp => kpoints%kp_env(1)%kpoint_env
1296 0 : CALL get_mo_set(kp%mos(1, 1), mo_coeff=mo_coeff)
1297 0 : CALL cp_fm_get_info(mo_coeff, matrix_struct=mo_struct)
1298 0 : 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 0 : matrix_ks(1, 1)%matrix)
1301 :
1302 0 : IF (do_diis) THEN
1303 0 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
1304 0 : cdiis_ncall = kpoints%scf_diis_buffer%ncall
1305 0 : 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 0 : 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 0 : !$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 0 : 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 0 : scf_section, kpoints%para_env_inter_kp)
1325 0 : diis_weight = 1.0_dp
1326 0 : IF (use_adiis) THEN
1327 0 : diis_weight = scf_env%scf_subspace_buffer%diis_weight
1328 0 : IF (.NOT. ieee_is_finite(diis_error) .OR. diis_error >= scf_control%eps_diis) THEN
1329 0 : cdiis_step = .FALSE.
1330 0 : IF (cdiis_ncall < kpoints%scf_diis_buffer%nbuffer) THEN
1331 0 : kpoints%scf_diis_buffer%ncall = cdiis_ncall
1332 : ELSE
1333 0 : CALL qs_diis_b_clear_kp(kpoints%scf_diis_buffer)
1334 0 : diis_weight = 0.0_dp
1335 : END IF
1336 : END IF
1337 0 : IF (cdiis_step) THEN
1338 0 : 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 0 : diis_weight = 0.0_dp
1342 : END IF
1343 0 : scf_env%scf_subspace_buffer%diis_weight = diis_weight
1344 0 : diis_step = cdiis_step .AND. diis_weight >= 1.0_dp
1345 0 : IF (.NOT. diis_step) diis_error = scf_env%scf_subspace_buffer%last_old_fock_weight
1346 : ELSE
1347 0 : 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 0 : !$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 0 : DEALLOCATE (coeffs)
1362 : ELSE
1363 0 : 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 0 : !$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 0 : IF (update_p) THEN
1377 0 : IF (PRESENT(probe)) THEN
1378 0 : scf_control%smear%do_smear = .FALSE.
1379 0 : CALL kpoint_set_mo_occupation(kpoints, scf_control%smear, probe=probe)
1380 : ELSE
1381 0 : CALL kpoint_set_mo_occupation(kpoints, scf_control%smear)
1382 : END IF
1383 0 : 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 0 : overlap_rs=matrix_s)
1387 : END IF
1388 :
1389 0 : CALL kp_transform_plan_release(transform_plan)
1390 0 : CALL kp_diag_omp_workspaces_release(workspace)
1391 0 : CALL timestop(handle)
1392 :
1393 0 : 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 0 : 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 0 : ALLOCATE (workspace(nworkers))
1412 0 : DO i = 1, nworkers
1413 0 : CALL cp_cfm_create(workspace(i)%cksmat, matrix_struct)
1414 0 : CALL cp_cfm_create(workspace(i)%csmat, matrix_struct)
1415 0 : CALL cp_cfm_create(workspace(i)%csmat_base, matrix_struct)
1416 0 : CALL cp_cfm_create(workspace(i)%cwork, matrix_struct)
1417 0 : CALL cp_cfm_create(workspace(i)%cmos, mo_struct)
1418 : END DO
1419 :
1420 0 : END SUBROUTINE kp_diag_omp_workspaces_create
1421 :
1422 : ! **************************************************************************************************
1423 : !> \brief Release OpenMP k-point workspaces.
1424 : !> \param workspace ...
1425 : ! **************************************************************************************************
1426 0 : 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 0 : DO i = 1, SIZE(workspace)
1434 0 : CALL cp_cfm_release(workspace(i)%cksmat)
1435 0 : CALL cp_cfm_release(workspace(i)%csmat)
1436 0 : CALL cp_cfm_release(workspace(i)%csmat_base)
1437 0 : CALL cp_cfm_release(workspace(i)%cwork)
1438 0 : CALL cp_cfm_release(workspace(i)%cmos)
1439 : END DO
1440 0 : DEALLOCATE (workspace)
1441 :
1442 0 : 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 0 : 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 0 : COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: full
1463 : INTEGER :: col_offset, i, row_offset
1464 : LOGICAL :: found
1465 : REAL(KIND=dp) :: arg
1466 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
1467 :
1468 0 : full => target%local_data
1469 0 : full = z_zero
1470 0 : DO i = 1, plan%nentry
1471 : CALL dbcsr_get_readonly_block_p(source(ispin, plan%image(i))%matrix, &
1472 0 : plan%row(i), plan%col(i), block, found)
1473 0 : IF (.NOT. found) CYCLE
1474 0 : row_offset = plan%row_offset(i)
1475 0 : col_offset = plan%col_offset(i)
1476 0 : arg = DOT_PRODUCT(REAL(plan%cell(:, i), KIND=dp), xkp)
1477 0 : 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 0 : col_offset:col_offset + SIZE(block, 2) - 1) + phase*block
1482 0 : 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 0 : row_offset:row_offset + SIZE(block, 1) - 1) + CONJG(phase)*TRANSPOSE(block)
1487 : END IF
1488 : END DO
1489 :
1490 0 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1520 : TYPE(cp_fm_type), POINTER :: imos, rmos
1521 : TYPE(kpoint_env_type), POINTER :: kp
1522 :
1523 0 : kp => kpoints%kp_env(ikp)%kpoint_env
1524 0 : CALL kp_build_cfm(matrix_s, 1, xkp(1:3, ikp), transform_plan, workspace%csmat_base)
1525 0 : DO ispin = 1, nspin
1526 0 : CALL kp_build_cfm(matrix_ks, ispin, xkp(1:3, ikp), transform_plan, workspace%cksmat)
1527 0 : CALL cp_cfm_to_cfm(workspace%csmat_base, workspace%csmat)
1528 0 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
1529 0 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1530 0 : IF (scf_env%cholesky_method == cholesky_off) THEN
1531 : CALL cp_cfm_geeig_canon_local(workspace%cksmat, workspace%csmat, workspace%cmos, &
1532 0 : eigenvalues, workspace%cwork, scf_control%eps_eigval)
1533 : ELSE
1534 0 : CALL cp_cfm_geeig_local(workspace%cksmat, workspace%csmat, workspace%cmos, eigenvalues)
1535 : END IF
1536 0 : kp%mos(2, ispin)%eigenvalues = eigenvalues
1537 0 : CALL cp_cfm_to_fm(workspace%cmos, rmos, imos)
1538 : END DO
1539 :
1540 0 : 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 0 : 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 0 : kp => kpoints%kp_env(ikp)%kpoint_env
1572 : CALL kp_build_cfm(matrix_s, 1, xkp(1:3, ikp), transform_plan, &
1573 0 : workspace%csmat_base)
1574 0 : DO ispin = 1, nspin
1575 : CALL kp_build_cfm(matrix_ks, ispin, xkp(1:3, ikp), transform_plan, &
1576 0 : workspace%cksmat)
1577 0 : CALL cp_cfm_to_cfm(workspace%csmat_base, workspace%csmat)
1578 0 : CALL kp_diis_calc_err_local(kpoints, kp, workspace, ib, ispin, ikp)
1579 : END DO
1580 :
1581 0 : 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 0 : 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 0 : 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 0 : 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 0 : mo_coeff=rmos, maxocc=maxocc)
1608 0 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1609 0 : CALL cp_fm_get_info(rmos, local_data=cr)
1610 0 : CALL cp_fm_get_info(imos, local_data=ci)
1611 :
1612 0 : c => workspace%cmos%local_data
1613 0 : c(1:nao, 1:nmo) = CMPLX(cr(1:nao, 1:nmo), ci(1:nao, 1:nmo), KIND=dp)
1614 :
1615 0 : h => kpoints%scf_diis_buffer%param(ib, ispin, ikp)%local_data
1616 0 : s => kpoints%scf_diis_buffer%smat(ikp)%local_data
1617 0 : h(:, :) = workspace%cksmat%local_data(:, :)
1618 0 : s(:, :) = workspace%csmat%local_data(:, :)
1619 :
1620 0 : kc => workspace%cksmat%local_data
1621 0 : sc => workspace%csmat%local_data
1622 : CALL zgemm('N', 'N', nao, homo, nao, CMPLX(maxocc, 0.0_dp, KIND=dp), &
1623 0 : 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 0 : s(1, 1), SIZE(s, 1), c(1, 1), SIZE(c, 1), z_zero, sc(1, 1), SIZE(sc, 1))
1626 :
1627 0 : 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 0 : 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 0 : sc(1, 1), SIZE(sc, 1), -z_one, error(1, 1), SIZE(error, 1))
1632 :
1633 0 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1670 : TYPE(cp_fm_type), POINTER :: imos, rmos
1671 : TYPE(kpoint_env_type), POINTER :: kp
1672 :
1673 0 : kp => kpoints%kp_env(ikp)%kpoint_env
1674 0 : DO ispin = 1, nspin
1675 0 : CALL cp_cfm_to_cfm(kpoints%scf_diis_buffer%smat(ikp), workspace%csmat)
1676 0 : CALL cp_cfm_set_all(workspace%cksmat, z_zero)
1677 0 : DO jb = 1, nb
1678 : CALL cp_cfm_scale_and_add(z_one, workspace%cksmat, coeffs(jb), &
1679 0 : kpoints%scf_diis_buffer%param(jb, ispin, ikp))
1680 : END DO
1681 0 : IF (use_adiis .AND. diis_weight < 1.0_dp) THEN
1682 0 : 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 0 : (1.0_dp - diis_weight)*workspace%csmat_base%local_data(:, :)
1686 : END IF
1687 0 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
1688 0 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1689 0 : IF (scf_env%cholesky_method == cholesky_off) THEN
1690 : CALL cp_cfm_geeig_canon_local(workspace%cksmat, workspace%csmat, workspace%cmos, &
1691 0 : eigenvalues, workspace%cwork, scf_control%eps_eigval)
1692 : ELSE
1693 0 : CALL cp_cfm_geeig_local(workspace%cksmat, workspace%csmat, workspace%cmos, eigenvalues)
1694 : END IF
1695 0 : kp%mos(2, ispin)%eigenvalues = eigenvalues
1696 0 : CALL cp_cfm_to_fm(workspace%cmos, rmos, imos)
1697 : END DO
1698 :
1699 0 : 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 22 : 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 22 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1726 22 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1727 : LOGICAL :: my_kpgrp, use_real_wfn
1728 22 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1729 22 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1730 22 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
1731 : TYPE(cp_cfm_type) :: cksmat, cmos, csmat, cwork
1732 22 : 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 22 : 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 22 : POINTER :: sab_nl
1742 : TYPE(qs_matrix_pools_type), POINTER :: mpools
1743 :
1744 22 : CALL timeset(routineN, handle)
1745 :
1746 22 : 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 22 : cell_to_index=cell_to_index)
1750 22 : CPASSERT(ASSOCIATED(sab_nl))
1751 22 : kplocal = kp_range(2) - kp_range(1) + 1
1752 :
1753 : ! use as template
1754 22 : tempmat => matrix_ks(1, 1)%matrix
1755 :
1756 : ! allocate some work matrices
1757 22 : ALLOCATE (rmatrix, cmatrix, tmpmat)
1758 22 : CALL dbcsr_create(rmatrix, template=tempmat, matrix_type=dbcsr_type_symmetric)
1759 22 : CALL dbcsr_create(cmatrix, template=tempmat, matrix_type=dbcsr_type_antisymmetric)
1760 22 : CALL dbcsr_create(tmpmat, template=tempmat, matrix_type=dbcsr_type_no_symmetry)
1761 22 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
1762 22 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
1763 :
1764 : ! fm pools to be used within a kpoint group
1765 22 : CALL get_kpoint_info(kpoints, mpools=mpools)
1766 22 : CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
1767 :
1768 22 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
1769 22 : CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
1770 :
1771 22 : IF (use_real_wfn) THEN
1772 0 : CALL cp_fm_create(rksmat, matrix_struct)
1773 0 : CALL cp_fm_create(rsmat, matrix_struct)
1774 : ELSE
1775 22 : CALL cp_cfm_create(cksmat, matrix_struct)
1776 22 : CALL cp_cfm_create(csmat, matrix_struct)
1777 22 : CALL cp_cfm_create(cwork, matrix_struct)
1778 22 : kp => kpoints%kp_env(1)%kpoint_env
1779 22 : CALL get_mo_set(kp%mos(1, 1), mo_coeff=mo_coeff)
1780 22 : CALL cp_fm_get_info(mo_coeff, matrix_struct=mo_struct)
1781 22 : CALL cp_cfm_create(cmos, mo_struct)
1782 : END IF
1783 :
1784 22 : para_env => kpoints%blacs_env_all%para_env
1785 22 : nspin = SIZE(matrix_ks, 1)
1786 22 : CALL kpoint_operator_context_create(op_ctx, kpoints, fmwork(1), nspin)
1787 : ! The real-wavefunction path keeps its own transfer bookkeeping.
1788 22 : IF (use_real_wfn) ALLOCATE (info(kplocal*nspin*nkp_groups, 2))
1789 :
1790 : ! Setup and start all the communication
1791 22 : indx = 0
1792 94 : DO ikp = 1, kplocal
1793 166 : DO ispin = 1, nspin
1794 144 : IF (use_real_wfn) THEN
1795 0 : DO igroup = 1, nkp_groups
1796 : ! number of current kpoint
1797 0 : ik = kp_dist(1, igroup) + ikp - 1
1798 0 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
1799 0 : indx = indx + 1
1800 : ! FT of matrices KS and S, then transfer to FM type
1801 0 : CALL dbcsr_set(rmatrix, 0.0_dp)
1802 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_ks, ispin=ispin, &
1803 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1804 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1805 0 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
1806 : ! s matrix is not spin dependent
1807 0 : CALL dbcsr_set(rmatrix, 0.0_dp)
1808 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_s, ispin=1, &
1809 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1810 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1811 0 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(3))
1812 : ! transfer to the kpoint group with a new blacs environment
1813 0 : IF (my_kpgrp) THEN
1814 0 : CALL cp_fm_start_copy_general(fmwork(1), rksmat, para_env, info(indx, 1))
1815 0 : CALL cp_fm_start_copy_general(fmwork(3), rsmat, para_env, info(indx, 2))
1816 : ELSE
1817 0 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
1818 0 : 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 72 : use_grid=.FALSE., matrix_row=ispin)
1826 : CALL kpoint_operator_start(op_ctx, ikp, ispin, matrix_s, kpoint_slot_s, &
1827 72 : 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 22 : indx = 0
1834 94 : DO ikp = 1, kplocal
1835 166 : DO ispin = 1, nspin
1836 72 : IF (use_real_wfn) THEN
1837 0 : DO igroup = 1, nkp_groups
1838 : ! number of current kpoint
1839 0 : ik = kp_dist(1, igroup) + ikp - 1
1840 0 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
1841 0 : indx = indx + 1
1842 0 : IF (my_kpgrp) THEN
1843 0 : CALL cp_fm_finish_copy_general(rksmat, info(indx, 1))
1844 0 : CALL cp_fm_finish_copy_general(rsmat, info(indx, 2))
1845 : END IF
1846 : END DO
1847 : ELSE
1848 72 : CALL kpoint_operator_finish(op_ctx, ikp, ispin, kpoint_slot_ks, cksmat)
1849 72 : 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 72 : kp => kpoints%kp_env(ikp)%kpoint_env
1855 144 : IF (use_real_wfn) THEN
1856 0 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=mo_coeff, eigenvalues=eigenvalues)
1857 0 : CALL cp_fm_geeig(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal)
1858 : ELSE
1859 72 : CALL get_mo_set(kp%mos(1, ispin), eigenvalues=eigenvalues)
1860 72 : 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 72 : 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 22 : IF (use_real_wfn) THEN
1871 0 : indx = 0
1872 0 : DO ikp = 1, kplocal
1873 0 : DO ispin = 1, nspin
1874 0 : DO igroup = 1, nkp_groups
1875 : ! number of current kpoint
1876 0 : ik = kp_dist(1, igroup) + ikp - 1
1877 0 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
1878 0 : indx = indx + 1
1879 0 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
1880 0 : CALL cp_fm_cleanup_copy_general(info(indx, 2))
1881 : END DO
1882 : END DO
1883 : END DO
1884 : END IF
1885 :
1886 22 : CALL kpoint_operator_context_release(op_ctx)
1887 :
1888 : ! All done
1889 22 : IF (ALLOCATED(info)) DEALLOCATE (info)
1890 :
1891 22 : CALL dbcsr_deallocate_matrix(rmatrix)
1892 22 : CALL dbcsr_deallocate_matrix(cmatrix)
1893 22 : CALL dbcsr_deallocate_matrix(tmpmat)
1894 :
1895 22 : IF (use_real_wfn) THEN
1896 0 : CALL cp_fm_release(rksmat)
1897 0 : CALL cp_fm_release(rsmat)
1898 : ELSE
1899 22 : CALL cp_cfm_release(cksmat)
1900 22 : CALL cp_cfm_release(csmat)
1901 22 : CALL cp_cfm_release(cwork)
1902 22 : CALL cp_cfm_release(cmos)
1903 : END IF
1904 22 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
1905 :
1906 22 : CALL timestop(handle)
1907 :
1908 44 : 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 10 : 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 10 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues, mo_occupations
1946 10 : 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 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
1951 10 : 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 10 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
1956 :
1957 10 : CALL timeset(routineN, handle)
1958 10 : NULLIFY (c0, chc, energy, evec, matrix_ks, mo_coeff, mo_eigenvalues, &
1959 10 : mo_occupations, dft_control, rho_ao, rho_ao_kp)
1960 :
1961 10 : logger => cp_get_default_logger()
1962 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIAG_SUB_SCF", &
1963 10 : extension=".scfLog")
1964 :
1965 : !Extra loop keeping mos unchanged and refining the subspace occupation
1966 10 : nspin = SIZE(mos)
1967 10 : CALL qs_rho_get(rho, rho_ao=rho_ao, rho_ao_kp=rho_ao_kp)
1968 :
1969 40 : ALLOCATE (eval_first(nspin))
1970 40 : ALLOCATE (occ_first(nspin))
1971 20 : DO ispin = 1, nspin
1972 : CALL get_mo_set(mo_set=mos(ispin), &
1973 : nmo=nmo, &
1974 : eigenvalues=mo_eigenvalues, &
1975 10 : occupation_numbers=mo_occupations)
1976 30 : ALLOCATE (eval_first(ispin)%array(nmo))
1977 20 : ALLOCATE (occ_first(ispin)%array(nmo))
1978 50 : eval_first(ispin)%array(1:nmo) = mo_eigenvalues(1:nmo)
1979 70 : occ_first(ispin)%array(1:nmo) = mo_occupations(1:nmo)
1980 : END DO
1981 :
1982 20 : DO ispin = 1, nspin
1983 : ! does not yet handle k-points
1984 10 : CALL dbcsr_copy(subspace_env%p_matrix_store(ispin)%matrix, rho_ao(ispin)%matrix)
1985 20 : CALL dbcsr_copy(rho_ao(ispin)%matrix, scf_env%p_mix_new(ispin, 1)%matrix)
1986 : END DO
1987 :
1988 10 : subspace_env%p_matrix_mix => scf_env%p_mix_new
1989 :
1990 10 : 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 10 : dft_control=dft_control)
1997 :
1998 : ! mixing storage allocation
1999 10 : 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 0 : scf_env%p_delta, nspin, subspace_env%mixing_store)
2002 0 : IF (dft_control%qs_control%gapw) THEN
2003 0 : 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 0 : para_env, rho_atom=rho_atom)
2006 0 : ELSE IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
2007 0 : CALL charge_mixing_init(subspace_env%mixing_store)
2008 0 : ELSE IF (dft_control%qs_control%semi_empirical) THEN
2009 0 : CPABORT('SE Code not possible')
2010 : ELSE
2011 0 : CALL mixing_init(subspace_env%mixing_method, rho, subspace_env%mixing_store, para_env)
2012 : END IF
2013 : END IF
2014 :
2015 10 : ene_old = 0.0_dp
2016 : ene_diff = 0.0_dp
2017 10 : IF (output_unit > 0) THEN
2018 0 : 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 0 : "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 10 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
2027 :
2028 22 : DO iloop = 1, subspace_env%max_iter
2029 20 : t1 = m_walltime()
2030 20 : converged = .FALSE.
2031 20 : ene_old = energy%total
2032 :
2033 20 : CALL qs_ks_did_change(ks_env, rho_changed=.TRUE.)
2034 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., &
2035 20 : just_energy=.FALSE., print_active=.FALSE.)
2036 :
2037 20 : max_val = 0.0_dp
2038 20 : sum_val = 0.0_dp
2039 20 : sum_band = 0.0_dp
2040 40 : 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 20 : mo_coeff=mo_coeff)
2047 :
2048 : !compute C'HC
2049 20 : chc => subspace_env%chc_mat(ispin)
2050 20 : evec => subspace_env%c_vec(ispin)
2051 20 : c0 => subspace_env%c0(ispin)
2052 20 : CALL cp_fm_to_fm(mo_coeff, c0)
2053 20 : CALL cp_fm_create(work, c0%matrix_struct)
2054 20 : CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, c0, work, nmo)
2055 20 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, rone, c0, work, rzero, chc)
2056 20 : CALL cp_fm_release(work)
2057 : !diagonalize C'HC
2058 20 : CALL choose_eigv_solver(chc, evec, mo_eigenvalues)
2059 :
2060 : !rotate the mos by the eigenvectors of C'HC
2061 20 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, rone, c0, evec, rzero, mo_coeff)
2062 :
2063 20 : IF (.NOT. scf_control%gce%do_gce) THEN
2064 : CALL set_mo_occupation(mo_set=mos(ispin), &
2065 20 : smear=scf_control%smear)
2066 : ELSE
2067 : CALL set_mo_occupation(mo_set=mos(ispin), &
2068 : smear=scf_control%smear, &
2069 0 : gce=scf_control%gce)
2070 : END IF
2071 :
2072 : ! does not yet handle k-points
2073 : CALL calculate_density_matrix(mos(ispin), &
2074 20 : subspace_env%p_matrix_mix(ispin, 1)%matrix)
2075 :
2076 160 : DO i = 1, nmo
2077 100 : sum_band = sum_band + mo_eigenvalues(i)*mo_occupations(i)
2078 : END DO
2079 :
2080 : !check for self consistency
2081 : END DO
2082 :
2083 20 : IF (subspace_env%mixing_method == direct_mixing_nr) THEN
2084 : CALL scf_env_density_mixing(subspace_env%p_matrix_mix, &
2085 20 : 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 0 : subspace_env%p_matrix_mix, delta=iter_delta)
2089 : END IF
2090 :
2091 40 : DO ispin = 1, nspin
2092 : ! does not yet handle k-points
2093 40 : CALL dbcsr_copy(rho_ao(ispin)%matrix, subspace_env%p_matrix_mix(ispin, 1)%matrix)
2094 : END DO
2095 : ! update of density
2096 20 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
2097 : ! Mixing in reciprocal space
2098 20 : IF (subspace_env%mixing_method >= gspace_mixing_nr) THEN
2099 : CALL gspace_mixing(qs_env, scf_env%mixing_method, subspace_env%mixing_store, &
2100 0 : rho, para_env, scf_env%iter_count)
2101 : END IF
2102 :
2103 20 : ene_diff = energy%total - ene_old
2104 : converged = (ABS(ene_diff) < subspace_env%eps_ene .AND. &
2105 20 : iter_delta < subspace_env%eps_adapt*scf_env%iter_delta)
2106 20 : t2 = m_walltime()
2107 20 : 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 0 : iloop, t2 - t1, iter_delta, sum_band, energy%total, ene_diff
2110 0 : CALL m_flush(output_unit)
2111 : END IF
2112 22 : IF (converged) THEN
2113 8 : IF (output_unit > 0) WRITE (output_unit, "(T10,A,I6,A,/)") &
2114 0 : " Reached convergence in ", iloop, " iterations "
2115 : EXIT
2116 : END IF
2117 :
2118 : END DO ! iloop
2119 :
2120 10 : NULLIFY (subspace_env%p_matrix_mix)
2121 20 : DO ispin = 1, nspin
2122 : ! does not yet handle k-points
2123 10 : CALL dbcsr_copy(scf_env%p_mix_new(ispin, 1)%matrix, rho_ao(ispin)%matrix)
2124 10 : CALL dbcsr_copy(rho_ao(ispin)%matrix, subspace_env%p_matrix_store(ispin)%matrix)
2125 :
2126 20 : DEALLOCATE (eval_first(ispin)%array, occ_first(ispin)%array)
2127 : END DO
2128 10 : DEALLOCATE (eval_first, occ_first)
2129 :
2130 10 : CALL timestop(handle)
2131 :
2132 10 : END SUBROUTINE do_scf_diag_subspace
2133 :
2134 : ! **************************************************************************************************
2135 : !> \brief ...
2136 : !> \param subspace_env ...
2137 : !> \param qs_env ...
2138 : !> \param mos ...
2139 : ! **************************************************************************************************
2140 2 : 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 2 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
2152 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2153 2 : POINTER :: sab_orb
2154 :
2155 2 : CALL timeset(routineN, handle)
2156 :
2157 2 : NULLIFY (sab_orb, matrix_s)
2158 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb, &
2159 2 : matrix_s=matrix_s)
2160 :
2161 2 : nspin = SIZE(mos)
2162 : ! *** allocate p_atrix_store ***
2163 2 : IF (.NOT. ASSOCIATED(subspace_env%p_matrix_store)) THEN
2164 2 : CALL dbcsr_allocate_matrix_set(subspace_env%p_matrix_store, nspin)
2165 :
2166 4 : DO i = 1, nspin
2167 2 : 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 2 : name="DENSITY_STORE", matrix_type=dbcsr_type_symmetric)
2170 : CALL cp_dbcsr_alloc_block_from_nbl(subspace_env%p_matrix_store(i)%matrix, &
2171 2 : sab_orb)
2172 4 : CALL dbcsr_set(subspace_env%p_matrix_store(i)%matrix, 0.0_dp)
2173 : END DO
2174 :
2175 : END IF
2176 :
2177 8 : ALLOCATE (subspace_env%chc_mat(nspin))
2178 6 : ALLOCATE (subspace_env%c_vec(nspin))
2179 6 : ALLOCATE (subspace_env%c0(nspin))
2180 :
2181 4 : DO ispin = 1, nspin
2182 2 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
2183 2 : CALL cp_fm_create(subspace_env%c0(ispin), mo_coeff%matrix_struct)
2184 2 : 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 2 : context=mo_coeff%matrix_struct%context)
2188 2 : CALL cp_fm_create(subspace_env%chc_mat(ispin), fm_struct_tmp, "chc")
2189 2 : CALL cp_fm_create(subspace_env%c_vec(ispin), fm_struct_tmp, "vec")
2190 6 : CALL cp_fm_struct_release(fm_struct_tmp)
2191 : END DO
2192 :
2193 2 : CALL timestop(handle)
2194 :
2195 2 : 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 17876 : 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 17876 : nspin = SIZE(matrix_ks)
2224 :
2225 36518 : DO ispin = 1, nspin
2226 36518 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, scf_env%scf_work1(ispin))
2227 : END DO
2228 17876 : 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 15282 : scf_section=scf_section)
2233 : ELSE
2234 2594 : diis_step = .FALSE.
2235 : END IF
2236 :
2237 17876 : IF ((scf_env%iter_count > 1) .AND. (scf_env%iter_delta < scf_control%diagonalization%eps_jacobi)) THEN
2238 18 : use_jacobi = .TRUE.
2239 : ELSE
2240 17858 : use_jacobi = .FALSE.
2241 : END IF
2242 :
2243 : do_level_shift = ((scf_control%level_shift /= 0.0_dp) .AND. &
2244 17876 : ((scf_control%density_guess == core_guess) .OR. (scf_env%iter_count > 1)))
2245 17876 : IF (diis_step) THEN
2246 11870 : scf_env%iter_param = diis_error
2247 11870 : IF (use_jacobi) THEN
2248 18 : scf_env%iter_method = "DIIS/Jacobi"
2249 : ELSE
2250 11852 : scf_env%iter_method = "DIIS/Diag."
2251 : END IF
2252 : ELSE
2253 6006 : IF (scf_env%mixing_method == 1) THEN
2254 6006 : scf_env%iter_param = scf_env%p_mix_alpha
2255 6006 : IF (use_jacobi) THEN
2256 0 : scf_env%iter_method = "P_Mix/Jacobi"
2257 : ELSE
2258 6006 : scf_env%iter_method = "P_Mix/Diag."
2259 : END IF
2260 0 : ELSE IF (scf_env%mixing_method > 1) THEN
2261 0 : scf_env%iter_param = scf_env%mixing_store%alpha
2262 0 : IF (use_jacobi) THEN
2263 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Jacobi"
2264 : ELSE
2265 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Diag."
2266 : END IF
2267 : END IF
2268 : END IF
2269 17876 : scf_env%iter_delta = 0.0_dp
2270 :
2271 36518 : 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 36518 : jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
2279 : END DO
2280 :
2281 17876 : IF (.NOT. scf_control%gce%do_gce) THEN
2282 : CALL set_mo_occupation(mo_array=mos, &
2283 17876 : smear=scf_control%smear)
2284 : ELSE
2285 : CALL set_mo_occupation(mo_array=mos, &
2286 : smear=scf_control%smear, &
2287 0 : gce=scf_control%gce)
2288 : END IF
2289 :
2290 36518 : DO ispin = 1, nspin
2291 : ! does not yet handle k-points
2292 : CALL calculate_density_matrix(mos(ispin), &
2293 36518 : scf_env%p_mix_new(ispin, 1)%matrix)
2294 : END DO
2295 :
2296 17876 : 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 202 : 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 202 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
2324 : TYPE(cp_fm_type), POINTER :: mo_coeff
2325 :
2326 202 : NULLIFY (eigenvalues)
2327 :
2328 202 : nspin = SIZE(matrix_ks)
2329 :
2330 448 : DO ispin = 1, nspin
2331 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, &
2332 448 : scf_env%scf_work1(ispin))
2333 : END DO
2334 :
2335 202 : 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 180 : scf_section=scf_section)
2341 : ELSE
2342 22 : diis_step = .FALSE.
2343 : END IF
2344 :
2345 202 : eps_iter = scf_control%diagonalization%eps_iter
2346 202 : IF (diis_step) THEN
2347 130 : scf_env%iter_param = diis_error
2348 130 : scf_env%iter_method = "DIIS/OTdiag"
2349 274 : DO ispin = 1, nspin
2350 : CALL copy_fm_to_dbcsr(scf_env%scf_work1(ispin), &
2351 274 : matrix_ks(ispin)%matrix, keep_sparsity=.TRUE.)
2352 : END DO
2353 130 : eps_iter = MAX(eps_iter, scf_control%diagonalization%eps_adapt*diis_error)
2354 : ELSE
2355 72 : IF (scf_env%mixing_method == 1) THEN
2356 72 : scf_env%iter_param = scf_env%p_mix_alpha
2357 72 : scf_env%iter_method = "P_Mix/OTdiag."
2358 0 : ELSE IF (scf_env%mixing_method > 1) THEN
2359 0 : scf_env%iter_param = scf_env%mixing_store%alpha
2360 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/OTdiag."
2361 : END IF
2362 : END IF
2363 :
2364 202 : scf_env%iter_delta = 0.0_dp
2365 :
2366 448 : DO ispin = 1, nspin
2367 : CALL get_mo_set(mos(ispin), &
2368 : mo_coeff=mo_coeff, &
2369 : eigenvalues=eigenvalues, &
2370 : nmo=nmo, &
2371 246 : 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 246 : ot_settings=scf_control%diagonalization%ot_settings)
2381 : CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin)%matrix, &
2382 : evals_arg=eigenvalues, &
2383 246 : do_rotation=.TRUE.)
2384 : CALL copy_fm_to_dbcsr(mos(ispin)%mo_coeff, &
2385 694 : mos(ispin)%mo_coeff_b)
2386 : !fm->dbcsr
2387 : END DO
2388 :
2389 : CALL set_mo_occupation(mo_array=mos, &
2390 202 : smear=scf_control%smear)
2391 :
2392 448 : DO ispin = 1, nspin
2393 : ! does not yet handle k-points
2394 : CALL calculate_density_matrix(mos(ispin), &
2395 448 : scf_env%p_mix_new(ispin, 1)%matrix)
2396 : END DO
2397 :
2398 202 : 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 58 : 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 58 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
2424 : TYPE(cp_fm_type) :: coeff_im_work
2425 58 : 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 58 : POINTER :: sab_nl
2432 :
2433 58 : NULLIFY (coeff_im, coeff_re, eigenvalues, fmwork, kp, matrix_h_im, matrix_h_re, &
2434 58 : matrix_s_im, matrix_s_re, matrix_t_im, matrix_t_re, sab_nl)
2435 58 : added_mos_auto_grow = .FALSE.
2436 58 : diis_step = .FALSE.
2437 58 : eps_iter = scf_control%diagonalization%eps_iter
2438 58 : nspin = SIZE(matrix_ks, 1)
2439 : store_kinetic = scf_control%diagonalization%ot_settings%preconditioner_type == &
2440 58 : ot_precond_full_kinetic
2441 58 : IF (store_kinetic) THEN
2442 0 : CPASSERT(ASSOCIATED(matrix_t))
2443 : END IF
2444 58 : CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn, sab_nl=sab_nl)
2445 58 : 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 58 : IF (scf_env%iter_count <= 1) THEN
2450 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, &
2451 6 : .FALSE., diis_step)
2452 : END IF
2453 58 : 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 0 : 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 58 : .FALSE., diis_step, store_ot_matrices=.TRUE., transfer_only=.TRUE.)
2460 : END IF
2461 :
2462 58 : scf_env%iter_delta = 0.0_dp
2463 136 : DO ikp = 1, SIZE(kpoints%kp_env)
2464 78 : kp => kpoints%kp_env(ikp)%kpoint_env
2465 78 : CPASSERT(ASSOCIATED(kp%ot_hmat) .AND. ASSOCIATED(kp%ot_smat))
2466 78 : 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 78 : 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 0 : matrix_s_re, matrix_s_im)
2474 0 : CALL dbcsr_set(matrix_s_im, 0.0_dp)
2475 : END IF
2476 78 : IF (store_kinetic) THEN
2477 0 : CPASSERT(ASSOCIATED(kp%ot_tmat))
2478 0 : 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 0 : 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 0 : matrix_t_re, matrix_t_im)
2486 0 : CALL dbcsr_set(matrix_t_im, 0.0_dp)
2487 : END IF
2488 : END IF
2489 156 : DO ispin = 1, nspin
2490 78 : ks_spin = MIN(ispin, SIZE(kp%ot_hmat, 2))
2491 78 : 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 78 : 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 0 : matrix_h_re, matrix_h_im)
2499 0 : CALL dbcsr_set(matrix_h_im, 0.0_dp)
2500 : END IF
2501 78 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=coeff_re, eigenvalues=eigenvalues)
2502 78 : IF (use_real_wfn) THEN
2503 : CALL cp_fm_create(coeff_im_work, coeff_re%matrix_struct, &
2504 0 : name='real k-point OT imaginary work')
2505 0 : 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 0 : mo_set=kp%mos(1, ispin))
2515 0 : CALL cp_fm_release(coeff_im_work)
2516 : ELSE
2517 78 : 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 78 : mo_set=kp%mos(1, ispin))
2527 1254 : kp%mos(2, ispin)%eigenvalues = eigenvalues
2528 : END IF
2529 78 : CALL dbcsr_release_p(matrix_h_re)
2530 156 : CALL dbcsr_release_p(matrix_h_im)
2531 : END DO
2532 78 : CALL dbcsr_release_p(matrix_s_re)
2533 78 : CALL dbcsr_release_p(matrix_s_im)
2534 78 : IF (ASSOCIATED(matrix_t_re)) CALL dbcsr_release_p(matrix_t_re)
2535 136 : 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 58 : added_mos_auto_grow=added_mos_auto_grow)
2541 58 : IF (added_mos_auto_grow) RETURN
2542 :
2543 58 : CALL kpoint_density_matrices(kpoints)
2544 58 : fmwork => scf_env%scf_work1
2545 : CALL kpoint_density_transform(kpoints, scf_env%p_mix_new, .FALSE., &
2546 58 : matrix_s(1, 1)%matrix, sab_nl, fmwork, overlap_rs=matrix_s)
2547 58 : IF (scf_env%mixing_method == 0) THEN
2548 0 : scf_env%iter_method = 'NoMix/OTdiag.'
2549 58 : ELSE IF (scf_env%mixing_method == 1) THEN
2550 20 : scf_env%iter_param = scf_env%p_mix_alpha
2551 20 : scf_env%iter_method = 'P_Mix/OTdiag.'
2552 : ELSE
2553 38 : scf_env%iter_param = scf_env%mixing_store%alpha
2554 38 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//'/OTdiag.'
2555 : END IF
2556 58 : 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 1132 : 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 1132 : 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 1132 : CALL timeset(routineN, handle)
2602 :
2603 1132 : IF (scf_env%cholesky_method == cholesky_inverse) THEN
2604 0 : ortho => scf_env%ortho_m1
2605 : ELSE
2606 1132 : ortho => scf_env%ortho
2607 : END IF
2608 1132 : work => scf_env%scf_work2
2609 :
2610 1132 : ksa => scf_env%scf_work1(1)
2611 1132 : ksb => scf_env%scf_work1(2)
2612 :
2613 1132 : CALL copy_dbcsr_to_fm(matrix_ks(1)%matrix, ksa)
2614 1132 : 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 1132 : 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 1132 : mo_coeff=mob)
2633 :
2634 : ! Define the amount of level-shifting
2635 :
2636 1132 : 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 20 : level_shift_loc = scf_control%level_shift
2641 : ELSE
2642 1112 : level_shift_loc = 0.0_dp
2643 : END IF
2644 :
2645 : IF ((scf_env%iter_count > 1) .OR. &
2646 1132 : (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 1024 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, moa, 0.0_dp, work)
2653 1024 : CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, moa, work, 0.0_dp, ksa)
2654 :
2655 1024 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksb, moa, 0.0_dp, work)
2656 1024 : 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 1024 : IF (scf_control%roks_scheme == general_roks) THEN
2662 : CALL combine_ks_matrices(ksa, ksb, occa, occb, scf_control%roks_f, &
2663 0 : nalpha, nbeta)
2664 1024 : ELSE IF (scf_control%roks_scheme == high_spin_roks) THEN
2665 1024 : CALL combine_ks_matrices(ksa, ksb, occa, occb, scf_control%roks_parameter)
2666 : ELSE
2667 0 : 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 1024 : IF (orthogonal_basis) THEN
2674 : ! Q = C
2675 464 : mo2ao => moa
2676 : ELSE
2677 : ! Q = S*C
2678 560 : 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 560 : 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 1024 : CALL parallel_gemm("N", "T", nao, nao, nao, 1.0_dp, ksa, mo2ao, 0.0_dp, work)
2687 1024 : 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 108 : CALL cp_fm_to_fm(ksb, ksa)
2696 :
2697 : END IF
2698 :
2699 : ! Update DIIS buffer and possibly perform DIIS extrapolation step
2700 :
2701 1132 : IF (scf_env%iter_count > 1) THEN
2702 1018 : 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 460 : roks=.TRUE.)
2713 460 : 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 558 : roks=.TRUE.)
2726 : END IF
2727 : END IF
2728 :
2729 1132 : IF (diis_step) THEN
2730 714 : scf_env%iter_param = diis_error
2731 714 : scf_env%iter_method = "DIIS/Diag."
2732 : ELSE
2733 418 : IF (scf_env%mixing_method == 1) THEN
2734 418 : scf_env%iter_param = scf_env%p_mix_alpha
2735 418 : scf_env%iter_method = "P_Mix/Diag."
2736 0 : ELSE IF (scf_env%mixing_method > 1) THEN
2737 0 : scf_env%iter_param = scf_env%mixing_store%alpha
2738 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Diag."
2739 : END IF
2740 : END IF
2741 :
2742 1132 : scf_env%iter_delta = 0.0_dp
2743 :
2744 1132 : 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 20 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, moa, 0.0_dp, work)
2750 20 : 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 60 : DO imo = homob + 1, homoa
2755 60 : CALL cp_fm_add_to_element(ksa, imo, imo, 0.5_dp*level_shift_loc)
2756 : END DO
2757 220 : DO imo = homoa + 1, nmo
2758 220 : CALL cp_fm_add_to_element(ksa, imo, imo, level_shift_loc)
2759 : END DO
2760 :
2761 1112 : ELSE IF (.NOT. orthogonal_basis) THEN
2762 :
2763 : ! Transform the current Kohn-Sham matrix to an orthogonal basis
2764 590 : SELECT CASE (scf_env%cholesky_method)
2765 : CASE (cholesky_reduce)
2766 0 : CALL cp_fm_cholesky_reduce(ksa, ortho)
2767 : CASE (cholesky_restore)
2768 526 : CALL cp_fm_uplo_to_full(ksa, work)
2769 : CALL cp_fm_cholesky_restore(ksa, nao, ortho, work, &
2770 526 : "SOLVE", pos="RIGHT")
2771 : CALL cp_fm_cholesky_restore(work, nao, ortho, ksa, &
2772 526 : "SOLVE", pos="LEFT", transa="T")
2773 : CASE (cholesky_inverse)
2774 0 : CALL cp_fm_uplo_to_full(ksa, work)
2775 : CALL cp_fm_cholesky_restore(ksa, nao, ortho, work, &
2776 0 : "MULTIPLY", pos="RIGHT")
2777 : CALL cp_fm_cholesky_restore(work, nao, ortho, ksa, &
2778 0 : "MULTIPLY", pos="LEFT", transa="T")
2779 : CASE (cholesky_off)
2780 64 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, ortho, 0.0_dp, work)
2781 654 : 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 1132 : CALL choose_eigv_solver(ksa, work, eiga)
2789 :
2790 : ! Back-transformation of the orthonormal eigenvectors if needed
2791 :
2792 1132 : IF (level_shift_loc /= 0.0_dp) THEN
2793 : ! Use old MO set for back-transformation if level-shifting was applied
2794 20 : CALL cp_fm_to_fm(moa, ortho)
2795 20 : CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, ortho, work, 0.0_dp, moa)
2796 : ELSE
2797 1112 : IF (orthogonal_basis) THEN
2798 522 : CALL cp_fm_to_fm(work, moa)
2799 : ELSE
2800 1116 : SELECT CASE (scf_env%cholesky_method)
2801 : CASE (cholesky_reduce, cholesky_restore)
2802 526 : CALL cp_fm_cholesky_restore(work, nmo, ortho, moa, "SOLVE")
2803 : CASE (cholesky_inverse)
2804 0 : CALL cp_fm_cholesky_restore(work, nmo, ortho, moa, "MULTIPLY")
2805 : CASE (cholesky_off)
2806 590 : 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 1132 : IF (level_shift_loc /= 0.0_dp) THEN
2814 60 : DO imo = homob + 1, homoa
2815 60 : eiga(imo) = eiga(imo) - 0.5_dp*level_shift_loc
2816 : END DO
2817 220 : DO imo = homoa + 1, nmo
2818 220 : eiga(imo) = eiga(imo) - level_shift_loc
2819 : END DO
2820 : END IF
2821 :
2822 : ! Update also the beta MO set
2823 :
2824 34952 : eigb(:) = eiga(:)
2825 1132 : 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 1132 : CALL calculate_density_matrix(mos(1), scf_env%p_mix_new(1, 1)%matrix)
2831 1132 : CALL calculate_density_matrix(mos(2), scf_env%p_mix_new(2, 1)%matrix)
2832 :
2833 1132 : CALL timestop(handle)
2834 :
2835 1132 : 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 38 : 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 38 : 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 76 : logger => cp_get_default_logger()
2873 38 : CALL timeset(routineN, handle)
2874 :
2875 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%LANCZOS", &
2876 38 : extension=".scfLog")
2877 :
2878 38 : my_check_moconv_only = .FALSE.
2879 38 : IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
2880 :
2881 38 : NULLIFY (mo_coeff, ortho, work, ks)
2882 38 : NULLIFY (mo_eigenvalues)
2883 38 : NULLIFY (c0, c1)
2884 :
2885 38 : IF (scf_env%cholesky_method == cholesky_inverse) THEN
2886 38 : ortho => scf_env%ortho_m1
2887 : ELSE
2888 0 : ortho => scf_env%ortho
2889 : END IF
2890 38 : work => scf_env%scf_work2
2891 :
2892 76 : DO ispin = 1, SIZE(matrix_ks)
2893 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, &
2894 76 : scf_env%scf_work1(ispin))
2895 : END DO
2896 :
2897 38 : IF (scf_env%mixing_method == 1) THEN
2898 0 : scf_env%iter_param = scf_env%p_mix_alpha
2899 0 : scf_env%iter_method = "P_Mix/Lanczos"
2900 : ELSE
2901 : ! scf_env%iter_param = scf_env%mixing_store%alpha
2902 38 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Lanc."
2903 : END IF
2904 :
2905 76 : DO ispin = 1, SIZE(matrix_ks)
2906 :
2907 38 : ks => scf_env%scf_work1(ispin)
2908 38 : 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 38 : mo_coeff=mo_coeff)
2916 :
2917 38 : NULLIFY (c0, c1)
2918 38 : c0 => scf_env%krylov_space%mo_conv(ispin)
2919 38 : c1 => scf_env%krylov_space%mo_refine(ispin)
2920 38 : SELECT CASE (scf_env%cholesky_method)
2921 : CASE (cholesky_reduce)
2922 0 : CALL cp_fm_cholesky_reduce(ks, ortho)
2923 0 : CALL cp_fm_uplo_to_full(ks, work)
2924 0 : 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 0 : "SOLVE", pos="RIGHT")
2928 : CALL cp_fm_cholesky_restore(work, nao, ortho, ks, &
2929 0 : "SOLVE", pos="LEFT", transa="T")
2930 0 : 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 38 : "MULTIPLY", pos="RIGHT")
2934 : CALL cp_fm_cholesky_restore(work, nao, ortho, ks, &
2935 38 : "MULTIPLY", pos="LEFT", transa="T")
2936 76 : CALL cp_fm_cholesky_restore(mo_coeff, nmo, ortho, c0, "SOLVE")
2937 : END SELECT
2938 :
2939 38 : scf_env%krylov_space%nmo_nc = nmo
2940 38 : scf_env%krylov_space%nmo_conv = 0
2941 :
2942 38 : t1 = m_walltime()
2943 38 : IF (output_unit > 0) THEN
2944 0 : 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 0 : " Spin ", " Cycle ", &
2947 0 : " conv. MOS ", " B2MAX ", " B2MIN ", " Time", REPEAT("-", 60)
2948 : END IF
2949 38 : eps_iter = MAX(scf_env%krylov_space%eps_conv, scf_env%krylov_space%eps_adapt*scf_env%iter_delta)
2950 38 : iter = 0
2951 38 : converged = .FALSE.
2952 : !Check convergence of MOS
2953 114 : IF (my_check_moconv_only) THEN
2954 :
2955 : CALL lanczos_refinement(scf_env%krylov_space, ks, c0, c1, mo_eigenvalues, &
2956 0 : nao, eps_iter, ispin, check_moconv_only=my_check_moconv_only)
2957 0 : t2 = m_walltime()
2958 0 : 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 0 : ispin, iter, scf_env%krylov_space%nmo_conv, &
2961 0 : 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 620 : DO iter = 1, scf_env%krylov_space%max_iter
2968 : CALL lanczos_refinement_2v(scf_env%krylov_space, ks, c0, c1, mo_eigenvalues, &
2969 592 : nao, eps_iter, ispin)
2970 592 : t2 = m_walltime()
2971 592 : 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 0 : ispin, iter, scf_env%krylov_space%nmo_conv, &
2974 0 : scf_env%krylov_space%max_res_norm, scf_env%krylov_space%min_res_norm, t2 - t1
2975 : END IF
2976 592 : t1 = m_walltime()
2977 620 : IF (scf_env%krylov_space%max_res_norm < eps_iter) THEN
2978 10 : converged = .TRUE.
2979 10 : IF (output_unit > 0) WRITE (output_unit, *) &
2980 0 : " Reached convergence in ", iter, " iterations "
2981 : EXIT
2982 : END IF
2983 : END DO
2984 :
2985 38 : IF (.NOT. converged .AND. output_unit > 0) THEN
2986 : WRITE (output_unit, "(T4, A)") " WARNING Lanczos refinement could "// &
2987 0 : "not converge all the mos:"
2988 0 : WRITE (output_unit, "(T40,A,T70,I10)") " number of not converged mos ", &
2989 0 : scf_env%krylov_space%nmo_nc
2990 0 : WRITE (output_unit, "(T40,A,T70,E10.2)") " max norm of the residual ", &
2991 0 : 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 38 : IF (scf_env%cholesky_method == cholesky_inverse) THEN
3011 38 : CALL cp_fm_cholesky_restore(c0, nmo, ortho, mo_coeff, "MULTIPLY")
3012 : ELSE
3013 0 : CALL cp_fm_cholesky_restore(c0, nmo, ortho, mo_coeff, "SOLVE")
3014 : END IF
3015 :
3016 38 : IF (.NOT. scf_control%gce%do_gce) THEN
3017 : CALL set_mo_occupation(mo_set=mos(ispin), &
3018 38 : smear=scf_control%smear)
3019 : ELSE
3020 : CALL set_mo_occupation(mo_set=mos(ispin), &
3021 : smear=scf_control%smear, &
3022 0 : gce=scf_control%gce)
3023 : END IF
3024 :
3025 : ! does not yet handle k-points
3026 : CALL calculate_density_matrix(mos(ispin), &
3027 38 : scf_env%p_mix_new(ispin, 1)%matrix)
3028 : END IF
3029 : END DO ! ispin
3030 :
3031 38 : IF (output_unit > 0) THEN
3032 0 : 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 38 : "PRINT%LANCZOS")
3037 :
3038 38 : CALL timestop(handle)
3039 :
3040 38 : 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 94 : 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 94 : logger => cp_get_default_logger()
3075 94 : CALL timeset(routineN, handle)
3076 :
3077 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DAVIDSON", &
3078 94 : extension=".scfLog")
3079 :
3080 94 : IF (output_unit > 0) THEN
3081 0 : WRITE (output_unit, "(/T15,A)") '<<<<<<<<< DAVIDSON ITERATIONS <<<<<<<<<<'
3082 : END IF
3083 :
3084 94 : IF (scf_env%mixing_method == 1) THEN
3085 0 : scf_env%iter_param = scf_env%p_mix_alpha
3086 0 : scf_env%iter_method = "P_Mix/Dav."
3087 : ELSE
3088 94 : scf_env%iter_param = scf_env%mixing_store%alpha
3089 94 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Dav."
3090 : END IF
3091 :
3092 94 : my_check_moconv_only = .FALSE.
3093 94 : IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
3094 94 : do_prec = .FALSE.
3095 94 : 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 84 : do_prec = .TRUE.
3098 : END IF
3099 :
3100 94 : nspins = SIZE(matrix_ks)
3101 :
3102 94 : 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 18 : 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 18 : full_mo_set=.TRUE.)
3112 : END IF
3113 :
3114 198 : DO ispin = 1, nspins
3115 198 : IF (scf_env%block_davidson_env(ispin)%use_sparse_mos) THEN
3116 64 : IF (.NOT. do_prec) THEN
3117 : CALL generate_extended_space_sparse(scf_env%block_davidson_env(ispin), mos(ispin), &
3118 8 : 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 56 : scf_env%ot_preconditioner(ispin)%preconditioner)
3123 : END IF
3124 :
3125 : ELSE
3126 40 : IF (.NOT. do_prec) THEN
3127 : CALL generate_extended_space(scf_env%block_davidson_env(ispin), mos(ispin), &
3128 4 : 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 36 : scf_env%ot_preconditioner(ispin)%preconditioner)
3133 : END IF
3134 : END IF
3135 : END DO !ispin
3136 :
3137 94 : IF (.NOT. scf_control%gce%do_gce) THEN
3138 : CALL set_mo_occupation(mo_array=mos, &
3139 94 : smear=scf_control%smear)
3140 : ELSE
3141 : CALL set_mo_occupation(mo_array=mos, &
3142 : smear=scf_control%smear, &
3143 0 : gce=scf_control%gce)
3144 : END IF
3145 :
3146 198 : DO ispin = 1, nspins
3147 : ! does not yet handle k-points
3148 : CALL calculate_density_matrix(mos(ispin), &
3149 198 : scf_env%p_mix_new(ispin, 1)%matrix)
3150 : END DO
3151 :
3152 94 : IF (output_unit > 0) THEN
3153 0 : 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 94 : "PRINT%DAVIDSON")
3158 :
3159 94 : CALL timestop(handle)
3160 :
3161 94 : 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 932 : 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 932 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
3193 932 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3194 : LOGICAL :: do_prec, use_real_wfn
3195 : REAL(KIND=dp) :: eps_iter, eps_iter_empty
3196 932 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
3197 : TYPE(cp_cfm_type) :: cksmat, csmat, csmat_base
3198 932 : 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 932 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
3201 : TYPE(cp_logger_type), POINTER :: logger
3202 : TYPE(kpoint_env_type), POINTER :: kp
3203 932 : 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 932 : 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 1864 : logger => cp_get_default_logger()
3221 932 : CALL timeset(routineN, handle)
3222 :
3223 932 : 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 932 : extension=".scfLog")
3226 932 : IF (output_unit > 0) THEN
3227 0 : 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 932 : cell_to_index=cell_to_index, mpools=kp_mpools)
3233 932 : IF (use_real_wfn) THEN
3234 0 : CPABORT('Davidson Kp requires complex K-point MOS')
3235 : END IF
3236 932 : para_env => kpoints%blacs_env_all%para_env
3237 932 : kplocal = kp_range(2) - kp_range(1) + 1
3238 932 : nspin = SIZE(matrix_ks_kp, 1)
3239 932 : nchannels = kplocal*nspin
3240 :
3241 932 : IF (scf_env%mixing_method == 1) THEN
3242 684 : scf_env%iter_param = scf_env%p_mix_alpha
3243 684 : scf_env%iter_method = "P_Mix/Dav."
3244 248 : ELSE IF (scf_env%mixing_method > 1) THEN
3245 248 : scf_env%iter_param = scf_env%mixing_store%alpha
3246 248 : 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 932 : first_prec = MAX(scf_env%block_davidson_env(1)%first_prec, 2)
3257 932 : do_prec = .FALSE.
3258 932 : IF (scf_env%block_davidson_env(1)%prec_type /= ot_precond_none .AND. &
3259 : scf_env%iter_count >= first_prec) THEN
3260 838 : 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 932 : eps_iter = scf_env%block_davidson_env(1)%eps_iter
3277 932 : eps_iter_empty = eps_iter
3278 932 : 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 172 : 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 2 : 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 172 : eps_iter_cap)), scf_env%block_davidson_env(1)%eps_iter_prev)
3301 172 : 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 172 : scf_env%block_davidson_env(1)%eps_iter_prev = eps_iter
3305 172 : 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 932 : 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 0 : iw = cp_logger_get_default_unit_nr()
3312 0 : 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 0 : 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 932 : CALL get_qs_env(qs_env, mpools=mpools)
3322 932 : CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
3323 4660 : ALLOCATE (fmwork(4))
3324 4660 : DO iwork = 1, 4
3325 4660 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmwork(iwork))
3326 : END DO
3327 932 : CALL kpoint_operator_context_create(op_ctx, kpoints, fmwork(1), nspin)
3328 932 : CALL mpools_get(kp_mpools, ao_ao_fm_pools=kp_ao_ao_fm_pools)
3329 932 : matrix_struct => fm_pool_get_el_struct(kp_ao_ao_fm_pools(1)%pool)
3330 932 : CALL cp_cfm_create(cksmat, matrix_struct, name="dav_hk")
3331 932 : CALL cp_cfm_create(csmat, matrix_struct, name="dav_sk")
3332 932 : CALL cp_cfm_create(csmat_base, matrix_struct, name="dav_sk_base")
3333 :
3334 932 : 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 58 : qs_env, scf_env, matrix_ks_kp, matrix_s_kp, kpoints, op_ctx)
3339 : END IF
3340 :
3341 2936 : DO local_kpoint = 1, kplocal
3342 2004 : 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 2004 : cmat=csmat_base, matrix_row=kpoint_spin_free)
3349 :
3350 5016 : DO ispin = 1, nspin
3351 : CALL kpoint_operator_get(op_ctx, local_kpoint, ispin, matrix_ks_kp, &
3352 2080 : cmat=cksmat, matrix_row=ispin)
3353 2080 : CALL cp_cfm_to_cfm(csmat_base, csmat)
3354 :
3355 4084 : 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 68 : CALL cold_start_mo_init(kp, ispin, cksmat, csmat)
3360 : ELSE
3361 2012 : channel = (local_kpoint - 1)*nspin + ispin
3362 2012 : 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 1884 : 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 128 : 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 932 : CALL kpoint_set_mo_occupation(kpoints, scf_control%smear)
3378 932 : CALL kpoint_density_matrices(kpoints)
3379 : CALL kpoint_density_transform(kpoints, scf_env%p_mix_new, .FALSE., &
3380 932 : matrix_s_kp(1, 1)%matrix, sab_nl, fmwork, overlap_rs=matrix_s_kp)
3381 :
3382 932 : CALL cp_cfm_release(cksmat)
3383 932 : CALL cp_cfm_release(csmat)
3384 932 : CALL cp_cfm_release(csmat_base)
3385 932 : CALL kpoint_operator_context_release(op_ctx)
3386 4660 : DO iwork = 1, 4
3387 4660 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmwork(iwork))
3388 : END DO
3389 932 : DEALLOCATE (fmwork)
3390 :
3391 932 : IF (output_unit > 0) THEN
3392 0 : WRITE (output_unit, "(T15,A/)") '<<<<<<<<< END DAVIDSON ITERATION <<<<<<<<<<'
3393 : END IF
3394 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
3395 932 : "PRINT%DAVIDSON")
3396 :
3397 932 : CALL timestop(handle)
3398 :
3399 2796 : 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 204 : 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 68 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
3420 : TYPE(cp_cfm_type) :: cmos, cwork
3421 : TYPE(cp_fm_type), POINTER :: imos, rmos
3422 :
3423 68 : CALL timeset(routineN, handle)
3424 :
3425 68 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
3426 68 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
3427 68 : CALL cp_cfm_create(cmos, rmos%matrix_struct, name="dav_cold_cmos")
3428 68 : CALL cp_cfm_create(cwork, cksmat%matrix_struct, name="dav_cold_cwork")
3429 68 : CALL cp_cfm_geeig(cksmat, csmat, cmos, eigenvalues, cwork)
3430 552 : kp%mos(2, ispin)%eigenvalues = eigenvalues
3431 68 : CALL cp_cfm_to_fm(cmos, rmos, imos)
3432 68 : CALL cp_cfm_release(cmos)
3433 68 : CALL cp_cfm_release(cwork)
3434 :
3435 68 : CALL timestop(handle)
3436 :
3437 68 : 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 58 : 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 58 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
3465 58 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3466 : REAL(KIND=dp) :: energy_gap
3467 58 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
3468 58 : 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 58 : 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 58 : POINTER :: sab_nl
3479 :
3480 58 : 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 58 : cell_to_index=cell_to_index)
3484 58 : 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 58 : blacs_env => kpoints%blacs_env
3488 58 : para_env => kpoints%para_env_kp
3489 58 : kplocal = kp_range(2) - kp_range(1) + 1
3490 58 : nspin = SIZE(matrix_ks_kp, 1)
3491 58 : nchannels = kplocal*nspin
3492 58 : prec_type = scf_env%block_davidson_env(1)%prec_type
3493 58 : solver_type = scf_env%block_davidson_env(1)%solver_type
3494 58 : 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 58 : IF (ASSOCIATED(scf_env%ot_preconditioner)) THEN
3501 128 : DO channel = 1, SIZE(scf_env%ot_preconditioner)
3502 128 : IF (ASSOCIATED(scf_env%ot_preconditioner(channel)%preconditioner)) THEN
3503 88 : CALL destroy_preconditioner(scf_env%ot_preconditioner(channel)%preconditioner)
3504 88 : DEALLOCATE (scf_env%ot_preconditioner(channel)%preconditioner)
3505 : END IF
3506 : END DO
3507 40 : DEALLOCATE (scf_env%ot_preconditioner)
3508 40 : NULLIFY (scf_env%ot_preconditioner)
3509 : END IF
3510 306 : ALLOCATE (scf_env%ot_preconditioner(nchannels))
3511 190 : DO channel = 1, nchannels
3512 660 : ALLOCATE (scf_env%ot_preconditioner(channel)%preconditioner)
3513 : CALL init_preconditioner(scf_env%ot_preconditioner(channel)%preconditioner, &
3514 190 : para_env, blacs_env)
3515 : END DO
3516 :
3517 58 : NULLIFY (eigenvalues, matrix_h_im, matrix_h_re, matrix_s_im, matrix_s_re, &
3518 58 : matrix_t_im, matrix_t_re, mo_coeff_im, mo_coeff_re)
3519 58 : IF (prec_type == ot_precond_full_kinetic) THEN
3520 28 : CPASSERT(ASSOCIATED(matrix_t_kp))
3521 : END IF
3522 :
3523 186 : DO local_kpoint = 1, kplocal
3524 128 : 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 128 : CALL dbcsr_init_p(matrix_s_re)
3531 128 : 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 128 : matrix_row=kpoint_spin_free)
3535 :
3536 128 : IF (prec_type == ot_precond_full_kinetic) THEN
3537 56 : CALL dbcsr_init_p(matrix_t_re)
3538 56 : 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 56 : matrix_row=kpoint_spin_free)
3542 : END IF
3543 :
3544 260 : DO ispin = 1, nspin
3545 132 : channel = (local_kpoint - 1)*nspin + ispin
3546 260 : IF (prec_type == ot_precond_full_all .OR. prec_type == ot_precond_full_single_inverse) THEN
3547 60 : CALL dbcsr_init_p(matrix_h_re)
3548 60 : CALL dbcsr_init_p(matrix_h_im)
3549 : CALL kpoint_operator_get(op_ctx, local_kpoint, ispin, matrix_ks_kp, &
3550 60 : 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 60 : nmo=nmo)
3559 60 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=mo_coeff_im)
3560 60 : CALL cp_cfm_create(matrix_c, mo_coeff_re%matrix_struct, name="dav pc C")
3561 60 : CALL cp_cfm_scale_and_add_fm(z_zero, matrix_c, z_one, mo_coeff_re)
3562 60 : 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 60 : '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 60 : 'dav pc S', matrix_s)
3569 60 : 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 44 : energy_gap)
3573 : ELSE
3574 : CALL make_complex_full_single_inverse( &
3575 : scf_env%ot_preconditioner(channel)%preconditioner, &
3576 16 : matrix_c, matrix_h, matrix_s, energy_gap)
3577 : END IF
3578 60 : CALL cp_cfm_release(matrix_c)
3579 60 : CALL cp_cfm_release(matrix_h)
3580 60 : CALL cp_cfm_release(matrix_s)
3581 60 : CALL dbcsr_release_p(matrix_h_re)
3582 60 : CALL dbcsr_release_p(matrix_h_im)
3583 72 : ELSE IF (prec_type == ot_precond_full_kinetic) THEN
3584 : CALL make_preconditioner_complex_full_kinetic( &
3585 : scf_env%ot_preconditioner(channel)%preconditioner, &
3586 56 : matrix_t_re, matrix_t_im, matrix_s_re, matrix_s_im, energy_gap, solver_type)
3587 16 : 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 16 : matrix_s_re, matrix_s_im, solver_type)
3591 : END IF
3592 : END DO
3593 :
3594 128 : CALL dbcsr_release_p(matrix_s_re)
3595 128 : CALL dbcsr_release_p(matrix_s_im)
3596 186 : IF (prec_type == ot_precond_full_kinetic) THEN
3597 56 : CALL dbcsr_release_p(matrix_t_re)
3598 56 : CALL dbcsr_release_p(matrix_t_im)
3599 : END IF
3600 : END DO
3601 :
3602 58 : CALL timestop(handle)
3603 :
3604 58 : 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 12 : 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 12 : 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 12 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
3630 12 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3631 : LOGICAL :: my_kpgrp, use_real_wfn
3632 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
3633 12 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
3634 12 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
3635 : TYPE(cp_cfm_type) :: csmat, cwork
3636 12 : 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 12 : 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 12 : POINTER :: sab_nl
3645 : TYPE(qs_matrix_pools_type), POINTER :: mpools
3646 :
3647 12 : CALL timeset(routineN, handle)
3648 :
3649 12 : 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 12 : cell_to_index=cell_to_index)
3653 12 : CPASSERT(ASSOCIATED(sab_nl))
3654 12 : kplocal = kp_range(2) - kp_range(1) + 1
3655 :
3656 : ! allocate some work matrices
3657 12 : ALLOCATE (rmatrix, cmatrix, tmpmat)
3658 : CALL dbcsr_create(rmatrix, template=matrix_s(1, 1)%matrix, &
3659 12 : matrix_type=dbcsr_type_symmetric)
3660 : CALL dbcsr_create(cmatrix, template=matrix_s(1, 1)%matrix, &
3661 12 : matrix_type=dbcsr_type_antisymmetric)
3662 : CALL dbcsr_create(tmpmat, template=matrix_s(1, 1)%matrix, &
3663 12 : matrix_type=dbcsr_type_no_symmetry)
3664 12 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
3665 12 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
3666 :
3667 : ! fm pools to be used within a kpoint group
3668 12 : CALL get_kpoint_info(kpoints, mpools=mpools)
3669 12 : CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
3670 :
3671 12 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
3672 12 : CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
3673 :
3674 12 : IF (use_real_wfn) THEN
3675 0 : CALL cp_fm_create(rsmat, matrix_struct)
3676 : ELSE
3677 12 : CALL cp_cfm_create(csmat, matrix_struct)
3678 12 : CALL cp_cfm_create(cwork, matrix_struct)
3679 : END IF
3680 :
3681 12 : CALL cp_fm_get_info(fmwork(1), nrow_global=nao)
3682 60 : ALLOCATE (eigenvalues(nao), ceig(nao))
3683 :
3684 12 : 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 12 : CALL kpoint_operator_context_create(op_ctx, kpoints, fmwork(1), 1)
3688 : ! The real-wavefunction path keeps its own transfer bookkeeping.
3689 12 : IF (use_real_wfn) ALLOCATE (info(kplocal*nkp_groups, 1))
3690 :
3691 : ! Setup and start all the communication
3692 12 : indx = 0
3693 268 : DO ikp = 1, kplocal
3694 268 : IF (use_real_wfn) THEN
3695 0 : DO igroup = 1, nkp_groups
3696 : ! number of current kpoint
3697 0 : ik = kp_dist(1, igroup) + ikp - 1
3698 0 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
3699 0 : indx = indx + 1
3700 0 : CALL dbcsr_set(rmatrix, 0.0_dp)
3701 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_s, ispin=1, &
3702 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
3703 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
3704 0 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
3705 : ! transfer to the kpoint group with a new blacs environment
3706 0 : IF (my_kpgrp) THEN
3707 0 : CALL cp_fm_start_copy_general(fmwork(1), rsmat, para_env, info(indx, 1))
3708 : ELSE
3709 0 : 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 256 : use_grid=.FALSE., matrix_row=kpoint_spin_free)
3717 : END IF
3718 : END DO
3719 :
3720 : ! Finish communication then diagonalise in each group
3721 12 : indx = 0
3722 268 : DO ikp = 1, kplocal
3723 256 : IF (use_real_wfn) THEN
3724 0 : DO igroup = 1, nkp_groups
3725 : ! number of current kpoint
3726 0 : ik = kp_dist(1, igroup) + ikp - 1
3727 0 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
3728 0 : indx = indx + 1
3729 0 : IF (my_kpgrp) THEN
3730 0 : CALL cp_fm_finish_copy_general(rsmat, info(indx, 1))
3731 : END IF
3732 : END DO
3733 : ELSE
3734 256 : 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 256 : kp => kpoints%kp_env(ikp)%kpoint_env
3740 256 : IF (use_real_wfn) THEN
3741 0 : CALL choose_eigv_solver(rsmat, fmlocal, eigenvalues)
3742 : ELSE
3743 256 : CALL cp_cfm_heevd(csmat, cwork, eigenvalues)
3744 : END IF
3745 2680 : CPASSERT(ALL(eigenvalues(1:nao) >= 0.0_dp))
3746 268 : IF (use_real_wfn) THEN
3747 0 : CALL cp_fm_release(kp%shalf)
3748 0 : CALL cp_fm_create(kp%shalf, matrix_struct)
3749 0 : eigenvalues(1:nao) = SQRT(eigenvalues(1:nao))
3750 0 : CALL cp_fm_to_fm(fmlocal, rsmat)
3751 0 : CALL cp_fm_column_scale(rsmat, eigenvalues)
3752 : CALL parallel_gemm("N", "T", nao, nao, nao, 1.0_dp, rsmat, fmlocal, &
3753 0 : 0.0_dp, kp%shalf)
3754 : ELSE
3755 256 : CALL cp_cfm_release(kp%cshalf)
3756 256 : CALL cp_cfm_create(kp%cshalf, matrix_struct)
3757 2680 : ceig(1:nao) = SQRT(eigenvalues(1:nao))
3758 256 : CALL cp_cfm_to_cfm(cwork, csmat)
3759 256 : CALL cp_cfm_column_scale(csmat, ceig)
3760 : CALL parallel_gemm("N", "C", nao, nao, nao, cone, csmat, cwork, &
3761 256 : 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 12 : IF (use_real_wfn) THEN
3769 0 : indx = 0
3770 0 : DO ikp = 1, kplocal
3771 0 : DO igroup = 1, nkp_groups
3772 : ! number of current kpoint
3773 0 : ik = kp_dist(1, igroup) + ikp - 1
3774 0 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
3775 0 : indx = indx + 1
3776 0 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
3777 : END DO
3778 : END DO
3779 : END IF
3780 :
3781 12 : CALL kpoint_operator_context_release(op_ctx)
3782 :
3783 : ! All done
3784 12 : IF (ALLOCATED(info)) DEALLOCATE (info)
3785 12 : DEALLOCATE (eigenvalues, ceig)
3786 :
3787 12 : CALL dbcsr_deallocate_matrix(rmatrix)
3788 12 : CALL dbcsr_deallocate_matrix(cmatrix)
3789 12 : CALL dbcsr_deallocate_matrix(tmpmat)
3790 :
3791 12 : IF (use_real_wfn) THEN
3792 0 : CALL cp_fm_release(rsmat)
3793 : ELSE
3794 12 : CALL cp_cfm_release(csmat)
3795 12 : CALL cp_cfm_release(cwork)
3796 : END IF
3797 12 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
3798 :
3799 12 : CALL timestop(handle)
3800 :
3801 48 : END SUBROUTINE diag_kp_smat
3802 :
3803 : END MODULE qs_scf_diagonalization
|
