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