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