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