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