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