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