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