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