1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief module that contains the algorithms to perform an itrative 8!> diagonalization by the block-Davidson approach 9!> P. Blaha, et al J. Comp. Physics, 229, (2010), 453-460 10!> Iterative diagonalization in augmented plane wave based 11!> methods in electronic structure calculations 12!> \par History 13!> 05.2011 created [MI] 14!> \author MI 15! ************************************************************************************************** 16MODULE qs_scf_block_davidson 17 18 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 19 copy_fm_to_dbcsr,& 20 cp_dbcsr_m_by_n_from_row_template,& 21 cp_dbcsr_m_by_n_from_template,& 22 cp_dbcsr_sm_fm_multiply 23 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,& 24 cp_fm_scale_and_add,& 25 cp_fm_symm,& 26 cp_fm_transpose,& 27 cp_fm_triangular_invert,& 28 cp_fm_upper_to_full 29 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,& 30 cp_fm_cholesky_restore 31 USE cp_fm_diag, ONLY: choose_eigv_solver,& 32 cp_fm_power 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_release,& 39 cp_fm_set_all,& 40 cp_fm_to_fm,& 41 cp_fm_to_fm_submat,& 42 cp_fm_type,& 43 cp_fm_vectorsnorm 44 USE cp_gemm_interface, ONLY: cp_gemm 45 USE dbcsr_api, ONLY: & 46 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_get_diag, dbcsr_get_info, dbcsr_init_p, & 47 dbcsr_multiply, dbcsr_norm, dbcsr_norm_column, dbcsr_release_p, dbcsr_scale_by_vector, & 48 dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_real_default, dbcsr_type_symmetric 49 USE kinds, ONLY: dp 50 USE machine, ONLY: m_walltime 51 USE message_passing, ONLY: mp_sum 52 USE preconditioner, ONLY: apply_preconditioner 53 USE preconditioner_types, ONLY: preconditioner_type 54 USE qs_block_davidson_types, ONLY: davidson_type 55 USE qs_mo_types, ONLY: get_mo_set,& 56 mo_set_type 57#include "./base/base_uses.f90" 58 59 IMPLICIT NONE 60 PRIVATE 61 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_block_davidson' 62 63 PUBLIC :: generate_extended_space, generate_extended_space_sparse 64 65CONTAINS 66 67! ************************************************************************************************** 68!> \brief ... 69!> \param bdav_env ... 70!> \param mo_set ... 71!> \param matrix_h ... 72!> \param matrix_s ... 73!> \param output_unit ... 74!> \param preconditioner ... 75! ************************************************************************************************** 76 SUBROUTINE generate_extended_space(bdav_env, mo_set, matrix_h, matrix_s, output_unit, & 77 preconditioner) 78 79 TYPE(davidson_type) :: bdav_env 80 TYPE(mo_set_type), POINTER :: mo_set 81 TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s 82 INTEGER, INTENT(IN) :: output_unit 83 TYPE(preconditioner_type), OPTIONAL, POINTER :: preconditioner 84 85 CHARACTER(len=*), PARAMETER :: routineN = 'generate_extended_space', & 86 routineP = moduleN//':'//routineN 87 88 INTEGER :: handle, homo, i_first, i_last, imo, iter, j, jj, max_iter, n, nao, nmat, nmat2, & 89 nmo, nmo_converged, nmo_not_converged, nset, nset_conv, nset_not_conv 90 INTEGER, ALLOCATABLE, DIMENSION(:) :: iconv, inotconv 91 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iconv_set, inotconv_set 92 LOGICAL :: converged, do_apply_preconditioner 93 REAL(dp) :: lambda, max_norm, min_norm, t1, t2 94 REAL(dp), ALLOCATABLE, DIMENSION(:) :: ritz_coeff, vnorm 95 REAL(dp), DIMENSION(:), POINTER :: eig_not_conv, eigenvalues, evals 96 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 97 TYPE(cp_fm_type), POINTER :: c_conv, c_notconv, c_out, c_pz, c_z, & 98 h_block, h_fm, m_hc, m_sc, m_tmp, & 99 mo_coeff, mt_tmp, s_block, s_fm, & 100 v_block, w_block 101 TYPE(dbcsr_type), POINTER :: mo_coeff_b 102 103 CALL timeset(routineN, handle) 104 105 NULLIFY (mo_coeff, mo_coeff_b, eigenvalues) 106 107 do_apply_preconditioner = .FALSE. 108 IF (PRESENT(preconditioner)) do_apply_preconditioner = .TRUE. 109 CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, eigenvalues=eigenvalues, & 110 nao=nao, nmo=nmo, homo=homo) 111 IF (do_apply_preconditioner) THEN 112 max_iter = bdav_env%max_iter 113 ELSE 114 max_iter = 1 115 END IF 116 117 NULLIFY (c_conv, c_notconv, c_out, c_z, c_pz, m_hc, m_sc, m_tmp, mt_tmp) 118 NULLIFY (h_block, h_fm, s_block, s_fm, v_block, w_block) 119 NULLIFY (evals, eig_not_conv) 120 t1 = m_walltime() 121 IF (output_unit > 0) THEN 122 WRITE (output_unit, "(T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") & 123 " Cycle ", " conv. MOS ", " B2MAX ", " B2MIN ", " Time", REPEAT("-", 60) 124 END IF 125 126 ALLOCATE (iconv(nmo)) 127 ALLOCATE (inotconv(nmo)) 128 ALLOCATE (ritz_coeff(nmo)) 129 ALLOCATE (vnorm(nmo)) 130 131 converged = .FALSE. 132 DO iter = 1, max_iter 133 134 ! compute Ritz values 135 ritz_coeff = 0.0_dp 136 CALL cp_fm_create(m_hc, mo_coeff%matrix_struct, name="hc") 137 CALL cp_dbcsr_sm_fm_multiply(matrix_h, mo_coeff, m_hc, nmo) 138 CALL cp_fm_create(m_sc, mo_coeff%matrix_struct, name="sc") 139 CALL cp_dbcsr_sm_fm_multiply(matrix_s, mo_coeff, m_sc, nmo) 140 141 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo, ncol_global=nmo, & 142 context=mo_coeff%matrix_struct%context, & 143 para_env=mo_coeff%matrix_struct%para_env) 144 CALL cp_fm_create(m_tmp, fm_struct_tmp, name="matrix_tmp") 145 CALL cp_fm_struct_release(fm_struct_tmp) 146 147 CALL cp_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, mo_coeff, m_hc, 0.0_dp, m_tmp) 148 CALL cp_fm_get_diag(m_tmp, ritz_coeff) 149 CALL cp_fm_release(m_tmp) 150 151 ! Check for converged eigenvectors 152! CALL cp_fm_create(c_z,mo_coeff%matrix_struct,name="tmp") 153 c_z => bdav_env%matrix_z 154 c_pz => bdav_env%matrix_pz 155 CALL cp_fm_to_fm(m_sc, c_z) 156 CALL cp_fm_column_scale(c_z, ritz_coeff) 157 CALL cp_fm_scale_and_add(-1.0_dp, c_z, 1.0_dp, m_hc) 158 CALL cp_fm_vectorsnorm(c_z, vnorm) 159 160 nmo_converged = 0 161 nmo_not_converged = 0 162 max_norm = 0.0_dp 163 min_norm = 1.e10_dp 164 DO imo = 1, nmo 165 max_norm = MAX(max_norm, vnorm(imo)) 166 min_norm = MIN(min_norm, vnorm(imo)) 167 END DO 168 iconv = 0 169 inotconv = 0 170 DO imo = 1, nmo 171 IF (vnorm(imo) <= bdav_env%eps_iter) THEN 172 nmo_converged = nmo_converged + 1 173 iconv(nmo_converged) = imo 174 ELSE 175 nmo_not_converged = nmo_not_converged + 1 176 inotconv(nmo_not_converged) = imo 177 END IF 178 END DO 179 180 IF (nmo_converged > 0) THEN 181 ALLOCATE (iconv_set(nmo_converged, 2)) 182 ALLOCATE (inotconv_set(nmo_not_converged, 2)) 183 i_last = iconv(1) 184 nset = 0 185 DO j = 1, nmo_converged 186 imo = iconv(j) 187 188 IF (imo == i_last + 1) THEN 189 i_last = imo 190 iconv_set(nset, 2) = imo 191 ELSE 192 i_last = imo 193 nset = nset + 1 194 iconv_set(nset, 1) = imo 195 iconv_set(nset, 2) = imo 196 END IF 197 END DO 198 nset_conv = nset 199 200 i_last = inotconv(1) 201 nset = 0 202 DO j = 1, nmo_not_converged 203 imo = inotconv(j) 204 205 IF (imo == i_last + 1) THEN 206 i_last = imo 207 inotconv_set(nset, 2) = imo 208 ELSE 209 i_last = imo 210 nset = nset + 1 211 inotconv_set(nset, 1) = imo 212 inotconv_set(nset, 2) = imo 213 END IF 214 END DO 215 nset_not_conv = nset 216 CALL cp_fm_release(m_sc) 217 CALL cp_fm_release(m_hc) 218 NULLIFY (c_z, c_pz) 219 END IF 220 221 IF (REAL(nmo_converged, dp)/REAL(nmo, dp) > bdav_env%conv_percent) THEN 222 converged = .TRUE. 223 DEALLOCATE (iconv_set) 224 DEALLOCATE (inotconv_set) 225 t2 = m_walltime() 226 IF (output_unit > 0) THEN 227 WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') & 228 iter, nmo_converged, max_norm, min_norm, t2 - t1 229 230 WRITE (output_unit, *) " Reached convergence in ", iter, & 231 " Davidson iterations" 232 END IF 233 234 EXIT 235 END IF 236 237 IF (nmo_converged > 0) THEN 238 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nao, & 239 context=mo_coeff%matrix_struct%context, & 240 para_env=mo_coeff%matrix_struct%para_env) 241 !allocate h_fm 242 CALL cp_fm_create(h_fm, fm_struct_tmp, name="matrix_tmp") 243 !allocate s_fm 244 CALL cp_fm_create(s_fm, fm_struct_tmp, name="matrix_tmp") 245 !copy matrix_h in h_fm 246 CALL copy_dbcsr_to_fm(matrix_h, h_fm) 247 CALL cp_fm_upper_to_full(h_fm, s_fm) 248 249 !copy matrix_s in s_fm 250! CALL cp_fm_set_all(s_fm,0.0_dp) 251 CALL copy_dbcsr_to_fm(matrix_s, s_fm) 252 253 !allocate c_out 254 CALL cp_fm_create(c_out, fm_struct_tmp, name="matrix_tmp") 255 ! set c_out to zero 256 CALL cp_fm_set_all(c_out, 0.0_dp) 257 CALL cp_fm_struct_release(fm_struct_tmp) 258 259 !allocate c_conv 260 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_converged, & 261 context=mo_coeff%matrix_struct%context, & 262 para_env=mo_coeff%matrix_struct%para_env) 263 CALL cp_fm_create(c_conv, fm_struct_tmp, name="c_conv") 264 CALL cp_fm_set_all(c_conv, 0.0_dp) 265 !allocate m_tmp 266 CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_nxmc") 267 CALL cp_fm_struct_release(fm_struct_tmp) 268 END IF 269 270 !allocate c_notconv 271 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_not_converged, & 272 context=mo_coeff%matrix_struct%context, & 273 para_env=mo_coeff%matrix_struct%para_env) 274 CALL cp_fm_create(c_notconv, fm_struct_tmp, name="c_notconv") 275 CALL cp_fm_set_all(c_notconv, 0.0_dp) 276 IF (nmo_converged > 0) THEN 277 CALL cp_fm_create(m_hc, fm_struct_tmp, name="m_hc") 278 CALL cp_fm_create(m_sc, fm_struct_tmp, name="m_sc") 279 !allocate c_z 280 CALL cp_fm_create(c_z, fm_struct_tmp, name="c_z") 281 CALL cp_fm_create(c_pz, fm_struct_tmp, name="c_pz") 282 CALL cp_fm_set_all(c_z, 0.0_dp) 283 284 ! sum contributions to c_out 285 jj = 1 286 DO j = 1, nset_conv 287 i_first = iconv_set(j, 1) 288 i_last = iconv_set(j, 2) 289 n = i_last - i_first + 1 290 CALL cp_fm_to_fm_submat(mo_coeff, c_conv, nao, n, 1, i_first, 1, jj) 291 jj = jj + n 292 END DO 293 CALL cp_fm_symm('L', 'U', nao, nmo_converged, 1.0_dp, s_fm, c_conv, 0.0_dp, m_tmp) 294 CALL cp_gemm('N', 'T', nao, nao, nmo_converged, 1.0_dp, m_tmp, m_tmp, 0.0_dp, c_out) 295 296 ! project c_out out of H 297 lambda = 100.0_dp*ABS(eigenvalues(homo)) 298 CALL cp_fm_scale_and_add(lambda, c_out, 1.0_dp, h_fm) 299 CALL cp_fm_release(m_tmp) 300 CALL cp_fm_release(h_fm) 301 302 END IF 303 304 !allocate m_tmp 305 CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_nxm") 306 CALL cp_fm_struct_release(fm_struct_tmp) 307 IF (nmo_converged > 0) THEN 308 ALLOCATE (eig_not_conv(nmo_not_converged)) 309 jj = 1 310 DO j = 1, nset_not_conv 311 i_first = inotconv_set(j, 1) 312 i_last = inotconv_set(j, 2) 313 n = i_last - i_first + 1 314 CALL cp_fm_to_fm_submat(mo_coeff, c_notconv, nao, n, 1, i_first, 1, jj) 315 eig_not_conv(jj:jj + n - 1) = ritz_coeff(i_first:i_last) 316 jj = jj + n 317 END DO 318 CALL cp_gemm('N', 'N', nao, nmo_not_converged, nao, 1.0_dp, c_out, c_notconv, 0.0_dp, m_hc) 319 CALL cp_fm_symm('L', 'U', nao, nmo_not_converged, 1.0_dp, s_fm, c_notconv, 0.0_dp, m_sc) 320 ! extend suspace using only the not converged vectors 321 CALL cp_fm_to_fm(m_sc, m_tmp) 322 CALL cp_fm_column_scale(m_tmp, eig_not_conv) 323 CALL cp_fm_scale_and_add(-1.0_dp, m_tmp, 1.0_dp, m_hc) 324 DEALLOCATE (eig_not_conv) 325 CALL cp_fm_to_fm(m_tmp, c_z) 326 ELSE 327 CALL cp_fm_to_fm(mo_coeff, c_notconv) 328 END IF 329 330 !preconditioner 331 IF (do_apply_preconditioner) THEN 332 IF (preconditioner%in_use /= 0) THEN 333 CALL apply_preconditioner(preconditioner, c_z, c_pz) 334 ELSE 335 CALL cp_fm_to_fm(c_z, c_pz) 336 ENDIF 337 ELSE 338 CALL cp_fm_to_fm(c_z, c_pz) 339 END IF 340 CALL cp_fm_release(m_tmp) 341 342 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo_not_converged, ncol_global=nmo_not_converged, & 343 context=mo_coeff%matrix_struct%context, & 344 para_env=mo_coeff%matrix_struct%para_env) 345 346 CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_mxm") 347 CALL cp_fm_create(mt_tmp, fm_struct_tmp, name="mt_tmp_mxm") 348 CALL cp_fm_struct_release(fm_struct_tmp) 349 350 nmat = nmo_not_converged 351 nmat2 = 2*nmo_not_converged 352 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmat2, ncol_global=nmat2, & 353 context=mo_coeff%matrix_struct%context, & 354 para_env=mo_coeff%matrix_struct%para_env) 355 356 CALL cp_fm_create(s_block, fm_struct_tmp, name="sb") 357 CALL cp_fm_create(h_block, fm_struct_tmp, name="hb") 358 CALL cp_fm_create(v_block, fm_struct_tmp, name="vb") 359 CALL cp_fm_create(w_block, fm_struct_tmp, name="wb") 360 ALLOCATE (evals(nmat2)) 361 362 CALL cp_fm_struct_release(fm_struct_tmp) 363 364 ! compute CSC 365 CALL cp_fm_set_all(s_block, 0.0_dp, 1.0_dp) 366 367 ! compute CHC 368 CALL cp_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_notconv, m_hc, 0.0_dp, m_tmp) 369 CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1, 1) 370 371 ! compute ZSC 372 CALL cp_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_sc, 0.0_dp, m_tmp) 373 CALL cp_fm_to_fm_submat(m_tmp, s_block, nmat, nmat, 1, 1, 1 + nmat, 1) 374 CALL cp_fm_transpose(m_tmp, mt_tmp) 375 CALL cp_fm_to_fm_submat(mt_tmp, s_block, nmat, nmat, 1, 1, 1, 1 + nmat) 376 ! compute ZHC 377 CALL cp_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_hc, 0.0_dp, m_tmp) 378 CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1 + nmat, 1) 379 CALL cp_fm_transpose(m_tmp, mt_tmp) 380 CALL cp_fm_to_fm_submat(mt_tmp, h_block, nmat, nmat, 1, 1, 1, 1 + nmat) 381 382 CALL cp_fm_release(mt_tmp) 383 384 ! reuse m_sc and m_hc to computr HZ and SZ 385 IF (nmo_converged > 0) THEN 386 CALL cp_gemm('N', 'N', nao, nmat, nao, 1.0_dp, c_out, c_pz, 0.0_dp, m_hc) 387 CALL cp_fm_symm('L', 'U', nao, nmo_not_converged, 1.0_dp, s_fm, c_pz, 0.0_dp, m_sc) 388 389 CALL cp_fm_release(c_out) 390 CALL cp_fm_release(c_conv) 391 CALL cp_fm_release(s_fm) 392 ELSE 393 CALL cp_dbcsr_sm_fm_multiply(matrix_h, c_pz, m_hc, nmo) 394 CALL cp_dbcsr_sm_fm_multiply(matrix_s, c_pz, m_sc, nmo) 395 END IF 396 397 ! compute ZSZ 398 CALL cp_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_sc, 0.0_dp, m_tmp) 399 CALL cp_fm_to_fm_submat(m_tmp, s_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat) 400 ! compute ZHZ 401 CALL cp_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_hc, 0.0_dp, m_tmp) 402 CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat) 403 404 CALL cp_fm_release(m_sc) 405 406 ! solution of the reduced eigenvalues problem 407 CALL reduce_extended_space(s_block, h_block, v_block, w_block, evals, nmat2) 408 409 ! extract egenvectors 410 CALL cp_fm_to_fm_submat(v_block, m_tmp, nmat, nmat, 1, 1, 1, 1) 411 CALL cp_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, c_notconv, m_tmp, 0.0_dp, m_hc) 412 CALL cp_fm_to_fm_submat(v_block, m_tmp, nmat, nmat, 1 + nmat, 1, 1, 1) 413 CALL cp_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, c_pz, m_tmp, 1.0_dp, m_hc) 414 415 CALL cp_fm_release(m_tmp) 416 417 CALL cp_fm_release(c_notconv) 418 CALL cp_fm_release(s_block) 419 CALL cp_fm_release(h_block) 420 CALL cp_fm_release(w_block) 421 CALL cp_fm_release(v_block) 422 423 IF (nmo_converged > 0) THEN 424 CALL cp_fm_release(c_z) 425 CALL cp_fm_release(c_pz) 426 jj = 1 427 DO j = 1, nset_not_conv 428 i_first = inotconv_set(j, 1) 429 i_last = inotconv_set(j, 2) 430 n = i_last - i_first + 1 431 CALL cp_fm_to_fm_submat(m_hc, mo_coeff, nao, n, 1, jj, 1, i_first) 432 eigenvalues(i_first:i_last) = evals(jj:jj + n - 1) 433 jj = jj + n 434 END DO 435 DEALLOCATE (iconv_set) 436 DEALLOCATE (inotconv_set) 437 ELSE 438 CALL cp_fm_to_fm(m_hc, mo_coeff) 439 eigenvalues(1:nmo) = evals(1:nmo) 440 END IF 441 DEALLOCATE (evals) 442 443 CALL cp_fm_release(m_hc) 444 445 CALL copy_fm_to_dbcsr(mo_coeff, mo_coeff_b) !fm->dbcsr 446 447 t2 = m_walltime() 448 IF (output_unit > 0) THEN 449 WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') & 450 iter, nmo_converged, max_norm, min_norm, t2 - t1 451 END IF 452 t1 = m_walltime() 453 454 END DO ! iter 455 456 DEALLOCATE (iconv) 457 DEALLOCATE (inotconv) 458 DEALLOCATE (ritz_coeff) 459 DEALLOCATE (vnorm) 460 461 CALL timestop(handle) 462 END SUBROUTINE generate_extended_space 463 464! ************************************************************************************************** 465!> \brief ... 466!> \param bdav_env ... 467!> \param mo_set ... 468!> \param matrix_h ... 469!> \param matrix_s ... 470!> \param output_unit ... 471!> \param preconditioner ... 472! ************************************************************************************************** 473 SUBROUTINE generate_extended_space_sparse(bdav_env, mo_set, matrix_h, matrix_s, output_unit, & 474 preconditioner) 475 476 TYPE(davidson_type) :: bdav_env 477 TYPE(mo_set_type), POINTER :: mo_set 478 TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s 479 INTEGER, INTENT(IN) :: output_unit 480 TYPE(preconditioner_type), OPTIONAL, POINTER :: preconditioner 481 482 CHARACTER(len=*), PARAMETER :: routineN = 'generate_extended_space_sparse', & 483 routineP = moduleN//':'//routineN 484 485 INTEGER :: handle, homo, i_first, i_last, imo, iter, j, jj, k, max_iter, n, nao, nmat, & 486 nmat2, nmo, nmo_converged, nmo_not_converged, nset, nset_conv, nset_not_conv 487 INTEGER, ALLOCATABLE, DIMENSION(:) :: iconv, inotconv 488 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iconv_set, inotconv_set 489 LOGICAL :: converged, do_apply_preconditioner 490 REAL(dp) :: lambda, max_norm, min_norm, t1, t2 491 REAL(dp), ALLOCATABLE, DIMENSION(:) :: eig_not_conv, evals, ritz_coeff, vnorm 492 REAL(dp), DIMENSION(:), POINTER :: eigenvalues 493 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 494 TYPE(cp_fm_type), POINTER :: h_block, matrix_mm_fm, matrix_mmt_fm, matrix_nm_fm, & 495 matrix_z_fm, mo_coeff, mo_conv_fm, mo_notconv_fm, s_block, v_block, w_block 496 TYPE(dbcsr_type), POINTER :: c_out, matrix_hc, matrix_mm, matrix_pz, & 497 matrix_sc, matrix_z, mo_coeff_b, & 498 mo_conv, mo_notconv, smo_conv 499 500 CALL timeset(routineN, handle) 501 502 do_apply_preconditioner = .FALSE. 503 IF (PRESENT(preconditioner)) do_apply_preconditioner = .TRUE. 504 505 NULLIFY (mo_coeff, mo_coeff_b, matrix_hc, matrix_sc, matrix_z, matrix_pz, matrix_mm) 506 NULLIFY (mo_conv_fm, mo_notconv_fm, mo_conv, mo_notconv, smo_conv, c_out) 507 NULLIFY (matrix_mm_fm, matrix_mmt_fm, mo_coeff, matrix_nm_fm, matrix_z_fm) 508 NULLIFY (h_block, s_block, v_block, w_block) 509 NULLIFY (fm_struct_tmp) 510 CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, & 511 eigenvalues=eigenvalues, homo=homo, nao=nao, nmo=nmo) 512 IF (do_apply_preconditioner) THEN 513 max_iter = bdav_env%max_iter 514 ELSE 515 max_iter = 1 516 END IF 517 518 t1 = m_walltime() 519 IF (output_unit > 0) THEN 520 WRITE (output_unit, "(T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") & 521 " Cycle ", " conv. MOS ", " B2MAX ", " B2MIN ", " Time", REPEAT("-", 60) 522 END IF 523 524 ! Allocate array for Ritz values 525 ALLOCATE (ritz_coeff(nmo)) 526 ALLOCATE (iconv(nmo)) 527 ALLOCATE (inotconv(nmo)) 528 ALLOCATE (vnorm(nmo)) 529 530 converged = .FALSE. 531 DO iter = 1, max_iter 532 NULLIFY (c_out, mo_conv, mo_notconv_fm, mo_notconv) 533 ! Prepare HC and SC, using mo_coeff_b (sparse), these are still sparse 534 CALL dbcsr_init_p(matrix_hc) 535 CALL dbcsr_create(matrix_hc, template=mo_coeff_b, & 536 name="matrix_hc", & 537 matrix_type=dbcsr_type_no_symmetry, & 538 nze=0, data_type=dbcsr_type_real_default) 539 CALL dbcsr_init_p(matrix_sc) 540 CALL dbcsr_create(matrix_sc, template=mo_coeff_b, & 541 name="matrix_sc", & 542 matrix_type=dbcsr_type_no_symmetry, & 543 nze=0, data_type=dbcsr_type_real_default) 544 545 CALL dbcsr_get_info(mo_coeff_b, nfullrows_total=n, nfullcols_total=k) 546 CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_h, mo_coeff_b, 0.0_dp, matrix_hc, last_column=k) 547 CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, mo_coeff_b, 0.0_dp, matrix_sc, last_column=k) 548 549 ! compute Ritz values 550 ritz_coeff = 0.0_dp 551 ! Allocate Sparse matrices: nmoxnmo 552 ! matrix_mm 553 554 CALL dbcsr_init_p(matrix_mm) 555 CALL cp_dbcsr_m_by_n_from_template(matrix_mm, template=matrix_s, m=nmo, n=nmo, & 556 sym=dbcsr_type_no_symmetry) 557 558 CALL dbcsr_multiply('t', 'n', 1.0_dp, mo_coeff_b, matrix_hc, 0.0_dp, matrix_mm, last_column=k) 559 CALL dbcsr_get_diag(matrix_mm, ritz_coeff) 560 CALL mp_sum(ritz_coeff, mo_coeff%matrix_struct%para_env%group) 561 562 ! extended subspace P Z = P [H - theta S]C this ia another matrix of type and size as mo_coeff_b 563 CALL dbcsr_init_p(matrix_z) 564 CALL dbcsr_create(matrix_z, template=mo_coeff_b, & 565 name="matrix_z", & 566 matrix_type=dbcsr_type_no_symmetry, & 567 nze=0, data_type=dbcsr_type_real_default) 568 CALL dbcsr_copy(matrix_z, matrix_sc) 569 CALL dbcsr_scale_by_vector(matrix_z, ritz_coeff, side='right') 570 CALL dbcsr_add(matrix_z, matrix_hc, -1.0_dp, 1.0_dp) 571 572 ! Check for converged eigenvectors 573 vnorm = 0.0_dp 574 CALL dbcsr_norm(matrix_z, which_norm=dbcsr_norm_column, norm_vector=vnorm) 575 nmo_converged = 0 576 nmo_not_converged = 0 577 max_norm = 0.0_dp 578 min_norm = 1.e10_dp 579 DO imo = 1, nmo 580 max_norm = MAX(max_norm, vnorm(imo)) 581 min_norm = MIN(min_norm, vnorm(imo)) 582 END DO 583 iconv = 0 584 inotconv = 0 585 586 DO imo = 1, nmo 587 IF (vnorm(imo) <= bdav_env%eps_iter) THEN 588 nmo_converged = nmo_converged + 1 589 iconv(nmo_converged) = imo 590 ELSE 591 nmo_not_converged = nmo_not_converged + 1 592 inotconv(nmo_not_converged) = imo 593 END IF 594 END DO 595 596 IF (nmo_converged > 0) THEN 597 ALLOCATE (iconv_set(nmo_converged, 2)) 598 ALLOCATE (inotconv_set(nmo_not_converged, 2)) 599 i_last = iconv(1) 600 nset = 0 601 DO j = 1, nmo_converged 602 imo = iconv(j) 603 604 IF (imo == i_last + 1) THEN 605 i_last = imo 606 iconv_set(nset, 2) = imo 607 ELSE 608 i_last = imo 609 nset = nset + 1 610 iconv_set(nset, 1) = imo 611 iconv_set(nset, 2) = imo 612 END IF 613 END DO 614 nset_conv = nset 615 616 i_last = inotconv(1) 617 nset = 0 618 DO j = 1, nmo_not_converged 619 imo = inotconv(j) 620 621 IF (imo == i_last + 1) THEN 622 i_last = imo 623 inotconv_set(nset, 2) = imo 624 ELSE 625 i_last = imo 626 nset = nset + 1 627 inotconv_set(nset, 1) = imo 628 inotconv_set(nset, 2) = imo 629 END IF 630 END DO 631 nset_not_conv = nset 632 633 CALL dbcsr_release_p(matrix_hc) 634 CALL dbcsr_release_p(matrix_sc) 635 CALL dbcsr_release_p(matrix_z) 636 CALL dbcsr_release_p(matrix_mm) 637 END IF 638 639 IF (REAL(nmo_converged, dp)/REAL(nmo, dp) > bdav_env%conv_percent) THEN 640 DEALLOCATE (iconv_set) 641 642 DEALLOCATE (inotconv_set) 643 644 converged = .TRUE. 645 t2 = m_walltime() 646 IF (output_unit > 0) THEN 647 WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') & 648 iter, nmo_converged, max_norm, min_norm, t2 - t1 649 650 WRITE (output_unit, *) " Reached convergence in ", iter, & 651 " Davidson iterations" 652 END IF 653 654 EXIT 655 END IF 656 657 IF (nmo_converged > 0) THEN 658 659 !allocate mo_conv_fm 660 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_converged, & 661 context=mo_coeff%matrix_struct%context, & 662 para_env=mo_coeff%matrix_struct%para_env) 663 CALL cp_fm_create(mo_conv_fm, fm_struct_tmp, name="mo_conv_fm") 664 665 CALL cp_fm_struct_release(fm_struct_tmp) 666 667 ! extract mo_conv from mo_coeff full matrix 668 jj = 1 669 DO j = 1, nset_conv 670 i_first = iconv_set(j, 1) 671 i_last = iconv_set(j, 2) 672 n = i_last - i_first + 1 673 CALL cp_fm_to_fm_submat(mo_coeff, mo_conv_fm, nao, n, 1, i_first, 1, jj) 674 jj = jj + n 675 END DO 676 677 ! allocate c_out sparse matrix, to project out the converged MOS 678 CALL dbcsr_init_p(c_out) 679 CALL dbcsr_create(c_out, template=matrix_s, & 680 name="c_out", & 681 matrix_type=dbcsr_type_symmetric, & 682 nze=0, data_type=dbcsr_type_real_default) 683 684 ! allocate mo_conv sparse 685 CALL dbcsr_init_p(mo_conv) 686 CALL cp_dbcsr_m_by_n_from_row_template(mo_conv, template=matrix_s, n=nmo_converged, & 687 sym=dbcsr_type_no_symmetry) 688 689 CALL dbcsr_init_p(smo_conv) 690 CALL cp_dbcsr_m_by_n_from_row_template(smo_conv, template=matrix_s, n=nmo_converged, & 691 sym=dbcsr_type_no_symmetry) 692 693 CALL copy_fm_to_dbcsr(mo_conv_fm, mo_conv) !fm->dbcsr 694 695 CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, mo_conv, 0.0_dp, smo_conv, last_column=nmo_converged) 696 CALL dbcsr_multiply('n', 't', 1.0_dp, smo_conv, smo_conv, 0.0_dp, c_out, last_column=nao) 697 ! project c_out out of H 698 lambda = 100.0_dp*ABS(eigenvalues(homo)) 699 CALL dbcsr_add(c_out, matrix_h, lambda, 1.0_dp) 700 701 CALL dbcsr_release_p(mo_conv) 702 CALL dbcsr_release_p(smo_conv) 703 CALL cp_fm_release(mo_conv_fm) 704 705 !allocate c_notconv_fm 706 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_not_converged, & 707 context=mo_coeff%matrix_struct%context, & 708 para_env=mo_coeff%matrix_struct%para_env) 709 CALL cp_fm_create(mo_notconv_fm, fm_struct_tmp, name="mo_notconv_fm") 710 CALL cp_fm_struct_release(fm_struct_tmp) 711 712 ! extract mo_notconv from mo_coeff full matrix 713 ALLOCATE (eig_not_conv(nmo_not_converged)) 714 715 jj = 1 716 DO j = 1, nset_not_conv 717 i_first = inotconv_set(j, 1) 718 i_last = inotconv_set(j, 2) 719 n = i_last - i_first + 1 720 CALL cp_fm_to_fm_submat(mo_coeff, mo_notconv_fm, nao, n, 1, i_first, 1, jj) 721 eig_not_conv(jj:jj + n - 1) = ritz_coeff(i_first:i_last) 722 jj = jj + n 723 END DO 724 725 ! allocate mo_conv sparse 726 CALL dbcsr_init_p(mo_notconv) 727 CALL cp_dbcsr_m_by_n_from_row_template(mo_notconv, template=matrix_s, n=nmo_not_converged, & 728 sym=dbcsr_type_no_symmetry) 729 730 CALL dbcsr_init_p(matrix_hc) 731 CALL cp_dbcsr_m_by_n_from_row_template(matrix_hc, template=matrix_s, n=nmo_not_converged, & 732 sym=dbcsr_type_no_symmetry) 733 734 CALL dbcsr_init_p(matrix_sc) 735 CALL cp_dbcsr_m_by_n_from_row_template(matrix_sc, template=matrix_s, n=nmo_not_converged, & 736 sym=dbcsr_type_no_symmetry) 737 738 CALL dbcsr_init_p(matrix_z) 739 CALL cp_dbcsr_m_by_n_from_row_template(matrix_z, template=matrix_s, n=nmo_not_converged, & 740 sym=dbcsr_type_no_symmetry) 741 742 CALL copy_fm_to_dbcsr(mo_notconv_fm, mo_notconv) !fm->dbcsr 743 744 CALL dbcsr_multiply('n', 'n', 1.0_dp, c_out, mo_notconv, 0.0_dp, matrix_hc, & 745 last_column=nmo_not_converged) 746 CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, mo_notconv, 0.0_dp, matrix_sc, & 747 last_column=nmo_not_converged) 748 749 CALL dbcsr_copy(matrix_z, matrix_sc) 750 CALL dbcsr_scale_by_vector(matrix_z, eig_not_conv, side='right') 751 CALL dbcsr_add(matrix_z, matrix_hc, -1.0_dp, 1.0_dp) 752 753 DEALLOCATE (eig_not_conv) 754 755 ! matrix_mm 756 CALL dbcsr_init_p(matrix_mm) 757 CALL cp_dbcsr_m_by_n_from_template(matrix_mm, template=matrix_s, m=nmo_not_converged, n=nmo_not_converged, & 758 sym=dbcsr_type_no_symmetry) 759 760 CALL dbcsr_multiply('t', 'n', 1.0_dp, mo_notconv, matrix_hc, 0.0_dp, matrix_mm, & 761 last_column=nmo_not_converged) 762 763 ELSE 764 mo_notconv => mo_coeff_b 765 mo_notconv_fm => mo_coeff 766 c_out => matrix_h 767 END IF 768 769 ! allocate matrix_pz using as template matrix_z 770 CALL dbcsr_init_p(matrix_pz) 771 CALL dbcsr_create(matrix_pz, template=matrix_z, & 772 name="matrix_pz", & 773 matrix_type=dbcsr_type_no_symmetry, & 774 nze=0, data_type=dbcsr_type_real_default) 775 776 IF (do_apply_preconditioner) THEN 777 IF (preconditioner%in_use /= 0) THEN 778 CALL apply_preconditioner(preconditioner, matrix_z, matrix_pz) 779 ELSE 780 CALL dbcsr_copy(matrix_pz, matrix_z) 781 END IF 782 ELSE 783 CALL dbcsr_copy(matrix_pz, matrix_z) 784 END IF 785 786 !allocate NMOxNMO full matrices 787 nmat = nmo_not_converged 788 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmat, ncol_global=nmat, & 789 context=mo_coeff%matrix_struct%context, & 790 para_env=mo_coeff%matrix_struct%para_env) 791 CALL cp_fm_create(matrix_mm_fm, fm_struct_tmp, name="m_tmp_mxm") 792 CALL cp_fm_create(matrix_mmt_fm, fm_struct_tmp, name="mt_tmp_mxm") 793 CALL cp_fm_struct_release(fm_struct_tmp) 794 795 !allocate 2NMOx2NMO full matrices 796 nmat2 = 2*nmo_not_converged 797 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmat2, ncol_global=nmat2, & 798 context=mo_coeff%matrix_struct%context, & 799 para_env=mo_coeff%matrix_struct%para_env) 800 801 CALL cp_fm_create(s_block, fm_struct_tmp, name="sb") 802 CALL cp_fm_create(h_block, fm_struct_tmp, name="hb") 803 CALL cp_fm_create(v_block, fm_struct_tmp, name="vb") 804 CALL cp_fm_create(w_block, fm_struct_tmp, name="wb") 805 ALLOCATE (evals(nmat2)) 806 CALL cp_fm_struct_release(fm_struct_tmp) 807 808 ! compute CSC 809 CALL cp_fm_set_all(s_block, 0.0_dp, 1.0_dp) 810 ! compute CHC 811 CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm) 812 CALL cp_fm_to_fm_submat(matrix_mm_fm, h_block, nmat, nmat, 1, 1, 1, 1) 813 814 ! compute the bottom left ZSC (top right is transpose) 815 CALL dbcsr_multiply('t', 'n', 1.0_dp, matrix_pz, matrix_sc, 0.0_dp, matrix_mm, last_column=nmat) 816 ! set the bottom left part of S[C,Z] block matrix ZSC 817 !copy sparse to full 818 CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm) 819 CALL cp_fm_to_fm_submat(matrix_mm_fm, s_block, nmat, nmat, 1, 1, 1 + nmat, 1) 820 CALL cp_fm_transpose(matrix_mm_fm, matrix_mmt_fm) 821 CALL cp_fm_to_fm_submat(matrix_mmt_fm, s_block, nmat, nmat, 1, 1, 1, 1 + nmat) 822 823 ! compute the bottom left ZHC (top right is transpose) 824 CALL dbcsr_multiply('t', 'n', 1.0_dp, matrix_pz, matrix_hc, 0.0_dp, matrix_mm, last_column=nmat) 825 ! set the bottom left part of S[C,Z] block matrix ZHC 826 CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm) 827 CALL cp_fm_to_fm_submat(matrix_mm_fm, h_block, nmat, nmat, 1, 1, 1 + nmat, 1) 828 CALL cp_fm_transpose(matrix_mm_fm, matrix_mmt_fm) 829 CALL cp_fm_to_fm_submat(matrix_mmt_fm, h_block, nmat, nmat, 1, 1, 1, 1 + nmat) 830 831 CALL cp_fm_release(matrix_mmt_fm) 832 833 ! (reuse matrix_sc and matrix_hc to computr HZ and SZ) 834 CALL dbcsr_get_info(matrix_pz, nfullrows_total=n, nfullcols_total=k) 835 CALL dbcsr_multiply('n', 'n', 1.0_dp, c_out, matrix_pz, 0.0_dp, matrix_hc, last_column=k) 836 CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, matrix_pz, 0.0_dp, matrix_sc, last_column=k) 837 838 ! compute the bottom right ZSZ 839 CALL dbcsr_multiply('t', 'n', 1.0_dp, matrix_pz, matrix_sc, 0.0_dp, matrix_mm, last_column=k) 840 ! set the bottom right part of S[C,Z] block matrix ZSZ 841 CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm) 842 CALL cp_fm_to_fm_submat(matrix_mm_fm, s_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat) 843 844 ! compute the bottom right ZHZ 845 CALL dbcsr_multiply('t', 'n', 1.0_dp, matrix_pz, matrix_hc, 0.0_dp, matrix_mm, last_column=k) 846 ! set the bottom right part of H[C,Z] block matrix ZHZ 847 CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm) 848 CALL cp_fm_to_fm_submat(matrix_mm_fm, h_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat) 849 850 CALL dbcsr_release_p(matrix_mm) 851 CALL dbcsr_release_p(matrix_sc) 852 CALL dbcsr_release_p(matrix_hc) 853 854 CALL reduce_extended_space(s_block, h_block, v_block, w_block, evals, nmat2) 855 856 ! allocate two (nao x nmat) full matrix 857 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmat, & 858 context=mo_coeff%matrix_struct%context, & 859 para_env=mo_coeff%matrix_struct%para_env) 860 CALL cp_fm_create(matrix_nm_fm, fm_struct_tmp, name="m_nxm") 861 CALL cp_fm_create(matrix_z_fm, fm_struct_tmp, name="m_nxm") 862 CALL cp_fm_struct_release(fm_struct_tmp) 863 864 CALL copy_dbcsr_to_fm(matrix_pz, matrix_z_fm) 865 ! extract egenvectors 866 CALL cp_fm_to_fm_submat(v_block, matrix_mm_fm, nmat, nmat, 1, 1, 1, 1) 867 CALL cp_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, mo_notconv_fm, matrix_mm_fm, 0.0_dp, matrix_nm_fm) 868 CALL cp_fm_to_fm_submat(v_block, matrix_mm_fm, nmat, nmat, 1 + nmat, 1, 1, 1) 869 CALL cp_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, matrix_z_fm, matrix_mm_fm, 1.0_dp, matrix_nm_fm) 870 871 CALL dbcsr_release_p(matrix_z) 872 CALL dbcsr_release_p(matrix_pz) 873 CALL cp_fm_release(matrix_z_fm) 874 CALL cp_fm_release(s_block) 875 CALL cp_fm_release(h_block) 876 CALL cp_fm_release(w_block) 877 CALL cp_fm_release(v_block) 878 CALL cp_fm_release(matrix_mm_fm) 879 880 ! in case some vector are already converged only a subset of vectors are copied in the MOS 881 IF (nmo_converged > 0) THEN 882 jj = 1 883 DO j = 1, nset_not_conv 884 i_first = inotconv_set(j, 1) 885 i_last = inotconv_set(j, 2) 886 n = i_last - i_first + 1 887 CALL cp_fm_to_fm_submat(matrix_nm_fm, mo_coeff, nao, n, 1, jj, 1, i_first) 888 eigenvalues(i_first:i_last) = evals(jj:jj + n - 1) 889 jj = jj + n 890 END DO 891 DEALLOCATE (iconv_set) 892 DEALLOCATE (inotconv_set) 893 894 CALL dbcsr_release_p(mo_notconv) 895 CALL dbcsr_release_p(c_out) 896 CALL cp_fm_release(mo_notconv_fm) 897 ELSE 898 CALL cp_fm_to_fm(matrix_nm_fm, mo_coeff) 899 eigenvalues(1:nmo) = evals(1:nmo) 900 END IF 901 DEALLOCATE (evals) 902 903 CALL cp_fm_release(matrix_nm_fm) 904 CALL copy_fm_to_dbcsr(mo_coeff, mo_coeff_b) !fm->dbcsr 905 906 t2 = m_walltime() 907 IF (output_unit > 0) THEN 908 WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') & 909 iter, nmo_converged, max_norm, min_norm, t2 - t1 910 END IF 911 t1 = m_walltime() 912 913 END DO ! iter 914 915 DEALLOCATE (ritz_coeff) 916 DEALLOCATE (iconv) 917 DEALLOCATE (inotconv) 918 DEALLOCATE (vnorm) 919 920 CALL timestop(handle) 921 922 END SUBROUTINE generate_extended_space_sparse 923 924! ************************************************************************************************** 925 926! ************************************************************************************************** 927!> \brief ... 928!> \param s_block ... 929!> \param h_block ... 930!> \param v_block ... 931!> \param w_block ... 932!> \param evals ... 933!> \param ndim ... 934! ************************************************************************************************** 935 SUBROUTINE reduce_extended_space(s_block, h_block, v_block, w_block, evals, ndim) 936 937 TYPE(cp_fm_type), POINTER :: s_block, h_block, v_block, w_block 938 REAL(dp), DIMENSION(:) :: evals 939 INTEGER :: ndim 940 941 CHARACTER(len=*), PARAMETER :: routineN = 'reduce_extended_space', & 942 routineP = moduleN//':'//routineN 943 944 INTEGER :: handle, info 945 946 CALL timeset(routineN, handle) 947 948 CALL cp_fm_to_fm(s_block, w_block) 949 CALL cp_fm_cholesky_decompose(s_block, info_out=info) 950 IF (info == 0) THEN 951 CALL cp_fm_triangular_invert(s_block) 952 CALL cp_fm_cholesky_restore(H_block, ndim, S_block, w_block, "MULTIPLY", pos="RIGHT") 953 CALL cp_fm_cholesky_restore(w_block, ndim, S_block, H_block, "MULTIPLY", pos="LEFT", transa="T") 954 CALL choose_eigv_solver(H_block, w_block, evals) 955 CALL cp_fm_cholesky_restore(w_block, ndim, S_block, v_block, "MULTIPLY") 956 ELSE 957! S^(-1/2) 958 CALL cp_fm_power(w_block, s_block, -0.5_dp, 1.0E-5_dp, info) 959 CALL cp_fm_to_fm(w_block, s_block) 960 CALL cp_gemm('N', 'N', ndim, ndim, ndim, 1.0_dp, H_block, s_block, 0.0_dp, w_block) 961 CALL cp_gemm('N', 'N', ndim, ndim, ndim, 1.0_dp, s_block, w_block, 0.0_dp, H_block) 962 CALL choose_eigv_solver(H_block, w_block, evals) 963 CALL cp_gemm('N', 'N', ndim, ndim, ndim, 1.0_dp, s_block, w_block, 0.0_dp, v_block) 964 END IF 965 966 CALL timestop(handle) 967 968 END SUBROUTINE reduce_extended_space 969 970END MODULE qs_scf_block_davidson 971