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