1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Apply the direct inversion in the iterative subspace (DIIS) of Pulay 8!> in the framework of an SCF iteration for convergence acceleration 9!> \par Literature 10!> - P. Pulay, Chem. Phys. Lett. 73, 393 (1980) 11!> - P. Pulay, J. Comput. Chem. 3, 556 (1982) 12!> \par History 13!> - Changed to BLACS matrix usage (08.06.2001,MK) 14!> - rewritten to include LSD (1st attempt) (01.2003, Joost VandeVondele) 15!> - DIIS for ROKS (05.04.06,MK) 16!> \author Matthias Krack (28.06.2000) 17! ************************************************************************************************** 18MODULE qs_diis 19 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm 20 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,& 21 cp_fm_scale_and_add,& 22 cp_fm_symm,& 23 cp_fm_trace 24 USE cp_fm_struct, ONLY: cp_fm_struct_type 25 USE cp_fm_types, ONLY: cp_fm_create,& 26 cp_fm_get_info,& 27 cp_fm_maxabsval,& 28 cp_fm_p_type,& 29 cp_fm_set_all,& 30 cp_fm_to_fm,& 31 cp_fm_type 32 USE cp_gemm_interface, ONLY: cp_gemm 33 USE cp_log_handling, ONLY: cp_get_default_logger,& 34 cp_logger_type,& 35 cp_to_string 36 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 37 cp_print_key_unit_nr 38 USE cp_para_types, ONLY: cp_para_env_type 39 USE dbcsr_api, ONLY: & 40 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_dot, dbcsr_maxabs, dbcsr_multiply, & 41 dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_transposed, dbcsr_type 42 USE dm_ls_scf_types, ONLY: ls_scf_env_type 43 USE input_section_types, ONLY: section_vals_type 44 USE kinds, ONLY: default_string_length,& 45 dp 46 USE mathlib, ONLY: diamat_all 47 USE qs_diis_types, ONLY: qs_diis_buffer_type,& 48 qs_diis_buffer_type_sparse 49 USE qs_environment_types, ONLY: get_qs_env,& 50 qs_environment_type 51 USE qs_mo_types, ONLY: get_mo_set,& 52 mo_set_p_type 53 USE string_utilities, ONLY: compress 54#include "./base/base_uses.f90" 55 56 IMPLICIT NONE 57 58 PRIVATE 59 60 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_diis' 61 INTEGER, SAVE :: last_diis_b_id = 0 62 63 ! Public subroutines 64 65 PUBLIC :: qs_diis_b_clear, & 66 qs_diis_b_create, & 67 qs_diis_b_step 68 PUBLIC :: qs_diis_b_clear_sparse, & 69 qs_diis_b_create_sparse, & 70 qs_diis_b_step_4lscf 71 72CONTAINS 73 74! ************************************************************************************************** 75!> \brief Allocates an SCF DIIS buffer 76!> \param diis_buffer the buffer to create 77!> \param nbuffer ... 78!> \par History 79!> 02.2003 created [fawzi] 80!> \author fawzi 81! ************************************************************************************************** 82 SUBROUTINE qs_diis_b_create(diis_buffer, nbuffer) 83 84 TYPE(qs_diis_buffer_type), POINTER :: diis_buffer 85 INTEGER, INTENT(in) :: nbuffer 86 87 CHARACTER(len=*), PARAMETER :: routineN = 'qs_diis_b_create', & 88 routineP = moduleN//':'//routineN 89 90 INTEGER :: handle 91 92! ------------------------------------------------------------------------- 93 94 CALL timeset(routineN, handle) 95 96 ALLOCATE (diis_buffer) 97 98 NULLIFY (diis_buffer%b_matrix) 99 NULLIFY (diis_buffer%error) 100 NULLIFY (diis_buffer%parameter) 101 diis_buffer%nbuffer = nbuffer 102 diis_buffer%ncall = 0 103 last_diis_b_id = last_diis_b_id + 1 104 diis_buffer%id_nr = last_diis_b_id 105 diis_buffer%ref_count = 1 106 107 CALL timestop(handle) 108 109 END SUBROUTINE qs_diis_b_create 110 111! ************************************************************************************************** 112!> \brief Allocate and initialize a DIIS buffer for nao*nao parameter 113!> variables and with a buffer size of nbuffer. 114!> \param diis_buffer the buffer to initialize 115!> \param matrix_struct the structure for the matrix of the buffer 116!> \param nspin ... 117!> \param scf_section ... 118!> \par History 119!> - Creation (07.05.2001, Matthias Krack) 120!> - Changed to BLACS matrix usage (08.06.2001,MK) 121!> - DIIS for ROKS (05.04.06,MK) 122!> \author Matthias Krack 123!> \note 124!> check to allocate matrixes only when needed, using a linked list? 125! ************************************************************************************************** 126 SUBROUTINE qs_diis_b_check_i_alloc(diis_buffer, matrix_struct, nspin, & 127 scf_section) 128 129 TYPE(qs_diis_buffer_type), POINTER :: diis_buffer 130 TYPE(cp_fm_struct_type), POINTER :: matrix_struct 131 INTEGER, INTENT(IN) :: nspin 132 TYPE(section_vals_type), POINTER :: scf_section 133 134 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_check_i_alloc', & 135 routineP = moduleN//':'//routineN 136 137 INTEGER :: handle, ibuffer, ispin, nbuffer, & 138 output_unit 139 TYPE(cp_logger_type), POINTER :: logger 140 141! ------------------------------------------------------------------------- 142 143 CALL timeset(routineN, handle) 144 145 logger => cp_get_default_logger() 146 147 CPASSERT(ASSOCIATED(diis_buffer)) 148 CPASSERT(diis_buffer%ref_count > 0) 149 150 nbuffer = diis_buffer%nbuffer 151 152 IF (.NOT. ASSOCIATED(diis_buffer%error)) THEN 153 ALLOCATE (diis_buffer%error(nbuffer, nspin)) 154 155 DO ispin = 1, nspin 156 DO ibuffer = 1, nbuffer 157 NULLIFY (diis_buffer%error(ibuffer, ispin)%matrix) 158 CALL cp_fm_create(diis_buffer%error(ibuffer, ispin)%matrix, & 159 name="qs_diis_b"// & 160 TRIM(ADJUSTL(cp_to_string(diis_buffer%id_nr)))// & 161 "%error("// & 162 TRIM(ADJUSTL(cp_to_string(ibuffer)))//","// & 163 TRIM(ADJUSTL(cp_to_string(ibuffer)))//")", & 164 matrix_struct=matrix_struct) 165 END DO 166 END DO 167 END IF 168 169 IF (.NOT. ASSOCIATED(diis_buffer%parameter)) THEN 170 ALLOCATE (diis_buffer%parameter(nbuffer, nspin)) 171 172 DO ispin = 1, nspin 173 DO ibuffer = 1, nbuffer 174 NULLIFY (diis_buffer%parameter(ibuffer, ispin)%matrix) 175 CALL cp_fm_create(diis_buffer%parameter(ibuffer, ispin)%matrix, & 176 name="qs_diis_b"// & 177 TRIM(ADJUSTL(cp_to_string(diis_buffer%id_nr)))// & 178 "%parameter("// & 179 TRIM(ADJUSTL(cp_to_string(ibuffer)))//","// & 180 TRIM(ADJUSTL(cp_to_string(ibuffer)))//")", & 181 matrix_struct=matrix_struct) 182 END DO 183 END DO 184 END IF 185 186 IF (.NOT. ASSOCIATED(diis_buffer%b_matrix)) THEN 187 ALLOCATE (diis_buffer%b_matrix(nbuffer + 1, nbuffer + 1)) 188 diis_buffer%b_matrix = 0.0_dp 189 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", & 190 extension=".scfLog") 191 IF (output_unit > 0) THEN 192 WRITE (UNIT=output_unit, FMT="(/,T9,A)") & 193 "DIIS | The SCF DIIS buffer was allocated and initialized" 194 END IF 195 CALL cp_print_key_finished_output(output_unit, logger, scf_section, & 196 "PRINT%DIIS_INFO") 197 END IF 198 199 CALL timestop(handle) 200 201 END SUBROUTINE qs_diis_b_check_i_alloc 202 203! ************************************************************************************************** 204!> \brief Update the SCF DIIS buffer, and if appropriate does a diis step. 205!> \param diis_buffer ... 206!> \param mo_array ... 207!> \param kc ... 208!> \param sc ... 209!> \param delta ... 210!> \param error_max ... 211!> \param diis_step ... 212!> \param eps_diis ... 213!> \param nmixing ... 214!> \param s_matrix ... 215!> \param scf_section ... 216!> \param roks ... 217!> \par History 218!> - Creation (07.05.2001, Matthias Krack) 219!> - Changed to BLACS matrix usage (08.06.2001, MK) 220!> - 03.2003 rewamped [fawzi] 221!> - Adapted for high-spin ROKS (08.04.06,MK) 222!> \author Matthias Krack 223! ************************************************************************************************** 224 SUBROUTINE qs_diis_b_step(diis_buffer, mo_array, kc, sc, delta, error_max, & 225 diis_step, eps_diis, nmixing, s_matrix, scf_section, roks) 226 227 TYPE(qs_diis_buffer_type), POINTER :: diis_buffer 228 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mo_array 229 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: kc 230 TYPE(cp_fm_type), POINTER :: sc 231 REAL(KIND=dp), INTENT(IN) :: delta 232 REAL(KIND=dp), INTENT(OUT) :: error_max 233 LOGICAL, INTENT(OUT) :: diis_step 234 REAL(KIND=dp), INTENT(IN) :: eps_diis 235 INTEGER, INTENT(IN), OPTIONAL :: nmixing 236 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, & 237 POINTER :: s_matrix 238 TYPE(section_vals_type), POINTER :: scf_section 239 LOGICAL, INTENT(IN), OPTIONAL :: roks 240 241 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_step', routineP = moduleN//':'//routineN 242 REAL(KIND=dp), PARAMETER :: eigenvalue_threshold = 1.0E-12_dp 243 244 CHARACTER(LEN=2*default_string_length) :: message 245 INTEGER :: handle, homo, ib, imo, ispin, jb, & 246 my_nmixing, nao, nb, nb1, nmo, nspin, & 247 output_unit 248 LOGICAL :: eigenvectors_discarded, mo_uocc, my_roks 249 REAL(KIND=dp) :: maxocc, tmp 250 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ev, occ 251 REAL(KIND=dp), DIMENSION(:), POINTER :: occa, occb 252 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b 253 TYPE(cp_fm_struct_type), POINTER :: matrix_struct 254 TYPE(cp_fm_type), POINTER :: c, new_errors, old_errors, parameters 255 TYPE(cp_logger_type), POINTER :: logger 256 257! ------------------------------------------------------------------------- 258 259 CALL timeset(routineN, handle) 260 261 nspin = SIZE(mo_array) 262 diis_step = .FALSE. 263 264 IF (PRESENT(roks)) THEN 265 my_roks = .TRUE. 266 nspin = 1 267 ELSE 268 my_roks = .FALSE. 269 END IF 270 271 my_nmixing = 2 272 IF (PRESENT(nmixing)) my_nmixing = nmixing 273 274 NULLIFY (c, new_errors, old_errors, parameters, matrix_struct, a, b, occa, occb) 275 logger => cp_get_default_logger() 276 277 ! Quick return, if no DIIS is requested 278 279 IF (diis_buffer%nbuffer < 1) THEN 280 CALL timestop(handle) 281 RETURN 282 END IF 283 284 CALL cp_fm_get_info(kc(1)%matrix, & 285 matrix_struct=matrix_struct) 286 CALL qs_diis_b_check_i_alloc(diis_buffer, & 287 matrix_struct=matrix_struct, & 288 nspin=nspin, & 289 scf_section=scf_section) 290 291 error_max = 0.0_dp 292 293 ib = MODULO(diis_buffer%ncall, diis_buffer%nbuffer) + 1 294 diis_buffer%ncall = diis_buffer%ncall + 1 295 nb = MIN(diis_buffer%ncall, diis_buffer%nbuffer) 296 297 DO ispin = 1, nspin 298 299 CALL get_mo_set(mo_set=mo_array(ispin)%mo_set, & 300 nao=nao, & 301 nmo=nmo, & 302 homo=homo, & 303 mo_coeff=c, & 304 occupation_numbers=occa, & 305 uniform_occupation=mo_uocc, & 306 maxocc=maxocc) 307 308 new_errors => diis_buffer%error(ib, ispin)%matrix 309 parameters => diis_buffer%parameter(ib, ispin)%matrix 310 311 ! Copy the Kohn-Sham matrix K to the DIIS buffer 312 313 CALL cp_fm_to_fm(kc(ispin)%matrix, parameters) 314 315 IF (my_roks) THEN 316 317 ALLOCATE (occ(nmo)) 318 319 CALL get_mo_set(mo_set=mo_array(2)%mo_set, & 320 occupation_numbers=occb) 321 322 DO imo = 1, nmo 323 occ(imo) = SQRT(occa(imo) + occb(imo)) 324 END DO 325 326 CALL cp_fm_to_fm(c, sc) 327 CALL cp_fm_column_scale(sc, occ(1:homo)) 328 329 ! KC <- K*C 330 CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, parameters, sc, 0.0_dp, kc(ispin)%matrix) 331 332 IF (PRESENT(s_matrix)) THEN 333 CALL copy_dbcsr_to_fm(s_matrix(1)%matrix, new_errors) 334 ! SC <- S*C 335 CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, new_errors, c, 0.0_dp, sc) 336 CALL cp_fm_column_scale(sc, occ(1:homo)) 337 END IF 338 339 ! new_errors <- KC*(SC)^T - (SC)*(KC)^T = K*P*S - S*P*K 340 ! or for an orthogonal basis 341 ! new_errors <- KC*C^T - C*(KC)^T = K*P - P*K with S = I 342 CALL cp_gemm("N", "T", nao, nao, homo, 1.0_dp, sc, kc(ispin)%matrix, 0.0_dp, new_errors) 343 CALL cp_gemm("N", "T", nao, nao, homo, 1.0_dp, kc(ispin)%matrix, sc, -1.0_dp, new_errors) 344 345 DEALLOCATE (occ) 346 347 ELSE 348 349 ! KC <- K*C 350 CALL cp_fm_symm("L", "U", nao, homo, maxocc, parameters, c, 0.0_dp, kc(ispin)%matrix) 351 352 IF (PRESENT(s_matrix)) THEN 353 ! I guess that this copy can be avoided for LSD 354 CALL copy_dbcsr_to_fm(s_matrix(1)%matrix, new_errors) 355 ! sc <- S*C 356 CALL cp_fm_symm("L", "U", nao, homo, 2.0_dp, new_errors, c, 0.0_dp, sc) 357 ! new_errors <- KC*(SC)^T - (SC)*(KC)^T = K*P*S - S*P*K 358 CALL cp_gemm("N", "T", nao, nao, homo, 1.0_dp, sc, kc(ispin)%matrix, 0.0_dp, new_errors) 359 CALL cp_gemm("N", "T", nao, nao, homo, 1.0_dp, kc(ispin)%matrix, sc, -1.0_dp, new_errors) 360 ELSE 361 ! new_errors <- KC*(C)^T - C*(KC)^T = K*P - P*K 362 CALL cp_gemm("N", "T", nao, nao, homo, 1.0_dp, c, kc(ispin)%matrix, 0.0_dp, new_errors) 363 CALL cp_gemm("N", "T", nao, nao, homo, 1.0_dp, kc(ispin)%matrix, c, -1.0_dp, new_errors) 364 END IF 365 366 END IF 367 368 CALL cp_fm_maxabsval(new_errors, tmp) 369 error_max = MAX(error_max, tmp) 370 371 END DO 372 373 ! Check, if a DIIS step is appropiate 374 375 diis_step = ((diis_buffer%ncall >= my_nmixing) .AND. (delta < eps_diis)) 376 377 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", & 378 extension=".scfLog") 379 IF (output_unit > 0) THEN 380 WRITE (UNIT=output_unit, FMT="(/,T9,A,I4,/,(T9,A,ES12.3))") & 381 "DIIS | Current SCF DIIS buffer size: ", nb, & 382 "DIIS | Maximum SCF DIIS error vector element:", error_max, & 383 "DIIS | Current SCF convergence: ", delta, & 384 "DIIS | Threshold value for a DIIS step: ", eps_diis 385 IF (error_max < eps_diis) THEN 386 WRITE (UNIT=output_unit, FMT="(T9,A)") & 387 "DIIS | => The SCF DIIS buffer will be updated" 388 ELSE 389 WRITE (UNIT=output_unit, FMT="(T9,A)") & 390 "DIIS | => No update of the SCF DIIS buffer" 391 END IF 392 IF (diis_step .AND. (error_max < eps_diis)) THEN 393 WRITE (UNIT=output_unit, FMT="(T9,A,/)") & 394 "DIIS | => A SCF DIIS step will be performed" 395 ELSE 396 WRITE (UNIT=output_unit, FMT="(T9,A,/)") & 397 "DIIS | => No SCF DIIS step will be performed" 398 END IF 399 END IF 400 401 ! Update the SCF DIIS buffer 402 403 IF (error_max < eps_diis) THEN 404 405 b => diis_buffer%b_matrix 406 407 DO jb = 1, nb 408 b(jb, ib) = 0.0_dp 409 DO ispin = 1, nspin 410 old_errors => diis_buffer%error(jb, ispin)%matrix 411 new_errors => diis_buffer%error(ib, ispin)%matrix 412 CALL cp_fm_trace(old_errors, new_errors, tmp) 413 b(jb, ib) = b(jb, ib) + tmp 414 END DO 415 b(ib, jb) = b(jb, ib) 416 END DO 417 418 ELSE 419 420 diis_step = .FALSE. 421 422 END IF 423 424 ! Perform DIIS step 425 426 IF (diis_step) THEN 427 428 nb1 = nb + 1 429 430 ALLOCATE (a(nb1, nb1)) 431 ALLOCATE (b(nb1, nb1)) 432 ALLOCATE (ev(nb1)) 433 434 ! Set up the linear DIIS equation system 435 436 b(1:nb, 1:nb) = diis_buffer%b_matrix(1:nb, 1:nb) 437 438 b(1:nb, nb1) = -1.0_dp 439 b(nb1, 1:nb) = -1.0_dp 440 b(nb1, nb1) = 0.0_dp 441 442 ! Solve the linear DIIS equation system 443 444 ev(1:nb1) = 0.0_dp 445 CALL diamat_all(b(1:nb1, 1:nb1), ev(1:nb1)) 446 447 a(1:nb1, 1:nb1) = b(1:nb1, 1:nb1) 448 449 eigenvectors_discarded = .FALSE. 450 451 DO jb = 1, nb1 452 IF (ABS(ev(jb)) < eigenvalue_threshold) THEN 453 IF (output_unit > 0) THEN 454 IF (.NOT. eigenvectors_discarded) THEN 455 WRITE (UNIT=output_unit, FMT="(T9,A)") & 456 "DIIS | Checking eigenvalues of the DIIS error matrix" 457 END IF 458 WRITE (UNIT=message, FMT="(T9,A,I6,A,ES10.1,A,ES10.1)") & 459 "DIIS | Eigenvalue ", jb, " = ", ev(jb), " is smaller than "// & 460 "threshold ", eigenvalue_threshold 461 CALL compress(message) 462 WRITE (UNIT=output_unit, FMT="(T9,A)") TRIM(message) 463 eigenvectors_discarded = .TRUE. 464 END IF 465 a(1:nb1, jb) = 0.0_dp 466 ELSE 467 a(1:nb1, jb) = a(1:nb1, jb)/ev(jb) 468 END IF 469 END DO 470 471 IF ((output_unit > 0) .AND. eigenvectors_discarded) THEN 472 WRITE (UNIT=output_unit, FMT="(T9,A,/)") & 473 "DIIS | The corresponding eigenvectors were discarded" 474 END IF 475 476 ev(1:nb) = MATMUL(a(1:nb, 1:nb1), b(nb1, 1:nb1)) 477 478 ! Update Kohn-Sham matrix 479 480 DO ispin = 1, nspin 481 CALL cp_fm_set_all(kc(ispin)%matrix, 0.0_dp) 482 DO jb = 1, nb 483 parameters => diis_buffer%parameter(jb, ispin)%matrix 484 CALL cp_fm_scale_and_add(1.0_dp, kc(ispin)%matrix, -ev(jb), parameters) 485 END DO 486 END DO 487 488 DEALLOCATE (a) 489 DEALLOCATE (b) 490 DEALLOCATE (ev) 491 492 ELSE 493 494 DO ispin = 1, nspin 495 parameters => diis_buffer%parameter(ib, ispin)%matrix 496 CALL cp_fm_to_fm(parameters, kc(ispin)%matrix) 497 END DO 498 499 END IF 500 501 CALL cp_print_key_finished_output(output_unit, logger, scf_section, & 502 "PRINT%DIIS_INFO") 503 504 CALL timestop(handle) 505 506 END SUBROUTINE qs_diis_b_step 507 508! ************************************************************************************************** 509!> \brief clears the buffer 510!> \param diis_buffer the buffer to clear 511!> \par History 512!> 02.2003 created [fawzi] 513!> \author fawzi 514! ************************************************************************************************** 515 SUBROUTINE qs_diis_b_clear(diis_buffer) 516 517 TYPE(qs_diis_buffer_type), POINTER :: diis_buffer 518 519 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_clear', & 520 routineP = moduleN//':'//routineN 521 522 INTEGER :: handle 523 524! ------------------------------------------------------------------------- 525 526 CALL timeset(routineN, handle) 527 528 CPASSERT(ASSOCIATED(diis_buffer)) 529 CPASSERT(diis_buffer%ref_count > 0) 530 531 diis_buffer%ncall = 0 532 533 CALL timestop(handle) 534 535 END SUBROUTINE qs_diis_b_clear 536 537! ************************************************************************************************** 538!> \brief Update the SCF DIIS buffer in linear scaling SCF (LS-SCF), 539!> and if appropriate does a diis step. 540!> \param diis_buffer ... 541!> \param qs_env ... 542!> \param ls_scf_env ... 543!> \param unit_nr ... 544!> \param iscf ... 545!> \param diis_step ... 546!> \param eps_diis ... 547!> \param nmixing ... 548!> \param s_matrix ... 549!> \param threshold ... 550!> \par History 551!> - Adapted for LS-SCF (10-11-14) from qs_diis_b_step 552!> \author Fredy W. Aquino 553! ************************************************************************************************** 554 555 SUBROUTINE qs_diis_b_step_4lscf(diis_buffer, qs_env, ls_scf_env, unit_nr, iscf, & 556 diis_step, eps_diis, nmixing, s_matrix, threshold) 557! Note.- Input: ls_scf_env%matrix_p(ispin) , Density Matrix 558! matrix_ks (from qs_env) , Kohn-Sham Matrix (IN/OUT) 559 560 TYPE(qs_diis_buffer_type_sparse), POINTER :: diis_buffer 561 TYPE(qs_environment_type), POINTER :: qs_env 562 TYPE(ls_scf_env_type) :: ls_scf_env 563 INTEGER, INTENT(IN) :: unit_nr, iscf 564 LOGICAL, INTENT(OUT) :: diis_step 565 REAL(KIND=dp), INTENT(IN) :: eps_diis 566 INTEGER, INTENT(IN), OPTIONAL :: nmixing 567 TYPE(dbcsr_type), OPTIONAL :: s_matrix 568 REAL(KIND=dp), INTENT(IN) :: threshold 569 570 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_step_4lscf', & 571 routineP = moduleN//':'//routineN 572 REAL(KIND=dp), PARAMETER :: eigenvalue_threshold = 1.0E-12_dp 573 574 INTEGER :: handle, ib, ispin, jb, my_nmixing, nb, & 575 nb1, nspin 576 REAL(KIND=dp) :: error_max, tmp 577 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ev 578 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b 579 TYPE(cp_logger_type), POINTER :: logger 580 TYPE(cp_para_env_type), POINTER :: para_env 581 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks 582 TYPE(dbcsr_type) :: matrix_KSerr_t, matrix_tmp 583 TYPE(dbcsr_type), POINTER :: new_errors, old_errors, parameters 584 585 CALL timeset(routineN, handle) 586 nspin = ls_scf_env%nspins 587 diis_step = .FALSE. 588 my_nmixing = 2 589 IF (PRESENT(nmixing)) my_nmixing = nmixing 590 NULLIFY (new_errors, old_errors, parameters, a, b) 591 logger => cp_get_default_logger() 592 ! Quick return, if no DIIS is requested 593 IF (diis_buffer%nbuffer < 1) THEN 594 CALL timestop(handle) 595 RETURN 596 END IF 597 598! Getting current Kohn-Sham matrix from qs_env 599 CALL get_qs_env(qs_env, & 600 para_env=para_env, & 601 matrix_ks=matrix_ks) 602 CALL qs_diis_b_check_i_alloc_sparse( & 603 diis_buffer, & 604 ls_scf_env, & 605 nspin) 606 error_max = 0.0_dp 607 608 ib = MODULO(diis_buffer%ncall, diis_buffer%nbuffer) + 1 609 diis_buffer%ncall = diis_buffer%ncall + 1 610 nb = MIN(diis_buffer%ncall, diis_buffer%nbuffer) 611! Create scratch arrays 612 CALL dbcsr_create(matrix_tmp, & 613 template=ls_scf_env%matrix_ks(1), & 614 matrix_type='N') 615 CALL dbcsr_set(matrix_tmp, 0.0_dp) ! reset matrix 616 CALL dbcsr_create(matrix_KSerr_t, & 617 template=ls_scf_env%matrix_ks(1), & 618 matrix_type='N') 619 CALL dbcsr_set(matrix_KSerr_t, 0.0_dp) ! reset matrix 620 621 DO ispin = 1, nspin ! ------ Loop-ispin----START 622 623 new_errors => diis_buffer%error(ib, ispin)%matrix 624 parameters => diis_buffer%parameter(ib, ispin)%matrix 625 ! Copy the Kohn-Sham matrix K to the DIIS buffer 626 CALL dbcsr_copy(parameters, & ! out 627 matrix_ks(ispin)%matrix) ! in 628 629 IF (PRESENT(s_matrix)) THEN ! if-s_matrix ---------- START 630! Calculate Kohn-Sham error (non-orthogonal)= K*P*S-(K*P*S)^T 631! matrix_tmp = P*S 632 CALL dbcsr_multiply("N", "N", & 633 1.0_dp, ls_scf_env%matrix_p(ispin), & 634 s_matrix, & 635 0.0_dp, matrix_tmp, & 636 filter_eps=threshold) 637! new_errors= K*P*S 638 CALL dbcsr_multiply("N", "N", & 639 1.0_dp, matrix_ks(ispin)%matrix, & 640 matrix_tmp, & 641 0.0_dp, new_errors, & 642 filter_eps=threshold) 643! matrix_KSerr_t= transpose(K*P*S) 644 CALL dbcsr_transposed(matrix_KSerr_t, & 645 new_errors) 646! new_errors=K*P*S-transpose(K*P*S) 647 CALL dbcsr_add(new_errors, & 648 matrix_KSerr_t, & 649 1.0_dp, -1.0_dp) 650 ELSE ! if-s_matrix ---------- MID 651! Calculate Kohn-Sham error (orthogonal)= K*P - P*K 652! new_errors=K*P 653 CALL dbcsr_multiply("N", "N", & 654 1.0_dp, matrix_ks(ispin)%matrix, & 655 ls_scf_env%matrix_p(ispin), & 656 0.0_dp, new_errors, & 657 filter_eps=threshold) 658! matrix_KSerr_t= transpose(K*P) 659 CALL dbcsr_transposed(matrix_KSerr_t, & 660 new_errors) 661! new_errors=K*P-transpose(K*P) 662 CALL dbcsr_add(new_errors, & 663 matrix_KSerr_t, & 664 1.0_dp, -1.0_dp) 665 END IF ! if-s_matrix ---------- END 666 667 tmp = dbcsr_maxabs(new_errors) 668 error_max = MAX(error_max, tmp) 669 670 END DO ! ------ Loop-ispin----END 671 672 ! Check, if a DIIS step is appropiate 673 674 diis_step = (diis_buffer%ncall >= my_nmixing) 675 676 IF (unit_nr > 0) THEN 677 WRITE (unit_nr, '(A29,I3,A3,4(I3,A1))') & 678 "DIIS: (ncall,nbuffer,ib,nb)=(", iscf, ")=(", & 679 diis_buffer%ncall, ",", diis_buffer%nbuffer, ",", ib, ",", nb, ")" 680 WRITE (unit_nr, '(A57,I3,A3,L1,A1,F10.8,A1,F4.2,A1,L1,A1)') & 681 "DIIS: (diis_step,error_max,eps_diis,error_max<eps_diis)=(", & 682 iscf, ")=(", diis_step, ",", error_max, ",", eps_diis, ",", & 683 (error_max < eps_diis), ")" 684 WRITE (unit_nr, '(A75)') & 685 "DIIS: diis_step=T : Perform DIIS error_max<eps_diis=T : Update DIIS buffer" 686 ENDIF 687 688 ! Update the SCF DIIS buffer 689 IF (error_max < eps_diis) THEN 690 b => diis_buffer%b_matrix 691 DO jb = 1, nb 692 b(jb, ib) = 0.0_dp 693 DO ispin = 1, nspin 694 old_errors => diis_buffer%error(jb, ispin)%matrix 695 new_errors => diis_buffer%error(ib, ispin)%matrix 696 CALL dbcsr_dot(old_errors, & 697 new_errors, & 698 tmp) ! out : < f_i | f_j > 699 b(jb, ib) = b(jb, ib) + tmp 700 END DO ! end-loop-ispin 701 b(ib, jb) = b(jb, ib) 702 END DO ! end-loop-jb 703 ELSE 704 diis_step = .FALSE. 705 END IF 706 707 ! Perform DIIS step 708 IF (diis_step) THEN 709 nb1 = nb + 1 710 ALLOCATE (a(nb1, nb1)) 711 ALLOCATE (b(nb1, nb1)) 712 ALLOCATE (ev(nb1)) 713 ! Set up the linear DIIS equation system 714 b(1:nb, 1:nb) = diis_buffer%b_matrix(1:nb, 1:nb) 715 b(1:nb, nb1) = -1.0_dp 716 b(nb1, 1:nb) = -1.0_dp 717 b(nb1, nb1) = 0.0_dp 718 ! Solve the linear DIIS equation system 719 CALL diamat_all(b(1:nb1, 1:nb1), ev(1:nb1)) 720 a(1:nb1, 1:nb1) = b(1:nb1, 1:nb1) 721 DO jb = 1, nb1 722 IF (ABS(ev(jb)) < eigenvalue_threshold) THEN 723 a(1:nb1, jb) = 0.0_dp 724 ELSE 725 a(1:nb1, jb) = a(1:nb1, jb)/ev(jb) 726 END IF 727 END DO ! end-loop-jb 728 729 ev(1:nb) = MATMUL(a(1:nb, 1:nb1), b(nb1, 1:nb1)) 730 731 ! Update Kohn-Sham matrix 732 IF (iscf .GE. ls_scf_env%iter_ini_diis) THEN ! if-iscf-to-updateKS------ START 733 734 IF (unit_nr > 0) THEN 735 WRITE (unit_nr, '(A40,I3)') 'DIIS: Updating Kohn-Sham matrix at iscf=', iscf 736 ENDIF 737 738 DO ispin = 1, nspin 739 CALL dbcsr_set(matrix_ks(ispin)%matrix, & ! reset matrix 740 0.0_dp) 741 DO jb = 1, nb 742 parameters => diis_buffer%parameter(jb, ispin)%matrix 743 CALL dbcsr_add(matrix_ks(ispin)%matrix, parameters, & 744 1.0_dp, -ev(jb)) 745 END DO ! end-loop-jb 746 END DO ! end-loop-ispin 747 ENDIF ! if-iscf-to-updateKS------ END 748 749 DEALLOCATE (a) 750 DEALLOCATE (b) 751 DEALLOCATE (ev) 752 753 ELSE 754 DO ispin = 1, nspin 755 parameters => diis_buffer%parameter(ib, ispin)%matrix 756 CALL dbcsr_copy(parameters, & ! out 757 matrix_ks(ispin)%matrix) ! in 758 ENDDO ! end-loop-ispin 759 END IF 760 CALL dbcsr_release(matrix_tmp) 761 CALL dbcsr_release(matrix_KSerr_t) 762 CALL timestop(handle) 763 764 END SUBROUTINE qs_diis_b_step_4lscf 765 766! ************************************************************************************************** 767!> \brief Allocate and initialize a DIIS buffer with a buffer size of nbuffer. 768!> \param diis_buffer the buffer to initialize 769!> \param ls_scf_env ... 770!> \param nspin ... 771!> \par History 772!> - Adapted from qs_diis_b_check_i_alloc for sparse matrices and 773!> used in LS-SCF module (ls_scf_main) (10-11-14) 774!> \author Fredy W. Aquino 775!> \note 776!> check to allocate matrices only when needed 777! ************************************************************************************************** 778 779 SUBROUTINE qs_diis_b_check_i_alloc_sparse(diis_buffer, ls_scf_env, & 780 nspin) 781 782 TYPE(qs_diis_buffer_type_sparse), POINTER :: diis_buffer 783 TYPE(ls_scf_env_type) :: ls_scf_env 784 INTEGER, INTENT(IN) :: nspin 785 786 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_check_i_alloc_sparse', & 787 routineP = moduleN//':'//routineN 788 789 INTEGER :: handle, ibuffer, ispin, nbuffer 790 TYPE(cp_logger_type), POINTER :: logger 791 792! ------------------------------------------------------------------------- 793 794 CALL timeset(routineN, handle) 795 796 logger => cp_get_default_logger() 797 798 CPASSERT(ASSOCIATED(diis_buffer)) 799 CPASSERT(diis_buffer%ref_count > 0) 800 801 nbuffer = diis_buffer%nbuffer 802 803 IF (.NOT. ASSOCIATED(diis_buffer%error)) THEN 804 ALLOCATE (diis_buffer%error(nbuffer, nspin)) 805 806 DO ispin = 1, nspin 807 DO ibuffer = 1, nbuffer 808 ALLOCATE (diis_buffer%error(ibuffer, ispin)%matrix) 809 810 CALL dbcsr_create(diis_buffer%error(ibuffer, ispin)%matrix, & 811 template=ls_scf_env%matrix_ks(1), & 812 matrix_type='N') 813 END DO 814 END DO 815 END IF 816 817 IF (.NOT. ASSOCIATED(diis_buffer%parameter)) THEN 818 ALLOCATE (diis_buffer%parameter(nbuffer, nspin)) 819 820 DO ispin = 1, nspin 821 DO ibuffer = 1, nbuffer 822 ALLOCATE (diis_buffer%parameter(ibuffer, ispin)%matrix) 823 CALL dbcsr_create(diis_buffer%parameter(ibuffer, ispin)%matrix, & 824 template=ls_scf_env%matrix_ks(1), & 825 matrix_type='N') 826 END DO 827 END DO 828 END IF 829 830 IF (.NOT. ASSOCIATED(diis_buffer%b_matrix)) THEN 831 ALLOCATE (diis_buffer%b_matrix(nbuffer + 1, nbuffer + 1)) 832 833 diis_buffer%b_matrix = 0.0_dp 834 END IF 835 836 CALL timestop(handle) 837 838 END SUBROUTINE qs_diis_b_check_i_alloc_sparse 839 840! ************************************************************************************************** 841!> \brief clears the DIIS buffer in LS-SCF calculation 842!> \param diis_buffer the buffer to clear 843!> \par History 844!> 10-11-14 created [FA] modified from qs_diis_b_clear 845!> \author Fredy W. Aquino 846! ************************************************************************************************** 847 848 SUBROUTINE qs_diis_b_clear_sparse(diis_buffer) 849 850 TYPE(qs_diis_buffer_type_sparse), POINTER :: diis_buffer 851 852 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_clear_sparse', & 853 routineP = moduleN//':'//routineN 854 855 INTEGER :: handle 856 857! ------------------------------------------------------------------------- 858 859 CALL timeset(routineN, handle) 860 861 CPASSERT(ASSOCIATED(diis_buffer)) 862 CPASSERT(diis_buffer%ref_count > 0) 863 864 diis_buffer%ncall = 0 865 866 CALL timestop(handle) 867 868 END SUBROUTINE qs_diis_b_clear_sparse 869 870! ************************************************************************************************** 871!> \brief Allocates an SCF DIIS buffer for LS-SCF calculation 872!> \param diis_buffer the buffer to create 873!> \param nbuffer ... 874!> \par History 875!> 10-11-14 created [FA] modified from qs_diis_b_create 876!> \author Fredy W. Aquino 877! ************************************************************************************************** 878 SUBROUTINE qs_diis_b_create_sparse(diis_buffer, nbuffer) 879 880 TYPE(qs_diis_buffer_type_sparse), POINTER :: diis_buffer 881 INTEGER, INTENT(in) :: nbuffer 882 883 CHARACTER(len=*), PARAMETER :: routineN = 'qs_diis_b_create_sparse', & 884 routineP = moduleN//':'//routineN 885 886 INTEGER :: handle 887 888! ------------------------------------------------------------------------- 889 890 CALL timeset(routineN, handle) 891 892 ALLOCATE (diis_buffer) 893 894 NULLIFY (diis_buffer%b_matrix) 895 NULLIFY (diis_buffer%error) 896 NULLIFY (diis_buffer%parameter) 897 diis_buffer%nbuffer = nbuffer 898 diis_buffer%ncall = 0 899 last_diis_b_id = last_diis_b_id + 1 900 diis_buffer%id_nr = last_diis_b_id 901 diis_buffer%ref_count = 1 902 903 CALL timestop(handle) 904 905 END SUBROUTINE qs_diis_b_create_sparse 906 907END MODULE qs_diis 908