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