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