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 24459 : 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 24459 : diis_step = .FALSE.
101 24459 : energy_only = .FALSE.
102 24459 : 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 24459 : scf_env%iter_count = 0
108 24459 : exit_inner_loop = .FALSE.
109 :
110 24459 : 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 192939 : 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 192939 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
140 : TYPE(dft_control_type), POINTER :: dft_control
141 192939 : 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 192939 : CALL timeset(routineN, handle)
147 :
148 192939 : 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 192939 : has_unit_metric=has_unit_metric)
156 192939 : scf_env%iter_param = 0.0_dp
157 : disable_diis = dft_control%qs_control%xtb_control%do_tblite .AND. &
158 192939 : tb_native_scc_mixer_active(dft_control)
159 : IF (disable_diis) THEN
160 13012 : saved_eps_diis = scf_control%eps_diis
161 13012 : 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 192939 : 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 : ELSE IF (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 : ELSE IF (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 192939 : 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 96415 : 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 95805 : 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 95481 : 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 95467 : diis_step)
231 : END IF
232 : END IF
233 95805 : 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 18398 : 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 522 : has_unit_metric)
248 : ELSE
249 : CALL do_special_diag(scf_env, mos, matrix_ks, &
250 : scf_control, scf_section, &
251 17876 : 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 192939 : matrix_s, energy_only=energy_only, has_unit_metric=has_unit_metric)
288 : END SELECT
289 192939 : IF (disable_diis) scf_control%eps_diis = saved_eps_diis
290 :
291 192939 : energy%kTS = 0.0_dp
292 192939 : energy%efermi = 0.0_dp
293 192939 : CALL get_qs_env(qs_env, mos=mos)
294 417345 : DO ispin = 1, SIZE(mos)
295 224406 : energy%kTS = energy%kTS + mos(ispin)%kTS
296 417345 : energy%efermi = energy%efermi + mos(ispin)%mu
297 : END DO
298 192939 : energy%efermi = energy%efermi/REAL(SIZE(mos), KIND=dp)
299 :
300 192939 : CALL timestop(handle)
301 :
302 192939 : 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 33934 : 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 33934 : 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 33934 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos
330 : TYPE(qs_energy_type), POINTER :: energy
331 :
332 33934 : CALL timeset(routineN, handle)
333 :
334 33934 : NULLIFY (dft_control, kpoints, matrix_ks, matrix_s)
335 :
336 33934 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, kpoints=kpoints)
337 33934 : scf_env%iter_param = 0.0_dp
338 : disable_diis = dft_control%qs_control%xtb_control%do_tblite .AND. &
339 33934 : tb_native_scc_mixer_active(dft_control)
340 : IF (disable_diis) THEN
341 10638 : saved_eps_diis = scf_control%eps_diis
342 10638 : scf_control%eps_diis = 0.0_dp
343 : END IF
344 :
345 33934 : IF (dft_control%roks) THEN
346 0 : CPABORT("KP code: ROKS method not available: ")
347 : END IF
348 :
349 33934 : SELECT CASE (scf_env%method)
350 : CASE DEFAULT
351 : CALL cp_abort(__LOCATION__, &
352 : "KP code: Unknown scf method: "// &
353 0 : cp_to_string(scf_env%method))
354 : CASE (general_diag_method_nr)
355 : ! Diagonlization in non orthonormal case
356 33934 : CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s)
357 33934 : IF (dft_control%hairy_probes .EQV. .TRUE.) THEN
358 0 : scf_control%smear%do_smear = .FALSE.
359 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, .TRUE., &
360 0 : diis_step, diis_error, qs_env, probe)
361 : ELSE
362 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, .TRUE., &
363 33934 : diis_step, diis_error, qs_env)
364 : END IF
365 33934 : IF (diis_step) THEN
366 11130 : scf_env%iter_param = diis_error
367 11130 : scf_env%iter_method = "DIIS/Diag."
368 : ELSE
369 22804 : IF (scf_env%mixing_method == 0) THEN
370 0 : scf_env%iter_method = "NoMix/Diag."
371 22804 : ELSE IF (scf_env%mixing_method == 1) THEN
372 20294 : scf_env%iter_param = scf_env%p_mix_alpha
373 20294 : scf_env%iter_method = "P_Mix/Diag."
374 2510 : ELSE IF (scf_env%mixing_method > 1) THEN
375 2510 : scf_env%iter_param = scf_env%mixing_store%alpha
376 2510 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Diag."
377 : END IF
378 : END IF
379 : CASE (special_diag_method_nr)
380 0 : CALL get_qs_env(qs_env=qs_env, has_unit_metric=has_unit_metric)
381 0 : CPASSERT(has_unit_metric)
382 : ! Diagonlization in orthonormal case
383 : CALL cp_abort(__LOCATION__, &
384 : "KP code: Scf method not available: "// &
385 0 : cp_to_string(scf_env%method))
386 : CASE (ot_diag_method_nr, &
387 : block_krylov_diag_method_nr, &
388 : block_davidson_diag_method_nr, &
389 : ot_method_nr)
390 : CALL cp_abort(__LOCATION__, &
391 : "KP code: Scf method not available: "// &
392 0 : cp_to_string(scf_env%method))
393 : CASE (smeagol_method_nr)
394 : ! SMEAGOL interface
395 0 : diis_step = .FALSE.
396 0 : IF (scf_env%mixing_method == 0) THEN
397 0 : scf_env%iter_method = "NoMix/SMGL"
398 0 : ELSE IF (scf_env%mixing_method == 1) THEN
399 0 : scf_env%iter_param = scf_env%p_mix_alpha
400 0 : scf_env%iter_method = "P_Mix/SMGL"
401 0 : ELSE IF (scf_env%mixing_method > 1) THEN
402 0 : scf_env%iter_param = scf_env%mixing_store%alpha
403 0 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/SMGL"
404 : END IF
405 33934 : CALL run_smeagol_emtrans(qs_env, last=.FALSE., iter=scf_env%iter_count, rho_ao_kp=scf_env%p_mix_new)
406 : END SELECT
407 33934 : IF (disable_diis) scf_control%eps_diis = saved_eps_diis
408 :
409 33934 : CALL get_qs_env(qs_env=qs_env, energy=energy)
410 33934 : energy%kTS = 0.0_dp
411 33934 : energy%efermi = 0.0_dp
412 33934 : mos => kpoints%kp_env(1)%kpoint_env%mos
413 69678 : DO ispin = 1, SIZE(mos, 2)
414 35744 : energy%kTS = energy%kTS + mos(1, ispin)%kTS
415 69678 : energy%efermi = energy%efermi + mos(1, ispin)%mu
416 : END DO
417 33934 : energy%efermi = energy%efermi/REAL(SIZE(mos, 2), KIND=dp)
418 :
419 33934 : CALL timestop(handle)
420 :
421 33934 : END SUBROUTINE qs_scf_new_mos_kp
422 :
423 : ! **************************************************************************************************
424 : !> \brief the inner loop of scf, specific to using to the orbital transformation method
425 : !> basically, in goes the ks matrix out goes a new p matrix
426 : !> \param qs_env ...
427 : !> \param scf_env ...
428 : !> \param smear ...
429 : !> \param mos ...
430 : !> \param rho ...
431 : !> \param mo_derivs ...
432 : !> \param total_energy ...
433 : !> \param matrix_s ...
434 : !> \param energy_only ...
435 : !> \param has_unit_metric ...
436 : !> \par History
437 : !> 03.2006 created [Joost VandeVondele]
438 : !> 2013 moved from qs_scf [Florian Schiffmann]
439 : ! **************************************************************************************************
440 77848 : SUBROUTINE qs_scf_loop_do_ot(qs_env, scf_env, smear, mos, rho, mo_derivs, total_energy, &
441 : matrix_s, energy_only, has_unit_metric)
442 :
443 : TYPE(qs_environment_type), POINTER :: qs_env
444 : TYPE(qs_scf_env_type), POINTER :: scf_env
445 : TYPE(smear_type), POINTER :: smear
446 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
447 : TYPE(qs_rho_type), POINTER :: rho
448 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs
449 : REAL(KIND=dp), INTENT(IN) :: total_energy
450 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
451 : LOGICAL, INTENT(INOUT) :: energy_only
452 : LOGICAL, INTENT(IN) :: has_unit_metric
453 :
454 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_loop_do_ot'
455 :
456 : INTEGER :: handle, ispin
457 77848 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
458 : TYPE(dbcsr_type), POINTER :: orthogonality_metric
459 :
460 77848 : CALL timeset(routineN, handle)
461 77848 : NULLIFY (rho_ao)
462 :
463 77848 : CALL qs_rho_get(rho, rho_ao=rho_ao)
464 :
465 77848 : IF (has_unit_metric) THEN
466 18424 : NULLIFY (orthogonality_metric)
467 : ELSE
468 59424 : orthogonality_metric => matrix_s(1)%matrix
469 : END IF
470 :
471 : ! in case of LSD the first spin qs_ot_env will drive the minimization
472 : ! in the case of a restricted calculation, it will make sure the spin orbitals are equal
473 :
474 : CALL ot_scf_mini(mos, mo_derivs, smear, orthogonality_metric, &
475 : total_energy, energy_only, scf_env%iter_delta, &
476 77848 : scf_env%qs_ot_env)
477 :
478 168973 : DO ispin = 1, SIZE(mos)
479 168973 : CALL set_mo_occupation(mo_set=mos(ispin), smear=smear)
480 : END DO
481 :
482 168973 : DO ispin = 1, SIZE(mos)
483 : CALL calculate_density_matrix(mos(ispin), &
484 : rho_ao(ispin)%matrix, &
485 168973 : use_dbcsr=.TRUE.)
486 : END DO
487 :
488 77848 : scf_env%iter_method = scf_env%qs_ot_env(1)%OT_METHOD_FULL
489 77848 : scf_env%iter_param = scf_env%qs_ot_env(1)%ds_min
490 77848 : qs_env%broyden_adaptive_sigma = scf_env%qs_ot_env(1)%broyden_adaptive_sigma
491 :
492 77848 : CALL timestop(handle)
493 :
494 77848 : END SUBROUTINE qs_scf_loop_do_ot
495 :
496 : ! **************************************************************************************************
497 : !> \brief Performs the requested density mixing if any needed
498 : !> \param scf_env Holds SCF environment information
499 : !> \param rho All data for the electron density
500 : !> \param para_env Parallel environment
501 : !> \param diis_step Did we do a DIIS step?
502 : ! **************************************************************************************************
503 224043 : SUBROUTINE qs_scf_density_mixing(scf_env, rho, para_env, diis_step)
504 : TYPE(qs_scf_env_type), POINTER :: scf_env
505 : TYPE(qs_rho_type), POINTER :: rho
506 : TYPE(mp_para_env_type), POINTER :: para_env
507 : LOGICAL :: diis_step
508 :
509 224043 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
510 :
511 224043 : NULLIFY (rho_ao_kp)
512 :
513 224043 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
514 :
515 364502 : SELECT CASE (scf_env%mixing_method)
516 : CASE (direct_mixing_nr)
517 : CALL scf_env_density_mixing(scf_env%p_mix_new, &
518 : scf_env%mixing_store, rho_ao_kp, para_env, scf_env%iter_delta, scf_env%iter_count, &
519 140459 : diis=diis_step)
520 : CASE (gspace_mixing_nr, pulay_mixing_nr, broyden_mixing_nr, modified_broyden_mixing_nr, &
521 : multisecant_mixing_nr, new_pulay_mixing_nr)
522 : ! Compute the difference p_out-p_in
523 : CALL self_consistency_check(rho_ao_kp, scf_env%p_delta, para_env, scf_env%p_mix_new, &
524 5736 : delta=scf_env%iter_delta)
525 : CASE (no_mixing_nr)
526 : CASE DEFAULT
527 : CALL cp_abort(__LOCATION__, &
528 : "unknown scf mixing method: "// &
529 224043 : cp_to_string(scf_env%mixing_method))
530 : END SELECT
531 :
532 224043 : END SUBROUTINE qs_scf_density_mixing
533 :
534 : ! **************************************************************************************************
535 : !> \brief checks whether exit conditions for outer loop are satisfied
536 : !> \param qs_env ...
537 : !> \param scf_env ...
538 : !> \param scf_control ...
539 : !> \param should_stop ...
540 : !> \param outer_loop_converged ...
541 : !> \param exit_outer_loop ...
542 : ! **************************************************************************************************
543 25027 : SUBROUTINE qs_scf_check_outer_exit(qs_env, scf_env, scf_control, should_stop, &
544 : outer_loop_converged, exit_outer_loop)
545 : TYPE(qs_environment_type), POINTER :: qs_env
546 : TYPE(qs_scf_env_type), POINTER :: scf_env
547 : TYPE(scf_control_type), POINTER :: scf_control
548 : LOGICAL :: should_stop, outer_loop_converged, &
549 : exit_outer_loop
550 :
551 : REAL(KIND=dp) :: outer_loop_eps
552 :
553 25027 : outer_loop_converged = .TRUE.
554 25027 : IF (scf_control%outer_scf%have_scf) THEN
555 : ! We have an outer SCF loop...
556 6127 : scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count + 1
557 6127 : outer_loop_converged = .FALSE.
558 :
559 6127 : CALL outer_loop_gradient(qs_env, scf_env)
560 : ! Multiple constraints: get largest deviation
561 12340 : outer_loop_eps = SQRT(MAXVAL(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count)**2))
562 :
563 6127 : IF (outer_loop_eps < scf_control%outer_scf%eps_scf) outer_loop_converged = .TRUE.
564 : END IF
565 :
566 : exit_outer_loop = should_stop .OR. outer_loop_converged .OR. &
567 25027 : scf_env%outer_scf%iter_count > scf_control%outer_scf%max_scf
568 :
569 25027 : END SUBROUTINE qs_scf_check_outer_exit
570 :
571 : ! **************************************************************************************************
572 : !> \brief checks whether exit conditions for inner loop are satisfied
573 : !> \param qs_env ...
574 : !> \param scf_env ...
575 : !> \param scf_control ...
576 : !> \param should_stop ...
577 : !> \param just_energy ...
578 : !> \param exit_inner_loop ...
579 : !> \param inner_loop_converged ...
580 : !> \param output_unit ...
581 : ! **************************************************************************************************
582 224043 : SUBROUTINE qs_scf_check_inner_exit(qs_env, scf_env, scf_control, should_stop, just_energy, &
583 : exit_inner_loop, inner_loop_converged, output_unit)
584 : TYPE(qs_environment_type), POINTER :: qs_env
585 : TYPE(qs_scf_env_type), POINTER :: scf_env
586 : TYPE(scf_control_type), POINTER :: scf_control
587 : LOGICAL :: should_stop, just_energy, &
588 : exit_inner_loop, inner_loop_converged
589 : INTEGER :: output_unit
590 :
591 224043 : inner_loop_converged = .FALSE.
592 224043 : exit_inner_loop = .FALSE.
593 :
594 : CALL external_control(should_stop, "SCF", target_time=qs_env%target_time, &
595 224043 : start_time=qs_env%start_time)
596 224043 : IF (scf_env%iter_delta < scf_control%eps_scf) THEN
597 21251 : IF (output_unit > 0) THEN
598 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I5,A/)") &
599 10802 : "*** SCF run converged in ", scf_env%iter_count, " steps ***"
600 : END IF
601 21251 : inner_loop_converged = .TRUE.
602 21251 : exit_inner_loop = .TRUE.
603 202792 : ELSE IF (should_stop .OR. scf_env%iter_count >= scf_control%max_scf) THEN
604 4094 : inner_loop_converged = .FALSE.
605 4094 : IF (just_energy) THEN
606 886 : exit_inner_loop = .FALSE.
607 : ELSE
608 3208 : exit_inner_loop = .TRUE.
609 3208 : IF (output_unit > 0) THEN
610 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I5,A/)") &
611 1610 : "Leaving inner SCF loop after reaching ", scf_env%iter_count, " steps."
612 : END IF
613 : END IF
614 : END IF
615 :
616 224043 : END SUBROUTINE qs_scf_check_inner_exit
617 :
618 : ! **************************************************************************************************
619 : !> \brief undoing density mixing. Important upon convergence
620 : !> \param scf_env ...
621 : !> \param rho ...
622 : !> \param dft_control ...
623 : !> \param para_env ...
624 : !> \param diis_step ...
625 : ! **************************************************************************************************
626 24459 : SUBROUTINE qs_scf_undo_mixing(scf_env, rho, dft_control, para_env, diis_step)
627 : TYPE(qs_scf_env_type), POINTER :: scf_env
628 : TYPE(qs_rho_type), POINTER :: rho
629 : TYPE(dft_control_type), POINTER :: dft_control
630 : TYPE(mp_para_env_type), POINTER :: para_env
631 : LOGICAL :: diis_step
632 :
633 : CHARACTER(len=default_string_length) :: name
634 : INTEGER :: ic, ispin, nc
635 24459 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
636 :
637 24459 : NULLIFY (rho_ao_kp)
638 :
639 24459 : IF (scf_env%mixing_method > 0) THEN
640 17074 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
641 17074 : nc = SIZE(scf_env%p_mix_new, 2)
642 33460 : SELECT CASE (scf_env%mixing_method)
643 : CASE (direct_mixing_nr)
644 : CALL scf_env_density_mixing(scf_env%p_mix_new, scf_env%mixing_store, &
645 : rho_ao_kp, para_env, scf_env%iter_delta, &
646 : scf_env%iter_count, diis=diis_step, &
647 16386 : invert=.TRUE.)
648 159434 : DO ic = 1, nc
649 319408 : DO ispin = 1, dft_control%nspins
650 159974 : CALL dbcsr_get_info(rho_ao_kp(ispin, ic)%matrix, name=name) ! keep the name
651 303022 : CALL dbcsr_copy(rho_ao_kp(ispin, ic)%matrix, scf_env%p_mix_new(ispin, ic)%matrix, name=name)
652 : END DO
653 : END DO
654 : CASE (gspace_mixing_nr, pulay_mixing_nr, broyden_mixing_nr, modified_broyden_mixing_nr, &
655 : multisecant_mixing_nr, new_pulay_mixing_nr)
656 60718 : DO ic = 1, nc
657 87292 : DO ispin = 1, dft_control%nspins
658 43648 : CALL dbcsr_get_info(rho_ao_kp(ispin, ic)%matrix, name=name) ! keep the name
659 86604 : CALL dbcsr_copy(rho_ao_kp(ispin, ic)%matrix, scf_env%p_mix_new(ispin, ic)%matrix, name=name)
660 : END DO
661 : END DO
662 : END SELECT
663 : END IF
664 24459 : END SUBROUTINE qs_scf_undo_mixing
665 :
666 : ! **************************************************************************************************
667 : !> \brief Performs the updates rho (takes care of mixing as well)
668 : !> \param rho ...
669 : !> \param qs_env ...
670 : !> \param scf_env ...
671 : !> \param ks_env ...
672 : !> \param mix_rho ...
673 : ! **************************************************************************************************
674 224043 : SUBROUTINE qs_scf_rho_update(rho, qs_env, scf_env, ks_env, mix_rho)
675 : TYPE(qs_rho_type), POINTER :: rho
676 : TYPE(qs_environment_type), POINTER :: qs_env
677 : TYPE(qs_scf_env_type), POINTER :: scf_env
678 : TYPE(qs_ks_env_type), POINTER :: ks_env
679 : LOGICAL, INTENT(IN) :: mix_rho
680 :
681 : TYPE(mp_para_env_type), POINTER :: para_env
682 :
683 224043 : NULLIFY (para_env)
684 224043 : CALL get_qs_env(qs_env, para_env=para_env)
685 : ! ** update qs_env%rho
686 224043 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
687 : ! ** Density mixing through density matrix or on the reciprocal space grid (exclusive)
688 224043 : IF (mix_rho) THEN
689 : CALL gspace_mixing(qs_env, scf_env%mixing_method, scf_env%mixing_store, rho, &
690 5048 : para_env, scf_env%iter_count)
691 :
692 : END IF
693 224043 : CALL qs_ks_did_change(ks_env, rho_changed=.TRUE.)
694 :
695 224043 : END SUBROUTINE qs_scf_rho_update
696 :
697 : ! **************************************************************************************************
698 : !> \brief Performs the necessary steps before leaving innner scf loop
699 : !> \param scf_env ...
700 : !> \param qs_env ...
701 : !> \param diis_step ...
702 : !> \param output_unit ...
703 : ! **************************************************************************************************
704 24459 : SUBROUTINE qs_scf_inner_finalize(scf_env, qs_env, diis_step, output_unit)
705 : TYPE(qs_scf_env_type), POINTER :: scf_env
706 : TYPE(qs_environment_type), POINTER :: qs_env
707 : LOGICAL :: diis_step
708 : INTEGER, INTENT(IN) :: output_unit
709 :
710 : LOGICAL :: do_kpoints
711 : TYPE(dft_control_type), POINTER :: dft_control
712 : TYPE(mp_para_env_type), POINTER :: para_env
713 : TYPE(qs_energy_type), POINTER :: energy
714 : TYPE(qs_ks_env_type), POINTER :: ks_env
715 : TYPE(qs_rho_type), POINTER :: rho
716 :
717 24459 : NULLIFY (energy, rho, dft_control, ks_env)
718 :
719 : CALL get_qs_env(qs_env=qs_env, energy=energy, ks_env=ks_env, &
720 : rho=rho, dft_control=dft_control, para_env=para_env, &
721 24459 : do_kpoints=do_kpoints)
722 :
723 24459 : CALL cleanup_scf_loop(scf_env)
724 :
725 : ! now, print out energies and charges corresponding to the obtained wfn
726 : ! (this actually is not 100% consistent at this point)!
727 24459 : CALL qs_scf_print_summary(output_unit, qs_env)
728 :
729 24459 : CALL qs_scf_undo_mixing(scf_env, rho, dft_control, para_env, diis_step)
730 :
731 : ! *** update rspace rho since the mo changed
732 : ! *** this might not always be needed (i.e. no post calculation / no forces )
733 : ! *** but guarantees that rho and wfn are consistent at this point
734 24459 : CALL qs_scf_rho_update(rho, qs_env, scf_env, ks_env, mix_rho=.FALSE.)
735 :
736 24459 : END SUBROUTINE qs_scf_inner_finalize
737 :
738 : ! **************************************************************************************************
739 : !> \brief perform cleanup operations at the end of an scf loop
740 : !> \param scf_env ...
741 : !> \par History
742 : !> 03.2006 created [Joost VandeVondele]
743 : ! **************************************************************************************************
744 24459 : SUBROUTINE cleanup_scf_loop(scf_env)
745 : TYPE(qs_scf_env_type), INTENT(INOUT) :: scf_env
746 :
747 : CHARACTER(len=*), PARAMETER :: routineN = 'cleanup_scf_loop'
748 :
749 : INTEGER :: handle, ispin
750 :
751 24459 : CALL timeset(routineN, handle)
752 :
753 31844 : SELECT CASE (scf_env%method)
754 : CASE (ot_method_nr)
755 16079 : DO ispin = 1, SIZE(scf_env%qs_ot_env)
756 16079 : CALL ot_scf_destroy(scf_env%qs_ot_env(ispin))
757 : END DO
758 7385 : DEALLOCATE (scf_env%qs_ot_env)
759 : CASE (ot_diag_method_nr)
760 : !
761 : CASE (general_diag_method_nr)
762 : !
763 : CASE (special_diag_method_nr)
764 : !
765 : CASE (block_krylov_diag_method_nr, block_davidson_diag_method_nr)
766 : !
767 : CASE (filter_matrix_diag_method_nr)
768 : !
769 : CASE (smeagol_method_nr)
770 : !
771 : CASE DEFAULT
772 : CALL cp_abort(__LOCATION__, &
773 : "unknown scf method method:"// &
774 24459 : cp_to_string(scf_env%method))
775 : END SELECT
776 :
777 24459 : CALL timestop(handle)
778 :
779 24459 : END SUBROUTINE cleanup_scf_loop
780 :
781 : END MODULE qs_scf_loop_utils
|
