Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines for the Quickstep SCF run.
10 : !> \par History
11 : !> - Joost VandeVondele (02.2002)
12 : !> added code for: incremental (pab and gvg) update
13 : !> initialisation (init_cube, l_info)
14 : !> - Joost VandeVondele (02.2002)
15 : !> called the poisson code of the classical part
16 : !> this takes into account the spherical cutoff and allows for
17 : !> isolated systems
18 : !> - Joost VandeVondele (02.2002)
19 : !> added multiple grid feature
20 : !> changed to spherical cutoff consistently (?)
21 : !> therefore removed the gradient correct functionals
22 : !> - updated with the new QS data structures (10.04.02,MK)
23 : !> - copy_matrix replaced by transfer_matrix (11.04.02,MK)
24 : !> - nrebuild_rho and nrebuild_gvg unified (12.04.02,MK)
25 : !> - set_mo_occupation for smearing of the MO occupation numbers
26 : !> (17.04.02,MK)
27 : !> - MO level shifting added (22.04.02,MK)
28 : !> - Usage of TYPE mo_set_p_type
29 : !> - Joost VandeVondele (05.2002)
30 : !> added cholesky based diagonalisation
31 : !> - 05.2002 added pao method [fawzi]
32 : !> - parallel FFT (JGH 22.05.2002)
33 : !> - 06.2002 moved KS matrix construction to qs_build_KS_matrix.F [fawzi]
34 : !> - started to include more LSD (01.2003,Joost VandeVondele)
35 : !> - 02.2003 scf_env [fawzi]
36 : !> - got rid of nrebuild (01.2004, Joost VandeVondele)
37 : !> - 10.2004 removed pao [fawzi]
38 : !> - 03.2006 large cleaning action [Joost VandeVondele]
39 : !> - High-spin ROKS added (05.04.06,MK)
40 : !> - Mandes (10.2013)
41 : !> intermediate energy communication with external communicator added
42 : !> - kpoints (08.2014, JGH)
43 : !> - unified k-point and gamma-point code (2014.11) [Ole Schuett]
44 : !> - added extra SCF loop for CDFT constraints (12.2015) [Nico Holmberg]
45 : !> \author Matthias Krack (30.04.2001)
46 : ! **************************************************************************************************
47 : MODULE qs_scf
48 : USE atomic_kind_types, ONLY: atomic_kind_type
49 : USE cp_control_types, ONLY: dft_control_type
50 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
51 : dbcsr_deallocate_matrix,&
52 : dbcsr_get_info,&
53 : dbcsr_init_p,&
54 : dbcsr_p_type,&
55 : dbcsr_set,&
56 : dbcsr_type
57 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
58 : dbcsr_deallocate_matrix_set
59 : USE cp_files, ONLY: close_file
60 : USE cp_fm_types, ONLY: cp_fm_create,&
61 : cp_fm_release,&
62 : cp_fm_to_fm,&
63 : cp_fm_type
64 : USE cp_log_handling, ONLY: cp_add_default_logger,&
65 : cp_get_default_logger,&
66 : cp_logger_release,&
67 : cp_logger_type,&
68 : cp_rm_default_logger,&
69 : cp_to_string
70 : USE cp_output_handling, ONLY: cp_add_iter_level,&
71 : cp_iterate,&
72 : cp_p_file,&
73 : cp_print_key_should_output,&
74 : cp_print_key_unit_nr,&
75 : cp_rm_iter_level
76 : USE cp_result_methods, ONLY: get_results,&
77 : test_for_result
78 : USE cp_result_types, ONLY: cp_result_type
79 : USE ec_env_types, ONLY: energy_correction_type
80 : USE input_constants, ONLY: &
81 : broyden_type_1, broyden_type_1_explicit, broyden_type_1_explicit_ls, broyden_type_1_ls, &
82 : broyden_type_2, broyden_type_2_explicit, broyden_type_2_explicit_ls, broyden_type_2_ls, &
83 : cdft2ot, history_guess, ot2cdft, ot_precond_full_all, ot_precond_full_single, &
84 : ot_precond_full_single_inverse, ot_precond_none, ot_precond_s_inverse, &
85 : outer_scf_becke_constraint, outer_scf_hirshfeld_constraint, outer_scf_optimizer_broyden, &
86 : outer_scf_optimizer_newton_ls
87 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
88 : section_vals_type
89 : USE kinds, ONLY: default_path_length,&
90 : default_string_length,&
91 : dp
92 : USE kpoint_io, ONLY: write_kpoints_restart
93 : USE kpoint_types, ONLY: kpoint_type
94 : USE machine, ONLY: m_flush,&
95 : m_walltime
96 : USE mathlib, ONLY: invert_matrix
97 : USE message_passing, ONLY: mp_comm_type,&
98 : mp_para_env_type
99 : USE particle_types, ONLY: particle_type
100 : USE preconditioner, ONLY: prepare_preconditioner,&
101 : restart_preconditioner
102 : USE pw_env_types, ONLY: pw_env_get,&
103 : pw_env_type
104 : USE pw_pool_types, ONLY: pw_pool_type
105 : USE qs_block_davidson_types, ONLY: block_davidson_deallocate
106 : USE qs_cdft_scf_utils, ONLY: build_diagonal_jacobian,&
107 : create_tmp_logger,&
108 : initialize_inverse_jacobian,&
109 : prepare_jacobian_stencil,&
110 : print_inverse_jacobian,&
111 : restart_inverse_jacobian
112 : USE qs_cdft_types, ONLY: cdft_control_type
113 : USE qs_charges_types, ONLY: qs_charges_type
114 : USE qs_density_matrices, ONLY: calculate_density_matrix
115 : USE qs_density_mixing_types, ONLY: gspace_mixing_nr
116 : USE qs_diis, ONLY: qs_diis_b_clear,&
117 : qs_diis_b_clear_kp,&
118 : qs_diis_b_create,&
119 : qs_diis_b_create_kp
120 : USE qs_energy_types, ONLY: qs_energy_type
121 : USE qs_environment_types, ONLY: get_qs_env,&
122 : qs_environment_type,&
123 : set_qs_env
124 : USE qs_integrate_potential, ONLY: integrate_v_rspace
125 : USE qs_kind_types, ONLY: qs_kind_type
126 : USE qs_ks_methods, ONLY: evaluate_core_matrix_traces,&
127 : qs_ks_update_qs_env
128 : USE qs_ks_types, ONLY: get_ks_env,&
129 : qs_ks_did_change,&
130 : qs_ks_env_type
131 : USE qs_mo_io, ONLY: write_mo_set_to_restart
132 : USE qs_mo_methods, ONLY: make_basis_simple,&
133 : make_basis_sm
134 : USE qs_mo_occupation, ONLY: set_mo_occupation
135 : USE qs_mo_types, ONLY: deallocate_mo_set,&
136 : duplicate_mo_set,&
137 : get_mo_set,&
138 : mo_set_type,&
139 : reassign_allocated_mos
140 : USE qs_ot, ONLY: qs_ot_new_preconditioner
141 : USE qs_ot_scf, ONLY: ot_scf_init,&
142 : ot_scf_read_input
143 : USE qs_outer_scf, ONLY: outer_loop_gradient,&
144 : outer_loop_optimize,&
145 : outer_loop_purge_history,&
146 : outer_loop_switch,&
147 : outer_loop_update_qs_env
148 : USE qs_rho_methods, ONLY: qs_rho_update_rho
149 : USE qs_rho_types, ONLY: qs_rho_get,&
150 : qs_rho_type
151 : USE qs_scf_initialization, ONLY: qs_scf_env_initialize
152 : USE qs_scf_loop_utils, ONLY: qs_scf_check_inner_exit,&
153 : qs_scf_check_outer_exit,&
154 : qs_scf_density_mixing,&
155 : qs_scf_inner_finalize,&
156 : qs_scf_new_mos,&
157 : qs_scf_new_mos_kp,&
158 : qs_scf_rho_update,&
159 : qs_scf_set_loop_flags
160 : USE qs_scf_output, ONLY: qs_scf_cdft_info,&
161 : qs_scf_cdft_initial_info,&
162 : qs_scf_loop_info,&
163 : qs_scf_loop_print,&
164 : qs_scf_outer_loop_info,&
165 : qs_scf_write_mos
166 : USE qs_scf_post_scf, ONLY: qs_scf_compute_properties
167 : USE qs_scf_types, ONLY: &
168 : block_davidson_diag_method_nr, block_krylov_diag_method_nr, filter_matrix_diag_method_nr, &
169 : general_diag_method_nr, ot_diag_method_nr, ot_method_nr, qs_scf_env_type, &
170 : smeagol_method_nr, special_diag_method_nr
171 : USE qs_wf_history_methods, ONLY: wfi_purge_history,&
172 : wfi_update
173 : USE scf_control_types, ONLY: scf_control_type
174 : USE smeagol_interface, ONLY: run_smeagol_bulktrans,&
175 : run_smeagol_emtrans
176 : USE tblite_interface, ONLY: tb_get_energy,&
177 : tb_update_charges
178 : #include "./base/base_uses.f90"
179 :
180 : IMPLICIT NONE
181 :
182 : PRIVATE
183 :
184 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf'
185 : LOGICAL, PRIVATE :: reuse_precond = .FALSE.
186 : LOGICAL, PRIVATE :: used_history = .FALSE.
187 :
188 : PUBLIC :: scf, scf_env_cleanup, scf_env_do_scf, cdft_scf, init_scf_loop
189 :
190 : CONTAINS
191 :
192 : ! **************************************************************************************************
193 : !> \brief perform an scf procedure in the given qs_env
194 : !> \param qs_env the qs_environment where to perform the scf procedure
195 : !> \param has_converged ...
196 : !> \param total_scf_steps ...
197 : !> \par History
198 : !> 02.2003 introduced scf_env, moved real work to scf_env_do_scf [fawzi]
199 : !> \author fawzi
200 : !> \note
201 : ! **************************************************************************************************
202 19625 : SUBROUTINE scf(qs_env, has_converged, total_scf_steps)
203 : TYPE(qs_environment_type), POINTER :: qs_env
204 : LOGICAL, INTENT(OUT), OPTIONAL :: has_converged
205 : INTEGER, INTENT(OUT), OPTIONAL :: total_scf_steps
206 :
207 : INTEGER :: ihistory, max_scf_tmp, tsteps
208 : LOGICAL :: converged, outer_scf_loop, should_stop
209 : LOGICAL, SAVE :: first_step_flag = .TRUE.
210 19625 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient_history, variable_history
211 : TYPE(cp_logger_type), POINTER :: logger
212 : TYPE(dft_control_type), POINTER :: dft_control
213 : TYPE(qs_scf_env_type), POINTER :: scf_env
214 : TYPE(scf_control_type), POINTER :: scf_control
215 : TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
216 :
217 19625 : NULLIFY (scf_env)
218 19625 : logger => cp_get_default_logger()
219 19625 : CPASSERT(ASSOCIATED(qs_env))
220 19625 : IF (PRESENT(has_converged)) THEN
221 0 : has_converged = .FALSE.
222 : END IF
223 19625 : IF (PRESENT(total_scf_steps)) THEN
224 0 : total_scf_steps = 0
225 : END IF
226 : CALL get_qs_env(qs_env, scf_env=scf_env, input=input, &
227 19625 : dft_control=dft_control, scf_control=scf_control)
228 19625 : IF (scf_control%max_scf > 0) THEN
229 :
230 18983 : dft_section => section_vals_get_subs_vals(input, "DFT")
231 18983 : scf_section => section_vals_get_subs_vals(dft_section, "SCF")
232 :
233 18983 : IF (.NOT. ASSOCIATED(scf_env)) THEN
234 5793 : CALL qs_scf_env_initialize(qs_env, scf_env)
235 : ! Moved here from qs_scf_env_initialize to be able to have more scf_env
236 5793 : CALL set_qs_env(qs_env, scf_env=scf_env)
237 : ELSE
238 13190 : CALL qs_scf_env_initialize(qs_env, scf_env)
239 : END IF
240 :
241 18983 : IF ((scf_control%density_guess == history_guess) .AND. (first_step_flag)) THEN
242 2 : max_scf_tmp = scf_control%max_scf
243 2 : scf_control%max_scf = 1
244 2 : outer_scf_loop = scf_control%outer_scf%have_scf
245 2 : scf_control%outer_scf%have_scf = .FALSE.
246 : END IF
247 :
248 18983 : IF (.NOT. dft_control%qs_control%cdft) THEN
249 : CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
250 18657 : converged=converged, should_stop=should_stop, total_scf_steps=tsteps)
251 : ELSE
252 : ! Third SCF loop needed for CDFT with OT to properly restart OT inner loop
253 326 : CALL cdft_scf(qs_env=qs_env, should_stop=should_stop)
254 : END IF
255 :
256 : ! If SCF has not converged, then we should not start MP2
257 18983 : IF (ASSOCIATED(qs_env%mp2_env)) qs_env%mp2_env%hf_fail = .NOT. converged
258 :
259 : ! Add the converged outer_scf SCF gradient(s)/variable(s) to history
260 18983 : IF (scf_control%outer_scf%have_scf) THEN
261 4135 : ihistory = scf_env%outer_scf%iter_count
262 : CALL get_qs_env(qs_env, gradient_history=gradient_history, &
263 4135 : variable_history=variable_history)
264 : ! We only store the latest two values
265 8300 : gradient_history(:, 1) = gradient_history(:, 2)
266 16600 : gradient_history(:, 2) = scf_env%outer_scf%gradient(:, ihistory)
267 8300 : variable_history(:, 1) = variable_history(:, 2)
268 16600 : variable_history(:, 2) = scf_env%outer_scf%variables(:, ihistory)
269 : ! Reset flag
270 4135 : IF (used_history) used_history = .FALSE.
271 : ! Update a counter and check if the Jacobian should be deallocated
272 4135 : IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
273 64 : scf_control%outer_scf%cdft_opt_control%ijacobian(2) = scf_control%outer_scf%cdft_opt_control%ijacobian(2) + 1
274 : IF (scf_control%outer_scf%cdft_opt_control%ijacobian(2) >= &
275 64 : scf_control%outer_scf%cdft_opt_control%jacobian_freq(2) .AND. &
276 : scf_control%outer_scf%cdft_opt_control%jacobian_freq(2) > 0) &
277 50 : scf_env%outer_scf%deallocate_jacobian = .TRUE.
278 : END IF
279 : END IF
280 : ! *** add the converged wavefunction to the wavefunction history
281 18983 : IF ((ASSOCIATED(qs_env%wf_history)) .AND. &
282 : ((scf_control%density_guess /= history_guess) .OR. &
283 : (.NOT. first_step_flag))) THEN
284 18981 : IF (.NOT. dft_control%qs_control%cdft) THEN
285 18655 : CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
286 : ELSE
287 326 : IF (dft_control%qs_control%cdft_control%should_purge) THEN
288 0 : CALL wfi_purge_history(qs_env)
289 0 : CALL outer_loop_purge_history(qs_env)
290 0 : dft_control%qs_control%cdft_control%should_purge = .FALSE.
291 : ELSE
292 326 : CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
293 : END IF
294 : END IF
295 2 : ELSE IF ((scf_control%density_guess == history_guess) .AND. &
296 : (first_step_flag)) THEN
297 2 : scf_control%max_scf = max_scf_tmp
298 2 : scf_control%outer_scf%have_scf = outer_scf_loop
299 2 : first_step_flag = .FALSE.
300 : END IF
301 :
302 : ! *** compute properties that depend on the converged wavefunction
303 18983 : IF (.NOT. (should_stop)) CALL qs_scf_compute_properties(qs_env)
304 :
305 : ! *** SMEAGOL interface ***
306 18983 : IF (.NOT. (should_stop)) THEN
307 : ! compute properties that depend on the converged wavefunction ..
308 18983 : CALL run_smeagol_emtrans(qs_env, last=.TRUE., iter=0)
309 : ! .. or save matrices related to bulk leads
310 18983 : CALL run_smeagol_bulktrans(qs_env)
311 : END IF
312 :
313 : ! *** cleanup
314 18983 : CALL scf_env_cleanup(scf_env)
315 18983 : IF (dft_control%qs_control%cdft) &
316 326 : CALL cdft_control_cleanup(dft_control%qs_control%cdft_control)
317 :
318 18983 : IF (PRESENT(has_converged)) THEN
319 0 : has_converged = converged
320 : END IF
321 18983 : IF (PRESENT(total_scf_steps)) THEN
322 0 : total_scf_steps = tsteps
323 : END IF
324 :
325 : END IF
326 :
327 19625 : END SUBROUTINE scf
328 :
329 : ! **************************************************************************************************
330 : !> \brief perform an scf loop
331 : !> \param scf_env the scf_env where to perform the scf procedure
332 : !> \param scf_control ...
333 : !> \param qs_env the qs_env, the scf_env lives in
334 : !> \param converged will be true / false if converged is reached
335 : !> \param should_stop ...
336 : !> \param total_scf_steps ...
337 : !> \par History
338 : !> long history, see cvs and qs_scf module history
339 : !> 02.2003 introduced scf_env [fawzi]
340 : !> 09.2005 Frozen density approximation [TdK]
341 : !> 06.2007 Check for SCF iteration count early [jgh]
342 : !> 10.2019 switch_surf_dip [SGh]
343 : !> \author Matthias Krack
344 : !> \note
345 : ! **************************************************************************************************
346 19285 : SUBROUTINE scf_env_do_scf(scf_env, scf_control, qs_env, converged, should_stop, total_scf_steps)
347 :
348 : TYPE(qs_scf_env_type), POINTER :: scf_env
349 : TYPE(scf_control_type), POINTER :: scf_control
350 : TYPE(qs_environment_type), POINTER :: qs_env
351 : LOGICAL, INTENT(OUT) :: converged, should_stop
352 : INTEGER, INTENT(OUT) :: total_scf_steps
353 :
354 : CHARACTER(LEN=*), PARAMETER :: routineN = 'scf_env_do_scf'
355 :
356 : CHARACTER(LEN=default_string_length) :: description, name
357 : INTEGER :: ext_master_id, handle, handle2, i_tmp, &
358 : ic, ispin, iter_count, output_unit, &
359 : scf_energy_message_tag, total_steps
360 : LOGICAL :: diis_step, do_kpoints, energy_only, exit_inner_loop, exit_outer_loop, &
361 : inner_loop_converged, just_energy, outer_loop_converged
362 : REAL(KIND=dp) :: t1, t2
363 : REAL(KIND=dp), DIMENSION(3) :: res_val_3
364 19285 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
365 : TYPE(cp_logger_type), POINTER :: logger
366 : TYPE(cp_result_type), POINTER :: results
367 19285 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
368 19285 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
369 : TYPE(dft_control_type), POINTER :: dft_control
370 : TYPE(energy_correction_type), POINTER :: ec_env
371 : TYPE(kpoint_type), POINTER :: kpoints
372 19285 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_last_converged
373 : TYPE(mp_comm_type) :: external_comm
374 : TYPE(mp_para_env_type), POINTER :: para_env
375 19285 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
376 : TYPE(pw_env_type), POINTER :: pw_env
377 : TYPE(qs_charges_type), POINTER :: qs_charges
378 : TYPE(qs_energy_type), POINTER :: energy
379 19285 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
380 : TYPE(qs_ks_env_type), POINTER :: ks_env
381 : TYPE(qs_rho_type), POINTER :: rho
382 : TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
383 :
384 19285 : CALL timeset(routineN, handle)
385 :
386 19285 : NULLIFY (dft_control, rho, energy, &
387 19285 : logger, qs_charges, ks_env, mos, atomic_kind_set, qs_kind_set, &
388 19285 : particle_set, dft_section, input, &
389 19285 : scf_section, para_env, results, kpoints, pw_env, rho_ao_kp, mos_last_converged)
390 :
391 19285 : CPASSERT(ASSOCIATED(scf_env))
392 19285 : CPASSERT(ASSOCIATED(qs_env))
393 :
394 19285 : logger => cp_get_default_logger()
395 19285 : t1 = m_walltime()
396 :
397 : CALL get_qs_env(qs_env=qs_env, &
398 : energy=energy, &
399 : particle_set=particle_set, &
400 : qs_charges=qs_charges, &
401 : ks_env=ks_env, &
402 : atomic_kind_set=atomic_kind_set, &
403 : qs_kind_set=qs_kind_set, &
404 : rho=rho, &
405 : mos=mos, &
406 : input=input, &
407 : dft_control=dft_control, &
408 : do_kpoints=do_kpoints, &
409 : kpoints=kpoints, &
410 : results=results, &
411 : pw_env=pw_env, &
412 19285 : para_env=para_env)
413 :
414 19285 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
415 :
416 19285 : dft_section => section_vals_get_subs_vals(input, "DFT")
417 19285 : scf_section => section_vals_get_subs_vals(dft_section, "SCF")
418 :
419 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", &
420 19285 : extension=".scfLog")
421 :
422 19285 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
423 9825 : "SCF WAVEFUNCTION OPTIMIZATION"
424 :
425 : ! when switch_surf_dip is switched on, indicate storing mos from the last converged step
426 19285 : IF (dft_control%switch_surf_dip) THEN
427 2 : CALL get_qs_env(qs_env, mos_last_converged=mos_last_converged)
428 4 : DO ispin = 1, dft_control%nspins
429 4 : CALL reassign_allocated_mos(mos(ispin), mos_last_converged(ispin))
430 : END DO
431 2 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
432 1 : "COPIED mos_last_converged ---> mos"
433 : END IF
434 :
435 19285 : IF ((output_unit > 0) .AND. (.NOT. scf_control%use_ot)) THEN
436 : WRITE (UNIT=output_unit, &
437 : FMT="(/,T3,A,T12,A,T31,A,T39,A,T59,A,T75,A,/,T3,A)") &
438 6561 : "Step", "Update method", "Time", "Convergence", "Total energy", "Change", &
439 13122 : REPEAT("-", 78)
440 : END IF
441 19285 : CALL cp_add_iter_level(logger%iter_info, "QS_SCF")
442 :
443 : ! check for external communicator and if the intermediate energy should be sent
444 77140 : res_val_3(:) = -1.0_dp
445 19285 : description = "[EXT_SCF_ENER_COMM]"
446 19285 : IF (test_for_result(results, description=description)) THEN
447 : CALL get_results(results, description=description, &
448 0 : values=res_val_3, n_entries=i_tmp)
449 0 : CPASSERT(i_tmp == 3)
450 0 : IF (ALL(res_val_3(:) <= 0.0)) &
451 : CALL cp_abort(__LOCATION__, &
452 : " Trying to access result ("//TRIM(description)// &
453 0 : ") which is not correctly stored.")
454 0 : CALL external_comm%set_handle(NINT(res_val_3(1)))
455 : END IF
456 19285 : ext_master_id = NINT(res_val_3(2))
457 19285 : scf_energy_message_tag = NINT(res_val_3(3))
458 :
459 : ! *** outer loop of the scf, can treat other variables,
460 : ! *** such as lagrangian multipliers
461 19285 : scf_env%outer_scf%iter_count = 0
462 19285 : iter_count = 0
463 19285 : total_steps = 0
464 19285 : energy%tot_old = 0.0_dp
465 :
466 856 : scf_outer_loop: DO
467 :
468 : CALL init_scf_loop(scf_env=scf_env, qs_env=qs_env, &
469 20141 : scf_section=scf_section)
470 :
471 : CALL qs_scf_set_loop_flags(scf_env, diis_step, &
472 20141 : energy_only, just_energy, exit_inner_loop)
473 :
474 : ! decide whether to switch off dipole correction for convergence purposes
475 20141 : dft_control%surf_dip_correct_switch = dft_control%correct_surf_dip
476 20141 : IF ((dft_control%correct_surf_dip) .AND. (scf_control%outer_scf%have_scf) .AND. &
477 : (scf_env%outer_scf%iter_count > FLOOR(scf_control%outer_scf%max_scf/2.0_dp))) THEN
478 0 : IF (dft_control%switch_surf_dip) THEN
479 0 : dft_control%surf_dip_correct_switch = .FALSE.
480 0 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
481 0 : "SURFACE DIPOLE CORRECTION switched off"
482 : END IF
483 : END IF
484 :
485 174495 : scf_loop: DO
486 :
487 174495 : CALL timeset(routineN//"_inner_loop", handle2)
488 :
489 174495 : IF (.NOT. just_energy) scf_env%iter_count = scf_env%iter_count + 1
490 174495 : iter_count = iter_count + 1
491 174495 : CALL cp_iterate(logger%iter_info, last=.FALSE., iter_nr=iter_count)
492 :
493 174495 : IF (output_unit > 0) CALL m_flush(output_unit)
494 :
495 174495 : total_steps = total_steps + 1
496 174495 : just_energy = energy_only
497 :
498 : CALL qs_ks_update_qs_env(qs_env, just_energy=just_energy, &
499 174495 : calculate_forces=.FALSE.)
500 :
501 : ! print 'heavy weight' or relatively expensive quantities
502 174495 : CALL qs_scf_loop_print(qs_env, scf_env, para_env)
503 :
504 174495 : IF (do_kpoints) THEN
505 : ! kpoints
506 5428 : IF (dft_control%hairy_probes .EQV. .TRUE.) THEN
507 0 : scf_control%smear%do_smear = .FALSE.
508 0 : CALL qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step, dft_control%probe)
509 : ELSE
510 5428 : CALL qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step)
511 : END IF
512 : ELSE
513 : ! Gamma points only
514 169067 : IF (dft_control%hairy_probes .EQV. .TRUE.) THEN
515 14 : scf_control%smear%do_smear = .FALSE.
516 : CALL qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, energy_only, &
517 14 : dft_control%probe)
518 : ELSE
519 169053 : CALL qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, energy_only)
520 : END IF
521 : END IF
522 :
523 : ! Print requested MO information (can be computationally expensive with OT)
524 174495 : CALL qs_scf_write_mos(qs_env, scf_env, final_mos=.FALSE.)
525 :
526 174495 : IF (dft_control%qs_control%xtb_control%do_tblite) THEN
527 9674 : CPASSERT(scf_env%mixing_method > 0)
528 9674 : CALL tb_update_charges(qs_env, dft_control, qs_env%tb_tblite, .FALSE., .FALSE.)
529 9674 : CALL evaluate_core_matrix_traces(qs_env, rho_ao_ext=scf_env%p_mix_new)
530 9674 : CALL tb_get_energy(qs_env, qs_env%tb_tblite, energy)
531 : END IF
532 :
533 174495 : CALL qs_scf_density_mixing(scf_env, rho, para_env, diis_step)
534 :
535 174495 : t2 = m_walltime()
536 :
537 174495 : CALL qs_scf_loop_info(scf_env, output_unit, just_energy, t1, t2, energy)
538 :
539 174495 : IF (.NOT. just_energy) energy%tot_old = energy%total
540 :
541 : ! check for external communicator and if the intermediate energy should be sent
542 174495 : IF (scf_energy_message_tag > 0) THEN
543 0 : CALL external_comm%send(energy%total, ext_master_id, scf_energy_message_tag)
544 : END IF
545 :
546 : CALL qs_scf_check_inner_exit(qs_env, scf_env, scf_control, should_stop, just_energy, &
547 174495 : exit_inner_loop, inner_loop_converged, output_unit)
548 :
549 : ! In case we decide to exit we perform few more check to see if this one
550 : ! is really the last SCF step
551 174495 : IF (exit_inner_loop) THEN
552 :
553 20141 : CALL qs_scf_inner_finalize(scf_env, qs_env, diis_step, output_unit)
554 :
555 : CALL qs_scf_check_outer_exit(qs_env, scf_env, scf_control, should_stop, &
556 20141 : outer_loop_converged, exit_outer_loop)
557 :
558 : ! Let's tag the last SCF cycle so we can print informations only of the last step
559 20141 : IF (exit_outer_loop) CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=iter_count)
560 :
561 : END IF
562 :
563 174495 : IF (do_kpoints) THEN
564 5428 : CALL write_kpoints_restart(rho_ao_kp, kpoints, scf_env, dft_section, particle_set, qs_kind_set)
565 : ELSE
566 : ! Write wavefunction restart file
567 169067 : IF (scf_env%method == ot_method_nr) THEN
568 : ! With OT: provide the Kohn-Sham matrix for the calculation of the MO eigenvalues
569 75232 : CALL get_ks_env(ks_env=ks_env, matrix_ks=matrix_ks)
570 : CALL write_mo_set_to_restart(mos, particle_set, dft_section=dft_section, qs_kind_set=qs_kind_set, &
571 75232 : matrix_ks=matrix_ks)
572 : ELSE
573 93835 : CALL write_mo_set_to_restart(mos, particle_set, dft_section=dft_section, qs_kind_set=qs_kind_set)
574 : END IF
575 :
576 : END IF
577 :
578 : ! Exit if we have finished with the SCF inner loop
579 174495 : IF (exit_inner_loop) THEN
580 20141 : CALL timestop(handle2)
581 : EXIT scf_loop
582 : END IF
583 :
584 154354 : IF (.NOT. BTEST(cp_print_key_should_output(logger%iter_info, &
585 : scf_section, "PRINT%ITERATION_INFO/TIME_CUMUL"), cp_p_file)) &
586 154354 : t1 = m_walltime()
587 :
588 : ! mixing methods have the new density matrix in p_mix_new
589 154354 : IF (scf_env%mixing_method > 0) THEN
590 530658 : DO ic = 1, SIZE(rho_ao_kp, 2)
591 1025695 : DO ispin = 1, dft_control%nspins
592 495037 : CALL dbcsr_get_info(rho_ao_kp(ispin, ic)%matrix, name=name) ! keep the name
593 939408 : CALL dbcsr_copy(rho_ao_kp(ispin, ic)%matrix, scf_env%p_mix_new(ispin, ic)%matrix, name=name)
594 : END DO
595 : END DO
596 : END IF
597 :
598 : CALL qs_scf_rho_update(rho, qs_env, scf_env, ks_env, &
599 154354 : mix_rho=scf_env%mixing_method >= gspace_mixing_nr)
600 :
601 154354 : CALL timestop(handle2)
602 :
603 : END DO scf_loop
604 :
605 20141 : IF (.NOT. scf_control%outer_scf%have_scf) EXIT scf_outer_loop
606 :
607 : ! In case we use the OUTER SCF loop let's print some info..
608 : CALL qs_scf_outer_loop_info(output_unit, scf_control, scf_env, &
609 5291 : energy, total_steps, should_stop, outer_loop_converged)
610 :
611 : ! Save MOs to converged MOs if outer_loop_converged and surf_dip_correct_switch is true
612 5291 : IF (exit_outer_loop) THEN
613 4435 : IF ((dft_control%switch_surf_dip) .AND. (outer_loop_converged) .AND. &
614 : (dft_control%surf_dip_correct_switch)) THEN
615 4 : DO ispin = 1, dft_control%nspins
616 4 : CALL reassign_allocated_mos(mos_last_converged(ispin), mos(ispin))
617 : END DO
618 2 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
619 1 : "COPIED mos ---> mos_last_converged"
620 : END IF
621 : END IF
622 :
623 5291 : IF (exit_outer_loop) EXIT scf_outer_loop
624 :
625 : !
626 856 : CALL outer_loop_optimize(scf_env, scf_control)
627 856 : CALL outer_loop_update_qs_env(qs_env, scf_env)
628 20141 : CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.)
629 :
630 : END DO scf_outer_loop
631 :
632 19285 : converged = inner_loop_converged .AND. outer_loop_converged
633 19285 : total_scf_steps = total_steps
634 :
635 19285 : IF (dft_control%qs_control%cdft) &
636 : dft_control%qs_control%cdft_control%total_steps = &
637 626 : dft_control%qs_control%cdft_control%total_steps + total_steps
638 :
639 19285 : IF (.NOT. converged) THEN
640 2170 : IF (scf_control%ignore_convergence_failure .OR. should_stop) THEN
641 2170 : CALL cp_warn(__LOCATION__, "SCF run NOT converged")
642 : ELSE
643 : CALL cp_abort(__LOCATION__, &
644 : "SCF run NOT converged. To continue the calculation "// &
645 0 : "regardless, please set the keyword IGNORE_CONVERGENCE_FAILURE.")
646 : END IF
647 : END IF
648 :
649 : ! Skip Harris functional calculation if ground-state is NOT converged
650 19285 : IF (qs_env%energy_correction) THEN
651 656 : CALL get_qs_env(qs_env, ec_env=ec_env)
652 656 : ec_env%do_skip = .FALSE.
653 656 : IF (ec_env%skip_ec .AND. .NOT. converged) ec_env%do_skip = .TRUE.
654 : END IF
655 :
656 : ! if needed copy mo_coeff dbcsr->fm for later use in post_scf!fm->dbcsr
657 41130 : DO ispin = 1, SIZE(mos) !fm -> dbcsr
658 41130 : IF (mos(ispin)%use_mo_coeff_b) THEN !fm->dbcsr
659 7511 : IF (.NOT. ASSOCIATED(mos(ispin)%mo_coeff_b)) & !fm->dbcsr
660 0 : CPABORT("mo_coeff_b is not allocated") !fm->dbcsr
661 : CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, & !fm->dbcsr
662 7511 : mos(ispin)%mo_coeff) !fm -> dbcsr
663 : END IF !fm->dbcsr
664 : END DO !fm -> dbcsr
665 :
666 19285 : CALL cp_rm_iter_level(logger%iter_info, level_name="QS_SCF")
667 19285 : CALL timestop(handle)
668 :
669 19285 : END SUBROUTINE scf_env_do_scf
670 :
671 : ! **************************************************************************************************
672 : !> \brief inits those objects needed if you want to restart the scf with, say
673 : !> only a new initial guess, or different density functional or ...
674 : !> this will happen just before the scf loop starts
675 : !> \param scf_env ...
676 : !> \param qs_env ...
677 : !> \param scf_section ...
678 : !> \par History
679 : !> 03.2006 created [Joost VandeVondele]
680 : ! **************************************************************************************************
681 22281 : SUBROUTINE init_scf_loop(scf_env, qs_env, scf_section)
682 :
683 : TYPE(qs_scf_env_type), POINTER :: scf_env
684 : TYPE(qs_environment_type), POINTER :: qs_env
685 : TYPE(section_vals_type), POINTER :: scf_section
686 :
687 : CHARACTER(LEN=*), PARAMETER :: routineN = 'init_scf_loop'
688 :
689 : INTEGER :: handle, ispin, nmo, number_of_OT_envs
690 : LOGICAL :: do_kpoints, do_rotation, &
691 : has_unit_metric, is_full_all
692 : TYPE(cp_fm_type), POINTER :: mo_coeff
693 22281 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
694 : TYPE(dbcsr_type), POINTER :: orthogonality_metric
695 : TYPE(dft_control_type), POINTER :: dft_control
696 : TYPE(kpoint_type), POINTER :: kpoints
697 22281 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
698 : TYPE(scf_control_type), POINTER :: scf_control
699 :
700 22281 : CALL timeset(routineN, handle)
701 :
702 22281 : NULLIFY (scf_control, matrix_s, matrix_ks, dft_control, mos, mo_coeff, kpoints)
703 :
704 22281 : CPASSERT(ASSOCIATED(scf_env))
705 22281 : CPASSERT(ASSOCIATED(qs_env))
706 :
707 : CALL get_qs_env(qs_env=qs_env, &
708 : scf_control=scf_control, &
709 : dft_control=dft_control, &
710 : do_kpoints=do_kpoints, &
711 : kpoints=kpoints, &
712 22281 : mos=mos)
713 :
714 : ! if using mo_coeff_b then copy to fm
715 47525 : DO ispin = 1, SIZE(mos) !fm->dbcsr
716 47525 : IF (mos(1)%use_mo_coeff_b) THEN !fm->dbcsr
717 8542 : CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, mos(ispin)%mo_coeff) !fm->dbcsr
718 : END IF !fm->dbcsr
719 : END DO !fm->dbcsr
720 :
721 : ! this just guarantees that all mo_occupations match the eigenvalues, if smear
722 47525 : DO ispin = 1, dft_control%nspins
723 : ! do not reset mo_occupations if the maximum overlap method is in use
724 47525 : IF (.NOT. scf_control%diagonalization%mom) THEN
725 : !if the hair probes section is present, this sends hairy_probes to set_mo_occupation subroutine
726 : !and switches off the standard smearing
727 25200 : IF (dft_control%hairy_probes .EQV. .TRUE.) THEN
728 4 : IF (scf_env%outer_scf%iter_count > 0) THEN
729 0 : scf_control%smear%do_smear = .FALSE.
730 : CALL set_mo_occupation(mo_set=mos(ispin), &
731 : smear=scf_control%smear, &
732 0 : probe=dft_control%probe)
733 : END IF
734 : ELSE
735 : CALL set_mo_occupation(mo_set=mos(ispin), &
736 25196 : smear=scf_control%smear)
737 : END IF
738 : END IF
739 : END DO
740 :
741 22281 : SELECT CASE (scf_env%method)
742 : CASE DEFAULT
743 :
744 0 : CPABORT("unknown scf method method:"//cp_to_string(scf_env%method))
745 :
746 : CASE (filter_matrix_diag_method_nr)
747 :
748 10 : IF (.NOT. scf_env%skip_diis) THEN
749 0 : IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
750 0 : ALLOCATE (scf_env%scf_diis_buffer)
751 0 : CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis)
752 : END IF
753 0 : CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
754 : END IF
755 :
756 : CASE (general_diag_method_nr, special_diag_method_nr, block_krylov_diag_method_nr, smeagol_method_nr)
757 15082 : IF (.NOT. scf_env%skip_diis) THEN
758 14800 : IF (do_kpoints) THEN
759 878 : IF (.NOT. ASSOCIATED(kpoints%scf_diis_buffer)) THEN
760 140 : ALLOCATE (kpoints%scf_diis_buffer)
761 140 : CALL qs_diis_b_create_kp(kpoints%scf_diis_buffer, nbuffer=scf_control%max_diis)
762 : END IF
763 878 : CALL qs_diis_b_clear_kp(kpoints%scf_diis_buffer)
764 : ELSE
765 13922 : IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
766 4016 : ALLOCATE (scf_env%scf_diis_buffer)
767 4016 : CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis)
768 : END IF
769 13922 : CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
770 : END IF
771 : END IF
772 :
773 : CASE (ot_diag_method_nr)
774 8 : CALL get_qs_env(qs_env, matrix_ks=matrix_ks, matrix_s=matrix_s)
775 :
776 8 : IF (.NOT. scf_env%skip_diis) THEN
777 6 : IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
778 6 : ALLOCATE (scf_env%scf_diis_buffer)
779 6 : CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis)
780 : END IF
781 6 : CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
782 : END IF
783 :
784 : ! disable DFTB and SE for now
785 : IF (dft_control%qs_control%dftb .OR. &
786 8 : dft_control%qs_control%xtb .OR. &
787 : dft_control%qs_control%semi_empirical) THEN
788 0 : CPABORT("DFTB and SE not available with OT/DIAG")
789 : END IF
790 :
791 : ! if an old preconditioner is still around (i.e. outer SCF is active),
792 : ! remove it if this could be worthwhile
793 : CALL restart_preconditioner(qs_env, scf_env%ot_preconditioner, &
794 : scf_control%diagonalization%ot_settings%preconditioner_type, &
795 8 : dft_control%nspins)
796 :
797 : CALL prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, scf_env%ot_preconditioner, &
798 : scf_control%diagonalization%ot_settings%preconditioner_type, &
799 : scf_control%diagonalization%ot_settings%precond_solver_type, &
800 8 : scf_control%diagonalization%ot_settings%energy_gap, dft_control%nspins)
801 :
802 : CASE (block_davidson_diag_method_nr)
803 : ! Preconditioner initialized within the loop, when required
804 : CASE (ot_method_nr)
805 : CALL get_qs_env(qs_env, &
806 : has_unit_metric=has_unit_metric, &
807 : matrix_s=matrix_s, &
808 7165 : matrix_ks=matrix_ks)
809 :
810 : ! reortho the wavefunctions if we are having an outer scf and
811 : ! this is not the first iteration
812 : ! this is useful to avoid the build-up of numerical noise
813 : ! however, we can not play this trick if restricted (don't mix non-equivalent orbs)
814 7165 : IF (scf_control%do_outer_scf_reortho) THEN
815 6595 : IF (scf_control%outer_scf%have_scf .AND. .NOT. dft_control%restricted) THEN
816 4541 : IF (scf_env%outer_scf%iter_count > 0) THEN
817 1829 : DO ispin = 1, dft_control%nspins
818 993 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
819 1829 : IF (has_unit_metric) THEN
820 108 : CALL make_basis_simple(mo_coeff, nmo)
821 : ELSE
822 885 : CALL make_basis_sm(mo_coeff, nmo, matrix_s(1)%matrix)
823 : END IF
824 : END DO
825 : END IF
826 : END IF
827 : ELSE
828 : ! dont need any dirty trick for the numerically stable irac algorithm.
829 : END IF
830 :
831 7165 : IF (.NOT. ASSOCIATED(scf_env%qs_ot_env)) THEN
832 :
833 : ! restricted calculations require just one set of OT orbitals
834 7165 : number_of_OT_envs = dft_control%nspins
835 7165 : IF (dft_control%restricted) number_of_OT_envs = 1
836 :
837 1183480 : ALLOCATE (scf_env%qs_ot_env(number_of_OT_envs))
838 :
839 : ! XXX Joost XXX should disentangle reading input from this part
840 7165 : IF (scf_env%outer_scf%iter_count > 0) THEN
841 856 : IF (scf_env%iter_delta < scf_control%eps_diis) THEN
842 4 : scf_env%qs_ot_env(1)%settings%ot_state = 1
843 : END IF
844 : END IF
845 : !
846 7165 : CALL ot_scf_read_input(scf_env%qs_ot_env, scf_section)
847 : !
848 7165 : IF (scf_env%outer_scf%iter_count > 0) THEN
849 856 : IF (scf_env%qs_ot_env(1)%settings%ot_state == 1) THEN
850 : scf_control%max_scf = MAX(scf_env%qs_ot_env(1)%settings%max_scf_diis, &
851 4 : scf_control%max_scf)
852 : END IF
853 : END IF
854 :
855 : ! keep a note that we are restricted
856 7165 : IF (dft_control%restricted) THEN
857 92 : scf_env%qs_ot_env(1)%restricted = .TRUE.
858 : ! requires rotation
859 92 : IF (.NOT. scf_env%qs_ot_env(1)%settings%do_rotation) &
860 : CALL cp_abort(__LOCATION__, &
861 : "Restricted calculation with OT requires orbital rotation. Please "// &
862 0 : "activate the OT%ROTATION keyword!")
863 : ELSE
864 15401 : scf_env%qs_ot_env(:)%restricted = .FALSE.
865 : END IF
866 :
867 : ! this will rotate the MOs to be eigen states, which is not compatible with rotation
868 : ! e.g. mo_derivs here do not yet include potentially different occupations numbers
869 7165 : do_rotation = scf_env%qs_ot_env(1)%settings%do_rotation
870 : ! only full all needs rotation
871 7165 : is_full_all = scf_env%qs_ot_env(1)%settings%preconditioner_type == ot_precond_full_all
872 7165 : IF (do_rotation .AND. is_full_all) &
873 0 : CPABORT('PRECONDITIONER FULL_ALL is not compatible with ROTATION.')
874 :
875 : ! might need the KS matrix to init properly
876 : CALL qs_ks_update_qs_env(qs_env, just_energy=.FALSE., &
877 7165 : calculate_forces=.FALSE.)
878 :
879 : ! if an old preconditioner is still around (i.e. outer SCF is active),
880 : ! remove it if this could be worthwhile
881 7165 : IF (.NOT. reuse_precond) &
882 : CALL restart_preconditioner(qs_env, scf_env%ot_preconditioner, &
883 : scf_env%qs_ot_env(1)%settings%preconditioner_type, &
884 7165 : dft_control%nspins)
885 :
886 : !
887 : ! preconditioning still needs to be done correctly with has_unit_metric
888 : ! notice that a big part of the preconditioning (S^-1) is fine anyhow
889 : !
890 7165 : IF (has_unit_metric) THEN
891 1154 : NULLIFY (orthogonality_metric)
892 : ELSE
893 6011 : orthogonality_metric => matrix_s(1)%matrix
894 : END IF
895 :
896 7165 : IF (.NOT. reuse_precond) &
897 : CALL prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, scf_env%ot_preconditioner, &
898 : scf_env%qs_ot_env(1)%settings%preconditioner_type, &
899 : scf_env%qs_ot_env(1)%settings%precond_solver_type, &
900 : scf_env%qs_ot_env(1)%settings%energy_gap, dft_control%nspins, &
901 : has_unit_metric=has_unit_metric, &
902 7165 : chol_type=scf_env%qs_ot_env(1)%settings%cholesky_type)
903 7165 : IF (reuse_precond) reuse_precond = .FALSE.
904 :
905 : CALL ot_scf_init(mo_array=mos, matrix_s=orthogonality_metric, &
906 : broyden_adaptive_sigma=qs_env%broyden_adaptive_sigma, &
907 7165 : qs_ot_env=scf_env%qs_ot_env, matrix_ks=matrix_ks(1)%matrix)
908 :
909 12325 : SELECT CASE (scf_env%qs_ot_env(1)%settings%preconditioner_type)
910 : CASE (ot_precond_none)
911 : CASE (ot_precond_full_all, ot_precond_full_single_inverse)
912 11289 : DO ispin = 1, SIZE(scf_env%qs_ot_env)
913 : CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin), &
914 11289 : scf_env%ot_preconditioner(ispin)%preconditioner)
915 : END DO
916 : CASE (ot_precond_s_inverse, ot_precond_full_single)
917 152 : DO ispin = 1, SIZE(scf_env%qs_ot_env)
918 : CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin), &
919 152 : scf_env%ot_preconditioner(1)%preconditioner)
920 : END DO
921 : CASE DEFAULT
922 8616 : DO ispin = 1, SIZE(scf_env%qs_ot_env)
923 : CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin), &
924 2666 : scf_env%ot_preconditioner(1)%preconditioner)
925 : END DO
926 : END SELECT
927 : END IF
928 :
929 : ! if we have non-uniform occupations we should be using rotation
930 7165 : do_rotation = scf_env%qs_ot_env(1)%settings%do_rotation
931 37958 : DO ispin = 1, SIZE(mos)
932 15677 : IF (.NOT. mos(ispin)%uniform_occupation) THEN
933 0 : CPASSERT(do_rotation)
934 : END IF
935 : END DO
936 : END SELECT
937 :
938 : ! another safety check
939 22281 : IF (dft_control%low_spin_roks) THEN
940 24 : CPASSERT(scf_env%method == ot_method_nr)
941 24 : do_rotation = scf_env%qs_ot_env(1)%settings%do_rotation
942 24 : CPASSERT(do_rotation)
943 : END IF
944 :
945 22281 : CALL timestop(handle)
946 :
947 22281 : END SUBROUTINE init_scf_loop
948 :
949 : ! **************************************************************************************************
950 : !> \brief perform cleanup operations (like releasing temporary storage)
951 : !> at the end of the scf
952 : !> \param scf_env ...
953 : !> \par History
954 : !> 02.2003 created [fawzi]
955 : !> \author fawzi
956 : ! **************************************************************************************************
957 19029 : SUBROUTINE scf_env_cleanup(scf_env)
958 : TYPE(qs_scf_env_type), INTENT(INOUT) :: scf_env
959 :
960 : CHARACTER(len=*), PARAMETER :: routineN = 'scf_env_cleanup'
961 :
962 : INTEGER :: handle
963 :
964 19029 : CALL timeset(routineN, handle)
965 :
966 : ! Release SCF work storage
967 19029 : CALL cp_fm_release(scf_env%scf_work1)
968 :
969 19029 : IF (ASSOCIATED(scf_env%scf_work1_red)) THEN
970 48 : CALL cp_fm_release(scf_env%scf_work1_red)
971 : END IF
972 19029 : IF (ASSOCIATED(scf_env%scf_work2)) THEN
973 12896 : CALL cp_fm_release(scf_env%scf_work2)
974 12896 : DEALLOCATE (scf_env%scf_work2)
975 : NULLIFY (scf_env%scf_work2)
976 : END IF
977 19029 : IF (ASSOCIATED(scf_env%scf_work2_red)) THEN
978 48 : CALL cp_fm_release(scf_env%scf_work2_red)
979 48 : DEALLOCATE (scf_env%scf_work2_red)
980 : NULLIFY (scf_env%scf_work2_red)
981 : END IF
982 19029 : IF (ASSOCIATED(scf_env%ortho)) THEN
983 10232 : CALL cp_fm_release(scf_env%ortho)
984 10232 : DEALLOCATE (scf_env%ortho)
985 : NULLIFY (scf_env%ortho)
986 : END IF
987 19029 : IF (ASSOCIATED(scf_env%ortho_red)) THEN
988 48 : CALL cp_fm_release(scf_env%ortho_red)
989 48 : DEALLOCATE (scf_env%ortho_red)
990 : NULLIFY (scf_env%ortho_red)
991 : END IF
992 19029 : IF (ASSOCIATED(scf_env%ortho_m1)) THEN
993 48 : CALL cp_fm_release(scf_env%ortho_m1)
994 48 : DEALLOCATE (scf_env%ortho_m1)
995 : NULLIFY (scf_env%ortho_m1)
996 : END IF
997 19029 : IF (ASSOCIATED(scf_env%ortho_m1_red)) THEN
998 6 : CALL cp_fm_release(scf_env%ortho_m1_red)
999 6 : DEALLOCATE (scf_env%ortho_m1_red)
1000 : NULLIFY (scf_env%ortho_m1_red)
1001 : END IF
1002 :
1003 19029 : IF (ASSOCIATED(scf_env%ortho_dbcsr)) THEN
1004 58 : CALL dbcsr_deallocate_matrix(scf_env%ortho_dbcsr)
1005 : END IF
1006 19029 : IF (ASSOCIATED(scf_env%buf1_dbcsr)) THEN
1007 58 : CALL dbcsr_deallocate_matrix(scf_env%buf1_dbcsr)
1008 : END IF
1009 19029 : IF (ASSOCIATED(scf_env%buf2_dbcsr)) THEN
1010 58 : CALL dbcsr_deallocate_matrix(scf_env%buf2_dbcsr)
1011 : END IF
1012 :
1013 19029 : IF (ASSOCIATED(scf_env%p_mix_new)) THEN
1014 12912 : CALL dbcsr_deallocate_matrix_set(scf_env%p_mix_new)
1015 : END IF
1016 :
1017 19029 : IF (ASSOCIATED(scf_env%p_delta)) THEN
1018 276 : CALL dbcsr_deallocate_matrix_set(scf_env%p_delta)
1019 : END IF
1020 :
1021 : ! Method dependent cleanup
1022 19045 : SELECT CASE (scf_env%method)
1023 : CASE (ot_method_nr)
1024 : !
1025 : CASE (ot_diag_method_nr)
1026 : !
1027 : CASE (general_diag_method_nr)
1028 : !
1029 : CASE (special_diag_method_nr)
1030 : !
1031 : CASE (block_krylov_diag_method_nr)
1032 : CASE (block_davidson_diag_method_nr)
1033 16 : CALL block_davidson_deallocate(scf_env%block_davidson_env)
1034 : CASE (filter_matrix_diag_method_nr)
1035 : !
1036 : CASE (smeagol_method_nr)
1037 : !
1038 : CASE DEFAULT
1039 19029 : CPABORT("unknown scf method method:"//cp_to_string(scf_env%method))
1040 : END SELECT
1041 :
1042 19029 : IF (ASSOCIATED(scf_env%outer_scf%variables)) THEN
1043 4139 : DEALLOCATE (scf_env%outer_scf%variables)
1044 : END IF
1045 19029 : IF (ASSOCIATED(scf_env%outer_scf%count)) THEN
1046 4139 : DEALLOCATE (scf_env%outer_scf%count)
1047 : END IF
1048 19029 : IF (ASSOCIATED(scf_env%outer_scf%gradient)) THEN
1049 4139 : DEALLOCATE (scf_env%outer_scf%gradient)
1050 : END IF
1051 19029 : IF (ASSOCIATED(scf_env%outer_scf%energy)) THEN
1052 4139 : DEALLOCATE (scf_env%outer_scf%energy)
1053 : END IF
1054 19029 : IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian) .AND. &
1055 : scf_env%outer_scf%deallocate_jacobian) THEN
1056 50 : DEALLOCATE (scf_env%outer_scf%inv_jacobian)
1057 : END IF
1058 :
1059 19029 : CALL timestop(handle)
1060 :
1061 19029 : END SUBROUTINE scf_env_cleanup
1062 :
1063 : ! **************************************************************************************************
1064 : !> \brief perform a CDFT scf procedure in the given qs_env
1065 : !> \param qs_env the qs_environment where to perform the scf procedure
1066 : !> \param should_stop flag determining if calculation should stop
1067 : !> \par History
1068 : !> 12.2015 Created
1069 : !> \author Nico Holmberg
1070 : ! **************************************************************************************************
1071 652 : SUBROUTINE cdft_scf(qs_env, should_stop)
1072 : TYPE(qs_environment_type), POINTER :: qs_env
1073 : LOGICAL, INTENT(OUT) :: should_stop
1074 :
1075 : CHARACTER(len=*), PARAMETER :: routineN = 'cdft_scf'
1076 :
1077 : INTEGER :: handle, iatom, ispin, ivar, nmo, nvar, &
1078 : output_unit, tsteps
1079 : LOGICAL :: cdft_loop_converged, converged, &
1080 : exit_cdft_loop, first_iteration, &
1081 : my_uocc, uniform_occupation
1082 326 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_occupations
1083 : TYPE(cdft_control_type), POINTER :: cdft_control
1084 : TYPE(cp_logger_type), POINTER :: logger
1085 326 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
1086 : TYPE(dft_control_type), POINTER :: dft_control
1087 326 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1088 : TYPE(pw_env_type), POINTER :: pw_env
1089 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1090 : TYPE(qs_energy_type), POINTER :: energy
1091 : TYPE(qs_ks_env_type), POINTER :: ks_env
1092 : TYPE(qs_rho_type), POINTER :: rho
1093 : TYPE(qs_scf_env_type), POINTER :: scf_env
1094 : TYPE(scf_control_type), POINTER :: scf_control
1095 : TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
1096 :
1097 326 : NULLIFY (scf_env, ks_env, energy, rho, matrix_s, rho_ao, cdft_control, logger, &
1098 326 : dft_control, pw_env, auxbas_pw_pool, energy, ks_env, scf_env, dft_section, &
1099 326 : input, scf_section, scf_control, mos, mo_occupations)
1100 652 : logger => cp_get_default_logger()
1101 :
1102 326 : CPASSERT(ASSOCIATED(qs_env))
1103 : CALL get_qs_env(qs_env, scf_env=scf_env, energy=energy, &
1104 : dft_control=dft_control, scf_control=scf_control, &
1105 326 : ks_env=ks_env, input=input)
1106 :
1107 326 : CALL timeset(routineN//"_loop", handle)
1108 326 : dft_section => section_vals_get_subs_vals(input, "DFT")
1109 326 : scf_section => section_vals_get_subs_vals(dft_section, "SCF")
1110 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", &
1111 326 : extension=".scfLog")
1112 326 : first_iteration = .TRUE.
1113 :
1114 326 : cdft_control => dft_control%qs_control%cdft_control
1115 :
1116 326 : scf_env%outer_scf%iter_count = 0
1117 326 : cdft_control%total_steps = 0
1118 :
1119 : ! Write some info about the CDFT calculation
1120 326 : IF (output_unit > 0) THEN
1121 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
1122 181 : "CDFT EXTERNAL SCF WAVEFUNCTION OPTIMIZATION"
1123 181 : CALL qs_scf_cdft_initial_info(output_unit, cdft_control)
1124 : END IF
1125 326 : IF (cdft_control%reuse_precond) THEN
1126 0 : reuse_precond = .FALSE.
1127 0 : cdft_control%nreused = 0
1128 : END IF
1129 512 : cdft_outer_loop: DO
1130 : ! Change outer_scf settings to OT settings
1131 512 : CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot)
1132 : ! Solve electronic structure with fixed value of constraint
1133 : CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
1134 512 : converged=converged, should_stop=should_stop, total_scf_steps=tsteps)
1135 : ! Decide whether to reuse the preconditioner on the next iteration
1136 512 : IF (cdft_control%reuse_precond) THEN
1137 : ! For convergence in exactly one step, the preconditioner is always reused (assuming max_reuse > 0)
1138 : ! usually this means that the electronic structure has already converged to the correct state
1139 : ! but the constraint optimizer keeps jumping over the optimal solution
1140 : IF (scf_env%outer_scf%iter_count == 1 .AND. scf_env%iter_count == 1 &
1141 0 : .AND. cdft_control%total_steps /= 1) &
1142 0 : cdft_control%nreused = cdft_control%nreused - 1
1143 : ! SCF converged in less than precond_freq steps
1144 : IF (scf_env%outer_scf%iter_count == 1 .AND. scf_env%iter_count <= cdft_control%precond_freq .AND. &
1145 0 : cdft_control%total_steps /= 1 .AND. cdft_control%nreused < cdft_control%max_reuse) THEN
1146 0 : reuse_precond = .TRUE.
1147 0 : cdft_control%nreused = cdft_control%nreused + 1
1148 : ELSE
1149 0 : reuse_precond = .FALSE.
1150 0 : cdft_control%nreused = 0
1151 : END IF
1152 : END IF
1153 : ! Update history purging counters
1154 512 : IF (first_iteration .AND. cdft_control%purge_history) THEN
1155 0 : cdft_control%istep = cdft_control%istep + 1
1156 0 : IF (scf_env%outer_scf%iter_count > 1) THEN
1157 0 : cdft_control%nbad_conv = cdft_control%nbad_conv + 1
1158 0 : IF (cdft_control%nbad_conv >= cdft_control%purge_freq .AND. &
1159 : cdft_control%istep >= cdft_control%purge_offset) THEN
1160 0 : cdft_control%nbad_conv = 0
1161 0 : cdft_control%istep = 0
1162 0 : cdft_control%should_purge = .TRUE.
1163 : END IF
1164 : END IF
1165 : END IF
1166 512 : first_iteration = .FALSE.
1167 : ! Change outer_scf settings to CDFT settings
1168 512 : CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft)
1169 : CALL qs_scf_check_outer_exit(qs_env, scf_env, scf_control, should_stop, &
1170 512 : cdft_loop_converged, exit_cdft_loop)
1171 : CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
1172 : energy, cdft_control%total_steps, &
1173 512 : should_stop, cdft_loop_converged, cdft_loop=.TRUE.)
1174 512 : IF (exit_cdft_loop) EXIT cdft_outer_loop
1175 : ! Check if the inverse Jacobian needs to be calculated
1176 186 : CALL qs_calculate_inverse_jacobian(qs_env)
1177 : ! Check if a line search should be performed to find an optimal step size for the optimizer
1178 186 : CALL qs_cdft_line_search(qs_env)
1179 : ! Optimize constraint
1180 186 : CALL outer_loop_optimize(scf_env, scf_control)
1181 186 : CALL outer_loop_update_qs_env(qs_env, scf_env)
1182 512 : CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.)
1183 : END DO cdft_outer_loop
1184 :
1185 326 : cdft_control%ienergy = cdft_control%ienergy + 1
1186 :
1187 : ! Store needed arrays for ET coupling calculation
1188 326 : IF (cdft_control%do_et) THEN
1189 176 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, mos=mos)
1190 176 : nvar = SIZE(cdft_control%target)
1191 : ! Matrix representation of weight function
1192 708 : ALLOCATE (cdft_control%wmat(nvar))
1193 356 : DO ivar = 1, nvar
1194 180 : CALL dbcsr_init_p(cdft_control%wmat(ivar)%matrix)
1195 : CALL dbcsr_copy(cdft_control%wmat(ivar)%matrix, matrix_s(1)%matrix, &
1196 180 : name="ET_RESTRAINT_MATRIX")
1197 180 : CALL dbcsr_set(cdft_control%wmat(ivar)%matrix, 0.0_dp)
1198 : CALL integrate_v_rspace(cdft_control%group(ivar)%weight, &
1199 : hmat=cdft_control%wmat(ivar), qs_env=qs_env, &
1200 : calculate_forces=.FALSE., &
1201 356 : gapw=dft_control%qs_control%gapw)
1202 : END DO
1203 : ! Overlap matrix
1204 176 : CALL dbcsr_init_p(cdft_control%matrix_s%matrix)
1205 : CALL dbcsr_copy(cdft_control%matrix_s%matrix, matrix_s(1)%matrix, &
1206 176 : name="OVERLAP")
1207 : ! Molecular orbital coefficients
1208 176 : NULLIFY (cdft_control%mo_coeff)
1209 880 : ALLOCATE (cdft_control%mo_coeff(dft_control%nspins))
1210 528 : DO ispin = 1, dft_control%nspins
1211 : CALL cp_fm_create(matrix=cdft_control%mo_coeff(ispin), &
1212 : matrix_struct=qs_env%mos(ispin)%mo_coeff%matrix_struct, &
1213 352 : name="MO_COEFF_A"//TRIM(ADJUSTL(cp_to_string(ispin)))//"MATRIX")
1214 : CALL cp_fm_to_fm(qs_env%mos(ispin)%mo_coeff, &
1215 528 : cdft_control%mo_coeff(ispin))
1216 : END DO
1217 : ! Density matrix
1218 176 : IF (cdft_control%calculate_metric) THEN
1219 24 : CALL get_qs_env(qs_env, rho=rho)
1220 24 : CALL qs_rho_get(rho, rho_ao=rho_ao)
1221 120 : ALLOCATE (cdft_control%matrix_p(dft_control%nspins))
1222 72 : DO ispin = 1, dft_control%nspins
1223 48 : NULLIFY (cdft_control%matrix_p(ispin)%matrix)
1224 48 : CALL dbcsr_init_p(cdft_control%matrix_p(ispin)%matrix)
1225 : CALL dbcsr_copy(cdft_control%matrix_p(ispin)%matrix, rho_ao(ispin)%matrix, &
1226 72 : name="DENSITY MATRIX")
1227 : END DO
1228 : END IF
1229 : ! Copy occupation numbers if non-uniform occupation
1230 176 : uniform_occupation = .TRUE.
1231 528 : DO ispin = 1, dft_control%nspins
1232 352 : CALL get_mo_set(mo_set=mos(ispin), uniform_occupation=my_uocc)
1233 584 : uniform_occupation = uniform_occupation .AND. my_uocc
1234 : END DO
1235 176 : IF (.NOT. uniform_occupation) THEN
1236 140 : ALLOCATE (cdft_control%occupations(dft_control%nspins))
1237 84 : DO ispin = 1, dft_control%nspins
1238 : CALL get_mo_set(mo_set=mos(ispin), &
1239 : nmo=nmo, &
1240 56 : occupation_numbers=mo_occupations)
1241 168 : ALLOCATE (cdft_control%occupations(ispin)%array(nmo))
1242 588 : cdft_control%occupations(ispin)%array(1:nmo) = mo_occupations(1:nmo)
1243 : END DO
1244 : END IF
1245 : END IF
1246 :
1247 : ! Deallocate constraint storage if forces are not needed
1248 : ! In case of a simulation with multiple force_evals,
1249 : ! deallocate only if weight function should not be copied to different force_evals
1250 326 : IF (.NOT. (cdft_control%save_pot .OR. cdft_control%transfer_pot)) THEN
1251 148 : CALL get_qs_env(qs_env, pw_env=pw_env)
1252 148 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1253 308 : DO iatom = 1, SIZE(cdft_control%group)
1254 160 : CALL auxbas_pw_pool%give_back_pw(cdft_control%group(iatom)%weight)
1255 308 : DEALLOCATE (cdft_control%group(iatom)%weight)
1256 : END DO
1257 148 : IF (cdft_control%atomic_charges) THEN
1258 256 : DO iatom = 1, cdft_control%natoms
1259 256 : CALL auxbas_pw_pool%give_back_pw(cdft_control%charge(iatom))
1260 : END DO
1261 84 : DEALLOCATE (cdft_control%charge)
1262 : END IF
1263 148 : IF (cdft_control%type == outer_scf_becke_constraint .AND. &
1264 : cdft_control%becke_control%cavity_confine) THEN
1265 120 : IF (.NOT. ASSOCIATED(cdft_control%becke_control%cavity_mat)) THEN
1266 110 : CALL auxbas_pw_pool%give_back_pw(cdft_control%becke_control%cavity)
1267 : ELSE
1268 10 : DEALLOCATE (cdft_control%becke_control%cavity_mat)
1269 : END IF
1270 28 : ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1271 20 : IF (ASSOCIATED(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)) THEN
1272 0 : CALL auxbas_pw_pool%give_back_pw(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)
1273 : END IF
1274 : END IF
1275 148 : IF (ASSOCIATED(cdft_control%charges_fragment)) DEALLOCATE (cdft_control%charges_fragment)
1276 148 : cdft_control%need_pot = .TRUE.
1277 148 : cdft_control%external_control = .FALSE.
1278 : END IF
1279 :
1280 326 : CALL timestop(handle)
1281 :
1282 326 : END SUBROUTINE cdft_scf
1283 :
1284 : ! **************************************************************************************************
1285 : !> \brief perform cleanup operations for cdft_control
1286 : !> \param cdft_control container for the external CDFT SCF loop variables
1287 : !> \par History
1288 : !> 12.2015 created [Nico Holmberg]
1289 : !> \author Nico Holmberg
1290 : ! **************************************************************************************************
1291 326 : SUBROUTINE cdft_control_cleanup(cdft_control)
1292 : TYPE(cdft_control_type), POINTER :: cdft_control
1293 :
1294 326 : IF (ASSOCIATED(cdft_control%constraint%variables)) &
1295 326 : DEALLOCATE (cdft_control%constraint%variables)
1296 326 : IF (ASSOCIATED(cdft_control%constraint%count)) &
1297 326 : DEALLOCATE (cdft_control%constraint%count)
1298 326 : IF (ASSOCIATED(cdft_control%constraint%gradient)) &
1299 326 : DEALLOCATE (cdft_control%constraint%gradient)
1300 326 : IF (ASSOCIATED(cdft_control%constraint%energy)) &
1301 326 : DEALLOCATE (cdft_control%constraint%energy)
1302 326 : IF (ASSOCIATED(cdft_control%constraint%inv_jacobian) .AND. &
1303 : cdft_control%constraint%deallocate_jacobian) &
1304 4 : DEALLOCATE (cdft_control%constraint%inv_jacobian)
1305 :
1306 326 : END SUBROUTINE cdft_control_cleanup
1307 :
1308 : ! **************************************************************************************************
1309 : !> \brief Calculates the finite difference inverse Jacobian
1310 : !> \param qs_env the qs_environment_type where to compute the Jacobian
1311 : !> \par History
1312 : !> 01.2017 created [Nico Holmberg]
1313 : ! **************************************************************************************************
1314 186 : SUBROUTINE qs_calculate_inverse_jacobian(qs_env)
1315 : TYPE(qs_environment_type), POINTER :: qs_env
1316 :
1317 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_calculate_inverse_jacobian'
1318 :
1319 : CHARACTER(len=default_path_length) :: project_name
1320 : INTEGER :: counter, handle, i, ispin, iter_count, &
1321 : iwork, j, max_scf, nspins, nsteps, &
1322 : nvar, nwork, output_unit, pwork, &
1323 : tsteps, twork
1324 : LOGICAL :: converged, explicit_jacobian, &
1325 : should_build, should_stop, &
1326 : use_md_history
1327 : REAL(KIND=dp) :: inv_error, step_size
1328 186 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coeff, dh, step_multiplier
1329 186 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: jacobian
1330 186 : REAL(KIND=dp), DIMENSION(:), POINTER :: energy
1331 186 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient, inv_jacobian
1332 : TYPE(cdft_control_type), POINTER :: cdft_control
1333 : TYPE(cp_logger_type), POINTER :: logger, tmp_logger
1334 186 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv
1335 186 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
1336 : TYPE(dft_control_type), POINTER :: dft_control
1337 186 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mos_stashed
1338 186 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1339 : TYPE(mp_para_env_type), POINTER :: para_env
1340 : TYPE(qs_energy_type), POINTER :: energy_qs
1341 : TYPE(qs_ks_env_type), POINTER :: ks_env
1342 : TYPE(qs_rho_type), POINTER :: rho
1343 : TYPE(qs_scf_env_type), POINTER :: scf_env
1344 : TYPE(scf_control_type), POINTER :: scf_control
1345 :
1346 186 : NULLIFY (energy, gradient, p_rmpv, rho_ao_kp, mos, rho, &
1347 186 : ks_env, scf_env, scf_control, dft_control, cdft_control, &
1348 186 : inv_jacobian, para_env, tmp_logger, energy_qs)
1349 372 : logger => cp_get_default_logger()
1350 :
1351 186 : CPASSERT(ASSOCIATED(qs_env))
1352 : CALL get_qs_env(qs_env, scf_env=scf_env, ks_env=ks_env, &
1353 : scf_control=scf_control, mos=mos, rho=rho, &
1354 : dft_control=dft_control, &
1355 186 : para_env=para_env, energy=energy_qs)
1356 186 : explicit_jacobian = .FALSE.
1357 186 : should_build = .FALSE.
1358 186 : use_md_history = .FALSE.
1359 186 : iter_count = scf_env%outer_scf%iter_count
1360 : ! Quick exit if optimizer does not require Jacobian
1361 186 : IF (.NOT. ASSOCIATED(scf_control%outer_scf%cdft_opt_control)) RETURN
1362 : ! Check if Jacobian should be calculated and initialize
1363 118 : CALL timeset(routineN, handle)
1364 118 : CALL initialize_inverse_jacobian(scf_control, scf_env, explicit_jacobian, should_build, used_history)
1365 118 : IF (scf_control%outer_scf%cdft_opt_control%jacobian_restart) THEN
1366 : ! Restart from previously calculated inverse Jacobian
1367 6 : should_build = .FALSE.
1368 6 : CALL restart_inverse_jacobian(qs_env)
1369 : END IF
1370 118 : IF (should_build) THEN
1371 78 : scf_env%outer_scf%deallocate_jacobian = .FALSE.
1372 : ! Actually need to (re)build the Jacobian
1373 78 : IF (explicit_jacobian) THEN
1374 : ! Build Jacobian with finite differences
1375 62 : cdft_control => dft_control%qs_control%cdft_control
1376 62 : IF (.NOT. ASSOCIATED(cdft_control)) &
1377 : CALL cp_abort(__LOCATION__, &
1378 : "Optimizers that need the explicit Jacobian can"// &
1379 0 : " only be used together with a valid CDFT constraint.")
1380 : ! Redirect output from Jacobian calculation to a new file by creating a temporary logger
1381 62 : project_name = logger%iter_info%project_name
1382 62 : CALL create_tmp_logger(para_env, project_name, "-JacobianInfo.out", output_unit, tmp_logger)
1383 : ! Save last converged state so we can roll back to it (mo_coeff and some outer_loop variables)
1384 62 : nspins = dft_control%nspins
1385 310 : ALLOCATE (mos_stashed(nspins))
1386 186 : DO ispin = 1, nspins
1387 186 : CALL duplicate_mo_set(mos_stashed(ispin), mos(ispin))
1388 : END DO
1389 62 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1390 62 : p_rmpv => rho_ao_kp(:, 1)
1391 : ! Allocate work
1392 62 : nvar = SIZE(scf_env%outer_scf%variables, 1)
1393 62 : max_scf = scf_control%outer_scf%max_scf + 1
1394 248 : ALLOCATE (gradient(nvar, max_scf))
1395 1310 : gradient = scf_env%outer_scf%gradient
1396 186 : ALLOCATE (energy(max_scf))
1397 594 : energy = scf_env%outer_scf%energy
1398 248 : ALLOCATE (jacobian(nvar, nvar))
1399 282 : jacobian = 0.0_dp
1400 62 : nsteps = cdft_control%total_steps
1401 : ! Setup finite difference scheme
1402 62 : CALL prepare_jacobian_stencil(qs_env, output_unit, nwork, pwork, coeff, step_multiplier, dh)
1403 62 : twork = pwork - nwork
1404 148 : DO i = 1, nvar
1405 282 : jacobian(i, :) = coeff(0)*scf_env%outer_scf%gradient(i, iter_count)
1406 : END DO
1407 : ! Calculate the Jacobian by perturbing each Lagrangian and recalculating the energy self-consistently
1408 62 : CALL cp_add_default_logger(tmp_logger)
1409 148 : DO i = 1, nvar
1410 86 : IF (output_unit > 0) THEN
1411 43 : WRITE (output_unit, FMT="(A)") " "
1412 43 : WRITE (output_unit, FMT="(A)") " #####################################"
1413 : WRITE (output_unit, '(A,I3,A,I3,A)') &
1414 43 : " ### Constraint ", i, " of ", nvar, " ###"
1415 43 : WRITE (output_unit, FMT="(A)") " #####################################"
1416 : END IF
1417 86 : counter = 0
1418 332 : DO iwork = nwork, pwork
1419 184 : IF (iwork == 0) CYCLE
1420 98 : counter = counter + 1
1421 98 : IF (output_unit > 0) THEN
1422 49 : WRITE (output_unit, FMT="(A)") " #####################################"
1423 : WRITE (output_unit, '(A,I3,A,I3,A)') &
1424 49 : " ### Energy evaluation ", counter, " of ", twork, " ###"
1425 49 : WRITE (output_unit, FMT="(A)") " #####################################"
1426 : END IF
1427 98 : IF (SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step) == 1) THEN
1428 90 : step_size = scf_control%outer_scf%cdft_opt_control%jacobian_step(1)
1429 : ELSE
1430 8 : step_size = scf_control%outer_scf%cdft_opt_control%jacobian_step(i)
1431 : END IF
1432 244 : scf_env%outer_scf%variables(:, iter_count + 1) = scf_env%outer_scf%variables(:, iter_count)
1433 : scf_env%outer_scf%variables(i, iter_count + 1) = scf_env%outer_scf%variables(i, iter_count) + &
1434 98 : step_multiplier(iwork)*step_size
1435 98 : CALL outer_loop_update_qs_env(qs_env, scf_env)
1436 98 : CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.)
1437 98 : CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot)
1438 : CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
1439 98 : converged=converged, should_stop=should_stop, total_scf_steps=tsteps)
1440 98 : CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft)
1441 : ! Update (iter_count + 1) element of gradient and print constraint info
1442 98 : scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count + 1
1443 98 : CALL outer_loop_gradient(qs_env, scf_env)
1444 : CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
1445 : energy_qs, cdft_control%total_steps, &
1446 98 : should_stop=.FALSE., outer_loop_converged=.FALSE., cdft_loop=.FALSE.)
1447 98 : scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count - 1
1448 : ! Update Jacobian
1449 244 : DO j = 1, nvar
1450 244 : jacobian(j, i) = jacobian(j, i) + coeff(iwork)*scf_env%outer_scf%gradient(j, iter_count + 1)
1451 : END DO
1452 : ! Reset everything to last converged state
1453 244 : scf_env%outer_scf%variables(:, iter_count + 1) = 0.0_dp
1454 2026 : scf_env%outer_scf%gradient = gradient
1455 878 : scf_env%outer_scf%energy = energy
1456 98 : cdft_control%total_steps = nsteps
1457 294 : DO ispin = 1, nspins
1458 196 : CALL deallocate_mo_set(mos(ispin))
1459 196 : CALL duplicate_mo_set(mos(ispin), mos_stashed(ispin))
1460 : CALL calculate_density_matrix(mos(ispin), &
1461 294 : p_rmpv(ispin)%matrix)
1462 : END DO
1463 98 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1464 368 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1465 : END DO
1466 : END DO
1467 62 : CALL cp_rm_default_logger()
1468 62 : CALL cp_logger_release(tmp_logger)
1469 : ! Finalize and invert Jacobian
1470 148 : DO j = 1, nvar
1471 282 : DO i = 1, nvar
1472 220 : jacobian(i, j) = jacobian(i, j)/dh(j)
1473 : END DO
1474 : END DO
1475 62 : IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
1476 102 : ALLOCATE (scf_env%outer_scf%inv_jacobian(nvar, nvar))
1477 62 : inv_jacobian => scf_env%outer_scf%inv_jacobian
1478 62 : CALL invert_matrix(jacobian, inv_jacobian, inv_error)
1479 62 : scf_control%outer_scf%cdft_opt_control%broyden_update = .FALSE.
1480 : ! Release temporary storage
1481 186 : DO ispin = 1, nspins
1482 186 : CALL deallocate_mo_set(mos_stashed(ispin))
1483 : END DO
1484 62 : DEALLOCATE (mos_stashed, jacobian, gradient, energy, coeff, step_multiplier, dh)
1485 186 : IF (output_unit > 0) THEN
1486 : WRITE (output_unit, FMT="(/,A)") &
1487 31 : " ================================== JACOBIAN CALCULATED =================================="
1488 31 : CALL close_file(unit_number=output_unit)
1489 : END IF
1490 : ELSE
1491 : ! Build a strictly diagonal Jacobian from history and invert it
1492 16 : CALL build_diagonal_jacobian(qs_env, used_history)
1493 : END IF
1494 : END IF
1495 118 : IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian) .AND. para_env%is_source()) THEN
1496 : ! Write restart file for inverse Jacobian
1497 55 : CALL print_inverse_jacobian(logger, scf_env%outer_scf%inv_jacobian, iter_count)
1498 : END IF
1499 : ! Update counter
1500 118 : scf_control%outer_scf%cdft_opt_control%ijacobian(1) = scf_control%outer_scf%cdft_opt_control%ijacobian(1) + 1
1501 118 : CALL timestop(handle)
1502 :
1503 372 : END SUBROUTINE qs_calculate_inverse_jacobian
1504 :
1505 : ! **************************************************************************************************
1506 : !> \brief Perform backtracking line search to find the optimal step size for the CDFT constraint
1507 : !> optimizer. Assumes that the CDFT gradient function is a smooth function of the constraint
1508 : !> variables.
1509 : !> \param qs_env the qs_environment_type where to perform the line search
1510 : !> \par History
1511 : !> 02.2017 created [Nico Holmberg]
1512 : ! **************************************************************************************************
1513 186 : SUBROUTINE qs_cdft_line_search(qs_env)
1514 : TYPE(qs_environment_type), POINTER :: qs_env
1515 :
1516 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_cdft_line_search'
1517 :
1518 : CHARACTER(len=default_path_length) :: project_name
1519 : INTEGER :: handle, i, ispin, iter_count, &
1520 : max_linesearch, max_scf, nspins, &
1521 : nsteps, nvar, output_unit, tsteps
1522 : LOGICAL :: continue_ls, continue_ls_exit, converged, do_linesearch, found_solution, &
1523 : reached_maxls, should_exit, should_stop, sign_changed
1524 186 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: positive_sign
1525 : REAL(KIND=dp) :: alpha, alpha_ls, factor, norm_ls
1526 186 : REAL(KIND=dp), DIMENSION(:), POINTER :: energy
1527 186 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient, inv_jacobian
1528 : REAL(KIND=dp), EXTERNAL :: dnrm2
1529 : TYPE(cdft_control_type), POINTER :: cdft_control
1530 : TYPE(cp_logger_type), POINTER :: logger, tmp_logger
1531 186 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv
1532 186 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
1533 : TYPE(dft_control_type), POINTER :: dft_control
1534 186 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1535 : TYPE(mp_para_env_type), POINTER :: para_env
1536 : TYPE(qs_energy_type), POINTER :: energy_qs
1537 : TYPE(qs_ks_env_type), POINTER :: ks_env
1538 : TYPE(qs_rho_type), POINTER :: rho
1539 : TYPE(qs_scf_env_type), POINTER :: scf_env
1540 : TYPE(scf_control_type), POINTER :: scf_control
1541 :
1542 186 : CALL timeset(routineN, handle)
1543 :
1544 186 : NULLIFY (energy, gradient, p_rmpv, rho_ao_kp, mos, rho, &
1545 186 : ks_env, scf_env, scf_control, dft_control, &
1546 186 : cdft_control, inv_jacobian, para_env, &
1547 186 : tmp_logger, energy_qs)
1548 186 : logger => cp_get_default_logger()
1549 :
1550 186 : CPASSERT(ASSOCIATED(qs_env))
1551 : CALL get_qs_env(qs_env, scf_env=scf_env, ks_env=ks_env, &
1552 : scf_control=scf_control, mos=mos, rho=rho, &
1553 : dft_control=dft_control, &
1554 186 : para_env=para_env, energy=energy_qs)
1555 186 : do_linesearch = .FALSE.
1556 186 : SELECT CASE (scf_control%outer_scf%optimizer)
1557 : CASE DEFAULT
1558 : do_linesearch = .FALSE.
1559 : CASE (outer_scf_optimizer_newton_ls)
1560 24 : do_linesearch = .TRUE.
1561 : CASE (outer_scf_optimizer_broyden)
1562 186 : SELECT CASE (scf_control%outer_scf%cdft_opt_control%broyden_type)
1563 : CASE (broyden_type_1, broyden_type_2, broyden_type_1_explicit, broyden_type_2_explicit)
1564 0 : do_linesearch = .FALSE.
1565 : CASE (broyden_type_1_ls, broyden_type_1_explicit_ls, broyden_type_2_ls, broyden_type_2_explicit_ls)
1566 0 : cdft_control => dft_control%qs_control%cdft_control
1567 0 : IF (.NOT. ASSOCIATED(cdft_control)) &
1568 : CALL cp_abort(__LOCATION__, &
1569 : "Optimizers that perform a line search can"// &
1570 0 : " only be used together with a valid CDFT constraint")
1571 0 : IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
1572 24 : do_linesearch = .TRUE.
1573 : END SELECT
1574 : END SELECT
1575 : IF (do_linesearch) THEN
1576 8 : BLOCK
1577 8 : TYPE(mo_set_type), DIMENSION(:), ALLOCATABLE :: mos_ls, mos_stashed
1578 8 : cdft_control => dft_control%qs_control%cdft_control
1579 8 : IF (.NOT. ASSOCIATED(cdft_control)) &
1580 : CALL cp_abort(__LOCATION__, &
1581 : "Optimizers that perform a line search can"// &
1582 0 : " only be used together with a valid CDFT constraint")
1583 8 : CPASSERT(ASSOCIATED(scf_env%outer_scf%inv_jacobian))
1584 8 : CPASSERT(ASSOCIATED(scf_control%outer_scf%cdft_opt_control))
1585 8 : alpha = scf_control%outer_scf%cdft_opt_control%newton_step_save
1586 8 : iter_count = scf_env%outer_scf%iter_count
1587 : ! Redirect output from line search procedure to a new file by creating a temporary logger
1588 8 : project_name = logger%iter_info%project_name
1589 8 : CALL create_tmp_logger(para_env, project_name, "-LineSearch.out", output_unit, tmp_logger)
1590 : ! Save last converged state so we can roll back to it (mo_coeff and some outer_loop variables)
1591 8 : nspins = dft_control%nspins
1592 40 : ALLOCATE (mos_stashed(nspins))
1593 24 : DO ispin = 1, nspins
1594 24 : CALL duplicate_mo_set(mos_stashed(ispin), mos(ispin))
1595 : END DO
1596 8 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1597 8 : p_rmpv => rho_ao_kp(:, 1)
1598 8 : nsteps = cdft_control%total_steps
1599 : ! Allocate work
1600 8 : nvar = SIZE(scf_env%outer_scf%variables, 1)
1601 8 : max_scf = scf_control%outer_scf%max_scf + 1
1602 8 : max_linesearch = scf_control%outer_scf%cdft_opt_control%max_ls
1603 8 : continue_ls = scf_control%outer_scf%cdft_opt_control%continue_ls
1604 8 : factor = scf_control%outer_scf%cdft_opt_control%factor_ls
1605 8 : continue_ls_exit = .FALSE.
1606 8 : found_solution = .FALSE.
1607 32 : ALLOCATE (gradient(nvar, max_scf))
1608 104 : gradient = scf_env%outer_scf%gradient
1609 24 : ALLOCATE (energy(max_scf))
1610 56 : energy = scf_env%outer_scf%energy
1611 8 : reached_maxls = .FALSE.
1612 : ! Broyden optimizers: perform update of inv_jacobian if necessary
1613 8 : IF (scf_control%outer_scf%cdft_opt_control%broyden_update) THEN
1614 0 : CALL outer_loop_optimize(scf_env, scf_control)
1615 : ! Reset the variables and prevent a reupdate of inv_jacobian
1616 0 : scf_env%outer_scf%variables(:, iter_count + 1) = 0
1617 0 : scf_control%outer_scf%cdft_opt_control%broyden_update = .FALSE.
1618 : END IF
1619 : ! Print some info
1620 8 : IF (output_unit > 0) THEN
1621 : WRITE (output_unit, FMT="(/,A)") &
1622 4 : " ================================== LINE SEARCH STARTED =================================="
1623 : WRITE (output_unit, FMT="(A,I5,A)") &
1624 4 : " Evaluating optimal step size for optimizer using a maximum of", max_linesearch, " steps"
1625 4 : IF (continue_ls) THEN
1626 : WRITE (output_unit, FMT="(A)") &
1627 2 : " Line search continues until best step size is found or max steps are reached"
1628 : END IF
1629 : WRITE (output_unit, '(/,A,F5.3)') &
1630 4 : " Initial step size: ", alpha
1631 : WRITE (output_unit, '(/,A,F5.3)') &
1632 4 : " Step size update factor: ", factor
1633 : WRITE (output_unit, '(/,A,I10,A,I10)') &
1634 4 : " Energy evaluation: ", cdft_control%ienergy, ", CDFT SCF iteration: ", iter_count
1635 : END IF
1636 : ! Perform backtracking line search
1637 8 : CALL cp_add_default_logger(tmp_logger)
1638 16 : DO i = 1, max_linesearch
1639 16 : IF (output_unit > 0) THEN
1640 8 : WRITE (output_unit, FMT="(A)") " "
1641 8 : WRITE (output_unit, FMT="(A)") " #####################################"
1642 : WRITE (output_unit, '(A,I10,A)') &
1643 8 : " ### Line search step: ", i, " ###"
1644 8 : WRITE (output_unit, FMT="(A)") " #####################################"
1645 : END IF
1646 16 : inv_jacobian => scf_env%outer_scf%inv_jacobian
1647 : ! Newton update of CDFT variables with a step size of alpha
1648 : scf_env%outer_scf%variables(:, iter_count + 1) = scf_env%outer_scf%variables(:, iter_count) - alpha* &
1649 144 : MATMUL(inv_jacobian, scf_env%outer_scf%gradient(:, iter_count))
1650 : ! With updated CDFT variables, perform SCF
1651 16 : CALL outer_loop_update_qs_env(qs_env, scf_env)
1652 16 : CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.)
1653 16 : CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot)
1654 : CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
1655 16 : converged=converged, should_stop=should_stop, total_scf_steps=tsteps)
1656 16 : CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft)
1657 : ! Update (iter_count + 1) element of gradient and print constraint info
1658 16 : scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count + 1
1659 16 : CALL outer_loop_gradient(qs_env, scf_env)
1660 : CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
1661 : energy_qs, cdft_control%total_steps, &
1662 16 : should_stop=.FALSE., outer_loop_converged=.FALSE., cdft_loop=.FALSE.)
1663 16 : scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count - 1
1664 : ! Store sign of initial gradient for each variable for continue_ls
1665 16 : IF (continue_ls .AND. .NOT. ALLOCATED(positive_sign)) THEN
1666 12 : ALLOCATE (positive_sign(nvar))
1667 8 : DO ispin = 1, nvar
1668 8 : positive_sign(ispin) = scf_env%outer_scf%gradient(ispin, iter_count + 1) >= 0.0_dp
1669 : END DO
1670 : END IF
1671 : ! Check if the L2 norm of the gradient decreased
1672 16 : inv_jacobian => scf_env%outer_scf%inv_jacobian
1673 16 : IF (dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count + 1), 1) < &
1674 : dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count), 1)) THEN
1675 : ! Optimal step size found
1676 14 : IF (.NOT. continue_ls) THEN
1677 : should_exit = .TRUE.
1678 : ELSE
1679 : ! But line search continues for at least one more iteration in an attempt to find a better solution
1680 : ! if max number of steps is not exceeded
1681 10 : IF (found_solution) THEN
1682 : ! Check if the norm also decreased w.r.t. to previously found solution
1683 6 : IF (dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count + 1), 1) > norm_ls) THEN
1684 : ! Norm increased => accept previous solution and exit
1685 : continue_ls_exit = .TRUE.
1686 : END IF
1687 : END IF
1688 : ! Store current state and the value of alpha
1689 10 : IF (.NOT. continue_ls_exit) THEN
1690 10 : should_exit = .FALSE.
1691 10 : alpha_ls = alpha
1692 10 : found_solution = .TRUE.
1693 10 : norm_ls = dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count + 1), 1)
1694 : ! Check if the sign of the gradient has changed for all variables (w.r.t initial gradient)
1695 : ! In this case we should exit because further line search steps will just increase the norm
1696 10 : sign_changed = .TRUE.
1697 20 : DO ispin = 1, nvar
1698 : sign_changed = sign_changed .AND. (positive_sign(ispin) .NEQV. &
1699 28 : scf_env%outer_scf%gradient(ispin, iter_count + 1) >= 0.0_dp)
1700 : END DO
1701 10 : IF (.NOT. ALLOCATED(mos_ls)) THEN
1702 16 : ALLOCATE (mos_ls(nspins))
1703 : ELSE
1704 18 : DO ispin = 1, nspins
1705 18 : CALL deallocate_mo_set(mos_ls(ispin))
1706 : END DO
1707 : END IF
1708 30 : DO ispin = 1, nspins
1709 30 : CALL duplicate_mo_set(mos_ls(ispin), mos(ispin))
1710 : END DO
1711 10 : alpha = alpha*factor
1712 : ! Exit on last iteration
1713 10 : IF (i == max_linesearch) continue_ls_exit = .TRUE.
1714 : ! Exit if constraint target is satisfied to requested tolerance
1715 30 : IF (SQRT(MAXVAL(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count + 1)**2)) < &
1716 : scf_control%outer_scf%eps_scf) &
1717 2 : continue_ls_exit = .TRUE.
1718 : ! Exit if line search jumped over the optimal step length
1719 10 : IF (sign_changed) continue_ls_exit = .TRUE.
1720 : END IF
1721 : END IF
1722 : ELSE
1723 : ! Gradient increased => alpha is too large (if the gradient function is smooth)
1724 2 : should_exit = .FALSE.
1725 : ! Update alpha using Armijo's scheme
1726 2 : alpha = alpha*factor
1727 : END IF
1728 16 : IF (continue_ls_exit) THEN
1729 : ! Continuation of line search did not yield a better alpha, use previously located solution and exit
1730 4 : alpha = alpha_ls
1731 12 : DO ispin = 1, nspins
1732 8 : CALL deallocate_mo_set(mos(ispin))
1733 8 : CALL duplicate_mo_set(mos(ispin), mos_ls(ispin))
1734 : CALL calculate_density_matrix(mos(ispin), &
1735 8 : p_rmpv(ispin)%matrix)
1736 12 : CALL deallocate_mo_set(mos_ls(ispin))
1737 : END DO
1738 4 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1739 4 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1740 4 : DEALLOCATE (mos_ls)
1741 : should_exit = .TRUE.
1742 : END IF
1743 : ! Reached max steps and SCF converged: continue with last iterated step size
1744 12 : IF (.NOT. should_exit .AND. &
1745 : (i == max_linesearch .AND. converged .AND. .NOT. found_solution)) THEN
1746 0 : should_exit = .TRUE.
1747 0 : reached_maxls = .TRUE.
1748 0 : alpha = alpha*(1.0_dp/factor)
1749 : END IF
1750 : ! Reset outer SCF environment to last converged state
1751 32 : scf_env%outer_scf%variables(:, iter_count + 1) = 0.0_dp
1752 208 : scf_env%outer_scf%gradient = gradient
1753 112 : scf_env%outer_scf%energy = energy
1754 : ! Exit line search if a suitable step size was found
1755 16 : IF (should_exit) EXIT
1756 : ! Reset the electronic structure
1757 8 : cdft_control%total_steps = nsteps
1758 24 : DO ispin = 1, nspins
1759 16 : CALL deallocate_mo_set(mos(ispin))
1760 16 : CALL duplicate_mo_set(mos(ispin), mos_stashed(ispin))
1761 : CALL calculate_density_matrix(mos(ispin), &
1762 24 : p_rmpv(ispin)%matrix)
1763 : END DO
1764 8 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1765 24 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1766 : END DO
1767 8 : scf_control%outer_scf%cdft_opt_control%newton_step = alpha
1768 8 : IF (.NOT. should_exit) THEN
1769 : CALL cp_warn(__LOCATION__, &
1770 0 : "Line search did not converge. CDFT SCF proceeds with fixed step size.")
1771 0 : scf_control%outer_scf%cdft_opt_control%newton_step = scf_control%outer_scf%cdft_opt_control%newton_step_save
1772 : END IF
1773 8 : IF (reached_maxls) &
1774 : CALL cp_warn(__LOCATION__, &
1775 0 : "Line search did not converge. CDFT SCF proceeds with lasted iterated step size.")
1776 8 : CALL cp_rm_default_logger()
1777 8 : CALL cp_logger_release(tmp_logger)
1778 : ! Release temporary storage
1779 24 : DO ispin = 1, nspins
1780 24 : CALL deallocate_mo_set(mos_stashed(ispin))
1781 : END DO
1782 8 : DEALLOCATE (mos_stashed, gradient, energy)
1783 8 : IF (ALLOCATED(positive_sign)) DEALLOCATE (positive_sign)
1784 20 : IF (output_unit > 0) THEN
1785 : WRITE (output_unit, FMT="(/,A)") &
1786 4 : " ================================== LINE SEARCH COMPLETE =================================="
1787 4 : CALL close_file(unit_number=output_unit)
1788 : END IF
1789 : END BLOCK
1790 : END IF
1791 :
1792 186 : CALL timestop(handle)
1793 :
1794 186 : END SUBROUTINE qs_cdft_line_search
1795 :
1796 16 : END MODULE qs_scf
|
