1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7MODULE qs_tddfpt_eigensolver 8 USE cp_blacs_env, ONLY: cp_blacs_env_type 9 USE cp_control_types, ONLY: dft_control_type,& 10 tddfpt_control_type 11 USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t,& 12 cp_dbcsr_sm_fm_multiply 13 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,& 14 cp_fm_symm,& 15 cp_fm_trace 16 USE cp_fm_diag, ONLY: cp_fm_syevd 17 USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,& 18 fm_pools_create_fm_vect,& 19 fm_pools_give_back_fm_vect 20 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 21 cp_fm_struct_p_type,& 22 cp_fm_struct_release,& 23 cp_fm_struct_type 24 USE cp_fm_types, ONLY: cp_fm_create,& 25 cp_fm_get_element,& 26 cp_fm_p_type,& 27 cp_fm_release,& 28 cp_fm_set_all,& 29 cp_fm_set_element,& 30 cp_fm_to_fm 31 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit,& 32 cp_to_string 33 USE cp_para_types, ONLY: cp_para_env_type 34 USE dbcsr_api, ONLY: dbcsr_p_type,& 35 dbcsr_set 36 USE input_constants, ONLY: tddfpt_davidson,& 37 tddfpt_lanczos 38 USE kinds, ONLY: default_string_length,& 39 dp 40 USE physcon, ONLY: evolt 41 USE qs_environment_types, ONLY: get_qs_env,& 42 qs_environment_type 43 USE qs_matrix_pools, ONLY: mpools_get 44 USE qs_p_env_methods, ONLY: p_op_l1,& 45 p_op_l2,& 46 p_postortho,& 47 p_preortho 48 USE qs_p_env_types, ONLY: qs_p_env_type 49 USE qs_tddfpt_types, ONLY: tddfpt_env_type 50 USE qs_tddfpt_utils, ONLY: co_initial_guess,& 51 normalize,& 52 reorthogonalize 53#include "./base/base_uses.f90" 54 55 IMPLICIT NONE 56 57 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt_eigensolver' 58 59 PRIVATE 60 61 PUBLIC :: eigensolver 62 63CONTAINS 64 65! ************************************************************************************************** 66!> \brief ... 67!> \param p_env ... 68!> \param qs_env ... 69!> \param t_env ... 70! ************************************************************************************************** 71 SUBROUTINE eigensolver(p_env, qs_env, t_env) 72 73 TYPE(qs_p_env_type), POINTER :: p_env 74 TYPE(qs_environment_type), POINTER :: qs_env 75 TYPE(tddfpt_env_type), INTENT(INOUT) :: t_env 76 77 CHARACTER(len=*), PARAMETER :: routineN = 'eigensolver', routineP = moduleN//':'//routineN 78 79 INTEGER :: handle, n_ev, nspins, output_unit, & 80 restarts 81 LOGICAL :: do_kernel_save 82 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ievals 83 TYPE(dft_control_type), POINTER :: dft_control 84 85 CALL timeset(routineN, handle) 86 87 NULLIFY (dft_control) 88 89 output_unit = cp_logger_get_default_io_unit() 90 91 CALL get_qs_env(qs_env, dft_control=dft_control) 92 n_ev = dft_control%tddfpt_control%n_ev 93 nspins = dft_control%nspins 94 95 ALLOCATE (ievals(n_ev)) 96 97 !---------------! 98 ! initial guess ! 99 !---------------! 100 do_kernel_save = dft_control%tddfpt_control%do_kernel 101 dft_control%tddfpt_control%do_kernel = .FALSE. 102 IF (output_unit > 0) THEN 103 WRITE (output_unit, *) " Generating initial guess" 104 WRITE (output_unit, *) 105 END IF 106 IF (ASSOCIATED(dft_control%tddfpt_control%lumos)) THEN 107 CALL co_initial_guess(t_env%evecs, ievals, n_ev, qs_env) 108 ELSE 109 IF (output_unit > 0) WRITE (output_unit, *) "LUMOS are needed in TDDFPT!" 110 CPABORT("") 111 END IF 112 DO restarts = 1, dft_control%tddfpt_control%n_restarts 113 IF (iterative_solver(ievals, t_env, p_env, qs_env, ievals)) EXIT 114 IF (output_unit > 0) THEN 115 WRITE (output_unit, *) " Restarting" 116 WRITE (output_unit, *) 117 END IF 118 END DO 119 dft_control%tddfpt_control%do_kernel = do_kernel_save 120 121 !-----------------! 122 ! call the solver ! 123 !-----------------! 124 IF (output_unit > 0) THEN 125 WRITE (output_unit, *) 126 WRITE (output_unit, *) " Doing TDDFPT calculation" 127 WRITE (output_unit, *) 128 END IF 129 DO restarts = 1, dft_control%tddfpt_control%n_restarts 130 IF (iterative_solver(ievals, t_env, p_env, qs_env, t_env%evals)) EXIT 131 IF (output_unit > 0) THEN 132 WRITE (output_unit, *) " Restarting" 133 WRITE (output_unit, *) 134 END IF 135 END DO 136 137 !---------! 138 ! cleanup ! 139 !---------! 140 DEALLOCATE (ievals) 141 142 CALL timestop(handle) 143 144 END SUBROUTINE eigensolver 145 146 ! in_evals : approximations to the eigenvalues for the preconditioner 147 ! t_env : TD-DFT environment values 148 ! p_env : perturbation environment values 149 ! qs_env : general Quickstep environment values 150 ! out_evals : the resulting eigenvalues 151 ! error : used for error handling 152 ! 153 ! res : the function will return wheter the eigenvalues are converged or not 154 155! ************************************************************************************************** 156!> \brief ... 157!> \param in_evals ... 158!> \param t_env ... 159!> \param p_env ... 160!> \param qs_env ... 161!> \param out_evals ... 162!> \return ... 163! ************************************************************************************************** 164 FUNCTION iterative_solver(in_evals, & 165 t_env, p_env, qs_env, & 166 out_evals) RESULT(res) 167 168 REAL(KIND=dp), DIMENSION(:) :: in_evals 169 TYPE(tddfpt_env_type), INTENT(INOUT) :: t_env 170 TYPE(qs_p_env_type), POINTER :: p_env 171 TYPE(qs_environment_type), POINTER :: qs_env 172 REAL(kind=dp), DIMENSION(:), OPTIONAL :: out_evals 173 LOGICAL :: res 174 175 CHARACTER(len=*), PARAMETER :: routineN = 'iterative_solver', & 176 routineP = moduleN//':'//routineN 177 178 CHARACTER :: mode 179 INTEGER :: col, handle, i, iev, iter, j, k, & 180 max_krylovspace_dim, max_kv, n_ev, & 181 n_kv, nspins, output_unit, row, spin 182 INTEGER, ALLOCATABLE, DIMENSION(:) :: must_improve 183 REAL(dp) :: Atilde_ij, convergence, tmp, tmp2 184 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals_difference, evals_tmp 185 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: evals 186 TYPE(cp_blacs_env_type), POINTER :: blacs_env 187 TYPE(cp_fm_p_type) :: Atilde, Us 188 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: R, X 189 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: Ab, b, Sb 190 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools 191 TYPE(cp_fm_struct_p_type), DIMENSION(:), POINTER :: kv_fm_struct 192 TYPE(cp_fm_struct_type), POINTER :: tilde_fm_struct 193 TYPE(cp_para_env_type), POINTER :: para_env 194 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 195 TYPE(dft_control_type), POINTER :: dft_control 196 TYPE(tddfpt_control_type), POINTER :: tddfpt_control 197 198 res = .FALSE. 199 200 CALL timeset(routineN, handle) 201 202 NULLIFY (ao_mo_fm_pools, X, R, b, Ab, Sb, tddfpt_control, & 203 tilde_fm_struct, Atilde%matrix, Us%matrix, matrix_s, dft_control, & 204 para_env, blacs_env) 205 206 CALL get_qs_env(qs_env, & 207 matrix_s=matrix_s, & 208 dft_control=dft_control, & 209 para_env=para_env, & 210 blacs_env=blacs_env) 211 212 tddfpt_control => dft_control%tddfpt_control 213 output_unit = cp_logger_get_default_io_unit() 214 n_ev = tddfpt_control%n_ev 215 nspins = dft_control%nspins 216 217 IF (dft_control%tddfpt_control%diag_method == tddfpt_lanczos) THEN 218 mode = 'L' 219 ELSE IF (dft_control%tddfpt_control%diag_method == tddfpt_davidson) THEN 220 mode = 'D' 221 END IF 222 223 !-----------------------------------------! 224 ! determine the size of the problem ! 225 ! and how many krylov space vetors to use ! 226 !-----------------------------------------! 227 max_krylovspace_dim = SUM(p_env%n_ao(1:nspins)*p_env%n_mo(1:nspins)) 228 max_kv = tddfpt_control%max_kv 229 IF (max_krylovspace_dim <= max_kv) THEN 230 max_kv = max_krylovspace_dim 231 IF (output_unit > 0) THEN 232 WRITE (output_unit, *) " Setting the maximum number of krylov vectors to ", max_kv, "!!" 233 END IF 234 END IF 235! CPPostcondition(max_krylovspace_dim>=max_kv,cp_failure_level,routineP,failure) 236 237 !----------------------! 238 ! allocate the vectors ! 239 !----------------------! 240 CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools) 241 CALL fm_pools_create_fm_vect(ao_mo_fm_pools, X, name=routineP//":X") 242 CALL fm_pools_create_fm_vect(ao_mo_fm_pools, R, name=routineP//":R") 243 244 ALLOCATE (evals_difference(n_ev)) 245 246 ALLOCATE (must_improve(n_ev)) 247 248 ALLOCATE (evals(max_kv, 0:max_kv)) 249 ALLOCATE (evals_tmp(max_kv)) 250 251 ALLOCATE (b(max_kv, nspins), Ab(max_kv, nspins), & 252 Sb(max_kv, nspins)) 253 254 ALLOCATE (kv_fm_struct(nspins)) 255 256 DO spin = 1, nspins 257 NULLIFY (kv_fm_struct(spin)%struct) 258 CALL cp_fm_struct_create(kv_fm_struct(spin)%struct, para_env, blacs_env, & 259 p_env%n_ao(spin), p_env%n_mo(spin)) 260 END DO 261 DO spin = 1, nspins 262 DO i = 1, max_kv 263 NULLIFY (b(i, spin)%matrix, Ab(i, spin)%matrix, Sb(i, spin)%matrix) 264 END DO 265 END DO 266 267 IF (output_unit > 0) THEN 268 WRITE (output_unit, '(2X,A,T69,A)') & 269 "nvec", "Convergence" 270 WRITE (output_unit, '(2X,A)') & 271 "-----------------------------------------------------------------------------" 272 END IF 273 274 iter = 1 275 k = 0 276 n_kv = n_ev 277 iteration: DO 278 279 CALL allocate_krylov_vectors(b, "b-", k + 1, n_kv, nspins, kv_fm_struct) 280 CALL allocate_krylov_vectors(Ab, "Ab-", k + 1, n_kv, nspins, kv_fm_struct) 281 CALL allocate_krylov_vectors(Sb, "Sb-", k + 1, n_kv, nspins, kv_fm_struct) 282 283 DO i = 1, n_kv 284 k = k + 1 285 286 IF (k <= SIZE(t_env%evecs, 1)) THEN ! the first iteration 287 288 ! take the initial guess 289 DO spin = 1, nspins 290 CALL cp_fm_to_fm(t_env%evecs(k, spin)%matrix, b(k, spin)%matrix) 291 END DO 292 293 ELSE 294 295 ! create a new vector 296 IF (mode == 'L') THEN 297 298 DO spin = 1, nspins 299 IF (tddfpt_control%invert_S) THEN 300 CALL cp_fm_symm('L', 'U', p_env%n_ao(spin), p_env%n_mo(spin), & 301 1.0_dp, t_env%invS(spin)%matrix, Ab(k - 1, spin)%matrix, & 302 0.0_dp, b(k, spin)%matrix) 303 ELSE 304 CALL cp_fm_to_fm(Ab(k - 1, spin)%matrix, b(k, spin)%matrix) 305 END IF 306 END DO 307 308 ELSE IF (mode == 'D') THEN 309 310 iev = must_improve(i) 311 ! create the new davidson vector 312 DO spin = 1, nspins 313 314 CALL cp_fm_set_all(R(spin)%matrix, 0.0_dp) 315 DO j = 1, k - i 316 CALL cp_fm_to_fm(Ab(j, spin)%matrix, X(spin)%matrix) 317 CALL cp_fm_scale_and_add(1.0_dp, X(spin)%matrix, & 318 -evals(iev, iter - 1), Sb(j, spin)%matrix) 319 CALL cp_fm_get_element(Us%matrix, j, iev, tmp) 320 CALL cp_fm_scale_and_add(1.0_dp, R(spin)%matrix, & 321 tmp, X(spin)%matrix) 322 END DO 323 324 IF (tddfpt_control%invert_S) THEN 325 CALL cp_fm_symm('L', 'U', p_env%n_ao(spin), p_env%n_mo(spin), & 326 1.0_dp, t_env%invS(spin)%matrix, R(spin)%matrix, & 327 0.0_dp, X(spin)%matrix) 328 ELSE 329 CALL cp_fm_to_fm(R(spin)%matrix, X(spin)%matrix) 330 END IF 331 332 !----------------! 333 ! preconditioner ! 334 !----------------! 335 IF (dft_control%tddfpt_control%precond) THEN 336 DO col = 1, p_env%n_mo(spin) 337 IF (col <= n_ev) THEN 338 tmp2 = ABS(evals(iev, iter - 1) - in_evals(col)) 339 ELSE 340 tmp2 = ABS(evals(iev, iter - 1) - (in_evals(n_ev) + 10.0_dp)) 341 END IF 342 ! protect against division by 0 by a introducing a cutoff. 343 tmp2 = MAX(tmp2, 100*EPSILON(1.0_dp)) 344 DO row = 1, p_env%n_ao(spin) 345 CALL cp_fm_get_element(X(spin)%matrix, row, col, tmp) 346 CALL cp_fm_set_element(b(k, spin)%matrix, row, col, tmp/tmp2) 347 END DO 348 END DO 349 ELSE 350 CALL cp_fm_to_fm(X(spin)%matrix, b(k, spin)%matrix) 351 END IF 352 353 END DO 354 355 ELSE 356 IF (output_unit > 0) WRITE (output_unit, *) "unknown mode" 357 CPABORT("") 358 END IF 359 360 END IF 361 362 CALL p_preortho(p_env, qs_env, b(k, :)) 363 DO j = 1, tddfpt_control%n_reortho 364 CALL reorthogonalize(b(k, :), b, Sb, R, k - 1) ! R is temp 365 END DO 366 CALL normalize(b(k, :), R, matrix_s) ! R is temp 367 DO spin = 1, nspins 368 CALL cp_fm_to_fm(b(k, spin)%matrix, X(spin)%matrix) 369 END DO 370 CALL apply_op(X, Ab(k, :), p_env, qs_env, & 371 dft_control%tddfpt_control%do_kernel) 372 CALL p_postortho(p_env, qs_env, Ab(k, :)) 373 DO spin = 1, nspins 374 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, & 375 b(k, spin)%matrix, & 376 Sb(k, spin)%matrix, & 377 p_env%n_mo(spin)) 378 END DO 379 END DO 380 381 !--------------------------------------------! 382 ! deallocate memory for the reduced matrices ! 383 !--------------------------------------------! 384 IF (ASSOCIATED(Atilde%matrix)) CALL cp_fm_release(Atilde%matrix) 385 IF (ASSOCIATED(Us%matrix)) CALL cp_fm_release(Us%matrix) 386 IF (ASSOCIATED(tilde_fm_struct)) CALL cp_fm_struct_release(tilde_fm_struct) 387 388 !------------------------------------------! 389 ! allocate memory for the reduced matrices ! 390 !------------------------------------------! 391 CALL cp_fm_struct_create(tilde_fm_struct, para_env, blacs_env, k, k) 392 CALL cp_fm_create(Atilde%matrix, & 393 tilde_fm_struct, & 394 routineP//"Atilde") 395 CALL cp_fm_create(Us%matrix, & 396 tilde_fm_struct, & 397 routineP//"Us") 398 399 !---------------------------------------! 400 ! calc the matrix Atilde = transp(b)*Ab ! 401 !---------------------------------------! 402 DO i = 1, k 403 DO j = 1, k 404 Atilde_ij = 0.0_dp 405 DO spin = 1, nspins 406 CALL cp_fm_trace(b(i, spin)%matrix, Ab(j, spin)%matrix, tmp) 407 Atilde_ij = Atilde_ij + tmp 408 END DO 409 CALL cp_fm_set_element(Atilde%matrix, i, j, Atilde_ij) 410 END DO 411 END DO 412 413 !--------------------! 414 ! diagonalize Atilde ! 415 !--------------------! 416 evals_tmp(:) = evals(:, iter) 417 CALL cp_fm_syevd(Atilde%matrix, Us%matrix, evals_tmp(:)) 418 evals(:, iter) = evals_tmp(:) 419 420 !-------------------! 421 ! check convergence ! 422 !-------------------! 423 evals_difference = 1.0_dp 424 IF (iter /= 1) THEN 425 426 evals_difference(:) = ABS((evals(1:n_ev, iter - 1) - evals(1:n_ev, iter))) 427 ! For debugging 428 IF (output_unit > 0) THEN 429 WRITE (output_unit, *) 430 DO i = 1, n_ev 431 WRITE (output_unit, '(2X,F10.7,T69,ES11.4)') evals(i, iter)*evolt, evals_difference(i) 432 END DO 433 WRITE (output_unit, *) 434 END IF 435 436 convergence = MAXVAL(evals_difference) 437 IF (output_unit > 0) WRITE (output_unit, '(2X,I4,T69,ES11.4)') k, convergence 438 439 IF (convergence < tddfpt_control%tolerance) THEN 440 res = .TRUE. 441 EXIT iteration 442 END IF 443 END IF 444 445 IF (mode == 'L') THEN 446 n_kv = 1 447 ELSE 448 must_improve = 0 449 DO i = 1, n_ev 450 IF (evals_difference(i) > tddfpt_control%tolerance) must_improve(i) = 1 451 END DO 452!! Set must_improve to 1 if all the vectors should 453!! be updated in one iteration. 454!! must_improve = 1 455 n_kv = SUM(must_improve) 456 j = 1 457 DO i = 1, n_ev 458 IF (must_improve(i) == 1) THEN 459 must_improve(j) = i 460 j = j + 1 461 END IF 462 END DO 463 END IF 464 465 IF (k + n_kv > max_kv) EXIT iteration 466 467 iter = iter + 1 468 469 END DO iteration 470 471 IF (PRESENT(out_evals)) THEN 472 out_evals(1:n_ev) = evals(1:n_ev, iter) 473 END IF 474 475 DO spin = 1, nspins 476 DO j = 1, n_ev 477 CALL cp_fm_set_all(t_env%evecs(j, spin)%matrix, 0.0_dp) 478 DO i = 1, k 479 CALL cp_fm_get_element(Us%matrix, i, j, tmp) 480 CALL cp_fm_scale_and_add(1.0_dp, t_env%evecs(j, spin)%matrix, & 481 tmp, b(i, spin)%matrix) 482 END DO 483 END DO 484 END DO 485 486 DO spin = 1, nspins 487 DO i = 1, max_kv 488 IF (ASSOCIATED(b(i, spin)%matrix)) & 489 CALL cp_fm_release(b(i, spin)%matrix) 490 IF (ASSOCIATED(Ab(i, spin)%matrix)) & 491 CALL cp_fm_release(Ab(i, spin)%matrix) 492 IF (ASSOCIATED(Sb(i, spin)%matrix)) & 493 CALL cp_fm_release(Sb(i, spin)%matrix) 494 END DO 495 END DO 496 497 !----------! 498 ! clean up ! 499 !----------! 500 IF (ASSOCIATED(Atilde%matrix)) CALL cp_fm_release(Atilde%matrix) 501 IF (ASSOCIATED(Us%matrix)) CALL cp_fm_release(Us%matrix) 502 IF (ASSOCIATED(tilde_fm_struct)) CALL cp_fm_struct_release(tilde_fm_struct) 503 CALL fm_pools_give_back_fm_vect(ao_mo_fm_pools, X) 504 CALL fm_pools_give_back_fm_vect(ao_mo_fm_pools, R) 505 NULLIFY (X, R) 506 DO spin = 1, nspins 507 CALL cp_fm_struct_release(kv_fm_struct(spin)%struct) 508 DO i = 1, max_kv 509 IF (ASSOCIATED(b(i, spin)%matrix)) & 510 CALL cp_fm_release(b(i, spin)%matrix) 511 IF (ASSOCIATED(Ab(i, spin)%matrix)) & 512 CALL cp_fm_release(Ab(i, spin)%matrix) 513 IF (ASSOCIATED(Sb(i, spin)%matrix)) & 514 CALL cp_fm_release(Sb(i, spin)%matrix) 515 END DO 516 END DO 517 DEALLOCATE (b, Ab, Sb, evals, evals_tmp, evals_difference, must_improve, kv_fm_struct) 518 519 CALL timestop(handle) 520 521 END FUNCTION iterative_solver 522 523 ! X : the vector on which to apply the op 524 ! R : the result 525 ! t_env : td-dft environment (mainly control information) 526 ! p_env : perturbation environment (variables) 527 ! both of these carry info for the tddfpt calculation 528 ! qs_env : info about a quickstep ground state calculation 529 530! ************************************************************************************************** 531!> \brief ... 532!> \param X ... 533!> \param R ... 534!> \param p_env ... 535!> \param qs_env ... 536!> \param do_kernel ... 537! ************************************************************************************************** 538 SUBROUTINE apply_op(X, R, p_env, qs_env, do_kernel) 539 540 TYPE(cp_fm_p_type), DIMENSION(:) :: X, R 541 TYPE(qs_p_env_type), POINTER :: p_env 542 TYPE(qs_environment_type), POINTER :: qs_env 543 LOGICAL, INTENT(IN) :: do_kernel 544 545 CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_op', routineP = moduleN//':'//routineN 546 547 INTEGER :: handle, nspins, spin 548 INTEGER, SAVE :: counter = 0 549 TYPE(dft_control_type), POINTER :: dft_control 550 551 NULLIFY (dft_control) 552 553 CALL timeset(routineN, handle) 554 555 counter = counter + 1 556 CALL get_qs_env(qs_env, dft_control=dft_control) 557 nspins = dft_control%nspins 558 559 !------------! 560 ! R = HX-SXL ! 561 !------------! 562 CALL p_op_l1(p_env, qs_env, X, R) ! acts on both spins, result in R 563 564 !-----------------! 565 ! calc P1 and ! 566 ! R = R + K(P1)*C ! 567 !-----------------! 568 IF (do_kernel) THEN 569 DO spin = 1, nspins 570 CALL dbcsr_set(p_env%p1(spin)%matrix, 0.0_dp) ! optimize? 571 CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(spin)%matrix, & 572 matrix_v=p_env%psi0d(spin)%matrix, & 573 matrix_g=X(spin)%matrix, & 574 ncol=p_env%n_mo(spin), & 575 alpha=0.5_dp) 576 CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(spin)%matrix, & 577 matrix_v=X(spin)%matrix, & 578 matrix_g=p_env%psi0d(spin)%matrix, & 579 ncol=p_env%n_mo(spin), & 580 alpha=0.5_dp) 581 END DO 582 DO spin = 1, nspins 583 CALL cp_fm_set_all(X(spin)%matrix, 0.0_dp) 584 END DO 585 CALL p_op_l2(p_env, qs_env, p_env%p1, X, & 586 alpha=1.0_dp, beta=0.0_dp) ! X = beta*X + alpha*K(P1)*C 587 DO spin = 1, nspins 588 CALL cp_fm_scale_and_add(1.0_dp, R(spin)%matrix, & 589 1.0_dp, X(spin)%matrix) ! add X to R 590 END DO 591 END IF 592 593 CALL timestop(handle) 594 595 END SUBROUTINE apply_op 596 597! ************************************************************************************************** 598!> \brief ... 599!> \param vectors ... 600!> \param vectors_name ... 601!> \param startv ... 602!> \param n_v ... 603!> \param nspins ... 604!> \param fm_struct ... 605! ************************************************************************************************** 606 SUBROUTINE allocate_krylov_vectors(vectors, vectors_name, & 607 startv, n_v, nspins, fm_struct) 608 609 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: vectors 610 CHARACTER(LEN=*), INTENT(IN) :: vectors_name 611 INTEGER, INTENT(IN) :: startv, n_v, nspins 612 TYPE(cp_fm_struct_p_type), DIMENSION(:), POINTER :: fm_struct 613 614 CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_krylov_vectors', & 615 routineP = moduleN//':'//routineN 616 617 CHARACTER(LEN=default_string_length) :: mat_name 618 INTEGER :: index, spin 619 620 DO spin = 1, nspins 621 DO index = startv, startv + n_v - 1 622 NULLIFY (vectors(index, spin)%matrix) 623 mat_name = routineP//vectors_name//TRIM(cp_to_string(index)) & 624 //","//TRIM(cp_to_string(spin)) 625 CALL cp_fm_create(vectors(index, spin)%matrix, & 626 fm_struct(spin)%struct, mat_name) 627 IF (.NOT. ASSOCIATED(vectors(index, spin)%matrix)) & 628 CPABORT("Could not allocate "//TRIM(mat_name)//".") 629 END DO 630 END DO 631 632 END SUBROUTINE allocate_krylov_vectors 633 634END MODULE qs_tddfpt_eigensolver 635