1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6MODULE qs_tddfpt2_eigensolver 7 USE cp_blacs_env, ONLY: cp_blacs_env_type 8 USE cp_control_types, ONLY: tddfpt2_control_type 9 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply 10 USE cp_fm_basic_linalg, ONLY: cp_fm_contracted_trace,& 11 cp_fm_scale,& 12 cp_fm_scale_and_add,& 13 cp_fm_trace 14 USE cp_fm_diag, ONLY: choose_eigv_solver 15 USE cp_fm_pool_types, ONLY: fm_pool_create_fm,& 16 fm_pool_give_back_fm 17 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 18 cp_fm_struct_release,& 19 cp_fm_struct_type 20 USE cp_fm_types, ONLY: & 21 cp_fm_copy_general, cp_fm_create, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_maxabsval, & 22 cp_fm_p_type, cp_fm_release, cp_fm_set_all, cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_type 23 USE cp_gemm_interface, ONLY: cp_gemm 24 USE cp_log_handling, ONLY: cp_logger_type 25 USE cp_output_handling, ONLY: cp_iterate 26 USE cp_para_types, ONLY: cp_para_env_type 27 USE dbcsr_api, ONLY: dbcsr_get_info,& 28 dbcsr_p_type,& 29 dbcsr_type 30 USE input_constants, ONLY: tddfpt_kernel_full,& 31 tddfpt_kernel_stda 32 USE input_section_types, ONLY: section_vals_type 33 USE kinds, ONLY: dp,& 34 int_8 35 USE machine, ONLY: m_flush,& 36 m_walltime 37 USE memory_utilities, ONLY: reallocate 38 USE physcon, ONLY: evolt 39 USE qs_environment_types, ONLY: get_qs_env,& 40 qs_environment_type 41 USE qs_scf_methods, ONLY: eigensolver 42 USE qs_tddfpt2_fhxc, ONLY: fhxc_kernel,& 43 stda_kernel 44 USE qs_tddfpt2_operators, ONLY: tddfpt_apply_energy_diff,& 45 tddfpt_apply_hfx 46 USE qs_tddfpt2_restart, ONLY: tddfpt_write_restart 47 USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type 48 USE qs_tddfpt2_types, ONLY: full_kernel_env_type,& 49 tddfpt_ground_state_mos,& 50 tddfpt_kernel_env_type,& 51 tddfpt_work_matrices 52 USE qs_tddfpt2_utils, ONLY: tddfpt_total_number_of_states 53#include "./base/base_uses.f90" 54 55 IMPLICIT NONE 56 57 PRIVATE 58 59 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_eigensolver' 60 61 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. 62 ! number of first derivative components (3: d/dx, d/dy, d/dz) 63 INTEGER, PARAMETER, PRIVATE :: nderivs = 3 64 INTEGER, PARAMETER, PRIVATE :: maxspins = 2 65 66 PUBLIC :: tddfpt_davidson_solver, tddfpt_orthogonalize_psi1_psi0, & 67 tddfpt_orthonormalize_psi1_psi1 68 69! ************************************************************************************************** 70 71CONTAINS 72 73! ************************************************************************************************** 74!> \brief Make TDDFPT trial vectors orthogonal to all occupied molecular orbitals. 75!> \param evects trial vectors distributed across all processors (modified on exit) 76!> \param S_C0_C0T matrix product S * C_0 * C_0^T, where C_0 is the ground state 77!> wave function for each spin expressed in atomic basis set, 78!> and S is the corresponding overlap matrix 79!> \par History 80!> * 05.2016 created [Sergey Chulkov] 81!> * 05.2019 use a temporary work matrix [JHU] 82!> \note Based on the subroutine p_preortho() which was created by Thomas Chassaing on 09.2002. 83!> Should be useless when ground state MOs are computed with extremely high accuracy, 84!> as all virtual orbitals are already orthogonal to the occupied ones by design. 85!> However, when the norm of residual vectors is relatively small (e.g. less then SCF_EPS), 86!> new Krylov's vectors seem to be random and should be orthogonalised even with respect to 87!> the occupied MOs. 88! ************************************************************************************************** 89 SUBROUTINE tddfpt_orthogonalize_psi1_psi0(evects, S_C0_C0T) 90 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: evects 91 TYPE(cp_fm_p_type), DIMENSION(:), INTENT(in) :: S_C0_C0T 92 93 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_orthogonalize_psi1_psi0', & 94 routineP = moduleN//':'//routineN 95 96 INTEGER :: handle, ispin, ivect, nactive, nao, & 97 nspins, nvects 98 TYPE(cp_fm_struct_type), POINTER :: matrix_struct 99 TYPE(cp_fm_type), POINTER :: evortho 100 101 CALL timeset(routineN, handle) 102 103 nspins = SIZE(evects, 1) 104 nvects = SIZE(evects, 2) 105 106 IF (nvects > 0) THEN 107 DO ispin = 1, nspins 108 CALL cp_fm_get_info(matrix=evects(ispin, 1)%matrix, matrix_struct=matrix_struct, & 109 nrow_global=nao, ncol_global=nactive) 110 NULLIFY (evortho) 111 CALL cp_fm_create(evortho, matrix_struct) 112 DO ivect = 1, nvects 113 ! evortho: C0 * C0^T * S * C1 == (S * C0 * C0^T)^T * C1 114 CALL cp_gemm('T', 'N', nao, nactive, nao, 1.0_dp, S_C0_C0T(ispin)%matrix, & 115 evects(ispin, ivect)%matrix, 0.0_dp, evortho) 116 CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, ivect)%matrix, -1.0_dp, evortho) 117 END DO 118 CALL cp_fm_release(evortho) 119 END DO 120 END IF 121 122 CALL timestop(handle) 123 124 END SUBROUTINE tddfpt_orthogonalize_psi1_psi0 125 126! ************************************************************************************************** 127!> \brief Check that orthogonalised TDDFPT trial vectors remain orthogonal to 128!> occupied molecular orbitals. 129!> \param evects trial vectors 130!> \param S_C0 matrix product S * C_0, where C_0 is the ground state wave function 131!> for each spin in atomic basis set, and S is the corresponding overlap matrix 132!> \param max_norm the largest possible overlap between the ground state and 133!> excited state wave functions 134!> \return true if trial vectors are non-orthogonal to occupied molecular orbitals 135!> \par History 136!> * 07.2016 created [Sergey Chulkov] 137!> * 05.2019 use temporary work matrices [JHU] 138! ************************************************************************************************** 139 FUNCTION tddfpt_is_nonorthogonal_psi1_psi0(evects, S_C0, max_norm) & 140 RESULT(is_nonortho) 141 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: evects 142 TYPE(cp_fm_p_type), DIMENSION(:), INTENT(in) :: S_C0 143 REAL(kind=dp), INTENT(in) :: max_norm 144 LOGICAL :: is_nonortho 145 146 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_is_nonorthogonal_psi1_psi0', & 147 routineP = moduleN//':'//routineN 148 149 INTEGER :: handle, ispin, ivect, nactive, nao, & 150 nocc, nspins, nvects 151 REAL(kind=dp) :: maxabs_val 152 TYPE(cp_fm_struct_type), POINTER :: matrix_struct, matrix_struct_tmp 153 TYPE(cp_fm_type), POINTER :: aortho 154 155 CALL timeset(routineN, handle) 156 157 nspins = SIZE(evects, 1) 158 nvects = SIZE(evects, 2) 159 160 is_nonortho = .FALSE. 161 162 loop: DO ispin = 1, nspins 163 CALL cp_fm_get_info(matrix=S_C0(ispin)%matrix, ncol_global=nocc) 164 CALL cp_fm_get_info(matrix=evects(ispin, 1)%matrix, matrix_struct=matrix_struct, & 165 nrow_global=nao, ncol_global=nactive) 166 CALL cp_fm_struct_create(matrix_struct_tmp, nrow_global=nocc, & 167 ncol_global=nactive, template_fmstruct=matrix_struct) 168 NULLIFY (aortho) 169 CALL cp_fm_create(aortho, matrix_struct_tmp) 170 CALL cp_fm_struct_release(matrix_struct_tmp) 171 DO ivect = 1, nvects 172 ! aortho = S_0^T * S * C_1 173 CALL cp_gemm('T', 'N', nocc, nactive, nao, 1.0_dp, S_C0(ispin)%matrix, & 174 evects(ispin, ivect)%matrix, 0.0_dp, aortho) 175 CALL cp_fm_maxabsval(aortho, maxabs_val) 176 is_nonortho = maxabs_val > max_norm 177 IF (is_nonortho) THEN 178 CALL cp_fm_release(aortho) 179 EXIT loop 180 END IF 181 END DO 182 CALL cp_fm_release(aortho) 183 END DO loop 184 185 CALL timestop(handle) 186 187 END FUNCTION tddfpt_is_nonorthogonal_psi1_psi0 188 189! ************************************************************************************************** 190!> \brief Make new TDDFPT trial vectors orthonormal to all previous trial vectors. 191!> \param evects trial vectors (modified on exit) 192!> \param nvects_new number of new trial vectors to orthogonalise 193!> \param S_evects set of matrices to store matrix product S * evects (modified on exit) 194!> \param matrix_s overlap matrix 195!> \par History 196!> * 05.2016 created [Sergey Chulkov] 197!> * 02.2017 caching the matrix product S * evects [Sergey Chulkov] 198!> \note \parblock 199!> Based on the subroutines reorthogonalize() and normalize() which were originally created 200!> by Thomas Chassaing on 03.2003. 201!> 202!> In order to orthogonalise a trial vector C3 = evects(:,3) with respect to previously 203!> orthogonalised vectors C1 = evects(:,1) and C2 = evects(:,2) we need to compute the 204!> quantity C3'' using the following formulae: 205!> C3' = C3 - Tr(C3^T * S * C1) * C1, 206!> C3'' = C3' - Tr(C3'^T * S * C2) * C2, 207!> which can be expanded as: 208!> C3'' = C3 - Tr(C3^T * S * C1) * C1 - Tr(C3^T * S * C2) * C2 + 209!> Tr(C3^T * S * C1) * Tr(C2^T * S * C1) * C2 . 210!> In case of unlimited float-point precision, the last term in above expression is exactly 0, 211!> due to orthogonality condition between C1 and C2. In this case the expression could be 212!> simplified as (taking into account the identity: Tr(A * S * B) = Tr(B * S * A)): 213!> C3'' = C3 - Tr(C1^T * S * C3) * C1 - Tr(C2^T * S * C3) * C2 , 214!> which means we do not need the variable S_evects to keep the matrix products S * Ci . 215!> 216!> In reality, however, we deal with limited float-point precision arithmetic meaning that 217!> the trace Tr(C2^T * S * C1) is close to 0 but does not equal to 0 exactly. The term 218!> Tr(C3^T * S * C1) * Tr(C2^T * S * C1) * C2 219!> can not be ignored anymore. Ignorance of this term will lead to numerical instability 220!> when the trace Tr(C3^T * S * C1) is large enough. 221!> \endparblock 222! ************************************************************************************************** 223 SUBROUTINE tddfpt_orthonormalize_psi1_psi1(evects, nvects_new, S_evects, matrix_s) 224 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: evects 225 INTEGER, INTENT(in) :: nvects_new 226 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: S_evects 227 TYPE(dbcsr_type), POINTER :: matrix_s 228 229 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_orthonormalize_psi1_psi1', & 230 routineP = moduleN//':'//routineN 231 232 INTEGER :: handle, ispin, ivect, jvect, nspins, & 233 nvects_old, nvects_total 234 INTEGER, DIMENSION(maxspins) :: nactive 235 REAL(kind=dp) :: norm 236 REAL(kind=dp), DIMENSION(maxspins) :: weights 237 238 CALL timeset(routineN, handle) 239 240 nspins = SIZE(evects, 1) 241 nvects_total = SIZE(evects, 2) 242 nvects_old = nvects_total - nvects_new 243 244 IF (debug_this_module) THEN 245 CPASSERT(SIZE(S_evects, 1) == nspins) 246 CPASSERT(SIZE(S_evects, 2) == nvects_total) 247 CPASSERT(nvects_old >= 0) 248 END IF 249 250 DO ispin = 1, nspins 251 CALL cp_fm_get_info(matrix=evects(ispin, 1)%matrix, ncol_global=nactive(ispin)) 252 END DO 253 254 DO jvect = nvects_old + 1, nvects_total 255 ! <psi1_i | psi1_j> 256 DO ivect = 1, jvect - 1 257 CALL cp_fm_trace(evects(:, jvect), S_evects(:, ivect), weights(1:nspins), accurate=.FALSE.) 258 norm = SUM(weights(1:nspins)) 259 260 DO ispin = 1, nspins 261 CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, jvect)%matrix, -norm, evects(ispin, ivect)%matrix) 262 END DO 263 END DO 264 265 ! <psi1_j | psi1_j> 266 DO ispin = 1, nspins 267 CALL cp_dbcsr_sm_fm_multiply(matrix_s, evects(ispin, jvect)%matrix, S_evects(ispin, jvect)%matrix, & 268 ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp) 269 END DO 270 271 CALL cp_fm_trace(evects(:, jvect), S_evects(:, jvect), weights(1:nspins), accurate=.FALSE.) 272 273 norm = SUM(weights(1:nspins)) 274 norm = 1.0_dp/SQRT(norm) 275 276 DO ispin = 1, nspins 277 CALL cp_fm_scale(norm, evects(ispin, jvect)%matrix) 278 CALL cp_fm_scale(norm, S_evects(ispin, jvect)%matrix) 279 END DO 280 END DO 281 282 CALL timestop(handle) 283 284 END SUBROUTINE tddfpt_orthonormalize_psi1_psi1 285 286! ************************************************************************************************** 287!> \brief Compute action matrix-vector products. 288!> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit) 289!> \param evects TDDFPT trial vectors 290!> \param S_evects cached matrix product S * evects where S is the overlap matrix 291!> in primary basis set 292!> \param gs_mos molecular orbitals optimised for the ground state 293!> \param tddfpt_control control section for tddfpt 294!> \param do_hfx flag that activates computation of exact-exchange terms 295!> \param matrix_ks Kohn-Sham matrix 296!> \param qs_env Quickstep environment 297!> \param kernel_env kernel environment 298!> \param sub_env parallel (sub)group environment 299!> \param work_matrices collection of work matrices (modified on exit) 300!> \par History 301!> * 06.2016 created [Sergey Chulkov] 302!> * 03.2017 refactored [Sergey Chulkov] 303! ************************************************************************************************** 304 SUBROUTINE tddfpt_compute_Aop_evects(Aop_evects, evects, S_evects, gs_mos, tddfpt_control, & 305 do_hfx, matrix_ks, qs_env, kernel_env, & 306 sub_env, work_matrices) 307 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: Aop_evects, evects, S_evects 308 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 309 INTENT(in) :: gs_mos 310 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control 311 LOGICAL, INTENT(in) :: do_hfx 312 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_ks 313 TYPE(qs_environment_type), POINTER :: qs_env 314 TYPE(tddfpt_kernel_env_type), INTENT(in) :: kernel_env 315 TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env 316 TYPE(tddfpt_work_matrices), INTENT(inout) :: work_matrices 317 318 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_compute_Aop_evects', & 319 routineP = moduleN//':'//routineN 320 321 INTEGER :: handle, ispin, ivect, nspins, nvects 322 INTEGER, DIMENSION(maxspins) :: nmo_occ 323 LOGICAL :: do_admm, is_rks_triplets, is_stda 324 TYPE(cp_para_env_type), POINTER :: para_env 325 TYPE(full_kernel_env_type), POINTER :: full_kernel_env, kernel_env_admm_aux 326 327 CALL timeset(routineN, handle) 328 329 full_kernel_env => kernel_env%full_kernel 330 331 nspins = SIZE(evects, 1) 332 nvects = SIZE(evects, 2) 333 do_admm = ASSOCIATED(sub_env%admm_A) 334 IF (do_admm) THEN 335 kernel_env_admm_aux => kernel_env%admm_kernel 336 ELSE 337 NULLIFY (kernel_env_admm_aux) 338 END IF 339 is_rks_triplets = tddfpt_control%rks_triplets 340 341 IF (debug_this_module) THEN 342 CPASSERT(nspins > 0) 343 CPASSERT(SIZE(Aop_evects, 1) == nspins) 344 CPASSERT(SIZE(Aop_evects, 2) == nvects) 345 CPASSERT(SIZE(S_evects, 1) == nspins) 346 CPASSERT(SIZE(S_evects, 2) == nvects) 347 CPASSERT(SIZE(gs_mos) == nspins) 348 END IF 349 350 DO ispin = 1, nspins 351 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ) 352 END DO 353 354 IF (nvects > 0) THEN 355 CALL cp_fm_get_info(evects(1, 1)%matrix, para_env=para_env) 356 IF (ALLOCATED(work_matrices%evects_sub)) THEN 357 DO ivect = 1, nvects 358 DO ispin = 1, nspins 359 CALL cp_fm_copy_general(evects(ispin, ivect)%matrix, & 360 work_matrices%evects_sub(ispin, ivect)%matrix, para_env) 361 END DO 362 END DO 363 END IF 364 365 IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN 366 ! full TDDFPT kernel 367 CALL fhxc_kernel(Aop_evects, evects, is_rks_triplets, do_hfx, qs_env, & 368 full_kernel_env, kernel_env_admm_aux, sub_env, work_matrices) 369 is_stda = .FALSE. 370 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN 371 ! sTDA kernel 372 CALL stda_kernel(Aop_evects, evects, is_rks_triplets, qs_env, tddfpt_control%stda_control, & 373 kernel_env%stda_kernel, sub_env, work_matrices) 374 is_stda = .TRUE. 375 ELSE 376 CPABORT("Kernel type not implemented") 377 END IF 378 379 IF (ALLOCATED(work_matrices%evects_sub)) THEN 380 DO ivect = 1, nvects 381 DO ispin = 1, nspins 382 CALL cp_fm_copy_general(work_matrices%Aop_evects_sub(ispin, ivect)%matrix, & 383 Aop_evects(ispin, ivect)%matrix, para_env) 384 END DO 385 END DO 386 END IF 387 388 ! orbital energy difference term 389 CALL tddfpt_apply_energy_diff(Aop_evects=Aop_evects, evects=evects, S_evects=S_evects, & 390 gs_mos=gs_mos, matrix_ks=matrix_ks) 391 392 IF (do_hfx .AND. .NOT. is_stda) THEN 393 CALL tddfpt_apply_hfx(Aop_evects=Aop_evects, evects=evects, gs_mos=gs_mos, do_admm=do_admm, & 394 qs_env=qs_env, work_rho_ia_ao=work_matrices%hfx_rho_ao, & 395 work_hmat=work_matrices%hfx_hmat, wfm_rho_orb=work_matrices%hfx_fm_ao_ao) 396 END IF 397 END IF 398 399 CALL timestop(handle) 400 401 END SUBROUTINE tddfpt_compute_Aop_evects 402 403! ************************************************************************************************** 404!> \brief Solve eigenproblem for the reduced action matrix and find new Ritz eigenvectors and 405!> eigenvalues. 406!> \param ritz_vects Ritz eigenvectors (initialised on exit) 407!> \param Aop_ritz approximate action of TDDFPT operator on Ritz vectors (initialised on exit) 408!> \param evals Ritz eigenvalues (initialised on exit) 409!> \param krylov_vects Krylov's vectors 410!> \param Aop_krylov action of TDDFPT operator on Krylov's vectors 411!> \param Atilde ... 412!> \param nkvo ... 413!> \param nkvn ... 414!> \par History 415!> * 06.2016 created [Sergey Chulkov] 416!> * 03.2017 altered prototype, OpenMP parallelisation [Sergey Chulkov] 417! ************************************************************************************************** 418 SUBROUTINE tddfpt_compute_ritz_vects(ritz_vects, Aop_ritz, evals, krylov_vects, Aop_krylov, & 419 Atilde, nkvo, nkvn) 420 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: ritz_vects, Aop_ritz 421 REAL(kind=dp), DIMENSION(:), INTENT(out) :: evals 422 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: krylov_vects, Aop_krylov 423 REAL(kind=dp), DIMENSION(:, :), POINTER :: Atilde 424 INTEGER, INTENT(IN) :: nkvo, nkvn 425 426 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_compute_ritz_vects', & 427 routineP = moduleN//':'//routineN 428 429 INTEGER :: handle, ikv, irv, ispin, nkvs, nrvs, & 430 nspins 431 REAL(kind=dp) :: act 432 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: at12, at21, at22, evects_Atilde 433 TYPE(cp_blacs_env_type), POINTER :: blacs_env_global 434 TYPE(cp_fm_struct_type), POINTER :: fm_struct 435 TYPE(cp_fm_type), POINTER :: Atilde_fm, evects_Atilde_fm 436 437 CALL timeset(routineN, handle) 438 439 nspins = SIZE(krylov_vects, 1) 440 nkvs = SIZE(krylov_vects, 2) 441 nrvs = SIZE(ritz_vects, 2) 442 CPASSERT(nkvs == nkvo + nkvn) 443 444 CALL cp_fm_get_info(krylov_vects(1, 1)%matrix, context=blacs_env_global) 445 446 CALL cp_fm_struct_create(fm_struct, nrow_global=nkvs, ncol_global=nkvs, context=blacs_env_global) 447 CALL cp_fm_create(Atilde_fm, fm_struct) 448 CALL cp_fm_create(evects_Atilde_fm, fm_struct) 449 CALL cp_fm_struct_release(fm_struct) 450 451 ! *** compute upper-diagonal reduced action matrix *** 452 CALL reallocate(Atilde, 1, nkvs, 1, nkvs) 453 ! TO DO: the subroutine 'cp_fm_contracted_trace' will compute all elements of 454 ! the matrix 'Atilde', however only upper-triangular elements are actually needed 455 ! 456 IF (nkvo == 0) THEN 457 CALL cp_fm_contracted_trace(Aop_krylov(:, 1:nkvs), krylov_vects(:, 1:nkvs), & 458 Atilde(1:nkvs, 1:nkvs), accurate=.FALSE.) 459 ELSE 460 ALLOCATE (at12(nkvn, nkvo), at21(nkvo, nkvn), at22(nkvn, nkvn)) 461 CALL cp_fm_contracted_trace(Aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, 1:nkvo), & 462 at12, accurate=.FALSE.) 463 Atilde(nkvo + 1:nkvs, 1:nkvo) = at12(1:nkvn, 1:nkvo) 464 CALL cp_fm_contracted_trace(Aop_krylov(:, 1:nkvo), krylov_vects(:, nkvo + 1:nkvs), & 465 at21, accurate=.FALSE.) 466 Atilde(1:nkvo, nkvo + 1:nkvs) = at21(1:nkvo, 1:nkvn) 467 CALL cp_fm_contracted_trace(Aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, nkvo + 1:nkvs), & 468 at22, accurate=.FALSE.) 469 Atilde(nkvo + 1:nkvs, nkvo + 1:nkvs) = at22(1:nkvn, 1:nkvn) 470 DEALLOCATE (at12, at21, at22) 471 END IF 472 Atilde = 0.5_dp*(Atilde + TRANSPOSE(Atilde)) 473 CALL cp_fm_set_submatrix(Atilde_fm, Atilde) 474 475 ! *** solve an eigenproblem for the reduced matrix *** 476 CALL choose_eigv_solver(Atilde_fm, evects_Atilde_fm, evals(1:nkvs)) 477 478 ALLOCATE (evects_Atilde(nkvs, nrvs)) 479 CALL cp_fm_get_submatrix(evects_Atilde_fm, evects_Atilde, start_row=1, start_col=1, n_rows=nkvs, n_cols=nrvs) 480 CALL cp_fm_release(evects_Atilde_fm) 481 CALL cp_fm_release(Atilde_fm) 482 483!$OMP PARALLEL DO DEFAULT(NONE), & 484!$OMP PRIVATE(act, ikv, irv, ispin), & 485!$OMP SHARED(Aop_krylov, Aop_ritz, krylov_vects, evects_Atilde, nkvs, nrvs, nspins, ritz_vects) 486 DO irv = 1, nrvs 487 DO ispin = 1, nspins 488 CALL cp_fm_set_all(ritz_vects(ispin, irv)%matrix, 0.0_dp) 489 CALL cp_fm_set_all(Aop_ritz(ispin, irv)%matrix, 0.0_dp) 490 END DO 491 492 DO ikv = 1, nkvs 493 act = evects_Atilde(ikv, irv) 494 DO ispin = 1, nspins 495 CALL cp_fm_scale_and_add(1.0_dp, ritz_vects(ispin, irv)%matrix, & 496 act, krylov_vects(ispin, ikv)%matrix) 497 CALL cp_fm_scale_and_add(1.0_dp, Aop_ritz(ispin, irv)%matrix, & 498 act, Aop_krylov(ispin, ikv)%matrix) 499 END DO 500 END DO 501 END DO 502!$OMP END PARALLEL DO 503 504 DEALLOCATE (evects_Atilde) 505 506 CALL timestop(handle) 507 508 END SUBROUTINE tddfpt_compute_ritz_vects 509 510! ************************************************************************************************** 511!> \brief Expand Krylov space by computing residual vectors. 512!> \param residual_vects residual vectors (modified on exit) 513!> \param evals Ritz eigenvalues (modified on exit) 514!> \param ritz_vects Ritz eigenvectors 515!> \param Aop_ritz approximate action of TDDFPT operator on Ritz vectors 516!> \param gs_mos molecular orbitals optimised for the ground state 517!> \param matrix_s overlap matrix 518!> \par History 519!> * 06.2016 created [Sergey Chulkov] 520!> * 03.2017 refactored to achieve significant performance gain [Sergey Chulkov] 521! ************************************************************************************************** 522 SUBROUTINE tddfpt_compute_residual_vects(residual_vects, evals, ritz_vects, Aop_ritz, gs_mos, & 523 matrix_s) 524 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: residual_vects 525 REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals 526 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: ritz_vects, Aop_ritz 527 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 528 INTENT(in) :: gs_mos 529 TYPE(dbcsr_type), POINTER :: matrix_s 530 531 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_compute_residual_vects', & 532 routineP = moduleN//':'//routineN 533 REAL(kind=dp), PARAMETER :: eref_scale = 0.99_dp, threshold = 16.0_dp*EPSILON(1.0_dp) 534 535 INTEGER :: handle, icol_local, irow_local, irv, & 536 ispin, nao, ncols_local, nrows_local, & 537 nrvs, nspins 538 INTEGER, DIMENSION(:), POINTER :: col_indices_local, row_indices_local 539 INTEGER, DIMENSION(maxspins) :: nactive, nmo_virt 540 REAL(kind=dp) :: e_occ_plus_lambda, eref, lambda 541 REAL(kind=dp), DIMENSION(:, :), POINTER :: weights_ldata 542 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: awork, vomat 543 TYPE(cp_fm_struct_type), POINTER :: ao_mo_struct, virt_mo_struct 544 545 CALL timeset(routineN, handle) 546 547 nspins = SIZE(residual_vects, 1) 548 nrvs = SIZE(residual_vects, 2) 549 550 IF (nrvs > 0) THEN 551 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao) 552 ALLOCATE (awork(nspins), vomat(nspins)) 553 DO ispin = 1, nspins 554 nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt) 555 ! 556 CALL cp_fm_get_info(matrix=ritz_vects(ispin, 1)%matrix, matrix_struct=ao_mo_struct, & 557 ncol_global=nactive(ispin)) 558 NULLIFY (awork(ispin)%matrix) 559 CALL cp_fm_create(awork(ispin)%matrix, ao_mo_struct) 560 ! 561 CALL cp_fm_struct_create(virt_mo_struct, nrow_global=nmo_virt(ispin), & 562 ncol_global=nactive(ispin), template_fmstruct=ao_mo_struct) 563 NULLIFY (vomat(ispin)%matrix) 564 CALL cp_fm_create(vomat(ispin)%matrix, virt_mo_struct) 565 CALL cp_fm_struct_release(virt_mo_struct) 566 END DO 567 568 ! *** actually compute residual vectors *** 569 DO irv = 1, nrvs 570 lambda = evals(irv) 571 572 DO ispin = 1, nspins 573 CALL cp_fm_get_info(vomat(ispin)%matrix, nrow_local=nrows_local, & 574 ncol_local=ncols_local, row_indices=row_indices_local, & 575 col_indices=col_indices_local, local_data=weights_ldata) 576 577 ! awork := Ab(ispin, irv) - evals(irv) b(ispin, irv), where 'b' is a Ritz vector 578 CALL cp_dbcsr_sm_fm_multiply(matrix_s, ritz_vects(ispin, irv)%matrix, awork(ispin)%matrix, & 579 ncol=nactive(ispin), alpha=-lambda, beta=0.0_dp) 580 CALL cp_fm_scale_and_add(1.0_dp, awork(ispin)%matrix, 1.0_dp, Aop_ritz(ispin, irv)%matrix) 581 ! 582 CALL cp_gemm('T', 'N', nmo_virt(ispin), nactive(ispin), nao, 1.0_dp, gs_mos(ispin)%mos_virt, & 583 awork(ispin)%matrix, 0.0_dp, vomat(ispin)%matrix) 584 585 DO icol_local = 1, ncols_local 586 e_occ_plus_lambda = gs_mos(ispin)%evals_occ(col_indices_local(icol_local)) + lambda 587 588 DO irow_local = 1, nrows_local 589 eref = gs_mos(ispin)%evals_virt(row_indices_local(irow_local)) - e_occ_plus_lambda 590 591 ! eref = e_virt - e_occ - lambda = e_virt - e_occ - (eref_scale*lambda + (1-eref_scale)*lambda); 592 ! eref_new = e_virt - e_occ - eref_scale*lambda = eref + (1 - eref_scale)*lambda 593 IF (ABS(eref) < threshold) & 594 eref = eref + (1.0_dp - eref_scale)*lambda 595 596 weights_ldata(irow_local, icol_local) = weights_ldata(irow_local, icol_local)/eref 597 END DO 598 END DO 599 600 CALL cp_gemm('N', 'N', nao, nactive(ispin), nmo_virt(ispin), 1.0_dp, gs_mos(ispin)%mos_virt, & 601 vomat(ispin)%matrix, 0.0_dp, residual_vects(ispin, irv)%matrix) 602 END DO 603 END DO 604 ! 605 DO ispin = 1, nspins 606 CALL cp_fm_release(awork(ispin)%matrix) 607 CALL cp_fm_release(vomat(ispin)%matrix) 608 END DO 609 DEALLOCATE (awork, vomat) 610 END IF 611 612 CALL timestop(handle) 613 614 END SUBROUTINE tddfpt_compute_residual_vects 615 616! ************************************************************************************************** 617!> \brief Perform Davidson iterations. 618!> \param evects TDDFPT trial vectors (modified on exit) 619!> \param evals TDDFPT eigenvalues (modified on exit) 620!> \param S_evects cached matrix product S * evects (modified on exit) 621!> \param gs_mos molecular orbitals optimised for the ground state 622!> \param do_hfx flag that activates computation of exact-exchange terms 623!> \param tddfpt_control TDDFPT control parameters 624!> \param matrix_ks Kohn-Sham matrix 625!> \param qs_env Quickstep environment 626!> \param kernel_env kernel environment 627!> \param sub_env parallel (sub)group environment 628!> \param logger CP2K logger 629!> \param iter_unit I/O unit to write basic iteration information 630!> \param energy_unit I/O unit to write detailed energy information 631!> \param tddfpt_print_section TDDFPT print input section (need to write TDDFPT restart files) 632!> \param work_matrices collection of work matrices (modified on exit) 633!> \return energy convergence achieved (in Hartree) 634!> \par History 635!> * 03.2017 code related to Davidson eigensolver has been moved here from the main subroutine 636!> tddfpt() [Sergey Chulkov] 637!> \note Based on the subroutines apply_op() and iterative_solver() originally created by 638!> Thomas Chassaing in 2002. 639! ************************************************************************************************** 640 FUNCTION tddfpt_davidson_solver(evects, evals, S_evects, gs_mos, do_hfx, tddfpt_control, & 641 matrix_ks, qs_env, kernel_env, & 642 sub_env, logger, iter_unit, energy_unit, & 643 tddfpt_print_section, work_matrices) RESULT(conv) 644 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(inout) :: evects 645 REAL(kind=dp), DIMENSION(:), INTENT(inout) :: evals 646 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(inout) :: S_evects 647 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 648 INTENT(in) :: gs_mos 649 LOGICAL, INTENT(in) :: do_hfx 650 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control 651 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks 652 TYPE(qs_environment_type), POINTER :: qs_env 653 TYPE(tddfpt_kernel_env_type), INTENT(in) :: kernel_env 654 TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env 655 TYPE(cp_logger_type), POINTER :: logger 656 INTEGER, INTENT(in) :: iter_unit, energy_unit 657 TYPE(section_vals_type), POINTER :: tddfpt_print_section 658 TYPE(tddfpt_work_matrices), INTENT(inout) :: work_matrices 659 REAL(kind=dp) :: conv 660 661 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_davidson_solver', & 662 routineP = moduleN//':'//routineN 663 664 INTEGER :: handle, ispin, istate, iter, & 665 max_krylov_vects, nspins, nstates, & 666 nstates_conv, nvects_exist, nvects_new 667 INTEGER(kind=int_8) :: nstates_total 668 LOGICAL :: is_nonortho 669 REAL(kind=dp) :: t1, t2 670 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals_last 671 REAL(kind=dp), DIMENSION(:, :), POINTER :: Atilde 672 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:, :) :: Aop_krylov, Aop_ritz, krylov_vects, & 673 S_krylov 674 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 675 676 CALL timeset(routineN, handle) 677 678 nspins = SIZE(gs_mos) 679 nstates = tddfpt_control%nstates 680 nstates_total = tddfpt_total_number_of_states(gs_mos) 681 682 IF (debug_this_module) THEN 683 CPASSERT(SIZE(evects, 1) == nspins) 684 CPASSERT(SIZE(evects, 2) == nstates) 685 CPASSERT(SIZE(evals) == nstates) 686 END IF 687 688 CALL get_qs_env(qs_env, matrix_s=matrix_s) 689 690 ! adjust the number of Krylov vectors 691 max_krylov_vects = tddfpt_control%nkvs 692 IF (max_krylov_vects < nstates) max_krylov_vects = nstates 693 IF (INT(max_krylov_vects, kind=int_8) > nstates_total) max_krylov_vects = INT(nstates_total) 694 695 ALLOCATE (Aop_ritz(nspins, nstates)) 696 DO istate = 1, nstates 697 DO ispin = 1, nspins 698 NULLIFY (Aop_ritz(ispin, istate)%matrix) 699 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, Aop_ritz(ispin, istate)%matrix) 700 END DO 701 END DO 702 703 ALLOCATE (evals_last(max_krylov_vects)) 704 ALLOCATE (Aop_krylov(nspins, max_krylov_vects), krylov_vects(nspins, max_krylov_vects), & 705 S_krylov(nspins, max_krylov_vects)) 706 DO istate = 1, max_krylov_vects 707 DO ispin = 1, nspins 708 NULLIFY (Aop_krylov(ispin, istate)%matrix, krylov_vects(ispin, istate)%matrix, & 709 S_krylov(ispin, istate)%matrix) 710 END DO 711 END DO 712 713 DO istate = 1, nstates 714 DO ispin = 1, nspins 715 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, krylov_vects(ispin, istate)%matrix) 716 CALL cp_fm_to_fm(evects(ispin, istate)%matrix, krylov_vects(ispin, istate)%matrix) 717 718 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, S_krylov(ispin, istate)%matrix) 719 CALL cp_fm_to_fm(S_evects(ispin, istate)%matrix, S_krylov(ispin, istate)%matrix) 720 721 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, Aop_krylov(ispin, istate)%matrix) 722 END DO 723 END DO 724 725 nvects_exist = 0 726 nvects_new = nstates 727 728 t1 = m_walltime() 729 730 ALLOCATE (Atilde(1, 1)) 731 732 DO 733 ! davidson iteration 734 CALL cp_iterate(logger%iter_info, iter_nr_out=iter) 735 736 CALL tddfpt_compute_Aop_evects(Aop_evects=Aop_krylov(:, nvects_exist + 1:nvects_exist + nvects_new), & 737 evects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), & 738 S_evects=S_krylov(:, nvects_exist + 1:nvects_exist + nvects_new), & 739 gs_mos=gs_mos, tddfpt_control=tddfpt_control, & 740 do_hfx=do_hfx, matrix_ks=matrix_ks, & 741 qs_env=qs_env, kernel_env=kernel_env, & 742 sub_env=sub_env, & 743 work_matrices=work_matrices) 744 745 CALL tddfpt_compute_ritz_vects(ritz_vects=evects, Aop_ritz=Aop_ritz, & 746 evals=evals_last(1:nvects_exist + nvects_new), & 747 krylov_vects=krylov_vects(:, 1:nvects_exist + nvects_new), & 748 Aop_krylov=Aop_krylov(:, 1:nvects_exist + nvects_new), & 749 Atilde=Atilde, nkvo=nvects_exist, nkvn=nvects_new) 750 751 CALL tddfpt_write_restart(evects=evects, evals=evals_last(1:nstates), gs_mos=gs_mos, & 752 logger=logger, tddfpt_print_section=tddfpt_print_section) 753 754 conv = MAXVAL(ABS(evals_last(1:nstates) - evals(1:nstates))) 755 756 nvects_exist = nvects_exist + nvects_new 757 IF (nvects_exist + nvects_new > max_krylov_vects) & 758 nvects_new = max_krylov_vects - nvects_exist 759 IF (iter >= tddfpt_control%niters) nvects_new = 0 760 761 IF (conv > tddfpt_control%conv .AND. nvects_new > 0) THEN 762 ! compute residual vectors for the next iteration 763 DO istate = 1, nvects_new 764 DO ispin = 1, nspins 765 NULLIFY (Aop_krylov(ispin, nvects_exist + istate)%matrix, & 766 krylov_vects(ispin, nvects_exist + istate)%matrix, & 767 S_krylov(ispin, nvects_exist + istate)%matrix) 768 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, & 769 krylov_vects(ispin, nvects_exist + istate)%matrix) 770 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, & 771 S_krylov(ispin, nvects_exist + istate)%matrix) 772 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, & 773 Aop_krylov(ispin, nvects_exist + istate)%matrix) 774 END DO 775 END DO 776 777 CALL tddfpt_compute_residual_vects(residual_vects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), & 778 evals=evals_last(1:nvects_new), & 779 ritz_vects=evects(:, 1:nvects_new), Aop_ritz=Aop_ritz(:, 1:nvects_new), & 780 gs_mos=gs_mos, matrix_s=matrix_s(1)%matrix) 781 782 CALL tddfpt_orthogonalize_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), & 783 work_matrices%S_C0_C0T) 784 785 CALL tddfpt_orthonormalize_psi1_psi1(krylov_vects(:, 1:nvects_exist + nvects_new), nvects_new, & 786 S_krylov(:, 1:nvects_exist + nvects_new), matrix_s(1)%matrix) 787 788 is_nonortho = tddfpt_is_nonorthogonal_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), & 789 work_matrices%S_C0, tddfpt_control%orthogonal_eps) 790 ELSE 791 ! convergence or the maximum number of Krylov vectors have been achieved 792 nvects_new = 0 793 is_nonortho = .FALSE. 794 END IF 795 796 t2 = m_walltime() 797 IF (energy_unit > 0) THEN 798 WRITE (energy_unit, '(/,4X,A,T14,A,T36,A)') "State", "Exc. energy (eV)", "Convergence (eV)" 799 DO istate = 1, nstates 800 WRITE (energy_unit, '(1X,I8,T12,F14.7,T38,ES11.4)') istate, & 801 evals_last(istate)*evolt, (evals_last(istate) - evals(istate))*evolt 802 END DO 803 WRITE (energy_unit, *) 804 CALL m_flush(energy_unit) 805 END IF 806 807 IF (iter_unit > 0) THEN 808 nstates_conv = 0 809 DO istate = 1, nstates 810 IF (ABS(evals_last(istate) - evals(istate)) <= tddfpt_control%conv) & 811 nstates_conv = nstates_conv + 1 812 END DO 813 814 WRITE (iter_unit, '(T7,I8,T24,F7.1,T40,ES11.4,T66,I8)') iter, t2 - t1, conv, nstates_conv 815 CALL m_flush(iter_unit) 816 END IF 817 818 t1 = t2 819 evals(1:nstates) = evals_last(1:nstates) 820 821 ! nvects_new == 0 if iter >= tddfpt_control%niters 822 IF (nvects_new == 0 .OR. is_nonortho) THEN 823 ! restart Davidson iterations 824 CALL tddfpt_orthogonalize_psi1_psi0(evects, work_matrices%S_C0_C0T) 825 CALL tddfpt_orthonormalize_psi1_psi1(evects, nstates, S_evects, matrix_s(1)%matrix) 826 827 EXIT 828 END IF 829 END DO 830 831 DEALLOCATE (Atilde) 832 833 DO istate = nvects_exist + nvects_new, 1, -1 834 DO ispin = nspins, 1, -1 835 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, Aop_krylov(ispin, istate)%matrix) 836 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, S_krylov(ispin, istate)%matrix) 837 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, krylov_vects(ispin, istate)%matrix) 838 END DO 839 END DO 840 DEALLOCATE (Aop_krylov, krylov_vects, S_krylov) 841 DEALLOCATE (evals_last) 842 843 DO istate = nstates, 1, -1 844 DO ispin = nspins, 1, -1 845 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, Aop_ritz(ispin, istate)%matrix) 846 END DO 847 END DO 848 DEALLOCATE (Aop_ritz) 849 850 CALL timestop(handle) 851 852 END FUNCTION tddfpt_davidson_solver 853 854END MODULE qs_tddfpt2_eigensolver 855