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