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