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