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