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