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