1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines for the Quickstep SCF run. 8!> \par History 9!> - Joost VandeVondele (02.2002) 10!> added code for: incremental (pab and gvg) update 11!> initialisation (init_cube, l_info) 12!> - Joost VandeVondele (02.2002) 13!> called the poisson code of the classical part 14!> this takes into account the spherical cutoff and allows for 15!> isolated systems 16!> - Joost VandeVondele (02.2002) 17!> added multiple grid feature 18!> changed to spherical cutoff consistently (?) 19!> therefore removed the gradient correct functionals 20!> - updated with the new QS data structures (10.04.02,MK) 21!> - copy_matrix replaced by transfer_matrix (11.04.02,MK) 22!> - nrebuild_rho and nrebuild_gvg unified (12.04.02,MK) 23!> - set_mo_occupation for smearing of the MO occupation numbers 24!> (17.04.02,MK) 25!> - MO level shifting added (22.04.02,MK) 26!> - Usage of TYPE mo_set_p_type 27!> - Joost VandeVondele (05.2002) 28!> added cholesky based diagonalisation 29!> - 05.2002 added pao method [fawzi] 30!> - parallel FFT (JGH 22.05.2002) 31!> - 06.2002 moved KS matrix construction to qs_build_KS_matrix.F [fawzi] 32!> - started to include more LSD (01.2003,Joost VandeVondele) 33!> - 02.2003 scf_env [fawzi] 34!> - got rid of nrebuild (01.2004, Joost VandeVondele) 35!> - 10.2004 removed pao [fawzi] 36!> - 03.2006 large cleaning action [Joost VandeVondele] 37!> - High-spin ROKS added (05.04.06,MK) 38!> - Mandes (10.2013) 39!> intermediate energy communication with external communicator added 40!> - kpoints (08.2014, JGH) 41!> - unified k-point and gamma-point code (2014.11) [Ole Schuett] 42!> - added extra SCF loop for CDFT constraints (12.2015) [Nico Holmberg] 43!> \author Matthias Krack (30.04.2001) 44! ************************************************************************************************** 45MODULE qs_scf 46 USE atomic_kind_types, ONLY: atomic_kind_type 47 USE cp_control_types, ONLY: dft_control_type 48 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 49 dbcsr_deallocate_matrix_set 50 USE cp_files, ONLY: close_file 51 USE cp_fm_types, ONLY: cp_fm_create,& 52 cp_fm_release,& 53 cp_fm_to_fm,& 54 cp_fm_type 55 USE cp_log_handling, ONLY: cp_add_default_logger,& 56 cp_get_default_logger,& 57 cp_logger_release,& 58 cp_logger_type,& 59 cp_rm_default_logger,& 60 cp_to_string 61 USE cp_output_handling, ONLY: cp_add_iter_level,& 62 cp_iterate,& 63 cp_p_file,& 64 cp_print_key_should_output,& 65 cp_print_key_unit_nr,& 66 cp_rm_iter_level 67 USE cp_para_types, ONLY: cp_para_env_type 68 USE cp_result_methods, ONLY: get_results,& 69 test_for_result 70 USE cp_result_types, ONLY: cp_result_type 71 USE dbcsr_api, ONLY: dbcsr_copy,& 72 dbcsr_deallocate_matrix,& 73 dbcsr_get_info,& 74 dbcsr_init_p,& 75 dbcsr_p_type,& 76 dbcsr_set,& 77 dbcsr_type 78 USE input_constants, ONLY: & 79 broyden_type_1, broyden_type_1_explicit, broyden_type_1_explicit_ls, broyden_type_1_ls, & 80 broyden_type_2, broyden_type_2_explicit, broyden_type_2_explicit_ls, broyden_type_2_ls, & 81 cdft2ot, history_guess, ot2cdft, ot_precond_full_all, ot_precond_full_single, & 82 ot_precond_full_single_inverse, ot_precond_none, ot_precond_s_inverse, & 83 outer_scf_becke_constraint, outer_scf_hirshfeld_constraint, outer_scf_optimizer_broyden, & 84 outer_scf_optimizer_newton_ls 85 USE input_section_types, ONLY: section_vals_get_subs_vals,& 86 section_vals_type 87 USE kinds, ONLY: default_path_length,& 88 default_string_length,& 89 dp 90 USE kpoint_io, ONLY: write_kpoints_restart 91 USE kpoint_types, ONLY: kpoint_type 92 USE machine, ONLY: m_flush,& 93 m_walltime 94 USE mathlib, ONLY: invert_matrix 95 USE message_passing, ONLY: mp_send 96 USE particle_types, ONLY: particle_type 97 USE preconditioner, ONLY: prepare_preconditioner,& 98 restart_preconditioner 99 USE pw_env_types, ONLY: pw_env_get,& 100 pw_env_type 101 USE pw_pool_types, ONLY: pw_pool_give_back_pw,& 102 pw_pool_type 103 USE qs_block_davidson_types, ONLY: block_davidson_deallocate 104 USE qs_cdft_scf_utils, ONLY: build_diagonal_jacobian,& 105 create_tmp_logger,& 106 initialize_inverse_jacobian,& 107 prepare_jacobian_stencil,& 108 print_inverse_jacobian,& 109 restart_inverse_jacobian 110 USE qs_cdft_types, ONLY: cdft_control_type 111 USE qs_charges_types, ONLY: qs_charges_type 112 USE qs_density_mixing_types, ONLY: gspace_mixing_nr 113 USE qs_diis, ONLY: qs_diis_b_clear,& 114 qs_diis_b_create 115 USE qs_energy_types, ONLY: qs_energy_type 116 USE qs_environment_types, ONLY: get_qs_env,& 117 qs_environment_type,& 118 set_qs_env 119 USE qs_integrate_potential, ONLY: integrate_v_rspace 120 USE qs_kind_types, ONLY: qs_kind_type 121 USE qs_ks_methods, ONLY: qs_ks_update_qs_env 122 USE qs_ks_types, ONLY: qs_ks_did_change,& 123 qs_ks_env_type 124 USE qs_mo_io, ONLY: write_mo_set 125 USE qs_mo_methods, ONLY: calculate_density_matrix,& 126 make_basis_simple,& 127 make_basis_sm 128 USE qs_mo_occupation, ONLY: set_mo_occupation 129 USE qs_mo_types, ONLY: deallocate_mo_set,& 130 duplicate_mo_set,& 131 get_mo_set,& 132 mo_set_p_type 133 USE qs_ot, ONLY: qs_ot_new_preconditioner 134 USE qs_ot_scf, ONLY: ot_scf_init,& 135 ot_scf_read_input 136 USE qs_outer_scf, ONLY: outer_loop_gradient,& 137 outer_loop_optimize,& 138 outer_loop_purge_history,& 139 outer_loop_switch,& 140 outer_loop_update_qs_env 141 USE qs_rho_methods, ONLY: qs_rho_update_rho 142 USE qs_rho_types, ONLY: qs_rho_get,& 143 qs_rho_type 144 USE qs_scf_initialization, ONLY: qs_scf_env_initialize 145 USE qs_scf_loop_utils, ONLY: qs_scf_check_inner_exit,& 146 qs_scf_check_outer_exit,& 147 qs_scf_density_mixing,& 148 qs_scf_inner_finalize,& 149 qs_scf_new_mos,& 150 qs_scf_new_mos_kp,& 151 qs_scf_rho_update,& 152 qs_scf_set_loop_flags 153 USE qs_scf_output, ONLY: qs_scf_cdft_info,& 154 qs_scf_cdft_initial_info,& 155 qs_scf_loop_info,& 156 qs_scf_loop_print,& 157 qs_scf_outer_loop_info,& 158 qs_scf_write_mos 159 USE qs_scf_post_scf, ONLY: qs_scf_compute_properties 160 USE qs_scf_types, ONLY: & 161 block_davidson_diag_method_nr, block_krylov_diag_method_nr, filter_matrix_diag_method_nr, & 162 general_diag_method_nr, ot_diag_method_nr, ot_method_nr, qs_scf_env_type, scf_env_release, & 163 special_diag_method_nr 164 USE qs_wf_history_methods, ONLY: wfi_purge_history,& 165 wfi_update 166 USE scf_control_types, ONLY: scf_control_type 167#include "./base/base_uses.f90" 168 169 IMPLICIT NONE 170 171 PRIVATE 172 173 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf' 174 LOGICAL, PRIVATE :: reuse_precond = .FALSE. 175 LOGICAL, PRIVATE :: used_history = .FALSE. 176 177 PUBLIC :: scf, scf_env_cleanup, scf_env_do_scf, cdft_scf 178 179CONTAINS 180 181! ************************************************************************************************** 182!> \brief perform an scf procedure in the given qs_env 183!> \param qs_env the qs_environment where to perform the scf procedure 184!> \par History 185!> 02.2003 introduced scf_env, moved real work to scf_env_do_scf [fawzi] 186!> \author fawzi 187!> \note 188! ************************************************************************************************** 189 SUBROUTINE scf(qs_env) 190 TYPE(qs_environment_type), POINTER :: qs_env 191 192 CHARACTER(len=*), PARAMETER :: routineN = 'scf', routineP = moduleN//':'//routineN 193 194 INTEGER :: ihistory, max_scf_tmp 195 LOGICAL :: converged, outer_scf_loop, should_stop 196 LOGICAL, SAVE :: first_step_flag = .TRUE. 197 REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient_history, variable_history 198 TYPE(cp_logger_type), POINTER :: logger 199 TYPE(dft_control_type), POINTER :: dft_control 200 TYPE(qs_scf_env_type), POINTER :: scf_env 201 TYPE(scf_control_type), POINTER :: scf_control 202 TYPE(section_vals_type), POINTER :: dft_section, input, scf_section 203 204 NULLIFY (scf_env) 205 logger => cp_get_default_logger() 206 CPASSERT(ASSOCIATED(qs_env)) 207 CALL get_qs_env(qs_env, scf_env=scf_env, input=input, & 208 dft_control=dft_control, scf_control=scf_control) 209 IF (scf_control%max_scf > 0) THEN 210 211 dft_section => section_vals_get_subs_vals(input, "DFT") 212 scf_section => section_vals_get_subs_vals(dft_section, "SCF") 213 214 IF (.NOT. ASSOCIATED(scf_env)) THEN 215 CALL qs_scf_env_initialize(qs_env, scf_env) 216 ! Moved here from qs_scf_env_initialize to be able to have more scf_env 217 CALL set_qs_env(qs_env, scf_env=scf_env) 218 CALL scf_env_release(scf_env) 219 CALL get_qs_env(qs_env=qs_env, scf_env=scf_env) 220 ELSE 221 CALL qs_scf_env_initialize(qs_env, scf_env) 222 ENDIF 223 224 IF ((scf_control%density_guess .EQ. history_guess) .AND. (first_step_flag)) THEN 225 max_scf_tmp = scf_control%max_scf 226 scf_control%max_scf = 1 227 outer_scf_loop = scf_control%outer_scf%have_scf 228 scf_control%outer_scf%have_scf = .FALSE. 229 END IF 230 231 IF (.NOT. dft_control%qs_control%cdft) THEN 232 CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, & 233 converged=converged, should_stop=should_stop) 234 ELSE 235 ! Third SCF loop needed for CDFT with OT to properly restart OT inner loop 236 CALL cdft_scf(qs_env=qs_env, should_stop=should_stop) 237 END IF 238 239 ! If SCF has not converged, then we should not start MP2 240 IF (ASSOCIATED(qs_env%mp2_env)) qs_env%mp2_env%hf_fail = .NOT. converged 241 242 ! Add the converged outer_scf SCF gradient(s)/variable(s) to history 243 IF (scf_control%outer_scf%have_scf) THEN 244 ihistory = scf_env%outer_scf%iter_count 245 CALL get_qs_env(qs_env, gradient_history=gradient_history, & 246 variable_history=variable_history) 247 ! We only store the latest two values 248 gradient_history(:, 1) = gradient_history(:, 2) 249 gradient_history(:, 2) = scf_env%outer_scf%gradient(:, ihistory) 250 variable_history(:, 1) = variable_history(:, 2) 251 variable_history(:, 2) = scf_env%outer_scf%variables(:, ihistory) 252 ! Reset flag 253 IF (used_history) used_history = .FALSE. 254 ! Update a counter and check if the Jacobian should be deallocated 255 IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN 256 scf_control%outer_scf%cdft_opt_control%ijacobian(2) = scf_control%outer_scf%cdft_opt_control%ijacobian(2) + 1 257 IF (scf_control%outer_scf%cdft_opt_control%ijacobian(2) .GE. & 258 scf_control%outer_scf%cdft_opt_control%jacobian_freq(2) .AND. & 259 scf_control%outer_scf%cdft_opt_control%jacobian_freq(2) > 0) & 260 scf_env%outer_scf%deallocate_jacobian = .TRUE. 261 END IF 262 END IF 263 ! *** add the converged wavefunction to the wavefunction history 264 IF ((ASSOCIATED(qs_env%wf_history)) .AND. & 265 ((scf_control%density_guess .NE. history_guess) .OR. & 266 (.NOT. first_step_flag))) THEN 267 IF (.NOT. dft_control%qs_control%cdft) THEN 268 CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp) 269 ELSE 270 IF (dft_control%qs_control%cdft_control%should_purge) THEN 271 CALL wfi_purge_history(qs_env) 272 CALL outer_loop_purge_history(qs_env) 273 dft_control%qs_control%cdft_control%should_purge = .FALSE. 274 ELSE 275 CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp) 276 END IF 277 END IF 278 ELSE IF ((scf_control%density_guess .EQ. history_guess) .AND. & 279 (first_step_flag)) THEN 280 scf_control%max_scf = max_scf_tmp 281 scf_control%outer_scf%have_scf = outer_scf_loop 282 first_step_flag = .FALSE. 283 END IF 284 285 ! *** compute properties that depend on the converged wavefunction 286 IF (.NOT. (should_stop)) CALL qs_scf_compute_properties(qs_env) 287 288 ! *** cleanup 289 CALL scf_env_cleanup(scf_env) 290 IF (dft_control%qs_control%cdft) & 291 CALL cdft_control_cleanup(dft_control%qs_control%cdft_control) 292 293 END IF 294 295 END SUBROUTINE scf 296 297! ************************************************************************************************** 298!> \brief perform an scf loop 299!> \param scf_env the scf_env where to perform the scf procedure 300!> \param scf_control ... 301!> \param qs_env the qs_env, the scf_env lives in 302!> \param converged will be true / false if converged is reached 303!> \param should_stop ... 304!> \par History 305!> long history, see cvs and qs_scf module history 306!> 02.2003 introduced scf_env [fawzi] 307!> 09.2005 Frozen density approximation [TdK] 308!> 06.2007 Check for SCF iteration count early [jgh] 309!> \author Matthias Krack 310!> \note 311! ************************************************************************************************** 312 SUBROUTINE scf_env_do_scf(scf_env, scf_control, qs_env, converged, should_stop) 313 314 TYPE(qs_scf_env_type), POINTER :: scf_env 315 TYPE(scf_control_type), POINTER :: scf_control 316 TYPE(qs_environment_type), POINTER :: qs_env 317 LOGICAL, INTENT(OUT) :: converged, should_stop 318 319 CHARACTER(LEN=*), PARAMETER :: routineN = 'scf_env_do_scf', routineP = moduleN//':'//routineN 320 321 CHARACTER(LEN=default_string_length) :: description, name 322 INTEGER :: ext_master_id, external_comm, handle, handle2, i_tmp, ic, ispin, iter_count, & 323 output_unit, scf_energy_message_tag, total_steps 324 LOGICAL :: diis_step, do_kpoints, energy_only, exit_inner_loop, exit_outer_loop, & 325 inner_loop_converged, just_energy, outer_loop_converged 326 REAL(KIND=dp) :: t1, t2 327 REAL(KIND=dp), DIMENSION(3) :: res_val_3 328 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 329 TYPE(cp_logger_type), POINTER :: logger 330 TYPE(cp_para_env_type), POINTER :: para_env 331 TYPE(cp_result_type), POINTER :: results 332 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp 333 TYPE(dft_control_type), POINTER :: dft_control 334 TYPE(kpoint_type), POINTER :: kpoints 335 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 336 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 337 TYPE(pw_env_type), POINTER :: pw_env 338 TYPE(qs_charges_type), POINTER :: qs_charges 339 TYPE(qs_energy_type), POINTER :: energy 340 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 341 TYPE(qs_ks_env_type), POINTER :: ks_env 342 TYPE(qs_rho_type), POINTER :: rho 343 TYPE(section_vals_type), POINTER :: dft_section, input, scf_section 344 345 CALL timeset(routineN, handle) 346 347 NULLIFY (dft_control, rho, energy, & 348 logger, qs_charges, ks_env, mos, atomic_kind_set, qs_kind_set, & 349 particle_set, dft_section, input, & 350 scf_section, para_env, results, kpoints, pw_env, rho_ao_kp) 351 352 CPASSERT(ASSOCIATED(scf_env)) 353 CPASSERT(scf_env%ref_count > 0) 354 CPASSERT(ASSOCIATED(qs_env)) 355 CPASSERT(qs_env%ref_count > 0) 356 357 logger => cp_get_default_logger() 358 t1 = m_walltime() 359 360 CALL get_qs_env(qs_env=qs_env, & 361 energy=energy, & 362 particle_set=particle_set, & 363 qs_charges=qs_charges, & 364 ks_env=ks_env, & 365 atomic_kind_set=atomic_kind_set, & 366 qs_kind_set=qs_kind_set, & 367 rho=rho, & 368 mos=mos, & 369 input=input, & 370 dft_control=dft_control, & 371 do_kpoints=do_kpoints, & 372 kpoints=kpoints, & 373 results=results, & 374 pw_env=pw_env, & 375 para_env=para_env) 376 377 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp) 378 379 dft_section => section_vals_get_subs_vals(input, "DFT") 380 scf_section => section_vals_get_subs_vals(dft_section, "SCF") 381 382 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", & 383 extension=".scfLog") 384 385 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") & 386 "SCF WAVEFUNCTION OPTIMIZATION" 387 388 IF ((output_unit > 0) .AND. (.NOT. scf_control%use_ot)) THEN 389 WRITE (UNIT=output_unit, & 390 FMT="(/,T3,A,T12,A,T31,A,T39,A,T59,A,T75,A,/,T3,A)") & 391 "Step", "Update method", "Time", "Convergence", "Total energy", "Change", & 392 REPEAT("-", 78) 393 END IF 394 CALL cp_add_iter_level(logger%iter_info, "QS_SCF") 395 396 ! check for external communicator and if the intermediate energy should be sent 397 res_val_3(:) = -1.0_dp 398 description = "[EXT_SCF_ENER_COMM]" 399 IF (test_for_result(results, description=description)) THEN 400 CALL get_results(results, description=description, & 401 values=res_val_3, n_entries=i_tmp) 402 CPASSERT(i_tmp .EQ. 3) 403 IF (ALL(res_val_3(:) .LE. 0.0)) & 404 CALL cp_abort(__LOCATION__, & 405 " Trying to access result ("//TRIM(description)// & 406 ") which is not correctly stored.") 407 END IF 408 external_comm = NINT(res_val_3(1)) 409 ext_master_id = NINT(res_val_3(2)) 410 scf_energy_message_tag = NINT(res_val_3(3)) 411 412 ! *** outer loop of the scf, can treat other variables, 413 ! *** such as lagrangian multipliers 414 scf_env%outer_scf%iter_count = 0 415 iter_count = 0 416 total_steps = 0 417 energy%tot_old = 0.0_dp 418 419 scf_outer_loop: DO 420 421 CALL init_scf_loop(scf_env=scf_env, qs_env=qs_env, & 422 scf_section=scf_section) 423 424 CALL qs_scf_set_loop_flags(scf_env, diis_step, & 425 energy_only, just_energy, exit_inner_loop) 426 427 scf_loop: DO 428 429 CALL timeset(routineN//"_inner_loop", handle2) 430 431 scf_env%iter_count = scf_env%iter_count + 1 432 iter_count = iter_count + 1 433 CALL cp_iterate(logger%iter_info, last=.FALSE., iter_nr=iter_count) 434 435 IF (output_unit > 0) CALL m_flush(output_unit) 436 437 total_steps = total_steps + 1 438 just_energy = energy_only 439 440 CALL qs_ks_update_qs_env(qs_env, just_energy=just_energy, & 441 calculate_forces=.FALSE.) 442 443 ! print 'heavy weight' or relatively expensive quantities 444 CALL qs_scf_loop_print(qs_env, scf_env, para_env) 445 446 IF (do_kpoints) THEN 447 ! kpoints 448 CALL qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step) 449 ELSE 450 ! Gamma points only 451 CALL qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, energy_only) 452 END IF 453 454 ! another heavy weight print object, print controlled by dft_section 455 CALL qs_scf_write_mos(mos, atomic_kind_set, qs_kind_set, particle_set, dft_section) 456 457 CALL qs_scf_density_mixing(scf_env, rho, para_env, diis_step) 458 459 t2 = m_walltime() 460 461 CALL qs_scf_loop_info(scf_env, output_unit, just_energy, t1, t2, energy) 462 463 IF (.NOT. just_energy) energy%tot_old = energy%total 464 465 ! check for external communicator and if the intermediate energy should be sent 466 IF (scf_energy_message_tag .GT. 0) THEN 467 CALL mp_send(energy%total, ext_master_id, scf_energy_message_tag, external_comm) 468 END IF 469 470 CALL qs_scf_check_inner_exit(qs_env, scf_env, scf_control, should_stop, exit_inner_loop, & 471 inner_loop_converged, output_unit) 472 473 ! In case we decide to exit we perform few more check to see if this one 474 ! is really the last SCF step 475 IF (exit_inner_loop) THEN 476 477 CALL qs_scf_inner_finalize(scf_env, qs_env, diis_step, output_unit) 478 479 CALL qs_scf_check_outer_exit(qs_env, scf_env, scf_control, should_stop, & 480 outer_loop_converged, exit_outer_loop) 481 482 ! Let's tag the last SCF cycle so we can print informations only of the last step 483 IF (exit_outer_loop) CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=iter_count) 484 485 END IF 486 487 IF (do_kpoints) THEN 488 CALL write_kpoints_restart(rho_ao_kp, kpoints, scf_env, dft_section, particle_set, qs_kind_set) 489 ELSE 490 ! Write Wavefunction restart file 491 CALL write_mo_set(mos, particle_set, dft_section=dft_section, qs_kind_set=qs_kind_set) 492 END IF 493 494 ! Exit if we have finished with the SCF inner loop 495 IF (exit_inner_loop) THEN 496 CALL timestop(handle2) 497 EXIT scf_loop 498 END IF 499 500 IF (.NOT. BTEST(cp_print_key_should_output(logger%iter_info, & 501 scf_section, "PRINT%ITERATION_INFO/TIME_CUMUL"), cp_p_file)) & 502 t1 = m_walltime() 503 504 ! mixing methods have the new density matrix in p_mix_new 505 IF (scf_env%mixing_method > 0) THEN 506 DO ic = 1, SIZE(rho_ao_kp, 2) 507 DO ispin = 1, dft_control%nspins 508 CALL dbcsr_get_info(rho_ao_kp(ispin, ic)%matrix, name=name) ! keep the name 509 CALL dbcsr_copy(rho_ao_kp(ispin, ic)%matrix, scf_env%p_mix_new(ispin, ic)%matrix, name=name) 510 END DO 511 END DO 512 END IF 513 514 CALL qs_scf_rho_update(rho, qs_env, scf_env, ks_env, & 515 mix_rho=scf_env%mixing_method >= gspace_mixing_nr) 516 517 CALL timestop(handle2) 518 519 END DO scf_loop 520 521 IF (.NOT. scf_control%outer_scf%have_scf) EXIT scf_outer_loop 522 523 ! In case we use the OUTER SCF loop let's print some info.. 524 CALL qs_scf_outer_loop_info(output_unit, scf_control, scf_env, & 525 energy, total_steps, should_stop, outer_loop_converged) 526 527 IF (exit_outer_loop) EXIT scf_outer_loop 528 529 ! 530 CALL outer_loop_optimize(scf_env, scf_control) 531 CALL outer_loop_update_qs_env(qs_env, scf_env) 532 CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.) 533 534 END DO scf_outer_loop 535 536 converged = inner_loop_converged .AND. outer_loop_converged 537 538 IF (dft_control%qs_control%cdft) & 539 dft_control%qs_control%cdft_control%total_steps = & 540 dft_control%qs_control%cdft_control%total_steps + total_steps 541 542 IF (.NOT. converged) CPWARN("SCF run NOT converged") 543 544 IF (.NOT. converged .AND. scf_control%stop_higher_iter_level) THEN 545 logger%iter_info%last_iter(logger%iter_info%n_rlevel - 1) = .TRUE. 546 CPWARN("MD iteration also stops") 547 END IF 548 549 ! if needed copy mo_coeff dbcsr->fm for later use in post_scf!fm->dbcsr 550 DO ispin = 1, SIZE(mos) !fm -> dbcsr 551 IF (mos(ispin)%mo_set%use_mo_coeff_b) THEN !fm->dbcsr 552 IF (.NOT. ASSOCIATED(mos(ispin)%mo_set%mo_coeff_b)) & !fm->dbcsr 553 CPABORT("mo_coeff_b is not allocated") !fm->dbcsr 554 CALL copy_dbcsr_to_fm(mos(ispin)%mo_set%mo_coeff_b, & !fm->dbcsr 555 mos(ispin)%mo_set%mo_coeff) !fm -> dbcsr 556 ENDIF !fm->dbcsr 557 ENDDO !fm -> dbcsr 558 559 CALL cp_rm_iter_level(logger%iter_info, level_name="QS_SCF") 560 CALL timestop(handle) 561 562 END SUBROUTINE scf_env_do_scf 563 564! ************************************************************************************************** 565!> \brief inits those objects needed if you want to restart the scf with, say 566!> only a new initial guess, or different density functional or ... 567!> this will happen just before the scf loop starts 568!> \param scf_env ... 569!> \param qs_env ... 570!> \param scf_section ... 571!> \par History 572!> 03.2006 created [Joost VandeVondele] 573! ************************************************************************************************** 574 SUBROUTINE init_scf_loop(scf_env, qs_env, scf_section) 575 576 TYPE(qs_scf_env_type), POINTER :: scf_env 577 TYPE(qs_environment_type), POINTER :: qs_env 578 TYPE(section_vals_type), POINTER :: scf_section 579 580 CHARACTER(LEN=*), PARAMETER :: routineN = 'init_scf_loop', routineP = moduleN//':'//routineN 581 582 INTEGER :: handle, ispin, nmo, number_of_OT_envs 583 LOGICAL :: do_rotation, has_unit_metric, is_full_all 584 TYPE(cp_fm_type), POINTER :: mo_coeff 585 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s 586 TYPE(dbcsr_type), POINTER :: orthogonality_metric 587 TYPE(dft_control_type), POINTER :: dft_control 588 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 589 TYPE(scf_control_type), POINTER :: scf_control 590 591 CALL timeset(routineN, handle) 592 593 NULLIFY (scf_control, matrix_s, matrix_ks, dft_control, mos, mo_coeff) 594 595 CPASSERT(ASSOCIATED(scf_env)) 596 CPASSERT(scf_env%ref_count > 0) 597 CPASSERT(ASSOCIATED(qs_env)) 598 CPASSERT(qs_env%ref_count > 0) 599 600 CALL get_qs_env(qs_env=qs_env, & 601 scf_control=scf_control, & 602 dft_control=dft_control, & 603 mos=mos) 604 605 ! if using mo_coeff_b then copy to fm 606 DO ispin = 1, SIZE(mos) !fm->dbcsr 607 IF (mos(1)%mo_set%use_mo_coeff_b) THEN !fm->dbcsr 608 CALL copy_dbcsr_to_fm(mos(ispin)%mo_set%mo_coeff_b, mos(ispin)%mo_set%mo_coeff) !fm->dbcsr 609 ENDIF !fm->dbcsr 610 ENDDO !fm->dbcsr 611 612 ! this just guarantees that all mo_occupations match the eigenvalues, if smear 613 DO ispin = 1, dft_control%nspins 614 ! do not reset mo_occupations if the maximum overlap method is in use 615 IF (.NOT. scf_control%diagonalization%mom) & 616 CALL set_mo_occupation(mo_set=mos(ispin)%mo_set, & 617 smear=scf_control%smear) 618 ENDDO 619 620 SELECT CASE (scf_env%method) 621 CASE DEFAULT 622 623 CPABORT("unknown scf method method:"//cp_to_string(scf_env%method)) 624 625 CASE (filter_matrix_diag_method_nr) 626 627 IF (.NOT. scf_env%skip_diis) THEN 628 IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN 629 CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis) 630 END IF 631 CALL qs_diis_b_clear(scf_env%scf_diis_buffer) 632 END IF 633 634 CASE (general_diag_method_nr, special_diag_method_nr, block_krylov_diag_method_nr) 635 IF (.NOT. scf_env%skip_diis) THEN 636 IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN 637 CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis) 638 END IF 639 CALL qs_diis_b_clear(scf_env%scf_diis_buffer) 640 END IF 641 642 CASE (ot_diag_method_nr) 643 CALL get_qs_env(qs_env, matrix_ks=matrix_ks, matrix_s=matrix_s) 644 645 IF (.NOT. scf_env%skip_diis) THEN 646 IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN 647 CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis) 648 END IF 649 CALL qs_diis_b_clear(scf_env%scf_diis_buffer) 650 END IF 651 652 ! disable DFTB and SE for now 653 IF (dft_control%qs_control%dftb .OR. & 654 dft_control%qs_control%xtb .OR. & 655 dft_control%qs_control%semi_empirical) THEN 656 CPABORT("DFTB and SE not available with OT/DIAG") 657 END IF 658 659 ! if an old preconditioner is still around (i.e. outer SCF is active), 660 ! remove it if this could be worthwhile 661 CALL restart_preconditioner(qs_env, scf_env%ot_preconditioner, & 662 scf_control%diagonalization%ot_settings%preconditioner_type, & 663 dft_control%nspins) 664 665 CALL prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, scf_env%ot_preconditioner, & 666 scf_control%diagonalization%ot_settings%preconditioner_type, & 667 scf_control%diagonalization%ot_settings%precond_solver_type, & 668 scf_control%diagonalization%ot_settings%energy_gap, dft_control%nspins) 669 670 CASE (block_davidson_diag_method_nr) 671 ! Preconditioner initialized within the loop, when required 672 CASE (ot_method_nr) 673 CALL get_qs_env(qs_env, & 674 has_unit_metric=has_unit_metric, & 675 matrix_s=matrix_s, & 676 matrix_ks=matrix_ks) 677 678 ! reortho the wavefunctions if we are having an outer scf and 679 ! this is not the first iteration 680 ! this is useful to avoid the build-up of numerical noise 681 ! however, we can not play this trick if restricted (don't mix non-equivalent orbs) 682 IF (scf_control%do_outer_scf_reortho) THEN 683 IF (scf_control%outer_scf%have_scf .AND. .NOT. dft_control%restricted) THEN 684 IF (scf_env%outer_scf%iter_count > 0) THEN 685 DO ispin = 1, dft_control%nspins 686 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, nmo=nmo) 687 IF (has_unit_metric) THEN 688 CALL make_basis_simple(mo_coeff, nmo) 689 ELSE 690 CALL make_basis_sm(mo_coeff, nmo, matrix_s(1)%matrix) 691 ENDIF 692 ENDDO 693 ENDIF 694 ENDIF 695 ELSE 696 ! dont need any dirty trick for the numerically stable irac algorithm. 697 ENDIF 698 699 IF (.NOT. ASSOCIATED(scf_env%qs_ot_env)) THEN 700 701 ! restricted calculations require just one set of OT orbitals 702 number_of_OT_envs = dft_control%nspins 703 IF (dft_control%restricted) number_of_OT_envs = 1 704 705 ALLOCATE (scf_env%qs_ot_env(number_of_OT_envs)) 706 707 ! XXX Joost XXX should disentangle reading input from this part 708 CALL ot_scf_read_input(scf_env%qs_ot_env, scf_section) 709 710 ! keep a note that we are restricted 711 IF (dft_control%restricted) THEN 712 scf_env%qs_ot_env(1)%restricted = .TRUE. 713 ! requires rotation 714 IF (.NOT. scf_env%qs_ot_env(1)%settings%do_rotation) & 715 CALL cp_abort(__LOCATION__, & 716 "Restricted calculation with OT requires orbital rotation. Please "// & 717 "activate the OT%ROTATION keyword!") 718 ELSE 719 scf_env%qs_ot_env(:)%restricted = .FALSE. 720 ENDIF 721 722 ! might need the KS matrix to init properly 723 CALL qs_ks_update_qs_env(qs_env, just_energy=.FALSE., & 724 calculate_forces=.FALSE.) 725 726 ! if an old preconditioner is still around (i.e. outer SCF is active), 727 ! remove it if this could be worthwhile 728 IF (.NOT. reuse_precond) & 729 CALL restart_preconditioner(qs_env, scf_env%ot_preconditioner, & 730 scf_env%qs_ot_env(1)%settings%preconditioner_type, & 731 dft_control%nspins) 732 733 ! 734 ! preconditioning still needs to be done correctly with has_unit_metric 735 ! notice that a big part of the preconditioning (S^-1) is fine anyhow 736 ! 737 IF (has_unit_metric) THEN 738 NULLIFY (orthogonality_metric) 739 ELSE 740 orthogonality_metric => matrix_s(1)%matrix 741 ENDIF 742 743 IF (.NOT. reuse_precond) & 744 CALL prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, scf_env%ot_preconditioner, & 745 scf_env%qs_ot_env(1)%settings%preconditioner_type, & 746 scf_env%qs_ot_env(1)%settings%precond_solver_type, & 747 scf_env%qs_ot_env(1)%settings%energy_gap, dft_control%nspins, & 748 has_unit_metric=has_unit_metric, & 749 chol_type=scf_env%qs_ot_env(1)%settings%cholesky_type) 750 IF (reuse_precond) reuse_precond = .FALSE. 751 752 CALL ot_scf_init(mo_array=mos, matrix_s=orthogonality_metric, & 753 broyden_adaptive_sigma=qs_env%broyden_adaptive_sigma, & 754 qs_ot_env=scf_env%qs_ot_env, matrix_ks=matrix_ks(1)%matrix) 755 756 SELECT CASE (scf_env%qs_ot_env(1)%settings%preconditioner_type) 757 CASE (ot_precond_none) 758 CASE (ot_precond_full_all, ot_precond_full_single_inverse) 759 ! this will rotate the MOs to be eigen states, which is not compatible with rotation 760 ! e.g. mo_derivs here do not yet include potentially different occupations numbers 761 do_rotation = scf_env%qs_ot_env(1)%settings%do_rotation 762 ! only full all needs rotation 763 is_full_all = scf_env%qs_ot_env(1)%settings%preconditioner_type == ot_precond_full_all 764 CPASSERT(.NOT. (do_rotation .AND. is_full_all)) 765 DO ispin = 1, SIZE(scf_env%qs_ot_env) 766 CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin), & 767 scf_env%ot_preconditioner(ispin)%preconditioner) 768 ENDDO 769 CASE (ot_precond_s_inverse, ot_precond_full_single) 770 DO ispin = 1, SIZE(scf_env%qs_ot_env) 771 CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin), & 772 scf_env%ot_preconditioner(1)%preconditioner) 773 ENDDO 774 CASE DEFAULT 775 DO ispin = 1, SIZE(scf_env%qs_ot_env) 776 CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin), & 777 scf_env%ot_preconditioner(1)%preconditioner) 778 ENDDO 779 END SELECT 780 ENDIF 781 782 ! if we have non-uniform occupations we should be using rotation 783 do_rotation = scf_env%qs_ot_env(1)%settings%do_rotation 784 DO ispin = 1, SIZE(mos) 785 IF (.NOT. mos(ispin)%mo_set%uniform_occupation) THEN 786 CPASSERT(do_rotation) 787 ENDIF 788 ENDDO 789 END SELECT 790 791 ! another safety check 792 IF (dft_control%low_spin_roks) THEN 793 CPASSERT(scf_env%method == ot_method_nr) 794 do_rotation = scf_env%qs_ot_env(1)%settings%do_rotation 795 CPASSERT(do_rotation) 796 ENDIF 797 798 CALL timestop(handle) 799 800 END SUBROUTINE init_scf_loop 801 802! ************************************************************************************************** 803!> \brief perform cleanup operations (like releasing temporary storage) 804!> at the end of the scf 805!> \param scf_env ... 806!> \par History 807!> 02.2003 created [fawzi] 808!> \author fawzi 809! ************************************************************************************************** 810 SUBROUTINE scf_env_cleanup(scf_env) 811 TYPE(qs_scf_env_type), POINTER :: scf_env 812 813 CHARACTER(len=*), PARAMETER :: routineN = 'scf_env_cleanup', & 814 routineP = moduleN//':'//routineN 815 816 INTEGER :: handle, ispin 817 818 CALL timeset(routineN, handle) 819 820 CPASSERT(ASSOCIATED(scf_env)) 821 CPASSERT(scf_env%ref_count > 0) 822 823! *** Release SCF work storage *** 824 825 IF (ASSOCIATED(scf_env%scf_work1)) THEN 826 DO ispin = 1, SIZE(scf_env%scf_work1) 827 CALL cp_fm_release(scf_env%scf_work1(ispin)%matrix) 828 ENDDO 829 DEALLOCATE (scf_env%scf_work1) 830 ENDIF 831 IF (ASSOCIATED(scf_env%scf_work2)) CALL cp_fm_release(scf_env%scf_work2) 832 IF (ASSOCIATED(scf_env%ortho)) CALL cp_fm_release(scf_env%ortho) 833 IF (ASSOCIATED(scf_env%ortho_m1)) CALL cp_fm_release(scf_env%ortho_m1) 834 835 IF (ASSOCIATED(scf_env%ortho_dbcsr)) THEN 836 CALL dbcsr_deallocate_matrix(scf_env%ortho_dbcsr) 837 END IF 838 IF (ASSOCIATED(scf_env%buf1_dbcsr)) THEN 839 CALL dbcsr_deallocate_matrix(scf_env%buf1_dbcsr) 840 END IF 841 IF (ASSOCIATED(scf_env%buf2_dbcsr)) THEN 842 CALL dbcsr_deallocate_matrix(scf_env%buf2_dbcsr) 843 END IF 844 845 IF (ASSOCIATED(scf_env%p_mix_new)) THEN 846 CALL dbcsr_deallocate_matrix_set(scf_env%p_mix_new) 847 END IF 848 849 IF (ASSOCIATED(scf_env%p_delta)) THEN 850 CALL dbcsr_deallocate_matrix_set(scf_env%p_delta) 851 END IF 852 853! *** method dependent cleanup 854 SELECT CASE (scf_env%method) 855 CASE (ot_method_nr) 856 ! 857 CASE (ot_diag_method_nr) 858 ! 859 CASE (general_diag_method_nr) 860 ! 861 CASE (special_diag_method_nr) 862 ! 863 CASE (block_krylov_diag_method_nr) 864 CASE (block_davidson_diag_method_nr) 865 CALL block_davidson_deallocate(scf_env%block_davidson_env) 866 CASE (filter_matrix_diag_method_nr) 867 ! 868 CASE DEFAULT 869 CPABORT("unknown scf method method:"//cp_to_string(scf_env%method)) 870 END SELECT 871 872 IF (ASSOCIATED(scf_env%outer_scf%variables)) THEN 873 DEALLOCATE (scf_env%outer_scf%variables) 874 ENDIF 875 IF (ASSOCIATED(scf_env%outer_scf%count)) THEN 876 DEALLOCATE (scf_env%outer_scf%count) 877 ENDIF 878 IF (ASSOCIATED(scf_env%outer_scf%gradient)) THEN 879 DEALLOCATE (scf_env%outer_scf%gradient) 880 ENDIF 881 IF (ASSOCIATED(scf_env%outer_scf%energy)) THEN 882 DEALLOCATE (scf_env%outer_scf%energy) 883 ENDIF 884 IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian) .AND. & 885 scf_env%outer_scf%deallocate_jacobian) THEN 886 DEALLOCATE (scf_env%outer_scf%inv_jacobian) 887 ENDIF 888 889 CALL timestop(handle) 890 891 END SUBROUTINE scf_env_cleanup 892 893! ************************************************************************************************** 894!> \brief perform a CDFT scf procedure in the given qs_env 895!> \param qs_env the qs_environment where to perform the scf procedure 896!> \param should_stop flag determining if calculation should stop 897!> \par History 898!> 12.2015 Created 899!> \author Nico Holmberg 900! ************************************************************************************************** 901 SUBROUTINE cdft_scf(qs_env, should_stop) 902 TYPE(qs_environment_type), POINTER :: qs_env 903 LOGICAL, INTENT(OUT) :: should_stop 904 905 CHARACTER(len=*), PARAMETER :: routineN = 'cdft_scf', routineP = moduleN//':'//routineN 906 907 INTEGER :: handle, iatom, ispin, ivar, nmo, nvar, & 908 output_unit 909 LOGICAL :: cdft_loop_converged, converged, & 910 exit_cdft_loop, first_iteration, & 911 my_uocc, uniform_occupation 912 REAL(KIND=dp), DIMENSION(:), POINTER :: mo_occupations 913 TYPE(cdft_control_type), POINTER :: cdft_control 914 TYPE(cp_logger_type), POINTER :: logger 915 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao 916 TYPE(dft_control_type), POINTER :: dft_control 917 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 918 TYPE(pw_env_type), POINTER :: pw_env 919 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 920 TYPE(qs_energy_type), POINTER :: energy 921 TYPE(qs_ks_env_type), POINTER :: ks_env 922 TYPE(qs_rho_type), POINTER :: rho 923 TYPE(qs_scf_env_type), POINTER :: scf_env 924 TYPE(scf_control_type), POINTER :: scf_control 925 TYPE(section_vals_type), POINTER :: dft_section, input, scf_section 926 927 NULLIFY (scf_env, ks_env, energy, rho, matrix_s, rho_ao, cdft_control, logger, & 928 dft_control, pw_env, auxbas_pw_pool, energy, ks_env, scf_env, dft_section, & 929 input, scf_section, scf_control, mos, mo_occupations) 930 logger => cp_get_default_logger() 931 932 CPASSERT(ASSOCIATED(qs_env)) 933 CALL get_qs_env(qs_env, scf_env=scf_env, energy=energy, & 934 dft_control=dft_control, scf_control=scf_control, & 935 ks_env=ks_env, input=input) 936 937 CALL timeset(routineN//"_loop", handle) 938 dft_section => section_vals_get_subs_vals(input, "DFT") 939 scf_section => section_vals_get_subs_vals(dft_section, "SCF") 940 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", & 941 extension=".scfLog") 942 first_iteration = .TRUE. 943 944 cdft_control => dft_control%qs_control%cdft_control 945 946 scf_env%outer_scf%iter_count = 0 947 cdft_control%total_steps = 0 948 949 ! Write some info about the CDFT calculation 950 IF (output_unit > 0) THEN 951 WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") & 952 "CDFT EXTERNAL SCF WAVEFUNCTION OPTIMIZATION" 953 CALL qs_scf_cdft_initial_info(output_unit, cdft_control) 954 END IF 955 IF (cdft_control%reuse_precond) THEN 956 reuse_precond = .FALSE. 957 cdft_control%nreused = 0 958 END IF 959 cdft_outer_loop: DO 960 ! Change outer_scf settings to OT settings 961 CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot) 962 ! Solve electronic structure with fixed value of constraint 963 CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, & 964 converged=converged, should_stop=should_stop) 965 ! Decide whether to reuse the preconditioner on the next iteration 966 IF (cdft_control%reuse_precond) THEN 967 ! For convergence in exactly one step, the preconditioner is always reused (assuming max_reuse > 0) 968 ! usually this means that the electronic structure has already converged to the correct state 969 ! but the constraint optimizer keeps jumping over the optimal solution 970 IF (scf_env%outer_scf%iter_count == 1 .AND. scf_env%iter_count == 1 & 971 .AND. cdft_control%total_steps /= 1) & 972 cdft_control%nreused = cdft_control%nreused - 1 973 ! SCF converged in less than precond_freq steps 974 IF (scf_env%outer_scf%iter_count == 1 .AND. scf_env%iter_count .LE. cdft_control%precond_freq .AND. & 975 cdft_control%total_steps /= 1 .AND. cdft_control%nreused .LT. cdft_control%max_reuse) THEN 976 reuse_precond = .TRUE. 977 cdft_control%nreused = cdft_control%nreused + 1 978 ELSE 979 reuse_precond = .FALSE. 980 cdft_control%nreused = 0 981 END IF 982 END IF 983 ! Update history purging counters 984 IF (first_iteration .AND. cdft_control%purge_history) THEN 985 cdft_control%istep = cdft_control%istep + 1 986 IF (scf_env%outer_scf%iter_count .GT. 1) THEN 987 cdft_control%nbad_conv = cdft_control%nbad_conv + 1 988 IF (cdft_control%nbad_conv .GE. cdft_control%purge_freq .AND. & 989 cdft_control%istep .GE. cdft_control%purge_offset) THEN 990 cdft_control%nbad_conv = 0 991 cdft_control%istep = 0 992 cdft_control%should_purge = .TRUE. 993 END IF 994 END IF 995 END IF 996 first_iteration = .FALSE. 997 ! Change outer_scf settings to CDFT settings 998 CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft) 999 CALL qs_scf_check_outer_exit(qs_env, scf_env, scf_control, should_stop, & 1000 cdft_loop_converged, exit_cdft_loop) 1001 CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, & 1002 energy, cdft_control%total_steps, & 1003 should_stop, cdft_loop_converged, cdft_loop=.TRUE.) 1004 IF (exit_cdft_loop) EXIT cdft_outer_loop 1005 ! Check if the inverse Jacobian needs to be calculated 1006 CALL qs_calculate_inverse_jacobian(qs_env) 1007 ! Check if a line search should be performed to find an optimal step size for the optimizer 1008 CALL qs_cdft_line_search(qs_env) 1009 ! Optimize constraint 1010 CALL outer_loop_optimize(scf_env, scf_control) 1011 CALL outer_loop_update_qs_env(qs_env, scf_env) 1012 CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.) 1013 END DO cdft_outer_loop 1014 1015 cdft_control%ienergy = cdft_control%ienergy + 1 1016 1017 ! Store needed arrays for ET coupling calculation 1018 IF (cdft_control%do_et) THEN 1019 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, mos=mos) 1020 nvar = SIZE(cdft_control%target) 1021 ! Matrix representation of weight function 1022 ALLOCATE (cdft_control%wmat(nvar)) 1023 DO ivar = 1, nvar 1024 CALL dbcsr_init_p(cdft_control%wmat(ivar)%matrix) 1025 CALL dbcsr_copy(cdft_control%wmat(ivar)%matrix, matrix_s(1)%matrix, & 1026 name="ET_RESTRAINT_MATRIX") 1027 CALL dbcsr_set(cdft_control%wmat(ivar)%matrix, 0.0_dp) 1028 CALL integrate_v_rspace(cdft_control%group(ivar)%weight, & 1029 hmat=cdft_control%wmat(ivar), qs_env=qs_env, & 1030 calculate_forces=.FALSE., & 1031 gapw=dft_control%qs_control%gapw) 1032 END DO 1033 ! Overlap matrix 1034 CALL dbcsr_init_p(cdft_control%matrix_s%matrix) 1035 CALL dbcsr_copy(cdft_control%matrix_s%matrix, matrix_s(1)%matrix, & 1036 name="OVERLAP") 1037 ! Molecular orbital coefficients 1038 NULLIFY (cdft_control%mo_coeff) 1039 ALLOCATE (cdft_control%mo_coeff(dft_control%nspins)) 1040 DO ispin = 1, dft_control%nspins 1041 NULLIFY (cdft_control%mo_coeff(ispin)%matrix) 1042 CALL cp_fm_create(matrix=cdft_control%mo_coeff(ispin)%matrix, & 1043 matrix_struct=qs_env%mos(ispin)%mo_set%mo_coeff%matrix_struct, & 1044 name="MO_COEFF_A"//TRIM(ADJUSTL(cp_to_string(ispin)))//"MATRIX") 1045 CALL cp_fm_to_fm(qs_env%mos(ispin)%mo_set%mo_coeff, & 1046 cdft_control%mo_coeff(ispin)%matrix) 1047 END DO 1048 ! Density matrix 1049 IF (cdft_control%calculate_metric) THEN 1050 CALL get_qs_env(qs_env, rho=rho) 1051 CALL qs_rho_get(rho, rho_ao=rho_ao) 1052 ALLOCATE (cdft_control%matrix_p(dft_control%nspins)) 1053 DO ispin = 1, dft_control%nspins 1054 NULLIFY (cdft_control%matrix_p(ispin)%matrix) 1055 CALL dbcsr_init_p(cdft_control%matrix_p(ispin)%matrix) 1056 CALL dbcsr_copy(cdft_control%matrix_p(ispin)%matrix, rho_ao(ispin)%matrix, & 1057 name="DENSITY MATRIX") 1058 END DO 1059 END IF 1060 ! Copy occupation numbers if non-uniform occupation 1061 uniform_occupation = .TRUE. 1062 DO ispin = 1, dft_control%nspins 1063 CALL get_mo_set(mo_set=mos(ispin)%mo_set, uniform_occupation=my_uocc) 1064 uniform_occupation = uniform_occupation .AND. my_uocc 1065 END DO 1066 IF (.NOT. uniform_occupation) THEN 1067 ALLOCATE (cdft_control%occupations(dft_control%nspins)) 1068 DO ispin = 1, dft_control%nspins 1069 CALL get_mo_set(mo_set=mos(ispin)%mo_set, & 1070 nmo=nmo, & 1071 occupation_numbers=mo_occupations) 1072 ALLOCATE (cdft_control%occupations(ispin)%array(nmo)) 1073 cdft_control%occupations(ispin)%array(1:nmo) = mo_occupations(1:nmo) 1074 END DO 1075 END IF 1076 END IF 1077 1078 ! Deallocate constraint storage if forces are not needed 1079 ! In case of a simulation with multiple force_evals, 1080 ! deallocate only if weight function should not be copied to different force_evals 1081 IF (.NOT. (cdft_control%save_pot .OR. cdft_control%transfer_pot)) THEN 1082 CALL get_qs_env(qs_env, pw_env=pw_env) 1083 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 1084 DO iatom = 1, SIZE(cdft_control%group) 1085 CALL pw_pool_give_back_pw(auxbas_pw_pool, cdft_control%group(iatom)%weight%pw) 1086 END DO 1087 IF (cdft_control%atomic_charges) THEN 1088 DO iatom = 1, cdft_control%natoms 1089 CALL pw_pool_give_back_pw(auxbas_pw_pool, cdft_control%charge(iatom)%pw) 1090 END DO 1091 DEALLOCATE (cdft_control%charge) 1092 END IF 1093 IF (cdft_control%type == outer_scf_becke_constraint .AND. & 1094 cdft_control%becke_control%cavity_confine) THEN 1095 IF (.NOT. ASSOCIATED(cdft_control%becke_control%cavity_mat)) THEN 1096 CALL pw_pool_give_back_pw(auxbas_pw_pool, cdft_control%becke_control%cavity%pw) 1097 ELSE 1098 DEALLOCATE (cdft_control%becke_control%cavity_mat) 1099 END IF 1100 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN 1101 IF (ASSOCIATED(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)) & 1102 CALL pw_pool_give_back_pw(auxbas_pw_pool, cdft_control%hirshfeld_control%hirshfeld_env%fnorm%pw) 1103 END IF 1104 IF (ASSOCIATED(cdft_control%charges_fragment)) DEALLOCATE (cdft_control%charges_fragment) 1105 cdft_control%need_pot = .TRUE. 1106 cdft_control%external_control = .FALSE. 1107 END IF 1108 1109 CALL timestop(handle) 1110 1111 END SUBROUTINE cdft_scf 1112 1113! ************************************************************************************************** 1114!> \brief perform cleanup operations for cdft_control 1115!> \param cdft_control container for the external CDFT SCF loop variables 1116!> \par History 1117!> 12.2015 created [Nico Holmberg] 1118!> \author Nico Holmberg 1119! ************************************************************************************************** 1120 SUBROUTINE cdft_control_cleanup(cdft_control) 1121 TYPE(cdft_control_type), POINTER :: cdft_control 1122 1123 CHARACTER(len=*), PARAMETER :: routineN = 'cdft_control_cleanup', & 1124 routineP = moduleN//':'//routineN 1125 1126 IF (ASSOCIATED(cdft_control%constraint%variables)) & 1127 DEALLOCATE (cdft_control%constraint%variables) 1128 IF (ASSOCIATED(cdft_control%constraint%count)) & 1129 DEALLOCATE (cdft_control%constraint%count) 1130 IF (ASSOCIATED(cdft_control%constraint%gradient)) & 1131 DEALLOCATE (cdft_control%constraint%gradient) 1132 IF (ASSOCIATED(cdft_control%constraint%energy)) & 1133 DEALLOCATE (cdft_control%constraint%energy) 1134 IF (ASSOCIATED(cdft_control%constraint%inv_jacobian) .AND. & 1135 cdft_control%constraint%deallocate_jacobian) & 1136 DEALLOCATE (cdft_control%constraint%inv_jacobian) 1137 1138 END SUBROUTINE cdft_control_cleanup 1139 1140! ************************************************************************************************** 1141!> \brief Calculates the finite difference inverse Jacobian 1142!> \param qs_env the qs_environment_type where to compute the Jacobian 1143!> \par History 1144!> 01.2017 created [Nico Holmberg] 1145! ************************************************************************************************** 1146 SUBROUTINE qs_calculate_inverse_jacobian(qs_env) 1147 TYPE(qs_environment_type), POINTER :: qs_env 1148 1149 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_calculate_inverse_jacobian', & 1150 routineP = moduleN//':'//routineN 1151 1152 CHARACTER(len=default_path_length) :: project_name 1153 INTEGER :: counter, handle, i, ispin, iter_count, & 1154 iwork, j, max_scf, nspins, nsteps, & 1155 nvar, nwork, output_unit, pwork, twork 1156 LOGICAL :: converged, explicit_jacobian, & 1157 should_build, should_stop, & 1158 use_md_history 1159 REAL(KIND=dp) :: inv_error, step_size 1160 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coeff, dh, step_multiplier 1161 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: jacobian 1162 REAL(KIND=dp), DIMENSION(:), POINTER :: energy 1163 REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient, inv_jacobian 1164 TYPE(cdft_control_type), POINTER :: cdft_control 1165 TYPE(cp_logger_type), POINTER :: logger, tmp_logger 1166 TYPE(cp_para_env_type), POINTER :: para_env 1167 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv 1168 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp 1169 TYPE(dft_control_type), POINTER :: dft_control 1170 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos, mos_stashed 1171 TYPE(qs_energy_type), POINTER :: energy_qs 1172 TYPE(qs_ks_env_type), POINTER :: ks_env 1173 TYPE(qs_rho_type), POINTER :: rho 1174 TYPE(qs_scf_env_type), POINTER :: scf_env 1175 TYPE(scf_control_type), POINTER :: scf_control 1176 1177 NULLIFY (energy, gradient, p_rmpv, rho_ao_kp, mos, rho, mos_stashed, & 1178 ks_env, scf_env, scf_control, dft_control, cdft_control, & 1179 inv_jacobian, para_env, tmp_logger, energy_qs) 1180 logger => cp_get_default_logger() 1181 1182 CPASSERT(ASSOCIATED(qs_env)) 1183 CALL get_qs_env(qs_env, scf_env=scf_env, ks_env=ks_env, & 1184 scf_control=scf_control, mos=mos, rho=rho, & 1185 dft_control=dft_control, & 1186 para_env=para_env, energy=energy_qs) 1187 explicit_jacobian = .FALSE. 1188 should_build = .FALSE. 1189 use_md_history = .FALSE. 1190 iter_count = scf_env%outer_scf%iter_count 1191 ! Quick exit if optimizer does not require Jacobian 1192 IF (.NOT. ASSOCIATED(scf_control%outer_scf%cdft_opt_control)) RETURN 1193 ! Check if Jacobian should be calculated and initialize 1194 CALL timeset(routineN, handle) 1195 CALL initialize_inverse_jacobian(scf_control, scf_env, explicit_jacobian, should_build, used_history) 1196 IF (scf_control%outer_scf%cdft_opt_control%jacobian_restart) THEN 1197 ! Restart from previously calculated inverse Jacobian 1198 should_build = .FALSE. 1199 CALL restart_inverse_jacobian(qs_env) 1200 END IF 1201 IF (should_build) THEN 1202 scf_env%outer_scf%deallocate_jacobian = .FALSE. 1203 ! Actually need to (re)build the Jacobian 1204 IF (explicit_jacobian) THEN 1205 ! Build Jacobian with finite differences 1206 cdft_control => dft_control%qs_control%cdft_control 1207 IF (.NOT. ASSOCIATED(cdft_control)) & 1208 CALL cp_abort(__LOCATION__, & 1209 "Optimizers that need the explicit Jacobian can"// & 1210 " only be used together with a valid CDFT constraint.") 1211 ! Redirect output from Jacobian calculation to a new file by creating a temporary logger 1212 project_name = logger%iter_info%project_name 1213 CALL create_tmp_logger(para_env, project_name, "-JacobianInfo.out", output_unit, tmp_logger) 1214 ! Save last converged state so we can roll back to it (mo_coeff and some outer_loop variables) 1215 nspins = dft_control%nspins 1216 ALLOCATE (mos_stashed(nspins)) 1217 DO ispin = 1, nspins 1218 CALL duplicate_mo_set(mos_stashed(ispin)%mo_set, mos(ispin)%mo_set) 1219 END DO 1220 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp) 1221 p_rmpv => rho_ao_kp(:, 1) 1222 ! Allocate work 1223 nvar = SIZE(scf_env%outer_scf%variables, 1) 1224 max_scf = scf_control%outer_scf%max_scf + 1 1225 ALLOCATE (gradient(nvar, max_scf)) 1226 gradient = scf_env%outer_scf%gradient 1227 ALLOCATE (energy(max_scf)) 1228 energy = scf_env%outer_scf%energy 1229 ALLOCATE (jacobian(nvar, nvar)) 1230 jacobian = 0.0_dp 1231 nsteps = cdft_control%total_steps 1232 ! Setup finite difference scheme 1233 CALL prepare_jacobian_stencil(qs_env, output_unit, nwork, pwork, coeff, step_multiplier, dh) 1234 twork = pwork - nwork 1235 DO i = 1, nvar 1236 jacobian(i, :) = coeff(0)*scf_env%outer_scf%gradient(i, iter_count) 1237 END DO 1238 ! Calculate the Jacobian by perturbing each Lagrangian and recalculating the energy self-consistently 1239 CALL cp_add_default_logger(tmp_logger) 1240 DO i = 1, nvar 1241 IF (output_unit > 0) THEN 1242 WRITE (output_unit, FMT="(A)") " " 1243 WRITE (output_unit, FMT="(A)") " #####################################" 1244 WRITE (output_unit, '(A,I3,A,I3,A)') & 1245 " ### Constraint ", i, " of ", nvar, " ###" 1246 WRITE (output_unit, FMT="(A)") " #####################################" 1247 END IF 1248 counter = 0 1249 DO iwork = nwork, pwork 1250 IF (iwork == 0) CYCLE 1251 counter = counter + 1 1252 IF (output_unit > 0) THEN 1253 WRITE (output_unit, FMT="(A)") " #####################################" 1254 WRITE (output_unit, '(A,I3,A,I3,A)') & 1255 " ### Energy evaluation ", counter, " of ", twork, " ###" 1256 WRITE (output_unit, FMT="(A)") " #####################################" 1257 END IF 1258 IF (SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step) == 1) THEN 1259 step_size = scf_control%outer_scf%cdft_opt_control%jacobian_step(1) 1260 ELSE 1261 step_size = scf_control%outer_scf%cdft_opt_control%jacobian_step(i) 1262 END IF 1263 scf_env%outer_scf%variables(:, iter_count + 1) = scf_env%outer_scf%variables(:, iter_count) 1264 scf_env%outer_scf%variables(i, iter_count + 1) = scf_env%outer_scf%variables(i, iter_count) + & 1265 step_multiplier(iwork)*step_size 1266 CALL outer_loop_update_qs_env(qs_env, scf_env) 1267 CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.) 1268 CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot) 1269 CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, & 1270 converged=converged, should_stop=should_stop) 1271 CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft) 1272 ! Update (iter_count + 1) element of gradient and print constraint info 1273 scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count + 1 1274 CALL outer_loop_gradient(qs_env, scf_env) 1275 CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, & 1276 energy_qs, cdft_control%total_steps, & 1277 should_stop=.FALSE., outer_loop_converged=.FALSE., cdft_loop=.FALSE.) 1278 scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count - 1 1279 ! Update Jacobian 1280 DO j = 1, nvar 1281 jacobian(j, i) = jacobian(j, i) + coeff(iwork)*scf_env%outer_scf%gradient(j, iter_count + 1) 1282 END DO 1283 ! Reset everything to last converged state 1284 scf_env%outer_scf%variables(:, iter_count + 1) = 0.0_dp 1285 scf_env%outer_scf%gradient = gradient 1286 scf_env%outer_scf%energy = energy 1287 cdft_control%total_steps = nsteps 1288 DO ispin = 1, nspins 1289 CALL deallocate_mo_set(mos(ispin)%mo_set) 1290 CALL duplicate_mo_set(mos(ispin)%mo_set, mos_stashed(ispin)%mo_set) 1291 CALL calculate_density_matrix(mos(ispin)%mo_set, & 1292 p_rmpv(ispin)%matrix) 1293 END DO 1294 CALL qs_rho_update_rho(rho, qs_env=qs_env) 1295 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) 1296 END DO 1297 END DO 1298 CALL cp_rm_default_logger() 1299 CALL cp_logger_release(tmp_logger) 1300 ! Finalize and invert Jacobian 1301 DO j = 1, nvar 1302 DO i = 1, nvar 1303 jacobian(i, j) = jacobian(i, j)/dh(j) 1304 END DO 1305 END DO 1306 IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) & 1307 ALLOCATE (scf_env%outer_scf%inv_jacobian(nvar, nvar)) 1308 inv_jacobian => scf_env%outer_scf%inv_jacobian 1309 CALL invert_matrix(jacobian, inv_jacobian, inv_error) 1310 scf_control%outer_scf%cdft_opt_control%broyden_update = .FALSE. 1311 ! Release temporary storage 1312 DO ispin = 1, nspins 1313 CALL deallocate_mo_set(mos_stashed(ispin)%mo_set) 1314 END DO 1315 DEALLOCATE (mos_stashed, jacobian, gradient, energy, coeff, step_multiplier, dh) 1316 IF (output_unit > 0) THEN 1317 WRITE (output_unit, FMT="(/,A)") & 1318 " ================================== JACOBIAN CALCULATED ==================================" 1319 CALL close_file(unit_number=output_unit) 1320 END IF 1321 ELSE 1322 ! Build a strictly diagonal Jacobian from history and invert it 1323 CALL build_diagonal_jacobian(qs_env, used_history) 1324 END IF 1325 END IF 1326 IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian) .AND. para_env%ionode) THEN 1327 ! Write restart file for inverse Jacobian 1328 CALL print_inverse_jacobian(logger, scf_env%outer_scf%inv_jacobian, iter_count) 1329 END IF 1330 ! Update counter 1331 scf_control%outer_scf%cdft_opt_control%ijacobian(1) = scf_control%outer_scf%cdft_opt_control%ijacobian(1) + 1 1332 CALL timestop(handle) 1333 1334 END SUBROUTINE qs_calculate_inverse_jacobian 1335 1336! ************************************************************************************************** 1337!> \brief Perform backtracking line search to find the optimal step size for the CDFT constraint 1338!> optimizer. Assumes that the CDFT gradient function is a smooth function of the constraint 1339!> variables. 1340!> \param qs_env the qs_environment_type where to perform the line search 1341!> \par History 1342!> 02.2017 created [Nico Holmberg] 1343! ************************************************************************************************** 1344 SUBROUTINE qs_cdft_line_search(qs_env) 1345 TYPE(qs_environment_type), POINTER :: qs_env 1346 1347 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_cdft_line_search', & 1348 routineP = moduleN//':'//routineN 1349 1350 CHARACTER(len=default_path_length) :: project_name 1351 INTEGER :: handle, i, ispin, iter_count, & 1352 max_linesearch, max_scf, nspins, & 1353 nsteps, nvar, output_unit 1354 LOGICAL :: continue_ls, continue_ls_exit, converged, do_linesearch, found_solution, & 1355 reached_maxls, should_exit, should_stop, sign_changed 1356 LOGICAL, ALLOCATABLE, DIMENSION(:) :: positive_sign 1357 REAL(KIND=dp) :: alpha, alpha_ls, factor, norm_ls 1358 REAL(KIND=dp), DIMENSION(:), POINTER :: energy 1359 REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient, inv_jacobian 1360 REAL(KIND=dp), EXTERNAL :: dnrm2 1361 TYPE(cdft_control_type), POINTER :: cdft_control 1362 TYPE(cp_logger_type), POINTER :: logger, tmp_logger 1363 TYPE(cp_para_env_type), POINTER :: para_env 1364 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv 1365 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp 1366 TYPE(dft_control_type), POINTER :: dft_control 1367 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos, mos_ls, mos_stashed 1368 TYPE(qs_energy_type), POINTER :: energy_qs 1369 TYPE(qs_ks_env_type), POINTER :: ks_env 1370 TYPE(qs_rho_type), POINTER :: rho 1371 TYPE(qs_scf_env_type), POINTER :: scf_env 1372 TYPE(scf_control_type), POINTER :: scf_control 1373 1374 CALL timeset(routineN, handle) 1375 1376 NULLIFY (energy, gradient, p_rmpv, rho_ao_kp, mos, rho, & 1377 mos_stashed, ks_env, scf_env, scf_control, dft_control, & 1378 cdft_control, inv_jacobian, para_env, & 1379 tmp_logger, energy_qs, mos_ls) 1380 logger => cp_get_default_logger() 1381 1382 CPASSERT(ASSOCIATED(qs_env)) 1383 CALL get_qs_env(qs_env, scf_env=scf_env, ks_env=ks_env, & 1384 scf_control=scf_control, mos=mos, rho=rho, & 1385 dft_control=dft_control, & 1386 para_env=para_env, energy=energy_qs) 1387 do_linesearch = .FALSE. 1388 SELECT CASE (scf_control%outer_scf%optimizer) 1389 CASE DEFAULT 1390 do_linesearch = .FALSE. 1391 CASE (outer_scf_optimizer_newton_ls) 1392 do_linesearch = .TRUE. 1393 CASE (outer_scf_optimizer_broyden) 1394 SELECT CASE (scf_control%outer_scf%cdft_opt_control%broyden_type) 1395 CASE (broyden_type_1, broyden_type_2, broyden_type_1_explicit, broyden_type_2_explicit) 1396 do_linesearch = .FALSE. 1397 CASE (broyden_type_1_ls, broyden_type_1_explicit_ls, broyden_type_2_ls, broyden_type_2_explicit_ls) 1398 cdft_control => dft_control%qs_control%cdft_control 1399 IF (.NOT. ASSOCIATED(cdft_control)) & 1400 CALL cp_abort(__LOCATION__, & 1401 "Optimizers that perform a line search can"// & 1402 " only be used together with a valid CDFT constraint") 1403 IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) & 1404 do_linesearch = .TRUE. 1405 END SELECT 1406 END SELECT 1407 IF (do_linesearch) THEN 1408 cdft_control => dft_control%qs_control%cdft_control 1409 IF (.NOT. ASSOCIATED(cdft_control)) & 1410 CALL cp_abort(__LOCATION__, & 1411 "Optimizers that perform a line search can"// & 1412 " only be used together with a valid CDFT constraint") 1413 CPASSERT(ASSOCIATED(scf_env%outer_scf%inv_jacobian)) 1414 CPASSERT(ASSOCIATED(scf_control%outer_scf%cdft_opt_control)) 1415 alpha = scf_control%outer_scf%cdft_opt_control%newton_step_save 1416 iter_count = scf_env%outer_scf%iter_count 1417 ! Redirect output from line search procedure to a new file by creating a temporary logger 1418 project_name = logger%iter_info%project_name 1419 CALL create_tmp_logger(para_env, project_name, "-LineSearch.out", output_unit, tmp_logger) 1420 ! Save last converged state so we can roll back to it (mo_coeff and some outer_loop variables) 1421 nspins = dft_control%nspins 1422 ALLOCATE (mos_stashed(nspins)) 1423 DO ispin = 1, nspins 1424 CALL duplicate_mo_set(mos_stashed(ispin)%mo_set, mos(ispin)%mo_set) 1425 END DO 1426 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp) 1427 p_rmpv => rho_ao_kp(:, 1) 1428 nsteps = cdft_control%total_steps 1429 ! Allocate work 1430 nvar = SIZE(scf_env%outer_scf%variables, 1) 1431 max_scf = scf_control%outer_scf%max_scf + 1 1432 max_linesearch = scf_control%outer_scf%cdft_opt_control%max_ls 1433 continue_ls = scf_control%outer_scf%cdft_opt_control%continue_ls 1434 factor = scf_control%outer_scf%cdft_opt_control%factor_ls 1435 continue_ls_exit = .FALSE. 1436 found_solution = .FALSE. 1437 ALLOCATE (gradient(nvar, max_scf)) 1438 gradient = scf_env%outer_scf%gradient 1439 ALLOCATE (energy(max_scf)) 1440 energy = scf_env%outer_scf%energy 1441 reached_maxls = .FALSE. 1442 ! Broyden optimizers: perform update of inv_jacobian if necessary 1443 IF (scf_control%outer_scf%cdft_opt_control%broyden_update) THEN 1444 CALL outer_loop_optimize(scf_env, scf_control) 1445 ! Reset the variables and prevent a reupdate of inv_jacobian 1446 scf_env%outer_scf%variables(:, iter_count + 1) = 0 1447 scf_control%outer_scf%cdft_opt_control%broyden_update = .FALSE. 1448 END IF 1449 ! Print some info 1450 IF (output_unit > 0) THEN 1451 WRITE (output_unit, FMT="(/,A)") & 1452 " ================================== LINE SEARCH STARTED ==================================" 1453 WRITE (output_unit, FMT="(A,I5,A)") & 1454 " Evaluating optimal step size for optimizer using a maximum of", max_linesearch, " steps" 1455 IF (continue_ls) THEN 1456 WRITE (output_unit, FMT="(A)") & 1457 " Line search continues until best step size is found or max steps are reached" 1458 END IF 1459 WRITE (output_unit, '(/,A,F5.3)') & 1460 " Initial step size: ", alpha 1461 WRITE (output_unit, '(/,A,F5.3)') & 1462 " Step size update factor: ", factor 1463 WRITE (output_unit, '(/,A,I10,A,I10)') & 1464 " Energy evaluation: ", cdft_control%ienergy, ", CDFT SCF iteration: ", iter_count 1465 END IF 1466 ! Perform backtracking line search 1467 CALL cp_add_default_logger(tmp_logger) 1468 DO i = 1, max_linesearch 1469 IF (output_unit > 0) THEN 1470 WRITE (output_unit, FMT="(A)") " " 1471 WRITE (output_unit, FMT="(A)") " #####################################" 1472 WRITE (output_unit, '(A,I10,A)') & 1473 " ### Line search step: ", i, " ###" 1474 WRITE (output_unit, FMT="(A)") " #####################################" 1475 END IF 1476 inv_jacobian => scf_env%outer_scf%inv_jacobian 1477 ! Newton update of CDFT variables with a step size of alpha 1478 scf_env%outer_scf%variables(:, iter_count + 1) = scf_env%outer_scf%variables(:, iter_count) - alpha* & 1479 MATMUL(inv_jacobian, scf_env%outer_scf%gradient(:, iter_count)) 1480 ! With updated CDFT variables, perform SCF 1481 CALL outer_loop_update_qs_env(qs_env, scf_env) 1482 CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.) 1483 CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot) 1484 CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, & 1485 converged=converged, should_stop=should_stop) 1486 CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft) 1487 ! Update (iter_count + 1) element of gradient and print constraint info 1488 scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count + 1 1489 CALL outer_loop_gradient(qs_env, scf_env) 1490 CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, & 1491 energy_qs, cdft_control%total_steps, & 1492 should_stop=.FALSE., outer_loop_converged=.FALSE., cdft_loop=.FALSE.) 1493 scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count - 1 1494 ! Store sign of initial gradient for each variable for continue_ls 1495 IF (continue_ls .AND. .NOT. ALLOCATED(positive_sign)) THEN 1496 ALLOCATE (positive_sign(nvar)) 1497 DO ispin = 1, nvar 1498 positive_sign(ispin) = scf_env%outer_scf%gradient(ispin, iter_count + 1) .GE. 0.0_dp 1499 END DO 1500 END IF 1501 ! Check if the L2 norm of the gradient decreased 1502 inv_jacobian => scf_env%outer_scf%inv_jacobian 1503 IF (dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count + 1), 1) .LT. & 1504 dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count), 1)) THEN 1505 ! Optimal step size found 1506 IF (.NOT. continue_ls) THEN 1507 should_exit = .TRUE. 1508 ELSE 1509 ! But line search continues for at least one more iteration in an attempt to find a better solution 1510 ! if max number of steps is not exceeded 1511 IF (found_solution) THEN 1512 ! Check if the norm also decreased w.r.t. to previously found solution 1513 IF (dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count + 1), 1) .GT. norm_ls) THEN 1514 ! Norm increased => accept previous solution and exit 1515 continue_ls_exit = .TRUE. 1516 END IF 1517 END IF 1518 ! Store current state and the value of alpha 1519 IF (.NOT. continue_ls_exit) THEN 1520 should_exit = .FALSE. 1521 alpha_ls = alpha 1522 found_solution = .TRUE. 1523 norm_ls = dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count + 1), 1) 1524 ! Check if the sign of the gradient has changed for all variables (w.r.t initial gradient) 1525 ! In this case we should exit because further line search steps will just increase the norm 1526 sign_changed = .TRUE. 1527 DO ispin = 1, nvar 1528 sign_changed = sign_changed .AND. (positive_sign(ispin) .NEQV. & 1529 scf_env%outer_scf%gradient(ispin, iter_count + 1) .GE. 0.0_dp) 1530 END DO 1531 IF (.NOT. ASSOCIATED(mos_ls)) THEN 1532 ALLOCATE (mos_ls(nspins)) 1533 ELSE 1534 DO ispin = 1, nspins 1535 CALL deallocate_mo_set(mos_ls(ispin)%mo_set) 1536 END DO 1537 END IF 1538 DO ispin = 1, nspins 1539 CALL duplicate_mo_set(mos_ls(ispin)%mo_set, mos(ispin)%mo_set) 1540 END DO 1541 alpha = alpha*factor 1542 ! Exit on last iteration 1543 IF (i == max_linesearch) continue_ls_exit = .TRUE. 1544 ! Exit if constraint target is satisfied to requested tolerance 1545 IF (SQRT(MAXVAL(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count + 1)**2)) .LT. & 1546 scf_control%outer_scf%eps_scf) & 1547 continue_ls_exit = .TRUE. 1548 ! Exit if line search jumped over the optimal step length 1549 IF (sign_changed) continue_ls_exit = .TRUE. 1550 END IF 1551 END IF 1552 ELSE 1553 ! Gradient increased => alpha is too large (if the gradient function is smooth) 1554 should_exit = .FALSE. 1555 ! Update alpha using Armijo's scheme 1556 alpha = alpha*factor 1557 END IF 1558 IF (continue_ls_exit) THEN 1559 ! Continuation of line search did not yield a better alpha, use previously located solution and exit 1560 alpha = alpha_ls 1561 DO ispin = 1, nspins 1562 CALL deallocate_mo_set(mos(ispin)%mo_set) 1563 CALL duplicate_mo_set(mos(ispin)%mo_set, mos_ls(ispin)%mo_set) 1564 CALL calculate_density_matrix(mos(ispin)%mo_set, & 1565 p_rmpv(ispin)%matrix) 1566 CALL deallocate_mo_set(mos_ls(ispin)%mo_set) 1567 END DO 1568 CALL qs_rho_update_rho(rho, qs_env=qs_env) 1569 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) 1570 DEALLOCATE (mos_ls) 1571 should_exit = .TRUE. 1572 END IF 1573 ! Reached max steps and SCF converged: continue with last iterated step size 1574 IF (.NOT. should_exit .AND. & 1575 (i == max_linesearch .AND. converged .AND. .NOT. found_solution)) THEN 1576 should_exit = .TRUE. 1577 reached_maxls = .TRUE. 1578 alpha = alpha*(1.0_dp/factor) 1579 END IF 1580 ! Reset outer SCF environment to last converged state 1581 scf_env%outer_scf%variables(:, iter_count + 1) = 0.0_dp 1582 scf_env%outer_scf%gradient = gradient 1583 scf_env%outer_scf%energy = energy 1584 ! Exit line search if a suitable step size was found 1585 IF (should_exit) EXIT 1586 ! Reset the electronic structure 1587 cdft_control%total_steps = nsteps 1588 DO ispin = 1, nspins 1589 CALL deallocate_mo_set(mos(ispin)%mo_set) 1590 CALL duplicate_mo_set(mos(ispin)%mo_set, mos_stashed(ispin)%mo_set) 1591 CALL calculate_density_matrix(mos(ispin)%mo_set, & 1592 p_rmpv(ispin)%matrix) 1593 END DO 1594 CALL qs_rho_update_rho(rho, qs_env=qs_env) 1595 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) 1596 END DO 1597 scf_control%outer_scf%cdft_opt_control%newton_step = alpha 1598 IF (.NOT. should_exit) THEN 1599 CALL cp_warn(__LOCATION__, & 1600 "Line search did not converge. CDFT SCF proceeds with fixed step size.") 1601 scf_control%outer_scf%cdft_opt_control%newton_step = scf_control%outer_scf%cdft_opt_control%newton_step_save 1602 END IF 1603 IF (reached_maxls) & 1604 CALL cp_warn(__LOCATION__, & 1605 "Line search did not converge. CDFT SCF proceeds with lasted iterated step size.") 1606 CALL cp_rm_default_logger() 1607 CALL cp_logger_release(tmp_logger) 1608 ! Release temporary storage 1609 DO ispin = 1, nspins 1610 CALL deallocate_mo_set(mos_stashed(ispin)%mo_set) 1611 END DO 1612 DEALLOCATE (mos_stashed, gradient, energy) 1613 IF (ALLOCATED(positive_sign)) DEALLOCATE (positive_sign) 1614 IF (output_unit > 0) THEN 1615 WRITE (output_unit, FMT="(/,A)") & 1616 " ================================== LINE SEARCH COMPLETE ==================================" 1617 CALL close_file(unit_number=output_unit) 1618 END IF 1619 END IF 1620 1621 CALL timestop(handle) 1622 1623 END SUBROUTINE qs_cdft_line_search 1624 1625END MODULE qs_scf 1626