1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6MODULE qs_tddfpt2_methods 7 USE admm_types, ONLY: admm_type 8 USE bibliography, ONLY: Iannuzzi2005,& 9 cite_reference 10 USE cell_types, ONLY: cell_type 11 USE cp_array_utils, ONLY: cp_1d_r_p_type 12 USE cp_blacs_env, ONLY: cp_blacs_env_type,& 13 get_blacs_info 14 USE cp_cfm_basic_linalg, ONLY: cp_cfm_solve 15 USE cp_cfm_types, ONLY: cp_cfm_create,& 16 cp_cfm_p_type,& 17 cp_cfm_release,& 18 cp_cfm_set_all,& 19 cp_cfm_to_fm,& 20 cp_fm_to_cfm 21 USE cp_control_types, ONLY: dft_control_type,& 22 tddfpt2_control_type 23 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 24 copy_fm_to_dbcsr,& 25 cp_dbcsr_sm_fm_multiply,& 26 dbcsr_allocate_matrix_set,& 27 dbcsr_deallocate_matrix_set 28 USE cp_files, ONLY: close_file,& 29 open_file 30 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,& 31 cp_fm_contracted_trace,& 32 cp_fm_scale,& 33 cp_fm_scale_and_add,& 34 cp_fm_trace,& 35 cp_fm_triangular_invert 36 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose 37 USE cp_fm_diag, ONLY: choose_eigv_solver 38 USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,& 39 cp_fm_pool_type,& 40 fm_pool_create,& 41 fm_pool_create_fm,& 42 fm_pool_give_back_fm,& 43 fm_pool_release 44 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 45 cp_fm_struct_p_type,& 46 cp_fm_struct_release,& 47 cp_fm_struct_type 48 USE cp_fm_types, ONLY: & 49 cp_fm_copy_general, cp_fm_create, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_maxabsval, & 50 cp_fm_p_type, cp_fm_read_unformatted, cp_fm_release, cp_fm_set_all, cp_fm_set_submatrix, & 51 cp_fm_to_fm, cp_fm_type, cp_fm_write_unformatted 52 USE cp_gemm_interface, ONLY: cp_gemm 53 USE cp_log_handling, ONLY: cp_get_default_logger,& 54 cp_logger_get_default_io_unit,& 55 cp_logger_type 56 USE cp_output_handling, ONLY: cp_add_iter_level,& 57 cp_iterate,& 58 cp_p_file,& 59 cp_print_key_finished_output,& 60 cp_print_key_generate_filename,& 61 cp_print_key_should_output,& 62 cp_print_key_unit_nr,& 63 cp_rm_iter_level 64 USE cp_para_types, ONLY: cp_para_env_type 65 USE dbcsr_api, ONLY: & 66 dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_distribution_type, dbcsr_get_info, & 67 dbcsr_init_p, dbcsr_p_type, dbcsr_scale, dbcsr_set, dbcsr_type 68 USE hfx_admm_utils, ONLY: tddft_hfx_matrix 69 USE input_constants, ONLY: cholesky_dbcsr,& 70 cholesky_inverse,& 71 cholesky_off,& 72 cholesky_restore,& 73 tddfpt_dipole_berry,& 74 tddfpt_dipole_length,& 75 tddfpt_dipole_velocity 76 USE input_section_types, ONLY: section_get_ival,& 77 section_get_rval,& 78 section_vals_get,& 79 section_vals_get_subs_vals,& 80 section_vals_type,& 81 section_vals_val_get 82 USE kahan_sum, ONLY: accurate_dot_product,& 83 accurate_sum 84 USE kinds, ONLY: default_path_length,& 85 dp,& 86 int_8 87 USE machine, ONLY: m_flush,& 88 m_walltime 89 USE mathconstants, ONLY: twopi,& 90 z_one,& 91 z_zero 92 USE message_passing, ONLY: mp_bcast,& 93 mp_irecv,& 94 mp_isend,& 95 mp_min,& 96 mp_sum,& 97 mp_wait 98 USE moments_utils, ONLY: get_reference_point 99 USE physcon, ONLY: evolt 100 USE pw_env_types, ONLY: pw_env_get,& 101 pw_env_type 102 USE pw_methods, ONLY: pw_axpy,& 103 pw_scale,& 104 pw_transfer,& 105 pw_zero 106 USE pw_poisson_methods, ONLY: pw_poisson_solve 107 USE pw_poisson_types, ONLY: pw_poisson_type 108 USE pw_pool_types, ONLY: pw_pool_create_pw,& 109 pw_pool_give_back_pw,& 110 pw_pool_type 111 USE pw_types, ONLY: COMPLEXDATA1D,& 112 REALDATA3D,& 113 REALSPACE,& 114 RECIPROCALSPACE,& 115 pw_p_type,& 116 pw_type 117 USE qs_collocate_density, ONLY: calculate_rho_elec 118 USE qs_environment_types, ONLY: get_qs_env,& 119 qs_environment_type 120 USE qs_integrate_potential, ONLY: integrate_v_rspace 121 USE qs_ks_types, ONLY: qs_ks_env_type 122 USE qs_mo_types, ONLY: allocate_mo_set,& 123 deallocate_mo_set,& 124 get_mo_set,& 125 init_mo_set,& 126 mo_set_p_type,& 127 mo_set_type 128 USE qs_moments, ONLY: build_berry_moment_matrix 129 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 130 USE qs_operators_ao, ONLY: rRc_xyz_ao 131 USE qs_rho_methods, ONLY: qs_rho_rebuild,& 132 qs_rho_update_rho 133 USE qs_rho_types, ONLY: qs_rho_create,& 134 qs_rho_get,& 135 qs_rho_release,& 136 qs_rho_set,& 137 qs_rho_type 138 USE qs_scf_methods, ONLY: eigensolver 139 USE qs_scf_post_gpw, ONLY: make_lumo 140 USE qs_scf_types, ONLY: ot_method_nr,& 141 qs_scf_env_type 142 USE qs_tddfpt2_subgroups, ONLY: tddfpt_dbcsr_create_by_dist,& 143 tddfpt_sub_env_init,& 144 tddfpt_sub_env_release,& 145 tddfpt_subgroup_env_type 146 USE string_utilities, ONLY: integer_to_string 147 USE util, ONLY: sort 148 USE xc, ONLY: xc_calc_2nd_deriv,& 149 xc_prep_2nd_deriv 150 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 151 xc_dset_release 152 USE xc_derivatives, ONLY: xc_functionals_get_needs 153 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 154 USE xc_rho_set_types, ONLY: xc_rho_set_create,& 155 xc_rho_set_release,& 156 xc_rho_set_type,& 157 xc_rho_set_update 158#include "./base/base_uses.f90" 159 160 IMPLICIT NONE 161 162 PRIVATE 163 164 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_methods' 165 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. 166 ! number of first derivative components (3: d/dx, d/dy, d/dz) 167 INTEGER, PARAMETER, PRIVATE :: nderivs = 3 168 INTEGER, PARAMETER, PRIVATE :: maxspins = 2 169 170 PUBLIC :: tddfpt 171 172! ************************************************************************************************** 173!> \brief Ground state molecular orbitals. 174!> \par History 175!> * 06.2016 created [Sergey Chulkov] 176! ************************************************************************************************** 177 TYPE tddfpt_ground_state_mos 178 !> occupied MOs stored in a matrix form [nao x nmo_occ] 179 TYPE(cp_fm_type), POINTER :: mos_occ 180 !> virtual MOs stored in a matrix form [nao x nmo_virt] 181 TYPE(cp_fm_type), POINTER :: mos_virt 182 !> occupied orbital energies 183 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals_occ 184 !> virtual orbital energies 185 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals_virt 186 !> phase of occupied MOs; +1.0 -- positive, -1.0 -- negative; 187 !> it is mainly needed to make the restart file transferable 188 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: phases_occ 189 END TYPE tddfpt_ground_state_mos 190 191! ************************************************************************************************** 192!> \brief Set of temporary ("work") matrices. 193!> \par History 194!> * 01.2017 created [Sergey Chulkov] 195! ************************************************************************************************** 196 TYPE tddfpt_work_matrices 197 ! 198 ! *** globally distributed dense matrices *** 199 ! 200 !> pool of dense [nao x nmo_occ(spin)] matrices; 201 !> used mainly to dynamically expand the list of trial vectors 202 TYPE(cp_fm_pool_p_type), ALLOCATABLE, DIMENSION(:) :: fm_pool_ao_mo_occ 203 !> S * mos_occ(spin) 204 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: S_C0 205 !> S * \rho_0(spin) 206 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: S_C0_C0T 207 !> globally distributed dense matrices with shape [nao x nmo_occ(spin)] 208 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: wfm_ao_mo_occ 209 !> globally distributed dense matrices with shape [nmo_occ(spin) x nmo_occ(spin)] 210 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: wfm_mo_occ_mo_occ 211 !> globally distributed dense matrices with shape [nmo_virt(spin) x nmo_occ(spin)] 212 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: wfm_mo_virt_mo_occ 213 ! 214 ! *** dense matrices distributed across parallel (sub)groups *** 215 ! 216 !> evects_sub(1:nspins, 1:nstates): a copy of the last 'nstates' trial vectors distributed 217 !> across parallel (sub)groups. Here 'nstates' is the number of requested excited states which 218 !> is typically much smaller than the total number of Krylov's vectors. Allocated only if 219 !> the number of parallel groups > 1, otherwise we use the original globally distributed vectors. 220 !> evects_sub(spin, state) == null() means that the trial vector is assigned to a different (sub)group 221 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:, :) :: evects_sub 222 !> action of TDDFPT operator on trial vectors distributed across parallel (sub)groups 223 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:, :) :: Aop_evects_sub 224 !> electron density expressed in terms of atomic orbitals using primary basis set 225 TYPE(cp_fm_type), POINTER :: rho_ao_orb_fm_sub 226 ! 227 ! NOTE: we do not need the next 2 matrices in case of a sparse matrix 'tddfpt_subgroup_env_type%admm_A' 228 ! 229 !> electron density expressed in terms of atomic orbitals using auxiliary basis set; 230 !> can be seen as a group-specific version of the matrix 'admm_type%work_aux_aux' 231 TYPE(cp_fm_type), POINTER :: rho_ao_aux_fit_fm_sub 232 !> group-specific version of the matrix 'admm_type%work_aux_orb' with shape [nao_aux x nao] 233 TYPE(cp_fm_type), POINTER :: wfm_aux_orb_sub 234 ! 235 ! *** sparse matrices distributed across parallel (sub)groups *** 236 ! 237 !> sparse matrix with shape [nao x nao] distributed across subgroups; 238 !> Aop_evects_sub(spin,:) = A_ia_munu_sub(spin) * mos_occ(spin) 239 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: A_ia_munu_sub 240 ! 241 ! *** structures to store electron densities distributed across parallel (sub)groups *** 242 ! 243 !> electron density in terms of primary basis set 244 TYPE(qs_rho_type), POINTER :: rho_orb_struct_sub 245 !> electron density in terms of auxiliary basis set 246 TYPE(qs_rho_type), POINTER :: rho_aux_fit_struct_sub 247 !> group-specific copy of a Coulomb/xc-potential on a real-space grid 248 TYPE(pw_p_type), DIMENSION(:), POINTER :: A_ia_rspace_sub 249 !> group-specific copy of a reciprocal-space grid 250 TYPE(pw_p_type), DIMENSION(:), POINTER :: wpw_gspace_sub 251 !> group-specific copy of a real-space grid 252 TYPE(pw_p_type), DIMENSION(:), POINTER :: wpw_rspace_sub 253 ! 254 ! *** globally distributed matrices required to compute exact exchange terms *** 255 ! 256 !> globally distributed version of the matrix 'rho_ao_orb_fm_sub' to store the electron density 257 TYPE(cp_fm_type), POINTER :: hfx_fm_ao_ao 258 !> sparse matrix to store the electron density in terms of auxiliary (ADMM calculation) 259 !> or primary (regular calculation) basis set 260 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: hfx_rho_ao 261 !> exact exchange expressed in terms of auxiliary or primary basis set 262 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: hfx_hmat 263 END TYPE tddfpt_work_matrices 264 265! ************************************************************************************************** 266!> \brief Collection of variables required to evaluate adiabatic TDDFPT kernel. 267!> \par History 268!> * 12.2016 created [Sergey Chulkov] 269! ************************************************************************************************** 270 TYPE tddfpt_kernel_env_type 271 ! ground state electron density 272 TYPE(xc_rho_set_type), POINTER :: xc_rho_set 273 ! response density 274 TYPE(xc_rho_set_type), POINTER :: xc_rho1_set 275 !> first and second derivatives of exchange-correlation functional 276 TYPE(xc_derivative_set_type), POINTER :: xc_deriv_set 277 !> XC input section 278 TYPE(section_vals_type), POINTER :: xc_section 279 !> flags which indicate required components of the exchange-correlation functional 280 !> (density, gradient, etc) 281 TYPE(xc_rho_cflags_type) :: xc_rho1_cflags 282 !> the method used to compute derivatives 283 INTEGER :: deriv_method_id 284 !> the density smoothing method 285 INTEGER :: rho_smooth_id 286 !> scaling coefficients in the linear combination: 287 !> K = alpha * K_{\alpha,\alpha} + beta * K_{\alpha,\beta} 288 REAL(kind=dp) :: alpha, beta 289 END TYPE tddfpt_kernel_env_type 290 291CONTAINS 292 293! ************************************************************************************************** 294!> \brief Perform TDDFPT calculation. 295!> \param qs_env Quickstep environment 296!> \par History 297!> * 05.2016 created [Sergey Chulkov] 298!> * 06.2016 refactored to be used with Davidson eigensolver [Sergey Chulkov] 299!> * 03.2017 cleaned and refactored [Sergey Chulkov] 300!> \note Based on the subroutines tddfpt_env_init(), and tddfpt_env_deallocate(). 301! ************************************************************************************************** 302 SUBROUTINE tddfpt(qs_env) 303 TYPE(qs_environment_type), POINTER :: qs_env 304 305 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt', routineP = moduleN//':'//routineN 306 307 CHARACTER(len=20) :: nstates_str 308 INTEGER :: energy_unit, handle, ispin, istate, iter, log_unit, mult, niters, nmo_avail, & 309 nmo_occ, nmo_virt, nspins, nstates, nstates_read 310 LOGICAL :: do_admm, do_hfx, explicit_xc, & 311 is_restarted 312 REAL(kind=dp) :: C_hf, conv 313 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals 314 REAL(kind=dp), DIMENSION(:), POINTER :: evals_virt_spin 315 TYPE(admm_type), POINTER :: admm_env 316 TYPE(cell_type), POINTER :: cell 317 TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: evals_virt 318 TYPE(cp_blacs_env_type), POINTER :: blacs_env 319 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:, :) :: dipole_op_mos_occ, evects, S_evects 320 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mos_virt 321 TYPE(cp_fm_type), POINTER :: mos_virt_spin 322 TYPE(cp_logger_type), POINTER :: logger 323 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s 324 TYPE(dft_control_type), POINTER :: dft_control 325 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 326 TYPE(qs_scf_env_type), POINTER :: scf_env 327 TYPE(section_vals_type), POINTER :: hfx_section, input, & 328 tddfpt_print_section, tddfpt_section, & 329 xc_section 330 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control 331 TYPE(tddfpt_ground_state_mos), ALLOCATABLE, & 332 DIMENSION(:) :: gs_mos 333 TYPE(tddfpt_kernel_env_type) :: kernel_env, kernel_env_admm_aux 334 TYPE(tddfpt_subgroup_env_type) :: sub_env 335 TYPE(tddfpt_work_matrices) :: work_matrices 336 337 CALL timeset(routineN, handle) 338 logger => cp_get_default_logger() 339 340 CALL cite_reference(Iannuzzi2005) 341 342 NULLIFY (blacs_env, cell, dft_control, input, matrix_ks, matrix_s, mos) 343 CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, dft_control=dft_control, & 344 input=input, matrix_ks=matrix_ks, matrix_s=matrix_s, mos=mos, scf_env=scf_env) 345 tddfpt_control => dft_control%tddfpt2_control 346 347 nspins = dft_control%nspins 348 CPASSERT(dft_control%nimages <= 1) 349 350 IF (tddfpt_control%nstates <= 0) THEN 351 CALL integer_to_string(tddfpt_control%nstates, nstates_str) 352 CALL cp_warn(__LOCATION__, "TDDFPT calculation was requested for "// & 353 TRIM(nstates_str)//" excited states: nothing to do.") 354 CALL timestop(handle) 355 RETURN 356 END IF 357 358 tddfpt_section => section_vals_get_subs_vals(input, "PROPERTIES%TDDFPT") 359 tddfpt_print_section => section_vals_get_subs_vals(tddfpt_section, "PRINT") 360 361 xc_section => section_vals_get_subs_vals(tddfpt_section, "XC%XC_FUNCTIONAL") 362 CALL section_vals_get(xc_section, explicit=explicit_xc) 363 IF (explicit_xc) THEN 364 xc_section => section_vals_get_subs_vals(tddfpt_section, "XC") 365 ELSE 366 xc_section => section_vals_get_subs_vals(input, "DFT%XC") 367 END IF 368 hfx_section => section_vals_get_subs_vals(xc_section, "HF") 369 370 CALL section_vals_get(hfx_section, explicit=do_hfx) 371 IF (do_hfx) THEN 372 CALL section_vals_val_get(hfx_section, "FRACTION", r_val=C_hf) 373 do_hfx = (C_hf /= 0.0_dp) 374 END IF 375 376 do_admm = do_hfx .AND. dft_control%do_admm 377 IF (do_admm) THEN 378 IF (explicit_xc) THEN 379 ! 'admm_env%xc_section_primary' and 'admm_env%xc_section_aux' need to be redefined 380 CALL cp_abort(__LOCATION__, & 381 "ADMM is not implemented for a TDDFT kernel XC-functional which is different from "// & 382 "the one used for the ground-state calculation. A ground-state 'admm_env' cannot be reused.") 383 END IF 384 385 CALL get_qs_env(qs_env, admm_env=admm_env) 386 END IF 387 388 ! reset rks_triplets if UKS is in use 389 IF (tddfpt_control%rks_triplets .AND. nspins > 1) THEN 390 tddfpt_control%rks_triplets = .FALSE. 391 CALL cp_warn(__LOCATION__, "Keyword RKS_TRIPLETS has been ignored for spin-polarised calculations") 392 END IF 393 394 ! obtain occupied and virtual (unoccupied) ground-state Kohn-Sham orbitals 395 ALLOCATE (gs_mos(nspins)) 396 397 ! when the number of unoccupied orbitals is limited and OT has been used for the ground-state DFT calculation, 398 ! compute the missing unoccupied orbitals using OT as well. 399 NULLIFY (evals_virt, evals_virt_spin, mos_virt, mos_virt_spin) 400 IF (ASSOCIATED(scf_env)) THEN 401 IF (scf_env%method == ot_method_nr .AND. tddfpt_control%nlumo > 0) THEN 402 ! As OT with ADDED_MOS/=0 is currently not implemented, the following block is equivalent to: 403 ! nmo_virt = tddfpt_control%nlumo 404 405 ! number of already computed unoccupied orbitals (added_mos) . 406 nmo_virt = HUGE(nmo_virt) 407 DO ispin = 1, nspins 408 CALL get_mo_set(mos(ispin)%mo_set, nmo=nmo_avail, homo=nmo_occ) 409 nmo_virt = MIN(nmo_virt, nmo_avail - nmo_occ) 410 END DO 411 ! number of unoccupied orbitals to compute 412 nmo_virt = tddfpt_control%nlumo - nmo_virt 413 414 IF (nmo_virt > 0) THEN 415 ALLOCATE (evals_virt(nspins), mos_virt(nspins)) 416 ! the number of actually computed unoccupied orbitals will be stored as 'nmo_avail' 417 CALL make_lumo(qs_env, scf_env, mos_virt, evals_virt, nmo_virt, nmo_avail) 418 END IF 419 END IF 420 END IF 421 422 DO ispin = 1, nspins 423 IF (ASSOCIATED(evals_virt)) & 424 evals_virt_spin => evals_virt(ispin)%array 425 IF (ASSOCIATED(mos_virt)) & 426 mos_virt_spin => mos_virt(ispin)%matrix 427 428 CALL tddfpt_init_ground_state_mos(gs_mos=gs_mos(ispin), mo_set=mos(ispin)%mo_set, nlumo=tddfpt_control%nlumo, & 429 blacs_env=blacs_env, cholesky_method=cholesky_restore, & 430 matrix_ks=matrix_ks(ispin)%matrix, matrix_s=matrix_s(1)%matrix, & 431 mos_virt=mos_virt_spin, evals_virt=evals_virt_spin) 432 END DO 433 434 IF (ASSOCIATED(evals_virt)) THEN 435 DO ispin = 1, SIZE(evals_virt) 436 IF (ASSOCIATED(evals_virt(ispin)%array)) DEALLOCATE (evals_virt(ispin)%array) 437 END DO 438 DEALLOCATE (evals_virt) 439 END IF 440 441 IF (ASSOCIATED(mos_virt)) THEN 442 DO ispin = 1, SIZE(mos_virt) 443 IF (ASSOCIATED(mos_virt(ispin)%matrix)) CALL cp_fm_release(mos_virt(ispin)%matrix) 444 END DO 445 DEALLOCATE (mos_virt) 446 END IF 447 448 ! components of the dipole operator 449 CALL tddfpt_dipole_operator(dipole_op_mos_occ, tddfpt_control, gs_mos, qs_env) 450 451 ! multiplicity of molecular system 452 IF (nspins > 1) THEN 453 mult = ABS(SIZE(gs_mos(1)%evals_occ) - SIZE(gs_mos(2)%evals_occ)) + 1 454 IF (mult > 2) & 455 CALL cp_warn(__LOCATION__, "There is a convergence issue for multiplicity >= 3") 456 ELSE 457 IF (tddfpt_control%rks_triplets) THEN 458 mult = 3 459 ELSE 460 mult = 1 461 END IF 462 END IF 463 464 ! split mpi communicator 465 ALLOCATE (evects(nspins, 1)) 466 DO ispin = 1, nspins 467 evects(ispin, 1)%matrix => gs_mos(ispin)%mos_occ 468 END DO 469 CALL tddfpt_sub_env_init(sub_env, qs_env, mos_occ=evects(:, 1)) 470 DEALLOCATE (evects) 471 472 ! allocate pools and work matrices 473 nstates = tddfpt_control%nstates 474 CALL tddfpt_create_work_matrices(work_matrices, gs_mos, nstates, do_hfx, qs_env, sub_env) 475 476 ! create kernel environment 477 CALL tddfpt_construct_ground_state_orb_density(rho_orb_struct=work_matrices%rho_orb_struct_sub, & 478 is_rks_triplets=tddfpt_control%rks_triplets, & 479 qs_env=qs_env, sub_env=sub_env, & 480 wfm_rho_orb=work_matrices%rho_ao_orb_fm_sub) 481 482 IF (do_admm) THEN 483 CALL tddfpt_create_kernel_env(kernel_env=kernel_env, & 484 rho_struct_sub=work_matrices%rho_orb_struct_sub, & 485 xc_section=admm_env%xc_section_primary, & 486 is_rks_triplets=tddfpt_control%rks_triplets, sub_env=sub_env) 487 488 CALL tddfpt_construct_aux_fit_density(rho_orb_struct=work_matrices%rho_orb_struct_sub, & 489 rho_aux_fit_struct=work_matrices%rho_aux_fit_struct_sub, & 490 qs_env=qs_env, sub_env=sub_env, & 491 wfm_rho_orb=work_matrices%rho_ao_orb_fm_sub, & 492 wfm_rho_aux_fit=work_matrices%rho_ao_aux_fit_fm_sub, & 493 wfm_aux_orb=work_matrices%wfm_aux_orb_sub) 494 495 CALL tddfpt_create_kernel_env(kernel_env=kernel_env_admm_aux, & 496 rho_struct_sub=work_matrices%rho_aux_fit_struct_sub, & 497 xc_section=admm_env%xc_section_aux, & 498 is_rks_triplets=tddfpt_control%rks_triplets, sub_env=sub_env) 499 ELSE 500 CALL tddfpt_create_kernel_env(kernel_env=kernel_env, & 501 rho_struct_sub=work_matrices%rho_orb_struct_sub, & 502 xc_section=xc_section, & 503 is_rks_triplets=tddfpt_control%rks_triplets, sub_env=sub_env) 504 END IF 505 506 ALLOCATE (evals(nstates)) 507 ALLOCATE (evects(nspins, nstates), S_evects(nspins, nstates)) 508 DO istate = 1, nstates 509 DO ispin = 1, nspins 510 NULLIFY (evects(ispin, istate)%matrix, S_evects(ispin, istate)%matrix) 511 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, S_evects(ispin, istate)%matrix) 512 END DO 513 END DO 514 515 ! reuse Ritz vectors from the previous calculation if available 516 IF (tddfpt_control%is_restart) THEN 517 nstates_read = tddfpt_read_restart(evects=evects, evals=evals, gs_mos=gs_mos, & 518 logger=logger, tddfpt_section=tddfpt_section, & 519 tddfpt_print_section=tddfpt_print_section, & 520 fm_pool_ao_mo_occ=work_matrices%fm_pool_ao_mo_occ, & 521 blacs_env_global=blacs_env) 522 ELSE 523 nstates_read = 0 524 END IF 525 526 is_restarted = nstates_read >= nstates 527 528 ! build the list of missed singly excited states and sort them in ascending order according to their excitation energies 529 log_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, "GUESS_VECTORS", extension=".tddfptLog") 530 CALL tddfpt_guess_vectors(evects=evects, evals=evals, gs_mos=gs_mos, log_unit=log_unit) 531 CALL cp_print_key_finished_output(log_unit, logger, tddfpt_print_section, "GUESS_VECTORS") 532 533 CALL tddfpt_orthogonalize_psi1_psi0(evects, work_matrices%S_C0_C0T, work_matrices%wfm_ao_mo_occ) 534 CALL tddfpt_orthonormalize_psi1_psi1(evects, nstates, S_evects, matrix_s(1)%matrix) 535 536 niters = tddfpt_control%niters 537 IF (niters > 0) THEN 538 log_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, "ITERATION_INFO", extension=".tddfptLog") 539 energy_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, "DETAILED_ENERGY", extension=".tddfptLog") 540 541 IF (log_unit > 0) THEN 542 WRITE (log_unit, '(/,1X,A,/)') 'TDDFPT WAVEFUNCTION OPTIMIZATION' 543 WRITE (log_unit, '(5X,A,T15,A,T24,A,T40,A)') "Step", "Time", "Convergence", "Conv. states" 544 WRITE (log_unit, '(1X,50("-"))') 545 END IF 546 547 CALL cp_add_iter_level(logger%iter_info, "TDDFT_SCF") 548 549 DO 550 ! *** perform Davidson iterations *** 551 conv = tddfpt_davidson_solver(evects=evects, evals=evals, S_evects=S_evects, gs_mos=gs_mos, & 552 do_hfx=do_hfx, tddfpt_control=tddfpt_control, qs_env=qs_env, & 553 kernel_env=kernel_env, kernel_env_admm_aux=kernel_env_admm_aux, & 554 sub_env=sub_env, logger=logger, & 555 iter_unit=log_unit, energy_unit=energy_unit, & 556 tddfpt_print_section=tddfpt_print_section, work_matrices=work_matrices) 557 558 ! at this point at least one of the following conditions are met: 559 ! a) convergence criteria has been achieved; 560 ! b) maximum number of iterations has been reached; 561 ! c) Davidson iterations must be restarted due to lack of Krylov vectors or numerical instability 562 563 CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter) 564 ! terminate the loop if either (a) or (b) is true ... 565 IF ((conv <= tddfpt_control%conv .AND. is_restarted) .OR. iter >= niters) EXIT 566 567 ! ... otherwise restart Davidson iterations 568 is_restarted = .TRUE. 569 IF (log_unit > 0) THEN 570 WRITE (log_unit, '(1X,10("-"),1X,A,1X,11("-"))') "Restart Davidson iterations" 571 CALL m_flush(log_unit) 572 END IF 573 END DO 574 575 ! write TDDFPT restart file at the last iteration if requested to do so 576 CALL cp_iterate(logger%iter_info, increment=0, last=.TRUE.) 577 CALL tddfpt_write_restart(evects=evects, evals=evals, gs_mos=gs_mos, & 578 logger=logger, tddfpt_print_section=tddfpt_print_section) 579 580 CALL cp_rm_iter_level(logger%iter_info, "TDDFT_SCF") 581 582 ! print convergence summary 583 IF (log_unit > 0) THEN 584 CALL integer_to_string(iter, nstates_str) 585 IF (conv <= tddfpt_control%conv) THEN 586 WRITE (log_unit, '(/,1X,A)') "*** TDDFPT run converged in "//TRIM(nstates_str)//" iteration(s) ***" 587 ELSE 588 WRITE (log_unit, '(/,1X,A)') "*** TDDFPT run did NOT converge after "//TRIM(nstates_str)//" iteration(s) ***" 589 END IF 590 END IF 591 592 CALL cp_print_key_finished_output(energy_unit, logger, tddfpt_print_section, "DETAILED_ENERGY") 593 CALL cp_print_key_finished_output(log_unit, logger, tddfpt_print_section, "ITERATION_INFO") 594 ELSE 595 CALL cp_warn(__LOCATION__, "Skipping TDDFPT wavefunction optimization") 596 END IF 597 598 ! *** print summary information *** 599 log_unit = cp_logger_get_default_io_unit() 600 601 CALL tddfpt_print_summary(log_unit, evects, evals, mult, dipole_op_mos_occ) 602 CALL tddfpt_print_excitation_analysis(log_unit, evects, gs_mos, matrix_s(1)%matrix, & 603 min_amplitude=tddfpt_control%min_excitation_amplitude) 604 605 ! -- clean up all useless stuff 606 DO istate = SIZE(evects, 2), 1, -1 607 DO ispin = nspins, 1, -1 608 CALL cp_fm_release(evects(ispin, istate)%matrix) 609 CALL cp_fm_release(S_evects(ispin, istate)%matrix) 610 END DO 611 END DO 612 DEALLOCATE (evects, S_evects, evals) 613 614 IF (do_admm) THEN 615 CALL tddfpt_release_kernel_env(kernel_env_admm_aux) 616 END IF 617 CALL tddfpt_release_kernel_env(kernel_env) 618 CALL tddfpt_release_work_matrices(work_matrices, sub_env) 619 CALL tddfpt_sub_env_release(sub_env) 620 621 IF (ALLOCATED(dipole_op_mos_occ)) THEN 622 DO ispin = nspins, 1, -1 623 DO istate = SIZE(dipole_op_mos_occ, 1), 1, -1 624 CALL cp_fm_release(dipole_op_mos_occ(istate, ispin)%matrix) 625 END DO 626 END DO 627 DEALLOCATE (dipole_op_mos_occ) 628 END IF 629 630 DO ispin = nspins, 1, -1 631 CALL tddfpt_release_ground_state_mos(gs_mos(ispin)) 632 END DO 633 DEALLOCATE (gs_mos) 634 635 CALL timestop(handle) 636 END SUBROUTINE tddfpt 637 638! ************************************************************************************************** 639!> \brief Allocate work matrices. 640!> \param work_matrices work matrices (allocated on exit) 641!> \param gs_mos occupied and virtual molecular orbitals optimised for the ground state 642!> \param nstates number of excited states to converge 643!> \param do_hfx flag that requested to allocate work matrices required for computation 644!> of exact-exchange terms 645!> \param qs_env Quickstep environment 646!> \param sub_env parallel group environment 647!> \par History 648!> * 02.2017 created [Sergey Chulkov] 649! ************************************************************************************************** 650 SUBROUTINE tddfpt_create_work_matrices(work_matrices, gs_mos, nstates, do_hfx, qs_env, sub_env) 651 TYPE(tddfpt_work_matrices), INTENT(out) :: work_matrices 652 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 653 INTENT(in) :: gs_mos 654 INTEGER, INTENT(in) :: nstates 655 LOGICAL, INTENT(in) :: do_hfx 656 TYPE(qs_environment_type), POINTER :: qs_env 657 TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env 658 659 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_create_work_matrices', & 660 routineP = moduleN//':'//routineN 661 662 INTEGER :: handle, igroup, ispin, istate, nao, & 663 nao_aux, ngroups, nspins 664 INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_virt 665 LOGICAL :: do_admm 666 TYPE(cp_blacs_env_type), POINTER :: blacs_env 667 TYPE(cp_fm_struct_p_type), DIMENSION(maxspins) :: fm_struct_evects 668 TYPE(cp_fm_struct_type), POINTER :: fm_struct 669 TYPE(cp_para_env_type), POINTER :: para_env 670 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist 671 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_s_aux_fit, rho_ia_ao 672 TYPE(dbcsr_type), POINTER :: dbcsr_template_hfx 673 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 674 POINTER :: sab_hfx 675 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 676 677 CALL timeset(routineN, handle) 678 679 nspins = SIZE(gs_mos) 680 CALL get_qs_env(qs_env, blacs_env=blacs_env, matrix_s=matrix_s) 681 CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao) 682 683 DO ispin = 1, nspins 684 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ) 685 nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt) 686 END DO 687 688 do_admm = do_hfx .AND. ASSOCIATED(sub_env%admm_A) 689 IF (do_admm) THEN 690 CALL get_qs_env(qs_env, matrix_s_aux_fit=matrix_s_aux_fit) 691 CALL dbcsr_get_info(matrix_s_aux_fit(1)%matrix, nfullrows_total=nao_aux) 692 END IF 693 694 NULLIFY (fm_struct) 695 ALLOCATE (work_matrices%fm_pool_ao_mo_occ(nspins)) 696 DO ispin = 1, nspins 697 NULLIFY (work_matrices%fm_pool_ao_mo_occ(ispin)%pool) 698 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nmo_occ(ispin), context=blacs_env) 699 CALL fm_pool_create(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, fm_struct) 700 CALL cp_fm_struct_release(fm_struct) 701 END DO 702 703 ALLOCATE (work_matrices%S_C0_C0T(nspins)) 704 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env) 705 DO ispin = 1, nspins 706 NULLIFY (work_matrices%S_C0_C0T(ispin)%matrix) 707 CALL cp_fm_create(work_matrices%S_C0_C0T(ispin)%matrix, fm_struct) 708 END DO 709 CALL cp_fm_struct_release(fm_struct) 710 711 ALLOCATE (work_matrices%S_C0(nspins), work_matrices%wfm_ao_mo_occ(nspins)) 712 ALLOCATE (work_matrices%wfm_mo_occ_mo_occ(nspins), work_matrices%wfm_mo_virt_mo_occ(nspins)) 713 DO ispin = 1, nspins 714 NULLIFY (work_matrices%S_C0(ispin)%matrix, work_matrices%wfm_ao_mo_occ(ispin)%matrix) 715 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, work_matrices%S_C0(ispin)%matrix) 716 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, work_matrices%wfm_ao_mo_occ(ispin)%matrix) 717 718 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, gs_mos(ispin)%mos_occ, work_matrices%S_C0(ispin)%matrix, & 719 ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp) 720 CALL cp_gemm('N', 'T', nao, nao, nmo_occ(ispin), 1.0_dp, work_matrices%S_C0(ispin)%matrix, & 721 gs_mos(ispin)%mos_occ, 0.0_dp, work_matrices%S_C0_C0T(ispin)%matrix) 722 723 NULLIFY (work_matrices%wfm_mo_occ_mo_occ(ispin)%matrix) 724 CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_occ(ispin), ncol_global=nmo_occ(ispin), context=blacs_env) 725 CALL cp_fm_create(work_matrices%wfm_mo_occ_mo_occ(ispin)%matrix, fm_struct) 726 CALL cp_fm_struct_release(fm_struct) 727 728 NULLIFY (work_matrices%wfm_mo_virt_mo_occ(ispin)%matrix) 729 CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(ispin), ncol_global=nmo_occ(ispin), context=blacs_env) 730 CALL cp_fm_create(work_matrices%wfm_mo_virt_mo_occ(ispin)%matrix, fm_struct) 731 CALL cp_fm_struct_release(fm_struct) 732 END DO 733 734 IF (sub_env%is_split) THEN 735 DO ispin = 1, nspins 736 NULLIFY (fm_struct_evects(ispin)%struct) 737 CALL cp_fm_struct_create(fm_struct_evects(ispin)%struct, nrow_global=nao, & 738 ncol_global=nmo_occ(ispin), context=sub_env%blacs_env) 739 END DO 740 741 ALLOCATE (work_matrices%evects_sub(nspins, nstates), work_matrices%Aop_evects_sub(nspins, nstates)) 742 DO istate = 1, nstates 743 DO ispin = 1, nspins 744 NULLIFY (work_matrices%evects_sub(ispin, istate)%matrix) 745 NULLIFY (work_matrices%Aop_evects_sub(ispin, istate)%matrix) 746 END DO 747 END DO 748 749 CALL get_blacs_info(blacs_env, para_env=para_env) 750 igroup = sub_env%group_distribution(para_env%mepos) 751 ngroups = sub_env%ngroups 752 753 DO istate = ngroups - igroup, nstates, ngroups 754 DO ispin = 1, nspins 755 CALL cp_fm_create(work_matrices%evects_sub(ispin, istate)%matrix, fm_struct_evects(ispin)%struct) 756 CALL cp_fm_create(work_matrices%Aop_evects_sub(ispin, istate)%matrix, fm_struct_evects(ispin)%struct) 757 END DO 758 END DO 759 760 DO ispin = nspins, 1, -1 761 CALL cp_fm_struct_release(fm_struct_evects(ispin)%struct) 762 END DO 763 END IF 764 765 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=sub_env%blacs_env) 766 NULLIFY (work_matrices%rho_ao_orb_fm_sub) 767 CALL cp_fm_create(work_matrices%rho_ao_orb_fm_sub, fm_struct) 768 CALL cp_fm_struct_release(fm_struct) 769 770 NULLIFY (work_matrices%rho_ao_aux_fit_fm_sub, work_matrices%wfm_aux_orb_sub) 771 IF (do_admm) THEN 772 CALL cp_fm_struct_create(fm_struct, nrow_global=nao_aux, ncol_global=nao_aux, context=sub_env%blacs_env) 773 CALL cp_fm_create(work_matrices%rho_ao_aux_fit_fm_sub, fm_struct) 774 CALL cp_fm_struct_release(fm_struct) 775 776 CALL cp_fm_struct_create(fm_struct, nrow_global=nao_aux, ncol_global=nao, context=sub_env%blacs_env) 777 CALL cp_fm_create(work_matrices%wfm_aux_orb_sub, fm_struct) 778 CALL cp_fm_struct_release(fm_struct) 779 END IF 780 781 ! group-specific dbcsr matrices 782 NULLIFY (work_matrices%A_ia_munu_sub) 783 CALL dbcsr_allocate_matrix_set(work_matrices%A_ia_munu_sub, nspins) 784 DO ispin = 1, nspins 785 CALL dbcsr_init_p(work_matrices%A_ia_munu_sub(ispin)%matrix) 786 CALL tddfpt_dbcsr_create_by_dist(work_matrices%A_ia_munu_sub(ispin)%matrix, template=matrix_s(1)%matrix, & 787 dbcsr_dist=sub_env%dbcsr_dist, sab=sub_env%sab_orb) 788 END DO 789 790 ! group-specific response density 791 NULLIFY (rho_ia_ao) 792 CALL dbcsr_allocate_matrix_set(rho_ia_ao, nspins) 793 DO ispin = 1, nspins 794 CALL dbcsr_init_p(rho_ia_ao(ispin)%matrix) 795 CALL tddfpt_dbcsr_create_by_dist(rho_ia_ao(ispin)%matrix, template=matrix_s(1)%matrix, & 796 dbcsr_dist=sub_env%dbcsr_dist, sab=sub_env%sab_orb) 797 END DO 798 799 NULLIFY (work_matrices%rho_orb_struct_sub) 800 CALL qs_rho_create(work_matrices%rho_orb_struct_sub) 801 CALL qs_rho_set(work_matrices%rho_orb_struct_sub, rho_ao=rho_ia_ao) 802 CALL qs_rho_rebuild(work_matrices%rho_orb_struct_sub, qs_env, rebuild_ao=.FALSE., & 803 rebuild_grids=.TRUE., pw_env_external=sub_env%pw_env) 804 805 NULLIFY (work_matrices%rho_aux_fit_struct_sub) 806 IF (do_admm) THEN 807 NULLIFY (rho_ia_ao) 808 CALL dbcsr_allocate_matrix_set(rho_ia_ao, nspins) 809 DO ispin = 1, nspins 810 CALL dbcsr_init_p(rho_ia_ao(ispin)%matrix) 811 CALL tddfpt_dbcsr_create_by_dist(rho_ia_ao(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix, & 812 dbcsr_dist=sub_env%dbcsr_dist, sab=sub_env%sab_aux_fit) 813 END DO 814 815 CALL qs_rho_create(work_matrices%rho_aux_fit_struct_sub) 816 CALL qs_rho_set(work_matrices%rho_aux_fit_struct_sub, rho_ao=rho_ia_ao) 817 CALL qs_rho_rebuild(work_matrices%rho_aux_fit_struct_sub, qs_env, rebuild_ao=.FALSE., & 818 rebuild_grids=.TRUE., pw_env_external=sub_env%pw_env) 819 END IF 820 821 ! work plain-wave grids 822 CALL pw_env_get(sub_env%pw_env, auxbas_pw_pool=auxbas_pw_pool) 823 ALLOCATE (work_matrices%A_ia_rspace_sub(nspins)) 824 ALLOCATE (work_matrices%wpw_gspace_sub(nspins), work_matrices%wpw_rspace_sub(nspins)) 825 DO ispin = 1, nspins 826 NULLIFY (work_matrices%A_ia_rspace_sub(ispin)%pw) 827 CALL pw_pool_create_pw(auxbas_pw_pool, work_matrices%A_ia_rspace_sub(ispin)%pw, & 828 use_data=REALDATA3D, in_space=REALSPACE) 829 830 NULLIFY (work_matrices%wpw_gspace_sub(ispin)%pw, work_matrices%wpw_rspace_sub(ispin)%pw) 831 CALL pw_pool_create_pw(auxbas_pw_pool, work_matrices%wpw_gspace_sub(ispin)%pw, & 832 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 833 CALL pw_pool_create_pw(auxbas_pw_pool, work_matrices%wpw_rspace_sub(ispin)%pw, & 834 use_data=REALDATA3D, in_space=REALSPACE) 835 END DO 836 837 ! HFX-related globally distributed matrices 838 NULLIFY (work_matrices%hfx_fm_ao_ao, work_matrices%hfx_rho_ao, work_matrices%hfx_hmat) 839 IF (do_hfx) THEN 840 IF (do_admm) THEN 841 CALL get_qs_env(qs_env, dbcsr_dist=dbcsr_dist, sab_aux_fit=sab_hfx) 842 dbcsr_template_hfx => matrix_s_aux_fit(1)%matrix 843 ELSE 844 CALL get_qs_env(qs_env, dbcsr_dist=dbcsr_dist, sab_orb=sab_hfx) 845 dbcsr_template_hfx => matrix_s(1)%matrix 846 END IF 847 848 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env) 849 CALL cp_fm_create(work_matrices%hfx_fm_ao_ao, fm_struct) 850 CALL cp_fm_struct_release(fm_struct) 851 852 CALL dbcsr_allocate_matrix_set(work_matrices%hfx_rho_ao, nspins) 853 DO ispin = 1, nspins 854 CALL dbcsr_init_p(work_matrices%hfx_rho_ao(ispin)%matrix) 855 CALL tddfpt_dbcsr_create_by_dist(work_matrices%hfx_rho_ao(ispin)%matrix, & 856 template=dbcsr_template_hfx, dbcsr_dist=dbcsr_dist, sab=sab_hfx) 857 END DO 858 859 CALL dbcsr_allocate_matrix_set(work_matrices%hfx_hmat, nspins) 860 DO ispin = 1, nspins 861 CALL dbcsr_init_p(work_matrices%hfx_hmat(ispin)%matrix) 862 CALL tddfpt_dbcsr_create_by_dist(work_matrices%hfx_hmat(ispin)%matrix, & 863 template=dbcsr_template_hfx, dbcsr_dist=dbcsr_dist, sab=sab_hfx) 864 END DO 865 END IF 866 867 CALL timestop(handle) 868 END SUBROUTINE tddfpt_create_work_matrices 869 870! ************************************************************************************************** 871!> \brief Release work matrices. 872!> \param work_matrices work matrices (destroyed on exit) 873!> \param sub_env parallel group environment 874!> \par History 875!> * 02.2017 created [Sergey Chulkov] 876! ************************************************************************************************** 877 SUBROUTINE tddfpt_release_work_matrices(work_matrices, sub_env) 878 TYPE(tddfpt_work_matrices), INTENT(inout) :: work_matrices 879 TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env 880 881 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_release_work_matrices', & 882 routineP = moduleN//':'//routineN 883 884 INTEGER :: handle, ispin, istate 885 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 886 887 CALL timeset(routineN, handle) 888 889 ! HFX-ralated matrices 890 IF (ASSOCIATED(work_matrices%hfx_hmat)) THEN 891 DO ispin = SIZE(work_matrices%hfx_hmat), 1, -1 892 CALL dbcsr_deallocate_matrix(work_matrices%hfx_hmat(ispin)%matrix) 893 END DO 894 DEALLOCATE (work_matrices%hfx_hmat) 895 END IF 896 897 IF (ASSOCIATED(work_matrices%hfx_rho_ao)) THEN 898 DO ispin = SIZE(work_matrices%hfx_rho_ao), 1, -1 899 CALL dbcsr_deallocate_matrix(work_matrices%hfx_rho_ao(ispin)%matrix) 900 END DO 901 DEALLOCATE (work_matrices%hfx_rho_ao) 902 END IF 903 904 IF (ASSOCIATED(work_matrices%hfx_fm_ao_ao)) & 905 CALL cp_fm_release(work_matrices%hfx_fm_ao_ao) 906 907 ! real-space and reciprocal-space grids 908 CALL pw_env_get(sub_env%pw_env, auxbas_pw_pool=auxbas_pw_pool) 909 DO ispin = SIZE(work_matrices%wpw_rspace_sub), 1, -1 910 CALL pw_pool_give_back_pw(auxbas_pw_pool, work_matrices%wpw_rspace_sub(ispin)%pw) 911 CALL pw_pool_give_back_pw(auxbas_pw_pool, work_matrices%wpw_gspace_sub(ispin)%pw) 912 CALL pw_pool_give_back_pw(auxbas_pw_pool, work_matrices%A_ia_rspace_sub(ispin)%pw) 913 END DO 914 DEALLOCATE (work_matrices%A_ia_rspace_sub, work_matrices%wpw_gspace_sub, work_matrices%wpw_rspace_sub) 915 916 IF (ASSOCIATED(work_matrices%rho_aux_fit_struct_sub)) & 917 CALL qs_rho_release(work_matrices%rho_aux_fit_struct_sub) 918 CALL qs_rho_release(work_matrices%rho_orb_struct_sub) 919 920 DO ispin = SIZE(work_matrices%A_ia_munu_sub), 1, -1 921 CALL dbcsr_deallocate_matrix(work_matrices%A_ia_munu_sub(ispin)%matrix) 922 END DO 923 DEALLOCATE (work_matrices%A_ia_munu_sub) 924 925 IF (ASSOCIATED(work_matrices%wfm_aux_orb_sub)) & 926 CALL cp_fm_release(work_matrices%wfm_aux_orb_sub) 927 IF (ASSOCIATED(work_matrices%rho_ao_aux_fit_fm_sub)) & 928 CALL cp_fm_release(work_matrices%rho_ao_aux_fit_fm_sub) 929 CALL cp_fm_release(work_matrices%rho_ao_orb_fm_sub) 930 931 IF (ALLOCATED(work_matrices%evects_sub)) THEN 932 DO istate = SIZE(work_matrices%evects_sub, 2), 1, -1 933 DO ispin = SIZE(work_matrices%evects_sub, 1), 1, -1 934 CALL cp_fm_release(work_matrices%Aop_evects_sub(ispin, istate)%matrix) 935 CALL cp_fm_release(work_matrices%evects_sub(ispin, istate)%matrix) 936 END DO 937 END DO 938 DEALLOCATE (work_matrices%Aop_evects_sub, work_matrices%evects_sub) 939 END IF 940 941 DO ispin = SIZE(work_matrices%fm_pool_ao_mo_occ), 1, -1 942 CALL cp_fm_release(work_matrices%wfm_mo_virt_mo_occ(ispin)%matrix) 943 CALL cp_fm_release(work_matrices%wfm_mo_occ_mo_occ(ispin)%matrix) 944 CALL cp_fm_release(work_matrices%wfm_ao_mo_occ(ispin)%matrix) 945 CALL cp_fm_release(work_matrices%S_C0(ispin)%matrix) 946 CALL cp_fm_release(work_matrices%S_C0_C0T(ispin)%matrix) 947 END DO 948 DEALLOCATE (work_matrices%S_C0, work_matrices%S_C0_C0T, work_matrices%wfm_ao_mo_occ) 949 DEALLOCATE (work_matrices%wfm_mo_occ_mo_occ, work_matrices%wfm_mo_virt_mo_occ) 950 951 DO ispin = SIZE(work_matrices%fm_pool_ao_mo_occ), 1, -1 952 CALL fm_pool_release(work_matrices%fm_pool_ao_mo_occ(ispin)%pool) 953 END DO 954 DEALLOCATE (work_matrices%fm_pool_ao_mo_occ) 955 956 CALL timestop(handle) 957 END SUBROUTINE tddfpt_release_work_matrices 958 959! ************************************************************************************************** 960!> \brief Create kernel environment. 961!> \param kernel_env kernel environment (allocated and initialised on exit) 962!> \param rho_struct_sub ground state charge density 963!> \param xc_section input section which defines an exchange-correlation functional 964!> \param is_rks_triplets indicates that the triplet excited states calculation using 965!> spin-unpolarised molecular orbitals has been requested 966!> \param sub_env parallel group environment 967!> \par History 968!> * 02.2017 created [Sergey Chulkov] 969!> * 06.2018 the charge density needs to be provided via a dummy argument [Sergey Chulkov] 970! ************************************************************************************************** 971 SUBROUTINE tddfpt_create_kernel_env(kernel_env, rho_struct_sub, xc_section, is_rks_triplets, sub_env) 972 TYPE(tddfpt_kernel_env_type), INTENT(out) :: kernel_env 973 TYPE(qs_rho_type), POINTER :: rho_struct_sub 974 TYPE(section_vals_type), POINTER :: xc_section 975 LOGICAL, INTENT(in) :: is_rks_triplets 976 TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env 977 978 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_create_kernel_env', & 979 routineP = moduleN//':'//routineN 980 981 INTEGER :: handle, ispin, nao, nspins 982 INTEGER, DIMENSION(maxspins) :: nmo_occ 983 LOGICAL :: lsd 984 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_ij_r, rho_ij_r2, tau_ij_r, tau_ij_r2 985 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 986 TYPE(section_vals_type), POINTER :: xc_fun_section 987 988 CALL timeset(routineN, handle) 989 990 nspins = SIZE(sub_env%mos_occ) 991 lsd = (nspins > 1) .OR. is_rks_triplets 992 993 DO ispin = 1, nspins 994 CALL cp_fm_get_info(sub_env%mos_occ(ispin)%matrix, nrow_global=nao, ncol_global=nmo_occ(ispin)) 995 END DO 996 997 CALL pw_env_get(sub_env%pw_env, auxbas_pw_pool=auxbas_pw_pool) 998 999 CALL qs_rho_get(rho_struct_sub, rho_r=rho_ij_r, tau_r=tau_ij_r) 1000 1001 NULLIFY (kernel_env%xc_rho_set, kernel_env%xc_rho1_set, kernel_env%xc_deriv_set) 1002 1003 IF (is_rks_triplets) THEN 1004 ! we are about to compute triplet states using spin-restricted reference MOs; 1005 ! we still need the beta-spin density component in order to compute the TDDFT kernel 1006 ALLOCATE (rho_ij_r2(2)) 1007 rho_ij_r2(1)%pw => rho_ij_r(1)%pw 1008 rho_ij_r2(2)%pw => rho_ij_r(1)%pw 1009 1010 IF (ASSOCIATED(tau_ij_r)) THEN 1011 ALLOCATE (tau_ij_r2(2)) 1012 tau_ij_r2(1)%pw => tau_ij_r(1)%pw 1013 tau_ij_r2(2)%pw => tau_ij_r(1)%pw 1014 END IF 1015 1016 CALL xc_prep_2nd_deriv(kernel_env%xc_deriv_set, kernel_env%xc_rho_set, rho_ij_r2, & 1017 auxbas_pw_pool, xc_section=xc_section, tau_r=tau_ij_r2) 1018 1019 IF (ASSOCIATED(tau_ij_r)) & 1020 DEALLOCATE (tau_ij_r2) 1021 1022 DEALLOCATE (rho_ij_r2) 1023 ELSE 1024 CALL xc_prep_2nd_deriv(kernel_env%xc_deriv_set, kernel_env%xc_rho_set, rho_ij_r, & 1025 auxbas_pw_pool, xc_section=xc_section, tau_r=tau_ij_r) 1026 END IF 1027 1028 ! ++ allocate structure for response density 1029 kernel_env%xc_section => xc_section 1030 kernel_env%deriv_method_id = section_get_ival(xc_section, "XC_GRID%XC_DERIV") 1031 kernel_env%rho_smooth_id = section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO") 1032 1033 xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL") 1034 kernel_env%xc_rho1_cflags = xc_functionals_get_needs(functionals=xc_fun_section, lsd=lsd, add_basic_components=.TRUE.) 1035 1036 CALL xc_rho_set_create(kernel_env%xc_rho1_set, auxbas_pw_pool%pw_grid%bounds_local, & 1037 rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), & 1038 drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), & 1039 tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF")) 1040 1041 kernel_env%alpha = 1.0_dp 1042 kernel_env%beta = 0.0_dp 1043 1044 ! kernel_env%beta is taken into account in spin-restricted case only 1045 IF (nspins == 1) THEN 1046 IF (is_rks_triplets) THEN 1047 ! K_{triplets} = K_{alpha,alpha} - K_{alpha,beta} 1048 kernel_env%beta = -1.0_dp 1049 ELSE 1050 ! alpha beta 1051 ! K_{singlets} = K_{alpha,alpha} + K_{alpha,beta} = 2 * K_{alpha,alpha} + 0 * K_{alpha,beta}, 1052 ! due to the following relation : K_{alpha,alpha,singlets} == K_{alpha,beta,singlets} 1053 kernel_env%alpha = 2.0_dp 1054 END IF 1055 END IF 1056 1057 CALL timestop(handle) 1058 END SUBROUTINE tddfpt_create_kernel_env 1059 1060! ************************************************************************************************** 1061!> \brief Release kernel environment. 1062!> \param kernel_env kernel environment (destroyed on exit) 1063!> \par History 1064!> * 02.2017 created [Sergey Chulkov] 1065! ************************************************************************************************** 1066 SUBROUTINE tddfpt_release_kernel_env(kernel_env) 1067 TYPE(tddfpt_kernel_env_type), INTENT(inout) :: kernel_env 1068 1069 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_release_kernel_env', & 1070 routineP = moduleN//':'//routineN 1071 1072 INTEGER :: handle 1073 1074 CALL timeset(routineN, handle) 1075 1076 CALL xc_rho_set_release(kernel_env%xc_rho1_set) 1077 CALL xc_dset_release(kernel_env%xc_deriv_set) 1078 CALL xc_rho_set_release(kernel_env%xc_rho_set) 1079 1080 CALL timestop(handle) 1081 END SUBROUTINE tddfpt_release_kernel_env 1082 1083! ************************************************************************************************** 1084!> \brief Generate all virtual molecular orbitals for a given spin by diagonalising 1085!> the corresponding Kohn-Sham matrix. 1086!> \param gs_mos structure to store occupied and virtual molecular orbitals 1087!> (allocated and initialised on exit) 1088!> \param mo_set ground state molecular orbitals for a given spin 1089!> \param nlumo number of unoccupied states to consider (-1 means all states) 1090!> \param blacs_env BLACS parallel environment 1091!> \param cholesky_method Cholesky method to compute the inverse overlap matrix 1092!> \param matrix_ks Kohn-Sham matrix for a given spin 1093!> \param matrix_s overlap matrix 1094!> \param mos_virt precomputed (OT) expansion coefficients of virtual molecular orbitals 1095!> (in addition to the ADDED_MOS, if present) 1096!> \param evals_virt orbital energies of precomputed (OT) virtual molecular orbitals 1097!> \par History 1098!> * 05.2016 created as tddfpt_lumos() [Sergey Chulkov] 1099!> * 06.2016 renamed, altered prototype [Sergey Chulkov] 1100! ************************************************************************************************** 1101 SUBROUTINE tddfpt_init_ground_state_mos(gs_mos, mo_set, nlumo, blacs_env, cholesky_method, matrix_ks, matrix_s, & 1102 mos_virt, evals_virt) 1103 TYPE(tddfpt_ground_state_mos), INTENT(out) :: gs_mos 1104 TYPE(mo_set_type), POINTER :: mo_set 1105 INTEGER, INTENT(in) :: nlumo 1106 TYPE(cp_blacs_env_type), POINTER :: blacs_env 1107 INTEGER, INTENT(in) :: cholesky_method 1108 TYPE(dbcsr_type), POINTER :: matrix_ks, matrix_s 1109 TYPE(cp_fm_type), POINTER :: mos_virt 1110 REAL(kind=dp), DIMENSION(:), POINTER :: evals_virt 1111 1112 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_init_ground_state_mos', & 1113 routineP = moduleN//':'//routineN 1114 REAL(kind=dp), PARAMETER :: eps_dp = EPSILON(0.0_dp) 1115 1116 INTEGER :: cholesky_method_inout, handle, icol_global, icol_local, imo, irow_global, & 1117 irow_local, nao, ncol_local, nelectrons, nmo_occ, nmo_scf, nmo_virt, nrow_local, sign_int 1118 INTEGER, ALLOCATABLE, DIMENSION(:) :: minrow_neg_array, minrow_pos_array, & 1119 sum_sign_array 1120 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 1121 LOGICAL :: do_eigen 1122 REAL(kind=dp) :: element, maxocc 1123 REAL(kind=dp), DIMENSION(:), POINTER :: mo_evals_extended, mo_occ_extended, & 1124 mo_occ_scf 1125 REAL(KIND=dp), DIMENSION(:, :), POINTER :: my_block 1126 TYPE(cp_fm_pool_type), POINTER :: wfn_fm_pool 1127 TYPE(cp_fm_struct_type), POINTER :: ao_ao_fm_struct, ao_mo_occ_fm_struct, & 1128 ao_mo_virt_fm_struct, wfn_fm_struct 1129 TYPE(cp_fm_type), POINTER :: matrix_ks_fm, mo_coeff_extended, & 1130 ortho_fm, work_fm 1131 TYPE(cp_para_env_type), POINTER :: para_env 1132 TYPE(mo_set_type), POINTER :: mos_extended 1133 1134 CALL timeset(routineN, handle) 1135 1136 CALL get_blacs_info(blacs_env, para_env=para_env) 1137 1138 CALL get_mo_set(mo_set, nao=nao, nmo=nmo_scf, homo=nmo_occ, maxocc=maxocc, & 1139 nelectron=nelectrons, occupation_numbers=mo_occ_scf) 1140 1141 nmo_virt = nao - nmo_occ 1142 IF (nlumo >= 0) & 1143 nmo_virt = MIN(nmo_virt, nlumo) 1144 1145 IF (nmo_virt <= 0) & 1146 CALL cp_abort(__LOCATION__, & 1147 'At least one unoccupied molecular orbital is required to calculate excited states.') 1148 1149 do_eigen = .FALSE. 1150 ! diagonalise the Kohn-Sham matrix one more time if the number of available unoccupied states are too small 1151 IF (ASSOCIATED(evals_virt)) THEN 1152 CPASSERT(ASSOCIATED(mos_virt)) 1153 IF (nmo_virt > nmo_scf - nmo_occ + SIZE(evals_virt)) do_eigen = .TRUE. 1154 ELSE 1155 IF (nmo_virt > nmo_scf - nmo_occ) do_eigen = .TRUE. 1156 END IF 1157 1158 ! ++ allocate storage space for gs_mos 1159 NULLIFY (ao_mo_occ_fm_struct, ao_mo_virt_fm_struct) 1160 CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env) 1161 CALL cp_fm_struct_create(ao_mo_virt_fm_struct, nrow_global=nao, ncol_global=nmo_virt, context=blacs_env) 1162 1163 NULLIFY (gs_mos%mos_occ, gs_mos%mos_virt) 1164 CALL cp_fm_create(gs_mos%mos_occ, ao_mo_occ_fm_struct) 1165 CALL cp_fm_create(gs_mos%mos_virt, ao_mo_virt_fm_struct) 1166 1167 ALLOCATE (gs_mos%evals_occ(nmo_occ)) 1168 ALLOCATE (gs_mos%evals_virt(nmo_virt)) 1169 ALLOCATE (gs_mos%phases_occ(nmo_occ)) 1170 1171 CALL cp_fm_struct_release(ao_mo_virt_fm_struct) 1172 CALL cp_fm_struct_release(ao_mo_occ_fm_struct) 1173 1174 ! ++ nullify pointers 1175 NULLIFY (ao_ao_fm_struct, wfn_fm_struct, wfn_fm_pool) 1176 NULLIFY (mos_extended, mo_coeff_extended, mo_evals_extended, mo_occ_extended) 1177 CALL cp_fm_struct_create(ao_ao_fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env) 1178 1179 IF (do_eigen) THEN 1180 ! ++ set of molecular orbitals 1181 CALL cp_fm_struct_create(wfn_fm_struct, nrow_global=nao, ncol_global=nmo_occ + nmo_virt, context=blacs_env) 1182 CALL fm_pool_create(wfn_fm_pool, wfn_fm_struct) 1183 1184 CALL allocate_mo_set(mos_extended, nao, nmo_occ + nmo_virt, nelectrons, & 1185 REAL(nelectrons, dp), maxocc, flexible_electron_count=0.0_dp) 1186 CALL init_mo_set(mos_extended, fm_pool=wfn_fm_pool, name="mos-extended") 1187 CALL fm_pool_release(wfn_fm_pool) 1188 CALL cp_fm_struct_release(wfn_fm_struct) 1189 CALL get_mo_set(mos_extended, mo_coeff=mo_coeff_extended, & 1190 eigenvalues=mo_evals_extended, occupation_numbers=mo_occ_extended) 1191 1192 ! use the explicit loop in order to avoid temporary arrays. 1193 ! 1194 ! The assignment statement : mo_occ_extended(1:nmo_scf) = mo_occ_scf(1:nmo_scf) 1195 ! implies temporary arrays as a compiler does not know in advance that the pointers 1196 ! on both sides of the statement point to non-overlapped memory regions 1197 DO imo = 1, nmo_scf 1198 mo_occ_extended(imo) = mo_occ_scf(imo) 1199 END DO 1200 mo_occ_extended(nmo_scf + 1:) = 0.0_dp 1201 1202 ! ++ allocate temporary matrices 1203 NULLIFY (matrix_ks_fm, ortho_fm, work_fm) 1204 CALL cp_fm_create(matrix_ks_fm, ao_ao_fm_struct) 1205 CALL cp_fm_create(ortho_fm, ao_ao_fm_struct) 1206 CALL cp_fm_create(work_fm, ao_ao_fm_struct) 1207 1208 CALL cp_fm_struct_release(ao_ao_fm_struct) 1209 1210 ! some stuff from the subroutine general_eigenproblem() 1211 CALL copy_dbcsr_to_fm(matrix_s, ortho_fm) 1212 CALL copy_dbcsr_to_fm(matrix_ks, matrix_ks_fm) 1213 1214 IF (cholesky_method == cholesky_dbcsr) THEN 1215 CPABORT('CHOLESKY DBCSR_INVERSE is not implemented in TDDFT.') 1216 ELSE IF (cholesky_method == cholesky_off) THEN 1217 CPABORT('CHOLESKY OFF is not implemented in TDDFT.') 1218 ELSE 1219 CALL cp_fm_cholesky_decompose(ortho_fm) 1220 IF (cholesky_method == cholesky_inverse) THEN 1221 CALL cp_fm_triangular_invert(ortho_fm) 1222 END IF 1223 1224 ! need to store 'cholesky_method' in a temporary variable, as the subroutine eigensolver() will update this variable 1225 cholesky_method_inout = cholesky_method 1226 CALL eigensolver(matrix_ks_fm=matrix_ks_fm, mo_set=mos_extended, ortho=ortho_fm, & 1227 work=work_fm, cholesky_method=cholesky_method_inout, & 1228 do_level_shift=.FALSE., level_shift=0.0_dp, matrix_u_fm=null(), use_jacobi=.FALSE.) 1229 END IF 1230 1231 ! -- clean up needless matrices 1232 CALL cp_fm_release(work_fm) 1233 CALL cp_fm_release(ortho_fm) 1234 CALL cp_fm_release(matrix_ks_fm) 1235 ELSE 1236 CALL get_mo_set(mo_set, mo_coeff=mo_coeff_extended, & 1237 eigenvalues=mo_evals_extended, occupation_numbers=mo_occ_extended) 1238 END IF 1239 1240 ! compute the phase of molecular orbitals; 1241 ! matrix work_fm holds occupied molecular orbital coefficients distributed among all the processors 1242 CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env) 1243 CALL cp_fm_create(work_fm, ao_mo_occ_fm_struct) 1244 CALL cp_fm_struct_release(ao_mo_occ_fm_struct) 1245 1246 CALL cp_fm_to_fm(mo_coeff_extended, work_fm, ncol=nmo_occ, source_start=1, target_start=1) 1247 CALL cp_fm_get_info(work_fm, nrow_local=nrow_local, ncol_local=ncol_local, & 1248 row_indices=row_indices, col_indices=col_indices, local_data=my_block) 1249 1250 ALLOCATE (minrow_neg_array(nmo_occ), minrow_pos_array(nmo_occ), sum_sign_array(nmo_occ)) 1251 minrow_neg_array(:) = nao 1252 minrow_pos_array(:) = nao 1253 sum_sign_array(:) = 0 1254 DO icol_local = 1, ncol_local 1255 icol_global = col_indices(icol_local) 1256 1257 DO irow_local = 1, nrow_local 1258 element = my_block(irow_local, icol_local) 1259 1260 sign_int = 0 1261 IF (element >= eps_dp) THEN 1262 sign_int = 1 1263 ELSE IF (element <= -eps_dp) THEN 1264 sign_int = -1 1265 END IF 1266 1267 sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int 1268 1269 irow_global = row_indices(irow_local) 1270 IF (sign_int > 0) THEN 1271 IF (minrow_pos_array(icol_global) > irow_global) & 1272 minrow_pos_array(icol_global) = irow_global 1273 ELSE IF (sign_int < 0) THEN 1274 IF (minrow_neg_array(icol_global) > irow_global) & 1275 minrow_neg_array(icol_global) = irow_global 1276 END IF 1277 END DO 1278 END DO 1279 1280 CALL mp_sum(sum_sign_array, para_env%group) 1281 CALL mp_min(minrow_neg_array, para_env%group) 1282 CALL mp_min(minrow_pos_array, para_env%group) 1283 1284 DO icol_local = 1, nmo_occ 1285 IF (sum_sign_array(icol_local) > 0) THEN 1286 ! most of the expansion coefficients are positive => MO's phase = +1 1287 gs_mos%phases_occ(icol_local) = 1.0_dp 1288 ELSE IF (sum_sign_array(icol_local) < 0) THEN 1289 ! most of the expansion coefficients are negative => MO's phase = -1 1290 gs_mos%phases_occ(icol_local) = -1.0_dp 1291 ELSE 1292 ! equal number of positive and negative expansion coefficients 1293 IF (minrow_pos_array(icol_local) <= minrow_neg_array(icol_local)) THEN 1294 ! the first positive expansion coefficient has a lower index then 1295 ! the first negative expansion coefficient; MO's phase = +1 1296 gs_mos%phases_occ(icol_local) = 1.0_dp 1297 ELSE 1298 ! MO's phase = -1 1299 gs_mos%phases_occ(icol_local) = -1.0_dp 1300 END IF 1301 END IF 1302 END DO 1303 1304 DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array) 1305 CALL cp_fm_release(work_fm) 1306 1307 ! return the requested occupied and virtual molecular orbitals and corresponding orbital energies 1308 CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_occ, ncol=nmo_occ, source_start=1, target_start=1) 1309 gs_mos%evals_occ(1:nmo_occ) = mo_evals_extended(1:nmo_occ) 1310 1311 IF (ASSOCIATED(evals_virt) .AND. (.NOT. do_eigen) .AND. nmo_virt > nmo_scf - nmo_occ) THEN 1312 CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_scf - nmo_occ, & 1313 source_start=nmo_occ + 1, target_start=1) 1314 CALL cp_fm_to_fm(mos_virt, gs_mos%mos_virt, ncol=nmo_virt - (nmo_scf - nmo_occ), & 1315 source_start=1, target_start=nmo_scf - nmo_occ + 1) 1316 gs_mos%evals_virt(1:nmo_scf - nmo_occ) = evals_virt(nmo_occ + 1:nmo_occ + nmo_scf) 1317 gs_mos%evals_virt(nmo_scf - nmo_occ + 1:nmo_virt) = evals_virt(1:nmo_virt - (nmo_scf - nmo_occ)) 1318 ELSE 1319 CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_virt, source_start=nmo_occ + 1, target_start=1) 1320 gs_mos%evals_virt(1:nmo_virt) = mo_evals_extended(nmo_occ + 1:nmo_occ + nmo_virt) 1321 END IF 1322 1323 IF (do_eigen) & 1324 CALL deallocate_mo_set(mos_extended) 1325 CALL timestop(handle) 1326 END SUBROUTINE tddfpt_init_ground_state_mos 1327 1328! ************************************************************************************************** 1329!> \brief Release molecular orbitals. 1330!> \param gs_mos structure that holds occupied and virtual molecular orbitals 1331!> \par History 1332!> * 06.2016 created [Sergey Chulkov] 1333! ************************************************************************************************** 1334 SUBROUTINE tddfpt_release_ground_state_mos(gs_mos) 1335 TYPE(tddfpt_ground_state_mos), INTENT(inout) :: gs_mos 1336 1337 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_release_ground_state_mos', & 1338 routineP = moduleN//':'//routineN 1339 1340 INTEGER :: handle 1341 1342 CALL timeset(routineN, handle) 1343 1344 IF (ALLOCATED(gs_mos%phases_occ)) & 1345 DEALLOCATE (gs_mos%phases_occ) 1346 1347 IF (ALLOCATED(gs_mos%evals_virt)) & 1348 DEALLOCATE (gs_mos%evals_virt) 1349 1350 IF (ALLOCATED(gs_mos%evals_occ)) & 1351 DEALLOCATE (gs_mos%evals_occ) 1352 1353 IF (ASSOCIATED(gs_mos%mos_virt)) & 1354 CALL cp_fm_release(gs_mos%mos_virt) 1355 1356 IF (ASSOCIATED(gs_mos%mos_occ)) & 1357 CALL cp_fm_release(gs_mos%mos_occ) 1358 1359 CALL timestop(handle) 1360 END SUBROUTINE tddfpt_release_ground_state_mos 1361 1362! ************************************************************************************************** 1363!> \brief Compute the number of possible singly excited states (occ -> virt) 1364!> \param gs_mos occupied and virtual molecular orbitals optimised for the ground state 1365!> \return the number of possible single excitations 1366!> \par History 1367!> * 01.2017 created [Sergey Chulkov] 1368! ************************************************************************************************** 1369 PURE FUNCTION tddfpt_total_number_of_states(gs_mos) RESULT(nstates_total) 1370 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 1371 INTENT(in) :: gs_mos 1372 INTEGER(kind=int_8) :: nstates_total 1373 1374 INTEGER :: ispin, nspins 1375 1376 nstates_total = 0 1377 nspins = SIZE(gs_mos) 1378 1379 DO ispin = 1, nspins 1380 nstates_total = nstates_total + & 1381 SIZE(gs_mos(ispin)%evals_occ, kind=int_8)* & 1382 SIZE(gs_mos(ispin)%evals_virt, kind=int_8) 1383 END DO 1384 END FUNCTION tddfpt_total_number_of_states 1385 1386! ************************************************************************************************** 1387!> \brief Generate missed guess vectors. 1388!> \param evects guess vectors distributed across all processors (initialised on exit) 1389!> \param evals guessed transition energies (initialised on exit) 1390!> \param gs_mos occupied and virtual molecular orbitals optimised for the ground state 1391!> \param log_unit output unit 1392!> \par History 1393!> * 05.2016 created as tddfpt_guess() [Sergey Chulkov] 1394!> * 06.2016 renamed, altered prototype, supports spin-polarised density [Sergey Chulkov] 1395!> * 01.2017 simplified prototype, do not compute all possible singly-excited states 1396!> [Sergey Chulkov] 1397!> \note \parblock 1398!> Based on the subroutine co_initial_guess() which was originally created by 1399!> Thomas Chassaing on 06.2003. 1400!> 1401!> Only not associated guess vectors 'evects(spin, state)%matrix' are allocated and 1402!> initialised; associated vectors assumed to be initialised elsewhere (e.g. using 1403!> a restart file). 1404!> \endparblock 1405! ************************************************************************************************** 1406 SUBROUTINE tddfpt_guess_vectors(evects, evals, gs_mos, log_unit) 1407 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(inout) :: evects 1408 REAL(kind=dp), DIMENSION(:), INTENT(inout) :: evals 1409 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 1410 INTENT(in) :: gs_mos 1411 INTEGER, INTENT(in) :: log_unit 1412 1413 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_guess_vectors', & 1414 routineP = moduleN//':'//routineN 1415 1416 CHARACTER(len=5) :: spin_label 1417 INTEGER :: handle, imo_occ, imo_virt, ind, ispin, & 1418 istate, jspin, nspins, nstates, & 1419 nstates_occ_virt_alpha, & 1420 nstates_selected 1421 INTEGER, ALLOCATABLE, DIMENSION(:) :: inds 1422 INTEGER, DIMENSION(maxspins) :: nmo_occ_avail, nmo_occ_selected, & 1423 nmo_virt_selected 1424 REAL(kind=dp) :: e_occ 1425 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: e_virt_minus_occ 1426 TYPE(cp_fm_struct_p_type), DIMENSION(maxspins) :: fm_struct_evects 1427 1428 CALL timeset(routineN, handle) 1429 1430 nspins = SIZE(evects, 1) 1431 nstates = SIZE(evects, 2) 1432 1433 IF (debug_this_module) THEN 1434 CPASSERT(nstates > 0) 1435 CPASSERT(nspins == 1 .OR. nspins == 2) 1436 CPASSERT(SIZE(gs_mos) == nspins) 1437 END IF 1438 1439 DO ispin = 1, nspins 1440 ! number of occupied orbitals for each spin component 1441 nmo_occ_avail(ispin) = SIZE(gs_mos(ispin)%evals_occ) 1442 ! number of occupied and virtual orbitals which can potentially 1443 ! contribute to the excited states in question. 1444 nmo_occ_selected(ispin) = MIN(nmo_occ_avail(ispin), nstates) 1445 nmo_virt_selected(ispin) = MIN(SIZE(gs_mos(ispin)%evals_virt), nstates) 1446 1447 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct_evects(ispin)%struct) 1448 END DO 1449 1450 ! TO DO: the variable 'nstates_selected' should probably be declared as INTEGER(kind=int_8), 1451 ! however we need a special version of the subroutine sort() in order to do so 1452 nstates_selected = DOT_PRODUCT(nmo_occ_selected(1:nspins), nmo_virt_selected(1:nspins)) 1453 1454 ALLOCATE (inds(nstates_selected)) 1455 ALLOCATE (e_virt_minus_occ(nstates_selected)) 1456 1457 istate = 0 1458 DO ispin = 1, nspins 1459 DO imo_occ = 1, nmo_occ_selected(ispin) 1460 ! Here imo_occ enumerate Occupied orbitals in inverse order (from the last to the first element) 1461 e_occ = gs_mos(ispin)%evals_occ(nmo_occ_avail(ispin) - imo_occ + 1) 1462 1463 DO imo_virt = 1, nmo_virt_selected(ispin) 1464 istate = istate + 1 1465 e_virt_minus_occ(istate) = gs_mos(ispin)%evals_virt(imo_virt) - e_occ 1466 END DO 1467 END DO 1468 END DO 1469 1470 IF (debug_this_module) THEN 1471 CPASSERT(istate == nstates_selected) 1472 END IF 1473 1474 CALL sort(e_virt_minus_occ, nstates_selected, inds) 1475 1476 IF (nspins == 1) THEN 1477 ispin = 1 1478 spin_label = ' ' 1479 END IF 1480 1481 nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(1) 1482 IF (log_unit > 0) THEN 1483 WRITE (log_unit, '(/,21X,A)') "TDDFPT initial guess" 1484 WRITE (log_unit, '(1X,60("-"))') 1485 WRITE (log_unit, '(5X,A)') "State Occupied -> Virtual Excitation" 1486 WRITE (log_unit, '(5X,A)') "number orbital orbital energy (eV)" 1487 WRITE (log_unit, '(1X,60("-"))') 1488 END IF 1489 1490 DO istate = 1, nstates 1491 IF (ASSOCIATED(evects(1, istate)%matrix)) THEN 1492 IF (log_unit > 0) & 1493 WRITE (log_unit, '(1X,I8,11X,A19,8X,F14.5)') & 1494 istate, "*** restarted ***", evals(istate)*evolt 1495 ELSE 1496 ind = inds(istate) - 1 1497 IF (nspins > 1) THEN 1498 IF (ind < nstates_occ_virt_alpha) THEN 1499 ispin = 1 1500 spin_label = '(alp)' 1501 ELSE 1502 ispin = 2 1503 ind = ind - nstates_occ_virt_alpha 1504 spin_label = '(bet)' 1505 END IF 1506 END IF 1507 1508 imo_occ = nmo_occ_avail(ispin) - ind/nmo_virt_selected(ispin) 1509 imo_virt = MOD(ind, nmo_virt_selected(ispin)) + 1 1510 evals(istate) = e_virt_minus_occ(istate) 1511 1512 IF (log_unit > 0) & 1513 WRITE (log_unit, '(1X,I8,5X,I8,1X,A5,3X,I8,1X,A5,2X,F14.5)') & 1514 istate, imo_occ, spin_label, nmo_occ_avail(ispin) + imo_virt, spin_label, e_virt_minus_occ(istate)*evolt 1515 1516 DO jspin = 1, nspins 1517 ! .NOT. ASSOCIATED(evects(jspin, istate)%matrix)) 1518 CALL cp_fm_create(evects(jspin, istate)%matrix, fm_struct_evects(jspin)%struct) 1519 CALL cp_fm_set_all(evects(jspin, istate)%matrix, 0.0_dp) 1520 1521 IF (jspin == ispin) & 1522 CALL cp_fm_to_fm(gs_mos(ispin)%mos_virt, evects(ispin, istate)%matrix, & 1523 ncol=1, source_start=imo_virt, target_start=imo_occ) 1524 END DO 1525 END IF 1526 END DO 1527 1528 IF (log_unit > 0) & 1529 WRITE (log_unit, '(/,1X,A,T30,I24)') 'Number of active states:', tddfpt_total_number_of_states(gs_mos) 1530 1531 DEALLOCATE (e_virt_minus_occ) 1532 DEALLOCATE (inds) 1533 CALL timestop(handle) 1534 END SUBROUTINE tddfpt_guess_vectors 1535 1536! ************************************************************************************************** 1537!> \brief Make TDDFPT trial vectors orthogonal to all occupied molecular orbitals. 1538!> \param evects trial vectors distributed across all processors (modified on exit) 1539!> \param S_C0_C0T matrix product S * C_0 * C_0^T, where C_0 is the ground state 1540!> wave function for each spin expressed in atomic basis set, 1541!> and S is the corresponding overlap matrix 1542!> \param work_fm_ao_mo_occ work matrices with shape [nao x nmo_occ(spin)] to store matrix products 1543!> C_0 * C_0^T * S * evects (modified on exit) 1544!> \par History 1545!> * 05.2016 created [Sergey Chulkov] 1546!> \note Based on the subroutine p_preortho() which was created by Thomas Chassaing on 09.2002. 1547!> Should be useless when ground state MOs are computed with extremely high accuracy, 1548!> as all virtual orbitals are already orthogonal to the occupied ones by design. 1549!> However, when the norm of residual vectors is relatively small (e.g. less then SCF_EPS), 1550!> new Krylov's vectors seem to be random and should be orthogonalised even with respect to 1551!> the occupied MOs. 1552! ************************************************************************************************** 1553 SUBROUTINE tddfpt_orthogonalize_psi1_psi0(evects, S_C0_C0T, work_fm_ao_mo_occ) 1554 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: evects 1555 TYPE(cp_fm_p_type), DIMENSION(:), INTENT(in) :: S_C0_C0T, work_fm_ao_mo_occ 1556 1557 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_orthogonalize_psi1_psi0', & 1558 routineP = moduleN//':'//routineN 1559 1560 INTEGER :: handle, ispin, ivect, nao, nspins, nvects 1561 INTEGER, DIMENSION(maxspins) :: nmo_occ 1562 1563 CALL timeset(routineN, handle) 1564 1565 nspins = SIZE(evects, 1) 1566 nvects = SIZE(evects, 2) 1567 1568 IF (debug_this_module) THEN 1569 CPASSERT(nspins > 0 .AND. nspins <= maxspins) 1570 CPASSERT(SIZE(S_C0_C0T) == nspins) 1571 CPASSERT(SIZE(work_fm_ao_mo_occ) == nspins) 1572 END IF 1573 1574 CALL cp_fm_get_info(matrix=work_fm_ao_mo_occ(1)%matrix, nrow_global=nao, ncol_global=nmo_occ(1)) 1575 DO ispin = 2, nspins 1576 CALL cp_fm_get_info(matrix=work_fm_ao_mo_occ(ispin)%matrix, ncol_global=nmo_occ(ispin)) 1577 END DO 1578 1579 IF (nvects > 0) THEN 1580 1581 DO ivect = 1, nvects 1582 DO ispin = 1, nspins 1583 ! work_fm_ao_mo_occ: C0 * C0^T * S * C1 == (S * C0 * C0^T)^T * C1 1584 CALL cp_gemm('T', 'N', nao, nmo_occ(ispin), nao, 1.0_dp, S_C0_C0T(ispin)%matrix, & 1585 evects(ispin, ivect)%matrix, 0.0_dp, work_fm_ao_mo_occ(ispin)%matrix) 1586 1587 CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, ivect)%matrix, -1.0_dp, work_fm_ao_mo_occ(ispin)%matrix) 1588 END DO 1589 END DO 1590 END IF 1591 1592 CALL timestop(handle) 1593 END SUBROUTINE tddfpt_orthogonalize_psi1_psi0 1594 1595! ************************************************************************************************** 1596!> \brief Check that orthogonalised TDDFPT trial vectors remain orthogonal to 1597!> occupied molecular orbitals. 1598!> \param evects trial vectors 1599!> \param S_C0 matrix product S * C_0, where C_0 is the ground state wave function 1600!> for each spin in atomic basis set, and S is the corresponding overlap matrix 1601!> \param max_norm the largest possible overlap between the ground state and 1602!> excited state wave functions 1603!> \param work_fm_mo_occ_mo_occ work matrices with shape [nmo_occ(spin) x nmo_occ(spin)] 1604!> to store matrix products S_0^T * S * evects (modified on exit) 1605!> \return true if trial vectors are non-orthogonal to occupied molecular orbitals 1606!> \par History 1607!> * 07.2016 created [Sergey Chulkov] 1608! ************************************************************************************************** 1609 FUNCTION tddfpt_is_nonorthogonal_psi1_psi0(evects, S_C0, max_norm, work_fm_mo_occ_mo_occ) RESULT(is_nonortho) 1610 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: evects 1611 TYPE(cp_fm_p_type), DIMENSION(:), INTENT(in) :: S_C0 1612 REAL(kind=dp), INTENT(in) :: max_norm 1613 TYPE(cp_fm_p_type), DIMENSION(:), INTENT(in) :: work_fm_mo_occ_mo_occ 1614 LOGICAL :: is_nonortho 1615 1616 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_is_nonorthogonal_psi1_psi0', & 1617 routineP = moduleN//':'//routineN 1618 1619 INTEGER :: handle, ispin, ivect, nao, nspins, nvects 1620 INTEGER, DIMENSION(maxspins) :: nmo_occ 1621 REAL(kind=dp) :: maxabs_val 1622 1623 CALL timeset(routineN, handle) 1624 1625 nspins = SIZE(evects, 1) 1626 nvects = SIZE(evects, 2) 1627 1628 IF (debug_this_module) THEN 1629 CPASSERT(nspins > 0 .AND. nspins <= maxspins) 1630 CPASSERT(SIZE(S_C0) == nspins) 1631 CPASSERT(SIZE(work_fm_mo_occ_mo_occ) == nspins) 1632 END IF 1633 1634 CALL cp_fm_get_info(matrix=S_C0(1)%matrix, nrow_global=nao, ncol_global=nmo_occ(1)) 1635 DO ispin = 2, nspins 1636 CALL cp_fm_get_info(matrix=S_C0(ispin)%matrix, ncol_global=nmo_occ(ispin)) 1637 END DO 1638 1639 is_nonortho = .FALSE. 1640 1641 loop: DO ivect = 1, nvects 1642 DO ispin = 1, nspins 1643 ! work_fm_mo_occ_mo_occ = S_0^T * S * C_1 1644 CALL cp_gemm('T', 'N', nmo_occ(ispin), nmo_occ(ispin), nao, 1.0_dp, S_C0(ispin)%matrix, & 1645 evects(ispin, ivect)%matrix, 0.0_dp, work_fm_mo_occ_mo_occ(ispin)%matrix) 1646 1647 CALL cp_fm_maxabsval(work_fm_mo_occ_mo_occ(ispin)%matrix, maxabs_val) 1648 is_nonortho = maxabs_val > max_norm 1649 IF (is_nonortho) EXIT loop 1650 END DO 1651 END DO loop 1652 1653 CALL timestop(handle) 1654 END FUNCTION tddfpt_is_nonorthogonal_psi1_psi0 1655 1656! ************************************************************************************************** 1657!> \brief Make new TDDFPT trial vectors orthonormal to all previous trial vectors. 1658!> \param evects trial vectors (modified on exit) 1659!> \param nvects_new number of new trial vectors to orthogonalise 1660!> \param S_evects set of matrices to store matrix product S * evects (modified on exit) 1661!> \param matrix_s overlap matrix 1662!> \par History 1663!> * 05.2016 created [Sergey Chulkov] 1664!> * 02.2017 caching the matrix product S * evects [Sergey Chulkov] 1665!> \note \parblock 1666!> Based on the subroutines reorthogonalize() and normalize() which were originally created 1667!> by Thomas Chassaing on 03.2003. 1668!> 1669!> In order to orthogonalise a trial vector C3 = evects(:,3) with respect to previously 1670!> orthogonalised vectors C1 = evects(:,1) and C2 = evects(:,2) we need to compute the 1671!> quantity C3'' using the following formulae: 1672!> C3' = C3 - Tr(C3^T * S * C1) * C1, 1673!> C3'' = C3' - Tr(C3'^T * S * C2) * C2, 1674!> which can be expanded as: 1675!> C3'' = C3 - Tr(C3^T * S * C1) * C1 - Tr(C3^T * S * C2) * C2 + 1676!> Tr(C3^T * S * C1) * Tr(C2^T * S * C1) * C2 . 1677!> In case of unlimited float-point precision, the last term in above expression is exactly 0, 1678!> due to orthogonality condition between C1 and C2. In this case the expression could be 1679!> simplified as (taking into account the identity: Tr(A * S * B) = Tr(B * S * A)): 1680!> C3'' = C3 - Tr(C1^T * S * C3) * C1 - Tr(C2^T * S * C3) * C2 , 1681!> which means we do not need the variable S_evects to keep the matrix products S * Ci . 1682!> 1683!> In reality, however, we deal with limited float-point precision arithmetic meaning that 1684!> the trace Tr(C2^T * S * C1) is close to 0 but does not equal to 0 exactly. The term 1685!> Tr(C3^T * S * C1) * Tr(C2^T * S * C1) * C2 1686!> can not be ignored anymore. Ignorance of this term will lead to numerical instability 1687!> when the trace Tr(C3^T * S * C1) is large enough. 1688!> \endparblock 1689! ************************************************************************************************** 1690 SUBROUTINE tddfpt_orthonormalize_psi1_psi1(evects, nvects_new, S_evects, matrix_s) 1691 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: evects 1692 INTEGER, INTENT(in) :: nvects_new 1693 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: S_evects 1694 TYPE(dbcsr_type), POINTER :: matrix_s 1695 1696 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_orthonormalize_psi1_psi1', & 1697 routineP = moduleN//':'//routineN 1698 1699 INTEGER :: handle, ispin, ivect, jvect, nspins, & 1700 nvects_old, nvects_total 1701 INTEGER, DIMENSION(maxspins) :: nmo_occ 1702 REAL(kind=dp) :: norm 1703 REAL(kind=dp), DIMENSION(maxspins) :: weights 1704 1705 CALL timeset(routineN, handle) 1706 1707 nspins = SIZE(evects, 1) 1708 nvects_total = SIZE(evects, 2) 1709 nvects_old = nvects_total - nvects_new 1710 1711 IF (debug_this_module) THEN 1712 CPASSERT(SIZE(S_evects, 1) == nspins) 1713 CPASSERT(SIZE(S_evects, 2) == nvects_total) 1714 CPASSERT(nvects_old >= 0) 1715 END IF 1716 1717 DO ispin = 1, nspins 1718 CALL cp_fm_get_info(matrix=evects(ispin, 1)%matrix, ncol_global=nmo_occ(ispin)) 1719 END DO 1720 1721 DO jvect = nvects_old + 1, nvects_total 1722 ! <psi1_i | psi1_j> 1723 DO ivect = 1, jvect - 1 1724 CALL cp_fm_trace(evects(:, jvect), S_evects(:, ivect), weights(1:nspins)) 1725 norm = accurate_sum(weights(1:nspins)) 1726 1727 DO ispin = 1, nspins 1728 CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, jvect)%matrix, -norm, evects(ispin, ivect)%matrix) 1729 END DO 1730 END DO 1731 1732 ! <psi1_j | psi1_j> 1733 DO ispin = 1, nspins 1734 CALL cp_dbcsr_sm_fm_multiply(matrix_s, evects(ispin, jvect)%matrix, S_evects(ispin, jvect)%matrix, & 1735 ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp) 1736 END DO 1737 1738 CALL cp_fm_trace(evects(:, jvect), S_evects(:, jvect), weights(1:nspins)) 1739 1740 norm = accurate_sum(weights(1:nspins)) 1741 norm = 1.0_dp/SQRT(norm) 1742 1743 DO ispin = 1, nspins 1744 CALL cp_fm_scale(norm, evects(ispin, jvect)%matrix) 1745 CALL cp_fm_scale(norm, S_evects(ispin, jvect)%matrix) 1746 END DO 1747 END DO 1748 1749 CALL timestop(handle) 1750 END SUBROUTINE tddfpt_orthonormalize_psi1_psi1 1751 1752! ************************************************************************************************** 1753!> \brief Apply orbital energy difference term: 1754!> Aop_evects(spin,state) += KS(spin) * evects(spin,state) - 1755!> S * evects(spin,state) * diag(evals_occ(spin)) 1756!> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit) 1757!> \param evects trial vectors C_{1,i} 1758!> \param S_evects S * C_{1,i} 1759!> \param gs_mos molecular orbitals optimised for the ground state (only occupied orbital 1760!> energies [component %evals_occ] are needed) 1761!> \param matrix_ks Kohn-Sham matrix 1762!> \param work_fm_ao_mo_occ work matrices with shape [nao x nmo_occ(spin)] to store matrix 1763!> products S * evects * diag(evals_occ). 1764!> \par History 1765!> * 05.2016 initialise all matrix elements in one go [Sergey Chulkov] 1766!> * 03.2017 renamed from tddfpt_init_energy_diff(), altered prototype [Sergey Chulkov] 1767!> \note Based on the subroutine p_op_l1() which was originally created by 1768!> Thomas Chassaing on 08.2002. 1769! ************************************************************************************************** 1770 SUBROUTINE tddfpt_apply_energy_diff(Aop_evects, evects, S_evects, gs_mos, matrix_ks, work_fm_ao_mo_occ) 1771 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: Aop_evects, evects, S_evects 1772 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 1773 INTENT(in) :: gs_mos 1774 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_ks 1775 TYPE(cp_fm_p_type), DIMENSION(:), INTENT(in) :: work_fm_ao_mo_occ 1776 1777 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_energy_diff', & 1778 routineP = moduleN//':'//routineN 1779 1780 INTEGER :: handle, ispin, ivect, nspins, nvects 1781 INTEGER, DIMENSION(maxspins) :: nmo_occ 1782 1783 CALL timeset(routineN, handle) 1784 1785 nspins = SIZE(evects, 1) 1786 nvects = SIZE(evects, 2) 1787 1788 DO ispin = 1, nspins 1789 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ) 1790 END DO 1791 1792 DO ivect = 1, nvects 1793 DO ispin = 1, nspins 1794 CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, evects(ispin, ivect)%matrix, & 1795 Aop_evects(ispin, ivect)%matrix, ncol=nmo_occ(ispin), & 1796 alpha=1.0_dp, beta=1.0_dp) 1797 1798 CALL cp_fm_to_fm(S_evects(ispin, ivect)%matrix, work_fm_ao_mo_occ(ispin)%matrix) 1799 CALL cp_fm_column_scale(work_fm_ao_mo_occ(ispin)%matrix, gs_mos(ispin)%evals_occ) 1800 1801 ! KS * C1 - S * C1 * occupied_orbital_energies 1802 CALL cp_fm_scale_and_add(1.0_dp, Aop_evects(ispin, ivect)%matrix, -1.0_dp, work_fm_ao_mo_occ(ispin)%matrix) 1803 END DO 1804 END DO 1805 1806 CALL timestop(handle) 1807 END SUBROUTINE tddfpt_apply_energy_diff 1808 1809! ************************************************************************************************** 1810!> \brief Update v_rspace by adding coulomb term. 1811!> \param A_ia_rspace action of TDDFPT operator on the trial vector expressed in a plane wave 1812!> representation (modified on exit) 1813!> \param rho_ia_g response density in reciprocal space for the given trial vector 1814!> \param pw_env plain wave environment 1815!> \param work_v_gspace work reciprocal-space grid to store Coulomb potential (modified on exit) 1816!> \param work_v_rspace work real-space grid to store Coulomb potential (modified on exit) 1817!> \par History 1818!> * 05.2016 compute all coulomb terms in one go [Sergey Chulkov] 1819!> * 03.2017 proceed excited states sequentially; minimise the number of conversions between 1820!> DBCSR and FM matrices [Sergey Chulkov] 1821!> * 06.2018 return the action expressed in the plane wave representation instead of the one 1822!> in the atomic basis set representation 1823!> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by 1824!> Mohamed Fawzi on 10.2002. 1825! ************************************************************************************************** 1826 SUBROUTINE tddfpt_apply_coulomb(A_ia_rspace, rho_ia_g, pw_env, work_v_gspace, work_v_rspace) 1827 TYPE(pw_p_type), DIMENSION(:), POINTER :: A_ia_rspace 1828 TYPE(pw_type), POINTER :: rho_ia_g 1829 TYPE(pw_env_type), POINTER :: pw_env 1830 TYPE(pw_p_type), INTENT(inout) :: work_v_gspace, work_v_rspace 1831 1832 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_coulomb', & 1833 routineP = moduleN//':'//routineN 1834 1835 INTEGER :: handle, ispin, nspins 1836 REAL(kind=dp) :: alpha, pair_energy 1837 TYPE(pw_poisson_type), POINTER :: poisson_env 1838 1839 CALL timeset(routineN, handle) 1840 1841 nspins = SIZE(A_ia_rspace) 1842 CALL pw_env_get(pw_env, poisson_env=poisson_env) 1843 1844 IF (nspins > 1) THEN 1845 alpha = 1.0_dp 1846 ELSE 1847 ! spin-restricted case: alpha == 2 due to singlet state. 1848 ! In case of triplet states alpha == 0, so we should not call this subroutine at all. 1849 alpha = 2.0_dp 1850 END IF 1851 1852 CALL pw_poisson_solve(poisson_env, rho_ia_g, pair_energy, work_v_gspace%pw) 1853 CALL pw_transfer(work_v_gspace%pw, work_v_rspace%pw) 1854 1855 ! (i a || j b) = ( i_alpha a_alpha + i_beta a_beta || j_alpha b_alpha + j_beta b_beta) = 1856 ! tr (Cj_alpha^T * [J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu] * Cb_alpha) + 1857 ! tr (Cj_beta^T * [J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu] * Cb_beta) 1858 DO ispin = 1, nspins 1859 CALL pw_axpy(work_v_rspace%pw, A_ia_rspace(ispin)%pw, alpha) 1860 END DO 1861 1862 ! TO DO: tau component 1863 1864 CALL timestop(handle) 1865 END SUBROUTINE tddfpt_apply_coulomb 1866 1867! ************************************************************************************************** 1868!> \brief Update A_ia_munu by adding exchange-correlation term. 1869!> \param A_ia_rspace action of TDDFPT operator on trial vectors expressed in a plane wave 1870!> representation (modified on exit) 1871!> \param kernel_env kernel environment 1872!> \param rho_ia_struct response density for the given trial vector 1873!> \param is_rks_triplets indicates that the triplet excited states calculation using 1874!> spin-unpolarised molecular orbitals has been requested 1875!> \param pw_env plain wave environment 1876!> \param work_v_xc work real-space grid to store the gradient of the exchange-correlation 1877!> potential with respect to the response density (modified on exit) 1878!> \par History 1879!> * 05.2016 compute all kernel terms in one go [Sergey Chulkov] 1880!> * 03.2017 proceed excited states sequentially; minimise the number of conversions between 1881!> DBCSR and FM matrices [Sergey Chulkov] 1882!> * 06.2018 return the action expressed in the plane wave representation instead of the one 1883!> in the atomic basis set representation 1884!> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by 1885!> Mohamed Fawzi on 10.2002. 1886! ************************************************************************************************** 1887 SUBROUTINE tddfpt_apply_xc(A_ia_rspace, kernel_env, rho_ia_struct, is_rks_triplets, pw_env, work_v_xc) 1888 TYPE(pw_p_type), DIMENSION(:), POINTER :: A_ia_rspace 1889 TYPE(tddfpt_kernel_env_type), INTENT(in) :: kernel_env 1890 TYPE(qs_rho_type), POINTER :: rho_ia_struct 1891 LOGICAL, INTENT(in) :: is_rks_triplets 1892 TYPE(pw_env_type), POINTER :: pw_env 1893 TYPE(pw_p_type), DIMENSION(:), POINTER :: work_v_xc 1894 1895 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_xc', & 1896 routineP = moduleN//':'//routineN 1897 1898 INTEGER :: handle, ispin, nspins 1899 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_ia_g, rho_ia_g2, rho_ia_r, & 1900 rho_ia_r2, tau_ia_r, tau_ia_r2 1901 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 1902 1903 CALL timeset(routineN, handle) 1904 1905 nspins = SIZE(A_ia_rspace) 1906 CALL qs_rho_get(rho_ia_struct, rho_g=rho_ia_g, rho_r=rho_ia_r, tau_r=tau_ia_r) 1907 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 1908 1909 IF (debug_this_module) THEN 1910 CPASSERT(SIZE(rho_ia_g) == nspins) 1911 CPASSERT(SIZE(rho_ia_r) == nspins) 1912 CPASSERT((.NOT. ASSOCIATED(tau_ia_r)) .OR. SIZE(tau_ia_r) == nspins) 1913 CPASSERT((.NOT. is_rks_triplets) .OR. nspins == 1) 1914 END IF 1915 1916 NULLIFY (tau_ia_r2) 1917 IF (is_rks_triplets) THEN 1918 ALLOCATE (rho_ia_r2(2)) 1919 ALLOCATE (rho_ia_g2(2)) 1920 rho_ia_r2(1)%pw => rho_ia_r(1)%pw 1921 rho_ia_r2(2)%pw => rho_ia_r(1)%pw 1922 rho_ia_g2(1)%pw => rho_ia_g(1)%pw 1923 rho_ia_g2(2)%pw => rho_ia_g(1)%pw 1924 1925 IF (ASSOCIATED(tau_ia_r)) THEN 1926 ALLOCATE (tau_ia_r2(2)) 1927 tau_ia_r2(1)%pw => tau_ia_r(1)%pw 1928 tau_ia_r2(2)%pw => tau_ia_r(1)%pw 1929 END IF 1930 ELSE 1931 ALLOCATE (rho_ia_r2(nspins)) 1932 ALLOCATE (rho_ia_g2(nspins)) 1933 DO ispin = 1, nspins 1934 rho_ia_r2(ispin)%pw => rho_ia_r(ispin)%pw 1935 rho_ia_g2(ispin)%pw => rho_ia_g(ispin)%pw 1936 END DO 1937 1938 IF (ASSOCIATED(tau_ia_r)) THEN 1939 ALLOCATE (tau_ia_r2(nspins)) 1940 DO ispin = 1, nspins 1941 tau_ia_r2(ispin)%pw => tau_ia_r(ispin)%pw 1942 END DO 1943 END IF 1944 END IF 1945 1946 DO ispin = 1, nspins 1947 CALL pw_zero(work_v_xc(ispin)%pw) 1948 END DO 1949 1950 CALL xc_rho_set_update(rho_set=kernel_env%xc_rho1_set, rho_r=rho_ia_r2, rho_g=rho_ia_g2, tau=tau_ia_r2, & 1951 needs=kernel_env%xc_rho1_cflags, xc_deriv_method_id=kernel_env%deriv_method_id, & 1952 xc_rho_smooth_id=kernel_env%rho_smooth_id, pw_pool=auxbas_pw_pool) 1953 1954 CALL xc_calc_2nd_deriv(v_xc=work_v_xc, deriv_set=kernel_env%xc_deriv_set, rho_set=kernel_env%xc_rho_set, & 1955 rho1_set=kernel_env%xc_rho1_set, pw_pool=auxbas_pw_pool, & 1956 xc_section=kernel_env%xc_section, gapw=.FALSE., tddfpt_fac=kernel_env%beta) 1957 1958 DO ispin = 1, nspins 1959 ! pw2 = pw2 + alpha * pw1 1960 CALL pw_axpy(work_v_xc(ispin)%pw, A_ia_rspace(ispin)%pw, kernel_env%alpha) 1961 END DO 1962 1963 DEALLOCATE (rho_ia_g2, rho_ia_r2) 1964 1965 CALL timestop(handle) 1966 END SUBROUTINE tddfpt_apply_xc 1967 1968! ************************************************************************************************** 1969!> \brief Compute the ground-state charge density expressed in primary basis set. 1970!> \param rho_orb_struct ground-state density in primary basis set 1971!> \param is_rks_triplets indicates that the triplet excited states calculation using 1972!> spin-unpolarised molecular orbitals has been requested 1973!> \param qs_env Quickstep environment 1974!> \param sub_env parallel (sub)group environment 1975!> \param wfm_rho_orb work dense matrix with shape [nao x nao] distributed among 1976!> processors of the given parallel group (modified on exit) 1977!> \par History 1978!> * 06.2018 created by splitting the subroutine tddfpt_apply_admm_correction() in two 1979!> subroutines tddfpt_construct_ground_state_orb_density() and 1980!> tddfpt_construct_aux_fit_density [Sergey Chulkov] 1981! ************************************************************************************************** 1982 SUBROUTINE tddfpt_construct_ground_state_orb_density(rho_orb_struct, is_rks_triplets, qs_env, sub_env, wfm_rho_orb) 1983 TYPE(qs_rho_type), POINTER :: rho_orb_struct 1984 LOGICAL, INTENT(in) :: is_rks_triplets 1985 TYPE(qs_environment_type), POINTER :: qs_env 1986 TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env 1987 TYPE(cp_fm_type), POINTER :: wfm_rho_orb 1988 1989 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_construct_ground_state_orb_density', & 1990 routineP = moduleN//':'//routineN 1991 1992 INTEGER :: handle, ispin, nao, nspins 1993 INTEGER, DIMENSION(maxspins) :: nmo_occ 1994 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ij_ao 1995 1996 CALL timeset(routineN, handle) 1997 1998 nspins = SIZE(sub_env%mos_occ) 1999 DO ispin = 1, nspins 2000 CALL cp_fm_get_info(sub_env%mos_occ(ispin)%matrix, nrow_global=nao, ncol_global=nmo_occ(ispin)) 2001 END DO 2002 2003 CALL qs_rho_get(rho_orb_struct, rho_ao=rho_ij_ao) 2004 DO ispin = 1, nspins 2005 CALL cp_gemm('N', 'T', nao, nao, nmo_occ(ispin), 1.0_dp, & 2006 sub_env%mos_occ(ispin)%matrix, sub_env%mos_occ(ispin)%matrix, & 2007 0.0_dp, wfm_rho_orb) 2008 2009 CALL copy_fm_to_dbcsr(wfm_rho_orb, rho_ij_ao(ispin)%matrix, keep_sparsity=.TRUE.) 2010 END DO 2011 2012 ! take into account that all MOs are doubly occupied in spin-restricted case 2013 IF (nspins == 1 .AND. (.NOT. is_rks_triplets)) & 2014 CALL dbcsr_scale(rho_ij_ao(1)%matrix, 2.0_dp) 2015 2016 CALL qs_rho_update_rho(rho_orb_struct, qs_env, & 2017 pw_env_external=sub_env%pw_env, & 2018 task_list_external=sub_env%task_list_orb) 2019 2020 CALL timestop(handle) 2021 END SUBROUTINE tddfpt_construct_ground_state_orb_density 2022 2023! ************************************************************************************************** 2024!> \brief Project a charge density expressed in primary basis set into the auxiliary basis set. 2025!> \param rho_orb_struct response density in primary basis set 2026!> \param rho_aux_fit_struct response density in auxiliary basis set (modified on exit) 2027!> \param qs_env Quickstep environment 2028!> \param sub_env parallel (sub)group environment 2029!> \param wfm_rho_orb work dense matrix with shape [nao x nao] distributed among 2030!> processors of the given parallel group (modified on exit) 2031!> \param wfm_rho_aux_fit work dense matrix with shape [nao_aux x nao_aux] distributed among 2032!> processors of the given parallel group (modified on exit) 2033!> \param wfm_aux_orb work dense matrix with shape [nao_aux x nao] distributed among 2034!> processors of the given parallel group (modified on exit) 2035!> \par History 2036!> * 03.2017 the subroutine tddfpt_apply_admm_correction() was originally created by splitting 2037!> the subroutine tddfpt_apply_hfx() in two parts [Sergey Chulkov] 2038!> * 06.2018 created by splitting the subroutine tddfpt_apply_admm_correction() in two subroutines 2039!> tddfpt_construct_ground_state_orb_density() and tddfpt_construct_aux_fit_density() 2040!> in order to avoid code duplication [Sergey Chulkov] 2041! ************************************************************************************************** 2042 SUBROUTINE tddfpt_construct_aux_fit_density(rho_orb_struct, rho_aux_fit_struct, qs_env, sub_env, & 2043 wfm_rho_orb, wfm_rho_aux_fit, wfm_aux_orb) 2044 TYPE(qs_rho_type), POINTER :: rho_orb_struct, rho_aux_fit_struct 2045 TYPE(qs_environment_type), POINTER :: qs_env 2046 TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env 2047 TYPE(cp_fm_type), POINTER :: wfm_rho_orb, wfm_rho_aux_fit, wfm_aux_orb 2048 2049 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_construct_aux_fit_density', & 2050 routineP = moduleN//':'//routineN 2051 2052 INTEGER :: handle, ispin, nao, nao_aux, nspins 2053 REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_aux_fit_r 2054 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_aux_fit, rho_ao_orb 2055 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_aux_fit_g, rho_aux_fit_r 2056 TYPE(qs_ks_env_type), POINTER :: ks_env 2057 2058 CALL timeset(routineN, handle) 2059 2060 CPASSERT(ASSOCIATED(sub_env%admm_A)) 2061 2062 CALL get_qs_env(qs_env, ks_env=ks_env) 2063 CALL qs_rho_get(rho_orb_struct, rho_ao=rho_ao_orb) 2064 CALL qs_rho_get(rho_aux_fit_struct, rho_ao=rho_ao_aux_fit, rho_g=rho_aux_fit_g, & 2065 rho_r=rho_aux_fit_r, tot_rho_r=tot_rho_aux_fit_r) 2066 2067 nspins = SIZE(rho_ao_orb) 2068 2069 CALL cp_fm_get_info(sub_env%admm_A, nrow_global=nao_aux, ncol_global=nao) 2070 DO ispin = 1, nspins 2071 ! TO DO: consider sub_env%admm_A to be a DBCSR matrix 2072 CALL copy_dbcsr_to_fm(rho_ao_orb(ispin)%matrix, wfm_rho_orb) 2073 CALL cp_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, sub_env%admm_A, & 2074 wfm_rho_orb, 0.0_dp, wfm_aux_orb) 2075 CALL cp_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, sub_env%admm_A, wfm_aux_orb, & 2076 0.0_dp, wfm_rho_aux_fit) 2077 CALL copy_fm_to_dbcsr(wfm_rho_aux_fit, rho_ao_aux_fit(ispin)%matrix, keep_sparsity=.TRUE.) 2078 2079 CALL calculate_rho_elec(matrix_p=rho_ao_aux_fit(ispin)%matrix, & 2080 rho=rho_aux_fit_r(ispin), rho_gspace=rho_aux_fit_g(ispin), & 2081 total_rho=tot_rho_aux_fit_r(ispin), ks_env=ks_env, & 2082 soft_valid=.FALSE., basis_type="AUX_FIT", & 2083 pw_env_external=sub_env%pw_env, task_list_external=sub_env%task_list_aux_fit) 2084 END DO 2085 2086 CALL timestop(handle) 2087 END SUBROUTINE tddfpt_construct_aux_fit_density 2088 2089! ************************************************************************************************** 2090!> \brief Update action of TDDFPT operator on trial vectors by adding exact-exchange term. 2091!> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit) 2092!> \param evects trial vectors 2093!> \param gs_mos molecular orbitals optimised for the ground state (only occupied 2094!> molecular orbitals [component %mos_occ] are needed) 2095!> \param do_admm perform auxiliary density matrix method calculations 2096!> \param qs_env Quickstep environment 2097!> \param work_rho_ia_ao work sparse matrix with shape [nao x nao] distributed globally 2098!> to store response density (modified on exit) 2099!> \param work_hmat work sparse matrix with shape [nao x nao] distributed globally 2100!> (modified on exit) 2101!> \param wfm_rho_orb work dense matrix with shape [nao x nao] distributed globally 2102!> (modified on exit) 2103!> \par History 2104!> * 05.2016 compute all exact-exchange terms in one go [Sergey Chulkov] 2105!> * 03.2017 code related to ADMM correction is now moved to tddfpt_apply_admm_correction() 2106!> in order to compute this correction within parallel groups [Sergey Chulkov] 2107!> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by 2108!> Mohamed Fawzi on 10.2002. 2109! ************************************************************************************************** 2110 SUBROUTINE tddfpt_apply_hfx(Aop_evects, evects, gs_mos, do_admm, qs_env, & 2111 work_rho_ia_ao, work_hmat, wfm_rho_orb) 2112 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: Aop_evects, evects 2113 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 2114 INTENT(in) :: gs_mos 2115 LOGICAL, INTENT(in) :: do_admm 2116 TYPE(qs_environment_type), POINTER :: qs_env 2117 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: work_rho_ia_ao, work_hmat 2118 TYPE(cp_fm_type), POINTER :: wfm_rho_orb 2119 2120 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_hfx', & 2121 routineP = moduleN//':'//routineN 2122 2123 INTEGER :: handle, ispin, ivect, nao, nao_aux, & 2124 nspins, nvects 2125 INTEGER, DIMENSION(maxspins) :: nmo_occ 2126 REAL(kind=dp) :: alpha 2127 TYPE(admm_type), POINTER :: admm_env 2128 2129 CALL timeset(routineN, handle) 2130 2131 nspins = SIZE(evects, 1) 2132 nvects = SIZE(evects, 2) 2133 2134 IF (nspins > 1) THEN 2135 alpha = 2.0_dp 2136 ELSE 2137 alpha = 4.0_dp 2138 END IF 2139 2140 CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao) 2141 DO ispin = 1, nspins 2142 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ) 2143 END DO 2144 2145 IF (do_admm) THEN 2146 CALL get_qs_env(qs_env, admm_env=admm_env) 2147 CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux) 2148 END IF 2149 2150 ! some stuff from qs_ks_build_kohn_sham_matrix 2151 ! TO DO: add SIC support 2152 DO ivect = 1, nvects 2153 DO ispin = 1, nspins 2154 CALL cp_gemm('N', 'T', nao, nao, nmo_occ(ispin), 0.5_dp, gs_mos(ispin)%mos_occ, & 2155 evects(ispin, ivect)%matrix, 0.0_dp, wfm_rho_orb) 2156 CALL cp_gemm('N', 'T', nao, nao, nmo_occ(ispin), 0.5_dp, evects(ispin, ivect)%matrix, & 2157 gs_mos(ispin)%mos_occ, 1.0_dp, wfm_rho_orb) 2158 2159 CALL dbcsr_set(work_hmat(ispin)%matrix, 0.0_dp) 2160 IF (do_admm) THEN 2161 CALL cp_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, & 2162 wfm_rho_orb, 0.0_dp, admm_env%work_aux_orb) 2163 CALL cp_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%A, admm_env%work_aux_orb, & 2164 0.0_dp, admm_env%work_aux_aux) 2165 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_rho_ia_ao(ispin)%matrix, keep_sparsity=.TRUE.) 2166 ELSE 2167 CALL copy_fm_to_dbcsr(wfm_rho_orb, work_rho_ia_ao(ispin)%matrix, keep_sparsity=.TRUE.) 2168 END IF 2169 END DO 2170 2171 CALL tddft_hfx_matrix(work_hmat, work_rho_ia_ao, qs_env) 2172 2173 IF (do_admm) THEN 2174 DO ispin = 1, nspins 2175 CALL cp_dbcsr_sm_fm_multiply(work_hmat(ispin)%matrix, admm_env%A, admm_env%work_aux_orb, & 2176 ncol=nao, alpha=1.0_dp, beta=0.0_dp) 2177 2178 CALL cp_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, & 2179 admm_env%work_aux_orb, 0.0_dp, wfm_rho_orb) 2180 2181 CALL cp_gemm('N', 'N', nao, nmo_occ(ispin), nao, alpha, wfm_rho_orb, & 2182 gs_mos(ispin)%mos_occ, 1.0_dp, Aop_evects(ispin, ivect)%matrix) 2183 END DO 2184 ELSE 2185 DO ispin = 1, nspins 2186 CALL cp_dbcsr_sm_fm_multiply(work_hmat(ispin)%matrix, gs_mos(ispin)%mos_occ, & 2187 Aop_evects(ispin, ivect)%matrix, ncol=nmo_occ(ispin), & 2188 alpha=alpha, beta=1.0_dp) 2189 END DO 2190 END IF 2191 END DO 2192 2193 CALL timestop(handle) 2194 END SUBROUTINE tddfpt_apply_hfx 2195 2196! ************************************************************************************************** 2197!> \brief Compute action matrix-vector products. 2198!> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit) 2199!> \param evects TDDFPT trial vectors 2200!> \param S_evects cached matrix product S * evects where S is the overlap matrix 2201!> in primary basis set 2202!> \param gs_mos molecular orbitals optimised for the ground state 2203!> \param is_rks_triplets indicates that a triplet excited states calculation using 2204!> spin-unpolarised molecular orbitals has been requested 2205!> \param do_hfx flag that activates computation of exact-exchange terms 2206!> \param matrix_ks Kohn-Sham matrix 2207!> \param qs_env Quickstep environment 2208!> \param kernel_env kernel environment 2209!> \param kernel_env_admm_aux kernel environment for ADMM correction 2210!> \param sub_env parallel (sub)group environment 2211!> \param work_matrices collection of work matrices (modified on exit) 2212!> \par History 2213!> * 06.2016 created [Sergey Chulkov] 2214!> * 03.2017 refactored [Sergey Chulkov] 2215! ************************************************************************************************** 2216 SUBROUTINE tddfpt_compute_Aop_evects(Aop_evects, evects, S_evects, gs_mos, is_rks_triplets, & 2217 do_hfx, matrix_ks, qs_env, kernel_env, kernel_env_admm_aux, & 2218 sub_env, work_matrices) 2219 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: Aop_evects, evects, S_evects 2220 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 2221 INTENT(in) :: gs_mos 2222 LOGICAL, INTENT(in) :: is_rks_triplets, do_hfx 2223 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_ks 2224 TYPE(qs_environment_type), POINTER :: qs_env 2225 TYPE(tddfpt_kernel_env_type), INTENT(in) :: kernel_env, kernel_env_admm_aux 2226 TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env 2227 TYPE(tddfpt_work_matrices), INTENT(inout) :: work_matrices 2228 2229 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_compute_Aop_evects', & 2230 routineP = moduleN//':'//routineN 2231 2232 INTEGER :: handle, ispin, ivect, nao, nspins, nvects 2233 INTEGER, DIMENSION(maxspins) :: nmo_occ 2234 LOGICAL :: do_admm 2235 TYPE(cp_para_env_type), POINTER :: para_env 2236 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ia_ao, rho_ia_ao_aux_fit 2237 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_ia_g, rho_ia_g_aux_fit, rho_ia_r, & 2238 rho_ia_r_aux_fit, tau_ia_r, & 2239 tau_ia_r_aux_fit 2240 2241 CALL timeset(routineN, handle) 2242 2243 nspins = SIZE(evects, 1) 2244 nvects = SIZE(evects, 2) 2245 do_admm = ASSOCIATED(sub_env%admm_A) 2246 2247 IF (debug_this_module) THEN 2248 CPASSERT(nspins > 0) 2249 CPASSERT(SIZE(Aop_evects, 1) == nspins) 2250 CPASSERT(SIZE(Aop_evects, 2) == nvects) 2251 CPASSERT(SIZE(S_evects, 1) == nspins) 2252 CPASSERT(SIZE(S_evects, 2) == nvects) 2253 CPASSERT(SIZE(gs_mos) == nspins) 2254 END IF 2255 2256 CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao) 2257 DO ispin = 1, nspins 2258 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ) 2259 END DO 2260 2261 IF (nvects > 0) THEN 2262 CALL cp_fm_get_info(evects(1, 1)%matrix, para_env=para_env) 2263 CALL qs_rho_get(work_matrices%rho_orb_struct_sub, rho_ao=rho_ia_ao, rho_g=rho_ia_g, rho_r=rho_ia_r, tau_r=tau_ia_r) 2264 IF (do_hfx .AND. do_admm) THEN 2265 CALL qs_rho_get(work_matrices%rho_aux_fit_struct_sub, rho_ao=rho_ia_ao_aux_fit, rho_g=rho_ia_g_aux_fit, & 2266 rho_r=rho_ia_r_aux_fit, tau_r=tau_ia_r_aux_fit) 2267 END IF 2268 2269 IF (ALLOCATED(work_matrices%evects_sub)) THEN 2270 DO ivect = 1, nvects 2271 DO ispin = 1, nspins 2272 CALL cp_fm_copy_general(evects(ispin, ivect)%matrix, work_matrices%evects_sub(ispin, ivect)%matrix, para_env) 2273 END DO 2274 END DO 2275 END IF 2276 2277 DO ivect = 1, nvects 2278 IF (ALLOCATED(work_matrices%evects_sub)) THEN 2279 IF (ASSOCIATED(work_matrices%evects_sub(1, ivect)%matrix)) THEN 2280 DO ispin = 1, nspins 2281 CALL cp_gemm('N', 'T', nao, nao, nmo_occ(ispin), & 2282 0.5_dp, sub_env%mos_occ(ispin)%matrix, & 2283 work_matrices%evects_sub(ispin, ivect)%matrix, & 2284 0.0_dp, work_matrices%rho_ao_orb_fm_sub) 2285 CALL cp_gemm('N', 'T', nao, nao, nmo_occ(ispin), & 2286 0.5_dp, work_matrices%evects_sub(ispin, ivect)%matrix, & 2287 sub_env%mos_occ(ispin)%matrix, & 2288 1.0_dp, work_matrices%rho_ao_orb_fm_sub) 2289 2290 CALL copy_fm_to_dbcsr(work_matrices%rho_ao_orb_fm_sub, & 2291 rho_ia_ao(ispin)%matrix, keep_sparsity=.TRUE.) 2292 END DO 2293 ELSE 2294 ! skip trial vectors which are assigned to different parallel groups 2295 CYCLE 2296 END IF 2297 ELSE 2298 DO ispin = 1, nspins 2299 CALL cp_gemm('N', 'T', nao, nao, nmo_occ(ispin), 0.5_dp, sub_env%mos_occ(ispin)%matrix, & 2300 evects(ispin, ivect)%matrix, 0.0_dp, work_matrices%rho_ao_orb_fm_sub) 2301 CALL cp_gemm('N', 'T', nao, nao, nmo_occ(ispin), 0.5_dp, evects(ispin, ivect)%matrix, & 2302 sub_env%mos_occ(ispin)%matrix, 1.0_dp, work_matrices%rho_ao_orb_fm_sub) 2303 2304 CALL copy_fm_to_dbcsr(work_matrices%rho_ao_orb_fm_sub, & 2305 rho_ia_ao(ispin)%matrix, keep_sparsity=.TRUE.) 2306 END DO 2307 END IF 2308 2309 CALL qs_rho_update_rho(work_matrices%rho_orb_struct_sub, qs_env, & 2310 pw_env_external=sub_env%pw_env, task_list_external=sub_env%task_list_orb) 2311 2312 DO ispin = 1, nspins 2313 CALL dbcsr_set(work_matrices%A_ia_munu_sub(ispin)%matrix, 0.0_dp) 2314 END DO 2315 2316 ! electron-hole exchange-correlation interaction 2317 DO ispin = 1, nspins 2318 CALL pw_zero(work_matrices%A_ia_rspace_sub(ispin)%pw) 2319 END DO 2320 2321 ! C_x d^{2}E_{x}^{DFT}[\rho] / d\rho^2 2322 ! + C_{HF} d^{2}E_{x, ADMM}^{DFT}[\rho] / d\rho^2 in case of ADMM calculation 2323 CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, & 2324 rho_ia_struct=work_matrices%rho_orb_struct_sub, & 2325 is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, & 2326 work_v_xc=work_matrices%wpw_rspace_sub) 2327 2328 ! ADMM correction 2329 IF (do_admm) THEN 2330 CALL tddfpt_construct_aux_fit_density(rho_orb_struct=work_matrices%rho_orb_struct_sub, & 2331 rho_aux_fit_struct=work_matrices%rho_aux_fit_struct_sub, & 2332 qs_env=qs_env, sub_env=sub_env, & 2333 wfm_rho_orb=work_matrices%rho_ao_orb_fm_sub, & 2334 wfm_rho_aux_fit=work_matrices%rho_ao_aux_fit_fm_sub, & 2335 wfm_aux_orb=work_matrices%wfm_aux_orb_sub) 2336 2337 ! - C_{HF} d^{2}E_{x, ADMM}^{DFT}[\hat{\rho}] / d\hat{\rho}^2 2338 CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, & 2339 kernel_env=kernel_env_admm_aux, & 2340 rho_ia_struct=work_matrices%rho_aux_fit_struct_sub, & 2341 is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, & 2342 work_v_xc=work_matrices%wpw_rspace_sub) 2343 END IF 2344 2345 ! electron-hole Coulomb interaction 2346 IF (.NOT. is_rks_triplets) THEN 2347 ! a sum J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu can be computed by solving 2348 ! the Poisson equation for combined density (rho_{ia,alpha} + rho_{ia,beta}) . 2349 ! The following action will destroy reciprocal-space grid in spin-unrestricted case. 2350 DO ispin = 2, nspins 2351 CALL pw_axpy(rho_ia_g(ispin)%pw, rho_ia_g(1)%pw) 2352 END DO 2353 2354 CALL tddfpt_apply_coulomb(A_ia_rspace=work_matrices%A_ia_rspace_sub, & 2355 rho_ia_g=rho_ia_g(1)%pw, pw_env=sub_env%pw_env, & 2356 work_v_gspace=work_matrices%wpw_gspace_sub(1), & 2357 work_v_rspace=work_matrices%wpw_rspace_sub(1)) 2358 END IF 2359 2360 ! convert from the plane-wave representation into the Gaussian basis set representation 2361 DO ispin = 1, nspins 2362 CALL pw_scale(work_matrices%A_ia_rspace_sub(ispin)%pw, work_matrices%A_ia_rspace_sub(ispin)%pw%pw_grid%dvol) 2363 2364 CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), & 2365 hmat=work_matrices%A_ia_munu_sub(ispin), & 2366 qs_env=qs_env, calculate_forces=.FALSE., gapw=.FALSE., & 2367 pw_env_external=sub_env%pw_env, task_list_external=sub_env%task_list_orb) 2368 END DO 2369 2370 IF (ALLOCATED(work_matrices%evects_sub)) THEN 2371 DO ispin = 1, nspins 2372 CALL cp_dbcsr_sm_fm_multiply(work_matrices%A_ia_munu_sub(ispin)%matrix, sub_env%mos_occ(ispin)%matrix, & 2373 work_matrices%Aop_evects_sub(ispin, ivect)%matrix, & 2374 ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp) 2375 END DO 2376 ELSE 2377 DO ispin = 1, nspins 2378 CALL cp_dbcsr_sm_fm_multiply(work_matrices%A_ia_munu_sub(ispin)%matrix, sub_env%mos_occ(ispin)%matrix, & 2379 Aop_evects(ispin, ivect)%matrix, ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp) 2380 END DO 2381 END IF 2382 END DO 2383 2384 IF (ALLOCATED(work_matrices%evects_sub)) THEN 2385 DO ivect = 1, nvects 2386 DO ispin = 1, nspins 2387 CALL cp_fm_copy_general(work_matrices%Aop_evects_sub(ispin, ivect)%matrix, & 2388 Aop_evects(ispin, ivect)%matrix, para_env) 2389 END DO 2390 END DO 2391 END IF 2392 2393 ! orbital energy difference term 2394 CALL tddfpt_apply_energy_diff(Aop_evects=Aop_evects, evects=evects, S_evects=S_evects, gs_mos=gs_mos, & 2395 matrix_ks=matrix_ks, work_fm_ao_mo_occ=work_matrices%wfm_ao_mo_occ) 2396 2397 IF (do_hfx) THEN 2398 CALL tddfpt_apply_hfx(Aop_evects=Aop_evects, evects=evects, gs_mos=gs_mos, do_admm=do_admm, & 2399 qs_env=qs_env, work_rho_ia_ao=work_matrices%hfx_rho_ao, & 2400 work_hmat=work_matrices%hfx_hmat, wfm_rho_orb=work_matrices%hfx_fm_ao_ao) 2401 END IF 2402 END IF 2403 2404 CALL timestop(handle) 2405 END SUBROUTINE tddfpt_compute_Aop_evects 2406 2407! ************************************************************************************************** 2408!> \brief Solve eigenproblem for the reduced action matrix and find new Ritz eigenvectors and 2409!> eigenvalues. 2410!> \param ritz_vects Ritz eigenvectors (initialised on exit) 2411!> \param Aop_ritz approximate action of TDDFPT operator on Ritz vectors (initialised on exit) 2412!> \param evals Ritz eigenvalues (initialised on exit) 2413!> \param krylov_vects Krylov's vectors 2414!> \param Aop_krylov action of TDDFPT operator on Krylov's vectors 2415!> \par History 2416!> * 06.2016 created [Sergey Chulkov] 2417!> * 03.2017 altered prototype, OpenMP parallelisation [Sergey Chulkov] 2418! ************************************************************************************************** 2419 SUBROUTINE tddfpt_compute_ritz_vects(ritz_vects, Aop_ritz, evals, krylov_vects, Aop_krylov) 2420 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: ritz_vects, Aop_ritz 2421 REAL(kind=dp), DIMENSION(:), INTENT(out) :: evals 2422 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: krylov_vects, Aop_krylov 2423 2424 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_compute_ritz_vects', & 2425 routineP = moduleN//':'//routineN 2426 2427 INTEGER :: handle, ikv, irv, ispin, nkvs, nrvs, & 2428 nspins 2429 REAL(kind=dp) :: act 2430 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: Atilde, evects_Atilde 2431 TYPE(cp_blacs_env_type), POINTER :: blacs_env_global 2432 TYPE(cp_fm_struct_type), POINTER :: fm_struct 2433 TYPE(cp_fm_type), POINTER :: Atilde_fm, evects_Atilde_fm 2434 2435 CALL timeset(routineN, handle) 2436 2437 nspins = SIZE(krylov_vects, 1) 2438 nkvs = SIZE(krylov_vects, 2) 2439 nrvs = SIZE(ritz_vects, 2) 2440 2441 CALL cp_fm_get_info(krylov_vects(1, 1)%matrix, context=blacs_env_global) 2442 2443 CALL cp_fm_struct_create(fm_struct, nrow_global=nkvs, ncol_global=nkvs, context=blacs_env_global) 2444 CALL cp_fm_create(Atilde_fm, fm_struct) 2445 CALL cp_fm_create(evects_Atilde_fm, fm_struct) 2446 CALL cp_fm_struct_release(fm_struct) 2447 2448 ! *** compute upper-diagonal reduced action matrix *** 2449 ALLOCATE (Atilde(nkvs, nkvs)) 2450 ! TO DO: the subroutine 'cp_fm_contracted_trace' will compute all elements of 2451 ! the matrix 'Atilde', however only upper-triangular elements are actually needed 2452 CALL cp_fm_contracted_trace(Aop_krylov, krylov_vects, Atilde) 2453 CALL cp_fm_set_submatrix(Atilde_fm, Atilde) 2454 DEALLOCATE (Atilde) 2455 2456 ! *** solve an eigenproblem for the reduced matrix *** 2457 CALL choose_eigv_solver(Atilde_fm, evects_Atilde_fm, evals(1:nkvs)) 2458 2459 ALLOCATE (evects_Atilde(nkvs, nrvs)) 2460 CALL cp_fm_get_submatrix(evects_Atilde_fm, evects_Atilde, start_row=1, start_col=1, n_rows=nkvs, n_cols=nrvs) 2461 CALL cp_fm_release(evects_Atilde_fm) 2462 CALL cp_fm_release(Atilde_fm) 2463 2464!$OMP PARALLEL DO DEFAULT(NONE), & 2465!$OMP PRIVATE(act, ikv, irv, ispin), & 2466!$OMP SHARED(Aop_krylov, Aop_ritz, krylov_vects, evects_Atilde, nkvs, nrvs, nspins, ritz_vects) 2467 DO irv = 1, nrvs 2468 DO ispin = 1, nspins 2469 CALL cp_fm_set_all(ritz_vects(ispin, irv)%matrix, 0.0_dp) 2470 CALL cp_fm_set_all(Aop_ritz(ispin, irv)%matrix, 0.0_dp) 2471 END DO 2472 2473 DO ikv = 1, nkvs 2474 act = evects_Atilde(ikv, irv) 2475 DO ispin = 1, nspins 2476 CALL cp_fm_scale_and_add(1.0_dp, ritz_vects(ispin, irv)%matrix, & 2477 act, krylov_vects(ispin, ikv)%matrix) 2478 CALL cp_fm_scale_and_add(1.0_dp, Aop_ritz(ispin, irv)%matrix, & 2479 act, Aop_krylov(ispin, ikv)%matrix) 2480 END DO 2481 END DO 2482 END DO 2483!$OMP END PARALLEL DO 2484 2485 DEALLOCATE (evects_Atilde) 2486 2487 CALL timestop(handle) 2488 END SUBROUTINE tddfpt_compute_ritz_vects 2489 2490! ************************************************************************************************** 2491!> \brief Expand Krylov space by computing residual vectors. 2492!> \param residual_vects residual vectors (modified on exit) 2493!> \param evals Ritz eigenvalues (modified on exit) 2494!> \param ritz_vects Ritz eigenvectors 2495!> \param Aop_ritz approximate action of TDDFPT operator on Ritz vectors 2496!> \param gs_mos molecular orbitals optimised for the ground state 2497!> \param matrix_s overlap matrix 2498!> \param work_fm_ao_mo_occ work dense matrices with shape [nao x nmo_occ(spin)] 2499!> (modified on exit) 2500!> \param work_fm_mo_virt_mo_occ work dense matrices with shape [nmo_virt(spin) x nmo_occ(spin)] 2501!> (modified on exit) 2502!> \par History 2503!> * 06.2016 created [Sergey Chulkov] 2504!> * 03.2017 refactored to achieve significant performance gain [Sergey Chulkov] 2505! ************************************************************************************************** 2506 SUBROUTINE tddfpt_compute_residual_vects(residual_vects, evals, ritz_vects, Aop_ritz, gs_mos, & 2507 matrix_s, work_fm_ao_mo_occ, work_fm_mo_virt_mo_occ) 2508 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: residual_vects 2509 REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals 2510 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: ritz_vects, Aop_ritz 2511 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 2512 INTENT(in) :: gs_mos 2513 TYPE(dbcsr_type), POINTER :: matrix_s 2514 TYPE(cp_fm_p_type), DIMENSION(:), INTENT(in) :: work_fm_ao_mo_occ, work_fm_mo_virt_mo_occ 2515 2516 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_compute_residual_vects', & 2517 routineP = moduleN//':'//routineN 2518 REAL(kind=dp), PARAMETER :: eref_scale = 0.99_dp, threshold = 16.0_dp*EPSILON(1.0_dp) 2519 2520 INTEGER :: handle, icol_local, irow_local, irv, & 2521 ispin, nao, ncols_local, nrows_local, & 2522 nrvs, nspins 2523 INTEGER, DIMENSION(:), POINTER :: col_indices_local, row_indices_local 2524 INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_virt 2525 REAL(kind=dp) :: e_occ_plus_lambda, eref, lambda 2526 REAL(kind=dp), DIMENSION(:, :), POINTER :: weights_ldata 2527 2528 CALL timeset(routineN, handle) 2529 2530 nspins = SIZE(residual_vects, 1) 2531 nrvs = SIZE(residual_vects, 2) 2532 2533 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao) 2534 DO ispin = 1, nspins 2535 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ) 2536 nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt) 2537 END DO 2538 2539 IF (nrvs > 0) THEN 2540 ! *** actually compute residual vectors *** 2541 DO irv = 1, nrvs 2542 lambda = evals(irv) 2543 2544 DO ispin = 1, nspins 2545 CALL cp_fm_get_info(work_fm_mo_virt_mo_occ(ispin)%matrix, nrow_local=nrows_local, ncol_local=ncols_local, & 2546 row_indices=row_indices_local, col_indices=col_indices_local, local_data=weights_ldata) 2547 2548 ! work_fm_ao_mo_occ := Ab(ispin, irv) - evals(irv) b(ispin, irv), where 'b' is a Ritz vector 2549 CALL cp_dbcsr_sm_fm_multiply(matrix_s, ritz_vects(ispin, irv)%matrix, work_fm_ao_mo_occ(ispin)%matrix, & 2550 ncol=nmo_occ(ispin), alpha=-lambda, beta=0.0_dp) 2551 CALL cp_fm_scale_and_add(1.0_dp, work_fm_ao_mo_occ(ispin)%matrix, 1.0_dp, Aop_ritz(ispin, irv)%matrix) 2552 2553 CALL cp_gemm('T', 'N', nmo_virt(ispin), nmo_occ(ispin), nao, 1.0_dp, gs_mos(ispin)%mos_virt, & 2554 work_fm_ao_mo_occ(ispin)%matrix, 0.0_dp, work_fm_mo_virt_mo_occ(ispin)%matrix) 2555 2556 DO icol_local = 1, ncols_local 2557 e_occ_plus_lambda = gs_mos(ispin)%evals_occ(col_indices_local(icol_local)) + lambda 2558 2559 DO irow_local = 1, nrows_local 2560 eref = gs_mos(ispin)%evals_virt(row_indices_local(irow_local)) - e_occ_plus_lambda 2561 2562 ! eref = e_virt - e_occ - lambda = e_virt - e_occ - (eref_scale*lambda + (1-eref_scale)*lambda); 2563 ! eref_new = e_virt - e_occ - eref_scale*lambda = eref + (1 - eref_scale)*lambda 2564 IF (ABS(eref) < threshold) & 2565 eref = eref + (1.0_dp - eref_scale)*lambda 2566 2567 weights_ldata(irow_local, icol_local) = weights_ldata(irow_local, icol_local)/eref 2568 END DO 2569 END DO 2570 2571 CALL cp_gemm('N', 'N', nao, nmo_occ(ispin), nmo_virt(ispin), 1.0_dp, gs_mos(ispin)%mos_virt, & 2572 work_fm_mo_virt_mo_occ(ispin)%matrix, 0.0_dp, residual_vects(ispin, irv)%matrix) 2573 END DO 2574 END DO 2575 END IF 2576 2577 CALL timestop(handle) 2578 END SUBROUTINE tddfpt_compute_residual_vects 2579 2580! ************************************************************************************************** 2581!> \brief Perform Davidson iterations. 2582!> \param evects TDDFPT trial vectors (modified on exit) 2583!> \param evals TDDFPT eigenvalues (modified on exit) 2584!> \param S_evects cached matrix product S * evects (modified on exit) 2585!> \param gs_mos molecular orbitals optimised for the ground state 2586!> \param do_hfx flag that activates computation of exact-exchange terms 2587!> \param tddfpt_control TDDFPT control parameters 2588!> \param qs_env Quickstep environment 2589!> \param kernel_env kernel environment 2590!> \param kernel_env_admm_aux kernel environment for ADMM correction 2591!> \param sub_env parallel (sub)group environment 2592!> \param logger CP2K logger 2593!> \param iter_unit I/O unit to write basic iteration information 2594!> \param energy_unit I/O unit to write detailed energy information 2595!> \param tddfpt_print_section TDDFPT print input section (need to write TDDFPT restart files) 2596!> \param work_matrices collection of work matrices (modified on exit) 2597!> \return energy convergence achieved (in Hartree) 2598!> \par History 2599!> * 03.2017 code related to Davidson eigensolver has been moved here from the main subroutine 2600!> tddfpt() [Sergey Chulkov] 2601!> \note Based on the subroutines apply_op() and iterative_solver() originally created by 2602!> Thomas Chassaing in 2002. 2603! ************************************************************************************************** 2604 FUNCTION tddfpt_davidson_solver(evects, evals, S_evects, gs_mos, do_hfx, tddfpt_control, & 2605 qs_env, kernel_env, kernel_env_admm_aux, & 2606 sub_env, logger, iter_unit, energy_unit, & 2607 tddfpt_print_section, work_matrices) RESULT(conv) 2608 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(inout) :: evects 2609 REAL(kind=dp), DIMENSION(:), INTENT(inout) :: evals 2610 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(inout) :: S_evects 2611 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 2612 INTENT(in) :: gs_mos 2613 LOGICAL, INTENT(in) :: do_hfx 2614 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control 2615 TYPE(qs_environment_type), POINTER :: qs_env 2616 TYPE(tddfpt_kernel_env_type), INTENT(in) :: kernel_env, kernel_env_admm_aux 2617 TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env 2618 TYPE(cp_logger_type), POINTER :: logger 2619 INTEGER, INTENT(in) :: iter_unit, energy_unit 2620 TYPE(section_vals_type), POINTER :: tddfpt_print_section 2621 TYPE(tddfpt_work_matrices), INTENT(inout) :: work_matrices 2622 REAL(kind=dp) :: conv 2623 2624 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_davidson_solver', & 2625 routineP = moduleN//':'//routineN 2626 2627 INTEGER :: handle, ispin, istate, iter, & 2628 max_krylov_vects, nspins, nstates, & 2629 nstates_conv, nvects_exist, nvects_new 2630 INTEGER(kind=int_8) :: nstates_total 2631 LOGICAL :: is_nonortho 2632 REAL(kind=dp) :: t1, t2 2633 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals_last 2634 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:, :) :: Aop_krylov, Aop_ritz, krylov_vects, & 2635 S_krylov 2636 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s 2637 2638 CALL timeset(routineN, handle) 2639 2640 nspins = SIZE(gs_mos) 2641 nstates = tddfpt_control%nstates 2642 nstates_total = tddfpt_total_number_of_states(gs_mos) 2643 2644 IF (debug_this_module) THEN 2645 CPASSERT(SIZE(evects, 1) == nspins) 2646 CPASSERT(SIZE(evects, 2) == nstates) 2647 CPASSERT(SIZE(evals) == nstates) 2648 END IF 2649 2650 CALL get_qs_env(qs_env, matrix_ks=matrix_ks, matrix_s=matrix_s) 2651 2652 ! adjust the number of Krylov vectors 2653 max_krylov_vects = tddfpt_control%nkvs 2654 IF (max_krylov_vects < nstates) max_krylov_vects = nstates 2655 IF (INT(max_krylov_vects, kind=int_8) > nstates_total) max_krylov_vects = INT(nstates_total) 2656 2657 ALLOCATE (Aop_ritz(nspins, nstates)) 2658 DO istate = 1, nstates 2659 DO ispin = 1, nspins 2660 NULLIFY (Aop_ritz(ispin, istate)%matrix) 2661 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, Aop_ritz(ispin, istate)%matrix) 2662 END DO 2663 END DO 2664 2665 ALLOCATE (evals_last(max_krylov_vects)) 2666 ALLOCATE (Aop_krylov(nspins, max_krylov_vects), krylov_vects(nspins, max_krylov_vects), S_krylov(nspins, max_krylov_vects)) 2667 DO istate = 1, max_krylov_vects 2668 DO ispin = 1, nspins 2669 NULLIFY (Aop_krylov(ispin, istate)%matrix, krylov_vects(ispin, istate)%matrix, S_krylov(ispin, istate)%matrix) 2670 END DO 2671 END DO 2672 2673 DO istate = 1, nstates 2674 DO ispin = 1, nspins 2675 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, krylov_vects(ispin, istate)%matrix) 2676 CALL cp_fm_to_fm(evects(ispin, istate)%matrix, krylov_vects(ispin, istate)%matrix) 2677 2678 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, S_krylov(ispin, istate)%matrix) 2679 CALL cp_fm_to_fm(S_evects(ispin, istate)%matrix, S_krylov(ispin, istate)%matrix) 2680 2681 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, Aop_krylov(ispin, istate)%matrix) 2682 END DO 2683 END DO 2684 2685 nvects_exist = 0 2686 nvects_new = nstates 2687 2688 t1 = m_walltime() 2689 2690 DO 2691 ! davidson iteration 2692 CALL cp_iterate(logger%iter_info, iter_nr_out=iter) 2693 2694 CALL tddfpt_compute_Aop_evects(Aop_evects=Aop_krylov(:, nvects_exist + 1:nvects_exist + nvects_new), & 2695 evects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), & 2696 S_evects=S_krylov(:, nvects_exist + 1:nvects_exist + nvects_new), & 2697 gs_mos=gs_mos, is_rks_triplets=tddfpt_control%rks_triplets, & 2698 do_hfx=do_hfx, matrix_ks=matrix_ks, & 2699 qs_env=qs_env, kernel_env=kernel_env, & 2700 kernel_env_admm_aux=kernel_env_admm_aux, & 2701 sub_env=sub_env, & 2702 work_matrices=work_matrices) 2703 2704 CALL tddfpt_compute_ritz_vects(ritz_vects=evects, Aop_ritz=Aop_ritz, & 2705 evals=evals_last(1:nvects_exist + nvects_new), & 2706 krylov_vects=krylov_vects(:, 1:nvects_exist + nvects_new), & 2707 Aop_krylov=Aop_krylov(:, 1:nvects_exist + nvects_new)) 2708 2709 CALL tddfpt_write_restart(evects=evects, evals=evals_last(1:nstates), gs_mos=gs_mos, & 2710 logger=logger, tddfpt_print_section=tddfpt_print_section) 2711 2712 conv = MAXVAL(ABS(evals_last(1:nstates) - evals(1:nstates))) 2713 2714 nvects_exist = nvects_exist + nvects_new 2715 IF (nvects_exist + nvects_new > max_krylov_vects) & 2716 nvects_new = max_krylov_vects - nvects_exist 2717 IF (iter >= tddfpt_control%niters) nvects_new = 0 2718 2719 IF (conv > tddfpt_control%conv .AND. nvects_new > 0) THEN 2720 ! compute residual vectors for the next iteration 2721 DO istate = 1, nvects_new 2722 DO ispin = 1, nspins 2723 NULLIFY (Aop_krylov(ispin, nvects_exist + istate)%matrix, & 2724 krylov_vects(ispin, nvects_exist + istate)%matrix, & 2725 S_krylov(ispin, nvects_exist + istate)%matrix) 2726 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, & 2727 krylov_vects(ispin, nvects_exist + istate)%matrix) 2728 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, & 2729 S_krylov(ispin, nvects_exist + istate)%matrix) 2730 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, & 2731 Aop_krylov(ispin, nvects_exist + istate)%matrix) 2732 END DO 2733 END DO 2734 2735 CALL tddfpt_compute_residual_vects(residual_vects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), & 2736 evals=evals_last(1:nvects_new), & 2737 ritz_vects=evects(:, 1:nvects_new), Aop_ritz=Aop_ritz(:, 1:nvects_new), & 2738 gs_mos=gs_mos, matrix_s=matrix_s(1)%matrix, & 2739 work_fm_ao_mo_occ=work_matrices%wfm_ao_mo_occ, & 2740 work_fm_mo_virt_mo_occ=work_matrices%wfm_mo_virt_mo_occ) 2741 2742 CALL tddfpt_orthogonalize_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), & 2743 work_matrices%S_C0_C0T, work_matrices%wfm_ao_mo_occ) 2744 2745 CALL tddfpt_orthonormalize_psi1_psi1(krylov_vects(:, 1:nvects_exist + nvects_new), nvects_new, & 2746 S_krylov(:, 1:nvects_exist + nvects_new), matrix_s(1)%matrix) 2747 2748 is_nonortho = tddfpt_is_nonorthogonal_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), & 2749 work_matrices%S_C0, tddfpt_control%orthogonal_eps, & 2750 work_matrices%wfm_mo_occ_mo_occ) 2751 ELSE 2752 ! convergence or the maximum number of Krylov vectors have been achieved 2753 nvects_new = 0 2754 is_nonortho = .FALSE. 2755 END IF 2756 2757 t2 = m_walltime() 2758 IF (energy_unit > 0) THEN 2759 WRITE (energy_unit, '(/,4X,A,T14,A,T36,A)') "State", "Exc. energy (eV)", "Convergence (eV)" 2760 DO istate = 1, nstates 2761 WRITE (energy_unit, '(1X,I8,T12,F14.7,T38,ES11.4)') istate, & 2762 evals_last(istate)*evolt, (evals_last(istate) - evals(istate))*evolt 2763 END DO 2764 WRITE (energy_unit, *) 2765 CALL m_flush(energy_unit) 2766 END IF 2767 2768 IF (iter_unit > 0) THEN 2769 nstates_conv = 0 2770 DO istate = 1, nstates 2771 IF (ABS(evals_last(istate) - evals(istate)) <= tddfpt_control%conv) & 2772 nstates_conv = nstates_conv + 1 2773 END DO 2774 2775 WRITE (iter_unit, '(1X,I8,T12,F7.1,T24,ES11.4,T42,I8)') iter, t2 - t1, conv, nstates_conv 2776 CALL m_flush(iter_unit) 2777 END IF 2778 2779 t1 = t2 2780 evals(1:nstates) = evals_last(1:nstates) 2781 2782 ! nvects_new == 0 if iter >= tddfpt_control%niters 2783 IF (nvects_new == 0 .OR. is_nonortho) THEN 2784 ! restart Davidson iterations 2785 CALL tddfpt_orthogonalize_psi1_psi0(evects, work_matrices%S_C0_C0T, work_matrices%wfm_ao_mo_occ) 2786 CALL tddfpt_orthonormalize_psi1_psi1(evects, nstates, S_evects, matrix_s(1)%matrix) 2787 2788 EXIT 2789 END IF 2790 END DO 2791 2792 DO istate = nvects_exist + nvects_new, 1, -1 2793 DO ispin = nspins, 1, -1 2794 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, Aop_krylov(ispin, istate)%matrix) 2795 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, S_krylov(ispin, istate)%matrix) 2796 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, krylov_vects(ispin, istate)%matrix) 2797 END DO 2798 END DO 2799 DEALLOCATE (Aop_krylov, krylov_vects, S_krylov) 2800 DEALLOCATE (evals_last) 2801 2802 DO istate = nstates, 1, -1 2803 DO ispin = nspins, 1, -1 2804 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, Aop_ritz(ispin, istate)%matrix) 2805 END DO 2806 END DO 2807 DEALLOCATE (Aop_ritz) 2808 2809 CALL timestop(handle) 2810 END FUNCTION tddfpt_davidson_solver 2811 2812! ************************************************************************************************** 2813!> \brief Compute the action of the dipole operator on the ground state wave function. 2814!> \param dipole_op_mos_occ 2-D array [x,y,z ; spin] of matrices where to put the computed quantity 2815!> (allocated and initialised on exit) 2816!> \param tddfpt_control TDDFPT control parameters 2817!> \param gs_mos molecular orbitals optimised for the ground state 2818!> \param qs_env Quickstep environment 2819!> \par History 2820!> * 05.2016 created as 'tddfpt_print_summary' [Sergey Chulkov] 2821!> * 06.2018 dipole operator based on the Berry-phase formula [Sergey Chulkov] 2822!> * 08.2018 splited of from 'tddfpt_print_summary' and merged with code from 'tddfpt' 2823!> [Sergey Chulkov] 2824!> \note \parblock 2825!> Adapted version of the subroutine find_contributions() which was originally created 2826!> by Thomas Chassaing on 02.2005. 2827!> 2828!> The relation between dipole integrals in velocity and length forms are the following: 2829!> \f[<\psi_i|\nabla|\psi_a> = <\psi_i|\vec{r}|\hat{H}\psi_a> - <\hat{H}\psi_i|\vec{r}|\psi_a> 2830!> = (\epsilon_a - \epsilon_i) <\psi_i|\vec{r}|\psi_a> .\f], 2831!> due to the commutation identity: 2832!> \f[\vec{r}\hat{H} - \hat{H}\vec{r} = [\vec{r},\hat{H}] = [\vec{r},-1/2 \nabla^2] = \nabla\f] . 2833!> \endparblock 2834! ************************************************************************************************** 2835 SUBROUTINE tddfpt_dipole_operator(dipole_op_mos_occ, tddfpt_control, gs_mos, qs_env) 2836 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:, :), & 2837 INTENT(inout) :: dipole_op_mos_occ 2838 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control 2839 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 2840 INTENT(in) :: gs_mos 2841 TYPE(qs_environment_type), POINTER :: qs_env 2842 2843 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_dipole_operator', & 2844 routineP = moduleN//':'//routineN 2845 2846 INTEGER :: handle, i_cos_sin, icol, ideriv, irow, & 2847 ispin, jderiv, nao, ncols_local, & 2848 ndim_periodic, nrows_local, nspins 2849 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 2850 INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_virt 2851 REAL(kind=dp) :: eval_occ 2852 REAL(kind=dp), DIMENSION(3) :: kvec, reference_point 2853 REAL(kind=dp), DIMENSION(:, :), POINTER :: local_data_ediff, local_data_wfm 2854 TYPE(cell_type), POINTER :: cell 2855 TYPE(cp_blacs_env_type), POINTER :: blacs_env 2856 TYPE(cp_cfm_p_type), ALLOCATABLE, DIMENSION(:) :: gamma_00, gamma_inv_00 2857 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: S_mos_virt 2858 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:, :) :: dBerry_mos_occ, gamma_real_imag, opvec 2859 TYPE(cp_fm_struct_type), POINTER :: fm_struct 2860 TYPE(cp_fm_type), POINTER :: ediff_inv, rRc_mos_occ, wfm_ao_ao, & 2861 wfm_mo_virt_mo_occ 2862 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: berry_cossin_xyz, matrix_s, rRc_xyz 2863 TYPE(pw_env_type), POINTER :: pw_env 2864 TYPE(pw_poisson_type), POINTER :: poisson_env 2865 2866 CALL timeset(routineN, handle) 2867 2868 NULLIFY (blacs_env, cell, matrix_s, pw_env) 2869 CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, matrix_s=matrix_s, pw_env=pw_env) 2870 2871 nspins = SIZE(gs_mos) 2872 CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao) 2873 DO ispin = 1, nspins 2874 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ) 2875 nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt) 2876 END DO 2877 2878 ! +++ allocate dipole operator matrices (must be deallocated elsewhere) 2879 ALLOCATE (dipole_op_mos_occ(nderivs, nspins)) 2880 DO ispin = 1, nspins 2881 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct) 2882 2883 DO ideriv = 1, nderivs 2884 NULLIFY (dipole_op_mos_occ(ideriv, ispin)%matrix) 2885 CALL cp_fm_create(dipole_op_mos_occ(ideriv, ispin)%matrix, fm_struct) 2886 END DO 2887 END DO 2888 2889 ! +++ allocate work matrices 2890 ALLOCATE (S_mos_virt(nspins)) 2891 DO ispin = 1, nspins 2892 NULLIFY (S_mos_virt(ispin)%matrix) 2893 CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct) 2894 CALL cp_fm_create(S_mos_virt(ispin)%matrix, fm_struct) 2895 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, & 2896 gs_mos(ispin)%mos_virt, & 2897 S_mos_virt(ispin)%matrix, & 2898 ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp) 2899 END DO 2900 2901 ! check that the chosen dipole operator is consistent with the periodic boundary conditions used 2902 CALL pw_env_get(pw_env, poisson_env=poisson_env) 2903 ndim_periodic = COUNT(poisson_env%parameters%periodic == 1) 2904 2905 SELECT CASE (tddfpt_control%dipole_form) 2906 CASE (tddfpt_dipole_berry) 2907 IF (ndim_periodic /= 3) THEN 2908 CALL cp_warn(__LOCATION__, & 2909 "Fully periodic Poisson solver (PERIODIC xyz) is needed "// & 2910 "for oscillator strengths based on the Berry phase formula") 2911 END IF 2912 2913 NULLIFY (berry_cossin_xyz) 2914 ! index: 1 = Re[exp(-i * G_t * t)], 2915 ! 2 = Im[exp(-i * G_t * t)]; 2916 ! t = x,y,z 2917 CALL dbcsr_allocate_matrix_set(berry_cossin_xyz, 2) 2918 2919 DO i_cos_sin = 1, 2 2920 CALL dbcsr_init_p(berry_cossin_xyz(i_cos_sin)%matrix) 2921 CALL dbcsr_copy(berry_cossin_xyz(i_cos_sin)%matrix, matrix_s(1)%matrix) 2922 END DO 2923 2924 ! +++ allocate berry-phase-related work matrices 2925 ALLOCATE (gamma_00(nspins), gamma_inv_00(nspins), gamma_real_imag(2, nspins), opvec(2, nspins)) 2926 ALLOCATE (dBerry_mos_occ(nderivs, nspins)) 2927 DO ispin = 1, nspins 2928 NULLIFY (fm_struct) 2929 CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_occ(ispin), ncol_global=nmo_occ(ispin), context=blacs_env) 2930 2931 NULLIFY (gamma_00(ispin)%matrix, gamma_inv_00(ispin)%matrix) 2932 CALL cp_cfm_create(gamma_00(ispin)%matrix, fm_struct) 2933 CALL cp_cfm_create(gamma_inv_00(ispin)%matrix, fm_struct) 2934 2935 DO i_cos_sin = 1, 2 2936 NULLIFY (gamma_real_imag(i_cos_sin, ispin)%matrix) 2937 CALL cp_fm_create(gamma_real_imag(i_cos_sin, ispin)%matrix, fm_struct) 2938 END DO 2939 CALL cp_fm_struct_release(fm_struct) 2940 2941 ! G_real C_0, G_imag C_0 2942 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct) 2943 DO i_cos_sin = 1, 2 2944 NULLIFY (opvec(i_cos_sin, ispin)%matrix) 2945 CALL cp_fm_create(opvec(i_cos_sin, ispin)%matrix, fm_struct) 2946 END DO 2947 2948 ! dBerry * C_0 2949 DO ideriv = 1, nderivs 2950 NULLIFY (dBerry_mos_occ(ideriv, ispin)%matrix) 2951 CALL cp_fm_create(dBerry_mos_occ(ideriv, ispin)%matrix, fm_struct) 2952 CALL cp_fm_set_all(dBerry_mos_occ(ideriv, ispin)%matrix, 0.0_dp) 2953 END DO 2954 END DO 2955 2956 DO ideriv = 1, nderivs 2957 kvec(:) = twopi*cell%h_inv(ideriv, :) 2958 DO i_cos_sin = 1, 2 2959 CALL dbcsr_set(berry_cossin_xyz(i_cos_sin)%matrix, 0.0_dp) 2960 END DO 2961 CALL build_berry_moment_matrix(qs_env, berry_cossin_xyz(1)%matrix, berry_cossin_xyz(2)%matrix, kvec) 2962 2963 DO ispin = 1, nspins 2964 ! i_cos_sin = 1: cos (real) component; opvec(1) = gamma_real C_0 2965 ! i_cos_sin = 2: sin (imaginary) component; opvec(2) = gamma_imag C_0 2966 DO i_cos_sin = 1, 2 2967 CALL cp_dbcsr_sm_fm_multiply(berry_cossin_xyz(i_cos_sin)%matrix, & 2968 gs_mos(ispin)%mos_occ, & 2969 opvec(i_cos_sin, ispin)%matrix, & 2970 ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp) 2971 END DO 2972 2973 CALL cp_gemm('T', 'N', nmo_occ(ispin), nmo_occ(ispin), nao, & 2974 1.0_dp, gs_mos(ispin)%mos_occ, opvec(1, ispin)%matrix, & 2975 0.0_dp, gamma_real_imag(1, ispin)%matrix) 2976 2977 CALL cp_gemm('T', 'N', nmo_occ(ispin), nmo_occ(ispin), nao, & 2978 -1.0_dp, gs_mos(ispin)%mos_occ, opvec(2, ispin)%matrix, & 2979 0.0_dp, gamma_real_imag(2, ispin)%matrix) 2980 2981 CALL cp_fm_to_cfm(msourcer=gamma_real_imag(1, ispin)%matrix, & 2982 msourcei=gamma_real_imag(2, ispin)%matrix, & 2983 mtarget=gamma_00(ispin)%matrix) 2984 2985 ! gamma_inv_00 = Q = [C_0^T (gamma_real - i gamma_imag) C_0] ^ {-1} 2986 CALL cp_cfm_set_all(gamma_inv_00(ispin)%matrix, z_zero, z_one) 2987 CALL cp_cfm_solve(gamma_00(ispin)%matrix, gamma_inv_00(ispin)%matrix) 2988 2989 CALL cp_cfm_to_fm(msource=gamma_inv_00(ispin)%matrix, & 2990 mtargetr=gamma_real_imag(1, ispin)%matrix, & 2991 mtargeti=gamma_real_imag(2, ispin)%matrix) 2992 2993 ! dBerry_mos_occ is identical to dBerry_psi0 from qs_linres_op % polar_operators() 2994 CALL cp_gemm("N", "N", nao, nmo_occ(ispin), nmo_occ(ispin), & 2995 1.0_dp, opvec(1, ispin)%matrix, gamma_real_imag(2, ispin)%matrix, & 2996 0.0_dp, dipole_op_mos_occ(1, ispin)%matrix) 2997 CALL cp_gemm("N", "N", nao, nmo_occ(ispin), nmo_occ(ispin), & 2998 -1.0_dp, opvec(2, ispin)%matrix, gamma_real_imag(1, ispin)%matrix, & 2999 1.0_dp, dipole_op_mos_occ(1, ispin)%matrix) 3000 3001 DO jderiv = 1, nderivs 3002 CALL cp_fm_scale_and_add(1.0_dp, dBerry_mos_occ(jderiv, ispin)%matrix, & 3003 cell%hmat(jderiv, ideriv), dipole_op_mos_occ(1, ispin)%matrix) 3004 END DO 3005 END DO 3006 END DO 3007 3008 ! --- release berry-phase-related work matrices 3009 DO ispin = nspins, 1, -1 3010 DO i_cos_sin = SIZE(opvec, 1), 1, -1 3011 CALL cp_fm_release(opvec(i_cos_sin, ispin)%matrix) 3012 END DO 3013 3014 DO i_cos_sin = SIZE(gamma_real_imag, 1), 1, -1 3015 CALL cp_fm_release(gamma_real_imag(i_cos_sin, ispin)%matrix) 3016 END DO 3017 3018 CALL cp_cfm_release(gamma_inv_00(ispin)%matrix) 3019 CALL cp_cfm_release(gamma_00(ispin)%matrix) 3020 END DO 3021 DEALLOCATE (gamma_00, gamma_inv_00, gamma_real_imag, opvec) 3022 CALL dbcsr_deallocate_matrix_set(berry_cossin_xyz) 3023 3024 NULLIFY (fm_struct, wfm_ao_ao) 3025 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env) 3026 CALL cp_fm_create(wfm_ao_ao, fm_struct) 3027 CALL cp_fm_struct_release(fm_struct) 3028 3029 ! trans_dipole = 2|e|/|G_mu| * Tr Imag(evects^T * (gamma_real - i gamma_imag) * C_0 * gamma_inv_00) + 3030 ! 2|e|/|G_mu| * Tr Imag(C_0^T * (gamma_real - i gamma_imag) * evects * gamma_inv_00) , 3031 ! 3032 ! Taking into account the symmetry of the matrices 'gamma_real' and 'gamma_imag' and the fact 3033 ! that the response wave-function is a real-valued function, the above expression can be simplified as 3034 ! trans_dipole = 4|e|/|G_mu| * Tr Imag(evects^T * (gamma_real - i gamma_imag) * C_0 * gamma_inv_00) 3035 ! 3036 ! 1/|G_mu| = |lattice_vector_mu| / (2*pi) . 3037 DO ispin = 1, nspins 3038 ! wfm_ao_ao = S * mos_virt * mos_virt^T 3039 CALL cp_gemm('N', 'T', nao, nao, nmo_virt(ispin), & 3040 1.0_dp/twopi, S_mos_virt(ispin)%matrix, gs_mos(ispin)%mos_virt, & 3041 0.0_dp, wfm_ao_ao) 3042 3043 DO ideriv = 1, nderivs 3044 CALL cp_gemm('N', 'N', nao, nmo_occ(ispin), nao, & 3045 1.0_dp, wfm_ao_ao, dBerry_mos_occ(ideriv, ispin)%matrix, & 3046 0.0_dp, dipole_op_mos_occ(ideriv, ispin)%matrix) 3047 END DO 3048 END DO 3049 3050 CALL cp_fm_release(wfm_ao_ao) 3051 3052 DO ispin = nspins, 1, -1 3053 DO ideriv = SIZE(dBerry_mos_occ, 1), 1, -1 3054 CALL cp_fm_release(dBerry_mos_occ(ideriv, ispin)%matrix) 3055 END DO 3056 END DO 3057 DEALLOCATE (dBerry_mos_occ) 3058 3059 CASE (tddfpt_dipole_length) 3060 IF (ndim_periodic /= 0) THEN 3061 CALL cp_warn(__LOCATION__, & 3062 "Non-periodic Poisson solver (PERIODIC none) is needed "// & 3063 "for oscillator strengths based on the length operator") 3064 END IF 3065 3066 ! compute components of the dipole operator in the length form 3067 NULLIFY (rRc_xyz) 3068 CALL dbcsr_allocate_matrix_set(rRc_xyz, nderivs) 3069 3070 DO ideriv = 1, nderivs 3071 CALL dbcsr_init_p(rRc_xyz(ideriv)%matrix) 3072 CALL dbcsr_copy(rRc_xyz(ideriv)%matrix, matrix_s(1)%matrix) 3073 END DO 3074 3075 CALL get_reference_point(reference_point, qs_env=qs_env, & 3076 reference=tddfpt_control%dipole_reference, & 3077 ref_point=tddfpt_control%dipole_ref_point) 3078 3079 CALL rRc_xyz_ao(op=rRc_xyz, qs_env=qs_env, rc=reference_point, order=1, minimum_image=.FALSE., soft=.FALSE.) 3080 3081 NULLIFY (fm_struct, wfm_ao_ao) 3082 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env) 3083 CALL cp_fm_create(wfm_ao_ao, fm_struct) 3084 CALL cp_fm_struct_release(fm_struct) 3085 3086 DO ispin = 1, nspins 3087 NULLIFY (rRc_mos_occ) 3088 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct) 3089 CALL cp_fm_create(rRc_mos_occ, fm_struct) 3090 3091 ! wfm_ao_ao = S * mos_virt * mos_virt^T 3092 CALL cp_gemm('N', 'T', nao, nao, nmo_virt(ispin), & 3093 1.0_dp, S_mos_virt(ispin)%matrix, gs_mos(ispin)%mos_virt, & 3094 0.0_dp, wfm_ao_ao) 3095 3096 DO ideriv = 1, nderivs 3097 CALL cp_dbcsr_sm_fm_multiply(rRc_xyz(ideriv)%matrix, & 3098 gs_mos(ispin)%mos_occ, & 3099 rRc_mos_occ, & 3100 ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp) 3101 3102 CALL cp_gemm('N', 'N', nao, nmo_occ(ispin), nao, & 3103 1.0_dp, wfm_ao_ao, rRc_mos_occ, & 3104 0.0_dp, dipole_op_mos_occ(ideriv, ispin)%matrix) 3105 END DO 3106 3107 CALL cp_fm_release(rRc_mos_occ) 3108 END DO 3109 3110 CALL cp_fm_release(wfm_ao_ao) 3111 CALL dbcsr_deallocate_matrix_set(rRc_xyz) 3112 3113 CASE (tddfpt_dipole_velocity) 3114 IF (SIZE(matrix_s) < nderivs + 1) THEN 3115 CPABORT("Where is the derivative?") 3116 END IF 3117 3118 DO ispin = 1, nspins 3119 NULLIFY (fm_struct, ediff_inv, wfm_mo_virt_mo_occ) 3120 CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(ispin), ncol_global=nmo_occ(ispin), context=blacs_env) 3121 CALL cp_fm_create(ediff_inv, fm_struct) 3122 CALL cp_fm_create(wfm_mo_virt_mo_occ, fm_struct) 3123 CALL cp_fm_struct_release(fm_struct) 3124 3125 CALL cp_fm_get_info(ediff_inv, nrow_local=nrows_local, ncol_local=ncols_local, & 3126 row_indices=row_indices, col_indices=col_indices, local_data=local_data_ediff) 3127 CALL cp_fm_get_info(wfm_mo_virt_mo_occ, local_data=local_data_wfm) 3128 3129!$OMP PARALLEL DO DEFAULT(NONE), & 3130!$OMP PRIVATE(eval_occ, icol, irow), & 3131!$OMP SHARED(col_indices, gs_mos, ispin, local_data_ediff, ncols_local, nrows_local, row_indices) 3132 DO icol = 1, ncols_local 3133 ! E_occ_i ; imo_occ = col_indices(icol) 3134 eval_occ = gs_mos(ispin)%evals_occ(col_indices(icol)) 3135 3136 DO irow = 1, nrows_local 3137 ! ediff_inv_weights(a, i) = 1.0 / (E_virt_a - E_occ_i) 3138 ! imo_virt = row_indices(irow) 3139 local_data_ediff(irow, icol) = 1.0_dp/(gs_mos(ispin)%evals_virt(row_indices(irow)) - eval_occ) 3140 END DO 3141 END DO 3142!$OMP END PARALLEL DO 3143 3144 DO ideriv = 1, nderivs 3145 CALL cp_dbcsr_sm_fm_multiply(matrix_s(ideriv + 1)%matrix, & 3146 gs_mos(ispin)%mos_occ, & 3147 dipole_op_mos_occ(ideriv, ispin)%matrix, & 3148 ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp) 3149 3150 CALL cp_gemm('T', 'N', nmo_virt(ispin), nmo_occ(ispin), nao, & 3151 1.0_dp, gs_mos(ispin)%mos_virt, dipole_op_mos_occ(ideriv, ispin)%matrix, & 3152 0.0_dp, wfm_mo_virt_mo_occ) 3153 3154 ! in-place element-wise (Schur) product; 3155 ! avoid allocation of a temporary [nmo_virt x nmo_occ] matrix which is needed for cp_fm_schur_product() subroutine call 3156 3157!$OMP PARALLEL DO DEFAULT(NONE), & 3158!$OMP PRIVATE(icol, irow), & 3159!$OMP SHARED(ispin, local_data_ediff, local_data_wfm, ncols_local, nrows_local) 3160 DO icol = 1, ncols_local 3161 DO irow = 1, nrows_local 3162 local_data_wfm(irow, icol) = local_data_wfm(irow, icol)*local_data_ediff(irow, icol) 3163 END DO 3164 END DO 3165!$OMP END PARALLEL DO 3166 3167 CALL cp_gemm('N', 'N', nao, nmo_occ(ispin), nmo_virt(ispin), & 3168 1.0_dp, S_mos_virt(ispin)%matrix, wfm_mo_virt_mo_occ, & 3169 0.0_dp, dipole_op_mos_occ(ideriv, ispin)%matrix) 3170 END DO 3171 3172 CALL cp_fm_release(wfm_mo_virt_mo_occ) 3173 CALL cp_fm_release(ediff_inv) 3174 END DO 3175 3176 CASE DEFAULT 3177 CPABORT("Unimplemented form of the dipole operator") 3178 END SELECT 3179 3180 ! --- release work matrices 3181 DO ispin = nspins, 1, -1 3182 CALL cp_fm_release(S_mos_virt(ispin)%matrix) 3183 END DO 3184 DEALLOCATE (S_mos_virt) 3185 3186 CALL timestop(handle) 3187 END SUBROUTINE tddfpt_dipole_operator 3188 3189! ************************************************************************************************** 3190!> \brief Print final TDDFPT excitation energies and oscillator strengths. 3191!> \param log_unit output unit 3192!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins; 3193!> SIZE(evects,2) -- number of excited states to print) 3194!> \param evals TDDFPT eigenvalues 3195!> \param mult multiplicity 3196!> \param dipole_op_mos_occ action of the dipole operator on the ground state wave function 3197!> [x,y,z ; spin] 3198!> \par History 3199!> * 05.2016 created [Sergey Chulkov] 3200!> * 06.2016 transition dipole moments and oscillator strengths [Sergey Chulkov] 3201!> * 07.2016 spin-unpolarised electron density [Sergey Chulkov] 3202!> * 08.2018 compute 'dipole_op_mos_occ' in a separate subroutine [Sergey Chulkov] 3203!> \note \parblock 3204!> Adapted version of the subroutine find_contributions() which was originally created 3205!> by Thomas Chassaing on 02.2005. 3206!> 3207!> Transition dipole moment along direction 'd' is computed as following: 3208!> \f[ t_d(spin) = Tr[evects^T dipole\_op\_mos\_occ(d, spin)] .\f] 3209!> \endparblock 3210! ************************************************************************************************** 3211 SUBROUTINE tddfpt_print_summary(log_unit, evects, evals, mult, dipole_op_mos_occ) 3212 INTEGER, INTENT(in) :: log_unit 3213 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: evects 3214 REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals 3215 INTEGER, INTENT(in) :: mult 3216 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: dipole_op_mos_occ 3217 3218 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_print_summary', & 3219 routineP = moduleN//':'//routineN 3220 3221 CHARACTER(len=1) :: lsd_str 3222 CHARACTER(len=20) :: mult_str 3223 INTEGER :: handle, ideriv, ispin, istate, nspins, & 3224 nstates 3225 REAL(kind=dp) :: osc_strength 3226 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: trans_dipoles 3227 3228 CALL timeset(routineN, handle) 3229 3230 nspins = SIZE(evects, 1) 3231 nstates = SIZE(evects, 2) 3232 3233 IF (nspins > 1) THEN 3234 lsd_str = 'U' 3235 ELSE 3236 lsd_str = 'R' 3237 END IF 3238 3239 ! *** summary header *** 3240 IF (log_unit > 0) THEN 3241 CALL integer_to_string(mult, mult_str) 3242 WRITE (log_unit, '(/,1X,A1,A,1X,A,/)') lsd_str, "-TDDFPT states of multiplicity", TRIM(mult_str) 3243 3244 WRITE (log_unit, '(T10,A,T19,A,T37,A,T69,A)') "State", "Excitation", "Transition dipole (a.u.)", "Oscillator" 3245 WRITE (log_unit, '(T10,A,T19,A,T37,A,T49,A,T61,A,T67,A)') "number", "energy (eV)", "x", "y", "z", "strength (a.u.)" 3246 WRITE (log_unit, '(T10,72("-"))') 3247 END IF 3248 3249 ! transition dipole moment 3250 ALLOCATE (trans_dipoles(nstates, nderivs, nspins)) 3251 trans_dipoles(:, :, :) = 0.0_dp 3252 3253 ! nspins == 1 .AND. mult == 3 : spin-flip transitions are forbidden due to symmetry reasons 3254 IF (nspins > 1 .OR. mult == 1) THEN 3255 DO ispin = 1, nspins 3256 DO ideriv = 1, nderivs 3257 CALL cp_fm_trace(evects(ispin, :), dipole_op_mos_occ(ideriv, ispin)%matrix, trans_dipoles(:, ideriv, ispin)) 3258 END DO 3259 END DO 3260 3261 IF (nspins == 1) THEN 3262 trans_dipoles(:, :, 1) = SQRT(2.0_dp)*trans_dipoles(:, :, 1) 3263 ELSE 3264 trans_dipoles(:, :, 1) = SQRT(trans_dipoles(:, :, 1)**2 + trans_dipoles(:, :, 2)**2) 3265 END IF 3266 END IF 3267 3268 ! *** summary information *** 3269 DO istate = 1, nstates 3270 IF (log_unit > 0) THEN 3271 osc_strength = 2.0_dp/3.0_dp*evals(istate)* & 3272 accurate_dot_product(trans_dipoles(istate, :, 1), trans_dipoles(istate, :, 1)) 3273 3274 WRITE (log_unit, '(1X,A,T9,I7,T19,F11.5,T31,3(1X,ES11.4E2),T69,ES12.5E2)') & 3275 "TDDFPT|", istate, evals(istate)*evolt, trans_dipoles(istate, 1:nderivs, 1), osc_strength 3276 END IF 3277 END DO 3278 3279 DEALLOCATE (trans_dipoles) 3280 3281 CALL timestop(handle) 3282 END SUBROUTINE tddfpt_print_summary 3283 3284! ************************************************************************************************** 3285!> \brief Print excitation analysis. 3286!> \param log_unit output unit 3287!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins; 3288!> SIZE(evects,2) -- number of excited states to print) 3289!> \param gs_mos molecular orbitals optimised for the ground state 3290!> \param matrix_s overlap matrix 3291!> \param min_amplitude the smallest excitation amplitude to print 3292!> \par History 3293!> * 05.2016 created as 'tddfpt_print_summary' [Sergey Chulkov] 3294!> * 08.2018 splited of from 'tddfpt_print_summary' [Sergey Chulkov] 3295! ************************************************************************************************** 3296 SUBROUTINE tddfpt_print_excitation_analysis(log_unit, evects, gs_mos, matrix_s, min_amplitude) 3297 INTEGER, INTENT(in) :: log_unit 3298 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: evects 3299 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 3300 INTENT(in) :: gs_mos 3301 TYPE(dbcsr_type), POINTER :: matrix_s 3302 REAL(kind=dp), INTENT(in) :: min_amplitude 3303 3304 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_print_excitation_analysis', & 3305 routineP = moduleN//':'//routineN 3306 3307 CHARACTER(len=5) :: spin_label 3308 INTEGER :: handle, icol, iproc, irow, ispin, istate, nao, ncols_local, nrows_local, nspins, & 3309 nstates, send_handler, send_handler2, state_spin 3310 INTEGER(kind=int_8) :: iexc, imo_occ, imo_virt, ind, nexcs, & 3311 nexcs_local, nexcs_max_local, & 3312 nmo_virt_occ, nmo_virt_occ_alpha 3313 INTEGER(kind=int_8), ALLOCATABLE, DIMENSION(:) :: inds_local, inds_recv, nexcs_recv 3314 INTEGER(kind=int_8), DIMENSION(1) :: nexcs_send 3315 INTEGER(kind=int_8), DIMENSION(maxspins) :: nmo_occ8, nmo_virt8 3316 INTEGER, ALLOCATABLE, DIMENSION(:) :: inds, recv_handlers, recv_handlers2 3317 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 3318 INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_virt 3319 LOGICAL :: do_exc_analysis 3320 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: weights_local, weights_neg_abs_recv, & 3321 weights_recv 3322 REAL(kind=dp), DIMENSION(:, :), POINTER :: local_data 3323 TYPE(cp_blacs_env_type), POINTER :: blacs_env 3324 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: S_mos_virt, weights_fm 3325 TYPE(cp_fm_struct_type), POINTER :: fm_struct 3326 TYPE(cp_para_env_type), POINTER :: para_env 3327 3328 CALL timeset(routineN, handle) 3329 3330 nspins = SIZE(evects, 1) 3331 nstates = SIZE(evects, 2) 3332 do_exc_analysis = min_amplitude < 1.0_dp 3333 3334 CALL cp_fm_get_info(gs_mos(1)%mos_occ, context=blacs_env, para_env=para_env) 3335 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao) 3336 3337 DO ispin = 1, nspins 3338 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ) 3339 nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt) 3340 nmo_occ8(ispin) = SIZE(gs_mos(ispin)%evals_occ, kind=int_8) 3341 nmo_virt8(ispin) = SIZE(gs_mos(ispin)%evals_virt, kind=int_8) 3342 END DO 3343 3344 ! *** excitation analysis *** 3345 IF (do_exc_analysis) THEN 3346 CPASSERT(log_unit <= 0 .OR. para_env%ionode) 3347 nmo_virt_occ_alpha = INT(nmo_virt(1), int_8)*INT(nmo_occ(1), int_8) 3348 3349 IF (log_unit > 0) THEN 3350 WRITE (log_unit, '(/,1X,A,/)') "Excitation analysis" 3351 3352 WRITE (log_unit, '(3X,A,T17,A,T34,A,T49,A)') "State", "Occupied", "Virtual", "Excitation" 3353 WRITE (log_unit, '(3X,A,T18,A,T34,A,T49,A)') "number", "orbital", "orbital", "amplitude" 3354 WRITE (log_unit, '(1X,57("-"))') 3355 3356 IF (nspins == 1) THEN 3357 state_spin = 1 3358 spin_label = ' ' 3359 END IF 3360 END IF 3361 3362 ALLOCATE (S_mos_virt(nspins), weights_fm(nspins)) 3363 DO ispin = 1, nspins 3364 NULLIFY (S_mos_virt(ispin)%matrix) 3365 CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct) 3366 CALL cp_fm_create(S_mos_virt(ispin)%matrix, fm_struct) 3367 CALL cp_dbcsr_sm_fm_multiply(matrix_s, & 3368 gs_mos(ispin)%mos_virt, & 3369 S_mos_virt(ispin)%matrix, & 3370 ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp) 3371 3372 NULLIFY (fm_struct, weights_fm(ispin)%matrix) 3373 CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(ispin), ncol_global=nmo_occ(ispin), context=blacs_env) 3374 CALL cp_fm_create(weights_fm(ispin)%matrix, fm_struct) 3375 CALL cp_fm_struct_release(fm_struct) 3376 END DO 3377 3378 nexcs_max_local = 0 3379 DO ispin = 1, nspins 3380 CALL cp_fm_get_info(weights_fm(ispin)%matrix, nrow_local=nrows_local, ncol_local=ncols_local) 3381 nexcs_max_local = nexcs_max_local + INT(nrows_local, int_8)*INT(ncols_local, int_8) 3382 END DO 3383 3384 ALLOCATE (weights_local(nexcs_max_local), inds_local(nexcs_max_local)) 3385 3386 DO istate = 1, nstates 3387 nexcs_local = 0 3388 nmo_virt_occ = 0 3389 3390 ! analyse matrix elements locally and transfer only significant 3391 ! excitations to the master node for subsequent ordering 3392 DO ispin = 1, nspins 3393 ! compute excitation amplitudes 3394 CALL cp_gemm('T', 'N', nmo_virt(ispin), nmo_occ(ispin), nao, 1.0_dp, S_mos_virt(ispin)%matrix, & 3395 evects(ispin, istate)%matrix, 0.0_dp, weights_fm(ispin)%matrix) 3396 3397 CALL cp_fm_get_info(weights_fm(ispin)%matrix, nrow_local=nrows_local, ncol_local=ncols_local, & 3398 row_indices=row_indices, col_indices=col_indices, local_data=local_data) 3399 3400 ! locate single excitations with significant amplitudes (>= min_amplitude) 3401 DO icol = 1, ncols_local 3402 DO irow = 1, nrows_local 3403 IF (ABS(local_data(irow, icol)) >= min_amplitude) THEN 3404 ! number of non-negligible excitations 3405 nexcs_local = nexcs_local + 1 3406 ! excitation amplitude 3407 weights_local(nexcs_local) = local_data(irow, icol) 3408 ! index of single excitation (ivirt, iocc, ispin) in compressed form 3409 inds_local(nexcs_local) = nmo_virt_occ + INT(row_indices(irow), int_8) + & 3410 INT(col_indices(icol) - 1, int_8)*nmo_virt8(ispin) 3411 END IF 3412 END DO 3413 END DO 3414 3415 nmo_virt_occ = nmo_virt_occ + nmo_virt8(ispin)*nmo_occ8(ispin) 3416 END DO 3417 3418 IF (para_env%ionode) THEN 3419 ! master node 3420 ALLOCATE (nexcs_recv(para_env%num_pe), recv_handlers(para_env%num_pe), recv_handlers2(para_env%num_pe)) 3421 3422 ! collect number of non-negligible excitations from other nodes 3423 DO iproc = 1, para_env%num_pe 3424 IF (iproc - 1 /= para_env%mepos) THEN 3425 CALL mp_irecv(nexcs_recv(iproc:iproc), iproc - 1, para_env%group, recv_handlers(iproc), 0) 3426 ELSE 3427 nexcs_recv(iproc) = nexcs_local 3428 END IF 3429 END DO 3430 3431 DO iproc = 1, para_env%num_pe 3432 IF (iproc - 1 /= para_env%mepos) & 3433 CALL mp_wait(recv_handlers(iproc)) 3434 END DO 3435 3436 ! compute total number of non-negligible excitations 3437 nexcs = 0 3438 DO iproc = 1, para_env%num_pe 3439 nexcs = nexcs + nexcs_recv(iproc) 3440 END DO 3441 3442 ! receive indices and amplitudes of selected excitations 3443 ALLOCATE (weights_recv(nexcs), weights_neg_abs_recv(nexcs)) 3444 ALLOCATE (inds_recv(nexcs), inds(nexcs)) 3445 3446 nmo_virt_occ = 0 3447 DO iproc = 1, para_env%num_pe 3448 IF (nexcs_recv(iproc) > 0) THEN 3449 IF (iproc - 1 /= para_env%mepos) THEN 3450 ! excitation amplitudes 3451 CALL mp_irecv(weights_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)), & 3452 iproc - 1, para_env%group, recv_handlers(iproc), 1) 3453 ! compressed indices 3454 CALL mp_irecv(inds_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)), & 3455 iproc - 1, para_env%group, recv_handlers2(iproc), 2) 3456 ELSE 3457 ! data on master node 3458 weights_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)) = weights_local(1:nexcs_recv(iproc)) 3459 inds_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)) = inds_local(1:nexcs_recv(iproc)) 3460 END IF 3461 3462 nmo_virt_occ = nmo_virt_occ + nexcs_recv(iproc) 3463 END IF 3464 END DO 3465 3466 DO iproc = 1, para_env%num_pe 3467 IF (iproc - 1 /= para_env%mepos .AND. nexcs_recv(iproc) > 0) THEN 3468 CALL mp_wait(recv_handlers(iproc)) 3469 CALL mp_wait(recv_handlers2(iproc)) 3470 END IF 3471 END DO 3472 3473 DEALLOCATE (nexcs_recv, recv_handlers, recv_handlers2) 3474 ELSE 3475 ! working node: send the number of selected excited states to the master node 3476 nexcs_send(1) = nexcs_local 3477 CALL mp_isend(nexcs_send, para_env%source, para_env%group, send_handler, 0) 3478 CALL mp_wait(send_handler) 3479 3480 IF (nexcs_local > 0) THEN 3481 ! send excitation amplitudes 3482 CALL mp_isend(weights_local(1:nexcs_local), para_env%source, para_env%group, send_handler, 1) 3483 ! send compressed indices 3484 CALL mp_isend(inds_local(1:nexcs_local), para_env%source, para_env%group, send_handler2, 2) 3485 3486 CALL mp_wait(send_handler) 3487 CALL mp_wait(send_handler2) 3488 END IF 3489 END IF 3490 3491 ! sort non-negligible excitations on the master node according to their amplitudes, 3492 ! uncompress indices and print summary information 3493 IF (para_env%ionode .AND. log_unit > 0) THEN 3494 weights_neg_abs_recv(:) = -ABS(weights_recv) 3495 CALL sort(weights_neg_abs_recv, INT(nexcs), inds) 3496 3497 WRITE (log_unit, '(1X,I8)') istate 3498 3499 DO iexc = 1, nexcs 3500 ind = inds_recv(inds(iexc)) - 1 3501 IF (nspins > 1) THEN 3502 IF (ind < nmo_virt_occ_alpha) THEN 3503 state_spin = 1 3504 spin_label = '(alp)' 3505 ELSE 3506 state_spin = 2 3507 ind = ind - nmo_virt_occ_alpha 3508 spin_label = '(bet)' 3509 END IF 3510 END IF 3511 3512 imo_occ = ind/nmo_virt8(state_spin) + 1 3513 imo_virt = MOD(ind, nmo_virt8(state_spin)) + 1 3514 3515 WRITE (log_unit, '(T14,I8,1X,A5,T30,I8,1X,A5,T50,F9.6)') imo_occ, spin_label, & 3516 nmo_occ8(state_spin) + imo_virt, spin_label, weights_recv(inds(iexc)) 3517 END DO 3518 END IF 3519 3520 ! deallocate temporary arrays 3521 IF (para_env%ionode) & 3522 DEALLOCATE (weights_recv, weights_neg_abs_recv, inds_recv, inds) 3523 END DO 3524 3525 DEALLOCATE (weights_local, inds_local) 3526 END IF 3527 3528 DO ispin = nspins, 1, -1 3529 CALL cp_fm_release(weights_fm(ispin)%matrix) 3530 CALL cp_fm_release(S_mos_virt(ispin)%matrix) 3531 END DO 3532 DEALLOCATE (S_mos_virt, weights_fm) 3533 3534 CALL timestop(handle) 3535 END SUBROUTINE tddfpt_print_excitation_analysis 3536 3537! ************************************************************************************************** 3538!> \brief Write Ritz vectors to a binary restart file. 3539!> \param evects vectors to store 3540!> \param evals TDDFPT eigenvalues 3541!> \param gs_mos structure that holds ground state occupied and virtual 3542!> molecular orbitals 3543!> \param logger a logger object 3544!> \param tddfpt_print_section TDDFPT%PRINT input section 3545!> \par History 3546!> * 08.2016 created [Sergey Chulkov] 3547! ************************************************************************************************** 3548 SUBROUTINE tddfpt_write_restart(evects, evals, gs_mos, logger, tddfpt_print_section) 3549 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: evects 3550 REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals 3551 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 3552 INTENT(in) :: gs_mos 3553 TYPE(cp_logger_type), POINTER :: logger 3554 TYPE(section_vals_type), POINTER :: tddfpt_print_section 3555 3556 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_write_restart', & 3557 routineP = moduleN//':'//routineN 3558 3559 INTEGER :: handle, ispin, istate, nao, nspins, & 3560 nstates, ounit 3561 INTEGER, DIMENSION(maxspins) :: nmo_occ 3562 3563 IF (BTEST(cp_print_key_should_output(logger%iter_info, tddfpt_print_section, "RESTART"), cp_p_file)) THEN 3564 CALL timeset(routineN, handle) 3565 3566 nspins = SIZE(evects, 1) 3567 nstates = SIZE(evects, 2) 3568 3569 IF (debug_this_module) THEN 3570 CPASSERT(SIZE(evals) == nstates) 3571 CPASSERT(nspins > 0) 3572 CPASSERT(nstates > 0) 3573 END IF 3574 3575 CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao) 3576 3577 DO ispin = 1, nspins 3578 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ) 3579 END DO 3580 3581 ounit = cp_print_key_unit_nr(logger, tddfpt_print_section, "RESTART", & 3582 extension=".tdwfn", file_status="REPLACE", file_action="WRITE", & 3583 do_backup=.TRUE., file_form="UNFORMATTED") 3584 3585 IF (ounit > 0) THEN 3586 WRITE (ounit) nstates, nspins, nao 3587 WRITE (ounit) nmo_occ(1:nspins) 3588 WRITE (ounit) evals 3589 END IF 3590 3591 DO istate = 1, nstates 3592 DO ispin = 1, nspins 3593 ! TDDFPT wave function is actually stored as a linear combination of virtual MOs 3594 ! that replaces the corresponding deoccupied MO. Unfortunately, the phase 3595 ! of the occupied MOs varies depending on the eigensolver used as well as 3596 ! how eigenvectors are distributed across computational cores. The phase is important 3597 ! because TDDFPT wave functions are used to compute a response electron density 3598 ! \rho^{-} = 1/2 * [C_{0} * evect^T + evect * C_{0}^{-}], where C_{0} is the expansion 3599 ! coefficients of the reference ground-state wave function. To make the restart file 3600 ! transferable, TDDFPT wave functions are stored in assumption that all ground state 3601 ! MOs have a positive phase. 3602 CALL cp_fm_column_scale(evects(ispin, istate)%matrix, gs_mos(ispin)%phases_occ) 3603 3604 CALL cp_fm_write_unformatted(evects(ispin, istate)%matrix, ounit) 3605 3606 CALL cp_fm_column_scale(evects(ispin, istate)%matrix, gs_mos(ispin)%phases_occ) 3607 END DO 3608 END DO 3609 3610 CALL cp_print_key_finished_output(ounit, logger, tddfpt_print_section, "RESTART") 3611 3612 CALL timestop(handle) 3613 END IF 3614 END SUBROUTINE tddfpt_write_restart 3615 3616! ************************************************************************************************** 3617!> \brief Initialise initial guess vectors by reading (un-normalised) Ritz vectors 3618!> from a binary restart file. 3619!> \param evects vectors to initialise (initialised on exit) 3620!> \param evals TDDFPT eigenvalues (initialised on exit) 3621!> \param gs_mos structure that holds ground state occupied and virtual 3622!> molecular orbitals 3623!> \param logger a logger object 3624!> \param tddfpt_section TDDFPT input section 3625!> \param tddfpt_print_section TDDFPT%PRINT input section 3626!> \param fm_pool_ao_mo_occ pools of dense matrices with shape [nao x nmo_occ(spin)] 3627!> \param blacs_env_global BLACS parallel environment involving all the processor 3628!> \return the number of excited states found in the restart file 3629!> \par History 3630!> * 08.2016 created [Sergey Chulkov] 3631! ************************************************************************************************** 3632 FUNCTION tddfpt_read_restart(evects, evals, gs_mos, logger, tddfpt_section, tddfpt_print_section, & 3633 fm_pool_ao_mo_occ, blacs_env_global) RESULT(nstates_read) 3634 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: evects 3635 REAL(kind=dp), DIMENSION(:), INTENT(out) :: evals 3636 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 3637 INTENT(in) :: gs_mos 3638 TYPE(cp_logger_type), POINTER :: logger 3639 TYPE(section_vals_type), POINTER :: tddfpt_section, tddfpt_print_section 3640 TYPE(cp_fm_pool_p_type), DIMENSION(:), INTENT(in) :: fm_pool_ao_mo_occ 3641 TYPE(cp_blacs_env_type), POINTER :: blacs_env_global 3642 INTEGER :: nstates_read 3643 3644 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_read_restart', & 3645 routineP = moduleN//':'//routineN 3646 3647 CHARACTER(len=20) :: read_str, ref_str 3648 CHARACTER(LEN=default_path_length) :: filename 3649 INTEGER :: handle, ispin, istate, iunit, n_rep_val, & 3650 nao, nao_read, nspins, nspins_read, & 3651 nstates 3652 INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_occ_read 3653 LOGICAL :: file_exists 3654 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals_read 3655 TYPE(cp_para_env_type), POINTER :: para_env_global 3656 TYPE(section_vals_type), POINTER :: print_key 3657 3658 CALL timeset(routineN, handle) 3659 3660 ! generate restart file name 3661 CALL section_vals_val_get(tddfpt_section, "WFN_RESTART_FILE_NAME", n_rep_val=n_rep_val) 3662 IF (n_rep_val > 0) THEN 3663 CALL section_vals_val_get(tddfpt_section, "WFN_RESTART_FILE_NAME", c_val=filename) 3664 ELSE 3665 print_key => section_vals_get_subs_vals(tddfpt_print_section, "RESTART") 3666 filename = cp_print_key_generate_filename(logger, print_key, & 3667 extension=".tdwfn", my_local=.FALSE.) 3668 END IF 3669 3670 CALL get_blacs_info(blacs_env_global, para_env=para_env_global) 3671 3672 IF (para_env_global%ionode) THEN 3673 INQUIRE (FILE=filename, exist=file_exists) 3674 3675 IF (.NOT. file_exists) THEN 3676 nstates_read = 0 3677 CALL mp_bcast(nstates_read, para_env_global%source, para_env_global%group) 3678 3679 CALL cp_warn(__LOCATION__, & 3680 "User requested to restart the TDDFPT wave functions from the file '"//TRIM(filename)// & 3681 "' which does not exist. Guess wave functions will be constructed using Kohn-Sham orbitals.") 3682 CALL timestop(handle) 3683 RETURN 3684 END IF 3685 3686 CALL open_file(file_name=filename, file_action="READ", file_form="UNFORMATTED", file_status="OLD", unit_number=iunit) 3687 END IF 3688 3689 nspins = SIZE(evects, 1) 3690 nstates = SIZE(evects, 2) 3691 CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao) 3692 3693 DO ispin = 1, nspins 3694 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ) 3695 END DO 3696 3697 IF (para_env_global%ionode) THEN 3698 READ (iunit) nstates_read, nspins_read, nao_read 3699 3700 IF (nspins_read /= nspins) THEN 3701 CALL integer_to_string(nspins, ref_str) 3702 CALL integer_to_string(nspins_read, read_str) 3703 CALL cp_abort(__LOCATION__, & 3704 "Restarted TDDFPT wave function contains incompatible number of spin components ("// & 3705 TRIM(read_str)//" instead of "//TRIM(ref_str)//").") 3706 END IF 3707 3708 IF (nao_read /= nao) THEN 3709 CALL integer_to_string(nao, ref_str) 3710 CALL integer_to_string(nao_read, read_str) 3711 CALL cp_abort(__LOCATION__, & 3712 "Incompatible number of atomic orbitals ("//TRIM(read_str)//" instead of "//TRIM(ref_str)//").") 3713 END IF 3714 3715 READ (iunit) nmo_occ_read(1:nspins) 3716 3717 DO ispin = 1, nspins 3718 IF (nmo_occ_read(ispin) /= nmo_occ(ispin)) THEN 3719 CALL cp_abort(__LOCATION__, & 3720 "Incompatible number of electrons and/or multiplicity.") 3721 END IF 3722 END DO 3723 3724 IF (nstates_read /= nstates) THEN 3725 CALL integer_to_string(nstates, ref_str) 3726 CALL integer_to_string(nstates_read, read_str) 3727 CALL cp_warn(__LOCATION__, & 3728 "TDDFPT restart file contains "//TRIM(read_str)// & 3729 " wave function(s) however "//TRIM(ref_str)// & 3730 " excited states were requested.") 3731 END IF 3732 END IF 3733 CALL mp_bcast(nstates_read, para_env_global%source, para_env_global%group) 3734 3735 ! exit if restart file does not exist 3736 IF (nstates_read <= 0) THEN 3737 CALL timestop(handle) 3738 RETURN 3739 END IF 3740 3741 IF (para_env_global%ionode) THEN 3742 ALLOCATE (evals_read(nstates_read)) 3743 READ (iunit) evals_read 3744 IF (nstates_read <= nstates) THEN 3745 evals(1:nstates_read) = evals_read(1:nstates_read) 3746 ELSE 3747 evals(1:nstates) = evals_read(1:nstates) 3748 END IF 3749 DEALLOCATE (evals_read) 3750 END IF 3751 CALL mp_bcast(evals, para_env_global%source, para_env_global%group) 3752 3753 DO istate = 1, nstates_read 3754 DO ispin = 1, nspins 3755 IF (istate <= nstates) THEN 3756 CALL fm_pool_create_fm(fm_pool_ao_mo_occ(ispin)%pool, evects(ispin, istate)%matrix) 3757 3758 CALL cp_fm_read_unformatted(evects(ispin, istate)%matrix, iunit) 3759 3760 CALL cp_fm_column_scale(evects(ispin, istate)%matrix, gs_mos(ispin)%phases_occ) 3761 END IF 3762 END DO 3763 END DO 3764 3765 IF (para_env_global%ionode) & 3766 CALL close_file(unit_number=iunit) 3767 3768 CALL timestop(handle) 3769 END FUNCTION tddfpt_read_restart 3770END MODULE qs_tddfpt2_methods 3771