(git:9754b87)
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 scf_loop: DO
480
481 CALL timeset(routinen//"_inner_loop", handle2)
482
483 scf_env%iter_count = scf_env%iter_count + 1
484 iter_count = iter_count + 1
485 CALL cp_iterate(logger%iter_info, last=.false., iter_nr=iter_count)
486
487 IF (output_unit > 0) CALL m_flush(output_unit)
488
489 total_steps = total_steps + 1
490 just_energy = energy_only
491
492 CALL qs_ks_update_qs_env(qs_env, just_energy=just_energy, &
493 calculate_forces=.false.)
494
495 ! print 'heavy weight' or relatively expensive quantities
496 CALL qs_scf_loop_print(qs_env, scf_env, para_env)
497
498 IF (do_kpoints) THEN
499 ! kpoints
500 CALL qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step)
501 ELSE
502 ! Gamma points only
503 CALL qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, energy_only)
504 END IF
505
506 ! Print requested MO information (can be computationally expensive with OT)
507 CALL qs_scf_write_mos(qs_env, scf_env, final_mos=.false.)
508
509 CALL qs_scf_density_mixing(scf_env, rho, para_env, diis_step)
510
511 t2 = m_walltime()
512
513 CALL qs_scf_loop_info(scf_env, output_unit, just_energy, t1, t2, energy)
514
515 IF (.NOT. just_energy) energy%tot_old = energy%total
516
517 ! check for external communicator and if the intermediate energy should be sent
518 IF (scf_energy_message_tag .GT. 0) THEN
519 CALL external_comm%send(energy%total, ext_master_id, scf_energy_message_tag)
520 END IF
521
522 CALL qs_scf_check_inner_exit(qs_env, scf_env, scf_control, should_stop, exit_inner_loop, &
523 inner_loop_converged, output_unit)
524
525 ! In case we decide to exit we perform few more check to see if this one
526 ! is really the last SCF step
527 IF (exit_inner_loop) THEN
528
529 CALL qs_scf_inner_finalize(scf_env, qs_env, diis_step, output_unit)
530
531 CALL qs_scf_check_outer_exit(qs_env, scf_env, scf_control, should_stop, &
532 outer_loop_converged, exit_outer_loop)
533
534 ! Let's tag the last SCF cycle so we can print informations only of the last step
535 IF (exit_outer_loop) CALL cp_iterate(logger%iter_info, last=.true., iter_nr=iter_count)
536
537 END IF
538
539 IF (do_kpoints) THEN
540 CALL write_kpoints_restart(rho_ao_kp, kpoints, scf_env, dft_section, particle_set, qs_kind_set)
541 ELSE
542 ! Write Wavefunction restart file
543 CALL write_mo_set_to_restart(mos, particle_set, dft_section=dft_section, qs_kind_set=qs_kind_set)
544 END IF
545
546 ! Exit if we have finished with the SCF inner loop
547 IF (exit_inner_loop) THEN
548 CALL timestop(handle2)
549 EXIT scf_loop
550 END IF
551
552 IF (.NOT. btest(cp_print_key_should_output(logger%iter_info, &
553 scf_section, "PRINT%ITERATION_INFO/TIME_CUMUL"), cp_p_file)) &
554 t1 = m_walltime()
555
556 ! mixing methods have the new density matrix in p_mix_new
557 IF (scf_env%mixing_method > 0) THEN
558 DO ic = 1, SIZE(rho_ao_kp, 2)
559 DO ispin = 1, dft_control%nspins
560 CALL dbcsr_get_info(rho_ao_kp(ispin, ic)%matrix, name=name) ! keep the name
561 CALL dbcsr_copy(rho_ao_kp(ispin, ic)%matrix, scf_env%p_mix_new(ispin, ic)%matrix, name=name)
562 END DO
563 END DO
564 END IF
565
566 CALL qs_scf_rho_update(rho, qs_env, scf_env, ks_env, &
567 mix_rho=scf_env%mixing_method >= gspace_mixing_nr)
568
569 CALL timestop(handle2)
570
571 END DO scf_loop
572
573 IF (.NOT. scf_control%outer_scf%have_scf) EXIT scf_outer_loop
574
575 ! In case we use the OUTER SCF loop let's print some info..
576 CALL qs_scf_outer_loop_info(output_unit, scf_control, scf_env, &
577 energy, total_steps, should_stop, outer_loop_converged)
578
579! save mos to converged mos if outer_loop_converged and surf_dip_correct_switch is true
580 IF (exit_outer_loop) THEN
581 IF ((dft_control%switch_surf_dip) .AND. (outer_loop_converged) .AND. &
582 (dft_control%surf_dip_correct_switch)) THEN
583 DO ispin = 1, dft_control%nspins
584 CALL reassign_allocated_mos(mos_last_converged(ispin), mos(ispin))
585 END DO
586 IF (output_unit > 0) WRITE (unit=output_unit, fmt="(/,/,T2,A)") &
587 "COPIED mos ---> mos_last_converged"
588 END IF
589 END IF
590
591 IF (exit_outer_loop) EXIT scf_outer_loop
592
593 !
594 CALL outer_loop_optimize(scf_env, scf_control)
595 CALL outer_loop_update_qs_env(qs_env, scf_env)
596 CALL qs_ks_did_change(ks_env, potential_changed=.true.)
597
598 END DO scf_outer_loop
599
600 converged = inner_loop_converged .AND. outer_loop_converged
601 total_scf_steps = total_steps
602
603 IF (dft_control%qs_control%cdft) &
604 dft_control%qs_control%cdft_control%total_steps = &
605 dft_control%qs_control%cdft_control%total_steps + total_steps
606
607 IF (.NOT. converged) THEN
608 IF (scf_control%ignore_convergence_failure .OR. should_stop) THEN
609 CALL cp_warn(__location__, "SCF run NOT converged")
610 ELSE
611 CALL cp_abort(__location__, &
612 "SCF run NOT converged. To continue the calculation "// &
613 "regardless, please set the keyword IGNORE_CONVERGENCE_FAILURE.")
614 END IF
615 END IF
616
617 ! Skip Harris functional calculation if ground-state is NOT converged
618 IF (qs_env%energy_correction) THEN
619 CALL get_qs_env(qs_env, ec_env=ec_env)
620 ec_env%do_skip = .false.
621 IF (ec_env%skip_ec .AND. .NOT. converged) ec_env%do_skip = .true.
622 END IF
623
624 ! if needed copy mo_coeff dbcsr->fm for later use in post_scf!fm->dbcsr
625 DO ispin = 1, SIZE(mos) !fm -> dbcsr
626 IF (mos(ispin)%use_mo_coeff_b) THEN !fm->dbcsr
627 IF (.NOT. ASSOCIATED(mos(ispin)%mo_coeff_b)) & !fm->dbcsr
628 cpabort("mo_coeff_b is not allocated") !fm->dbcsr
629 CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, & !fm->dbcsr
630 mos(ispin)%mo_coeff) !fm -> dbcsr
631 END IF !fm->dbcsr
632 END DO !fm -> dbcsr
633
634 CALL cp_rm_iter_level(logger%iter_info, level_name="QS_SCF")
635 CALL timestop(handle)
636
637 END SUBROUTINE scf_env_do_scf
638
639! **************************************************************************************************
640!> \brief inits those objects needed if you want to restart the scf with, say
641!> only a new initial guess, or different density functional or ...
642!> this will happen just before the scf loop starts
643!> \param scf_env ...
644!> \param qs_env ...
645!> \param scf_section ...
646!> \par History
647!> 03.2006 created [Joost VandeVondele]
648! **************************************************************************************************
649 SUBROUTINE init_scf_loop(scf_env, qs_env, scf_section)
650
651 TYPE(qs_scf_env_type), POINTER :: scf_env
652 TYPE(qs_environment_type), POINTER :: qs_env
653 TYPE(section_vals_type), POINTER :: scf_section
654
655 CHARACTER(LEN=*), PARAMETER :: routinen = 'init_scf_loop'
656
657 INTEGER :: handle, ispin, nmo, number_of_ot_envs
658 LOGICAL :: do_kpoints, do_rotation, &
659 has_unit_metric, is_full_all
660 TYPE(cp_fm_type), POINTER :: mo_coeff
661 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
662 TYPE(dbcsr_type), POINTER :: orthogonality_metric
663 TYPE(dft_control_type), POINTER :: dft_control
664 TYPE(kpoint_type), POINTER :: kpoints
665 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
666 TYPE(scf_control_type), POINTER :: scf_control
667
668 CALL timeset(routinen, handle)
669
670 NULLIFY (scf_control, matrix_s, matrix_ks, dft_control, mos, mo_coeff, kpoints)
671
672 cpassert(ASSOCIATED(scf_env))
673 cpassert(ASSOCIATED(qs_env))
674
675 CALL get_qs_env(qs_env=qs_env, &
676 scf_control=scf_control, &
677 dft_control=dft_control, &
678 do_kpoints=do_kpoints, &
679 kpoints=kpoints, &
680 mos=mos)
681
682 ! if using mo_coeff_b then copy to fm
683 DO ispin = 1, SIZE(mos) !fm->dbcsr
684 IF (mos(1)%use_mo_coeff_b) THEN !fm->dbcsr
685 CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, mos(ispin)%mo_coeff) !fm->dbcsr
686 END IF !fm->dbcsr
687 END DO !fm->dbcsr
688
689 ! this just guarantees that all mo_occupations match the eigenvalues, if smear
690 DO ispin = 1, dft_control%nspins
691 ! do not reset mo_occupations if the maximum overlap method is in use
692 IF (.NOT. scf_control%diagonalization%mom) &
693 CALL set_mo_occupation(mo_set=mos(ispin), &
694 smear=scf_control%smear)
695 END DO
696
697 SELECT CASE (scf_env%method)
698 CASE DEFAULT
699
700 cpabort("unknown scf method method:"//cp_to_string(scf_env%method))
701
703
704 IF (.NOT. scf_env%skip_diis) THEN
705 IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
706 ALLOCATE (scf_env%scf_diis_buffer)
707 CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis)
708 END IF
709 CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
710 END IF
711
713 IF (.NOT. scf_env%skip_diis) THEN
714 IF (do_kpoints) THEN
715 IF (.NOT. ASSOCIATED(kpoints%scf_diis_buffer)) THEN
716 ALLOCATE (kpoints%scf_diis_buffer)
717 CALL qs_diis_b_create_kp(kpoints%scf_diis_buffer, nbuffer=scf_control%max_diis)
718 END IF
719 CALL qs_diis_b_clear_kp(kpoints%scf_diis_buffer)
720 ELSE
721 IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
722 ALLOCATE (scf_env%scf_diis_buffer)
723 CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis)
724 END IF
725 CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
726 END IF
727 END IF
728
729 CASE (ot_diag_method_nr)
730 CALL get_qs_env(qs_env, matrix_ks=matrix_ks, matrix_s=matrix_s)
731
732 IF (.NOT. scf_env%skip_diis) THEN
733 IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
734 ALLOCATE (scf_env%scf_diis_buffer)
735 CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis)
736 END IF
737 CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
738 END IF
739
740 ! disable DFTB and SE for now
741 IF (dft_control%qs_control%dftb .OR. &
742 dft_control%qs_control%xtb .OR. &
743 dft_control%qs_control%semi_empirical) THEN
744 cpabort("DFTB and SE not available with OT/DIAG")
745 END IF
746
747 ! if an old preconditioner is still around (i.e. outer SCF is active),
748 ! remove it if this could be worthwhile
749 CALL restart_preconditioner(qs_env, scf_env%ot_preconditioner, &
750 scf_control%diagonalization%ot_settings%preconditioner_type, &
751 dft_control%nspins)
752
753 CALL prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, scf_env%ot_preconditioner, &
754 scf_control%diagonalization%ot_settings%preconditioner_type, &
755 scf_control%diagonalization%ot_settings%precond_solver_type, &
756 scf_control%diagonalization%ot_settings%energy_gap, dft_control%nspins)
757
759 ! Preconditioner initialized within the loop, when required
760 CASE (ot_method_nr)
761 CALL get_qs_env(qs_env, &
762 has_unit_metric=has_unit_metric, &
763 matrix_s=matrix_s, &
764 matrix_ks=matrix_ks)
765
766 ! reortho the wavefunctions if we are having an outer scf and
767 ! this is not the first iteration
768 ! this is useful to avoid the build-up of numerical noise
769 ! however, we can not play this trick if restricted (don't mix non-equivalent orbs)
770 IF (scf_control%do_outer_scf_reortho) THEN
771 IF (scf_control%outer_scf%have_scf .AND. .NOT. dft_control%restricted) THEN
772 IF (scf_env%outer_scf%iter_count > 0) THEN
773 DO ispin = 1, dft_control%nspins
774 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
775 IF (has_unit_metric) THEN
776 CALL make_basis_simple(mo_coeff, nmo)
777 ELSE
778 CALL make_basis_sm(mo_coeff, nmo, matrix_s(1)%matrix)
779 END IF
780 END DO
781 END IF
782 END IF
783 ELSE
784 ! dont need any dirty trick for the numerically stable irac algorithm.
785 END IF
786
787 IF (.NOT. ASSOCIATED(scf_env%qs_ot_env)) THEN
788
789 ! restricted calculations require just one set of OT orbitals
790 number_of_ot_envs = dft_control%nspins
791 IF (dft_control%restricted) number_of_ot_envs = 1
792
793 ALLOCATE (scf_env%qs_ot_env(number_of_ot_envs))
794
795 ! XXX Joost XXX should disentangle reading input from this part
796 CALL ot_scf_read_input(scf_env%qs_ot_env, scf_section)
797
798 ! keep a note that we are restricted
799 IF (dft_control%restricted) THEN
800 scf_env%qs_ot_env(1)%restricted = .true.
801 ! requires rotation
802 IF (.NOT. scf_env%qs_ot_env(1)%settings%do_rotation) &
803 CALL cp_abort(__location__, &
804 "Restricted calculation with OT requires orbital rotation. Please "// &
805 "activate the OT%ROTATION keyword!")
806 ELSE
807 scf_env%qs_ot_env(:)%restricted = .false.
808 END IF
809
810 ! this will rotate the MOs to be eigen states, which is not compatible with rotation
811 ! e.g. mo_derivs here do not yet include potentially different occupations numbers
812 do_rotation = scf_env%qs_ot_env(1)%settings%do_rotation
813 ! only full all needs rotation
814 is_full_all = scf_env%qs_ot_env(1)%settings%preconditioner_type == ot_precond_full_all
815 IF (do_rotation .AND. is_full_all) &
816 cpabort('PRECONDITIONER FULL_ALL is not compatible with ROTATION.')
817
818 ! might need the KS matrix to init properly
819 CALL qs_ks_update_qs_env(qs_env, just_energy=.false., &
820 calculate_forces=.false.)
821
822 ! if an old preconditioner is still around (i.e. outer SCF is active),
823 ! remove it if this could be worthwhile
824 IF (.NOT. reuse_precond) &
825 CALL restart_preconditioner(qs_env, scf_env%ot_preconditioner, &
826 scf_env%qs_ot_env(1)%settings%preconditioner_type, &
827 dft_control%nspins)
828
829 !
830 ! preconditioning still needs to be done correctly with has_unit_metric
831 ! notice that a big part of the preconditioning (S^-1) is fine anyhow
832 !
833 IF (has_unit_metric) THEN
834 NULLIFY (orthogonality_metric)
835 ELSE
836 orthogonality_metric => matrix_s(1)%matrix
837 END IF
838
839 IF (.NOT. reuse_precond) &
840 CALL prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, scf_env%ot_preconditioner, &
841 scf_env%qs_ot_env(1)%settings%preconditioner_type, &
842 scf_env%qs_ot_env(1)%settings%precond_solver_type, &
843 scf_env%qs_ot_env(1)%settings%energy_gap, dft_control%nspins, &
844 has_unit_metric=has_unit_metric, &
845 chol_type=scf_env%qs_ot_env(1)%settings%cholesky_type)
846 IF (reuse_precond) reuse_precond = .false.
847
848 CALL ot_scf_init(mo_array=mos, matrix_s=orthogonality_metric, &
849 broyden_adaptive_sigma=qs_env%broyden_adaptive_sigma, &
850 qs_ot_env=scf_env%qs_ot_env, matrix_ks=matrix_ks(1)%matrix)
851
852 SELECT CASE (scf_env%qs_ot_env(1)%settings%preconditioner_type)
853 CASE (ot_precond_none)
855 DO ispin = 1, SIZE(scf_env%qs_ot_env)
856 CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin), &
857 scf_env%ot_preconditioner(ispin)%preconditioner)
858 END DO
860 DO ispin = 1, SIZE(scf_env%qs_ot_env)
861 CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin), &
862 scf_env%ot_preconditioner(1)%preconditioner)
863 END DO
864 CASE DEFAULT
865 DO ispin = 1, SIZE(scf_env%qs_ot_env)
866 CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin), &
867 scf_env%ot_preconditioner(1)%preconditioner)
868 END DO
869 END SELECT
870 END IF
871
872 ! if we have non-uniform occupations we should be using rotation
873 do_rotation = scf_env%qs_ot_env(1)%settings%do_rotation
874 DO ispin = 1, SIZE(mos)
875 IF (.NOT. mos(ispin)%uniform_occupation) THEN
876 cpassert(do_rotation)
877 END IF
878 END DO
879 END SELECT
880
881 ! another safety check
882 IF (dft_control%low_spin_roks) THEN
883 cpassert(scf_env%method == ot_method_nr)
884 do_rotation = scf_env%qs_ot_env(1)%settings%do_rotation
885 cpassert(do_rotation)
886 END IF
887
888 CALL timestop(handle)
889
890 END SUBROUTINE init_scf_loop
891
892! **************************************************************************************************
893!> \brief perform cleanup operations (like releasing temporary storage)
894!> at the end of the scf
895!> \param scf_env ...
896!> \par History
897!> 02.2003 created [fawzi]
898!> \author fawzi
899! **************************************************************************************************
900 SUBROUTINE scf_env_cleanup(scf_env)
901 TYPE(qs_scf_env_type), INTENT(INOUT) :: scf_env
902
903 CHARACTER(len=*), PARAMETER :: routinen = 'scf_env_cleanup'
904
905 INTEGER :: handle
906
907 CALL timeset(routinen, handle)
908
909 ! Release SCF work storage
910 CALL cp_fm_release(scf_env%scf_work1)
911
912 IF (ASSOCIATED(scf_env%scf_work1_red)) THEN
913 CALL cp_fm_release(scf_env%scf_work1_red)
914 END IF
915 IF (ASSOCIATED(scf_env%scf_work2)) THEN
916 CALL cp_fm_release(scf_env%scf_work2)
917 DEALLOCATE (scf_env%scf_work2)
918 NULLIFY (scf_env%scf_work2)
919 END IF
920 IF (ASSOCIATED(scf_env%scf_work2_red)) THEN
921 CALL cp_fm_release(scf_env%scf_work2_red)
922 DEALLOCATE (scf_env%scf_work2_red)
923 NULLIFY (scf_env%scf_work2_red)
924 END IF
925 IF (ASSOCIATED(scf_env%ortho)) THEN
926 CALL cp_fm_release(scf_env%ortho)
927 DEALLOCATE (scf_env%ortho)
928 NULLIFY (scf_env%ortho)
929 END IF
930 IF (ASSOCIATED(scf_env%ortho_red)) THEN
931 CALL cp_fm_release(scf_env%ortho_red)
932 DEALLOCATE (scf_env%ortho_red)
933 NULLIFY (scf_env%ortho_red)
934 END IF
935 IF (ASSOCIATED(scf_env%ortho_m1)) THEN
936 CALL cp_fm_release(scf_env%ortho_m1)
937 DEALLOCATE (scf_env%ortho_m1)
938 NULLIFY (scf_env%ortho_m1)
939 END IF
940 IF (ASSOCIATED(scf_env%ortho_m1_red)) THEN
941 CALL cp_fm_release(scf_env%ortho_m1_red)
942 DEALLOCATE (scf_env%ortho_m1_red)
943 NULLIFY (scf_env%ortho_m1_red)
944 END IF
945
946 IF (ASSOCIATED(scf_env%ortho_dbcsr)) THEN
947 CALL dbcsr_deallocate_matrix(scf_env%ortho_dbcsr)
948 END IF
949 IF (ASSOCIATED(scf_env%buf1_dbcsr)) THEN
950 CALL dbcsr_deallocate_matrix(scf_env%buf1_dbcsr)
951 END IF
952 IF (ASSOCIATED(scf_env%buf2_dbcsr)) THEN
953 CALL dbcsr_deallocate_matrix(scf_env%buf2_dbcsr)
954 END IF
955
956 IF (ASSOCIATED(scf_env%p_mix_new)) THEN
957 CALL dbcsr_deallocate_matrix_set(scf_env%p_mix_new)
958 END IF
959
960 IF (ASSOCIATED(scf_env%p_delta)) THEN
961 CALL dbcsr_deallocate_matrix_set(scf_env%p_delta)
962 END IF
963
964 ! Method dependent cleanup
965 SELECT CASE (scf_env%method)
966 CASE (ot_method_nr)
967 !
968 CASE (ot_diag_method_nr)
969 !
971 !
973 !
976 CALL block_davidson_deallocate(scf_env%block_davidson_env)
978 !
979 CASE (smeagol_method_nr)
980 !
981 CASE DEFAULT
982 cpabort("unknown scf method method:"//cp_to_string(scf_env%method))
983 END SELECT
984
985 IF (ASSOCIATED(scf_env%outer_scf%variables)) THEN
986 DEALLOCATE (scf_env%outer_scf%variables)
987 END IF
988 IF (ASSOCIATED(scf_env%outer_scf%count)) THEN
989 DEALLOCATE (scf_env%outer_scf%count)
990 END IF
991 IF (ASSOCIATED(scf_env%outer_scf%gradient)) THEN
992 DEALLOCATE (scf_env%outer_scf%gradient)
993 END IF
994 IF (ASSOCIATED(scf_env%outer_scf%energy)) THEN
995 DEALLOCATE (scf_env%outer_scf%energy)
996 END IF
997 IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian) .AND. &
998 scf_env%outer_scf%deallocate_jacobian) THEN
999 DEALLOCATE (scf_env%outer_scf%inv_jacobian)
1000 END IF
1001
1002 CALL timestop(handle)
1003
1004 END SUBROUTINE scf_env_cleanup
1005
1006! **************************************************************************************************
1007!> \brief perform a CDFT scf procedure in the given qs_env
1008!> \param qs_env the qs_environment where to perform the scf procedure
1009!> \param should_stop flag determining if calculation should stop
1010!> \par History
1011!> 12.2015 Created
1012!> \author Nico Holmberg
1013! **************************************************************************************************
1014 SUBROUTINE cdft_scf(qs_env, should_stop)
1015 TYPE(qs_environment_type), POINTER :: qs_env
1016 LOGICAL, INTENT(OUT) :: should_stop
1017
1018 CHARACTER(len=*), PARAMETER :: routinen = 'cdft_scf'
1019
1020 INTEGER :: handle, iatom, ispin, ivar, nmo, nvar, &
1021 output_unit, tsteps
1022 LOGICAL :: cdft_loop_converged, converged, &
1023 exit_cdft_loop, first_iteration, &
1024 my_uocc, uniform_occupation
1025 REAL(kind=dp), DIMENSION(:), POINTER :: mo_occupations
1026 TYPE(cdft_control_type), POINTER :: cdft_control
1027 TYPE(cp_logger_type), POINTER :: logger
1028 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
1029 TYPE(dft_control_type), POINTER :: dft_control
1030 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1031 TYPE(pw_env_type), POINTER :: pw_env
1032 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1033 TYPE(qs_energy_type), POINTER :: energy
1034 TYPE(qs_ks_env_type), POINTER :: ks_env
1035 TYPE(qs_rho_type), POINTER :: rho
1036 TYPE(qs_scf_env_type), POINTER :: scf_env
1037 TYPE(scf_control_type), POINTER :: scf_control
1038 TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
1039
1040 NULLIFY (scf_env, ks_env, energy, rho, matrix_s, rho_ao, cdft_control, logger, &
1041 dft_control, pw_env, auxbas_pw_pool, energy, ks_env, scf_env, dft_section, &
1042 input, scf_section, scf_control, mos, mo_occupations)
1043 logger => cp_get_default_logger()
1044
1045 cpassert(ASSOCIATED(qs_env))
1046 CALL get_qs_env(qs_env, scf_env=scf_env, energy=energy, &
1047 dft_control=dft_control, scf_control=scf_control, &
1048 ks_env=ks_env, input=input)
1049
1050 CALL timeset(routinen//"_loop", handle)
1051 dft_section => section_vals_get_subs_vals(input, "DFT")
1052 scf_section => section_vals_get_subs_vals(dft_section, "SCF")
1053 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", &
1054 extension=".scfLog")
1055 first_iteration = .true.
1056
1057 cdft_control => dft_control%qs_control%cdft_control
1058
1059 scf_env%outer_scf%iter_count = 0
1060 cdft_control%total_steps = 0
1061
1062 ! Write some info about the CDFT calculation
1063 IF (output_unit > 0) THEN
1064 WRITE (unit=output_unit, fmt="(/,/,T2,A)") &
1065 "CDFT EXTERNAL SCF WAVEFUNCTION OPTIMIZATION"
1066 CALL qs_scf_cdft_initial_info(output_unit, cdft_control)
1067 END IF
1068 IF (cdft_control%reuse_precond) THEN
1069 reuse_precond = .false.
1070 cdft_control%nreused = 0
1071 END IF
1072 cdft_outer_loop: DO
1073 ! Change outer_scf settings to OT settings
1074 CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot)
1075 ! Solve electronic structure with fixed value of constraint
1076 CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
1077 converged=converged, should_stop=should_stop, total_scf_steps=tsteps)
1078 ! Decide whether to reuse the preconditioner on the next iteration
1079 IF (cdft_control%reuse_precond) THEN
1080 ! For convergence in exactly one step, the preconditioner is always reused (assuming max_reuse > 0)
1081 ! usually this means that the electronic structure has already converged to the correct state
1082 ! but the constraint optimizer keeps jumping over the optimal solution
1083 IF (scf_env%outer_scf%iter_count == 1 .AND. scf_env%iter_count == 1 &
1084 .AND. cdft_control%total_steps /= 1) &
1085 cdft_control%nreused = cdft_control%nreused - 1
1086 ! SCF converged in less than precond_freq steps
1087 IF (scf_env%outer_scf%iter_count == 1 .AND. scf_env%iter_count .LE. cdft_control%precond_freq .AND. &
1088 cdft_control%total_steps /= 1 .AND. cdft_control%nreused .LT. cdft_control%max_reuse) THEN
1089 reuse_precond = .true.
1090 cdft_control%nreused = cdft_control%nreused + 1
1091 ELSE
1092 reuse_precond = .false.
1093 cdft_control%nreused = 0
1094 END IF
1095 END IF
1096 ! Update history purging counters
1097 IF (first_iteration .AND. cdft_control%purge_history) THEN
1098 cdft_control%istep = cdft_control%istep + 1
1099 IF (scf_env%outer_scf%iter_count .GT. 1) THEN
1100 cdft_control%nbad_conv = cdft_control%nbad_conv + 1
1101 IF (cdft_control%nbad_conv .GE. cdft_control%purge_freq .AND. &
1102 cdft_control%istep .GE. cdft_control%purge_offset) THEN
1103 cdft_control%nbad_conv = 0
1104 cdft_control%istep = 0
1105 cdft_control%should_purge = .true.
1106 END IF
1107 END IF
1108 END IF
1109 first_iteration = .false.
1110 ! Change outer_scf settings to CDFT settings
1111 CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft)
1112 CALL qs_scf_check_outer_exit(qs_env, scf_env, scf_control, should_stop, &
1113 cdft_loop_converged, exit_cdft_loop)
1114 CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
1115 energy, cdft_control%total_steps, &
1116 should_stop, cdft_loop_converged, cdft_loop=.true.)
1117 IF (exit_cdft_loop) EXIT cdft_outer_loop
1118 ! Check if the inverse Jacobian needs to be calculated
1119 CALL qs_calculate_inverse_jacobian(qs_env)
1120 ! Check if a line search should be performed to find an optimal step size for the optimizer
1121 CALL qs_cdft_line_search(qs_env)
1122 ! Optimize constraint
1123 CALL outer_loop_optimize(scf_env, scf_control)
1124 CALL outer_loop_update_qs_env(qs_env, scf_env)
1125 CALL qs_ks_did_change(ks_env, potential_changed=.true.)
1126 END DO cdft_outer_loop
1127
1128 cdft_control%ienergy = cdft_control%ienergy + 1
1129
1130 ! Store needed arrays for ET coupling calculation
1131 IF (cdft_control%do_et) THEN
1132 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, mos=mos)
1133 nvar = SIZE(cdft_control%target)
1134 ! Matrix representation of weight function
1135 ALLOCATE (cdft_control%wmat(nvar))
1136 DO ivar = 1, nvar
1137 CALL dbcsr_init_p(cdft_control%wmat(ivar)%matrix)
1138 CALL dbcsr_copy(cdft_control%wmat(ivar)%matrix, matrix_s(1)%matrix, &
1139 name="ET_RESTRAINT_MATRIX")
1140 CALL dbcsr_set(cdft_control%wmat(ivar)%matrix, 0.0_dp)
1141 CALL integrate_v_rspace(cdft_control%group(ivar)%weight, &
1142 hmat=cdft_control%wmat(ivar), qs_env=qs_env, &
1143 calculate_forces=.false., &
1144 gapw=dft_control%qs_control%gapw)
1145 END DO
1146 ! Overlap matrix
1147 CALL dbcsr_init_p(cdft_control%matrix_s%matrix)
1148 CALL dbcsr_copy(cdft_control%matrix_s%matrix, matrix_s(1)%matrix, &
1149 name="OVERLAP")
1150 ! Molecular orbital coefficients
1151 NULLIFY (cdft_control%mo_coeff)
1152 ALLOCATE (cdft_control%mo_coeff(dft_control%nspins))
1153 DO ispin = 1, dft_control%nspins
1154 CALL cp_fm_create(matrix=cdft_control%mo_coeff(ispin), &
1155 matrix_struct=qs_env%mos(ispin)%mo_coeff%matrix_struct, &
1156 name="MO_COEFF_A"//trim(adjustl(cp_to_string(ispin)))//"MATRIX")
1157 CALL cp_fm_to_fm(qs_env%mos(ispin)%mo_coeff, &
1158 cdft_control%mo_coeff(ispin))
1159 END DO
1160 ! Density matrix
1161 IF (cdft_control%calculate_metric) THEN
1162 CALL get_qs_env(qs_env, rho=rho)
1163 CALL qs_rho_get(rho, rho_ao=rho_ao)
1164 ALLOCATE (cdft_control%matrix_p(dft_control%nspins))
1165 DO ispin = 1, dft_control%nspins
1166 NULLIFY (cdft_control%matrix_p(ispin)%matrix)
1167 CALL dbcsr_init_p(cdft_control%matrix_p(ispin)%matrix)
1168 CALL dbcsr_copy(cdft_control%matrix_p(ispin)%matrix, rho_ao(ispin)%matrix, &
1169 name="DENSITY MATRIX")
1170 END DO
1171 END IF
1172 ! Copy occupation numbers if non-uniform occupation
1173 uniform_occupation = .true.
1174 DO ispin = 1, dft_control%nspins
1175 CALL get_mo_set(mo_set=mos(ispin), uniform_occupation=my_uocc)
1176 uniform_occupation = uniform_occupation .AND. my_uocc
1177 END DO
1178 IF (.NOT. uniform_occupation) THEN
1179 ALLOCATE (cdft_control%occupations(dft_control%nspins))
1180 DO ispin = 1, dft_control%nspins
1181 CALL get_mo_set(mo_set=mos(ispin), &
1182 nmo=nmo, &
1183 occupation_numbers=mo_occupations)
1184 ALLOCATE (cdft_control%occupations(ispin)%array(nmo))
1185 cdft_control%occupations(ispin)%array(1:nmo) = mo_occupations(1:nmo)
1186 END DO
1187 END IF
1188 END IF
1189
1190 ! Deallocate constraint storage if forces are not needed
1191 ! In case of a simulation with multiple force_evals,
1192 ! deallocate only if weight function should not be copied to different force_evals
1193 IF (.NOT. (cdft_control%save_pot .OR. cdft_control%transfer_pot)) THEN
1194 CALL get_qs_env(qs_env, pw_env=pw_env)
1195 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1196 DO iatom = 1, SIZE(cdft_control%group)
1197 CALL auxbas_pw_pool%give_back_pw(cdft_control%group(iatom)%weight)
1198 DEALLOCATE (cdft_control%group(iatom)%weight)
1199 END DO
1200 IF (cdft_control%atomic_charges) THEN
1201 DO iatom = 1, cdft_control%natoms
1202 CALL auxbas_pw_pool%give_back_pw(cdft_control%charge(iatom))
1203 END DO
1204 DEALLOCATE (cdft_control%charge)
1205 END IF
1206 IF (cdft_control%type == outer_scf_becke_constraint .AND. &
1207 cdft_control%becke_control%cavity_confine) THEN
1208 IF (.NOT. ASSOCIATED(cdft_control%becke_control%cavity_mat)) THEN
1209 CALL auxbas_pw_pool%give_back_pw(cdft_control%becke_control%cavity)
1210 ELSE
1211 DEALLOCATE (cdft_control%becke_control%cavity_mat)
1212 END IF
1213 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1214 IF (ASSOCIATED(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)) THEN
1215 CALL auxbas_pw_pool%give_back_pw(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)
1216 END IF
1217 END IF
1218 IF (ASSOCIATED(cdft_control%charges_fragment)) DEALLOCATE (cdft_control%charges_fragment)
1219 cdft_control%need_pot = .true.
1220 cdft_control%external_control = .false.
1221 END IF
1222
1223 CALL timestop(handle)
1224
1225 END SUBROUTINE cdft_scf
1226
1227! **************************************************************************************************
1228!> \brief perform cleanup operations for cdft_control
1229!> \param cdft_control container for the external CDFT SCF loop variables
1230!> \par History
1231!> 12.2015 created [Nico Holmberg]
1232!> \author Nico Holmberg
1233! **************************************************************************************************
1234 SUBROUTINE cdft_control_cleanup(cdft_control)
1235 TYPE(cdft_control_type), POINTER :: cdft_control
1236
1237 IF (ASSOCIATED(cdft_control%constraint%variables)) &
1238 DEALLOCATE (cdft_control%constraint%variables)
1239 IF (ASSOCIATED(cdft_control%constraint%count)) &
1240 DEALLOCATE (cdft_control%constraint%count)
1241 IF (ASSOCIATED(cdft_control%constraint%gradient)) &
1242 DEALLOCATE (cdft_control%constraint%gradient)
1243 IF (ASSOCIATED(cdft_control%constraint%energy)) &
1244 DEALLOCATE (cdft_control%constraint%energy)
1245 IF (ASSOCIATED(cdft_control%constraint%inv_jacobian) .AND. &
1246 cdft_control%constraint%deallocate_jacobian) &
1247 DEALLOCATE (cdft_control%constraint%inv_jacobian)
1248
1249 END SUBROUTINE cdft_control_cleanup
1250
1251! **************************************************************************************************
1252!> \brief Calculates the finite difference inverse Jacobian
1253!> \param qs_env the qs_environment_type where to compute the Jacobian
1254!> \par History
1255!> 01.2017 created [Nico Holmberg]
1256! **************************************************************************************************
1257 SUBROUTINE qs_calculate_inverse_jacobian(qs_env)
1258 TYPE(qs_environment_type), POINTER :: qs_env
1259
1260 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_calculate_inverse_jacobian'
1261
1262 CHARACTER(len=default_path_length) :: project_name
1263 INTEGER :: counter, handle, i, ispin, iter_count, &
1264 iwork, j, max_scf, nspins, nsteps, &
1265 nvar, nwork, output_unit, pwork, &
1266 tsteps, twork
1267 LOGICAL :: converged, explicit_jacobian, &
1268 should_build, should_stop, &
1269 use_md_history
1270 REAL(kind=dp) :: inv_error, step_size
1271 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coeff, dh, step_multiplier
1272 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: jacobian
1273 REAL(kind=dp), DIMENSION(:), POINTER :: energy
1274 REAL(kind=dp), DIMENSION(:, :), POINTER :: gradient, inv_jacobian
1275 TYPE(cdft_control_type), POINTER :: cdft_control
1276 TYPE(cp_logger_type), POINTER :: logger, tmp_logger
1277 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv
1278 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
1279 TYPE(dft_control_type), POINTER :: dft_control
1280 TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mos_stashed
1281 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1282 TYPE(mp_para_env_type), POINTER :: para_env
1283 TYPE(qs_energy_type), POINTER :: energy_qs
1284 TYPE(qs_ks_env_type), POINTER :: ks_env
1285 TYPE(qs_rho_type), POINTER :: rho
1286 TYPE(qs_scf_env_type), POINTER :: scf_env
1287 TYPE(scf_control_type), POINTER :: scf_control
1288
1289 NULLIFY (energy, gradient, p_rmpv, rho_ao_kp, mos, rho, &
1290 ks_env, scf_env, scf_control, dft_control, cdft_control, &
1291 inv_jacobian, para_env, tmp_logger, energy_qs)
1292 logger => cp_get_default_logger()
1293
1294 cpassert(ASSOCIATED(qs_env))
1295 CALL get_qs_env(qs_env, scf_env=scf_env, ks_env=ks_env, &
1296 scf_control=scf_control, mos=mos, rho=rho, &
1297 dft_control=dft_control, &
1298 para_env=para_env, energy=energy_qs)
1299 explicit_jacobian = .false.
1300 should_build = .false.
1301 use_md_history = .false.
1302 iter_count = scf_env%outer_scf%iter_count
1303 ! Quick exit if optimizer does not require Jacobian
1304 IF (.NOT. ASSOCIATED(scf_control%outer_scf%cdft_opt_control)) RETURN
1305 ! Check if Jacobian should be calculated and initialize
1306 CALL timeset(routinen, handle)
1307 CALL initialize_inverse_jacobian(scf_control, scf_env, explicit_jacobian, should_build, used_history)
1308 IF (scf_control%outer_scf%cdft_opt_control%jacobian_restart) THEN
1309 ! Restart from previously calculated inverse Jacobian
1310 should_build = .false.
1311 CALL restart_inverse_jacobian(qs_env)
1312 END IF
1313 IF (should_build) THEN
1314 scf_env%outer_scf%deallocate_jacobian = .false.
1315 ! Actually need to (re)build the Jacobian
1316 IF (explicit_jacobian) THEN
1317 ! Build Jacobian with finite differences
1318 cdft_control => dft_control%qs_control%cdft_control
1319 IF (.NOT. ASSOCIATED(cdft_control)) &
1320 CALL cp_abort(__location__, &
1321 "Optimizers that need the explicit Jacobian can"// &
1322 " only be used together with a valid CDFT constraint.")
1323 ! Redirect output from Jacobian calculation to a new file by creating a temporary logger
1324 project_name = logger%iter_info%project_name
1325 CALL create_tmp_logger(para_env, project_name, "-JacobianInfo.out", output_unit, tmp_logger)
1326 ! Save last converged state so we can roll back to it (mo_coeff and some outer_loop variables)
1327 nspins = dft_control%nspins
1328 ALLOCATE (mos_stashed(nspins))
1329 DO ispin = 1, nspins
1330 CALL duplicate_mo_set(mos_stashed(ispin), mos(ispin))
1331 END DO
1332 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1333 p_rmpv => rho_ao_kp(:, 1)
1334 ! Allocate work
1335 nvar = SIZE(scf_env%outer_scf%variables, 1)
1336 max_scf = scf_control%outer_scf%max_scf + 1
1337 ALLOCATE (gradient(nvar, max_scf))
1338 gradient = scf_env%outer_scf%gradient
1339 ALLOCATE (energy(max_scf))
1340 energy = scf_env%outer_scf%energy
1341 ALLOCATE (jacobian(nvar, nvar))
1342 jacobian = 0.0_dp
1343 nsteps = cdft_control%total_steps
1344 ! Setup finite difference scheme
1345 CALL prepare_jacobian_stencil(qs_env, output_unit, nwork, pwork, coeff, step_multiplier, dh)
1346 twork = pwork - nwork
1347 DO i = 1, nvar
1348 jacobian(i, :) = coeff(0)*scf_env%outer_scf%gradient(i, iter_count)
1349 END DO
1350 ! Calculate the Jacobian by perturbing each Lagrangian and recalculating the energy self-consistently
1351 CALL cp_add_default_logger(tmp_logger)
1352 DO i = 1, nvar
1353 IF (output_unit > 0) THEN
1354 WRITE (output_unit, fmt="(A)") " "
1355 WRITE (output_unit, fmt="(A)") " #####################################"
1356 WRITE (output_unit, '(A,I3,A,I3,A)') &
1357 " ### Constraint ", i, " of ", nvar, " ###"
1358 WRITE (output_unit, fmt="(A)") " #####################################"
1359 END IF
1360 counter = 0
1361 DO iwork = nwork, pwork
1362 IF (iwork == 0) cycle
1363 counter = counter + 1
1364 IF (output_unit > 0) THEN
1365 WRITE (output_unit, fmt="(A)") " #####################################"
1366 WRITE (output_unit, '(A,I3,A,I3,A)') &
1367 " ### Energy evaluation ", counter, " of ", twork, " ###"
1368 WRITE (output_unit, fmt="(A)") " #####################################"
1369 END IF
1370 IF (SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step) == 1) THEN
1371 step_size = scf_control%outer_scf%cdft_opt_control%jacobian_step(1)
1372 ELSE
1373 step_size = scf_control%outer_scf%cdft_opt_control%jacobian_step(i)
1374 END IF
1375 scf_env%outer_scf%variables(:, iter_count + 1) = scf_env%outer_scf%variables(:, iter_count)
1376 scf_env%outer_scf%variables(i, iter_count + 1) = scf_env%outer_scf%variables(i, iter_count) + &
1377 step_multiplier(iwork)*step_size
1378 CALL outer_loop_update_qs_env(qs_env, scf_env)
1379 CALL qs_ks_did_change(ks_env, potential_changed=.true.)
1380 CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot)
1381 CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
1382 converged=converged, should_stop=should_stop, total_scf_steps=tsteps)
1383 CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft)
1384 ! Update (iter_count + 1) element of gradient and print constraint info
1385 scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count + 1
1386 CALL outer_loop_gradient(qs_env, scf_env)
1387 CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
1388 energy_qs, cdft_control%total_steps, &
1389 should_stop=.false., outer_loop_converged=.false., cdft_loop=.false.)
1390 scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count - 1
1391 ! Update Jacobian
1392 DO j = 1, nvar
1393 jacobian(j, i) = jacobian(j, i) + coeff(iwork)*scf_env%outer_scf%gradient(j, iter_count + 1)
1394 END DO
1395 ! Reset everything to last converged state
1396 scf_env%outer_scf%variables(:, iter_count + 1) = 0.0_dp
1397 scf_env%outer_scf%gradient = gradient
1398 scf_env%outer_scf%energy = energy
1399 cdft_control%total_steps = nsteps
1400 DO ispin = 1, nspins
1401 CALL deallocate_mo_set(mos(ispin))
1402 CALL duplicate_mo_set(mos(ispin), mos_stashed(ispin))
1403 CALL calculate_density_matrix(mos(ispin), &
1404 p_rmpv(ispin)%matrix)
1405 END DO
1406 CALL qs_rho_update_rho(rho, qs_env=qs_env)
1407 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
1408 END DO
1409 END DO
1411 CALL cp_logger_release(tmp_logger)
1412 ! Finalize and invert Jacobian
1413 DO j = 1, nvar
1414 DO i = 1, nvar
1415 jacobian(i, j) = jacobian(i, j)/dh(j)
1416 END DO
1417 END DO
1418 IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
1419 ALLOCATE (scf_env%outer_scf%inv_jacobian(nvar, nvar))
1420 inv_jacobian => scf_env%outer_scf%inv_jacobian
1421 CALL invert_matrix(jacobian, inv_jacobian, inv_error)
1422 scf_control%outer_scf%cdft_opt_control%broyden_update = .false.
1423 ! Release temporary storage
1424 DO ispin = 1, nspins
1425 CALL deallocate_mo_set(mos_stashed(ispin))
1426 END DO
1427 DEALLOCATE (mos_stashed, jacobian, gradient, energy, coeff, step_multiplier, dh)
1428 IF (output_unit > 0) THEN
1429 WRITE (output_unit, fmt="(/,A)") &
1430 " ================================== JACOBIAN CALCULATED =================================="
1431 CALL close_file(unit_number=output_unit)
1432 END IF
1433 ELSE
1434 ! Build a strictly diagonal Jacobian from history and invert it
1435 CALL build_diagonal_jacobian(qs_env, used_history)
1436 END IF
1437 END IF
1438 IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian) .AND. para_env%is_source()) THEN
1439 ! Write restart file for inverse Jacobian
1440 CALL print_inverse_jacobian(logger, scf_env%outer_scf%inv_jacobian, iter_count)
1441 END IF
1442 ! Update counter
1443 scf_control%outer_scf%cdft_opt_control%ijacobian(1) = scf_control%outer_scf%cdft_opt_control%ijacobian(1) + 1
1444 CALL timestop(handle)
1445
1446 END SUBROUTINE qs_calculate_inverse_jacobian
1447
1448! **************************************************************************************************
1449!> \brief Perform backtracking line search to find the optimal step size for the CDFT constraint
1450!> optimizer. Assumes that the CDFT gradient function is a smooth function of the constraint
1451!> variables.
1452!> \param qs_env the qs_environment_type where to perform the line search
1453!> \par History
1454!> 02.2017 created [Nico Holmberg]
1455! **************************************************************************************************
1456 SUBROUTINE qs_cdft_line_search(qs_env)
1457 TYPE(qs_environment_type), POINTER :: qs_env
1458
1459 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_cdft_line_search'
1460
1461 CHARACTER(len=default_path_length) :: project_name
1462 INTEGER :: handle, i, ispin, iter_count, &
1463 max_linesearch, max_scf, nspins, &
1464 nsteps, nvar, output_unit, tsteps
1465 LOGICAL :: continue_ls, continue_ls_exit, converged, do_linesearch, found_solution, &
1466 reached_maxls, should_exit, should_stop, sign_changed
1467 LOGICAL, ALLOCATABLE, DIMENSION(:) :: positive_sign
1468 REAL(kind=dp) :: alpha, alpha_ls, factor, norm_ls
1469 REAL(kind=dp), DIMENSION(:), POINTER :: energy
1470 REAL(kind=dp), DIMENSION(:, :), POINTER :: gradient, inv_jacobian
1471 REAL(kind=dp), EXTERNAL :: dnrm2
1472 TYPE(cdft_control_type), POINTER :: cdft_control
1473 TYPE(cp_logger_type), POINTER :: logger, tmp_logger
1474 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv
1475 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
1476 TYPE(dft_control_type), POINTER :: dft_control
1477 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1478 TYPE(mp_para_env_type), POINTER :: para_env
1479 TYPE(qs_energy_type), POINTER :: energy_qs
1480 TYPE(qs_ks_env_type), POINTER :: ks_env
1481 TYPE(qs_rho_type), POINTER :: rho
1482 TYPE(qs_scf_env_type), POINTER :: scf_env
1483 TYPE(scf_control_type), POINTER :: scf_control
1484
1485 CALL timeset(routinen, handle)
1486
1487 NULLIFY (energy, gradient, p_rmpv, rho_ao_kp, mos, rho, &
1488 ks_env, scf_env, scf_control, dft_control, &
1489 cdft_control, inv_jacobian, para_env, &
1490 tmp_logger, energy_qs)
1491 logger => cp_get_default_logger()
1492
1493 cpassert(ASSOCIATED(qs_env))
1494 CALL get_qs_env(qs_env, scf_env=scf_env, ks_env=ks_env, &
1495 scf_control=scf_control, mos=mos, rho=rho, &
1496 dft_control=dft_control, &
1497 para_env=para_env, energy=energy_qs)
1498 do_linesearch = .false.
1499 SELECT CASE (scf_control%outer_scf%optimizer)
1500 CASE DEFAULT
1501 do_linesearch = .false.
1503 do_linesearch = .true.
1505 SELECT CASE (scf_control%outer_scf%cdft_opt_control%broyden_type)
1507 do_linesearch = .false.
1509 cdft_control => dft_control%qs_control%cdft_control
1510 IF (.NOT. ASSOCIATED(cdft_control)) &
1511 CALL cp_abort(__location__, &
1512 "Optimizers that perform a line search can"// &
1513 " only be used together with a valid CDFT constraint")
1514 IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
1515 do_linesearch = .true.
1516 END SELECT
1517 END SELECT
1518 IF (do_linesearch) THEN
1519 block
1520 TYPE(mo_set_type), DIMENSION(:), ALLOCATABLE :: mos_ls, mos_stashed
1521 cdft_control => dft_control%qs_control%cdft_control
1522 IF (.NOT. ASSOCIATED(cdft_control)) &
1523 CALL cp_abort(__location__, &
1524 "Optimizers that perform a line search can"// &
1525 " only be used together with a valid CDFT constraint")
1526 cpassert(ASSOCIATED(scf_env%outer_scf%inv_jacobian))
1527 cpassert(ASSOCIATED(scf_control%outer_scf%cdft_opt_control))
1528 alpha = scf_control%outer_scf%cdft_opt_control%newton_step_save
1529 iter_count = scf_env%outer_scf%iter_count
1530 ! Redirect output from line search procedure to a new file by creating a temporary logger
1531 project_name = logger%iter_info%project_name
1532 CALL create_tmp_logger(para_env, project_name, "-LineSearch.out", output_unit, tmp_logger)
1533 ! Save last converged state so we can roll back to it (mo_coeff and some outer_loop variables)
1534 nspins = dft_control%nspins
1535 ALLOCATE (mos_stashed(nspins))
1536 DO ispin = 1, nspins
1537 CALL duplicate_mo_set(mos_stashed(ispin), mos(ispin))
1538 END DO
1539 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1540 p_rmpv => rho_ao_kp(:, 1)
1541 nsteps = cdft_control%total_steps
1542 ! Allocate work
1543 nvar = SIZE(scf_env%outer_scf%variables, 1)
1544 max_scf = scf_control%outer_scf%max_scf + 1
1545 max_linesearch = scf_control%outer_scf%cdft_opt_control%max_ls
1546 continue_ls = scf_control%outer_scf%cdft_opt_control%continue_ls
1547 factor = scf_control%outer_scf%cdft_opt_control%factor_ls
1548 continue_ls_exit = .false.
1549 found_solution = .false.
1550 ALLOCATE (gradient(nvar, max_scf))
1551 gradient = scf_env%outer_scf%gradient
1552 ALLOCATE (energy(max_scf))
1553 energy = scf_env%outer_scf%energy
1554 reached_maxls = .false.
1555 ! Broyden optimizers: perform update of inv_jacobian if necessary
1556 IF (scf_control%outer_scf%cdft_opt_control%broyden_update) THEN
1557 CALL outer_loop_optimize(scf_env, scf_control)
1558 ! Reset the variables and prevent a reupdate of inv_jacobian
1559 scf_env%outer_scf%variables(:, iter_count + 1) = 0
1560 scf_control%outer_scf%cdft_opt_control%broyden_update = .false.
1561 END IF
1562 ! Print some info
1563 IF (output_unit > 0) THEN
1564 WRITE (output_unit, fmt="(/,A)") &
1565 " ================================== LINE SEARCH STARTED =================================="
1566 WRITE (output_unit, fmt="(A,I5,A)") &
1567 " Evaluating optimal step size for optimizer using a maximum of", max_linesearch, " steps"
1568 IF (continue_ls) THEN
1569 WRITE (output_unit, fmt="(A)") &
1570 " Line search continues until best step size is found or max steps are reached"
1571 END IF
1572 WRITE (output_unit, '(/,A,F5.3)') &
1573 " Initial step size: ", alpha
1574 WRITE (output_unit, '(/,A,F5.3)') &
1575 " Step size update factor: ", factor
1576 WRITE (output_unit, '(/,A,I10,A,I10)') &
1577 " Energy evaluation: ", cdft_control%ienergy, ", CDFT SCF iteration: ", iter_count
1578 END IF
1579 ! Perform backtracking line search
1580 CALL cp_add_default_logger(tmp_logger)
1581 DO i = 1, max_linesearch
1582 IF (output_unit > 0) THEN
1583 WRITE (output_unit, fmt="(A)") " "
1584 WRITE (output_unit, fmt="(A)") " #####################################"
1585 WRITE (output_unit, '(A,I10,A)') &
1586 " ### Line search step: ", i, " ###"
1587 WRITE (output_unit, fmt="(A)") " #####################################"
1588 END IF
1589 inv_jacobian => scf_env%outer_scf%inv_jacobian
1590 ! Newton update of CDFT variables with a step size of alpha
1591 scf_env%outer_scf%variables(:, iter_count + 1) = scf_env%outer_scf%variables(:, iter_count) - alpha* &
1592 matmul(inv_jacobian, scf_env%outer_scf%gradient(:, iter_count))
1593 ! With updated CDFT variables, perform SCF
1594 CALL outer_loop_update_qs_env(qs_env, scf_env)
1595 CALL qs_ks_did_change(ks_env, potential_changed=.true.)
1596 CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot)
1597 CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
1598 converged=converged, should_stop=should_stop, total_scf_steps=tsteps)
1599 CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft)
1600 ! Update (iter_count + 1) element of gradient and print constraint info
1601 scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count + 1
1602 CALL outer_loop_gradient(qs_env, scf_env)
1603 CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
1604 energy_qs, cdft_control%total_steps, &
1605 should_stop=.false., outer_loop_converged=.false., cdft_loop=.false.)
1606 scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count - 1
1607 ! Store sign of initial gradient for each variable for continue_ls
1608 IF (continue_ls .AND. .NOT. ALLOCATED(positive_sign)) THEN
1609 ALLOCATE (positive_sign(nvar))
1610 DO ispin = 1, nvar
1611 positive_sign(ispin) = scf_env%outer_scf%gradient(ispin, iter_count + 1) .GE. 0.0_dp
1612 END DO
1613 END IF
1614 ! Check if the L2 norm of the gradient decreased
1615 inv_jacobian => scf_env%outer_scf%inv_jacobian
1616 IF (dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count + 1), 1) .LT. &
1617 dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count), 1)) THEN
1618 ! Optimal step size found
1619 IF (.NOT. continue_ls) THEN
1620 should_exit = .true.
1621 ELSE
1622 ! But line search continues for at least one more iteration in an attempt to find a better solution
1623 ! if max number of steps is not exceeded
1624 IF (found_solution) THEN
1625 ! Check if the norm also decreased w.r.t. to previously found solution
1626 IF (dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count + 1), 1) .GT. norm_ls) THEN
1627 ! Norm increased => accept previous solution and exit
1628 continue_ls_exit = .true.
1629 END IF
1630 END IF
1631 ! Store current state and the value of alpha
1632 IF (.NOT. continue_ls_exit) THEN
1633 should_exit = .false.
1634 alpha_ls = alpha
1635 found_solution = .true.
1636 norm_ls = dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count + 1), 1)
1637 ! Check if the sign of the gradient has changed for all variables (w.r.t initial gradient)
1638 ! In this case we should exit because further line search steps will just increase the norm
1639 sign_changed = .true.
1640 DO ispin = 1, nvar
1641 sign_changed = sign_changed .AND. (positive_sign(ispin) .NEQV. &
1642 scf_env%outer_scf%gradient(ispin, iter_count + 1) .GE. 0.0_dp)
1643 END DO
1644 IF (.NOT. ALLOCATED(mos_ls)) THEN
1645 ALLOCATE (mos_ls(nspins))
1646 ELSE
1647 DO ispin = 1, nspins
1648 CALL deallocate_mo_set(mos_ls(ispin))
1649 END DO
1650 END IF
1651 DO ispin = 1, nspins
1652 CALL duplicate_mo_set(mos_ls(ispin), mos(ispin))
1653 END DO
1654 alpha = alpha*factor
1655 ! Exit on last iteration
1656 IF (i == max_linesearch) continue_ls_exit = .true.
1657 ! Exit if constraint target is satisfied to requested tolerance
1658 IF (sqrt(maxval(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count + 1)**2)) .LT. &
1659 scf_control%outer_scf%eps_scf) &
1660 continue_ls_exit = .true.
1661 ! Exit if line search jumped over the optimal step length
1662 IF (sign_changed) continue_ls_exit = .true.
1663 END IF
1664 END IF
1665 ELSE
1666 ! Gradient increased => alpha is too large (if the gradient function is smooth)
1667 should_exit = .false.
1668 ! Update alpha using Armijo's scheme
1669 alpha = alpha*factor
1670 END IF
1671 IF (continue_ls_exit) THEN
1672 ! Continuation of line search did not yield a better alpha, use previously located solution and exit
1673 alpha = alpha_ls
1674 DO ispin = 1, nspins
1675 CALL deallocate_mo_set(mos(ispin))
1676 CALL duplicate_mo_set(mos(ispin), mos_ls(ispin))
1677 CALL calculate_density_matrix(mos(ispin), &
1678 p_rmpv(ispin)%matrix)
1679 CALL deallocate_mo_set(mos_ls(ispin))
1680 END DO
1681 CALL qs_rho_update_rho(rho, qs_env=qs_env)
1682 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
1683 DEALLOCATE (mos_ls)
1684 should_exit = .true.
1685 END IF
1686 ! Reached max steps and SCF converged: continue with last iterated step size
1687 IF (.NOT. should_exit .AND. &
1688 (i == max_linesearch .AND. converged .AND. .NOT. found_solution)) THEN
1689 should_exit = .true.
1690 reached_maxls = .true.
1691 alpha = alpha*(1.0_dp/factor)
1692 END IF
1693 ! Reset outer SCF environment to last converged state
1694 scf_env%outer_scf%variables(:, iter_count + 1) = 0.0_dp
1695 scf_env%outer_scf%gradient = gradient
1696 scf_env%outer_scf%energy = energy
1697 ! Exit line search if a suitable step size was found
1698 IF (should_exit) EXIT
1699 ! Reset the electronic structure
1700 cdft_control%total_steps = nsteps
1701 DO ispin = 1, nspins
1702 CALL deallocate_mo_set(mos(ispin))
1703 CALL duplicate_mo_set(mos(ispin), mos_stashed(ispin))
1704 CALL calculate_density_matrix(mos(ispin), &
1705 p_rmpv(ispin)%matrix)
1706 END DO
1707 CALL qs_rho_update_rho(rho, qs_env=qs_env)
1708 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
1709 END DO
1710 scf_control%outer_scf%cdft_opt_control%newton_step = alpha
1711 IF (.NOT. should_exit) THEN
1712 CALL cp_warn(__location__, &
1713 "Line search did not converge. CDFT SCF proceeds with fixed step size.")
1714 scf_control%outer_scf%cdft_opt_control%newton_step = scf_control%outer_scf%cdft_opt_control%newton_step_save
1715 END IF
1716 IF (reached_maxls) &
1717 CALL cp_warn(__location__, &
1718 "Line search did not converge. CDFT SCF proceeds with lasted iterated step size.")
1720 CALL cp_logger_release(tmp_logger)
1721 ! Release temporary storage
1722 DO ispin = 1, nspins
1723 CALL deallocate_mo_set(mos_stashed(ispin))
1724 END DO
1725 DEALLOCATE (mos_stashed, gradient, energy)
1726 IF (ALLOCATED(positive_sign)) DEALLOCATE (positive_sign)
1727 IF (output_unit > 0) THEN
1728 WRITE (output_unit, fmt="(/,A)") &
1729 " ================================== LINE SEARCH COMPLETE =================================="
1730 CALL close_file(unit_number=output_unit)
1731 END IF
1732 END block
1733 END IF
1734
1735 CALL timestop(handle)
1736
1737 END SUBROUTINE qs_cdft_line_search
1738
1739END 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_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_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
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:650
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:901
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:1015
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.