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