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