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