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