Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines for the Quickstep SCF run.
10 : !> \par History
11 : !> - Joost VandeVondele (02.2002)
12 : !> added code for: incremental (pab and gvg) update
13 : !> initialisation (init_cube, l_info)
14 : !> - Joost VandeVondele (02.2002)
15 : !> called the poisson code of the classical part
16 : !> this takes into account the spherical cutoff and allows for
17 : !> isolated systems
18 : !> - Joost VandeVondele (02.2002)
19 : !> added multiple grid feature
20 : !> changed to spherical cutoff consistently (?)
21 : !> therefore removed the gradient correct functionals
22 : !> - updated with the new QS data structures (10.04.02,MK)
23 : !> - copy_matrix replaced by transfer_matrix (11.04.02,MK)
24 : !> - nrebuild_rho and nrebuild_gvg unified (12.04.02,MK)
25 : !> - set_mo_occupation for smearing of the MO occupation numbers
26 : !> (17.04.02,MK)
27 : !> - MO level shifting added (22.04.02,MK)
28 : !> - Usage of TYPE mo_set_p_type
29 : !> - Joost VandeVondele (05.2002)
30 : !> added cholesky based diagonalisation
31 : !> - 05.2002 added pao method [fawzi]
32 : !> - parallel FFT (JGH 22.05.2002)
33 : !> - 06.2002 moved KS matrix construction to qs_build_KS_matrix.F [fawzi]
34 : !> - started to include more LSD (01.2003,Joost VandeVondele)
35 : !> - 02.2003 scf_env [fawzi]
36 : !> - got rid of nrebuild (01.2004, Joost VandeVondele)
37 : !> - 10.2004 removed pao [fawzi]
38 : !> - 03.2006 large cleaning action [Joost VandeVondele]
39 : !> - High-spin ROKS added (05.04.06,MK)
40 : !> - Mandes (10.2013)
41 : !> intermediate energy communication with external communicator added
42 : !> - kpoints (08.2014, JGH)
43 : !> - unified k-point and gamma-point code (2014.11) [Ole Schuett]
44 : !> - added extra SCF loop for CDFT constraints (12.2015) [Nico Holmberg]
45 : !> \author Matthias Krack (30.04.2001)
46 : ! **************************************************************************************************
47 : MODULE qs_scf
48 : USE atomic_kind_types, ONLY: atomic_kind_type
49 : USE cp_blacs_env, ONLY: cp_blacs_env_type
50 : USE cp_control_types, ONLY: dft_control_type
51 : USE cp_dbcsr_api, ONLY: &
52 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_get_info, &
53 : dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, dbcsr_release_p, dbcsr_set, dbcsr_type, &
54 : dbcsr_type_no_symmetry
55 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
56 : copy_fm_to_dbcsr,&
57 : cp_dbcsr_m_by_n_from_row_template,&
58 : dbcsr_deallocate_matrix_set
59 : USE cp_files, ONLY: close_file
60 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type
61 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
62 : cp_fm_struct_release,&
63 : cp_fm_struct_type
64 : USE cp_fm_types, ONLY: cp_fm_create,&
65 : cp_fm_init_random,&
66 : cp_fm_release,&
67 : cp_fm_to_fm,&
68 : cp_fm_type
69 : USE cp_log_handling, ONLY: cp_add_default_logger,&
70 : cp_get_default_logger,&
71 : cp_logger_release,&
72 : cp_logger_type,&
73 : cp_rm_default_logger,&
74 : cp_to_string
75 : USE cp_output_handling, ONLY: cp_add_iter_level,&
76 : cp_iterate,&
77 : cp_p_file,&
78 : cp_print_key_should_output,&
79 : cp_print_key_unit_nr,&
80 : cp_rm_iter_level
81 : USE cp_result_methods, ONLY: get_results,&
82 : test_for_result
83 : USE cp_result_types, ONLY: cp_result_type
84 : USE ec_env_types, ONLY: energy_correction_type
85 : USE input_constants, ONLY: &
86 : broyden_type_1, broyden_type_1_explicit, broyden_type_1_explicit_ls, broyden_type_1_ls, &
87 : broyden_type_2, broyden_type_2_explicit, broyden_type_2_explicit_ls, broyden_type_2_ls, &
88 : cdft2ot, cholesky_dbcsr, diag_update_method_adiis, history_guess, ot2cdft, &
89 : ot_precond_fermi_low_rank, ot_precond_full_all, ot_precond_full_all_covariant, &
90 : ot_precond_full_kinetic, ot_precond_full_single, ot_precond_full_single_inverse, &
91 : ot_precond_none, ot_precond_s_inverse, outer_scf_becke_constraint, &
92 : outer_scf_hirshfeld_constraint, outer_scf_optimizer_broyden, &
93 : outer_scf_optimizer_newton_ls, sic_eo, smear_fermi_dirac, smear_gaussian, smear_mp, &
94 : smear_mv, tblite_scc_mixer_tblite
95 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
96 : section_vals_type
97 : USE kinds, ONLY: default_path_length,&
98 : default_string_length,&
99 : dp
100 : USE kpoint_io, ONLY: write_kpoints_restart
101 : USE kpoint_methods, ONLY: kpoint_initialize_mo_set,&
102 : kpoint_initialize_mos
103 : USE kpoint_types, ONLY: get_kpoint_info,&
104 : kpoint_env_type,&
105 : kpoint_type
106 : USE machine, ONLY: m_flush,&
107 : m_walltime
108 : USE mathlib, ONLY: invert_matrix
109 : USE message_passing, ONLY: mp_comm_type,&
110 : mp_para_env_type
111 : USE particle_types, ONLY: particle_type
112 : USE physcon, ONLY: evolt
113 : USE preconditioner, ONLY: &
114 : make_preconditioner_complex_fermi_low_rank, make_preconditioner_complex_full_all, &
115 : make_preconditioner_complex_full_all_covariant, make_preconditioner_complex_full_kinetic, &
116 : make_preconditioner_complex_full_s_inverse, make_preconditioner_complex_full_single, &
117 : make_preconditioner_complex_full_single_inverse, prepare_preconditioner, &
118 : restart_preconditioner
119 : USE preconditioner_types, ONLY: destroy_preconditioner,&
120 : init_preconditioner
121 : USE pw_env_types, ONLY: pw_env_get,&
122 : pw_env_type
123 : USE pw_pool_types, ONLY: pw_pool_type
124 : USE qs_block_davidson_types, ONLY: block_davidson_deallocate
125 : USE qs_cdft_scf_utils, ONLY: build_diagonal_jacobian,&
126 : create_tmp_logger,&
127 : initialize_inverse_jacobian,&
128 : prepare_jacobian_stencil,&
129 : print_inverse_jacobian,&
130 : restart_inverse_jacobian
131 : USE qs_cdft_types, ONLY: cdft_control_type
132 : USE qs_charge_mixing, ONLY: charge_mixing_scc_error
133 : USE qs_charges_types, ONLY: qs_charges_type
134 : USE qs_density_matrices, ONLY: calculate_density_matrix
135 : USE qs_density_mixing_types, ONLY: gspace_mixing_nr
136 : USE qs_diis, ONLY: qs_diis_b_clear,&
137 : qs_diis_b_clear_kp,&
138 : qs_diis_b_create,&
139 : qs_diis_b_create_kp
140 : USE qs_energy_types, ONLY: qs_energy_type
141 : USE qs_environment_types, ONLY: get_qs_env,&
142 : qs_environment_type,&
143 : set_qs_env
144 : USE qs_fod, ONLY: qs_fod_validate
145 : USE qs_integrate_potential, ONLY: integrate_v_rspace
146 : USE qs_kind_types, ONLY: get_qs_kind,&
147 : qs_kind_type
148 : USE qs_kpoint_operators, ONLY: kpoint_operator_get_local,&
149 : kpoint_operator_store
150 : USE qs_kpoint_state, ONLY: qs_kpoint_copy_spin_mos,&
151 : qs_kpoint_mos_initialized,&
152 : qs_kpoint_set_fixed_occupations,&
153 : qs_kpoint_state_canonicalize_fixed,&
154 : qs_kpoint_state_commit,&
155 : qs_kpoint_state_prepare_fixed_density
156 : USE qs_ks_atom, ONLY: update_ks_atom
157 : USE qs_ks_methods, ONLY: evaluate_core_matrix_traces,&
158 : qs_ks_update_qs_env
159 : USE qs_ks_types, ONLY: get_ks_env,&
160 : qs_ks_did_change,&
161 : qs_ks_env_type
162 : USE qs_local_rho_types, ONLY: local_rho_set_create,&
163 : local_rho_set_release,&
164 : local_rho_type
165 : USE qs_matrix_pools, ONLY: mpools_get,&
166 : mpools_rebuild_fm_pools,&
167 : mpools_release
168 : USE qs_mo_io, ONLY: write_mo_set_to_restart
169 : USE qs_mo_methods, ONLY: make_basis_simple,&
170 : make_basis_sm
171 : USE qs_mo_occupation, ONLY: set_mo_occupation
172 : USE qs_mo_types, ONLY: allocate_mo_set,&
173 : deallocate_mo_set,&
174 : duplicate_mo_set,&
175 : get_mo_set,&
176 : init_mo_set,&
177 : mo_set_restrict,&
178 : mo_set_type,&
179 : reassign_allocated_mos
180 : USE qs_ot, ONLY: qs_ot_get_orbitals_ref_complex,&
181 : qs_ot_get_p,&
182 : qs_ot_get_p_complex,&
183 : qs_ot_new_preconditioner
184 : USE qs_ot_scf, ONLY: ot_scf_destroy,&
185 : ot_scf_init,&
186 : ot_scf_read_input
187 : USE qs_ot_types, ONLY: &
188 : qs_ot_allocate, qs_ot_allocate_complex_state, qs_ot_channel_index, &
189 : qs_ot_check_channel_context, qs_ot_init, qs_ot_kpoint_preconditioner_solver_supported, &
190 : qs_ot_kpoint_preconditioner_supported, qs_ot_number_of_channels, qs_ot_set_context
191 : USE qs_outer_scf, ONLY: outer_loop_gradient,&
192 : outer_loop_optimize,&
193 : outer_loop_purge_history,&
194 : outer_loop_switch,&
195 : outer_loop_update_qs_env
196 : USE qs_rho_atom_methods, ONLY: allocate_rho_atom_internals
197 : USE qs_rho_atom_types, ONLY: zero_rho_atom_integrals
198 : USE qs_rho_methods, ONLY: qs_rho_update_rho
199 : USE qs_rho_types, ONLY: qs_rho_get,&
200 : qs_rho_type
201 : USE qs_scf_diagonalization, ONLY: do_general_diag_kp
202 : USE qs_scf_initialization, ONLY: qs_scf_env_initialize
203 : USE qs_scf_loop_utils, ONLY: &
204 : qs_scf_candidate_density_delta, qs_scf_check_inner_exit, qs_scf_check_outer_exit, &
205 : qs_scf_commit_density_candidate, qs_scf_density_mixing, qs_scf_inner_finalize, &
206 : qs_scf_kp_search_endpoint, qs_scf_new_mos, qs_scf_new_mos_kp, qs_scf_rho_update, &
207 : qs_scf_set_loop_flags
208 : USE qs_scf_output, ONLY: qs_scf_cdft_info,&
209 : qs_scf_cdft_initial_info,&
210 : qs_scf_gce_info,&
211 : qs_scf_loop_info,&
212 : qs_scf_loop_print,&
213 : qs_scf_outer_loop_info,&
214 : qs_scf_write_mos
215 : USE qs_scf_post_scf, ONLY: qs_scf_compute_properties
216 : USE qs_scf_subspace, ONLY: qs_scf_subspace_build,&
217 : qs_scf_subspace_push,&
218 : qs_scf_subspace_restart,&
219 : qs_scf_subspace_update_shift
220 : USE qs_scf_subspace_types, ONLY: qs_scf_subspace_buffer_clear,&
221 : qs_scf_subspace_buffer_create,&
222 : qs_scf_subspace_buffer_release
223 : USE qs_scf_types, ONLY: &
224 : block_davidson_diag_method_nr, block_krylov_diag_method_nr, filter_matrix_diag_method_nr, &
225 : general_diag_method_nr, ot_diag_method_nr, ot_method_nr, qs_scf_env_type, &
226 : smeagol_method_nr, special_diag_method_nr
227 : USE qs_vxc_atom, ONLY: gapw_cdft_one_center
228 : USE qs_wf_history_methods, ONLY: wfi_purge_history,&
229 : wfi_update
230 : USE scf_control_types, ONLY: scf_control_type
231 : USE smeagol_interface, ONLY: run_smeagol_bulktrans,&
232 : run_smeagol_emtrans
233 : USE tblite_interface, ONLY: tb_get_energy,&
234 : tb_native_scc_mixer_active,&
235 : tb_scf_mixer_error,&
236 : tb_update_charges
237 : #include "./base/base_uses.f90"
238 :
239 : IMPLICIT NONE
240 :
241 : PRIVATE
242 :
243 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf'
244 : INTEGER, PARAMETER, PRIVATE :: kp_ot_entry_initial = 1, &
245 : kp_ot_entry_reuse = 2, &
246 : kp_ot_entry_resized = 3, &
247 : kp_ot_entry_refresh = 4
248 : LOGICAL, PRIVATE :: reuse_precond = .FALSE.
249 : LOGICAL, PRIVATE :: used_history = .FALSE.
250 :
251 : PUBLIC :: scf, scf_env_cleanup, scf_env_do_scf, cdft_scf, init_scf_loop
252 :
253 : CONTAINS
254 :
255 : ! **************************************************************************************************
256 : !> \brief perform an scf procedure in the given qs_env
257 : !> \param qs_env the qs_environment where to perform the scf procedure
258 : !> \param has_converged ...
259 : !> \param total_scf_steps ...
260 : !> \par History
261 : !> 02.2003 introduced scf_env, moved real work to scf_env_do_scf [fawzi]
262 : !> \author fawzi
263 : !> \note
264 : ! **************************************************************************************************
265 24791 : SUBROUTINE scf(qs_env, has_converged, total_scf_steps)
266 : TYPE(qs_environment_type), POINTER :: qs_env
267 : LOGICAL, INTENT(OUT), OPTIONAL :: has_converged
268 : INTEGER, INTENT(OUT), OPTIONAL :: total_scf_steps
269 :
270 : INTEGER :: ihistory, max_scf_tmp, tsteps
271 : LOGICAL :: converged, outer_scf_loop, should_stop
272 : LOGICAL, SAVE :: first_step_flag = .TRUE.
273 24791 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient_history, variable_history
274 : TYPE(cp_logger_type), POINTER :: logger
275 : TYPE(dft_control_type), POINTER :: dft_control
276 : TYPE(qs_scf_env_type), POINTER :: scf_env
277 : TYPE(scf_control_type), POINTER :: scf_control
278 : TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
279 :
280 24791 : NULLIFY (scf_env)
281 49582 : logger => cp_get_default_logger()
282 24791 : CPASSERT(ASSOCIATED(qs_env))
283 24791 : IF (PRESENT(has_converged)) THEN
284 0 : has_converged = .FALSE.
285 : END IF
286 24791 : IF (PRESENT(total_scf_steps)) THEN
287 0 : total_scf_steps = 0
288 : END IF
289 : CALL get_qs_env(qs_env, scf_env=scf_env, input=input, &
290 24791 : dft_control=dft_control, scf_control=scf_control)
291 24791 : CALL qs_fod_validate(input, logger, qs_env)
292 24791 : qs_env%scf_convergence_available = .FALSE.
293 24791 : qs_env%scf_converged = .FALSE.
294 24791 : IF (scf_control%max_scf > 0) THEN
295 :
296 24149 : dft_section => section_vals_get_subs_vals(input, "DFT")
297 24149 : scf_section => section_vals_get_subs_vals(dft_section, "SCF")
298 :
299 24149 : IF (.NOT. ASSOCIATED(scf_env)) THEN
300 7143 : CALL qs_scf_env_initialize(qs_env, scf_env)
301 : ! Moved here from qs_scf_env_initialize to be able to have more scf_env
302 7143 : CALL set_qs_env(qs_env, scf_env=scf_env)
303 : ELSE
304 17006 : CALL qs_scf_env_initialize(qs_env, scf_env)
305 : END IF
306 :
307 24149 : IF ((scf_control%density_guess == history_guess) .AND. (first_step_flag)) THEN
308 2 : max_scf_tmp = scf_control%max_scf
309 2 : scf_control%max_scf = 1
310 2 : outer_scf_loop = scf_control%outer_scf%have_scf
311 2 : scf_control%outer_scf%have_scf = .FALSE.
312 : END IF
313 :
314 24149 : IF (.NOT. dft_control%qs_control%cdft) THEN
315 : CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
316 23789 : converged=converged, should_stop=should_stop, total_scf_steps=tsteps)
317 : ELSE
318 : ! Third SCF loop needed for CDFT with OT to properly restart OT inner loop
319 : CALL cdft_scf(qs_env=qs_env, should_stop=should_stop, &
320 360 : has_converged=converged, total_scf_steps=tsteps)
321 : END IF
322 :
323 24149 : qs_env%scf_convergence_available = .TRUE.
324 24149 : qs_env%scf_converged = converged
325 :
326 : ! If SCF has not converged, then we should not start MP2
327 24149 : IF (ASSOCIATED(qs_env%mp2_env)) qs_env%mp2_env%hf_fail = .NOT. converged
328 :
329 : ! Add the converged outer_scf SCF gradient(s)/variable(s) to history
330 24149 : IF (scf_control%outer_scf%have_scf) THEN
331 4347 : ihistory = scf_env%outer_scf%iter_count
332 : CALL get_qs_env(qs_env, gradient_history=gradient_history, &
333 4347 : variable_history=variable_history)
334 : ! We only store the latest two values
335 8726 : gradient_history(:, 1) = gradient_history(:, 2)
336 17452 : gradient_history(:, 2) = scf_env%outer_scf%gradient(:, ihistory)
337 8726 : variable_history(:, 1) = variable_history(:, 2)
338 17452 : variable_history(:, 2) = scf_env%outer_scf%variables(:, ihistory)
339 : ! Reset flag
340 4347 : IF (used_history) used_history = .FALSE.
341 : ! Update a counter and check if the Jacobian should be deallocated
342 4347 : IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
343 64 : scf_control%outer_scf%cdft_opt_control%ijacobian(2) = scf_control%outer_scf%cdft_opt_control%ijacobian(2) + 1
344 : IF (scf_control%outer_scf%cdft_opt_control%ijacobian(2) >= &
345 64 : scf_control%outer_scf%cdft_opt_control%jacobian_freq(2) .AND. &
346 : scf_control%outer_scf%cdft_opt_control%jacobian_freq(2) > 0) THEN
347 50 : scf_env%outer_scf%deallocate_jacobian = .TRUE.
348 : END IF
349 : END IF
350 : END IF
351 : ! *** add the converged wavefunction to the wavefunction history
352 24149 : IF ((ASSOCIATED(qs_env%wf_history)) .AND. &
353 : ((scf_control%density_guess /= history_guess) .OR. &
354 : (.NOT. first_step_flag))) THEN
355 24147 : IF (.NOT. dft_control%qs_control%cdft) THEN
356 : ! No next geometry step in a standalone single-point calculation.
357 23787 : IF (.NOT. qs_env%skip_wf_history) THEN
358 20675 : CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
359 : END IF
360 : ELSE
361 360 : IF (dft_control%qs_control%cdft_control%should_purge) THEN
362 0 : CALL wfi_purge_history(qs_env)
363 0 : CALL outer_loop_purge_history(qs_env)
364 0 : dft_control%qs_control%cdft_control%should_purge = .FALSE.
365 : ELSE
366 360 : CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
367 : END IF
368 : END IF
369 2 : ELSE IF ((scf_control%density_guess == history_guess) .AND. &
370 : (first_step_flag)) THEN
371 2 : scf_control%max_scf = max_scf_tmp
372 2 : scf_control%outer_scf%have_scf = outer_scf_loop
373 2 : first_step_flag = .FALSE.
374 : END IF
375 :
376 : ! *** compute properties that depend on the converged wavefunction
377 24149 : IF (.NOT. (should_stop)) CALL qs_scf_compute_properties(qs_env)
378 :
379 : ! *** SMEAGOL interface ***
380 24149 : IF (.NOT. (should_stop)) THEN
381 : ! compute properties that depend on the converged wavefunction ..
382 24149 : CALL run_smeagol_emtrans(qs_env, last=.TRUE., iter=0)
383 : ! .. or save matrices related to bulk leads
384 24149 : CALL run_smeagol_bulktrans(qs_env)
385 : END IF
386 :
387 : ! *** cleanup
388 24149 : CALL scf_env_cleanup(scf_env)
389 24149 : IF (dft_control%qs_control%cdft) THEN
390 360 : CALL cdft_control_cleanup(dft_control%qs_control%cdft_control)
391 : END IF
392 :
393 24149 : IF (PRESENT(has_converged)) THEN
394 0 : has_converged = converged
395 : END IF
396 24149 : IF (PRESENT(total_scf_steps)) THEN
397 0 : total_scf_steps = tsteps
398 : END IF
399 :
400 : END IF
401 :
402 24791 : END SUBROUTINE scf
403 :
404 : ! **************************************************************************************************
405 : !> \brief perform an scf loop
406 : !> \param scf_env the scf_env where to perform the scf procedure
407 : !> \param scf_control ...
408 : !> \param qs_env the qs_env, the scf_env lives in
409 : !> \param converged will be true / false if converged is reached
410 : !> \param should_stop ...
411 : !> \param total_scf_steps ...
412 : !> \par History
413 : !> long history, see cvs and qs_scf module history
414 : !> 02.2003 introduced scf_env [fawzi]
415 : !> 09.2005 Frozen density approximation [TdK]
416 : !> 06.2007 Check for SCF iteration count early [jgh]
417 : !> 10.2019 switch_surf_dip [SGh]
418 : !> \author Matthias Krack
419 : !> \note
420 : ! **************************************************************************************************
421 24501 : SUBROUTINE scf_env_do_scf(scf_env, scf_control, qs_env, converged, should_stop, total_scf_steps)
422 :
423 : TYPE(qs_scf_env_type), POINTER :: scf_env
424 : TYPE(scf_control_type), POINTER :: scf_control
425 : TYPE(qs_environment_type), POINTER :: qs_env
426 : LOGICAL, INTENT(OUT) :: converged, should_stop
427 : INTEGER, INTENT(OUT) :: total_scf_steps
428 :
429 : CHARACTER(LEN=*), PARAMETER :: routineN = 'scf_env_do_scf'
430 : INTEGER, PARAMETER :: max_ot_kp_subspace_refreshes = 4
431 : REAL(KIND=dp), PARAMETER :: adiis_stagnation_ratio = 1.0E-2_dp
432 :
433 : CHARACTER(LEN=default_string_length) :: description, name
434 : INTEGER :: accepted_ot_kp_searches, ext_master_id, handle, handle2, i_tmp, ic, ispin, &
435 : iter_count, kp_ot_entry_reason, ot_kp_subspace_refresh_count, &
436 : ot_kp_subspace_refresh_iter_count, output_unit, scf_energy_message_tag, total_steps
437 : LOGICAL :: added_mos_auto_grow, adiis_pushed, adiis_restarted, adiis_stagnated, &
438 : adiis_validation, density_full_step, diis_step, do_kpoints, energy_only, exit_inner_loop, &
439 : exit_outer_loop, inner_loop_converged, internal_tblite_density_full_step, &
440 : internal_tblite_mixer, just_energy, ot_kp_subspace_refresh, &
441 : ot_kp_subspace_refresh_pending, outer_loop_converged, tblite_native_mixer
442 : REAL(KIND=dp) :: adiis_raw_delta, adiis_step_delta, t1, t2
443 : REAL(KIND=dp), DIMENSION(3) :: res_val_3
444 24501 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
445 : TYPE(cp_logger_type), POINTER :: logger
446 : TYPE(cp_result_type), POINTER :: results
447 24501 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
448 24501 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, rho_ao_kp
449 : TYPE(dft_control_type), POINTER :: dft_control
450 : TYPE(energy_correction_type), POINTER :: ec_env
451 : TYPE(kpoint_type), POINTER :: kpoints
452 24501 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_last_converged
453 : TYPE(mp_comm_type) :: external_comm
454 : TYPE(mp_para_env_type), POINTER :: para_env
455 24501 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
456 : TYPE(pw_env_type), POINTER :: pw_env
457 : TYPE(qs_charges_type), POINTER :: qs_charges
458 : TYPE(qs_energy_type), POINTER :: energy
459 24501 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
460 : TYPE(qs_ks_env_type), POINTER :: ks_env
461 : TYPE(qs_rho_type), POINTER :: rho
462 : TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
463 :
464 : ! Weak orbital metrics can need several bounded REF rebuilds to close the canonical density tail.
465 :
466 24501 : CALL timeset(routineN, handle)
467 :
468 24501 : NULLIFY (dft_control, rho, energy, &
469 24501 : logger, qs_charges, ks_env, mos, atomic_kind_set, qs_kind_set, &
470 24501 : particle_set, dft_section, input, &
471 24501 : scf_section, para_env, results, kpoints, pw_env, matrix_ks, &
472 24501 : matrix_ks_kp, rho_ao_kp, mos_last_converged)
473 :
474 24501 : CPASSERT(ASSOCIATED(scf_env))
475 24501 : CPASSERT(ASSOCIATED(qs_env))
476 :
477 24501 : logger => cp_get_default_logger()
478 24501 : t1 = m_walltime()
479 :
480 : CALL get_qs_env(qs_env=qs_env, &
481 : energy=energy, &
482 : particle_set=particle_set, &
483 : qs_charges=qs_charges, &
484 : ks_env=ks_env, &
485 : atomic_kind_set=atomic_kind_set, &
486 : qs_kind_set=qs_kind_set, &
487 : rho=rho, &
488 : mos=mos, &
489 : matrix_ks_kp=matrix_ks_kp, &
490 : input=input, &
491 : dft_control=dft_control, &
492 : do_kpoints=do_kpoints, &
493 : kpoints=kpoints, &
494 : results=results, &
495 : pw_env=pw_env, &
496 24501 : para_env=para_env)
497 : tblite_native_mixer = dft_control%qs_control%xtb_control%do_tblite .AND. &
498 : scf_env%method /= ot_method_nr .AND. &
499 24501 : tb_native_scc_mixer_active(dft_control)
500 : internal_tblite_mixer = (dft_control%qs_control%dftb .AND. &
501 : dft_control%qs_control%dftb_control%tblite_scc_mixer == tblite_scc_mixer_tblite) .OR. &
502 : (dft_control%qs_control%xtb .AND. &
503 : .NOT. dft_control%qs_control%xtb_control%do_tblite .AND. &
504 24501 : dft_control%qs_control%xtb_control%tblite_scc_mixer == tblite_scc_mixer_tblite)
505 : internal_tblite_density_full_step = dft_control%qs_control%xtb .AND. &
506 : .NOT. dft_control%qs_control%xtb_control%do_tblite .AND. &
507 24501 : dft_control%qs_control%xtb_control%tblite_scc_mixer == tblite_scc_mixer_tblite
508 :
509 24501 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
510 :
511 24501 : dft_section => section_vals_get_subs_vals(input, "DFT")
512 24501 : scf_section => section_vals_get_subs_vals(dft_section, "SCF")
513 :
514 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", &
515 24501 : extension=".scfLog")
516 :
517 24501 : IF (scf_control%gce%do_gce .AND. output_unit > 0) THEN
518 1 : WRITE (UNIT=output_unit, FMT="(/,T2,78('-'))")
519 : WRITE (UNIT=output_unit, FMT="(T31,A)") &
520 1 : "GRAND-CANONICAL SCF"
521 : WRITE (UNIT=output_unit, FMT="(T20,A,F12.6,A)") &
522 1 : "Target work function (TWF):", &
523 2 : evolt*scf_control%gce%target_workfunction, " eV"
524 1 : WRITE (UNIT=output_unit, FMT="(T2,78('-'))")
525 : END IF
526 :
527 24501 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
528 12414 : "SCF WAVEFUNCTION OPTIMIZATION"
529 :
530 : ! when switch_surf_dip is switched on, indicate storing mos from the last converged step
531 24501 : IF (dft_control%switch_surf_dip) THEN
532 2 : CALL get_qs_env(qs_env, mos_last_converged=mos_last_converged)
533 4 : DO ispin = 1, dft_control%nspins
534 4 : CALL reassign_allocated_mos(mos(ispin), mos_last_converged(ispin))
535 : END DO
536 2 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
537 1 : "COPIED mos_last_converged ---> mos"
538 : END IF
539 :
540 24501 : IF ((output_unit > 0) .AND. (.NOT. scf_control%use_ot)) THEN
541 : WRITE (UNIT=output_unit, &
542 : FMT="(/,T3,A,T12,A,T31,A,T39,A,T59,A,T75,A,/,T3,A)") &
543 8905 : "Step", "Update method", "Time", "Convergence", "Total energy", "Change", &
544 17810 : REPEAT("-", 78)
545 : END IF
546 24501 : CALL cp_add_iter_level(logger%iter_info, "QS_SCF")
547 :
548 : ! check for external communicator and if the intermediate energy should be sent
549 98004 : res_val_3(:) = -1.0_dp
550 24501 : description = "[EXT_SCF_ENER_COMM]"
551 24501 : IF (test_for_result(results, description=description)) THEN
552 : CALL get_results(results, description=description, &
553 0 : values=res_val_3, n_entries=i_tmp)
554 0 : CPASSERT(i_tmp == 3)
555 0 : IF (ALL(res_val_3(:) <= 0.0_dp)) THEN
556 : CALL cp_abort(__LOCATION__, &
557 : " Trying to access result ("//TRIM(description)// &
558 0 : ") which is not correctly stored.")
559 : END IF
560 0 : CALL external_comm%set_handle(NINT(res_val_3(1)))
561 : END IF
562 24501 : ext_master_id = NINT(res_val_3(2))
563 24501 : scf_energy_message_tag = NINT(res_val_3(3))
564 :
565 : ! *** outer loop of the scf, can treat other variables,
566 : ! *** such as lagrangian multipliers
567 24501 : scf_env%outer_scf%iter_count = 0
568 24501 : accepted_ot_kp_searches = 0
569 24501 : iter_count = 0
570 24501 : kp_ot_entry_reason = kp_ot_entry_initial
571 24501 : ot_kp_subspace_refresh_count = 0
572 24501 : ot_kp_subspace_refresh_iter_count = 0
573 24501 : ot_kp_subspace_refresh_pending = .FALSE.
574 24501 : total_steps = 0
575 24501 : energy%tot_old = 0.0_dp
576 :
577 : scf_outer_loop: DO
578 :
579 : CALL init_scf_loop(scf_env=scf_env, qs_env=qs_env, &
580 : scf_section=scf_section, &
581 25439 : kp_ot_entry_reason=kp_ot_entry_reason)
582 25439 : kp_ot_entry_reason = kp_ot_entry_reuse
583 25439 : scf_env%adiis_shift = scf_control%diagonalization%adiis_shift
584 :
585 : CALL qs_scf_set_loop_flags(scf_env, diis_step, &
586 25439 : energy_only, just_energy, exit_inner_loop)
587 25439 : IF (ot_kp_subspace_refresh_pending) THEN
588 24 : scf_env%iter_count = ot_kp_subspace_refresh_iter_count
589 24 : ot_kp_subspace_refresh_pending = .FALSE.
590 : END IF
591 :
592 : ! decide whether to switch off dipole correction for convergence purposes
593 25439 : dft_control%surf_dip_correct_switch = dft_control%correct_surf_dip
594 25439 : IF ((dft_control%correct_surf_dip) .AND. (scf_control%outer_scf%have_scf) .AND. &
595 : (scf_env%outer_scf%iter_count > FLOOR(scf_control%outer_scf%max_scf/2.0_dp))) THEN
596 0 : IF (dft_control%switch_surf_dip) THEN
597 0 : dft_control%surf_dip_correct_switch = .FALSE.
598 0 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
599 0 : "SURFACE DIPOLE CORRECTION switched off"
600 : END IF
601 : END IF
602 :
603 235719 : scf_loop: DO
604 :
605 235719 : CALL timeset(routineN//"_inner_loop", handle2)
606 :
607 235719 : IF (.NOT. just_energy) scf_env%iter_count = scf_env%iter_count + 1
608 235719 : iter_count = iter_count + 1
609 235719 : CALL cp_iterate(logger%iter_info, last=.FALSE., iter_nr=iter_count)
610 :
611 235719 : IF (output_unit > 0) CALL m_flush(output_unit)
612 :
613 235719 : total_steps = total_steps + 1
614 235719 : just_energy = energy_only
615 :
616 : CALL qs_ks_update_qs_env(qs_env, just_energy=just_energy, &
617 235719 : calculate_forces=.FALSE.)
618 :
619 235719 : scf_env%raw_map_delta = 0.0_dp
620 235719 : scf_env%raw_map_delta_valid = .FALSE.
621 : ! print 'heavy weight' or relatively expensive quantities
622 235719 : CALL qs_scf_loop_print(qs_env, scf_env, para_env)
623 :
624 235719 : added_mos_auto_grow = .FALSE.
625 235719 : ot_kp_subspace_refresh = .FALSE.
626 : adiis_validation = scf_control%diagonalization%update_method == diag_update_method_adiis .AND. &
627 235719 : scf_env%adiis_check_next
628 235719 : adiis_step_delta = scf_env%step_norm
629 235719 : scf_env%adiis_validated = .FALSE.
630 :
631 235719 : IF (scf_control%diagonalization%update_method == diag_update_method_adiis) THEN
632 : ! ADIIS is a Fock-space SCF method and never uses density mixing. A small
633 : ! accelerated step is checked against the unmodified F[P] map before
634 : ! convergence is accepted. The check detects an artificial subspace zero;
635 : ! it is not a second, equally strict convergence criterion.
636 148 : CPASSERT(ASSOCIATED(scf_env%scf_subspace_buffer))
637 148 : scf_env%scf_subspace_buffer%last_restart = .FALSE.
638 148 : scf_env%scf_subspace_buffer%use_combined_fock = .FALSE.
639 148 : scf_env%scf_subspace_buffer%last_old_fock_weight = 0.0_dp
640 148 : IF (.NOT. adiis_validation .AND. scf_env%iter_count > 1) THEN
641 : ! Only accepted densities whose raw F[P] has just been evaluated enter history.
642 : ! The initial guess and raw-validation trial endpoint are deliberately excluded.
643 : ! A shifted candidate needs a current ADIIS fallback even if CDIIS was
644 : ! active last time: its acceptance is decided later, during diagonalization.
645 : IF (scf_control%diagonalization%adiis_shift > 0.0_dp .OR. &
646 136 : scf_env%scf_subspace_buffer%diis_weight < 1.0_dp .OR. &
647 : scf_env%iter_delta >= scf_control%eps_diis) THEN
648 : CALL qs_scf_subspace_push(scf_env%scf_subspace_buffer, matrix_ks_kp, rho_ao_kp, &
649 110 : energy%total, adiis_pushed)
650 110 : IF (.NOT. adiis_pushed) CPABORT("Failed to append the accepted ADIIS SCF state")
651 110 : CALL qs_scf_subspace_build(scf_env%scf_subspace_buffer)
652 110 : IF (scf_control%diagonalization%adiis_shift > 0.0_dp) THEN
653 50 : CALL qs_scf_subspace_update_shift(scf_env%scf_subspace_buffer, scf_env%adiis_shift)
654 : END IF
655 : END IF
656 : END IF
657 148 : IF (do_kpoints) THEN
658 : CALL qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step, &
659 80 : added_mos_auto_grow=added_mos_auto_grow)
660 : ELSE
661 68 : CALL qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, energy_only)
662 : END IF
663 :
664 148 : IF (adiis_validation) THEN
665 0 : CALL qs_scf_candidate_density_delta(scf_env, rho, para_env, adiis_raw_delta)
666 0 : scf_env%raw_map_delta = adiis_raw_delta
667 0 : scf_env%raw_map_delta_valid = .TRUE.
668 0 : IF (adiis_raw_delta < scf_control%eps_scf) THEN
669 0 : scf_env%iter_delta = adiis_raw_delta
670 0 : scf_env%iter_param = 0.0_dp
671 0 : scf_env%iter_method = "ADIIS/Chk."
672 0 : scf_env%adiis_validated = .TRUE.
673 0 : scf_env%adiis_check_next = .FALSE.
674 : ELSE
675 0 : adiis_stagnated = adiis_step_delta <= adiis_stagnation_ratio*adiis_raw_delta
676 0 : IF (adiis_stagnated) THEN
677 : ! A nearly vanishing ADIIS step can be a fixed point of the current
678 : ! subspace without being a fixed point of the raw SCF map. Discard
679 : ! the stale subspace, retain the current paired P,F[P] state, and
680 : ! accept the already computed raw candidate as the restart step.
681 : ! This avoids duplicate history entries and a second diagonalization.
682 : CALL qs_scf_subspace_restart(scf_env%scf_subspace_buffer, matrix_ks_kp, rho_ao_kp, &
683 0 : energy%total, adiis_restarted)
684 0 : IF (.NOT. adiis_restarted) THEN
685 0 : CPABORT("Failed to restart the ADIIS history from the current SCF state")
686 : END IF
687 0 : scf_env%iter_delta = adiis_raw_delta
688 0 : scf_env%iter_param = 0.0_dp
689 0 : scf_env%iter_method = "ADIIS/Rst."
690 0 : scf_env%adiis_check_next = .FALSE.
691 : ELSE
692 : ! The raw map is not converged, but the accelerated step is not an
693 : ! artificial subspace zero. Preserve the useful history, discard the
694 : ! raw trial, and repeat the normal ADIIS step in this iteration.
695 : CALL qs_scf_subspace_push(scf_env%scf_subspace_buffer, matrix_ks_kp, rho_ao_kp, &
696 0 : energy%total, adiis_pushed)
697 0 : IF (.NOT. adiis_pushed) CPABORT("Failed to append the accepted ADIIS SCF state")
698 0 : CALL qs_scf_subspace_build(scf_env%scf_subspace_buffer)
699 0 : IF (do_kpoints) THEN
700 : CALL qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step, &
701 0 : added_mos_auto_grow=added_mos_auto_grow)
702 : ELSE
703 0 : CALL qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, energy_only)
704 : END IF
705 0 : CALL qs_scf_candidate_density_delta(scf_env, rho, para_env, scf_env%iter_delta)
706 0 : scf_env%adiis_validated = scf_env%iter_delta < scf_control%eps_scf
707 0 : IF (scf_env%adiis_validated .AND. .NOT. diis_step) THEN
708 0 : scf_env%iter_param = 0.0_dp
709 0 : scf_env%iter_method = "ADIIS/Chk."
710 : END IF
711 : ! A successful sanity check must not recursively demand another raw
712 : ! endpoint check. A later, independently small step can trigger one.
713 0 : scf_env%adiis_check_next = .FALSE.
714 : END IF
715 : END IF
716 : END IF
717 :
718 : IF (.NOT. adiis_validation) THEN
719 148 : CALL qs_scf_candidate_density_delta(scf_env, rho, para_env, scf_env%iter_delta)
720 : scf_env%adiis_check_next = .NOT. diis_step .AND. &
721 296 : scf_env%iter_delta < scf_control%eps_scf
722 : END IF
723 : ELSE
724 : ! Keep the existing density-mixing SCF path unchanged.
725 235571 : IF (do_kpoints) THEN
726 : ! kpoints
727 36660 : IF (dft_control%hairy_probes .EQV. .TRUE.) THEN
728 0 : scf_control%smear%do_smear = .FALSE.
729 : CALL qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step, dft_control%probe, &
730 0 : energy_only=energy_only)
731 : ELSE
732 : CALL qs_scf_new_mos_kp( &
733 : qs_env, scf_env, scf_control, diis_step, &
734 : ot_kp_subspace_refresh=ot_kp_subspace_refresh, &
735 : allow_ot_kp_subspace_refresh= &
736 : ot_kp_subspace_refresh_count == 0, &
737 : allow_ot_kp_exit_refresh= &
738 : ot_kp_subspace_refresh_count > 0 .AND. &
739 : ot_kp_subspace_refresh_count < max_ot_kp_subspace_refreshes, &
740 : accepted_ot_kp_searches=accepted_ot_kp_searches, &
741 36660 : added_mos_auto_grow=added_mos_auto_grow, energy_only=energy_only)
742 : END IF
743 : ELSE
744 : ! Gamma points only
745 198911 : IF (dft_control%hairy_probes .EQV. .TRUE.) THEN
746 14 : scf_control%smear%do_smear = .FALSE.
747 : CALL qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, energy_only, &
748 14 : dft_control%probe)
749 : ELSE
750 198897 : CALL qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, energy_only)
751 : END IF
752 : END IF
753 : END IF
754 :
755 235719 : IF (added_mos_auto_grow) THEN
756 10 : CALL qs_scf_grow_added_mos_auto_kp(qs_env, scf_env, scf_control, output_unit)
757 10 : kp_ot_entry_reason = kp_ot_entry_resized
758 10 : total_steps = MAX(0, total_steps - 1)
759 10 : iter_count = MAX(0, iter_count - 1)
760 10 : IF (.NOT. just_energy) scf_env%iter_count = MAX(0, scf_env%iter_count - 1)
761 10 : CALL timestop(handle2)
762 10 : CYCLE scf_outer_loop
763 : END IF
764 :
765 235709 : IF (ot_kp_subspace_refresh) THEN
766 24 : IF (output_unit > 0) THEN
767 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
768 12 : "K-point OT: rebuilding physical virtual subspace by full KS diagonalization."
769 : END IF
770 24 : CALL qs_kpoint_state_commit(qs_env, update_occupations=.FALSE.)
771 24 : IF (ASSOCIATED(scf_env%qs_ot_env)) THEN
772 54 : DO ispin = 1, SIZE(scf_env%qs_ot_env)
773 54 : CALL ot_scf_destroy(scf_env%qs_ot_env(ispin))
774 : END DO
775 24 : DEALLOCATE (scf_env%qs_ot_env)
776 24 : NULLIFY (scf_env%qs_ot_env)
777 : END IF
778 24 : IF (ASSOCIATED(kpoints%scf_diis_buffer)) THEN
779 0 : CALL qs_diis_b_clear_kp(kpoints%scf_diis_buffer)
780 : END IF
781 24 : reuse_precond = .FALSE.
782 24 : accepted_ot_kp_searches = 0
783 24 : ot_kp_subspace_refresh_count = ot_kp_subspace_refresh_count + 1
784 24 : total_steps = MAX(0, total_steps - 1)
785 24 : iter_count = MAX(0, iter_count - 1)
786 24 : IF (.NOT. just_energy) scf_env%iter_count = MAX(0, scf_env%iter_count - 1)
787 24 : ot_kp_subspace_refresh_iter_count = scf_env%iter_count
788 24 : ot_kp_subspace_refresh_pending = .TRUE.
789 24 : kp_ot_entry_reason = kp_ot_entry_refresh
790 24 : CALL timestop(handle2)
791 24 : CYCLE scf_outer_loop
792 : END IF
793 :
794 235685 : IF (do_kpoints .AND. scf_env%method == ot_method_nr .AND. &
795 : qs_scf_kp_search_endpoint(scf_env%iter_method)) THEN
796 1440 : accepted_ot_kp_searches = accepted_ot_kp_searches + 1
797 : END IF
798 :
799 : ! Print requested MO information (can be computationally expensive with OT)
800 235685 : CALL qs_scf_write_mos(qs_env, scf_env, final_mos=.FALSE.)
801 :
802 235685 : IF (dft_control%qs_control%xtb_control%do_tblite) THEN
803 25400 : IF (scf_env%method == ot_method_nr) THEN
804 834 : CALL tb_update_charges(qs_env, dft_control, qs_env%tb_tblite, .FALSE., .TRUE.)
805 834 : CALL evaluate_core_matrix_traces(qs_env)
806 : ELSE
807 24566 : CPASSERT(scf_env%mixing_method > 0)
808 24566 : CALL tb_update_charges(qs_env, dft_control, qs_env%tb_tblite, .FALSE., .FALSE.)
809 24566 : CALL evaluate_core_matrix_traces(qs_env, rho_ao_ext=scf_env%p_mix_new)
810 : END IF
811 25400 : CALL tb_get_energy(qs_env, qs_env%tb_tblite, energy)
812 : END IF
813 :
814 235685 : IF (scf_control%diagonalization%update_method == diag_update_method_adiis) THEN
815 148 : CALL qs_scf_commit_density_candidate(scf_env, rho)
816 148 : density_full_step = .TRUE.
817 : ELSE
818 235537 : density_full_step = diis_step .OR. tblite_native_mixer .OR. internal_tblite_density_full_step
819 235537 : CALL qs_scf_density_mixing(scf_env, rho, para_env, density_full_step)
820 : END IF
821 235685 : IF (dft_control%qs_control%xtb_control%do_tblite .AND. &
822 : .NOT. (dft_control%qs_control%do_ls_scf .OR. scf_control%use_ot)) THEN
823 : scf_env%iter_delta = MAX(scf_env%iter_delta, &
824 : tb_scf_mixer_error(dft_control, qs_env%tb_tblite, &
825 24566 : scf_control%eps_scf))
826 : END IF
827 235685 : IF (dft_control%qs_control%dftb .OR. &
828 : (dft_control%qs_control%xtb .AND. .NOT. dft_control%qs_control%xtb_control%do_tblite)) THEN
829 : scf_env%iter_delta = MAX(scf_env%iter_delta, &
830 57094 : charge_mixing_scc_error(scf_env%mixing_store, scf_control%eps_scf))
831 : END IF
832 235685 : IF (tblite_native_mixer) THEN
833 22408 : scf_env%iter_param = dft_control%qs_control%xtb_control%tblite_mixer_damping
834 22408 : scf_env%iter_method = "TBLite/Diag"
835 213277 : ELSE IF (internal_tblite_mixer) THEN
836 30 : scf_env%iter_method = "TBLite/Diag"
837 30 : IF (dft_control%qs_control%dftb) THEN
838 18 : scf_env%iter_param = dft_control%qs_control%dftb_control%tblite_mixer_damping
839 : ELSE
840 12 : scf_env%iter_param = dft_control%qs_control%xtb_control%tblite_mixer_damping
841 : END IF
842 : END IF
843 235685 : scf_env%step_norm = scf_env%iter_delta
844 :
845 235685 : t2 = m_walltime()
846 :
847 : CALL qs_scf_loop_info(scf_env, output_unit, just_energy, t1, t2, energy, &
848 235685 : scf_control%diagonalization%adiis_verbose)
849 :
850 235685 : IF (scf_control%gce%do_gce) THEN
851 64 : CALL qs_scf_gce_info(output_unit, qs_env, just_energy)
852 : END IF
853 :
854 235685 : IF (.NOT. just_energy) energy%tot_old = energy%total
855 :
856 : ! check for external communicator and if the intermediate energy should be sent
857 235685 : IF (scf_energy_message_tag > 0) THEN
858 0 : CALL external_comm%send(energy%total, ext_master_id, scf_energy_message_tag)
859 : END IF
860 :
861 : CALL qs_scf_check_inner_exit(qs_env, scf_env, scf_control, should_stop, just_energy, &
862 235685 : exit_inner_loop, inner_loop_converged, output_unit)
863 :
864 : ! In case we decide to exit we perform few more check to see if this one
865 : ! is really the last SCF step
866 235685 : IF (exit_inner_loop) THEN
867 :
868 25405 : CALL qs_scf_inner_finalize(scf_env, qs_env, density_full_step, output_unit)
869 :
870 : CALL qs_scf_check_outer_exit(qs_env, scf_env, scf_control, should_stop, &
871 25405 : outer_loop_converged, exit_outer_loop)
872 :
873 : ! Let's tag the last SCF cycle so we can print informations only of the last step
874 25405 : IF (exit_outer_loop) CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=iter_count)
875 :
876 : END IF
877 :
878 235685 : IF (do_kpoints) THEN
879 36706 : CALL write_kpoints_restart(rho_ao_kp, kpoints, scf_env, dft_section, particle_set, qs_kind_set)
880 : ELSE
881 198979 : IF (.NOT. dft_control%mtlr_dft_with_perturbation) THEN
882 : ! Write wavefunction restart file
883 197135 : IF (scf_env%method == ot_method_nr) THEN
884 : ! With OT: provide the Kohn-Sham matrix for the calculation of the MO eigenvalues
885 80552 : CALL get_ks_env(ks_env=ks_env, matrix_ks=matrix_ks)
886 : CALL write_mo_set_to_restart(mos, particle_set, dft_section=dft_section, qs_kind_set=qs_kind_set, &
887 80552 : matrix_ks=matrix_ks)
888 : ELSE
889 116583 : CALL write_mo_set_to_restart(mos, particle_set, dft_section=dft_section, qs_kind_set=qs_kind_set)
890 : END IF
891 : END IF
892 : END IF
893 :
894 : ! Exit if we have finished with the SCF inner loop
895 235685 : IF (exit_inner_loop) THEN
896 25405 : CALL timestop(handle2)
897 : EXIT scf_loop
898 : END IF
899 :
900 210280 : IF (.NOT. BTEST(cp_print_key_should_output(logger%iter_info, &
901 : scf_section, "PRINT%ITERATION_INFO/TIME_CUMUL"), cp_p_file)) THEN
902 210280 : t1 = m_walltime()
903 : END IF
904 :
905 : ! mixing methods have the new density matrix in p_mix_new
906 210280 : IF (scf_env%mixing_method > 0) THEN
907 1577012 : DO ic = 1, SIZE(rho_ao_kp, 2)
908 3146403 : DO ispin = 1, dft_control%nspins
909 1569391 : CALL dbcsr_get_info(rho_ao_kp(ispin, ic)%matrix, name=name) ! keep the name
910 3012838 : CALL dbcsr_copy(rho_ao_kp(ispin, ic)%matrix, scf_env%p_mix_new(ispin, ic)%matrix, name=name)
911 : END DO
912 : END DO
913 : END IF
914 :
915 : CALL qs_scf_rho_update(rho, qs_env, scf_env, ks_env, &
916 210280 : mix_rho=scf_env%mixing_method >= gspace_mixing_nr)
917 :
918 210280 : CALL timestop(handle2)
919 :
920 : END DO scf_loop
921 :
922 25405 : IF (.NOT. scf_control%outer_scf%have_scf) EXIT scf_outer_loop
923 :
924 : ! In case we use the OUTER SCF loop let's print some info..
925 : CALL qs_scf_outer_loop_info(output_unit, scf_control, scf_env, &
926 5601 : energy, total_steps, should_stop, outer_loop_converged)
927 :
928 : ! Save MOs to converged MOs if outer_loop_converged and surf_dip_correct_switch is true
929 5601 : IF (exit_outer_loop) THEN
930 4697 : IF ((dft_control%switch_surf_dip) .AND. (outer_loop_converged) .AND. &
931 : (dft_control%surf_dip_correct_switch)) THEN
932 4 : DO ispin = 1, dft_control%nspins
933 4 : CALL reassign_allocated_mos(mos_last_converged(ispin), mos(ispin))
934 : END DO
935 2 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
936 1 : "COPIED mos ---> mos_last_converged"
937 : END IF
938 : END IF
939 :
940 5601 : IF (exit_outer_loop) EXIT scf_outer_loop
941 :
942 : !
943 904 : CALL outer_loop_optimize(scf_env, scf_control)
944 904 : CALL outer_loop_update_qs_env(qs_env, scf_env)
945 25405 : CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.)
946 :
947 : END DO scf_outer_loop
948 :
949 24501 : converged = inner_loop_converged .AND. outer_loop_converged
950 24501 : total_scf_steps = total_steps
951 :
952 24501 : IF (dft_control%qs_control%cdft) THEN
953 : dft_control%qs_control%cdft_control%total_steps = &
954 710 : dft_control%qs_control%cdft_control%total_steps + total_steps
955 : END IF
956 :
957 24501 : IF (.NOT. converged) THEN
958 2344 : IF (scf_control%ignore_convergence_failure .OR. should_stop) THEN
959 2344 : CALL cp_warn(__LOCATION__, "SCF run NOT converged")
960 : ELSE
961 : CALL cp_abort(__LOCATION__, &
962 : "SCF run NOT converged. To continue the calculation "// &
963 0 : "regardless, please set the keyword IGNORE_CONVERGENCE_FAILURE.")
964 : END IF
965 : END IF
966 :
967 : ! Skip Harris functional calculation if ground-state is NOT converged
968 24501 : IF (qs_env%energy_correction) THEN
969 678 : CALL get_qs_env(qs_env, ec_env=ec_env)
970 678 : ec_env%do_skip = .FALSE.
971 678 : IF (ec_env%skip_ec .AND. .NOT. converged) ec_env%do_skip = .TRUE.
972 : END IF
973 :
974 : ! if needed copy mo_coeff dbcsr->fm for later use in post_scf!fm->dbcsr
975 52782 : DO ispin = 1, SIZE(mos) !fm -> dbcsr
976 52782 : IF (mos(ispin)%use_mo_coeff_b) THEN !fm->dbcsr
977 8299 : IF (.NOT. ASSOCIATED(mos(ispin)%mo_coeff_b)) THEN
978 : !fm->dbcsr
979 0 : CPABORT("mo_coeff_b is not allocated")
980 : END IF !fm->dbcsr
981 : CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, & !fm->dbcsr
982 8299 : mos(ispin)%mo_coeff) !fm -> dbcsr
983 : END IF !fm->dbcsr
984 : END DO !fm -> dbcsr
985 :
986 24501 : CALL cp_rm_iter_level(logger%iter_info, level_name="QS_SCF")
987 24501 : CALL timestop(handle)
988 :
989 24501 : END SUBROUTINE scf_env_do_scf
990 :
991 : ! **************************************************************************************************
992 : !> \brief grow an automatic k-point smearing virtual space and rebuild dimensioned state
993 : !> \param qs_env ...
994 : !> \param scf_env ...
995 : !> \param scf_control ...
996 : !> \param output_unit ...
997 : ! **************************************************************************************************
998 10 : SUBROUTINE qs_scf_grow_added_mos_auto_kp(qs_env, scf_env, scf_control, output_unit)
999 :
1000 : TYPE(qs_environment_type), POINTER :: qs_env
1001 : TYPE(qs_scf_env_type), POINTER :: scf_env
1002 : TYPE(scf_control_type), POINTER :: scf_control
1003 : INTEGER, INTENT(IN) :: output_unit
1004 :
1005 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_grow_added_mos_auto_kp'
1006 :
1007 : INTEGER :: base_nmo, common_base_nmo, current_added, current_target_nmo, grow_by, handle, &
1008 : ic, ik, ispin, nmo_mat, nspins, old_nmo, target_nmo
1009 : INTEGER, DIMENSION(2) :: base_nmo_spin, new_added, new_nmo
1010 : LOGICAL :: diag_step, has_unit_metric, need_resize, &
1011 : shared_spin_auto
1012 : REAL(KIND=dp) :: energy_step, flexible_electron_count, &
1013 : maxocc, n_el_f
1014 10 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers, &
1015 10 : old_eigenvalues, old_occupation_numbers
1016 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1017 10 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
1018 : TYPE(cp_fm_type), POINTER :: mo_coeff
1019 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs
1020 10 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
1021 : TYPE(dft_control_type), POINTER :: dft_control
1022 : TYPE(kpoint_env_type), POINTER :: kp
1023 : TYPE(kpoint_type), POINTER :: kpoints
1024 10 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: old_mos
1025 10 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1026 : TYPE(mp_para_env_type), POINTER :: para_env
1027 :
1028 10 : CALL timeset(routineN, handle)
1029 :
1030 10 : NULLIFY (ao_mo_fm_pools, blacs_env, dft_control, eigenvalues, kpoints, matrix_ks, matrix_s, &
1031 10 : mo_coeff, mo_derivs, mos, occupation_numbers, old_eigenvalues, &
1032 10 : old_occupation_numbers, para_env)
1033 :
1034 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env, dft_control=dft_control, &
1035 : kpoints=kpoints, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s, &
1036 : mo_derivs=mo_derivs, mos=mos, para_env=para_env, &
1037 10 : has_unit_metric=has_unit_metric)
1038 :
1039 10 : CPASSERT(ASSOCIATED(mos))
1040 10 : CPASSERT(ASSOCIATED(kpoints))
1041 10 : CPASSERT(ASSOCIATED(dft_control))
1042 10 : CPASSERT(ASSOCIATED(scf_control))
1043 10 : CPASSERT(ASSOCIATED(matrix_ks))
1044 10 : CPASSERT(ASSOCIATED(matrix_s))
1045 :
1046 10 : nspins = dft_control%nspins
1047 10 : CPASSERT(nspins >= 1 .AND. nspins <= SIZE(new_nmo))
1048 30 : new_added = scf_control%added_mos
1049 10 : new_nmo(:) = 0
1050 10 : need_resize = .FALSE.
1051 30 : shared_spin_auto = nspins == 2 .AND. ALL(scf_control%added_mos_auto(1:2))
1052 :
1053 10 : IF (shared_spin_auto) THEN
1054 12 : DO ispin = 1, nspins
1055 8 : current_added = MAX(0, scf_control%added_mos(ispin))
1056 12 : base_nmo_spin(ispin) = MAX(0, mos(ispin)%nmo - current_added)
1057 : END DO
1058 12 : common_base_nmo = MAXVAL(base_nmo_spin(1:nspins))
1059 12 : current_target_nmo = MAXVAL(mos(1:nspins)%nmo)
1060 4 : current_added = MAX(0, current_target_nmo - common_base_nmo)
1061 4 : grow_by = MAX(4, current_added)
1062 12 : target_nmo = MIN(MINVAL(mos(1:nspins)%nao), current_target_nmo + grow_by)
1063 12 : DO ispin = 1, nspins
1064 8 : new_nmo(ispin) = target_nmo
1065 8 : new_added(ispin) = target_nmo - base_nmo_spin(ispin)
1066 12 : need_resize = need_resize .OR. new_nmo(ispin) > mos(ispin)%nmo
1067 : END DO
1068 : ELSE
1069 12 : DO ispin = 1, nspins
1070 6 : old_nmo = mos(ispin)%nmo
1071 6 : new_nmo(ispin) = old_nmo
1072 6 : IF (.NOT. scf_control%added_mos_auto(ispin)) CYCLE
1073 6 : current_added = MAX(0, scf_control%added_mos(ispin))
1074 6 : base_nmo = MAX(0, old_nmo - current_added)
1075 6 : grow_by = MAX(4, current_added)
1076 6 : new_added(ispin) = MIN(mos(ispin)%nao - base_nmo, current_added + grow_by)
1077 6 : new_nmo(ispin) = base_nmo + new_added(ispin)
1078 12 : need_resize = need_resize .OR. new_nmo(ispin) > old_nmo
1079 : END DO
1080 :
1081 6 : IF (nspins == 2) THEN
1082 0 : target_nmo = MAXVAL(new_nmo(1:nspins))
1083 0 : DO ispin = 1, nspins
1084 0 : IF (target_nmo > mos(ispin)%nao) THEN
1085 : CALL cp_abort(__LOCATION__, &
1086 0 : "K-point ADDED_MOS AUTO exhausted the AO basis while matching spin bands.")
1087 : END IF
1088 0 : IF (target_nmo > mos(ispin)%nmo .AND. .NOT. scf_control%added_mos_auto(ispin)) THEN
1089 : CALL cp_abort(__LOCATION__, &
1090 0 : "K-point ADDED_MOS AUTO needs to grow a spin channel with explicit ADDED_MOS.")
1091 : END IF
1092 0 : base_nmo = MAX(0, mos(ispin)%nmo - MAX(0, scf_control%added_mos(ispin)))
1093 0 : new_nmo(ispin) = target_nmo
1094 0 : new_added(ispin) = target_nmo - base_nmo
1095 0 : need_resize = need_resize .OR. new_nmo(ispin) > mos(ispin)%nmo
1096 : END DO
1097 : END IF
1098 : END IF
1099 :
1100 10 : IF (.NOT. need_resize) THEN
1101 : CALL cp_abort(__LOCATION__, &
1102 0 : "K-point ADDED_MOS AUTO cannot grow further although the highest band is still occupied.")
1103 : END IF
1104 :
1105 24 : scf_control%added_mos(1:nspins) = new_added(1:nspins)
1106 10 : IF (output_unit > 0) THEN
1107 5 : IF (nspins == 2) THEN
1108 : WRITE (UNIT=output_unit, FMT="(T2,A,2I5)") &
1109 2 : "K-point ADDED_MOS AUTO: growing virtual-space buffer to:", &
1110 8 : scf_control%added_mos(1:nspins)
1111 : ELSE
1112 : WRITE (UNIT=output_unit, FMT="(T2,A,I0)") &
1113 3 : "K-point ADDED_MOS AUTO: growing virtual-space buffer to: ", &
1114 6 : scf_control%added_mos(1)
1115 : END IF
1116 : END IF
1117 :
1118 54 : ALLOCATE (old_mos(nspins))
1119 24 : DO ispin = 1, nspins
1120 14 : CALL duplicate_mo_set(old_mos(ispin), mos(ispin))
1121 : CALL get_mo_set(old_mos(ispin), maxocc=maxocc, n_el_f=n_el_f, &
1122 14 : flexible_electron_count=flexible_electron_count)
1123 14 : CALL deallocate_mo_set(mos(ispin))
1124 : CALL allocate_mo_set(mo_set=mos(ispin), nao=old_mos(ispin)%nao, nmo=new_nmo(ispin), &
1125 : nelectron=old_mos(ispin)%nelectron, n_el_f=n_el_f, maxocc=maxocc, &
1126 14 : flexible_electron_count=flexible_electron_count)
1127 38 : mos(ispin)%use_mo_coeff_b = old_mos(ispin)%use_mo_coeff_b
1128 : END DO
1129 :
1130 10 : CALL mpools_rebuild_fm_pools(qs_env%mpools, mos=mos, blacs_env=blacs_env, para_env=para_env)
1131 10 : CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
1132 :
1133 24 : DO ispin = 1, nspins
1134 : CALL init_mo_set(mos(ispin), fm_pool=ao_mo_fm_pools(ispin)%pool, &
1135 14 : name="qs_env%mo"//TRIM(ADJUSTL(cp_to_string(ispin))))
1136 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, eigenvalues=eigenvalues, &
1137 14 : occupation_numbers=occupation_numbers)
1138 : CALL get_mo_set(old_mos(ispin), eigenvalues=old_eigenvalues, &
1139 14 : occupation_numbers=old_occupation_numbers)
1140 :
1141 14 : old_nmo = old_mos(ispin)%nmo
1142 14 : CALL cp_fm_init_random(mo_coeff, mos(ispin)%nmo)
1143 14 : IF (has_unit_metric) THEN
1144 0 : CALL make_basis_simple(mo_coeff, mos(ispin)%nmo)
1145 : ELSE
1146 14 : CALL make_basis_sm(mo_coeff, mos(ispin)%nmo, matrix_s(1, 1)%matrix)
1147 : END IF
1148 :
1149 244 : eigenvalues(1:old_nmo) = old_eigenvalues(1:old_nmo)
1150 186 : occupation_numbers(:) = 0.0_dp
1151 244 : occupation_numbers(1:old_nmo) = old_occupation_numbers(1:old_nmo)
1152 14 : IF (mos(ispin)%nmo > old_nmo) THEN
1153 14 : energy_step = MAX(scf_control%smear%electronic_temperature, 1.0E-3_dp)
1154 78 : DO ic = old_nmo + 1, mos(ispin)%nmo
1155 78 : eigenvalues(ic) = eigenvalues(old_nmo) + energy_step*REAL(ic - old_nmo, KIND=dp)
1156 : END DO
1157 : END IF
1158 14 : mos(ispin)%homo = old_mos(ispin)%homo
1159 14 : mos(ispin)%lfomo = old_mos(ispin)%lfomo
1160 14 : mos(ispin)%kTS = old_mos(ispin)%kTS
1161 14 : mos(ispin)%mu = old_mos(ispin)%mu
1162 24 : mos(ispin)%uniform_occupation = old_mos(ispin)%uniform_occupation
1163 : END DO
1164 :
1165 10 : IF (dft_control%restricted) CALL mo_set_restrict(mos)
1166 :
1167 24 : DO ispin = 1, nspins
1168 14 : IF (.NOT. mos(ispin)%use_mo_coeff_b) CYCLE
1169 10 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
1170 10 : CALL dbcsr_init_p(mos(ispin)%mo_coeff_b)
1171 : CALL cp_dbcsr_m_by_n_from_row_template(mos(ispin)%mo_coeff_b, template=matrix_s(1, 1)%matrix, &
1172 10 : n=mos(ispin)%nmo, sym=dbcsr_type_no_symmetry)
1173 24 : CALL copy_fm_to_dbcsr(mo_coeff, mos(ispin)%mo_coeff_b)
1174 : END DO
1175 :
1176 10 : IF (ASSOCIATED(mo_derivs)) THEN
1177 12 : DO ispin = 1, SIZE(mo_derivs)
1178 12 : IF (ASSOCIATED(mo_derivs(ispin)%matrix)) CALL dbcsr_release_p(mo_derivs(ispin)%matrix)
1179 : END DO
1180 6 : DEALLOCATE (mo_derivs)
1181 : NULLIFY (mo_derivs)
1182 6 : CALL set_qs_env(qs_env, mo_derivs=mo_derivs)
1183 : END IF
1184 10 : IF (qs_env%requires_mo_derivs) THEN
1185 0 : nmo_mat = MERGE(1, nspins, dft_control%restricted)
1186 0 : ALLOCATE (mo_derivs(nmo_mat))
1187 0 : DO ispin = 1, nmo_mat
1188 0 : NULLIFY (mo_derivs(ispin)%matrix)
1189 0 : CALL dbcsr_init_p(mo_derivs(ispin)%matrix)
1190 : CALL dbcsr_create(mo_derivs(ispin)%matrix, template=mos(ispin)%mo_coeff_b, &
1191 0 : name="mo_derivs", matrix_type=dbcsr_type_no_symmetry)
1192 : END DO
1193 0 : CALL set_qs_env(qs_env, mo_derivs=mo_derivs)
1194 : END IF
1195 :
1196 22 : DO ik = 1, SIZE(kpoints%kp_env)
1197 12 : kp => kpoints%kp_env(ik)%kpoint_env
1198 12 : IF (ASSOCIATED(kp%mos)) THEN
1199 28 : DO ispin = 1, SIZE(kp%mos, 2)
1200 60 : DO ic = 1, SIZE(kp%mos, 1)
1201 48 : CALL deallocate_mo_set(kp%mos(ic, ispin))
1202 : END DO
1203 : END DO
1204 12 : DEALLOCATE (kp%mos)
1205 : NULLIFY (kp%mos)
1206 : END IF
1207 12 : CALL cp_fm_release(kp%pmat)
1208 22 : CALL cp_fm_release(kp%wmat)
1209 : END DO
1210 10 : CALL mpools_release(kpoints%mpools)
1211 10 : CALL kpoint_initialize_mos(kpoints, mos)
1212 10 : CALL kpoint_initialize_mo_set(kpoints)
1213 :
1214 10 : diag_step = .FALSE.
1215 10 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, .FALSE., diag_step)
1216 10 : IF (dft_control%restricted) CALL qs_kpoint_copy_spin_mos(kpoints, nspins)
1217 : CALL qs_kpoint_state_commit( &
1218 : qs_env, update_occupations=.TRUE., &
1219 : separate_spin_occupations=dft_control%restricted, &
1220 20 : fixed_occupations=.NOT. (dft_control%smear .OR. scf_control%smear%do_smear))
1221 :
1222 10 : IF (ASSOCIATED(kpoints%scf_diis_buffer)) CALL qs_diis_b_clear_kp(kpoints%scf_diis_buffer)
1223 10 : IF (ASSOCIATED(scf_env%qs_ot_env)) THEN
1224 14 : DO ispin = 1, SIZE(scf_env%qs_ot_env)
1225 14 : CALL ot_scf_destroy(scf_env%qs_ot_env(ispin))
1226 : END DO
1227 6 : DEALLOCATE (scf_env%qs_ot_env)
1228 6 : NULLIFY (scf_env%qs_ot_env)
1229 : END IF
1230 10 : reuse_precond = .FALSE.
1231 :
1232 24 : DO ispin = 1, nspins
1233 24 : CALL deallocate_mo_set(old_mos(ispin))
1234 : END DO
1235 10 : DEALLOCATE (old_mos)
1236 :
1237 10 : CALL timestop(handle)
1238 :
1239 10 : END SUBROUTINE qs_scf_grow_added_mos_auto_kp
1240 :
1241 : ! **************************************************************************************************
1242 : !> \brief inits those objects needed if you want to restart the scf with, say
1243 : !> only a new initial guess, or different density functional or ...
1244 : !> this will happen just before the scf loop starts
1245 : !> \param scf_env ...
1246 : !> \param qs_env ...
1247 : !> \param scf_section ...
1248 : !> \param kp_ot_entry_reason reason for constructing a new k-point OT workspace
1249 : !> \par History
1250 : !> 03.2006 created [Joost VandeVondele]
1251 : ! **************************************************************************************************
1252 28275 : SUBROUTINE init_scf_loop(scf_env, qs_env, scf_section, kp_ot_entry_reason)
1253 :
1254 : TYPE(qs_scf_env_type), POINTER :: scf_env
1255 : TYPE(qs_environment_type), POINTER :: qs_env
1256 : TYPE(section_vals_type), POINTER :: scf_section
1257 : INTEGER, INTENT(IN), OPTIONAL :: kp_ot_entry_reason
1258 :
1259 : CHARACTER(LEN=*), PARAMETER :: routineN = 'init_scf_loop'
1260 :
1261 : INTEGER :: entry_reason, handle, ikind, ispin, &
1262 : nkpoint, nmo, nspin_ot
1263 : INTEGER, DIMENSION(2) :: kp_range
1264 : LOGICAL :: do_adiis, do_kpoints, do_rotation, fixed_density_prepared, has_unit_metric, &
1265 : kp_diis_step, kpoint_mos_initialized
1266 : REAL(KIND=dp) :: u_ramping
1267 28275 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
1268 : TYPE(cp_fm_type), POINTER :: mo_coeff
1269 28275 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
1270 28275 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp, matrix_t_kp, &
1271 28275 : rho_ao_kp
1272 : TYPE(dbcsr_type), POINTER :: orthogonality_metric
1273 : TYPE(dft_control_type), POINTER :: dft_control
1274 : TYPE(kpoint_type), POINTER :: kpoints
1275 28275 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1276 28275 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1277 : TYPE(qs_rho_type), POINTER :: rho
1278 : TYPE(scf_control_type), POINTER :: scf_control
1279 :
1280 28275 : CALL timeset(routineN, handle)
1281 :
1282 28275 : NULLIFY (scf_control, matrix_ks_kp, matrix_s, matrix_s_kp, matrix_t_kp, matrix_ks, &
1283 28275 : rho_ao_kp, dft_control, mos, mo_coeff, kpoints, qs_kind_set, rho, wkp)
1284 :
1285 28275 : CPASSERT(ASSOCIATED(scf_env))
1286 28275 : CPASSERT(ASSOCIATED(qs_env))
1287 28275 : entry_reason = kp_ot_entry_initial
1288 28275 : IF (PRESENT(kp_ot_entry_reason)) entry_reason = kp_ot_entry_reason
1289 28275 : fixed_density_prepared = .FALSE.
1290 :
1291 : CALL get_qs_env(qs_env=qs_env, &
1292 : scf_control=scf_control, &
1293 : dft_control=dft_control, &
1294 : do_kpoints=do_kpoints, &
1295 : kpoints=kpoints, &
1296 : mos=mos, &
1297 : rho=rho, &
1298 28275 : qs_kind_set=qs_kind_set)
1299 :
1300 28275 : nkpoint = 1
1301 28275 : kp_range = 0
1302 28275 : nspin_ot = dft_control%nspins
1303 28275 : IF (dft_control%restricted) nspin_ot = 1
1304 28275 : IF (do_kpoints) THEN
1305 3702 : CALL get_kpoint_info(kpoints, nkp=nkpoint, kp_range=kp_range, wkp=wkp)
1306 : END IF
1307 :
1308 : ! if using mo_coeff_b then copy to fm
1309 60903 : DO ispin = 1, SIZE(mos) !fm->dbcsr
1310 60903 : IF (mos(1)%use_mo_coeff_b) THEN !fm->dbcsr
1311 9434 : CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, mos(ispin)%mo_coeff) !fm->dbcsr
1312 : END IF !fm->dbcsr
1313 : END DO !fm->dbcsr
1314 :
1315 : ! this just guarantees that all mo_occupations match the eigenvalues, if smear
1316 60903 : DO ispin = 1, dft_control%nspins
1317 : ! do not reset mo_occupations if the maximum overlap method is in use
1318 60903 : IF (.NOT. scf_control%diagonalization%mom) THEN
1319 : !if the hair probes section is present, this sends hairy_probes to set_mo_occupation subroutine
1320 : !and switches off the standard smearing
1321 32584 : IF (dft_control%hairy_probes .EQV. .TRUE.) THEN
1322 4 : IF (scf_env%outer_scf%iter_count > 0) THEN
1323 0 : scf_control%smear%do_smear = .FALSE.
1324 : CALL set_mo_occupation(mo_set=mos(ispin), &
1325 : smear=scf_control%smear, &
1326 0 : probe=dft_control%probe)
1327 : END IF
1328 : ELSE
1329 32580 : IF (.NOT. scf_control%gce%do_gce) THEN
1330 : CALL set_mo_occupation(mo_set=mos(ispin), &
1331 : smear=scf_control%smear, &
1332 32576 : emit_warnings=.NOT. do_kpoints)
1333 : ELSE
1334 : CALL set_mo_occupation(mo_set=mos(ispin), &
1335 : smear=scf_control%smear, &
1336 : gce=scf_control%gce, &
1337 4 : emit_warnings=.NOT. do_kpoints)
1338 : END IF
1339 : END IF
1340 : END IF
1341 : END DO
1342 :
1343 28275 : do_adiis = scf_control%diagonalization%update_method == diag_update_method_adiis
1344 28275 : IF (do_adiis) THEN
1345 12 : IF (scf_env%method /= general_diag_method_nr) THEN
1346 0 : CPABORT("ADIIS currently requires SCF%DIAGONALIZATION ALGORITHM STANDARD")
1347 : END IF
1348 12 : IF (scf_env%cholesky_method == cholesky_dbcsr) THEN
1349 0 : CPABORT("ADIIS is not yet compatible with CHOLESKY INVERSE_DBCSR")
1350 : END IF
1351 12 : IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb .OR. &
1352 : dft_control%qs_control%semi_empirical) THEN
1353 0 : CPABORT("ADIIS currently requires a Quickstep DFT KS matrix")
1354 : END IF
1355 12 : IF (dft_control%do_admm_mo) THEN
1356 0 : CPABORT("ADIIS is not yet compatible with ADMM-MO")
1357 : END IF
1358 12 : IF (dft_control%sic_method_id == sic_eo) THEN
1359 0 : CPABORT("ADIIS is not yet compatible with explicit-orbital SIC")
1360 : END IF
1361 12 : IF (dft_control%apply_period_efield) THEN
1362 0 : CPABORT("ADIIS is not yet compatible with PERIODIC_EFIELD")
1363 : END IF
1364 12 : IF (dft_control%dft_plus_u) THEN
1365 8 : CPASSERT(ASSOCIATED(qs_kind_set))
1366 24 : DO ikind = 1, SIZE(qs_kind_set)
1367 16 : CALL get_qs_kind(qs_kind_set(ikind), u_ramping=u_ramping)
1368 24 : IF (u_ramping > 0.0_dp) THEN
1369 0 : CPABORT("ADIIS is not yet compatible with DFT+U U_RAMPING")
1370 : END IF
1371 : END DO
1372 : END IF
1373 12 : IF (dft_control%smear .OR. scf_control%smear%do_smear) THEN
1374 0 : CPABORT("ADIIS finite-temperature smearing support is not implemented yet")
1375 : END IF
1376 : IF (dft_control%roks .OR. scf_control%diagonalization%mom .OR. &
1377 12 : dft_control%hairy_probes .OR. scf_control%gce%do_gce) THEN
1378 0 : CPABORT("ADIIS currently supports ordinary RKS/UKS occupations only")
1379 : END IF
1380 12 : IF (scf_control%do_diag_sub) THEN
1381 0 : CPABORT("ADIIS is not yet compatible with DIAG_SUB_SCF")
1382 : END IF
1383 12 : IF (scf_control%diagonalization%max_history < 1) THEN
1384 0 : CPABORT("ADIIS requires SCF%ADIIS%MAX_HISTORY >= 1")
1385 : END IF
1386 12 : IF (.NOT. ASSOCIATED(scf_env%scf_subspace_buffer)) THEN
1387 12 : ALLOCATE (scf_env%scf_subspace_buffer)
1388 : CALL qs_scf_subspace_buffer_create(scf_env%scf_subspace_buffer, &
1389 12 : scf_control%diagonalization%max_history)
1390 0 : ELSE IF (scf_env%scf_subspace_buffer%nbuffer /= scf_control%diagonalization%max_history) THEN
1391 0 : CALL qs_scf_subspace_buffer_release(scf_env%scf_subspace_buffer)
1392 : CALL qs_scf_subspace_buffer_create(scf_env%scf_subspace_buffer, &
1393 0 : scf_control%diagonalization%max_history)
1394 : END IF
1395 12 : CALL qs_scf_subspace_buffer_clear(scf_env%scf_subspace_buffer)
1396 : END IF
1397 :
1398 28275 : SELECT CASE (scf_env%method)
1399 : CASE DEFAULT
1400 :
1401 : CALL cp_abort(__LOCATION__, &
1402 0 : "Unknown SCF method <"//TRIM(cp_to_string(scf_env%method))//"> found. Check the code!")
1403 :
1404 : CASE (filter_matrix_diag_method_nr)
1405 :
1406 10 : IF (.NOT. scf_env%skip_diis) THEN
1407 0 : IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
1408 0 : ALLOCATE (scf_env%scf_diis_buffer)
1409 0 : CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis)
1410 : END IF
1411 0 : CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
1412 : END IF
1413 :
1414 : CASE (general_diag_method_nr, special_diag_method_nr, block_krylov_diag_method_nr, smeagol_method_nr)
1415 20454 : IF (.NOT. scf_env%skip_diis) THEN
1416 19582 : IF (do_kpoints) THEN
1417 3206 : IF (.NOT. ASSOCIATED(kpoints%scf_diis_buffer)) THEN
1418 2556 : ALLOCATE (kpoints%scf_diis_buffer)
1419 2556 : CALL qs_diis_b_create_kp(kpoints%scf_diis_buffer, nbuffer=scf_control%max_diis)
1420 : END IF
1421 3206 : CALL qs_diis_b_clear_kp(kpoints%scf_diis_buffer)
1422 : ELSE
1423 16376 : IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
1424 4492 : ALLOCATE (scf_env%scf_diis_buffer)
1425 4492 : CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis)
1426 : END IF
1427 16376 : CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
1428 : END IF
1429 : END IF
1430 :
1431 : CASE (ot_diag_method_nr)
1432 : IF (scf_control%diagonalization%ot_settings%preconditioner_type == &
1433 24 : ot_precond_full_kinetic .AND. &
1434 : (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb .OR. &
1435 : dft_control%qs_control%semi_empirical)) THEN
1436 : CALL cp_warn(__LOCATION__, &
1437 : "FULL_KINETIC is unavailable for semi-empirical DIAGONALIZATION%OT; "// &
1438 4 : "falling back to FULL_S_INVERSE.")
1439 4 : scf_control%diagonalization%ot_settings%preconditioner_type = ot_precond_s_inverse
1440 4 : scf_control%diagonalization%ot_settings%preconditioner_name = "FULL_S_INVERSE"
1441 : END IF
1442 24 : IF (do_kpoints) THEN
1443 10 : IF (.NOT. scf_env%skip_diis) THEN
1444 6 : IF (.NOT. ASSOCIATED(kpoints%scf_diis_buffer)) THEN
1445 6 : ALLOCATE (kpoints%scf_diis_buffer)
1446 6 : CALL qs_diis_b_create_kp(kpoints%scf_diis_buffer, nbuffer=scf_control%max_diis)
1447 : END IF
1448 6 : CALL qs_diis_b_clear_kp(kpoints%scf_diis_buffer)
1449 : END IF
1450 10 : CALL timestop(handle)
1451 164 : RETURN
1452 : END IF
1453 14 : CALL get_qs_env(qs_env, matrix_ks=matrix_ks, matrix_s=matrix_s)
1454 :
1455 14 : IF (.NOT. scf_env%skip_diis) THEN
1456 12 : IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
1457 12 : ALLOCATE (scf_env%scf_diis_buffer)
1458 12 : CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis)
1459 : END IF
1460 12 : CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
1461 : END IF
1462 :
1463 : ! if an old preconditioner is still around (i.e. outer SCF is active),
1464 : ! remove it if this could be worthwhile
1465 : CALL restart_preconditioner(qs_env, scf_env%ot_preconditioner, &
1466 : scf_control%diagonalization%ot_settings%preconditioner_type, &
1467 14 : dft_control%nspins)
1468 :
1469 : CALL prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, scf_env%ot_preconditioner, &
1470 : scf_control%diagonalization%ot_settings%preconditioner_type, &
1471 : scf_control%diagonalization%ot_settings%precond_solver_type, &
1472 : scf_control%diagonalization%ot_settings%energy_gap, dft_control%nspins, &
1473 : chebyshev_degree=scf_control%diagonalization%ot_settings%chebyshev_degree, &
1474 : low_rank_base=scf_control%diagonalization%ot_settings%low_rank_base, &
1475 : fermi_low_rank_max_rank= &
1476 : scf_control%diagonalization%ot_settings%fermi_low_rank_max_rank, &
1477 : lattice_fft=scf_control%diagonalization%ot_settings%lattice_fft, &
1478 : lattice_fft_local_cells= &
1479 14 : scf_control%diagonalization%ot_settings%lattice_fft_local_cells)
1480 :
1481 : CASE (block_davidson_diag_method_nr)
1482 : ! Preconditioner initialized within the loop, when required
1483 : CASE (ot_method_nr)
1484 7733 : IF (do_kpoints) THEN
1485 154 : IF (.NOT. ASSOCIATED(scf_env%qs_ot_env)) THEN
1486 : CALL get_qs_env(qs_env, matrix_s_kp=matrix_s_kp, &
1487 154 : kinetic_kp=matrix_t_kp)
1488 154 : CPASSERT(ASSOCIATED(matrix_s_kp))
1489 : CALL allocate_qs_ot_envs(scf_env=scf_env, &
1490 : scf_control=scf_control, &
1491 : dft_control=dft_control, &
1492 : scf_section=scf_section, &
1493 : do_kpoints=do_kpoints, &
1494 : nkpoint=nkpoint, &
1495 : nspin_ot=nspin_ot, &
1496 : kp_range=kp_range, &
1497 154 : wkp=wkp)
1498 :
1499 : ! A previous zero-width filling may have changed the occupied rank at each k point.
1500 : ! Restore the fixed OT ranks before validating or reconstructing the active space.
1501 154 : IF (.NOT. scf_env%qs_ot_env(1)%settings%do_ener) THEN
1502 66 : CALL qs_kpoint_set_fixed_occupations(kpoints)
1503 : END IF
1504 : kpoint_mos_initialized = qs_kpoint_mos_initialized( &
1505 154 : kpoints, require_full_space=scf_env%qs_ot_env(1)%settings%do_ener)
1506 106 : SELECT CASE (entry_reason)
1507 : CASE (kp_ot_entry_initial)
1508 106 : IF (kpoint_mos_initialized) THEN
1509 6 : IF (dft_control%restricted) THEN
1510 0 : CALL qs_kpoint_copy_spin_mos(kpoints, dft_control%nspins)
1511 : END IF
1512 : CALL qs_kpoint_state_commit( &
1513 : qs_env, update_occupations=.TRUE., &
1514 : separate_spin_occupations=dft_control%restricted, &
1515 8 : fixed_occupations=.NOT. (dft_control%smear .OR. scf_control%smear%do_smear))
1516 : END IF
1517 : CASE (kp_ot_entry_reuse, kp_ot_entry_resized, kp_ot_entry_refresh)
1518 48 : CPASSERT(kpoint_mos_initialized)
1519 : CASE DEFAULT
1520 154 : CPABORT("Invalid k-point OT workspace entry reason.")
1521 : END SELECT
1522 :
1523 : ! For fixed occupations, reconstruct the occupied subspace directly from the
1524 : ! density-only guess. This preserves its physical content instead of replacing it
1525 : ! with the potentially remote ground-state projector of H[P].
1526 154 : IF (entry_reason == kp_ot_entry_initial .AND. .NOT. kpoint_mos_initialized .AND. &
1527 : .NOT. scf_env%qs_ot_env(1)%settings%do_ener) THEN
1528 46 : CPASSERT(ASSOCIATED(rho))
1529 46 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1530 46 : CPASSERT(ASSOCIATED(rho_ao_kp))
1531 46 : CALL kpoint_operator_store(kpoints, scf_env%scf_work1(1), rho_ao_kp, matrix_s_kp)
1532 46 : CALL qs_kpoint_state_prepare_fixed_density(kpoints)
1533 46 : IF (SIZE(rho_ao_kp, 1) < nspin_ot) THEN
1534 0 : CPASSERT(SIZE(rho_ao_kp, 1) == 1)
1535 0 : CALL qs_kpoint_copy_spin_mos(kpoints, nspin_ot)
1536 : END IF
1537 46 : IF (dft_control%restricted) THEN
1538 2 : CALL qs_kpoint_copy_spin_mos(kpoints, dft_control%nspins)
1539 : END IF
1540 : CALL qs_kpoint_state_commit( &
1541 : qs_env, update_occupations=.TRUE., &
1542 46 : separate_spin_occupations=dft_control%restricted, fixed_occupations=.TRUE.)
1543 :
1544 : ! The reconstructed projector replaces the density-only input. Rebuild H and E
1545 : ! so the first minimizer call sees one coherent physical electronic state.
1546 46 : CALL qs_ks_update_qs_env(qs_env, just_energy=.FALSE., calculate_forces=.FALSE.)
1547 46 : fixed_density_prepared = .TRUE.
1548 :
1549 : ELSE
1550 : ! Build the physical Hamiltonian from the committed state or a Mermin
1551 : ! density-only guess. Ordinary workspace re-entry only needs this rebuild.
1552 108 : CALL qs_ks_update_qs_env(qs_env, just_energy=.FALSE., calculate_forces=.FALSE.)
1553 108 : CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks_kp)
1554 108 : CPASSERT(ASSOCIATED(matrix_ks_kp))
1555 :
1556 : ! Mermin initialization and deliberate REF refreshes require a full orbital
1557 : ! space, so use the physical Hamiltonian replacement for those cases only.
1558 108 : IF ((entry_reason == kp_ot_entry_initial .AND. .NOT. kpoint_mos_initialized) .OR. &
1559 : entry_reason == kp_ot_entry_refresh) THEN
1560 78 : kp_diis_step = .FALSE.
1561 : CALL do_general_diag_kp(matrix_ks_kp, matrix_s_kp, kpoints, scf_env, &
1562 78 : scf_control, .FALSE., kp_diis_step)
1563 78 : IF (SIZE(matrix_ks_kp, 1) < nspin_ot) THEN
1564 0 : CPASSERT(SIZE(matrix_ks_kp, 1) == 1)
1565 0 : CALL qs_kpoint_copy_spin_mos(kpoints, nspin_ot)
1566 : END IF
1567 78 : IF (dft_control%restricted) THEN
1568 4 : CALL qs_kpoint_copy_spin_mos(kpoints, dft_control%nspins)
1569 : END IF
1570 : CALL qs_kpoint_state_commit( &
1571 : qs_env, update_occupations=.TRUE., &
1572 : separate_spin_occupations=dft_control%restricted, &
1573 156 : fixed_occupations=.NOT. (dft_control%smear .OR. scf_control%smear%do_smear))
1574 :
1575 : ! The accepted projector replaced the density-only input. Rebuild H and E so
1576 : ! the first minimizer call sees one coherent physical electronic state.
1577 : CALL qs_ks_update_qs_env(qs_env, just_energy=.FALSE., &
1578 78 : calculate_forces=.FALSE.)
1579 : END IF
1580 : END IF
1581 :
1582 154 : CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks_kp)
1583 154 : CPASSERT(ASSOCIATED(matrix_ks_kp))
1584 :
1585 : ! Redistribute the physical operators to their owning k-point groups. OT must not
1586 : ! call the collective real-space transform independently with group-local k values.
1587 154 : CALL kpoint_operator_store(kpoints, scf_env%scf_work1(1), matrix_ks_kp, matrix_s_kp, matrix_t_kp)
1588 :
1589 : ! Natural orbitals define the fixed occupied projector but not its internal gauge.
1590 : ! Canonicalize that subspace with the matching physical Hamiltonian and populate
1591 : ! the MO energy labels before constructing H-dependent OT preconditioners.
1592 154 : IF (fixed_density_prepared) CALL qs_kpoint_state_canonicalize_fixed(kpoints)
1593 :
1594 : CALL allocate_qs_ot_kpoint_state(qs_env=qs_env, &
1595 : scf_env=scf_env, &
1596 : dft_control=dft_control, &
1597 : kpoints=kpoints, &
1598 : matrix_ks_kp=matrix_ks_kp, &
1599 : matrix_s_kp=matrix_s_kp, &
1600 : matrix_t_kp=matrix_t_kp, &
1601 154 : nspin_ot=nspin_ot)
1602 : END IF
1603 :
1604 154 : CALL timestop(handle)
1605 154 : RETURN
1606 : END IF
1607 :
1608 : CALL get_qs_env(qs_env, &
1609 : has_unit_metric=has_unit_metric, &
1610 : matrix_s=matrix_s, &
1611 7579 : matrix_ks=matrix_ks)
1612 :
1613 : ! reortho the wavefunctions if we are having an outer scf and
1614 : ! this is not the first iteration
1615 : ! this is useful to avoid the build-up of numerical noise
1616 : ! however, we can not play this trick if restricted (don't mix non-equivalent orbs)
1617 7579 : IF (scf_control%do_outer_scf_reortho) THEN
1618 7005 : IF (scf_control%outer_scf%have_scf .AND. .NOT. dft_control%restricted) THEN
1619 4755 : IF (scf_env%outer_scf%iter_count > 0) THEN
1620 1909 : DO ispin = 1, dft_control%nspins
1621 1043 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1622 1909 : IF (has_unit_metric) THEN
1623 108 : CALL make_basis_simple(mo_coeff, nmo)
1624 : ELSE
1625 935 : CALL make_basis_sm(mo_coeff, nmo, matrix_s(1)%matrix)
1626 : END IF
1627 : END DO
1628 : END IF
1629 : END IF
1630 : ELSE
1631 : ! dont need any dirty trick for the numerically stable irac algorithm.
1632 : END IF
1633 :
1634 7579 : IF (.NOT. ASSOCIATED(scf_env%qs_ot_env)) THEN
1635 :
1636 : CALL allocate_qs_ot_envs(scf_env=scf_env, &
1637 : scf_control=scf_control, &
1638 : dft_control=dft_control, &
1639 : scf_section=scf_section, &
1640 : do_kpoints=do_kpoints, &
1641 : nkpoint=nkpoint, &
1642 : nspin_ot=nspin_ot, &
1643 : kp_range=kp_range, &
1644 7579 : wkp=wkp)
1645 :
1646 7579 : IF (scf_env%qs_ot_env(1)%settings%preconditioner_type == &
1647 : ot_precond_full_all_covariant) THEN
1648 8 : IF (.NOT. do_kpoints) THEN
1649 16 : DO ispin = 1, SIZE(mos)
1650 16 : IF (.NOT. mos(ispin)%uniform_occupation) THEN
1651 : CALL cp_abort(__LOCATION__, &
1652 0 : 'PRECONDITIONER FULL_ALL_COVARIANT currently requires uniform occupations.')
1653 : END IF
1654 : END DO
1655 : END IF
1656 : END IF
1657 :
1658 : ! might need the KS matrix to init properly
1659 : CALL qs_ks_update_qs_env(qs_env, just_energy=.FALSE., &
1660 7579 : calculate_forces=.FALSE.)
1661 :
1662 : ! if an old preconditioner is still around (i.e. outer SCF is active),
1663 : ! remove it if this could be worthwhile
1664 7579 : IF (.NOT. reuse_precond) THEN
1665 : CALL restart_preconditioner(qs_env, scf_env%ot_preconditioner, &
1666 : scf_env%qs_ot_env(1)%settings%preconditioner_type, &
1667 7579 : dft_control%nspins)
1668 : END IF
1669 :
1670 : !
1671 : ! preconditioning still needs to be done correctly with has_unit_metric
1672 : ! notice that a big part of the preconditioning (S^-1) is fine anyhow
1673 : !
1674 7579 : IF (has_unit_metric) THEN
1675 1154 : NULLIFY (orthogonality_metric)
1676 : ELSE
1677 6425 : orthogonality_metric => matrix_s(1)%matrix
1678 : END IF
1679 :
1680 7579 : IF (.NOT. reuse_precond) THEN
1681 : CALL prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, scf_env%ot_preconditioner, &
1682 : scf_env%qs_ot_env(1)%settings%preconditioner_type, &
1683 : scf_env%qs_ot_env(1)%settings%precond_solver_type, &
1684 : scf_env%qs_ot_env(1)%settings%energy_gap, dft_control%nspins, &
1685 : has_unit_metric=has_unit_metric, &
1686 : chol_type=scf_env%qs_ot_env(1)%settings%cholesky_type, &
1687 : chebyshev_degree=scf_env%qs_ot_env(1)%settings%chebyshev_degree, &
1688 : low_rank_base=scf_env%qs_ot_env(1)%settings%low_rank_base, &
1689 : fermi_low_rank_max_rank= &
1690 : scf_env%qs_ot_env(1)%settings%fermi_low_rank_max_rank, &
1691 : lattice_fft=scf_env%qs_ot_env(1)%settings%lattice_fft, &
1692 : lattice_fft_local_cells= &
1693 7579 : scf_env%qs_ot_env(1)%settings%lattice_fft_local_cells)
1694 : END IF
1695 7579 : IF (reuse_precond) reuse_precond = .FALSE.
1696 :
1697 : CALL ot_scf_init(mo_array=mos, matrix_s=orthogonality_metric, &
1698 : broyden_adaptive_sigma=qs_env%broyden_adaptive_sigma, &
1699 7579 : qs_ot_env=scf_env%qs_ot_env, matrix_ks=matrix_ks(1)%matrix)
1700 :
1701 13127 : SELECT CASE (scf_env%qs_ot_env(1)%settings%preconditioner_type)
1702 : CASE (ot_precond_none)
1703 : CASE (ot_precond_full_all, ot_precond_fermi_low_rank, ot_precond_full_all_covariant, &
1704 : ot_precond_full_single_inverse)
1705 12321 : DO ispin = 1, SIZE(scf_env%qs_ot_env)
1706 : CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin), &
1707 12321 : scf_env%ot_preconditioner(ispin)%preconditioner)
1708 : END DO
1709 : CASE (ot_precond_s_inverse, ot_precond_full_single)
1710 160 : DO ispin = 1, SIZE(scf_env%qs_ot_env)
1711 : CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin), &
1712 160 : scf_env%ot_preconditioner(1)%preconditioner)
1713 : END DO
1714 : CASE DEFAULT
1715 9062 : DO ispin = 1, SIZE(scf_env%qs_ot_env)
1716 : CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin), &
1717 2720 : scf_env%ot_preconditioner(1)%preconditioner)
1718 : END DO
1719 : END SELECT
1720 : END IF
1721 :
1722 : ! if we have non-uniform occupations we should be using rotation
1723 7579 : do_rotation = scf_env%qs_ot_env(1)%settings%do_rotation
1724 45046 : DO ispin = 1, SIZE(mos)
1725 16771 : IF (.NOT. mos(ispin)%uniform_occupation) THEN
1726 0 : CPASSERT(do_rotation)
1727 : END IF
1728 : END DO
1729 : END SELECT
1730 :
1731 : ! another safety check
1732 28111 : IF (dft_control%low_spin_roks) THEN
1733 24 : CPASSERT(scf_env%method == ot_method_nr)
1734 24 : do_rotation = scf_env%qs_ot_env(1)%settings%do_rotation
1735 24 : CPASSERT(do_rotation)
1736 : END IF
1737 :
1738 28111 : CALL timestop(handle)
1739 :
1740 28275 : END SUBROUTINE init_scf_loop
1741 :
1742 : ! **************************************************************************************************
1743 : !> \brief allocate OT environments and label their spin/k-point channel identity
1744 : !> \param scf_env ...
1745 : !> \param scf_control ...
1746 : !> \param dft_control ...
1747 : !> \param scf_section ...
1748 : !> \param do_kpoints ...
1749 : !> \param nkpoint ...
1750 : !> \param nspin_ot ...
1751 : !> \param kp_range ...
1752 : !> \param wkp ...
1753 : ! **************************************************************************************************
1754 7733 : SUBROUTINE allocate_qs_ot_envs(scf_env, scf_control, dft_control, scf_section, &
1755 : do_kpoints, nkpoint, nspin_ot, kp_range, wkp)
1756 : TYPE(qs_scf_env_type), POINTER :: scf_env
1757 : TYPE(scf_control_type), POINTER :: scf_control
1758 : TYPE(dft_control_type), POINTER :: dft_control
1759 : TYPE(section_vals_type), POINTER :: scf_section
1760 : LOGICAL, INTENT(IN) :: do_kpoints
1761 : INTEGER, INTENT(IN) :: nkpoint, nspin_ot
1762 : INTEGER, DIMENSION(2), INTENT(IN) :: kp_range
1763 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
1764 :
1765 : CHARACTER(len=*), PARAMETER :: kpoint_precond_error = &
1766 : "Complex K-point OT supports PRECONDITIONER NONE, FERMI_LOW_RANK, "// &
1767 : "FULL_S_INVERSE, FULL_KINETIC, FULL_SINGLE, FULL_SINGLE_INVERSE, FULL_ALL, or "// &
1768 : "FULL_ALL_COVARIANT.", kpoint_solver_error = &
1769 : "Complex K-point FERMI_LOW_RANK, FULL_SINGLE, FULL_ALL, and FULL_ALL_COVARIANT "// &
1770 : "support PRECOND_SOLVER DEFAULT; FULL_S_INVERSE, FULL_KINETIC, and "// &
1771 : "FULL_SINGLE_INVERSE support DEFAULT or INVERSE_CHOLESKY."
1772 :
1773 : INTEGER :: ikpoint, ispin, local_kpoint, &
1774 : number_of_OT_envs, ot_channel
1775 : LOGICAL :: do_rotation, do_smear, is_full_all
1776 : REAL(KIND=dp) :: kpoint_weight
1777 :
1778 7733 : CPASSERT(.NOT. ASSOCIATED(scf_env%qs_ot_env))
1779 :
1780 : ! Restricted calculations require just one set of OT orbitals per k-point.
1781 : number_of_OT_envs = qs_ot_number_of_channels(dft_control%nspins, &
1782 : nkpoint=nkpoint, &
1783 7733 : restricted=dft_control%restricted)
1784 :
1785 1296538 : ALLOCATE (scf_env%qs_ot_env(number_of_OT_envs))
1786 :
1787 : ! XXX Joost XXX should disentangle reading input from this part
1788 7733 : IF (scf_env%outer_scf%iter_count > 0) THEN
1789 904 : IF (scf_env%iter_delta < scf_control%eps_diis) THEN
1790 6 : scf_env%qs_ot_env(1)%settings%ot_state = 1
1791 : END IF
1792 : END IF
1793 :
1794 7733 : CALL ot_scf_read_input(scf_env%qs_ot_env, scf_section, do_kpoints)
1795 :
1796 7733 : IF (dft_control%restricted .AND. .NOT. do_kpoints .AND. &
1797 : scf_env%qs_ot_env(1)%settings%preconditioner_type == &
1798 : ot_precond_full_all_covariant) THEN
1799 0 : CPABORT('Gamma-point PRECONDITIONER FULL_ALL_COVARIANT does not currently support ROKS.')
1800 : END IF
1801 :
1802 7733 : IF (do_kpoints) THEN
1803 154 : do_smear = dft_control%smear .OR. scf_control%smear%do_smear
1804 : IF (do_smear) THEN
1805 88 : IF (.NOT. scf_env%qs_ot_env(1)%settings%do_ener) THEN
1806 0 : CPABORT("K-point OT smearing requires OT%ENERGIES and OT%ROTATION.")
1807 : END IF
1808 88 : IF (.NOT. scf_env%qs_ot_env(1)%settings%do_rotation) THEN
1809 0 : CPABORT("K-point OT smearing requires OT%ROTATION.")
1810 : END IF
1811 88 : SELECT CASE (scf_control%smear%method)
1812 : CASE (smear_fermi_dirac, smear_gaussian, smear_mp, smear_mv)
1813 0 : CONTINUE
1814 : CASE DEFAULT
1815 88 : CPABORT("K-point Mermin OT does not support the selected smearing method.")
1816 : END SELECT
1817 88 : IF (scf_env%qs_ot_env(1)%settings%ot_algorithm /= "TOD" .AND. &
1818 : scf_env%qs_ot_env(1)%settings%ot_algorithm /= "REF") THEN
1819 0 : CPABORT("K-point Mermin OT currently supports OT%ALGORITHM STRICT or IRAC.")
1820 : END IF
1821 : IF (scf_env%qs_ot_env(1)%settings%OT_METHOD /= "CG" .AND. &
1822 : scf_env%qs_ot_env(1)%settings%OT_METHOD /= "SD" .AND. &
1823 : scf_env%qs_ot_env(1)%settings%OT_METHOD /= "DIIS" .AND. &
1824 88 : scf_env%qs_ot_env(1)%settings%OT_METHOD /= "BROY" .AND. &
1825 : scf_env%qs_ot_env(1)%settings%OT_METHOD /= "LBFG") THEN
1826 0 : CPABORT("K-point Mermin OT supports OT%MINIMIZER CG, SD, DIIS, BROYDEN, or LBFGS.")
1827 : END IF
1828 66 : ELSE IF (scf_env%qs_ot_env(1)%settings%do_ener) THEN
1829 0 : CPABORT("OT%ENERGIES requires smearing in the complex K-point path.")
1830 : END IF
1831 154 : IF (scf_env%qs_ot_env(1)%settings%ot_algorithm /= "TOD" .AND. &
1832 : scf_env%qs_ot_env(1)%settings%ot_algorithm /= "REF") THEN
1833 0 : CPABORT("K-point OT currently supports OT%ALGORITHM STRICT or IRAC.")
1834 : END IF
1835 : IF (scf_env%qs_ot_env(1)%settings%OT_METHOD /= "CG" .AND. &
1836 : scf_env%qs_ot_env(1)%settings%OT_METHOD /= "SD" .AND. &
1837 : scf_env%qs_ot_env(1)%settings%OT_METHOD /= "DIIS" .AND. &
1838 154 : scf_env%qs_ot_env(1)%settings%OT_METHOD /= "BROY" .AND. &
1839 : scf_env%qs_ot_env(1)%settings%OT_METHOD /= "LBFG") THEN
1840 0 : CPABORT("K-point OT currently supports OT%MINIMIZER CG, SD, DIIS, BROYDEN, or LBFGS.")
1841 : END IF
1842 154 : IF (.NOT. qs_ot_kpoint_preconditioner_supported( &
1843 : scf_env%qs_ot_env(1)%settings%preconditioner_type, .FALSE.)) THEN
1844 0 : CPABORT(kpoint_precond_error)
1845 : END IF
1846 154 : IF (scf_env%qs_ot_env(1)%settings%preconditioner_type == ot_precond_full_kinetic .AND. &
1847 : (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
1848 : dft_control%qs_control%xtb)) THEN
1849 0 : CPABORT("FULL_KINETIC is unavailable for semi-empirical methods")
1850 : END IF
1851 154 : IF (.NOT. qs_ot_kpoint_preconditioner_solver_supported( &
1852 : scf_env%qs_ot_env(1)%settings%preconditioner_type, &
1853 : scf_env%qs_ot_env(1)%settings%precond_solver_type, .FALSE.)) THEN
1854 0 : CPABORT(kpoint_solver_error)
1855 : END IF
1856 154 : IF (scf_env%qs_ot_env(1)%settings%occupation_preconditioner .AND. .NOT. do_smear) THEN
1857 0 : CPABORT("K-point OCCUPATION_PRECONDITIONER requires smearing.")
1858 : END IF
1859 : END IF
1860 :
1861 7733 : IF (scf_env%outer_scf%iter_count > 0) THEN
1862 904 : IF (scf_env%qs_ot_env(1)%settings%ot_state == 1) THEN
1863 : scf_control%max_scf = MAX(scf_env%qs_ot_env(1)%settings%max_scf_diis, &
1864 6 : scf_control%max_scf)
1865 : END IF
1866 : END IF
1867 :
1868 : ! Keep a note that we are restricted.
1869 7733 : IF (dft_control%restricted) THEN
1870 208 : scf_env%qs_ot_env(:)%restricted = .TRUE.
1871 : ! requires rotation
1872 102 : IF (.NOT. scf_env%qs_ot_env(1)%settings%do_rotation) THEN
1873 : CALL cp_abort(__LOCATION__, &
1874 : "Restricted calculation with OT requires orbital rotation. Please "// &
1875 0 : "activate the OT%ROTATION keyword!")
1876 : END IF
1877 : ELSE
1878 16967 : scf_env%qs_ot_env(:)%restricted = .FALSE.
1879 : END IF
1880 :
1881 : ! This will rotate the MOs to be eigen states, which is not compatible with rotation.
1882 : ! e.g. mo_derivs here do not yet include potentially different occupation numbers.
1883 7733 : do_rotation = scf_env%qs_ot_env(1)%settings%do_rotation
1884 : ! Only full all needs rotation.
1885 7733 : is_full_all = scf_env%qs_ot_env(1)%settings%preconditioner_type == ot_precond_full_all
1886 7733 : IF (do_rotation .AND. is_full_all .AND. .NOT. do_kpoints) THEN
1887 0 : CPABORT('PRECONDITIONER FULL_ALL is not compatible with ROTATION.')
1888 : END IF
1889 :
1890 15618 : DO ikpoint = 1, nkpoint
1891 7885 : local_kpoint = 0
1892 7885 : IF (do_kpoints) THEN
1893 306 : IF (ikpoint >= kp_range(1) .AND. ikpoint <= kp_range(2)) THEN
1894 249 : local_kpoint = ikpoint - kp_range(1) + 1
1895 : END IF
1896 : END IF
1897 25060 : DO ispin = 1, nspin_ot
1898 9442 : ot_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot)
1899 17327 : IF (do_kpoints) THEN
1900 342 : kpoint_weight = 1.0_dp
1901 342 : IF (ASSOCIATED(wkp)) kpoint_weight = wkp(ikpoint)
1902 : CALL qs_ot_set_context(scf_env%qs_ot_env(ot_channel), &
1903 : spin_index=ispin, &
1904 : kpoint_index=ikpoint, &
1905 : local_kpoint_index=local_kpoint, &
1906 342 : kpoint_weight=kpoint_weight)
1907 : ELSE
1908 9100 : CALL qs_ot_set_context(scf_env%qs_ot_env(ot_channel), spin_index=ispin)
1909 : END IF
1910 : END DO
1911 : END DO
1912 :
1913 : CALL qs_ot_check_channel_context(scf_env%qs_ot_env, dft_control%nspins, &
1914 : nkpoint=nkpoint, &
1915 : restricted=dft_control%restricted, &
1916 : require_kpoint=do_kpoints, &
1917 : kp_range=kp_range, &
1918 7733 : wkp=wkp)
1919 :
1920 7733 : END SUBROUTINE allocate_qs_ot_envs
1921 :
1922 : ! **************************************************************************************************
1923 : !> \brief allocate local k-point OT minimizer state for labelled spin/k-point channels
1924 : !> \param qs_env ...
1925 : !> \param scf_env ...
1926 : !> \param dft_control ...
1927 : !> \param kpoints ...
1928 : !> \param matrix_ks_kp ...
1929 : !> \param matrix_s_kp ...
1930 : !> \param matrix_t_kp ...
1931 : !> \param nspin_ot ...
1932 : ! **************************************************************************************************
1933 154 : SUBROUTINE allocate_qs_ot_kpoint_state(qs_env, scf_env, dft_control, &
1934 : kpoints, matrix_ks_kp, matrix_s_kp, matrix_t_kp, nspin_ot)
1935 : TYPE(qs_environment_type), POINTER :: qs_env
1936 : TYPE(qs_scf_env_type), POINTER :: scf_env
1937 : TYPE(dft_control_type), POINTER :: dft_control
1938 : TYPE(kpoint_type), POINTER :: kpoints
1939 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp, matrix_t_kp
1940 : INTEGER, INTENT(IN) :: nspin_ot
1941 :
1942 : CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_qs_ot_kpoint_state'
1943 :
1944 : INTEGER :: energy_dimension, energy_spin, &
1945 : energy_start, handle, ikpoint, &
1946 : ikpoint_local, ispin, nao, nmo, nocc, &
1947 : ot_channel
1948 : LOGICAL :: use_real_wfn
1949 154 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1950 : TYPE(cp_fm_struct_type), POINTER :: active_mo_struct
1951 : TYPE(cp_fm_type) :: active_mo_coeff, active_mo_coeff_im
1952 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_im
1953 : TYPE(dbcsr_type), POINTER :: matrix_sk_im, matrix_sk_re
1954 : TYPE(kpoint_env_type), POINTER :: kp
1955 : TYPE(mo_set_type), POINTER :: mo_set, mo_set_im
1956 :
1957 154 : CALL timeset(routineN, handle)
1958 :
1959 154 : NULLIFY (active_mo_struct, eigenvalues, kp, matrix_sk_im, matrix_sk_re, &
1960 154 : mo_coeff, mo_coeff_im, mo_set, mo_set_im)
1961 :
1962 154 : CPASSERT(ASSOCIATED(qs_env))
1963 154 : CPASSERT(ASSOCIATED(scf_env%qs_ot_env))
1964 154 : CPASSERT(ASSOCIATED(kpoints))
1965 154 : CPASSERT(ASSOCIATED(kpoints%kp_env))
1966 154 : CPASSERT(ASSOCIATED(matrix_ks_kp))
1967 154 : CPASSERT(ASSOCIATED(matrix_s_kp))
1968 154 : CPASSERT(ASSOCIATED(matrix_s_kp(1, 1)%matrix))
1969 154 : CPASSERT(.NOT. dft_control%restricted .OR. nspin_ot == 1)
1970 :
1971 154 : CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn)
1972 154 : IF (use_real_wfn) THEN
1973 0 : CPABORT("K-point OT currently requires complex K-point wavefunctions.")
1974 : END IF
1975 :
1976 403 : DO ikpoint_local = 1, SIZE(kpoints%kp_env)
1977 249 : kp => kpoints%kp_env(ikpoint_local)%kpoint_env
1978 249 : CPASSERT(ASSOCIATED(kp))
1979 249 : CPASSERT(ASSOCIATED(kp%mos))
1980 249 : CPASSERT(ASSOCIATED(kp%ot_smat))
1981 249 : CPASSERT(nspin_ot <= SIZE(kp%mos, 2))
1982 249 : ikpoint = kp%nkpoint
1983 249 : CPASSERT(ikpoint >= 1)
1984 249 : CPASSERT(SIZE(kp%ot_smat) >= 2)
1985 : CALL kpoint_operator_get_local(matrix_s_kp, kpoints, kp, 1, &
1986 : kp%ot_smat(1), kp%ot_smat(2), &
1987 249 : matrix_sk_re, matrix_sk_im)
1988 528 : DO ispin = 1, nspin_ot
1989 279 : ot_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot)
1990 279 : CPASSERT(.NOT. scf_env%qs_ot_env(ot_channel)%state_allocated)
1991 :
1992 279 : NULLIFY (mo_coeff, mo_coeff_im, mo_set_im)
1993 279 : mo_set => kp%mos(1, ispin)
1994 : CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, eigenvalues=eigenvalues, &
1995 279 : homo=nocc, nao=nao, nmo=nmo)
1996 279 : CPASSERT(ASSOCIATED(mo_coeff))
1997 279 : CPASSERT(ASSOCIATED(eigenvalues))
1998 279 : CPASSERT(nao > 0)
1999 279 : CPASSERT(nmo > 0)
2000 279 : IF (scf_env%qs_ot_env(ot_channel)%settings%do_ener) nocc = nmo
2001 279 : IF (nocc < 1 .OR. nocc > nmo) THEN
2002 : CALL cp_abort(__LOCATION__, &
2003 0 : "K-point OT requires at least one occupied orbital in every spin channel.")
2004 : END IF
2005 :
2006 : CALL cp_fm_struct_create(active_mo_struct, template_fmstruct=mo_coeff%matrix_struct, &
2007 279 : ncol_global=nocc)
2008 279 : CALL cp_fm_create(active_mo_coeff, active_mo_struct)
2009 279 : CALL cp_fm_to_fm(mo_coeff, active_mo_coeff, nocc)
2010 :
2011 279 : energy_dimension = nocc
2012 279 : IF (dft_control%restricted .AND. &
2013 : scf_env%qs_ot_env(ot_channel)%settings%do_ener) THEN
2014 8 : energy_dimension = nocc*SIZE(kp%mos, 2)
2015 : END IF
2016 : CALL qs_ot_allocate(scf_env%qs_ot_env(ot_channel), &
2017 : matrix_sk_re, &
2018 279 : active_mo_struct, energy_dimension=energy_dimension)
2019 :
2020 279 : CALL dbcsr_set(scf_env%qs_ot_env(ot_channel)%matrix_c0, 0.0_dp)
2021 279 : CALL copy_fm_to_dbcsr(active_mo_coeff, scf_env%qs_ot_env(ot_channel)%matrix_c0)
2022 :
2023 279 : IF (.NOT. use_real_wfn) THEN
2024 279 : CPASSERT(SIZE(kp%mos, 1) >= 2)
2025 279 : mo_set_im => kp%mos(2, ispin)
2026 279 : CALL get_mo_set(mo_set=mo_set_im, mo_coeff=mo_coeff_im)
2027 279 : CPASSERT(ASSOCIATED(mo_coeff_im))
2028 279 : CALL cp_fm_create(active_mo_coeff_im, active_mo_struct)
2029 279 : CALL cp_fm_to_fm(mo_coeff_im, active_mo_coeff_im, nocc)
2030 279 : CALL qs_ot_allocate_complex_state(scf_env%qs_ot_env(ot_channel), matrix_sk_re)
2031 279 : CALL dbcsr_set(scf_env%qs_ot_env(ot_channel)%matrix_c0_im, 0.0_dp)
2032 279 : CALL dbcsr_set(scf_env%qs_ot_env(ot_channel)%matrix_sc0_im, 0.0_dp)
2033 279 : CALL copy_fm_to_dbcsr(active_mo_coeff_im, scf_env%qs_ot_env(ot_channel)%matrix_c0_im)
2034 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_sk_re, &
2035 : scf_env%qs_ot_env(ot_channel)%matrix_c0, 0.0_dp, &
2036 279 : scf_env%qs_ot_env(ot_channel)%matrix_sc0)
2037 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_sk_im, &
2038 : scf_env%qs_ot_env(ot_channel)%matrix_c0_im, 0.0_dp, &
2039 279 : scf_env%qs_ot_env(ot_channel)%matrix_x)
2040 : CALL dbcsr_add(scf_env%qs_ot_env(ot_channel)%matrix_sc0, &
2041 : scf_env%qs_ot_env(ot_channel)%matrix_x, &
2042 279 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2043 :
2044 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_sk_re, &
2045 : scf_env%qs_ot_env(ot_channel)%matrix_c0_im, 0.0_dp, &
2046 279 : scf_env%qs_ot_env(ot_channel)%matrix_sc0_im)
2047 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_sk_im, &
2048 : scf_env%qs_ot_env(ot_channel)%matrix_c0, 0.0_dp, &
2049 279 : scf_env%qs_ot_env(ot_channel)%matrix_x)
2050 : CALL dbcsr_add(scf_env%qs_ot_env(ot_channel)%matrix_sc0_im, &
2051 : scf_env%qs_ot_env(ot_channel)%matrix_x, &
2052 279 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2053 : END IF
2054 :
2055 279 : CALL qs_ot_init(scf_env%qs_ot_env(ot_channel))
2056 279 : IF (scf_env%qs_ot_env(ot_channel)%settings%do_ener) THEN
2057 120 : IF (dft_control%restricted) THEN
2058 24 : DO energy_spin = 1, SIZE(kp%mos, 2)
2059 16 : CALL get_mo_set(kp%mos(1, energy_spin), eigenvalues=eigenvalues)
2060 16 : energy_start = (energy_spin - 1)*nocc + 1
2061 : scf_env%qs_ot_env(ot_channel)%ener_x(energy_start:energy_start + nocc - 1) = &
2062 280 : eigenvalues(1:nocc)
2063 : END DO
2064 : ELSE
2065 1588 : scf_env%qs_ot_env(ot_channel)%ener_x(:) = eigenvalues(1:nocc)
2066 : END IF
2067 : END IF
2068 279 : scf_env%qs_ot_env(ot_channel)%broyden_adaptive_sigma = qs_env%broyden_adaptive_sigma
2069 :
2070 279 : CALL dbcsr_set(scf_env%qs_ot_env(ot_channel)%matrix_x, 0.0_dp)
2071 279 : CALL dbcsr_set(scf_env%qs_ot_env(ot_channel)%matrix_sx, 0.0_dp)
2072 279 : IF (scf_env%qs_ot_env(ot_channel)%has_complex_kpoint_state) THEN
2073 279 : CALL dbcsr_set(scf_env%qs_ot_env(ot_channel)%matrix_x_im, 0.0_dp)
2074 279 : CALL dbcsr_set(scf_env%qs_ot_env(ot_channel)%matrix_sx_im, 0.0_dp)
2075 : END IF
2076 :
2077 377 : SELECT CASE (scf_env%qs_ot_env(ot_channel)%settings%ot_algorithm)
2078 : CASE ("TOD")
2079 98 : IF (scf_env%qs_ot_env(ot_channel)%has_complex_kpoint_state) THEN
2080 : CALL qs_ot_get_p_complex(scf_env%qs_ot_env(ot_channel)%matrix_x, &
2081 : scf_env%qs_ot_env(ot_channel)%matrix_x_im, &
2082 : scf_env%qs_ot_env(ot_channel)%matrix_sx, &
2083 : scf_env%qs_ot_env(ot_channel)%matrix_sx_im, &
2084 98 : scf_env%qs_ot_env(ot_channel))
2085 : ELSE
2086 : CALL qs_ot_get_p(scf_env%qs_ot_env(ot_channel)%matrix_x, &
2087 : scf_env%qs_ot_env(ot_channel)%matrix_sx, &
2088 0 : scf_env%qs_ot_env(ot_channel))
2089 : END IF
2090 : CASE ("REF")
2091 : CALL dbcsr_copy(scf_env%qs_ot_env(ot_channel)%matrix_x, &
2092 181 : scf_env%qs_ot_env(ot_channel)%matrix_c0)
2093 : CALL dbcsr_copy(scf_env%qs_ot_env(ot_channel)%matrix_sx, &
2094 181 : scf_env%qs_ot_env(ot_channel)%matrix_sc0)
2095 181 : IF (scf_env%qs_ot_env(ot_channel)%has_complex_kpoint_state) THEN
2096 : CALL dbcsr_copy(scf_env%qs_ot_env(ot_channel)%matrix_x_im, &
2097 181 : scf_env%qs_ot_env(ot_channel)%matrix_c0_im)
2098 : CALL dbcsr_copy(scf_env%qs_ot_env(ot_channel)%matrix_sx_im, &
2099 181 : scf_env%qs_ot_env(ot_channel)%matrix_sc0_im)
2100 : CALL qs_ot_get_orbitals_ref_complex(scf_env%qs_ot_env(ot_channel)%matrix_c0, &
2101 : scf_env%qs_ot_env(ot_channel)%matrix_c0_im, &
2102 : matrix_sk_re, matrix_sk_im, &
2103 181 : scf_env%qs_ot_env(ot_channel))
2104 : END IF
2105 : CASE DEFAULT
2106 279 : CPABORT("Algorithm not yet implemented")
2107 : END SELECT
2108 :
2109 279 : CPASSERT(scf_env%qs_ot_env(ot_channel)%state_allocated)
2110 279 : IF (.NOT. use_real_wfn) THEN
2111 279 : CPASSERT(scf_env%qs_ot_env(ot_channel)%has_complex_kpoint_state)
2112 : END IF
2113 279 : CALL cp_fm_release(active_mo_coeff)
2114 279 : IF (.NOT. use_real_wfn) THEN
2115 279 : CALL cp_fm_release(active_mo_coeff_im)
2116 : END IF
2117 1086 : CALL cp_fm_struct_release(active_mo_struct)
2118 : END DO
2119 249 : CALL dbcsr_release_p(matrix_sk_re)
2120 403 : CALL dbcsr_release_p(matrix_sk_im)
2121 : END DO
2122 :
2123 : CALL qs_ot_check_channel_context(scf_env%qs_ot_env, dft_control%nspins, &
2124 : nkpoint=kpoints%nkp, &
2125 : restricted=ANY(scf_env%qs_ot_env(:)%restricted), &
2126 : require_kpoint=.TRUE., &
2127 : kp_range=kpoints%kp_range, &
2128 : require_local_state=.TRUE., &
2129 482 : require_complex_state=.NOT. use_real_wfn)
2130 :
2131 : CALL prepare_qs_ot_kpoint_preconditioners(scf_env, kpoints, matrix_ks_kp, &
2132 154 : matrix_s_kp, matrix_t_kp, nspin_ot)
2133 :
2134 154 : CALL timestop(handle)
2135 :
2136 154 : END SUBROUTINE allocate_qs_ot_kpoint_state
2137 :
2138 : ! **************************************************************************************************
2139 : !> \brief prepare and attach orbital preconditioners for complex k-point OT channels
2140 : !> \param scf_env ...
2141 : !> \param kpoints ...
2142 : !> \param matrix_ks ...
2143 : !> \param matrix_s ...
2144 : !> \param matrix_t ...
2145 : !> \param nspin_ot ...
2146 : ! **************************************************************************************************
2147 154 : SUBROUTINE prepare_qs_ot_kpoint_preconditioners(scf_env, kpoints, matrix_ks, matrix_s, matrix_t, nspin_ot)
2148 : TYPE(qs_scf_env_type), POINTER :: scf_env
2149 : TYPE(kpoint_type), POINTER :: kpoints
2150 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s, matrix_t
2151 : INTEGER, INTENT(IN) :: nspin_ot
2152 :
2153 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_qs_ot_kpoint_preconditioners'
2154 :
2155 : INTEGER :: handle, ikpoint, ispin, ks_spin, local_kpoint, n_ot_channels, nocc, &
2156 : occupation_spin, ot_channel, prec_type, source_channel
2157 154 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: occupation_signature
2158 154 : REAL(KIND=dp), DIMENSION(:), POINTER :: occupation_numbers
2159 : TYPE(dbcsr_type), POINTER :: matrix_h_im, matrix_h_re, matrix_s_im, &
2160 : matrix_s_re, matrix_t_im, matrix_t_re
2161 : TYPE(kpoint_env_type), POINTER :: kp
2162 :
2163 154 : CALL timeset(routineN, handle)
2164 :
2165 154 : NULLIFY (kp, matrix_h_im, matrix_h_re, matrix_s_im, matrix_s_re, matrix_t_im, matrix_t_re, &
2166 154 : occupation_numbers)
2167 154 : CPASSERT(ASSOCIATED(scf_env%qs_ot_env))
2168 154 : CPASSERT(ASSOCIATED(kpoints))
2169 154 : CPASSERT(ASSOCIATED(matrix_ks))
2170 154 : CPASSERT(ASSOCIATED(matrix_s))
2171 :
2172 154 : prec_type = scf_env%qs_ot_env(1)%settings%preconditioner_type
2173 154 : IF (prec_type == ot_precond_none) THEN
2174 8 : CALL timestop(handle)
2175 : RETURN
2176 : END IF
2177 0 : SELECT CASE (prec_type)
2178 : CASE (ot_precond_fermi_low_rank, ot_precond_full_all, ot_precond_full_all_covariant, &
2179 : ot_precond_full_single, ot_precond_full_single_inverse, &
2180 : ot_precond_full_kinetic, ot_precond_s_inverse)
2181 : CASE DEFAULT
2182 146 : CPABORT("Unsupported complex K-point OT preconditioner")
2183 : END SELECT
2184 146 : n_ot_channels = SIZE(scf_env%qs_ot_env)
2185 146 : IF (ASSOCIATED(scf_env%ot_preconditioner)) THEN
2186 174 : DO ot_channel = 1, SIZE(scf_env%ot_preconditioner)
2187 174 : IF (ASSOCIATED(scf_env%ot_preconditioner(ot_channel)%preconditioner)) THEN
2188 78 : CALL destroy_preconditioner(scf_env%ot_preconditioner(ot_channel)%preconditioner)
2189 78 : DEALLOCATE (scf_env%ot_preconditioner(ot_channel)%preconditioner)
2190 : END IF
2191 : END DO
2192 58 : DEALLOCATE (scf_env%ot_preconditioner)
2193 58 : NULLIFY (scf_env%ot_preconditioner)
2194 : END IF
2195 764 : ALLOCATE (scf_env%ot_preconditioner(n_ot_channels))
2196 387 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2197 241 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2198 241 : ikpoint = kp%nkpoint
2199 650 : DO ispin = 1, nspin_ot
2200 263 : ot_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot)
2201 263 : CPASSERT(scf_env%qs_ot_env(ot_channel)%state_allocated)
2202 1052 : ALLOCATE (scf_env%ot_preconditioner(ot_channel)%preconditioner)
2203 : CALL init_preconditioner(scf_env%ot_preconditioner(ot_channel)%preconditioner, &
2204 : para_env=scf_env%qs_ot_env(ot_channel)%para_env, &
2205 504 : blacs_env=scf_env%qs_ot_env(ot_channel)%blacs_env)
2206 : END DO
2207 : END DO
2208 :
2209 387 : DO local_kpoint = 1, SIZE(kpoints%kp_env)
2210 241 : kp => kpoints%kp_env(local_kpoint)%kpoint_env
2211 241 : CPASSERT(ASSOCIATED(kp))
2212 241 : CPASSERT(ASSOCIATED(kp%ot_hmat))
2213 241 : CPASSERT(ASSOCIATED(kp%ot_smat))
2214 241 : CPASSERT(SIZE(kp%ot_smat) >= 2)
2215 241 : ikpoint = kp%nkpoint
2216 : CALL kpoint_operator_get_local(matrix_s, kpoints, kp, 1, &
2217 : kp%ot_smat(1), kp%ot_smat(2), &
2218 241 : matrix_s_re, matrix_s_im)
2219 241 : IF (prec_type == ot_precond_full_kinetic) THEN
2220 64 : CPASSERT(ASSOCIATED(matrix_t))
2221 64 : CPASSERT(ASSOCIATED(kp%ot_tmat))
2222 64 : CPASSERT(SIZE(kp%ot_tmat) >= 2)
2223 : CALL kpoint_operator_get_local(matrix_t, kpoints, kp, 1, &
2224 : kp%ot_tmat(1), kp%ot_tmat(2), &
2225 64 : matrix_t_re, matrix_t_im)
2226 : END IF
2227 504 : DO ispin = 1, nspin_ot
2228 263 : ot_channel = qs_ot_channel_index(ispin, ikpoint, nspin_ot)
2229 263 : IF (ispin > 1 .AND. &
2230 : (prec_type == ot_precond_full_kinetic .OR. prec_type == ot_precond_s_inverse)) THEN
2231 10 : source_channel = qs_ot_channel_index(1, ikpoint, nspin_ot)
2232 10 : IF (.NOT. (ASSOCIATED(scf_env%ot_preconditioner(source_channel)%preconditioner%complex_fm))) THEN
2233 0 : CALL cp_abort(__LOCATION__, "Missing shared complex OT preconditioner")
2234 : END IF
2235 : scf_env%ot_preconditioner(ot_channel)%preconditioner%complex_fm => &
2236 10 : scf_env%ot_preconditioner(source_channel)%preconditioner%complex_fm
2237 10 : scf_env%ot_preconditioner(ot_channel)%preconditioner%owns_complex_fm = .FALSE.
2238 : scf_env%ot_preconditioner(ot_channel)%preconditioner%energy_gap = &
2239 10 : scf_env%ot_preconditioner(source_channel)%preconditioner%energy_gap
2240 : scf_env%ot_preconditioner(ot_channel)%preconditioner%in_use = &
2241 10 : scf_env%ot_preconditioner(source_channel)%preconditioner%in_use
2242 : scf_env%ot_preconditioner(ot_channel)%preconditioner%solver = &
2243 10 : scf_env%ot_preconditioner(source_channel)%preconditioner%solver
2244 : CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ot_channel), &
2245 10 : scf_env%ot_preconditioner(ot_channel)%preconditioner)
2246 10 : CYCLE
2247 : END IF
2248 : IF (prec_type == ot_precond_fermi_low_rank .OR. &
2249 : prec_type == ot_precond_full_all .OR. &
2250 : prec_type == ot_precond_full_all_covariant .OR. &
2251 253 : prec_type == ot_precond_full_single .OR. &
2252 : prec_type == ot_precond_full_single_inverse) THEN
2253 139 : ks_spin = MIN(ispin, SIZE(kp%ot_hmat, 2))
2254 139 : CPASSERT(SIZE(kp%ot_hmat, 1) >= 2)
2255 : CALL kpoint_operator_get_local(matrix_ks, kpoints, kp, ks_spin, &
2256 : kp%ot_hmat(1, ks_spin), kp%ot_hmat(2, ks_spin), &
2257 139 : matrix_h_re, matrix_h_im)
2258 : END IF
2259 6 : SELECT CASE (prec_type)
2260 : CASE (ot_precond_fermi_low_rank)
2261 : CALL make_preconditioner_complex_fermi_low_rank( &
2262 : scf_env%ot_preconditioner(ot_channel)%preconditioner, &
2263 : scf_env%qs_ot_env(ot_channel)%matrix_c0, &
2264 : scf_env%qs_ot_env(ot_channel)%matrix_c0_im, &
2265 : matrix_h_re, matrix_h_im, matrix_s_re, matrix_s_im, &
2266 : scf_env%qs_ot_env(ot_channel)%settings%energy_gap, &
2267 : scf_env%qs_ot_env(ot_channel)%settings%fermi_low_rank_max_rank, &
2268 6 : scf_env%qs_ot_env(ot_channel)%settings%precond_solver_type)
2269 : CASE (ot_precond_full_all)
2270 : CALL make_preconditioner_complex_full_all( &
2271 : scf_env%ot_preconditioner(ot_channel)%preconditioner, &
2272 : scf_env%qs_ot_env(ot_channel)%matrix_c0, &
2273 : scf_env%qs_ot_env(ot_channel)%matrix_c0_im, &
2274 : matrix_h_re, matrix_h_im, matrix_s_re, matrix_s_im, &
2275 : kp%mos(1, ispin), &
2276 : scf_env%qs_ot_env(ot_channel)%settings%energy_gap, &
2277 113 : scf_env%qs_ot_env(ot_channel)%settings%precond_solver_type)
2278 : CASE (ot_precond_full_single)
2279 : CALL make_preconditioner_complex_full_single( &
2280 : scf_env%ot_preconditioner(ot_channel)%preconditioner, &
2281 : matrix_h_re, matrix_h_im, matrix_s_re, matrix_s_im, &
2282 : kp%mos(1, ispin), &
2283 : scf_env%qs_ot_env(ot_channel)%settings%energy_gap, &
2284 4 : scf_env%qs_ot_env(ot_channel)%settings%precond_solver_type)
2285 : CASE (ot_precond_full_all_covariant)
2286 12 : IF (nspin_ot == 1 .AND. SIZE(kp%mos, 2) > 1) THEN
2287 : CALL dbcsr_get_info(scf_env%qs_ot_env(ot_channel)%matrix_c0, &
2288 8 : nfullcols_total=nocc)
2289 32 : ALLOCATE (occupation_signature(nocc, SIZE(kp%mos, 2)))
2290 24 : DO occupation_spin = 1, SIZE(kp%mos, 2)
2291 : CALL get_mo_set(kp%mos(1, occupation_spin), &
2292 16 : occupation_numbers=occupation_numbers)
2293 16 : CPASSERT(ASSOCIATED(occupation_numbers))
2294 16 : CPASSERT(SIZE(occupation_numbers) >= nocc)
2295 152 : occupation_signature(:, occupation_spin) = occupation_numbers(1:nocc)
2296 : END DO
2297 : CALL make_preconditioner_complex_full_all_covariant( &
2298 : scf_env%ot_preconditioner(ot_channel)%preconditioner, &
2299 : scf_env%qs_ot_env(ot_channel)%matrix_c0, &
2300 : scf_env%qs_ot_env(ot_channel)%matrix_c0_im, &
2301 : matrix_h_re, matrix_h_im, matrix_s_re, matrix_s_im, &
2302 : scf_env%qs_ot_env(ot_channel)%settings%energy_gap, &
2303 : scf_env%qs_ot_env(ot_channel)%settings%precond_solver_type, &
2304 8 : occupation_signature)
2305 16 : DEALLOCATE (occupation_signature)
2306 : ELSE
2307 : CALL make_preconditioner_complex_full_all_covariant( &
2308 : scf_env%ot_preconditioner(ot_channel)%preconditioner, &
2309 : scf_env%qs_ot_env(ot_channel)%matrix_c0, &
2310 : scf_env%qs_ot_env(ot_channel)%matrix_c0_im, &
2311 : matrix_h_re, matrix_h_im, matrix_s_re, matrix_s_im, &
2312 : scf_env%qs_ot_env(ot_channel)%settings%energy_gap, &
2313 4 : scf_env%qs_ot_env(ot_channel)%settings%precond_solver_type)
2314 : END IF
2315 : CASE (ot_precond_full_single_inverse)
2316 : CALL make_preconditioner_complex_full_single_inverse( &
2317 : scf_env%ot_preconditioner(ot_channel)%preconditioner, &
2318 : scf_env%qs_ot_env(ot_channel)%matrix_c0, &
2319 : scf_env%qs_ot_env(ot_channel)%matrix_c0_im, &
2320 : matrix_h_re, matrix_h_im, matrix_s_re, matrix_s_im, &
2321 : scf_env%qs_ot_env(ot_channel)%settings%energy_gap, &
2322 4 : scf_env%qs_ot_env(ot_channel)%settings%precond_solver_type)
2323 : CASE (ot_precond_full_kinetic)
2324 : CALL make_preconditioner_complex_full_kinetic( &
2325 : scf_env%ot_preconditioner(ot_channel)%preconditioner, &
2326 : matrix_t_re, matrix_t_im, matrix_s_re, matrix_s_im, &
2327 : scf_env%qs_ot_env(ot_channel)%settings%energy_gap, &
2328 64 : scf_env%qs_ot_env(ot_channel)%settings%precond_solver_type)
2329 : CASE (ot_precond_s_inverse)
2330 : CALL make_preconditioner_complex_full_s_inverse( &
2331 : scf_env%ot_preconditioner(ot_channel)%preconditioner, &
2332 : matrix_s_re, matrix_s_im, &
2333 253 : scf_env%qs_ot_env(ot_channel)%settings%precond_solver_type)
2334 : END SELECT
2335 : CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ot_channel), &
2336 253 : scf_env%ot_preconditioner(ot_channel)%preconditioner)
2337 253 : IF (ASSOCIATED(matrix_h_re)) CALL dbcsr_release_p(matrix_h_re)
2338 494 : IF (ASSOCIATED(matrix_h_im)) CALL dbcsr_release_p(matrix_h_im)
2339 : END DO
2340 241 : CALL dbcsr_release_p(matrix_s_re)
2341 241 : CALL dbcsr_release_p(matrix_s_im)
2342 241 : IF (ASSOCIATED(matrix_t_re)) CALL dbcsr_release_p(matrix_t_re)
2343 387 : IF (ASSOCIATED(matrix_t_im)) CALL dbcsr_release_p(matrix_t_im)
2344 : END DO
2345 :
2346 146 : CALL timestop(handle)
2347 :
2348 300 : END SUBROUTINE prepare_qs_ot_kpoint_preconditioners
2349 :
2350 : ! **************************************************************************************************
2351 : !> \brief perform cleanup operations (like releasing temporary storage)
2352 : !> at the end of the scf
2353 : !> \param scf_env ...
2354 : !> \par History
2355 : !> 02.2003 created [fawzi]
2356 : !> \author fawzi
2357 : ! **************************************************************************************************
2358 24195 : SUBROUTINE scf_env_cleanup(scf_env)
2359 : TYPE(qs_scf_env_type), INTENT(INOUT) :: scf_env
2360 :
2361 : CHARACTER(len=*), PARAMETER :: routineN = 'scf_env_cleanup'
2362 :
2363 : INTEGER :: handle
2364 :
2365 24195 : CALL timeset(routineN, handle)
2366 :
2367 : ! Release SCF work storage
2368 24195 : CALL cp_fm_release(scf_env%scf_work1)
2369 :
2370 24195 : IF (ASSOCIATED(scf_env%scf_work1_red)) THEN
2371 48 : CALL cp_fm_release(scf_env%scf_work1_red)
2372 : END IF
2373 24195 : IF (ASSOCIATED(scf_env%scf_work2)) THEN
2374 17538 : CALL cp_fm_release(scf_env%scf_work2)
2375 17538 : DEALLOCATE (scf_env%scf_work2)
2376 : NULLIFY (scf_env%scf_work2)
2377 : END IF
2378 24195 : IF (ASSOCIATED(scf_env%scf_work2_red)) THEN
2379 48 : CALL cp_fm_release(scf_env%scf_work2_red)
2380 48 : DEALLOCATE (scf_env%scf_work2_red)
2381 : NULLIFY (scf_env%scf_work2_red)
2382 : END IF
2383 24195 : IF (ASSOCIATED(scf_env%ortho)) THEN
2384 14868 : CALL cp_fm_release(scf_env%ortho)
2385 14868 : DEALLOCATE (scf_env%ortho)
2386 : NULLIFY (scf_env%ortho)
2387 : END IF
2388 24195 : IF (ASSOCIATED(scf_env%ortho_red)) THEN
2389 48 : CALL cp_fm_release(scf_env%ortho_red)
2390 48 : DEALLOCATE (scf_env%ortho_red)
2391 : NULLIFY (scf_env%ortho_red)
2392 : END IF
2393 24195 : IF (ASSOCIATED(scf_env%ortho_m1)) THEN
2394 56 : CALL cp_fm_release(scf_env%ortho_m1)
2395 56 : DEALLOCATE (scf_env%ortho_m1)
2396 : NULLIFY (scf_env%ortho_m1)
2397 : END IF
2398 24195 : IF (ASSOCIATED(scf_env%ortho_m1_red)) THEN
2399 6 : CALL cp_fm_release(scf_env%ortho_m1_red)
2400 6 : DEALLOCATE (scf_env%ortho_m1_red)
2401 : NULLIFY (scf_env%ortho_m1_red)
2402 : END IF
2403 :
2404 24195 : IF (ASSOCIATED(scf_env%ortho_dbcsr)) THEN
2405 58 : CALL dbcsr_deallocate_matrix(scf_env%ortho_dbcsr)
2406 : END IF
2407 24195 : IF (ASSOCIATED(scf_env%buf1_dbcsr)) THEN
2408 58 : CALL dbcsr_deallocate_matrix(scf_env%buf1_dbcsr)
2409 : END IF
2410 24195 : IF (ASSOCIATED(scf_env%buf2_dbcsr)) THEN
2411 58 : CALL dbcsr_deallocate_matrix(scf_env%buf2_dbcsr)
2412 : END IF
2413 :
2414 24195 : IF (ASSOCIATED(scf_env%p_mix_new)) THEN
2415 17588 : CALL dbcsr_deallocate_matrix_set(scf_env%p_mix_new)
2416 : END IF
2417 :
2418 24195 : IF (ASSOCIATED(scf_env%p_delta)) THEN
2419 838 : CALL dbcsr_deallocate_matrix_set(scf_env%p_delta)
2420 : END IF
2421 :
2422 : ! Method dependent cleanup
2423 24245 : SELECT CASE (scf_env%method)
2424 : CASE (ot_method_nr)
2425 : !
2426 : CASE (ot_diag_method_nr)
2427 : !
2428 : CASE (general_diag_method_nr)
2429 : !
2430 : CASE (special_diag_method_nr)
2431 : !
2432 : CASE (block_krylov_diag_method_nr)
2433 : CASE (block_davidson_diag_method_nr)
2434 50 : CALL block_davidson_deallocate(scf_env%block_davidson_env)
2435 : CASE (filter_matrix_diag_method_nr)
2436 : !
2437 : CASE (smeagol_method_nr)
2438 : !
2439 : CASE DEFAULT
2440 24195 : CPABORT("unknown scf method method:"//cp_to_string(scf_env%method))
2441 : END SELECT
2442 :
2443 24195 : IF (ASSOCIATED(scf_env%outer_scf%variables)) THEN
2444 4351 : DEALLOCATE (scf_env%outer_scf%variables)
2445 : END IF
2446 24195 : IF (ASSOCIATED(scf_env%outer_scf%count)) THEN
2447 4351 : DEALLOCATE (scf_env%outer_scf%count)
2448 : END IF
2449 24195 : IF (ASSOCIATED(scf_env%outer_scf%gradient)) THEN
2450 4351 : DEALLOCATE (scf_env%outer_scf%gradient)
2451 : END IF
2452 24195 : IF (ASSOCIATED(scf_env%outer_scf%energy)) THEN
2453 4351 : DEALLOCATE (scf_env%outer_scf%energy)
2454 : END IF
2455 24195 : IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian) .AND. &
2456 : scf_env%outer_scf%deallocate_jacobian) THEN
2457 50 : DEALLOCATE (scf_env%outer_scf%inv_jacobian)
2458 : END IF
2459 :
2460 24195 : CALL timestop(handle)
2461 :
2462 24195 : END SUBROUTINE scf_env_cleanup
2463 :
2464 : ! **************************************************************************************************
2465 : !> \brief perform a CDFT scf procedure in the given qs_env
2466 : !> \param qs_env the qs_environment where to perform the scf procedure
2467 : !> \param should_stop flag determining if calculation should stop
2468 : !> \param has_converged both the electronic and constraint loops converged
2469 : !> \param total_scf_steps number of electronic SCF steps over the constraint loop
2470 : !> \par History
2471 : !> 12.2015 Created
2472 : !> \author Nico Holmberg
2473 : ! **************************************************************************************************
2474 360 : SUBROUTINE cdft_scf(qs_env, should_stop, has_converged, total_scf_steps)
2475 : TYPE(qs_environment_type), POINTER :: qs_env
2476 : LOGICAL, INTENT(OUT) :: should_stop, has_converged
2477 : INTEGER, INTENT(OUT) :: total_scf_steps
2478 :
2479 : CHARACTER(len=*), PARAMETER :: routineN = 'cdft_scf'
2480 :
2481 : INTEGER :: handle, iatom, iimage, ispin, ivar, nmo, &
2482 : nvar, output_unit, tsteps
2483 : LOGICAL :: cdft_loop_converged, converged, &
2484 : exit_cdft_loop, first_iteration, &
2485 : my_uocc, uniform_occupation
2486 360 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: gapw_cdft_values
2487 360 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_occupations
2488 360 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2489 : TYPE(cdft_control_type), POINTER :: cdft_control
2490 : TYPE(cp_logger_type), POINTER :: logger
2491 360 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: gapw_wmat, matrix_s, rho_ao
2492 : TYPE(dft_control_type), POINTER :: dft_control
2493 : TYPE(local_rho_type), POINTER :: gapw_operator_rho
2494 360 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2495 : TYPE(mp_para_env_type), POINTER :: para_env
2496 : TYPE(pw_env_type), POINTER :: pw_env
2497 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2498 : TYPE(qs_energy_type), POINTER :: energy
2499 360 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2500 : TYPE(qs_ks_env_type), POINTER :: ks_env
2501 : TYPE(qs_rho_type), POINTER :: rho
2502 : TYPE(qs_scf_env_type), POINTER :: scf_env
2503 : TYPE(scf_control_type), POINTER :: scf_control
2504 : TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
2505 :
2506 360 : NULLIFY (atomic_kind_set, gapw_operator_rho, gapw_wmat, para_env, qs_kind_set, &
2507 360 : scf_env, ks_env, energy, rho, matrix_s, rho_ao, cdft_control, logger, &
2508 360 : dft_control, pw_env, auxbas_pw_pool, energy, ks_env, scf_env, dft_section, &
2509 360 : input, scf_section, scf_control, mos, mo_occupations)
2510 720 : logger => cp_get_default_logger()
2511 :
2512 360 : CPASSERT(ASSOCIATED(qs_env))
2513 : CALL get_qs_env(qs_env, scf_env=scf_env, energy=energy, &
2514 : dft_control=dft_control, scf_control=scf_control, &
2515 360 : ks_env=ks_env, input=input)
2516 :
2517 360 : CALL timeset(routineN//"_loop", handle)
2518 360 : dft_section => section_vals_get_subs_vals(input, "DFT")
2519 360 : scf_section => section_vals_get_subs_vals(dft_section, "SCF")
2520 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", &
2521 360 : extension=".scfLog")
2522 360 : first_iteration = .TRUE.
2523 :
2524 360 : cdft_control => dft_control%qs_control%cdft_control
2525 :
2526 360 : scf_env%outer_scf%iter_count = 0
2527 360 : cdft_control%total_steps = 0
2528 360 : total_scf_steps = 0
2529 :
2530 : ! Write some info about the CDFT calculation
2531 360 : IF (output_unit > 0) THEN
2532 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
2533 199 : "CDFT EXTERNAL SCF WAVEFUNCTION OPTIMIZATION"
2534 199 : CALL qs_scf_cdft_initial_info(output_unit, cdft_control)
2535 : END IF
2536 360 : IF (cdft_control%reuse_precond) THEN
2537 0 : reuse_precond = .FALSE.
2538 0 : cdft_control%nreused = 0
2539 : END IF
2540 596 : cdft_outer_loop: DO
2541 : ! Change outer_scf settings to OT settings
2542 596 : CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot)
2543 : ! Solve electronic structure with fixed value of constraint
2544 : CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
2545 596 : converged=converged, should_stop=should_stop, total_scf_steps=tsteps)
2546 596 : total_scf_steps = total_scf_steps + tsteps
2547 : ! Decide whether to reuse the preconditioner on the next iteration
2548 596 : IF (cdft_control%reuse_precond) THEN
2549 : ! For convergence in exactly one step, the preconditioner is always reused (assuming max_reuse > 0)
2550 : ! usually this means that the electronic structure has already converged to the correct state
2551 : ! but the constraint optimizer keeps jumping over the optimal solution
2552 : IF (scf_env%outer_scf%iter_count == 1 .AND. scf_env%iter_count == 1 &
2553 0 : .AND. cdft_control%total_steps /= 1) THEN
2554 0 : cdft_control%nreused = cdft_control%nreused - 1
2555 : END IF
2556 : ! SCF converged in less than precond_freq steps
2557 : IF (scf_env%outer_scf%iter_count == 1 .AND. scf_env%iter_count <= cdft_control%precond_freq .AND. &
2558 0 : cdft_control%total_steps /= 1 .AND. cdft_control%nreused < cdft_control%max_reuse) THEN
2559 0 : reuse_precond = .TRUE.
2560 0 : cdft_control%nreused = cdft_control%nreused + 1
2561 : ELSE
2562 0 : reuse_precond = .FALSE.
2563 0 : cdft_control%nreused = 0
2564 : END IF
2565 : END IF
2566 : ! Update history purging counters
2567 596 : IF (first_iteration .AND. cdft_control%purge_history) THEN
2568 0 : cdft_control%istep = cdft_control%istep + 1
2569 0 : IF (scf_env%outer_scf%iter_count > 1) THEN
2570 0 : cdft_control%nbad_conv = cdft_control%nbad_conv + 1
2571 0 : IF (cdft_control%nbad_conv >= cdft_control%purge_freq .AND. &
2572 : cdft_control%istep >= cdft_control%purge_offset) THEN
2573 0 : cdft_control%nbad_conv = 0
2574 0 : cdft_control%istep = 0
2575 0 : cdft_control%should_purge = .TRUE.
2576 : END IF
2577 : END IF
2578 : END IF
2579 596 : first_iteration = .FALSE.
2580 : ! Change outer_scf settings to CDFT settings
2581 596 : CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft)
2582 : CALL qs_scf_check_outer_exit(qs_env, scf_env, scf_control, should_stop, &
2583 596 : cdft_loop_converged, exit_cdft_loop)
2584 : CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
2585 : energy, cdft_control%total_steps, &
2586 596 : should_stop, cdft_loop_converged, cdft_loop=.TRUE.)
2587 596 : IF (exit_cdft_loop) EXIT cdft_outer_loop
2588 : ! Check if the inverse Jacobian needs to be calculated
2589 236 : CALL qs_calculate_inverse_jacobian(qs_env)
2590 : ! Check if a line search should be performed to find an optimal step size for the optimizer
2591 236 : CALL qs_cdft_line_search(qs_env)
2592 : ! Optimize constraint
2593 236 : CALL outer_loop_optimize(scf_env, scf_control)
2594 236 : CALL outer_loop_update_qs_env(qs_env, scf_env)
2595 596 : CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.)
2596 : END DO cdft_outer_loop
2597 :
2598 360 : has_converged = converged .AND. cdft_loop_converged
2599 360 : cdft_control%ienergy = cdft_control%ienergy + 1
2600 :
2601 : ! Store needed arrays for ET coupling calculation
2602 360 : IF (cdft_control%do_et) THEN
2603 190 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, mos=mos)
2604 190 : nvar = SIZE(cdft_control%target)
2605 190 : IF (dft_control%qs_control%gapw) THEN
2606 4 : IF (dft_control%nimages /= 1) THEN
2607 : CALL cp_abort(__LOCATION__, &
2608 0 : "GAPW CDFT-CI currently requires a Gamma-point calculation.")
2609 : END IF
2610 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, para_env=para_env, &
2611 4 : qs_kind_set=qs_kind_set, rho=rho)
2612 4 : CALL qs_rho_get(rho, rho_ao=rho_ao)
2613 4 : CALL local_rho_set_create(gapw_operator_rho)
2614 : CALL allocate_rho_atom_internals(gapw_operator_rho%rho_atom_set, atomic_kind_set, &
2615 4 : qs_kind_set, dft_control, para_env)
2616 24 : ALLOCATE (gapw_cdft_values(nvar), gapw_wmat(dft_control%nspins*dft_control%nimages))
2617 8 : DO iimage = 1, dft_control%nimages
2618 12 : DO ispin = 1, dft_control%nspins
2619 4 : CALL dbcsr_init_p(gapw_wmat(dft_control%nspins*(iimage - 1) + ispin)%matrix)
2620 : CALL dbcsr_copy(gapw_wmat(dft_control%nspins*(iimage - 1) + ispin)%matrix, &
2621 8 : matrix_s(iimage)%matrix, name="GAPW CDFT WEIGHT MATRIX")
2622 : END DO
2623 : END DO
2624 : END IF
2625 : ! Matrix representation of weight function
2626 764 : ALLOCATE (cdft_control%wmat(nvar))
2627 384 : DO ivar = 1, nvar
2628 194 : CALL dbcsr_init_p(cdft_control%wmat(ivar)%matrix)
2629 : CALL dbcsr_copy(cdft_control%wmat(ivar)%matrix, matrix_s(1)%matrix, &
2630 194 : name="ET_RESTRAINT_MATRIX")
2631 194 : CALL dbcsr_set(cdft_control%wmat(ivar)%matrix, 0.0_dp)
2632 : CALL integrate_v_rspace(cdft_control%group(ivar)%weight, &
2633 : hmat=cdft_control%wmat(ivar), qs_env=qs_env, &
2634 : calculate_forces=.FALSE., &
2635 194 : gapw=dft_control%qs_control%gapw)
2636 384 : IF (dft_control%qs_control%gapw) THEN
2637 8 : DO ispin = 1, SIZE(gapw_wmat)
2638 8 : CALL dbcsr_set(gapw_wmat(ispin)%matrix, 0.0_dp)
2639 : END DO
2640 4 : CALL zero_rho_atom_integrals(gapw_operator_rho%rho_atom_set)
2641 : CALL gapw_cdft_one_center(qs_env, energy_only=.FALSE., calculate_forces=.FALSE., &
2642 : values=gapw_cdft_values, operator_group=ivar, &
2643 4 : rho_atom_operator_set=gapw_operator_rho%rho_atom_set)
2644 : CALL update_ks_atom(qs_env, gapw_wmat, rho_ao, forces=.FALSE., &
2645 4 : rho_atom_external=gapw_operator_rho%rho_atom_set)
2646 4 : CALL dbcsr_add(cdft_control%wmat(ivar)%matrix, gapw_wmat(1)%matrix, 1.0_dp, 1.0_dp)
2647 : END IF
2648 : END DO
2649 190 : IF (dft_control%qs_control%gapw) THEN
2650 4 : CALL dbcsr_deallocate_matrix_set(gapw_wmat)
2651 4 : CALL local_rho_set_release(gapw_operator_rho)
2652 4 : DEALLOCATE (gapw_cdft_values)
2653 : END IF
2654 : ! Overlap matrix
2655 190 : CALL dbcsr_init_p(cdft_control%matrix_s%matrix)
2656 : CALL dbcsr_copy(cdft_control%matrix_s%matrix, matrix_s(1)%matrix, &
2657 190 : name="OVERLAP")
2658 : ! Molecular orbital coefficients
2659 190 : NULLIFY (cdft_control%mo_coeff)
2660 936 : ALLOCATE (cdft_control%mo_coeff(dft_control%nspins))
2661 556 : DO ispin = 1, dft_control%nspins
2662 : CALL cp_fm_create(matrix=cdft_control%mo_coeff(ispin), &
2663 : matrix_struct=qs_env%mos(ispin)%mo_coeff%matrix_struct, &
2664 366 : name="MO_COEFF_A"//TRIM(ADJUSTL(cp_to_string(ispin)))//"MATRIX")
2665 : CALL cp_fm_to_fm(qs_env%mos(ispin)%mo_coeff, &
2666 556 : cdft_control%mo_coeff(ispin))
2667 : END DO
2668 : ! Density matrix
2669 190 : IF (cdft_control%calculate_metric) THEN
2670 24 : CALL get_qs_env(qs_env, rho=rho)
2671 24 : CALL qs_rho_get(rho, rho_ao=rho_ao)
2672 120 : ALLOCATE (cdft_control%matrix_p(dft_control%nspins))
2673 72 : DO ispin = 1, dft_control%nspins
2674 48 : NULLIFY (cdft_control%matrix_p(ispin)%matrix)
2675 48 : CALL dbcsr_init_p(cdft_control%matrix_p(ispin)%matrix)
2676 : CALL dbcsr_copy(cdft_control%matrix_p(ispin)%matrix, rho_ao(ispin)%matrix, &
2677 72 : name="DENSITY MATRIX")
2678 : END DO
2679 : END IF
2680 : ! Copy occupation numbers if non-uniform occupation
2681 190 : uniform_occupation = .TRUE.
2682 556 : DO ispin = 1, dft_control%nspins
2683 366 : CALL get_mo_set(mo_set=mos(ispin), uniform_occupation=my_uocc)
2684 612 : uniform_occupation = uniform_occupation .AND. my_uocc
2685 : END DO
2686 190 : IF (.NOT. uniform_occupation) THEN
2687 140 : ALLOCATE (cdft_control%occupations(dft_control%nspins))
2688 84 : DO ispin = 1, dft_control%nspins
2689 : CALL get_mo_set(mo_set=mos(ispin), &
2690 : nmo=nmo, &
2691 56 : occupation_numbers=mo_occupations)
2692 168 : ALLOCATE (cdft_control%occupations(ispin)%array(nmo))
2693 588 : cdft_control%occupations(ispin)%array(1:nmo) = mo_occupations(1:nmo)
2694 : END DO
2695 : END IF
2696 : END IF
2697 :
2698 : ! Deallocate constraint storage if forces are not needed
2699 : ! In case of a simulation with multiple force_evals,
2700 : ! deallocate only if weight function should not be copied to different force_evals
2701 360 : IF (.NOT. (cdft_control%save_pot .OR. cdft_control%transfer_pot)) THEN
2702 168 : CALL get_qs_env(qs_env, pw_env=pw_env)
2703 168 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
2704 348 : DO iatom = 1, SIZE(cdft_control%group)
2705 180 : CALL auxbas_pw_pool%give_back_pw(cdft_control%group(iatom)%weight)
2706 348 : DEALLOCATE (cdft_control%group(iatom)%weight)
2707 : END DO
2708 168 : IF (cdft_control%atomic_charges) THEN
2709 266 : DO iatom = 1, cdft_control%natoms
2710 266 : CALL auxbas_pw_pool%give_back_pw(cdft_control%charge(iatom))
2711 : END DO
2712 88 : DEALLOCATE (cdft_control%charge)
2713 : END IF
2714 168 : IF (cdft_control%type == outer_scf_becke_constraint .AND. &
2715 : cdft_control%becke_control%cavity_confine) THEN
2716 122 : IF (.NOT. ASSOCIATED(cdft_control%becke_control%cavity_mat)) THEN
2717 112 : CALL auxbas_pw_pool%give_back_pw(cdft_control%becke_control%cavity)
2718 : ELSE
2719 10 : DEALLOCATE (cdft_control%becke_control%cavity_mat)
2720 : END IF
2721 46 : ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
2722 22 : IF (ASSOCIATED(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)) THEN
2723 0 : CALL auxbas_pw_pool%give_back_pw(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)
2724 : END IF
2725 : END IF
2726 168 : IF (ASSOCIATED(cdft_control%charges_fragment)) DEALLOCATE (cdft_control%charges_fragment)
2727 168 : cdft_control%need_pot = .TRUE.
2728 168 : cdft_control%external_control = .FALSE.
2729 : END IF
2730 :
2731 360 : CALL timestop(handle)
2732 :
2733 720 : END SUBROUTINE cdft_scf
2734 :
2735 : ! **************************************************************************************************
2736 : !> \brief perform cleanup operations for cdft_control
2737 : !> \param cdft_control container for the external CDFT SCF loop variables
2738 : !> \par History
2739 : !> 12.2015 created [Nico Holmberg]
2740 : !> \author Nico Holmberg
2741 : ! **************************************************************************************************
2742 360 : SUBROUTINE cdft_control_cleanup(cdft_control)
2743 : TYPE(cdft_control_type), POINTER :: cdft_control
2744 :
2745 360 : IF (ASSOCIATED(cdft_control%constraint%variables)) THEN
2746 360 : DEALLOCATE (cdft_control%constraint%variables)
2747 : END IF
2748 360 : IF (ASSOCIATED(cdft_control%constraint%count)) THEN
2749 360 : DEALLOCATE (cdft_control%constraint%count)
2750 : END IF
2751 360 : IF (ASSOCIATED(cdft_control%constraint%gradient)) THEN
2752 360 : DEALLOCATE (cdft_control%constraint%gradient)
2753 : END IF
2754 360 : IF (ASSOCIATED(cdft_control%constraint%energy)) THEN
2755 360 : DEALLOCATE (cdft_control%constraint%energy)
2756 : END IF
2757 360 : IF (ASSOCIATED(cdft_control%constraint%inv_jacobian) .AND. &
2758 : cdft_control%constraint%deallocate_jacobian) THEN
2759 4 : DEALLOCATE (cdft_control%constraint%inv_jacobian)
2760 : END IF
2761 :
2762 360 : END SUBROUTINE cdft_control_cleanup
2763 :
2764 : ! **************************************************************************************************
2765 : !> \brief Calculates the finite difference inverse Jacobian
2766 : !> \param qs_env the qs_environment_type where to compute the Jacobian
2767 : !> \par History
2768 : !> 01.2017 created [Nico Holmberg]
2769 : ! **************************************************************************************************
2770 236 : SUBROUTINE qs_calculate_inverse_jacobian(qs_env)
2771 : TYPE(qs_environment_type), POINTER :: qs_env
2772 :
2773 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_calculate_inverse_jacobian'
2774 :
2775 : CHARACTER(len=default_path_length) :: project_name
2776 : INTEGER :: counter, handle, i, ispin, iter_count, &
2777 : iwork, j, max_scf, nspins, nsteps, &
2778 : nvar, nwork, output_unit, pwork, &
2779 : tsteps, twork
2780 : LOGICAL :: converged, explicit_jacobian, &
2781 : should_build, should_stop, &
2782 : use_md_history
2783 : REAL(KIND=dp) :: inv_error, step_size
2784 236 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coeff, dh, step_multiplier
2785 236 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: jacobian
2786 236 : REAL(KIND=dp), DIMENSION(:), POINTER :: energy
2787 236 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient, inv_jacobian
2788 : TYPE(cdft_control_type), POINTER :: cdft_control
2789 : TYPE(cp_logger_type), POINTER :: logger, tmp_logger
2790 236 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv
2791 236 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
2792 : TYPE(dft_control_type), POINTER :: dft_control
2793 236 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mos_stashed
2794 236 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2795 : TYPE(mp_para_env_type), POINTER :: para_env
2796 : TYPE(qs_energy_type), POINTER :: energy_qs
2797 : TYPE(qs_ks_env_type), POINTER :: ks_env
2798 : TYPE(qs_rho_type), POINTER :: rho
2799 : TYPE(qs_scf_env_type), POINTER :: scf_env
2800 : TYPE(scf_control_type), POINTER :: scf_control
2801 :
2802 236 : NULLIFY (energy, gradient, p_rmpv, rho_ao_kp, mos, rho, &
2803 236 : ks_env, scf_env, scf_control, dft_control, cdft_control, &
2804 236 : inv_jacobian, para_env, tmp_logger, energy_qs)
2805 472 : logger => cp_get_default_logger()
2806 :
2807 236 : CPASSERT(ASSOCIATED(qs_env))
2808 : CALL get_qs_env(qs_env, scf_env=scf_env, ks_env=ks_env, &
2809 : scf_control=scf_control, mos=mos, rho=rho, &
2810 : dft_control=dft_control, &
2811 236 : para_env=para_env, energy=energy_qs)
2812 236 : explicit_jacobian = .FALSE.
2813 236 : should_build = .FALSE.
2814 236 : use_md_history = .FALSE.
2815 236 : iter_count = scf_env%outer_scf%iter_count
2816 : ! Quick exit if optimizer does not require Jacobian
2817 236 : IF (.NOT. ASSOCIATED(scf_control%outer_scf%cdft_opt_control)) RETURN
2818 : ! Check if Jacobian should be calculated and initialize
2819 118 : CALL timeset(routineN, handle)
2820 118 : CALL initialize_inverse_jacobian(scf_control, scf_env, explicit_jacobian, should_build, used_history)
2821 118 : IF (scf_control%outer_scf%cdft_opt_control%jacobian_restart) THEN
2822 : ! Restart from previously calculated inverse Jacobian
2823 6 : should_build = .FALSE.
2824 6 : CALL restart_inverse_jacobian(qs_env)
2825 : END IF
2826 118 : IF (should_build) THEN
2827 78 : scf_env%outer_scf%deallocate_jacobian = .FALSE.
2828 : ! Actually need to (re)build the Jacobian
2829 78 : IF (explicit_jacobian) THEN
2830 : ! Build Jacobian with finite differences
2831 62 : cdft_control => dft_control%qs_control%cdft_control
2832 62 : IF (.NOT. ASSOCIATED(cdft_control)) THEN
2833 : CALL cp_abort(__LOCATION__, &
2834 : "Optimizers that need the explicit Jacobian can"// &
2835 0 : " only be used together with a valid CDFT constraint.")
2836 : END IF
2837 : ! Redirect output from Jacobian calculation to a new file by creating a temporary logger
2838 62 : project_name = logger%iter_info%project_name
2839 62 : CALL create_tmp_logger(para_env, project_name, "-JacobianInfo.out", output_unit, tmp_logger)
2840 : ! Save last converged state so we can roll back to it (mo_coeff and some outer_loop variables)
2841 62 : nspins = dft_control%nspins
2842 310 : ALLOCATE (mos_stashed(nspins))
2843 186 : DO ispin = 1, nspins
2844 186 : CALL duplicate_mo_set(mos_stashed(ispin), mos(ispin))
2845 : END DO
2846 62 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
2847 62 : p_rmpv => rho_ao_kp(:, 1)
2848 : ! Allocate work
2849 62 : nvar = SIZE(scf_env%outer_scf%variables, 1)
2850 62 : max_scf = scf_control%outer_scf%max_scf + 1
2851 248 : ALLOCATE (gradient(nvar, max_scf))
2852 1310 : gradient = scf_env%outer_scf%gradient
2853 186 : ALLOCATE (energy(max_scf))
2854 594 : energy = scf_env%outer_scf%energy
2855 248 : ALLOCATE (jacobian(nvar, nvar))
2856 62 : jacobian = 0.0_dp
2857 62 : nsteps = cdft_control%total_steps
2858 : ! Setup finite difference scheme
2859 62 : CALL prepare_jacobian_stencil(qs_env, output_unit, nwork, pwork, coeff, step_multiplier, dh)
2860 62 : twork = pwork - nwork
2861 148 : DO i = 1, nvar
2862 282 : jacobian(i, :) = coeff(0)*scf_env%outer_scf%gradient(i, iter_count)
2863 : END DO
2864 : ! Calculate the Jacobian by perturbing each Lagrangian and recalculating the energy self-consistently
2865 62 : CALL cp_add_default_logger(tmp_logger)
2866 148 : DO i = 1, nvar
2867 86 : IF (output_unit > 0) THEN
2868 43 : WRITE (output_unit, FMT="(A)") " "
2869 43 : WRITE (output_unit, FMT="(A)") " #####################################"
2870 : WRITE (output_unit, '(A,I3,A,I3,A)') &
2871 43 : " ### Constraint ", i, " of ", nvar, " ###"
2872 43 : WRITE (output_unit, FMT="(A)") " #####################################"
2873 : END IF
2874 86 : counter = 0
2875 332 : DO iwork = nwork, pwork
2876 184 : IF (iwork == 0) CYCLE
2877 98 : counter = counter + 1
2878 98 : IF (output_unit > 0) THEN
2879 49 : WRITE (output_unit, FMT="(A)") " #####################################"
2880 : WRITE (output_unit, '(A,I3,A,I3,A)') &
2881 49 : " ### Energy evaluation ", counter, " of ", twork, " ###"
2882 49 : WRITE (output_unit, FMT="(A)") " #####################################"
2883 : END IF
2884 98 : IF (SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step) == 1) THEN
2885 90 : step_size = scf_control%outer_scf%cdft_opt_control%jacobian_step(1)
2886 : ELSE
2887 8 : step_size = scf_control%outer_scf%cdft_opt_control%jacobian_step(i)
2888 : END IF
2889 244 : scf_env%outer_scf%variables(:, iter_count + 1) = scf_env%outer_scf%variables(:, iter_count)
2890 : scf_env%outer_scf%variables(i, iter_count + 1) = scf_env%outer_scf%variables(i, iter_count) + &
2891 98 : step_multiplier(iwork)*step_size
2892 98 : CALL outer_loop_update_qs_env(qs_env, scf_env)
2893 98 : CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.)
2894 98 : CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot)
2895 : CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
2896 98 : converged=converged, should_stop=should_stop, total_scf_steps=tsteps)
2897 98 : CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft)
2898 : ! Update (iter_count + 1) element of gradient and print constraint info
2899 98 : scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count + 1
2900 98 : CALL outer_loop_gradient(qs_env, scf_env)
2901 : CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
2902 : energy_qs, cdft_control%total_steps, &
2903 98 : should_stop=.FALSE., outer_loop_converged=.FALSE., cdft_loop=.FALSE.)
2904 98 : scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count - 1
2905 : ! Update Jacobian
2906 244 : DO j = 1, nvar
2907 244 : jacobian(j, i) = jacobian(j, i) + coeff(iwork)*scf_env%outer_scf%gradient(j, iter_count + 1)
2908 : END DO
2909 : ! Reset everything to last converged state
2910 244 : scf_env%outer_scf%variables(:, iter_count + 1) = 0.0_dp
2911 2026 : scf_env%outer_scf%gradient = gradient
2912 878 : scf_env%outer_scf%energy = energy
2913 98 : cdft_control%total_steps = nsteps
2914 294 : DO ispin = 1, nspins
2915 196 : CALL deallocate_mo_set(mos(ispin))
2916 196 : CALL duplicate_mo_set(mos(ispin), mos_stashed(ispin))
2917 : CALL calculate_density_matrix(mos(ispin), &
2918 294 : p_rmpv(ispin)%matrix)
2919 : END DO
2920 98 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
2921 368 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
2922 : END DO
2923 : END DO
2924 62 : CALL cp_rm_default_logger()
2925 62 : CALL cp_logger_release(tmp_logger)
2926 : ! Finalize and invert Jacobian
2927 148 : DO j = 1, nvar
2928 282 : DO i = 1, nvar
2929 220 : jacobian(i, j) = jacobian(i, j)/dh(j)
2930 : END DO
2931 : END DO
2932 62 : IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
2933 102 : ALLOCATE (scf_env%outer_scf%inv_jacobian(nvar, nvar))
2934 : END IF
2935 62 : inv_jacobian => scf_env%outer_scf%inv_jacobian
2936 62 : CALL invert_matrix(jacobian, inv_jacobian, inv_error)
2937 62 : scf_control%outer_scf%cdft_opt_control%broyden_update = .FALSE.
2938 : ! Release temporary storage
2939 186 : DO ispin = 1, nspins
2940 186 : CALL deallocate_mo_set(mos_stashed(ispin))
2941 : END DO
2942 62 : DEALLOCATE (mos_stashed, jacobian, gradient, energy, coeff, step_multiplier, dh)
2943 186 : IF (output_unit > 0) THEN
2944 : WRITE (output_unit, FMT="(/,A)") &
2945 31 : " ================================== JACOBIAN CALCULATED =================================="
2946 31 : CALL close_file(unit_number=output_unit)
2947 : END IF
2948 : ELSE
2949 : ! Build a strictly diagonal Jacobian from history and invert it
2950 16 : CALL build_diagonal_jacobian(qs_env, used_history)
2951 : END IF
2952 : END IF
2953 118 : IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian) .AND. para_env%is_source()) THEN
2954 : ! Write restart file for inverse Jacobian
2955 55 : CALL print_inverse_jacobian(logger, scf_env%outer_scf%inv_jacobian, iter_count)
2956 : END IF
2957 : ! Update counter
2958 118 : scf_control%outer_scf%cdft_opt_control%ijacobian(1) = scf_control%outer_scf%cdft_opt_control%ijacobian(1) + 1
2959 118 : CALL timestop(handle)
2960 :
2961 472 : END SUBROUTINE qs_calculate_inverse_jacobian
2962 :
2963 : ! **************************************************************************************************
2964 : !> \brief Perform backtracking line search to find the optimal step size for the CDFT constraint
2965 : !> optimizer. Assumes that the CDFT gradient function is a smooth function of the constraint
2966 : !> variables.
2967 : !> \param qs_env the qs_environment_type where to perform the line search
2968 : !> \par History
2969 : !> 02.2017 created [Nico Holmberg]
2970 : ! **************************************************************************************************
2971 236 : SUBROUTINE qs_cdft_line_search(qs_env)
2972 : TYPE(qs_environment_type), POINTER :: qs_env
2973 :
2974 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_cdft_line_search'
2975 :
2976 : CHARACTER(len=default_path_length) :: project_name
2977 : INTEGER :: handle, i, ispin, iter_count, &
2978 : max_linesearch, max_scf, nspins, &
2979 : nsteps, nvar, output_unit, tsteps
2980 : LOGICAL :: continue_ls, continue_ls_exit, converged, do_linesearch, found_solution, &
2981 : reached_maxls, should_exit, should_stop, sign_changed
2982 236 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: positive_sign
2983 : REAL(KIND=dp) :: alpha, alpha_ls, factor, norm_ls
2984 236 : REAL(KIND=dp), DIMENSION(:), POINTER :: energy
2985 236 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient, inv_jacobian
2986 : REAL(KIND=dp), EXTERNAL :: dnrm2
2987 : TYPE(cdft_control_type), POINTER :: cdft_control
2988 : TYPE(cp_logger_type), POINTER :: logger, tmp_logger
2989 236 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv
2990 236 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
2991 : TYPE(dft_control_type), POINTER :: dft_control
2992 236 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2993 : TYPE(mp_para_env_type), POINTER :: para_env
2994 : TYPE(qs_energy_type), POINTER :: energy_qs
2995 : TYPE(qs_ks_env_type), POINTER :: ks_env
2996 : TYPE(qs_rho_type), POINTER :: rho
2997 : TYPE(qs_scf_env_type), POINTER :: scf_env
2998 : TYPE(scf_control_type), POINTER :: scf_control
2999 :
3000 236 : CALL timeset(routineN, handle)
3001 :
3002 236 : NULLIFY (energy, gradient, p_rmpv, rho_ao_kp, mos, rho, &
3003 236 : ks_env, scf_env, scf_control, dft_control, &
3004 236 : cdft_control, inv_jacobian, para_env, &
3005 236 : tmp_logger, energy_qs)
3006 236 : logger => cp_get_default_logger()
3007 :
3008 236 : CPASSERT(ASSOCIATED(qs_env))
3009 : CALL get_qs_env(qs_env, scf_env=scf_env, ks_env=ks_env, &
3010 : scf_control=scf_control, mos=mos, rho=rho, &
3011 : dft_control=dft_control, &
3012 236 : para_env=para_env, energy=energy_qs)
3013 236 : do_linesearch = .FALSE.
3014 236 : SELECT CASE (scf_control%outer_scf%optimizer)
3015 : CASE DEFAULT
3016 : do_linesearch = .FALSE.
3017 : CASE (outer_scf_optimizer_newton_ls)
3018 24 : do_linesearch = .TRUE.
3019 : CASE (outer_scf_optimizer_broyden)
3020 236 : SELECT CASE (scf_control%outer_scf%cdft_opt_control%broyden_type)
3021 : CASE (broyden_type_1, broyden_type_2, broyden_type_1_explicit, broyden_type_2_explicit)
3022 0 : do_linesearch = .FALSE.
3023 : CASE (broyden_type_1_ls, broyden_type_1_explicit_ls, broyden_type_2_ls, broyden_type_2_explicit_ls)
3024 0 : cdft_control => dft_control%qs_control%cdft_control
3025 0 : IF (.NOT. ASSOCIATED(cdft_control)) THEN
3026 : CALL cp_abort(__LOCATION__, &
3027 : "Optimizers that perform a line search can"// &
3028 0 : " only be used together with a valid CDFT constraint")
3029 : END IF
3030 24 : IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
3031 : do_linesearch = .TRUE.
3032 : END IF
3033 : END SELECT
3034 : END SELECT
3035 : IF (do_linesearch) THEN
3036 8 : BLOCK
3037 8 : TYPE(mo_set_type), DIMENSION(:), ALLOCATABLE :: mos_ls, mos_stashed
3038 8 : cdft_control => dft_control%qs_control%cdft_control
3039 8 : IF (.NOT. ASSOCIATED(cdft_control)) THEN
3040 : CALL cp_abort(__LOCATION__, &
3041 : "Optimizers that perform a line search can"// &
3042 0 : " only be used together with a valid CDFT constraint")
3043 : END IF
3044 8 : CPASSERT(ASSOCIATED(scf_env%outer_scf%inv_jacobian))
3045 8 : CPASSERT(ASSOCIATED(scf_control%outer_scf%cdft_opt_control))
3046 8 : alpha = scf_control%outer_scf%cdft_opt_control%newton_step_save
3047 8 : iter_count = scf_env%outer_scf%iter_count
3048 : ! Redirect output from line search procedure to a new file by creating a temporary logger
3049 8 : project_name = logger%iter_info%project_name
3050 8 : CALL create_tmp_logger(para_env, project_name, "-LineSearch.out", output_unit, tmp_logger)
3051 : ! Save last converged state so we can roll back to it (mo_coeff and some outer_loop variables)
3052 8 : nspins = dft_control%nspins
3053 40 : ALLOCATE (mos_stashed(nspins))
3054 24 : DO ispin = 1, nspins
3055 24 : CALL duplicate_mo_set(mos_stashed(ispin), mos(ispin))
3056 : END DO
3057 8 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
3058 8 : p_rmpv => rho_ao_kp(:, 1)
3059 8 : nsteps = cdft_control%total_steps
3060 : ! Allocate work
3061 8 : nvar = SIZE(scf_env%outer_scf%variables, 1)
3062 8 : max_scf = scf_control%outer_scf%max_scf + 1
3063 8 : max_linesearch = scf_control%outer_scf%cdft_opt_control%max_ls
3064 8 : continue_ls = scf_control%outer_scf%cdft_opt_control%continue_ls
3065 8 : factor = scf_control%outer_scf%cdft_opt_control%factor_ls
3066 8 : continue_ls_exit = .FALSE.
3067 8 : found_solution = .FALSE.
3068 32 : ALLOCATE (gradient(nvar, max_scf))
3069 104 : gradient = scf_env%outer_scf%gradient
3070 24 : ALLOCATE (energy(max_scf))
3071 56 : energy = scf_env%outer_scf%energy
3072 8 : reached_maxls = .FALSE.
3073 : ! Broyden optimizers: perform update of inv_jacobian if necessary
3074 8 : IF (scf_control%outer_scf%cdft_opt_control%broyden_update) THEN
3075 0 : CALL outer_loop_optimize(scf_env, scf_control)
3076 : ! Reset the variables and prevent a reupdate of inv_jacobian
3077 0 : scf_env%outer_scf%variables(:, iter_count + 1) = 0
3078 0 : scf_control%outer_scf%cdft_opt_control%broyden_update = .FALSE.
3079 : END IF
3080 : ! Print some info
3081 8 : IF (output_unit > 0) THEN
3082 : WRITE (output_unit, FMT="(/,A)") &
3083 4 : " ================================== LINE SEARCH STARTED =================================="
3084 : WRITE (output_unit, FMT="(A,I5,A)") &
3085 4 : " Evaluating optimal step size for optimizer using a maximum of", max_linesearch, " steps"
3086 4 : IF (continue_ls) THEN
3087 : WRITE (output_unit, FMT="(A)") &
3088 2 : " Line search continues until best step size is found or max steps are reached"
3089 : END IF
3090 : WRITE (output_unit, '(/,A,F5.3)') &
3091 4 : " Initial step size: ", alpha
3092 : WRITE (output_unit, '(/,A,F5.3)') &
3093 4 : " Step size update factor: ", factor
3094 : WRITE (output_unit, '(/,A,I10,A,I10)') &
3095 4 : " Energy evaluation: ", cdft_control%ienergy, ", CDFT SCF iteration: ", iter_count
3096 : END IF
3097 : ! Perform backtracking line search
3098 8 : CALL cp_add_default_logger(tmp_logger)
3099 16 : DO i = 1, max_linesearch
3100 16 : IF (output_unit > 0) THEN
3101 8 : WRITE (output_unit, FMT="(A)") " "
3102 8 : WRITE (output_unit, FMT="(A)") " #####################################"
3103 : WRITE (output_unit, '(A,I10,A)') &
3104 8 : " ### Line search step: ", i, " ###"
3105 8 : WRITE (output_unit, FMT="(A)") " #####################################"
3106 : END IF
3107 16 : inv_jacobian => scf_env%outer_scf%inv_jacobian
3108 : ! Newton update of CDFT variables with a step size of alpha
3109 : scf_env%outer_scf%variables(:, iter_count + 1) = scf_env%outer_scf%variables(:, iter_count) - alpha* &
3110 128 : MATMUL(inv_jacobian, scf_env%outer_scf%gradient(:, iter_count))
3111 : ! With updated CDFT variables, perform SCF
3112 16 : CALL outer_loop_update_qs_env(qs_env, scf_env)
3113 16 : CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.)
3114 16 : CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot)
3115 : CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
3116 16 : converged=converged, should_stop=should_stop, total_scf_steps=tsteps)
3117 16 : CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft)
3118 : ! Update (iter_count + 1) element of gradient and print constraint info
3119 16 : scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count + 1
3120 16 : CALL outer_loop_gradient(qs_env, scf_env)
3121 : CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
3122 : energy_qs, cdft_control%total_steps, &
3123 16 : should_stop=.FALSE., outer_loop_converged=.FALSE., cdft_loop=.FALSE.)
3124 16 : scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count - 1
3125 : ! Store sign of initial gradient for each variable for continue_ls
3126 16 : IF (continue_ls .AND. .NOT. ALLOCATED(positive_sign)) THEN
3127 12 : ALLOCATE (positive_sign(nvar))
3128 8 : DO ispin = 1, nvar
3129 8 : positive_sign(ispin) = scf_env%outer_scf%gradient(ispin, iter_count + 1) >= 0.0_dp
3130 : END DO
3131 : END IF
3132 : ! Check if the L2 norm of the gradient decreased
3133 16 : inv_jacobian => scf_env%outer_scf%inv_jacobian
3134 16 : IF (dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count + 1), 1) < &
3135 : dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count), 1)) THEN
3136 : ! Optimal step size found
3137 14 : IF (.NOT. continue_ls) THEN
3138 : should_exit = .TRUE.
3139 : ELSE
3140 : ! But line search continues for at least one more iteration in an attempt to find a better solution
3141 : ! if max number of steps is not exceeded
3142 10 : IF (found_solution) THEN
3143 : ! Check if the norm also decreased w.r.t. to previously found solution
3144 6 : IF (dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count + 1), 1) > norm_ls) THEN
3145 : ! Norm increased => accept previous solution and exit
3146 : continue_ls_exit = .TRUE.
3147 : END IF
3148 : END IF
3149 : ! Store current state and the value of alpha
3150 10 : IF (.NOT. continue_ls_exit) THEN
3151 10 : should_exit = .FALSE.
3152 10 : alpha_ls = alpha
3153 10 : found_solution = .TRUE.
3154 10 : norm_ls = dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count + 1), 1)
3155 : ! Check if the sign of the gradient has changed for all variables (w.r.t initial gradient)
3156 : ! In this case we should exit because further line search steps will just increase the norm
3157 10 : sign_changed = .TRUE.
3158 20 : DO ispin = 1, nvar
3159 : sign_changed = sign_changed .AND. (positive_sign(ispin) .NEQV. &
3160 28 : scf_env%outer_scf%gradient(ispin, iter_count + 1) >= 0.0_dp)
3161 : END DO
3162 10 : IF (.NOT. ALLOCATED(mos_ls)) THEN
3163 16 : ALLOCATE (mos_ls(nspins))
3164 : ELSE
3165 18 : DO ispin = 1, nspins
3166 18 : CALL deallocate_mo_set(mos_ls(ispin))
3167 : END DO
3168 : END IF
3169 30 : DO ispin = 1, nspins
3170 30 : CALL duplicate_mo_set(mos_ls(ispin), mos(ispin))
3171 : END DO
3172 10 : alpha = alpha*factor
3173 : ! Exit on last iteration
3174 10 : IF (i == max_linesearch) continue_ls_exit = .TRUE.
3175 : ! Exit if constraint target is satisfied to requested tolerance
3176 20 : IF (SQRT(MAXVAL(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count + 1)**2)) < &
3177 : scf_control%outer_scf%eps_scf) THEN
3178 2 : continue_ls_exit = .TRUE.
3179 : END IF
3180 : ! Exit if line search jumped over the optimal step length
3181 10 : IF (sign_changed) continue_ls_exit = .TRUE.
3182 : END IF
3183 : END IF
3184 : ELSE
3185 : ! Gradient increased => alpha is too large (if the gradient function is smooth)
3186 2 : should_exit = .FALSE.
3187 : ! Update alpha using Armijo's scheme
3188 2 : alpha = alpha*factor
3189 : END IF
3190 14 : IF (continue_ls_exit) THEN
3191 : ! Continuation of line search did not yield a better alpha, use previously located solution and exit
3192 4 : alpha = alpha_ls
3193 12 : DO ispin = 1, nspins
3194 8 : CALL deallocate_mo_set(mos(ispin))
3195 8 : CALL duplicate_mo_set(mos(ispin), mos_ls(ispin))
3196 : CALL calculate_density_matrix(mos(ispin), &
3197 8 : p_rmpv(ispin)%matrix)
3198 12 : CALL deallocate_mo_set(mos_ls(ispin))
3199 : END DO
3200 4 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
3201 4 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
3202 4 : DEALLOCATE (mos_ls)
3203 : should_exit = .TRUE.
3204 : END IF
3205 : ! Reached max steps and SCF converged: continue with last iterated step size
3206 12 : IF (.NOT. should_exit .AND. &
3207 : (i == max_linesearch .AND. converged .AND. .NOT. found_solution)) THEN
3208 0 : should_exit = .TRUE.
3209 0 : reached_maxls = .TRUE.
3210 0 : alpha = alpha*(1.0_dp/factor)
3211 : END IF
3212 : ! Reset outer SCF environment to last converged state
3213 32 : scf_env%outer_scf%variables(:, iter_count + 1) = 0.0_dp
3214 208 : scf_env%outer_scf%gradient = gradient
3215 112 : scf_env%outer_scf%energy = energy
3216 : ! Exit line search if a suitable step size was found
3217 16 : IF (should_exit) EXIT
3218 : ! Reset the electronic structure
3219 8 : cdft_control%total_steps = nsteps
3220 24 : DO ispin = 1, nspins
3221 16 : CALL deallocate_mo_set(mos(ispin))
3222 16 : CALL duplicate_mo_set(mos(ispin), mos_stashed(ispin))
3223 : CALL calculate_density_matrix(mos(ispin), &
3224 24 : p_rmpv(ispin)%matrix)
3225 : END DO
3226 8 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
3227 24 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
3228 : END DO
3229 8 : scf_control%outer_scf%cdft_opt_control%newton_step = alpha
3230 8 : IF (.NOT. should_exit) THEN
3231 : CALL cp_warn(__LOCATION__, &
3232 0 : "Line search did not converge. CDFT SCF proceeds with fixed step size.")
3233 0 : scf_control%outer_scf%cdft_opt_control%newton_step = scf_control%outer_scf%cdft_opt_control%newton_step_save
3234 : END IF
3235 8 : IF (reached_maxls) THEN
3236 : CALL cp_warn(__LOCATION__, &
3237 0 : "Line search did not converge. CDFT SCF proceeds with lasted iterated step size.")
3238 : END IF
3239 8 : CALL cp_rm_default_logger()
3240 8 : CALL cp_logger_release(tmp_logger)
3241 : ! Release temporary storage
3242 24 : DO ispin = 1, nspins
3243 24 : CALL deallocate_mo_set(mos_stashed(ispin))
3244 : END DO
3245 8 : DEALLOCATE (mos_stashed, gradient, energy)
3246 8 : IF (ALLOCATED(positive_sign)) DEALLOCATE (positive_sign)
3247 20 : IF (output_unit > 0) THEN
3248 : WRITE (output_unit, FMT="(/,A)") &
3249 4 : " ================================== LINE SEARCH COMPLETE =================================="
3250 4 : CALL close_file(unit_number=output_unit)
3251 : END IF
3252 : END BLOCK
3253 : END IF
3254 :
3255 236 : CALL timestop(handle)
3256 :
3257 236 : END SUBROUTINE qs_cdft_line_search
3258 :
3259 16 : END MODULE qs_scf
|
