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