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 : !> \brief Utility routines for qs_scf
9 : ! **************************************************************************************************
10 : MODULE qs_scf_loop_utils
11 : USE cp_control_types, ONLY: dft_control_type,&
12 : hairy_probes_type
13 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
14 : dbcsr_get_info,&
15 : dbcsr_p_type,&
16 : dbcsr_type
17 : USE cp_external_control, ONLY: external_control
18 : USE cp_log_handling, ONLY: cp_to_string
19 : USE input_section_types, ONLY: section_vals_type
20 : USE kinds, ONLY: default_string_length,&
21 : dp
22 : USE kpoint_types, ONLY: kpoint_type
23 : USE message_passing, ONLY: mp_para_env_type
24 : USE qs_density_matrices, ONLY: calculate_density_matrix
25 : USE qs_density_mixing_types, ONLY: broyden_mixing_nr,&
26 : direct_mixing_nr,&
27 : gspace_mixing_nr,&
28 : multisecant_mixing_nr,&
29 : no_mixing_nr,&
30 : pulay_mixing_nr
31 : USE qs_energy_types, ONLY: qs_energy_type
32 : USE qs_environment_types, ONLY: get_qs_env,&
33 : qs_environment_type
34 : USE qs_fb_env_methods, ONLY: fb_env_do_diag
35 : USE qs_gspace_mixing, ONLY: gspace_mixing
36 : USE qs_ks_types, ONLY: qs_ks_did_change,&
37 : qs_ks_env_type
38 : USE qs_mixing_utils, ONLY: self_consistency_check
39 : USE qs_mo_occupation, ONLY: set_mo_occupation
40 : USE qs_mo_types, ONLY: mo_set_type
41 : USE qs_mom_methods, ONLY: do_mom_diag
42 : USE qs_ot_scf, ONLY: ot_scf_destroy,&
43 : ot_scf_mini
44 : USE qs_outer_scf, ONLY: outer_loop_gradient
45 : USE qs_rho_methods, ONLY: qs_rho_update_rho
46 : USE qs_rho_types, ONLY: qs_rho_get,&
47 : qs_rho_type
48 : USE qs_scf_diagonalization, ONLY: do_block_davidson_diag,&
49 : do_block_krylov_diag,&
50 : do_general_diag,&
51 : do_general_diag_kp,&
52 : do_ot_diag,&
53 : do_roks_diag,&
54 : do_scf_diag_subspace,&
55 : do_special_diag
56 : USE qs_scf_methods, ONLY: scf_env_density_mixing
57 : USE qs_scf_output, ONLY: qs_scf_print_summary
58 : USE qs_scf_types, ONLY: &
59 : block_davidson_diag_method_nr, block_krylov_diag_method_nr, filter_matrix_diag_method_nr, &
60 : general_diag_method_nr, ot_diag_method_nr, ot_method_nr, qs_scf_env_type, &
61 : smeagol_method_nr, special_diag_method_nr
62 : USE scf_control_types, ONLY: scf_control_type,&
63 : smear_type
64 : USE smeagol_interface, ONLY: run_smeagol_emtrans
65 : USE tblite_interface, ONLY: tb_native_scc_mixer_active
66 : #include "./base/base_uses.f90"
67 :
68 : IMPLICIT NONE
69 :
70 : PRIVATE
71 :
72 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_loop_utils'
73 :
74 : PUBLIC :: qs_scf_set_loop_flags, &
75 : qs_scf_new_mos, qs_scf_new_mos_kp, &
76 : qs_scf_density_mixing, qs_scf_check_inner_exit, &
77 : qs_scf_check_outer_exit, qs_scf_inner_finalize, qs_scf_rho_update
78 :
79 : CONTAINS
80 :
81 : ! **************************************************************************************************
82 : !> \brief computes properties for a given hamiltonian using the current wfn
83 : !> \param scf_env ...
84 : !> \param diis_step ...
85 : !> \param energy_only ...
86 : !> \param just_energy ...
87 : !> \param exit_inner_loop ...
88 : ! **************************************************************************************************
89 21304 : SUBROUTINE qs_scf_set_loop_flags(scf_env, diis_step, &
90 : energy_only, just_energy, exit_inner_loop)
91 :
92 : TYPE(qs_scf_env_type), POINTER :: scf_env
93 : LOGICAL :: diis_step, energy_only, just_energy, &
94 : exit_inner_loop
95 :
96 : ! Some flags needed to be set at the beginning of the loop
97 :
98 21304 : diis_step = .FALSE.
99 21304 : energy_only = .FALSE.
100 21304 : just_energy = .FALSE.
101 :
102 : ! SCF loop, optimisation of the wfn coefficients
103 : ! qs_env%rho%rho_r and qs_env%rho%rho_g should be up to date here
104 :
105 21304 : scf_env%iter_count = 0
106 21304 : exit_inner_loop = .FALSE.
107 :
108 21304 : END SUBROUTINE qs_scf_set_loop_flags
109 :
110 : ! **************************************************************************************************
111 : !> \brief takes known energy and derivatives and produces new wfns
112 : !> and or density matrix
113 : !> \param qs_env ...
114 : !> \param scf_env ...
115 : !> \param scf_control ...
116 : !> \param scf_section ...
117 : !> \param diis_step ...
118 : !> \param energy_only ...
119 : !> \param probe ...
120 : ! **************************************************************************************************
121 176991 : SUBROUTINE qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, &
122 : energy_only, probe)
123 : TYPE(qs_environment_type), POINTER :: qs_env
124 : TYPE(qs_scf_env_type), POINTER :: scf_env
125 : TYPE(scf_control_type), POINTER :: scf_control
126 : TYPE(section_vals_type), POINTER :: scf_section
127 : LOGICAL :: diis_step, energy_only
128 : TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
129 : POINTER :: probe
130 :
131 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_new_mos'
132 :
133 : INTEGER :: handle, ispin
134 : LOGICAL :: disable_diis, has_unit_metric, &
135 : skip_diag_sub
136 : REAL(KIND=dp) :: saved_eps_diis
137 176991 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
138 : TYPE(dft_control_type), POINTER :: dft_control
139 176991 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
140 : TYPE(qs_energy_type), POINTER :: energy
141 : TYPE(qs_ks_env_type), POINTER :: ks_env
142 : TYPE(qs_rho_type), POINTER :: rho
143 :
144 176991 : CALL timeset(routineN, handle)
145 :
146 176991 : NULLIFY (energy, ks_env, matrix_ks, matrix_s, rho, mos, dft_control)
147 :
148 : CALL get_qs_env(qs_env=qs_env, &
149 : matrix_s=matrix_s, energy=energy, &
150 : ks_env=ks_env, &
151 : matrix_ks=matrix_ks, rho=rho, mos=mos, &
152 : dft_control=dft_control, &
153 176991 : has_unit_metric=has_unit_metric)
154 176991 : scf_env%iter_param = 0.0_dp
155 : disable_diis = dft_control%qs_control%xtb_control%do_tblite .AND. &
156 176991 : tb_native_scc_mixer_active(dft_control)
157 : IF (disable_diis) THEN
158 10008 : saved_eps_diis = scf_control%eps_diis
159 10008 : scf_control%eps_diis = 0.0_dp
160 : END IF
161 :
162 : ! transfer total_zeff_corr from qs_env to scf_env only if
163 : ! correct_el_density_dip is switched on [SGh]
164 176991 : IF (dft_control%correct_el_density_dip) THEN
165 40 : scf_env%sum_zeff_corr = qs_env%total_zeff_corr
166 40 : IF (ABS(qs_env%total_zeff_corr) > 0.0_dp) THEN
167 40 : IF (scf_env%method /= general_diag_method_nr) THEN
168 : CALL cp_abort(__LOCATION__, &
169 : "Please use ALGORITHM STANDARD in "// &
170 : "SCF%DIAGONALIZATION if "// &
171 : "CORE_CORRECTION /= 0.0 and "// &
172 0 : "SURFACE_DIPOLE_CORRECTION TRUE ")
173 40 : ELSEIF (dft_control%roks) THEN
174 : CALL cp_abort(__LOCATION__, &
175 : "Combination of "// &
176 : "CORE_CORRECTION /= 0.0 and "// &
177 : "SURFACE_DIPOLE_CORRECTION TRUE "// &
178 0 : "is not implemented with ROKS")
179 40 : ELSEIF (scf_control%diagonalization%mom) THEN
180 : CALL cp_abort(__LOCATION__, &
181 : "Combination of "// &
182 : "CORE_CORRECTION /= 0.0 and "// &
183 : "SURFACE_DIPOLE_CORRECTION TRUE "// &
184 0 : "is not implemented with SCF%MOM")
185 : END IF
186 : END IF
187 : END IF
188 :
189 176991 : SELECT CASE (scf_env%method)
190 : CASE DEFAULT
191 : CALL cp_abort(__LOCATION__, &
192 : "unknown scf method: "// &
193 0 : cp_to_string(scf_env%method))
194 :
195 : ! *************************************************************************
196 : ! Filter matrix diagonalisation: ugly implementation at this point of time
197 : ! *************************************************************************
198 : CASE (filter_matrix_diag_method_nr)
199 :
200 80 : IF (ABS(qs_env%total_zeff_corr) > 0.0_dp) THEN
201 : CALL cp_abort(__LOCATION__, &
202 : "CORE_CORRECTION /= 0.0 plus SURFACE_DIPOLE_CORRECTION TRUE "// &
203 0 : "requires SCF%DIAGONALIZATION: ALGORITHM STANDARD")
204 : END IF
205 : CALL fb_env_do_diag(scf_env%filter_matrix_env, qs_env, &
206 80 : matrix_ks, matrix_s, scf_section, diis_step)
207 :
208 : ! Diagonlization in non orthonormal case
209 : CASE (general_diag_method_nr)
210 82997 : IF (dft_control%roks) THEN
211 : CALL do_roks_diag(scf_env, mos, matrix_ks, matrix_s, &
212 : scf_control, scf_section, diis_step, &
213 592 : has_unit_metric)
214 : ELSE
215 82405 : IF (scf_control%diagonalization%mom) THEN
216 : CALL do_mom_diag(scf_env, mos, matrix_ks, &
217 : matrix_s, scf_control, scf_section, &
218 324 : diis_step)
219 : ELSE
220 82081 : IF (dft_control%hairy_probes .EQV. .TRUE.) THEN
221 : CALL do_general_diag(scf_env, mos, matrix_ks, &
222 : matrix_s, scf_control, scf_section, &
223 : diis_step, &
224 14 : probe)
225 : ELSE
226 : CALL do_general_diag(scf_env, mos, matrix_ks, &
227 : matrix_s, scf_control, scf_section, &
228 82067 : diis_step)
229 : END IF
230 : END IF
231 82405 : IF (scf_control%do_diag_sub) THEN
232 : skip_diag_sub = (scf_env%subspace_env%eps_diag_sub > 0.0_dp) .AND. &
233 10 : (scf_env%iter_count == 1 .OR. scf_env%iter_delta > scf_env%subspace_env%eps_diag_sub)
234 : IF (.NOT. skip_diag_sub) THEN
235 : CALL do_scf_diag_subspace(qs_env, scf_env, scf_env%subspace_env, mos, rho, &
236 10 : ks_env, scf_section, scf_control)
237 : END IF
238 : END IF
239 : END IF
240 : ! Diagonlization in orthonormal case
241 : CASE (special_diag_method_nr)
242 18410 : IF (dft_control%roks) THEN
243 : CALL do_roks_diag(scf_env, mos, matrix_ks, matrix_s, &
244 : scf_control, scf_section, diis_step, &
245 512 : has_unit_metric)
246 : ELSE
247 : CALL do_special_diag(scf_env, mos, matrix_ks, &
248 : scf_control, scf_section, &
249 17898 : diis_step)
250 : END IF
251 : ! OT diagonalization
252 : CASE (ot_diag_method_nr)
253 : CALL do_ot_diag(scf_env, mos, matrix_ks, matrix_s, &
254 64 : scf_control, scf_section, diis_step)
255 : ! Block Krylov diagonlization
256 : CASE (block_krylov_diag_method_nr)
257 40 : IF ((scf_env%krylov_space%eps_std_diag > 0.0_dp) .AND. &
258 : (scf_env%iter_count == 1 .OR. scf_env%iter_delta > scf_env%krylov_space%eps_std_diag)) THEN
259 2 : IF (scf_env%krylov_space%always_check_conv) THEN
260 : CALL do_block_krylov_diag(scf_env, mos, matrix_ks, &
261 0 : scf_control, scf_section, check_moconv_only=.TRUE.)
262 : END IF
263 : CALL do_general_diag(scf_env, mos, matrix_ks, &
264 2 : matrix_s, scf_control, scf_section, diis_step)
265 : ELSE
266 : CALL do_block_krylov_diag(scf_env, mos, matrix_ks, &
267 38 : scf_control, scf_section)
268 : END IF
269 40 : IF (scf_control%do_diag_sub) THEN
270 : skip_diag_sub = (scf_env%subspace_env%eps_diag_sub > 0.0_dp) .AND. &
271 0 : (scf_env%iter_count == 1 .OR. scf_env%iter_delta > scf_env%subspace_env%eps_diag_sub)
272 : IF (.NOT. skip_diag_sub) THEN
273 : CALL do_scf_diag_subspace(qs_env, scf_env, scf_env%subspace_env, mos, rho, &
274 0 : ks_env, scf_section, scf_control)
275 : END IF
276 : END IF
277 : ! Block Davidson diagonlization
278 : CASE (block_davidson_diag_method_nr)
279 : CALL do_block_davidson_diag(qs_env, scf_env, mos, matrix_ks, matrix_s, scf_control, &
280 84 : scf_section, .FALSE.)
281 : ! OT without diagonlization. Needs special treatment for SCP runs
282 : CASE (ot_method_nr)
283 : CALL qs_scf_loop_do_ot(qs_env, scf_env, scf_control%smear, mos, rho, &
284 : qs_env%mo_derivs, energy%total, &
285 176991 : matrix_s, energy_only=energy_only, has_unit_metric=has_unit_metric)
286 : END SELECT
287 176991 : IF (disable_diis) scf_control%eps_diis = saved_eps_diis
288 :
289 176991 : energy%kTS = 0.0_dp
290 176991 : energy%efermi = 0.0_dp
291 176991 : CALL get_qs_env(qs_env, mos=mos)
292 378119 : DO ispin = 1, SIZE(mos)
293 201128 : energy%kTS = energy%kTS + mos(ispin)%kTS
294 378119 : energy%efermi = energy%efermi + mos(ispin)%mu
295 : END DO
296 176991 : energy%efermi = energy%efermi/REAL(SIZE(mos), KIND=dp)
297 :
298 176991 : CALL timestop(handle)
299 :
300 176991 : END SUBROUTINE qs_scf_new_mos
301 :
302 : ! **************************************************************************************************
303 : !> \brief Updates MOs and density matrix using diagonalization
304 : !> Kpoint code
305 : !> \param qs_env ...
306 : !> \param scf_env ...
307 : !> \param scf_control ...
308 : !> \param diis_step ...
309 : !> \param probe ...
310 : ! **************************************************************************************************
311 12396 : SUBROUTINE qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step, probe)
312 : TYPE(qs_environment_type), POINTER :: qs_env
313 : TYPE(qs_scf_env_type), POINTER :: scf_env
314 : TYPE(scf_control_type), POINTER :: scf_control
315 : LOGICAL :: diis_step
316 : TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
317 : POINTER :: probe
318 :
319 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_new_mos_kp'
320 :
321 : INTEGER :: handle, ispin
322 : LOGICAL :: disable_diis, has_unit_metric
323 : REAL(dp) :: diis_error, saved_eps_diis
324 12396 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
325 : TYPE(dft_control_type), POINTER :: dft_control
326 : TYPE(kpoint_type), POINTER :: kpoints
327 12396 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos
328 : TYPE(qs_energy_type), POINTER :: energy
329 :
330 12396 : CALL timeset(routineN, handle)
331 :
332 12396 : NULLIFY (dft_control, kpoints, matrix_ks, matrix_s)
333 :
334 12396 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, kpoints=kpoints)
335 12396 : scf_env%iter_param = 0.0_dp
336 : disable_diis = dft_control%qs_control%xtb_control%do_tblite .AND. &
337 12396 : tb_native_scc_mixer_active(dft_control)
338 : IF (disable_diis) THEN
339 3736 : saved_eps_diis = scf_control%eps_diis
340 3736 : scf_control%eps_diis = 0.0_dp
341 : END IF
342 :
343 12396 : IF (dft_control%roks) &
344 0 : CPABORT("KP code: ROKS method not available: ")
345 :
346 12396 : SELECT CASE (scf_env%method)
347 : CASE DEFAULT
348 : CALL cp_abort(__LOCATION__, &
349 : "KP code: Unknown scf method: "// &
350 0 : cp_to_string(scf_env%method))
351 : CASE (general_diag_method_nr)
352 : ! Diagonlization in non orthonormal case
353 12396 : CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s)
354 12396 : IF (dft_control%hairy_probes .EQV. .TRUE.) THEN
355 0 : scf_control%smear%do_smear = .FALSE.
356 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, .TRUE., &
357 0 : diis_step, diis_error, qs_env, probe)
358 : ELSE
359 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, .TRUE., &
360 12396 : diis_step, diis_error, qs_env)
361 : END IF
362 12396 : IF (diis_step) THEN
363 3740 : scf_env%iter_param = diis_error
364 3740 : scf_env%iter_method = "DIIS/Diag."
365 : ELSE
366 8656 : IF (scf_env%mixing_method == 0) THEN
367 0 : scf_env%iter_method = "NoMix/Diag."
368 8656 : ELSE IF (scf_env%mixing_method == 1) THEN
369 7338 : scf_env%iter_param = scf_env%p_mix_alpha
370 7338 : scf_env%iter_method = "P_Mix/Diag."
371 1318 : ELSEIF (scf_env%mixing_method > 1) THEN
372 1318 : scf_env%iter_param = scf_env%mixing_store%alpha
373 1318 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Diag."
374 : END IF
375 : END IF
376 : CASE (special_diag_method_nr)
377 0 : CALL get_qs_env(qs_env=qs_env, has_unit_metric=has_unit_metric)
378 0 : CPASSERT(has_unit_metric)
379 : ! Diagonlization in orthonormal case
380 : CALL cp_abort(__LOCATION__, &
381 : "KP code: Scf method not available: "// &
382 0 : cp_to_string(scf_env%method))
383 : CASE (ot_diag_method_nr, &
384 : block_krylov_diag_method_nr, &
385 : block_davidson_diag_method_nr, &
386 : ot_method_nr)
387 : CALL cp_abort(__LOCATION__, &
388 : "KP code: Scf method not available: "// &
389 0 : cp_to_string(scf_env%method))
390 : CASE (smeagol_method_nr)
391 : ! SMEAGOL interface
392 0 : diis_step = .FALSE.
393 0 : IF (scf_env%mixing_method == 0) THEN
394 0 : scf_env%iter_method = "NoMix/SMGL"
395 0 : ELSE IF (scf_env%mixing_method == 1) THEN
396 0 : scf_env%iter_param = scf_env%p_mix_alpha
397 0 : scf_env%iter_method = "P_Mix/SMGL"
398 0 : ELSE IF (scf_env%mixing_method > 1) THEN
399 0 : scf_env%iter_param = scf_env%mixing_store%alpha
400 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/SMGL"
401 : END IF
402 12396 : CALL run_smeagol_emtrans(qs_env, last=.FALSE., iter=scf_env%iter_count, rho_ao_kp=scf_env%p_mix_new)
403 : END SELECT
404 12396 : IF (disable_diis) scf_control%eps_diis = saved_eps_diis
405 :
406 12396 : CALL get_qs_env(qs_env=qs_env, energy=energy)
407 12396 : energy%kTS = 0.0_dp
408 12396 : energy%efermi = 0.0_dp
409 12396 : mos => kpoints%kp_env(1)%kpoint_env%mos
410 25454 : DO ispin = 1, SIZE(mos, 2)
411 13058 : energy%kTS = energy%kTS + mos(1, ispin)%kTS
412 25454 : energy%efermi = energy%efermi + mos(1, ispin)%mu
413 : END DO
414 12396 : energy%efermi = energy%efermi/REAL(SIZE(mos, 2), KIND=dp)
415 :
416 12396 : CALL timestop(handle)
417 :
418 12396 : END SUBROUTINE qs_scf_new_mos_kp
419 :
420 : ! **************************************************************************************************
421 : !> \brief the inner loop of scf, specific to using to the orbital transformation method
422 : !> basically, in goes the ks matrix out goes a new p matrix
423 : !> \param qs_env ...
424 : !> \param scf_env ...
425 : !> \param smear ...
426 : !> \param mos ...
427 : !> \param rho ...
428 : !> \param mo_derivs ...
429 : !> \param total_energy ...
430 : !> \param matrix_s ...
431 : !> \param energy_only ...
432 : !> \param has_unit_metric ...
433 : !> \par History
434 : !> 03.2006 created [Joost VandeVondele]
435 : !> 2013 moved from qs_scf [Florian Schiffmann]
436 : ! **************************************************************************************************
437 75316 : SUBROUTINE qs_scf_loop_do_ot(qs_env, scf_env, smear, mos, rho, mo_derivs, total_energy, &
438 : matrix_s, energy_only, has_unit_metric)
439 :
440 : TYPE(qs_environment_type), POINTER :: qs_env
441 : TYPE(qs_scf_env_type), POINTER :: scf_env
442 : TYPE(smear_type), POINTER :: smear
443 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
444 : TYPE(qs_rho_type), POINTER :: rho
445 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs
446 : REAL(KIND=dp), INTENT(IN) :: total_energy
447 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
448 : LOGICAL, INTENT(INOUT) :: energy_only
449 : LOGICAL, INTENT(IN) :: has_unit_metric
450 :
451 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_loop_do_ot'
452 :
453 : INTEGER :: handle, ispin
454 75316 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
455 : TYPE(dbcsr_type), POINTER :: orthogonality_metric
456 :
457 75316 : CALL timeset(routineN, handle)
458 75316 : NULLIFY (rho_ao)
459 :
460 75316 : CALL qs_rho_get(rho, rho_ao=rho_ao)
461 :
462 75316 : IF (has_unit_metric) THEN
463 18290 : NULLIFY (orthogonality_metric)
464 : ELSE
465 57026 : orthogonality_metric => matrix_s(1)%matrix
466 : END IF
467 :
468 : ! in case of LSD the first spin qs_ot_env will drive the minimization
469 : ! in the case of a restricted calculation, it will make sure the spin orbitals are equal
470 :
471 : CALL ot_scf_mini(mos, mo_derivs, smear, orthogonality_metric, &
472 : total_energy, energy_only, scf_env%iter_delta, &
473 75316 : scf_env%qs_ot_env)
474 :
475 163361 : DO ispin = 1, SIZE(mos)
476 163361 : CALL set_mo_occupation(mo_set=mos(ispin), smear=smear)
477 : END DO
478 :
479 163361 : DO ispin = 1, SIZE(mos)
480 : CALL calculate_density_matrix(mos(ispin), &
481 : rho_ao(ispin)%matrix, &
482 163361 : use_dbcsr=.TRUE.)
483 : END DO
484 :
485 75316 : scf_env%iter_method = scf_env%qs_ot_env(1)%OT_METHOD_FULL
486 75316 : scf_env%iter_param = scf_env%qs_ot_env(1)%ds_min
487 75316 : qs_env%broyden_adaptive_sigma = scf_env%qs_ot_env(1)%broyden_adaptive_sigma
488 :
489 75316 : CALL timestop(handle)
490 :
491 75316 : END SUBROUTINE qs_scf_loop_do_ot
492 :
493 : ! **************************************************************************************************
494 : !> \brief Performs the requested density mixing if any needed
495 : !> \param scf_env Holds SCF environment information
496 : !> \param rho All data for the electron density
497 : !> \param para_env Parallel environment
498 : !> \param diis_step Did we do a DIIS step?
499 : ! **************************************************************************************************
500 187113 : SUBROUTINE qs_scf_density_mixing(scf_env, rho, para_env, diis_step)
501 : TYPE(qs_scf_env_type), POINTER :: scf_env
502 : TYPE(qs_rho_type), POINTER :: rho
503 : TYPE(mp_para_env_type), POINTER :: para_env
504 : LOGICAL :: diis_step
505 :
506 187113 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
507 :
508 187113 : NULLIFY (rho_ao_kp)
509 :
510 187113 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
511 :
512 295670 : SELECT CASE (scf_env%mixing_method)
513 : CASE (direct_mixing_nr)
514 : CALL scf_env_density_mixing(scf_env%p_mix_new, &
515 : scf_env%mixing_store, rho_ao_kp, para_env, scf_env%iter_delta, scf_env%iter_count, &
516 108557 : diis=diis_step)
517 : CASE (gspace_mixing_nr, pulay_mixing_nr, broyden_mixing_nr, &
518 : multisecant_mixing_nr)
519 : ! Compute the difference p_out-p_in
520 : CALL self_consistency_check(rho_ao_kp, scf_env%p_delta, para_env, scf_env%p_mix_new, &
521 3240 : delta=scf_env%iter_delta)
522 : CASE (no_mixing_nr)
523 : CASE DEFAULT
524 : CALL cp_abort(__LOCATION__, &
525 : "unknown scf mixing method: "// &
526 187113 : cp_to_string(scf_env%mixing_method))
527 : END SELECT
528 :
529 187113 : END SUBROUTINE qs_scf_density_mixing
530 :
531 : ! **************************************************************************************************
532 : !> \brief checks whether exit conditions for outer loop are satisfied
533 : !> \param qs_env ...
534 : !> \param scf_env ...
535 : !> \param scf_control ...
536 : !> \param should_stop ...
537 : !> \param outer_loop_converged ...
538 : !> \param exit_outer_loop ...
539 : ! **************************************************************************************************
540 21815 : SUBROUTINE qs_scf_check_outer_exit(qs_env, scf_env, scf_control, should_stop, &
541 : outer_loop_converged, exit_outer_loop)
542 : TYPE(qs_environment_type), POINTER :: qs_env
543 : TYPE(qs_scf_env_type), POINTER :: scf_env
544 : TYPE(scf_control_type), POINTER :: scf_control
545 : LOGICAL :: should_stop, outer_loop_converged, &
546 : exit_outer_loop
547 :
548 : REAL(KIND=dp) :: outer_loop_eps
549 :
550 21815 : outer_loop_converged = .TRUE.
551 21815 : IF (scf_control%outer_scf%have_scf) THEN
552 : ! We have an outer SCF loop...
553 5811 : scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count + 1
554 5811 : outer_loop_converged = .FALSE.
555 :
556 5811 : CALL outer_loop_gradient(qs_env, scf_env)
557 : ! Multiple constraints: get largest deviation
558 17519 : outer_loop_eps = SQRT(MAXVAL(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count)**2))
559 :
560 5811 : IF (outer_loop_eps < scf_control%outer_scf%eps_scf) outer_loop_converged = .TRUE.
561 : END IF
562 :
563 : exit_outer_loop = should_stop .OR. outer_loop_converged .OR. &
564 21815 : scf_env%outer_scf%iter_count > scf_control%outer_scf%max_scf
565 :
566 21815 : END SUBROUTINE qs_scf_check_outer_exit
567 :
568 : ! **************************************************************************************************
569 : !> \brief checks whether exit conditions for inner loop are satisfied
570 : !> \param qs_env ...
571 : !> \param scf_env ...
572 : !> \param scf_control ...
573 : !> \param should_stop ...
574 : !> \param just_energy ...
575 : !> \param exit_inner_loop ...
576 : !> \param inner_loop_converged ...
577 : !> \param output_unit ...
578 : ! **************************************************************************************************
579 187113 : SUBROUTINE qs_scf_check_inner_exit(qs_env, scf_env, scf_control, should_stop, just_energy, &
580 : exit_inner_loop, inner_loop_converged, output_unit)
581 : TYPE(qs_environment_type), POINTER :: qs_env
582 : TYPE(qs_scf_env_type), POINTER :: scf_env
583 : TYPE(scf_control_type), POINTER :: scf_control
584 : LOGICAL :: should_stop, just_energy, &
585 : exit_inner_loop, inner_loop_converged
586 : INTEGER :: output_unit
587 :
588 187113 : inner_loop_converged = .FALSE.
589 187113 : exit_inner_loop = .FALSE.
590 :
591 : CALL external_control(should_stop, "SCF", target_time=qs_env%target_time, &
592 187113 : start_time=qs_env%start_time)
593 187113 : IF (scf_env%iter_delta < scf_control%eps_scf) THEN
594 18331 : IF (output_unit > 0) THEN
595 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I5,A/)") &
596 9348 : "*** SCF run converged in ", scf_env%iter_count, " steps ***"
597 : END IF
598 18331 : inner_loop_converged = .TRUE.
599 18331 : exit_inner_loop = .TRUE.
600 168782 : ELSE IF (should_stop .OR. scf_env%iter_count >= scf_control%max_scf) THEN
601 3894 : inner_loop_converged = .FALSE.
602 3894 : IF (just_energy) THEN
603 922 : exit_inner_loop = .FALSE.
604 : ELSE
605 2972 : exit_inner_loop = .TRUE.
606 2972 : IF (output_unit > 0) THEN
607 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I5,A/)") &
608 1492 : "Leaving inner SCF loop after reaching ", scf_env%iter_count, " steps."
609 : END IF
610 : END IF
611 : END IF
612 :
613 187113 : END SUBROUTINE qs_scf_check_inner_exit
614 :
615 : ! **************************************************************************************************
616 : !> \brief undoing density mixing. Important upon convergence
617 : !> \param scf_env ...
618 : !> \param rho ...
619 : !> \param dft_control ...
620 : !> \param para_env ...
621 : !> \param diis_step ...
622 : ! **************************************************************************************************
623 21303 : SUBROUTINE qs_scf_undo_mixing(scf_env, rho, dft_control, para_env, diis_step)
624 : TYPE(qs_scf_env_type), POINTER :: scf_env
625 : TYPE(qs_rho_type), POINTER :: rho
626 : TYPE(dft_control_type), POINTER :: dft_control
627 : TYPE(mp_para_env_type), POINTER :: para_env
628 : LOGICAL :: diis_step
629 :
630 : CHARACTER(len=default_string_length) :: name
631 : INTEGER :: ic, ispin, nc
632 21303 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
633 :
634 21303 : NULLIFY (rho_ao_kp)
635 :
636 21303 : IF (scf_env%mixing_method > 0) THEN
637 14126 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
638 14126 : nc = SIZE(scf_env%p_mix_new, 2)
639 27878 : SELECT CASE (scf_env%mixing_method)
640 : CASE (direct_mixing_nr)
641 : CALL scf_env_density_mixing(scf_env%p_mix_new, scf_env%mixing_store, &
642 : rho_ao_kp, para_env, scf_env%iter_delta, &
643 : scf_env%iter_count, diis=diis_step, &
644 13752 : invert=.TRUE.)
645 104086 : DO ic = 1, nc
646 205464 : DO ispin = 1, dft_control%nspins
647 101378 : CALL dbcsr_get_info(rho_ao_kp(ispin, ic)%matrix, name=name) ! keep the name
648 191712 : CALL dbcsr_copy(rho_ao_kp(ispin, ic)%matrix, scf_env%p_mix_new(ispin, ic)%matrix, name=name)
649 : END DO
650 : END DO
651 : CASE (gspace_mixing_nr, pulay_mixing_nr, broyden_mixing_nr, &
652 : multisecant_mixing_nr)
653 42298 : DO ic = 1, nc
654 56580 : DO ispin = 1, dft_control%nspins
655 28408 : CALL dbcsr_get_info(rho_ao_kp(ispin, ic)%matrix, name=name) ! keep the name
656 56206 : CALL dbcsr_copy(rho_ao_kp(ispin, ic)%matrix, scf_env%p_mix_new(ispin, ic)%matrix, name=name)
657 : END DO
658 : END DO
659 : END SELECT
660 : END IF
661 21303 : END SUBROUTINE qs_scf_undo_mixing
662 :
663 : ! **************************************************************************************************
664 : !> \brief Performs the updates rho (takes care of mixing as well)
665 : !> \param rho ...
666 : !> \param qs_env ...
667 : !> \param scf_env ...
668 : !> \param ks_env ...
669 : !> \param mix_rho ...
670 : ! **************************************************************************************************
671 187113 : SUBROUTINE qs_scf_rho_update(rho, qs_env, scf_env, ks_env, mix_rho)
672 : TYPE(qs_rho_type), POINTER :: rho
673 : TYPE(qs_environment_type), POINTER :: qs_env
674 : TYPE(qs_scf_env_type), POINTER :: scf_env
675 : TYPE(qs_ks_env_type), POINTER :: ks_env
676 : LOGICAL, INTENT(IN) :: mix_rho
677 :
678 : TYPE(mp_para_env_type), POINTER :: para_env
679 :
680 187113 : NULLIFY (para_env)
681 187113 : CALL get_qs_env(qs_env, para_env=para_env)
682 : ! ** update qs_env%rho
683 187113 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
684 : ! ** Density mixing through density matrix or on the reciprocal space grid (exclusive)
685 187113 : IF (mix_rho) THEN
686 : CALL gspace_mixing(qs_env, scf_env%mixing_method, scf_env%mixing_store, rho, &
687 2866 : para_env, scf_env%iter_count)
688 :
689 : END IF
690 187113 : CALL qs_ks_did_change(ks_env, rho_changed=.TRUE.)
691 :
692 187113 : END SUBROUTINE qs_scf_rho_update
693 :
694 : ! **************************************************************************************************
695 : !> \brief Performs the necessary steps before leaving innner scf loop
696 : !> \param scf_env ...
697 : !> \param qs_env ...
698 : !> \param diis_step ...
699 : !> \param output_unit ...
700 : ! **************************************************************************************************
701 21303 : SUBROUTINE qs_scf_inner_finalize(scf_env, qs_env, diis_step, output_unit)
702 : TYPE(qs_scf_env_type), POINTER :: scf_env
703 : TYPE(qs_environment_type), POINTER :: qs_env
704 : LOGICAL :: diis_step
705 : INTEGER, INTENT(IN) :: output_unit
706 :
707 : LOGICAL :: do_kpoints
708 : TYPE(dft_control_type), POINTER :: dft_control
709 : TYPE(mp_para_env_type), POINTER :: para_env
710 : TYPE(qs_energy_type), POINTER :: energy
711 : TYPE(qs_ks_env_type), POINTER :: ks_env
712 : TYPE(qs_rho_type), POINTER :: rho
713 :
714 21303 : NULLIFY (energy, rho, dft_control, ks_env)
715 :
716 : CALL get_qs_env(qs_env=qs_env, energy=energy, ks_env=ks_env, &
717 : rho=rho, dft_control=dft_control, para_env=para_env, &
718 21303 : do_kpoints=do_kpoints)
719 :
720 21303 : CALL cleanup_scf_loop(scf_env)
721 :
722 : ! now, print out energies and charges corresponding to the obtained wfn
723 : ! (this actually is not 100% consistent at this point)!
724 21303 : CALL qs_scf_print_summary(output_unit, qs_env)
725 :
726 21303 : CALL qs_scf_undo_mixing(scf_env, rho, dft_control, para_env, diis_step)
727 :
728 : ! *** update rspace rho since the mo changed
729 : ! *** this might not always be needed (i.e. no post calculation / no forces )
730 : ! *** but guarantees that rho and wfn are consistent at this point
731 21303 : CALL qs_scf_rho_update(rho, qs_env, scf_env, ks_env, mix_rho=.FALSE.)
732 :
733 21303 : END SUBROUTINE qs_scf_inner_finalize
734 :
735 : ! **************************************************************************************************
736 : !> \brief perform cleanup operations at the end of an scf loop
737 : !> \param scf_env ...
738 : !> \par History
739 : !> 03.2006 created [Joost VandeVondele]
740 : ! **************************************************************************************************
741 21303 : SUBROUTINE cleanup_scf_loop(scf_env)
742 : TYPE(qs_scf_env_type), INTENT(INOUT) :: scf_env
743 :
744 : CHARACTER(len=*), PARAMETER :: routineN = 'cleanup_scf_loop'
745 :
746 : INTEGER :: handle, ispin
747 :
748 21303 : CALL timeset(routineN, handle)
749 :
750 28480 : SELECT CASE (scf_env%method)
751 : CASE (ot_method_nr)
752 15609 : DO ispin = 1, SIZE(scf_env%qs_ot_env)
753 15609 : CALL ot_scf_destroy(scf_env%qs_ot_env(ispin))
754 : END DO
755 7177 : DEALLOCATE (scf_env%qs_ot_env)
756 : CASE (ot_diag_method_nr)
757 : !
758 : CASE (general_diag_method_nr)
759 : !
760 : CASE (special_diag_method_nr)
761 : !
762 : CASE (block_krylov_diag_method_nr, block_davidson_diag_method_nr)
763 : !
764 : CASE (filter_matrix_diag_method_nr)
765 : !
766 : CASE (smeagol_method_nr)
767 : !
768 : CASE DEFAULT
769 : CALL cp_abort(__LOCATION__, &
770 : "unknown scf method method:"// &
771 21303 : cp_to_string(scf_env%method))
772 : END SELECT
773 :
774 21303 : CALL timestop(handle)
775 :
776 21303 : END SUBROUTINE cleanup_scf_loop
777 :
778 : END MODULE qs_scf_loop_utils
|
