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