(git:1f285aa)
qs_scf.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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_operations, ONLY: copy_dbcsr_to_fm,&
51  dbcsr_deallocate_matrix_set
52  USE cp_files, ONLY: close_file
53  USE cp_fm_types, ONLY: cp_fm_create,&
54  cp_fm_release,&
55  cp_fm_to_fm,&
56  cp_fm_type
57  USE cp_log_handling, ONLY: cp_add_default_logger,&
58  cp_get_default_logger,&
59  cp_logger_release,&
60  cp_logger_type,&
61  cp_rm_default_logger,&
62  cp_to_string
63  USE cp_output_handling, ONLY: cp_add_iter_level,&
64  cp_iterate,&
65  cp_p_file,&
66  cp_print_key_should_output,&
67  cp_print_key_unit_nr,&
68  cp_rm_iter_level
69  USE cp_result_methods, ONLY: get_results,&
70  test_for_result
71  USE cp_result_types, ONLY: cp_result_type
72  USE dbcsr_api, ONLY: dbcsr_copy,&
73  dbcsr_deallocate_matrix,&
74  dbcsr_get_info,&
75  dbcsr_init_p,&
76  dbcsr_p_type,&
77  dbcsr_set,&
78  dbcsr_type
79  USE ec_env_types, ONLY: energy_correction_type
80  USE input_constants, ONLY: &
81  broyden_type_1, broyden_type_1_explicit, broyden_type_1_explicit_ls, broyden_type_1_ls, &
82  broyden_type_2, broyden_type_2_explicit, broyden_type_2_explicit_ls, broyden_type_2_ls, &
83  cdft2ot, history_guess, ot2cdft, ot_precond_full_all, ot_precond_full_single, &
84  ot_precond_full_single_inverse, ot_precond_none, ot_precond_s_inverse, &
85  outer_scf_becke_constraint, outer_scf_hirshfeld_constraint, outer_scf_optimizer_broyden, &
86  outer_scf_optimizer_newton_ls
87  USE input_section_types, ONLY: section_vals_get_subs_vals,&
88  section_vals_type
89  USE kinds, ONLY: default_path_length,&
90  default_string_length,&
91  dp
92  USE kpoint_io, ONLY: write_kpoints_restart
93  USE kpoint_types, ONLY: kpoint_type
94  USE machine, ONLY: m_flush,&
95  m_walltime
96  USE mathlib, ONLY: invert_matrix
97  USE message_passing, ONLY: mp_comm_type,&
98  mp_para_env_type
99  USE particle_types, ONLY: particle_type
100  USE preconditioner, ONLY: prepare_preconditioner,&
101  restart_preconditioner
102  USE pw_env_types, ONLY: pw_env_get,&
103  pw_env_type
104  USE pw_pool_types, ONLY: pw_pool_type
105  USE qs_block_davidson_types, ONLY: block_davidson_deallocate
106  USE qs_cdft_scf_utils, ONLY: build_diagonal_jacobian,&
107  create_tmp_logger,&
108  initialize_inverse_jacobian,&
109  prepare_jacobian_stencil,&
110  print_inverse_jacobian,&
111  restart_inverse_jacobian
112  USE qs_cdft_types, ONLY: cdft_control_type
113  USE qs_charges_types, ONLY: qs_charges_type
114  USE qs_density_matrices, ONLY: calculate_density_matrix
115  USE qs_density_mixing_types, ONLY: gspace_mixing_nr
116  USE qs_diis, ONLY: qs_diis_b_clear,&
117  qs_diis_b_clear_kp,&
118  qs_diis_b_create,&
119  qs_diis_b_create_kp
120  USE qs_energy_types, ONLY: qs_energy_type
121  USE qs_environment_types, ONLY: get_qs_env,&
122  qs_environment_type,&
123  set_qs_env
124  USE qs_integrate_potential, ONLY: integrate_v_rspace
125  USE qs_kind_types, ONLY: qs_kind_type
126  USE qs_ks_methods, ONLY: qs_ks_update_qs_env
127  USE qs_ks_types, ONLY: qs_ks_did_change,&
128  qs_ks_env_type
129  USE qs_mo_io, ONLY: write_mo_set_to_restart
130  USE qs_mo_methods, ONLY: make_basis_simple,&
131  make_basis_sm
132  USE qs_mo_occupation, ONLY: set_mo_occupation
133  USE qs_mo_types, ONLY: deallocate_mo_set,&
134  duplicate_mo_set,&
135  get_mo_set,&
136  mo_set_type,&
137  reassign_allocated_mos
138  USE qs_ot, ONLY: qs_ot_new_preconditioner
139  USE qs_ot_scf, ONLY: ot_scf_init,&
140  ot_scf_read_input
141  USE qs_outer_scf, ONLY: outer_loop_gradient,&
142  outer_loop_optimize,&
143  outer_loop_purge_history,&
144  outer_loop_switch,&
145  outer_loop_update_qs_env
146  USE qs_rho_methods, ONLY: qs_rho_update_rho
147  USE qs_rho_types, ONLY: qs_rho_get,&
148  qs_rho_type
149  USE qs_scf_initialization, ONLY: qs_scf_env_initialize
150  USE qs_scf_loop_utils, ONLY: qs_scf_check_inner_exit,&
151  qs_scf_check_outer_exit,&
152  qs_scf_density_mixing,&
153  qs_scf_inner_finalize,&
154  qs_scf_new_mos,&
155  qs_scf_new_mos_kp,&
156  qs_scf_rho_update,&
157  qs_scf_set_loop_flags
158  USE qs_scf_output, ONLY: qs_scf_cdft_info,&
159  qs_scf_cdft_initial_info,&
160  qs_scf_loop_info,&
161  qs_scf_loop_print,&
162  qs_scf_outer_loop_info,&
163  qs_scf_write_mos
164  USE qs_scf_post_scf, ONLY: qs_scf_compute_properties
165  USE qs_scf_types, ONLY: block_davidson_diag_method_nr,&
166  block_krylov_diag_method_nr,&
167  filter_matrix_diag_method_nr,&
168  general_diag_method_nr,&
169  ot_diag_method_nr,&
170  ot_method_nr,&
171  qs_scf_env_type,&
172  special_diag_method_nr
173  USE qs_wf_history_methods, ONLY: wfi_purge_history,&
174  wfi_update
175  USE scf_control_types, ONLY: scf_control_type
176 #include "./base/base_uses.f90"
177 
178  IMPLICIT NONE
179 
180  PRIVATE
181 
182  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf'
183  LOGICAL, PRIVATE :: reuse_precond = .false.
184  LOGICAL, PRIVATE :: used_history = .false.
185 
186  PUBLIC :: scf, scf_env_cleanup, scf_env_do_scf, cdft_scf
187 
188 CONTAINS
189 
190 ! **************************************************************************************************
191 !> \brief perform an scf procedure in the given qs_env
192 !> \param qs_env the qs_environment where to perform the scf procedure
193 !> \param has_converged ...
194 !> \param total_scf_steps ...
195 !> \par History
196 !> 02.2003 introduced scf_env, moved real work to scf_env_do_scf [fawzi]
197 !> \author fawzi
198 !> \note
199 ! **************************************************************************************************
200  SUBROUTINE scf(qs_env, has_converged, total_scf_steps)
201  TYPE(qs_environment_type), POINTER :: qs_env
202  LOGICAL, INTENT(OUT), OPTIONAL :: has_converged
203  INTEGER, INTENT(OUT), OPTIONAL :: total_scf_steps
204 
205  INTEGER :: ihistory, max_scf_tmp, tsteps
206  LOGICAL :: converged, outer_scf_loop, should_stop
207  LOGICAL, SAVE :: first_step_flag = .true.
208  REAL(kind=dp), DIMENSION(:, :), POINTER :: gradient_history, variable_history
209  TYPE(cp_logger_type), POINTER :: logger
210  TYPE(dft_control_type), POINTER :: dft_control
211  TYPE(qs_scf_env_type), POINTER :: scf_env
212  TYPE(scf_control_type), POINTER :: scf_control
213  TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
214 
215  NULLIFY (scf_env)
216  logger => cp_get_default_logger()
217  cpassert(ASSOCIATED(qs_env))
218  IF (PRESENT(has_converged)) THEN
219  has_converged = .false.
220  END IF
221  IF (PRESENT(total_scf_steps)) THEN
222  total_scf_steps = 0
223  END IF
224  CALL get_qs_env(qs_env, scf_env=scf_env, input=input, &
225  dft_control=dft_control, scf_control=scf_control)
226  IF (scf_control%max_scf > 0) THEN
227 
228  dft_section => section_vals_get_subs_vals(input, "DFT")
229  scf_section => section_vals_get_subs_vals(dft_section, "SCF")
230 
231  IF (.NOT. ASSOCIATED(scf_env)) THEN
232  CALL qs_scf_env_initialize(qs_env, scf_env)
233  ! Moved here from qs_scf_env_initialize to be able to have more scf_env
234  CALL set_qs_env(qs_env, scf_env=scf_env)
235  ELSE
236  CALL qs_scf_env_initialize(qs_env, scf_env)
237  END IF
238 
239  IF ((scf_control%density_guess .EQ. history_guess) .AND. (first_step_flag)) THEN
240  max_scf_tmp = scf_control%max_scf
241  scf_control%max_scf = 1
242  outer_scf_loop = scf_control%outer_scf%have_scf
243  scf_control%outer_scf%have_scf = .false.
244  END IF
245 
246  IF (.NOT. dft_control%qs_control%cdft) THEN
247  CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
248  converged=converged, should_stop=should_stop, total_scf_steps=tsteps)
249  ELSE
250  ! Third SCF loop needed for CDFT with OT to properly restart OT inner loop
251  CALL cdft_scf(qs_env=qs_env, should_stop=should_stop)
252  END IF
253 
254  ! If SCF has not converged, then we should not start MP2
255  IF (ASSOCIATED(qs_env%mp2_env)) qs_env%mp2_env%hf_fail = .NOT. converged
256 
257  ! Add the converged outer_scf SCF gradient(s)/variable(s) to history
258  IF (scf_control%outer_scf%have_scf) THEN
259  ihistory = scf_env%outer_scf%iter_count
260  CALL get_qs_env(qs_env, gradient_history=gradient_history, &
261  variable_history=variable_history)
262  ! We only store the latest two values
263  gradient_history(:, 1) = gradient_history(:, 2)
264  gradient_history(:, 2) = scf_env%outer_scf%gradient(:, ihistory)
265  variable_history(:, 1) = variable_history(:, 2)
266  variable_history(:, 2) = scf_env%outer_scf%variables(:, ihistory)
267  ! Reset flag
268  IF (used_history) used_history = .false.
269  ! Update a counter and check if the Jacobian should be deallocated
270  IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
271  scf_control%outer_scf%cdft_opt_control%ijacobian(2) = scf_control%outer_scf%cdft_opt_control%ijacobian(2) + 1
272  IF (scf_control%outer_scf%cdft_opt_control%ijacobian(2) .GE. &
273  scf_control%outer_scf%cdft_opt_control%jacobian_freq(2) .AND. &
274  scf_control%outer_scf%cdft_opt_control%jacobian_freq(2) > 0) &
275  scf_env%outer_scf%deallocate_jacobian = .true.
276  END IF
277  END IF
278  ! *** add the converged wavefunction to the wavefunction history
279  IF ((ASSOCIATED(qs_env%wf_history)) .AND. &
280  ((scf_control%density_guess .NE. history_guess) .OR. &
281  (.NOT. first_step_flag))) THEN
282  IF (.NOT. dft_control%qs_control%cdft) THEN
283  CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
284  ELSE
285  IF (dft_control%qs_control%cdft_control%should_purge) THEN
286  CALL wfi_purge_history(qs_env)
287  CALL outer_loop_purge_history(qs_env)
288  dft_control%qs_control%cdft_control%should_purge = .false.
289  ELSE
290  CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
291  END IF
292  END IF
293  ELSE IF ((scf_control%density_guess .EQ. history_guess) .AND. &
294  (first_step_flag)) THEN
295  scf_control%max_scf = max_scf_tmp
296  scf_control%outer_scf%have_scf = outer_scf_loop
297  first_step_flag = .false.
298  END IF
299 
300  ! *** compute properties that depend on the converged wavefunction
301  IF (.NOT. (should_stop)) CALL qs_scf_compute_properties(qs_env)
302 
303  ! *** cleanup
304  CALL scf_env_cleanup(scf_env)
305  IF (dft_control%qs_control%cdft) &
306  CALL cdft_control_cleanup(dft_control%qs_control%cdft_control)
307 
308  IF (PRESENT(has_converged)) THEN
309  has_converged = converged
310  END IF
311  IF (PRESENT(total_scf_steps)) THEN
312  total_scf_steps = tsteps
313  END IF
314 
315  END IF
316 
317  END SUBROUTINE scf
318 
319 ! **************************************************************************************************
320 !> \brief perform an scf loop
321 !> \param scf_env the scf_env where to perform the scf procedure
322 !> \param scf_control ...
323 !> \param qs_env the qs_env, the scf_env lives in
324 !> \param converged will be true / false if converged is reached
325 !> \param should_stop ...
326 !> \param total_scf_steps ...
327 !> \par History
328 !> long history, see cvs and qs_scf module history
329 !> 02.2003 introduced scf_env [fawzi]
330 !> 09.2005 Frozen density approximation [TdK]
331 !> 06.2007 Check for SCF iteration count early [jgh]
332 !> 10.2019 switch_surf_dip [SGh]
333 !> \author Matthias Krack
334 !> \note
335 ! **************************************************************************************************
336  SUBROUTINE scf_env_do_scf(scf_env, scf_control, qs_env, converged, should_stop, total_scf_steps)
337 
338  TYPE(qs_scf_env_type), POINTER :: scf_env
339  TYPE(scf_control_type), POINTER :: scf_control
340  TYPE(qs_environment_type), POINTER :: qs_env
341  LOGICAL, INTENT(OUT) :: converged, should_stop
342  INTEGER, INTENT(OUT) :: total_scf_steps
343 
344  CHARACTER(LEN=*), PARAMETER :: routinen = 'scf_env_do_scf'
345 
346  CHARACTER(LEN=default_string_length) :: description, name
347  INTEGER :: ext_master_id, handle, handle2, i_tmp, &
348  ic, ispin, iter_count, output_unit, &
349  scf_energy_message_tag, total_steps
350  LOGICAL :: diis_step, do_kpoints, energy_only, exit_inner_loop, exit_outer_loop, &
351  inner_loop_converged, just_energy, outer_loop_converged
352  REAL(kind=dp) :: t1, t2
353  REAL(kind=dp), DIMENSION(3) :: res_val_3
354  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
355  TYPE(cp_logger_type), POINTER :: logger
356  TYPE(cp_result_type), POINTER :: results
357  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
358  TYPE(dft_control_type), POINTER :: dft_control
359  TYPE(energy_correction_type), POINTER :: ec_env
360  TYPE(kpoint_type), POINTER :: kpoints
361  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_last_converged
362  TYPE(mp_comm_type) :: external_comm
363  TYPE(mp_para_env_type), POINTER :: para_env
364  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
365  TYPE(pw_env_type), POINTER :: pw_env
366  TYPE(qs_charges_type), POINTER :: qs_charges
367  TYPE(qs_energy_type), POINTER :: energy
368  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
369  TYPE(qs_ks_env_type), POINTER :: ks_env
370  TYPE(qs_rho_type), POINTER :: rho
371  TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
372 
373  CALL timeset(routinen, handle)
374 
375  NULLIFY (dft_control, rho, energy, &
376  logger, qs_charges, ks_env, mos, atomic_kind_set, qs_kind_set, &
377  particle_set, dft_section, input, &
378  scf_section, para_env, results, kpoints, pw_env, rho_ao_kp, mos_last_converged)
379 
380  cpassert(ASSOCIATED(scf_env))
381  cpassert(ASSOCIATED(qs_env))
382 
383  logger => cp_get_default_logger()
384  t1 = m_walltime()
385 
386  CALL get_qs_env(qs_env=qs_env, &
387  energy=energy, &
388  particle_set=particle_set, &
389  qs_charges=qs_charges, &
390  ks_env=ks_env, &
391  atomic_kind_set=atomic_kind_set, &
392  qs_kind_set=qs_kind_set, &
393  rho=rho, &
394  mos=mos, &
395  input=input, &
396  dft_control=dft_control, &
397  do_kpoints=do_kpoints, &
398  kpoints=kpoints, &
399  results=results, &
400  pw_env=pw_env, &
401  para_env=para_env)
402 
403  CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
404 
405  dft_section => section_vals_get_subs_vals(input, "DFT")
406  scf_section => section_vals_get_subs_vals(dft_section, "SCF")
407 
408  output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", &
409  extension=".scfLog")
410 
411  IF (output_unit > 0) WRITE (unit=output_unit, fmt="(/,/,T2,A)") &
412  "SCF WAVEFUNCTION OPTIMIZATION"
413 
414 ! when switch_surf_dip is switched on, indicate storing mos from the last converged step
415  IF (dft_control%switch_surf_dip) THEN
416  CALL get_qs_env(qs_env, mos_last_converged=mos_last_converged)
417  DO ispin = 1, dft_control%nspins
418  CALL reassign_allocated_mos(mos(ispin), mos_last_converged(ispin))
419  END DO
420  IF (output_unit > 0) WRITE (unit=output_unit, fmt="(/,/,T2,A)") &
421  "COPIED mos_last_converged ---> mos"
422  END IF
423 
424  IF ((output_unit > 0) .AND. (.NOT. scf_control%use_ot)) THEN
425  WRITE (unit=output_unit, &
426  fmt="(/,T3,A,T12,A,T31,A,T39,A,T59,A,T75,A,/,T3,A)") &
427  "Step", "Update method", "Time", "Convergence", "Total energy", "Change", &
428  repeat("-", 78)
429  END IF
430  CALL cp_add_iter_level(logger%iter_info, "QS_SCF")
431 
432  ! check for external communicator and if the intermediate energy should be sent
433  res_val_3(:) = -1.0_dp
434  description = "[EXT_SCF_ENER_COMM]"
435  IF (test_for_result(results, description=description)) THEN
436  CALL get_results(results, description=description, &
437  values=res_val_3, n_entries=i_tmp)
438  cpassert(i_tmp .EQ. 3)
439  IF (all(res_val_3(:) .LE. 0.0)) &
440  CALL cp_abort(__location__, &
441  " Trying to access result ("//trim(description)// &
442  ") which is not correctly stored.")
443  CALL external_comm%set_handle(nint(res_val_3(1)))
444  END IF
445  ext_master_id = nint(res_val_3(2))
446  scf_energy_message_tag = nint(res_val_3(3))
447 
448  ! *** outer loop of the scf, can treat other variables,
449  ! *** such as lagrangian multipliers
450  scf_env%outer_scf%iter_count = 0
451  iter_count = 0
452  total_steps = 0
453  energy%tot_old = 0.0_dp
454 
455  scf_outer_loop: DO
456 
457  CALL init_scf_loop(scf_env=scf_env, qs_env=qs_env, &
458  scf_section=scf_section)
459 
460  CALL qs_scf_set_loop_flags(scf_env, diis_step, &
461  energy_only, just_energy, exit_inner_loop)
462 
463 ! decide whether to switch off dipole correction for convergence purposes
464  dft_control%surf_dip_correct_switch = dft_control%correct_surf_dip
465  IF ((dft_control%correct_surf_dip) .AND. (scf_control%outer_scf%have_scf) .AND. &
466  (scf_env%outer_scf%iter_count > floor(scf_control%outer_scf%max_scf/2.0_dp))) THEN
467  IF (dft_control%switch_surf_dip) THEN
468  dft_control%surf_dip_correct_switch = .false.
469  IF (output_unit > 0) WRITE (unit=output_unit, fmt="(/,/,T2,A)") &
470  "SURFACE DIPOLE CORRECTION switched off"
471  END IF
472  END IF
473  scf_loop: DO
474 
475  CALL timeset(routinen//"_inner_loop", handle2)
476 
477  scf_env%iter_count = scf_env%iter_count + 1
478  iter_count = iter_count + 1
479  CALL cp_iterate(logger%iter_info, last=.false., iter_nr=iter_count)
480 
481  IF (output_unit > 0) CALL m_flush(output_unit)
482 
483  total_steps = total_steps + 1
484  just_energy = energy_only
485 
486  CALL qs_ks_update_qs_env(qs_env, just_energy=just_energy, &
487  calculate_forces=.false.)
488 
489  ! print 'heavy weight' or relatively expensive quantities
490  CALL qs_scf_loop_print(qs_env, scf_env, para_env)
491 
492  IF (do_kpoints) THEN
493  ! kpoints
494  CALL qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step)
495  ELSE
496  ! Gamma points only
497  CALL qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, energy_only)
498  END IF
499 
500  ! Print requested MO information (can be computationally expensive with OT)
501  CALL qs_scf_write_mos(qs_env, scf_env, final_mos=.false.)
502 
503  CALL qs_scf_density_mixing(scf_env, rho, para_env, diis_step)
504 
505  t2 = m_walltime()
506 
507  CALL qs_scf_loop_info(scf_env, output_unit, just_energy, t1, t2, energy)
508 
509  IF (.NOT. just_energy) energy%tot_old = energy%total
510 
511  ! check for external communicator and if the intermediate energy should be sent
512  IF (scf_energy_message_tag .GT. 0) THEN
513  CALL external_comm%send(energy%total, ext_master_id, scf_energy_message_tag)
514  END IF
515 
516  CALL qs_scf_check_inner_exit(qs_env, scf_env, scf_control, should_stop, exit_inner_loop, &
517  inner_loop_converged, output_unit)
518 
519  ! In case we decide to exit we perform few more check to see if this one
520  ! is really the last SCF step
521  IF (exit_inner_loop) THEN
522 
523  CALL qs_scf_inner_finalize(scf_env, qs_env, diis_step, output_unit)
524 
525  CALL qs_scf_check_outer_exit(qs_env, scf_env, scf_control, should_stop, &
526  outer_loop_converged, exit_outer_loop)
527 
528  ! Let's tag the last SCF cycle so we can print informations only of the last step
529  IF (exit_outer_loop) CALL cp_iterate(logger%iter_info, last=.true., iter_nr=iter_count)
530 
531  END IF
532 
533  IF (do_kpoints) THEN
534  CALL write_kpoints_restart(rho_ao_kp, kpoints, scf_env, dft_section, particle_set, qs_kind_set)
535  ELSE
536  ! Write Wavefunction restart file
537  CALL write_mo_set_to_restart(mos, particle_set, dft_section=dft_section, qs_kind_set=qs_kind_set)
538  END IF
539 
540  ! Exit if we have finished with the SCF inner loop
541  IF (exit_inner_loop) THEN
542  CALL timestop(handle2)
543  EXIT scf_loop
544  END IF
545 
546  IF (.NOT. btest(cp_print_key_should_output(logger%iter_info, &
547  scf_section, "PRINT%ITERATION_INFO/TIME_CUMUL"), cp_p_file)) &
548  t1 = m_walltime()
549 
550  ! mixing methods have the new density matrix in p_mix_new
551  IF (scf_env%mixing_method > 0) THEN
552  DO ic = 1, SIZE(rho_ao_kp, 2)
553  DO ispin = 1, dft_control%nspins
554  CALL dbcsr_get_info(rho_ao_kp(ispin, ic)%matrix, name=name) ! keep the name
555  CALL dbcsr_copy(rho_ao_kp(ispin, ic)%matrix, scf_env%p_mix_new(ispin, ic)%matrix, name=name)
556  END DO
557  END DO
558  END IF
559 
560  CALL qs_scf_rho_update(rho, qs_env, scf_env, ks_env, &
561  mix_rho=scf_env%mixing_method >= gspace_mixing_nr)
562 
563  CALL timestop(handle2)
564 
565  END DO scf_loop
566 
567  IF (.NOT. scf_control%outer_scf%have_scf) EXIT scf_outer_loop
568 
569  ! In case we use the OUTER SCF loop let's print some info..
570  CALL qs_scf_outer_loop_info(output_unit, scf_control, scf_env, &
571  energy, total_steps, should_stop, outer_loop_converged)
572 
573 ! save mos to converged mos if outer_loop_converged and surf_dip_correct_switch is true
574  IF (exit_outer_loop) THEN
575  IF ((dft_control%switch_surf_dip) .AND. (outer_loop_converged) .AND. &
576  (dft_control%surf_dip_correct_switch)) THEN
577  DO ispin = 1, dft_control%nspins
578  CALL reassign_allocated_mos(mos_last_converged(ispin), mos(ispin))
579  END DO
580  IF (output_unit > 0) WRITE (unit=output_unit, fmt="(/,/,T2,A)") &
581  "COPIED mos ---> mos_last_converged"
582  END IF
583  END IF
584 
585  IF (exit_outer_loop) EXIT scf_outer_loop
586 
587  !
588  CALL outer_loop_optimize(scf_env, scf_control)
589  CALL outer_loop_update_qs_env(qs_env, scf_env)
590  CALL qs_ks_did_change(ks_env, potential_changed=.true.)
591 
592  END DO scf_outer_loop
593 
594  converged = inner_loop_converged .AND. outer_loop_converged
595  total_scf_steps = total_steps
596 
597  IF (dft_control%qs_control%cdft) &
598  dft_control%qs_control%cdft_control%total_steps = &
599  dft_control%qs_control%cdft_control%total_steps + total_steps
600 
601  IF (.NOT. converged) THEN
602  IF (scf_control%ignore_convergence_failure .OR. should_stop) THEN
603  CALL cp_warn(__location__, "SCF run NOT converged")
604  ELSE
605  CALL cp_abort(__location__, &
606  "SCF run NOT converged. To continue the calculation "// &
607  "regardless, please set the keyword IGNORE_CONVERGENCE_FAILURE.")
608  END IF
609  END IF
610 
611  ! Skip Harris functional calculation if ground-state is NOT converged
612  IF (qs_env%energy_correction) THEN
613  CALL get_qs_env(qs_env, ec_env=ec_env)
614  ec_env%do_skip = .false.
615  IF (ec_env%skip_ec .AND. .NOT. converged) ec_env%do_skip = .true.
616  END IF
617 
618  ! if needed copy mo_coeff dbcsr->fm for later use in post_scf!fm->dbcsr
619  DO ispin = 1, SIZE(mos) !fm -> dbcsr
620  IF (mos(ispin)%use_mo_coeff_b) THEN !fm->dbcsr
621  IF (.NOT. ASSOCIATED(mos(ispin)%mo_coeff_b)) & !fm->dbcsr
622  cpabort("mo_coeff_b is not allocated") !fm->dbcsr
623  CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, & !fm->dbcsr
624  mos(ispin)%mo_coeff) !fm -> dbcsr
625  END IF !fm->dbcsr
626  END DO !fm -> dbcsr
627 
628  CALL cp_rm_iter_level(logger%iter_info, level_name="QS_SCF")
629  CALL timestop(handle)
630 
631  END SUBROUTINE scf_env_do_scf
632 
633 ! **************************************************************************************************
634 !> \brief inits those objects needed if you want to restart the scf with, say
635 !> only a new initial guess, or different density functional or ...
636 !> this will happen just before the scf loop starts
637 !> \param scf_env ...
638 !> \param qs_env ...
639 !> \param scf_section ...
640 !> \par History
641 !> 03.2006 created [Joost VandeVondele]
642 ! **************************************************************************************************
643  SUBROUTINE init_scf_loop(scf_env, qs_env, scf_section)
644 
645  TYPE(qs_scf_env_type), POINTER :: scf_env
646  TYPE(qs_environment_type), POINTER :: qs_env
647  TYPE(section_vals_type), POINTER :: scf_section
648 
649  CHARACTER(LEN=*), PARAMETER :: routinen = 'init_scf_loop'
650 
651  INTEGER :: handle, ispin, nmo, number_of_ot_envs
652  LOGICAL :: do_kpoints, do_rotation, &
653  has_unit_metric, is_full_all
654  TYPE(cp_fm_type), POINTER :: mo_coeff
655  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
656  TYPE(dbcsr_type), POINTER :: orthogonality_metric
657  TYPE(dft_control_type), POINTER :: dft_control
658  TYPE(kpoint_type), POINTER :: kpoints
659  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
660  TYPE(scf_control_type), POINTER :: scf_control
661 
662  CALL timeset(routinen, handle)
663 
664  NULLIFY (scf_control, matrix_s, matrix_ks, dft_control, mos, mo_coeff, kpoints)
665 
666  cpassert(ASSOCIATED(scf_env))
667  cpassert(ASSOCIATED(qs_env))
668 
669  CALL get_qs_env(qs_env=qs_env, &
670  scf_control=scf_control, &
671  dft_control=dft_control, &
672  do_kpoints=do_kpoints, &
673  kpoints=kpoints, &
674  mos=mos)
675 
676  ! if using mo_coeff_b then copy to fm
677  DO ispin = 1, SIZE(mos) !fm->dbcsr
678  IF (mos(1)%use_mo_coeff_b) THEN !fm->dbcsr
679  CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, mos(ispin)%mo_coeff) !fm->dbcsr
680  END IF !fm->dbcsr
681  END DO !fm->dbcsr
682 
683  ! this just guarantees that all mo_occupations match the eigenvalues, if smear
684  DO ispin = 1, dft_control%nspins
685  ! do not reset mo_occupations if the maximum overlap method is in use
686  IF (.NOT. scf_control%diagonalization%mom) &
687  CALL set_mo_occupation(mo_set=mos(ispin), &
688  smear=scf_control%smear)
689  END DO
690 
691  SELECT CASE (scf_env%method)
692  CASE DEFAULT
693 
694  cpabort("unknown scf method method:"//cp_to_string(scf_env%method))
695 
696  CASE (filter_matrix_diag_method_nr)
697 
698  IF (.NOT. scf_env%skip_diis) THEN
699  IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
700  ALLOCATE (scf_env%scf_diis_buffer)
701  CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis)
702  END IF
703  CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
704  END IF
705 
706  CASE (general_diag_method_nr, special_diag_method_nr, block_krylov_diag_method_nr)
707  IF (.NOT. scf_env%skip_diis) THEN
708  IF (do_kpoints) THEN
709  IF (.NOT. ASSOCIATED(kpoints%scf_diis_buffer)) THEN
710  ALLOCATE (kpoints%scf_diis_buffer)
711  CALL qs_diis_b_create_kp(kpoints%scf_diis_buffer, nbuffer=scf_control%max_diis)
712  END IF
713  CALL qs_diis_b_clear_kp(kpoints%scf_diis_buffer)
714  ELSE
715  IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
716  ALLOCATE (scf_env%scf_diis_buffer)
717  CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis)
718  END IF
719  CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
720  END IF
721  END IF
722 
723  CASE (ot_diag_method_nr)
724  CALL get_qs_env(qs_env, matrix_ks=matrix_ks, matrix_s=matrix_s)
725 
726  IF (.NOT. scf_env%skip_diis) THEN
727  IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
728  ALLOCATE (scf_env%scf_diis_buffer)
729  CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis)
730  END IF
731  CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
732  END IF
733 
734  ! disable DFTB and SE for now
735  IF (dft_control%qs_control%dftb .OR. &
736  dft_control%qs_control%xtb .OR. &
737  dft_control%qs_control%semi_empirical) THEN
738  cpabort("DFTB and SE not available with OT/DIAG")
739  END IF
740 
741  ! if an old preconditioner is still around (i.e. outer SCF is active),
742  ! remove it if this could be worthwhile
743  CALL restart_preconditioner(qs_env, scf_env%ot_preconditioner, &
744  scf_control%diagonalization%ot_settings%preconditioner_type, &
745  dft_control%nspins)
746 
747  CALL prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, scf_env%ot_preconditioner, &
748  scf_control%diagonalization%ot_settings%preconditioner_type, &
749  scf_control%diagonalization%ot_settings%precond_solver_type, &
750  scf_control%diagonalization%ot_settings%energy_gap, dft_control%nspins)
751 
752  CASE (block_davidson_diag_method_nr)
753  ! Preconditioner initialized within the loop, when required
754  CASE (ot_method_nr)
755  CALL get_qs_env(qs_env, &
756  has_unit_metric=has_unit_metric, &
757  matrix_s=matrix_s, &
758  matrix_ks=matrix_ks)
759 
760  ! reortho the wavefunctions if we are having an outer scf and
761  ! this is not the first iteration
762  ! this is useful to avoid the build-up of numerical noise
763  ! however, we can not play this trick if restricted (don't mix non-equivalent orbs)
764  IF (scf_control%do_outer_scf_reortho) THEN
765  IF (scf_control%outer_scf%have_scf .AND. .NOT. dft_control%restricted) THEN
766  IF (scf_env%outer_scf%iter_count > 0) THEN
767  DO ispin = 1, dft_control%nspins
768  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
769  IF (has_unit_metric) THEN
770  CALL make_basis_simple(mo_coeff, nmo)
771  ELSE
772  CALL make_basis_sm(mo_coeff, nmo, matrix_s(1)%matrix)
773  END IF
774  END DO
775  END IF
776  END IF
777  ELSE
778  ! dont need any dirty trick for the numerically stable irac algorithm.
779  END IF
780 
781  IF (.NOT. ASSOCIATED(scf_env%qs_ot_env)) THEN
782 
783  ! restricted calculations require just one set of OT orbitals
784  number_of_ot_envs = dft_control%nspins
785  IF (dft_control%restricted) number_of_ot_envs = 1
786 
787  ALLOCATE (scf_env%qs_ot_env(number_of_ot_envs))
788 
789  ! XXX Joost XXX should disentangle reading input from this part
790  CALL ot_scf_read_input(scf_env%qs_ot_env, scf_section)
791 
792  ! keep a note that we are restricted
793  IF (dft_control%restricted) THEN
794  scf_env%qs_ot_env(1)%restricted = .true.
795  ! requires rotation
796  IF (.NOT. scf_env%qs_ot_env(1)%settings%do_rotation) &
797  CALL cp_abort(__location__, &
798  "Restricted calculation with OT requires orbital rotation. Please "// &
799  "activate the OT%ROTATION keyword!")
800  ELSE
801  scf_env%qs_ot_env(:)%restricted = .false.
802  END IF
803 
804  ! this will rotate the MOs to be eigen states, which is not compatible with rotation
805  ! e.g. mo_derivs here do not yet include potentially different occupations numbers
806  do_rotation = scf_env%qs_ot_env(1)%settings%do_rotation
807  ! only full all needs rotation
808  is_full_all = scf_env%qs_ot_env(1)%settings%preconditioner_type == ot_precond_full_all
809  IF (do_rotation .AND. is_full_all) &
810  cpabort('PRECONDITIONER FULL_ALL is not compatible with ROTATION.')
811 
812  ! might need the KS matrix to init properly
813  CALL qs_ks_update_qs_env(qs_env, just_energy=.false., &
814  calculate_forces=.false.)
815 
816  ! if an old preconditioner is still around (i.e. outer SCF is active),
817  ! remove it if this could be worthwhile
818  IF (.NOT. reuse_precond) &
819  CALL restart_preconditioner(qs_env, scf_env%ot_preconditioner, &
820  scf_env%qs_ot_env(1)%settings%preconditioner_type, &
821  dft_control%nspins)
822 
823  !
824  ! preconditioning still needs to be done correctly with has_unit_metric
825  ! notice that a big part of the preconditioning (S^-1) is fine anyhow
826  !
827  IF (has_unit_metric) THEN
828  NULLIFY (orthogonality_metric)
829  ELSE
830  orthogonality_metric => matrix_s(1)%matrix
831  END IF
832 
833  IF (.NOT. reuse_precond) &
834  CALL prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, scf_env%ot_preconditioner, &
835  scf_env%qs_ot_env(1)%settings%preconditioner_type, &
836  scf_env%qs_ot_env(1)%settings%precond_solver_type, &
837  scf_env%qs_ot_env(1)%settings%energy_gap, dft_control%nspins, &
838  has_unit_metric=has_unit_metric, &
839  chol_type=scf_env%qs_ot_env(1)%settings%cholesky_type)
840  IF (reuse_precond) reuse_precond = .false.
841 
842  CALL ot_scf_init(mo_array=mos, matrix_s=orthogonality_metric, &
843  broyden_adaptive_sigma=qs_env%broyden_adaptive_sigma, &
844  qs_ot_env=scf_env%qs_ot_env, matrix_ks=matrix_ks(1)%matrix)
845 
846  SELECT CASE (scf_env%qs_ot_env(1)%settings%preconditioner_type)
847  CASE (ot_precond_none)
848  CASE (ot_precond_full_all, ot_precond_full_single_inverse)
849  DO ispin = 1, SIZE(scf_env%qs_ot_env)
850  CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin), &
851  scf_env%ot_preconditioner(ispin)%preconditioner)
852  END DO
853  CASE (ot_precond_s_inverse, ot_precond_full_single)
854  DO ispin = 1, SIZE(scf_env%qs_ot_env)
855  CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin), &
856  scf_env%ot_preconditioner(1)%preconditioner)
857  END DO
858  CASE DEFAULT
859  DO ispin = 1, SIZE(scf_env%qs_ot_env)
860  CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin), &
861  scf_env%ot_preconditioner(1)%preconditioner)
862  END DO
863  END SELECT
864  END IF
865 
866  ! if we have non-uniform occupations we should be using rotation
867  do_rotation = scf_env%qs_ot_env(1)%settings%do_rotation
868  DO ispin = 1, SIZE(mos)
869  IF (.NOT. mos(ispin)%uniform_occupation) THEN
870  cpassert(do_rotation)
871  END IF
872  END DO
873  END SELECT
874 
875  ! another safety check
876  IF (dft_control%low_spin_roks) THEN
877  cpassert(scf_env%method == ot_method_nr)
878  do_rotation = scf_env%qs_ot_env(1)%settings%do_rotation
879  cpassert(do_rotation)
880  END IF
881 
882  CALL timestop(handle)
883 
884  END SUBROUTINE init_scf_loop
885 
886 ! **************************************************************************************************
887 !> \brief perform cleanup operations (like releasing temporary storage)
888 !> at the end of the scf
889 !> \param scf_env ...
890 !> \par History
891 !> 02.2003 created [fawzi]
892 !> \author fawzi
893 ! **************************************************************************************************
894  SUBROUTINE scf_env_cleanup(scf_env)
895  TYPE(qs_scf_env_type), INTENT(INOUT) :: scf_env
896 
897  CHARACTER(len=*), PARAMETER :: routinen = 'scf_env_cleanup'
898 
899  INTEGER :: handle
900 
901  CALL timeset(routinen, handle)
902 
903  ! Release SCF work storage
904  CALL cp_fm_release(scf_env%scf_work1)
905 
906  IF (ASSOCIATED(scf_env%scf_work2)) THEN
907  CALL cp_fm_release(scf_env%scf_work2)
908  DEALLOCATE (scf_env%scf_work2)
909  END IF
910  IF (ASSOCIATED(scf_env%ortho)) THEN
911  CALL cp_fm_release(scf_env%ortho)
912  DEALLOCATE (scf_env%ortho)
913  END IF
914  IF (ASSOCIATED(scf_env%ortho_m1)) THEN
915  CALL cp_fm_release(scf_env%ortho_m1)
916  DEALLOCATE (scf_env%ortho_m1)
917  END IF
918 
919  IF (ASSOCIATED(scf_env%ortho_dbcsr)) THEN
920  CALL dbcsr_deallocate_matrix(scf_env%ortho_dbcsr)
921  END IF
922  IF (ASSOCIATED(scf_env%buf1_dbcsr)) THEN
923  CALL dbcsr_deallocate_matrix(scf_env%buf1_dbcsr)
924  END IF
925  IF (ASSOCIATED(scf_env%buf2_dbcsr)) THEN
926  CALL dbcsr_deallocate_matrix(scf_env%buf2_dbcsr)
927  END IF
928 
929  IF (ASSOCIATED(scf_env%p_mix_new)) THEN
930  CALL dbcsr_deallocate_matrix_set(scf_env%p_mix_new)
931  END IF
932 
933  IF (ASSOCIATED(scf_env%p_delta)) THEN
934  CALL dbcsr_deallocate_matrix_set(scf_env%p_delta)
935  END IF
936 
937  ! Method dependent cleanup
938  SELECT CASE (scf_env%method)
939  CASE (ot_method_nr)
940  !
941  CASE (ot_diag_method_nr)
942  !
943  CASE (general_diag_method_nr)
944  !
945  CASE (special_diag_method_nr)
946  !
947  CASE (block_krylov_diag_method_nr)
948  CASE (block_davidson_diag_method_nr)
949  CALL block_davidson_deallocate(scf_env%block_davidson_env)
950  CASE (filter_matrix_diag_method_nr)
951  !
952  CASE DEFAULT
953  cpabort("unknown scf method method:"//cp_to_string(scf_env%method))
954  END SELECT
955 
956  IF (ASSOCIATED(scf_env%outer_scf%variables)) THEN
957  DEALLOCATE (scf_env%outer_scf%variables)
958  END IF
959  IF (ASSOCIATED(scf_env%outer_scf%count)) THEN
960  DEALLOCATE (scf_env%outer_scf%count)
961  END IF
962  IF (ASSOCIATED(scf_env%outer_scf%gradient)) THEN
963  DEALLOCATE (scf_env%outer_scf%gradient)
964  END IF
965  IF (ASSOCIATED(scf_env%outer_scf%energy)) THEN
966  DEALLOCATE (scf_env%outer_scf%energy)
967  END IF
968  IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian) .AND. &
969  scf_env%outer_scf%deallocate_jacobian) THEN
970  DEALLOCATE (scf_env%outer_scf%inv_jacobian)
971  END IF
972 
973  CALL timestop(handle)
974 
975  END SUBROUTINE scf_env_cleanup
976 
977 ! **************************************************************************************************
978 !> \brief perform a CDFT scf procedure in the given qs_env
979 !> \param qs_env the qs_environment where to perform the scf procedure
980 !> \param should_stop flag determining if calculation should stop
981 !> \par History
982 !> 12.2015 Created
983 !> \author Nico Holmberg
984 ! **************************************************************************************************
985  SUBROUTINE cdft_scf(qs_env, should_stop)
986  TYPE(qs_environment_type), POINTER :: qs_env
987  LOGICAL, INTENT(OUT) :: should_stop
988 
989  CHARACTER(len=*), PARAMETER :: routinen = 'cdft_scf'
990 
991  INTEGER :: handle, iatom, ispin, ivar, nmo, nvar, &
992  output_unit, tsteps
993  LOGICAL :: cdft_loop_converged, converged, &
994  exit_cdft_loop, first_iteration, &
995  my_uocc, uniform_occupation
996  REAL(kind=dp), DIMENSION(:), POINTER :: mo_occupations
997  TYPE(cdft_control_type), POINTER :: cdft_control
998  TYPE(cp_logger_type), POINTER :: logger
999  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
1000  TYPE(dft_control_type), POINTER :: dft_control
1001  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1002  TYPE(pw_env_type), POINTER :: pw_env
1003  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1004  TYPE(qs_energy_type), POINTER :: energy
1005  TYPE(qs_ks_env_type), POINTER :: ks_env
1006  TYPE(qs_rho_type), POINTER :: rho
1007  TYPE(qs_scf_env_type), POINTER :: scf_env
1008  TYPE(scf_control_type), POINTER :: scf_control
1009  TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
1010 
1011  NULLIFY (scf_env, ks_env, energy, rho, matrix_s, rho_ao, cdft_control, logger, &
1012  dft_control, pw_env, auxbas_pw_pool, energy, ks_env, scf_env, dft_section, &
1013  input, scf_section, scf_control, mos, mo_occupations)
1014  logger => cp_get_default_logger()
1015 
1016  cpassert(ASSOCIATED(qs_env))
1017  CALL get_qs_env(qs_env, scf_env=scf_env, energy=energy, &
1018  dft_control=dft_control, scf_control=scf_control, &
1019  ks_env=ks_env, input=input)
1020 
1021  CALL timeset(routinen//"_loop", handle)
1022  dft_section => section_vals_get_subs_vals(input, "DFT")
1023  scf_section => section_vals_get_subs_vals(dft_section, "SCF")
1024  output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", &
1025  extension=".scfLog")
1026  first_iteration = .true.
1027 
1028  cdft_control => dft_control%qs_control%cdft_control
1029 
1030  scf_env%outer_scf%iter_count = 0
1031  cdft_control%total_steps = 0
1032 
1033  ! Write some info about the CDFT calculation
1034  IF (output_unit > 0) THEN
1035  WRITE (unit=output_unit, fmt="(/,/,T2,A)") &
1036  "CDFT EXTERNAL SCF WAVEFUNCTION OPTIMIZATION"
1037  CALL qs_scf_cdft_initial_info(output_unit, cdft_control)
1038  END IF
1039  IF (cdft_control%reuse_precond) THEN
1040  reuse_precond = .false.
1041  cdft_control%nreused = 0
1042  END IF
1043  cdft_outer_loop: DO
1044  ! Change outer_scf settings to OT settings
1045  CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot)
1046  ! Solve electronic structure with fixed value of constraint
1047  CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
1048  converged=converged, should_stop=should_stop, total_scf_steps=tsteps)
1049  ! Decide whether to reuse the preconditioner on the next iteration
1050  IF (cdft_control%reuse_precond) THEN
1051  ! For convergence in exactly one step, the preconditioner is always reused (assuming max_reuse > 0)
1052  ! usually this means that the electronic structure has already converged to the correct state
1053  ! but the constraint optimizer keeps jumping over the optimal solution
1054  IF (scf_env%outer_scf%iter_count == 1 .AND. scf_env%iter_count == 1 &
1055  .AND. cdft_control%total_steps /= 1) &
1056  cdft_control%nreused = cdft_control%nreused - 1
1057  ! SCF converged in less than precond_freq steps
1058  IF (scf_env%outer_scf%iter_count == 1 .AND. scf_env%iter_count .LE. cdft_control%precond_freq .AND. &
1059  cdft_control%total_steps /= 1 .AND. cdft_control%nreused .LT. cdft_control%max_reuse) THEN
1060  reuse_precond = .true.
1061  cdft_control%nreused = cdft_control%nreused + 1
1062  ELSE
1063  reuse_precond = .false.
1064  cdft_control%nreused = 0
1065  END IF
1066  END IF
1067  ! Update history purging counters
1068  IF (first_iteration .AND. cdft_control%purge_history) THEN
1069  cdft_control%istep = cdft_control%istep + 1
1070  IF (scf_env%outer_scf%iter_count .GT. 1) THEN
1071  cdft_control%nbad_conv = cdft_control%nbad_conv + 1
1072  IF (cdft_control%nbad_conv .GE. cdft_control%purge_freq .AND. &
1073  cdft_control%istep .GE. cdft_control%purge_offset) THEN
1074  cdft_control%nbad_conv = 0
1075  cdft_control%istep = 0
1076  cdft_control%should_purge = .true.
1077  END IF
1078  END IF
1079  END IF
1080  first_iteration = .false.
1081  ! Change outer_scf settings to CDFT settings
1082  CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft)
1083  CALL qs_scf_check_outer_exit(qs_env, scf_env, scf_control, should_stop, &
1084  cdft_loop_converged, exit_cdft_loop)
1085  CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
1086  energy, cdft_control%total_steps, &
1087  should_stop, cdft_loop_converged, cdft_loop=.true.)
1088  IF (exit_cdft_loop) EXIT cdft_outer_loop
1089  ! Check if the inverse Jacobian needs to be calculated
1090  CALL qs_calculate_inverse_jacobian(qs_env)
1091  ! Check if a line search should be performed to find an optimal step size for the optimizer
1092  CALL qs_cdft_line_search(qs_env)
1093  ! Optimize constraint
1094  CALL outer_loop_optimize(scf_env, scf_control)
1095  CALL outer_loop_update_qs_env(qs_env, scf_env)
1096  CALL qs_ks_did_change(ks_env, potential_changed=.true.)
1097  END DO cdft_outer_loop
1098 
1099  cdft_control%ienergy = cdft_control%ienergy + 1
1100 
1101  ! Store needed arrays for ET coupling calculation
1102  IF (cdft_control%do_et) THEN
1103  CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, mos=mos)
1104  nvar = SIZE(cdft_control%target)
1105  ! Matrix representation of weight function
1106  ALLOCATE (cdft_control%wmat(nvar))
1107  DO ivar = 1, nvar
1108  CALL dbcsr_init_p(cdft_control%wmat(ivar)%matrix)
1109  CALL dbcsr_copy(cdft_control%wmat(ivar)%matrix, matrix_s(1)%matrix, &
1110  name="ET_RESTRAINT_MATRIX")
1111  CALL dbcsr_set(cdft_control%wmat(ivar)%matrix, 0.0_dp)
1112  CALL integrate_v_rspace(cdft_control%group(ivar)%weight, &
1113  hmat=cdft_control%wmat(ivar), qs_env=qs_env, &
1114  calculate_forces=.false., &
1115  gapw=dft_control%qs_control%gapw)
1116  END DO
1117  ! Overlap matrix
1118  CALL dbcsr_init_p(cdft_control%matrix_s%matrix)
1119  CALL dbcsr_copy(cdft_control%matrix_s%matrix, matrix_s(1)%matrix, &
1120  name="OVERLAP")
1121  ! Molecular orbital coefficients
1122  NULLIFY (cdft_control%mo_coeff)
1123  ALLOCATE (cdft_control%mo_coeff(dft_control%nspins))
1124  DO ispin = 1, dft_control%nspins
1125  CALL cp_fm_create(matrix=cdft_control%mo_coeff(ispin), &
1126  matrix_struct=qs_env%mos(ispin)%mo_coeff%matrix_struct, &
1127  name="MO_COEFF_A"//trim(adjustl(cp_to_string(ispin)))//"MATRIX")
1128  CALL cp_fm_to_fm(qs_env%mos(ispin)%mo_coeff, &
1129  cdft_control%mo_coeff(ispin))
1130  END DO
1131  ! Density matrix
1132  IF (cdft_control%calculate_metric) THEN
1133  CALL get_qs_env(qs_env, rho=rho)
1134  CALL qs_rho_get(rho, rho_ao=rho_ao)
1135  ALLOCATE (cdft_control%matrix_p(dft_control%nspins))
1136  DO ispin = 1, dft_control%nspins
1137  NULLIFY (cdft_control%matrix_p(ispin)%matrix)
1138  CALL dbcsr_init_p(cdft_control%matrix_p(ispin)%matrix)
1139  CALL dbcsr_copy(cdft_control%matrix_p(ispin)%matrix, rho_ao(ispin)%matrix, &
1140  name="DENSITY MATRIX")
1141  END DO
1142  END IF
1143  ! Copy occupation numbers if non-uniform occupation
1144  uniform_occupation = .true.
1145  DO ispin = 1, dft_control%nspins
1146  CALL get_mo_set(mo_set=mos(ispin), uniform_occupation=my_uocc)
1147  uniform_occupation = uniform_occupation .AND. my_uocc
1148  END DO
1149  IF (.NOT. uniform_occupation) THEN
1150  ALLOCATE (cdft_control%occupations(dft_control%nspins))
1151  DO ispin = 1, dft_control%nspins
1152  CALL get_mo_set(mo_set=mos(ispin), &
1153  nmo=nmo, &
1154  occupation_numbers=mo_occupations)
1155  ALLOCATE (cdft_control%occupations(ispin)%array(nmo))
1156  cdft_control%occupations(ispin)%array(1:nmo) = mo_occupations(1:nmo)
1157  END DO
1158  END IF
1159  END IF
1160 
1161  ! Deallocate constraint storage if forces are not needed
1162  ! In case of a simulation with multiple force_evals,
1163  ! deallocate only if weight function should not be copied to different force_evals
1164  IF (.NOT. (cdft_control%save_pot .OR. cdft_control%transfer_pot)) THEN
1165  CALL get_qs_env(qs_env, pw_env=pw_env)
1166  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1167  DO iatom = 1, SIZE(cdft_control%group)
1168  CALL auxbas_pw_pool%give_back_pw(cdft_control%group(iatom)%weight)
1169  DEALLOCATE (cdft_control%group(iatom)%weight)
1170  END DO
1171  IF (cdft_control%atomic_charges) THEN
1172  DO iatom = 1, cdft_control%natoms
1173  CALL auxbas_pw_pool%give_back_pw(cdft_control%charge(iatom))
1174  END DO
1175  DEALLOCATE (cdft_control%charge)
1176  END IF
1177  IF (cdft_control%type == outer_scf_becke_constraint .AND. &
1178  cdft_control%becke_control%cavity_confine) THEN
1179  IF (.NOT. ASSOCIATED(cdft_control%becke_control%cavity_mat)) THEN
1180  CALL auxbas_pw_pool%give_back_pw(cdft_control%becke_control%cavity)
1181  ELSE
1182  DEALLOCATE (cdft_control%becke_control%cavity_mat)
1183  END IF
1184  ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1185  IF (ASSOCIATED(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)) THEN
1186  CALL auxbas_pw_pool%give_back_pw(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)
1187  END IF
1188  END IF
1189  IF (ASSOCIATED(cdft_control%charges_fragment)) DEALLOCATE (cdft_control%charges_fragment)
1190  cdft_control%need_pot = .true.
1191  cdft_control%external_control = .false.
1192  END IF
1193 
1194  CALL timestop(handle)
1195 
1196  END SUBROUTINE cdft_scf
1197 
1198 ! **************************************************************************************************
1199 !> \brief perform cleanup operations for cdft_control
1200 !> \param cdft_control container for the external CDFT SCF loop variables
1201 !> \par History
1202 !> 12.2015 created [Nico Holmberg]
1203 !> \author Nico Holmberg
1204 ! **************************************************************************************************
1205  SUBROUTINE cdft_control_cleanup(cdft_control)
1206  TYPE(cdft_control_type), POINTER :: cdft_control
1207 
1208  IF (ASSOCIATED(cdft_control%constraint%variables)) &
1209  DEALLOCATE (cdft_control%constraint%variables)
1210  IF (ASSOCIATED(cdft_control%constraint%count)) &
1211  DEALLOCATE (cdft_control%constraint%count)
1212  IF (ASSOCIATED(cdft_control%constraint%gradient)) &
1213  DEALLOCATE (cdft_control%constraint%gradient)
1214  IF (ASSOCIATED(cdft_control%constraint%energy)) &
1215  DEALLOCATE (cdft_control%constraint%energy)
1216  IF (ASSOCIATED(cdft_control%constraint%inv_jacobian) .AND. &
1217  cdft_control%constraint%deallocate_jacobian) &
1218  DEALLOCATE (cdft_control%constraint%inv_jacobian)
1219 
1220  END SUBROUTINE cdft_control_cleanup
1221 
1222 ! **************************************************************************************************
1223 !> \brief Calculates the finite difference inverse Jacobian
1224 !> \param qs_env the qs_environment_type where to compute the Jacobian
1225 !> \par History
1226 !> 01.2017 created [Nico Holmberg]
1227 ! **************************************************************************************************
1228  SUBROUTINE qs_calculate_inverse_jacobian(qs_env)
1229  TYPE(qs_environment_type), POINTER :: qs_env
1230 
1231  CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_calculate_inverse_jacobian'
1232 
1233  CHARACTER(len=default_path_length) :: project_name
1234  INTEGER :: counter, handle, i, ispin, iter_count, &
1235  iwork, j, max_scf, nspins, nsteps, &
1236  nvar, nwork, output_unit, pwork, &
1237  tsteps, twork
1238  LOGICAL :: converged, explicit_jacobian, &
1239  should_build, should_stop, &
1240  use_md_history
1241  REAL(kind=dp) :: inv_error, step_size
1242  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coeff, dh, step_multiplier
1243  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: jacobian
1244  REAL(kind=dp), DIMENSION(:), POINTER :: energy
1245  REAL(kind=dp), DIMENSION(:, :), POINTER :: gradient, inv_jacobian
1246  TYPE(cdft_control_type), POINTER :: cdft_control
1247  TYPE(cp_logger_type), POINTER :: logger, tmp_logger
1248  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv
1249  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
1250  TYPE(dft_control_type), POINTER :: dft_control
1251  TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mos_stashed
1252  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1253  TYPE(mp_para_env_type), POINTER :: para_env
1254  TYPE(qs_energy_type), POINTER :: energy_qs
1255  TYPE(qs_ks_env_type), POINTER :: ks_env
1256  TYPE(qs_rho_type), POINTER :: rho
1257  TYPE(qs_scf_env_type), POINTER :: scf_env
1258  TYPE(scf_control_type), POINTER :: scf_control
1259 
1260  NULLIFY (energy, gradient, p_rmpv, rho_ao_kp, mos, rho, &
1261  ks_env, scf_env, scf_control, dft_control, cdft_control, &
1262  inv_jacobian, para_env, tmp_logger, energy_qs)
1263  logger => cp_get_default_logger()
1264 
1265  cpassert(ASSOCIATED(qs_env))
1266  CALL get_qs_env(qs_env, scf_env=scf_env, ks_env=ks_env, &
1267  scf_control=scf_control, mos=mos, rho=rho, &
1268  dft_control=dft_control, &
1269  para_env=para_env, energy=energy_qs)
1270  explicit_jacobian = .false.
1271  should_build = .false.
1272  use_md_history = .false.
1273  iter_count = scf_env%outer_scf%iter_count
1274  ! Quick exit if optimizer does not require Jacobian
1275  IF (.NOT. ASSOCIATED(scf_control%outer_scf%cdft_opt_control)) RETURN
1276  ! Check if Jacobian should be calculated and initialize
1277  CALL timeset(routinen, handle)
1278  CALL initialize_inverse_jacobian(scf_control, scf_env, explicit_jacobian, should_build, used_history)
1279  IF (scf_control%outer_scf%cdft_opt_control%jacobian_restart) THEN
1280  ! Restart from previously calculated inverse Jacobian
1281  should_build = .false.
1282  CALL restart_inverse_jacobian(qs_env)
1283  END IF
1284  IF (should_build) THEN
1285  scf_env%outer_scf%deallocate_jacobian = .false.
1286  ! Actually need to (re)build the Jacobian
1287  IF (explicit_jacobian) THEN
1288  ! Build Jacobian with finite differences
1289  cdft_control => dft_control%qs_control%cdft_control
1290  IF (.NOT. ASSOCIATED(cdft_control)) &
1291  CALL cp_abort(__location__, &
1292  "Optimizers that need the explicit Jacobian can"// &
1293  " only be used together with a valid CDFT constraint.")
1294  ! Redirect output from Jacobian calculation to a new file by creating a temporary logger
1295  project_name = logger%iter_info%project_name
1296  CALL create_tmp_logger(para_env, project_name, "-JacobianInfo.out", output_unit, tmp_logger)
1297  ! Save last converged state so we can roll back to it (mo_coeff and some outer_loop variables)
1298  nspins = dft_control%nspins
1299  ALLOCATE (mos_stashed(nspins))
1300  DO ispin = 1, nspins
1301  CALL duplicate_mo_set(mos_stashed(ispin), mos(ispin))
1302  END DO
1303  CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1304  p_rmpv => rho_ao_kp(:, 1)
1305  ! Allocate work
1306  nvar = SIZE(scf_env%outer_scf%variables, 1)
1307  max_scf = scf_control%outer_scf%max_scf + 1
1308  ALLOCATE (gradient(nvar, max_scf))
1309  gradient = scf_env%outer_scf%gradient
1310  ALLOCATE (energy(max_scf))
1311  energy = scf_env%outer_scf%energy
1312  ALLOCATE (jacobian(nvar, nvar))
1313  jacobian = 0.0_dp
1314  nsteps = cdft_control%total_steps
1315  ! Setup finite difference scheme
1316  CALL prepare_jacobian_stencil(qs_env, output_unit, nwork, pwork, coeff, step_multiplier, dh)
1317  twork = pwork - nwork
1318  DO i = 1, nvar
1319  jacobian(i, :) = coeff(0)*scf_env%outer_scf%gradient(i, iter_count)
1320  END DO
1321  ! Calculate the Jacobian by perturbing each Lagrangian and recalculating the energy self-consistently
1322  CALL cp_add_default_logger(tmp_logger)
1323  DO i = 1, nvar
1324  IF (output_unit > 0) THEN
1325  WRITE (output_unit, fmt="(A)") " "
1326  WRITE (output_unit, fmt="(A)") " #####################################"
1327  WRITE (output_unit, '(A,I3,A,I3,A)') &
1328  " ### Constraint ", i, " of ", nvar, " ###"
1329  WRITE (output_unit, fmt="(A)") " #####################################"
1330  END IF
1331  counter = 0
1332  DO iwork = nwork, pwork
1333  IF (iwork == 0) cycle
1334  counter = counter + 1
1335  IF (output_unit > 0) THEN
1336  WRITE (output_unit, fmt="(A)") " #####################################"
1337  WRITE (output_unit, '(A,I3,A,I3,A)') &
1338  " ### Energy evaluation ", counter, " of ", twork, " ###"
1339  WRITE (output_unit, fmt="(A)") " #####################################"
1340  END IF
1341  IF (SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step) == 1) THEN
1342  step_size = scf_control%outer_scf%cdft_opt_control%jacobian_step(1)
1343  ELSE
1344  step_size = scf_control%outer_scf%cdft_opt_control%jacobian_step(i)
1345  END IF
1346  scf_env%outer_scf%variables(:, iter_count + 1) = scf_env%outer_scf%variables(:, iter_count)
1347  scf_env%outer_scf%variables(i, iter_count + 1) = scf_env%outer_scf%variables(i, iter_count) + &
1348  step_multiplier(iwork)*step_size
1349  CALL outer_loop_update_qs_env(qs_env, scf_env)
1350  CALL qs_ks_did_change(ks_env, potential_changed=.true.)
1351  CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot)
1352  CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
1353  converged=converged, should_stop=should_stop, total_scf_steps=tsteps)
1354  CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft)
1355  ! Update (iter_count + 1) element of gradient and print constraint info
1356  scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count + 1
1357  CALL outer_loop_gradient(qs_env, scf_env)
1358  CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
1359  energy_qs, cdft_control%total_steps, &
1360  should_stop=.false., outer_loop_converged=.false., cdft_loop=.false.)
1361  scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count - 1
1362  ! Update Jacobian
1363  DO j = 1, nvar
1364  jacobian(j, i) = jacobian(j, i) + coeff(iwork)*scf_env%outer_scf%gradient(j, iter_count + 1)
1365  END DO
1366  ! Reset everything to last converged state
1367  scf_env%outer_scf%variables(:, iter_count + 1) = 0.0_dp
1368  scf_env%outer_scf%gradient = gradient
1369  scf_env%outer_scf%energy = energy
1370  cdft_control%total_steps = nsteps
1371  DO ispin = 1, nspins
1372  CALL deallocate_mo_set(mos(ispin))
1373  CALL duplicate_mo_set(mos(ispin), mos_stashed(ispin))
1374  CALL calculate_density_matrix(mos(ispin), &
1375  p_rmpv(ispin)%matrix)
1376  END DO
1377  CALL qs_rho_update_rho(rho, qs_env=qs_env)
1378  CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
1379  END DO
1380  END DO
1381  CALL cp_rm_default_logger()
1382  CALL cp_logger_release(tmp_logger)
1383  ! Finalize and invert Jacobian
1384  DO j = 1, nvar
1385  DO i = 1, nvar
1386  jacobian(i, j) = jacobian(i, j)/dh(j)
1387  END DO
1388  END DO
1389  IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
1390  ALLOCATE (scf_env%outer_scf%inv_jacobian(nvar, nvar))
1391  inv_jacobian => scf_env%outer_scf%inv_jacobian
1392  CALL invert_matrix(jacobian, inv_jacobian, inv_error)
1393  scf_control%outer_scf%cdft_opt_control%broyden_update = .false.
1394  ! Release temporary storage
1395  DO ispin = 1, nspins
1396  CALL deallocate_mo_set(mos_stashed(ispin))
1397  END DO
1398  DEALLOCATE (mos_stashed, jacobian, gradient, energy, coeff, step_multiplier, dh)
1399  IF (output_unit > 0) THEN
1400  WRITE (output_unit, fmt="(/,A)") &
1401  " ================================== JACOBIAN CALCULATED =================================="
1402  CALL close_file(unit_number=output_unit)
1403  END IF
1404  ELSE
1405  ! Build a strictly diagonal Jacobian from history and invert it
1406  CALL build_diagonal_jacobian(qs_env, used_history)
1407  END IF
1408  END IF
1409  IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian) .AND. para_env%is_source()) THEN
1410  ! Write restart file for inverse Jacobian
1411  CALL print_inverse_jacobian(logger, scf_env%outer_scf%inv_jacobian, iter_count)
1412  END IF
1413  ! Update counter
1414  scf_control%outer_scf%cdft_opt_control%ijacobian(1) = scf_control%outer_scf%cdft_opt_control%ijacobian(1) + 1
1415  CALL timestop(handle)
1416 
1417  END SUBROUTINE qs_calculate_inverse_jacobian
1418 
1419 ! **************************************************************************************************
1420 !> \brief Perform backtracking line search to find the optimal step size for the CDFT constraint
1421 !> optimizer. Assumes that the CDFT gradient function is a smooth function of the constraint
1422 !> variables.
1423 !> \param qs_env the qs_environment_type where to perform the line search
1424 !> \par History
1425 !> 02.2017 created [Nico Holmberg]
1426 ! **************************************************************************************************
1427  SUBROUTINE qs_cdft_line_search(qs_env)
1428  TYPE(qs_environment_type), POINTER :: qs_env
1429 
1430  CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_cdft_line_search'
1431 
1432  CHARACTER(len=default_path_length) :: project_name
1433  INTEGER :: handle, i, ispin, iter_count, &
1434  max_linesearch, max_scf, nspins, &
1435  nsteps, nvar, output_unit, tsteps
1436  LOGICAL :: continue_ls, continue_ls_exit, converged, do_linesearch, found_solution, &
1437  reached_maxls, should_exit, should_stop, sign_changed
1438  LOGICAL, ALLOCATABLE, DIMENSION(:) :: positive_sign
1439  REAL(kind=dp) :: alpha, alpha_ls, factor, norm_ls
1440  REAL(kind=dp), DIMENSION(:), POINTER :: energy
1441  REAL(kind=dp), DIMENSION(:, :), POINTER :: gradient, inv_jacobian
1442  REAL(kind=dp), EXTERNAL :: dnrm2
1443  TYPE(cdft_control_type), POINTER :: cdft_control
1444  TYPE(cp_logger_type), POINTER :: logger, tmp_logger
1445  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv
1446  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
1447  TYPE(dft_control_type), POINTER :: dft_control
1448  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1449  TYPE(mp_para_env_type), POINTER :: para_env
1450  TYPE(qs_energy_type), POINTER :: energy_qs
1451  TYPE(qs_ks_env_type), POINTER :: ks_env
1452  TYPE(qs_rho_type), POINTER :: rho
1453  TYPE(qs_scf_env_type), POINTER :: scf_env
1454  TYPE(scf_control_type), POINTER :: scf_control
1455 
1456  CALL timeset(routinen, handle)
1457 
1458  NULLIFY (energy, gradient, p_rmpv, rho_ao_kp, mos, rho, &
1459  ks_env, scf_env, scf_control, dft_control, &
1460  cdft_control, inv_jacobian, para_env, &
1461  tmp_logger, energy_qs)
1462  logger => cp_get_default_logger()
1463 
1464  cpassert(ASSOCIATED(qs_env))
1465  CALL get_qs_env(qs_env, scf_env=scf_env, ks_env=ks_env, &
1466  scf_control=scf_control, mos=mos, rho=rho, &
1467  dft_control=dft_control, &
1468  para_env=para_env, energy=energy_qs)
1469  do_linesearch = .false.
1470  SELECT CASE (scf_control%outer_scf%optimizer)
1471  CASE DEFAULT
1472  do_linesearch = .false.
1473  CASE (outer_scf_optimizer_newton_ls)
1474  do_linesearch = .true.
1475  CASE (outer_scf_optimizer_broyden)
1476  SELECT CASE (scf_control%outer_scf%cdft_opt_control%broyden_type)
1477  CASE (broyden_type_1, broyden_type_2, broyden_type_1_explicit, broyden_type_2_explicit)
1478  do_linesearch = .false.
1479  CASE (broyden_type_1_ls, broyden_type_1_explicit_ls, broyden_type_2_ls, broyden_type_2_explicit_ls)
1480  cdft_control => dft_control%qs_control%cdft_control
1481  IF (.NOT. ASSOCIATED(cdft_control)) &
1482  CALL cp_abort(__location__, &
1483  "Optimizers that perform a line search can"// &
1484  " only be used together with a valid CDFT constraint")
1485  IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
1486  do_linesearch = .true.
1487  END SELECT
1488  END SELECT
1489  IF (do_linesearch) THEN
1490  block
1491  TYPE(mo_set_type), DIMENSION(:), ALLOCATABLE :: mos_ls, mos_stashed
1492  cdft_control => dft_control%qs_control%cdft_control
1493  IF (.NOT. ASSOCIATED(cdft_control)) &
1494  CALL cp_abort(__location__, &
1495  "Optimizers that perform a line search can"// &
1496  " only be used together with a valid CDFT constraint")
1497  cpassert(ASSOCIATED(scf_env%outer_scf%inv_jacobian))
1498  cpassert(ASSOCIATED(scf_control%outer_scf%cdft_opt_control))
1499  alpha = scf_control%outer_scf%cdft_opt_control%newton_step_save
1500  iter_count = scf_env%outer_scf%iter_count
1501  ! Redirect output from line search procedure to a new file by creating a temporary logger
1502  project_name = logger%iter_info%project_name
1503  CALL create_tmp_logger(para_env, project_name, "-LineSearch.out", output_unit, tmp_logger)
1504  ! Save last converged state so we can roll back to it (mo_coeff and some outer_loop variables)
1505  nspins = dft_control%nspins
1506  ALLOCATE (mos_stashed(nspins))
1507  DO ispin = 1, nspins
1508  CALL duplicate_mo_set(mos_stashed(ispin), mos(ispin))
1509  END DO
1510  CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1511  p_rmpv => rho_ao_kp(:, 1)
1512  nsteps = cdft_control%total_steps
1513  ! Allocate work
1514  nvar = SIZE(scf_env%outer_scf%variables, 1)
1515  max_scf = scf_control%outer_scf%max_scf + 1
1516  max_linesearch = scf_control%outer_scf%cdft_opt_control%max_ls
1517  continue_ls = scf_control%outer_scf%cdft_opt_control%continue_ls
1518  factor = scf_control%outer_scf%cdft_opt_control%factor_ls
1519  continue_ls_exit = .false.
1520  found_solution = .false.
1521  ALLOCATE (gradient(nvar, max_scf))
1522  gradient = scf_env%outer_scf%gradient
1523  ALLOCATE (energy(max_scf))
1524  energy = scf_env%outer_scf%energy
1525  reached_maxls = .false.
1526  ! Broyden optimizers: perform update of inv_jacobian if necessary
1527  IF (scf_control%outer_scf%cdft_opt_control%broyden_update) THEN
1528  CALL outer_loop_optimize(scf_env, scf_control)
1529  ! Reset the variables and prevent a reupdate of inv_jacobian
1530  scf_env%outer_scf%variables(:, iter_count + 1) = 0
1531  scf_control%outer_scf%cdft_opt_control%broyden_update = .false.
1532  END IF
1533  ! Print some info
1534  IF (output_unit > 0) THEN
1535  WRITE (output_unit, fmt="(/,A)") &
1536  " ================================== LINE SEARCH STARTED =================================="
1537  WRITE (output_unit, fmt="(A,I5,A)") &
1538  " Evaluating optimal step size for optimizer using a maximum of", max_linesearch, " steps"
1539  IF (continue_ls) THEN
1540  WRITE (output_unit, fmt="(A)") &
1541  " Line search continues until best step size is found or max steps are reached"
1542  END IF
1543  WRITE (output_unit, '(/,A,F5.3)') &
1544  " Initial step size: ", alpha
1545  WRITE (output_unit, '(/,A,F5.3)') &
1546  " Step size update factor: ", factor
1547  WRITE (output_unit, '(/,A,I10,A,I10)') &
1548  " Energy evaluation: ", cdft_control%ienergy, ", CDFT SCF iteration: ", iter_count
1549  END IF
1550  ! Perform backtracking line search
1551  CALL cp_add_default_logger(tmp_logger)
1552  DO i = 1, max_linesearch
1553  IF (output_unit > 0) THEN
1554  WRITE (output_unit, fmt="(A)") " "
1555  WRITE (output_unit, fmt="(A)") " #####################################"
1556  WRITE (output_unit, '(A,I10,A)') &
1557  " ### Line search step: ", i, " ###"
1558  WRITE (output_unit, fmt="(A)") " #####################################"
1559  END IF
1560  inv_jacobian => scf_env%outer_scf%inv_jacobian
1561  ! Newton update of CDFT variables with a step size of alpha
1562  scf_env%outer_scf%variables(:, iter_count + 1) = scf_env%outer_scf%variables(:, iter_count) - alpha* &
1563  matmul(inv_jacobian, scf_env%outer_scf%gradient(:, iter_count))
1564  ! With updated CDFT variables, perform SCF
1565  CALL outer_loop_update_qs_env(qs_env, scf_env)
1566  CALL qs_ks_did_change(ks_env, potential_changed=.true.)
1567  CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot)
1568  CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
1569  converged=converged, should_stop=should_stop, total_scf_steps=tsteps)
1570  CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft)
1571  ! Update (iter_count + 1) element of gradient and print constraint info
1572  scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count + 1
1573  CALL outer_loop_gradient(qs_env, scf_env)
1574  CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
1575  energy_qs, cdft_control%total_steps, &
1576  should_stop=.false., outer_loop_converged=.false., cdft_loop=.false.)
1577  scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count - 1
1578  ! Store sign of initial gradient for each variable for continue_ls
1579  IF (continue_ls .AND. .NOT. ALLOCATED(positive_sign)) THEN
1580  ALLOCATE (positive_sign(nvar))
1581  DO ispin = 1, nvar
1582  positive_sign(ispin) = scf_env%outer_scf%gradient(ispin, iter_count + 1) .GE. 0.0_dp
1583  END DO
1584  END IF
1585  ! Check if the L2 norm of the gradient decreased
1586  inv_jacobian => scf_env%outer_scf%inv_jacobian
1587  IF (dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count + 1), 1) .LT. &
1588  dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count), 1)) THEN
1589  ! Optimal step size found
1590  IF (.NOT. continue_ls) THEN
1591  should_exit = .true.
1592  ELSE
1593  ! But line search continues for at least one more iteration in an attempt to find a better solution
1594  ! if max number of steps is not exceeded
1595  IF (found_solution) THEN
1596  ! Check if the norm also decreased w.r.t. to previously found solution
1597  IF (dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count + 1), 1) .GT. norm_ls) THEN
1598  ! Norm increased => accept previous solution and exit
1599  continue_ls_exit = .true.
1600  END IF
1601  END IF
1602  ! Store current state and the value of alpha
1603  IF (.NOT. continue_ls_exit) THEN
1604  should_exit = .false.
1605  alpha_ls = alpha
1606  found_solution = .true.
1607  norm_ls = dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count + 1), 1)
1608  ! Check if the sign of the gradient has changed for all variables (w.r.t initial gradient)
1609  ! In this case we should exit because further line search steps will just increase the norm
1610  sign_changed = .true.
1611  DO ispin = 1, nvar
1612  sign_changed = sign_changed .AND. (positive_sign(ispin) .NEQV. &
1613  scf_env%outer_scf%gradient(ispin, iter_count + 1) .GE. 0.0_dp)
1614  END DO
1615  IF (.NOT. ALLOCATED(mos_ls)) THEN
1616  ALLOCATE (mos_ls(nspins))
1617  ELSE
1618  DO ispin = 1, nspins
1619  CALL deallocate_mo_set(mos_ls(ispin))
1620  END DO
1621  END IF
1622  DO ispin = 1, nspins
1623  CALL duplicate_mo_set(mos_ls(ispin), mos(ispin))
1624  END DO
1625  alpha = alpha*factor
1626  ! Exit on last iteration
1627  IF (i == max_linesearch) continue_ls_exit = .true.
1628  ! Exit if constraint target is satisfied to requested tolerance
1629  IF (sqrt(maxval(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count + 1)**2)) .LT. &
1630  scf_control%outer_scf%eps_scf) &
1631  continue_ls_exit = .true.
1632  ! Exit if line search jumped over the optimal step length
1633  IF (sign_changed) continue_ls_exit = .true.
1634  END IF
1635  END IF
1636  ELSE
1637  ! Gradient increased => alpha is too large (if the gradient function is smooth)
1638  should_exit = .false.
1639  ! Update alpha using Armijo's scheme
1640  alpha = alpha*factor
1641  END IF
1642  IF (continue_ls_exit) THEN
1643  ! Continuation of line search did not yield a better alpha, use previously located solution and exit
1644  alpha = alpha_ls
1645  DO ispin = 1, nspins
1646  CALL deallocate_mo_set(mos(ispin))
1647  CALL duplicate_mo_set(mos(ispin), mos_ls(ispin))
1648  CALL calculate_density_matrix(mos(ispin), &
1649  p_rmpv(ispin)%matrix)
1650  CALL deallocate_mo_set(mos_ls(ispin))
1651  END DO
1652  CALL qs_rho_update_rho(rho, qs_env=qs_env)
1653  CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
1654  DEALLOCATE (mos_ls)
1655  should_exit = .true.
1656  END IF
1657  ! Reached max steps and SCF converged: continue with last iterated step size
1658  IF (.NOT. should_exit .AND. &
1659  (i == max_linesearch .AND. converged .AND. .NOT. found_solution)) THEN
1660  should_exit = .true.
1661  reached_maxls = .true.
1662  alpha = alpha*(1.0_dp/factor)
1663  END IF
1664  ! Reset outer SCF environment to last converged state
1665  scf_env%outer_scf%variables(:, iter_count + 1) = 0.0_dp
1666  scf_env%outer_scf%gradient = gradient
1667  scf_env%outer_scf%energy = energy
1668  ! Exit line search if a suitable step size was found
1669  IF (should_exit) EXIT
1670  ! Reset the electronic structure
1671  cdft_control%total_steps = nsteps
1672  DO ispin = 1, nspins
1673  CALL deallocate_mo_set(mos(ispin))
1674  CALL duplicate_mo_set(mos(ispin), mos_stashed(ispin))
1675  CALL calculate_density_matrix(mos(ispin), &
1676  p_rmpv(ispin)%matrix)
1677  END DO
1678  CALL qs_rho_update_rho(rho, qs_env=qs_env)
1679  CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
1680  END DO
1681  scf_control%outer_scf%cdft_opt_control%newton_step = alpha
1682  IF (.NOT. should_exit) THEN
1683  CALL cp_warn(__location__, &
1684  "Line search did not converge. CDFT SCF proceeds with fixed step size.")
1685  scf_control%outer_scf%cdft_opt_control%newton_step = scf_control%outer_scf%cdft_opt_control%newton_step_save
1686  END IF
1687  IF (reached_maxls) &
1688  CALL cp_warn(__location__, &
1689  "Line search did not converge. CDFT SCF proceeds with lasted iterated step size.")
1690  CALL cp_rm_default_logger()
1691  CALL cp_logger_release(tmp_logger)
1692  ! Release temporary storage
1693  DO ispin = 1, nspins
1694  CALL deallocate_mo_set(mos_stashed(ispin))
1695  END DO
1696  DEALLOCATE (mos_stashed, gradient, energy)
1697  IF (ALLOCATED(positive_sign)) DEALLOCATE (positive_sign)
1698  IF (output_unit > 0) THEN
1699  WRITE (output_unit, fmt="(/,A)") &
1700  " ================================== LINE SEARCH COMPLETE =================================="
1701  CALL close_file(unit_number=output_unit)
1702  END IF
1703  END block
1704  END IF
1705 
1706  CALL timestop(handle)
1707 
1708  END SUBROUTINE qs_cdft_line_search
1709 
1710 END MODULE qs_scf
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition: cp_dbcsr_operations.F:208
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition: cp_files.F:119
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_rm_default_logger
subroutine, public cp_rm_default_logger()
the cousin of cp_add_default_logger, decrements the stack, so that the default logger is what it has ...
Definition: cp_log_handling.F:213
cp_log_handling::cp_logger_release
subroutine, public cp_logger_release(logger)
releases this logger
Definition: cp_log_handling.F:423
cp_log_handling::cp_add_default_logger
subroutine, public cp_add_default_logger(logger)
adds a default logger. MUST be called before logging occours
Definition: cp_log_handling.F:188
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition: cp_output_handling.F:103
cp_output_handling::cp_iterate
subroutine, public cp_iterate(iteration_info, last, iter_nr, increment, iter_nr_out)
adds one to the actual iteration
Definition: cp_output_handling.F:574
cp_output_handling::cp_rm_iter_level
subroutine, public cp_rm_iter_level(iteration_info, level_name, n_rlevel_att)
Removes an iteration level.
Definition: cp_output_handling.F:655
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition: cp_output_handling.F:345
cp_output_handling::cp_add_iter_level
subroutine, public cp_add_iter_level(iteration_info, level_name, n_rlevel_new)
Adds an iteration level.
Definition: cp_output_handling.F:612
cp_result_methods
set of type/routines to handle the storage of results in force_envs
Definition: cp_result_methods.F:17
cp_result_methods::test_for_result
logical function, public test_for_result(results, description)
test for a certain result in the result_list
Definition: cp_result_methods.F:127
cp_result_types
set of type/routines to handle the storage of results in force_envs
Definition: cp_result_types.F:18
ec_env_types
Types needed for a for a Energy Correction.
Definition: ec_env_types.F:14
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::broyden_type_2_explicit_ls
integer, parameter, public broyden_type_2_explicit_ls
Definition: input_constants.F:737
input_constants::broyden_type_1_explicit
integer, parameter, public broyden_type_1_explicit
Definition: input_constants.F:737
input_constants::broyden_type_2_ls
integer, parameter, public broyden_type_2_ls
Definition: input_constants.F:737
input_constants::broyden_type_1
integer, parameter, public broyden_type_1
Definition: input_constants.F:737
input_constants::outer_scf_optimizer_broyden
integer, parameter, public outer_scf_optimizer_broyden
Definition: input_constants.F:727
input_constants::broyden_type_1_explicit_ls
integer, parameter, public broyden_type_1_explicit_ls
Definition: input_constants.F:737
input_constants::broyden_type_2_explicit
integer, parameter, public broyden_type_2_explicit
Definition: input_constants.F:737
input_constants::history_guess
integer, parameter, public history_guess
Definition: input_constants.F:201
input_constants::broyden_type_2
integer, parameter, public broyden_type_2
Definition: input_constants.F:737
input_constants::cdft2ot
integer, parameter, public cdft2ot
Definition: input_constants.F:762
input_constants::outer_scf_becke_constraint
integer, parameter, public outer_scf_becke_constraint
Definition: input_constants.F:718
input_constants::ot_precond_full_single
integer, parameter, public ot_precond_full_single
Definition: input_constants.F:513
input_constants::outer_scf_hirshfeld_constraint
integer, parameter, public outer_scf_hirshfeld_constraint
Definition: input_constants.F:718
input_constants::ot_precond_none
integer, parameter, public ot_precond_none
Definition: input_constants.F:513
input_constants::ot_precond_full_single_inverse
integer, parameter, public ot_precond_full_single_inverse
Definition: input_constants.F:513
input_constants::outer_scf_optimizer_newton_ls
integer, parameter, public outer_scf_optimizer_newton_ls
Definition: input_constants.F:727
input_constants::broyden_type_1_ls
integer, parameter, public broyden_type_1_ls
Definition: input_constants.F:737
input_constants::ot_precond_s_inverse
integer, parameter, public ot_precond_s_inverse
Definition: input_constants.F:513
input_constants::ot_precond_full_all
integer, parameter, public ot_precond_full_all
Definition: input_constants.F:513
input_constants::ot2cdft
integer, parameter, public ot2cdft
Definition: input_constants.F:762
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition: kinds.F:58
kpoint_io
Restart file for k point calculations.
Definition: kpoint_io.F:13
kpoint_io::write_kpoints_restart
subroutine, public write_kpoints_restart(denmat, kpoints, scf_env, dft_section, particle_set, qs_kind_set)
...
Definition: kpoint_io.F:75
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
machine
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition: machine.F:106
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition: machine.F:123
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
preconditioner
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition: preconditioner.F:15
preconditioner::restart_preconditioner
subroutine, public restart_preconditioner(qs_env, preconditioner, prec_type, nspins)
Allows for a restart of the preconditioner depending on the method it purges all arrays or keeps them...
Definition: preconditioner.F:244
preconditioner::prepare_preconditioner
subroutine, public prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, ot_preconditioner, prec_type, solver_type, energy_gap, nspins, has_unit_metric, convert_to_dbcsr, chol_type, full_mo_set)
...
Definition: preconditioner.F:309
pw_env_types
container for various plainwaves related things
Definition: pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
qs_block_davidson_types
module that contains the algorithms to perform an itrative diagonalization by the block-Davidson appr...
Definition: qs_block_davidson_types.F:18
qs_block_davidson_types::block_davidson_deallocate
subroutine, public block_davidson_deallocate(bdav_env)
...
Definition: qs_block_davidson_types.F:138
qs_cdft_scf_utils
Auxiliary routines for performing a constrained DFT SCF run with Quickstep.
Definition: qs_cdft_scf_utils.F:14
qs_cdft_scf_utils::prepare_jacobian_stencil
subroutine, public prepare_jacobian_stencil(qs_env, output_unit, nwork, pwork, coeff, step_multiplier, dh)
Prepares the finite difference stencil for computing the Jacobian. The constraints are re-evaluated b...
Definition: qs_cdft_scf_utils.F:65
qs_cdft_scf_utils::build_diagonal_jacobian
subroutine, public build_diagonal_jacobian(qs_env, used_history)
Builds a strictly diagonal inverse Jacobian from MD/SCF history.
Definition: qs_cdft_scf_utils.F:205
qs_cdft_scf_utils::print_inverse_jacobian
subroutine, public print_inverse_jacobian(logger, inv_jacobian, iter_count)
Prints the finite difference inverse Jacobian to file.
Definition: qs_cdft_scf_utils.F:326
qs_cdft_scf_utils::create_tmp_logger
subroutine, public create_tmp_logger(para_env, project_name, suffix, output_unit, tmp_logger)
Creates a temporary logger for redirecting output to a new file.
Definition: qs_cdft_scf_utils.F:364
qs_cdft_scf_utils::restart_inverse_jacobian
subroutine, public restart_inverse_jacobian(qs_env)
Restarts the finite difference inverse Jacobian.
Definition: qs_cdft_scf_utils.F:284
qs_cdft_scf_utils::initialize_inverse_jacobian
subroutine, public initialize_inverse_jacobian(scf_control, scf_env, explicit_jacobian, should_build, used_history)
Checks if the inverse Jacobian should be calculated and initializes the calculation.
Definition: qs_cdft_scf_utils.F:399
qs_cdft_types
Defines CDFT control structures.
Definition: qs_cdft_types.F:14
qs_charges_types
container for information about total charges on the grids
Definition: qs_charges_types.F:14
qs_density_matrices
collects routines that calculate density matrices
Definition: qs_density_matrices.F:14
qs_density_mixing_types
module that contains the definitions of the scf types
Definition: qs_density_mixing_types.F:14
qs_density_mixing_types::gspace_mixing_nr
integer, parameter, public gspace_mixing_nr
Definition: qs_density_mixing_types.F:45
qs_diis
Apply the direct inversion in the iterative subspace (DIIS) of Pulay in the framework of an SCF itera...
Definition: qs_diis.F:21
qs_diis::qs_diis_b_clear
pure subroutine, public qs_diis_b_clear(diis_buffer)
clears the buffer
Definition: qs_diis.F:517
qs_diis::qs_diis_b_create
subroutine, public qs_diis_b_create(diis_buffer, nbuffer)
Allocates an SCF DIIS buffer.
Definition: qs_diis.F:101
qs_diis::qs_diis_b_create_kp
subroutine, public qs_diis_b_create_kp(diis_buffer, nbuffer)
Allocates an SCF DIIS buffer for k-points.
Definition: qs_diis.F:865
qs_diis::qs_diis_b_clear_kp
pure subroutine, public qs_diis_b_clear_kp(diis_buffer)
clears the buffer
Definition: qs_diis.F:973
qs_energy_types
Definition: qs_energy_types.F:14
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Set the QUICKSTEP environment.
Definition: qs_environment_types.F:1035
qs_integrate_potential
Integrate single or product functions over a potential on a RS grid.
Definition: qs_integrate_potential.F:22
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_ks_methods
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition: qs_ks_methods.F:22
qs_ks_methods::qs_ks_update_qs_env
subroutine, public qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, print_active)
updates the Kohn Sham matrix of the given qs_env (facility method)
Definition: qs_ks_methods.F:1057
qs_ks_types
Definition: qs_ks_types.F:15
qs_ks_types::qs_ks_did_change
subroutine, public qs_ks_did_change(ks_env, s_mstruct_changed, rho_changed, potential_changed, full_reset)
tells that some of the things relevant to the ks calculation did change. has to be called when change...
Definition: qs_ks_types.F:919
qs_mo_io
Definition and initialisation of the mo data type.
Definition: qs_mo_io.F:21
qs_mo_io::write_mo_set_to_restart
subroutine, public write_mo_set_to_restart(mo_array, particle_set, dft_section, qs_kind_set)
...
Definition: qs_mo_io.F:107
qs_mo_methods
collects routines that perform operations directly related to MOs
Definition: qs_mo_methods.F:14
qs_mo_methods::make_basis_simple
subroutine, public make_basis_simple(vmatrix, ncol)
given a set of vectors, return an orthogonal (C^T C == 1) set spanning the same space (notice,...
Definition: qs_mo_methods.F:352
qs_mo_methods::make_basis_sm
subroutine, public make_basis_sm(vmatrix, ncol, matrix_s)
returns an S-orthonormal basis v (v^T S v ==1)
Definition: qs_mo_methods.F:87
qs_mo_occupation
Set occupation of molecular orbitals.
Definition: qs_mo_occupation.F:15
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::duplicate_mo_set
subroutine, public duplicate_mo_set(mo_set_new, mo_set_old)
allocate a new mo_set, and copy the old data
Definition: qs_mo_types.F:149
qs_mo_types::deallocate_mo_set
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
Definition: qs_mo_types.F:352
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
qs_mo_types::reassign_allocated_mos
subroutine, public reassign_allocated_mos(mo_set_new, mo_set_old)
reassign an already allocated mo_set
Definition: qs_mo_types.F:109
qs_ot_scf
basic functionality for using ot in the scf routines.
Definition: qs_ot_scf.F:14
qs_ot_scf::ot_scf_read_input
subroutine, public ot_scf_read_input(qs_ot_env, scf_section)
...
Definition: qs_ot_scf.F:68
qs_ot_scf::ot_scf_init
subroutine, public ot_scf_init(mo_array, matrix_s, qs_ot_env, matrix_ks, broyden_adaptive_sigma)
...
Definition: qs_ot_scf.F:352
qs_ot
orbital transformations
Definition: qs_ot.F:15
qs_ot::qs_ot_new_preconditioner
subroutine, public qs_ot_new_preconditioner(qs_ot_env, preconditioner)
...
Definition: qs_ot.F:70
qs_outer_scf
Routines for performing an outer scf loop.
Definition: qs_outer_scf.F:14
qs_outer_scf::outer_loop_purge_history
subroutine, public outer_loop_purge_history(qs_env)
purges outer_scf_history zeroing everything except the latest value of the outer_scf variable stored ...
Definition: qs_outer_scf.F:651
qs_outer_scf::outer_loop_switch
subroutine, public outer_loop_switch(scf_env, scf_control, cdft_control, dir)
switch between two outer_scf envs stored in cdft_control
Definition: qs_outer_scf.F:539
qs_outer_scf::outer_loop_optimize
subroutine, public outer_loop_optimize(scf_env, scf_control)
optimizes the parameters of the outer_scf
Definition: qs_outer_scf.F:187
qs_outer_scf::outer_loop_update_qs_env
subroutine, public outer_loop_update_qs_env(qs_env, scf_env)
propagates the updated variables to wherever they need to be set in qs_env
Definition: qs_outer_scf.F:387
qs_outer_scf::outer_loop_gradient
subroutine, public outer_loop_gradient(qs_env, scf_env)
computes the gradient wrt to the outer loop variables
Definition: qs_outer_scf.F:103
qs_rho_methods
methods of the rho structure (defined in qs_rho_types)
Definition: qs_rho_methods.F:15
qs_rho_methods::qs_rho_update_rho
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
Definition: qs_rho_methods.F:375
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
qs_scf_initialization
Utility routines for qs_scf.
Definition: qs_scf_initialization.F:11
qs_scf_initialization::qs_scf_env_initialize
subroutine, public qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section)
initializes input parameters if needed or restores values from previous runs to fill scf_env with the...
Definition: qs_scf_initialization.F:133
qs_scf_loop_utils
Utility routines for qs_scf.
Definition: qs_scf_loop_utils.F:10
qs_scf_loop_utils::qs_scf_new_mos
subroutine, public qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, energy_only)
takes known energy and derivatives and produces new wfns and or density matrix
Definition: qs_scf_loop_utils.F:123
qs_scf_loop_utils::qs_scf_inner_finalize
subroutine, public qs_scf_inner_finalize(scf_env, qs_env, diis_step, output_unit)
Performs the necessary steps before leaving innner scf loop.
Definition: qs_scf_loop_utils.F:650
qs_scf_loop_utils::qs_scf_check_inner_exit
subroutine, public qs_scf_check_inner_exit(qs_env, scf_env, scf_control, should_stop, exit_inner_loop, inner_loop_converged, output_unit)
checks whether exit conditions for inner loop are satisfied
Definition: qs_scf_loop_utils.F:533
qs_scf_loop_utils::qs_scf_set_loop_flags
subroutine, public qs_scf_set_loop_flags(scf_env, diis_step, energy_only, just_energy, exit_inner_loop)
computes properties for a given hamiltonian using the current wfn
Definition: qs_scf_loop_utils.F:92
qs_scf_loop_utils::qs_scf_rho_update
subroutine, public qs_scf_rho_update(rho, qs_env, scf_env, ks_env, mix_rho)
Performs the updates rho (takes care of mixing as well)
Definition: qs_scf_loop_utils.F:620
qs_scf_loop_utils::qs_scf_new_mos_kp
subroutine, public qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step)
Updates MOs and density matrix using diagonalization Kpoint code.
Definition: qs_scf_loop_utils.F:293
qs_scf_loop_utils::qs_scf_check_outer_exit
subroutine, public qs_scf_check_outer_exit(qs_env, scf_env, scf_control, should_stop, outer_loop_converged, exit_outer_loop)
checks whether exit conditions for outer loop are satisfied
Definition: qs_scf_loop_utils.F:495
qs_scf_loop_utils::qs_scf_density_mixing
subroutine, public qs_scf_density_mixing(scf_env, rho, para_env, diis_step)
Performs the requested density mixing if any needed.
Definition: qs_scf_loop_utils.F:454
qs_scf_output
Definition: qs_scf_output.F:8
qs_scf_output::qs_scf_loop_print
subroutine, public qs_scf_loop_print(qs_env, scf_env, para_env)
collects the 'heavy duty' printing tasks out of the SCF loop
Definition: qs_scf_output.F:809
qs_scf_output::qs_scf_outer_loop_info
subroutine, public qs_scf_outer_loop_info(output_unit, scf_control, scf_env, energy, total_steps, should_stop, outer_loop_converged)
writes basic information obtained in a scf outer loop step
Definition: qs_scf_output.F:491
qs_scf_output::qs_scf_write_mos
subroutine, public qs_scf_write_mos(qs_env, scf_env, final_mos)
Write the MO eigenvector, eigenvalues, and occupation numbers to the output unit.
Definition: qs_scf_output.F:178
qs_scf_output::qs_scf_loop_info
subroutine, public qs_scf_loop_info(scf_env, output_unit, just_energy, t1, t2, energy)
writes basic information obtained in a scf step
Definition: qs_scf_output.F:528
qs_scf_output::qs_scf_cdft_info
subroutine, public qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, energy, total_steps, should_stop, outer_loop_converged, cdft_loop)
writes CDFT constraint information and optionally CDFT scf loop info
Definition: qs_scf_output.F:950
qs_scf_output::qs_scf_cdft_initial_info
subroutine, public qs_scf_cdft_initial_info(output_unit, cdft_control)
writes information about the CDFT env
Definition: qs_scf_output.F:990
qs_scf_post_scf
Utility routines for qs_scf.
Definition: qs_scf_post_scf.F:11
qs_scf_post_scf::qs_scf_compute_properties
subroutine, public qs_scf_compute_properties(qs_env, wf_type, do_mp2)
computes properties for a given hamilonian using the current wfn
Definition: qs_scf_post_scf.F:45
qs_scf_types
module that contains the definitions of the scf types
Definition: qs_scf_types.F:14
qs_scf_types::ot_diag_method_nr
integer, parameter, public ot_diag_method_nr
Definition: qs_scf_types.F:51
qs_scf_types::filter_matrix_diag_method_nr
integer, parameter, public filter_matrix_diag_method_nr
Definition: qs_scf_types.F:51
qs_scf_types::block_davidson_diag_method_nr
integer, parameter, public block_davidson_diag_method_nr
Definition: qs_scf_types.F:51
qs_scf_types::ot_method_nr
integer, parameter, public ot_method_nr
Definition: qs_scf_types.F:51
qs_scf_types::special_diag_method_nr
integer, parameter, public special_diag_method_nr
Definition: qs_scf_types.F:51
qs_scf_types::block_krylov_diag_method_nr
integer, parameter, public block_krylov_diag_method_nr
Definition: qs_scf_types.F:51
qs_scf_types::general_diag_method_nr
integer, parameter, public general_diag_method_nr
Definition: qs_scf_types.F:51
qs_scf
Routines for the Quickstep SCF run.
Definition: qs_scf.F:47
qs_scf::scf
subroutine, public scf(qs_env, has_converged, total_scf_steps)
perform an scf procedure in the given qs_env
Definition: qs_scf.F:201
qs_scf::scf_env_cleanup
subroutine, public scf_env_cleanup(scf_env)
perform cleanup operations (like releasing temporary storage) at the end of the scf
Definition: qs_scf.F:895
qs_scf::scf_env_do_scf
subroutine, public scf_env_do_scf(scf_env, scf_control, qs_env, converged, should_stop, total_scf_steps)
perform an scf loop
Definition: qs_scf.F:337
qs_scf::cdft_scf
subroutine, public cdft_scf(qs_env, should_stop)
perform a CDFT scf procedure in the given qs_env
Definition: qs_scf.F:986
qs_wf_history_methods
Storage of past states of the qs_env. Methods to interpolate (or actually normally extrapolate) the n...
Definition: qs_wf_history_methods.F:21
qs_wf_history_methods::wfi_purge_history
subroutine, public wfi_purge_history(qs_env)
purges wf_history retaining only the latest snapshot
Definition: qs_wf_history_methods.F:1063
qs_wf_history_methods::wfi_update
subroutine, public wfi_update(wf_history, qs_env, dt)
updates the snapshot buffer, taking a new snapshot
Definition: qs_wf_history_methods.F:968
scf_control_types
parameters that control an scf iteration
Definition: scf_control_types.F:17