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_blacs_env, ONLY: cp_blacs_env_type
18 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,&
19 : cp_cfm_scale_and_add,&
20 : cp_cfm_scale_and_add_fm
21 : USE cp_cfm_diag, ONLY: cp_cfm_geeig,&
22 : cp_cfm_geeig_canon,&
23 : cp_cfm_geeig_canon_local,&
24 : cp_cfm_geeig_local,&
25 : cp_cfm_heevd
26 : USE cp_cfm_types, ONLY: cp_cfm_create,&
27 : cp_cfm_release,&
28 : cp_cfm_set_all,&
29 : cp_cfm_to_cfm,&
30 : cp_cfm_to_fm,&
31 : cp_cfm_type
32 : USE cp_control_types, ONLY: dft_control_type,&
33 : hairy_probes_type
34 : USE cp_dbcsr_api, ONLY: &
35 : dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
36 : dbcsr_get_readonly_block_p, dbcsr_init_p, dbcsr_p_type, dbcsr_release_p, dbcsr_set, &
37 : dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
38 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
39 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
40 : copy_fm_to_dbcsr,&
41 : copy_fm_to_dbcsr_bc,&
42 : cp_dbcsr_sm_fm_multiply,&
43 : dbcsr_allocate_matrix_set
44 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
45 : cp_fm_scale_and_add,&
46 : cp_fm_symm,&
47 : cp_fm_uplo_to_full
48 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_reduce,&
49 : cp_fm_cholesky_restore
50 : USE cp_fm_diag, ONLY: FM_DIAG_TYPE_CUSOLVER,&
51 : FM_DIAG_TYPE_DLAF,&
52 : choose_eigv_solver,&
53 : cp_fm_geeig,&
54 : cp_fm_geeig_canon,&
55 : diag_check_requested,&
56 : diag_type,&
57 : direct_generalized_diagonalization
58 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
59 : fm_pool_create_fm,&
60 : fm_pool_give_back_fm
61 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
62 : cp_fm_struct_release,&
63 : cp_fm_struct_type
64 : USE cp_fm_types, ONLY: &
65 : copy_info_type, cp_fm_add_to_element, cp_fm_cleanup_copy_general, cp_fm_create, &
66 : cp_fm_finish_copy_general, cp_fm_get_info, cp_fm_release, cp_fm_set_all, &
67 : cp_fm_start_copy_general, cp_fm_to_fm, cp_fm_type
68 : USE cp_log_handling, ONLY: cp_get_default_logger,&
69 : cp_logger_get_default_unit_nr,&
70 : cp_logger_type
71 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
72 : cp_print_key_unit_nr,&
73 : debug_print_level
74 : USE ieee_arithmetic, ONLY: ieee_is_finite
75 : USE input_constants, ONLY: &
76 : cholesky_dbcsr, cholesky_inverse, cholesky_off, cholesky_reduce, cholesky_restore, &
77 : core_guess, diag_update_method_adiis, general_roks, high_spin_roks, ot_precond_full_all, &
78 : ot_precond_full_kinetic, ot_precond_full_single_inverse, ot_precond_none, &
79 : ot_precond_s_inverse, restart_guess
80 : USE input_cp2k_kpoints, ONLY: lattice_fft_auto,&
81 : lattice_fft_off,&
82 : lattice_fft_on
83 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
84 : section_vals_type
85 : USE kinds, ONLY: dp,&
86 : int_8
87 : USE kpoint_methods, ONLY: &
88 : kp_transform_plan_create, kp_transform_plan_release, kp_transform_plan_type, &
89 : kpoint_density_matrices, kpoint_density_transform, kpoint_set_mo_occupation, &
90 : rskp_grid_type, rskp_transform, rskp_transform_grid_extract, rskp_transform_grid_prepare, &
91 : rskp_transform_grid_release
92 : USE kpoint_types, ONLY: get_kpoint_info,&
93 : kpoint_env_type,&
94 : kpoint_type
95 : USE machine, ONLY: m_flush,&
96 : m_walltime
97 : USE mathconstants, ONLY: gaussi,&
98 : twopi,&
99 : z_one,&
100 : z_zero
101 : USE message_passing, ONLY: mp_para_env_type
102 : USE parallel_gemm_api, ONLY: parallel_gemm
103 : USE preconditioner, ONLY: dbcsr_pair_to_cfm,&
104 : make_preconditioner_complex_full_kinetic,&
105 : make_preconditioner_complex_full_s_inverse,&
106 : prepare_preconditioner,&
107 : restart_preconditioner
108 : USE preconditioner_makes, ONLY: make_complex_full_all,&
109 : make_complex_full_single_inverse
110 : USE preconditioner_types, ONLY: destroy_preconditioner,&
111 : init_preconditioner
112 : USE qs_density_matrices, ONLY: calculate_density_matrix
113 : USE qs_density_mixing_types, ONLY: direct_mixing_nr,&
114 : gspace_mixing_nr
115 : USE qs_diis, ONLY: qs_diis_b_calc_err_kp,&
116 : qs_diis_b_check_i_alloc_kp,&
117 : qs_diis_b_clear,&
118 : qs_diis_b_clear_kp,&
119 : qs_diis_b_info_kp,&
120 : qs_diis_b_step,&
121 : qs_diis_b_step_kp
122 : USE qs_energy_types, ONLY: qs_energy_type
123 : USE qs_environment_types, ONLY: get_qs_env,&
124 : qs_environment_type
125 : USE qs_gspace_mixing, ONLY: gspace_mixing
126 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
127 : USE qs_ks_types, ONLY: qs_ks_did_change,&
128 : qs_ks_env_type
129 : USE qs_matrix_pools, ONLY: mpools_get,&
130 : qs_matrix_pools_type
131 : USE qs_mixing_utils, ONLY: charge_mixing_init,&
132 : mixing_allocate,&
133 : mixing_init,&
134 : self_consistency_check
135 : USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
136 : USE qs_mo_occupation, ONLY: set_mo_occupation
137 : USE qs_mo_types, ONLY: get_mo_set,&
138 : mo_set_type
139 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
140 : USE qs_ot_eigensolver, ONLY: ot_eigensolver,&
141 : ot_eigensolver_complex
142 : USE qs_rho_atom_types, ONLY: rho_atom_type
143 : USE qs_rho_methods, ONLY: qs_rho_update_rho
144 : USE qs_rho_types, ONLY: qs_rho_get,&
145 : qs_rho_type
146 : USE qs_scf_block_davidson, ONLY: generate_extended_space,&
147 : generate_extended_space_c,&
148 : generate_extended_space_sparse
149 : USE qs_scf_lanczos, ONLY: lanczos_refinement,&
150 : lanczos_refinement_2v
151 : USE qs_scf_methods, ONLY: combine_ks_matrices,&
152 : eigensolver,&
153 : eigensolver_dbcsr,&
154 : eigensolver_generalized,&
155 : eigensolver_simple,&
156 : eigensolver_symm,&
157 : scf_env_density_mixing
158 : USE qs_scf_types, ONLY: qs_scf_env_type,&
159 : subspace_env_type
160 : USE scf_control_types, ONLY: scf_control_type
161 :
162 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, &
163 : !$ omp_get_thread_num
164 : #include "./base/base_uses.f90"
165 :
166 : IMPLICIT NONE
167 :
168 : PRIVATE
169 :
170 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_diagonalization'
171 :
172 : TYPE, PRIVATE :: kp_diag_omp_workspace_type
173 : TYPE(cp_cfm_type) :: cksmat, cmos, csmat, csmat_base, cwork
174 : END TYPE kp_diag_omp_workspace_type
175 :
176 : PUBLIC :: do_general_diag, do_general_diag_kp, do_roks_diag, &
177 : do_special_diag, do_ot_diag, do_ot_diag_kp, do_block_davidson_diag, &
178 : do_block_davidson_diag_kp, &
179 : do_block_krylov_diag, do_scf_diag_subspace, diag_subspace_allocate, &
180 : general_eigenproblem, diag_kp_smat, diag_kp_basic, &
181 : qs_scf_distribute_kpoint_operator, qs_scf_get_ot_kpoint_operator
182 :
183 : CONTAINS
184 :
185 : ! **************************************************************************************************
186 : !> \brief Transform and redistribute a real-space operator to all local complex K-point groups.
187 : !> \param matrix_rs real-space image matrices of the operator
188 : !> \param kpoints K-point environment and group distribution
189 : !> \param fmwork full-matrix workspaces on the global communicator
190 : !> \param operator_re local real K-point operators
191 : !> \param operator_im local imaginary K-point operators
192 : ! **************************************************************************************************
193 56 : SUBROUTINE qs_scf_distribute_kpoint_operator( &
194 56 : matrix_rs, kpoints, fmwork, operator_re, operator_im)
195 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_rs
196 : TYPE(kpoint_type), POINTER :: kpoints
197 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: fmwork
198 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: operator_re, operator_im
199 :
200 : INTEGER :: ispin, local_kpoint, nlocal_kpoints, &
201 : nspin
202 56 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
203 56 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
204 56 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
205 : TYPE(dbcsr_type), POINTER :: matrix_im, matrix_re, matrix_work
206 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
207 56 : POINTER :: sab_nl
208 :
209 56 : NULLIFY (cell_to_index, kp_dist, matrix_im, matrix_re, matrix_work, sab_nl, xkp)
210 0 : CPASSERT(ASSOCIATED(matrix_rs))
211 56 : CPASSERT(ASSOCIATED(kpoints))
212 56 : CPASSERT(SIZE(fmwork) >= 2)
213 56 : nlocal_kpoints = SIZE(kpoints%kp_env)
214 56 : nspin = SIZE(matrix_rs, 1)
215 168 : CPASSERT(ALL(SHAPE(operator_re) == [nlocal_kpoints, nspin]))
216 168 : CPASSERT(ALL(SHAPE(operator_im) == [nlocal_kpoints, nspin]))
217 : CALL get_kpoint_info(kpoints, kp_dist=kp_dist, xkp=xkp, &
218 56 : cell_to_index=cell_to_index, sab_nl=sab_nl)
219 :
220 56 : ALLOCATE (matrix_re, matrix_im, matrix_work)
221 56 : CALL dbcsr_create(matrix_re, template=matrix_rs(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
222 56 : CALL dbcsr_create(matrix_im, template=matrix_rs(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
223 56 : CALL dbcsr_create(matrix_work, template=matrix_rs(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
224 56 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_re, sab_nl)
225 56 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_im, sab_nl)
226 :
227 112 : DO local_kpoint = 1, nlocal_kpoints
228 168 : DO ispin = 1, nspin
229 : CALL distribute_kpoint_operator_channel( &
230 : matrix_rs, kpoints, local_kpoint, ispin, kp_dist, xkp, cell_to_index, sab_nl, &
231 : matrix_re, matrix_im, matrix_work, fmwork, &
232 112 : operator_re(local_kpoint, ispin), operator_im(local_kpoint, ispin))
233 : END DO
234 : END DO
235 :
236 56 : CALL dbcsr_deallocate_matrix(matrix_re)
237 56 : CALL dbcsr_deallocate_matrix(matrix_im)
238 56 : CALL dbcsr_deallocate_matrix(matrix_work)
239 :
240 56 : END SUBROUTINE qs_scf_distribute_kpoint_operator
241 :
242 : ! **************************************************************************************************
243 : !> \brief Transform and redistribute one local k-point channel while all k-point
244 : !> groups participate in the same global communication sequence.
245 : !> \param matrix_rs real-space image matrices of the operator
246 : !> \param kpoints K-point environment and group distribution
247 : !> \param local_kpoint local k-point index, shared by all groups
248 : !> \param ispin spin component of matrix_rs
249 : !> \param kp_dist global k-point ranges of the groups
250 : !> \param xkp K-point coordinates
251 : !> \param cell_to_index ...
252 : !> \param sab_nl ...
253 : !> \param matrix_re symmetric real-part scratch on the global DBCSR distribution
254 : !> \param matrix_im antisymmetric imaginary-part scratch on the global DBCSR distribution
255 : !> \param matrix_work no-symmetry desymmetrization scratch
256 : !> \param fmwork full-matrix workspaces on the global communicator
257 : !> \param operator_re local real K-point operator on the owning group
258 : !> \param operator_im local imaginary K-point operator on the owning group
259 : ! **************************************************************************************************
260 4384 : SUBROUTINE distribute_kpoint_operator_channel( &
261 : matrix_rs, kpoints, local_kpoint, ispin, kp_dist, xkp, cell_to_index, sab_nl, &
262 4384 : matrix_re, matrix_im, matrix_work, fmwork, operator_re, operator_im)
263 :
264 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_rs
265 : TYPE(kpoint_type), POINTER :: kpoints
266 : INTEGER, INTENT(IN) :: local_kpoint, ispin
267 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
268 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
269 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
270 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
271 : POINTER :: sab_nl
272 : TYPE(dbcsr_type), POINTER :: matrix_re, matrix_im, matrix_work
273 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: fmwork
274 : TYPE(cp_fm_type), INTENT(INOUT) :: operator_re, operator_im
275 :
276 : INTEGER :: igroup, ikpoint, nkp_groups
277 : LOGICAL :: my_kpoint_group
278 74528 : TYPE(copy_info_type) :: info_im, info_re
279 : TYPE(cp_fm_type) :: fmdummy
280 : TYPE(mp_para_env_type), POINTER :: para_env
281 :
282 4384 : CPASSERT(ASSOCIATED(matrix_rs))
283 4384 : CPASSERT(ASSOCIATED(kpoints))
284 4384 : CPASSERT(SIZE(fmwork) >= 2)
285 4384 : CPASSERT(local_kpoint >= 1 .AND. local_kpoint <= SIZE(kpoints%kp_env))
286 4384 : CPASSERT(ispin >= 1 .AND. ispin <= SIZE(matrix_rs, 1))
287 4384 : nkp_groups = SIZE(kp_dist, 2)
288 4384 : para_env => kpoints%blacs_env_all%para_env
289 :
290 10952 : DO igroup = 1, nkp_groups
291 6568 : ikpoint = kp_dist(1, igroup) + local_kpoint - 1
292 : my_kpoint_group = ikpoint >= kpoints%kp_range(1) .AND. &
293 6568 : ikpoint <= kpoints%kp_range(2)
294 6568 : CALL dbcsr_set(matrix_re, 0.0_dp)
295 6568 : CALL dbcsr_set(matrix_im, 0.0_dp)
296 : CALL rskp_transform( &
297 : rmatrix=matrix_re, cmatrix=matrix_im, rsmat=matrix_rs, ispin=ispin, &
298 6568 : xkp=xkp(1:3, ikpoint), cell_to_index=cell_to_index, sab_nl=sab_nl)
299 6568 : CALL dbcsr_desymmetrize(matrix_re, matrix_work)
300 6568 : CALL copy_dbcsr_to_fm(matrix_work, fmwork(1))
301 6568 : CALL dbcsr_desymmetrize(matrix_im, matrix_work)
302 6568 : CALL copy_dbcsr_to_fm(matrix_work, fmwork(2))
303 6568 : IF (my_kpoint_group) THEN
304 4384 : CALL cp_fm_start_copy_general(fmwork(1), operator_re, para_env, info_re)
305 4384 : CALL cp_fm_start_copy_general(fmwork(2), operator_im, para_env, info_im)
306 4384 : CALL cp_fm_finish_copy_general(operator_re, info_re)
307 4384 : CALL cp_fm_finish_copy_general(operator_im, info_im)
308 : ELSE
309 2184 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info_re)
310 2184 : CALL cp_fm_start_copy_general(fmwork(2), fmdummy, para_env, info_im)
311 : END IF
312 6568 : CALL cp_fm_cleanup_copy_general(info_re)
313 10952 : CALL cp_fm_cleanup_copy_general(info_im)
314 : END DO
315 :
316 8768 : END SUBROUTINE distribute_kpoint_operator_channel
317 :
318 : ! **************************************************************************************************
319 : !> \brief Build one complex k-point operator for OT in the local k-point-group context.
320 : !> A single group keeps the established direct real-space transform. Multiple groups use
321 : !> the redistributed full-matrix cache because their communicators do not own the complete
322 : !> real-space DBCSR matrix.
323 : !> \param matrix_rs global real-space operator
324 : !> \param kpoints k-point environment
325 : !> \param kp local k-point environment
326 : !> \param ispin spin component of matrix_rs
327 : !> \param cache_re local real part for the multiple-group path
328 : !> \param cache_im local imaginary part for the multiple-group path
329 : !> \param matrix_re local real DBCSR operator
330 : !> \param matrix_im local imaginary DBCSR operator
331 : ! **************************************************************************************************
332 9608 : SUBROUTINE qs_scf_get_ot_kpoint_operator(matrix_rs, kpoints, kp, ispin, cache_re, cache_im, &
333 : matrix_re, matrix_im)
334 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_rs
335 : TYPE(kpoint_type), POINTER :: kpoints
336 : TYPE(kpoint_env_type), POINTER :: kp
337 : INTEGER, INTENT(IN) :: ispin
338 : TYPE(cp_fm_type), INTENT(IN) :: cache_re, cache_im
339 : TYPE(dbcsr_type), POINTER :: matrix_re, matrix_im
340 :
341 : INTEGER :: nkp_groups
342 9608 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
343 9608 : REAL(dp), DIMENSION(:, :), POINTER :: xkp
344 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
345 9608 : POINTER :: sab_nl
346 :
347 9608 : NULLIFY (cell_to_index, sab_nl, xkp)
348 0 : CPASSERT(ASSOCIATED(matrix_rs))
349 9608 : CPASSERT(ispin >= 1 .AND. ispin <= SIZE(matrix_rs, 1))
350 9608 : CPASSERT(ASSOCIATED(kpoints))
351 9608 : CPASSERT(ASSOCIATED(kp))
352 :
353 9608 : CALL dbcsr_init_p(matrix_re)
354 9608 : CALL dbcsr_init_p(matrix_im)
355 : CALL get_kpoint_info(kpoints, nkp_groups=nkp_groups, xkp=xkp, &
356 9608 : cell_to_index=cell_to_index, sab_nl=sab_nl)
357 9608 : IF (nkp_groups == 1) THEN
358 6748 : CPASSERT(ASSOCIATED(xkp))
359 6748 : CPASSERT(ASSOCIATED(cell_to_index))
360 6748 : CPASSERT(ASSOCIATED(sab_nl))
361 : CALL dbcsr_create(matrix_re, template=matrix_rs(1, 1)%matrix, &
362 6748 : matrix_type=dbcsr_type_symmetric)
363 : CALL dbcsr_create(matrix_im, template=matrix_rs(1, 1)%matrix, &
364 6748 : matrix_type=dbcsr_type_antisymmetric)
365 6748 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_re, sab_nl)
366 6748 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_im, sab_nl)
367 6748 : CALL dbcsr_set(matrix_re, 0.0_dp)
368 6748 : CALL dbcsr_set(matrix_im, 0.0_dp)
369 : CALL rskp_transform(rmatrix=matrix_re, cmatrix=matrix_im, rsmat=matrix_rs, ispin=ispin, &
370 6748 : xkp=xkp(1:3, kp%nkpoint), cell_to_index=cell_to_index, sab_nl=sab_nl)
371 : ELSE
372 2860 : CALL copy_fm_to_dbcsr_bc(cache_re, matrix_re)
373 2860 : CALL copy_fm_to_dbcsr_bc(cache_im, matrix_im)
374 : END IF
375 9608 : END SUBROUTINE qs_scf_get_ot_kpoint_operator
376 :
377 : ! **************************************************************************************************
378 : !> \brief the inner loop of scf, specific to diagonalization with S matrix
379 : !> basically, in goes the ks matrix out goes a new p matrix
380 : !> \param scf_env ...
381 : !> \param mos ...
382 : !> \param matrix_ks ...
383 : !> \param matrix_s ...
384 : !> \param scf_control ...
385 : !> \param scf_section ...
386 : !> \param diis_step ...
387 : !> \par History
388 : !> 03.2006 created [Joost VandeVondele]
389 : ! **************************************************************************************************
390 :
391 99629 : SUBROUTINE general_eigenproblem(scf_env, mos, matrix_ks, &
392 : matrix_s, scf_control, scf_section, &
393 : diis_step)
394 :
395 : TYPE(qs_scf_env_type), POINTER :: scf_env
396 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
397 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
398 : TYPE(scf_control_type), POINTER :: scf_control
399 : TYPE(section_vals_type), POINTER :: scf_section
400 : LOGICAL, INTENT(INOUT) :: diis_step
401 :
402 : INTEGER :: cdiis_cholesky_method, cdiis_ncall, &
403 : ispin, nspin
404 : LOGICAL :: cdiis_step, do_level_shift, owns_ortho, &
405 : use_adiis, use_cdiis, use_jacobi, &
406 : use_subspace_fock
407 : REAL(KIND=dp) :: cdiis_error, diis_error, diis_weight, &
408 : eps_diis
409 : TYPE(cp_fm_type), POINTER :: ortho
410 : TYPE(dbcsr_type), POINTER :: effective_ks_dbcsr, ortho_dbcsr
411 :
412 99629 : nspin = SIZE(matrix_ks)
413 99629 : NULLIFY (effective_ks_dbcsr, ortho, ortho_dbcsr)
414 :
415 99629 : use_adiis = scf_control%diagonalization%update_method == diag_update_method_adiis
416 99629 : use_subspace_fock = .FALSE.
417 99629 : IF (use_adiis) THEN
418 98 : IF (.NOT. ASSOCIATED(scf_env%scf_subspace_buffer)) THEN
419 0 : CPABORT("ADIIS is only available in the main Quickstep SCF workflow")
420 : END IF
421 98 : use_subspace_fock = scf_env%scf_subspace_buffer%use_combined_fock
422 : END IF
423 :
424 99629 : eps_diis = scf_control%eps_diis
425 99629 : use_cdiis = .FALSE.
426 99629 : IF (use_adiis .AND. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
427 : IF (.NOT. scf_env%skip_diis .AND. scf_env%iter_count > 1 .AND. &
428 : .NOT. scf_env%adiis_check_next .AND. &
429 98 : scf_env%scf_subspace_buffer%diis_state_valid .AND. &
430 : scf_env%scf_subspace_buffer%nstored > 0) THEN
431 70 : use_cdiis = scf_env%iter_delta < eps_diis .OR. scf_env%scf_diis_buffer%ncall > 0
432 : END IF
433 : END IF
434 :
435 216598 : DO ispin = 1, nspin
436 216598 : IF (use_cdiis) THEN
437 96 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, scf_env%scf_work1(ispin))
438 116873 : ELSE IF (use_subspace_fock) THEN
439 : CALL copy_dbcsr_to_fm(scf_env%scf_subspace_buffer%combined_fock(ispin, 1)%matrix, &
440 32 : scf_env%scf_work1(ispin))
441 : ELSE
442 116841 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, scf_env%scf_work1(ispin))
443 : END IF
444 : END DO
445 :
446 99629 : IF (use_adiis) THEN
447 98 : diis_step = .FALSE.
448 98 : diis_error = scf_env%scf_subspace_buffer%last_old_fock_weight
449 98 : diis_weight = scf_env%scf_subspace_buffer%diis_weight
450 98 : IF (use_cdiis) THEN
451 48 : cdiis_ncall = scf_env%scf_diis_buffer%ncall
452 : CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
453 : scf_env%scf_work2, scf_env%iter_delta, cdiis_error, cdiis_step, &
454 48 : eps_diis, 2, s_matrix=matrix_s, scf_section=scf_section)
455 48 : IF (.NOT. ieee_is_finite(cdiis_error) .OR. cdiis_error >= eps_diis) THEN
456 0 : cdiis_step = .FALSE.
457 0 : IF (cdiis_ncall < scf_env%scf_diis_buffer%nbuffer) THEN
458 : ! The rejected state has not overwritten a full circular buffer.
459 0 : scf_env%scf_diis_buffer%ncall = cdiis_ncall
460 : ELSE
461 0 : CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
462 0 : diis_weight = 0.0_dp
463 : END IF
464 : END IF
465 48 : IF (cdiis_step) THEN
466 : ! Use one half-weight handover step; MAX_DIIS only controls the buffer capacity.
467 34 : diis_weight = MIN(1.0_dp, diis_weight + 0.5_dp)
468 : ELSE
469 : ! A rejected or deferred CDIIS step must not leave a partial raw-Fock blend active.
470 14 : diis_weight = 0.0_dp
471 : END IF
472 48 : IF (diis_weight < 1.0_dp) THEN
473 60 : DO ispin = 1, nspin
474 40 : IF (use_subspace_fock) THEN
475 : CALL copy_dbcsr_to_fm( &
476 24 : scf_env%scf_subspace_buffer%combined_fock(ispin, 1)%matrix, scf_env%scf_work2)
477 : ELSE
478 16 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, scf_env%scf_work2)
479 : END IF
480 : CALL cp_fm_scale_and_add(diis_weight, scf_env%scf_work1(ispin), &
481 60 : 1.0_dp - diis_weight, scf_env%scf_work2)
482 : END DO
483 : END IF
484 48 : diis_step = cdiis_step .AND. diis_weight >= 1.0_dp
485 48 : IF (diis_step) diis_error = cdiis_error
486 50 : ELSE IF (ASSOCIATED(scf_env%scf_diis_buffer)) THEN
487 50 : IF (scf_env%scf_diis_buffer%ncall > 0) CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
488 50 : diis_weight = 0.0_dp
489 : END IF
490 98 : scf_env%scf_subspace_buffer%diis_weight = diis_weight
491 : ELSE
492 99531 : IF (scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis) THEN
493 : CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
494 : scf_env%scf_work2, scf_env%iter_delta, diis_error, diis_step, &
495 : eps_diis, scf_control%nmixing, &
496 : s_matrix=matrix_s, &
497 80951 : scf_section=scf_section)
498 : ELSE
499 18580 : diis_step = .FALSE.
500 : END IF
501 : END IF
502 :
503 : do_level_shift = ((scf_control%level_shift /= 0.0_dp) .AND. &
504 : ((scf_control%density_guess == core_guess) .OR. &
505 99629 : (scf_env%iter_count > 1)))
506 :
507 99629 : IF (.NOT. use_adiis .AND. (scf_env%iter_count > 1) .AND. &
508 : (scf_env%iter_delta < scf_control%diagonalization%eps_jacobi)) THEN
509 0 : use_jacobi = .TRUE.
510 : ELSE
511 99531 : use_jacobi = .FALSE.
512 : END IF
513 :
514 98 : IF (use_adiis .AND. .NOT. diis_step) THEN
515 70 : scf_env%iter_param = diis_error
516 70 : IF (scf_env%scf_subspace_buffer%last_restart) THEN
517 0 : scf_env%iter_method = "ADIIS/Rst."
518 : ELSE
519 70 : scf_env%iter_method = "ADIIS/Diag."
520 : END IF
521 99559 : ELSE IF (diis_step) THEN
522 49719 : scf_env%iter_param = diis_error
523 49719 : IF (use_jacobi) THEN
524 0 : scf_env%iter_method = "DIIS/Jacobi"
525 : ELSE
526 49719 : scf_env%iter_method = "DIIS/Diag."
527 : END IF
528 : ELSE
529 49840 : IF (scf_env%mixing_method == 0) THEN
530 0 : scf_env%iter_method = "NoMix/Diag."
531 49840 : ELSE IF (scf_env%mixing_method == 1) THEN
532 45016 : scf_env%iter_param = scf_env%p_mix_alpha
533 45016 : IF (use_jacobi) THEN
534 0 : scf_env%iter_method = "P_Mix/Jacobi"
535 : ELSE
536 45016 : scf_env%iter_method = "P_Mix/Diag."
537 : END IF
538 4824 : ELSE IF (scf_env%mixing_method > 1) THEN
539 4824 : scf_env%iter_param = scf_env%mixing_store%alpha
540 4824 : IF (use_jacobi) THEN
541 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Jacobi"
542 : ELSE
543 4824 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Diag."
544 : END IF
545 : END IF
546 : END IF
547 :
548 99629 : IF (scf_env%cholesky_method == cholesky_dbcsr .AND. .NOT. use_cdiis) THEN
549 1072 : ortho_dbcsr => scf_env%ortho_dbcsr
550 3198 : DO ispin = 1, nspin
551 2126 : effective_ks_dbcsr => matrix_ks(ispin)%matrix
552 2126 : IF (use_subspace_fock) THEN
553 0 : effective_ks_dbcsr => scf_env%scf_subspace_buffer%combined_fock(ispin, 1)%matrix
554 : END IF
555 : CALL eigensolver_dbcsr(matrix_ks=effective_ks_dbcsr, matrix_ks_fm=scf_env%scf_work1(ispin), &
556 : mo_set=mos(ispin), &
557 : ortho_dbcsr=ortho_dbcsr, &
558 3198 : ksbuf1=scf_env%buf1_dbcsr, ksbuf2=scf_env%buf2_dbcsr)
559 : END DO
560 :
561 98557 : ELSE IF (scf_env%cholesky_method > cholesky_off) THEN
562 98123 : cdiis_cholesky_method = scf_env%cholesky_method
563 98123 : IF (cdiis_cholesky_method == cholesky_dbcsr) THEN
564 0 : cdiis_cholesky_method = cholesky_inverse
565 : END IF
566 98123 : IF (scf_env%cholesky_method == cholesky_inverse) THEN
567 264 : ortho => scf_env%ortho_m1
568 : ELSE
569 97859 : ortho => scf_env%ortho
570 : END IF
571 :
572 98123 : owns_ortho = .FALSE.
573 98123 : IF (.NOT. ASSOCIATED(ortho)) THEN
574 0 : ALLOCATE (ortho)
575 0 : owns_ortho = .TRUE.
576 : END IF
577 :
578 212450 : DO ispin = 1, nspin
579 212450 : IF (direct_generalized_diagonalization .AND. .NOT. do_level_shift) THEN
580 : CALL eigensolver_generalized(matrix_ks_fm=scf_env%scf_work1(ispin), &
581 : matrix_s=matrix_s(ispin)%matrix, &
582 : mo_set=mos(ispin), &
583 0 : work=scf_env%scf_work2)
584 : ELSE
585 114327 : IF (do_level_shift) THEN
586 : CALL eigensolver(matrix_ks_fm=scf_env%scf_work1(ispin), &
587 : mo_set=mos(ispin), &
588 : ortho=ortho, &
589 : work=scf_env%scf_work2, &
590 : cholesky_method=cdiis_cholesky_method, &
591 : do_level_shift=do_level_shift, &
592 : level_shift=scf_control%level_shift, &
593 : matrix_u_fm=scf_env%ortho, &
594 124 : use_jacobi=use_jacobi)
595 : ELSE
596 : CALL eigensolver(matrix_ks_fm=scf_env%scf_work1(ispin), &
597 : mo_set=mos(ispin), &
598 : ortho=ortho, &
599 : work=scf_env%scf_work2, &
600 : cholesky_method=cdiis_cholesky_method, &
601 : do_level_shift=do_level_shift, &
602 : level_shift=scf_control%level_shift, &
603 114203 : use_jacobi=use_jacobi)
604 : END IF
605 : END IF
606 : END DO
607 :
608 98123 : IF (scf_env%cholesky_method /= cholesky_dbcsr) THEN
609 98123 : scf_env%cholesky_method = cdiis_cholesky_method
610 : END IF
611 :
612 98123 : IF (owns_ortho) DEALLOCATE (ortho)
613 : ELSE
614 434 : ortho => scf_env%ortho
615 :
616 434 : owns_ortho = .FALSE.
617 434 : IF (.NOT. ASSOCIATED(ortho)) THEN
618 0 : ALLOCATE (ortho)
619 0 : owns_ortho = .TRUE.
620 : END IF
621 :
622 434 : IF (do_level_shift) THEN
623 172 : DO ispin = 1, nspin
624 : IF (ASSOCIATED(scf_env%scf_work1_red) .AND. ASSOCIATED(scf_env%scf_work2_red) &
625 172 : .AND. ASSOCIATED(scf_env%ortho_red) .AND. ASSOCIATED(scf_env%ortho_m1_red)) THEN
626 : CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
627 : mo_set=mos(ispin), &
628 : ortho=ortho, &
629 : work=scf_env%scf_work2, &
630 : do_level_shift=do_level_shift, &
631 : level_shift=scf_control%level_shift, &
632 : matrix_u_fm=scf_env%ortho_m1, &
633 : use_jacobi=use_jacobi, &
634 : jacobi_threshold=scf_control%diagonalization%jacobi_threshold, &
635 : matrix_ks_fm_red=scf_env%scf_work1_red(ispin), &
636 : ortho_red=scf_env%ortho_red, &
637 : work_red=scf_env%scf_work2_red, &
638 86 : matrix_u_fm_red=scf_env%ortho_m1_red)
639 : ELSE
640 : CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
641 : mo_set=mos(ispin), &
642 : ortho=ortho, &
643 : work=scf_env%scf_work2, &
644 : do_level_shift=do_level_shift, &
645 : level_shift=scf_control%level_shift, &
646 : matrix_u_fm=scf_env%ortho_m1, &
647 : use_jacobi=use_jacobi, &
648 0 : jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
649 : END IF
650 : END DO
651 : ELSE
652 778 : DO ispin = 1, nspin
653 : IF (ASSOCIATED(scf_env%scf_work1_red) .AND. ASSOCIATED(scf_env%scf_work2_red) &
654 778 : .AND. ASSOCIATED(scf_env%ortho_red)) THEN
655 : CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
656 : mo_set=mos(ispin), &
657 : ortho=ortho, &
658 : work=scf_env%scf_work2, &
659 : do_level_shift=do_level_shift, &
660 : level_shift=scf_control%level_shift, &
661 : use_jacobi=use_jacobi, &
662 : jacobi_threshold=scf_control%diagonalization%jacobi_threshold, &
663 : matrix_ks_fm_red=scf_env%scf_work1_red(ispin), &
664 : ortho_red=scf_env%ortho_red, &
665 430 : work_red=scf_env%scf_work2_red)
666 : ELSE
667 : CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
668 : mo_set=mos(ispin), &
669 : ortho=ortho, &
670 : work=scf_env%scf_work2, &
671 : do_level_shift=do_level_shift, &
672 : level_shift=scf_control%level_shift, &
673 : use_jacobi=use_jacobi, &
674 0 : jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
675 : END IF
676 : END DO
677 : END IF
678 :
679 434 : IF (owns_ortho) DEALLOCATE (ortho)
680 : END IF
681 :
682 99629 : END SUBROUTINE general_eigenproblem
683 :
684 : ! **************************************************************************************************
685 : !> \brief ...
686 : !> \param scf_env ...
687 : !> \param mos ...
688 : !> \param matrix_ks ...
689 : !> \param matrix_s ...
690 : !> \param scf_control ...
691 : !> \param scf_section ...
692 : !> \param diis_step ...
693 : !> \param probe ...
694 : ! **************************************************************************************************
695 98589 : SUBROUTINE do_general_diag(scf_env, mos, matrix_ks, &
696 : matrix_s, scf_control, scf_section, &
697 : diis_step, probe)
698 :
699 : TYPE(qs_scf_env_type), POINTER :: scf_env
700 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
701 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
702 : TYPE(scf_control_type), POINTER :: scf_control
703 : TYPE(section_vals_type), POINTER :: scf_section
704 : LOGICAL, INTENT(INOUT) :: diis_step
705 : TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
706 : POINTER :: probe
707 :
708 : INTEGER :: ispin, nspin
709 : REAL(KIND=dp) :: total_zeff_corr
710 :
711 98589 : nspin = SIZE(matrix_ks)
712 :
713 : CALL general_eigenproblem(scf_env, mos, matrix_ks, &
714 98589 : matrix_s, scf_control, scf_section, diis_step)
715 :
716 : total_zeff_corr = 0.0_dp
717 98589 : total_zeff_corr = scf_env%sum_zeff_corr
718 :
719 98589 : IF (ABS(total_zeff_corr) > 0.0_dp) THEN
720 : CALL set_mo_occupation(mo_array=mos, &
721 40 : smear=scf_control%smear, tot_zeff_corr=total_zeff_corr)
722 : ELSE
723 98549 : IF (PRESENT(probe) .EQV. .TRUE.) THEN
724 14 : scf_control%smear%do_smear = .FALSE.
725 : CALL set_mo_occupation(mo_array=mos, &
726 : smear=scf_control%smear, &
727 14 : probe=probe)
728 : ELSE
729 98535 : IF (.NOT. scf_control%gce%do_gce) THEN
730 : CALL set_mo_occupation(mo_array=mos, &
731 98471 : smear=scf_control%smear)
732 : ELSE
733 : CALL set_mo_occupation(mo_array=mos, &
734 : smear=scf_control%smear, &
735 64 : gce=scf_control%gce)
736 : END IF
737 : END IF
738 : END IF
739 :
740 213478 : DO ispin = 1, nspin
741 : CALL calculate_density_matrix(mos(ispin), &
742 213478 : scf_env%p_mix_new(ispin, 1)%matrix)
743 : END DO
744 :
745 98589 : END SUBROUTINE do_general_diag
746 :
747 : ! **************************************************************************************************
748 : !> \brief Kpoint diagonalization routine
749 : !> Transforms matrices to kpoint, distributes kpoint groups, performs
750 : !> general diagonalization (no storgae of overlap decomposition), stores
751 : !> MOs, calculates occupation numbers, calculates density matrices
752 : !> in kpoint representation, transforms density matrices to real space
753 : !> \param matrix_ks Kohn-sham matrices (RS indices, global)
754 : !> \param matrix_s Overlap matrices (RS indices, global)
755 : !> \param kpoints Kpoint environment
756 : !> \param scf_env SCF environment
757 : !> \param scf_control SCF control variables
758 : !> \param update_p ...
759 : !> \param diis_step ...
760 : !> \param diis_error ...
761 : !> \param qs_env ...
762 : !> \param probe ...
763 : !> \param matrix_t optional kinetic operator retained for k-point OT preconditioning
764 : !> \param store_ot_matrices retain H(k) and S(k) in the local k-point-group context for OT
765 : !> \param transfer_only redistribute and retain OT operators without diagonalization or density update
766 : !> \param added_mos_auto_grow true when the adaptive virtual-space buffer must grow
767 : !> \par History
768 : !> 08.2014 created [JGH]
769 : ! **************************************************************************************************
770 34262 : SUBROUTINE do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, update_p, &
771 : diis_step, diis_error, qs_env, probe, matrix_t, &
772 : store_ot_matrices, transfer_only, added_mos_auto_grow)
773 :
774 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
775 : TYPE(kpoint_type), POINTER :: kpoints
776 : TYPE(qs_scf_env_type), POINTER :: scf_env
777 : TYPE(scf_control_type), POINTER :: scf_control
778 : LOGICAL, INTENT(IN) :: update_p
779 : LOGICAL, INTENT(INOUT) :: diis_step
780 : REAL(dp), INTENT(INOUT), OPTIONAL :: diis_error
781 : TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
782 : TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
783 : POINTER :: probe
784 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
785 : POINTER :: matrix_t
786 : LOGICAL, INTENT(IN), OPTIONAL :: store_ot_matrices, transfer_only
787 : LOGICAL, INTENT(OUT), OPTIONAL :: added_mos_auto_grow
788 :
789 : CHARACTER(len=*), PARAMETER :: routineN = 'do_general_diag_kp'
790 :
791 34262 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coeffs
792 : INTEGER :: cdiis_ncall, handle, ib, ic, igroup, ik, ikp, indx, ispin, jb, kplocal, &
793 : lattice_fft, nb, ninfo, nkp, nkp_groups, nspin, omp_threads, ot_components
794 : INTEGER, DIMENSION(2) :: kp_range
795 : INTEGER, DIMENSION(3) :: nkp_grid
796 34262 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
797 34262 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
798 : LOGICAL :: adiis_cdiis, cdiis_step, do_diis, my_added_mos_auto_grow, my_kpgrp, &
799 : my_store_ot_kinetic, my_store_ot_matrices, my_transfer_only, spin_grid_fft, use_adiis, &
800 : use_grid_fft, use_omp_kpoints, use_real_wfn
801 : REAL(KIND=dp) :: diis_weight
802 34262 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
803 34262 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
804 34262 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
805 411144 : TYPE(copy_info_type), DIMENSION(2) :: blend_info
806 : TYPE(cp_cfm_type) :: cksmat, cmos, csmat, ctmat, cwork
807 34262 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools
808 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct, mo_struct
809 : TYPE(cp_fm_type) :: fmdummy, fmlocal, rksmat, rsmat, rtmat
810 34262 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
811 : TYPE(cp_fm_type), POINTER :: imos, mo_coeff, rmos
812 68524 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: adiis_ks, effective_ks
813 : TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix, tmpmat
814 : TYPE(kpoint_env_type), POINTER :: kp
815 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_global
816 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
817 34262 : POINTER :: sab_nl
818 : TYPE(qs_matrix_pools_type), POINTER :: mpools
819 34262 : TYPE(rskp_grid_type) :: s_grid
820 34262 : TYPE(rskp_grid_type), ALLOCATABLE, DIMENSION(:) :: ks_grid
821 : TYPE(section_vals_type), POINTER :: scf_section
822 :
823 34262 : CALL timeset(routineN, handle)
824 :
825 34262 : my_store_ot_matrices = .FALSE.
826 34262 : IF (PRESENT(store_ot_matrices)) my_store_ot_matrices = store_ot_matrices
827 34262 : my_transfer_only = .FALSE.
828 34262 : IF (PRESENT(transfer_only)) my_transfer_only = transfer_only
829 1110 : IF (my_transfer_only) my_store_ot_matrices = .TRUE.
830 34262 : my_store_ot_kinetic = PRESENT(matrix_t)
831 34262 : IF (my_store_ot_kinetic) my_store_ot_kinetic = ASSOCIATED(matrix_t)
832 34262 : CPASSERT(.NOT. my_transfer_only .OR. .NOT. update_p)
833 34262 : CPASSERT(.NOT. my_store_ot_kinetic .OR. my_transfer_only)
834 34262 : my_added_mos_auto_grow = .FALSE.
835 34262 : IF (PRESENT(added_mos_auto_grow)) added_mos_auto_grow = .FALSE.
836 :
837 34262 : NULLIFY (adiis_ks, effective_ks, sab_nl)
838 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
839 : nkp_groups=nkp_groups, kp_dist=kp_dist, sab_nl=sab_nl, &
840 34262 : cell_to_index=cell_to_index, nkp_grid=nkp_grid, lattice_fft=lattice_fft)
841 34262 : CPASSERT(ASSOCIATED(sab_nl))
842 34262 : kplocal = kp_range(2) - kp_range(1) + 1
843 :
844 : !Whether we use DIIS for k-points
845 34262 : do_diis = .FALSE.
846 34262 : use_adiis = update_p .AND. scf_control%diagonalization%update_method == diag_update_method_adiis
847 : IF (use_adiis) THEN
848 122 : CPASSERT(ASSOCIATED(scf_env%scf_subspace_buffer))
849 : END IF
850 122 : adiis_cdiis = .FALSE.
851 122 : IF (use_adiis .AND. ASSOCIATED(kpoints%scf_diis_buffer)) THEN
852 : IF (scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis .AND. &
853 : .NOT. scf_env%adiis_check_next .AND. &
854 : scf_env%scf_subspace_buffer%diis_state_valid .AND. &
855 : scf_env%scf_subspace_buffer%nstored > 0 .AND. &
856 122 : PRESENT(diis_error) .AND. PRESENT(qs_env)) THEN
857 : adiis_cdiis = scf_env%iter_delta < scf_control%eps_diis .OR. &
858 108 : kpoints%scf_diis_buffer%ncall > 0
859 : END IF
860 : END IF
861 : IF ((.NOT. use_adiis .AND. scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis .AND. &
862 34262 : .NOT. use_real_wfn .AND. PRESENT(diis_error) .AND. PRESENT(qs_env)) .OR. &
863 27126 : adiis_cdiis) do_diis = .TRUE.
864 34262 : IF (my_transfer_only) do_diis = .FALSE.
865 :
866 34262 : adiis_ks => matrix_ks
867 34262 : IF (use_adiis) THEN
868 122 : IF (scf_env%scf_subspace_buffer%use_combined_fock) THEN
869 50 : adiis_ks => scf_env%scf_subspace_buffer%combined_fock
870 : END IF
871 : END IF
872 34262 : effective_ks => matrix_ks
873 34262 : IF (use_adiis .AND. .NOT. adiis_cdiis) THEN
874 50 : effective_ks => adiis_ks
875 : END IF
876 34262 : IF (use_adiis .AND. .NOT. do_diis) THEN
877 50 : scf_env%scf_subspace_buffer%diis_weight = 0.0_dp
878 50 : IF (ASSOCIATED(kpoints%scf_diis_buffer)) THEN
879 50 : IF (kpoints%scf_diis_buffer%ncall > 0) CALL qs_diis_b_clear_kp(kpoints%scf_diis_buffer)
880 : END IF
881 : END IF
882 :
883 34262 : para_env => kpoints%blacs_env_all%para_env
884 34262 : nspin = SIZE(matrix_ks, 1)
885 34262 : omp_threads = 1
886 34262 : !$ omp_threads = omp_get_max_threads()
887 : ! The local executor operates on the complex k-point representation.
888 : use_omp_kpoints = para_env%num_pe == 1 .AND. omp_threads > 1 .AND. kplocal > 1 .AND. &
889 34262 : .NOT. use_real_wfn .AND. .NOT. my_store_ot_matrices
890 : ! Enter the threaded executor when the ADIIS/CDIIS handover starts, before
891 : ! its DIIS buffer is allocated with an executor-specific matrix structure.
892 34262 : IF (use_adiis) THEN
893 122 : use_omp_kpoints = use_omp_kpoints .AND. adiis_cdiis
894 : END IF
895 : ! Keep accelerator and distributed eigensolver selections on their existing paths.
896 34262 : IF (diag_type == FM_DIAG_TYPE_CUSOLVER .OR. diag_type == FM_DIAG_TYPE_DLAF) THEN
897 0 : use_omp_kpoints = .FALSE.
898 : END IF
899 34262 : IF (diag_check_requested()) use_omp_kpoints = .FALSE.
900 :
901 34208 : IF (use_omp_kpoints) THEN
902 : CALL do_general_diag_kp_omp(effective_ks, matrix_s, adiis_ks, kpoints, scf_env, scf_control, update_p, &
903 : diis_step, do_diis, qs_env, diis_error, probe, sab_nl, &
904 0 : cell_to_index, xkp)
905 0 : CALL timestop(handle)
906 0 : RETURN
907 : END IF
908 :
909 : ! allocate some work matrices
910 34262 : ALLOCATE (rmatrix, cmatrix, tmpmat)
911 : CALL dbcsr_create(rmatrix, template=matrix_ks(1, 1)%matrix, &
912 34262 : matrix_type=dbcsr_type_symmetric)
913 : CALL dbcsr_create(cmatrix, template=matrix_ks(1, 1)%matrix, &
914 34262 : matrix_type=dbcsr_type_antisymmetric)
915 : CALL dbcsr_create(tmpmat, template=matrix_ks(1, 1)%matrix, &
916 34262 : matrix_type=dbcsr_type_no_symmetry)
917 34262 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
918 34262 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
919 :
920 34262 : fmwork => scf_env%scf_work1
921 :
922 : ! fm pools to be used within a kpoint group
923 34262 : CALL get_kpoint_info(kpoints, mpools=mpools)
924 34262 : CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
925 :
926 34262 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
927 34262 : CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
928 :
929 34262 : IF (use_real_wfn) THEN
930 130 : CALL cp_fm_create(rksmat, matrix_struct)
931 130 : CALL cp_fm_create(rsmat, matrix_struct)
932 130 : IF (my_store_ot_kinetic) CALL cp_fm_create(rtmat, matrix_struct)
933 : END IF
934 34262 : IF (.NOT. use_real_wfn .OR. do_diis) THEN
935 34132 : CALL cp_cfm_create(cksmat, matrix_struct)
936 34132 : CALL cp_cfm_create(csmat, matrix_struct)
937 34132 : CALL cp_cfm_create(cwork, matrix_struct)
938 34132 : IF (my_store_ot_kinetic) CALL cp_cfm_create(ctmat, matrix_struct)
939 34132 : IF (.NOT. use_real_wfn) THEN
940 34132 : kp => kpoints%kp_env(1)%kpoint_env
941 34132 : CALL get_mo_set(kp%mos(1, 1), mo_coeff=mo_coeff)
942 34132 : CALL cp_fm_get_info(mo_coeff, matrix_struct=mo_struct)
943 34132 : CALL cp_cfm_create(cmos, mo_struct)
944 : END IF
945 : END IF
946 :
947 34262 : para_env => kpoints%blacs_env_all%para_env
948 34262 : nspin = SIZE(matrix_ks, 1)
949 :
950 : ! For a complete Monkhorst-Pack grid, transform every MPI-local DBCSR block over all k points
951 : ! in one lattice FFT. The established per-k phase sum remains the exact memory/irregular-grid
952 : ! fallback. Divide the conservative memory budget between H spins and S.
953 138956 : ALLOCATE (ks_grid(nspin))
954 34254 : SELECT CASE (lattice_fft)
955 : CASE (lattice_fft_auto)
956 34254 : use_grid_fft = nkp >= 27
957 : CASE (lattice_fft_on)
958 8 : use_grid_fft = .TRUE.
959 : CASE (lattice_fft_off)
960 0 : use_grid_fft = .FALSE.
961 : CASE DEFAULT
962 34262 : CPABORT("Unknown lattice FFT selection mode")
963 : END SELECT
964 34262 : IF (use_grid_fft) THEN
965 : CALL rskp_transform_grid_prepare(s_grid, rmatrix, matrix_s, 1, xkp, nkp_grid, &
966 : cell_to_index, sab_nl, use_grid_fft, &
967 212 : max_storage_bytes=512_int_8*1024_int_8**2/INT(nspin + 1, int_8))
968 : END IF
969 34262 : IF (use_grid_fft) THEN
970 8 : DO ispin = 1, nspin
971 : CALL rskp_transform_grid_prepare(ks_grid(ispin), rmatrix, effective_ks, ispin, xkp, nkp_grid, &
972 : cell_to_index, sab_nl, spin_grid_fft, &
973 4 : max_storage_bytes=512_int_8*1024_int_8**2/INT(nspin + 1, int_8))
974 8 : use_grid_fft = use_grid_fft .AND. spin_grid_fft
975 : END DO
976 : END IF
977 34262 : IF (.NOT. use_grid_fft) THEN
978 34258 : CALL rskp_transform_grid_release(s_grid)
979 70424 : DO ispin = 1, nspin
980 70424 : CALL rskp_transform_grid_release(ks_grid(ispin))
981 : END DO
982 : END IF
983 :
984 34262 : IF (my_store_ot_matrices) THEN
985 1110 : ot_components = MERGE(1, 2, use_real_wfn)
986 2850 : DO ikp = 1, kplocal
987 1740 : kp => kpoints%kp_env(ikp)%kpoint_env
988 1740 : CALL cp_fm_release(kp%ot_hmat)
989 1740 : CALL cp_fm_release(kp%ot_smat)
990 17580 : ALLOCATE (kp%ot_hmat(ot_components, nspin), kp%ot_smat(ot_components))
991 5220 : DO ic = 1, ot_components
992 3480 : CALL cp_fm_create(kp%ot_smat(ic), matrix_struct)
993 8820 : DO ispin = 1, nspin
994 7080 : CALL cp_fm_create(kp%ot_hmat(ic, ispin), matrix_struct)
995 : END DO
996 : END DO
997 2850 : IF (my_store_ot_kinetic) THEN
998 192 : CALL cp_fm_release(kp%ot_tmat)
999 768 : ALLOCATE (kp%ot_tmat(ot_components))
1000 576 : DO ic = 1, ot_components
1001 576 : CALL cp_fm_create(kp%ot_tmat(ic), matrix_struct)
1002 : END DO
1003 : END IF
1004 : END DO
1005 : END IF
1006 :
1007 34262 : ninfo = 4
1008 34262 : IF (my_store_ot_kinetic) ninfo = 6
1009 921452 : ALLOCATE (info(kplocal*nspin*nkp_groups, ninfo))
1010 :
1011 : ! Setup and start all the communication
1012 34262 : indx = 0
1013 97400 : DO ikp = 1, kplocal
1014 166572 : DO ispin = 1, nspin
1015 242578 : DO igroup = 1, nkp_groups
1016 : ! number of current kpoint
1017 110268 : ik = kp_dist(1, igroup) + ikp - 1
1018 110268 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
1019 110268 : indx = indx + 1
1020 110268 : IF (use_real_wfn) THEN
1021 : ! FT of matrices KS and S, then transfer to FM type
1022 156 : IF (use_grid_fft) THEN
1023 0 : CALL rskp_transform_grid_extract(ks_grid(ispin), ik, rmatrix)
1024 : ELSE
1025 156 : CALL dbcsr_set(rmatrix, 0.0_dp)
1026 : CALL rskp_transform(rmatrix=rmatrix, rsmat=effective_ks, ispin=ispin, &
1027 156 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1028 : END IF
1029 156 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1030 156 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
1031 : ! s matrix is not spin dependent
1032 156 : IF (use_grid_fft) THEN
1033 0 : CALL rskp_transform_grid_extract(s_grid, ik, rmatrix)
1034 : ELSE
1035 156 : CALL dbcsr_set(rmatrix, 0.0_dp)
1036 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_s, ispin=1, &
1037 156 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1038 : END IF
1039 156 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1040 156 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(3))
1041 : ELSE
1042 : ! FT of matrices KS and S, then transfer to FM type
1043 110112 : IF (use_grid_fft) THEN
1044 256 : CALL rskp_transform_grid_extract(ks_grid(ispin), ik, rmatrix, cmatrix)
1045 : ELSE
1046 109856 : CALL dbcsr_set(rmatrix, 0.0_dp)
1047 109856 : CALL dbcsr_set(cmatrix, 0.0_dp)
1048 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=effective_ks, ispin=ispin, &
1049 109856 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1050 : END IF
1051 110112 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1052 110112 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
1053 110112 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
1054 110112 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(2))
1055 : ! s matrix is not spin dependent, double the work
1056 110112 : IF (use_grid_fft) THEN
1057 256 : CALL rskp_transform_grid_extract(s_grid, ik, rmatrix, cmatrix)
1058 : ELSE
1059 109856 : CALL dbcsr_set(rmatrix, 0.0_dp)
1060 109856 : CALL dbcsr_set(cmatrix, 0.0_dp)
1061 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_s, ispin=1, &
1062 109856 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1063 : END IF
1064 110112 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1065 110112 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(3))
1066 110112 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
1067 110112 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(4))
1068 : END IF
1069 : ! transfer to kpoint group
1070 : ! redistribution of matrices, new blacs environment
1071 : ! fmwork -> fmlocal -> rksmat/cksmat
1072 : ! fmwork -> fmlocal -> rsmat/csmat
1073 110268 : IF (use_real_wfn) THEN
1074 156 : IF (my_kpgrp) THEN
1075 156 : CALL cp_fm_start_copy_general(fmwork(1), rksmat, para_env, info(indx, 1))
1076 156 : CALL cp_fm_start_copy_general(fmwork(3), rsmat, para_env, info(indx, 2))
1077 : ELSE
1078 0 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
1079 0 : CALL cp_fm_start_copy_general(fmwork(3), fmdummy, para_env, info(indx, 2))
1080 : END IF
1081 : ELSE
1082 110112 : IF (my_kpgrp) THEN
1083 69016 : CALL cp_fm_start_copy_general(fmwork(1), fmlocal, para_env, info(indx, 1))
1084 69016 : CALL cp_fm_start_copy_general(fmwork(2), fmlocal, para_env, info(indx, 2))
1085 69016 : CALL cp_fm_start_copy_general(fmwork(3), fmlocal, para_env, info(indx, 3))
1086 69016 : CALL cp_fm_start_copy_general(fmwork(4), fmlocal, para_env, info(indx, 4))
1087 : ELSE
1088 41096 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
1089 41096 : CALL cp_fm_start_copy_general(fmwork(2), fmdummy, para_env, info(indx, 2))
1090 41096 : CALL cp_fm_start_copy_general(fmwork(3), fmdummy, para_env, info(indx, 3))
1091 41096 : CALL cp_fm_start_copy_general(fmwork(4), fmdummy, para_env, info(indx, 4))
1092 : END IF
1093 : END IF
1094 179440 : IF (my_store_ot_kinetic) THEN
1095 240 : CALL dbcsr_set(rmatrix, 0.0_dp)
1096 240 : IF (use_real_wfn) THEN
1097 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_t, ispin=1, &
1098 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1099 : ELSE
1100 240 : CALL dbcsr_set(cmatrix, 0.0_dp)
1101 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_t, ispin=1, &
1102 240 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1103 : END IF
1104 240 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1105 240 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(3))
1106 240 : IF (use_real_wfn) THEN
1107 0 : IF (my_kpgrp) THEN
1108 0 : CALL cp_fm_start_copy_general(fmwork(3), rtmat, para_env, info(indx, 3))
1109 : ELSE
1110 0 : CALL cp_fm_start_copy_general(fmwork(3), fmdummy, para_env, info(indx, 3))
1111 : END IF
1112 : ELSE
1113 240 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
1114 240 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(4))
1115 240 : IF (my_kpgrp) THEN
1116 232 : CALL cp_fm_start_copy_general(fmwork(3), fmlocal, para_env, info(indx, 5))
1117 232 : CALL cp_fm_start_copy_general(fmwork(4), fmlocal, para_env, info(indx, 6))
1118 : ELSE
1119 8 : CALL cp_fm_start_copy_general(fmwork(3), fmdummy, para_env, info(indx, 5))
1120 8 : CALL cp_fm_start_copy_general(fmwork(4), fmdummy, para_env, info(indx, 6))
1121 : END IF
1122 : END IF
1123 : END IF
1124 : END DO
1125 : END DO
1126 : END DO
1127 :
1128 : ! Finish communication then diagonalise in each group
1129 34262 : IF (do_diis) THEN
1130 27126 : CALL get_qs_env(qs_env, para_env=para_env_global)
1131 27126 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
1132 27126 : cdiis_ncall = kpoints%scf_diis_buffer%ncall
1133 27126 : CALL qs_diis_b_info_kp(kpoints%scf_diis_buffer, ib, nb)
1134 27126 : indx = 0
1135 71066 : DO ikp = 1, kplocal
1136 120054 : DO ispin = 1, nspin
1137 129132 : DO igroup = 1, nkp_groups
1138 : ! number of current kpoint
1139 80144 : ik = kp_dist(1, igroup) + ikp - 1
1140 80144 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
1141 31156 : indx = indx + 1
1142 48988 : IF (my_kpgrp) THEN
1143 48988 : IF (use_real_wfn) THEN
1144 0 : CALL cp_fm_finish_copy_general(rksmat, info(indx, 1))
1145 0 : CALL cp_fm_finish_copy_general(rsmat, info(indx, 2))
1146 0 : IF (my_store_ot_kinetic) THEN
1147 0 : CALL cp_fm_finish_copy_general(rtmat, info(indx, 3))
1148 : END IF
1149 0 : CALL cp_cfm_scale_and_add_fm(z_zero, cksmat, z_one, rksmat)
1150 0 : CALL cp_cfm_scale_and_add_fm(z_zero, csmat, z_one, rsmat)
1151 : ELSE
1152 48988 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
1153 48988 : CALL cp_cfm_scale_and_add_fm(z_zero, cksmat, z_one, fmlocal)
1154 48988 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
1155 48988 : CALL cp_cfm_scale_and_add_fm(z_one, cksmat, gaussi, fmlocal)
1156 48988 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 3))
1157 48988 : CALL cp_cfm_scale_and_add_fm(z_zero, csmat, z_one, fmlocal)
1158 48988 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 4))
1159 48988 : CALL cp_cfm_scale_and_add_fm(z_one, csmat, gaussi, fmlocal)
1160 48988 : IF (my_store_ot_kinetic) THEN
1161 0 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 5))
1162 0 : CALL cp_cfm_scale_and_add_fm(z_zero, ctmat, z_one, fmlocal)
1163 0 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 6))
1164 0 : CALL cp_cfm_scale_and_add_fm(z_one, ctmat, gaussi, fmlocal)
1165 : END IF
1166 : END IF
1167 : END IF
1168 : END DO !igroup
1169 :
1170 48988 : kp => kpoints%kp_env(ikp)%kpoint_env
1171 : CALL qs_diis_b_calc_err_kp(kpoints%scf_diis_buffer, ib, kp%mos, cksmat, csmat, &
1172 92928 : ispin, ikp, kplocal, scf_section, real_wfn=use_real_wfn)
1173 :
1174 : END DO !ispin
1175 : END DO !ikp
1176 :
1177 81378 : ALLOCATE (coeffs(nb))
1178 : CALL qs_diis_b_step_kp(kpoints%scf_diis_buffer, coeffs, ib, nb, scf_env%iter_delta, diis_error, &
1179 : cdiis_step, scf_control%eps_diis, nspin, nkp, kplocal, &
1180 : MERGE(2, scf_control%nmixing, use_adiis), &
1181 54180 : scf_section, para_env_global)
1182 27126 : diis_weight = 1.0_dp
1183 27126 : IF (use_adiis) THEN
1184 72 : diis_weight = scf_env%scf_subspace_buffer%diis_weight
1185 72 : IF (.NOT. ieee_is_finite(diis_error) .OR. diis_error >= scf_control%eps_diis) THEN
1186 8 : cdiis_step = .FALSE.
1187 8 : IF (cdiis_ncall < kpoints%scf_diis_buffer%nbuffer) THEN
1188 : ! The rejected state has not overwritten a full circular buffer.
1189 8 : kpoints%scf_diis_buffer%ncall = cdiis_ncall
1190 : ELSE
1191 0 : CALL qs_diis_b_clear_kp(kpoints%scf_diis_buffer)
1192 0 : diis_weight = 0.0_dp
1193 : END IF
1194 : END IF
1195 72 : IF (cdiis_step) THEN
1196 : ! Use one half-weight handover step; MAX_DIIS only controls the buffer capacity.
1197 56 : diis_weight = MIN(1.0_dp, diis_weight + 0.5_dp)
1198 : ELSE
1199 : ! A rejected or deferred CDIIS step must not leave a partial raw-Fock blend active.
1200 : diis_weight = 0.0_dp
1201 : END IF
1202 72 : scf_env%scf_subspace_buffer%diis_weight = diis_weight
1203 72 : diis_step = cdiis_step .AND. diis_weight >= 1.0_dp
1204 72 : IF (.NOT. diis_step) diis_error = scf_env%scf_subspace_buffer%last_old_fock_weight
1205 : ELSE
1206 27054 : diis_step = cdiis_step
1207 : END IF
1208 :
1209 : !build the ks matrices and diagonalize
1210 98192 : DO ikp = 1, kplocal
1211 120054 : DO ispin = 1, nspin
1212 48988 : IF (use_adiis) THEN
1213 120 : IF (diis_weight < 1.0_dp) THEN
1214 108 : DO igroup = 1, nkp_groups
1215 72 : ik = kp_dist(1, igroup) + ikp - 1
1216 72 : my_kpgrp = ik >= kp_range(1) .AND. ik <= kp_range(2)
1217 72 : CALL dbcsr_set(rmatrix, 0.0_dp)
1218 72 : CALL dbcsr_set(cmatrix, 0.0_dp)
1219 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, &
1220 : rsmat=adiis_ks, ispin=ispin, &
1221 72 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1222 72 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1223 72 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
1224 72 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
1225 72 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(2))
1226 72 : IF (my_kpgrp) THEN
1227 36 : CALL cp_fm_start_copy_general(fmwork(1), fmlocal, para_env, blend_info(1))
1228 36 : CALL cp_fm_start_copy_general(fmwork(2), fmlocal, para_env, blend_info(2))
1229 36 : CALL cp_fm_finish_copy_general(fmlocal, blend_info(1))
1230 36 : CALL cp_cfm_scale_and_add_fm(z_zero, cwork, z_one, fmlocal)
1231 36 : CALL cp_fm_finish_copy_general(fmlocal, blend_info(2))
1232 36 : CALL cp_cfm_scale_and_add_fm(z_one, cwork, gaussi, fmlocal)
1233 : ELSE
1234 36 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, blend_info(1))
1235 36 : CALL cp_fm_start_copy_general(fmwork(2), fmdummy, para_env, blend_info(2))
1236 : END IF
1237 72 : CALL cp_fm_cleanup_copy_general(blend_info(1))
1238 108 : CALL cp_fm_cleanup_copy_general(blend_info(2))
1239 : END DO
1240 : END IF
1241 : END IF
1242 48988 : kp => kpoints%kp_env(ikp)%kpoint_env
1243 48988 : CALL cp_cfm_to_cfm(kpoints%scf_diis_buffer%smat(ikp), csmat)
1244 :
1245 48988 : CALL cp_cfm_set_all(cksmat, z_zero)
1246 212544 : DO jb = 1, nb
1247 212544 : CALL cp_cfm_scale_and_add(z_one, cksmat, coeffs(jb), kpoints%scf_diis_buffer%param(jb, ispin, ikp))
1248 : END DO
1249 48988 : IF (use_adiis) THEN
1250 120 : IF (diis_weight < 1.0_dp) THEN
1251 : CALL cp_cfm_scale_and_add(CMPLX(diis_weight, KIND=dp), cksmat, &
1252 36 : CMPLX(1.0_dp - diis_weight, KIND=dp), cwork)
1253 : END IF
1254 : END IF
1255 :
1256 92928 : IF (use_real_wfn) THEN
1257 0 : CALL cp_cfm_to_fm(cksmat, rksmat, fmlocal)
1258 0 : CALL cp_cfm_to_fm(csmat, rsmat, fmlocal)
1259 0 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=mo_coeff, eigenvalues=eigenvalues)
1260 0 : IF (scf_env%cholesky_method == cholesky_off) THEN
1261 : CALL cp_fm_geeig_canon(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal, &
1262 0 : scf_control%eps_eigval)
1263 : ELSE
1264 0 : CALL cp_fm_geeig(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal)
1265 : END IF
1266 : ELSE
1267 48988 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
1268 48988 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1269 48988 : IF (scf_env%cholesky_method == cholesky_off) THEN
1270 : CALL cp_cfm_geeig_canon(cksmat, csmat, cmos, eigenvalues, cwork, &
1271 50 : scf_control%eps_eigval)
1272 : ELSE
1273 48938 : CALL cp_cfm_geeig(cksmat, csmat, cmos, eigenvalues, cwork)
1274 : END IF
1275 : ! copy eigenvalues to imag set (keep them in sync)
1276 1704216 : kp%mos(2, ispin)%eigenvalues = eigenvalues
1277 : ! split real and imaginary part of mos
1278 48988 : CALL cp_cfm_to_fm(cmos, rmos, imos)
1279 : END IF
1280 : END DO
1281 : END DO
1282 :
1283 : ELSE !no DIIS
1284 7136 : IF (.NOT. update_p .OR. use_adiis) diis_step = .FALSE.
1285 7136 : indx = 0
1286 26334 : DO ikp = 1, kplocal
1287 46518 : DO ispin = 1, nspin
1288 50308 : DO igroup = 1, nkp_groups
1289 : ! number of current kpoint
1290 30124 : ik = kp_dist(1, igroup) + ikp - 1
1291 30124 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
1292 9940 : indx = indx + 1
1293 20184 : IF (my_kpgrp) THEN
1294 20184 : IF (use_real_wfn) THEN
1295 156 : CALL cp_fm_finish_copy_general(rksmat, info(indx, 1))
1296 156 : CALL cp_fm_finish_copy_general(rsmat, info(indx, 2))
1297 156 : IF (my_store_ot_kinetic) THEN
1298 0 : CALL cp_fm_finish_copy_general(rtmat, info(indx, 3))
1299 : END IF
1300 : ELSE
1301 20028 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
1302 20028 : CALL cp_cfm_scale_and_add_fm(z_zero, cksmat, z_one, fmlocal)
1303 20028 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
1304 20028 : CALL cp_cfm_scale_and_add_fm(z_one, cksmat, gaussi, fmlocal)
1305 20028 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 3))
1306 20028 : CALL cp_cfm_scale_and_add_fm(z_zero, csmat, z_one, fmlocal)
1307 20028 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 4))
1308 20028 : CALL cp_cfm_scale_and_add_fm(z_one, csmat, gaussi, fmlocal)
1309 20028 : IF (my_store_ot_kinetic) THEN
1310 232 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 5))
1311 232 : CALL cp_cfm_scale_and_add_fm(z_zero, ctmat, z_one, fmlocal)
1312 232 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 6))
1313 232 : CALL cp_cfm_scale_and_add_fm(z_one, ctmat, gaussi, fmlocal)
1314 : END IF
1315 : END IF
1316 : END IF
1317 : END DO
1318 :
1319 : ! Each kpoint group has now information on a kpoint to be diagonalized
1320 : ! General eigensolver Hermite or Symmetric
1321 20184 : kp => kpoints%kp_env(ikp)%kpoint_env
1322 39382 : IF (use_real_wfn) THEN
1323 156 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=mo_coeff, eigenvalues=eigenvalues)
1324 156 : IF (my_store_ot_matrices) THEN
1325 0 : CALL cp_fm_to_fm(rksmat, kp%ot_hmat(1, ispin))
1326 0 : IF (ispin == 1) CALL cp_fm_to_fm(rsmat, kp%ot_smat(1))
1327 0 : IF (ispin == 1 .AND. my_store_ot_kinetic) THEN
1328 0 : CALL cp_fm_to_fm(rtmat, kp%ot_tmat(1))
1329 : END IF
1330 : END IF
1331 156 : IF (.NOT. my_transfer_only) THEN
1332 156 : IF (scf_env%cholesky_method == cholesky_off) THEN
1333 : CALL cp_fm_geeig_canon(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal, &
1334 40 : scf_control%eps_eigval)
1335 : ELSE
1336 116 : CALL cp_fm_geeig(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal)
1337 : END IF
1338 : END IF
1339 : ELSE
1340 20028 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
1341 20028 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1342 20028 : IF (my_store_ot_matrices) THEN
1343 1800 : CALL cp_cfm_to_fm(cksmat, kp%ot_hmat(1, ispin), kp%ot_hmat(2, ispin))
1344 1800 : IF (ispin == 1) THEN
1345 1740 : CALL cp_cfm_to_fm(csmat, kp%ot_smat(1), kp%ot_smat(2))
1346 1740 : IF (my_store_ot_kinetic) THEN
1347 192 : CALL cp_cfm_to_fm(ctmat, kp%ot_tmat(1), kp%ot_tmat(2))
1348 : END IF
1349 : END IF
1350 : END IF
1351 20028 : IF (.NOT. my_transfer_only) THEN
1352 18228 : IF (scf_env%cholesky_method == cholesky_off) THEN
1353 : CALL cp_cfm_geeig_canon(cksmat, csmat, cmos, eigenvalues, cwork, &
1354 3914 : scf_control%eps_eigval)
1355 : ELSE
1356 14314 : CALL cp_cfm_geeig(cksmat, csmat, cmos, eigenvalues, cwork)
1357 : END IF
1358 : ! copy eigenvalues to imag set (keep them in sync)
1359 710576 : kp%mos(2, ispin)%eigenvalues = eigenvalues
1360 : ! split real and imaginary part of mos
1361 18228 : CALL cp_cfm_to_fm(cmos, rmos, imos)
1362 : END IF
1363 : END IF
1364 : END DO
1365 : END DO
1366 : END IF
1367 :
1368 : ! Clean up communication
1369 34262 : indx = 0
1370 97400 : DO ikp = 1, kplocal
1371 166572 : DO ispin = 1, nspin
1372 242578 : DO igroup = 1, nkp_groups
1373 : ! number of current kpoint
1374 110268 : ik = kp_dist(1, igroup) + ikp - 1
1375 110268 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
1376 110268 : indx = indx + 1
1377 179440 : IF (use_real_wfn) THEN
1378 156 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
1379 156 : CALL cp_fm_cleanup_copy_general(info(indx, 2))
1380 156 : IF (my_store_ot_kinetic) THEN
1381 0 : CALL cp_fm_cleanup_copy_general(info(indx, 3))
1382 : END IF
1383 : ELSE
1384 110112 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
1385 110112 : CALL cp_fm_cleanup_copy_general(info(indx, 2))
1386 110112 : CALL cp_fm_cleanup_copy_general(info(indx, 3))
1387 110112 : CALL cp_fm_cleanup_copy_general(info(indx, 4))
1388 110112 : IF (my_store_ot_kinetic) THEN
1389 240 : CALL cp_fm_cleanup_copy_general(info(indx, 5))
1390 240 : CALL cp_fm_cleanup_copy_general(info(indx, 6))
1391 : END IF
1392 : END IF
1393 : END DO
1394 : END DO
1395 : END DO
1396 :
1397 : ! All done
1398 475814 : DEALLOCATE (info)
1399 :
1400 34262 : IF (update_p .AND. .NOT. my_transfer_only) THEN
1401 : ! MO occupations
1402 33008 : IF (PRESENT(probe) .EQV. .TRUE.) THEN
1403 0 : scf_control%smear%do_smear = .FALSE.
1404 : CALL kpoint_set_mo_occupation(kpoints, scf_control%smear, &
1405 0 : probe=probe)
1406 : ELSE
1407 : CALL kpoint_set_mo_occupation(kpoints, scf_control%smear, &
1408 : added_mos_auto=scf_control%added_mos_auto, &
1409 33008 : added_mos_auto_grow=my_added_mos_auto_grow)
1410 : END IF
1411 33008 : IF (my_added_mos_auto_grow) THEN
1412 0 : IF (PRESENT(added_mos_auto_grow)) added_mos_auto_grow = .TRUE.
1413 : ELSE
1414 : ! density matrices
1415 33008 : CALL kpoint_density_matrices(kpoints)
1416 : ! density matrices in real space
1417 : CALL kpoint_density_transform(kpoints, scf_env%p_mix_new, .FALSE., &
1418 33008 : matrix_s(1, 1)%matrix, sab_nl, fmwork, overlap_rs=matrix_s)
1419 : END IF
1420 : END IF
1421 :
1422 34262 : CALL dbcsr_deallocate_matrix(rmatrix)
1423 34262 : CALL dbcsr_deallocate_matrix(cmatrix)
1424 34262 : CALL dbcsr_deallocate_matrix(tmpmat)
1425 34262 : CALL rskp_transform_grid_release(s_grid)
1426 70432 : DO ispin = 1, nspin
1427 70432 : CALL rskp_transform_grid_release(ks_grid(ispin))
1428 : END DO
1429 70432 : DEALLOCATE (ks_grid)
1430 :
1431 34262 : IF (use_real_wfn) THEN
1432 130 : CALL cp_fm_release(rksmat)
1433 130 : CALL cp_fm_release(rsmat)
1434 130 : IF (my_store_ot_kinetic) CALL cp_fm_release(rtmat)
1435 : END IF
1436 34262 : IF (.NOT. use_real_wfn .OR. do_diis) THEN
1437 34132 : CALL cp_cfm_release(cksmat)
1438 34132 : CALL cp_cfm_release(csmat)
1439 34132 : CALL cp_cfm_release(cwork)
1440 34132 : IF (my_store_ot_kinetic) CALL cp_cfm_release(ctmat)
1441 : END IF
1442 34262 : IF (.NOT. use_real_wfn) THEN
1443 34132 : CALL cp_cfm_release(cmos)
1444 : END IF
1445 34262 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
1446 :
1447 34262 : CALL timestop(handle)
1448 :
1449 137048 : END SUBROUTINE do_general_diag_kp
1450 :
1451 : ! **************************************************************************************************
1452 : !> \brief OpenMP executor for k-points whose BLACS grid contains one MPI rank.
1453 : !> Each worker owns all dense temporary matrices and reads the shared
1454 : !> DBCSR inputs without mutation. MPI-backed full-matrix operations are
1455 : !> replaced by local LAPACK calls.
1456 : !> \param matrix_ks ...
1457 : !> \param matrix_s ...
1458 : !> \param adiis_ks ...
1459 : !> \param kpoints ...
1460 : !> \param scf_env ...
1461 : !> \param scf_control ...
1462 : !> \param update_p ...
1463 : !> \param diis_step ...
1464 : !> \param do_diis ...
1465 : !> \param qs_env ...
1466 : !> \param diis_error ...
1467 : !> \param probe ...
1468 : !> \param sab_nl ...
1469 : !> \param cell_to_index ...
1470 : !> \param xkp ...
1471 : ! **************************************************************************************************
1472 0 : SUBROUTINE do_general_diag_kp_omp(matrix_ks, matrix_s, adiis_ks, kpoints, scf_env, scf_control, update_p, &
1473 : diis_step, do_diis, qs_env, diis_error, probe, sab_nl, &
1474 : cell_to_index, xkp)
1475 :
1476 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s, adiis_ks
1477 : TYPE(kpoint_type), POINTER :: kpoints
1478 : TYPE(qs_scf_env_type), POINTER :: scf_env
1479 : TYPE(scf_control_type), POINTER :: scf_control
1480 : LOGICAL, INTENT(IN) :: update_p
1481 : LOGICAL, INTENT(INOUT) :: diis_step
1482 : LOGICAL, INTENT(IN) :: do_diis
1483 : TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
1484 : REAL(dp), INTENT(INOUT), OPTIONAL :: diis_error
1485 : TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
1486 : POINTER :: probe
1487 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1488 : POINTER :: sab_nl
1489 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1490 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1491 :
1492 : CHARACTER(len=*), PARAMETER :: routineN = 'do_general_diag_kp_omp'
1493 :
1494 0 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coeffs
1495 : INTEGER :: cdiis_ncall, handle, ib, ikp, kplocal, &
1496 : nb, nspin, nworkers, thread
1497 : LOGICAL :: cdiis_step, use_adiis
1498 : REAL(KIND=dp) :: diis_weight
1499 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct, mo_struct
1500 : TYPE(cp_fm_type), POINTER :: mo_coeff
1501 : TYPE(kp_diag_omp_workspace_type), ALLOCATABLE, &
1502 0 : DIMENSION(:) :: workspace
1503 0 : TYPE(kp_transform_plan_type) :: transform_plan
1504 : TYPE(kpoint_env_type), POINTER :: kp
1505 : TYPE(mp_para_env_type), POINTER :: para_env_global
1506 : TYPE(section_vals_type), POINTER :: scf_section
1507 :
1508 0 : CALL timeset(routineN, handle)
1509 0 : NULLIFY (matrix_struct, mo_struct, mo_coeff, kp, para_env_global, scf_section)
1510 :
1511 0 : kplocal = SIZE(xkp, 2)
1512 0 : nspin = SIZE(matrix_ks, 1)
1513 : use_adiis = update_p .AND. &
1514 0 : scf_control%diagonalization%update_method == diag_update_method_adiis
1515 : nworkers = 1
1516 0 : !$ nworkers = MIN(omp_get_max_threads(), kplocal)
1517 0 : CALL cp_fm_get_info(scf_env%scf_work1(1), matrix_struct=matrix_struct)
1518 0 : kp => kpoints%kp_env(1)%kpoint_env
1519 0 : CALL get_mo_set(kp%mos(1, 1), mo_coeff=mo_coeff)
1520 0 : CALL cp_fm_get_info(mo_coeff, matrix_struct=mo_struct)
1521 0 : CALL kp_diag_omp_workspaces_create(workspace, nworkers, matrix_struct, mo_struct)
1522 : CALL kp_transform_plan_create(transform_plan, sab_nl, cell_to_index, SIZE(matrix_ks, 2), &
1523 0 : matrix_ks(1, 1)%matrix)
1524 :
1525 0 : IF (do_diis) THEN
1526 0 : CALL get_qs_env(qs_env, para_env=para_env_global)
1527 0 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
1528 0 : cdiis_ncall = kpoints%scf_diis_buffer%ncall
1529 0 : CALL qs_diis_b_info_kp(kpoints%scf_diis_buffer, ib, nb)
1530 : CALL qs_diis_b_check_i_alloc_kp(kpoints%scf_diis_buffer, matrix_struct, nspin, kplocal, &
1531 0 : scf_section)
1532 :
1533 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(DYNAMIC, 1) NUM_THREADS(nworkers) &
1534 : !$OMP SHARED(workspace,matrix_ks,matrix_s,kpoints,ib,transform_plan,xkp,kplocal,nspin,nworkers) &
1535 0 : !$OMP PRIVATE(ikp,thread)
1536 : DO ikp = 1, kplocal
1537 : thread = 1
1538 : !$ thread = omp_get_thread_num() + 1
1539 : CALL kp_diis_error_job_local(ikp, workspace(thread), matrix_ks, matrix_s, kpoints, &
1540 : ib, transform_plan, xkp, nspin)
1541 : END DO
1542 : !$OMP END PARALLEL DO
1543 :
1544 0 : ALLOCATE (coeffs(nb))
1545 : CALL qs_diis_b_step_kp(kpoints%scf_diis_buffer, coeffs, ib, nb, scf_env%iter_delta, diis_error, &
1546 : cdiis_step, scf_control%eps_diis, nspin, kplocal, kplocal, &
1547 : MERGE(2, scf_control%nmixing, use_adiis), &
1548 0 : scf_section, para_env_global)
1549 0 : diis_weight = 1.0_dp
1550 0 : IF (use_adiis) THEN
1551 0 : diis_weight = scf_env%scf_subspace_buffer%diis_weight
1552 0 : IF (.NOT. ieee_is_finite(diis_error) .OR. diis_error >= scf_control%eps_diis) THEN
1553 0 : cdiis_step = .FALSE.
1554 0 : IF (cdiis_ncall < kpoints%scf_diis_buffer%nbuffer) THEN
1555 0 : kpoints%scf_diis_buffer%ncall = cdiis_ncall
1556 : ELSE
1557 0 : CALL qs_diis_b_clear_kp(kpoints%scf_diis_buffer)
1558 0 : diis_weight = 0.0_dp
1559 : END IF
1560 : END IF
1561 0 : IF (cdiis_step) THEN
1562 0 : diis_weight = MIN(1.0_dp, diis_weight + 0.5_dp)
1563 : ELSE
1564 : ! A rejected or deferred CDIIS step must not leave a partial raw-Fock blend active.
1565 0 : diis_weight = 0.0_dp
1566 : END IF
1567 0 : scf_env%scf_subspace_buffer%diis_weight = diis_weight
1568 0 : diis_step = cdiis_step .AND. diis_weight >= 1.0_dp
1569 0 : IF (.NOT. diis_step) diis_error = scf_env%scf_subspace_buffer%last_old_fock_weight
1570 : ELSE
1571 0 : diis_step = cdiis_step
1572 : END IF
1573 :
1574 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(DYNAMIC, 1) NUM_THREADS(nworkers) &
1575 : !$OMP SHARED(workspace,kpoints,coeffs,nb,scf_env,scf_control,kplocal,nspin,nworkers, &
1576 0 : !$OMP adiis_ks,transform_plan,xkp,use_adiis,diis_weight) PRIVATE(ikp,thread)
1577 : DO ikp = 1, kplocal
1578 : thread = 1
1579 : !$ thread = omp_get_thread_num() + 1
1580 : CALL kp_diis_diag_job(ikp, workspace(thread), kpoints, coeffs, nb, scf_env, &
1581 : scf_control, nspin, adiis_ks, transform_plan, xkp, &
1582 : use_adiis, diis_weight)
1583 : END DO
1584 : !$OMP END PARALLEL DO
1585 0 : DEALLOCATE (coeffs)
1586 : ELSE
1587 0 : diis_step = .FALSE.
1588 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(DYNAMIC, 1) NUM_THREADS(nworkers) &
1589 : !$OMP SHARED(workspace,matrix_ks,matrix_s,kpoints,scf_env,scf_control,transform_plan,xkp,kplocal,nspin,nworkers) &
1590 0 : !$OMP PRIVATE(ikp,thread)
1591 : DO ikp = 1, kplocal
1592 : thread = 1
1593 : !$ thread = omp_get_thread_num() + 1
1594 : CALL kp_diag_job(ikp, workspace(thread), matrix_ks, matrix_s, kpoints, scf_env, scf_control, &
1595 : transform_plan, xkp, nspin)
1596 : END DO
1597 : !$OMP END PARALLEL DO
1598 : END IF
1599 :
1600 0 : IF (update_p) THEN
1601 0 : IF (PRESENT(probe)) THEN
1602 0 : scf_control%smear%do_smear = .FALSE.
1603 0 : CALL kpoint_set_mo_occupation(kpoints, scf_control%smear, probe=probe)
1604 : ELSE
1605 0 : CALL kpoint_set_mo_occupation(kpoints, scf_control%smear)
1606 : END IF
1607 0 : CALL kpoint_density_matrices(kpoints)
1608 : CALL kpoint_density_transform(kpoints, scf_env%p_mix_new, .FALSE., &
1609 : matrix_s(1, 1)%matrix, sab_nl, scf_env%scf_work1, &
1610 0 : overlap_rs=matrix_s)
1611 : END IF
1612 :
1613 0 : CALL kp_transform_plan_release(transform_plan)
1614 0 : CALL kp_diag_omp_workspaces_release(workspace)
1615 0 : CALL timestop(handle)
1616 :
1617 0 : END SUBROUTINE do_general_diag_kp_omp
1618 :
1619 : ! **************************************************************************************************
1620 : !> \brief Create one independent dense workspace per OpenMP worker.
1621 : !> \param workspace ...
1622 : !> \param nworkers ...
1623 : !> \param matrix_struct ...
1624 : !> \param mo_struct ...
1625 : ! **************************************************************************************************
1626 0 : SUBROUTINE kp_diag_omp_workspaces_create(workspace, nworkers, matrix_struct, mo_struct)
1627 :
1628 : TYPE(kp_diag_omp_workspace_type), ALLOCATABLE, &
1629 : DIMENSION(:), INTENT(OUT) :: workspace
1630 : INTEGER, INTENT(IN) :: nworkers
1631 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct, mo_struct
1632 :
1633 : INTEGER :: i
1634 :
1635 0 : ALLOCATE (workspace(nworkers))
1636 0 : DO i = 1, nworkers
1637 0 : CALL cp_cfm_create(workspace(i)%cksmat, matrix_struct)
1638 0 : CALL cp_cfm_create(workspace(i)%csmat, matrix_struct)
1639 0 : CALL cp_cfm_create(workspace(i)%csmat_base, matrix_struct)
1640 0 : CALL cp_cfm_create(workspace(i)%cwork, matrix_struct)
1641 0 : CALL cp_cfm_create(workspace(i)%cmos, mo_struct)
1642 : END DO
1643 :
1644 0 : END SUBROUTINE kp_diag_omp_workspaces_create
1645 :
1646 : ! **************************************************************************************************
1647 : !> \brief Release OpenMP k-point workspaces.
1648 : !> \param workspace ...
1649 : ! **************************************************************************************************
1650 0 : SUBROUTINE kp_diag_omp_workspaces_release(workspace)
1651 :
1652 : TYPE(kp_diag_omp_workspace_type), ALLOCATABLE, &
1653 : DIMENSION(:), INTENT(INOUT) :: workspace
1654 :
1655 : INTEGER :: i
1656 :
1657 0 : DO i = 1, SIZE(workspace)
1658 0 : CALL cp_cfm_release(workspace(i)%cksmat)
1659 0 : CALL cp_cfm_release(workspace(i)%csmat)
1660 0 : CALL cp_cfm_release(workspace(i)%csmat_base)
1661 0 : CALL cp_cfm_release(workspace(i)%cwork)
1662 0 : CALL cp_cfm_release(workspace(i)%cmos)
1663 : END DO
1664 0 : DEALLOCATE (workspace)
1665 :
1666 0 : END SUBROUTINE kp_diag_omp_workspaces_release
1667 :
1668 : ! **************************************************************************************************
1669 : !> \brief Fuse the real-space Fourier transform and DBCSR-to-dense expansion.
1670 : !> The target is private to one worker and the source blocks are read-only.
1671 : !> \param source ...
1672 : !> \param ispin ...
1673 : !> \param xkp ...
1674 : !> \param plan ...
1675 : !> \param TARGET ...
1676 : ! **************************************************************************************************
1677 0 : SUBROUTINE kp_build_cfm(source, ispin, xkp, plan, TARGET)
1678 :
1679 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: source
1680 : INTEGER, INTENT(IN) :: ispin
1681 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
1682 : TYPE(kp_transform_plan_type), INTENT(IN) :: plan
1683 : TYPE(cp_cfm_type), INTENT(IN) :: target
1684 :
1685 : COMPLEX(KIND=dp) :: phase
1686 0 : COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: full
1687 : INTEGER :: col_offset, i, row_offset
1688 : LOGICAL :: found
1689 : REAL(KIND=dp) :: arg
1690 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
1691 :
1692 0 : full => target%local_data
1693 0 : full = z_zero
1694 0 : DO i = 1, plan%nentry
1695 : CALL dbcsr_get_readonly_block_p(source(ispin, plan%image(i))%matrix, &
1696 0 : plan%row(i), plan%col(i), block, found)
1697 0 : IF (.NOT. found) CYCLE
1698 0 : row_offset = plan%row_offset(i)
1699 0 : col_offset = plan%col_offset(i)
1700 0 : arg = DOT_PRODUCT(REAL(plan%cell(:, i), KIND=dp), xkp)
1701 0 : phase = CMPLX(COS(twopi*arg), plan%symmetry_sign(i)*SIN(twopi*arg), KIND=dp)
1702 : full(row_offset:row_offset + SIZE(block, 1) - 1, &
1703 : col_offset:col_offset + SIZE(block, 2) - 1) = &
1704 : full(row_offset:row_offset + SIZE(block, 1) - 1, &
1705 0 : col_offset:col_offset + SIZE(block, 2) - 1) + phase*block
1706 0 : IF (plan%symmetric .AND. plan%row(i) /= plan%col(i)) THEN
1707 : full(col_offset:col_offset + SIZE(block, 2) - 1, &
1708 : row_offset:row_offset + SIZE(block, 1) - 1) = &
1709 : full(col_offset:col_offset + SIZE(block, 2) - 1, &
1710 0 : row_offset:row_offset + SIZE(block, 1) - 1) + CONJG(phase)*TRANSPOSE(block)
1711 : END IF
1712 : END DO
1713 :
1714 0 : END SUBROUTINE kp_build_cfm
1715 :
1716 : ! **************************************************************************************************
1717 : !> \brief Transform and diagonalize one k-point without DIIS.
1718 : !> \param ikp ...
1719 : !> \param workspace ...
1720 : !> \param matrix_ks ...
1721 : !> \param matrix_s ...
1722 : !> \param kpoints ...
1723 : !> \param scf_env ...
1724 : !> \param scf_control ...
1725 : !> \param transform_plan ...
1726 : !> \param xkp ...
1727 : !> \param nspin ...
1728 : ! **************************************************************************************************
1729 0 : SUBROUTINE kp_diag_job(ikp, workspace, matrix_ks, matrix_s, kpoints, scf_env, scf_control, &
1730 : transform_plan, xkp, nspin)
1731 :
1732 : INTEGER, INTENT(IN) :: ikp
1733 : TYPE(kp_diag_omp_workspace_type), INTENT(INOUT) :: workspace
1734 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
1735 : TYPE(kpoint_type), POINTER :: kpoints
1736 : TYPE(qs_scf_env_type), POINTER :: scf_env
1737 : TYPE(scf_control_type), POINTER :: scf_control
1738 : TYPE(kp_transform_plan_type), INTENT(IN) :: transform_plan
1739 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1740 : INTEGER, INTENT(IN) :: nspin
1741 :
1742 : INTEGER :: ispin
1743 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1744 : TYPE(cp_fm_type), POINTER :: imos, rmos
1745 : TYPE(kpoint_env_type), POINTER :: kp
1746 :
1747 0 : kp => kpoints%kp_env(ikp)%kpoint_env
1748 0 : CALL kp_build_cfm(matrix_s, 1, xkp(1:3, ikp), transform_plan, workspace%csmat_base)
1749 0 : DO ispin = 1, nspin
1750 0 : CALL kp_build_cfm(matrix_ks, ispin, xkp(1:3, ikp), transform_plan, workspace%cksmat)
1751 0 : CALL cp_cfm_to_cfm(workspace%csmat_base, workspace%csmat)
1752 0 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
1753 0 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1754 0 : IF (scf_env%cholesky_method == cholesky_off) THEN
1755 : CALL cp_cfm_geeig_canon_local(workspace%cksmat, workspace%csmat, workspace%cmos, &
1756 0 : eigenvalues, workspace%cwork, scf_control%eps_eigval)
1757 : ELSE
1758 0 : CALL cp_cfm_geeig_local(workspace%cksmat, workspace%csmat, workspace%cmos, eigenvalues)
1759 : END IF
1760 0 : kp%mos(2, ispin)%eigenvalues = eigenvalues
1761 0 : CALL cp_cfm_to_fm(workspace%cmos, rmos, imos)
1762 : END DO
1763 :
1764 0 : END SUBROUTINE kp_diag_job
1765 :
1766 : ! **************************************************************************************************
1767 : !> \brief Build the DIIS error for one k-point after the shared buffer has
1768 : !> already been allocated. All products use local BLAS, so worker
1769 : !> threads never enter ScaLAPACK or MPI.
1770 : !> \param ikp ...
1771 : !> \param workspace ...
1772 : !> \param matrix_ks ...
1773 : !> \param matrix_s ...
1774 : !> \param kpoints ...
1775 : !> \param ib ...
1776 : !> \param transform_plan ...
1777 : !> \param xkp ...
1778 : !> \param nspin ...
1779 : ! **************************************************************************************************
1780 0 : SUBROUTINE kp_diis_error_job_local(ikp, workspace, matrix_ks, matrix_s, kpoints, ib, &
1781 : transform_plan, xkp, nspin)
1782 :
1783 : INTEGER, INTENT(IN) :: ikp
1784 : TYPE(kp_diag_omp_workspace_type), INTENT(INOUT) :: workspace
1785 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
1786 : TYPE(kpoint_type), POINTER :: kpoints
1787 : INTEGER, INTENT(IN) :: ib
1788 : TYPE(kp_transform_plan_type), INTENT(IN) :: transform_plan
1789 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1790 : INTEGER, INTENT(IN) :: nspin
1791 :
1792 : INTEGER :: ispin
1793 : TYPE(kpoint_env_type), POINTER :: kp
1794 :
1795 0 : kp => kpoints%kp_env(ikp)%kpoint_env
1796 : CALL kp_build_cfm(matrix_s, 1, xkp(1:3, ikp), transform_plan, &
1797 0 : workspace%csmat_base)
1798 0 : DO ispin = 1, nspin
1799 : CALL kp_build_cfm(matrix_ks, ispin, xkp(1:3, ikp), transform_plan, &
1800 0 : workspace%cksmat)
1801 0 : CALL cp_cfm_to_cfm(workspace%csmat_base, workspace%csmat)
1802 0 : CALL kp_diis_calc_err_local(kpoints, kp, workspace, ib, ispin, ikp)
1803 : END DO
1804 :
1805 0 : END SUBROUTINE kp_diis_error_job_local
1806 :
1807 : ! **************************************************************************************************
1808 : !> \brief Local-BLAS equivalent of qs_diis_b_calc_err_kp for an allocated buffer.
1809 : !> \param kpoints ...
1810 : !> \param kp ...
1811 : !> \param workspace ...
1812 : !> \param ib ...
1813 : !> \param ispin ...
1814 : !> \param ikp ...
1815 : ! **************************************************************************************************
1816 0 : SUBROUTINE kp_diis_calc_err_local(kpoints, kp, workspace, ib, ispin, ikp)
1817 :
1818 : TYPE(kpoint_type), POINTER :: kpoints
1819 : TYPE(kpoint_env_type), POINTER :: kp
1820 : TYPE(kp_diag_omp_workspace_type), INTENT(INOUT) :: workspace
1821 : INTEGER, INTENT(IN) :: ib, ispin, ikp
1822 :
1823 0 : COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: c, error, h, kc, s, sc
1824 : INTEGER :: homo, nao, nmo
1825 : REAL(KIND=dp) :: maxocc
1826 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1827 0 : POINTER :: ci, cr
1828 : TYPE(cp_fm_type), POINTER :: imos, rmos
1829 :
1830 : CALL get_mo_set(kp%mos(1, ispin), nao=nao, nmo=nmo, homo=homo, &
1831 0 : mo_coeff=rmos, maxocc=maxocc)
1832 0 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1833 0 : CALL cp_fm_get_info(rmos, local_data=cr)
1834 0 : CALL cp_fm_get_info(imos, local_data=ci)
1835 :
1836 0 : c => workspace%cmos%local_data
1837 0 : c(1:nao, 1:nmo) = CMPLX(cr(1:nao, 1:nmo), ci(1:nao, 1:nmo), KIND=dp)
1838 :
1839 0 : h => kpoints%scf_diis_buffer%param(ib, ispin, ikp)%local_data
1840 0 : s => kpoints%scf_diis_buffer%smat(ikp)%local_data
1841 0 : h(:, :) = workspace%cksmat%local_data(:, :)
1842 0 : s(:, :) = workspace%csmat%local_data(:, :)
1843 :
1844 0 : kc => workspace%cksmat%local_data
1845 0 : sc => workspace%csmat%local_data
1846 : CALL zgemm('N', 'N', nao, homo, nao, CMPLX(maxocc, 0.0_dp, KIND=dp), &
1847 0 : h(1, 1), SIZE(h, 1), c(1, 1), SIZE(c, 1), z_zero, kc(1, 1), SIZE(kc, 1))
1848 : CALL zgemm('N', 'N', nao, homo, nao, CMPLX(2.0_dp, 0.0_dp, KIND=dp), &
1849 0 : s(1, 1), SIZE(s, 1), c(1, 1), SIZE(c, 1), z_zero, sc(1, 1), SIZE(sc, 1))
1850 :
1851 0 : error => kpoints%scf_diis_buffer%error(ib, ispin, ikp)%local_data
1852 : CALL zgemm('N', 'T', nao, nao, homo, z_one, sc(1, 1), SIZE(sc, 1), &
1853 0 : kc(1, 1), SIZE(kc, 1), z_zero, error(1, 1), SIZE(error, 1))
1854 : CALL zgemm('N', 'T', nao, nao, homo, z_one, kc(1, 1), SIZE(kc, 1), &
1855 0 : sc(1, 1), SIZE(sc, 1), -z_one, error(1, 1), SIZE(error, 1))
1856 :
1857 0 : END SUBROUTINE kp_diis_calc_err_local
1858 :
1859 : ! **************************************************************************************************
1860 : !> \brief Form the DIIS Hamiltonian and diagonalize one k-point.
1861 : !> \param ikp ...
1862 : !> \param workspace ...
1863 : !> \param kpoints ...
1864 : !> \param coeffs ...
1865 : !> \param nb ...
1866 : !> \param scf_env ...
1867 : !> \param scf_control ...
1868 : !> \param nspin ...
1869 : !> \param adiis_ks ...
1870 : !> \param transform_plan ...
1871 : !> \param xkp ...
1872 : !> \param use_adiis ...
1873 : !> \param diis_weight ...
1874 : ! **************************************************************************************************
1875 0 : SUBROUTINE kp_diis_diag_job(ikp, workspace, kpoints, coeffs, nb, scf_env, scf_control, nspin, &
1876 : adiis_ks, transform_plan, xkp, use_adiis, diis_weight)
1877 :
1878 : INTEGER, INTENT(IN) :: ikp
1879 : TYPE(kp_diag_omp_workspace_type), INTENT(INOUT) :: workspace
1880 : TYPE(kpoint_type), POINTER :: kpoints
1881 : COMPLEX(KIND=dp), DIMENSION(:), INTENT(IN) :: coeffs
1882 : INTEGER, INTENT(IN) :: nb
1883 : TYPE(qs_scf_env_type), POINTER :: scf_env
1884 : TYPE(scf_control_type), POINTER :: scf_control
1885 : INTEGER, INTENT(IN) :: nspin
1886 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: adiis_ks
1887 : TYPE(kp_transform_plan_type), INTENT(IN) :: transform_plan
1888 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1889 : LOGICAL, INTENT(IN) :: use_adiis
1890 : REAL(KIND=dp), INTENT(IN) :: diis_weight
1891 :
1892 : INTEGER :: ispin, jb
1893 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1894 : TYPE(cp_fm_type), POINTER :: imos, rmos
1895 : TYPE(kpoint_env_type), POINTER :: kp
1896 :
1897 0 : kp => kpoints%kp_env(ikp)%kpoint_env
1898 0 : DO ispin = 1, nspin
1899 0 : CALL cp_cfm_to_cfm(kpoints%scf_diis_buffer%smat(ikp), workspace%csmat)
1900 0 : CALL cp_cfm_set_all(workspace%cksmat, z_zero)
1901 0 : DO jb = 1, nb
1902 : CALL cp_cfm_scale_and_add(z_one, workspace%cksmat, coeffs(jb), &
1903 0 : kpoints%scf_diis_buffer%param(jb, ispin, ikp))
1904 : END DO
1905 0 : IF (use_adiis .AND. diis_weight < 1.0_dp) THEN
1906 0 : CALL kp_build_cfm(adiis_ks, ispin, xkp(1:3, ikp), transform_plan, workspace%csmat_base)
1907 : workspace%cksmat%local_data(:, :) = &
1908 : diis_weight*workspace%cksmat%local_data(:, :) + &
1909 0 : (1.0_dp - diis_weight)*workspace%csmat_base%local_data(:, :)
1910 : END IF
1911 0 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
1912 0 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1913 0 : IF (scf_env%cholesky_method == cholesky_off) THEN
1914 : CALL cp_cfm_geeig_canon_local(workspace%cksmat, workspace%csmat, workspace%cmos, &
1915 0 : eigenvalues, workspace%cwork, scf_control%eps_eigval)
1916 : ELSE
1917 0 : CALL cp_cfm_geeig_local(workspace%cksmat, workspace%csmat, workspace%cmos, eigenvalues)
1918 : END IF
1919 0 : kp%mos(2, ispin)%eigenvalues = eigenvalues
1920 0 : CALL cp_cfm_to_fm(workspace%cmos, rmos, imos)
1921 : END DO
1922 :
1923 0 : END SUBROUTINE kp_diis_diag_job
1924 :
1925 : ! **************************************************************************************************
1926 : !> \brief Kpoint diagonalization routine
1927 : !> Transforms matrices to kpoint, distributes kpoint groups, performs
1928 : !> general diagonalization (no storgae of overlap decomposition), stores
1929 : !> MOs, calculates occupation numbers, calculates density matrices
1930 : !> in kpoint representation, transforms density matrices to real space
1931 : !> \param matrix_ks Kohn-sham matrices (RS indices, global)
1932 : !> \param matrix_s Overlap matrices (RS indices, global)
1933 : !> \param kpoints Kpoint environment
1934 : !> \param fmwork FM work matrices [at least dimension 4] in full para_env
1935 : !> \par History
1936 : !> 08.2014 created [JGH]
1937 : ! **************************************************************************************************
1938 22 : SUBROUTINE diag_kp_basic(matrix_ks, matrix_s, kpoints, fmwork)
1939 :
1940 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
1941 : TYPE(kpoint_type), POINTER :: kpoints
1942 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
1943 :
1944 : CHARACTER(len=*), PARAMETER :: routineN = 'diag_kp_basic'
1945 :
1946 : INTEGER :: handle, igroup, ik, ikp, indx, ispin, &
1947 : kplocal, nkp, nkp_groups, nspin
1948 : INTEGER, DIMENSION(2) :: kp_range
1949 22 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1950 22 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1951 : LOGICAL :: my_kpgrp, use_real_wfn
1952 22 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1953 22 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1954 22 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
1955 : TYPE(cp_cfm_type) :: cksmat, cmos, csmat, cwork
1956 22 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools
1957 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct, mo_struct
1958 : TYPE(cp_fm_type) :: fmdummy, fmlocal, rksmat, rsmat
1959 : TYPE(cp_fm_type), POINTER :: imos, mo_coeff, rmos
1960 : TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix, tempmat, tmpmat
1961 : TYPE(kpoint_env_type), POINTER :: kp
1962 : TYPE(mp_para_env_type), POINTER :: para_env
1963 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1964 22 : POINTER :: sab_nl
1965 : TYPE(qs_matrix_pools_type), POINTER :: mpools
1966 :
1967 22 : CALL timeset(routineN, handle)
1968 :
1969 22 : NULLIFY (sab_nl)
1970 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
1971 : nkp_groups=nkp_groups, kp_dist=kp_dist, sab_nl=sab_nl, &
1972 22 : cell_to_index=cell_to_index)
1973 22 : CPASSERT(ASSOCIATED(sab_nl))
1974 22 : kplocal = kp_range(2) - kp_range(1) + 1
1975 :
1976 : ! use as template
1977 22 : tempmat => matrix_ks(1, 1)%matrix
1978 :
1979 : ! allocate some work matrices
1980 22 : ALLOCATE (rmatrix, cmatrix, tmpmat)
1981 22 : CALL dbcsr_create(rmatrix, template=tempmat, matrix_type=dbcsr_type_symmetric)
1982 22 : CALL dbcsr_create(cmatrix, template=tempmat, matrix_type=dbcsr_type_antisymmetric)
1983 22 : CALL dbcsr_create(tmpmat, template=tempmat, matrix_type=dbcsr_type_no_symmetry)
1984 22 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
1985 22 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
1986 :
1987 : ! fm pools to be used within a kpoint group
1988 22 : CALL get_kpoint_info(kpoints, mpools=mpools)
1989 22 : CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
1990 :
1991 22 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
1992 22 : CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
1993 :
1994 22 : IF (use_real_wfn) THEN
1995 0 : CALL cp_fm_create(rksmat, matrix_struct)
1996 0 : CALL cp_fm_create(rsmat, matrix_struct)
1997 : ELSE
1998 22 : CALL cp_cfm_create(cksmat, matrix_struct)
1999 22 : CALL cp_cfm_create(csmat, matrix_struct)
2000 22 : CALL cp_cfm_create(cwork, matrix_struct)
2001 22 : kp => kpoints%kp_env(1)%kpoint_env
2002 22 : CALL get_mo_set(kp%mos(1, 1), mo_coeff=mo_coeff)
2003 22 : CALL cp_fm_get_info(mo_coeff, matrix_struct=mo_struct)
2004 22 : CALL cp_cfm_create(cmos, mo_struct)
2005 : END IF
2006 :
2007 22 : para_env => kpoints%blacs_env_all%para_env
2008 22 : nspin = SIZE(matrix_ks, 1)
2009 634 : ALLOCATE (info(kplocal*nspin*nkp_groups, 4))
2010 :
2011 : ! Setup and start all the communication
2012 22 : indx = 0
2013 94 : DO ikp = 1, kplocal
2014 166 : DO ispin = 1, nspin
2015 220 : DO igroup = 1, nkp_groups
2016 : ! number of current kpoint
2017 76 : ik = kp_dist(1, igroup) + ikp - 1
2018 76 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
2019 76 : indx = indx + 1
2020 76 : IF (use_real_wfn) THEN
2021 : ! FT of matrices KS and S, then transfer to FM type
2022 0 : CALL dbcsr_set(rmatrix, 0.0_dp)
2023 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_ks, ispin=ispin, &
2024 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
2025 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
2026 0 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
2027 : ! s matrix is not spin dependent
2028 0 : CALL dbcsr_set(rmatrix, 0.0_dp)
2029 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_s, ispin=1, &
2030 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
2031 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
2032 0 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(3))
2033 : ELSE
2034 : ! FT of matrices KS and S, then transfer to FM type
2035 76 : CALL dbcsr_set(rmatrix, 0.0_dp)
2036 76 : CALL dbcsr_set(cmatrix, 0.0_dp)
2037 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_ks, ispin=ispin, &
2038 76 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
2039 76 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
2040 76 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
2041 76 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
2042 76 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(2))
2043 : ! s matrix is not spin dependent, double the work
2044 76 : CALL dbcsr_set(rmatrix, 0.0_dp)
2045 76 : CALL dbcsr_set(cmatrix, 0.0_dp)
2046 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_s, ispin=1, &
2047 76 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
2048 76 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
2049 76 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(3))
2050 76 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
2051 76 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(4))
2052 : END IF
2053 : ! transfer to kpoint group
2054 : ! redistribution of matrices, new blacs environment
2055 : ! fmwork -> fmlocal -> rksmat/cksmat
2056 : ! fmwork -> fmlocal -> rsmat/csmat
2057 148 : IF (use_real_wfn) THEN
2058 0 : IF (my_kpgrp) THEN
2059 0 : CALL cp_fm_start_copy_general(fmwork(1), rksmat, para_env, info(indx, 1))
2060 0 : CALL cp_fm_start_copy_general(fmwork(3), rsmat, para_env, info(indx, 2))
2061 : ELSE
2062 0 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
2063 0 : CALL cp_fm_start_copy_general(fmwork(3), fmdummy, para_env, info(indx, 2))
2064 : END IF
2065 : ELSE
2066 76 : IF (my_kpgrp) THEN
2067 72 : CALL cp_fm_start_copy_general(fmwork(1), fmlocal, para_env, info(indx, 1))
2068 72 : CALL cp_fm_start_copy_general(fmwork(2), fmlocal, para_env, info(indx, 2))
2069 72 : CALL cp_fm_start_copy_general(fmwork(3), fmlocal, para_env, info(indx, 3))
2070 72 : CALL cp_fm_start_copy_general(fmwork(4), fmlocal, para_env, info(indx, 4))
2071 : ELSE
2072 4 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
2073 4 : CALL cp_fm_start_copy_general(fmwork(2), fmdummy, para_env, info(indx, 2))
2074 4 : CALL cp_fm_start_copy_general(fmwork(3), fmdummy, para_env, info(indx, 3))
2075 4 : CALL cp_fm_start_copy_general(fmwork(4), fmdummy, para_env, info(indx, 4))
2076 : END IF
2077 : END IF
2078 : END DO
2079 : END DO
2080 : END DO
2081 :
2082 : ! Finish communication then diagonalise in each group
2083 : indx = 0
2084 94 : DO ikp = 1, kplocal
2085 166 : DO ispin = 1, nspin
2086 148 : DO igroup = 1, nkp_groups
2087 : ! number of current kpoint
2088 76 : ik = kp_dist(1, igroup) + ikp - 1
2089 76 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
2090 4 : indx = indx + 1
2091 72 : IF (my_kpgrp) THEN
2092 72 : IF (use_real_wfn) THEN
2093 0 : CALL cp_fm_finish_copy_general(rksmat, info(indx, 1))
2094 0 : CALL cp_fm_finish_copy_general(rsmat, info(indx, 2))
2095 : ELSE
2096 72 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
2097 72 : CALL cp_cfm_scale_and_add_fm(z_zero, cksmat, z_one, fmlocal)
2098 72 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
2099 72 : CALL cp_cfm_scale_and_add_fm(z_one, cksmat, gaussi, fmlocal)
2100 72 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 3))
2101 72 : CALL cp_cfm_scale_and_add_fm(z_zero, csmat, z_one, fmlocal)
2102 72 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 4))
2103 72 : CALL cp_cfm_scale_and_add_fm(z_one, csmat, gaussi, fmlocal)
2104 : END IF
2105 : END IF
2106 : END DO
2107 :
2108 : ! Each kpoint group has now information on a kpoint to be diagonalized
2109 : ! General eigensolver Hermite or Symmetric
2110 72 : kp => kpoints%kp_env(ikp)%kpoint_env
2111 144 : IF (use_real_wfn) THEN
2112 0 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=mo_coeff, eigenvalues=eigenvalues)
2113 0 : CALL cp_fm_geeig(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal)
2114 : ELSE
2115 72 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
2116 72 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
2117 72 : CALL cp_cfm_geeig(cksmat, csmat, cmos, eigenvalues, cwork)
2118 : ! copy eigenvalues to imag set (keep them in sync)
2119 4752 : kp%mos(2, ispin)%eigenvalues = eigenvalues
2120 : ! split real and imaginary part of mos
2121 72 : CALL cp_cfm_to_fm(cmos, rmos, imos)
2122 : END IF
2123 : END DO
2124 : END DO
2125 :
2126 : ! Clean up communication
2127 : indx = 0
2128 94 : DO ikp = 1, kplocal
2129 166 : DO ispin = 1, nspin
2130 220 : DO igroup = 1, nkp_groups
2131 : ! number of current kpoint
2132 76 : ik = kp_dist(1, igroup) + ikp - 1
2133 76 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
2134 76 : indx = indx + 1
2135 148 : IF (use_real_wfn) THEN
2136 0 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
2137 0 : CALL cp_fm_cleanup_copy_general(info(indx, 2))
2138 : ELSE
2139 76 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
2140 76 : CALL cp_fm_cleanup_copy_general(info(indx, 2))
2141 76 : CALL cp_fm_cleanup_copy_general(info(indx, 3))
2142 76 : CALL cp_fm_cleanup_copy_general(info(indx, 4))
2143 : END IF
2144 : END DO
2145 : END DO
2146 : END DO
2147 :
2148 : ! All done
2149 348 : DEALLOCATE (info)
2150 :
2151 22 : CALL dbcsr_deallocate_matrix(rmatrix)
2152 22 : CALL dbcsr_deallocate_matrix(cmatrix)
2153 22 : CALL dbcsr_deallocate_matrix(tmpmat)
2154 :
2155 22 : IF (use_real_wfn) THEN
2156 0 : CALL cp_fm_release(rksmat)
2157 0 : CALL cp_fm_release(rsmat)
2158 : ELSE
2159 22 : CALL cp_cfm_release(cksmat)
2160 22 : CALL cp_cfm_release(csmat)
2161 22 : CALL cp_cfm_release(cwork)
2162 22 : CALL cp_cfm_release(cmos)
2163 : END IF
2164 22 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
2165 :
2166 22 : CALL timestop(handle)
2167 :
2168 44 : END SUBROUTINE diag_kp_basic
2169 :
2170 : ! **************************************************************************************************
2171 : !> \brief inner loop within MOS subspace, to refine occupation and density,
2172 : !> before next diagonalization of the Hamiltonian
2173 : !> \param qs_env ...
2174 : !> \param scf_env ...
2175 : !> \param subspace_env ...
2176 : !> \param mos ...
2177 : !> \param rho ...
2178 : !> \param ks_env ...
2179 : !> \param scf_section ...
2180 : !> \param scf_control ...
2181 : !> \par History
2182 : !> 09.2009 created [MI]
2183 : !> \note it is assumed that when diagonalization is used, also some mixing procedure is active
2184 : ! **************************************************************************************************
2185 10 : SUBROUTINE do_scf_diag_subspace(qs_env, scf_env, subspace_env, mos, rho, &
2186 : ks_env, scf_section, scf_control)
2187 :
2188 : TYPE(qs_environment_type), POINTER :: qs_env
2189 : TYPE(qs_scf_env_type), POINTER :: scf_env
2190 : TYPE(subspace_env_type), POINTER :: subspace_env
2191 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
2192 : TYPE(qs_rho_type), POINTER :: rho
2193 : TYPE(qs_ks_env_type), POINTER :: ks_env
2194 : TYPE(section_vals_type), POINTER :: scf_section
2195 : TYPE(scf_control_type), POINTER :: scf_control
2196 :
2197 : CHARACTER(LEN=*), PARAMETER :: routineN = 'do_scf_diag_subspace'
2198 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
2199 :
2200 : INTEGER :: handle, i, iloop, ispin, nao, nmo, &
2201 : nspin, output_unit
2202 : LOGICAL :: converged
2203 : REAL(dp) :: ene_diff, ene_old, iter_delta, max_val, &
2204 : sum_band, sum_val, t1, t2
2205 10 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues, mo_occupations
2206 10 : TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: eval_first, occ_first
2207 : TYPE(cp_fm_type) :: work
2208 : TYPE(cp_fm_type), POINTER :: c0, chc, evec, mo_coeff
2209 : TYPE(cp_logger_type), POINTER :: logger
2210 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
2211 10 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
2212 : TYPE(dft_control_type), POINTER :: dft_control
2213 : TYPE(mp_para_env_type), POINTER :: para_env
2214 : TYPE(qs_energy_type), POINTER :: energy
2215 10 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
2216 :
2217 10 : CALL timeset(routineN, handle)
2218 10 : NULLIFY (c0, chc, energy, evec, matrix_ks, mo_coeff, mo_eigenvalues, &
2219 10 : mo_occupations, dft_control, rho_ao, rho_ao_kp)
2220 :
2221 10 : logger => cp_get_default_logger()
2222 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIAG_SUB_SCF", &
2223 10 : extension=".scfLog")
2224 :
2225 : !Extra loop keeping mos unchanged and refining the subspace occupation
2226 10 : nspin = SIZE(mos)
2227 10 : CALL qs_rho_get(rho, rho_ao=rho_ao, rho_ao_kp=rho_ao_kp)
2228 :
2229 40 : ALLOCATE (eval_first(nspin))
2230 40 : ALLOCATE (occ_first(nspin))
2231 20 : DO ispin = 1, nspin
2232 : CALL get_mo_set(mo_set=mos(ispin), &
2233 : nmo=nmo, &
2234 : eigenvalues=mo_eigenvalues, &
2235 10 : occupation_numbers=mo_occupations)
2236 30 : ALLOCATE (eval_first(ispin)%array(nmo))
2237 20 : ALLOCATE (occ_first(ispin)%array(nmo))
2238 50 : eval_first(ispin)%array(1:nmo) = mo_eigenvalues(1:nmo)
2239 70 : occ_first(ispin)%array(1:nmo) = mo_occupations(1:nmo)
2240 : END DO
2241 :
2242 20 : DO ispin = 1, nspin
2243 : ! does not yet handle k-points
2244 10 : CALL dbcsr_copy(subspace_env%p_matrix_store(ispin)%matrix, rho_ao(ispin)%matrix)
2245 20 : CALL dbcsr_copy(rho_ao(ispin)%matrix, scf_env%p_mix_new(ispin, 1)%matrix)
2246 : END DO
2247 :
2248 10 : subspace_env%p_matrix_mix => scf_env%p_mix_new
2249 :
2250 10 : NULLIFY (matrix_ks, energy, para_env, matrix_s)
2251 : CALL get_qs_env(qs_env, &
2252 : matrix_ks=matrix_ks, &
2253 : energy=energy, &
2254 : matrix_s=matrix_s, &
2255 : para_env=para_env, &
2256 10 : dft_control=dft_control)
2257 :
2258 : ! mixing storage allocation
2259 10 : IF (subspace_env%mixing_method >= gspace_mixing_nr) THEN
2260 : CALL mixing_allocate(qs_env, subspace_env%mixing_method, scf_env%p_mix_new, &
2261 0 : scf_env%p_delta, nspin, subspace_env%mixing_store)
2262 0 : IF (dft_control%qs_control%gapw) THEN
2263 0 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
2264 : CALL mixing_init(subspace_env%mixing_method, rho, subspace_env%mixing_store, &
2265 0 : para_env, rho_atom=rho_atom)
2266 0 : ELSE IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
2267 0 : CALL charge_mixing_init(subspace_env%mixing_store)
2268 0 : ELSE IF (dft_control%qs_control%semi_empirical) THEN
2269 0 : CPABORT('SE Code not possible')
2270 : ELSE
2271 0 : CALL mixing_init(subspace_env%mixing_method, rho, subspace_env%mixing_store, para_env)
2272 : END IF
2273 : END IF
2274 :
2275 10 : ene_old = 0.0_dp
2276 : ene_diff = 0.0_dp
2277 10 : IF (output_unit > 0) THEN
2278 0 : WRITE (output_unit, "(/T19,A)") '<<<<<<<<< SUBSPACE ROTATION <<<<<<<<<<'
2279 : WRITE (output_unit, "(T4,A,T13,A,T21,A,T38,A,T51,A,T65,A/,T4,A)") &
2280 0 : "In-step", "Time", "Convergence", "Band ene.", "Total ene.", "Energy diff.", REPEAT("-", 74)
2281 : END IF
2282 :
2283 : ! recalculate density matrix here
2284 :
2285 : ! update of density
2286 10 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
2287 :
2288 22 : DO iloop = 1, subspace_env%max_iter
2289 20 : t1 = m_walltime()
2290 20 : converged = .FALSE.
2291 20 : ene_old = energy%total
2292 :
2293 20 : CALL qs_ks_did_change(ks_env, rho_changed=.TRUE.)
2294 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., &
2295 20 : just_energy=.FALSE., print_active=.FALSE.)
2296 :
2297 20 : max_val = 0.0_dp
2298 20 : sum_val = 0.0_dp
2299 20 : sum_band = 0.0_dp
2300 40 : DO ispin = 1, SIZE(matrix_ks)
2301 : CALL get_mo_set(mo_set=mos(ispin), &
2302 : nao=nao, &
2303 : nmo=nmo, &
2304 : eigenvalues=mo_eigenvalues, &
2305 : occupation_numbers=mo_occupations, &
2306 20 : mo_coeff=mo_coeff)
2307 :
2308 : !compute C'HC
2309 20 : chc => subspace_env%chc_mat(ispin)
2310 20 : evec => subspace_env%c_vec(ispin)
2311 20 : c0 => subspace_env%c0(ispin)
2312 20 : CALL cp_fm_to_fm(mo_coeff, c0)
2313 20 : CALL cp_fm_create(work, c0%matrix_struct)
2314 20 : CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, c0, work, nmo)
2315 20 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, rone, c0, work, rzero, chc)
2316 20 : CALL cp_fm_release(work)
2317 : !diagonalize C'HC
2318 20 : CALL choose_eigv_solver(chc, evec, mo_eigenvalues)
2319 :
2320 : !rotate the mos by the eigenvectors of C'HC
2321 20 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, rone, c0, evec, rzero, mo_coeff)
2322 :
2323 20 : IF (.NOT. scf_control%gce%do_gce) THEN
2324 : CALL set_mo_occupation(mo_set=mos(ispin), &
2325 20 : smear=scf_control%smear)
2326 : ELSE
2327 : CALL set_mo_occupation(mo_set=mos(ispin), &
2328 : smear=scf_control%smear, &
2329 0 : gce=scf_control%gce)
2330 : END IF
2331 :
2332 : ! does not yet handle k-points
2333 : CALL calculate_density_matrix(mos(ispin), &
2334 20 : subspace_env%p_matrix_mix(ispin, 1)%matrix)
2335 :
2336 160 : DO i = 1, nmo
2337 100 : sum_band = sum_band + mo_eigenvalues(i)*mo_occupations(i)
2338 : END DO
2339 :
2340 : !check for self consistency
2341 : END DO
2342 :
2343 20 : IF (subspace_env%mixing_method == direct_mixing_nr) THEN
2344 : CALL scf_env_density_mixing(subspace_env%p_matrix_mix, &
2345 20 : scf_env%mixing_store, rho_ao_kp, para_env, iter_delta, iloop)
2346 : ELSE
2347 : CALL self_consistency_check(rho_ao_kp, scf_env%p_delta, para_env, &
2348 0 : subspace_env%p_matrix_mix, delta=iter_delta)
2349 : END IF
2350 :
2351 40 : DO ispin = 1, nspin
2352 : ! does not yet handle k-points
2353 40 : CALL dbcsr_copy(rho_ao(ispin)%matrix, subspace_env%p_matrix_mix(ispin, 1)%matrix)
2354 : END DO
2355 : ! update of density
2356 20 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
2357 : ! Mixing in reciprocal space
2358 20 : IF (subspace_env%mixing_method >= gspace_mixing_nr) THEN
2359 : CALL gspace_mixing(qs_env, scf_env%mixing_method, subspace_env%mixing_store, &
2360 0 : rho, para_env, scf_env%iter_count)
2361 : END IF
2362 :
2363 20 : ene_diff = energy%total - ene_old
2364 : converged = (ABS(ene_diff) < subspace_env%eps_ene .AND. &
2365 20 : iter_delta < subspace_env%eps_adapt*scf_env%iter_delta)
2366 20 : t2 = m_walltime()
2367 20 : IF (output_unit > 0) THEN
2368 : WRITE (output_unit, "(T4,I5,T11,F8.3,T18,E14.4,T34,F12.5,T46,F16.8,T62,E14.4)") &
2369 0 : iloop, t2 - t1, iter_delta, sum_band, energy%total, ene_diff
2370 0 : CALL m_flush(output_unit)
2371 : END IF
2372 22 : IF (converged) THEN
2373 8 : IF (output_unit > 0) WRITE (output_unit, "(T10,A,I6,A,/)") &
2374 0 : " Reached convergence in ", iloop, " iterations "
2375 : EXIT
2376 : END IF
2377 :
2378 : END DO ! iloop
2379 :
2380 10 : NULLIFY (subspace_env%p_matrix_mix)
2381 20 : DO ispin = 1, nspin
2382 : ! does not yet handle k-points
2383 10 : CALL dbcsr_copy(scf_env%p_mix_new(ispin, 1)%matrix, rho_ao(ispin)%matrix)
2384 10 : CALL dbcsr_copy(rho_ao(ispin)%matrix, subspace_env%p_matrix_store(ispin)%matrix)
2385 :
2386 20 : DEALLOCATE (eval_first(ispin)%array, occ_first(ispin)%array)
2387 : END DO
2388 10 : DEALLOCATE (eval_first, occ_first)
2389 :
2390 10 : CALL timestop(handle)
2391 :
2392 10 : END SUBROUTINE do_scf_diag_subspace
2393 :
2394 : ! **************************************************************************************************
2395 : !> \brief ...
2396 : !> \param subspace_env ...
2397 : !> \param qs_env ...
2398 : !> \param mos ...
2399 : ! **************************************************************************************************
2400 2 : SUBROUTINE diag_subspace_allocate(subspace_env, qs_env, mos)
2401 :
2402 : TYPE(subspace_env_type), POINTER :: subspace_env
2403 : TYPE(qs_environment_type), POINTER :: qs_env
2404 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
2405 :
2406 : CHARACTER(LEN=*), PARAMETER :: routineN = 'diag_subspace_allocate'
2407 :
2408 : INTEGER :: handle, i, ispin, nmo, nspin
2409 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
2410 : TYPE(cp_fm_type), POINTER :: mo_coeff
2411 2 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
2412 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2413 2 : POINTER :: sab_orb
2414 :
2415 2 : CALL timeset(routineN, handle)
2416 :
2417 2 : NULLIFY (sab_orb, matrix_s)
2418 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb, &
2419 2 : matrix_s=matrix_s)
2420 :
2421 2 : nspin = SIZE(mos)
2422 : ! *** allocate p_atrix_store ***
2423 2 : IF (.NOT. ASSOCIATED(subspace_env%p_matrix_store)) THEN
2424 2 : CALL dbcsr_allocate_matrix_set(subspace_env%p_matrix_store, nspin)
2425 :
2426 4 : DO i = 1, nspin
2427 2 : ALLOCATE (subspace_env%p_matrix_store(i)%matrix)
2428 : CALL dbcsr_create(matrix=subspace_env%p_matrix_store(i)%matrix, template=matrix_s(1)%matrix, &
2429 2 : name="DENSITY_STORE", matrix_type=dbcsr_type_symmetric)
2430 : CALL cp_dbcsr_alloc_block_from_nbl(subspace_env%p_matrix_store(i)%matrix, &
2431 2 : sab_orb)
2432 4 : CALL dbcsr_set(subspace_env%p_matrix_store(i)%matrix, 0.0_dp)
2433 : END DO
2434 :
2435 : END IF
2436 :
2437 8 : ALLOCATE (subspace_env%chc_mat(nspin))
2438 6 : ALLOCATE (subspace_env%c_vec(nspin))
2439 6 : ALLOCATE (subspace_env%c0(nspin))
2440 :
2441 4 : DO ispin = 1, nspin
2442 2 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
2443 2 : CALL cp_fm_create(subspace_env%c0(ispin), mo_coeff%matrix_struct)
2444 2 : NULLIFY (fm_struct_tmp)
2445 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo, ncol_global=nmo, &
2446 : para_env=mo_coeff%matrix_struct%para_env, &
2447 2 : context=mo_coeff%matrix_struct%context)
2448 2 : CALL cp_fm_create(subspace_env%chc_mat(ispin), fm_struct_tmp, "chc")
2449 2 : CALL cp_fm_create(subspace_env%c_vec(ispin), fm_struct_tmp, "vec")
2450 6 : CALL cp_fm_struct_release(fm_struct_tmp)
2451 : END DO
2452 :
2453 2 : CALL timestop(handle)
2454 :
2455 2 : END SUBROUTINE diag_subspace_allocate
2456 :
2457 : ! **************************************************************************************************
2458 : !> \brief the inner loop of scf, specific to diagonalization without S matrix
2459 : !> basically, in goes the ks matrix out goes a new p matrix
2460 : !> \param scf_env ...
2461 : !> \param mos ...
2462 : !> \param matrix_ks ...
2463 : !> \param scf_control ...
2464 : !> \param scf_section ...
2465 : !> \param diis_step ...
2466 : !> \par History
2467 : !> 03.2006 created [Joost VandeVondele]
2468 : ! **************************************************************************************************
2469 17876 : SUBROUTINE do_special_diag(scf_env, mos, matrix_ks, scf_control, &
2470 : scf_section, diis_step)
2471 :
2472 : TYPE(qs_scf_env_type), POINTER :: scf_env
2473 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
2474 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
2475 : TYPE(scf_control_type), POINTER :: scf_control
2476 : TYPE(section_vals_type), POINTER :: scf_section
2477 : LOGICAL, INTENT(INOUT) :: diis_step
2478 :
2479 : INTEGER :: ispin, nspin
2480 : LOGICAL :: do_level_shift, use_jacobi
2481 : REAL(KIND=dp) :: diis_error
2482 :
2483 17876 : nspin = SIZE(matrix_ks)
2484 :
2485 36518 : DO ispin = 1, nspin
2486 36518 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, scf_env%scf_work1(ispin))
2487 : END DO
2488 17876 : IF (scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis) THEN
2489 : CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
2490 : scf_env%scf_work2, scf_env%iter_delta, diis_error, diis_step, &
2491 : scf_control%eps_diis, scf_control%nmixing, &
2492 15282 : scf_section=scf_section)
2493 : ELSE
2494 2594 : diis_step = .FALSE.
2495 : END IF
2496 :
2497 17876 : IF ((scf_env%iter_count > 1) .AND. (scf_env%iter_delta < scf_control%diagonalization%eps_jacobi)) THEN
2498 18 : use_jacobi = .TRUE.
2499 : ELSE
2500 17858 : use_jacobi = .FALSE.
2501 : END IF
2502 :
2503 : do_level_shift = ((scf_control%level_shift /= 0.0_dp) .AND. &
2504 17876 : ((scf_control%density_guess == core_guess) .OR. (scf_env%iter_count > 1)))
2505 17876 : IF (diis_step) THEN
2506 11870 : scf_env%iter_param = diis_error
2507 11870 : IF (use_jacobi) THEN
2508 18 : scf_env%iter_method = "DIIS/Jacobi"
2509 : ELSE
2510 11852 : scf_env%iter_method = "DIIS/Diag."
2511 : END IF
2512 : ELSE
2513 6006 : IF (scf_env%mixing_method == 1) THEN
2514 6006 : scf_env%iter_param = scf_env%p_mix_alpha
2515 6006 : IF (use_jacobi) THEN
2516 0 : scf_env%iter_method = "P_Mix/Jacobi"
2517 : ELSE
2518 6006 : scf_env%iter_method = "P_Mix/Diag."
2519 : END IF
2520 0 : ELSE IF (scf_env%mixing_method > 1) THEN
2521 0 : scf_env%iter_param = scf_env%mixing_store%alpha
2522 0 : IF (use_jacobi) THEN
2523 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Jacobi"
2524 : ELSE
2525 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Diag."
2526 : END IF
2527 : END IF
2528 : END IF
2529 17876 : scf_env%iter_delta = 0.0_dp
2530 :
2531 36518 : DO ispin = 1, nspin
2532 : CALL eigensolver_simple(matrix_ks=scf_env%scf_work1(ispin), &
2533 : mo_set=mos(ispin), &
2534 : work=scf_env%scf_work2, &
2535 : do_level_shift=do_level_shift, &
2536 : level_shift=scf_control%level_shift, &
2537 : use_jacobi=use_jacobi, &
2538 36518 : jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
2539 : END DO
2540 :
2541 17876 : IF (.NOT. scf_control%gce%do_gce) THEN
2542 : CALL set_mo_occupation(mo_array=mos, &
2543 17876 : smear=scf_control%smear)
2544 : ELSE
2545 : CALL set_mo_occupation(mo_array=mos, &
2546 : smear=scf_control%smear, &
2547 0 : gce=scf_control%gce)
2548 : END IF
2549 :
2550 36518 : DO ispin = 1, nspin
2551 : ! does not yet handle k-points
2552 : CALL calculate_density_matrix(mos(ispin), &
2553 36518 : scf_env%p_mix_new(ispin, 1)%matrix)
2554 : END DO
2555 :
2556 17876 : END SUBROUTINE do_special_diag
2557 :
2558 : ! **************************************************************************************************
2559 : !> \brief the inner loop of scf, specific to iterative diagonalization using OT
2560 : !> with S matrix; basically, in goes the ks matrix out goes a new p matrix
2561 : !> \param scf_env ...
2562 : !> \param mos ...
2563 : !> \param matrix_ks ...
2564 : !> \param matrix_s ...
2565 : !> \param scf_control ...
2566 : !> \param scf_section ...
2567 : !> \param diis_step ...
2568 : !> \par History
2569 : !> 10.2008 created [JGH]
2570 : ! **************************************************************************************************
2571 202 : SUBROUTINE do_ot_diag(scf_env, mos, matrix_ks, matrix_s, &
2572 : scf_control, scf_section, diis_step)
2573 :
2574 : TYPE(qs_scf_env_type), POINTER :: scf_env
2575 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
2576 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
2577 : TYPE(scf_control_type), POINTER :: scf_control
2578 : TYPE(section_vals_type), POINTER :: scf_section
2579 : LOGICAL, INTENT(INOUT) :: diis_step
2580 :
2581 : INTEGER :: homo, ispin, nmo, nspin
2582 : REAL(KIND=dp) :: diis_error, eps_iter
2583 202 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
2584 : TYPE(cp_fm_type), POINTER :: mo_coeff
2585 :
2586 202 : NULLIFY (eigenvalues)
2587 :
2588 202 : nspin = SIZE(matrix_ks)
2589 :
2590 448 : DO ispin = 1, nspin
2591 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, &
2592 448 : scf_env%scf_work1(ispin))
2593 : END DO
2594 :
2595 202 : IF ((scf_env%iter_count > 1) .AND. (.NOT. scf_env%skip_diis)) THEN
2596 : CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
2597 : scf_env%scf_work2, scf_env%iter_delta, diis_error, diis_step, &
2598 : scf_control%eps_diis, scf_control%nmixing, &
2599 : s_matrix=matrix_s, &
2600 180 : scf_section=scf_section)
2601 : ELSE
2602 22 : diis_step = .FALSE.
2603 : END IF
2604 :
2605 202 : eps_iter = scf_control%diagonalization%eps_iter
2606 202 : IF (diis_step) THEN
2607 130 : scf_env%iter_param = diis_error
2608 130 : scf_env%iter_method = "DIIS/OTdiag"
2609 274 : DO ispin = 1, nspin
2610 : CALL copy_fm_to_dbcsr(scf_env%scf_work1(ispin), &
2611 274 : matrix_ks(ispin)%matrix, keep_sparsity=.TRUE.)
2612 : END DO
2613 130 : eps_iter = MAX(eps_iter, scf_control%diagonalization%eps_adapt*diis_error)
2614 : ELSE
2615 72 : IF (scf_env%mixing_method == 1) THEN
2616 72 : scf_env%iter_param = scf_env%p_mix_alpha
2617 72 : scf_env%iter_method = "P_Mix/OTdiag."
2618 0 : ELSE IF (scf_env%mixing_method > 1) THEN
2619 0 : scf_env%iter_param = scf_env%mixing_store%alpha
2620 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/OTdiag."
2621 : END IF
2622 : END IF
2623 :
2624 202 : scf_env%iter_delta = 0.0_dp
2625 :
2626 448 : DO ispin = 1, nspin
2627 : CALL get_mo_set(mos(ispin), &
2628 : mo_coeff=mo_coeff, &
2629 : eigenvalues=eigenvalues, &
2630 : nmo=nmo, &
2631 246 : homo=homo)
2632 : CALL ot_eigensolver(matrix_h=matrix_ks(ispin)%matrix, &
2633 : matrix_s=matrix_s(1)%matrix, &
2634 : matrix_c_fm=mo_coeff, &
2635 : preconditioner=scf_env%ot_preconditioner( &
2636 : MIN(ispin, SIZE(scf_env%ot_preconditioner)))%preconditioner, &
2637 : eps_gradient=eps_iter, &
2638 : iter_max=scf_control%diagonalization%max_iter, &
2639 : silent=.TRUE., &
2640 246 : ot_settings=scf_control%diagonalization%ot_settings)
2641 : CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin)%matrix, &
2642 : evals_arg=eigenvalues, &
2643 246 : do_rotation=.TRUE.)
2644 : CALL copy_fm_to_dbcsr(mos(ispin)%mo_coeff, &
2645 694 : mos(ispin)%mo_coeff_b)
2646 : !fm->dbcsr
2647 : END DO
2648 :
2649 : CALL set_mo_occupation(mo_array=mos, &
2650 202 : smear=scf_control%smear)
2651 :
2652 448 : DO ispin = 1, nspin
2653 : ! does not yet handle k-points
2654 : CALL calculate_density_matrix(mos(ispin), &
2655 448 : scf_env%p_mix_new(ispin, 1)%matrix)
2656 : END DO
2657 :
2658 202 : END SUBROUTINE do_ot_diag
2659 :
2660 : ! **************************************************************************************************
2661 : !> \brief fixed-H OT diagonalization for real or complex K-point channels
2662 : !> \param matrix_ks real-space Kohn-Sham matrices
2663 : !> \param matrix_s real-space overlap matrices
2664 : !> \param matrix_t ...
2665 : !> \param kpoints K-point environment
2666 : !> \param scf_env SCF environment
2667 : !> \param scf_control SCF controls
2668 : !> \param diis_step set false; density mixing remains the outer accelerator
2669 : !> \param added_mos_auto_grow request a larger adaptive virtual-space buffer
2670 : ! **************************************************************************************************
2671 58 : SUBROUTINE do_ot_diag_kp(matrix_ks, matrix_s, matrix_t, kpoints, scf_env, scf_control, &
2672 : diis_step, added_mos_auto_grow)
2673 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s, matrix_t
2674 : TYPE(kpoint_type), POINTER :: kpoints
2675 : TYPE(qs_scf_env_type), POINTER :: scf_env
2676 : TYPE(scf_control_type), POINTER :: scf_control
2677 : LOGICAL, INTENT(INOUT) :: diis_step
2678 : LOGICAL, INTENT(OUT) :: added_mos_auto_grow
2679 :
2680 : INTEGER :: ikp, ispin, ks_spin, nspin
2681 : LOGICAL :: store_kinetic, use_real_wfn
2682 : REAL(KIND=dp) :: eps_iter
2683 58 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
2684 : TYPE(cp_fm_type) :: coeff_im_work
2685 58 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
2686 : TYPE(cp_fm_type), POINTER :: coeff_im, coeff_re
2687 : TYPE(dbcsr_type), POINTER :: matrix_h_im, matrix_h_re, matrix_s_im, &
2688 : matrix_s_re, matrix_t_im, matrix_t_re
2689 : TYPE(kpoint_env_type), POINTER :: kp
2690 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2691 58 : POINTER :: sab_nl
2692 :
2693 58 : NULLIFY (coeff_im, coeff_re, eigenvalues, fmwork, kp, matrix_h_im, matrix_h_re, &
2694 58 : matrix_s_im, matrix_s_re, matrix_t_im, matrix_t_re, sab_nl)
2695 58 : added_mos_auto_grow = .FALSE.
2696 58 : diis_step = .FALSE.
2697 58 : eps_iter = scf_control%diagonalization%eps_iter
2698 58 : nspin = SIZE(matrix_ks, 1)
2699 : store_kinetic = scf_control%diagonalization%ot_settings%preconditioner_type == &
2700 58 : ot_precond_full_kinetic
2701 58 : IF (store_kinetic) THEN
2702 0 : CPASSERT(ASSOCIATED(matrix_t))
2703 : END IF
2704 58 : CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn, sab_nl=sab_nl)
2705 58 : CPASSERT(ASSOCIATED(sab_nl))
2706 :
2707 : ! A full diagonalization supplies an orthonormal start after initialization or an
2708 : ! ADDED_MOS AUTO resize. Subsequent SCF steps only redistribute H(k) and S(k).
2709 58 : IF (scf_env%iter_count <= 1) THEN
2710 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, &
2711 6 : .FALSE., diis_step)
2712 : END IF
2713 58 : IF (store_kinetic) THEN
2714 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, &
2715 : .FALSE., diis_step, matrix_t=matrix_t, &
2716 0 : store_ot_matrices=.TRUE., transfer_only=.TRUE.)
2717 : ELSE
2718 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, &
2719 58 : .FALSE., diis_step, store_ot_matrices=.TRUE., transfer_only=.TRUE.)
2720 : END IF
2721 :
2722 58 : scf_env%iter_delta = 0.0_dp
2723 136 : DO ikp = 1, SIZE(kpoints%kp_env)
2724 78 : kp => kpoints%kp_env(ikp)%kpoint_env
2725 78 : CPASSERT(ASSOCIATED(kp%ot_hmat) .AND. ASSOCIATED(kp%ot_smat))
2726 78 : IF (SIZE(kp%ot_smat) >= 2) THEN
2727 : CALL qs_scf_get_ot_kpoint_operator(matrix_s, kpoints, kp, 1, &
2728 : kp%ot_smat(1), kp%ot_smat(2), &
2729 78 : matrix_s_re, matrix_s_im)
2730 : ELSE
2731 : CALL qs_scf_get_ot_kpoint_operator(matrix_s, kpoints, kp, 1, &
2732 : kp%ot_smat(1), kp%ot_smat(1), &
2733 0 : matrix_s_re, matrix_s_im)
2734 0 : CALL dbcsr_set(matrix_s_im, 0.0_dp)
2735 : END IF
2736 78 : IF (store_kinetic) THEN
2737 0 : CPASSERT(ASSOCIATED(kp%ot_tmat))
2738 0 : IF (SIZE(kp%ot_tmat) >= 2) THEN
2739 : CALL qs_scf_get_ot_kpoint_operator(matrix_t, kpoints, kp, 1, &
2740 : kp%ot_tmat(1), kp%ot_tmat(2), &
2741 0 : matrix_t_re, matrix_t_im)
2742 : ELSE
2743 : CALL qs_scf_get_ot_kpoint_operator(matrix_t, kpoints, kp, 1, &
2744 : kp%ot_tmat(1), kp%ot_tmat(1), &
2745 0 : matrix_t_re, matrix_t_im)
2746 0 : CALL dbcsr_set(matrix_t_im, 0.0_dp)
2747 : END IF
2748 : END IF
2749 156 : DO ispin = 1, nspin
2750 78 : ks_spin = MIN(ispin, SIZE(kp%ot_hmat, 2))
2751 78 : IF (SIZE(kp%ot_hmat, 1) >= 2) THEN
2752 : CALL qs_scf_get_ot_kpoint_operator(matrix_ks, kpoints, kp, ks_spin, &
2753 : kp%ot_hmat(1, ks_spin), kp%ot_hmat(2, ks_spin), &
2754 78 : matrix_h_re, matrix_h_im)
2755 : ELSE
2756 : CALL qs_scf_get_ot_kpoint_operator(matrix_ks, kpoints, kp, ks_spin, &
2757 : kp%ot_hmat(1, ks_spin), kp%ot_hmat(1, ks_spin), &
2758 0 : matrix_h_re, matrix_h_im)
2759 0 : CALL dbcsr_set(matrix_h_im, 0.0_dp)
2760 : END IF
2761 78 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=coeff_re, eigenvalues=eigenvalues)
2762 78 : IF (use_real_wfn) THEN
2763 : CALL cp_fm_create(coeff_im_work, coeff_re%matrix_struct, &
2764 0 : name='real k-point OT imaginary work')
2765 0 : CALL cp_fm_set_all(coeff_im_work, 0.0_dp)
2766 : CALL ot_eigensolver_complex(matrix_h=matrix_h_re, matrix_h_im=matrix_h_im, &
2767 : matrix_s=matrix_s_re, matrix_s_im=matrix_s_im, &
2768 : matrix_c_fm=coeff_re, matrix_c_fm_im=coeff_im_work, &
2769 : eps_gradient=eps_iter, &
2770 : iter_max=scf_control%diagonalization%max_iter, &
2771 : eigenvalues=eigenvalues, silent=.TRUE., &
2772 : ot_settings=scf_control%diagonalization%ot_settings, &
2773 : matrix_t=matrix_t_re, matrix_t_im=matrix_t_im, &
2774 0 : mo_set=kp%mos(1, ispin))
2775 0 : CALL cp_fm_release(coeff_im_work)
2776 : ELSE
2777 78 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=coeff_im)
2778 : CALL ot_eigensolver_complex(matrix_h=matrix_h_re, matrix_h_im=matrix_h_im, &
2779 : matrix_s=matrix_s_re, matrix_s_im=matrix_s_im, &
2780 : matrix_c_fm=coeff_re, matrix_c_fm_im=coeff_im, &
2781 : eps_gradient=eps_iter, &
2782 : iter_max=scf_control%diagonalization%max_iter, &
2783 : eigenvalues=eigenvalues, silent=.TRUE., &
2784 : ot_settings=scf_control%diagonalization%ot_settings, &
2785 : matrix_t=matrix_t_re, matrix_t_im=matrix_t_im, &
2786 78 : mo_set=kp%mos(1, ispin))
2787 1254 : kp%mos(2, ispin)%eigenvalues = eigenvalues
2788 : END IF
2789 78 : CALL dbcsr_release_p(matrix_h_re)
2790 156 : CALL dbcsr_release_p(matrix_h_im)
2791 : END DO
2792 78 : CALL dbcsr_release_p(matrix_s_re)
2793 78 : CALL dbcsr_release_p(matrix_s_im)
2794 78 : IF (ASSOCIATED(matrix_t_re)) CALL dbcsr_release_p(matrix_t_re)
2795 136 : IF (ASSOCIATED(matrix_t_im)) CALL dbcsr_release_p(matrix_t_im)
2796 : END DO
2797 :
2798 : CALL kpoint_set_mo_occupation(kpoints, scf_control%smear, &
2799 : added_mos_auto=scf_control%added_mos_auto, &
2800 58 : added_mos_auto_grow=added_mos_auto_grow)
2801 58 : IF (added_mos_auto_grow) RETURN
2802 :
2803 58 : CALL kpoint_density_matrices(kpoints)
2804 58 : fmwork => scf_env%scf_work1
2805 : CALL kpoint_density_transform(kpoints, scf_env%p_mix_new, .FALSE., &
2806 58 : matrix_s(1, 1)%matrix, sab_nl, fmwork, overlap_rs=matrix_s)
2807 58 : IF (scf_env%mixing_method == 0) THEN
2808 0 : scf_env%iter_method = 'NoMix/OTdiag.'
2809 58 : ELSE IF (scf_env%mixing_method == 1) THEN
2810 20 : scf_env%iter_param = scf_env%p_mix_alpha
2811 20 : scf_env%iter_method = 'P_Mix/OTdiag.'
2812 : ELSE
2813 38 : scf_env%iter_param = scf_env%mixing_store%alpha
2814 38 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//'/OTdiag.'
2815 : END IF
2816 58 : END SUBROUTINE do_ot_diag_kp
2817 :
2818 : ! **************************************************************************************************
2819 : !> \brief Solve a set restricted open Kohn-Sham (ROKS) equations based on the
2820 : !> alpha and beta Kohn-Sham matrices from unrestricted Kohn-Sham.
2821 : !> \param scf_env ...
2822 : !> \param mos ...
2823 : !> \param matrix_ks ...
2824 : !> \param matrix_s ...
2825 : !> \param scf_control ...
2826 : !> \param scf_section ...
2827 : !> \param diis_step ...
2828 : !> \param orthogonal_basis ...
2829 : !> \par History
2830 : !> 04.2006 created [MK]
2831 : !> Revised (01.05.06,MK)
2832 : !> \note
2833 : !> this is only a high-spin ROKS.
2834 : ! **************************************************************************************************
2835 1132 : SUBROUTINE do_roks_diag(scf_env, mos, matrix_ks, matrix_s, &
2836 : scf_control, scf_section, diis_step, &
2837 : orthogonal_basis)
2838 :
2839 : ! Literature: - C. C. J. Roothaan, Rev. Mod. Phys. 32, 179 (1960)
2840 : ! - M. F. Guest and V. R. Saunders, Mol. Phys. 28(3), 819 (1974)
2841 : ! - M. Filatov and S. Shaik, Chem. Phys. Lett. 288, 689 (1998)
2842 :
2843 : TYPE(qs_scf_env_type), POINTER :: scf_env
2844 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
2845 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
2846 : TYPE(scf_control_type), POINTER :: scf_control
2847 : TYPE(section_vals_type), POINTER :: scf_section
2848 : LOGICAL, INTENT(INOUT) :: diis_step
2849 : LOGICAL, INTENT(IN) :: orthogonal_basis
2850 :
2851 : CHARACTER(LEN=*), PARAMETER :: routineN = 'do_roks_diag'
2852 :
2853 : INTEGER :: handle, homoa, homob, imo, nalpha, nao, &
2854 : nbeta, nmo
2855 : REAL(KIND=dp) :: diis_error, level_shift_loc
2856 1132 : REAL(KIND=dp), DIMENSION(:), POINTER :: eiga, eigb, occa, occb
2857 : TYPE(cp_fm_type), POINTER :: ksa, ksb, mo2ao, moa, mob, ortho, work
2858 :
2859 : ! -------------------------------------------------------------------------
2860 :
2861 1132 : CALL timeset(routineN, handle)
2862 :
2863 1132 : IF (scf_env%cholesky_method == cholesky_inverse) THEN
2864 0 : ortho => scf_env%ortho_m1
2865 : ELSE
2866 1132 : ortho => scf_env%ortho
2867 : END IF
2868 1132 : work => scf_env%scf_work2
2869 :
2870 1132 : ksa => scf_env%scf_work1(1)
2871 1132 : ksb => scf_env%scf_work1(2)
2872 :
2873 1132 : CALL copy_dbcsr_to_fm(matrix_ks(1)%matrix, ksa)
2874 1132 : CALL copy_dbcsr_to_fm(matrix_ks(2)%matrix, ksb)
2875 :
2876 : ! Get MO information
2877 :
2878 : CALL get_mo_set(mo_set=mos(1), &
2879 : nao=nao, &
2880 : nmo=nmo, &
2881 : nelectron=nalpha, &
2882 : homo=homoa, &
2883 : eigenvalues=eiga, &
2884 : occupation_numbers=occa, &
2885 1132 : mo_coeff=moa)
2886 :
2887 : CALL get_mo_set(mo_set=mos(2), &
2888 : nelectron=nbeta, &
2889 : homo=homob, &
2890 : eigenvalues=eigb, &
2891 : occupation_numbers=occb, &
2892 1132 : mo_coeff=mob)
2893 :
2894 : ! Define the amount of level-shifting
2895 :
2896 1132 : IF ((scf_control%level_shift /= 0.0_dp) .AND. &
2897 : ((scf_control%density_guess == core_guess) .OR. &
2898 : (scf_control%density_guess == restart_guess) .OR. &
2899 : (scf_env%iter_count > 1))) THEN
2900 20 : level_shift_loc = scf_control%level_shift
2901 : ELSE
2902 1112 : level_shift_loc = 0.0_dp
2903 : END IF
2904 :
2905 : IF ((scf_env%iter_count > 1) .OR. &
2906 1132 : (scf_control%density_guess == core_guess) .OR. &
2907 : (scf_control%density_guess == restart_guess)) THEN
2908 :
2909 : ! Transform the spin unrestricted alpha and beta Kohn-Sham matrices
2910 : ! from AO basis to MO basis: K(MO) = C(T)*K(AO)*C
2911 :
2912 1024 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, moa, 0.0_dp, work)
2913 1024 : CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, moa, work, 0.0_dp, ksa)
2914 :
2915 1024 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksb, moa, 0.0_dp, work)
2916 1024 : CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, moa, work, 0.0_dp, ksb)
2917 :
2918 : ! Combine the spin unrestricted alpha and beta Kohn-Sham matrices
2919 : ! in the MO basis
2920 :
2921 1024 : IF (scf_control%roks_scheme == general_roks) THEN
2922 : CALL combine_ks_matrices(ksa, ksb, occa, occb, scf_control%roks_f, &
2923 0 : nalpha, nbeta)
2924 1024 : ELSE IF (scf_control%roks_scheme == high_spin_roks) THEN
2925 1024 : CALL combine_ks_matrices(ksa, ksb, occa, occb, scf_control%roks_parameter)
2926 : ELSE
2927 0 : CPABORT("Unknown ROKS scheme requested")
2928 : END IF
2929 :
2930 : ! Back-transform the restricted open Kohn-Sham matrix from MO basis
2931 : ! to AO basis
2932 :
2933 1024 : IF (orthogonal_basis) THEN
2934 : ! Q = C
2935 464 : mo2ao => moa
2936 : ELSE
2937 : ! Q = S*C
2938 560 : mo2ao => mob
2939 : !MK CALL copy_sm_to_fm(matrix_s(1)%matrix,work)
2940 : !MK CALL cp_fm_symm("L", "U",nao, nao, 1.0_dp, work, moa, 0.0_dp, mo2ao)
2941 560 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, moa, mo2ao, nao)
2942 : END IF
2943 :
2944 : ! K(AO) = Q*K(MO)*Q(T)
2945 :
2946 1024 : CALL parallel_gemm("N", "T", nao, nao, nao, 1.0_dp, ksa, mo2ao, 0.0_dp, work)
2947 1024 : CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, mo2ao, work, 0.0_dp, ksa)
2948 :
2949 : ELSE
2950 :
2951 : ! No transformation matrix available, yet. The closed shell part,
2952 : ! i.e. the beta Kohn-Sham matrix in AO basis, is taken.
2953 : ! There might be better choices, anyhow.
2954 :
2955 108 : CALL cp_fm_to_fm(ksb, ksa)
2956 :
2957 : END IF
2958 :
2959 : ! Update DIIS buffer and possibly perform DIIS extrapolation step
2960 :
2961 1132 : IF (scf_env%iter_count > 1) THEN
2962 1018 : IF (orthogonal_basis) THEN
2963 : CALL qs_diis_b_step(diis_buffer=scf_env%scf_diis_buffer, &
2964 : mo_array=mos, &
2965 : kc=scf_env%scf_work1, &
2966 : sc=work, &
2967 : delta=scf_env%iter_delta, &
2968 : error_max=diis_error, &
2969 : diis_step=diis_step, &
2970 : eps_diis=scf_control%eps_diis, &
2971 : scf_section=scf_section, &
2972 460 : roks=.TRUE.)
2973 460 : CPASSERT(scf_env%iter_delta == scf_env%iter_delta)
2974 : ELSE
2975 : CALL qs_diis_b_step(diis_buffer=scf_env%scf_diis_buffer, &
2976 : mo_array=mos, &
2977 : kc=scf_env%scf_work1, &
2978 : sc=work, &
2979 : delta=scf_env%iter_delta, &
2980 : error_max=diis_error, &
2981 : diis_step=diis_step, &
2982 : eps_diis=scf_control%eps_diis, &
2983 : scf_section=scf_section, &
2984 : s_matrix=matrix_s, &
2985 558 : roks=.TRUE.)
2986 : END IF
2987 : END IF
2988 :
2989 1132 : IF (diis_step) THEN
2990 714 : scf_env%iter_param = diis_error
2991 714 : scf_env%iter_method = "DIIS/Diag."
2992 : ELSE
2993 418 : IF (scf_env%mixing_method == 1) THEN
2994 418 : scf_env%iter_param = scf_env%p_mix_alpha
2995 418 : scf_env%iter_method = "P_Mix/Diag."
2996 0 : ELSE IF (scf_env%mixing_method > 1) THEN
2997 0 : scf_env%iter_param = scf_env%mixing_store%alpha
2998 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Diag."
2999 : END IF
3000 : END IF
3001 :
3002 1132 : scf_env%iter_delta = 0.0_dp
3003 :
3004 1132 : IF (level_shift_loc /= 0.0_dp) THEN
3005 :
3006 : ! Transform the current Kohn-Sham matrix from AO to MO basis
3007 : ! for level-shifting using the current MO set
3008 :
3009 20 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, moa, 0.0_dp, work)
3010 20 : CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, moa, work, 0.0_dp, ksa)
3011 :
3012 : ! Apply level-shifting using 50:50 split of the shift (could be relaxed)
3013 :
3014 60 : DO imo = homob + 1, homoa
3015 60 : CALL cp_fm_add_to_element(ksa, imo, imo, 0.5_dp*level_shift_loc)
3016 : END DO
3017 220 : DO imo = homoa + 1, nmo
3018 220 : CALL cp_fm_add_to_element(ksa, imo, imo, level_shift_loc)
3019 : END DO
3020 :
3021 1112 : ELSE IF (.NOT. orthogonal_basis) THEN
3022 :
3023 : ! Transform the current Kohn-Sham matrix to an orthogonal basis
3024 590 : SELECT CASE (scf_env%cholesky_method)
3025 : CASE (cholesky_reduce)
3026 0 : CALL cp_fm_cholesky_reduce(ksa, ortho)
3027 : CASE (cholesky_restore)
3028 526 : CALL cp_fm_uplo_to_full(ksa, work)
3029 : CALL cp_fm_cholesky_restore(ksa, nao, ortho, work, &
3030 526 : "SOLVE", pos="RIGHT")
3031 : CALL cp_fm_cholesky_restore(work, nao, ortho, ksa, &
3032 526 : "SOLVE", pos="LEFT", transa="T")
3033 : CASE (cholesky_inverse)
3034 0 : CALL cp_fm_uplo_to_full(ksa, work)
3035 : CALL cp_fm_cholesky_restore(ksa, nao, ortho, work, &
3036 0 : "MULTIPLY", pos="RIGHT")
3037 : CALL cp_fm_cholesky_restore(work, nao, ortho, ksa, &
3038 0 : "MULTIPLY", pos="LEFT", transa="T")
3039 : CASE (cholesky_off)
3040 64 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, ortho, 0.0_dp, work)
3041 654 : CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, ortho, work, 0.0_dp, ksa)
3042 : END SELECT
3043 :
3044 : END IF
3045 :
3046 : ! Diagonalization of the ROKS operator matrix
3047 :
3048 1132 : CALL choose_eigv_solver(ksa, work, eiga)
3049 :
3050 : ! Back-transformation of the orthonormal eigenvectors if needed
3051 :
3052 1132 : IF (level_shift_loc /= 0.0_dp) THEN
3053 : ! Use old MO set for back-transformation if level-shifting was applied
3054 20 : CALL cp_fm_to_fm(moa, ortho)
3055 20 : CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, ortho, work, 0.0_dp, moa)
3056 : ELSE
3057 1112 : IF (orthogonal_basis) THEN
3058 522 : CALL cp_fm_to_fm(work, moa)
3059 : ELSE
3060 1116 : SELECT CASE (scf_env%cholesky_method)
3061 : CASE (cholesky_reduce, cholesky_restore)
3062 526 : CALL cp_fm_cholesky_restore(work, nmo, ortho, moa, "SOLVE")
3063 : CASE (cholesky_inverse)
3064 0 : CALL cp_fm_cholesky_restore(work, nmo, ortho, moa, "MULTIPLY")
3065 : CASE (cholesky_off)
3066 590 : CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, ortho, work, 0.0_dp, moa)
3067 : END SELECT
3068 : END IF
3069 : END IF
3070 :
3071 : ! Correct MO eigenvalues, if level-shifting was applied
3072 :
3073 1132 : IF (level_shift_loc /= 0.0_dp) THEN
3074 60 : DO imo = homob + 1, homoa
3075 60 : eiga(imo) = eiga(imo) - 0.5_dp*level_shift_loc
3076 : END DO
3077 220 : DO imo = homoa + 1, nmo
3078 220 : eiga(imo) = eiga(imo) - level_shift_loc
3079 : END DO
3080 : END IF
3081 :
3082 : ! Update also the beta MO set
3083 :
3084 34952 : eigb(:) = eiga(:)
3085 1132 : CALL cp_fm_to_fm(moa, mob)
3086 :
3087 : ! Calculate the new alpha and beta density matrix
3088 :
3089 : ! does not yet handle k-points
3090 1132 : CALL calculate_density_matrix(mos(1), scf_env%p_mix_new(1, 1)%matrix)
3091 1132 : CALL calculate_density_matrix(mos(2), scf_env%p_mix_new(2, 1)%matrix)
3092 :
3093 1132 : CALL timestop(handle)
3094 :
3095 1132 : END SUBROUTINE do_roks_diag
3096 :
3097 : ! **************************************************************************************************
3098 : !> \brief iterative diagonalization using the block Krylov-space approach
3099 : !> \param scf_env ...
3100 : !> \param mos ...
3101 : !> \param matrix_ks ...
3102 : !> \param scf_control ...
3103 : !> \param scf_section ...
3104 : !> \param check_moconv_only ...
3105 : !> \param
3106 : !> \par History
3107 : !> 05.2009 created [MI]
3108 : ! **************************************************************************************************
3109 :
3110 38 : SUBROUTINE do_block_krylov_diag(scf_env, mos, matrix_ks, &
3111 : scf_control, scf_section, check_moconv_only)
3112 :
3113 : TYPE(qs_scf_env_type), POINTER :: scf_env
3114 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
3115 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
3116 : TYPE(scf_control_type), POINTER :: scf_control
3117 : TYPE(section_vals_type), POINTER :: scf_section
3118 : LOGICAL, INTENT(IN), OPTIONAL :: check_moconv_only
3119 :
3120 : CHARACTER(LEN=*), PARAMETER :: routineN = 'do_block_krylov_diag'
3121 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
3122 :
3123 : INTEGER :: handle, homo, ispin, iter, nao, nmo, &
3124 : output_unit
3125 : LOGICAL :: converged, my_check_moconv_only
3126 : REAL(dp) :: eps_iter, t1, t2
3127 38 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
3128 : TYPE(cp_fm_type), POINTER :: c0, c1, chc, evec, ks, mo_coeff, ortho, &
3129 : work
3130 : TYPE(cp_logger_type), POINTER :: logger
3131 :
3132 76 : logger => cp_get_default_logger()
3133 38 : CALL timeset(routineN, handle)
3134 :
3135 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%LANCZOS", &
3136 38 : extension=".scfLog")
3137 :
3138 38 : my_check_moconv_only = .FALSE.
3139 38 : IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
3140 :
3141 38 : NULLIFY (mo_coeff, ortho, work, ks)
3142 38 : NULLIFY (mo_eigenvalues)
3143 38 : NULLIFY (c0, c1)
3144 :
3145 38 : IF (scf_env%cholesky_method == cholesky_inverse) THEN
3146 38 : ortho => scf_env%ortho_m1
3147 : ELSE
3148 0 : ortho => scf_env%ortho
3149 : END IF
3150 38 : work => scf_env%scf_work2
3151 :
3152 76 : DO ispin = 1, SIZE(matrix_ks)
3153 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, &
3154 76 : scf_env%scf_work1(ispin))
3155 : END DO
3156 :
3157 38 : IF (scf_env%mixing_method == 1) THEN
3158 0 : scf_env%iter_param = scf_env%p_mix_alpha
3159 0 : scf_env%iter_method = "P_Mix/Lanczos"
3160 : ELSE
3161 : ! scf_env%iter_param = scf_env%mixing_store%alpha
3162 38 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Lanc."
3163 : END IF
3164 :
3165 76 : DO ispin = 1, SIZE(matrix_ks)
3166 :
3167 38 : ks => scf_env%scf_work1(ispin)
3168 38 : CALL cp_fm_uplo_to_full(ks, work)
3169 :
3170 : CALL get_mo_set(mo_set=mos(ispin), &
3171 : nao=nao, &
3172 : nmo=nmo, &
3173 : homo=homo, &
3174 : eigenvalues=mo_eigenvalues, &
3175 38 : mo_coeff=mo_coeff)
3176 :
3177 38 : NULLIFY (c0, c1)
3178 38 : c0 => scf_env%krylov_space%mo_conv(ispin)
3179 38 : c1 => scf_env%krylov_space%mo_refine(ispin)
3180 38 : SELECT CASE (scf_env%cholesky_method)
3181 : CASE (cholesky_reduce)
3182 0 : CALL cp_fm_cholesky_reduce(ks, ortho)
3183 0 : CALL cp_fm_uplo_to_full(ks, work)
3184 0 : CALL cp_fm_cholesky_restore(mo_coeff, nmo, ortho, c0, "MULTIPLY")
3185 : CASE (cholesky_restore)
3186 : CALL cp_fm_cholesky_restore(ks, nao, ortho, work, &
3187 0 : "SOLVE", pos="RIGHT")
3188 : CALL cp_fm_cholesky_restore(work, nao, ortho, ks, &
3189 0 : "SOLVE", pos="LEFT", transa="T")
3190 0 : CALL cp_fm_cholesky_restore(mo_coeff, nmo, ortho, c0, "MULTIPLY")
3191 : CASE (cholesky_inverse)
3192 : CALL cp_fm_cholesky_restore(ks, nao, ortho, work, &
3193 38 : "MULTIPLY", pos="RIGHT")
3194 : CALL cp_fm_cholesky_restore(work, nao, ortho, ks, &
3195 38 : "MULTIPLY", pos="LEFT", transa="T")
3196 76 : CALL cp_fm_cholesky_restore(mo_coeff, nmo, ortho, c0, "SOLVE")
3197 : END SELECT
3198 :
3199 38 : scf_env%krylov_space%nmo_nc = nmo
3200 38 : scf_env%krylov_space%nmo_conv = 0
3201 :
3202 38 : t1 = m_walltime()
3203 38 : IF (output_unit > 0) THEN
3204 0 : WRITE (output_unit, "(/T15,A)") '<<<<<<<<< LANCZOS REFINEMENT <<<<<<<<<<'
3205 : WRITE (output_unit, "(T8,A,T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") &
3206 0 : " Spin ", " Cycle ", &
3207 0 : " conv. MOS ", " B2MAX ", " B2MIN ", " Time", REPEAT("-", 60)
3208 : END IF
3209 38 : eps_iter = MAX(scf_env%krylov_space%eps_conv, scf_env%krylov_space%eps_adapt*scf_env%iter_delta)
3210 38 : iter = 0
3211 38 : converged = .FALSE.
3212 : !Check convergence of MOS
3213 114 : IF (my_check_moconv_only) THEN
3214 :
3215 : CALL lanczos_refinement(scf_env%krylov_space, ks, c0, c1, mo_eigenvalues, &
3216 0 : nao, eps_iter, ispin, check_moconv_only=my_check_moconv_only)
3217 0 : t2 = m_walltime()
3218 0 : IF (output_unit > 0) THEN
3219 : WRITE (output_unit, '(T8,I3,T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
3220 0 : ispin, iter, scf_env%krylov_space%nmo_conv, &
3221 0 : scf_env%krylov_space%max_res_norm, scf_env%krylov_space%min_res_norm, t2 - t1
3222 : END IF
3223 :
3224 : CYCLE
3225 : ELSE
3226 : !Block Lanczos refinement
3227 620 : DO iter = 1, scf_env%krylov_space%max_iter
3228 : CALL lanczos_refinement_2v(scf_env%krylov_space, ks, c0, c1, mo_eigenvalues, &
3229 592 : nao, eps_iter, ispin)
3230 592 : t2 = m_walltime()
3231 592 : IF (output_unit > 0) THEN
3232 : WRITE (output_unit, '(T8,I3,T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
3233 0 : ispin, iter, scf_env%krylov_space%nmo_conv, &
3234 0 : scf_env%krylov_space%max_res_norm, scf_env%krylov_space%min_res_norm, t2 - t1
3235 : END IF
3236 592 : t1 = m_walltime()
3237 620 : IF (scf_env%krylov_space%max_res_norm < eps_iter) THEN
3238 10 : converged = .TRUE.
3239 10 : IF (output_unit > 0) WRITE (output_unit, *) &
3240 0 : " Reached convergence in ", iter, " iterations "
3241 : EXIT
3242 : END IF
3243 : END DO
3244 :
3245 38 : IF (.NOT. converged .AND. output_unit > 0) THEN
3246 : WRITE (output_unit, "(T4, A)") " WARNING Lanczos refinement could "// &
3247 0 : "not converge all the mos:"
3248 0 : WRITE (output_unit, "(T40,A,T70,I10)") " number of not converged mos ", &
3249 0 : scf_env%krylov_space%nmo_nc
3250 0 : WRITE (output_unit, "(T40,A,T70,E10.2)") " max norm of the residual ", &
3251 0 : scf_env%krylov_space%max_res_norm
3252 :
3253 : END IF
3254 :
3255 : ! For the moment skip the re-orthogonalization
3256 : IF (.FALSE.) THEN
3257 : !Re-orthogonalization
3258 : NULLIFY (chc, evec)
3259 : chc => scf_env%krylov_space%chc_mat(ispin)
3260 : evec => scf_env%krylov_space%c_vec(ispin)
3261 : CALL parallel_gemm('N', 'N', nao, nmo, nao, rone, ks, c0, rzero, work)
3262 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, rone, c0, work, rzero, chc)
3263 : !Diagonalize (C^t)HC
3264 : CALL choose_eigv_solver(chc, evec, mo_eigenvalues)
3265 : !Rotate the C vectors
3266 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, rone, c0, evec, rzero, c1)
3267 : c0 => scf_env%krylov_space%mo_refine(ispin)
3268 : END IF
3269 :
3270 38 : IF (scf_env%cholesky_method == cholesky_inverse) THEN
3271 38 : CALL cp_fm_cholesky_restore(c0, nmo, ortho, mo_coeff, "MULTIPLY")
3272 : ELSE
3273 0 : CALL cp_fm_cholesky_restore(c0, nmo, ortho, mo_coeff, "SOLVE")
3274 : END IF
3275 :
3276 38 : IF (.NOT. scf_control%gce%do_gce) THEN
3277 : CALL set_mo_occupation(mo_set=mos(ispin), &
3278 38 : smear=scf_control%smear)
3279 : ELSE
3280 : CALL set_mo_occupation(mo_set=mos(ispin), &
3281 : smear=scf_control%smear, &
3282 0 : gce=scf_control%gce)
3283 : END IF
3284 :
3285 : ! does not yet handle k-points
3286 : CALL calculate_density_matrix(mos(ispin), &
3287 38 : scf_env%p_mix_new(ispin, 1)%matrix)
3288 : END IF
3289 : END DO ! ispin
3290 :
3291 38 : IF (output_unit > 0) THEN
3292 0 : WRITE (output_unit, "(T15,A/)") '<<<<<<<<< END LANCZOS REFINEMENT <<<<<<<<<<'
3293 : END IF
3294 :
3295 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
3296 38 : "PRINT%LANCZOS")
3297 :
3298 38 : CALL timestop(handle)
3299 :
3300 38 : END SUBROUTINE do_block_krylov_diag
3301 :
3302 : ! **************************************************************************************************
3303 : !> \brief iterative diagonalization using the block davidson space approach
3304 : !> \param qs_env ...
3305 : !> \param scf_env ...
3306 : !> \param mos ...
3307 : !> \param matrix_ks ...
3308 : !> \param matrix_s ...
3309 : !> \param scf_control ...
3310 : !> \param scf_section ...
3311 : !> \param check_moconv_only ...
3312 : !> \param
3313 : !> \par History
3314 : !> 05.2011 created [MI]
3315 : ! **************************************************************************************************
3316 :
3317 94 : SUBROUTINE do_block_davidson_diag(qs_env, scf_env, mos, matrix_ks, matrix_s, &
3318 : scf_control, scf_section, check_moconv_only)
3319 :
3320 : TYPE(qs_environment_type), POINTER :: qs_env
3321 : TYPE(qs_scf_env_type), POINTER :: scf_env
3322 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
3323 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
3324 : TYPE(scf_control_type), POINTER :: scf_control
3325 : TYPE(section_vals_type), POINTER :: scf_section
3326 : LOGICAL, INTENT(IN), OPTIONAL :: check_moconv_only
3327 :
3328 : CHARACTER(LEN=*), PARAMETER :: routineN = 'do_block_davidson_diag'
3329 :
3330 : INTEGER :: handle, ispin, nspins, output_unit
3331 : LOGICAL :: do_prec, my_check_moconv_only
3332 : TYPE(cp_logger_type), POINTER :: logger
3333 :
3334 94 : logger => cp_get_default_logger()
3335 94 : CALL timeset(routineN, handle)
3336 :
3337 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DAVIDSON", &
3338 94 : extension=".scfLog")
3339 :
3340 94 : IF (output_unit > 0) THEN
3341 0 : WRITE (output_unit, "(/T15,A)") '<<<<<<<<< DAVIDSON ITERATIONS <<<<<<<<<<'
3342 : END IF
3343 :
3344 94 : IF (scf_env%mixing_method == 1) THEN
3345 0 : scf_env%iter_param = scf_env%p_mix_alpha
3346 0 : scf_env%iter_method = "P_Mix/Dav."
3347 : ELSE
3348 94 : scf_env%iter_param = scf_env%mixing_store%alpha
3349 94 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Dav."
3350 : END IF
3351 :
3352 94 : my_check_moconv_only = .FALSE.
3353 94 : IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
3354 94 : do_prec = .FALSE.
3355 94 : IF (scf_env%block_davidson_env(1)%prec_type /= 0 .AND. &
3356 : scf_env%iter_count >= scf_env%block_davidson_env(1)%first_prec) THEN
3357 84 : do_prec = .TRUE.
3358 : END IF
3359 :
3360 94 : nspins = SIZE(matrix_ks)
3361 :
3362 94 : IF (do_prec .AND. (scf_env%iter_count == scf_env%block_davidson_env(1)%first_prec .OR. &
3363 : MODULO(scf_env%iter_count, scf_env%block_davidson_env(1)%niter_new_prec) == 0)) THEN
3364 : CALL restart_preconditioner(qs_env, scf_env%ot_preconditioner, &
3365 18 : prec_type=scf_env%block_davidson_env(1)%prec_type, nspins=nspins)
3366 : CALL prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, scf_env%ot_preconditioner, &
3367 : scf_env%block_davidson_env(1)%prec_type, &
3368 : scf_env%block_davidson_env(1)%solver_type, &
3369 : scf_env%block_davidson_env(1)%energy_gap, nspins, &
3370 : convert_to_dbcsr=scf_env%block_davidson_env(1)%use_sparse_mos, &
3371 18 : full_mo_set=.TRUE.)
3372 : END IF
3373 :
3374 198 : DO ispin = 1, nspins
3375 198 : IF (scf_env%block_davidson_env(ispin)%use_sparse_mos) THEN
3376 64 : IF (.NOT. do_prec) THEN
3377 : CALL generate_extended_space_sparse(scf_env%block_davidson_env(ispin), mos(ispin), &
3378 8 : matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit)
3379 : ELSE
3380 : CALL generate_extended_space_sparse(scf_env%block_davidson_env(ispin), mos(ispin), &
3381 : matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit, &
3382 56 : scf_env%ot_preconditioner(ispin)%preconditioner)
3383 : END IF
3384 :
3385 : ELSE
3386 40 : IF (.NOT. do_prec) THEN
3387 : CALL generate_extended_space(scf_env%block_davidson_env(ispin), mos(ispin), &
3388 4 : matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit)
3389 : ELSE
3390 : CALL generate_extended_space(scf_env%block_davidson_env(ispin), mos(ispin), &
3391 : matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit, &
3392 36 : scf_env%ot_preconditioner(ispin)%preconditioner)
3393 : END IF
3394 : END IF
3395 : END DO !ispin
3396 :
3397 94 : IF (.NOT. scf_control%gce%do_gce) THEN
3398 : CALL set_mo_occupation(mo_array=mos, &
3399 94 : smear=scf_control%smear)
3400 : ELSE
3401 : CALL set_mo_occupation(mo_array=mos, &
3402 : smear=scf_control%smear, &
3403 0 : gce=scf_control%gce)
3404 : END IF
3405 :
3406 198 : DO ispin = 1, nspins
3407 : ! does not yet handle k-points
3408 : CALL calculate_density_matrix(mos(ispin), &
3409 198 : scf_env%p_mix_new(ispin, 1)%matrix)
3410 : END DO
3411 :
3412 94 : IF (output_unit > 0) THEN
3413 0 : WRITE (output_unit, "(T15,A/)") '<<<<<<<<< END DAVIDSON ITERATION <<<<<<<<<<'
3414 : END IF
3415 :
3416 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
3417 94 : "PRINT%DAVIDSON")
3418 :
3419 94 : CALL timestop(handle)
3420 :
3421 94 : END SUBROUTINE do_block_davidson_diag
3422 :
3423 : ! **************************************************************************************************
3424 : !> \brief block-Davidson SCF step for complex K points, one channel per
3425 : !> (local kpoint, spin); mirrors do_block_davidson_diag for the
3426 : !> preconditioner scheduling and do_general_diag_kp for the H(k)/S(k)
3427 : !> assembly and the density exit chain
3428 : !> \param qs_env ...
3429 : !> \param scf_env ...
3430 : !> \param matrix_ks_kp real-space Kohn-Sham matrices, dbcsr image set
3431 : !> \param matrix_s_kp real-space overlap matrix, dbcsr image set
3432 : !> \param kpoints ...
3433 : !> \param scf_control ...
3434 : ! **************************************************************************************************
3435 932 : SUBROUTINE do_block_davidson_diag_kp(qs_env, scf_env, matrix_ks_kp, matrix_s_kp, kpoints, &
3436 : scf_control)
3437 :
3438 : TYPE(qs_environment_type), POINTER :: qs_env
3439 : TYPE(qs_scf_env_type), POINTER :: scf_env
3440 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
3441 : TYPE(kpoint_type), POINTER :: kpoints
3442 : TYPE(scf_control_type), POINTER :: scf_control
3443 :
3444 : CHARACTER(LEN=*), PARAMETER :: routineN = 'do_block_davidson_diag_kp'
3445 : REAL(KIND=dp), PARAMETER :: empty_factor = 5.0_dp, &
3446 : eps_iter_cap = 1.0E-4_dp
3447 :
3448 : INTEGER :: channel, first_prec, handle, ispin, iw, &
3449 : iwork, kplocal, local_kpoint, &
3450 : nchannels, nspin, output_unit
3451 : INTEGER, DIMENSION(2) :: kp_range
3452 932 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
3453 932 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3454 : LOGICAL :: do_prec, use_real_wfn
3455 : REAL(KIND=dp) :: eps_iter, eps_iter_empty
3456 932 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
3457 : TYPE(cp_cfm_type) :: cksmat, csmat, csmat_base
3458 932 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools, kp_ao_ao_fm_pools
3459 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
3460 : TYPE(cp_fm_type) :: fmlocal_im, fmlocal_re
3461 932 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
3462 : TYPE(cp_logger_type), POINTER :: logger
3463 : TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix, tmpmat
3464 : TYPE(kpoint_env_type), POINTER :: kp
3465 : TYPE(mp_para_env_type), POINTER :: para_env
3466 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3467 932 : POINTER :: sab_nl
3468 : TYPE(qs_matrix_pools_type), POINTER :: kp_mpools, mpools
3469 : TYPE(section_vals_type), POINTER :: scf_section
3470 :
3471 : ! bounds of the adaptive inner tolerance. The loose phase is capped for
3472 : ! the residual-norm semantics of EPS_ITER. A residual of order 1/k would
3473 : ! leave an eigenvector error of order 10 percent at a typical metallic
3474 : ! gap. Unoccupied MOS converge empty_factor times looser. The tracking
3475 : ! divides the mixing delta by noise_scale: the SCF convergence measure
3476 : ! is the LARGEST ELEMENT of the density residual, in which the
3477 : ! per-column residuals of a whole block accumulate. The measured
3478 : ! amplification is ~20 on a smeared Bi2S3/MoS2 slab, so tracking at
3479 : ! delta/100 keeps the injected noise below the outer convergence level.
3480 :
3481 1864 : logger => cp_get_default_logger()
3482 932 : CALL timeset(routineN, handle)
3483 :
3484 932 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
3485 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DAVIDSON", &
3486 932 : extension=".scfLog")
3487 932 : IF (output_unit > 0) THEN
3488 0 : WRITE (output_unit, "(/T15,A)") '<<<<<<<<< DAVIDSON ITERATIONS <<<<<<<<<<'
3489 : END IF
3490 :
3491 : CALL get_kpoint_info(kpoints, kp_range=kp_range, kp_dist=kp_dist, xkp=xkp, &
3492 : use_real_wfn=use_real_wfn, sab_nl=sab_nl, &
3493 932 : cell_to_index=cell_to_index, mpools=kp_mpools)
3494 932 : IF (use_real_wfn) THEN
3495 0 : CPABORT('Davidson Kp requires complex K-point MOS')
3496 : END IF
3497 932 : para_env => kpoints%blacs_env_all%para_env
3498 932 : kplocal = kp_range(2) - kp_range(1) + 1
3499 932 : nspin = SIZE(matrix_ks_kp, 1)
3500 932 : nchannels = kplocal*nspin
3501 :
3502 932 : IF (scf_env%mixing_method == 1) THEN
3503 684 : scf_env%iter_param = scf_env%p_mix_alpha
3504 684 : scf_env%iter_method = "P_Mix/Dav."
3505 248 : ELSE IF (scf_env%mixing_method > 1) THEN
3506 248 : scf_env%iter_param = scf_env%mixing_store%alpha
3507 248 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Dav."
3508 : END IF
3509 :
3510 : ! preconditioner scheduling, mirrors do_block_davidson_diag. The first
3511 : ! K-point SCF step cold-starts unless extrapolation prefilled the MOS
3512 : ! (mos_prefilled). The FULL_ALL maker consumes the MO coefficients and
3513 : ! eigenvalues, so building before the first solve completed would feed it
3514 : ! uninitialized data and stall the Davidson convergence. Keep first_prec
3515 : ! clamped past that step in both cases. Warming it up under prefilled MOS
3516 : ! would mix two variables at once.
3517 932 : first_prec = MAX(scf_env%block_davidson_env(1)%first_prec, 2)
3518 932 : do_prec = .FALSE.
3519 932 : IF (scf_env%block_davidson_env(1)%prec_type /= ot_precond_none .AND. &
3520 : scf_env%iter_count >= first_prec) THEN
3521 838 : do_prec = .TRUE.
3522 : END IF
3523 : ! adaptive inner tolerance from EPS_ADAPT: loose while the density is
3524 : ! still far from self-consistency, tightening with the outer error
3525 : ! scaled by noise_scale (mirrors the block Krylov scheduling in
3526 : ! do_block_krylov_diag). The
3527 : ! ratchet is monotone across SCF steps: the threshold may tighten but
3528 : ! never relax again. A rebounding outer error would otherwise re-loosen
3529 : ! the solver, and the run could stop on a near-zero mixing delta of a
3530 : ! pair of identical, sloppily converged densities. The cap bounds the
3531 : ! early loose phase. Unoccupied MOS converge empty_factor times
3532 : ! looser (the occupancy decides which columns qualify, see
3533 : ! generate_extended_space_c). There is no restoration of the
3534 : ! user threshold at the end. The tolerance stays a margin below the
3535 : ! outer error. The MOS the run exits with satisfy that same margin
3536 : ! rather than the nominal EPS_ITER.
3537 932 : eps_iter = scf_env%block_davidson_env(1)%eps_iter
3538 932 : eps_iter_empty = eps_iter
3539 932 : IF (scf_control%diagonalization%eps_adapt > 0.0_dp .AND. scf_env%iter_delta > 0.0_dp) THEN
3540 : ! stall escape: when the mixing delta sank below the tolerance the
3541 : ! last solve used, that solve's own sloppiness may be the noise floor
3542 : ! pinning the delta. Tighten one decade below what that solve used.
3543 : ! Tightening relative to the OPERATING tolerance instead of the
3544 : ! dipped delta bounds the reaction to transient mixing dips (a
3545 : ! momentary overshoot must not collapse the threshold by orders of
3546 : ! magnitude and destabilise the mixer), while a genuine stall moves
3547 : ! the threshold down decade by decade. Checking on the stall alone
3548 : ! keeps the loose phase as loose as EPS_ADAPT allows
3549 172 : IF (scf_env%iter_delta < scf_env%block_davidson_env(1)%eps_iter_used) THEN
3550 : scf_env%block_davidson_env(1)%eps_iter_prev = &
3551 : MIN(scf_env%block_davidson_env(1)%eps_iter_prev, &
3552 : MAX(scf_env%block_davidson_env(1)%eps_iter, &
3553 2 : 0.1_dp*scf_env%block_davidson_env(1)%eps_iter_used))
3554 : END IF
3555 : ! candidate: looser than the user value, capped; ratchet: only
3556 : ! tighter than the previous adaptive value, which lives in the
3557 : ! persistent channel settings (a local variable would reset on
3558 : ! every call into this driver)
3559 : eps_iter = MIN(MAX(eps_iter, MIN(scf_control%diagonalization%eps_adapt*scf_env%iter_delta/ &
3560 : scf_env%block_davidson_env(1)%noise_scale, &
3561 172 : eps_iter_cap)), scf_env%block_davidson_env(1)%eps_iter_prev)
3562 172 : eps_iter_empty = empty_factor*eps_iter
3563 : ! record what this solve uses, after every override above: the outer
3564 : ! loop compares the mixing delta against it (qs_scf_check_inner_exit)
3565 172 : scf_env%block_davidson_env(1)%eps_iter_prev = eps_iter
3566 172 : scf_env%block_davidson_env(1)%eps_iter_used = eps_iter
3567 : END IF
3568 : ! schedule trace of the adaptive tolerance, one line per SCF step at
3569 : ! DEBUG print level
3570 932 : IF (scf_control%diagonalization%eps_adapt > 0.0_dp .AND. para_env%is_source() .AND. &
3571 : logger%iter_info%print_level >= debug_print_level) THEN
3572 0 : iw = cp_logger_get_default_unit_nr()
3573 0 : IF (iw > 0) THEN
3574 : WRITE (UNIT=iw, FMT="(T3,'EPS_ADAPT step',I5,' delta=',ES9.2,' eps_occ=',ES9.2,' eps_empty=',ES9.2)") &
3575 0 : scf_env%iter_count, scf_env%iter_delta, eps_iter, eps_iter_empty
3576 : END IF
3577 : END IF
3578 :
3579 : ! work matrices: real-space Fourier scratch (dbcsr), global scratch fms,
3580 : ! two k-group fms as redistribution targets, and the complex H(k)/S(k)
3581 932 : ALLOCATE (rmatrix, cmatrix, tmpmat)
3582 : CALL dbcsr_create(rmatrix, template=matrix_ks_kp(1, 1)%matrix, &
3583 932 : matrix_type=dbcsr_type_symmetric)
3584 : CALL dbcsr_create(cmatrix, template=matrix_ks_kp(1, 1)%matrix, &
3585 932 : matrix_type=dbcsr_type_antisymmetric)
3586 : CALL dbcsr_create(tmpmat, template=matrix_ks_kp(1, 1)%matrix, &
3587 932 : matrix_type=dbcsr_type_no_symmetry)
3588 932 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
3589 932 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
3590 932 : CALL cp_dbcsr_alloc_block_from_nbl(tmpmat, sab_nl)
3591 :
3592 932 : CALL get_qs_env(qs_env, mpools=mpools)
3593 932 : CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
3594 4660 : ALLOCATE (fmwork(4))
3595 4660 : DO iwork = 1, 4
3596 4660 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmwork(iwork))
3597 : END DO
3598 932 : CALL mpools_get(kp_mpools, ao_ao_fm_pools=kp_ao_ao_fm_pools)
3599 932 : CALL fm_pool_create_fm(kp_ao_ao_fm_pools(1)%pool, fmlocal_re)
3600 932 : CALL fm_pool_create_fm(kp_ao_ao_fm_pools(1)%pool, fmlocal_im)
3601 932 : CALL cp_fm_get_info(fmlocal_re, matrix_struct=matrix_struct)
3602 932 : CALL cp_cfm_create(cksmat, matrix_struct, name="dav_hk")
3603 932 : CALL cp_cfm_create(csmat, matrix_struct, name="dav_sk")
3604 932 : CALL cp_cfm_create(csmat_base, matrix_struct, name="dav_sk_base")
3605 :
3606 932 : IF (do_prec .AND. (scf_env%iter_count == first_prec .OR. &
3607 : MODULO(scf_env%iter_count, &
3608 : scf_env%block_davidson_env(1)%niter_new_prec) == 0)) THEN
3609 : CALL davidson_kp_prepare_preconditioners( &
3610 : qs_env, scf_env, matrix_ks_kp, matrix_s_kp, kpoints, &
3611 58 : rmatrix, cmatrix, tmpmat, fmwork, fmlocal_re, fmlocal_im)
3612 : END IF
3613 :
3614 2936 : DO local_kpoint = 1, kplocal
3615 2004 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
3616 :
3617 : ! Every k-point group executes the same sequence of global transforms.
3618 : ! Only the group owning this local_kpoint receives the redistributed
3619 : ! matrix; the others participate through the dummy destination.
3620 : CALL distribute_kpoint_operator_channel( &
3621 : matrix_s_kp, kpoints, local_kpoint, 1, kp_dist, xkp, cell_to_index, sab_nl, &
3622 2004 : rmatrix, cmatrix, tmpmat, fmwork, fmlocal_re, fmlocal_im)
3623 2004 : CALL cp_cfm_scale_and_add_fm(z_zero, csmat_base, z_one, fmlocal_re)
3624 2004 : CALL cp_cfm_scale_and_add_fm(z_one, csmat_base, gaussi, fmlocal_im)
3625 :
3626 5016 : DO ispin = 1, nspin
3627 : CALL distribute_kpoint_operator_channel( &
3628 : matrix_ks_kp, kpoints, local_kpoint, ispin, kp_dist, xkp, cell_to_index, sab_nl, &
3629 2080 : rmatrix, cmatrix, tmpmat, fmwork, fmlocal_re, fmlocal_im)
3630 2080 : CALL cp_cfm_scale_and_add_fm(z_zero, cksmat, z_one, fmlocal_re)
3631 2080 : CALL cp_cfm_scale_and_add_fm(z_one, cksmat, gaussi, fmlocal_im)
3632 2080 : CALL cp_cfm_to_cfm(csmat_base, csmat)
3633 :
3634 4084 : IF (scf_env%iter_count == 1 .AND. .NOT. kp%mos_prefilled) THEN
3635 : ! K-point MOS enter the first SCF step allocated but uninitialized.
3636 : ! The stored S(k) copy is restored for every spin because
3637 : ! cp_cfm_geeig overwrites its matrix arguments.
3638 68 : CALL cold_start_mo_init(kp, ispin, cksmat, csmat)
3639 : ELSE
3640 2012 : channel = (local_kpoint - 1)*nspin + ispin
3641 2012 : IF (do_prec) THEN
3642 : CALL generate_extended_space_c(scf_env%block_davidson_env(channel), &
3643 : kp%mos(1:2, ispin), cksmat, csmat, output_unit, &
3644 : eps_iter, eps_iter_empty, &
3645 1884 : scf_env%ot_preconditioner(channel)%preconditioner)
3646 : ELSE
3647 : CALL generate_extended_space_c(scf_env%block_davidson_env(channel), &
3648 : kp%mos(1:2, ispin), cksmat, csmat, output_unit, &
3649 128 : eps_iter, eps_iter_empty)
3650 : END IF
3651 : END IF
3652 : END DO
3653 : END DO
3654 :
3655 : ! occupations and density, same exit chain as do_general_diag_kp
3656 932 : CALL kpoint_set_mo_occupation(kpoints, scf_control%smear)
3657 932 : CALL kpoint_density_matrices(kpoints)
3658 : CALL kpoint_density_transform(kpoints, scf_env%p_mix_new, .FALSE., &
3659 932 : matrix_s_kp(1, 1)%matrix, sab_nl, fmwork, overlap_rs=matrix_s_kp)
3660 :
3661 932 : CALL cp_cfm_release(cksmat)
3662 932 : CALL cp_cfm_release(csmat)
3663 932 : CALL cp_cfm_release(csmat_base)
3664 932 : CALL fm_pool_give_back_fm(kp_ao_ao_fm_pools(1)%pool, fmlocal_re)
3665 932 : CALL fm_pool_give_back_fm(kp_ao_ao_fm_pools(1)%pool, fmlocal_im)
3666 4660 : DO iwork = 1, 4
3667 4660 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmwork(iwork))
3668 : END DO
3669 932 : DEALLOCATE (fmwork)
3670 932 : CALL dbcsr_deallocate_matrix(rmatrix)
3671 932 : CALL dbcsr_deallocate_matrix(cmatrix)
3672 932 : CALL dbcsr_deallocate_matrix(tmpmat)
3673 :
3674 932 : IF (output_unit > 0) THEN
3675 0 : WRITE (output_unit, "(T15,A/)") '<<<<<<<<< END DAVIDSON ITERATION <<<<<<<<<<'
3676 : END IF
3677 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
3678 932 : "PRINT%DAVIDSON")
3679 :
3680 932 : CALL timestop(handle)
3681 :
3682 2796 : END SUBROUTINE do_block_davidson_diag_kp
3683 :
3684 : ! **************************************************************************************************
3685 : !> \brief first-SCF-step initialization of the K-point MOS of one spin by one
3686 : !> direct diagonalization of H(k); both arguments are assembled by the
3687 : !> caller because cp_cfm_geeig overwrites them
3688 : !> \param kp kpoint environment holding the MOS pair
3689 : !> \param ispin ...
3690 : !> \param cksmat complex H(k), assembled by the caller
3691 : !> \param csmat complex S(k), assembled by the caller
3692 : ! **************************************************************************************************
3693 204 : SUBROUTINE cold_start_mo_init(kp, ispin, cksmat, csmat)
3694 :
3695 : TYPE(kpoint_env_type), POINTER :: kp
3696 : INTEGER, INTENT(IN) :: ispin
3697 : TYPE(cp_cfm_type) :: cksmat, csmat
3698 :
3699 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cold_start_mo_init'
3700 :
3701 : INTEGER :: handle
3702 68 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
3703 : TYPE(cp_cfm_type) :: cmos, cwork
3704 : TYPE(cp_fm_type), POINTER :: imos, rmos
3705 :
3706 68 : CALL timeset(routineN, handle)
3707 :
3708 68 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
3709 68 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
3710 68 : CALL cp_cfm_create(cmos, rmos%matrix_struct, name="dav_cold_cmos")
3711 68 : CALL cp_cfm_create(cwork, cksmat%matrix_struct, name="dav_cold_cwork")
3712 68 : CALL cp_cfm_geeig(cksmat, csmat, cmos, eigenvalues, cwork)
3713 552 : kp%mos(2, ispin)%eigenvalues = eigenvalues
3714 68 : CALL cp_cfm_to_fm(cmos, rmos, imos)
3715 68 : CALL cp_cfm_release(cmos)
3716 68 : CALL cp_cfm_release(cwork)
3717 :
3718 68 : CALL timestop(handle)
3719 :
3720 68 : END SUBROUTINE cold_start_mo_init
3721 :
3722 : ! **************************************************************************************************
3723 : !> \brief builds one complex K-point preconditioner per (local kpoint, spin)
3724 : !> channel; mirrors prepare_qs_ot_kpoint_preconditioners
3725 : !> \param qs_env ...
3726 : !> \param scf_env ...
3727 : !> \param matrix_ks_kp ...
3728 : !> \param matrix_s_kp ...
3729 : !> \param kpoints ...
3730 : !> \param rmatrix ...
3731 : !> \param cmatrix ...
3732 : !> \param tmpmat ...
3733 : !> \param fmwork ...
3734 : !> \param fmlocal_re ...
3735 : !> \param fmlocal_im ...
3736 : ! **************************************************************************************************
3737 58 : SUBROUTINE davidson_kp_prepare_preconditioners( &
3738 : qs_env, scf_env, matrix_ks_kp, matrix_s_kp, kpoints, &
3739 : rmatrix, cmatrix, tmpmat, fmwork, fmlocal_re, fmlocal_im)
3740 :
3741 : TYPE(qs_environment_type), POINTER :: qs_env
3742 : TYPE(qs_scf_env_type), POINTER :: scf_env
3743 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
3744 : TYPE(kpoint_type), POINTER :: kpoints
3745 : TYPE(dbcsr_type), POINTER :: rmatrix, cmatrix, tmpmat
3746 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
3747 : TYPE(cp_fm_type) :: fmlocal_re, fmlocal_im
3748 :
3749 : CHARACTER(LEN=*), PARAMETER :: routineN = 'davidson_kp_prepare_preconditioners'
3750 :
3751 : INTEGER :: channel, handle, ispin, kplocal, &
3752 : local_kpoint, nchannels, nmo, nspin, &
3753 : prec_type, solver_type
3754 : INTEGER, DIMENSION(2) :: kp_range
3755 58 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
3756 58 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3757 : REAL(KIND=dp) :: energy_gap
3758 58 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
3759 58 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
3760 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
3761 : TYPE(cp_cfm_type) :: matrix_c, matrix_h, matrix_s
3762 : TYPE(cp_fm_type), POINTER :: mo_coeff_im, mo_coeff_re
3763 58 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_t_kp
3764 : TYPE(dbcsr_type), POINTER :: matrix_h_im, matrix_h_re, matrix_s_im, &
3765 : matrix_s_re, matrix_t_im, matrix_t_re
3766 : TYPE(kpoint_env_type), POINTER :: kp
3767 : TYPE(mp_para_env_type), POINTER :: para_env
3768 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3769 58 : POINTER :: sab_nl
3770 :
3771 58 : CALL timeset(routineN, handle)
3772 :
3773 : CALL get_kpoint_info(kpoints, kp_range=kp_range, kp_dist=kp_dist, xkp=xkp, sab_nl=sab_nl, &
3774 58 : cell_to_index=cell_to_index)
3775 58 : CALL get_qs_env(qs_env, kinetic_kp=matrix_t_kp)
3776 : ! the preconditioner fms must live on the kpoint-group BLACS context,
3777 : ! the context of the MOS they are applied with
3778 58 : blacs_env => kpoints%blacs_env
3779 58 : para_env => kpoints%para_env_kp
3780 58 : kplocal = kp_range(2) - kp_range(1) + 1
3781 58 : nspin = SIZE(matrix_ks_kp, 1)
3782 58 : nchannels = kplocal*nspin
3783 58 : prec_type = scf_env%block_davidson_env(1)%prec_type
3784 58 : solver_type = scf_env%block_davidson_env(1)%solver_type
3785 58 : energy_gap = scf_env%block_davidson_env(1)%energy_gap
3786 :
3787 : ! channel storage: every (kpoint, spin) channel carries its own
3788 : ! preconditioner, because S(k) and H(k) differ per kpoint. The Gamma
3789 : ! reuse pattern of restart_preconditioner (a single entry for the
3790 : ! spin-independent types) must not be used here.
3791 58 : IF (ASSOCIATED(scf_env%ot_preconditioner)) THEN
3792 128 : DO channel = 1, SIZE(scf_env%ot_preconditioner)
3793 128 : IF (ASSOCIATED(scf_env%ot_preconditioner(channel)%preconditioner)) THEN
3794 88 : CALL destroy_preconditioner(scf_env%ot_preconditioner(channel)%preconditioner)
3795 88 : DEALLOCATE (scf_env%ot_preconditioner(channel)%preconditioner)
3796 : END IF
3797 : END DO
3798 40 : DEALLOCATE (scf_env%ot_preconditioner)
3799 40 : NULLIFY (scf_env%ot_preconditioner)
3800 : END IF
3801 306 : ALLOCATE (scf_env%ot_preconditioner(nchannels))
3802 190 : DO channel = 1, nchannels
3803 132 : ALLOCATE (scf_env%ot_preconditioner(channel)%preconditioner)
3804 : CALL init_preconditioner(scf_env%ot_preconditioner(channel)%preconditioner, &
3805 190 : para_env, blacs_env)
3806 : END DO
3807 :
3808 58 : NULLIFY (eigenvalues, matrix_h_im, matrix_h_re, matrix_s_im, matrix_s_re, &
3809 58 : matrix_t_im, matrix_t_re, mo_coeff_im, mo_coeff_re)
3810 58 : IF (prec_type == ot_precond_full_kinetic) THEN
3811 28 : CPASSERT(ASSOCIATED(matrix_t_kp))
3812 : END IF
3813 :
3814 186 : DO local_kpoint = 1, kplocal
3815 128 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
3816 :
3817 : ! S(k) is spin independent. Build it with the same global-group
3818 : ! redistribution protocol as the Davidson operator itself, then make
3819 : ! a group-local block-cyclic DBCSR pair for the preconditioner makers.
3820 : CALL distribute_kpoint_operator_channel( &
3821 : matrix_s_kp, kpoints, local_kpoint, 1, kp_dist, xkp, cell_to_index, sab_nl, &
3822 128 : rmatrix, cmatrix, tmpmat, fmwork, fmlocal_re, fmlocal_im)
3823 128 : CALL dbcsr_init_p(matrix_s_re)
3824 128 : CALL dbcsr_init_p(matrix_s_im)
3825 128 : CALL copy_fm_to_dbcsr_bc(fmlocal_re, matrix_s_re)
3826 128 : CALL copy_fm_to_dbcsr_bc(fmlocal_im, matrix_s_im)
3827 :
3828 128 : IF (prec_type == ot_precond_full_kinetic) THEN
3829 : CALL distribute_kpoint_operator_channel( &
3830 : matrix_t_kp, kpoints, local_kpoint, 1, kp_dist, xkp, cell_to_index, sab_nl, &
3831 56 : rmatrix, cmatrix, tmpmat, fmwork, fmlocal_re, fmlocal_im)
3832 56 : CALL dbcsr_init_p(matrix_t_re)
3833 56 : CALL dbcsr_init_p(matrix_t_im)
3834 56 : CALL copy_fm_to_dbcsr_bc(fmlocal_re, matrix_t_re)
3835 56 : CALL copy_fm_to_dbcsr_bc(fmlocal_im, matrix_t_im)
3836 : END IF
3837 :
3838 260 : DO ispin = 1, nspin
3839 132 : channel = (local_kpoint - 1)*nspin + ispin
3840 260 : IF (prec_type == ot_precond_full_all .OR. prec_type == ot_precond_full_single_inverse) THEN
3841 : CALL distribute_kpoint_operator_channel( &
3842 : matrix_ks_kp, kpoints, local_kpoint, ispin, kp_dist, xkp, cell_to_index, sab_nl, &
3843 60 : rmatrix, cmatrix, tmpmat, fmwork, fmlocal_re, fmlocal_im)
3844 60 : CALL dbcsr_init_p(matrix_h_re)
3845 60 : CALL dbcsr_init_p(matrix_h_im)
3846 60 : CALL copy_fm_to_dbcsr_bc(fmlocal_re, matrix_h_re)
3847 60 : CALL copy_fm_to_dbcsr_bc(fmlocal_im, matrix_h_im)
3848 : ! Davidson iterates the full MO set, so its active space is all nmo
3849 : ! MOs of the channel: assemble the complex channel matrices on the
3850 : ! kpoint-group context and call the complex H-based cores directly
3851 : ! (mirrors the complex preconditioners of the K-point OT case; the
3852 : ! FULL_SINGLE_INVERSE variant trades the state selectivity of
3853 : ! FULL_ALL for one shifted factorization and a gemm application)
3854 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=mo_coeff_re, eigenvalues=eigenvalues, &
3855 60 : nmo=nmo)
3856 60 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=mo_coeff_im)
3857 60 : CALL cp_cfm_create(matrix_c, mo_coeff_re%matrix_struct, name="dav pc C")
3858 60 : CALL cp_cfm_scale_and_add_fm(z_zero, matrix_c, z_one, mo_coeff_re)
3859 60 : CALL cp_cfm_scale_and_add_fm(z_one, matrix_c, gaussi, mo_coeff_im)
3860 : CALL dbcsr_pair_to_cfm(matrix_h_re, matrix_h_im, &
3861 : scf_env%ot_preconditioner(channel)%preconditioner, &
3862 60 : 'dav pc H', matrix_h)
3863 : CALL dbcsr_pair_to_cfm(matrix_s_re, matrix_s_im, &
3864 : scf_env%ot_preconditioner(channel)%preconditioner, &
3865 60 : 'dav pc S', matrix_s)
3866 60 : IF (prec_type == ot_precond_full_all) THEN
3867 : CALL make_complex_full_all(scf_env%ot_preconditioner(channel)%preconditioner, &
3868 : matrix_c, matrix_h, matrix_s, eigenvalues(1:nmo), &
3869 44 : energy_gap)
3870 : ELSE
3871 : CALL make_complex_full_single_inverse( &
3872 : scf_env%ot_preconditioner(channel)%preconditioner, &
3873 16 : matrix_c, matrix_h, matrix_s, energy_gap)
3874 : END IF
3875 60 : CALL cp_cfm_release(matrix_c)
3876 60 : CALL cp_cfm_release(matrix_h)
3877 60 : CALL cp_cfm_release(matrix_s)
3878 60 : CALL dbcsr_release_p(matrix_h_re)
3879 60 : CALL dbcsr_release_p(matrix_h_im)
3880 72 : ELSE IF (prec_type == ot_precond_full_kinetic) THEN
3881 : CALL make_preconditioner_complex_full_kinetic( &
3882 : scf_env%ot_preconditioner(channel)%preconditioner, &
3883 56 : matrix_t_re, matrix_t_im, matrix_s_re, matrix_s_im, energy_gap, solver_type)
3884 16 : ELSE IF (prec_type == ot_precond_s_inverse) THEN
3885 : CALL make_preconditioner_complex_full_s_inverse( &
3886 : scf_env%ot_preconditioner(channel)%preconditioner, &
3887 16 : matrix_s_re, matrix_s_im, solver_type)
3888 : END IF
3889 : END DO
3890 :
3891 128 : CALL dbcsr_release_p(matrix_s_re)
3892 128 : CALL dbcsr_release_p(matrix_s_im)
3893 186 : IF (prec_type == ot_precond_full_kinetic) THEN
3894 56 : CALL dbcsr_release_p(matrix_t_re)
3895 56 : CALL dbcsr_release_p(matrix_t_im)
3896 : END IF
3897 : END DO
3898 :
3899 58 : CALL timestop(handle)
3900 :
3901 58 : END SUBROUTINE davidson_kp_prepare_preconditioners
3902 :
3903 : ! **************************************************************************************************
3904 : !> \brief Kpoint diagonalization routine
3905 : !> Transforms matrices to kpoint, distributes kpoint groups, performs diagonalization
3906 : !> \param matrix_s Overlap matrices (RS indices, global)
3907 : !> \param kpoints Kpoint environment
3908 : !> \param fmwork full matrices distributed over all groups
3909 : !> \par History
3910 : !> 02.2026 created [JGH]
3911 : ! **************************************************************************************************
3912 12 : SUBROUTINE diag_kp_smat(matrix_s, kpoints, fmwork)
3913 :
3914 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
3915 : TYPE(kpoint_type), POINTER :: kpoints
3916 : TYPE(cp_fm_type), DIMENSION(:) :: fmwork
3917 :
3918 : CHARACTER(len=*), PARAMETER :: routineN = 'diag_kp_smat'
3919 : COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp, 0.0_dp), &
3920 : czero = (0.0_dp, 0.0_dp)
3921 :
3922 12 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ceig
3923 : INTEGER :: handle, igroup, ik, ikp, indx, kplocal, &
3924 : nao, nkp, nkp_groups
3925 : INTEGER, DIMENSION(2) :: kp_range
3926 12 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
3927 12 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3928 : LOGICAL :: my_kpgrp, use_real_wfn
3929 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
3930 12 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
3931 12 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
3932 : TYPE(cp_cfm_type) :: csmat, cwork
3933 12 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools
3934 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
3935 : TYPE(cp_fm_type) :: fmdummy, fmlocal, rsmat
3936 : TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix, tmpmat
3937 : TYPE(kpoint_env_type), POINTER :: kp
3938 : TYPE(mp_para_env_type), POINTER :: para_env
3939 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3940 12 : POINTER :: sab_nl
3941 : TYPE(qs_matrix_pools_type), POINTER :: mpools
3942 :
3943 12 : CALL timeset(routineN, handle)
3944 :
3945 12 : NULLIFY (sab_nl)
3946 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
3947 : nkp_groups=nkp_groups, kp_dist=kp_dist, sab_nl=sab_nl, &
3948 12 : cell_to_index=cell_to_index)
3949 12 : CPASSERT(ASSOCIATED(sab_nl))
3950 12 : kplocal = kp_range(2) - kp_range(1) + 1
3951 :
3952 : ! allocate some work matrices
3953 12 : ALLOCATE (rmatrix, cmatrix, tmpmat)
3954 : CALL dbcsr_create(rmatrix, template=matrix_s(1, 1)%matrix, &
3955 12 : matrix_type=dbcsr_type_symmetric)
3956 : CALL dbcsr_create(cmatrix, template=matrix_s(1, 1)%matrix, &
3957 12 : matrix_type=dbcsr_type_antisymmetric)
3958 : CALL dbcsr_create(tmpmat, template=matrix_s(1, 1)%matrix, &
3959 12 : matrix_type=dbcsr_type_no_symmetry)
3960 12 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
3961 12 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
3962 :
3963 : ! fm pools to be used within a kpoint group
3964 12 : CALL get_kpoint_info(kpoints, mpools=mpools)
3965 12 : CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
3966 :
3967 12 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
3968 12 : CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
3969 :
3970 12 : IF (use_real_wfn) THEN
3971 0 : CALL cp_fm_create(rsmat, matrix_struct)
3972 : ELSE
3973 12 : CALL cp_cfm_create(csmat, matrix_struct)
3974 12 : CALL cp_cfm_create(cwork, matrix_struct)
3975 : END IF
3976 :
3977 12 : CALL cp_fm_get_info(fmwork(1), nrow_global=nao)
3978 60 : ALLOCATE (eigenvalues(nao), ceig(nao))
3979 :
3980 12 : para_env => kpoints%blacs_env_all%para_env
3981 644 : ALLOCATE (info(kplocal*nkp_groups, 2))
3982 :
3983 : ! Setup and start all the communication
3984 12 : indx = 0
3985 268 : DO ikp = 1, kplocal
3986 524 : DO igroup = 1, nkp_groups
3987 : ! number of current kpoint
3988 256 : ik = kp_dist(1, igroup) + ikp - 1
3989 256 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
3990 256 : indx = indx + 1
3991 256 : IF (use_real_wfn) THEN
3992 0 : CALL dbcsr_set(rmatrix, 0.0_dp)
3993 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_s, ispin=1, &
3994 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
3995 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
3996 0 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
3997 : ELSE
3998 256 : CALL dbcsr_set(rmatrix, 0.0_dp)
3999 256 : CALL dbcsr_set(cmatrix, 0.0_dp)
4000 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_s, ispin=1, &
4001 256 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
4002 256 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
4003 256 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
4004 256 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
4005 256 : CALL copy_dbcsr_to_fm(tmpmat, fmwork(2))
4006 : END IF
4007 : ! transfer to kpoint group
4008 : ! redistribution of matrices, new blacs environment
4009 : ! fmwork -> fmlocal -> rsmat/csmat
4010 512 : IF (use_real_wfn) THEN
4011 0 : IF (my_kpgrp) THEN
4012 0 : CALL cp_fm_start_copy_general(fmwork(1), rsmat, para_env, info(indx, 1))
4013 : ELSE
4014 0 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
4015 : END IF
4016 : ELSE
4017 256 : IF (my_kpgrp) THEN
4018 256 : CALL cp_fm_start_copy_general(fmwork(1), fmlocal, para_env, info(indx, 1))
4019 256 : CALL cp_fm_start_copy_general(fmwork(2), fmlocal, para_env, info(indx, 2))
4020 : ELSE
4021 0 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
4022 0 : CALL cp_fm_start_copy_general(fmwork(2), fmdummy, para_env, info(indx, 2))
4023 : END IF
4024 : END IF
4025 : END DO
4026 : END DO
4027 :
4028 : ! Finish communication then diagonalise in each group
4029 : indx = 0
4030 268 : DO ikp = 1, kplocal
4031 512 : DO igroup = 1, nkp_groups
4032 : ! number of current kpoint
4033 256 : ik = kp_dist(1, igroup) + ikp - 1
4034 256 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
4035 0 : indx = indx + 1
4036 256 : IF (my_kpgrp) THEN
4037 256 : IF (use_real_wfn) THEN
4038 0 : CALL cp_fm_finish_copy_general(rsmat, info(indx, 1))
4039 : ELSE
4040 256 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
4041 256 : CALL cp_cfm_scale_and_add_fm(z_zero, csmat, z_one, fmlocal)
4042 256 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
4043 256 : CALL cp_cfm_scale_and_add_fm(z_one, csmat, gaussi, fmlocal)
4044 : END IF
4045 : END IF
4046 : END DO
4047 :
4048 : ! Each kpoint group has now information on a kpoint to be diagonalized
4049 : ! Eigensolver Hermite or Symmetric
4050 256 : kp => kpoints%kp_env(ikp)%kpoint_env
4051 256 : IF (use_real_wfn) THEN
4052 0 : CALL choose_eigv_solver(rsmat, fmlocal, eigenvalues)
4053 : ELSE
4054 256 : CALL cp_cfm_heevd(csmat, cwork, eigenvalues)
4055 : END IF
4056 2680 : CPASSERT(ALL(eigenvalues(1:nao) >= 0.0_dp))
4057 268 : IF (use_real_wfn) THEN
4058 0 : CALL cp_fm_release(kp%shalf)
4059 0 : CALL cp_fm_create(kp%shalf, matrix_struct)
4060 0 : eigenvalues(1:nao) = SQRT(eigenvalues(1:nao))
4061 0 : CALL cp_fm_to_fm(fmlocal, rsmat)
4062 0 : CALL cp_fm_column_scale(rsmat, eigenvalues)
4063 : CALL parallel_gemm("N", "T", nao, nao, nao, 1.0_dp, rsmat, fmlocal, &
4064 0 : 0.0_dp, kp%shalf)
4065 : ELSE
4066 256 : CALL cp_cfm_release(kp%cshalf)
4067 256 : CALL cp_cfm_create(kp%cshalf, matrix_struct)
4068 2680 : ceig(1:nao) = SQRT(eigenvalues(1:nao))
4069 256 : CALL cp_cfm_to_cfm(cwork, csmat)
4070 256 : CALL cp_cfm_column_scale(csmat, ceig)
4071 : CALL parallel_gemm("N", "C", nao, nao, nao, cone, csmat, cwork, &
4072 256 : czero, kp%cshalf)
4073 : END IF
4074 : END DO
4075 :
4076 : ! Clean up communication
4077 : indx = 0
4078 268 : DO ikp = 1, kplocal
4079 524 : DO igroup = 1, nkp_groups
4080 : ! number of current kpoint
4081 256 : ik = kp_dist(1, igroup) + ikp - 1
4082 256 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
4083 256 : indx = indx + 1
4084 512 : IF (use_real_wfn) THEN
4085 0 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
4086 : ELSE
4087 256 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
4088 256 : CALL cp_fm_cleanup_copy_general(info(indx, 2))
4089 : END IF
4090 : END DO
4091 : END DO
4092 :
4093 : ! All done
4094 524 : DEALLOCATE (info)
4095 12 : DEALLOCATE (eigenvalues, ceig)
4096 :
4097 12 : CALL dbcsr_deallocate_matrix(rmatrix)
4098 12 : CALL dbcsr_deallocate_matrix(cmatrix)
4099 12 : CALL dbcsr_deallocate_matrix(tmpmat)
4100 :
4101 12 : IF (use_real_wfn) THEN
4102 0 : CALL cp_fm_release(rsmat)
4103 : ELSE
4104 12 : CALL cp_cfm_release(csmat)
4105 12 : CALL cp_cfm_release(cwork)
4106 : END IF
4107 12 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
4108 :
4109 12 : CALL timestop(handle)
4110 :
4111 36 : END SUBROUTINE diag_kp_smat
4112 :
4113 : END MODULE qs_scf_diagonalization
|
