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