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