1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief computes preconditioners, and implements methods to apply them 8!> currently used in qs_ot 9!> \par History 10!> - [UB] 2009-05-13 Adding stable approximate inverse (full and sparse) 11!> \author Joost VandeVondele (09.2002) 12! ************************************************************************************************** 13MODULE preconditioner_makes 14 USE arnoldi_api, ONLY: arnoldi_data_type,& 15 arnoldi_ev,& 16 deallocate_arnoldi_data,& 17 get_selected_ritz_val,& 18 get_selected_ritz_vector,& 19 set_arnoldi_initial_vector,& 20 setup_arnoldi_data 21 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 22 cp_dbcsr_m_by_n_from_template,& 23 cp_dbcsr_sm_fm_multiply,& 24 cp_fm_to_dbcsr_row_template 25 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,& 26 cp_fm_triangular_invert,& 27 cp_fm_triangular_multiply,& 28 cp_fm_upper_to_full 29 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,& 30 cp_fm_cholesky_reduce,& 31 cp_fm_cholesky_restore 32 USE cp_fm_diag, ONLY: choose_eigv_solver 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_diag,& 38 cp_fm_get_info,& 39 cp_fm_release,& 40 cp_fm_to_fm,& 41 cp_fm_type 42 USE cp_gemm_interface, ONLY: cp_gemm 43 USE dbcsr_api, ONLY: & 44 dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_get_info, dbcsr_multiply, & 45 dbcsr_p_type, dbcsr_release, dbcsr_type, dbcsr_type_symmetric 46 USE input_constants, ONLY: & 47 cholesky_inverse, cholesky_reduce, ot_precond_full_all, ot_precond_full_kinetic, & 48 ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_s_inverse, & 49 ot_precond_solver_default, ot_precond_solver_inv_chol 50 USE kinds, ONLY: dp 51 USE preconditioner_types, ONLY: preconditioner_type 52#include "./base/base_uses.f90" 53 54 IMPLICIT NONE 55 56 PRIVATE 57 58 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner_makes' 59 60 PUBLIC :: make_preconditioner_matrix 61 62CONTAINS 63 64! ************************************************************************************************** 65!> \brief ... 66!> \param preconditioner_env ... 67!> \param matrix_h ... 68!> \param matrix_s ... 69!> \param matrix_t ... 70!> \param mo_coeff ... 71!> \param energy_homo ... 72!> \param eigenvalues_ot ... 73!> \param energy_gap ... 74!> \param my_solver_type ... 75! ************************************************************************************************** 76 SUBROUTINE make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, & 77 energy_homo, eigenvalues_ot, energy_gap, & 78 my_solver_type) 79 TYPE(preconditioner_type) :: preconditioner_env 80 TYPE(dbcsr_type), POINTER :: matrix_h 81 TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s, matrix_t 82 TYPE(cp_fm_type), POINTER :: mo_coeff 83 REAL(KIND=dp) :: energy_homo 84 REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues_ot 85 REAL(KIND=dp) :: energy_gap 86 INTEGER :: my_solver_type 87 88 CHARACTER(len=*), PARAMETER :: routineN = 'make_preconditioner_matrix', & 89 routineP = moduleN//':'//routineN 90 91 INTEGER :: precon_type 92 93 precon_type = preconditioner_env%in_use 94 SELECT CASE (precon_type) 95 CASE (ot_precond_full_single) 96 IF (my_solver_type .NE. ot_precond_solver_default) & 97 CPABORT("Only PRECOND_SOLVER DEFAULT for the moment") 98 IF (PRESENT(matrix_s)) THEN 99 CALL make_full_single(preconditioner_env, preconditioner_env%fm, & 100 matrix_h, matrix_s, energy_homo, energy_gap) 101 ELSE 102 CALL make_full_single_ortho(preconditioner_env, preconditioner_env%fm, & 103 matrix_h, energy_homo, energy_gap) 104 END IF 105 106 CASE (ot_precond_s_inverse) 107 IF (my_solver_type .EQ. ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol 108 IF (.NOT. PRESENT(matrix_s)) & 109 CPABORT("Type for S=1 not implemented") 110 CALL make_full_s_inverse(preconditioner_env, matrix_s) 111 112 CASE (ot_precond_full_kinetic) 113 IF (my_solver_type .EQ. ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol 114 IF (.NOT. (PRESENT(matrix_s) .AND. PRESENT(matrix_t))) & 115 CPABORT("Type for S=1 not implemented") 116 CALL make_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap) 117 CASE (ot_precond_full_single_inverse) 118 IF (my_solver_type .EQ. ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol 119 CALL make_full_single_inverse(preconditioner_env, mo_coeff, matrix_h, energy_gap, & 120 matrix_s=matrix_s) 121 CASE (ot_precond_full_all) 122 IF (my_solver_type .NE. ot_precond_solver_default) THEN 123 CPABORT("Only PRECOND_SOLVER DEFAULT for the moment") 124 ENDIF 125 IF (PRESENT(matrix_s)) THEN 126 CALL make_full_all(preconditioner_env, mo_coeff, matrix_h, matrix_s, & 127 eigenvalues_ot, energy_gap) 128 ELSE 129 CALL make_full_all_ortho(preconditioner_env, mo_coeff, matrix_h, & 130 eigenvalues_ot, energy_gap) 131 END IF 132 133 CASE DEFAULT 134 CPABORT("Type not implemented") 135 END SELECT 136 137 END SUBROUTINE make_preconditioner_matrix 138 139! ************************************************************************************************** 140!> \brief Simply takes the overlap matrix as preconditioner 141!> \param preconditioner_env ... 142!> \param matrix_s ... 143! ************************************************************************************************** 144 SUBROUTINE make_full_s_inverse(preconditioner_env, matrix_s) 145 TYPE(preconditioner_type) :: preconditioner_env 146 TYPE(dbcsr_type), POINTER :: matrix_s 147 148 CHARACTER(len=*), PARAMETER :: routineN = 'make_full_s_inverse', & 149 routineP = moduleN//':'//routineN 150 151 INTEGER :: handle 152 153 CALL timeset(routineN, handle) 154 155 CPASSERT(ASSOCIATED(matrix_s)) 156 157 IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN 158 ALLOCATE (preconditioner_env%sparse_matrix) 159 END IF 160 CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_s, name="full_kinetic") 161 162 CALL timestop(handle) 163 164 END SUBROUTINE make_full_s_inverse 165 166! ************************************************************************************************** 167!> \brief kinetic matrix+shift*overlap as preconditioner. Cheap but could 168!> be better 169!> \param preconditioner_env ... 170!> \param matrix_t ... 171!> \param matrix_s ... 172!> \param energy_gap ... 173! ************************************************************************************************** 174 SUBROUTINE make_full_kinetic(preconditioner_env, matrix_t, matrix_s, & 175 energy_gap) 176 TYPE(preconditioner_type) :: preconditioner_env 177 TYPE(dbcsr_type), POINTER :: matrix_t, matrix_s 178 REAL(KIND=dp) :: energy_gap 179 180 CHARACTER(len=*), PARAMETER :: routineN = 'make_full_kinetic', & 181 routineP = moduleN//':'//routineN 182 183 INTEGER :: handle 184 REAL(KIND=dp) :: shift 185 186 CALL timeset(routineN, handle) 187 188 CPASSERT(ASSOCIATED(matrix_t)) 189 CPASSERT(ASSOCIATED(matrix_s)) 190 191 IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN 192 ALLOCATE (preconditioner_env%sparse_matrix) 193 END IF 194 CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_t, name="full_kinetic") 195 196 shift = MAX(0.0_dp, energy_gap) 197 198 CALL dbcsr_add(preconditioner_env%sparse_matrix, matrix_s, & 199 alpha_scalar=1.0_dp, beta_scalar=shift) 200 201 CALL timestop(handle) 202 203 END SUBROUTINE make_full_kinetic 204 205! ************************************************************************************************** 206!> \brief full_single_preconditioner 207!> \param preconditioner_env ... 208!> \param fm ... 209!> \param matrix_h ... 210!> \param matrix_s ... 211!> \param energy_homo ... 212!> \param energy_gap ... 213! ************************************************************************************************** 214 SUBROUTINE make_full_single(preconditioner_env, fm, matrix_h, matrix_s, & 215 energy_homo, energy_gap) 216 TYPE(preconditioner_type) :: preconditioner_env 217 TYPE(cp_fm_type), POINTER :: fm 218 TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s 219 REAL(KIND=dp) :: energy_homo, energy_gap 220 221 CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single', & 222 routineP = moduleN//':'//routineN 223 224 INTEGER :: handle, i, n 225 REAL(KIND=dp), DIMENSION(:), POINTER :: evals 226 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 227 TYPE(cp_fm_type), POINTER :: fm_h, fm_s 228 229 CALL timeset(routineN, handle) 230 231 NULLIFY (fm_h, fm_s, fm_struct_tmp, evals) 232 233 IF (ASSOCIATED(fm)) THEN 234 CALL cp_fm_release(fm) 235 ENDIF 236 CALL dbcsr_get_info(matrix_h, nfullrows_total=n) 237 ALLOCATE (evals(n)) 238 239 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, & 240 context=preconditioner_env%ctxt, & 241 para_env=preconditioner_env%para_env) 242 CALL cp_fm_create(fm, fm_struct_tmp, name="preconditioner") 243 CALL cp_fm_create(fm_h, fm_struct_tmp, name="fm_h") 244 CALL cp_fm_create(fm_s, fm_struct_tmp, name="fm_s") 245 CALL cp_fm_struct_release(fm_struct_tmp) 246 247 CALL copy_dbcsr_to_fm(matrix_h, fm_h) 248 CALL copy_dbcsr_to_fm(matrix_s, fm_s) 249 CALL cp_fm_cholesky_decompose(fm_s) 250 251 SELECT CASE (preconditioner_env%cholesky_use) 252 CASE (cholesky_inverse) 253! if cho inverse 254 CALL cp_fm_triangular_invert(fm_s) 255 CALL cp_fm_upper_to_full(fm_h, fm) 256 257 CALL cp_fm_triangular_multiply(fm_s, fm_h, side="R", transpose_tr=.FALSE., & 258 invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp) 259 CALL cp_fm_triangular_multiply(fm_s, fm_h, side="L", transpose_tr=.TRUE., & 260 invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp) 261 CASE (cholesky_reduce) 262 CALL cp_fm_cholesky_reduce(fm_h, fm_s) 263 CASE DEFAULT 264 CPABORT("cholesky type not implemented") 265 END SELECT 266 267 CALL choose_eigv_solver(fm_h, fm, evals) 268 269 SELECT CASE (preconditioner_env%cholesky_use) 270 CASE (cholesky_inverse) 271 CALL cp_fm_triangular_multiply(fm_s, fm, side="L", transpose_tr=.FALSE., & 272 invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp) 273 DO i = 1, n 274 evals(i) = 1.0_dp/MAX(evals(i) - energy_homo, energy_gap) 275 ENDDO 276 CALL cp_fm_to_fm(fm, fm_h) 277 CASE (cholesky_reduce) 278 CALL cp_fm_cholesky_restore(fm, n, fm_s, fm_h, "SOLVE") 279 DO i = 1, n 280 evals(i) = 1.0_dp/MAX(evals(i) - energy_homo, energy_gap) 281 ENDDO 282 CALL cp_fm_to_fm(fm_h, fm) 283 END SELECT 284 285 CALL cp_fm_column_scale(fm, evals) 286 CALL cp_gemm('N', 'T', n, n, n, 1.0_dp, fm, fm_h, 0.0_dp, fm_s) 287 CALL cp_fm_to_fm(fm_s, fm) 288 289 DEALLOCATE (evals) 290 CALL cp_fm_release(fm_h) 291 CALL cp_fm_release(fm_s) 292 293 CALL timestop(handle) 294 295 END SUBROUTINE make_full_single 296 297! ************************************************************************************************** 298!> \brief full single in the orthonormal basis 299!> \param preconditioner_env ... 300!> \param fm ... 301!> \param matrix_h ... 302!> \param energy_homo ... 303!> \param energy_gap ... 304! ************************************************************************************************** 305 SUBROUTINE make_full_single_ortho(preconditioner_env, fm, matrix_h, & 306 energy_homo, energy_gap) 307 TYPE(preconditioner_type) :: preconditioner_env 308 TYPE(cp_fm_type), POINTER :: fm 309 TYPE(dbcsr_type), POINTER :: matrix_h 310 REAL(KIND=dp) :: energy_homo, energy_gap 311 312 CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single_ortho', & 313 routineP = moduleN//':'//routineN 314 315 INTEGER :: handle, i, n 316 REAL(KIND=dp), DIMENSION(:), POINTER :: evals 317 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 318 TYPE(cp_fm_type), POINTER :: fm_h, fm_s 319 320 CALL timeset(routineN, handle) 321 NULLIFY (fm_h, fm_s, fm_struct_tmp, evals) 322 323 IF (ASSOCIATED(fm)) THEN 324 CALL cp_fm_release(fm) 325 ENDIF 326 CALL dbcsr_get_info(matrix_h, nfullrows_total=n) 327 ALLOCATE (evals(n)) 328 329 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, & 330 context=preconditioner_env%ctxt, & 331 para_env=preconditioner_env%para_env) 332 CALL cp_fm_create(fm, fm_struct_tmp, name="preconditioner") 333 CALL cp_fm_create(fm_h, fm_struct_tmp, name="fm_h") 334 CALL cp_fm_create(fm_s, fm_struct_tmp, name="fm_s") 335 CALL cp_fm_struct_release(fm_struct_tmp) 336 337 CALL copy_dbcsr_to_fm(matrix_h, fm_h) 338 339 CALL choose_eigv_solver(fm_h, fm, evals) 340 DO i = 1, n 341 evals(i) = 1.0_dp/MAX(evals(i) - energy_homo, energy_gap) 342 ENDDO 343 CALL cp_fm_to_fm(fm, fm_h) 344 CALL cp_fm_column_scale(fm, evals) 345 CALL cp_gemm('N', 'T', n, n, n, 1.0_dp, fm, fm_h, 0.0_dp, fm_s) 346 CALL cp_fm_to_fm(fm_s, fm) 347 348 DEALLOCATE (evals) 349 CALL cp_fm_release(fm_h) 350 CALL cp_fm_release(fm_s) 351 352 CALL timestop(handle) 353 354 END SUBROUTINE make_full_single_ortho 355 356! ************************************************************************************************** 357!> \brief generates a state by state preconditioner based on the full hamiltonian matrix 358!> \param preconditioner_env ... 359!> \param matrix_c0 ... 360!> \param matrix_h ... 361!> \param matrix_s ... 362!> \param c0_evals ... 363!> \param energy_gap should be a slight underestimate of the physical energy gap for almost all systems 364!> the c0 are already ritz states of (h,s) 365!> \par History 366!> 10.2006 made more stable [Joost VandeVondele] 367!> \note 368!> includes error estimate on the hamiltonian matrix to result in a stable preconditioner 369!> a preconditioner for each eigenstate i is generated by keeping the factorized form 370!> U diag( something i ) U^T. It is important to only precondition in the subspace orthogonal to c0. 371!> not only is it the only part that matters, it also simplifies the computation of 372!> the lagrangian multipliers in the OT minimization (i.e. if the c0 here is different 373!> from the c0 used in the OT setup, there will be a bug). 374! ************************************************************************************************** 375 SUBROUTINE make_full_all(preconditioner_env, matrix_c0, matrix_h, matrix_s, c0_evals, energy_gap) 376 TYPE(preconditioner_type) :: preconditioner_env 377 TYPE(cp_fm_type), POINTER :: matrix_c0 378 TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s 379 REAL(KIND=dp), DIMENSION(:), POINTER :: c0_evals 380 REAL(KIND=dp) :: energy_gap 381 382 CHARACTER(len=*), PARAMETER :: routineN = 'make_full_all', routineP = moduleN//':'//routineN 383 REAL(KIND=dp), PARAMETER :: fudge_factor = 0.25_dp, & 384 lambda_base = 10.0_dp 385 386 INTEGER :: handle, k, n 387 REAL(KIND=dp) :: error_estimate, lambda 388 REAL(KIND=dp), DIMENSION(:), POINTER :: diag, norms, shifted_evals 389 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 390 TYPE(cp_fm_type), POINTER :: matrix_hc0, matrix_left, matrix_pre, & 391 matrix_s1, matrix_s2, matrix_sc0, & 392 matrix_shc0, matrix_tmp, ortho 393 394 CALL timeset(routineN, handle) 395 396 IF (ASSOCIATED(preconditioner_env%fm)) CALL cp_fm_release(preconditioner_env%fm) 397 CALL cp_fm_get_info(matrix_c0, nrow_global=n, ncol_global=k) 398 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, & 399 context=preconditioner_env%ctxt, & 400 para_env=preconditioner_env%para_env) 401 CALL cp_fm_create(preconditioner_env%fm, fm_struct_tmp, name="preconditioner_env%fm") 402 matrix_pre => preconditioner_env%fm 403 CALL cp_fm_create(ortho, fm_struct_tmp, name="ortho") 404 CALL cp_fm_create(matrix_tmp, fm_struct_tmp, name="matrix_tmp") 405 CALL cp_fm_struct_release(fm_struct_tmp) 406 ALLOCATE (preconditioner_env%full_evals(n)) 407 ALLOCATE (preconditioner_env%occ_evals(k)) 408 409 ! 0) cholesky decompose the overlap matrix, if this fails the basis is singular, 410 ! more than EPS_DEFAULT 411 CALL copy_dbcsr_to_fm(matrix_s, ortho) 412 CALL cp_fm_cholesky_decompose(ortho) 413! if cho inverse 414 IF (preconditioner_env%cholesky_use == cholesky_inverse) THEN 415 CALL cp_fm_triangular_invert(ortho) 416 END IF 417 ! 1) Construct a new H matrix, which has the current C0 as eigenvectors, 418 ! possibly shifted by an amount lambda, 419 ! and the same spectrum as the original H matrix in the space orthogonal to the C0 420 ! with P=C0 C0 ^ T 421 ! (1 - PS)^T H (1-PS) + (PS)^T (H - lambda S ) (PS) 422 ! we exploit that the C0 are already the ritz states of H 423 CALL cp_fm_create(matrix_sc0, matrix_c0%matrix_struct, name="sc0") 424 CALL cp_dbcsr_sm_fm_multiply(matrix_s, matrix_c0, matrix_sc0, k) 425 CALL cp_fm_create(matrix_hc0, matrix_c0%matrix_struct, name="hc0") 426 CALL cp_dbcsr_sm_fm_multiply(matrix_h, matrix_c0, matrix_hc0, k) 427 428 ! An aside, try to estimate the error on the ritz values, we'll need it later on 429 CALL cp_fm_create(matrix_shc0, matrix_c0%matrix_struct, name="shc0") 430 431 SELECT CASE (preconditioner_env%cholesky_use) 432 CASE (cholesky_inverse) 433! if cho inverse 434 CALL cp_fm_to_fm(matrix_hc0, matrix_shc0) 435 CALL cp_fm_triangular_multiply(ortho, matrix_shc0, side="L", transpose_tr=.TRUE., & 436 invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=k, alpha=1.0_dp) 437 CASE (cholesky_reduce) 438 CALL cp_fm_cholesky_restore(matrix_hc0, k, ortho, matrix_shc0, "SOLVE", transa="T") 439 CASE DEFAULT 440 CPABORT("cholesky type not implemented") 441 END SELECT 442 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, & 443 context=preconditioner_env%ctxt, & 444 para_env=preconditioner_env%para_env) 445 CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1") 446 CALL cp_fm_struct_release(fm_struct_tmp) 447 ! since we only use diagonal elements this is a bit of a waste 448 CALL cp_gemm('T', 'N', k, k, n, 1.0_dp, matrix_shc0, matrix_shc0, 0.0_dp, matrix_s1) 449 ALLOCATE (diag(k)) 450 CALL cp_fm_get_diag(matrix_s1, diag) 451 error_estimate = MAXVAL(SQRT(ABS(diag - c0_evals**2))) 452 DEALLOCATE (diag) 453 CALL cp_fm_release(matrix_s1) 454 CALL cp_fm_release(matrix_shc0) 455 ! we'll only use the energy gap, if our estimate of the error on the eigenvalues 456 ! is small enough. A large error combined with a small energy gap would otherwise lead to 457 ! an aggressive but bad preconditioner. Only when the error is small (MD), we can precondition 458 ! aggressively 459 preconditioner_env%energy_gap = MAX(energy_gap, error_estimate*fudge_factor) 460 CALL copy_dbcsr_to_fm(matrix_h, matrix_tmp) 461 CALL cp_fm_upper_to_full(matrix_tmp, matrix_pre) 462 ! tmp = H ( 1 - PS ) 463 CALL cp_gemm('N', 'T', n, n, k, -1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp) 464 465 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=n, & 466 context=preconditioner_env%ctxt, & 467 para_env=preconditioner_env%para_env) 468 CALL cp_fm_create(matrix_left, fm_struct_tmp, name="matrix_left") 469 CALL cp_fm_struct_release(fm_struct_tmp) 470 CALL cp_gemm('T', 'N', k, n, n, 1.0_dp, matrix_c0, matrix_tmp, 0.0_dp, matrix_left) 471 ! tmp = (1 - PS)^T H (1-PS) 472 CALL cp_gemm('N', 'N', n, n, k, -1.0_dp, matrix_sc0, matrix_left, 1.0_dp, matrix_tmp) 473 CALL cp_fm_release(matrix_left) 474 475 ALLOCATE (shifted_evals(k)) 476 lambda = lambda_base + error_estimate 477 shifted_evals = c0_evals - lambda 478 CALL cp_fm_to_fm(matrix_sc0, matrix_hc0) 479 CALL cp_fm_column_scale(matrix_hc0, shifted_evals) 480 CALL cp_gemm('N', 'T', n, n, k, 1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp) 481 482 ! 2) diagonalize this operator 483 SELECT CASE (preconditioner_env%cholesky_use) 484 CASE (cholesky_inverse) 485 CALL cp_fm_triangular_multiply(ortho, matrix_tmp, side="R", transpose_tr=.FALSE., & 486 invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp) 487 CALL cp_fm_triangular_multiply(ortho, matrix_tmp, side="L", transpose_tr=.TRUE., & 488 invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp) 489 CASE (cholesky_reduce) 490 CALL cp_fm_cholesky_reduce(matrix_tmp, ortho) 491 END SELECT 492 CALL choose_eigv_solver(matrix_tmp, matrix_pre, preconditioner_env%full_evals) 493 SELECT CASE (preconditioner_env%cholesky_use) 494 CASE (cholesky_inverse) 495 CALL cp_fm_triangular_multiply(ortho, matrix_pre, side="L", transpose_tr=.FALSE., & 496 invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp) 497 CALL cp_fm_to_fm(matrix_pre, matrix_tmp) 498 CASE (cholesky_reduce) 499 CALL cp_fm_cholesky_restore(matrix_pre, n, ortho, matrix_tmp, "SOLVE") 500 CALL cp_fm_to_fm(matrix_tmp, matrix_pre) 501 END SELECT 502 503 ! test that the subspace remained conserved 504 IF (.FALSE.) THEN 505 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, & 506 context=preconditioner_env%ctxt, & 507 para_env=preconditioner_env%para_env) 508 CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1") 509 CALL cp_fm_create(matrix_s2, fm_struct_tmp, name="matrix_s2") 510 CALL cp_fm_struct_release(fm_struct_tmp) 511 ALLOCATE (norms(k)) 512 CALL cp_gemm('T', 'N', k, k, n, 1.0_dp, matrix_sc0, matrix_tmp, 0.0_dp, matrix_s1) 513 CALL choose_eigv_solver(matrix_s1, matrix_s2, norms) 514 WRITE (*, *) "matrix norm deviation (should be close to zero): ", MAXVAL(ABS(ABS(norms) - 1.0_dp)) 515 DEALLOCATE (norms) 516 CALL cp_fm_release(matrix_s1) 517 CALL cp_fm_release(matrix_s2) 518 ENDIF 519 520 ! 3) replace the lowest k evals and evecs with what they should be 521 preconditioner_env%occ_evals = c0_evals 522 ! notice, this choice causes the preconditioner to be constant when applied to sc0 (see apply_full_all) 523 preconditioner_env%full_evals(1:k) = c0_evals 524 CALL cp_fm_to_fm(matrix_c0, matrix_pre, k, 1, 1) 525 526 CALL cp_fm_release(matrix_sc0) 527 CALL cp_fm_release(matrix_hc0) 528 CALL cp_fm_release(ortho) 529 CALL cp_fm_release(matrix_tmp) 530 DEALLOCATE (shifted_evals) 531 CALL timestop(handle) 532 533 END SUBROUTINE make_full_all 534 535! ************************************************************************************************** 536!> \brief full all in the orthonormal basis 537!> \param preconditioner_env ... 538!> \param matrix_c0 ... 539!> \param matrix_h ... 540!> \param c0_evals ... 541!> \param energy_gap ... 542! ************************************************************************************************** 543 SUBROUTINE make_full_all_ortho(preconditioner_env, matrix_c0, matrix_h, c0_evals, energy_gap) 544 545 TYPE(preconditioner_type) :: preconditioner_env 546 TYPE(cp_fm_type), POINTER :: matrix_c0 547 TYPE(dbcsr_type), POINTER :: matrix_h 548 REAL(KIND=dp), DIMENSION(:), POINTER :: c0_evals 549 REAL(KIND=dp) :: energy_gap 550 551 CHARACTER(len=*), PARAMETER :: routineN = 'make_full_all_ortho', & 552 routineP = moduleN//':'//routineN 553 REAL(KIND=dp), PARAMETER :: fudge_factor = 0.25_dp, & 554 lambda_base = 10.0_dp 555 556 INTEGER :: handle, k, n 557 REAL(KIND=dp) :: error_estimate, lambda 558 REAL(KIND=dp), DIMENSION(:), POINTER :: diag, norms, shifted_evals 559 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 560 TYPE(cp_fm_type), POINTER :: matrix_hc0, matrix_left, matrix_pre, & 561 matrix_s1, matrix_s2, matrix_sc0, & 562 matrix_tmp 563 564 CALL timeset(routineN, handle) 565 566 IF (ASSOCIATED(preconditioner_env%fm)) CALL cp_fm_release(preconditioner_env%fm) 567 CALL cp_fm_get_info(matrix_c0, nrow_global=n, ncol_global=k) 568 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, & 569 context=preconditioner_env%ctxt, & 570 para_env=preconditioner_env%para_env) 571 CALL cp_fm_create(preconditioner_env%fm, fm_struct_tmp, name="preconditioner_env%fm") 572 matrix_pre => preconditioner_env%fm 573 CALL cp_fm_create(matrix_tmp, fm_struct_tmp, name="matrix_tmp") 574 CALL cp_fm_struct_release(fm_struct_tmp) 575 ALLOCATE (preconditioner_env%full_evals(n)) 576 ALLOCATE (preconditioner_env%occ_evals(k)) 577 578 ! 1) Construct a new H matrix, which has the current C0 as eigenvectors, 579 ! possibly shifted by an amount lambda, 580 ! and the same spectrum as the original H matrix in the space orthogonal to the C0 581 ! with P=C0 C0 ^ T 582 ! (1 - PS)^T H (1-PS) + (PS)^T (H - lambda S ) (PS) 583 ! we exploit that the C0 are already the ritz states of H 584 CALL cp_fm_create(matrix_sc0, matrix_c0%matrix_struct, name="sc0") 585 CALL cp_fm_to_fm(matrix_c0, matrix_sc0) 586 CALL cp_fm_create(matrix_hc0, matrix_c0%matrix_struct, name="hc0") 587 CALL cp_dbcsr_sm_fm_multiply(matrix_h, matrix_c0, matrix_hc0, k) 588 589 ! An aside, try to estimate the error on the ritz values, we'll need it later on 590 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, & 591 context=preconditioner_env%ctxt, & 592 para_env=preconditioner_env%para_env) 593 CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1") 594 CALL cp_fm_struct_release(fm_struct_tmp) 595 ! since we only use diagonal elements this is a bit of a waste 596 CALL cp_gemm('T', 'N', k, k, n, 1.0_dp, matrix_hc0, matrix_hc0, 0.0_dp, matrix_s1) 597 ALLOCATE (diag(k)) 598 CALL cp_fm_get_diag(matrix_s1, diag) 599 error_estimate = MAXVAL(SQRT(ABS(diag - c0_evals**2))) 600 DEALLOCATE (diag) 601 CALL cp_fm_release(matrix_s1) 602 ! we'll only use the energy gap, if our estimate of the error on the eigenvalues 603 ! is small enough. A large error combined with a small energy gap would otherwise lead to 604 ! an aggressive but bad preconditioner. Only when the error is small (MD), we can precondition 605 ! aggressively 606 preconditioner_env%energy_gap = MAX(energy_gap, error_estimate*fudge_factor) 607 608 CALL copy_dbcsr_to_fm(matrix_h, matrix_tmp) 609 CALL cp_fm_upper_to_full(matrix_tmp, matrix_pre) 610 ! tmp = H ( 1 - PS ) 611 CALL cp_gemm('N', 'T', n, n, k, -1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp) 612 613 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=n, & 614 context=preconditioner_env%ctxt, & 615 para_env=preconditioner_env%para_env) 616 CALL cp_fm_create(matrix_left, fm_struct_tmp, name="matrix_left") 617 CALL cp_fm_struct_release(fm_struct_tmp) 618 CALL cp_gemm('T', 'N', k, n, n, 1.0_dp, matrix_c0, matrix_tmp, 0.0_dp, matrix_left) 619 ! tmp = (1 - PS)^T H (1-PS) 620 CALL cp_gemm('N', 'N', n, n, k, -1.0_dp, matrix_sc0, matrix_left, 1.0_dp, matrix_tmp) 621 CALL cp_fm_release(matrix_left) 622 623 ALLOCATE (shifted_evals(k)) 624 lambda = lambda_base + error_estimate 625 shifted_evals = c0_evals - lambda 626 CALL cp_fm_to_fm(matrix_sc0, matrix_hc0) 627 CALL cp_fm_column_scale(matrix_hc0, shifted_evals) 628 CALL cp_gemm('N', 'T', n, n, k, 1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp) 629 630 ! 2) diagonalize this operator 631 CALL choose_eigv_solver(matrix_tmp, matrix_pre, preconditioner_env%full_evals) 632 633 ! test that the subspace remained conserved 634 IF (.FALSE.) THEN 635 CALL cp_fm_to_fm(matrix_pre, matrix_tmp) 636 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, & 637 context=preconditioner_env%ctxt, & 638 para_env=preconditioner_env%para_env) 639 CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1") 640 CALL cp_fm_create(matrix_s2, fm_struct_tmp, name="matrix_s2") 641 CALL cp_fm_struct_release(fm_struct_tmp) 642 ALLOCATE (norms(k)) 643 CALL cp_gemm('T', 'N', k, k, n, 1.0_dp, matrix_sc0, matrix_tmp, 0.0_dp, matrix_s1) 644 CALL choose_eigv_solver(matrix_s1, matrix_s2, norms) 645 646 WRITE (*, *) "matrix norm deviation (should be close to zero): ", MAXVAL(ABS(ABS(norms) - 1.0_dp)) 647 DEALLOCATE (norms) 648 CALL cp_fm_release(matrix_s1) 649 CALL cp_fm_release(matrix_s2) 650 ENDIF 651 652 ! 3) replace the lowest k evals and evecs with what they should be 653 preconditioner_env%occ_evals = c0_evals 654 ! notice, this choice causes the preconditioner to be constant when applied to sc0 (see apply_full_all) 655 preconditioner_env%full_evals(1:k) = c0_evals 656 CALL cp_fm_to_fm(matrix_c0, matrix_pre, k, 1, 1) 657 658 CALL cp_fm_release(matrix_sc0) 659 CALL cp_fm_release(matrix_hc0) 660 CALL cp_fm_release(matrix_tmp) 661 DEALLOCATE (shifted_evals) 662 663 CALL timestop(handle) 664 665 END SUBROUTINE make_full_all_ortho 666 667! ************************************************************************************************** 668!> \brief generates a preconditioner matrix H-lambda S+(SC)(2.0*CT*H*C+delta)(SC)^T 669!> for later inversion. 670!> H is the Kohn Sham matrix 671!> lambda*S shifts the spectrum of the generalized form up by lambda 672!> the last term only shifts the occupied space (reversing them in energy order) 673!> This form is implicitely multiplied from both sides by S^0.5 674!> This ensures we precondition the correct quantity 675!> Before this reads S^-0.5 H S^-0.5 + lambda + (S^0.5 C)shifts(S^0.5 C)T 676!> which might be a bit more obvious 677!> Replaced the old full_single_inverse at revision 14616 678!> \param preconditioner_env the preconditioner env 679!> \param matrix_c0 the MO coefficient matrix (fm) 680!> \param matrix_h Kohn-Sham matrix (dbcsr) 681!> \param energy_gap an additional shift in lambda=-E_homo+energy_gap 682!> \param matrix_s the overlap matrix if not orthonormal (dbcsr, optional) 683! ************************************************************************************************** 684 SUBROUTINE make_full_single_inverse(preconditioner_env, matrix_c0, matrix_h, energy_gap, matrix_s) 685 TYPE(preconditioner_type) :: preconditioner_env 686 TYPE(cp_fm_type), POINTER :: matrix_c0 687 TYPE(dbcsr_type), POINTER :: matrix_h 688 REAL(KIND=dp) :: energy_gap 689 TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s 690 691 CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single_inverse', & 692 routineP = moduleN//':'//routineN 693 694 INTEGER :: handle, k, n 695 REAL(KIND=dp) :: max_ev, min_ev, pre_shift 696 TYPE(arnoldi_data_type) :: my_arnoldi 697 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrices 698 TYPE(dbcsr_type), TARGET :: dbcsr_cThc, dbcsr_hc, dbcsr_sc, mo_dbcsr 699 700 CALL timeset(routineN, handle) 701 702 ! Allocate all working matrices needed 703 CALL cp_fm_get_info(matrix_c0, nrow_global=n, ncol_global=k) 704 ! copy the fm MO's to a sparse matrix, can be solved better if the sparse version is already present 705 ! but for the time beeing this will do 706 CALL cp_fm_to_dbcsr_row_template(mo_dbcsr, matrix_c0, matrix_h) 707 CALL dbcsr_create(dbcsr_sc, template=mo_dbcsr) 708 CALL dbcsr_create(dbcsr_hc, template=mo_dbcsr) 709 CALL cp_dbcsr_m_by_n_from_template(dbcsr_cThc, matrix_h, k, k, sym=dbcsr_type_symmetric) 710 711 ! Check whether the output matrix was already created, if not do it now 712 IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN 713 ALLOCATE (preconditioner_env%sparse_matrix) 714 END IF 715 716 ! Put the first term of the preconditioner (H) into the output matrix 717 CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_h) 718 719 ! Precompute some matrices 720 ! S*C, if orthonormal this will be simply C so a copy will do 721 IF (PRESENT(matrix_s)) THEN 722 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s, mo_dbcsr, 0.0_dp, dbcsr_sc) 723 ELSE 724 CALL dbcsr_copy(dbcsr_sc, mo_dbcsr) 725 END IF 726 727!----------------------------compute the occupied subspace and shift it ------------------------------------ 728 ! cT*H*C which will be used to shift the occupied states to 0 729 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_h, mo_dbcsr, 0.0_dp, dbcsr_hc) 730 CALL dbcsr_multiply("T", "N", 1.0_dp, mo_dbcsr, dbcsr_hc, 0.0_dp, dbcsr_cThc) 731 732 ! Compute the Energy of the HOMO. We will use this as a reference energy 733 ALLOCATE (matrices(1)) 734 matrices(1)%matrix => dbcsr_cThc 735 CALL setup_arnoldi_data(my_arnoldi, matrices, max_iter=20, threshold=1.0E-3_dp, selection_crit=2, & 736 nval_request=1, nrestarts=8, generalized_ev=.FALSE., iram=.FALSE.) 737 IF (ASSOCIATED(preconditioner_env%max_ev_vector)) & 738 CALL set_arnoldi_initial_vector(my_arnoldi, preconditioner_env%max_ev_vector) 739 CALL arnoldi_ev(matrices, my_arnoldi) 740 max_ev = REAL(get_selected_ritz_val(my_arnoldi, 1), dp) 741 742 ! save the ev as guess for the next time 743 IF (.NOT. ASSOCIATED(preconditioner_env%max_ev_vector)) ALLOCATE (preconditioner_env%max_ev_vector) 744 CALL get_selected_ritz_vector(my_arnoldi, 1, matrices(1)%matrix, preconditioner_env%max_ev_vector) 745 CALL deallocate_arnoldi_data(my_arnoldi) 746 DEALLOCATE (matrices) 747 748 ! Lets shift the occupied states a bit further up, -1.0 because we gonna subtract it from H 749 CALL dbcsr_add_on_diag(dbcsr_cThc, -0.5_dp) 750 ! Get the AO representation of the shift (see above why S is needed), W-matrix like object 751 CALL dbcsr_multiply("N", "N", 2.0_dp, dbcsr_sc, dbcsr_cThc, 0.0_dp, dbcsr_hc) 752 CALL dbcsr_multiply("N", "T", -1.0_dp, dbcsr_hc, dbcsr_sc, 1.0_dp, preconditioner_env%sparse_matrix) 753 754!-------------------------------------compute eigenvalues of H ---------------------------------------------- 755 ! Setup the arnoldi procedure to compute the lowest ev. if S is present this has to be the generalized ev 756 IF (PRESENT(matrix_s)) THEN 757 ALLOCATE (matrices(2)) 758 matrices(1)%matrix => preconditioner_env%sparse_matrix 759 matrices(2)%matrix => matrix_s 760 CALL setup_arnoldi_data(my_arnoldi, matrices, max_iter=20, threshold=2.0E-2_dp, selection_crit=3, & 761 nval_request=1, nrestarts=21, generalized_ev=.TRUE., iram=.FALSE.) 762 ELSE 763 ALLOCATE (matrices(1)) 764 matrices(1)%matrix => preconditioner_env%sparse_matrix 765 CALL setup_arnoldi_data(my_arnoldi, matrices, max_iter=20, threshold=2.0E-2_dp, selection_crit=3, & 766 nval_request=1, nrestarts=8, generalized_ev=.FALSE., iram=.FALSE.) 767 END IF 768 IF (ASSOCIATED(preconditioner_env%min_ev_vector)) & 769 CALL set_arnoldi_initial_vector(my_arnoldi, preconditioner_env%min_ev_vector) 770 771 ! compute the LUMO energy 772 CALL arnoldi_ev(matrices, my_arnoldi) 773 min_eV = REAL(get_selected_ritz_val(my_arnoldi, 1), dp) 774 775 ! save the lumo vector for restarting in the next step 776 IF (.NOT. ASSOCIATED(preconditioner_env%min_ev_vector)) ALLOCATE (preconditioner_env%min_ev_vector) 777 CALL get_selected_ritz_vector(my_arnoldi, 1, matrices(1)%matrix, preconditioner_env%min_ev_vector) 778 CALL deallocate_arnoldi_data(my_arnoldi) 779 DEALLOCATE (matrices) 780 781!-------------------------------------compute eigenvalues of H ---------------------------------------------- 782 ! Shift the Lumo to the 1.5*the computed energy_gap or the external energy gap value 783 ! The factor 1.5 is determined by trying. If the LUMO is positive, enough, just leave it alone 784 pre_shift = MAX(1.5_dp*(min_ev - max_ev), energy_gap) 785 IF (min_ev .LT. pre_shift) THEN 786 pre_shift = pre_shift - min_ev 787 ELSE 788 pre_shift = 0.0_dp 789 END IF 790 IF (PRESENT(matrix_s)) THEN 791 CALL dbcsr_add(preconditioner_env%sparse_matrix, matrix_s, 1.0_dp, pre_shift) 792 ELSE 793 CALL dbcsr_add_on_diag(preconditioner_env%sparse_matrix, pre_shift) 794 END IF 795 796 CALL dbcsr_release(mo_dbcsr) 797 CALL dbcsr_release(dbcsr_hc) 798 CALL dbcsr_release(dbcsr_sc) 799 CALL dbcsr_release(dbcsr_cThc) 800 801 CALL timestop(handle) 802 803 END SUBROUTINE make_full_single_inverse 804 805END MODULE preconditioner_makes 806 807