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