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_cfm_basic_linalg, ONLY: cp_cfm_column_scale,&
18 : cp_cfm_scale_and_add,&
19 : cp_cfm_scale_and_add_fm
20 : USE cp_cfm_diag, ONLY: cp_cfm_geeig,&
21 : cp_cfm_geeig_canon,&
22 : cp_cfm_heevd
23 : USE cp_cfm_types, ONLY: cp_cfm_create,&
24 : cp_cfm_release,&
25 : cp_cfm_set_all,&
26 : cp_cfm_to_cfm,&
27 : cp_cfm_to_fm,&
28 : cp_cfm_type
29 : USE cp_control_types, ONLY: dft_control_type,&
30 : hairy_probes_type
31 : USE cp_dbcsr_api, ONLY: &
32 : dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_p_type, &
33 : dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
34 : dbcsr_type_symmetric
35 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
36 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
37 : copy_fm_to_dbcsr,&
38 : cp_dbcsr_sm_fm_multiply,&
39 : dbcsr_allocate_matrix_set
40 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
41 : cp_fm_symm,&
42 : cp_fm_uplo_to_full
43 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_reduce,&
44 : cp_fm_cholesky_restore
45 : USE cp_fm_diag, ONLY: FM_DIAG_TYPE_CUSOLVER,&
46 : choose_eigv_solver,&
47 : cp_fm_geeig,&
48 : cp_fm_geeig_canon,&
49 : cusolver_generalized,&
50 : diag_type
51 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
52 : fm_pool_create_fm,&
53 : fm_pool_give_back_fm
54 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
55 : cp_fm_struct_release,&
56 : cp_fm_struct_type
57 : USE cp_fm_types, ONLY: &
58 : copy_info_type, cp_fm_add_to_element, cp_fm_cleanup_copy_general, cp_fm_create, &
59 : cp_fm_finish_copy_general, cp_fm_get_info, cp_fm_release, cp_fm_start_copy_general, &
60 : cp_fm_to_fm, cp_fm_type
61 : USE cp_log_handling, ONLY: cp_get_default_logger,&
62 : cp_logger_type
63 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
64 : cp_print_key_unit_nr
65 : USE input_constants, ONLY: &
66 : cholesky_dbcsr, cholesky_inverse, cholesky_off, cholesky_reduce, cholesky_restore, &
67 : core_guess, general_roks, high_spin_roks, restart_guess
68 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
69 : section_vals_type
70 : USE kinds, ONLY: dp
71 : USE kpoint_methods, ONLY: kpoint_density_matrices,&
72 : kpoint_density_transform,&
73 : kpoint_set_mo_occupation,&
74 : rskp_transform
75 : USE kpoint_types, ONLY: get_kpoint_info,&
76 : kpoint_env_type,&
77 : kpoint_type
78 : USE machine, ONLY: m_flush,&
79 : m_walltime
80 : USE mathconstants, ONLY: gaussi,&
81 : z_one,&
82 : z_zero
83 : USE message_passing, ONLY: mp_para_env_type
84 : USE parallel_gemm_api, ONLY: parallel_gemm
85 : USE preconditioner, ONLY: prepare_preconditioner,&
86 : restart_preconditioner
87 : USE qs_density_matrices, ONLY: calculate_density_matrix
88 : USE qs_density_mixing_types, ONLY: direct_mixing_nr,&
89 : gspace_mixing_nr
90 : USE qs_diis, ONLY: qs_diis_b_calc_err_kp,&
91 : qs_diis_b_info_kp,&
92 : qs_diis_b_step,&
93 : qs_diis_b_step_kp
94 : USE qs_energy_types, ONLY: qs_energy_type
95 : USE qs_environment_types, ONLY: get_qs_env,&
96 : qs_environment_type
97 : USE qs_gspace_mixing, ONLY: gspace_mixing
98 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
99 : USE qs_ks_types, ONLY: qs_ks_did_change,&
100 : qs_ks_env_type
101 : USE qs_matrix_pools, ONLY: mpools_get,&
102 : qs_matrix_pools_type
103 : USE qs_mixing_utils, ONLY: charge_mixing_init,&
104 : mixing_allocate,&
105 : mixing_init,&
106 : self_consistency_check
107 : USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
108 : USE qs_mo_occupation, ONLY: set_mo_occupation
109 : USE qs_mo_types, ONLY: get_mo_set,&
110 : mo_set_type
111 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
112 : USE qs_ot_eigensolver, ONLY: ot_eigensolver
113 : USE qs_rho_atom_types, ONLY: rho_atom_type
114 : USE qs_rho_methods, ONLY: qs_rho_update_rho
115 : USE qs_rho_types, ONLY: qs_rho_get,&
116 : qs_rho_type
117 : USE qs_scf_block_davidson, ONLY: generate_extended_space,&
118 : generate_extended_space_sparse
119 : USE qs_scf_lanczos, ONLY: lanczos_refinement,&
120 : lanczos_refinement_2v
121 : USE qs_scf_methods, ONLY: combine_ks_matrices,&
122 : eigensolver,&
123 : eigensolver_dbcsr,&
124 : eigensolver_generalized,&
125 : eigensolver_simple,&
126 : eigensolver_symm,&
127 : scf_env_density_mixing
128 : USE qs_scf_types, ONLY: qs_scf_env_type,&
129 : subspace_env_type
130 : USE scf_control_types, ONLY: scf_control_type
131 : #include "./base/base_uses.f90"
132 :
133 : IMPLICIT NONE
134 :
135 : PRIVATE
136 :
137 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_diagonalization'
138 :
139 : PUBLIC :: do_general_diag, do_general_diag_kp, do_roks_diag, &
140 : do_special_diag, do_ot_diag, do_block_davidson_diag, &
141 : do_block_krylov_diag, do_scf_diag_subspace, diag_subspace_allocate, &
142 : general_eigenproblem, diag_kp_smat
143 :
144 : CONTAINS
145 :
146 : ! **************************************************************************************************
147 : !> \brief the inner loop of scf, specific to diagonalization with S matrix
148 : !> basically, in goes the ks matrix out goes a new p matrix
149 : !> \param scf_env ...
150 : !> \param mos ...
151 : !> \param matrix_ks ...
152 : !> \param matrix_s ...
153 : !> \param scf_control ...
154 : !> \param scf_section ...
155 : !> \param diis_step ...
156 : !> \par History
157 : !> 03.2006 created [Joost VandeVondele]
158 : ! **************************************************************************************************
159 :
160 77173 : SUBROUTINE general_eigenproblem(scf_env, mos, matrix_ks, &
161 : matrix_s, scf_control, scf_section, &
162 : diis_step)
163 :
164 : TYPE(qs_scf_env_type), POINTER :: scf_env
165 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
166 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
167 : TYPE(scf_control_type), POINTER :: scf_control
168 : TYPE(section_vals_type), POINTER :: scf_section
169 : LOGICAL, INTENT(INOUT) :: diis_step
170 :
171 : INTEGER :: ispin, nspin
172 : LOGICAL :: do_level_shift, owns_ortho, use_jacobi
173 : REAL(KIND=dp) :: diis_error, eps_diis
174 : TYPE(cp_fm_type), POINTER :: ortho
175 : TYPE(dbcsr_type), POINTER :: ortho_dbcsr
176 :
177 77173 : nspin = SIZE(matrix_ks)
178 77173 : NULLIFY (ortho, ortho_dbcsr)
179 :
180 164266 : DO ispin = 1, nspin
181 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, &
182 164266 : scf_env%scf_work1(ispin))
183 : END DO
184 :
185 77173 : eps_diis = scf_control%eps_diis
186 :
187 77173 : IF (scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis) THEN
188 : CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
189 : scf_env%scf_work2, scf_env%iter_delta, diis_error, diis_step, &
190 : eps_diis, scf_control%nmixing, &
191 : s_matrix=matrix_s, &
192 64397 : scf_section=scf_section)
193 : ELSE
194 12776 : diis_step = .FALSE.
195 : END IF
196 :
197 : do_level_shift = ((scf_control%level_shift /= 0.0_dp) .AND. &
198 : ((scf_control%density_guess == core_guess) .OR. &
199 77173 : (scf_env%iter_count > 1)))
200 :
201 77173 : IF ((scf_env%iter_count > 1) .AND. &
202 : (scf_env%iter_delta < scf_control%diagonalization%eps_jacobi)) THEN
203 0 : use_jacobi = .TRUE.
204 : ELSE
205 77173 : use_jacobi = .FALSE.
206 : END IF
207 :
208 77173 : IF (diis_step) THEN
209 45595 : scf_env%iter_param = diis_error
210 45595 : IF (use_jacobi) THEN
211 0 : scf_env%iter_method = "DIIS/Jacobi"
212 : ELSE
213 45595 : scf_env%iter_method = "DIIS/Diag."
214 : END IF
215 : ELSE
216 31578 : IF (scf_env%mixing_method == 0) THEN
217 0 : scf_env%iter_method = "NoMix/Diag."
218 31578 : ELSE IF (scf_env%mixing_method == 1) THEN
219 30164 : scf_env%iter_param = scf_env%p_mix_alpha
220 30164 : IF (use_jacobi) THEN
221 0 : scf_env%iter_method = "P_Mix/Jacobi"
222 : ELSE
223 30164 : scf_env%iter_method = "P_Mix/Diag."
224 : END IF
225 1414 : ELSEIF (scf_env%mixing_method > 1) THEN
226 1414 : scf_env%iter_param = scf_env%mixing_store%alpha
227 1414 : IF (use_jacobi) THEN
228 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Jacobi"
229 : ELSE
230 1414 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Diag."
231 : END IF
232 : END IF
233 : END IF
234 :
235 77173 : IF (scf_env%cholesky_method == cholesky_dbcsr) THEN
236 1064 : ortho_dbcsr => scf_env%ortho_dbcsr
237 3182 : DO ispin = 1, nspin
238 : CALL eigensolver_dbcsr(matrix_ks=matrix_ks(ispin)%matrix, matrix_ks_fm=scf_env%scf_work1(ispin), &
239 : mo_set=mos(ispin), &
240 : ortho_dbcsr=ortho_dbcsr, &
241 3182 : ksbuf1=scf_env%buf1_dbcsr, ksbuf2=scf_env%buf2_dbcsr)
242 : END DO
243 :
244 76109 : ELSE IF (scf_env%cholesky_method > cholesky_off) THEN
245 75683 : IF (scf_env%cholesky_method == cholesky_inverse) THEN
246 144 : ortho => scf_env%ortho_m1
247 : ELSE
248 75539 : ortho => scf_env%ortho
249 : END IF
250 :
251 75683 : owns_ortho = .FALSE.
252 75683 : IF (.NOT. ASSOCIATED(ortho)) THEN
253 0 : ALLOCATE (ortho)
254 0 : owns_ortho = .TRUE.
255 : END IF
256 :
257 160150 : DO ispin = 1, nspin
258 160150 : IF (diag_type == FM_DIAG_TYPE_CUSOLVER .AND. cusolver_generalized .AND. .NOT. do_level_shift) THEN
259 : CALL eigensolver_generalized(matrix_ks_fm=scf_env%scf_work1(ispin), &
260 : matrix_s=matrix_s(ispin)%matrix, &
261 : mo_set=mos(ispin), &
262 0 : work=scf_env%scf_work2)
263 : ELSE
264 84467 : IF (do_level_shift) THEN
265 : CALL eigensolver(matrix_ks_fm=scf_env%scf_work1(ispin), &
266 : mo_set=mos(ispin), &
267 : ortho=ortho, &
268 : work=scf_env%scf_work2, &
269 : cholesky_method=scf_env%cholesky_method, &
270 : do_level_shift=do_level_shift, &
271 : level_shift=scf_control%level_shift, &
272 : matrix_u_fm=scf_env%ortho, &
273 124 : use_jacobi=use_jacobi)
274 : ELSE
275 : CALL eigensolver(matrix_ks_fm=scf_env%scf_work1(ispin), &
276 : mo_set=mos(ispin), &
277 : ortho=ortho, &
278 : work=scf_env%scf_work2, &
279 : cholesky_method=scf_env%cholesky_method, &
280 : do_level_shift=do_level_shift, &
281 : level_shift=scf_control%level_shift, &
282 84343 : use_jacobi=use_jacobi)
283 : END IF
284 : END IF
285 : END DO
286 :
287 75683 : IF (owns_ortho) DEALLOCATE (ortho)
288 : ELSE
289 426 : ortho => scf_env%ortho
290 :
291 426 : owns_ortho = .FALSE.
292 426 : IF (.NOT. ASSOCIATED(ortho)) THEN
293 0 : ALLOCATE (ortho)
294 0 : owns_ortho = .TRUE.
295 : END IF
296 :
297 426 : IF (do_level_shift) THEN
298 172 : DO ispin = 1, nspin
299 : IF (ASSOCIATED(scf_env%scf_work1_red) .AND. ASSOCIATED(scf_env%scf_work2_red) &
300 172 : .AND. ASSOCIATED(scf_env%ortho_red) .AND. ASSOCIATED(scf_env%ortho_m1_red)) THEN
301 : CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
302 : mo_set=mos(ispin), &
303 : ortho=ortho, &
304 : work=scf_env%scf_work2, &
305 : do_level_shift=do_level_shift, &
306 : level_shift=scf_control%level_shift, &
307 : matrix_u_fm=scf_env%ortho_m1, &
308 : use_jacobi=use_jacobi, &
309 : jacobi_threshold=scf_control%diagonalization%jacobi_threshold, &
310 : matrix_ks_fm_red=scf_env%scf_work1_red(ispin), &
311 : ortho_red=scf_env%ortho_red, &
312 : work_red=scf_env%scf_work2_red, &
313 86 : matrix_u_fm_red=scf_env%ortho_m1_red)
314 : ELSE
315 : CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
316 : mo_set=mos(ispin), &
317 : ortho=ortho, &
318 : work=scf_env%scf_work2, &
319 : do_level_shift=do_level_shift, &
320 : level_shift=scf_control%level_shift, &
321 : matrix_u_fm=scf_env%ortho_m1, &
322 : use_jacobi=use_jacobi, &
323 0 : jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
324 : END IF
325 : END DO
326 : ELSE
327 762 : DO ispin = 1, nspin
328 : IF (ASSOCIATED(scf_env%scf_work1_red) .AND. ASSOCIATED(scf_env%scf_work2_red) &
329 762 : .AND. ASSOCIATED(scf_env%ortho_red)) THEN
330 : CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
331 : mo_set=mos(ispin), &
332 : ortho=ortho, &
333 : work=scf_env%scf_work2, &
334 : do_level_shift=do_level_shift, &
335 : level_shift=scf_control%level_shift, &
336 : use_jacobi=use_jacobi, &
337 : jacobi_threshold=scf_control%diagonalization%jacobi_threshold, &
338 : matrix_ks_fm_red=scf_env%scf_work1_red(ispin), &
339 : ortho_red=scf_env%ortho_red, &
340 422 : work_red=scf_env%scf_work2_red)
341 : ELSE
342 : CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
343 : mo_set=mos(ispin), &
344 : ortho=ortho, &
345 : work=scf_env%scf_work2, &
346 : do_level_shift=do_level_shift, &
347 : level_shift=scf_control%level_shift, &
348 : use_jacobi=use_jacobi, &
349 0 : jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
350 : END IF
351 : END DO
352 : END IF
353 :
354 426 : IF (owns_ortho) DEALLOCATE (ortho)
355 : END IF
356 :
357 77173 : END SUBROUTINE general_eigenproblem
358 :
359 : ! **************************************************************************************************
360 : !> \brief ...
361 : !> \param scf_env ...
362 : !> \param mos ...
363 : !> \param matrix_ks ...
364 : !> \param matrix_s ...
365 : !> \param scf_control ...
366 : !> \param scf_section ...
367 : !> \param diis_step ...
368 : !> \param probe ...
369 : ! **************************************************************************************************
370 76145 : SUBROUTINE do_general_diag(scf_env, mos, matrix_ks, &
371 : matrix_s, scf_control, scf_section, &
372 : diis_step, probe)
373 :
374 : TYPE(qs_scf_env_type), POINTER :: scf_env
375 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
376 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
377 : TYPE(scf_control_type), POINTER :: scf_control
378 : TYPE(section_vals_type), POINTER :: scf_section
379 : LOGICAL, INTENT(INOUT) :: diis_step
380 : TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
381 : POINTER :: probe
382 :
383 : INTEGER :: ispin, nspin
384 : REAL(KIND=dp) :: total_zeff_corr
385 :
386 76145 : nspin = SIZE(matrix_ks)
387 :
388 : CALL general_eigenproblem(scf_env, mos, matrix_ks, &
389 76145 : matrix_s, scf_control, scf_section, diis_step)
390 :
391 : total_zeff_corr = 0.0_dp
392 76145 : total_zeff_corr = scf_env%sum_zeff_corr
393 :
394 76145 : IF (ABS(total_zeff_corr) > 0.0_dp) THEN
395 : CALL set_mo_occupation(mo_array=mos, &
396 40 : smear=scf_control%smear, tot_zeff_corr=total_zeff_corr)
397 : ELSE
398 76105 : IF (PRESENT(probe) .EQV. .TRUE.) THEN
399 14 : scf_control%smear%do_smear = .FALSE.
400 : CALL set_mo_occupation(mo_array=mos, &
401 : smear=scf_control%smear, &
402 14 : probe=probe)
403 : ELSE
404 : CALL set_mo_occupation(mo_array=mos, &
405 76091 : smear=scf_control%smear)
406 : END IF
407 : END IF
408 :
409 161182 : DO ispin = 1, nspin
410 : CALL calculate_density_matrix(mos(ispin), &
411 161182 : scf_env%p_mix_new(ispin, 1)%matrix)
412 : END DO
413 :
414 76145 : END SUBROUTINE do_general_diag
415 :
416 : ! **************************************************************************************************
417 : !> \brief Kpoint diagonalization routine
418 : !> Transforms matrices to kpoint, distributes kpoint groups, performs
419 : !> general diagonalization (no storgae of overlap decomposition), stores
420 : !> MOs, calculates occupation numbers, calculates density matrices
421 : !> in kpoint representation, transforms density matrices to real space
422 : !> \param matrix_ks Kohn-sham matrices (RS indices, global)
423 : !> \param matrix_s Overlap matrices (RS indices, global)
424 : !> \param kpoints Kpoint environment
425 : !> \param scf_env SCF environment
426 : !> \param scf_control SCF control variables
427 : !> \param update_p ...
428 : !> \param diis_step ...
429 : !> \param diis_error ...
430 : !> \param qs_env ...
431 : !> \param probe ...
432 : !> \par History
433 : !> 08.2014 created [JGH]
434 : ! **************************************************************************************************
435 5676 : SUBROUTINE do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, update_p, &
436 : diis_step, diis_error, qs_env, probe)
437 :
438 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
439 : TYPE(kpoint_type), POINTER :: kpoints
440 : TYPE(qs_scf_env_type), POINTER :: scf_env
441 : TYPE(scf_control_type), POINTER :: scf_control
442 : LOGICAL, INTENT(IN) :: update_p
443 : LOGICAL, INTENT(INOUT) :: diis_step
444 : REAL(dp), INTENT(INOUT), OPTIONAL :: diis_error
445 : TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
446 : TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
447 : POINTER :: probe
448 :
449 : CHARACTER(len=*), PARAMETER :: routineN = 'do_general_diag_kp'
450 :
451 5676 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coeffs
452 : INTEGER :: handle, ib, igroup, ik, ikp, indx, &
453 : ispin, jb, kplocal, nb, nkp, &
454 : nkp_groups, nspin
455 : INTEGER, DIMENSION(2) :: kp_range
456 5676 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
457 5676 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
458 : LOGICAL :: do_diis, my_kpgrp, use_real_wfn
459 5676 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
460 5676 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
461 5676 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
462 : TYPE(cp_cfm_type) :: cksmat, cmos, csmat, cwork
463 5676 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools
464 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct, mo_struct
465 : TYPE(cp_fm_type) :: fmdummy, fmlocal, rksmat, rsmat
466 5676 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
467 : TYPE(cp_fm_type), POINTER :: imos, mo_coeff, rmos
468 : TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix, tmpmat
469 : TYPE(kpoint_env_type), POINTER :: kp
470 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_global
471 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
472 5676 : POINTER :: sab_nl
473 : TYPE(qs_matrix_pools_type), POINTER :: mpools
474 : TYPE(section_vals_type), POINTER :: scf_section
475 :
476 5676 : CALL timeset(routineN, handle)
477 :
478 5676 : NULLIFY (sab_nl)
479 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
480 : nkp_groups=nkp_groups, kp_dist=kp_dist, sab_nl=sab_nl, &
481 5676 : cell_to_index=cell_to_index)
482 5676 : CPASSERT(ASSOCIATED(sab_nl))
483 5676 : kplocal = kp_range(2) - kp_range(1) + 1
484 :
485 : !Whether we use DIIS for k-points
486 5676 : do_diis = .FALSE.
487 : IF (scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis .AND. .NOT. use_real_wfn &
488 5676 : .AND. PRESENT(diis_error) .AND. PRESENT(qs_env)) do_diis = .TRUE.
489 :
490 : ! allocate some work matrices
491 5676 : ALLOCATE (rmatrix, cmatrix, tmpmat)
492 : CALL dbcsr_create(rmatrix, template=matrix_ks(1, 1)%matrix, &
493 5676 : matrix_type=dbcsr_type_symmetric)
494 : CALL dbcsr_create(cmatrix, template=matrix_ks(1, 1)%matrix, &
495 5676 : matrix_type=dbcsr_type_antisymmetric)
496 : CALL dbcsr_create(tmpmat, template=matrix_ks(1, 1)%matrix, &
497 5676 : matrix_type=dbcsr_type_no_symmetry)
498 5676 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
499 5676 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
500 :
501 5676 : fmwork => scf_env%scf_work1
502 :
503 : ! fm pools to be used within a kpoint group
504 5676 : CALL get_kpoint_info(kpoints, mpools=mpools)
505 5676 : CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
506 :
507 5676 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
508 5676 : CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
509 :
510 5676 : IF (use_real_wfn) THEN
511 52 : CALL cp_fm_create(rksmat, matrix_struct)
512 52 : CALL cp_fm_create(rsmat, matrix_struct)
513 : ELSE
514 5624 : CALL cp_cfm_create(cksmat, matrix_struct)
515 5624 : CALL cp_cfm_create(csmat, matrix_struct)
516 5624 : CALL cp_cfm_create(cwork, matrix_struct)
517 5624 : kp => kpoints%kp_env(1)%kpoint_env
518 5624 : CALL get_mo_set(kp%mos(1, 1), mo_coeff=mo_coeff)
519 5624 : CALL cp_fm_get_info(mo_coeff, matrix_struct=mo_struct)
520 5624 : CALL cp_cfm_create(cmos, mo_struct)
521 : END IF
522 :
523 5676 : para_env => kpoints%blacs_env_all%para_env
524 5676 : nspin = SIZE(matrix_ks, 1)
525 185404 : ALLOCATE (info(kplocal*nspin*nkp_groups, 4))
526 :
527 : ! Setup and start all the communication
528 5676 : indx = 0
529 21136 : DO ikp = 1, kplocal
530 38488 : DO ispin = 1, nspin
531 57878 : DO igroup = 1, nkp_groups
532 : ! number of current kpoint
533 25066 : ik = kp_dist(1, igroup) + ikp - 1
534 25066 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
535 25066 : indx = indx + 1
536 25066 : IF (use_real_wfn) THEN
537 : ! FT of matrices KS and S, then transfer to FM type
538 62 : CALL dbcsr_set(rmatrix, 0.0_dp)
539 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_ks, ispin=ispin, &
540 62 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
541 62 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
542 62 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
543 : ! s matrix is not spin dependent
544 62 : CALL dbcsr_set(rmatrix, 0.0_dp)
545 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_s, ispin=1, &
546 62 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
547 62 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
548 62 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(3))
549 : ELSE
550 : ! FT of matrices KS and S, then transfer to FM type
551 25004 : CALL dbcsr_set(rmatrix, 0.0_dp)
552 25004 : CALL dbcsr_set(cmatrix, 0.0_dp)
553 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_ks, ispin=ispin, &
554 25004 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
555 25004 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
556 25004 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
557 25004 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
558 25004 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(2))
559 : ! s matrix is not spin dependent, double the work
560 25004 : CALL dbcsr_set(rmatrix, 0.0_dp)
561 25004 : CALL dbcsr_set(cmatrix, 0.0_dp)
562 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_s, ispin=1, &
563 25004 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
564 25004 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
565 25004 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(3))
566 25004 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
567 25004 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(4))
568 : END IF
569 : ! transfer to kpoint group
570 : ! redistribution of matrices, new blacs environment
571 : ! fmwork -> fmlocal -> rksmat/cksmat
572 : ! fmwork -> fmlocal -> rsmat/csmat
573 42418 : IF (use_real_wfn) THEN
574 62 : IF (my_kpgrp) THEN
575 62 : CALL cp_fm_start_copy_general(fmwork(1), rksmat, para_env, info(indx, 1))
576 62 : CALL cp_fm_start_copy_general(fmwork(3), rsmat, para_env, info(indx, 2))
577 : ELSE
578 0 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
579 0 : CALL cp_fm_start_copy_general(fmwork(3), fmdummy, para_env, info(indx, 2))
580 : END IF
581 : ELSE
582 25004 : IF (my_kpgrp) THEN
583 17290 : CALL cp_fm_start_copy_general(fmwork(1), fmlocal, para_env, info(indx, 1))
584 17290 : CALL cp_fm_start_copy_general(fmwork(2), fmlocal, para_env, info(indx, 2))
585 17290 : CALL cp_fm_start_copy_general(fmwork(3), fmlocal, para_env, info(indx, 3))
586 17290 : CALL cp_fm_start_copy_general(fmwork(4), fmlocal, para_env, info(indx, 4))
587 : ELSE
588 7714 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
589 7714 : CALL cp_fm_start_copy_general(fmwork(2), fmdummy, para_env, info(indx, 2))
590 7714 : CALL cp_fm_start_copy_general(fmwork(3), fmdummy, para_env, info(indx, 3))
591 7714 : CALL cp_fm_start_copy_general(fmwork(4), fmdummy, para_env, info(indx, 4))
592 : END IF
593 : END IF
594 : END DO
595 : END DO
596 : END DO
597 :
598 : ! Finish communication then diagonalise in each group
599 5676 : IF (do_diis) THEN
600 4068 : CALL get_qs_env(qs_env, para_env=para_env_global)
601 4068 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
602 4068 : CALL qs_diis_b_info_kp(kpoints%scf_diis_buffer, ib, nb)
603 4068 : indx = 0
604 13484 : DO ikp = 1, kplocal
605 24302 : DO ispin = 1, nspin
606 25894 : DO igroup = 1, nkp_groups
607 : ! number of current kpoint
608 15076 : ik = kp_dist(1, igroup) + ikp - 1
609 15076 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
610 4258 : indx = indx + 1
611 10818 : IF (my_kpgrp) THEN
612 10818 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
613 10818 : CALL cp_cfm_scale_and_add_fm(z_zero, cksmat, z_one, fmlocal)
614 10818 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
615 10818 : CALL cp_cfm_scale_and_add_fm(z_one, cksmat, gaussi, fmlocal)
616 10818 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 3))
617 10818 : CALL cp_cfm_scale_and_add_fm(z_zero, csmat, z_one, fmlocal)
618 10818 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 4))
619 10818 : CALL cp_cfm_scale_and_add_fm(z_one, csmat, gaussi, fmlocal)
620 : END IF
621 : END DO !igroup
622 :
623 10818 : kp => kpoints%kp_env(ikp)%kpoint_env
624 : CALL qs_diis_b_calc_err_kp(kpoints%scf_diis_buffer, ib, kp%mos, cksmat, csmat, &
625 20234 : ispin, ikp, kplocal, scf_section)
626 :
627 : END DO !ispin
628 : END DO !ikp
629 :
630 12204 : ALLOCATE (coeffs(nb))
631 : CALL qs_diis_b_step_kp(kpoints%scf_diis_buffer, coeffs, ib, nb, scf_env%iter_delta, diis_error, &
632 : diis_step, scf_control%eps_diis, nspin, nkp, kplocal, scf_control%nmixing, &
633 4068 : scf_section, para_env_global)
634 :
635 : !build the ks matrices and idagonalize
636 17552 : DO ikp = 1, kplocal
637 24302 : DO ispin = 1, nspin
638 10818 : kp => kpoints%kp_env(ikp)%kpoint_env
639 10818 : CALL cp_cfm_to_cfm(kpoints%scf_diis_buffer%smat(ikp), csmat)
640 :
641 10818 : CALL cp_cfm_set_all(cksmat, z_zero)
642 42200 : DO jb = 1, nb
643 42200 : CALL cp_cfm_scale_and_add(z_one, cksmat, coeffs(jb), kpoints%scf_diis_buffer%param(jb, ispin, ikp))
644 : END DO
645 :
646 10818 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
647 10818 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
648 10818 : IF (scf_env%cholesky_method == cholesky_off) THEN
649 : CALL cp_cfm_geeig_canon(cksmat, csmat, cmos, eigenvalues, cwork, &
650 16 : scf_control%eps_eigval)
651 : ELSE
652 10802 : CALL cp_cfm_geeig(cksmat, csmat, cmos, eigenvalues, cwork)
653 : END IF
654 : ! copy eigenvalues to imag set (keep them in sync)
655 242748 : kp%mos(2, ispin)%eigenvalues = eigenvalues
656 : ! split real and imaginary part of mos
657 20234 : CALL cp_cfm_to_fm(cmos, rmos, imos)
658 : END DO
659 : END DO
660 :
661 : ELSE !no DIIS
662 1608 : diis_step = .FALSE.
663 1608 : indx = 0
664 7652 : DO ikp = 1, kplocal
665 14186 : DO ispin = 1, nspin
666 16524 : DO igroup = 1, nkp_groups
667 : ! number of current kpoint
668 9990 : ik = kp_dist(1, igroup) + ikp - 1
669 9990 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
670 3456 : indx = indx + 1
671 6534 : IF (my_kpgrp) THEN
672 6534 : IF (use_real_wfn) THEN
673 62 : CALL cp_fm_finish_copy_general(rksmat, info(indx, 1))
674 62 : CALL cp_fm_finish_copy_general(rsmat, info(indx, 2))
675 : ELSE
676 6472 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
677 6472 : CALL cp_cfm_scale_and_add_fm(z_zero, cksmat, z_one, fmlocal)
678 6472 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
679 6472 : CALL cp_cfm_scale_and_add_fm(z_one, cksmat, gaussi, fmlocal)
680 6472 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 3))
681 6472 : CALL cp_cfm_scale_and_add_fm(z_zero, csmat, z_one, fmlocal)
682 6472 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 4))
683 6472 : CALL cp_cfm_scale_and_add_fm(z_one, csmat, gaussi, fmlocal)
684 : END IF
685 : END IF
686 : END DO
687 :
688 : ! Each kpoint group has now information on a kpoint to be diagonalized
689 : ! General eigensolver Hermite or Symmetric
690 6534 : kp => kpoints%kp_env(ikp)%kpoint_env
691 12578 : IF (use_real_wfn) THEN
692 62 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=mo_coeff, eigenvalues=eigenvalues)
693 62 : IF (scf_env%cholesky_method == cholesky_off) THEN
694 : CALL cp_fm_geeig_canon(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal, &
695 40 : scf_control%eps_eigval)
696 : ELSE
697 22 : CALL cp_fm_geeig(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal)
698 : END IF
699 : ELSE
700 6472 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
701 6472 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
702 6472 : IF (scf_env%cholesky_method == cholesky_off) THEN
703 : CALL cp_cfm_geeig_canon(cksmat, csmat, cmos, eigenvalues, cwork, &
704 242 : scf_control%eps_eigval)
705 : ELSE
706 6230 : CALL cp_cfm_geeig(cksmat, csmat, cmos, eigenvalues, cwork)
707 : END IF
708 : ! copy eigenvalues to imag set (keep them in sync)
709 432284 : kp%mos(2, ispin)%eigenvalues = eigenvalues
710 : ! split real and imaginary part of mos
711 6472 : CALL cp_cfm_to_fm(cmos, rmos, imos)
712 : END IF
713 : END DO
714 : END DO
715 : END IF
716 :
717 : ! Clean up communication
718 5676 : indx = 0
719 21136 : DO ikp = 1, kplocal
720 38488 : DO ispin = 1, nspin
721 57878 : DO igroup = 1, nkp_groups
722 : ! number of current kpoint
723 25066 : ik = kp_dist(1, igroup) + ikp - 1
724 25066 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
725 25066 : indx = indx + 1
726 42418 : IF (use_real_wfn) THEN
727 62 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
728 62 : CALL cp_fm_cleanup_copy_general(info(indx, 2))
729 : ELSE
730 25004 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
731 25004 : CALL cp_fm_cleanup_copy_general(info(indx, 2))
732 25004 : CALL cp_fm_cleanup_copy_general(info(indx, 3))
733 25004 : CALL cp_fm_cleanup_copy_general(info(indx, 4))
734 : END IF
735 : END DO
736 : END DO
737 : END DO
738 :
739 : ! All done
740 111616 : DEALLOCATE (info)
741 :
742 5676 : IF (update_p) THEN
743 : ! MO occupations
744 5666 : IF (PRESENT(probe) .EQV. .TRUE.) THEN
745 0 : scf_control%smear%do_smear = .FALSE.
746 : CALL kpoint_set_mo_occupation(kpoints, scf_control%smear, &
747 0 : probe=probe)
748 : ELSE
749 5666 : CALL kpoint_set_mo_occupation(kpoints, scf_control%smear)
750 : END IF
751 : ! density matrices
752 5666 : CALL kpoint_density_matrices(kpoints)
753 : ! density matrices in real space
754 : CALL kpoint_density_transform(kpoints, scf_env%p_mix_new, .FALSE., &
755 5666 : matrix_s(1, 1)%matrix, sab_nl, fmwork)
756 : END IF
757 :
758 5676 : CALL dbcsr_deallocate_matrix(rmatrix)
759 5676 : CALL dbcsr_deallocate_matrix(cmatrix)
760 5676 : CALL dbcsr_deallocate_matrix(tmpmat)
761 :
762 5676 : IF (use_real_wfn) THEN
763 52 : CALL cp_fm_release(rksmat)
764 52 : CALL cp_fm_release(rsmat)
765 : ELSE
766 5624 : CALL cp_cfm_release(cksmat)
767 5624 : CALL cp_cfm_release(csmat)
768 5624 : CALL cp_cfm_release(cwork)
769 5624 : CALL cp_cfm_release(cmos)
770 : END IF
771 5676 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
772 :
773 5676 : CALL timestop(handle)
774 :
775 17028 : END SUBROUTINE do_general_diag_kp
776 :
777 : ! **************************************************************************************************
778 : !> \brief inner loop within MOS subspace, to refine occupation and density,
779 : !> before next diagonalization of the Hamiltonian
780 : !> \param qs_env ...
781 : !> \param scf_env ...
782 : !> \param subspace_env ...
783 : !> \param mos ...
784 : !> \param rho ...
785 : !> \param ks_env ...
786 : !> \param scf_section ...
787 : !> \param scf_control ...
788 : !> \par History
789 : !> 09.2009 created [MI]
790 : !> \note it is assumed that when diagonalization is used, also some mixing procedure is active
791 : ! **************************************************************************************************
792 10 : SUBROUTINE do_scf_diag_subspace(qs_env, scf_env, subspace_env, mos, rho, &
793 : ks_env, scf_section, scf_control)
794 :
795 : TYPE(qs_environment_type), POINTER :: qs_env
796 : TYPE(qs_scf_env_type), POINTER :: scf_env
797 : TYPE(subspace_env_type), POINTER :: subspace_env
798 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
799 : TYPE(qs_rho_type), POINTER :: rho
800 : TYPE(qs_ks_env_type), POINTER :: ks_env
801 : TYPE(section_vals_type), POINTER :: scf_section
802 : TYPE(scf_control_type), POINTER :: scf_control
803 :
804 : CHARACTER(LEN=*), PARAMETER :: routineN = 'do_scf_diag_subspace'
805 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
806 :
807 : INTEGER :: handle, i, iloop, ispin, nao, nmo, &
808 : nspin, output_unit
809 : LOGICAL :: converged
810 : REAL(dp) :: ene_diff, ene_old, iter_delta, max_val, &
811 : sum_band, sum_val, t1, t2
812 10 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues, mo_occupations
813 10 : TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: eval_first, occ_first
814 : TYPE(cp_fm_type) :: work
815 : TYPE(cp_fm_type), POINTER :: c0, chc, evec, mo_coeff
816 : TYPE(cp_logger_type), POINTER :: logger
817 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
818 10 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
819 : TYPE(dft_control_type), POINTER :: dft_control
820 : TYPE(mp_para_env_type), POINTER :: para_env
821 : TYPE(qs_energy_type), POINTER :: energy
822 10 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
823 :
824 10 : CALL timeset(routineN, handle)
825 10 : NULLIFY (c0, chc, energy, evec, matrix_ks, mo_coeff, mo_eigenvalues, &
826 10 : mo_occupations, dft_control, rho_ao, rho_ao_kp)
827 :
828 10 : logger => cp_get_default_logger()
829 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIAG_SUB_SCF", &
830 10 : extension=".scfLog")
831 :
832 : !Extra loop keeping mos unchanged and refining the subspace occupation
833 10 : nspin = SIZE(mos)
834 10 : CALL qs_rho_get(rho, rho_ao=rho_ao, rho_ao_kp=rho_ao_kp)
835 :
836 40 : ALLOCATE (eval_first(nspin))
837 40 : ALLOCATE (occ_first(nspin))
838 20 : DO ispin = 1, nspin
839 : CALL get_mo_set(mo_set=mos(ispin), &
840 : nmo=nmo, &
841 : eigenvalues=mo_eigenvalues, &
842 10 : occupation_numbers=mo_occupations)
843 30 : ALLOCATE (eval_first(ispin)%array(nmo))
844 20 : ALLOCATE (occ_first(ispin)%array(nmo))
845 50 : eval_first(ispin)%array(1:nmo) = mo_eigenvalues(1:nmo)
846 70 : occ_first(ispin)%array(1:nmo) = mo_occupations(1:nmo)
847 : END DO
848 :
849 20 : DO ispin = 1, nspin
850 : ! does not yet handle k-points
851 10 : CALL dbcsr_copy(subspace_env%p_matrix_store(ispin)%matrix, rho_ao(ispin)%matrix)
852 20 : CALL dbcsr_copy(rho_ao(ispin)%matrix, scf_env%p_mix_new(ispin, 1)%matrix)
853 : END DO
854 :
855 10 : subspace_env%p_matrix_mix => scf_env%p_mix_new
856 :
857 10 : NULLIFY (matrix_ks, energy, para_env, matrix_s)
858 : CALL get_qs_env(qs_env, &
859 : matrix_ks=matrix_ks, &
860 : energy=energy, &
861 : matrix_s=matrix_s, &
862 : para_env=para_env, &
863 10 : dft_control=dft_control)
864 :
865 : ! mixing storage allocation
866 10 : IF (subspace_env%mixing_method >= gspace_mixing_nr) THEN
867 : CALL mixing_allocate(qs_env, subspace_env%mixing_method, scf_env%p_mix_new, &
868 0 : scf_env%p_delta, nspin, subspace_env%mixing_store)
869 0 : IF (dft_control%qs_control%gapw) THEN
870 0 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
871 : CALL mixing_init(subspace_env%mixing_method, rho, subspace_env%mixing_store, &
872 0 : para_env, rho_atom=rho_atom)
873 0 : ELSEIF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
874 0 : CALL charge_mixing_init(subspace_env%mixing_store)
875 0 : ELSEIF (dft_control%qs_control%semi_empirical) THEN
876 0 : CPABORT('SE Code not possible')
877 : ELSE
878 0 : CALL mixing_init(subspace_env%mixing_method, rho, subspace_env%mixing_store, para_env)
879 : END IF
880 : END IF
881 :
882 10 : ene_old = 0.0_dp
883 10 : ene_diff = 0.0_dp
884 10 : IF (output_unit > 0) THEN
885 0 : WRITE (output_unit, "(/T19,A)") '<<<<<<<<< SUBSPACE ROTATION <<<<<<<<<<'
886 : WRITE (output_unit, "(T4,A,T13,A,T21,A,T38,A,T51,A,T65,A/,T4,A)") &
887 0 : "In-step", "Time", "Convergence", "Band ene.", "Total ene.", "Energy diff.", REPEAT("-", 74)
888 : END IF
889 :
890 : ! recalculate density matrix here
891 :
892 : ! update of density
893 10 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
894 :
895 22 : DO iloop = 1, subspace_env%max_iter
896 20 : t1 = m_walltime()
897 20 : converged = .FALSE.
898 20 : ene_old = energy%total
899 :
900 20 : CALL qs_ks_did_change(ks_env, rho_changed=.TRUE.)
901 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., &
902 20 : just_energy=.FALSE., print_active=.FALSE.)
903 :
904 20 : max_val = 0.0_dp
905 20 : sum_val = 0.0_dp
906 20 : sum_band = 0.0_dp
907 40 : DO ispin = 1, SIZE(matrix_ks)
908 : CALL get_mo_set(mo_set=mos(ispin), &
909 : nao=nao, &
910 : nmo=nmo, &
911 : eigenvalues=mo_eigenvalues, &
912 : occupation_numbers=mo_occupations, &
913 20 : mo_coeff=mo_coeff)
914 :
915 : !compute C'HC
916 20 : chc => subspace_env%chc_mat(ispin)
917 20 : evec => subspace_env%c_vec(ispin)
918 20 : c0 => subspace_env%c0(ispin)
919 20 : CALL cp_fm_to_fm(mo_coeff, c0)
920 20 : CALL cp_fm_create(work, c0%matrix_struct)
921 20 : CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, c0, work, nmo)
922 20 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, rone, c0, work, rzero, chc)
923 20 : CALL cp_fm_release(work)
924 : !diagonalize C'HC
925 20 : CALL choose_eigv_solver(chc, evec, mo_eigenvalues)
926 :
927 : !rotate the mos by the eigenvectors of C'HC
928 20 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, rone, c0, evec, rzero, mo_coeff)
929 :
930 : CALL set_mo_occupation(mo_set=mos(ispin), &
931 20 : smear=scf_control%smear)
932 :
933 : ! does not yet handle k-points
934 : CALL calculate_density_matrix(mos(ispin), &
935 20 : subspace_env%p_matrix_mix(ispin, 1)%matrix)
936 :
937 160 : DO i = 1, nmo
938 100 : sum_band = sum_band + mo_eigenvalues(i)*mo_occupations(i)
939 : END DO
940 :
941 : !check for self consistency
942 : END DO
943 :
944 20 : IF (subspace_env%mixing_method == direct_mixing_nr) THEN
945 : CALL scf_env_density_mixing(subspace_env%p_matrix_mix, &
946 20 : scf_env%mixing_store, rho_ao_kp, para_env, iter_delta, iloop)
947 : ELSE
948 : CALL self_consistency_check(rho_ao_kp, scf_env%p_delta, para_env, &
949 0 : subspace_env%p_matrix_mix, delta=iter_delta)
950 : END IF
951 :
952 40 : DO ispin = 1, nspin
953 : ! does not yet handle k-points
954 40 : CALL dbcsr_copy(rho_ao(ispin)%matrix, subspace_env%p_matrix_mix(ispin, 1)%matrix)
955 : END DO
956 : ! update of density
957 20 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
958 : ! Mixing in reciprocal space
959 20 : IF (subspace_env%mixing_method >= gspace_mixing_nr) THEN
960 : CALL gspace_mixing(qs_env, scf_env%mixing_method, subspace_env%mixing_store, &
961 0 : rho, para_env, scf_env%iter_count)
962 : END IF
963 :
964 20 : ene_diff = energy%total - ene_old
965 : converged = (ABS(ene_diff) < subspace_env%eps_ene .AND. &
966 20 : iter_delta < subspace_env%eps_adapt*scf_env%iter_delta)
967 20 : t2 = m_walltime()
968 20 : IF (output_unit > 0) THEN
969 : WRITE (output_unit, "(T4,I5,T11,F8.3,T18,E14.4,T34,F12.5,T46,F16.8,T62,E14.4)") &
970 0 : iloop, t2 - t1, iter_delta, sum_band, energy%total, ene_diff
971 0 : CALL m_flush(output_unit)
972 : END IF
973 22 : IF (converged) THEN
974 8 : IF (output_unit > 0) WRITE (output_unit, "(T10,A,I6,A,/)") &
975 0 : " Reached convergence in ", iloop, " iterations "
976 : EXIT
977 : END IF
978 :
979 : END DO ! iloop
980 :
981 10 : NULLIFY (subspace_env%p_matrix_mix)
982 20 : DO ispin = 1, nspin
983 : ! does not yet handle k-points
984 10 : CALL dbcsr_copy(scf_env%p_mix_new(ispin, 1)%matrix, rho_ao(ispin)%matrix)
985 10 : CALL dbcsr_copy(rho_ao(ispin)%matrix, subspace_env%p_matrix_store(ispin)%matrix)
986 :
987 20 : DEALLOCATE (eval_first(ispin)%array, occ_first(ispin)%array)
988 : END DO
989 10 : DEALLOCATE (eval_first, occ_first)
990 :
991 10 : CALL timestop(handle)
992 :
993 10 : END SUBROUTINE do_scf_diag_subspace
994 :
995 : ! **************************************************************************************************
996 : !> \brief ...
997 : !> \param subspace_env ...
998 : !> \param qs_env ...
999 : !> \param mos ...
1000 : ! **************************************************************************************************
1001 2 : SUBROUTINE diag_subspace_allocate(subspace_env, qs_env, mos)
1002 :
1003 : TYPE(subspace_env_type), POINTER :: subspace_env
1004 : TYPE(qs_environment_type), POINTER :: qs_env
1005 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1006 :
1007 : CHARACTER(LEN=*), PARAMETER :: routineN = 'diag_subspace_allocate'
1008 :
1009 : INTEGER :: handle, i, ispin, nmo, nspin
1010 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
1011 : TYPE(cp_fm_type), POINTER :: mo_coeff
1012 2 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1013 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1014 2 : POINTER :: sab_orb
1015 :
1016 2 : CALL timeset(routineN, handle)
1017 :
1018 2 : NULLIFY (sab_orb, matrix_s)
1019 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb, &
1020 2 : matrix_s=matrix_s)
1021 :
1022 2 : nspin = SIZE(mos)
1023 : ! *** allocate p_atrix_store ***
1024 2 : IF (.NOT. ASSOCIATED(subspace_env%p_matrix_store)) THEN
1025 2 : CALL dbcsr_allocate_matrix_set(subspace_env%p_matrix_store, nspin)
1026 :
1027 4 : DO i = 1, nspin
1028 2 : ALLOCATE (subspace_env%p_matrix_store(i)%matrix)
1029 : CALL dbcsr_create(matrix=subspace_env%p_matrix_store(i)%matrix, template=matrix_s(1)%matrix, &
1030 2 : name="DENSITY_STORE", matrix_type=dbcsr_type_symmetric)
1031 : CALL cp_dbcsr_alloc_block_from_nbl(subspace_env%p_matrix_store(i)%matrix, &
1032 2 : sab_orb)
1033 4 : CALL dbcsr_set(subspace_env%p_matrix_store(i)%matrix, 0.0_dp)
1034 : END DO
1035 :
1036 : END IF
1037 :
1038 8 : ALLOCATE (subspace_env%chc_mat(nspin))
1039 8 : ALLOCATE (subspace_env%c_vec(nspin))
1040 8 : ALLOCATE (subspace_env%c0(nspin))
1041 :
1042 4 : DO ispin = 1, nspin
1043 2 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1044 2 : CALL cp_fm_create(subspace_env%c0(ispin), mo_coeff%matrix_struct)
1045 2 : NULLIFY (fm_struct_tmp)
1046 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo, ncol_global=nmo, &
1047 : para_env=mo_coeff%matrix_struct%para_env, &
1048 2 : context=mo_coeff%matrix_struct%context)
1049 2 : CALL cp_fm_create(subspace_env%chc_mat(ispin), fm_struct_tmp, "chc")
1050 2 : CALL cp_fm_create(subspace_env%c_vec(ispin), fm_struct_tmp, "vec")
1051 6 : CALL cp_fm_struct_release(fm_struct_tmp)
1052 : END DO
1053 :
1054 2 : CALL timestop(handle)
1055 :
1056 2 : END SUBROUTINE diag_subspace_allocate
1057 :
1058 : ! **************************************************************************************************
1059 : !> \brief the inner loop of scf, specific to diagonalization without S matrix
1060 : !> basically, in goes the ks matrix out goes a new p matrix
1061 : !> \param scf_env ...
1062 : !> \param mos ...
1063 : !> \param matrix_ks ...
1064 : !> \param scf_control ...
1065 : !> \param scf_section ...
1066 : !> \param diis_step ...
1067 : !> \par History
1068 : !> 03.2006 created [Joost VandeVondele]
1069 : ! **************************************************************************************************
1070 17898 : SUBROUTINE do_special_diag(scf_env, mos, matrix_ks, scf_control, &
1071 : scf_section, diis_step)
1072 :
1073 : TYPE(qs_scf_env_type), POINTER :: scf_env
1074 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
1075 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1076 : TYPE(scf_control_type), POINTER :: scf_control
1077 : TYPE(section_vals_type), POINTER :: scf_section
1078 : LOGICAL, INTENT(INOUT) :: diis_step
1079 :
1080 : INTEGER :: ispin, nspin
1081 : LOGICAL :: do_level_shift, use_jacobi
1082 : REAL(KIND=dp) :: diis_error
1083 :
1084 17898 : nspin = SIZE(matrix_ks)
1085 :
1086 36570 : DO ispin = 1, nspin
1087 36570 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, scf_env%scf_work1(ispin))
1088 : END DO
1089 17898 : IF (scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis) THEN
1090 : CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
1091 : scf_env%scf_work2, scf_env%iter_delta, diis_error, diis_step, &
1092 : scf_control%eps_diis, scf_control%nmixing, &
1093 15304 : scf_section=scf_section)
1094 : ELSE
1095 2594 : diis_step = .FALSE.
1096 : END IF
1097 :
1098 17898 : IF ((scf_env%iter_count > 1) .AND. (scf_env%iter_delta < scf_control%diagonalization%eps_jacobi)) THEN
1099 18 : use_jacobi = .TRUE.
1100 : ELSE
1101 17880 : use_jacobi = .FALSE.
1102 : END IF
1103 :
1104 : do_level_shift = ((scf_control%level_shift /= 0.0_dp) .AND. &
1105 17898 : ((scf_control%density_guess == core_guess) .OR. (scf_env%iter_count > 1)))
1106 17898 : IF (diis_step) THEN
1107 11880 : scf_env%iter_param = diis_error
1108 11880 : IF (use_jacobi) THEN
1109 18 : scf_env%iter_method = "DIIS/Jacobi"
1110 : ELSE
1111 11862 : scf_env%iter_method = "DIIS/Diag."
1112 : END IF
1113 : ELSE
1114 6018 : IF (scf_env%mixing_method == 1) THEN
1115 6018 : scf_env%iter_param = scf_env%p_mix_alpha
1116 6018 : IF (use_jacobi) THEN
1117 0 : scf_env%iter_method = "P_Mix/Jacobi"
1118 : ELSE
1119 6018 : scf_env%iter_method = "P_Mix/Diag."
1120 : END IF
1121 0 : ELSEIF (scf_env%mixing_method > 1) THEN
1122 0 : scf_env%iter_param = scf_env%mixing_store%alpha
1123 0 : IF (use_jacobi) THEN
1124 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Jacobi"
1125 : ELSE
1126 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Diag."
1127 : END IF
1128 : END IF
1129 : END IF
1130 17898 : scf_env%iter_delta = 0.0_dp
1131 :
1132 36570 : DO ispin = 1, nspin
1133 : CALL eigensolver_simple(matrix_ks=scf_env%scf_work1(ispin), &
1134 : mo_set=mos(ispin), &
1135 : work=scf_env%scf_work2, &
1136 : do_level_shift=do_level_shift, &
1137 : level_shift=scf_control%level_shift, &
1138 : use_jacobi=use_jacobi, &
1139 36570 : jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
1140 : END DO
1141 :
1142 : CALL set_mo_occupation(mo_array=mos, &
1143 17898 : smear=scf_control%smear)
1144 :
1145 36570 : DO ispin = 1, nspin
1146 : ! does not yet handle k-points
1147 : CALL calculate_density_matrix(mos(ispin), &
1148 36570 : scf_env%p_mix_new(ispin, 1)%matrix)
1149 : END DO
1150 :
1151 17898 : END SUBROUTINE do_special_diag
1152 :
1153 : ! **************************************************************************************************
1154 : !> \brief the inner loop of scf, specific to iterative diagonalization using OT
1155 : !> with S matrix; basically, in goes the ks matrix out goes a new p matrix
1156 : !> \param scf_env ...
1157 : !> \param mos ...
1158 : !> \param matrix_ks ...
1159 : !> \param matrix_s ...
1160 : !> \param scf_control ...
1161 : !> \param scf_section ...
1162 : !> \param diis_step ...
1163 : !> \par History
1164 : !> 10.2008 created [JGH]
1165 : ! **************************************************************************************************
1166 64 : SUBROUTINE do_ot_diag(scf_env, mos, matrix_ks, matrix_s, &
1167 : scf_control, scf_section, diis_step)
1168 :
1169 : TYPE(qs_scf_env_type), POINTER :: scf_env
1170 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
1171 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
1172 : TYPE(scf_control_type), POINTER :: scf_control
1173 : TYPE(section_vals_type), POINTER :: scf_section
1174 : LOGICAL, INTENT(INOUT) :: diis_step
1175 :
1176 : INTEGER :: homo, ispin, nmo, nspin
1177 : REAL(KIND=dp) :: diis_error, eps_iter
1178 64 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1179 : TYPE(cp_fm_type), POINTER :: mo_coeff
1180 :
1181 64 : NULLIFY (eigenvalues)
1182 :
1183 64 : nspin = SIZE(matrix_ks)
1184 :
1185 172 : DO ispin = 1, nspin
1186 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, &
1187 172 : scf_env%scf_work1(ispin))
1188 : END DO
1189 :
1190 64 : IF ((scf_env%iter_count > 1) .AND. (.NOT. scf_env%skip_diis)) THEN
1191 : CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
1192 : scf_env%scf_work2, scf_env%iter_delta, diis_error, diis_step, &
1193 : scf_control%eps_diis, scf_control%nmixing, &
1194 : s_matrix=matrix_s, &
1195 48 : scf_section=scf_section)
1196 : ELSE
1197 16 : diis_step = .FALSE.
1198 : END IF
1199 :
1200 64 : eps_iter = scf_control%diagonalization%eps_iter
1201 64 : IF (diis_step) THEN
1202 20 : scf_env%iter_param = diis_error
1203 20 : scf_env%iter_method = "DIIS/OTdiag"
1204 54 : DO ispin = 1, nspin
1205 : CALL copy_fm_to_dbcsr(scf_env%scf_work1(ispin), &
1206 54 : matrix_ks(ispin)%matrix, keep_sparsity=.TRUE.)
1207 : END DO
1208 20 : eps_iter = MAX(eps_iter, scf_control%diagonalization%eps_adapt*diis_error)
1209 : ELSE
1210 44 : IF (scf_env%mixing_method == 1) THEN
1211 44 : scf_env%iter_param = scf_env%p_mix_alpha
1212 44 : scf_env%iter_method = "P_Mix/OTdiag."
1213 0 : ELSEIF (scf_env%mixing_method > 1) THEN
1214 0 : scf_env%iter_param = scf_env%mixing_store%alpha
1215 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/OTdiag."
1216 : END IF
1217 : END IF
1218 :
1219 64 : scf_env%iter_delta = 0.0_dp
1220 :
1221 172 : DO ispin = 1, nspin
1222 : CALL get_mo_set(mos(ispin), &
1223 : mo_coeff=mo_coeff, &
1224 : eigenvalues=eigenvalues, &
1225 : nmo=nmo, &
1226 108 : homo=homo)
1227 : CALL ot_eigensolver(matrix_h=matrix_ks(ispin)%matrix, &
1228 : matrix_s=matrix_s(1)%matrix, &
1229 : matrix_c_fm=mo_coeff, &
1230 : preconditioner=scf_env%ot_preconditioner(1)%preconditioner, &
1231 : eps_gradient=eps_iter, &
1232 : iter_max=scf_control%diagonalization%max_iter, &
1233 : silent=.TRUE., &
1234 108 : ot_settings=scf_control%diagonalization%ot_settings)
1235 : CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin)%matrix, &
1236 : evals_arg=eigenvalues, &
1237 108 : do_rotation=.TRUE.)
1238 : CALL copy_fm_to_dbcsr(mos(ispin)%mo_coeff, &
1239 280 : mos(ispin)%mo_coeff_b)
1240 : !fm->dbcsr
1241 : END DO
1242 :
1243 : CALL set_mo_occupation(mo_array=mos, &
1244 64 : smear=scf_control%smear)
1245 :
1246 172 : DO ispin = 1, nspin
1247 : ! does not yet handle k-points
1248 : CALL calculate_density_matrix(mos(ispin), &
1249 172 : scf_env%p_mix_new(ispin, 1)%matrix)
1250 : END DO
1251 :
1252 64 : END SUBROUTINE do_ot_diag
1253 :
1254 : ! **************************************************************************************************
1255 : !> \brief Solve a set restricted open Kohn-Sham (ROKS) equations based on the
1256 : !> alpha and beta Kohn-Sham matrices from unrestricted Kohn-Sham.
1257 : !> \param scf_env ...
1258 : !> \param mos ...
1259 : !> \param matrix_ks ...
1260 : !> \param matrix_s ...
1261 : !> \param scf_control ...
1262 : !> \param scf_section ...
1263 : !> \param diis_step ...
1264 : !> \param orthogonal_basis ...
1265 : !> \par History
1266 : !> 04.2006 created [MK]
1267 : !> Revised (01.05.06,MK)
1268 : !> \note
1269 : !> this is only a high-spin ROKS.
1270 : ! **************************************************************************************************
1271 1104 : SUBROUTINE do_roks_diag(scf_env, mos, matrix_ks, matrix_s, &
1272 : scf_control, scf_section, diis_step, &
1273 : orthogonal_basis)
1274 :
1275 : ! Literature: - C. C. J. Roothaan, Rev. Mod. Phys. 32, 179 (1960)
1276 : ! - M. F. Guest and V. R. Saunders, Mol. Phys. 28(3), 819 (1974)
1277 : ! - M. Filatov and S. Shaik, Chem. Phys. Lett. 288, 689 (1998)
1278 :
1279 : TYPE(qs_scf_env_type), POINTER :: scf_env
1280 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1281 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
1282 : TYPE(scf_control_type), POINTER :: scf_control
1283 : TYPE(section_vals_type), POINTER :: scf_section
1284 : LOGICAL, INTENT(INOUT) :: diis_step
1285 : LOGICAL, INTENT(IN) :: orthogonal_basis
1286 :
1287 : CHARACTER(LEN=*), PARAMETER :: routineN = 'do_roks_diag'
1288 :
1289 : INTEGER :: handle, homoa, homob, imo, nalpha, nao, &
1290 : nbeta, nmo
1291 : REAL(KIND=dp) :: diis_error, level_shift_loc
1292 1104 : REAL(KIND=dp), DIMENSION(:), POINTER :: eiga, eigb, occa, occb
1293 : TYPE(cp_fm_type), POINTER :: ksa, ksb, mo2ao, moa, mob, ortho, work
1294 :
1295 : ! -------------------------------------------------------------------------
1296 :
1297 1104 : CALL timeset(routineN, handle)
1298 :
1299 1104 : IF (scf_env%cholesky_method == cholesky_inverse) THEN
1300 0 : ortho => scf_env%ortho_m1
1301 : ELSE
1302 1104 : ortho => scf_env%ortho
1303 : END IF
1304 1104 : work => scf_env%scf_work2
1305 :
1306 1104 : ksa => scf_env%scf_work1(1)
1307 1104 : ksb => scf_env%scf_work1(2)
1308 :
1309 1104 : CALL copy_dbcsr_to_fm(matrix_ks(1)%matrix, ksa)
1310 1104 : CALL copy_dbcsr_to_fm(matrix_ks(2)%matrix, ksb)
1311 :
1312 : ! Get MO information
1313 :
1314 : CALL get_mo_set(mo_set=mos(1), &
1315 : nao=nao, &
1316 : nmo=nmo, &
1317 : nelectron=nalpha, &
1318 : homo=homoa, &
1319 : eigenvalues=eiga, &
1320 : occupation_numbers=occa, &
1321 1104 : mo_coeff=moa)
1322 :
1323 : CALL get_mo_set(mo_set=mos(2), &
1324 : nelectron=nbeta, &
1325 : homo=homob, &
1326 : eigenvalues=eigb, &
1327 : occupation_numbers=occb, &
1328 1104 : mo_coeff=mob)
1329 :
1330 : ! Define the amount of level-shifting
1331 :
1332 1104 : IF ((scf_control%level_shift /= 0.0_dp) .AND. &
1333 : ((scf_control%density_guess == core_guess) .OR. &
1334 : (scf_control%density_guess == restart_guess) .OR. &
1335 : (scf_env%iter_count > 1))) THEN
1336 20 : level_shift_loc = scf_control%level_shift
1337 : ELSE
1338 1084 : level_shift_loc = 0.0_dp
1339 : END IF
1340 :
1341 : IF ((scf_env%iter_count > 1) .OR. &
1342 1104 : (scf_control%density_guess == core_guess) .OR. &
1343 : (scf_control%density_guess == restart_guess)) THEN
1344 :
1345 : ! Transform the spin unrestricted alpha and beta Kohn-Sham matrices
1346 : ! from AO basis to MO basis: K(MO) = C(T)*K(AO)*C
1347 :
1348 998 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, moa, 0.0_dp, work)
1349 998 : CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, moa, work, 0.0_dp, ksa)
1350 :
1351 998 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksb, moa, 0.0_dp, work)
1352 998 : CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, moa, work, 0.0_dp, ksb)
1353 :
1354 : ! Combine the spin unrestricted alpha and beta Kohn-Sham matrices
1355 : ! in the MO basis
1356 :
1357 998 : IF (scf_control%roks_scheme == general_roks) THEN
1358 : CALL combine_ks_matrices(ksa, ksb, occa, occb, scf_control%roks_f, &
1359 0 : nalpha, nbeta)
1360 998 : ELSE IF (scf_control%roks_scheme == high_spin_roks) THEN
1361 998 : CALL combine_ks_matrices(ksa, ksb, occa, occb, scf_control%roks_parameter)
1362 : ELSE
1363 0 : CPABORT("Unknown ROKS scheme requested")
1364 : END IF
1365 :
1366 : ! Back-transform the restricted open Kohn-Sham matrix from MO basis
1367 : ! to AO basis
1368 :
1369 998 : IF (orthogonal_basis) THEN
1370 : ! Q = C
1371 454 : mo2ao => moa
1372 : ELSE
1373 : ! Q = S*C
1374 544 : mo2ao => mob
1375 : !MK CALL copy_sm_to_fm(matrix_s(1)%matrix,work)
1376 : !MK CALL cp_fm_symm("L", "U",nao, nao, 1.0_dp, work, moa, 0.0_dp, mo2ao)
1377 544 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, moa, mo2ao, nao)
1378 : END IF
1379 :
1380 : ! K(AO) = Q*K(MO)*Q(T)
1381 :
1382 998 : CALL parallel_gemm("N", "T", nao, nao, nao, 1.0_dp, ksa, mo2ao, 0.0_dp, work)
1383 998 : CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, mo2ao, work, 0.0_dp, ksa)
1384 :
1385 : ELSE
1386 :
1387 : ! No transformation matrix available, yet. The closed shell part,
1388 : ! i.e. the beta Kohn-Sham matrix in AO basis, is taken.
1389 : ! There might be better choices, anyhow.
1390 :
1391 106 : CALL cp_fm_to_fm(ksb, ksa)
1392 :
1393 : END IF
1394 :
1395 : ! Update DIIS buffer and possibly perform DIIS extrapolation step
1396 :
1397 1104 : IF (scf_env%iter_count > 1) THEN
1398 992 : IF (orthogonal_basis) THEN
1399 : CALL qs_diis_b_step(diis_buffer=scf_env%scf_diis_buffer, &
1400 : mo_array=mos, &
1401 : kc=scf_env%scf_work1, &
1402 : sc=work, &
1403 : delta=scf_env%iter_delta, &
1404 : error_max=diis_error, &
1405 : diis_step=diis_step, &
1406 : eps_diis=scf_control%eps_diis, &
1407 : scf_section=scf_section, &
1408 450 : roks=.TRUE.)
1409 450 : CPASSERT(scf_env%iter_delta == scf_env%iter_delta)
1410 : ELSE
1411 : CALL qs_diis_b_step(diis_buffer=scf_env%scf_diis_buffer, &
1412 : mo_array=mos, &
1413 : kc=scf_env%scf_work1, &
1414 : sc=work, &
1415 : delta=scf_env%iter_delta, &
1416 : error_max=diis_error, &
1417 : diis_step=diis_step, &
1418 : eps_diis=scf_control%eps_diis, &
1419 : scf_section=scf_section, &
1420 : s_matrix=matrix_s, &
1421 542 : roks=.TRUE.)
1422 : END IF
1423 : END IF
1424 :
1425 1104 : IF (diis_step) THEN
1426 692 : scf_env%iter_param = diis_error
1427 692 : scf_env%iter_method = "DIIS/Diag."
1428 : ELSE
1429 412 : IF (scf_env%mixing_method == 1) THEN
1430 412 : scf_env%iter_param = scf_env%p_mix_alpha
1431 412 : scf_env%iter_method = "P_Mix/Diag."
1432 0 : ELSEIF (scf_env%mixing_method > 1) THEN
1433 0 : scf_env%iter_param = scf_env%mixing_store%alpha
1434 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Diag."
1435 : END IF
1436 : END IF
1437 :
1438 1104 : scf_env%iter_delta = 0.0_dp
1439 :
1440 1104 : IF (level_shift_loc /= 0.0_dp) THEN
1441 :
1442 : ! Transform the current Kohn-Sham matrix from AO to MO basis
1443 : ! for level-shifting using the current MO set
1444 :
1445 20 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, moa, 0.0_dp, work)
1446 20 : CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, moa, work, 0.0_dp, ksa)
1447 :
1448 : ! Apply level-shifting using 50:50 split of the shift (could be relaxed)
1449 :
1450 60 : DO imo = homob + 1, homoa
1451 60 : CALL cp_fm_add_to_element(ksa, imo, imo, 0.5_dp*level_shift_loc)
1452 : END DO
1453 220 : DO imo = homoa + 1, nmo
1454 220 : CALL cp_fm_add_to_element(ksa, imo, imo, level_shift_loc)
1455 : END DO
1456 :
1457 1084 : ELSE IF (.NOT. orthogonal_basis) THEN
1458 :
1459 : ! Transform the current Kohn-Sham matrix to an orthogonal basis
1460 572 : SELECT CASE (scf_env%cholesky_method)
1461 : CASE (cholesky_reduce)
1462 0 : CALL cp_fm_cholesky_reduce(ksa, ortho)
1463 : CASE (cholesky_restore)
1464 508 : CALL cp_fm_uplo_to_full(ksa, work)
1465 : CALL cp_fm_cholesky_restore(ksa, nao, ortho, work, &
1466 508 : "SOLVE", pos="RIGHT")
1467 : CALL cp_fm_cholesky_restore(work, nao, ortho, ksa, &
1468 508 : "SOLVE", pos="LEFT", transa="T")
1469 : CASE (cholesky_inverse)
1470 0 : CALL cp_fm_uplo_to_full(ksa, work)
1471 : CALL cp_fm_cholesky_restore(ksa, nao, ortho, work, &
1472 0 : "MULTIPLY", pos="RIGHT")
1473 : CALL cp_fm_cholesky_restore(work, nao, ortho, ksa, &
1474 0 : "MULTIPLY", pos="LEFT", transa="T")
1475 : CASE (cholesky_off)
1476 64 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, ortho, 0.0_dp, work)
1477 636 : CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, ortho, work, 0.0_dp, ksa)
1478 : END SELECT
1479 :
1480 : END IF
1481 :
1482 : ! Diagonalization of the ROKS operator matrix
1483 :
1484 1104 : CALL choose_eigv_solver(ksa, work, eiga)
1485 :
1486 : ! Back-transformation of the orthonormal eigenvectors if needed
1487 :
1488 1104 : IF (level_shift_loc /= 0.0_dp) THEN
1489 : ! Use old MO set for back-transformation if level-shifting was applied
1490 20 : CALL cp_fm_to_fm(moa, ortho)
1491 20 : CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, ortho, work, 0.0_dp, moa)
1492 : ELSE
1493 1084 : IF (orthogonal_basis) THEN
1494 512 : CALL cp_fm_to_fm(work, moa)
1495 : ELSE
1496 1080 : SELECT CASE (scf_env%cholesky_method)
1497 : CASE (cholesky_reduce, cholesky_restore)
1498 508 : CALL cp_fm_cholesky_restore(work, nmo, ortho, moa, "SOLVE")
1499 : CASE (cholesky_inverse)
1500 0 : CALL cp_fm_cholesky_restore(work, nmo, ortho, moa, "MULTIPLY")
1501 : CASE (cholesky_off)
1502 572 : CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, ortho, work, 0.0_dp, moa)
1503 : END SELECT
1504 : END IF
1505 : END IF
1506 :
1507 : ! Correct MO eigenvalues, if level-shifting was applied
1508 :
1509 1104 : IF (level_shift_loc /= 0.0_dp) THEN
1510 60 : DO imo = homob + 1, homoa
1511 60 : eiga(imo) = eiga(imo) - 0.5_dp*level_shift_loc
1512 : END DO
1513 220 : DO imo = homoa + 1, nmo
1514 220 : eiga(imo) = eiga(imo) - level_shift_loc
1515 : END DO
1516 : END IF
1517 :
1518 : ! Update also the beta MO set
1519 :
1520 34556 : eigb(:) = eiga(:)
1521 1104 : CALL cp_fm_to_fm(moa, mob)
1522 :
1523 : ! Calculate the new alpha and beta density matrix
1524 :
1525 : ! does not yet handle k-points
1526 1104 : CALL calculate_density_matrix(mos(1), scf_env%p_mix_new(1, 1)%matrix)
1527 1104 : CALL calculate_density_matrix(mos(2), scf_env%p_mix_new(2, 1)%matrix)
1528 :
1529 1104 : CALL timestop(handle)
1530 :
1531 1104 : END SUBROUTINE do_roks_diag
1532 :
1533 : ! **************************************************************************************************
1534 : !> \brief iterative diagonalization using the block Krylov-space approach
1535 : !> \param scf_env ...
1536 : !> \param mos ...
1537 : !> \param matrix_ks ...
1538 : !> \param scf_control ...
1539 : !> \param scf_section ...
1540 : !> \param check_moconv_only ...
1541 : !> \param
1542 : !> \par History
1543 : !> 05.2009 created [MI]
1544 : ! **************************************************************************************************
1545 :
1546 38 : SUBROUTINE do_block_krylov_diag(scf_env, mos, matrix_ks, &
1547 : scf_control, scf_section, check_moconv_only)
1548 :
1549 : TYPE(qs_scf_env_type), POINTER :: scf_env
1550 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
1551 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1552 : TYPE(scf_control_type), POINTER :: scf_control
1553 : TYPE(section_vals_type), POINTER :: scf_section
1554 : LOGICAL, INTENT(IN), OPTIONAL :: check_moconv_only
1555 :
1556 : CHARACTER(LEN=*), PARAMETER :: routineN = 'do_block_krylov_diag'
1557 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1558 :
1559 : INTEGER :: handle, homo, ispin, iter, nao, nmo, &
1560 : output_unit
1561 : LOGICAL :: converged, my_check_moconv_only
1562 : REAL(dp) :: eps_iter, t1, t2
1563 38 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
1564 : TYPE(cp_fm_type), POINTER :: c0, c1, chc, evec, ks, mo_coeff, ortho, &
1565 : work
1566 : TYPE(cp_logger_type), POINTER :: logger
1567 :
1568 76 : logger => cp_get_default_logger()
1569 38 : CALL timeset(routineN, handle)
1570 :
1571 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%LANCZOS", &
1572 38 : extension=".scfLog")
1573 :
1574 38 : my_check_moconv_only = .FALSE.
1575 38 : IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
1576 :
1577 38 : NULLIFY (mo_coeff, ortho, work, ks)
1578 38 : NULLIFY (mo_eigenvalues)
1579 38 : NULLIFY (c0, c1)
1580 :
1581 38 : IF (scf_env%cholesky_method == cholesky_inverse) THEN
1582 38 : ortho => scf_env%ortho_m1
1583 : ELSE
1584 0 : ortho => scf_env%ortho
1585 : END IF
1586 38 : work => scf_env%scf_work2
1587 :
1588 76 : DO ispin = 1, SIZE(matrix_ks)
1589 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, &
1590 76 : scf_env%scf_work1(ispin))
1591 : END DO
1592 :
1593 38 : IF (scf_env%mixing_method == 1) THEN
1594 0 : scf_env%iter_param = scf_env%p_mix_alpha
1595 0 : scf_env%iter_method = "P_Mix/Lanczos"
1596 : ELSE
1597 : ! scf_env%iter_param = scf_env%mixing_store%alpha
1598 38 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Lanc."
1599 : END IF
1600 :
1601 76 : DO ispin = 1, SIZE(matrix_ks)
1602 :
1603 38 : ks => scf_env%scf_work1(ispin)
1604 38 : CALL cp_fm_uplo_to_full(ks, work)
1605 :
1606 : CALL get_mo_set(mo_set=mos(ispin), &
1607 : nao=nao, &
1608 : nmo=nmo, &
1609 : homo=homo, &
1610 : eigenvalues=mo_eigenvalues, &
1611 38 : mo_coeff=mo_coeff)
1612 :
1613 38 : NULLIFY (c0, c1)
1614 38 : c0 => scf_env%krylov_space%mo_conv(ispin)
1615 38 : c1 => scf_env%krylov_space%mo_refine(ispin)
1616 38 : SELECT CASE (scf_env%cholesky_method)
1617 : CASE (cholesky_reduce)
1618 0 : CALL cp_fm_cholesky_reduce(ks, ortho)
1619 0 : CALL cp_fm_uplo_to_full(ks, work)
1620 0 : CALL cp_fm_cholesky_restore(mo_coeff, nmo, ortho, c0, "MULTIPLY")
1621 : CASE (cholesky_restore)
1622 : CALL cp_fm_cholesky_restore(ks, nao, ortho, work, &
1623 0 : "SOLVE", pos="RIGHT")
1624 : CALL cp_fm_cholesky_restore(work, nao, ortho, ks, &
1625 0 : "SOLVE", pos="LEFT", transa="T")
1626 0 : CALL cp_fm_cholesky_restore(mo_coeff, nmo, ortho, c0, "MULTIPLY")
1627 : CASE (cholesky_inverse)
1628 : CALL cp_fm_cholesky_restore(ks, nao, ortho, work, &
1629 38 : "MULTIPLY", pos="RIGHT")
1630 : CALL cp_fm_cholesky_restore(work, nao, ortho, ks, &
1631 38 : "MULTIPLY", pos="LEFT", transa="T")
1632 76 : CALL cp_fm_cholesky_restore(mo_coeff, nmo, ortho, c0, "SOLVE")
1633 : END SELECT
1634 :
1635 38 : scf_env%krylov_space%nmo_nc = nmo
1636 38 : scf_env%krylov_space%nmo_conv = 0
1637 :
1638 38 : t1 = m_walltime()
1639 38 : IF (output_unit > 0) THEN
1640 0 : WRITE (output_unit, "(/T15,A)") '<<<<<<<<< LANCZOS REFINEMENT <<<<<<<<<<'
1641 : WRITE (output_unit, "(T8,A,T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") &
1642 0 : " Spin ", " Cycle ", &
1643 0 : " conv. MOS ", " B2MAX ", " B2MIN ", " Time", REPEAT("-", 60)
1644 : END IF
1645 38 : eps_iter = MAX(scf_env%krylov_space%eps_conv, scf_env%krylov_space%eps_adapt*scf_env%iter_delta)
1646 38 : iter = 0
1647 38 : converged = .FALSE.
1648 : !Check convergence of MOS
1649 114 : IF (my_check_moconv_only) THEN
1650 :
1651 : CALL lanczos_refinement(scf_env%krylov_space, ks, c0, c1, mo_eigenvalues, &
1652 0 : nao, eps_iter, ispin, check_moconv_only=my_check_moconv_only)
1653 0 : t2 = m_walltime()
1654 0 : IF (output_unit > 0) &
1655 : WRITE (output_unit, '(T8,I3,T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
1656 0 : ispin, iter, scf_env%krylov_space%nmo_conv, &
1657 0 : scf_env%krylov_space%max_res_norm, scf_env%krylov_space%min_res_norm, t2 - t1
1658 :
1659 : CYCLE
1660 : ELSE
1661 : !Block Lanczos refinement
1662 620 : DO iter = 1, scf_env%krylov_space%max_iter
1663 : CALL lanczos_refinement_2v(scf_env%krylov_space, ks, c0, c1, mo_eigenvalues, &
1664 592 : nao, eps_iter, ispin)
1665 592 : t2 = m_walltime()
1666 592 : IF (output_unit > 0) THEN
1667 : WRITE (output_unit, '(T8,I3,T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
1668 0 : ispin, iter, scf_env%krylov_space%nmo_conv, &
1669 0 : scf_env%krylov_space%max_res_norm, scf_env%krylov_space%min_res_norm, t2 - t1
1670 : END IF
1671 592 : t1 = m_walltime()
1672 620 : IF (scf_env%krylov_space%max_res_norm < eps_iter) THEN
1673 10 : converged = .TRUE.
1674 10 : IF (output_unit > 0) WRITE (output_unit, *) &
1675 0 : " Reached convergence in ", iter, " iterations "
1676 : EXIT
1677 : END IF
1678 : END DO
1679 :
1680 38 : IF (.NOT. converged .AND. output_unit > 0) THEN
1681 : WRITE (output_unit, "(T4, A)") " WARNING Lanczos refinement could "// &
1682 0 : "not converge all the mos:"
1683 0 : WRITE (output_unit, "(T40,A,T70,I10)") " number of not converged mos ", &
1684 0 : scf_env%krylov_space%nmo_nc
1685 0 : WRITE (output_unit, "(T40,A,T70,E10.2)") " max norm of the residual ", &
1686 0 : scf_env%krylov_space%max_res_norm
1687 :
1688 : END IF
1689 :
1690 : ! For the moment skip the re-orthogonalization
1691 : IF (.FALSE.) THEN
1692 : !Re-orthogonalization
1693 : NULLIFY (chc, evec)
1694 : chc => scf_env%krylov_space%chc_mat(ispin)
1695 : evec => scf_env%krylov_space%c_vec(ispin)
1696 : CALL parallel_gemm('N', 'N', nao, nmo, nao, rone, ks, c0, rzero, work)
1697 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, rone, c0, work, rzero, chc)
1698 : !Diagonalize (C^t)HC
1699 : CALL choose_eigv_solver(chc, evec, mo_eigenvalues)
1700 : !Rotate the C vectors
1701 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, rone, c0, evec, rzero, c1)
1702 : c0 => scf_env%krylov_space%mo_refine(ispin)
1703 : END IF
1704 :
1705 38 : IF (scf_env%cholesky_method == cholesky_inverse) THEN
1706 38 : CALL cp_fm_cholesky_restore(c0, nmo, ortho, mo_coeff, "MULTIPLY")
1707 : ELSE
1708 0 : CALL cp_fm_cholesky_restore(c0, nmo, ortho, mo_coeff, "SOLVE")
1709 : END IF
1710 :
1711 : CALL set_mo_occupation(mo_set=mos(ispin), &
1712 38 : smear=scf_control%smear)
1713 :
1714 : ! does not yet handle k-points
1715 : CALL calculate_density_matrix(mos(ispin), &
1716 38 : scf_env%p_mix_new(ispin, 1)%matrix)
1717 : END IF
1718 : END DO ! ispin
1719 :
1720 38 : IF (output_unit > 0) THEN
1721 0 : WRITE (output_unit, "(T15,A/)") '<<<<<<<<< END LANCZOS REFINEMENT <<<<<<<<<<'
1722 : END IF
1723 :
1724 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
1725 38 : "PRINT%LANCZOS")
1726 :
1727 38 : CALL timestop(handle)
1728 :
1729 38 : END SUBROUTINE do_block_krylov_diag
1730 :
1731 : ! **************************************************************************************************
1732 : !> \brief iterative diagonalization using the block davidson space approach
1733 : !> \param qs_env ...
1734 : !> \param scf_env ...
1735 : !> \param mos ...
1736 : !> \param matrix_ks ...
1737 : !> \param matrix_s ...
1738 : !> \param scf_control ...
1739 : !> \param scf_section ...
1740 : !> \param check_moconv_only ...
1741 : !> \param
1742 : !> \par History
1743 : !> 05.2011 created [MI]
1744 : ! **************************************************************************************************
1745 :
1746 84 : SUBROUTINE do_block_davidson_diag(qs_env, scf_env, mos, matrix_ks, matrix_s, &
1747 : scf_control, scf_section, check_moconv_only)
1748 :
1749 : TYPE(qs_environment_type), POINTER :: qs_env
1750 : TYPE(qs_scf_env_type), POINTER :: scf_env
1751 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
1752 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
1753 : TYPE(scf_control_type), POINTER :: scf_control
1754 : TYPE(section_vals_type), POINTER :: scf_section
1755 : LOGICAL, INTENT(IN), OPTIONAL :: check_moconv_only
1756 :
1757 : CHARACTER(LEN=*), PARAMETER :: routineN = 'do_block_davidson_diag'
1758 :
1759 : INTEGER :: handle, ispin, nspins, output_unit
1760 : LOGICAL :: do_prec, my_check_moconv_only
1761 : TYPE(cp_logger_type), POINTER :: logger
1762 :
1763 84 : logger => cp_get_default_logger()
1764 84 : CALL timeset(routineN, handle)
1765 :
1766 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DAVIDSON", &
1767 84 : extension=".scfLog")
1768 :
1769 84 : IF (output_unit > 0) &
1770 0 : WRITE (output_unit, "(/T15,A)") '<<<<<<<<< DAVIDSON ITERATIONS <<<<<<<<<<'
1771 :
1772 84 : IF (scf_env%mixing_method == 1) THEN
1773 0 : scf_env%iter_param = scf_env%p_mix_alpha
1774 0 : scf_env%iter_method = "P_Mix/Dav."
1775 : ELSE
1776 84 : scf_env%iter_param = scf_env%mixing_store%alpha
1777 84 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Dav."
1778 : END IF
1779 :
1780 84 : my_check_moconv_only = .FALSE.
1781 84 : IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
1782 84 : do_prec = .FALSE.
1783 84 : IF (scf_env%block_davidson_env(1)%prec_type /= 0 .AND. &
1784 : scf_env%iter_count >= scf_env%block_davidson_env(1)%first_prec) THEN
1785 76 : do_prec = .TRUE.
1786 : END IF
1787 :
1788 84 : nspins = SIZE(matrix_ks)
1789 :
1790 84 : IF (do_prec .AND. (scf_env%iter_count == scf_env%block_davidson_env(1)%first_prec .OR. &
1791 : MODULO(scf_env%iter_count, scf_env%block_davidson_env(1)%niter_new_prec) == 0)) THEN
1792 : CALL restart_preconditioner(qs_env, scf_env%ot_preconditioner, &
1793 16 : prec_type=scf_env%block_davidson_env(1)%prec_type, nspins=nspins)
1794 : CALL prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, scf_env%ot_preconditioner, &
1795 : scf_env%block_davidson_env(1)%prec_type, &
1796 : scf_env%block_davidson_env(1)%solver_type, &
1797 : scf_env%block_davidson_env(1)%energy_gap, nspins, &
1798 : convert_to_dbcsr=scf_env%block_davidson_env(1)%use_sparse_mos, &
1799 16 : full_mo_set=.TRUE.)
1800 : END IF
1801 :
1802 178 : DO ispin = 1, nspins
1803 178 : IF (scf_env%block_davidson_env(ispin)%use_sparse_mos) THEN
1804 64 : IF (.NOT. do_prec) THEN
1805 : CALL generate_extended_space_sparse(scf_env%block_davidson_env(ispin), mos(ispin), &
1806 8 : matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit)
1807 : ELSE
1808 : CALL generate_extended_space_sparse(scf_env%block_davidson_env(ispin), mos(ispin), &
1809 : matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit, &
1810 56 : scf_env%ot_preconditioner(ispin)%preconditioner)
1811 : END IF
1812 :
1813 : ELSE
1814 30 : IF (.NOT. do_prec) THEN
1815 : CALL generate_extended_space(scf_env%block_davidson_env(ispin), mos(ispin), &
1816 2 : matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit)
1817 : ELSE
1818 : CALL generate_extended_space(scf_env%block_davidson_env(ispin), mos(ispin), &
1819 : matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit, &
1820 28 : scf_env%ot_preconditioner(ispin)%preconditioner)
1821 : END IF
1822 : END IF
1823 : END DO !ispin
1824 :
1825 : CALL set_mo_occupation(mo_array=mos, &
1826 84 : smear=scf_control%smear)
1827 :
1828 178 : DO ispin = 1, nspins
1829 : ! does not yet handle k-points
1830 : CALL calculate_density_matrix(mos(ispin), &
1831 178 : scf_env%p_mix_new(ispin, 1)%matrix)
1832 : END DO
1833 :
1834 84 : IF (output_unit > 0) THEN
1835 0 : WRITE (output_unit, "(T15,A/)") '<<<<<<<<< END DAVIDSON ITERATION <<<<<<<<<<'
1836 : END IF
1837 :
1838 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
1839 84 : "PRINT%DAVIDSON")
1840 :
1841 84 : CALL timestop(handle)
1842 :
1843 84 : END SUBROUTINE do_block_davidson_diag
1844 :
1845 : ! **************************************************************************************************
1846 : !> \brief Kpoint diagonalization routine
1847 : !> Transforms matrices to kpoint, distributes kpoint groups, performs diagonalization
1848 : !> \param matrix_s Overlap matrices (RS indices, global)
1849 : !> \param kpoints Kpoint environment
1850 : !> \param fmwork full matrices distributed over all groups
1851 : !> \par History
1852 : !> 02.2026 created [JGH]
1853 : ! **************************************************************************************************
1854 0 : SUBROUTINE diag_kp_smat(matrix_s, kpoints, fmwork)
1855 :
1856 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
1857 : TYPE(kpoint_type), POINTER :: kpoints
1858 : TYPE(cp_fm_type), DIMENSION(:) :: fmwork
1859 :
1860 : CHARACTER(len=*), PARAMETER :: routineN = 'diag_kp_smat'
1861 : COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp, 0.0_dp), &
1862 : czero = (0.0_dp, 0.0_dp)
1863 :
1864 0 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ceig
1865 : INTEGER :: handle, igroup, ik, ikp, indx, kplocal, &
1866 : nao, nkp, nkp_groups
1867 : INTEGER, DIMENSION(2) :: kp_range
1868 0 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1869 0 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1870 : LOGICAL :: my_kpgrp, use_real_wfn
1871 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1872 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1873 0 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
1874 : TYPE(cp_cfm_type) :: csmat, cwork
1875 0 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools
1876 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1877 : TYPE(cp_fm_type) :: fmdummy, fmlocal, rsmat
1878 : TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix, tmpmat
1879 : TYPE(kpoint_env_type), POINTER :: kp
1880 : TYPE(mp_para_env_type), POINTER :: para_env
1881 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1882 0 : POINTER :: sab_nl
1883 : TYPE(qs_matrix_pools_type), POINTER :: mpools
1884 :
1885 0 : CALL timeset(routineN, handle)
1886 :
1887 0 : NULLIFY (sab_nl)
1888 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
1889 : nkp_groups=nkp_groups, kp_dist=kp_dist, sab_nl=sab_nl, &
1890 0 : cell_to_index=cell_to_index)
1891 0 : CPASSERT(ASSOCIATED(sab_nl))
1892 0 : kplocal = kp_range(2) - kp_range(1) + 1
1893 :
1894 : ! allocate some work matrices
1895 0 : ALLOCATE (rmatrix, cmatrix, tmpmat)
1896 : CALL dbcsr_create(rmatrix, template=matrix_s(1, 1)%matrix, &
1897 0 : matrix_type=dbcsr_type_symmetric)
1898 : CALL dbcsr_create(cmatrix, template=matrix_s(1, 1)%matrix, &
1899 0 : matrix_type=dbcsr_type_antisymmetric)
1900 : CALL dbcsr_create(tmpmat, template=matrix_s(1, 1)%matrix, &
1901 0 : matrix_type=dbcsr_type_no_symmetry)
1902 0 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
1903 0 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
1904 :
1905 : ! fm pools to be used within a kpoint group
1906 0 : CALL get_kpoint_info(kpoints, mpools=mpools)
1907 0 : CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
1908 :
1909 0 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
1910 0 : CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
1911 :
1912 0 : IF (use_real_wfn) THEN
1913 0 : CALL cp_fm_create(rsmat, matrix_struct)
1914 : ELSE
1915 0 : CALL cp_cfm_create(csmat, matrix_struct)
1916 0 : CALL cp_cfm_create(cwork, matrix_struct)
1917 : END IF
1918 :
1919 0 : CALL cp_fm_get_info(fmwork(1), nrow_global=nao)
1920 0 : ALLOCATE (eigenvalues(nao), ceig(nao))
1921 :
1922 0 : para_env => kpoints%blacs_env_all%para_env
1923 0 : ALLOCATE (info(kplocal*nkp_groups, 2))
1924 :
1925 : ! Setup and start all the communication
1926 0 : indx = 0
1927 0 : DO ikp = 1, kplocal
1928 0 : DO igroup = 1, nkp_groups
1929 : ! number of current kpoint
1930 0 : ik = kp_dist(1, igroup) + ikp - 1
1931 0 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
1932 0 : indx = indx + 1
1933 0 : IF (use_real_wfn) THEN
1934 0 : CALL dbcsr_set(rmatrix, 0.0_dp)
1935 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_s, ispin=1, &
1936 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1937 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1938 0 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
1939 : ELSE
1940 0 : CALL dbcsr_set(rmatrix, 0.0_dp)
1941 0 : CALL dbcsr_set(cmatrix, 0.0_dp)
1942 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_s, ispin=1, &
1943 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1944 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1945 0 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
1946 0 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
1947 0 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(2))
1948 : END IF
1949 : ! transfer to kpoint group
1950 : ! redistribution of matrices, new blacs environment
1951 : ! fmwork -> fmlocal -> rsmat/csmat
1952 0 : IF (use_real_wfn) THEN
1953 0 : IF (my_kpgrp) THEN
1954 0 : CALL cp_fm_start_copy_general(fmwork(1), rsmat, para_env, info(indx, 1))
1955 : ELSE
1956 0 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
1957 : END IF
1958 : ELSE
1959 0 : IF (my_kpgrp) THEN
1960 0 : CALL cp_fm_start_copy_general(fmwork(1), fmlocal, para_env, info(indx, 1))
1961 0 : CALL cp_fm_start_copy_general(fmwork(2), fmlocal, para_env, info(indx, 2))
1962 : ELSE
1963 0 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
1964 0 : CALL cp_fm_start_copy_general(fmwork(2), fmdummy, para_env, info(indx, 2))
1965 : END IF
1966 : END IF
1967 : END DO
1968 : END DO
1969 :
1970 : ! Finish communication then diagonalise in each group
1971 : indx = 0
1972 0 : DO ikp = 1, kplocal
1973 0 : DO igroup = 1, nkp_groups
1974 : ! number of current kpoint
1975 0 : ik = kp_dist(1, igroup) + ikp - 1
1976 0 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
1977 0 : indx = indx + 1
1978 0 : IF (my_kpgrp) THEN
1979 0 : IF (use_real_wfn) THEN
1980 0 : CALL cp_fm_finish_copy_general(rsmat, info(indx, 1))
1981 : ELSE
1982 0 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
1983 0 : CALL cp_cfm_scale_and_add_fm(z_zero, csmat, z_one, fmlocal)
1984 0 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
1985 0 : CALL cp_cfm_scale_and_add_fm(z_one, csmat, gaussi, fmlocal)
1986 : END IF
1987 : END IF
1988 : END DO
1989 :
1990 : ! Each kpoint group has now information on a kpoint to be diagonalized
1991 : ! Eigensolver Hermite or Symmetric
1992 0 : kp => kpoints%kp_env(ikp)%kpoint_env
1993 0 : IF (use_real_wfn) THEN
1994 0 : CALL choose_eigv_solver(rsmat, fmlocal, eigenvalues)
1995 : ELSE
1996 0 : CALL cp_cfm_heevd(csmat, cwork, eigenvalues)
1997 : END IF
1998 0 : CPASSERT(ALL(eigenvalues(1:nao) >= 0.0_dp))
1999 0 : IF (use_real_wfn) THEN
2000 0 : CALL cp_fm_release(kp%shalf)
2001 0 : CALL cp_fm_create(kp%shalf, matrix_struct)
2002 0 : eigenvalues(1:nao) = SQRT(eigenvalues(1:nao))
2003 0 : CALL cp_fm_to_fm(fmlocal, rsmat)
2004 0 : CALL cp_fm_column_scale(rsmat, eigenvalues)
2005 : CALL parallel_gemm("N", "T", nao, nao, nao, 1.0_dp, rsmat, fmlocal, &
2006 0 : 0.0_dp, kp%shalf)
2007 : ELSE
2008 0 : CALL cp_cfm_release(kp%cshalf)
2009 0 : CALL cp_cfm_create(kp%cshalf, matrix_struct)
2010 0 : ceig(1:nao) = SQRT(eigenvalues(1:nao))
2011 0 : CALL cp_cfm_to_cfm(cwork, csmat)
2012 0 : CALL cp_cfm_column_scale(csmat, ceig)
2013 : CALL parallel_gemm("N", "T", nao, nao, nao, cone, csmat, cwork, &
2014 0 : czero, kp%cshalf)
2015 : END IF
2016 : END DO
2017 :
2018 : ! Clean up communication
2019 : indx = 0
2020 0 : DO ikp = 1, kplocal
2021 0 : DO igroup = 1, nkp_groups
2022 : ! number of current kpoint
2023 0 : ik = kp_dist(1, igroup) + ikp - 1
2024 0 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
2025 0 : indx = indx + 1
2026 0 : IF (use_real_wfn) THEN
2027 0 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
2028 : ELSE
2029 0 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
2030 0 : CALL cp_fm_cleanup_copy_general(info(indx, 2))
2031 : END IF
2032 : END DO
2033 : END DO
2034 :
2035 : ! All done
2036 0 : DEALLOCATE (info)
2037 0 : DEALLOCATE (eigenvalues, ceig)
2038 :
2039 0 : CALL dbcsr_deallocate_matrix(rmatrix)
2040 0 : CALL dbcsr_deallocate_matrix(cmatrix)
2041 0 : CALL dbcsr_deallocate_matrix(tmpmat)
2042 :
2043 0 : IF (use_real_wfn) THEN
2044 0 : CALL cp_fm_release(rsmat)
2045 : ELSE
2046 0 : CALL cp_cfm_release(csmat)
2047 0 : CALL cp_cfm_release(cwork)
2048 : END IF
2049 0 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
2050 :
2051 0 : CALL timestop(handle)
2052 :
2053 0 : END SUBROUTINE diag_kp_smat
2054 :
2055 : END MODULE qs_scf_diagonalization
|
