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