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