1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Storage of past states of the qs_env. 8!> Methods to interpolate (or actually normally extrapolate) the 9!> new guess for density and wavefunctions. 10!> \note 11!> Most of the last snapshot should actually be in qs_env, but taking 12!> advantage of it would make the programming much convoluted 13!> \par History 14!> 02.2003 created [fawzi] 15!> 11.2003 Joost VandeVondele : Implemented Nth order PS extrapolation 16!> 02.2005 modified for KG_GPW [MI] 17!> \author fawzi 18! ************************************************************************************************** 19MODULE qs_wf_history_methods 20 USE bibliography, ONLY: Kolafa2004,& 21 VandeVondele2005a,& 22 cite_reference 23 USE cp_control_types, ONLY: dft_control_type 24 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,& 25 dbcsr_allocate_matrix_set,& 26 dbcsr_deallocate_matrix_set 27 USE cp_fm_basic_linalg, ONLY: cp_fm_scale,& 28 cp_fm_scale_and_add 29 USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,& 30 cp_fm_pool_type,& 31 fm_pools_create_fm_vect,& 32 fm_pools_give_back_fm_vect 33 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 34 cp_fm_struct_release,& 35 cp_fm_struct_type 36 USE cp_fm_types, ONLY: cp_fm_create,& 37 cp_fm_get_info,& 38 cp_fm_release,& 39 cp_fm_to_fm,& 40 cp_fm_type 41 USE cp_gemm_interface, ONLY: cp_gemm 42 USE cp_log_handling, ONLY: cp_get_default_logger,& 43 cp_logger_type,& 44 cp_to_string 45 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 46 cp_print_key_unit_nr,& 47 low_print_level 48 USE dbcsr_api, ONLY: dbcsr_add,& 49 dbcsr_copy,& 50 dbcsr_deallocate_matrix,& 51 dbcsr_p_type 52 USE input_constants, ONLY: & 53 wfi_aspc_nr, wfi_frozen_method_nr, wfi_linear_p_method_nr, wfi_linear_ps_method_nr, & 54 wfi_linear_wf_method_nr, wfi_ps_method_nr, wfi_use_guess_method_nr, & 55 wfi_use_prev_p_method_nr, wfi_use_prev_rho_r_method_nr, wfi_use_prev_wf_method_nr 56 USE kinds, ONLY: dp 57 USE mathlib, ONLY: binomial 58 USE pw_env_types, ONLY: pw_env_get,& 59 pw_env_type 60 USE pw_methods, ONLY: pw_copy 61 USE pw_pool_types, ONLY: pw_pool_create_pw,& 62 pw_pool_give_back_pw,& 63 pw_pool_type 64 USE pw_types, ONLY: COMPLEXDATA1D,& 65 REALDATA3D,& 66 REALSPACE,& 67 RECIPROCALSPACE,& 68 pw_p_type 69 USE qs_environment_types, ONLY: get_qs_env,& 70 qs_environment_type,& 71 set_qs_env 72 USE qs_ks_types, ONLY: qs_ks_did_change 73 USE qs_matrix_pools, ONLY: mpools_get,& 74 qs_matrix_pools_type 75 USE qs_mo_methods, ONLY: calculate_density_matrix,& 76 make_basis_cholesky,& 77 make_basis_lowdin,& 78 make_basis_simple,& 79 make_basis_sm 80 USE qs_mo_types, ONLY: get_mo_set,& 81 mo_set_p_type 82 USE qs_rho_methods, ONLY: qs_rho_update_rho 83 USE qs_rho_types, ONLY: qs_rho_get,& 84 qs_rho_set,& 85 qs_rho_type 86 USE qs_scf_types, ONLY: ot_method_nr,& 87 qs_scf_env_type 88 USE qs_wf_history_types, ONLY: qs_wf_history_type,& 89 qs_wf_snapshot_type,& 90 wfi_get_snapshot,& 91 wfi_release 92 USE scf_control_types, ONLY: scf_control_type 93#include "./base/base_uses.f90" 94 95 IMPLICIT NONE 96 PRIVATE 97 98 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE. 99 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_wf_history_methods' 100 INTEGER, SAVE, PRIVATE :: last_wfs_id = 0, last_wfi_id = 0 101 102 PUBLIC :: wfi_create, wfi_update, wfi_create_for_kp, & 103 wfi_extrapolate, wfi_get_method_label, & 104 reorthogonalize_vectors, wfi_purge_history 105 106CONTAINS 107 108! ************************************************************************************************** 109!> \brief allocates and initialize a wavefunction snapshot 110!> \param snapshot the snapshot to create 111!> \par History 112!> 02.2003 created [fawzi] 113!> 02.2005 added wf_mol [MI] 114!> \author fawzi 115! ************************************************************************************************** 116 SUBROUTINE wfs_create(snapshot) 117 TYPE(qs_wf_snapshot_type), POINTER :: snapshot 118 119 CHARACTER(len=*), PARAMETER :: routineN = 'wfs_create', routineP = moduleN//':'//routineN 120 121 ALLOCATE (snapshot) 122 last_wfs_id = last_wfs_id + 1 123 snapshot%id_nr = last_wfs_id 124 NULLIFY (snapshot%wf, snapshot%rho_r, & 125 snapshot%rho_g, snapshot%rho_ao, snapshot%rho_ao_kp, & 126 snapshot%overlap, snapshot%rho_frozen) 127 snapshot%dt = 1.0_dp 128 snapshot%ref_count = 1 129 END SUBROUTINE wfs_create 130 131! ************************************************************************************************** 132!> \brief updates the given snapshot 133!> \param snapshot the snapshot to be updated 134!> \param wf_history the history 135!> \param qs_env the qs_env that should be snapshotted 136!> \param dt the time of the snapshot (wrt. to the previous snapshot) 137!> \par History 138!> 02.2003 created [fawzi] 139!> 02.2005 added kg_fm_mol_set for KG_GPW [MI] 140!> \author fawzi 141! ************************************************************************************************** 142 SUBROUTINE wfs_update(snapshot, wf_history, qs_env, dt) 143 TYPE(qs_wf_snapshot_type), POINTER :: snapshot 144 TYPE(qs_wf_history_type), POINTER :: wf_history 145 TYPE(qs_environment_type), POINTER :: qs_env 146 REAL(KIND=dp), INTENT(in), OPTIONAL :: dt 147 148 CHARACTER(len=*), PARAMETER :: routineN = 'wfs_update', routineP = moduleN//':'//routineN 149 150 INTEGER :: handle, img, ispin, nimg, nspins 151 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_pools 152 TYPE(cp_fm_type), POINTER :: mo_coeff 153 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao 154 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp 155 TYPE(dft_control_type), POINTER :: dft_control 156 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 157 TYPE(pw_env_type), POINTER :: pw_env 158 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g, rho_r 159 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 160 TYPE(qs_rho_type), POINTER :: rho 161 162 CALL timeset(routineN, handle) 163 164 NULLIFY (pw_env, auxbas_pw_pool, ao_mo_pools, dft_control, mos, mo_coeff, & 165 rho, rho_r, rho_g, rho_ao, matrix_s) 166 CALL get_qs_env(qs_env, pw_env=pw_env, & 167 dft_control=dft_control, rho=rho) 168 CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_pools) 169 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 170 171 CPASSERT(ASSOCIATED(wf_history)) 172 CPASSERT(ASSOCIATED(dft_control)) 173 IF (.NOT. ASSOCIATED(snapshot)) THEN 174 CALL wfs_create(snapshot) 175 END IF 176 CPASSERT(wf_history%ref_count > 0) 177 CPASSERT(snapshot%ref_count > 0) 178 179 nspins = dft_control%nspins 180 snapshot%dt = 1.0_dp 181 IF (PRESENT(dt)) snapshot%dt = dt 182 IF (wf_history%store_wf) THEN 183 CALL get_qs_env(qs_env, mos=mos) 184 IF (.NOT. ASSOCIATED(snapshot%wf)) THEN 185 CALL fm_pools_create_fm_vect(ao_mo_pools, snapshot%wf, & 186 name="ws_snap"//TRIM(ADJUSTL(cp_to_string(snapshot%id_nr)))// & 187 "ws") 188 CPASSERT(nspins == SIZE(snapshot%wf)) 189 END IF 190 DO ispin = 1, nspins 191 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff) 192 CALL cp_fm_to_fm(mo_coeff, snapshot%wf(ispin)%matrix) 193 END DO 194 ELSE IF (ASSOCIATED(snapshot%wf)) THEN 195 CALL fm_pools_give_back_fm_vect(ao_mo_pools, snapshot%wf) 196 END IF 197 198 IF (wf_history%store_rho_r) THEN 199 CALL qs_rho_get(rho, rho_r=rho_r) 200 CPASSERT(ASSOCIATED(rho_r)) 201 IF (.NOT. ASSOCIATED(snapshot%rho_r)) THEN 202 ALLOCATE (snapshot%rho_r(nspins)) 203 DO ispin = 1, nspins 204 NULLIFY (snapshot%rho_r(ispin)%pw) 205 CALL pw_pool_create_pw(auxbas_pw_pool, snapshot%rho_r(ispin)%pw, & 206 in_space=REALSPACE, use_data=REALDATA3D) 207 END DO 208 END IF 209 DO ispin = 1, nspins 210 CALL pw_copy(rho_r(ispin)%pw, snapshot%rho_r(ispin)%pw) 211 END DO 212 ELSE IF (ASSOCIATED(snapshot%rho_r)) THEN 213 DO ispin = 1, SIZE(snapshot%rho_r) 214 CALL pw_pool_give_back_pw(auxbas_pw_pool, snapshot%rho_r(ispin)%pw) 215 END DO 216 DEALLOCATE (snapshot%rho_r) 217 END IF 218 219 IF (wf_history%store_rho_g) THEN 220 CALL qs_rho_get(rho, rho_g=rho_g) 221 CPASSERT(ASSOCIATED(rho_g)) 222 IF (.NOT. ASSOCIATED(snapshot%rho_g)) THEN 223 ALLOCATE (snapshot%rho_g(nspins)) 224 DO ispin = 1, nspins 225 NULLIFY (snapshot%rho_g(ispin)%pw) 226 CALL pw_pool_create_pw(auxbas_pw_pool, snapshot%rho_g(ispin)%pw, & 227 in_space=RECIPROCALSPACE, use_data=COMPLEXDATA1D) 228 END DO 229 END IF 230 DO ispin = 1, nspins 231 CALL pw_copy(rho_g(ispin)%pw, snapshot%rho_g(ispin)%pw) 232 END DO 233 ELSE IF (ASSOCIATED(snapshot%rho_g)) THEN 234 DO ispin = 1, SIZE(snapshot%rho_g) 235 CALL pw_pool_give_back_pw(auxbas_pw_pool, snapshot%rho_g(ispin)%pw) 236 END DO 237 DEALLOCATE (snapshot%rho_g) 238 END IF 239 240 IF (ASSOCIATED(snapshot%rho_ao)) THEN ! the sparsity might be different 241 ! (future struct:check) 242 CALL dbcsr_deallocate_matrix_set(snapshot%rho_ao) 243 END IF 244 IF (wf_history%store_rho_ao) THEN 245 CALL qs_rho_get(rho, rho_ao=rho_ao) 246 CPASSERT(ASSOCIATED(rho_ao)) 247 248 CALL dbcsr_allocate_matrix_set(snapshot%rho_ao, nspins) 249 DO ispin = 1, nspins 250 ALLOCATE (snapshot%rho_ao(ispin)%matrix) 251 CALL dbcsr_copy(snapshot%rho_ao(ispin)%matrix, rho_ao(ispin)%matrix) 252 END DO 253 END IF 254 255 IF (ASSOCIATED(snapshot%rho_ao_kp)) THEN ! the sparsity might be different 256 ! (future struct:check) 257 CALL dbcsr_deallocate_matrix_set(snapshot%rho_ao_kp) 258 END IF 259 IF (wf_history%store_rho_ao_kp) THEN 260 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp) 261 CPASSERT(ASSOCIATED(rho_ao_kp)) 262 263 nimg = dft_control%nimages 264 CALL dbcsr_allocate_matrix_set(snapshot%rho_ao_kp, nspins, nimg) 265 DO ispin = 1, nspins 266 DO img = 1, nimg 267 ALLOCATE (snapshot%rho_ao_kp(ispin, img)%matrix) 268 CALL dbcsr_copy(snapshot%rho_ao_kp(ispin, img)%matrix, & 269 rho_ao_kp(ispin, img)%matrix) 270 END DO 271 END DO 272 END IF 273 274 IF (ASSOCIATED(snapshot%overlap)) THEN ! the sparsity might be different 275 ! (future struct:check) 276 CALL dbcsr_deallocate_matrix(snapshot%overlap) 277 END IF 278 IF (wf_history%store_overlap) THEN 279 CALL get_qs_env(qs_env, matrix_s=matrix_s) 280 CPASSERT(ASSOCIATED(matrix_s)) 281 CPASSERT(ASSOCIATED(matrix_s(1)%matrix)) 282 ALLOCATE (snapshot%overlap) 283 CALL dbcsr_copy(snapshot%overlap, matrix_s(1)%matrix) 284 END IF 285 286 IF (wf_history%store_frozen_density) THEN 287 ! do nothing 288 ! CALL deallocate_matrix_set(snapshot%rho_frozen%rho_ao) 289 END IF 290 291 CALL timestop(handle) 292 293 END SUBROUTINE wfs_update 294 295! ************************************************************************************************** 296!> \brief ... 297!> \param wf_history ... 298!> \param interpolation_method_nr the tag of the method used for 299!> the extrapolation of the initial density for the next md step 300!> (see qs_wf_history_types:wfi_*_method_nr) 301!> \param extrapolation_order ... 302!> \param has_unit_metric ... 303!> \par History 304!> 02.2003 created [fawzi] 305!> \author fawzi 306! ************************************************************************************************** 307 SUBROUTINE wfi_create(wf_history, interpolation_method_nr, extrapolation_order, & 308 has_unit_metric) 309 TYPE(qs_wf_history_type), POINTER :: wf_history 310 INTEGER, INTENT(in) :: interpolation_method_nr, & 311 extrapolation_order 312 LOGICAL, INTENT(IN) :: has_unit_metric 313 314 CHARACTER(len=*), PARAMETER :: routineN = 'wfi_create', routineP = moduleN//':'//routineN 315 316 INTEGER :: handle, i 317 318 CALL timeset(routineN, handle) 319 320 ALLOCATE (wf_history) 321 last_wfi_id = last_wfi_id + 1 322 wf_history%id_nr = last_wfi_id 323 wf_history%ref_count = 1 324 wf_history%memory_depth = 0 325 wf_history%snapshot_count = 0 326 wf_history%last_state_index = 1 327 wf_history%store_wf = .FALSE. 328 wf_history%store_rho_r = .FALSE. 329 wf_history%store_rho_g = .FALSE. 330 wf_history%store_rho_ao = .FALSE. 331 wf_history%store_rho_ao_kp = .FALSE. 332 wf_history%store_overlap = .FALSE. 333 wf_history%store_frozen_density = .FALSE. 334 NULLIFY (wf_history%past_states) 335 336 wf_history%interpolation_method_nr = interpolation_method_nr 337 338 SELECT CASE (wf_history%interpolation_method_nr) 339 CASE (wfi_use_guess_method_nr) 340 wf_history%memory_depth = 0 341 CASE (wfi_use_prev_wf_method_nr) 342 wf_history%memory_depth = 0 343 CASE (wfi_use_prev_p_method_nr) 344 wf_history%memory_depth = 1 345 wf_history%store_rho_ao = .TRUE. 346 CASE (wfi_use_prev_rho_r_method_nr) 347 wf_history%memory_depth = 1 348 wf_history%store_rho_ao = .TRUE. 349 CASE (wfi_linear_wf_method_nr) 350 wf_history%memory_depth = 2 351 wf_history%store_wf = .TRUE. 352 CASE (wfi_linear_p_method_nr) 353 wf_history%memory_depth = 2 354 wf_history%store_rho_ao = .TRUE. 355 CASE (wfi_linear_ps_method_nr) 356 wf_history%memory_depth = 2 357 wf_history%store_wf = .TRUE. 358 IF (.NOT. has_unit_metric) wf_history%store_overlap = .TRUE. 359 CASE (wfi_ps_method_nr) 360 CALL cite_reference(VandeVondele2005a) 361 wf_history%memory_depth = extrapolation_order + 1 362 wf_history%store_wf = .TRUE. 363 IF (.NOT. has_unit_metric) wf_history%store_overlap = .TRUE. 364 CASE (wfi_frozen_method_nr) 365 wf_history%memory_depth = 1 366 wf_history%store_frozen_density = .TRUE. 367 CASE (wfi_aspc_nr) 368 wf_history%memory_depth = extrapolation_order + 2 369 wf_history%store_wf = .TRUE. 370 IF (.NOT. has_unit_metric) wf_history%store_overlap = .TRUE. 371 CASE default 372 CALL cp_abort(__LOCATION__, & 373 "Unknown interpolation method: "// & 374 TRIM(ADJUSTL(cp_to_string(interpolation_method_nr)))) 375 wf_history%interpolation_method_nr = wfi_use_prev_rho_r_method_nr 376 END SELECT 377 ALLOCATE (wf_history%past_states(wf_history%memory_depth)) 378 379 DO i = 1, SIZE(wf_history%past_states) 380 NULLIFY (wf_history%past_states(i)%snapshot) 381 END DO 382 383 CALL timestop(handle) 384 END SUBROUTINE wfi_create 385 386! ************************************************************************************************** 387!> \brief ... 388!> \param wf_history ... 389!> \par History 390!> 06.2015 created [jhu] 391!> \author fawzi 392! ************************************************************************************************** 393 SUBROUTINE wfi_create_for_kp(wf_history) 394 TYPE(qs_wf_history_type), POINTER :: wf_history 395 396 CHARACTER(len=*), PARAMETER :: routineN = 'wfi_create_for_kp', & 397 routineP = moduleN//':'//routineN 398 399 CPASSERT(ASSOCIATED(wf_history)) 400 IF (wf_history%store_rho_ao) THEN 401 wf_history%store_rho_ao_kp = .TRUE. 402 wf_history%store_rho_ao = .FALSE. 403 END IF 404 ! Check for KP compatible methods 405 IF (wf_history%store_wf) THEN 406 CPABORT("WFN based interpolation method not possible for kpoints.") 407 END IF 408 IF (wf_history%store_frozen_density) THEN 409 CPABORT("Frozen density initialization method not possible for kpoints.") 410 END IF 411 IF (wf_history%store_overlap) THEN 412 CPABORT("Inconsistent interpolation method for kpoints.") 413 END IF 414 415 END SUBROUTINE wfi_create_for_kp 416 417! ************************************************************************************************** 418!> \brief returns a string describing the interpolation method 419!> \param method_nr ... 420!> \return ... 421!> \par History 422!> 02.2003 created [fawzi] 423!> \author fawzi 424! ************************************************************************************************** 425 FUNCTION wfi_get_method_label(method_nr) RESULT(res) 426 INTEGER, INTENT(in) :: method_nr 427 CHARACTER(len=30) :: res 428 429 CHARACTER(len=*), PARAMETER :: routineN = 'wfi_get_method_label', & 430 routineP = moduleN//':'//routineN 431 432 res = "unknown" 433 SELECT CASE (method_nr) 434 CASE (wfi_use_prev_p_method_nr) 435 res = "previous_p" 436 CASE (wfi_use_prev_wf_method_nr) 437 res = "previous_wf" 438 CASE (wfi_use_prev_rho_r_method_nr) 439 res = "previous_rho_r" 440 CASE (wfi_use_guess_method_nr) 441 res = "initial_guess" 442 CASE (wfi_linear_wf_method_nr) 443 res = "mo linear" 444 CASE (wfi_linear_p_method_nr) 445 res = "P linear" 446 CASE (wfi_linear_ps_method_nr) 447 res = "PS linear" 448 CASE (wfi_ps_method_nr) 449 res = "PS Nth order" 450 CASE (wfi_frozen_method_nr) 451 res = "frozen density approximation" 452 CASE (wfi_aspc_nr) 453 res = "ASPC" 454 CASE default 455 CALL cp_abort(__LOCATION__, & 456 "Unknown interpolation method: "// & 457 TRIM(ADJUSTL(cp_to_string(method_nr)))) 458 END SELECT 459 END FUNCTION wfi_get_method_label 460 461! ************************************************************************************************** 462!> \brief calculates the new starting state for the scf for the next 463!> wf optimization 464!> \param wf_history the previous history needed to extrapolate 465!> \param qs_env the qs env with the latest result, and that will contain 466!> the new starting state 467!> \param dt the time at which to extrapolate (wrt. to the last snapshot) 468!> \param extrapolation_method_nr returns the extrapolation method used 469!> \param orthogonal_wf ... 470!> \par History 471!> 02.2003 created [fawzi] 472!> 11.2003 Joost VandeVondele : Implemented Nth order PS extrapolation 473!> \author fawzi 474! ************************************************************************************************** 475 SUBROUTINE wfi_extrapolate(wf_history, qs_env, dt, extrapolation_method_nr, & 476 orthogonal_wf) 477 TYPE(qs_wf_history_type), POINTER :: wf_history 478 TYPE(qs_environment_type), POINTER :: qs_env 479 REAL(KIND=dp), INTENT(IN) :: dt 480 INTEGER, INTENT(OUT), OPTIONAL :: extrapolation_method_nr 481 LOGICAL, INTENT(OUT), OPTIONAL :: orthogonal_wf 482 483 CHARACTER(len=*), PARAMETER :: routineN = 'wfi_extrapolate', & 484 routineP = moduleN//':'//routineN 485 486 INTEGER :: actual_extrapolation_method_nr, handle, & 487 i, img, io_unit, ispin, k, n, nmo, & 488 nvec, print_level 489 LOGICAL :: do_kpoints, my_orthogonal_wf, use_overlap 490 REAL(KIND=dp) :: alpha, t0, t1, t2 491 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools 492 TYPE(cp_fm_struct_type), POINTER :: matrix_struct, matrix_struct_new 493 TYPE(cp_fm_type), POINTER :: csc, fm_tmp, mo_coeff 494 TYPE(cp_logger_type), POINTER :: logger 495 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao, rho_frozen_ao 496 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp 497 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 498 TYPE(qs_rho_type), POINTER :: rho 499 TYPE(qs_wf_snapshot_type), POINTER :: t0_state, t1_state 500 501 NULLIFY (mos, ao_mo_fm_pools, t0_state, t1_state, mo_coeff, & 502 rho, rho_ao, rho_frozen_ao) 503 504 use_overlap = wf_history%store_overlap 505 506 CALL timeset(routineN, handle) 507 logger => cp_get_default_logger() 508 print_level = logger%iter_info%print_level 509 io_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%SCF%PRINT%PROGRAM_RUN_INFO", & 510 extension=".scfLog") 511 512 CPASSERT(ASSOCIATED(wf_history)) 513 CPASSERT(wf_history%ref_count > 0) 514 CPASSERT(ASSOCIATED(qs_env)) 515 CPASSERT(qs_env%ref_count > 0) 516 CALL get_qs_env(qs_env, mos=mos, rho=rho, do_kpoints=do_kpoints) 517 CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools) 518 ! chooses the method for this extrapolation 519 IF (wf_history%snapshot_count < 1) THEN 520 actual_extrapolation_method_nr = wfi_use_guess_method_nr 521 ELSE 522 actual_extrapolation_method_nr = wf_history%interpolation_method_nr 523 END IF 524 525 SELECT CASE (actual_extrapolation_method_nr) 526 CASE (wfi_linear_wf_method_nr) 527 IF (wf_history%snapshot_count < 2) THEN 528 actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr 529 END IF 530 CASE (wfi_linear_p_method_nr) 531 IF (wf_history%snapshot_count < 2) THEN 532 IF (do_kpoints) THEN 533 actual_extrapolation_method_nr = wfi_use_guess_method_nr 534 ELSE 535 actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr 536 ENDIF 537 END IF 538 CASE (wfi_linear_ps_method_nr) 539 IF (wf_history%snapshot_count < 2) THEN 540 actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr 541 END IF 542 END SELECT 543 544 IF (PRESENT(extrapolation_method_nr)) & 545 extrapolation_method_nr = actual_extrapolation_method_nr 546 my_orthogonal_wf = .FALSE. 547 548 SELECT CASE (actual_extrapolation_method_nr) 549 CASE (wfi_frozen_method_nr) 550 CPASSERT(.NOT. do_kpoints) 551 t0_state => wfi_get_snapshot(wf_history, wf_index=1) 552 CPASSERT(ASSOCIATED(t0_state%rho_frozen)) 553 554 nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count) 555 CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec) 556 557 CALL qs_rho_get(t0_state%rho_frozen, rho_ao=rho_frozen_ao) 558 CALL qs_rho_get(rho, rho_ao=rho_ao) 559 DO ispin = 1, SIZE(rho_frozen_ao) 560 CALL dbcsr_copy(rho_ao(ispin)%matrix, & 561 rho_frozen_ao(ispin)%matrix, & 562 keep_sparsity=.TRUE.) 563 END DO 564 !FM updating rho_ao directly with t0_state%rho_ao would have the 565 !FM wrong matrix structure 566 CALL qs_rho_update_rho(rho, qs_env=qs_env) 567 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) 568 569 my_orthogonal_wf = .FALSE. 570 CASE (wfi_use_prev_rho_r_method_nr) 571 t0_state => wfi_get_snapshot(wf_history, wf_index=1) 572 nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count) 573 CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec) 574 IF (do_kpoints) THEN 575 CPASSERT(ASSOCIATED(t0_state%rho_ao_kp)) 576 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp) 577 DO ispin = 1, SIZE(t0_state%rho_ao_kp, 1) 578 DO img = 1, SIZE(t0_state%rho_ao_kp, 2) 579 IF (img > SIZE(rho_ao_kp, 2)) THEN 580 CPWARN("Change in cell neighborlist: might affect quality of initial guess") 581 ELSE 582 CALL dbcsr_copy(rho_ao_kp(ispin, img)%matrix, & 583 t0_state%rho_ao_kp(ispin, img)%matrix, & 584 keep_sparsity=.TRUE.) 585 END IF 586 END DO 587 END DO 588 ELSE 589 CPASSERT(ASSOCIATED(t0_state%rho_ao)) 590 CALL qs_rho_get(rho, rho_ao=rho_ao) 591 DO ispin = 1, SIZE(t0_state%rho_ao) 592 CALL dbcsr_copy(rho_ao(ispin)%matrix, & 593 t0_state%rho_ao(ispin)%matrix, & 594 keep_sparsity=.TRUE.) 595 END DO 596 END IF 597 ! Why is rho_g valid at this point ? 598 CALL qs_rho_set(rho, rho_g_valid=.TRUE.) 599 600 ! does nothing 601 CASE (wfi_use_prev_wf_method_nr) 602 CPASSERT(.NOT. do_kpoints) 603 my_orthogonal_wf = .TRUE. 604 nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count) 605 CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec) 606 CALL qs_rho_get(rho, rho_ao=rho_ao) 607 DO ispin = 1, SIZE(mos) 608 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff, & 609 nmo=nmo) 610 CALL reorthogonalize_vectors(qs_env, & 611 v_matrix=mo_coeff, & 612 n_col=nmo) 613 CALL calculate_density_matrix(mo_set=mos(ispin)%mo_set, & 614 density_matrix=rho_ao(ispin)%matrix) 615 END DO 616 CALL qs_rho_update_rho(rho, qs_env=qs_env) 617 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) 618 619 CASE (wfi_use_prev_p_method_nr) 620 t0_state => wfi_get_snapshot(wf_history, wf_index=1) 621 nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count) 622 CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec) 623 IF (do_kpoints) THEN 624 CPASSERT(ASSOCIATED(t0_state%rho_ao_kp)) 625 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp) 626 DO ispin = 1, SIZE(t0_state%rho_ao_kp, 1) 627 DO img = 1, SIZE(t0_state%rho_ao_kp, 2) 628 IF (img > SIZE(rho_ao_kp, 2)) THEN 629 CPWARN("Change in cell neighborlist: might affect quality of initial guess") 630 ELSE 631 CALL dbcsr_copy(rho_ao_kp(ispin, img)%matrix, & 632 t0_state%rho_ao_kp(ispin, img)%matrix, & 633 keep_sparsity=.TRUE.) 634 END IF 635 END DO 636 END DO 637 ELSE 638 CPASSERT(ASSOCIATED(t0_state%rho_ao)) 639 CALL qs_rho_get(rho, rho_ao=rho_ao) 640 DO ispin = 1, SIZE(t0_state%rho_ao) 641 CALL dbcsr_copy(rho_ao(ispin)%matrix, & 642 t0_state%rho_ao(ispin)%matrix, & 643 keep_sparsity=.TRUE.) 644 END DO 645 END IF 646 !FM updating rho_ao directly with t0_state%rho_ao would have the 647 !FM wrong matrix structure 648 CALL qs_rho_update_rho(rho, qs_env=qs_env) 649 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) 650 651 CASE (wfi_use_guess_method_nr) 652 !FM more clean to do it here, but it 653 !FM might need to read a file (restart) and thus globenv 654 !FM I do not want globenv here, thus done by the caller 655 !FM (btw. it also needs the eigensolver, and unless you relocate it 656 !FM gives circular dependencies) 657 nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count) 658 CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec) 659 CASE (wfi_linear_wf_method_nr) 660 CPASSERT(.NOT. do_kpoints) 661 t0_state => wfi_get_snapshot(wf_history, wf_index=2) 662 t1_state => wfi_get_snapshot(wf_history, wf_index=1) 663 CPASSERT(ASSOCIATED(t0_state)) 664 CPASSERT(ASSOCIATED(t1_state)) 665 CPASSERT(ASSOCIATED(t0_state%wf)) 666 CPASSERT(ASSOCIATED(t1_state%wf)) 667 nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count) 668 CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec) 669 670 my_orthogonal_wf = .TRUE. 671 t0 = 0.0_dp 672 t1 = t1_state%dt 673 t2 = t1 + dt 674 CALL qs_rho_get(rho, rho_ao=rho_ao) 675 DO ispin = 1, SIZE(mos) 676 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff, & 677 nmo=nmo) 678 CALL cp_fm_scale_and_add(alpha=0.0_dp, & 679 matrix_a=mo_coeff, & 680 matrix_b=t1_state%wf(ispin)%matrix, & 681 beta=(t2 - t0)/(t1 - t0)) 682 ! this copy should be unnecessary 683 CALL cp_fm_scale_and_add(alpha=1.0_dp, & 684 matrix_a=mo_coeff, & 685 beta=(t1 - t2)/(t1 - t0), matrix_b=t0_state%wf(ispin)%matrix) 686 CALL reorthogonalize_vectors(qs_env, & 687 v_matrix=mo_coeff, & 688 n_col=nmo) 689 CALL calculate_density_matrix(mo_set=mos(ispin)%mo_set, & 690 density_matrix=rho_ao(ispin)%matrix) 691 END DO 692 CALL qs_rho_update_rho(rho, qs_env=qs_env) 693 694 CALL qs_ks_did_change(qs_env%ks_env, & 695 rho_changed=.TRUE.) 696 CASE (wfi_linear_p_method_nr) 697 t0_state => wfi_get_snapshot(wf_history, wf_index=2) 698 t1_state => wfi_get_snapshot(wf_history, wf_index=1) 699 CPASSERT(ASSOCIATED(t0_state)) 700 CPASSERT(ASSOCIATED(t1_state)) 701 IF (do_kpoints) THEN 702 CPASSERT(ASSOCIATED(t0_state%rho_ao_kp)) 703 CPASSERT(ASSOCIATED(t1_state%rho_ao_kp)) 704 ELSE 705 CPASSERT(ASSOCIATED(t0_state%rho_ao)) 706 CPASSERT(ASSOCIATED(t1_state%rho_ao)) 707 END IF 708 nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count) 709 CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec) 710 711 t0 = 0.0_dp 712 t1 = t1_state%dt 713 t2 = t1 + dt 714 IF (do_kpoints) THEN 715 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp) 716 DO ispin = 1, SIZE(rho_ao_kp, 1) 717 DO img = 1, SIZE(rho_ao_kp, 2) 718 IF (img > SIZE(t0_state%rho_ao_kp, 2) .OR. & 719 img > SIZE(t1_state%rho_ao_kp, 2)) THEN 720 CPWARN("Change in cell neighborlist: might affect quality of initial guess") 721 ELSE 722 CALL dbcsr_add(rho_ao_kp(ispin, img)%matrix, t1_state%rho_ao_kp(ispin, img)%matrix, & 723 alpha_scalar=0.0_dp, beta_scalar=(t2 - t0)/(t1 - t0)) ! this copy should be unnecessary 724 CALL dbcsr_add(rho_ao_kp(ispin, img)%matrix, t0_state%rho_ao_kp(ispin, img)%matrix, & 725 alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0)) 726 END IF 727 END DO 728 END DO 729 ELSE 730 CALL qs_rho_get(rho, rho_ao=rho_ao) 731 DO ispin = 1, SIZE(rho_ao) 732 CALL dbcsr_add(rho_ao(ispin)%matrix, t1_state%rho_ao(ispin)%matrix, & 733 alpha_scalar=0.0_dp, beta_scalar=(t2 - t0)/(t1 - t0)) ! this copy should be unnecessary 734 CALL dbcsr_add(rho_ao(ispin)%matrix, t0_state%rho_ao(ispin)%matrix, & 735 alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0)) 736 END DO 737 END IF 738 CALL qs_rho_update_rho(rho, qs_env=qs_env) 739 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) 740 CASE (wfi_linear_ps_method_nr) 741 ! wf not calculated, extract with PSC renormalized? 742 ! use wf_linear? 743 CPASSERT(.NOT. do_kpoints) 744 t0_state => wfi_get_snapshot(wf_history, wf_index=2) 745 t1_state => wfi_get_snapshot(wf_history, wf_index=1) 746 CPASSERT(ASSOCIATED(t0_state)) 747 CPASSERT(ASSOCIATED(t1_state)) 748 CPASSERT(ASSOCIATED(t0_state%wf)) 749 CPASSERT(ASSOCIATED(t1_state%wf)) 750 IF (wf_history%store_overlap) THEN 751 CPASSERT(ASSOCIATED(t0_state%overlap)) 752 CPASSERT(ASSOCIATED(t1_state%overlap)) 753 END IF 754 nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count) 755 IF (nvec >= wf_history%memory_depth) THEN 756 IF ((qs_env%scf_control%max_scf_hist .NE. 0) .AND. (qs_env%scf_control%eps_scf_hist .NE. 0)) THEN 757 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist 758 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist 759 qs_env%scf_control%outer_scf%have_scf = .FALSE. 760 ELSE IF (qs_env%scf_control%max_scf_hist .NE. 0) THEN 761 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist 762 qs_env%scf_control%outer_scf%have_scf = .FALSE. 763 ELSE IF (qs_env%scf_control%eps_scf_hist .NE. 0) THEN 764 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist 765 END IF 766 END IF 767 768 my_orthogonal_wf = .TRUE. 769 ! use PS_2=2 PS_1-PS_0 770 ! C_2 comes from using PS_2 as a projector acting on C_1 771 CALL qs_rho_get(rho, rho_ao=rho_ao) 772 DO ispin = 1, SIZE(mos) 773 NULLIFY (mo_coeff, matrix_struct, matrix_struct_new, csc) 774 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff) 775 CALL cp_fm_get_info(mo_coeff, nrow_global=n, ncol_global=k, & 776 matrix_struct=matrix_struct) 777 CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, & 778 nrow_global=k, ncol_global=k) 779 CALL cp_fm_create(csc, matrix_struct_new) 780 CALL cp_fm_struct_release(matrix_struct_new) 781 782 IF (use_overlap) THEN 783 CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin)%matrix, mo_coeff, k) 784 CALL cp_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin)%matrix, mo_coeff, 0.0_dp, csc) 785 ELSE 786 CALL cp_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin)%matrix, & 787 t1_state%wf(ispin)%matrix, 0.0_dp, csc) 788 END IF 789 CALL cp_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin)%matrix, csc, 0.0_dp, mo_coeff) 790 CALL cp_fm_release(csc) 791 CALL cp_fm_scale_and_add(-1.0_dp, mo_coeff, 2.0_dp, t1_state%wf(ispin)%matrix) 792 CALL reorthogonalize_vectors(qs_env, & 793 v_matrix=mo_coeff, & 794 n_col=k) 795 CALL calculate_density_matrix(mo_set=mos(ispin)%mo_set, & 796 density_matrix=rho_ao(ispin)%matrix) 797 END DO 798 CALL qs_rho_update_rho(rho, qs_env=qs_env) 799 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) 800 801 CASE (wfi_ps_method_nr) 802 CPASSERT(.NOT. do_kpoints) 803 ! figure out the actual number of vectors to use in the extrapolation: 804 nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count) 805 CPASSERT(nvec .GT. 0) 806 IF (nvec >= wf_history%memory_depth) THEN 807 IF ((qs_env%scf_control%max_scf_hist .NE. 0) .AND. (qs_env%scf_control%eps_scf_hist .NE. 0)) THEN 808 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist 809 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist 810 qs_env%scf_control%outer_scf%have_scf = .FALSE. 811 ELSE IF (qs_env%scf_control%max_scf_hist .NE. 0) THEN 812 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist 813 qs_env%scf_control%outer_scf%have_scf = .FALSE. 814 ELSE IF (qs_env%scf_control%eps_scf_hist .NE. 0) THEN 815 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist 816 END IF 817 END IF 818 819 my_orthogonal_wf = .TRUE. 820 DO ispin = 1, SIZE(mos) 821 NULLIFY (mo_coeff, matrix_struct, matrix_struct_new, csc, fm_tmp) 822 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff) 823 CALL cp_fm_get_info(mo_coeff, nrow_global=n, ncol_global=k, & 824 matrix_struct=matrix_struct) 825 CALL cp_fm_create(fm_tmp, matrix_struct) 826 CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, & 827 nrow_global=k, ncol_global=k) 828 CALL cp_fm_create(csc, matrix_struct_new) 829 CALL cp_fm_struct_release(matrix_struct_new) 830 ! first the most recent 831 t1_state => wfi_get_snapshot(wf_history, wf_index=1) 832 CALL cp_fm_to_fm(t1_state%wf(ispin)%matrix, mo_coeff) 833 alpha = nvec 834 CALL cp_fm_scale(alpha, mo_coeff) 835 CALL qs_rho_get(rho, rho_ao=rho_ao) 836 DO i = 2, nvec 837 t0_state => wfi_get_snapshot(wf_history, wf_index=i) 838 IF (use_overlap) THEN 839 CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin)%matrix, fm_tmp, k) 840 CALL cp_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin)%matrix, fm_tmp, 0.0_dp, csc) 841 ELSE 842 CALL cp_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin)%matrix, & 843 t1_state%wf(ispin)%matrix, 0.0_dp, csc) 844 END IF 845 CALL cp_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin)%matrix, csc, 0.0_dp, fm_tmp) 846 alpha = -1.0_dp*alpha*REAL(nvec - i + 1, dp)/REAL(i, dp) 847 CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp) 848 ENDDO 849 850 CALL cp_fm_release(csc) 851 CALL cp_fm_release(fm_tmp) 852 CALL reorthogonalize_vectors(qs_env, & 853 v_matrix=mo_coeff, & 854 n_col=k) 855 CALL calculate_density_matrix(mo_set=mos(ispin)%mo_set, & 856 density_matrix=rho_ao(ispin)%matrix) 857 END DO 858 CALL qs_rho_update_rho(rho, qs_env=qs_env) 859 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) 860 861 CASE (wfi_aspc_nr) 862 CPASSERT(.NOT. do_kpoints) 863 CALL cite_reference(Kolafa2004) 864 ! figure out the actual number of vectors to use in the extrapolation: 865 nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count) 866 CPASSERT(nvec .GT. 0) 867 IF (nvec >= wf_history%memory_depth) THEN 868 IF ((qs_env%scf_control%max_scf_hist .NE. 0) .AND. & 869 (qs_env%scf_control%eps_scf_hist .NE. 0)) THEN 870 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist 871 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist 872 qs_env%scf_control%outer_scf%have_scf = .FALSE. 873 ELSE IF (qs_env%scf_control%max_scf_hist .NE. 0) THEN 874 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist 875 qs_env%scf_control%outer_scf%have_scf = .FALSE. 876 ELSE IF (qs_env%scf_control%eps_scf_hist .NE. 0) THEN 877 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist 878 END IF 879 END IF 880 881 my_orthogonal_wf = .TRUE. 882 CALL qs_rho_get(rho, rho_ao=rho_ao) 883 DO ispin = 1, SIZE(mos) 884 NULLIFY (mo_coeff, matrix_struct, matrix_struct_new, csc, fm_tmp) 885 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff) 886 CALL cp_fm_get_info(mo_coeff, & 887 nrow_global=n, & 888 ncol_global=k, & 889 matrix_struct=matrix_struct) 890 CALL cp_fm_create(fm_tmp, matrix_struct) 891 CALL cp_fm_struct_create(matrix_struct_new, & 892 template_fmstruct=matrix_struct, & 893 nrow_global=k, & 894 ncol_global=k) 895 CALL cp_fm_create(csc, matrix_struct_new) 896 CALL cp_fm_struct_release(matrix_struct_new) 897 ! first the most recent 898 t1_state => wfi_get_snapshot(wf_history, & 899 wf_index=1) 900 CALL cp_fm_to_fm(t1_state%wf(ispin)%matrix, mo_coeff) 901 alpha = REAL(4*nvec - 2, KIND=dp)/REAL(nvec + 1, KIND=dp) 902 IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN 903 WRITE (UNIT=io_unit, FMT="(/,T2,A,/,/,T3,A,I0,/,/,T3,A2,I0,A4,F10.6)") & 904 "Parameters for the always stable predictor-corrector (ASPC) method:", & 905 "ASPC order: ", MAX(nvec - 2, 0), & 906 "B(", 1, ") = ", alpha 907 END IF 908 CALL cp_fm_scale(alpha, mo_coeff) 909 910 DO i = 2, nvec 911 t0_state => wfi_get_snapshot(wf_history, wf_index=i) 912 IF (use_overlap) THEN 913 CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin)%matrix, fm_tmp, k) 914 CALL cp_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin)%matrix, fm_tmp, 0.0_dp, csc) 915 ELSE 916 CALL cp_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin)%matrix, & 917 t1_state%wf(ispin)%matrix, 0.0_dp, csc) 918 END IF 919 CALL cp_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin)%matrix, csc, 0.0_dp, fm_tmp) 920 alpha = (-1.0_dp)**(i + 1)*REAL(i, KIND=dp)* & 921 binomial(2*nvec, nvec - i)/binomial(2*nvec - 2, nvec - 1) 922 IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN 923 WRITE (UNIT=io_unit, FMT="(T3,A2,I0,A4,F10.6)") & 924 "B(", i, ") = ", alpha 925 END IF 926 CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp) 927 END DO 928 CALL cp_fm_release(csc) 929 CALL cp_fm_release(fm_tmp) 930 CALL reorthogonalize_vectors(qs_env, & 931 v_matrix=mo_coeff, & 932 n_col=k) 933 CALL calculate_density_matrix(mo_set=mos(ispin)%mo_set, & 934 density_matrix=rho_ao(ispin)%matrix) 935 END DO 936 CALL qs_rho_update_rho(rho, qs_env=qs_env) 937 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) 938 939 CASE default 940 CALL cp_abort(__LOCATION__, & 941 "Unknown interpolation method: "// & 942 TRIM(ADJUSTL(cp_to_string(wf_history%interpolation_method_nr)))) 943 END SELECT 944 IF (PRESENT(orthogonal_wf)) orthogonal_wf = my_orthogonal_wf 945 CALL cp_print_key_finished_output(io_unit, logger, qs_env%input, & 946 "DFT%SCF%PRINT%PROGRAM_RUN_INFO") 947 CALL timestop(handle) 948 END SUBROUTINE wfi_extrapolate 949 950! ************************************************************************************************** 951!> \brief Decides if scf control variables has to changed due 952!> to using a WF extrapolation. 953!> \param qs_env The QS environment 954!> \param nvec ... 955!> \par History 956!> 11.2006 created [TdK] 957!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) 958! ************************************************************************************************** 959 SUBROUTINE wfi_set_history_variables(qs_env, nvec) 960 TYPE(qs_environment_type), POINTER :: qs_env 961 INTEGER, INTENT(IN) :: nvec 962 963 CHARACTER(len=*), PARAMETER :: routineN = 'wfi_set_history_variables', & 964 routineP = moduleN//':'//routineN 965 966 INTEGER :: handle 967 968 CALL timeset(routineN, handle) 969 970 CPASSERT(ASSOCIATED(qs_env)) 971 CPASSERT(qs_env%ref_count > 0) 972 973 IF (nvec >= qs_env%wf_history%memory_depth) THEN 974 IF ((qs_env%scf_control%max_scf_hist .NE. 0) .AND. (qs_env%scf_control%eps_scf_hist .NE. 0)) THEN 975 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist 976 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist 977 qs_env%scf_control%outer_scf%have_scf = .FALSE. 978 ELSE IF (qs_env%scf_control%max_scf_hist .NE. 0) THEN 979 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist 980 qs_env%scf_control%outer_scf%have_scf = .FALSE. 981 ELSE IF (qs_env%scf_control%eps_scf_hist .NE. 0) THEN 982 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist 983 qs_env%scf_control%outer_scf%eps_scf = qs_env%scf_control%eps_scf_hist 984 END IF 985 END IF 986 987 CALL timestop(handle) 988 989 END SUBROUTINE wfi_set_history_variables 990 991! ************************************************************************************************** 992!> \brief updates the snapshot buffer, taking a new snapshot 993!> \param wf_history the history buffer to update 994!> \param qs_env the qs_env we get the info from 995!> \param dt ... 996!> \par History 997!> 02.2003 created [fawzi] 998!> \author fawzi 999! ************************************************************************************************** 1000 SUBROUTINE wfi_update(wf_history, qs_env, dt) 1001 TYPE(qs_wf_history_type), POINTER :: wf_history 1002 TYPE(qs_environment_type), POINTER :: qs_env 1003 REAL(KIND=dp), INTENT(in) :: dt 1004 1005 CHARACTER(len=*), PARAMETER :: routineN = 'wfi_update', routineP = moduleN//':'//routineN 1006 1007 CPASSERT(ASSOCIATED(wf_history)) 1008 CPASSERT(wf_history%ref_count > 0) 1009 CPASSERT(ASSOCIATED(qs_env)) 1010 CPASSERT(qs_env%ref_count > 0) 1011 1012 wf_history%snapshot_count = wf_history%snapshot_count + 1 1013 IF (wf_history%memory_depth > 0) THEN 1014 wf_history%last_state_index = MODULO(wf_history%snapshot_count, & 1015 wf_history%memory_depth) + 1 1016 CALL wfs_update(snapshot=wf_history%past_states & 1017 (wf_history%last_state_index)%snapshot, wf_history=wf_history, & 1018 qs_env=qs_env, dt=dt) 1019 END IF 1020 END SUBROUTINE wfi_update 1021 1022! ************************************************************************************************** 1023!> \brief reorthogonalizes the mos 1024!> \param qs_env the qs_env in which to orthogonalize 1025!> \param v_matrix the vectors to orthogonalize 1026!> \param n_col number of column of v to orthogonalize 1027!> \par History 1028!> 04.2003 created [fawzi] 1029!> \author Fawzi Mohamed 1030! ************************************************************************************************** 1031 SUBROUTINE reorthogonalize_vectors(qs_env, v_matrix, n_col) 1032 TYPE(qs_environment_type), POINTER :: qs_env 1033 TYPE(cp_fm_type), POINTER :: v_matrix 1034 INTEGER, INTENT(in), OPTIONAL :: n_col 1035 1036 CHARACTER(len=*), PARAMETER :: routineN = 'reorthogonalize_vectors', & 1037 routineP = moduleN//':'//routineN 1038 1039 INTEGER :: handle, my_n_col 1040 LOGICAL :: has_unit_metric, & 1041 ortho_contains_cholesky, & 1042 smearing_is_used 1043 TYPE(cp_fm_pool_type), POINTER :: maxao_maxmo_fm_pool 1044 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 1045 TYPE(dft_control_type), POINTER :: dft_control 1046 TYPE(qs_matrix_pools_type), POINTER :: mpools 1047 TYPE(qs_scf_env_type), POINTER :: scf_env 1048 TYPE(scf_control_type), POINTER :: scf_control 1049 1050 NULLIFY (scf_env, scf_control, maxao_maxmo_fm_pool, matrix_s, mpools, dft_control) 1051 CALL timeset(routineN, handle) 1052 1053 CPASSERT(ASSOCIATED(qs_env)) 1054 CPASSERT(ASSOCIATED(v_matrix)) 1055 1056 CALL cp_fm_get_info(v_matrix, ncol_global=my_n_col) 1057 IF (PRESENT(n_col)) my_n_col = n_col 1058 CALL get_qs_env(qs_env, mpools=mpools, & 1059 scf_env=scf_env, & 1060 scf_control=scf_control, & 1061 matrix_s=matrix_s, & 1062 dft_control=dft_control) 1063 CALL mpools_get(mpools, maxao_maxmo_fm_pool=maxao_maxmo_fm_pool) 1064 IF (ASSOCIATED(scf_env)) THEN 1065 ortho_contains_cholesky = (scf_env%method /= ot_method_nr) .AND. & 1066 (scf_env%cholesky_method > 0) .AND. & 1067 ASSOCIATED(scf_env%ortho) 1068 ELSE 1069 ortho_contains_cholesky = .FALSE. 1070 END IF 1071 1072 CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric) 1073 smearing_is_used = .FALSE. 1074 IF (dft_control%smear) THEN 1075 smearing_is_used = .TRUE. 1076 END IF 1077 1078 IF (has_unit_metric) THEN 1079 CALL make_basis_simple(v_matrix, my_n_col) 1080 ELSE IF (smearing_is_used) THEN 1081 CALL make_basis_lowdin(vmatrix=v_matrix, ncol=my_n_col, & 1082 matrix_s=matrix_s(1)%matrix) 1083 ELSE IF (ortho_contains_cholesky) THEN 1084 CALL make_basis_cholesky(vmatrix=v_matrix, ncol=my_n_col, & 1085 ortho=scf_env%ortho) 1086 ELSE 1087 CALL make_basis_sm(v_matrix, my_n_col, matrix_s(1)%matrix) 1088 END IF 1089 CALL timestop(handle) 1090 END SUBROUTINE reorthogonalize_vectors 1091 1092! ************************************************************************************************** 1093!> \brief purges wf_history retaining only the latest snapshot 1094!> \param qs_env the qs env with the latest result, and that will contain 1095!> the purged wf_history 1096!> \par History 1097!> 05.2016 created [Nico Holmberg] 1098!> \author Nico Holmberg 1099! ************************************************************************************************** 1100 SUBROUTINE wfi_purge_history(qs_env) 1101 TYPE(qs_environment_type), POINTER :: qs_env 1102 1103 CHARACTER(len=*), PARAMETER :: routineN = 'wfi_purge_history', & 1104 routineP = moduleN//':'//routineN 1105 1106 INTEGER :: handle, io_unit, print_level 1107 TYPE(cp_logger_type), POINTER :: logger 1108 TYPE(dft_control_type), POINTER :: dft_control 1109 TYPE(qs_wf_history_type), POINTER :: wf_history 1110 1111 NULLIFY (dft_control, wf_history) 1112 1113 CALL timeset(routineN, handle) 1114 logger => cp_get_default_logger() 1115 print_level = logger%iter_info%print_level 1116 io_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%SCF%PRINT%PROGRAM_RUN_INFO", & 1117 extension=".scfLog") 1118 1119 CPASSERT(ASSOCIATED(qs_env)) 1120 CPASSERT(qs_env%ref_count > 0) 1121 CPASSERT(ASSOCIATED(qs_env%wf_history)) 1122 CPASSERT(qs_env%wf_history%ref_count > 0) 1123 CALL get_qs_env(qs_env, dft_control=dft_control) 1124 1125 SELECT CASE (qs_env%wf_history%interpolation_method_nr) 1126 CASE (wfi_use_guess_method_nr, wfi_use_prev_wf_method_nr, & 1127 wfi_use_prev_p_method_nr, wfi_use_prev_rho_r_method_nr, & 1128 wfi_frozen_method_nr) 1129 ! do nothing 1130 CASE (wfi_linear_wf_method_nr, wfi_linear_p_method_nr, & 1131 wfi_linear_ps_method_nr, wfi_ps_method_nr, & 1132 wfi_aspc_nr) 1133 IF (qs_env%wf_history%snapshot_count .GE. 2) THEN 1134 IF (debug_this_module .AND. io_unit > 0) & 1135 WRITE (io_unit, FMT="(T2,A)") "QS| Purging WFN history" 1136 CALL wfi_create(wf_history, interpolation_method_nr= & 1137 dft_control%qs_control%wf_interpolation_method_nr, & 1138 extrapolation_order=dft_control%qs_control%wf_extrapolation_order, & 1139 has_unit_metric=qs_env%has_unit_metric) 1140 CALL set_qs_env(qs_env=qs_env, & 1141 wf_history=wf_history) 1142 CALL wfi_release(wf_history) 1143 CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp) 1144 END IF 1145 CASE DEFAULT 1146 CPABORT("Unknown extrapolation method.") 1147 END SELECT 1148 CALL timestop(handle) 1149 1150 END SUBROUTINE wfi_purge_history 1151 1152END MODULE qs_wf_history_methods 1153