(git:eeadd9f)
Loading...
Searching...
No Matches
qs_scf_diagonalization.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Different diagonalization schemes that can be used
10!> for the iterative solution of the eigenvalue problem
11!> \par History
12!> started from routines previously located in the qs_scf module
13!> 05.2009
14! **************************************************************************************************
15MODULE qs_scf_diagonalization
16 USE cp_array_utils, ONLY: cp_1d_r_p_type
17 USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,&
18 cp_cfm_scale_and_add,&
19 cp_cfm_scale_and_add_fm
20 USE cp_cfm_diag, ONLY: cp_cfm_geeig,&
21 cp_cfm_geeig_canon,&
22 cp_cfm_heevd
23 USE cp_cfm_types, ONLY: cp_cfm_create,&
24 cp_cfm_release,&
25 cp_cfm_set_all,&
26 cp_cfm_to_cfm,&
27 cp_cfm_to_fm,&
28 cp_cfm_type
29 USE cp_control_types, ONLY: dft_control_type,&
30 hairy_probes_type
31 USE cp_dbcsr_api, ONLY: &
32 dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_p_type, &
33 dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
34 dbcsr_type_symmetric
35 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
36 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
37 copy_fm_to_dbcsr,&
38 cp_dbcsr_sm_fm_multiply,&
39 dbcsr_allocate_matrix_set
40 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
41 cp_fm_symm,&
42 cp_fm_uplo_to_full
43 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_reduce,&
44 cp_fm_cholesky_restore
45 USE cp_fm_diag, ONLY: fm_diag_type_cusolver,&
46 choose_eigv_solver,&
47 cp_fm_geeig,&
48 cp_fm_geeig_canon,&
49 cusolver_generalized,&
50 diag_type
51 USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
52 fm_pool_create_fm,&
53 fm_pool_give_back_fm
54 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
55 cp_fm_struct_release,&
56 cp_fm_struct_type
57 USE cp_fm_types, ONLY: &
58 copy_info_type, cp_fm_add_to_element, cp_fm_cleanup_copy_general, cp_fm_create, &
59 cp_fm_finish_copy_general, cp_fm_get_info, cp_fm_release, cp_fm_start_copy_general, &
60 cp_fm_to_fm, cp_fm_type
61 USE cp_log_handling, ONLY: cp_get_default_logger,&
62 cp_logger_type
63 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
64 cp_print_key_unit_nr
65 USE input_constants, ONLY: &
66 cholesky_dbcsr, cholesky_inverse, cholesky_off, cholesky_reduce, cholesky_restore, &
67 core_guess, general_roks, high_spin_roks, restart_guess
68 USE input_section_types, ONLY: section_vals_get_subs_vals,&
69 section_vals_type
70 USE kinds, ONLY: dp
71 USE kpoint_methods, ONLY: kpoint_density_matrices,&
72 kpoint_density_transform,&
73 kpoint_set_mo_occupation,&
74 rskp_transform
75 USE kpoint_types, ONLY: get_kpoint_info,&
76 kpoint_env_type,&
77 kpoint_type
78 USE machine, ONLY: m_flush,&
79 m_walltime
80 USE mathconstants, ONLY: gaussi,&
81 z_one,&
82 z_zero
83 USE message_passing, ONLY: mp_para_env_type
84 USE parallel_gemm_api, ONLY: parallel_gemm
85 USE preconditioner, ONLY: prepare_preconditioner,&
86 restart_preconditioner
87 USE qs_density_matrices, ONLY: calculate_density_matrix
88 USE qs_density_mixing_types, ONLY: direct_mixing_nr,&
89 gspace_mixing_nr
90 USE qs_diis, ONLY: qs_diis_b_calc_err_kp,&
91 qs_diis_b_info_kp,&
92 qs_diis_b_step,&
93 qs_diis_b_step_kp
94 USE qs_energy_types, ONLY: qs_energy_type
95 USE qs_environment_types, ONLY: get_qs_env,&
96 qs_environment_type
97 USE qs_gspace_mixing, ONLY: gspace_mixing
98 USE qs_ks_methods, ONLY: qs_ks_update_qs_env
99 USE qs_ks_types, ONLY: qs_ks_did_change,&
100 qs_ks_env_type
101 USE qs_matrix_pools, ONLY: mpools_get,&
102 qs_matrix_pools_type
103 USE qs_mixing_utils, ONLY: charge_mixing_init,&
104 mixing_allocate,&
105 mixing_init,&
106 self_consistency_check
107 USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
108 USE qs_mo_occupation, ONLY: set_mo_occupation
109 USE qs_mo_types, ONLY: get_mo_set,&
110 mo_set_type
111 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
112 USE qs_ot_eigensolver, ONLY: ot_eigensolver
113 USE qs_rho_atom_types, ONLY: rho_atom_type
114 USE qs_rho_methods, ONLY: qs_rho_update_rho
115 USE qs_rho_types, ONLY: qs_rho_get,&
116 qs_rho_type
117 USE qs_scf_block_davidson, ONLY: generate_extended_space,&
118 generate_extended_space_sparse
119 USE qs_scf_lanczos, ONLY: lanczos_refinement,&
120 lanczos_refinement_2v
121 USE qs_scf_methods, ONLY: combine_ks_matrices,&
122 eigensolver,&
123 eigensolver_dbcsr,&
124 eigensolver_generalized,&
125 eigensolver_simple,&
126 eigensolver_symm,&
127 scf_env_density_mixing
128 USE qs_scf_types, ONLY: qs_scf_env_type,&
129 subspace_env_type
130 USE scf_control_types, ONLY: scf_control_type
131#include "./base/base_uses.f90"
132
133 IMPLICIT NONE
134
135 PRIVATE
136
137 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_diagonalization'
138
139 PUBLIC :: do_general_diag, do_general_diag_kp, do_roks_diag, &
140 do_special_diag, do_ot_diag, do_block_davidson_diag, &
141 do_block_krylov_diag, do_scf_diag_subspace, diag_subspace_allocate, &
142 general_eigenproblem, diag_kp_smat
143
144CONTAINS
145
146! **************************************************************************************************
147!> \brief the inner loop of scf, specific to diagonalization with S matrix
148!> basically, in goes the ks matrix out goes a new p matrix
149!> \param scf_env ...
150!> \param mos ...
151!> \param matrix_ks ...
152!> \param matrix_s ...
153!> \param scf_control ...
154!> \param scf_section ...
155!> \param diis_step ...
156!> \par History
157!> 03.2006 created [Joost VandeVondele]
158! **************************************************************************************************
159
160 SUBROUTINE general_eigenproblem(scf_env, mos, matrix_ks, &
161 matrix_s, scf_control, scf_section, &
162 diis_step)
163
164 TYPE(qs_scf_env_type), POINTER :: scf_env
165 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
166 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
167 TYPE(scf_control_type), POINTER :: scf_control
168 TYPE(section_vals_type), POINTER :: scf_section
169 LOGICAL, INTENT(INOUT) :: diis_step
170
171 INTEGER :: ispin, nspin
172 LOGICAL :: do_level_shift, owns_ortho, use_jacobi
173 REAL(kind=dp) :: diis_error, eps_diis
174 TYPE(cp_fm_type), POINTER :: ortho
175 TYPE(dbcsr_type), POINTER :: ortho_dbcsr
176
177 nspin = SIZE(matrix_ks)
178 NULLIFY (ortho, ortho_dbcsr)
179
180 DO ispin = 1, nspin
181 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, &
182 scf_env%scf_work1(ispin))
183 END DO
184
185 eps_diis = scf_control%eps_diis
186
187 IF (scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis) THEN
188 CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
189 scf_env%scf_work2, scf_env%iter_delta, diis_error, diis_step, &
190 eps_diis, scf_control%nmixing, &
191 s_matrix=matrix_s, &
192 scf_section=scf_section)
193 ELSE
194 diis_step = .false.
195 END IF
196
197 do_level_shift = ((scf_control%level_shift /= 0.0_dp) .AND. &
198 ((scf_control%density_guess == core_guess) .OR. &
199 (scf_env%iter_count > 1)))
200
201 IF ((scf_env%iter_count > 1) .AND. &
202 (scf_env%iter_delta < scf_control%diagonalization%eps_jacobi)) THEN
203 use_jacobi = .true.
204 ELSE
205 use_jacobi = .false.
206 END IF
207
208 IF (diis_step) THEN
209 scf_env%iter_param = diis_error
210 IF (use_jacobi) THEN
211 scf_env%iter_method = "DIIS/Jacobi"
212 ELSE
213 scf_env%iter_method = "DIIS/Diag."
214 END IF
215 ELSE
216 IF (scf_env%mixing_method == 0) THEN
217 scf_env%iter_method = "NoMix/Diag."
218 ELSE IF (scf_env%mixing_method == 1) THEN
219 scf_env%iter_param = scf_env%p_mix_alpha
220 IF (use_jacobi) THEN
221 scf_env%iter_method = "P_Mix/Jacobi"
222 ELSE
223 scf_env%iter_method = "P_Mix/Diag."
224 END IF
225 ELSEIF (scf_env%mixing_method > 1) THEN
226 scf_env%iter_param = scf_env%mixing_store%alpha
227 IF (use_jacobi) THEN
228 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Jacobi"
229 ELSE
230 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Diag."
231 END IF
232 END IF
233 END IF
234
235 IF (scf_env%cholesky_method == cholesky_dbcsr) THEN
236 ortho_dbcsr => scf_env%ortho_dbcsr
237 DO ispin = 1, nspin
238 CALL eigensolver_dbcsr(matrix_ks=matrix_ks(ispin)%matrix, matrix_ks_fm=scf_env%scf_work1(ispin), &
239 mo_set=mos(ispin), &
240 ortho_dbcsr=ortho_dbcsr, &
241 ksbuf1=scf_env%buf1_dbcsr, ksbuf2=scf_env%buf2_dbcsr)
242 END DO
243
244 ELSE IF (scf_env%cholesky_method > cholesky_off) THEN
245 IF (scf_env%cholesky_method == cholesky_inverse) THEN
246 ortho => scf_env%ortho_m1
247 ELSE
248 ortho => scf_env%ortho
249 END IF
250
251 owns_ortho = .false.
252 IF (.NOT. ASSOCIATED(ortho)) THEN
253 ALLOCATE (ortho)
254 owns_ortho = .true.
255 END IF
256
257 DO ispin = 1, nspin
258 IF (diag_type == fm_diag_type_cusolver .AND. cusolver_generalized .AND. .NOT. do_level_shift) THEN
259 CALL eigensolver_generalized(matrix_ks_fm=scf_env%scf_work1(ispin), &
260 matrix_s=matrix_s(ispin)%matrix, &
261 mo_set=mos(ispin), &
262 work=scf_env%scf_work2)
263 ELSE
264 IF (do_level_shift) THEN
265 CALL eigensolver(matrix_ks_fm=scf_env%scf_work1(ispin), &
266 mo_set=mos(ispin), &
267 ortho=ortho, &
268 work=scf_env%scf_work2, &
269 cholesky_method=scf_env%cholesky_method, &
270 do_level_shift=do_level_shift, &
271 level_shift=scf_control%level_shift, &
272 matrix_u_fm=scf_env%ortho, &
273 use_jacobi=use_jacobi)
274 ELSE
275 CALL eigensolver(matrix_ks_fm=scf_env%scf_work1(ispin), &
276 mo_set=mos(ispin), &
277 ortho=ortho, &
278 work=scf_env%scf_work2, &
279 cholesky_method=scf_env%cholesky_method, &
280 do_level_shift=do_level_shift, &
281 level_shift=scf_control%level_shift, &
282 use_jacobi=use_jacobi)
283 END IF
284 END IF
285 END DO
286
287 IF (owns_ortho) DEALLOCATE (ortho)
288 ELSE
289 ortho => scf_env%ortho
290
291 owns_ortho = .false.
292 IF (.NOT. ASSOCIATED(ortho)) THEN
293 ALLOCATE (ortho)
294 owns_ortho = .true.
295 END IF
296
297 IF (do_level_shift) THEN
298 DO ispin = 1, nspin
299 IF (ASSOCIATED(scf_env%scf_work1_red) .AND. ASSOCIATED(scf_env%scf_work2_red) &
300 .AND. ASSOCIATED(scf_env%ortho_red) .AND. ASSOCIATED(scf_env%ortho_m1_red)) THEN
301 CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
302 mo_set=mos(ispin), &
303 ortho=ortho, &
304 work=scf_env%scf_work2, &
305 do_level_shift=do_level_shift, &
306 level_shift=scf_control%level_shift, &
307 matrix_u_fm=scf_env%ortho_m1, &
308 use_jacobi=use_jacobi, &
309 jacobi_threshold=scf_control%diagonalization%jacobi_threshold, &
310 matrix_ks_fm_red=scf_env%scf_work1_red(ispin), &
311 ortho_red=scf_env%ortho_red, &
312 work_red=scf_env%scf_work2_red, &
313 matrix_u_fm_red=scf_env%ortho_m1_red)
314 ELSE
315 CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
316 mo_set=mos(ispin), &
317 ortho=ortho, &
318 work=scf_env%scf_work2, &
319 do_level_shift=do_level_shift, &
320 level_shift=scf_control%level_shift, &
321 matrix_u_fm=scf_env%ortho_m1, &
322 use_jacobi=use_jacobi, &
323 jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
324 END IF
325 END DO
326 ELSE
327 DO ispin = 1, nspin
328 IF (ASSOCIATED(scf_env%scf_work1_red) .AND. ASSOCIATED(scf_env%scf_work2_red) &
329 .AND. ASSOCIATED(scf_env%ortho_red)) THEN
330 CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
331 mo_set=mos(ispin), &
332 ortho=ortho, &
333 work=scf_env%scf_work2, &
334 do_level_shift=do_level_shift, &
335 level_shift=scf_control%level_shift, &
336 use_jacobi=use_jacobi, &
337 jacobi_threshold=scf_control%diagonalization%jacobi_threshold, &
338 matrix_ks_fm_red=scf_env%scf_work1_red(ispin), &
339 ortho_red=scf_env%ortho_red, &
340 work_red=scf_env%scf_work2_red)
341 ELSE
342 CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
343 mo_set=mos(ispin), &
344 ortho=ortho, &
345 work=scf_env%scf_work2, &
346 do_level_shift=do_level_shift, &
347 level_shift=scf_control%level_shift, &
348 use_jacobi=use_jacobi, &
349 jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
350 END IF
351 END DO
352 END IF
353
354 IF (owns_ortho) DEALLOCATE (ortho)
355 END IF
356
357 END SUBROUTINE general_eigenproblem
358
359! **************************************************************************************************
360!> \brief ...
361!> \param scf_env ...
362!> \param mos ...
363!> \param matrix_ks ...
364!> \param matrix_s ...
365!> \param scf_control ...
366!> \param scf_section ...
367!> \param diis_step ...
368!> \param probe ...
369! **************************************************************************************************
370 SUBROUTINE do_general_diag(scf_env, mos, matrix_ks, &
371 matrix_s, scf_control, scf_section, &
372 diis_step, probe)
373
374 TYPE(qs_scf_env_type), POINTER :: scf_env
375 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
376 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
377 TYPE(scf_control_type), POINTER :: scf_control
378 TYPE(section_vals_type), POINTER :: scf_section
379 LOGICAL, INTENT(INOUT) :: diis_step
380 TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
381 POINTER :: probe
382
383 INTEGER :: ispin, nspin
384 REAL(kind=dp) :: total_zeff_corr
385
386 nspin = SIZE(matrix_ks)
387
388 CALL general_eigenproblem(scf_env, mos, matrix_ks, &
389 matrix_s, scf_control, scf_section, diis_step)
390
391 total_zeff_corr = 0.0_dp
392 total_zeff_corr = scf_env%sum_zeff_corr
393
394 IF (abs(total_zeff_corr) > 0.0_dp) THEN
395 CALL set_mo_occupation(mo_array=mos, &
396 smear=scf_control%smear, tot_zeff_corr=total_zeff_corr)
397 ELSE
398 IF (PRESENT(probe) .EQV. .true.) THEN
399 scf_control%smear%do_smear = .false.
400 CALL set_mo_occupation(mo_array=mos, &
401 smear=scf_control%smear, &
402 probe=probe)
403 ELSE
404 CALL set_mo_occupation(mo_array=mos, &
405 smear=scf_control%smear)
406 END IF
407 END IF
408
409 DO ispin = 1, nspin
410 CALL calculate_density_matrix(mos(ispin), &
411 scf_env%p_mix_new(ispin, 1)%matrix)
412 END DO
413
414 END SUBROUTINE do_general_diag
415
416! **************************************************************************************************
417!> \brief Kpoint diagonalization routine
418!> Transforms matrices to kpoint, distributes kpoint groups, performs
419!> general diagonalization (no storgae of overlap decomposition), stores
420!> MOs, calculates occupation numbers, calculates density matrices
421!> in kpoint representation, transforms density matrices to real space
422!> \param matrix_ks Kohn-sham matrices (RS indices, global)
423!> \param matrix_s Overlap matrices (RS indices, global)
424!> \param kpoints Kpoint environment
425!> \param scf_env SCF environment
426!> \param scf_control SCF control variables
427!> \param update_p ...
428!> \param diis_step ...
429!> \param diis_error ...
430!> \param qs_env ...
431!> \param probe ...
432!> \par History
433!> 08.2014 created [JGH]
434! **************************************************************************************************
435 SUBROUTINE do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, update_p, &
436 diis_step, diis_error, qs_env, probe)
437
438 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
439 TYPE(kpoint_type), POINTER :: kpoints
440 TYPE(qs_scf_env_type), POINTER :: scf_env
441 TYPE(scf_control_type), POINTER :: scf_control
442 LOGICAL, INTENT(IN) :: update_p
443 LOGICAL, INTENT(INOUT) :: diis_step
444 REAL(dp), INTENT(INOUT), OPTIONAL :: diis_error
445 TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
446 TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
447 POINTER :: probe
448
449 CHARACTER(len=*), PARAMETER :: routinen = 'do_general_diag_kp'
450
451 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coeffs
452 INTEGER :: handle, ib, igroup, ik, ikp, indx, &
453 ispin, jb, kplocal, nb, nkp, &
454 nkp_groups, nspin
455 INTEGER, DIMENSION(2) :: kp_range
456 INTEGER, DIMENSION(:, :), POINTER :: kp_dist
457 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
458 LOGICAL :: do_diis, my_kpgrp, use_real_wfn
459 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
460 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
461 TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
462 TYPE(cp_cfm_type) :: cksmat, cmos, csmat, cwork
463 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools
464 TYPE(cp_fm_struct_type), POINTER :: matrix_struct, mo_struct
465 TYPE(cp_fm_type) :: fmdummy, fmlocal, rksmat, rsmat
466 TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
467 TYPE(cp_fm_type), POINTER :: imos, mo_coeff, rmos
468 TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix, tmpmat
469 TYPE(kpoint_env_type), POINTER :: kp
470 TYPE(mp_para_env_type), POINTER :: para_env, para_env_global
471 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
472 POINTER :: sab_nl
473 TYPE(qs_matrix_pools_type), POINTER :: mpools
474 TYPE(section_vals_type), POINTER :: scf_section
475
476 CALL timeset(routinen, handle)
477
478 NULLIFY (sab_nl)
479 CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
480 nkp_groups=nkp_groups, kp_dist=kp_dist, sab_nl=sab_nl, &
481 cell_to_index=cell_to_index)
482 cpassert(ASSOCIATED(sab_nl))
483 kplocal = kp_range(2) - kp_range(1) + 1
484
485 !Whether we use DIIS for k-points
486 do_diis = .false.
487 IF (scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis .AND. .NOT. use_real_wfn &
488 .AND. PRESENT(diis_error) .AND. PRESENT(qs_env)) do_diis = .true.
489
490 ! allocate some work matrices
491 ALLOCATE (rmatrix, cmatrix, tmpmat)
492 CALL dbcsr_create(rmatrix, template=matrix_ks(1, 1)%matrix, &
493 matrix_type=dbcsr_type_symmetric)
494 CALL dbcsr_create(cmatrix, template=matrix_ks(1, 1)%matrix, &
495 matrix_type=dbcsr_type_antisymmetric)
496 CALL dbcsr_create(tmpmat, template=matrix_ks(1, 1)%matrix, &
497 matrix_type=dbcsr_type_no_symmetry)
498 CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
499 CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
500
501 fmwork => scf_env%scf_work1
502
503 ! fm pools to be used within a kpoint group
504 CALL get_kpoint_info(kpoints, mpools=mpools)
505 CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
506
507 CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
508 CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
509
510 IF (use_real_wfn) THEN
511 CALL cp_fm_create(rksmat, matrix_struct)
512 CALL cp_fm_create(rsmat, matrix_struct)
513 ELSE
514 CALL cp_cfm_create(cksmat, matrix_struct)
515 CALL cp_cfm_create(csmat, matrix_struct)
516 CALL cp_cfm_create(cwork, matrix_struct)
517 kp => kpoints%kp_env(1)%kpoint_env
518 CALL get_mo_set(kp%mos(1, 1), mo_coeff=mo_coeff)
519 CALL cp_fm_get_info(mo_coeff, matrix_struct=mo_struct)
520 CALL cp_cfm_create(cmos, mo_struct)
521 END IF
522
523 para_env => kpoints%blacs_env_all%para_env
524 nspin = SIZE(matrix_ks, 1)
525 ALLOCATE (info(kplocal*nspin*nkp_groups, 4))
526
527 ! Setup and start all the communication
528 indx = 0
529 DO ikp = 1, kplocal
530 DO ispin = 1, nspin
531 DO igroup = 1, nkp_groups
532 ! number of current kpoint
533 ik = kp_dist(1, igroup) + ikp - 1
534 my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
535 indx = indx + 1
536 IF (use_real_wfn) THEN
537 ! FT of matrices KS and S, then transfer to FM type
538 CALL dbcsr_set(rmatrix, 0.0_dp)
539 CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_ks, ispin=ispin, &
540 xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
541 CALL dbcsr_desymmetrize(rmatrix, tmpmat)
542 CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
543 ! s matrix is not spin dependent
544 CALL dbcsr_set(rmatrix, 0.0_dp)
545 CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_s, ispin=1, &
546 xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
547 CALL dbcsr_desymmetrize(rmatrix, tmpmat)
548 CALL copy_dbcsr_to_fm(tmpmat, fmwork(3))
549 ELSE
550 ! FT of matrices KS and S, then transfer to FM type
551 CALL dbcsr_set(rmatrix, 0.0_dp)
552 CALL dbcsr_set(cmatrix, 0.0_dp)
553 CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_ks, ispin=ispin, &
554 xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
555 CALL dbcsr_desymmetrize(rmatrix, tmpmat)
556 CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
557 CALL dbcsr_desymmetrize(cmatrix, tmpmat)
558 CALL copy_dbcsr_to_fm(tmpmat, fmwork(2))
559 ! s matrix is not spin dependent, double the work
560 CALL dbcsr_set(rmatrix, 0.0_dp)
561 CALL dbcsr_set(cmatrix, 0.0_dp)
562 CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_s, ispin=1, &
563 xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
564 CALL dbcsr_desymmetrize(rmatrix, tmpmat)
565 CALL copy_dbcsr_to_fm(tmpmat, fmwork(3))
566 CALL dbcsr_desymmetrize(cmatrix, tmpmat)
567 CALL copy_dbcsr_to_fm(tmpmat, fmwork(4))
568 END IF
569 ! transfer to kpoint group
570 ! redistribution of matrices, new blacs environment
571 ! fmwork -> fmlocal -> rksmat/cksmat
572 ! fmwork -> fmlocal -> rsmat/csmat
573 IF (use_real_wfn) THEN
574 IF (my_kpgrp) THEN
575 CALL cp_fm_start_copy_general(fmwork(1), rksmat, para_env, info(indx, 1))
576 CALL cp_fm_start_copy_general(fmwork(3), rsmat, para_env, info(indx, 2))
577 ELSE
578 CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
579 CALL cp_fm_start_copy_general(fmwork(3), fmdummy, para_env, info(indx, 2))
580 END IF
581 ELSE
582 IF (my_kpgrp) THEN
583 CALL cp_fm_start_copy_general(fmwork(1), fmlocal, para_env, info(indx, 1))
584 CALL cp_fm_start_copy_general(fmwork(2), fmlocal, para_env, info(indx, 2))
585 CALL cp_fm_start_copy_general(fmwork(3), fmlocal, para_env, info(indx, 3))
586 CALL cp_fm_start_copy_general(fmwork(4), fmlocal, para_env, info(indx, 4))
587 ELSE
588 CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
589 CALL cp_fm_start_copy_general(fmwork(2), fmdummy, para_env, info(indx, 2))
590 CALL cp_fm_start_copy_general(fmwork(3), fmdummy, para_env, info(indx, 3))
591 CALL cp_fm_start_copy_general(fmwork(4), fmdummy, para_env, info(indx, 4))
592 END IF
593 END IF
594 END DO
595 END DO
596 END DO
597
598 ! Finish communication then diagonalise in each group
599 IF (do_diis) THEN
600 CALL get_qs_env(qs_env, para_env=para_env_global)
601 scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
602 CALL qs_diis_b_info_kp(kpoints%scf_diis_buffer, ib, nb)
603 indx = 0
604 DO ikp = 1, kplocal
605 DO ispin = 1, nspin
606 DO igroup = 1, nkp_groups
607 ! number of current kpoint
608 ik = kp_dist(1, igroup) + ikp - 1
609 my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
610 indx = indx + 1
611 IF (my_kpgrp) THEN
612 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
613 CALL cp_cfm_scale_and_add_fm(z_zero, cksmat, z_one, fmlocal)
614 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
615 CALL cp_cfm_scale_and_add_fm(z_one, cksmat, gaussi, fmlocal)
616 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 3))
617 CALL cp_cfm_scale_and_add_fm(z_zero, csmat, z_one, fmlocal)
618 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 4))
619 CALL cp_cfm_scale_and_add_fm(z_one, csmat, gaussi, fmlocal)
620 END IF
621 END DO !igroup
622
623 kp => kpoints%kp_env(ikp)%kpoint_env
624 CALL qs_diis_b_calc_err_kp(kpoints%scf_diis_buffer, ib, kp%mos, cksmat, csmat, &
625 ispin, ikp, kplocal, scf_section)
626
627 END DO !ispin
628 END DO !ikp
629
630 ALLOCATE (coeffs(nb))
631 CALL qs_diis_b_step_kp(kpoints%scf_diis_buffer, coeffs, ib, nb, scf_env%iter_delta, diis_error, &
632 diis_step, scf_control%eps_diis, nspin, nkp, kplocal, scf_control%nmixing, &
633 scf_section, para_env_global)
634
635 !build the ks matrices and idagonalize
636 DO ikp = 1, kplocal
637 DO ispin = 1, nspin
638 kp => kpoints%kp_env(ikp)%kpoint_env
639 CALL cp_cfm_to_cfm(kpoints%scf_diis_buffer%smat(ikp), csmat)
640
641 CALL cp_cfm_set_all(cksmat, z_zero)
642 DO jb = 1, nb
643 CALL cp_cfm_scale_and_add(z_one, cksmat, coeffs(jb), kpoints%scf_diis_buffer%param(jb, ispin, ikp))
644 END DO
645
646 CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
647 CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
648 IF (scf_env%cholesky_method == cholesky_off) THEN
649 CALL cp_cfm_geeig_canon(cksmat, csmat, cmos, eigenvalues, cwork, &
650 scf_control%eps_eigval)
651 ELSE
652 CALL cp_cfm_geeig(cksmat, csmat, cmos, eigenvalues, cwork)
653 END IF
654 ! copy eigenvalues to imag set (keep them in sync)
655 kp%mos(2, ispin)%eigenvalues = eigenvalues
656 ! split real and imaginary part of mos
657 CALL cp_cfm_to_fm(cmos, rmos, imos)
658 END DO
659 END DO
660
661 ELSE !no DIIS
662 diis_step = .false.
663 indx = 0
664 DO ikp = 1, kplocal
665 DO ispin = 1, nspin
666 DO igroup = 1, nkp_groups
667 ! number of current kpoint
668 ik = kp_dist(1, igroup) + ikp - 1
669 my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
670 indx = indx + 1
671 IF (my_kpgrp) THEN
672 IF (use_real_wfn) THEN
673 CALL cp_fm_finish_copy_general(rksmat, info(indx, 1))
674 CALL cp_fm_finish_copy_general(rsmat, info(indx, 2))
675 ELSE
676 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
677 CALL cp_cfm_scale_and_add_fm(z_zero, cksmat, z_one, fmlocal)
678 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
679 CALL cp_cfm_scale_and_add_fm(z_one, cksmat, gaussi, fmlocal)
680 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 3))
681 CALL cp_cfm_scale_and_add_fm(z_zero, csmat, z_one, fmlocal)
682 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 4))
683 CALL cp_cfm_scale_and_add_fm(z_one, csmat, gaussi, fmlocal)
684 END IF
685 END IF
686 END DO
687
688 ! Each kpoint group has now information on a kpoint to be diagonalized
689 ! General eigensolver Hermite or Symmetric
690 kp => kpoints%kp_env(ikp)%kpoint_env
691 IF (use_real_wfn) THEN
692 CALL get_mo_set(kp%mos(1, ispin), mo_coeff=mo_coeff, eigenvalues=eigenvalues)
693 IF (scf_env%cholesky_method == cholesky_off) THEN
694 CALL cp_fm_geeig_canon(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal, &
695 scf_control%eps_eigval)
696 ELSE
697 CALL cp_fm_geeig(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal)
698 END IF
699 ELSE
700 CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
701 CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
702 IF (scf_env%cholesky_method == cholesky_off) THEN
703 CALL cp_cfm_geeig_canon(cksmat, csmat, cmos, eigenvalues, cwork, &
704 scf_control%eps_eigval)
705 ELSE
706 CALL cp_cfm_geeig(cksmat, csmat, cmos, eigenvalues, cwork)
707 END IF
708 ! copy eigenvalues to imag set (keep them in sync)
709 kp%mos(2, ispin)%eigenvalues = eigenvalues
710 ! split real and imaginary part of mos
711 CALL cp_cfm_to_fm(cmos, rmos, imos)
712 END IF
713 END DO
714 END DO
715 END IF
716
717 ! Clean up communication
718 indx = 0
719 DO ikp = 1, kplocal
720 DO ispin = 1, nspin
721 DO igroup = 1, nkp_groups
722 ! number of current kpoint
723 ik = kp_dist(1, igroup) + ikp - 1
724 my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
725 indx = indx + 1
726 IF (use_real_wfn) THEN
727 CALL cp_fm_cleanup_copy_general(info(indx, 1))
728 CALL cp_fm_cleanup_copy_general(info(indx, 2))
729 ELSE
730 CALL cp_fm_cleanup_copy_general(info(indx, 1))
731 CALL cp_fm_cleanup_copy_general(info(indx, 2))
732 CALL cp_fm_cleanup_copy_general(info(indx, 3))
733 CALL cp_fm_cleanup_copy_general(info(indx, 4))
734 END IF
735 END DO
736 END DO
737 END DO
738
739 ! All done
740 DEALLOCATE (info)
741
742 IF (update_p) THEN
743 ! MO occupations
744 IF (PRESENT(probe) .EQV. .true.) THEN
745 scf_control%smear%do_smear = .false.
746 CALL kpoint_set_mo_occupation(kpoints, scf_control%smear, &
747 probe=probe)
748 ELSE
749 CALL kpoint_set_mo_occupation(kpoints, scf_control%smear)
750 END IF
751 ! density matrices
752 CALL kpoint_density_matrices(kpoints)
753 ! density matrices in real space
754 CALL kpoint_density_transform(kpoints, scf_env%p_mix_new, .false., &
755 matrix_s(1, 1)%matrix, sab_nl, fmwork)
756 END IF
757
758 CALL dbcsr_deallocate_matrix(rmatrix)
759 CALL dbcsr_deallocate_matrix(cmatrix)
760 CALL dbcsr_deallocate_matrix(tmpmat)
761
762 IF (use_real_wfn) THEN
763 CALL cp_fm_release(rksmat)
764 CALL cp_fm_release(rsmat)
765 ELSE
766 CALL cp_cfm_release(cksmat)
767 CALL cp_cfm_release(csmat)
768 CALL cp_cfm_release(cwork)
769 CALL cp_cfm_release(cmos)
770 END IF
771 CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
772
773 CALL timestop(handle)
774
775 END SUBROUTINE do_general_diag_kp
776
777! **************************************************************************************************
778!> \brief inner loop within MOS subspace, to refine occupation and density,
779!> before next diagonalization of the Hamiltonian
780!> \param qs_env ...
781!> \param scf_env ...
782!> \param subspace_env ...
783!> \param mos ...
784!> \param rho ...
785!> \param ks_env ...
786!> \param scf_section ...
787!> \param scf_control ...
788!> \par History
789!> 09.2009 created [MI]
790!> \note it is assumed that when diagonalization is used, also some mixing procedure is active
791! **************************************************************************************************
792 SUBROUTINE do_scf_diag_subspace(qs_env, scf_env, subspace_env, mos, rho, &
793 ks_env, scf_section, scf_control)
794
795 TYPE(qs_environment_type), POINTER :: qs_env
796 TYPE(qs_scf_env_type), POINTER :: scf_env
797 TYPE(subspace_env_type), POINTER :: subspace_env
798 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
799 TYPE(qs_rho_type), POINTER :: rho
800 TYPE(qs_ks_env_type), POINTER :: ks_env
801 TYPE(section_vals_type), POINTER :: scf_section
802 TYPE(scf_control_type), POINTER :: scf_control
803
804 CHARACTER(LEN=*), PARAMETER :: routinen = 'do_scf_diag_subspace'
805 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
806
807 INTEGER :: handle, i, iloop, ispin, nao, nmo, &
808 nspin, output_unit
809 LOGICAL :: converged
810 REAL(dp) :: ene_diff, ene_old, iter_delta, max_val, &
811 sum_band, sum_val, t1, t2
812 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues, mo_occupations
813 TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: eval_first, occ_first
814 TYPE(cp_fm_type) :: work
815 TYPE(cp_fm_type), POINTER :: c0, chc, evec, mo_coeff
816 TYPE(cp_logger_type), POINTER :: logger
817 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
818 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
819 TYPE(dft_control_type), POINTER :: dft_control
820 TYPE(mp_para_env_type), POINTER :: para_env
821 TYPE(qs_energy_type), POINTER :: energy
822 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
823
824 CALL timeset(routinen, handle)
825 NULLIFY (c0, chc, energy, evec, matrix_ks, mo_coeff, mo_eigenvalues, &
826 mo_occupations, dft_control, rho_ao, rho_ao_kp)
827
828 logger => cp_get_default_logger()
829 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIAG_SUB_SCF", &
830 extension=".scfLog")
831
832 !Extra loop keeping mos unchanged and refining the subspace occupation
833 nspin = SIZE(mos)
834 CALL qs_rho_get(rho, rho_ao=rho_ao, rho_ao_kp=rho_ao_kp)
835
836 ALLOCATE (eval_first(nspin))
837 ALLOCATE (occ_first(nspin))
838 DO ispin = 1, nspin
839 CALL get_mo_set(mo_set=mos(ispin), &
840 nmo=nmo, &
841 eigenvalues=mo_eigenvalues, &
842 occupation_numbers=mo_occupations)
843 ALLOCATE (eval_first(ispin)%array(nmo))
844 ALLOCATE (occ_first(ispin)%array(nmo))
845 eval_first(ispin)%array(1:nmo) = mo_eigenvalues(1:nmo)
846 occ_first(ispin)%array(1:nmo) = mo_occupations(1:nmo)
847 END DO
848
849 DO ispin = 1, nspin
850 ! does not yet handle k-points
851 CALL dbcsr_copy(subspace_env%p_matrix_store(ispin)%matrix, rho_ao(ispin)%matrix)
852 CALL dbcsr_copy(rho_ao(ispin)%matrix, scf_env%p_mix_new(ispin, 1)%matrix)
853 END DO
854
855 subspace_env%p_matrix_mix => scf_env%p_mix_new
856
857 NULLIFY (matrix_ks, energy, para_env, matrix_s)
858 CALL get_qs_env(qs_env, &
859 matrix_ks=matrix_ks, &
860 energy=energy, &
861 matrix_s=matrix_s, &
862 para_env=para_env, &
863 dft_control=dft_control)
864
865 ! mixing storage allocation
866 IF (subspace_env%mixing_method >= gspace_mixing_nr) THEN
867 CALL mixing_allocate(qs_env, subspace_env%mixing_method, scf_env%p_mix_new, &
868 scf_env%p_delta, nspin, subspace_env%mixing_store)
869 IF (dft_control%qs_control%gapw) THEN
870 CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
871 CALL mixing_init(subspace_env%mixing_method, rho, subspace_env%mixing_store, &
872 para_env, rho_atom=rho_atom)
873 ELSEIF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
874 CALL charge_mixing_init(subspace_env%mixing_store)
875 ELSEIF (dft_control%qs_control%semi_empirical) THEN
876 cpabort('SE Code not possible')
877 ELSE
878 CALL mixing_init(subspace_env%mixing_method, rho, subspace_env%mixing_store, para_env)
879 END IF
880 END IF
881
882 ene_old = 0.0_dp
883 ene_diff = 0.0_dp
884 IF (output_unit > 0) THEN
885 WRITE (output_unit, "(/T19,A)") '<<<<<<<<< SUBSPACE ROTATION <<<<<<<<<<'
886 WRITE (output_unit, "(T4,A,T13,A,T21,A,T38,A,T51,A,T65,A/,T4,A)") &
887 "In-step", "Time", "Convergence", "Band ene.", "Total ene.", "Energy diff.", repeat("-", 74)
888 END IF
889
890 ! recalculate density matrix here
891
892 ! update of density
893 CALL qs_rho_update_rho(rho, qs_env=qs_env)
894
895 DO iloop = 1, subspace_env%max_iter
896 t1 = m_walltime()
897 converged = .false.
898 ene_old = energy%total
899
900 CALL qs_ks_did_change(ks_env, rho_changed=.true.)
901 CALL qs_ks_update_qs_env(qs_env, calculate_forces=.false., &
902 just_energy=.false., print_active=.false.)
903
904 max_val = 0.0_dp
905 sum_val = 0.0_dp
906 sum_band = 0.0_dp
907 DO ispin = 1, SIZE(matrix_ks)
908 CALL get_mo_set(mo_set=mos(ispin), &
909 nao=nao, &
910 nmo=nmo, &
911 eigenvalues=mo_eigenvalues, &
912 occupation_numbers=mo_occupations, &
913 mo_coeff=mo_coeff)
914
915 !compute C'HC
916 chc => subspace_env%chc_mat(ispin)
917 evec => subspace_env%c_vec(ispin)
918 c0 => subspace_env%c0(ispin)
919 CALL cp_fm_to_fm(mo_coeff, c0)
920 CALL cp_fm_create(work, c0%matrix_struct)
921 CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, c0, work, nmo)
922 CALL parallel_gemm('T', 'N', nmo, nmo, nao, rone, c0, work, rzero, chc)
923 CALL cp_fm_release(work)
924 !diagonalize C'HC
925 CALL choose_eigv_solver(chc, evec, mo_eigenvalues)
926
927 !rotate the mos by the eigenvectors of C'HC
928 CALL parallel_gemm('N', 'N', nao, nmo, nmo, rone, c0, evec, rzero, mo_coeff)
929
930 CALL set_mo_occupation(mo_set=mos(ispin), &
931 smear=scf_control%smear)
932
933 ! does not yet handle k-points
934 CALL calculate_density_matrix(mos(ispin), &
935 subspace_env%p_matrix_mix(ispin, 1)%matrix)
936
937 DO i = 1, nmo
938 sum_band = sum_band + mo_eigenvalues(i)*mo_occupations(i)
939 END DO
940
941 !check for self consistency
942 END DO
943
944 IF (subspace_env%mixing_method == direct_mixing_nr) THEN
945 CALL scf_env_density_mixing(subspace_env%p_matrix_mix, &
946 scf_env%mixing_store, rho_ao_kp, para_env, iter_delta, iloop)
947 ELSE
948 CALL self_consistency_check(rho_ao_kp, scf_env%p_delta, para_env, &
949 subspace_env%p_matrix_mix, delta=iter_delta)
950 END IF
951
952 DO ispin = 1, nspin
953 ! does not yet handle k-points
954 CALL dbcsr_copy(rho_ao(ispin)%matrix, subspace_env%p_matrix_mix(ispin, 1)%matrix)
955 END DO
956 ! update of density
957 CALL qs_rho_update_rho(rho, qs_env=qs_env)
958 ! Mixing in reciprocal space
959 IF (subspace_env%mixing_method >= gspace_mixing_nr) THEN
960 CALL gspace_mixing(qs_env, scf_env%mixing_method, subspace_env%mixing_store, &
961 rho, para_env, scf_env%iter_count)
962 END IF
963
964 ene_diff = energy%total - ene_old
965 converged = (abs(ene_diff) < subspace_env%eps_ene .AND. &
966 iter_delta < subspace_env%eps_adapt*scf_env%iter_delta)
967 t2 = m_walltime()
968 IF (output_unit > 0) THEN
969 WRITE (output_unit, "(T4,I5,T11,F8.3,T18,E14.4,T34,F12.5,T46,F16.8,T62,E14.4)") &
970 iloop, t2 - t1, iter_delta, sum_band, energy%total, ene_diff
971 CALL m_flush(output_unit)
972 END IF
973 IF (converged) THEN
974 IF (output_unit > 0) WRITE (output_unit, "(T10,A,I6,A,/)") &
975 " Reached convergence in ", iloop, " iterations "
976 EXIT
977 END IF
978
979 END DO ! iloop
980
981 NULLIFY (subspace_env%p_matrix_mix)
982 DO ispin = 1, nspin
983 ! does not yet handle k-points
984 CALL dbcsr_copy(scf_env%p_mix_new(ispin, 1)%matrix, rho_ao(ispin)%matrix)
985 CALL dbcsr_copy(rho_ao(ispin)%matrix, subspace_env%p_matrix_store(ispin)%matrix)
986
987 DEALLOCATE (eval_first(ispin)%array, occ_first(ispin)%array)
988 END DO
989 DEALLOCATE (eval_first, occ_first)
990
991 CALL timestop(handle)
992
993 END SUBROUTINE do_scf_diag_subspace
994
995! **************************************************************************************************
996!> \brief ...
997!> \param subspace_env ...
998!> \param qs_env ...
999!> \param mos ...
1000! **************************************************************************************************
1001 SUBROUTINE diag_subspace_allocate(subspace_env, qs_env, mos)
1002
1003 TYPE(subspace_env_type), POINTER :: subspace_env
1004 TYPE(qs_environment_type), POINTER :: qs_env
1005 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1006
1007 CHARACTER(LEN=*), PARAMETER :: routinen = 'diag_subspace_allocate'
1008
1009 INTEGER :: handle, i, ispin, nmo, nspin
1010 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
1011 TYPE(cp_fm_type), POINTER :: mo_coeff
1012 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1013 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1014 POINTER :: sab_orb
1015
1016 CALL timeset(routinen, handle)
1017
1018 NULLIFY (sab_orb, matrix_s)
1019 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb, &
1020 matrix_s=matrix_s)
1021
1022 nspin = SIZE(mos)
1023! *** allocate p_atrix_store ***
1024 IF (.NOT. ASSOCIATED(subspace_env%p_matrix_store)) THEN
1025 CALL dbcsr_allocate_matrix_set(subspace_env%p_matrix_store, nspin)
1026
1027 DO i = 1, nspin
1028 ALLOCATE (subspace_env%p_matrix_store(i)%matrix)
1029 CALL dbcsr_create(matrix=subspace_env%p_matrix_store(i)%matrix, template=matrix_s(1)%matrix, &
1030 name="DENSITY_STORE", matrix_type=dbcsr_type_symmetric)
1031 CALL cp_dbcsr_alloc_block_from_nbl(subspace_env%p_matrix_store(i)%matrix, &
1032 sab_orb)
1033 CALL dbcsr_set(subspace_env%p_matrix_store(i)%matrix, 0.0_dp)
1034 END DO
1035
1036 END IF
1037
1038 ALLOCATE (subspace_env%chc_mat(nspin))
1039 ALLOCATE (subspace_env%c_vec(nspin))
1040 ALLOCATE (subspace_env%c0(nspin))
1041
1042 DO ispin = 1, nspin
1043 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1044 CALL cp_fm_create(subspace_env%c0(ispin), mo_coeff%matrix_struct)
1045 NULLIFY (fm_struct_tmp)
1046 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo, ncol_global=nmo, &
1047 para_env=mo_coeff%matrix_struct%para_env, &
1048 context=mo_coeff%matrix_struct%context)
1049 CALL cp_fm_create(subspace_env%chc_mat(ispin), fm_struct_tmp, "chc")
1050 CALL cp_fm_create(subspace_env%c_vec(ispin), fm_struct_tmp, "vec")
1051 CALL cp_fm_struct_release(fm_struct_tmp)
1052 END DO
1053
1054 CALL timestop(handle)
1055
1056 END SUBROUTINE diag_subspace_allocate
1057
1058! **************************************************************************************************
1059!> \brief the inner loop of scf, specific to diagonalization without S matrix
1060!> basically, in goes the ks matrix out goes a new p matrix
1061!> \param scf_env ...
1062!> \param mos ...
1063!> \param matrix_ks ...
1064!> \param scf_control ...
1065!> \param scf_section ...
1066!> \param diis_step ...
1067!> \par History
1068!> 03.2006 created [Joost VandeVondele]
1069! **************************************************************************************************
1070 SUBROUTINE do_special_diag(scf_env, mos, matrix_ks, scf_control, &
1071 scf_section, diis_step)
1072
1073 TYPE(qs_scf_env_type), POINTER :: scf_env
1074 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
1075 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1076 TYPE(scf_control_type), POINTER :: scf_control
1077 TYPE(section_vals_type), POINTER :: scf_section
1078 LOGICAL, INTENT(INOUT) :: diis_step
1079
1080 INTEGER :: ispin, nspin
1081 LOGICAL :: do_level_shift, use_jacobi
1082 REAL(kind=dp) :: diis_error
1083
1084 nspin = SIZE(matrix_ks)
1085
1086 DO ispin = 1, nspin
1087 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, scf_env%scf_work1(ispin))
1088 END DO
1089 IF (scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis) THEN
1090 CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
1091 scf_env%scf_work2, scf_env%iter_delta, diis_error, diis_step, &
1092 scf_control%eps_diis, scf_control%nmixing, &
1093 scf_section=scf_section)
1094 ELSE
1095 diis_step = .false.
1096 END IF
1097
1098 IF ((scf_env%iter_count > 1) .AND. (scf_env%iter_delta < scf_control%diagonalization%eps_jacobi)) THEN
1099 use_jacobi = .true.
1100 ELSE
1101 use_jacobi = .false.
1102 END IF
1103
1104 do_level_shift = ((scf_control%level_shift /= 0.0_dp) .AND. &
1105 ((scf_control%density_guess == core_guess) .OR. (scf_env%iter_count > 1)))
1106 IF (diis_step) THEN
1107 scf_env%iter_param = diis_error
1108 IF (use_jacobi) THEN
1109 scf_env%iter_method = "DIIS/Jacobi"
1110 ELSE
1111 scf_env%iter_method = "DIIS/Diag."
1112 END IF
1113 ELSE
1114 IF (scf_env%mixing_method == 1) THEN
1115 scf_env%iter_param = scf_env%p_mix_alpha
1116 IF (use_jacobi) THEN
1117 scf_env%iter_method = "P_Mix/Jacobi"
1118 ELSE
1119 scf_env%iter_method = "P_Mix/Diag."
1120 END IF
1121 ELSEIF (scf_env%mixing_method > 1) THEN
1122 scf_env%iter_param = scf_env%mixing_store%alpha
1123 IF (use_jacobi) THEN
1124 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Jacobi"
1125 ELSE
1126 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Diag."
1127 END IF
1128 END IF
1129 END IF
1130 scf_env%iter_delta = 0.0_dp
1131
1132 DO ispin = 1, nspin
1133 CALL eigensolver_simple(matrix_ks=scf_env%scf_work1(ispin), &
1134 mo_set=mos(ispin), &
1135 work=scf_env%scf_work2, &
1136 do_level_shift=do_level_shift, &
1137 level_shift=scf_control%level_shift, &
1138 use_jacobi=use_jacobi, &
1139 jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
1140 END DO
1141
1142 CALL set_mo_occupation(mo_array=mos, &
1143 smear=scf_control%smear)
1144
1145 DO ispin = 1, nspin
1146 ! does not yet handle k-points
1147 CALL calculate_density_matrix(mos(ispin), &
1148 scf_env%p_mix_new(ispin, 1)%matrix)
1149 END DO
1150
1151 END SUBROUTINE do_special_diag
1152
1153! **************************************************************************************************
1154!> \brief the inner loop of scf, specific to iterative diagonalization using OT
1155!> with S matrix; basically, in goes the ks matrix out goes a new p matrix
1156!> \param scf_env ...
1157!> \param mos ...
1158!> \param matrix_ks ...
1159!> \param matrix_s ...
1160!> \param scf_control ...
1161!> \param scf_section ...
1162!> \param diis_step ...
1163!> \par History
1164!> 10.2008 created [JGH]
1165! **************************************************************************************************
1166 SUBROUTINE do_ot_diag(scf_env, mos, matrix_ks, matrix_s, &
1167 scf_control, scf_section, diis_step)
1168
1169 TYPE(qs_scf_env_type), POINTER :: scf_env
1170 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
1171 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
1172 TYPE(scf_control_type), POINTER :: scf_control
1173 TYPE(section_vals_type), POINTER :: scf_section
1174 LOGICAL, INTENT(INOUT) :: diis_step
1175
1176 INTEGER :: homo, ispin, nmo, nspin
1177 REAL(kind=dp) :: diis_error, eps_iter
1178 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
1179 TYPE(cp_fm_type), POINTER :: mo_coeff
1180
1181 NULLIFY (eigenvalues)
1182
1183 nspin = SIZE(matrix_ks)
1184
1185 DO ispin = 1, nspin
1186 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, &
1187 scf_env%scf_work1(ispin))
1188 END DO
1189
1190 IF ((scf_env%iter_count > 1) .AND. (.NOT. scf_env%skip_diis)) THEN
1191 CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
1192 scf_env%scf_work2, scf_env%iter_delta, diis_error, diis_step, &
1193 scf_control%eps_diis, scf_control%nmixing, &
1194 s_matrix=matrix_s, &
1195 scf_section=scf_section)
1196 ELSE
1197 diis_step = .false.
1198 END IF
1199
1200 eps_iter = scf_control%diagonalization%eps_iter
1201 IF (diis_step) THEN
1202 scf_env%iter_param = diis_error
1203 scf_env%iter_method = "DIIS/OTdiag"
1204 DO ispin = 1, nspin
1205 CALL copy_fm_to_dbcsr(scf_env%scf_work1(ispin), &
1206 matrix_ks(ispin)%matrix, keep_sparsity=.true.)
1207 END DO
1208 eps_iter = max(eps_iter, scf_control%diagonalization%eps_adapt*diis_error)
1209 ELSE
1210 IF (scf_env%mixing_method == 1) THEN
1211 scf_env%iter_param = scf_env%p_mix_alpha
1212 scf_env%iter_method = "P_Mix/OTdiag."
1213 ELSEIF (scf_env%mixing_method > 1) THEN
1214 scf_env%iter_param = scf_env%mixing_store%alpha
1215 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/OTdiag."
1216 END IF
1217 END IF
1218
1219 scf_env%iter_delta = 0.0_dp
1220
1221 DO ispin = 1, nspin
1222 CALL get_mo_set(mos(ispin), &
1223 mo_coeff=mo_coeff, &
1224 eigenvalues=eigenvalues, &
1225 nmo=nmo, &
1226 homo=homo)
1227 CALL ot_eigensolver(matrix_h=matrix_ks(ispin)%matrix, &
1228 matrix_s=matrix_s(1)%matrix, &
1229 matrix_c_fm=mo_coeff, &
1230 preconditioner=scf_env%ot_preconditioner(1)%preconditioner, &
1231 eps_gradient=eps_iter, &
1232 iter_max=scf_control%diagonalization%max_iter, &
1233 silent=.true., &
1234 ot_settings=scf_control%diagonalization%ot_settings)
1235 CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin)%matrix, &
1236 evals_arg=eigenvalues, &
1237 do_rotation=.true.)
1238 CALL copy_fm_to_dbcsr(mos(ispin)%mo_coeff, &
1239 mos(ispin)%mo_coeff_b)
1240 !fm->dbcsr
1241 END DO
1242
1243 CALL set_mo_occupation(mo_array=mos, &
1244 smear=scf_control%smear)
1245
1246 DO ispin = 1, nspin
1247 ! does not yet handle k-points
1248 CALL calculate_density_matrix(mos(ispin), &
1249 scf_env%p_mix_new(ispin, 1)%matrix)
1250 END DO
1251
1252 END SUBROUTINE do_ot_diag
1253
1254! **************************************************************************************************
1255!> \brief Solve a set restricted open Kohn-Sham (ROKS) equations based on the
1256!> alpha and beta Kohn-Sham matrices from unrestricted Kohn-Sham.
1257!> \param scf_env ...
1258!> \param mos ...
1259!> \param matrix_ks ...
1260!> \param matrix_s ...
1261!> \param scf_control ...
1262!> \param scf_section ...
1263!> \param diis_step ...
1264!> \param orthogonal_basis ...
1265!> \par History
1266!> 04.2006 created [MK]
1267!> Revised (01.05.06,MK)
1268!> \note
1269!> this is only a high-spin ROKS.
1270! **************************************************************************************************
1271 SUBROUTINE do_roks_diag(scf_env, mos, matrix_ks, matrix_s, &
1272 scf_control, scf_section, diis_step, &
1273 orthogonal_basis)
1274
1275 ! Literature: - C. C. J. Roothaan, Rev. Mod. Phys. 32, 179 (1960)
1276 ! - M. F. Guest and V. R. Saunders, Mol. Phys. 28(3), 819 (1974)
1277 ! - M. Filatov and S. Shaik, Chem. Phys. Lett. 288, 689 (1998)
1278
1279 TYPE(qs_scf_env_type), POINTER :: scf_env
1280 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1281 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
1282 TYPE(scf_control_type), POINTER :: scf_control
1283 TYPE(section_vals_type), POINTER :: scf_section
1284 LOGICAL, INTENT(INOUT) :: diis_step
1285 LOGICAL, INTENT(IN) :: orthogonal_basis
1286
1287 CHARACTER(LEN=*), PARAMETER :: routinen = 'do_roks_diag'
1288
1289 INTEGER :: handle, homoa, homob, imo, nalpha, nao, &
1290 nbeta, nmo
1291 REAL(kind=dp) :: diis_error, level_shift_loc
1292 REAL(kind=dp), DIMENSION(:), POINTER :: eiga, eigb, occa, occb
1293 TYPE(cp_fm_type), POINTER :: ksa, ksb, mo2ao, moa, mob, ortho, work
1294
1295! -------------------------------------------------------------------------
1296
1297 CALL timeset(routinen, handle)
1298
1299 IF (scf_env%cholesky_method == cholesky_inverse) THEN
1300 ortho => scf_env%ortho_m1
1301 ELSE
1302 ortho => scf_env%ortho
1303 END IF
1304 work => scf_env%scf_work2
1305
1306 ksa => scf_env%scf_work1(1)
1307 ksb => scf_env%scf_work1(2)
1308
1309 CALL copy_dbcsr_to_fm(matrix_ks(1)%matrix, ksa)
1310 CALL copy_dbcsr_to_fm(matrix_ks(2)%matrix, ksb)
1311
1312 ! Get MO information
1313
1314 CALL get_mo_set(mo_set=mos(1), &
1315 nao=nao, &
1316 nmo=nmo, &
1317 nelectron=nalpha, &
1318 homo=homoa, &
1319 eigenvalues=eiga, &
1320 occupation_numbers=occa, &
1321 mo_coeff=moa)
1322
1323 CALL get_mo_set(mo_set=mos(2), &
1324 nelectron=nbeta, &
1325 homo=homob, &
1326 eigenvalues=eigb, &
1327 occupation_numbers=occb, &
1328 mo_coeff=mob)
1329
1330 ! Define the amount of level-shifting
1331
1332 IF ((scf_control%level_shift /= 0.0_dp) .AND. &
1333 ((scf_control%density_guess == core_guess) .OR. &
1334 (scf_control%density_guess == restart_guess) .OR. &
1335 (scf_env%iter_count > 1))) THEN
1336 level_shift_loc = scf_control%level_shift
1337 ELSE
1338 level_shift_loc = 0.0_dp
1339 END IF
1340
1341 IF ((scf_env%iter_count > 1) .OR. &
1342 (scf_control%density_guess == core_guess) .OR. &
1343 (scf_control%density_guess == restart_guess)) THEN
1344
1345 ! Transform the spin unrestricted alpha and beta Kohn-Sham matrices
1346 ! from AO basis to MO basis: K(MO) = C(T)*K(AO)*C
1347
1348 CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, moa, 0.0_dp, work)
1349 CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, moa, work, 0.0_dp, ksa)
1350
1351 CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksb, moa, 0.0_dp, work)
1352 CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, moa, work, 0.0_dp, ksb)
1353
1354 ! Combine the spin unrestricted alpha and beta Kohn-Sham matrices
1355 ! in the MO basis
1356
1357 IF (scf_control%roks_scheme == general_roks) THEN
1358 CALL combine_ks_matrices(ksa, ksb, occa, occb, scf_control%roks_f, &
1359 nalpha, nbeta)
1360 ELSE IF (scf_control%roks_scheme == high_spin_roks) THEN
1361 CALL combine_ks_matrices(ksa, ksb, occa, occb, scf_control%roks_parameter)
1362 ELSE
1363 cpabort("Unknown ROKS scheme requested")
1364 END IF
1365
1366 ! Back-transform the restricted open Kohn-Sham matrix from MO basis
1367 ! to AO basis
1368
1369 IF (orthogonal_basis) THEN
1370 ! Q = C
1371 mo2ao => moa
1372 ELSE
1373 ! Q = S*C
1374 mo2ao => mob
1375!MK CALL copy_sm_to_fm(matrix_s(1)%matrix,work)
1376!MK CALL cp_fm_symm("L", "U",nao, nao, 1.0_dp, work, moa, 0.0_dp, mo2ao)
1377 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, moa, mo2ao, nao)
1378 END IF
1379
1380 ! K(AO) = Q*K(MO)*Q(T)
1381
1382 CALL parallel_gemm("N", "T", nao, nao, nao, 1.0_dp, ksa, mo2ao, 0.0_dp, work)
1383 CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, mo2ao, work, 0.0_dp, ksa)
1384
1385 ELSE
1386
1387 ! No transformation matrix available, yet. The closed shell part,
1388 ! i.e. the beta Kohn-Sham matrix in AO basis, is taken.
1389 ! There might be better choices, anyhow.
1390
1391 CALL cp_fm_to_fm(ksb, ksa)
1392
1393 END IF
1394
1395 ! Update DIIS buffer and possibly perform DIIS extrapolation step
1396
1397 IF (scf_env%iter_count > 1) THEN
1398 IF (orthogonal_basis) THEN
1399 CALL qs_diis_b_step(diis_buffer=scf_env%scf_diis_buffer, &
1400 mo_array=mos, &
1401 kc=scf_env%scf_work1, &
1402 sc=work, &
1403 delta=scf_env%iter_delta, &
1404 error_max=diis_error, &
1405 diis_step=diis_step, &
1406 eps_diis=scf_control%eps_diis, &
1407 scf_section=scf_section, &
1408 roks=.true.)
1409 cpassert(scf_env%iter_delta == scf_env%iter_delta)
1410 ELSE
1411 CALL qs_diis_b_step(diis_buffer=scf_env%scf_diis_buffer, &
1412 mo_array=mos, &
1413 kc=scf_env%scf_work1, &
1414 sc=work, &
1415 delta=scf_env%iter_delta, &
1416 error_max=diis_error, &
1417 diis_step=diis_step, &
1418 eps_diis=scf_control%eps_diis, &
1419 scf_section=scf_section, &
1420 s_matrix=matrix_s, &
1421 roks=.true.)
1422 END IF
1423 END IF
1424
1425 IF (diis_step) THEN
1426 scf_env%iter_param = diis_error
1427 scf_env%iter_method = "DIIS/Diag."
1428 ELSE
1429 IF (scf_env%mixing_method == 1) THEN
1430 scf_env%iter_param = scf_env%p_mix_alpha
1431 scf_env%iter_method = "P_Mix/Diag."
1432 ELSEIF (scf_env%mixing_method > 1) THEN
1433 scf_env%iter_param = scf_env%mixing_store%alpha
1434 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Diag."
1435 END IF
1436 END IF
1437
1438 scf_env%iter_delta = 0.0_dp
1439
1440 IF (level_shift_loc /= 0.0_dp) THEN
1441
1442 ! Transform the current Kohn-Sham matrix from AO to MO basis
1443 ! for level-shifting using the current MO set
1444
1445 CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, moa, 0.0_dp, work)
1446 CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, moa, work, 0.0_dp, ksa)
1447
1448 ! Apply level-shifting using 50:50 split of the shift (could be relaxed)
1449
1450 DO imo = homob + 1, homoa
1451 CALL cp_fm_add_to_element(ksa, imo, imo, 0.5_dp*level_shift_loc)
1452 END DO
1453 DO imo = homoa + 1, nmo
1454 CALL cp_fm_add_to_element(ksa, imo, imo, level_shift_loc)
1455 END DO
1456
1457 ELSE IF (.NOT. orthogonal_basis) THEN
1458
1459 ! Transform the current Kohn-Sham matrix to an orthogonal basis
1460 SELECT CASE (scf_env%cholesky_method)
1461 CASE (cholesky_reduce)
1462 CALL cp_fm_cholesky_reduce(ksa, ortho)
1463 CASE (cholesky_restore)
1464 CALL cp_fm_uplo_to_full(ksa, work)
1465 CALL cp_fm_cholesky_restore(ksa, nao, ortho, work, &
1466 "SOLVE", pos="RIGHT")
1467 CALL cp_fm_cholesky_restore(work, nao, ortho, ksa, &
1468 "SOLVE", pos="LEFT", transa="T")
1469 CASE (cholesky_inverse)
1470 CALL cp_fm_uplo_to_full(ksa, work)
1471 CALL cp_fm_cholesky_restore(ksa, nao, ortho, work, &
1472 "MULTIPLY", pos="RIGHT")
1473 CALL cp_fm_cholesky_restore(work, nao, ortho, ksa, &
1474 "MULTIPLY", pos="LEFT", transa="T")
1475 CASE (cholesky_off)
1476 CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, ortho, 0.0_dp, work)
1477 CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, ortho, work, 0.0_dp, ksa)
1478 END SELECT
1479
1480 END IF
1481
1482 ! Diagonalization of the ROKS operator matrix
1483
1484 CALL choose_eigv_solver(ksa, work, eiga)
1485
1486 ! Back-transformation of the orthonormal eigenvectors if needed
1487
1488 IF (level_shift_loc /= 0.0_dp) THEN
1489 ! Use old MO set for back-transformation if level-shifting was applied
1490 CALL cp_fm_to_fm(moa, ortho)
1491 CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, ortho, work, 0.0_dp, moa)
1492 ELSE
1493 IF (orthogonal_basis) THEN
1494 CALL cp_fm_to_fm(work, moa)
1495 ELSE
1496 SELECT CASE (scf_env%cholesky_method)
1497 CASE (cholesky_reduce, cholesky_restore)
1498 CALL cp_fm_cholesky_restore(work, nmo, ortho, moa, "SOLVE")
1499 CASE (cholesky_inverse)
1500 CALL cp_fm_cholesky_restore(work, nmo, ortho, moa, "MULTIPLY")
1501 CASE (cholesky_off)
1502 CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, ortho, work, 0.0_dp, moa)
1503 END SELECT
1504 END IF
1505 END IF
1506
1507 ! Correct MO eigenvalues, if level-shifting was applied
1508
1509 IF (level_shift_loc /= 0.0_dp) THEN
1510 DO imo = homob + 1, homoa
1511 eiga(imo) = eiga(imo) - 0.5_dp*level_shift_loc
1512 END DO
1513 DO imo = homoa + 1, nmo
1514 eiga(imo) = eiga(imo) - level_shift_loc
1515 END DO
1516 END IF
1517
1518 ! Update also the beta MO set
1519
1520 eigb(:) = eiga(:)
1521 CALL cp_fm_to_fm(moa, mob)
1522
1523 ! Calculate the new alpha and beta density matrix
1524
1525 ! does not yet handle k-points
1526 CALL calculate_density_matrix(mos(1), scf_env%p_mix_new(1, 1)%matrix)
1527 CALL calculate_density_matrix(mos(2), scf_env%p_mix_new(2, 1)%matrix)
1528
1529 CALL timestop(handle)
1530
1531 END SUBROUTINE do_roks_diag
1532
1533! **************************************************************************************************
1534!> \brief iterative diagonalization using the block Krylov-space approach
1535!> \param scf_env ...
1536!> \param mos ...
1537!> \param matrix_ks ...
1538!> \param scf_control ...
1539!> \param scf_section ...
1540!> \param check_moconv_only ...
1541!> \param
1542!> \par History
1543!> 05.2009 created [MI]
1544! **************************************************************************************************
1545
1546 SUBROUTINE do_block_krylov_diag(scf_env, mos, matrix_ks, &
1547 scf_control, scf_section, check_moconv_only)
1548
1549 TYPE(qs_scf_env_type), POINTER :: scf_env
1550 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
1551 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1552 TYPE(scf_control_type), POINTER :: scf_control
1553 TYPE(section_vals_type), POINTER :: scf_section
1554 LOGICAL, INTENT(IN), OPTIONAL :: check_moconv_only
1555
1556 CHARACTER(LEN=*), PARAMETER :: routinen = 'do_block_krylov_diag'
1557 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1558
1559 INTEGER :: handle, homo, ispin, iter, nao, nmo, &
1560 output_unit
1561 LOGICAL :: converged, my_check_moconv_only
1562 REAL(dp) :: eps_iter, t1, t2
1563 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues
1564 TYPE(cp_fm_type), POINTER :: c0, c1, chc, evec, ks, mo_coeff, ortho, &
1565 work
1566 TYPE(cp_logger_type), POINTER :: logger
1567
1568 logger => cp_get_default_logger()
1569 CALL timeset(routinen, handle)
1570
1571 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%LANCZOS", &
1572 extension=".scfLog")
1573
1574 my_check_moconv_only = .false.
1575 IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
1576
1577 NULLIFY (mo_coeff, ortho, work, ks)
1578 NULLIFY (mo_eigenvalues)
1579 NULLIFY (c0, c1)
1580
1581 IF (scf_env%cholesky_method == cholesky_inverse) THEN
1582 ortho => scf_env%ortho_m1
1583 ELSE
1584 ortho => scf_env%ortho
1585 END IF
1586 work => scf_env%scf_work2
1587
1588 DO ispin = 1, SIZE(matrix_ks)
1589 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, &
1590 scf_env%scf_work1(ispin))
1591 END DO
1592
1593 IF (scf_env%mixing_method == 1) THEN
1594 scf_env%iter_param = scf_env%p_mix_alpha
1595 scf_env%iter_method = "P_Mix/Lanczos"
1596 ELSE
1597! scf_env%iter_param = scf_env%mixing_store%alpha
1598 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Lanc."
1599 END IF
1600
1601 DO ispin = 1, SIZE(matrix_ks)
1602
1603 ks => scf_env%scf_work1(ispin)
1604 CALL cp_fm_uplo_to_full(ks, work)
1605
1606 CALL get_mo_set(mo_set=mos(ispin), &
1607 nao=nao, &
1608 nmo=nmo, &
1609 homo=homo, &
1610 eigenvalues=mo_eigenvalues, &
1611 mo_coeff=mo_coeff)
1612
1613 NULLIFY (c0, c1)
1614 c0 => scf_env%krylov_space%mo_conv(ispin)
1615 c1 => scf_env%krylov_space%mo_refine(ispin)
1616 SELECT CASE (scf_env%cholesky_method)
1617 CASE (cholesky_reduce)
1618 CALL cp_fm_cholesky_reduce(ks, ortho)
1619 CALL cp_fm_uplo_to_full(ks, work)
1620 CALL cp_fm_cholesky_restore(mo_coeff, nmo, ortho, c0, "MULTIPLY")
1621 CASE (cholesky_restore)
1622 CALL cp_fm_cholesky_restore(ks, nao, ortho, work, &
1623 "SOLVE", pos="RIGHT")
1624 CALL cp_fm_cholesky_restore(work, nao, ortho, ks, &
1625 "SOLVE", pos="LEFT", transa="T")
1626 CALL cp_fm_cholesky_restore(mo_coeff, nmo, ortho, c0, "MULTIPLY")
1627 CASE (cholesky_inverse)
1628 CALL cp_fm_cholesky_restore(ks, nao, ortho, work, &
1629 "MULTIPLY", pos="RIGHT")
1630 CALL cp_fm_cholesky_restore(work, nao, ortho, ks, &
1631 "MULTIPLY", pos="LEFT", transa="T")
1632 CALL cp_fm_cholesky_restore(mo_coeff, nmo, ortho, c0, "SOLVE")
1633 END SELECT
1634
1635 scf_env%krylov_space%nmo_nc = nmo
1636 scf_env%krylov_space%nmo_conv = 0
1637
1638 t1 = m_walltime()
1639 IF (output_unit > 0) THEN
1640 WRITE (output_unit, "(/T15,A)") '<<<<<<<<< LANCZOS REFINEMENT <<<<<<<<<<'
1641 WRITE (output_unit, "(T8,A,T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") &
1642 " Spin ", " Cycle ", &
1643 " conv. MOS ", " B2MAX ", " B2MIN ", " Time", repeat("-", 60)
1644 END IF
1645 eps_iter = max(scf_env%krylov_space%eps_conv, scf_env%krylov_space%eps_adapt*scf_env%iter_delta)
1646 iter = 0
1647 converged = .false.
1648 !Check convergence of MOS
1649 IF (my_check_moconv_only) THEN
1650
1651 CALL lanczos_refinement(scf_env%krylov_space, ks, c0, c1, mo_eigenvalues, &
1652 nao, eps_iter, ispin, check_moconv_only=my_check_moconv_only)
1653 t2 = m_walltime()
1654 IF (output_unit > 0) &
1655 WRITE (output_unit, '(T8,I3,T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
1656 ispin, iter, scf_env%krylov_space%nmo_conv, &
1657 scf_env%krylov_space%max_res_norm, scf_env%krylov_space%min_res_norm, t2 - t1
1658
1659 cycle
1660 ELSE
1661 !Block Lanczos refinement
1662 DO iter = 1, scf_env%krylov_space%max_iter
1663 CALL lanczos_refinement_2v(scf_env%krylov_space, ks, c0, c1, mo_eigenvalues, &
1664 nao, eps_iter, ispin)
1665 t2 = m_walltime()
1666 IF (output_unit > 0) THEN
1667 WRITE (output_unit, '(T8,I3,T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
1668 ispin, iter, scf_env%krylov_space%nmo_conv, &
1669 scf_env%krylov_space%max_res_norm, scf_env%krylov_space%min_res_norm, t2 - t1
1670 END IF
1671 t1 = m_walltime()
1672 IF (scf_env%krylov_space%max_res_norm < eps_iter) THEN
1673 converged = .true.
1674 IF (output_unit > 0) WRITE (output_unit, *) &
1675 " Reached convergence in ", iter, " iterations "
1676 EXIT
1677 END IF
1678 END DO
1679
1680 IF (.NOT. converged .AND. output_unit > 0) THEN
1681 WRITE (output_unit, "(T4, A)") " WARNING Lanczos refinement could "// &
1682 "not converge all the mos:"
1683 WRITE (output_unit, "(T40,A,T70,I10)") " number of not converged mos ", &
1684 scf_env%krylov_space%nmo_nc
1685 WRITE (output_unit, "(T40,A,T70,E10.2)") " max norm of the residual ", &
1686 scf_env%krylov_space%max_res_norm
1687
1688 END IF
1689
1690 ! For the moment skip the re-orthogonalization
1691 IF (.false.) THEN
1692 !Re-orthogonalization
1693 NULLIFY (chc, evec)
1694 chc => scf_env%krylov_space%chc_mat(ispin)
1695 evec => scf_env%krylov_space%c_vec(ispin)
1696 CALL parallel_gemm('N', 'N', nao, nmo, nao, rone, ks, c0, rzero, work)
1697 CALL parallel_gemm('T', 'N', nmo, nmo, nao, rone, c0, work, rzero, chc)
1698 !Diagonalize (C^t)HC
1699 CALL choose_eigv_solver(chc, evec, mo_eigenvalues)
1700 !Rotate the C vectors
1701 CALL parallel_gemm('N', 'N', nao, nmo, nmo, rone, c0, evec, rzero, c1)
1702 c0 => scf_env%krylov_space%mo_refine(ispin)
1703 END IF
1704
1705 IF (scf_env%cholesky_method == cholesky_inverse) THEN
1706 CALL cp_fm_cholesky_restore(c0, nmo, ortho, mo_coeff, "MULTIPLY")
1707 ELSE
1708 CALL cp_fm_cholesky_restore(c0, nmo, ortho, mo_coeff, "SOLVE")
1709 END IF
1710
1711 CALL set_mo_occupation(mo_set=mos(ispin), &
1712 smear=scf_control%smear)
1713
1714 ! does not yet handle k-points
1715 CALL calculate_density_matrix(mos(ispin), &
1716 scf_env%p_mix_new(ispin, 1)%matrix)
1717 END IF
1718 END DO ! ispin
1719
1720 IF (output_unit > 0) THEN
1721 WRITE (output_unit, "(T15,A/)") '<<<<<<<<< END LANCZOS REFINEMENT <<<<<<<<<<'
1722 END IF
1723
1724 CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
1725 "PRINT%LANCZOS")
1726
1727 CALL timestop(handle)
1728
1729 END SUBROUTINE do_block_krylov_diag
1730
1731! **************************************************************************************************
1732!> \brief iterative diagonalization using the block davidson space approach
1733!> \param qs_env ...
1734!> \param scf_env ...
1735!> \param mos ...
1736!> \param matrix_ks ...
1737!> \param matrix_s ...
1738!> \param scf_control ...
1739!> \param scf_section ...
1740!> \param check_moconv_only ...
1741!> \param
1742!> \par History
1743!> 05.2011 created [MI]
1744! **************************************************************************************************
1745
1746 SUBROUTINE do_block_davidson_diag(qs_env, scf_env, mos, matrix_ks, matrix_s, &
1747 scf_control, scf_section, check_moconv_only)
1748
1749 TYPE(qs_environment_type), POINTER :: qs_env
1750 TYPE(qs_scf_env_type), POINTER :: scf_env
1751 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
1752 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
1753 TYPE(scf_control_type), POINTER :: scf_control
1754 TYPE(section_vals_type), POINTER :: scf_section
1755 LOGICAL, INTENT(IN), OPTIONAL :: check_moconv_only
1756
1757 CHARACTER(LEN=*), PARAMETER :: routinen = 'do_block_davidson_diag'
1758
1759 INTEGER :: handle, ispin, nspins, output_unit
1760 LOGICAL :: do_prec, my_check_moconv_only
1761 TYPE(cp_logger_type), POINTER :: logger
1762
1763 logger => cp_get_default_logger()
1764 CALL timeset(routinen, handle)
1765
1766 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DAVIDSON", &
1767 extension=".scfLog")
1768
1769 IF (output_unit > 0) &
1770 WRITE (output_unit, "(/T15,A)") '<<<<<<<<< DAVIDSON ITERATIONS <<<<<<<<<<'
1771
1772 IF (scf_env%mixing_method == 1) THEN
1773 scf_env%iter_param = scf_env%p_mix_alpha
1774 scf_env%iter_method = "P_Mix/Dav."
1775 ELSE
1776 scf_env%iter_param = scf_env%mixing_store%alpha
1777 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Dav."
1778 END IF
1779
1780 my_check_moconv_only = .false.
1781 IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
1782 do_prec = .false.
1783 IF (scf_env%block_davidson_env(1)%prec_type /= 0 .AND. &
1784 scf_env%iter_count >= scf_env%block_davidson_env(1)%first_prec) THEN
1785 do_prec = .true.
1786 END IF
1787
1788 nspins = SIZE(matrix_ks)
1789
1790 IF (do_prec .AND. (scf_env%iter_count == scf_env%block_davidson_env(1)%first_prec .OR. &
1791 modulo(scf_env%iter_count, scf_env%block_davidson_env(1)%niter_new_prec) == 0)) THEN
1792 CALL restart_preconditioner(qs_env, scf_env%ot_preconditioner, &
1793 prec_type=scf_env%block_davidson_env(1)%prec_type, nspins=nspins)
1794 CALL prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, scf_env%ot_preconditioner, &
1795 scf_env%block_davidson_env(1)%prec_type, &
1796 scf_env%block_davidson_env(1)%solver_type, &
1797 scf_env%block_davidson_env(1)%energy_gap, nspins, &
1798 convert_to_dbcsr=scf_env%block_davidson_env(1)%use_sparse_mos, &
1799 full_mo_set=.true.)
1800 END IF
1801
1802 DO ispin = 1, nspins
1803 IF (scf_env%block_davidson_env(ispin)%use_sparse_mos) THEN
1804 IF (.NOT. do_prec) THEN
1805 CALL generate_extended_space_sparse(scf_env%block_davidson_env(ispin), mos(ispin), &
1806 matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit)
1807 ELSE
1808 CALL generate_extended_space_sparse(scf_env%block_davidson_env(ispin), mos(ispin), &
1809 matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit, &
1810 scf_env%ot_preconditioner(ispin)%preconditioner)
1811 END IF
1812
1813 ELSE
1814 IF (.NOT. do_prec) THEN
1815 CALL generate_extended_space(scf_env%block_davidson_env(ispin), mos(ispin), &
1816 matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit)
1817 ELSE
1818 CALL generate_extended_space(scf_env%block_davidson_env(ispin), mos(ispin), &
1819 matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit, &
1820 scf_env%ot_preconditioner(ispin)%preconditioner)
1821 END IF
1822 END IF
1823 END DO !ispin
1824
1825 CALL set_mo_occupation(mo_array=mos, &
1826 smear=scf_control%smear)
1827
1828 DO ispin = 1, nspins
1829 ! does not yet handle k-points
1830 CALL calculate_density_matrix(mos(ispin), &
1831 scf_env%p_mix_new(ispin, 1)%matrix)
1832 END DO
1833
1834 IF (output_unit > 0) THEN
1835 WRITE (output_unit, "(T15,A/)") '<<<<<<<<< END DAVIDSON ITERATION <<<<<<<<<<'
1836 END IF
1837
1838 CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
1839 "PRINT%DAVIDSON")
1840
1841 CALL timestop(handle)
1842
1843 END SUBROUTINE do_block_davidson_diag
1844
1845! **************************************************************************************************
1846!> \brief Kpoint diagonalization routine
1847!> Transforms matrices to kpoint, distributes kpoint groups, performs diagonalization
1848!> \param matrix_s Overlap matrices (RS indices, global)
1849!> \param kpoints Kpoint environment
1850!> \param fmwork full matrices distributed over all groups
1851!> \par History
1852!> 02.2026 created [JGH]
1853! **************************************************************************************************
1854 SUBROUTINE diag_kp_smat(matrix_s, kpoints, fmwork)
1855
1856 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
1857 TYPE(kpoint_type), POINTER :: kpoints
1858 TYPE(cp_fm_type), DIMENSION(:) :: fmwork
1859
1860 CHARACTER(len=*), PARAMETER :: routinen = 'diag_kp_smat'
1861 COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp, 0.0_dp), &
1862 czero = (0.0_dp, 0.0_dp)
1863
1864 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ceig
1865 INTEGER :: handle, igroup, ik, ikp, indx, kplocal, &
1866 nao, nkp, nkp_groups
1867 INTEGER, DIMENSION(2) :: kp_range
1868 INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1869 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1870 LOGICAL :: my_kpgrp, use_real_wfn
1871 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1872 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
1873 TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
1874 TYPE(cp_cfm_type) :: csmat, cwork
1875 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools
1876 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1877 TYPE(cp_fm_type) :: fmdummy, fmlocal, rsmat
1878 TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix, tmpmat
1879 TYPE(kpoint_env_type), POINTER :: kp
1880 TYPE(mp_para_env_type), POINTER :: para_env
1881 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1882 POINTER :: sab_nl
1883 TYPE(qs_matrix_pools_type), POINTER :: mpools
1884
1885 CALL timeset(routinen, handle)
1886
1887 NULLIFY (sab_nl)
1888 CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
1889 nkp_groups=nkp_groups, kp_dist=kp_dist, sab_nl=sab_nl, &
1890 cell_to_index=cell_to_index)
1891 cpassert(ASSOCIATED(sab_nl))
1892 kplocal = kp_range(2) - kp_range(1) + 1
1893
1894 ! allocate some work matrices
1895 ALLOCATE (rmatrix, cmatrix, tmpmat)
1896 CALL dbcsr_create(rmatrix, template=matrix_s(1, 1)%matrix, &
1897 matrix_type=dbcsr_type_symmetric)
1898 CALL dbcsr_create(cmatrix, template=matrix_s(1, 1)%matrix, &
1899 matrix_type=dbcsr_type_antisymmetric)
1900 CALL dbcsr_create(tmpmat, template=matrix_s(1, 1)%matrix, &
1901 matrix_type=dbcsr_type_no_symmetry)
1902 CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
1903 CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
1904
1905 ! fm pools to be used within a kpoint group
1906 CALL get_kpoint_info(kpoints, mpools=mpools)
1907 CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
1908
1909 CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
1910 CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
1911
1912 IF (use_real_wfn) THEN
1913 CALL cp_fm_create(rsmat, matrix_struct)
1914 ELSE
1915 CALL cp_cfm_create(csmat, matrix_struct)
1916 CALL cp_cfm_create(cwork, matrix_struct)
1917 END IF
1918
1919 CALL cp_fm_get_info(fmwork(1), nrow_global=nao)
1920 ALLOCATE (eigenvalues(nao), ceig(nao))
1921
1922 para_env => kpoints%blacs_env_all%para_env
1923 ALLOCATE (info(kplocal*nkp_groups, 2))
1924
1925 ! Setup and start all the communication
1926 indx = 0
1927 DO ikp = 1, kplocal
1928 DO igroup = 1, nkp_groups
1929 ! number of current kpoint
1930 ik = kp_dist(1, igroup) + ikp - 1
1931 my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
1932 indx = indx + 1
1933 IF (use_real_wfn) THEN
1934 CALL dbcsr_set(rmatrix, 0.0_dp)
1935 CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_s, ispin=1, &
1936 xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1937 CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1938 CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
1939 ELSE
1940 CALL dbcsr_set(rmatrix, 0.0_dp)
1941 CALL dbcsr_set(cmatrix, 0.0_dp)
1942 CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_s, ispin=1, &
1943 xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1944 CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1945 CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
1946 CALL dbcsr_desymmetrize(cmatrix, tmpmat)
1947 CALL copy_dbcsr_to_fm(tmpmat, fmwork(2))
1948 END IF
1949 ! transfer to kpoint group
1950 ! redistribution of matrices, new blacs environment
1951 ! fmwork -> fmlocal -> rsmat/csmat
1952 IF (use_real_wfn) THEN
1953 IF (my_kpgrp) THEN
1954 CALL cp_fm_start_copy_general(fmwork(1), rsmat, para_env, info(indx, 1))
1955 ELSE
1956 CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
1957 END IF
1958 ELSE
1959 IF (my_kpgrp) THEN
1960 CALL cp_fm_start_copy_general(fmwork(1), fmlocal, para_env, info(indx, 1))
1961 CALL cp_fm_start_copy_general(fmwork(2), fmlocal, para_env, info(indx, 2))
1962 ELSE
1963 CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
1964 CALL cp_fm_start_copy_general(fmwork(2), fmdummy, para_env, info(indx, 2))
1965 END IF
1966 END IF
1967 END DO
1968 END DO
1969
1970 ! Finish communication then diagonalise in each group
1971 indx = 0
1972 DO ikp = 1, kplocal
1973 DO igroup = 1, nkp_groups
1974 ! number of current kpoint
1975 ik = kp_dist(1, igroup) + ikp - 1
1976 my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
1977 indx = indx + 1
1978 IF (my_kpgrp) THEN
1979 IF (use_real_wfn) THEN
1980 CALL cp_fm_finish_copy_general(rsmat, info(indx, 1))
1981 ELSE
1982 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
1983 CALL cp_cfm_scale_and_add_fm(z_zero, csmat, z_one, fmlocal)
1984 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
1985 CALL cp_cfm_scale_and_add_fm(z_one, csmat, gaussi, fmlocal)
1986 END IF
1987 END IF
1988 END DO
1989
1990 ! Each kpoint group has now information on a kpoint to be diagonalized
1991 ! Eigensolver Hermite or Symmetric
1992 kp => kpoints%kp_env(ikp)%kpoint_env
1993 IF (use_real_wfn) THEN
1994 CALL choose_eigv_solver(rsmat, fmlocal, eigenvalues)
1995 ELSE
1996 CALL cp_cfm_heevd(csmat, cwork, eigenvalues)
1997 END IF
1998 cpassert(all(eigenvalues(1:nao) >= 0.0_dp))
1999 IF (use_real_wfn) THEN
2000 CALL cp_fm_release(kp%shalf)
2001 CALL cp_fm_create(kp%shalf, matrix_struct)
2002 eigenvalues(1:nao) = sqrt(eigenvalues(1:nao))
2003 CALL cp_fm_to_fm(fmlocal, rsmat)
2004 CALL cp_fm_column_scale(rsmat, eigenvalues)
2005 CALL parallel_gemm("N", "T", nao, nao, nao, 1.0_dp, rsmat, fmlocal, &
2006 0.0_dp, kp%shalf)
2007 ELSE
2008 CALL cp_cfm_release(kp%cshalf)
2009 CALL cp_cfm_create(kp%cshalf, matrix_struct)
2010 ceig(1:nao) = sqrt(eigenvalues(1:nao))
2011 CALL cp_cfm_to_cfm(cwork, csmat)
2012 CALL cp_cfm_column_scale(csmat, ceig)
2013 CALL parallel_gemm("N", "T", nao, nao, nao, cone, csmat, cwork, &
2014 czero, kp%cshalf)
2015 END IF
2016 END DO
2017
2018 ! Clean up communication
2019 indx = 0
2020 DO ikp = 1, kplocal
2021 DO igroup = 1, nkp_groups
2022 ! number of current kpoint
2023 ik = kp_dist(1, igroup) + ikp - 1
2024 my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
2025 indx = indx + 1
2026 IF (use_real_wfn) THEN
2027 CALL cp_fm_cleanup_copy_general(info(indx, 1))
2028 ELSE
2029 CALL cp_fm_cleanup_copy_general(info(indx, 1))
2030 CALL cp_fm_cleanup_copy_general(info(indx, 2))
2031 END IF
2032 END DO
2033 END DO
2034
2035 ! All done
2036 DEALLOCATE (info)
2037 DEALLOCATE (eigenvalues, ceig)
2038
2039 CALL dbcsr_deallocate_matrix(rmatrix)
2040 CALL dbcsr_deallocate_matrix(cmatrix)
2041 CALL dbcsr_deallocate_matrix(tmpmat)
2042
2043 IF (use_real_wfn) THEN
2044 CALL cp_fm_release(rsmat)
2045 ELSE
2046 CALL cp_cfm_release(csmat)
2047 CALL cp_cfm_release(cwork)
2048 END IF
2049 CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
2050
2051 CALL timestop(handle)
2052
2053 END SUBROUTINE diag_kp_smat
2054
2055END MODULE qs_scf_diagonalization
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:120
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:56
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:89
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:84
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
qs_density_matrices::calculate_density_matrix
Definition qs_density_matrices.F:58
qs_mo_methods::calculate_subspace_eigenvalues
Definition qs_mo_methods.F:67
qs_mo_occupation::set_mo_occupation
Definition qs_mo_occupation.F:45
qs_scf_methods::combine_ks_matrices
Definition qs_scf_methods.F:72
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition cp_array_utils.F:18
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_scale_and_add
subroutine, public cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b).
Definition cp_cfm_basic_linalg.F:240
cp_cfm_basic_linalg::cp_cfm_scale_and_add_fm
subroutine, public cp_cfm_scale_and_add_fm(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b). where b is a real matrix (adapted from cp_cf...
Definition cp_cfm_basic_linalg.F:351
cp_cfm_basic_linalg::cp_cfm_column_scale
subroutine, public cp_cfm_column_scale(matrix_a, scaling)
Scales columns of the full matrix by corresponding factors.
Definition cp_cfm_basic_linalg.F:633
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_geeig
subroutine, public cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition cp_cfm_diag.F:192
cp_cfm_diag::cp_cfm_heevd
subroutine, public cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
Perform a diagonalisation of a complex matrix.
Definition cp_cfm_diag.F:60
cp_cfm_diag::cp_cfm_geeig_canon
subroutine, public cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
General Eigenvalue Problem AX = BXE Use canonical orthogonalization.
Definition cp_cfm_diag.F:257
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:182
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:125
cp_cfm_types::cp_cfm_set_all
subroutine, public cp_cfm_set_all(matrix, alpha, beta)
Set all elements of the full matrix to alpha. Besides, set all diagonal matrix elements to beta (if g...
Definition cp_cfm_types.F:202
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition cp_cfm_types.F:788
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:524
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:370
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1194
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:552
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:214
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:111
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition cp_fm_basic_linalg.F:1948
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1812
cp_fm_basic_linalg::cp_fm_symm
subroutine, public cp_fm_symm(side, uplo, m, n, alpha, matrix_a, matrix_b, beta, matrix_c)
computes matrix_c = beta * matrix_c + alpha * matrix_a * matrix_b computes matrix_c = beta * matrix_c...
Definition cp_fm_basic_linalg.F:625
cp_fm_cholesky
various cholesky decomposition related routines
Definition cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_restore
subroutine, public cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, pos, transa)
apply Cholesky decomposition op can be "SOLVE" (out = U^-1 * in) or "MULTIPLY" (out = U * in) pos can...
Definition cp_fm_cholesky.F:287
cp_fm_cholesky::cp_fm_cholesky_reduce
subroutine, public cp_fm_cholesky_reduce(matrix, matrixb, itype)
reduce a matrix pencil A,B to normal form B has to be cholesky decomposed with cp_fm_cholesky_decompo...
Definition cp_fm_cholesky.F:229
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::fm_diag_type_cusolver
integer, parameter, public fm_diag_type_cusolver
Definition cp_fm_diag.F:102
cp_fm_diag::cp_fm_geeig
subroutine, public cp_fm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition cp_fm_diag.F:1360
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:239
cp_fm_diag::diag_type
integer, save, public diag_type
Definition cp_fm_diag.F:87
cp_fm_diag::cp_fm_geeig_canon
subroutine, public cp_fm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
General Eigenvalue Problem AX = BXE Use canonical diagonalization : U*s**(-1/2)
Definition cp_fm_diag.F:1415
cp_fm_diag::cusolver_generalized
logical, save, public cusolver_generalized
Definition cp_fm_diag.F:96
cp_fm_pool_types
pool for for elements that are retained and released
Definition cp_fm_pool_types.F:14
cp_fm_pool_types::fm_pool_create_fm
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
Definition cp_fm_pool_types.F:185
cp_fm_pool_types::fm_pool_give_back_fm
subroutine, public fm_pool_give_back_fm(pool, element)
returns the element to the pool
Definition cp_fm_pool_types.F:224
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_start_copy_general
subroutine, public cp_fm_start_copy_general(source, destination, para_env, info)
Initiates the copy operation: get distribution data, post MPI isend and irecvs.
Definition cp_fm_types.F:1645
cp_fm_types::cp_fm_cleanup_copy_general
subroutine, public cp_fm_cleanup_copy_general(info)
Completes the copy operation: wait for comms clean up MPI state.
Definition cp_fm_types.F:2023
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1087
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:169
cp_fm_types::cp_fm_add_to_element
subroutine, public cp_fm_add_to_element(matrix, irow_global, icol_global, alpha)
...
Definition cp_fm_types.F:2170
cp_fm_types::cp_fm_finish_copy_general
subroutine, public cp_fm_finish_copy_general(destination, info)
Completes the copy operation: wait for comms, unpack, clean up MPI state.
Definition cp_fm_types.F:1959
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::core_guess
integer, parameter, public core_guess
Definition input_constants.F:202
input_constants::cholesky_restore
integer, parameter, public cholesky_restore
Definition input_constants.F:196
input_constants::cholesky_dbcsr
integer, parameter, public cholesky_dbcsr
Definition input_constants.F:196
input_constants::cholesky_off
integer, parameter, public cholesky_off
Definition input_constants.F:196
input_constants::cholesky_reduce
integer, parameter, public cholesky_reduce
Definition input_constants.F:196
input_constants::high_spin_roks
integer, parameter, public high_spin_roks
Definition input_constants.F:789
input_constants::cholesky_inverse
integer, parameter, public cholesky_inverse
Definition input_constants.F:196
input_constants::general_roks
integer, parameter, public general_roks
Definition input_constants.F:789
input_constants::restart_guess
integer, parameter, public restart_guess
Definition input_constants.F:202
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kpoint_methods
Routines needed for kpoint calculation.
Definition kpoint_methods.F:15
kpoint_methods::rskp_transform
subroutine, public rskp_transform(rmatrix, cmatrix, rsmat, ispin, xkp, cell_to_index, sab_nl, is_complex, rs_sign)
Transformation of real space matrices to a kpoint.
Definition kpoint_methods.F:819
kpoint_methods::kpoint_density_transform
subroutine, public kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwork, for_aux_fit, pmat_ext)
generate real space density matrices in DBCSR format
Definition kpoint_methods.F:1456
kpoint_methods::kpoint_density_matrices
subroutine, public kpoint_density_matrices(kpoint, energy_weighted, for_aux_fit)
Calculate kpoint density matrices (rho(k), owned by kpoint groups)
Definition kpoint_methods.F:1136
kpoint_methods::kpoint_set_mo_occupation
subroutine, public kpoint_set_mo_occupation(kpoint, smear, probe)
Given the eigenvalues of all kpoints, calculates the occupation numbers.
Definition kpoint_methods.F:919
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:371
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:136
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:153
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:143
mathconstants::gaussi
complex(kind=dp), parameter, public gaussi
Definition mathconstants.F:142
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:143
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
preconditioner
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition preconditioner.F:15
preconditioner::restart_preconditioner
subroutine, public restart_preconditioner(qs_env, preconditioner, prec_type, nspins)
Allows for a restart of the preconditioner depending on the method it purges all arrays or keeps them...
Definition preconditioner.F:244
preconditioner::prepare_preconditioner
subroutine, public prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, ot_preconditioner, prec_type, solver_type, energy_gap, nspins, has_unit_metric, convert_to_dbcsr, chol_type, full_mo_set)
...
Definition preconditioner.F:309
qs_density_matrices
collects routines that calculate density matrices
Definition qs_density_matrices.F:14
qs_density_mixing_types
module that contains the definitions of the scf types
Definition qs_density_mixing_types.F:14
qs_density_mixing_types::direct_mixing_nr
integer, parameter, public direct_mixing_nr
Definition qs_density_mixing_types.F:45
qs_density_mixing_types::gspace_mixing_nr
integer, parameter, public gspace_mixing_nr
Definition qs_density_mixing_types.F:45
qs_diis
Apply the direct inversion in the iterative subspace (DIIS) of Pulay in the framework of an SCF itera...
Definition qs_diis.F:21
qs_diis::qs_diis_b_info_kp
subroutine, public qs_diis_b_info_kp(diis_buffer, ib, nb)
Update info about the current buffer step ib and the current number of buffers nb.
Definition qs_diis.F:989
qs_diis::qs_diis_b_calc_err_kp
subroutine, public qs_diis_b_calc_err_kp(diis_buffer, ib, mos, kc, sc, ispin, ikp, nkp_local, scf_section)
Calculate and store the error for a given k-point.
Definition qs_diis.F:1013
qs_diis::qs_diis_b_step_kp
subroutine, public qs_diis_b_step_kp(diis_buffer, coeffs, ib, nb, delta, error_max, diis_step, eps_diis, nspin, nkp, nkp_local, nmixing, scf_section, para_env)
Update the SCF DIIS buffer, and if appropriate does a diis step, for k-points.
Definition qs_diis.F:1100
qs_diis::qs_diis_b_step
subroutine, public qs_diis_b_step(diis_buffer, mo_array, kc, sc, delta, error_max, diis_step, eps_diis, nmixing, s_matrix, scf_section, roks)
Update the SCF DIIS buffer, and if appropriate does a diis step.
Definition qs_diis.F:230
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:530
qs_gspace_mixing
Definition qs_gspace_mixing.F:9
qs_gspace_mixing::gspace_mixing
subroutine, public gspace_mixing(qs_env, mixing_method, mixing_store, rho, para_env, iter_count)
Driver for the g-space mixing, calls the proper routine given the requested method.
Definition qs_gspace_mixing.F:66
qs_ks_methods
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition qs_ks_methods.F:22
qs_ks_methods::qs_ks_update_qs_env
subroutine, public qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, print_active)
updates the Kohn Sham matrix of the given qs_env (facility method)
Definition qs_ks_methods.F:1115
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::qs_ks_did_change
subroutine, public qs_ks_did_change(ks_env, s_mstruct_changed, rho_changed, potential_changed, full_reset)
tells that some of the things relevant to the ks calculation did change. has to be called when change...
Definition qs_ks_types.F:956
qs_matrix_pools
wrapper for the pools of matrixes
Definition qs_matrix_pools.F:14
qs_matrix_pools::mpools_get
subroutine, public mpools_get(mpools, ao_mo_fm_pools, ao_ao_fm_pools, mo_mo_fm_pools, ao_mosub_fm_pools, mosub_mosub_fm_pools, maxao_maxmo_fm_pool, maxao_maxao_fm_pool, maxmo_maxmo_fm_pool)
returns various attributes of the mpools (notably the pools contained in it)
Definition qs_matrix_pools.F:139
qs_mixing_utils
Definition qs_mixing_utils.F:9
qs_mixing_utils::charge_mixing_init
elemental subroutine, public charge_mixing_init(mixing_store)
initialiation needed when charge mixing is used
Definition qs_mixing_utils.F:643
qs_mixing_utils::mixing_init
subroutine, public mixing_init(mixing_method, rho, mixing_store, para_env, rho_atom)
initialiation needed when gspace mixing is used
Definition qs_mixing_utils.F:397
qs_mixing_utils::mixing_allocate
subroutine, public mixing_allocate(qs_env, mixing_method, p_mix_new, p_delta, nspins, mixing_store)
allocation needed when density mixing is used
Definition qs_mixing_utils.F:110
qs_mixing_utils::self_consistency_check
subroutine, public self_consistency_check(rho_ao, p_delta, para_env, p_out, delta)
...
Definition qs_mixing_utils.F:57
qs_mo_methods
collects routines that perform operations directly related to MOs
Definition qs_mo_methods.F:14
qs_mo_occupation
Set occupation of molecular orbitals.
Definition qs_mo_occupation.F:15
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:399
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_ot_eigensolver
an eigen-space solver for the generalised symmetric eigenvalue problem for sparse matrices,...
Definition qs_ot_eigensolver.F:13
qs_ot_eigensolver::ot_eigensolver
subroutine, public ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, matrix_c_fm, preconditioner, eps_gradient, iter_max, size_ortho_space, silent, ot_settings)
...
Definition qs_ot_eigensolver.F:74
qs_rho_atom_types
Definition qs_rho_atom_types.F:8
qs_rho_methods
methods of the rho structure (defined in qs_rho_types)
Definition qs_rho_methods.F:15
qs_rho_methods::qs_rho_update_rho
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
Definition qs_rho_methods.F:378
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_scf_block_davidson
module that contains the algorithms to perform an iterative diagonalization by the block-Davidson app...
Definition qs_scf_block_davidson.F:18
qs_scf_block_davidson::generate_extended_space_sparse
subroutine, public generate_extended_space_sparse(bdav_env, mo_set, matrix_h, matrix_s, output_unit, preconditioner)
...
Definition qs_scf_block_davidson.F:479
qs_scf_block_davidson::generate_extended_space
subroutine, public generate_extended_space(bdav_env, mo_set, matrix_h, matrix_s, output_unit, preconditioner)
...
Definition qs_scf_block_davidson.F:83
qs_scf_diagonalization
Different diagonalization schemes that can be used for the iterative solution of the eigenvalue probl...
Definition qs_scf_diagonalization.F:15
qs_scf_diagonalization::general_eigenproblem
subroutine, public general_eigenproblem(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step)
the inner loop of scf, specific to diagonalization with S matrix basically, in goes the ks matrix out...
Definition qs_scf_diagonalization.F:163
qs_scf_diagonalization::diag_subspace_allocate
subroutine, public diag_subspace_allocate(subspace_env, qs_env, mos)
...
Definition qs_scf_diagonalization.F:1002
qs_scf_diagonalization::do_ot_diag
subroutine, public do_ot_diag(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step)
the inner loop of scf, specific to iterative diagonalization using OT with S matrix; basically,...
Definition qs_scf_diagonalization.F:1168
qs_scf_diagonalization::do_block_davidson_diag
subroutine, public do_block_davidson_diag(qs_env, scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, check_moconv_only)
iterative diagonalization using the block davidson space approach
Definition qs_scf_diagonalization.F:1748
qs_scf_diagonalization::do_roks_diag
subroutine, public do_roks_diag(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step, orthogonal_basis)
Solve a set restricted open Kohn-Sham (ROKS) equations based on the alpha and beta Kohn-Sham matrices...
Definition qs_scf_diagonalization.F:1274
qs_scf_diagonalization::do_scf_diag_subspace
subroutine, public do_scf_diag_subspace(qs_env, scf_env, subspace_env, mos, rho, ks_env, scf_section, scf_control)
inner loop within MOS subspace, to refine occupation and density, before next diagonalization of the ...
Definition qs_scf_diagonalization.F:794
qs_scf_diagonalization::do_block_krylov_diag
subroutine, public do_block_krylov_diag(scf_env, mos, matrix_ks, scf_control, scf_section, check_moconv_only)
iterative diagonalization using the block Krylov-space approach
Definition qs_scf_diagonalization.F:1548
qs_scf_diagonalization::do_special_diag
subroutine, public do_special_diag(scf_env, mos, matrix_ks, scf_control, scf_section, diis_step)
the inner loop of scf, specific to diagonalization without S matrix basically, in goes the ks matrix ...
Definition qs_scf_diagonalization.F:1072
qs_scf_diagonalization::do_general_diag
subroutine, public do_general_diag(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step, probe)
...
Definition qs_scf_diagonalization.F:373
qs_scf_diagonalization::do_general_diag_kp
subroutine, public do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, update_p, diis_step, diis_error, qs_env, probe)
Kpoint diagonalization routine Transforms matrices to kpoint, distributes kpoint groups,...
Definition qs_scf_diagonalization.F:437
qs_scf_diagonalization::diag_kp_smat
subroutine, public diag_kp_smat(matrix_s, kpoints, fmwork)
Kpoint diagonalization routine Transforms matrices to kpoint, distributes kpoint groups,...
Definition qs_scf_diagonalization.F:1855
qs_scf_lanczos
module that contains the algorithms to perform an iterative diagonalization by the block-Lanczos appr...
Definition qs_scf_lanczos.F:15
qs_scf_lanczos::lanczos_refinement
subroutine, public lanczos_refinement(krylov_space, ks, c0, c1, eval, nao, eps_iter, ispin, check_moconv_only)
lanczos refinement by blocks of non-converged MOs
Definition qs_scf_lanczos.F:184
qs_scf_lanczos::lanczos_refinement_2v
subroutine, public lanczos_refinement_2v(krylov_space, ks, c0, c1, eval, nao, eps_iter, ispin, check_moconv_only)
...
Definition qs_scf_lanczos.F:519
qs_scf_methods
groups fairly general SCF methods, so that modules other than qs_scf can use them too split off from ...
Definition qs_scf_methods.F:19
qs_scf_methods::eigensolver_simple
subroutine, public eigensolver_simple(matrix_ks, mo_set, work, do_level_shift, level_shift, use_jacobi, jacobi_threshold)
...
Definition qs_scf_methods.F:474
qs_scf_methods::eigensolver_dbcsr
subroutine, public eigensolver_dbcsr(matrix_ks, matrix_ks_fm, mo_set, ortho_dbcsr, ksbuf1, ksbuf2)
...
Definition qs_scf_methods.F:312
qs_scf_methods::scf_env_density_mixing
subroutine, public scf_env_density_mixing(p_mix_new, mixing_store, rho_ao, para_env, iter_delta, iter_count, diis, invert)
perform (if requested) a density mixing
Definition qs_scf_methods.F:96
qs_scf_methods::eigensolver
subroutine, public eigensolver(matrix_ks_fm, mo_set, ortho, work, cholesky_method, do_level_shift, level_shift, matrix_u_fm, use_jacobi)
Diagonalise the Kohn-Sham matrix to get a new set of MO eigen- vectors and MO eigenvalues....
Definition qs_scf_methods.F:167
qs_scf_methods::eigensolver_symm
subroutine, public eigensolver_symm(matrix_ks_fm, mo_set, ortho, work, do_level_shift, level_shift, matrix_u_fm, use_jacobi, jacobi_threshold, ortho_red, work_red, matrix_ks_fm_red, matrix_u_fm_red)
...
Definition qs_scf_methods.F:378
qs_scf_methods::eigensolver_generalized
subroutine, public eigensolver_generalized(matrix_ks_fm, matrix_s, mo_set, work)
Solve the generalized eigenvalue problem using cusolverMpSygvd.
Definition qs_scf_methods.F:265
qs_scf_types
module that contains the definitions of the scf types
Definition qs_scf_types.F:14
scf_control_types
parameters that control an scf iteration
Definition scf_control_types.F:17
cp_array_utils::cp_1d_r_p_type
represent a pointer to a 1d array
Definition cp_array_utils.F:99
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:69
cp_control_types::dft_control_type
Definition cp_control_types.F:631
cp_control_types::hairy_probes_type
Definition cp_control_types.F:45
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_pool_types::cp_fm_pool_p_type
to create arrays of pools
Definition cp_fm_pool_types.F:70
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::copy_info_type
Stores the state of a copy between cp_fm_start_copy_general and cp_fm_finish_copy_general.
Definition cp_fm_types.F:141
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:115
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_env_type
Keeps information about a specific k-point.
Definition kpoint_types.F:77
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:156
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_matrix_pools::qs_matrix_pools_type
container for the pools of matrixes used by qs
Definition qs_matrix_pools.F:66
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_atom_types::rho_atom_type
Definition qs_rho_atom_types.F:23
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:97
qs_scf_types::subspace_env_type
Definition qs_scf_types.F:84
scf_control_types::scf_control_type
Definition scf_control_types.F:124